Jak korzystać z niewspieranej implementacji levinbers marquardt?

Question

staram się zminimalizować następującą funkcję próbkowania:

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

normalny sposób, aby zminimalizować taką FUNCT może być algorytm Levenberg-Marquardt. Chciałbym wykonać tę minimalizację w C++ i zrobić kilka początkowych testów z Eigenem, które dały oczekiwane rozwiązanie.

Moje pytanie jest następujące: Jestem przyzwyczajony do optymalizacji w python za pomocą np. scipy.optimize.fmin_powell. Tutaj parametrami funkcji wejściowej są (func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, direc=None). Więc mogę zdefiniować func(x0), dać wektor x0 i rozpocząć optymalizację. W razie potrzeby mogę zmienić parametry optymalizacji.

Teraz algorytm Levina-Marca działa w inny sposób. Muszę zdefiniować funkcję wektorową (dlaczego?) Ponadto nie mogę ustawić parametrów optymalizacji. Zgodnie z:
http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1LevenbergMarquardt.html Powinny być w stanie korzystać z setEpsilon() i innych funkcji zestawu.

Ale kiedy mam następujący kod:

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

Więc mam 2 pytania:

1. Dlaczego jest wektorem funkcji potrzebna i dlaczego nie skalar funkcja wystarczy?
   Referencje gdzie Rozglądałem się na odpowiedź:
   http://www.ultimatepp.org/reference$Eigen_demo$en-us.html
   http://www.alglib.net/optimization/levenbergmarquardt.php

2. Jak ustawić parametry optymalizacji stosując zbiór funkcji?

Answer 1

Sądzę, że znalazłem odpowiedzi.

1) Funkcja może pracować jako wektor funkcyjny i jako funkcja skalarna.
Jeśli istnieją parametry możliwe do rozwiązania, należy utworzyć macierz Jacobian m x m lub obliczać ją liczbowo. Aby wykonać mnożenie wektorów Matrix-Vector J(x[m]).transpose*f(x[m]), wektor funkcji f(x) powinien mieć elementy m. Może to być jedna z wielu różnych funkcji, ale możemy również nadać f1 kompletną funkcję i utworzyć inne elementy 0.

2) parametry można ustawić i odczytać za pomocą lm.parameters.maxfev = 2000;

Obie odpowiedzi były testowane w poniższym przykładzie kodu:

#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; 
} 

Answer 2

Wydaje się, że funkcja jest bardziej ogólne:

1. Załóżmy, że masz model m-parametrów.
2. Masz n obserwacji, do których chcesz dopasować model m-parametrów w sensie najmniejszych kwadratów.
3. Jakobin, jeśli przewidziano, będzie n razy m.

Musisz podać n wartości błędów w fvec. Ponadto istnieje nie ma potrzeby sortowania wartości f-, ponieważ domyślnie zakłada się, że całkowita funkcja błędu składa się z sumy kwadratów komponentów fvec.

Więc, jeśli się z poniższymi wskazówkami i zmienić kod do:

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

Będzie zbiegają się w śmiesznie małej liczby powtórzeń - jak mniejsza niż 5. Próbowałem też go na bardziej złożoną przykład - Benchmark Hahn1 przy http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml z m = 7 parametrami i n = 236 obserwacjami i zbiegnie się do znanego rozwiązania w zaledwie 11 iteracjach z obliczoną numerycznie liczbą Jakobian.

Answer 3

Jako alternatywę można po prostu stworzyć nową funktor takiego,

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

a następnie zainicjować LevenbergMarquardt instancję za pomocą tak,

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

Osobiście uważam, że to podejście czystsze bit .

Answer 4

Ta odpowiedź jest rozszerzeniem dwóch istniejących odpowiedzi: 1) Dostosowałem kod źródłowy dostarczony przez @Deepfreeze, aby uwzględnić dodatkowe komentarze i dwie różne funkcje testowe. 2) Używam sugestii z @ user3361661 do przepisania funkcji celu w poprawnej formie. Jak zasugerował, że zmniejsza liczbę iteracji na moim pierwszym problemem testowym od 67 do 4.

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

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

/***********************************************************************************************/ 

// Generic functor 
// See http://eigen.tuxfamily.org/index.php?title=Functors 
// C++ version of a function pointer that stores meta-data about the function 
template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic> 
struct Functor 
{ 

    // Information that tells the caller the numeric type (eg. double) and size (input/output dim) 
    typedef _Scalar Scalar; 
    enum { // Required by numerical differentiation module 
     InputsAtCompileTime = NX, 
     ValuesAtCompileTime = NY 
    }; 

    // Tell the caller the matrix sizes associated with the input, output, and jacobian 
    typedef Eigen::Matrix<Scalar,InputsAtCompileTime,1> InputType; 
    typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 
    typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 

    // Local copy of the number of inputs 
    int m_inputs, m_values; 

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

    // Get methods for users to determine function input and output dimensions 
    int inputs() const { return m_inputs; } 
    int values() const { return m_values; } 

}; 

/***********************************************************************************************/ 

// https://en.wikipedia.org/wiki/Test_functions_for_optimization 
// Booth Function 
// Implement f(x,y) = (x + 2*y -7)^2 + (2*x + y - 5)^2 
struct BoothFunctor : Functor<double> 
{ 
    // Simple constructor 
    BoothFunctor(): Functor<double>(2,2) {} 

    // Implementation of the objective function 
    int operator()(const Eigen::VectorXd &z, Eigen::VectorXd &fvec) const { 
    double x = z(0); double y = z(1); 
    /* 
    * Evaluate the Booth function. 
    * Important: LevenbergMarquardt is designed to work with objective functions that are a sum 
    * of squared terms. The algorithm takes this into account: do not do it yourself. 
    * In other words: objFun = sum(fvec(i)^2) 
    */ 
    fvec(0) = x + 2*y - 7; 
    fvec(1) = 2*x + y - 5; 
    return 0; 
    } 
}; 

/***********************************************************************************************/ 

// https://en.wikipedia.org/wiki/Test_functions_for_optimization 
// Himmelblau's Function 
// Implement f(x,y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2 
struct HimmelblauFunctor : Functor<double> 
{ 
    // Simple constructor 
    HimmelblauFunctor(): Functor<double>(2,2) {} 

    // Implementation of the objective function 
    int operator()(const Eigen::VectorXd &z, Eigen::VectorXd &fvec) const { 
    double x = z(0); double y = z(1); 
    /* 
    * Evaluate Himmelblau's function. 
    * Important: LevenbergMarquardt is designed to work with objective functions that are a sum 
    * of squared terms. The algorithm takes this into account: do not do it yourself. 
    * In other words: objFun = sum(fvec(i)^2) 
    */ 
    fvec(0) = x * x + y - 11; 
    fvec(1) = x + y * y - 7; 
    return 0; 
    } 
}; 

/***********************************************************************************************/ 

void testBoothFun() { 
    std::cout << "Testing the Booth function..." << std::endl; 
    Eigen::VectorXd zInit(2); zInit << 1.87, 2.032; 
    std::cout << "zInit: " << zInit.transpose() << std::endl; 
    Eigen::VectorXd zSoln(2); zSoln << 1.0, 3.0; 
    std::cout << "zSoln: " << zSoln.transpose() << std::endl; 

    BoothFunctor functor; 
    Eigen::NumericalDiff<BoothFunctor> numDiff(functor); 
    Eigen::LevenbergMarquardt<Eigen::NumericalDiff<BoothFunctor>,double> lm(numDiff); 
    lm.parameters.maxfev = 1000; 
    lm.parameters.xtol = 1.0e-10; 
    std::cout << "max fun eval: " << lm.parameters.maxfev << std::endl; 
    std::cout << "x tol: " << lm.parameters.xtol << std::endl; 

    Eigen::VectorXd z = zInit; 
    int ret = lm.minimize(z); 
    std::cout << "iter count: " << lm.iter << std::endl; 
    std::cout << "return status: " << ret << std::endl; 
    std::cout << "zSolver: " << z.transpose() << std::endl; 
    std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; 
} 

/***********************************************************************************************/ 

void testHimmelblauFun() { 
    std::cout << "Testing the Himmelblau function..." << std::endl; 
    // Eigen::VectorXd zInit(2); zInit << 0.0, 0.0; // soln 1 
    // Eigen::VectorXd zInit(2); zInit << -1, 1; // soln 2 
    // Eigen::VectorXd zInit(2); zInit << -1, -1; // soln 3 
    Eigen::VectorXd zInit(2); zInit << 1, -1; // soln 4 
    std::cout << "zInit: " << zInit.transpose() << std::endl; 
    std::cout << "soln 1: [3.0, 2.0]" << std::endl; 
    std::cout << "soln 2: [-2.805118, 3.131312]" << std::endl; 
    std::cout << "soln 3: [-3.77931, -3.28316]" << std::endl; 
    std::cout << "soln 4: [3.584428, -1.848126]" << std::endl; 

    HimmelblauFunctor functor; 
    Eigen::NumericalDiff<HimmelblauFunctor> numDiff(functor); 
    Eigen::LevenbergMarquardt<Eigen::NumericalDiff<HimmelblauFunctor>,double> lm(numDiff); 
    lm.parameters.maxfev = 1000; 
    lm.parameters.xtol = 1.0e-10; 
    std::cout << "max fun eval: " << lm.parameters.maxfev << std::endl; 
    std::cout << "x tol: " << lm.parameters.xtol << std::endl; 

    Eigen::VectorXd z = zInit; 
    int ret = lm.minimize(z); 
    std::cout << "iter count: " << lm.iter << std::endl; 
    std::cout << "return status: " << ret << std::endl; 
    std::cout << "zSolver: " << z.transpose() << std::endl; 
    std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; 
} 

/***********************************************************************************************/ 

int main(int argc, char *argv[]) 
{ 

std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; 
testBoothFun(); 
testHimmelblauFun(); 
return 0; 
} 

Wyjście w linii poleceń z uruchomieniem tego skryptu testowego jest:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
Testing the Booth function... 
zInit: 1.87 2.032 
zSoln: 1 3 
max fun eval: 1000 
x tol: 1e-10 
iter count: 4 
return status: 2 
zSolver: 1 3 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
Testing the Himmelblau function... 
zInit: 1 -1 
soln 1: [3.0, 2.0] 
soln 2: [-2.805118, 3.131312] 
soln 3: [-3.77931, -3.28316] 
soln 4: [3.584428, -1.848126] 
max fun eval: 1000 
x tol: 1e-10 
iter count: 8 
return status: 2 
zSolver: 3.58443 -1.84813 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~