c++ - 如何使用 Eigen 不受支持的 levenberg marquardt 实现？

Question

我正在尝试最小化以下示例函数:

F(x) = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])

最小化这种函数的一般方法是 Levenberg-Marquardt 算法。我想在 C++ 中执行这个最小化并做了一些初步测试使用 Eigen 得出预期的解决方案。

我的问题如下:我习惯于在 python 中使用 scipy.optimize.fmin_powell 进行优化。这里输入函数参数为 (func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, direc=None)。所以我可以定义一个 func(x0)，给出 x0 vector 并开始优化。如果需要我可以改变优化参数。

现在，Eigen Lev-Marq 算法以不同的方式工作。我需要定义一个函数 vector (为什么？)此外，我无法设置优化参数。根据:
http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1LevenbergMarquardt.html
我应该可以使用 setEpsilon() 和其他设置函数。

但是当我有以下代码时:

my_functor functor;
Eigen::NumericalDiff<my_functor> numDiff(functor);
Eigen::LevenbergMarquardt<Eigen::NumericalDiff<my_functor>,double> lm(numDiff);
lm.setEpsilon(); //doesn't exist!

所以我有两个问题:

1. 为什么需要函数 vector ，为什么函数标量不够？
   我搜索过答案的引用文献:
   http://www.ultimatepp.org/reference$Eigen_demo$en-us.html
   http://www.alglib.net/optimization/levenbergmarquardt.php

2. 如何使用set函数设置优化参数？

Answer 1

所以我相信我已经找到了答案。

1) 函数既可以作为函数 vector ，也可以作为函数标量。
如果有 m 个可解参数，则需要创建或计算一个 m x m 的雅可比矩阵。为了进行矩阵 vector 乘法 J(x[m]).transpose*f(x[m]) 函数 vector f(x) 应该有 m 个项目。这可以是 m 不同的功能，但我们也可以给 f1 完整的功能，使其他项目 0。

2) 可以使用lm.parameters.maxfev = 2000;

来设置和读取参数

两个答案都已在以下示例代码中进行了测试:

#include <iostream>
#include <Eigen/Dense>

#include <unsupported/Eigen/NonLinearOptimization>
#include <unsupported/Eigen/NumericalDiff>

// Generic functor
template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>
struct Functor
{
typedef _Scalar Scalar;
enum {
    InputsAtCompileTime = NX,
    ValuesAtCompileTime = NY
};
typedef Eigen::Matrix<Scalar,InputsAtCompileTime,1> InputType;
typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;

int m_inputs, m_values;

Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}

int inputs() const { return m_inputs; }
int values() const { return m_values; }

};

struct my_functor : Functor<double>
{
my_functor(void): Functor<double>(2,2) {}
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
{
    // Implement y = 10*(x0+3)^2 + (x1-5)^2
    fvec(0) = 10.0*pow(x(0)+3.0,2) +  pow(x(1)-5.0,2);
    fvec(1) = 0;

    return 0;
}
};


int main(int argc, char *argv[])
{
Eigen::VectorXd x(2);
x(0) = 2.0;
x(1) = 3.0;
std::cout << "x: " << x << std::endl;

my_functor functor;
Eigen::NumericalDiff<my_functor> numDiff(functor);
Eigen::LevenbergMarquardt<Eigen::NumericalDiff<my_functor>,double> lm(numDiff);
lm.parameters.maxfev = 2000;
lm.parameters.xtol = 1.0e-10;
std::cout << lm.parameters.maxfev << std::endl;

int ret = lm.minimize(x);
std::cout << lm.iter << std::endl;
std::cout << ret << std::endl;

std::cout << "x that minimizes the function: " << x << std::endl;

std::cout << "press [ENTER] to continue " << std::endl;
std::cin.get();
return 0;
}