ENH: Eigen-backed Levenberg-Marquardt; remove netlib MINPACK

Modules/Numerics/Optimizers/include/itkEigenLevenbergMarquardtEngine.h

@@ -0,0 +1,75 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  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
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkEigenLevenbergMarquardtEngine_h
+#define itkEigenLevenbergMarquardtEngine_h
+
+#include "ITKOptimizersExport.h"
+#include "vnl/vnl_vector.h"
+#include <functional>
+
+namespace itk::detail
+{
+/** \brief Tolerances and limits for the Eigen Levenberg-Marquardt engine.
+ *
+ * Fields mirror the MINPACK lmdif/lmder parameters that vnl_levenberg_marquardt
+ * also exposes, so the two backends accept identical settings.
+ */
+struct EigenLevenbergMarquardtOptions
+{
+  unsigned long maxFunctionEvaluations{ 2000 };
+  double        xTolerance{ 1e-8 };
+  double        gradientTolerance{ 1e-5 };
+  double        functionTolerance{ 1e-8 };
+  double        epsilonFunction{ 1e-11 };
+};
+
+/** \brief Result of an Eigen Levenberg-Marquardt minimization. */
+struct EigenLevenbergMarquardtResult
+{
+  vnl_vector<double> solution;
+  int                status{ 0 };
+  bool               converged{ false };
+  unsigned int       numberOfEvaluations{ 0 };
+  double             residualNorm{ 0.0 };
+};
+
+/** \brief Solve a nonlinear least-squares problem with Eigen's MINPACK-port
+ * Levenberg-Marquardt (unsupported/Eigen/NonLinearOptimization).
+ *
+ * This is the algorithm-identical counterpart to vnl_levenberg_marquardt
+ * (both are ports of MINPACK lmdif/lmder). Eigen is confined to the
+ * implementation; no Eigen type appears in this interface (see issue #6230).
+ *
+ * \param numberOfParameters  n, the size of x.
+ * \param numberOfResiduals   m, the size of the residual vector.
+ * \param residual            Fills m residuals from n parameters (raw pointers).
+ * \param jacobian            Fills the m-by-n Jacobian, row-major, from n
+ *                            parameters. If empty, forward differences are used.
+ * \param initialPosition     Starting x (size n).
+ * \param options             Tolerances / evaluation limit.
+ */
+ITKOptimizers_EXPORT EigenLevenbergMarquardtResult
+EigenLevenbergMarquardtSolve(unsigned int                                          numberOfParameters,
+                             unsigned int                                          numberOfResiduals,
+                             const std::function<void(const double *, double *)> & residual,
+                             const std::function<void(const double *, double *)> & jacobian,
+                             const vnl_vector<double> &                            initialPosition,
+                             const EigenLevenbergMarquardtOptions &                options);
+} // namespace itk::detail
+
+#endif

Modules/Numerics/Optimizers/include/itkLevenbergMarquardtOptimizer.h

@@ -19,9 +19,10 @@
 #define itkLevenbergMarquardtOptimizer_h
 
 #include "itkMultipleValuedNonLinearVnlOptimizer.h"
-#include "vnl/algo/vnl_levenberg_marquardt.h"
 #include "ITKOptimizersExport.h"
-#include <memory> // For unique_ptr.
+#include <string>
+
+class vnl_levenberg_marquardt;
 
 namespace itk
 {
@@ -55,12 +56,17 @@ class ITKOptimizers_EXPORT LevenbergMarquardtOptimizer : public MultipleValuedNo
   /** InternalParameters type alias. */
   using InternalParametersType = vnl_vector<double>;
 
-  /** Internal optimizer type. */
+#if !defined(ITK_LEGACY_REMOVE)
+  /** Internal optimizer type.
+   * \deprecated The minimization is performed by an Eigen-backed engine; this
+   * alias is retained only for source compatibility. */
   using InternalOptimizerType = vnl_levenberg_marquardt;
 
-  /** Method for getting access to the internal optimizer. */
-  vnl_levenberg_marquardt *
-  GetOptimizer() const;
+  /** \deprecated Access to the former internal vnl_levenberg_marquardt is no
+   * longer meaningful: the optimizer is backed by the Eigen MINPACK-port
+   * Levenberg-Marquardt. Returns nullptr. */
+  itkLegacyMacro(vnl_levenberg_marquardt * GetOptimizer() const);
+#endif
 
   /** Start optimization with an initial value. */
   void
@@ -96,12 +102,12 @@ class ITKOptimizers_EXPORT LevenbergMarquardtOptimizer : public MultipleValuedNo
   using CostFunctionAdaptorType = Superclass::CostFunctionAdaptorType;
 
 private:
-  bool                                   m_OptimizerInitialized{};
-  std::unique_ptr<InternalOptimizerType> m_VnlOptimizer;
-  unsigned int                           m_NumberOfIterations{};
-  double                                 m_ValueTolerance{};
-  double                                 m_GradientTolerance{};
-  double                                 m_EpsilonFunction{};
+  bool         m_OptimizerInitialized{};
+  unsigned int m_NumberOfIterations{};
+  double       m_ValueTolerance{};
+  double       m_GradientTolerance{};
+  double       m_EpsilonFunction{};
+  std::string  m_StopConditionDescription{};
 };
 } // end namespace itk
 

Modules/Numerics/Optimizers/itk-module.cmake

@@ -11,8 +11,11 @@ itk_module(
   ENABLE_SHARED
   DEPENDS
     ITKStatistics
+  COMPILE_DEPENDS
+    ITKEigen3
   TEST_DEPENDS
     ITKTransform
     ITKTestKernel
+    ITKGoogleTest
   DESCRIPTION "${DOCUMENTATION}"
 )

Modules/Numerics/Optimizers/src/CMakeLists.txt

@@ -4,6 +4,7 @@ set(
   itkConjugateGradientOptimizer.cxx
   itkCumulativeGaussianCostFunction.cxx
   itkCumulativeGaussianOptimizer.cxx
+  itkEigenLevenbergMarquardtEngine.cxx
   itkExhaustiveOptimizer.cxx
   itkFRPROptimizer.cxx
   itkGradientDescentOptimizer.cxx

Modules/Numerics/Optimizers/src/itkEigenLevenbergMarquardtEngine.cxx

@@ -0,0 +1,155 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  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
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "itkEigenLevenbergMarquardtEngine.h"
+
+#include "itk_eigen.h"
+#include "itkeigen/Eigen/Core"
+#include "itkeigen/unsupported/Eigen/NonLinearOptimization"
+#include "itkeigen/unsupported/Eigen/NumericalDiff"
+
+#include <vector>
+
+namespace itk::detail
+{
+namespace
+{
+/** Adapts the raw-pointer residual/Jacobian callbacks to an Eigen functor. */
+struct CallbackFunctor
+{
+  using Scalar = double;
+  using InputType = Eigen::VectorXd;
+  using ValueType = Eigen::VectorXd;
+  using JacobianType = Eigen::MatrixXd;
+  enum
+  {
+    InputsAtCompileTime = Eigen::Dynamic,
+    ValuesAtCompileTime = Eigen::Dynamic
+  };
+
+  int                                                   m_NumberOfParameters;
+  int                                                   m_NumberOfResiduals;
+  const std::function<void(const double *, double *)> & m_Residual;
+  const std::function<void(const double *, double *)> & m_Jacobian;
+
+  int
+  inputs() const
+  {
+    return m_NumberOfParameters;
+  }
+  int
+  values() const
+  {
+    return m_NumberOfResiduals;
+  }
+
+  int
+  operator()(const Eigen::VectorXd & x, Eigen::VectorXd & fvec) const
+  {
+    m_Residual(x.data(), fvec.data());
+    return 0;
+  }
+
+  // Eigen's Jacobian is column-major; the callback fills row-major (m-by-n).
+  int
+  df(const Eigen::VectorXd & x, Eigen::MatrixXd & fjac) const
+  {
+    std::vector<double> rowMajor(static_cast<size_t>(m_NumberOfResiduals) * m_NumberOfParameters);
+    m_Jacobian(x.data(), rowMajor.data());
+    for (int i = 0; i < m_NumberOfResiduals; ++i)
+    {
+      for (int j = 0; j < m_NumberOfParameters; ++j)
+      {
+        fjac(i, j) = rowMajor[static_cast<size_t>(i) * m_NumberOfParameters + j];
+      }
+    }
+    return 0;
+  }
+};
+
+bool
+StatusIsConverged(int status)
+{
+  // LevenbergMarquardtSpace::Status: 1..4 are the success codes
+  // (RelativeReductionTooSmall, RelativeErrorTooSmall,
+  //  RelativeErrorAndReductionTooSmall, CosinusTooSmall).
+  return status >= 1 && status <= 4;
+}
+} // namespace
+
+EigenLevenbergMarquardtResult
+EigenLevenbergMarquardtSolve(unsigned int                                          numberOfParameters,
+                             unsigned int                                          numberOfResiduals,
+                             const std::function<void(const double *, double *)> & residual,
+                             const std::function<void(const double *, double *)> & jacobian,
+                             const vnl_vector<double> &                            initialPosition,
+                             const EigenLevenbergMarquardtOptions &                options)
+{
+  Eigen::VectorXd x(numberOfParameters);
+  for (unsigned int i = 0; i < numberOfParameters; ++i)
+  {
+    x[static_cast<Eigen::Index>(i)] = initialPosition[i];
+  }
+
+  CallbackFunctor functor{
+    static_cast<int>(numberOfParameters), static_cast<int>(numberOfResiduals), residual, jacobian
+  };
+
+  int          status = 0;
+  unsigned int numberOfEvaluations = 0;
+  double       residualNorm = 0.0;
+
+  if (jacobian)
+  {
+    // Analytic Jacobian -> lmder.
+    Eigen::LevenbergMarquardt<CallbackFunctor> lm(functor);
+    lm.parameters.maxfev = static_cast<int>(options.maxFunctionEvaluations);
+    lm.parameters.xtol = options.xTolerance;
+    lm.parameters.ftol = options.functionTolerance;
+    lm.parameters.gtol = options.gradientTolerance;
+    status = lm.minimize(x);
+    numberOfEvaluations = static_cast<unsigned int>(lm.nfev);
+    residualNorm = lm.fvec.norm();
+  }
+  else
+  {
+    // Forward-difference Jacobian -> lmdif.
+    Eigen::NumericalDiff<CallbackFunctor>                            numericalDiff(functor, options.epsilonFunction);
+    Eigen::LevenbergMarquardt<Eigen::NumericalDiff<CallbackFunctor>> lm(numericalDiff);
+    lm.parameters.maxfev = static_cast<int>(options.maxFunctionEvaluations);
+    lm.parameters.xtol = options.xTolerance;
+    lm.parameters.ftol = options.functionTolerance;
+    lm.parameters.gtol = options.gradientTolerance;
+    status = lm.minimize(x);
+    numberOfEvaluations = static_cast<unsigned int>(lm.nfev);
+    residualNorm = lm.fvec.norm();
+  }
+
+  EigenLevenbergMarquardtResult result;
+  result.solution = vnl_vector<double>(numberOfParameters);
+  for (unsigned int i = 0; i < numberOfParameters; ++i)
+  {
+    result.solution[i] = x[static_cast<Eigen::Index>(i)];
+  }
+  result.status = status;
+  result.converged = StatusIsConverged(status);
+  result.numberOfEvaluations = numberOfEvaluations;
+  result.residualNorm = residualNorm;
+  return result;
+}
+} // namespace itk::detail