Làm cách nào để sử dụng triển khai thực hiện levenberg không được hỗ trợ bởi Eigen?

Question

Tôi đang cố gắng để giảm thiểu một hàm mẫu sau:

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

Một cách thông thường để giảm thiểu một funct như vậy có thể là thuật toán Levenberg-Marquardt. Tôi muốn thực hiện việc giảm thiểu này trong C++ và đã thực hiện một số thử nghiệm ban đầu với Eigen dẫn đến giải pháp dự kiến.

Câu hỏi của tôi là như sau: Tôi được sử dụng để tối ưu hóa trong trăn bằng cách sử dụng tức là scipy.optimize.fmin_powell. Ở đây các tham số hàm đầu vào là (func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, direc=None). Vì vậy, tôi có thể xác định một func(x0), cung cấp cho vector x0 và bắt đầu tối ưu hóa. Nếu cần, tôi có thể thay đổi thông số tối ưu hóa.

Bây giờ thuật toán Eigen Lev-Marq hoạt động theo một cách khác. Tôi cần phải xác định một chức năng vector (tại sao?) Hơn nữa tôi không thể quản lý để thiết lập các thông số tối ưu hóa. Theo:
http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1LevenbergMarquardt.html
Tôi có thể sử dụng setEpsilon() và các chức năng được đặt khác.

Nhưng khi tôi có đoạn mã sau:

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

Vì vậy, tôi có 2 câu hỏi:

1. Tại sao một vector chức năng cần thiết và tại sao không phải là một hàm vô hướng là đủ?
   Tài liệu tham khảo mà tôi đã tìm kiếm một câu trả lời:
   http://www.ultimatepp.org/reference$Eigen_demo$en-us.html
   http://www.alglib.net/optimization/levenbergmarquardt.php

2. Làm thế nào để thiết lập các thông số tối ưu hóa bằng cách sử dụng các chức năng thiết lập?

Answer 1

Vì vậy, tôi tin rằng tôi đã tìm thấy câu trả lời.

1) Chức năng có thể hoạt động như một vectơ chức năng và như một hàm vô hướng.
Nếu có m thông số có thể giải được, ma trận Jacobian của m x m cần được tạo hoặc tính toán bằng số. Để thực hiện phép nhân Matrix-Vector J(x[m]).transpose*f(x[m]), vector chức năng f(x) cần có m mục. Đây có thể là m các chức năng khác nhau, nhưng chúng tôi cũng có thể cung cấp chức năng hoàn chỉnh và làm cho các mục khác f10.

2) Các thông số có thể được thiết lập và đọc bằng lm.parameters.maxfev = 2000;

Cả hai câu trả lời đã được thử nghiệm trong các mã ví dụ sau:

#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

Dường như chức năng là tổng quát hơn:

1. Giả sử bạn có mô hình tham số m.
2. Bạn có n quan sát mà bạn muốn phù hợp với mô hình tham số m theo nghĩa tối thiểu.
3. Jacobian, nếu được cung cấp, sẽ là n lần m.

Bạn sẽ cần cung cấp n giá trị lỗi trong fvec. Ngoài ra, có không cần vuông giá trị f vì nó được giả định ngầm định rằng hàm lỗi tổng thể được tạo thành từ tổng các bình phương của các thành phần fvec.

Vì vậy, nếu bạn làm theo những hướng dẫn này và thay đổi mã để:

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

Nó sẽ hội tụ trong một số ridiculously nhỏ lặp đi lặp lại - như ít hơn 5. Tôi cũng đã thử nó trên một ví dụ phức tạp hơn - chuẩn Hahn1 tại http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml với m = 7 tham số và n = 236 quan sát và nó hội tụ với giải pháp đúng đã biết chỉ trong 11 lần lặp với Jacobian được tính toán bằng số.

Answer 3

Là một thay thế bạn có thể chỉ cần tạo một functor mới như thế này,

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

và sau đó khởi tạo instance LevenbergMarquardt sử dụng như thế này,

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

Cá nhân, tôi tìm cách tiếp cận này một chút bụi .

Answer 4

Câu trả lời này là phần mở rộng của hai câu trả lời hiện có: 1) Tôi đã điều chỉnh mã nguồn do @Deepfreeze cung cấp để bao gồm các nhận xét bổ sung và hai chức năng kiểm tra khác nhau. 2) Tôi sử dụng gợi ý từ @ user3361661 để viết lại hàm mục tiêu theo đúng mẫu. Như ông đề nghị, nó làm giảm số lần lặp về vấn đề thử nghiệm đầu tiên của tôi từ 67 đến 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; 
} 

Sản lượng tại dòng lệnh chạy kịch bản thử nghiệm này là:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
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 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~