summaryrefslogtreecommitdiff
path: root/kaldi_io/src/kaldi/matrix/optimization.h
blob: 66309acaad539b9fbc8575a7e479f02ce2a0ff27 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
// matrix/optimization.h

// Copyright 2012  Johns Hopkins University (author: Daniel Povey)
//
// See ../../COPYING for clarification regarding multiple authors
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//  http://www.apache.org/licenses/LICENSE-2.0
//
// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
// MERCHANTABLITY OR NON-INFRINGEMENT.
// See the Apache 2 License for the specific language governing permissions and
// limitations under the License.
//
// (*) incorporates, with permission, FFT code from his book
// "Signal Processing with Lapped Transforms", Artech, 1992.



#ifndef KALDI_MATRIX_OPTIMIZATION_H_
#define KALDI_MATRIX_OPTIMIZATION_H_

#include "matrix/kaldi-vector.h"
#include "matrix/kaldi-matrix.h"

namespace kaldi {


/// @addtogroup matrix_optimization
/// @{

struct LinearCgdOptions {
  int32 max_iters;  //  Maximum number of iters (if >= 0).
  BaseFloat max_error;  // Maximum 2-norm of the residual A x - b (convergence
                        // test)
  // Every time the residual 2-norm decreases by this recompute_residual_factor
  // since the last time it was computed from scratch, recompute it from
  // scratch.  This helps to keep the computed residual accurate even in the
  // presence of roundoff.
  BaseFloat recompute_residual_factor;
  
  LinearCgdOptions(): max_iters(-1),
                      max_error(0.0),
                      recompute_residual_factor(0.01) { }
};
  
/*
  This function uses linear conjugate gradient descent to approximately solve
  the system A x = b.  The value of x at entry corresponds to the initial guess
  of x.  The algorithm continues until the number of iterations equals b.Dim(),
  or until the 2-norm of (A x - b) is <= max_error, or until the number of
  iterations equals max_iter, whichever happens sooner.  It is a requirement
  that A be positive definite.
  It returns the number of iterations that were actually executed (this is
  useful for testing purposes).
*/
template<typename Real>
int32 LinearCgd(const LinearCgdOptions &opts,
                const SpMatrix<Real> &A, const VectorBase<Real> &b,
                VectorBase<Real> *x);






/**
   This is an implementation of L-BFGS.  It pushes responsibility for
   determining when to stop, onto the user.  There is no call-back here:
   everything is done via calls to the class itself (see the example in
   matrix-lib-test.cc).  This does not implement constrained L-BFGS, but it will
   handle constrained problems correctly as long as the function approaches
   +infinity (or -infinity for maximization problems) when it gets close to the
   bound of the constraint.  In these types of problems, you just let the
   function value be +infinity for minimization problems, or -infinity for
   maximization problems, outside these bounds).
*/

struct LbfgsOptions {
  bool minimize; // if true, we're minimizing, else maximizing.
  int m; // m is the number of stored vectors L-BFGS keeps.
  float first_step_learning_rate; // The very first step of L-BFGS is
  // like gradient descent.  If you want to configure the size of that step,
  // you can do it using this variable.
  float first_step_length; // If this variable is >0.0, it overrides
  // first_step_learning_rate; on the first step we choose an approximate
  // Hessian that is the multiple of the identity that would generate this
  // step-length, or 1.0 if the gradient is zero.
  float first_step_impr; // If this variable is >0.0, it overrides
  // first_step_learning_rate; on the first step we choose an approximate
  // Hessian that is the multiple of the identity that would generate this
  // amount of objective function improvement (assuming the "real" objf
  // was linear).
  float c1; // A constant in Armijo rule = Wolfe condition i)
  float c2; // A constant in Wolfe condition ii)
  float d; // An amount > 1.0 (default 2.0) that we initially multiply or
  // divide the step length by, in the line search.
  int max_line_search_iters; // after this many iters we restart L-BFGS.
  int avg_step_length; // number of iters to avg step length over, in
  // RecentStepLength().
  
  LbfgsOptions (bool minimize = true):
      minimize(minimize),
      m(10),
      first_step_learning_rate(1.0),
      first_step_length(0.0),
      first_step_impr(0.0),
      c1(1.0e-04),
      c2(0.9),
      d(2.0),
      max_line_search_iters(50),
      avg_step_length(4) { }
};

template<typename Real>
class OptimizeLbfgs {
 public:
  /// Initializer takes the starting value of x.
  OptimizeLbfgs(const VectorBase<Real> &x,
                const LbfgsOptions &opts);
  
  /// This returns the value of the variable x that has the best objective
  /// function so far, and the corresponding objective function value if
  /// requested.  This would typically be called only at the end.
  const VectorBase<Real>& GetValue(Real *objf_value = NULL) const;
  
  /// This returns the value at which the function wants us
  /// to compute the objective function and gradient.
  const VectorBase<Real>& GetProposedValue() const { return new_x_; }
  
  /// Returns the average magnitude of the last n steps (but not
  /// more than the number we have stored).  Before we have taken
  /// any steps, returns +infinity.  Note: if the most recent
  /// step length was 0, it returns 0, regardless of the other
  /// step lengths.  This makes it suitable as a convergence test
  /// (else we'd generate NaN's).
  Real RecentStepLength() const;
  
  /// The user calls this function to provide the class with the
  /// function and gradient info at the point GetProposedValue().
  /// If this point is outside the constraints you can set function_value
  /// to {+infinity,-infinity} for {minimization,maximization} problems.
  /// In this case the gradient, and also the second derivative (if you call
  /// the second overloaded version of this function) will be ignored.
  void DoStep(Real function_value,
              const VectorBase<Real> &gradient);
  
  /// The user can call this version of DoStep() if it is desired to set some
  /// kind of approximate Hessian on this iteration.  Note: it is a prerequisite
  /// that diag_approx_2nd_deriv must be strictly positive (minimizing), or
  /// negative (maximizing).
  void DoStep(Real function_value,
              const VectorBase<Real> &gradient,
              const VectorBase<Real> &diag_approx_2nd_deriv);
  
 private:
  KALDI_DISALLOW_COPY_AND_ASSIGN(OptimizeLbfgs);


  // The following variable says what stage of the computation we're at.
  // Refer to Algorithm 7.5 (L-BFGS) of Nodecdal & Wright, "Numerical
  // Optimization", 2nd edition.
  // kBeforeStep means we're about to do
  /// "compute p_k <-- - H_k \delta f_k" (i.e. Algorithm 7.4).
  // kWithinStep means we're at some point within line search; note
  // that line search is iterative so we can stay in this state more
  // than one time on each iteration.
  enum ComputationState {
    kBeforeStep,
    kWithinStep, // This means we're within the step-size computation, and
    // have not yet done the 1st function evaluation.
  };
  
  inline MatrixIndexT Dim() { return x_.Dim(); }
  inline MatrixIndexT M() { return opts_.m; }
  SubVector<Real> Y(MatrixIndexT i) {
    return SubVector<Real>(data_, (i % M()) * 2); // vector y_i
  }
  SubVector<Real> S(MatrixIndexT i) {
    return SubVector<Real>(data_, (i % M()) * 2 + 1); // vector s_i
  }
  // The following are subroutines within DoStep():
  bool AcceptStep(Real function_value,
                  const VectorBase<Real> &gradient);
  void Restart(const VectorBase<Real> &x,
               Real function_value,
               const VectorBase<Real> &gradient);
  void ComputeNewDirection(Real function_value,
                           const VectorBase<Real> &gradient);
  void ComputeHifNeeded(const VectorBase<Real> &gradient);
  void StepSizeIteration(Real function_value,
                         const VectorBase<Real> &gradient);
  void RecordStepLength(Real s);
  
  
  LbfgsOptions opts_;
  SignedMatrixIndexT k_; // Iteration number, starts from zero.  Gets set back to zero
  // when we restart.
  
  ComputationState computation_state_;
  bool H_was_set_; // True if the user specified H_; if false,
  // we'll use a heuristic to estimate it.


  Vector<Real> x_; // current x.
  Vector<Real> new_x_; // the x proposed in the line search.
  Vector<Real> best_x_; // the x with the best objective function so far
                        // (either the same as x_ or something in the current line search.)
  Vector<Real> deriv_; // The most recently evaluated derivative-- at x_k.
  Vector<Real> temp_;
  Real f_; // The function evaluated at x_k.
  Real best_f_; // the best objective function so far.
  Real d_; // a number d > 1.0, but during an iteration we may decrease this, when
  // we switch between armijo and wolfe failures.

  int num_wolfe_i_failures_; // the num times we decreased step size.
  int num_wolfe_ii_failures_; // the num times we increased step size.
  enum { kWolfeI, kWolfeII, kNone } last_failure_type_; // last type of step-search
  // failure on this iter.
  
  Vector<Real> H_; // Current inverse-Hessian estimate.  May be computed by this class itself,
  // or provided by user using 2nd form of SetGradientInfo().
  Matrix<Real> data_; // dimension (m*2) x dim.  Even rows store
  // gradients y_i, odd rows store steps s_i.
  Vector<Real> rho_; // dimension m; rho_(m) = 1/(y_m^T s_m), Eq. 7.17.

  std::vector<Real> step_lengths_; // The step sizes we took on the last
  // (up to m) iterations; these are not stored in a rotating buffer but
  // are shifted by one each time (this is more convenient when we
  // restart, as we keep this info past restarting).
  

};
  
/// @} 


} // end namespace kaldi



#endif