Do not implicitly negate the step in the TrustRegionMinimizer.

In the TrustRegionMinimizer, the step is currently implicitly negated.
This is done so that the linearized residual is |r - J*step|^2, which
corresponds to J*step = r, so neither J nor r have to be modified.
However, it leads to the rather unintuitive situation that the strategy
returns a step in positive gradient direction, which you would expect to
increase the function value. One way is to rename the "step" parameter in
the strategy to "negative_step" and document it.
This patch instead moves the negation inside the strategy, just around
the linear solver call, so that it is done in a local context and easier
to document.

Change-Id: Idb258149a01f61c64e22128ea221c5a30cd89c89

internal/ceres/dogleg_strategy.cc

diff --git a/internal/ceres/dogleg_strategy.cc b/internal/ceres/dogleg_strategy.cc
index 620bda7..44d484f 100644
--- a/internal/ceres/dogleg_strategy.cc
+++ b/internal/ceres/dogleg_strategy.cc
@@ -142,7 +142,7 @@
 }
 
 void DoglegStrategy::ComputeCauchyPoint(SparseMatrix* jacobian) {
-  // alpha * gradient is the Cauchy point.
+  // alpha * -gradient is the Cauchy point.
   Vector Jg(jacobian->num_rows());
   Jg.setZero();
   // The Jacobian is scaled implicitly by computing J * (D^-1 * (D^-1 * g))
@@ -173,7 +173,7 @@
   // the trust region. Rescale the Cauchy point to the trust region
   // and return.
   if  (gradient_norm * alpha_ >= radius_) {
-    dogleg_step = (radius_ / gradient_norm) * gradient_;
+    dogleg_step = -(radius_ / gradient_norm) * gradient_;
     dogleg_step_norm_ = radius_;
     dogleg_step.array() /= diagonal_.array();
     VLOG(3) << "Cauchy step size: " << dogleg_step_norm_
@@ -186,15 +186,15 @@
   // points and the point on it which intersects the trust region
   // boundary.
 
-  // a = alpha * gradient
+  // a = alpha * -gradient
   // b = gauss_newton_step
-  const double b_dot_a = alpha_ * gradient_.dot(gauss_newton_step_);
+  const double b_dot_a = -alpha_ * gradient_.dot(gauss_newton_step_);
   const double a_squared_norm = pow(alpha_ * gradient_norm, 2.0);
   const double b_minus_a_squared_norm =
       a_squared_norm - 2 * b_dot_a + pow(gauss_newton_norm, 2);
 
   // c = a' (b - a)
-  //   = alpha * gradient' gauss_newton_step - alpha^2 |gradient|^2
+  //   = alpha * -gradient' gauss_newton_step - alpha^2 |gradient|^2
   const double c = b_dot_a - a_squared_norm;
   const double d = sqrt(c * c + b_minus_a_squared_norm *
                         (pow(radius_, 2.0) - a_squared_norm));
@@ -203,7 +203,7 @@
       (c <= 0)
       ? (d - c) /  b_minus_a_squared_norm
       : (radius_ * radius_ - a_squared_norm) / (d + c);
-  dogleg_step = (alpha_ * (1.0 - beta)) * gradient_ + beta * gauss_newton_step_;
+  dogleg_step = (-alpha_ * (1.0 - beta)) * gradient_ + beta * gauss_newton_step_;
   dogleg_step_norm_ = dogleg_step.norm();
   dogleg_step.array() /= diagonal_.array();
   VLOG(3) << "Dogleg step size: " << dogleg_step_norm_
@@ -255,6 +255,9 @@
     lm_diagonal_ = diagonal_ * std::sqrt(mu_);
     solve_options.D = lm_diagonal_.data();
 
+    // As in the LevenbergMarquardtStrategy, solve Jy = r instead
+    // of Jx = -r and later set x = -y to avoid having to modify
+    // either jacobian or residuals.
     InvalidateArray(n, gauss_newton_step_.data());
     linear_solver_summary = linear_solver_->Solve(jacobian,
                                                   residuals,
@@ -277,7 +280,7 @@
   //   = - D (J^T J)^-1 D D^-1 g
   //   = D -(J^T J)^-1 g
   //
-  gauss_newton_step_.array() *= diagonal_.array();
+  gauss_newton_step_.array() *= -diagonal_.array();
 
   return linear_solver_summary;
 }

internal/ceres/levenberg_marquardt_strategy.cc

diff --git a/internal/ceres/levenberg_marquardt_strategy.cc b/internal/ceres/levenberg_marquardt_strategy.cc
index 298748f..576648a 100644
--- a/internal/ceres/levenberg_marquardt_strategy.cc
+++ b/internal/ceres/levenberg_marquardt_strategy.cc
@@ -98,12 +98,18 @@
   // to happen for the DENSE_QR and then DENSE_SCHUR solver when
   // the Jacobin is severly rank deficient and mu is too small.
   InvalidateArray(num_parameters, step);
+
+  // Instead of solving Jx = -r, solve Jy = r.
+  // Then x can be found as x = -y, but the inputs jacobian and residuals
+  // do not need to be modified.
   LinearSolver::Summary linear_solver_summary =
       linear_solver_->Solve(jacobian, residuals, solve_options, step);
   if (linear_solver_summary.termination_type == FAILURE ||
       !IsArrayValid(num_parameters, step)) {
     LOG(WARNING) << "Linear solver failure. Failed to compute a finite step.";
     linear_solver_summary.termination_type = FAILURE;
+  } else {
+    VectorRef(step, num_parameters) *= -1.0;
   }
 
   reuse_diagonal_ = true;

internal/ceres/trust_region_minimizer.cc

diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
index d08255f..476cc9d 100644
--- a/internal/ceres/trust_region_minimizer.cc
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -271,8 +271,8 @@
 
     double new_model_cost = 0.0;
     if (strategy_summary.termination_type != FAILURE) {
-      // new_model_cost = 1/2 |J * step - f|^2
-      model_residuals = -residuals;
+      // new_model_cost = 1/2 |f + J * step|^2
+      model_residuals = residuals;
       jacobian->RightMultiply(trust_region_step.data(), model_residuals.data());
       new_model_cost = model_residuals.squaredNorm() / 2.0;
 
@@ -328,7 +328,7 @@
       const double model_cost_change = max(kEpsilon, cost - new_model_cost);
 
       // Undo the Jacobian column scaling.
-      delta =  -(trust_region_step.array() * scale.array()).matrix();
+      delta = (trust_region_step.array() * scale.array()).matrix();
       iteration_summary.step_norm = delta.norm();
 
       // Convergence based on parameter_tolerance.