12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fixed maxIteration bug

release/4.3a0
Frank Dellaert 2010-02-14 07:14:42 +00:00
parent 1f165a9f85
commit 44094b494e
1 changed files with 16 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
cpp/iterative-inl.h
View File

@ -31,7 +31,7 @@ namespace gtsam {
k = 0; k = 0;
verbose = verb; verbose = verb;
steepest = steep; steepest = steep;
maxIterations == maxIt ? maxIt : dim(x) * (steepest ? 10 : 1); maxIterations = (maxIt > 0) ? maxIt : dim(x) * (steepest ? 10 : 1);
reset = (size_t) (sqrt(dim(x)) + 0.5); // when to reset reset = (size_t) (sqrt(dim(x)) + 0.5); // when to reset
// Start with g0 = A'*(A*x0-b), d0 = - g0 // Start with g0 = A'*(A*x0-b), d0 = - g0
@ -48,16 +48,18 @@ namespace gtsam {
} }
/** print */ /** print */
void print() { void print(const V& x) {
cout << "iteration = " << k << endl; cout << "iteration = " << k << endl;
cout << "dotg = " << gamma << endl; gtsam::print(x,"x");
gtsam::print(g, "g"); gtsam::print(g, "g");
cout << "dotg = " << gamma << endl;
gtsam::print(d, "d"); gtsam::print(d, "d");
gtsam::print(Ad, "Ad"); gtsam::print(Ad, "Ad");
} }
/** step the solution */ /** step the solution */
double takeOptimalStep(V& x) { double takeOptimalStep(V& x) {
// TODO: can we use gamma instead of dot(d,g) ????? Answer not trivial
double alpha = -dot(d, g) / dot(Ad, Ad); // calculate optimal step-size double alpha = -dot(d, g) / dot(Ad, Ad); // calculate optimal step-size
axpy(alpha, d, x); // // do step in new search direction, x += alpha*d axpy(alpha, d, x); // // do step in new search direction, x += alpha*d
return alpha; return alpha;
@ -68,8 +70,9 @@ namespace gtsam {
k += 1; // increase iteration number k += 1; // increase iteration number
double alpha = takeOptimalStep(x); double alpha = takeOptimalStep(x);
print(x);
if (k == maxIterations) return true; //----------------------------------> if (k >= maxIterations) return true; //---------------------------------->
// update gradient (or re-calculate at reset time) // update gradient (or re-calculate at reset time)
if (k % reset == 0) if (k % reset == 0)
@ -80,8 +83,8 @@ namespace gtsam {
// check for convergence // check for convergence
double new_gamma = dot(g, g); double new_gamma = dot(g, g);
if (verbose) cout << "iteration " << k << ": alpha = " << alpha << ", dotg = " << new_gamma << endl; if (verbose) cout << "iteration " << k << ": alpha = " << alpha
// print(); << ", dotg = " << new_gamma << endl;
if (new_gamma < threshold) return true; //----------------------------------> if (new_gamma < threshold) return true; //---------------------------------->
// calculate new search direction // calculate new search direction
@ -89,11 +92,11 @@ namespace gtsam {
d = g; d = g;
else { else {
double beta = new_gamma / gamma; double beta = new_gamma / gamma;
gamma = new_gamma;
// d = g + d*beta; // d = g + d*beta;
scal(beta, d); scal(beta, d);
axpy(1.0, g, d); axpy(1.0, g, d);
} }
gamma = new_gamma;
// In-place recalculation Ad <- A*d to avoid re-allocating Ad // In-place recalculation Ad <- A*d to avoid re-allocating Ad
Ab.multiplyInPlace(d, Ad); Ab.multiplyInPlace(d, Ad);
@ -109,11 +112,12 @@ namespace gtsam {
V conjugateGradients(const S& Ab, V x, bool verbose, double epsilon, V conjugateGradients(const S& Ab, V x, bool verbose, double epsilon,
double epsilon_abs, size_t maxIterations, bool steepest = false) { double epsilon_abs, size_t maxIterations, bool steepest = false) {
CGState<S, V, E> state(Ab, x, verbose, epsilon, epsilon_abs, maxIterations, steepest); CGState<S, V, E> state(Ab, x, verbose, epsilon, epsilon_abs, maxIterations,
steepest);
if (state.gamma < state.threshold) return x; if (state.gamma < state.threshold) return x;
if (verbose) cout << "CG: epsilon = " << epsilon << ", maxIterations = " if (verbose) cout << "CG: epsilon = " << epsilon << ", maxIterations = "
<< maxIterations << ", ||g0||^2 = " << state.gamma << state.maxIterations << ", ||g0||^2 = " << state.gamma
<< ", threshold = " << state.threshold << endl; << ", threshold = " << state.threshold << endl;
// loop maxIterations times // loop maxIterations times