TrustRegion.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
// This code was originally translated from Fortran code on
// https://ccpforge.cse.rl.ac.uk/gf/project/ral_nlls June 2016
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/RalNlls/TrustRegion.h"
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <string>
 
namespace Mantid::CurveFitting::NLLS {
 
const double EPSILON_MCH = std::numeric_limits<double>::epsilon();
 
void matmultInner(const DoubleFortranMatrix &J, DoubleFortranMatrix &A) {
  auto n = J.len2();
  A.allocate(n, n);
 
  A.mutator() = J.inspector().transpose() * J.inspector();
}
 
void getSvdJ(const DoubleFortranMatrix &J, double &s1, double &sn) {
  Eigen::BDCSVD<Eigen::MatrixXd> svd(J.inspector());
  auto S = svd.singularValues();
 
  s1 = S(0);
  sn = S(S.size() - 1);
}
 
double norm2(const DoubleFortranVector &v) {
  if (v.size() == 0)
    return 0.0;
  return v.norm();
}
 
void multJ(const DoubleFortranMatrix &J, const DoubleFortranVector &x, DoubleFortranVector &Jx) {
  // dgemv('N',m,n,alpha,J,m,x,1,beta,Jx,1);
  if (Jx.len() != J.len1()) {
    Jx.allocate(J.len1());
  }
 
  Jx.mutator() = J.inspector() * x.inspector();
}
 
void multJt(const DoubleFortranMatrix &J, const DoubleFortranVector &x, DoubleFortranVector &Jtx) {
  // dgemv('T',m,n,alpha,J,m,x,1,beta,Jtx,1)
  if (Jtx.len() != J.len2()) {
    Jtx.allocate(J.len2());
  }
 
  Jtx.mutator() = J.inspector().transpose() * x.inspector();
}
 
double dotProduct(const DoubleFortranVector &x, const DoubleFortranVector &y) { return x.dot(y); }
 
double evaluateModel(const DoubleFortranVector &f, const DoubleFortranMatrix &J, const DoubleFortranMatrix &hf,
                     const DoubleFortranVector &d, const nlls_options &options, evaluate_model_work &w) {
 
  // Jd = J*d
  multJ(J, d, w.Jd);
 
  // First, get the base
  // 0.5 (f^T f + f^T J d + d^T' J ^T J d )
  DoubleFortranVector temp = f;
  temp += w.Jd;
  w.md_gn = 0.5 * pow(norm2(temp), 2);
  double md = 0.0;
  switch (options.model) {
  case 1: // first-order (no Hessian)
    md = w.md_gn;
    break;
  default:
    // these have a dynamic H -- recalculate
    // H = J^T J + HF, HF is (an approx?) to the Hessian
    multJ(hf, d, w.Hd);
    md = w.md_gn + 0.5 * dotProduct(d, w.Hd);
  }
  return md;
}
 
double calculateRho(double normf, double normfnew, double md, const nlls_options &options) {
  UNUSED_ARG(options);
  auto actual_reduction = (0.5 * pow(normf, 2)) - (0.5 * pow(normfnew, 2));
  auto predicted_reduction = ((0.5 * pow(normf, 2)) - md);
  double rho = 0.0;
  if (fabs(actual_reduction) < 10 * EPSILON_MCH) {
    rho = ONE;
  } else if (fabs(predicted_reduction) < 10 * EPSILON_MCH) {
    rho = ONE;
  } else {
    rho = actual_reduction / predicted_reduction;
  }
  return rho;
}
 
void rankOneUpdate(DoubleFortranMatrix &hf, NLLS_workspace &w) {
 
  auto yts = dotProduct(w.d, w.y);
  if (fabs(yts) < sqrt(10 * EPSILON_MCH)) {
    return;
  }
 
  multJ(hf, w.d, w.Sks); // hfs = S_k * d
 
  w.ysharpSks = w.y_sharp;
  w.ysharpSks -= w.Sks;
 
  // now, let's scale hd (Nocedal and Wright, Section 10.2)
  auto dSks = fabs(dotProduct(w.d, w.Sks));
  auto alpha = fabs(dotProduct(w.d, w.y_sharp)) / dSks;
  alpha = std::min(ONE, alpha);
  hf *= alpha;
 
  // update S_k (again, as in N&W, Section 10.2)
 
  // hf = hf + (1/yts) (y# - Sk d)^T y:
  alpha = 1 / yts;
 
  w.hf.mutator() = alpha * w.ysharpSks.mutator() * w.y.mutator().transpose() + w.hf.mutator();
  w.hf.mutator() = alpha * w.y.mutator() * w.ysharpSks.mutator().transpose() + w.hf.mutator();
  alpha = -dotProduct(w.ysharpSks, w.d) / (pow(yts, 2));
  w.hf.mutator() = alpha * w.y.mutator() * w.y.mutator().transpose() + w.hf.mutator();
}
 
void updateTrustRegionRadius(double &rho, const nlls_options &options, NLLS_workspace &w) {
 
  switch (options.tr_update_strategy) {
  case 1: // default, step-function
    if (!std::isfinite(rho)) {
      w.Delta = std::max(options.radius_reduce, options.radius_reduce_max) * w.Delta;
      rho = -ONE; // set to be negative, so that the logic works....
    } else if (rho < options.eta_success_but_reduce) {
      // unsuccessful....reduce Delta
      w.Delta = std::max(options.radius_reduce, options.radius_reduce_max) * w.Delta;
    } else if (rho < options.eta_very_successful) {
      //  doing ok...retain status quo
    } else if (rho < options.eta_too_successful) {
      // more than very successful -- increase delta
      w.Delta = std::min(options.maximum_radius, options.radius_increase * w.normd);
      // increase based on normd = ||d||_D
      // if d is on the tr boundary, this is Delta
      // otherwise, point was within the tr, and there's no point
      // increasing the radius
    } else {
      // too successful....accept step, but don't change w.Delta
    }
    break;
  case 2: //  Continuous method
          //  Based on that proposed by Hans Bruun Nielsen, TR
          //  IMM-REP-1999-05
          //  http://www2.imm.dtu.dk/documents/ftp/tr99/tr05_99.pdf
    if (!std::isfinite(rho)) {
      w.Delta = std::max(options.radius_reduce, options.radius_reduce_max) * w.Delta;
      rho = -ONE; // set to be negative, so that the logic works....
    } else if (rho >= options.eta_too_successful) {
      // too successful....accept step, but don't change w.Delta
    } else if (rho > options.eta_successful) {
      w.Delta = w.Delta * std::min(options.radius_increase,
                                   std::max(options.radius_reduce,
                                            1 - ((options.radius_increase - 1) * (pow((1 - 2 * rho), w.tr_p)))));
      w.tr_nu = options.radius_reduce;
    } else {
      w.Delta = w.Delta * w.tr_nu;
      w.tr_nu = w.tr_nu * 0.5;
    }
    break;
  default:
    throw std::runtime_error("Bad strategy.");
  }
}
 
void testConvergence(double normF, double normJF, double normF0, double normJF0, const nlls_options &options,
                     nlls_inform &inform) {
 
  if (normF <= std::max(options.stop_g_absolute, options.stop_g_relative * normF0)) {
    inform.convergence_normf = 1;
    return;
  }
 
  if ((normJF / normF) <= std::max(options.stop_g_absolute, options.stop_g_relative * (normJF0 / normF0))) {
    inform.convergence_normg = 1;
  }
}
 
void applyScaling(const DoubleFortranMatrix &J, DoubleFortranMatrix &A, DoubleFortranVector &v,
                  DoubleFortranVector &scale, const nlls_options &options) {
  auto m = J.len1();
  auto n = J.len2();
  if (scale.len() != n) {
    scale.allocate(n);
  }
 
  switch (options.scale) {
  case 1:
  case 2:
    for (int ii = 1; ii <= n; ++ii) { // do ii = 1,n
      double temp = ZERO;
      if (options.scale == 1) {
        for (int jj = 1; jj <= m; ++jj) { // for_do(jj, 1,m)
          // get_element_of_matrix(J,m,jj,ii,Jij);
          temp = temp + pow(J(jj, ii), 2);
        }
      } else if (options.scale == 2) {
        for (int jj = 1; jj <= n; ++jj) { // for_do(jj, 1,n)
          temp = temp + pow(A(ii, jj), 2);
        }
      }
      if (temp < options.scale_min) {
        if (options.scale_trim_min) {
          temp = options.scale_min;
        } else {
          temp = ONE;
        }
      } else if (temp > options.scale_max) {
        if (options.scale_trim_max) {
          temp = options.scale_max;
        } else {
          temp = ONE;
        }
      }
      temp = sqrt(temp);
      if (options.scale_require_increase) {
        scale(ii) = std::max(temp, scale(ii));
      } else {
        scale(ii) = temp;
      }
    }
    break;
  default:
    throw std::runtime_error("Scaling error.");
  }
 
  // Now we have the w.diagonal scaling matrix, actually scale the
  // Hessian approximation and J^Tf
  for (int ii = 1; ii <= n; ++ii) { // for_do(ii, 1,n)
    double temp = scale(ii);
    v(ii) = v(ii) / temp;
    for (int jj = 1; jj <= n; ++jj) { // for_do(jj,1,n)
      A(ii, jj) = A(ii, jj) / temp;
      A(jj, ii) = A(jj, ii) / temp;
    }
  }
}
 
void allEigSymm(const DoubleFortranMatrix &A, DoubleFortranVector &ew, DoubleFortranMatrix &ev) {
  auto M = A;
  M.eigenSystem(ew, ev);
}
 
// This isn't used because we don't calculate second derivatives in Mantid
// If we start using them the method should be un-commented and used here
//
// void apply_second_order_info(int n, int m, const DoubleFortranVector &X,
//                             NLLS_workspace &w, eval_hf_type evalHF,
//                             params_base_type params,
//                             const nlls_options &options, nlls_inform &inform,
//                             const DoubleFortranVector &weights) {
//
//  if (options.exact_second_derivatives) {
//    DoubleFortranVector temp = w.f;
//    temp *= weights;
//    evalHF(inform.external_return, n, m, X, temp, w.hf, params);
//    inform.h_eval = inform.h_eval + 1;
//  } else {
//    // use the rank-one approximation...
//    rankOneUpdate(w.hf, w, n);
//  }
//}
 
} // namespace Mantid::CurveFitting::NLLS