Как использовать Eigen неподдерживаемую реализацию levenberg marquardt?

Ответ 1

Поэтому я считаю, что нашел ответы.

1) Функция может работать как вектор функции и как скаляр функции.
 Если существуют m решаемые параметры, необходимо создать или вычислить численную матрицу якобиана m x m. Чтобы сделать умножение матрицы-вектора J(x[m]).transpose*f(x[m]), вектор функции 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;
}

Ответ 2

В качестве альтернативы вы можете просто создать новый функтор, подобный этому,

struct my_functor_w_df : Eigen::NumericalDiff<my_functor> {};

а затем инициализировать экземпляр LevenbergMarquardt, используя это,

my_functor_w_df functor;
Eigen::LevenbergMarquardt<my_functor_w_df> lm(functor);

Лично я нахожу этот подход немного более чистым.

Ответ 3

Кажется, что функция более общая:

• Скажем, у вас есть модель параметров m.
• У вас есть n наблюдений, которым вы хотите соответствовать модели m-параметров в смысле наименьших квадратов.
• Якобиан, если он предусмотрен, будет n раз m.

Вам нужно будет указать n значений ошибки в fvec. Кроме того, не требуется квадрат значений f, потому что неявно предполагается, что общая функция ошибки состоит из суммы квадратов компонентов fvec.

Итак, если вы выполните следующие рекомендации и измените код на:

fvec(0) = sqrt(10.0)*(x(0)+3.0);
fvec(1) = x(1)-5.0;

Он будет сходиться в смехотворно небольшом числе итераций - как меньше 5. Я также пробовал его на более сложном примере - эталон Hahn1 в http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml с параметрами m = 7 и n = 236 наблюдений и он сходится к известному правильному решению только в 11 итерациях с численно рассчитанным якобианом.