I have been happily using Eigen in a project, and I am impressed! I now have need for a nonlinear solver, so I tried the HybridNonLinearSolver from the NonLinearOptimization library in Eigen unsupported. In my problem, I will need to numerically estimate the Jacobian. Unfortunately, the routine to compute the numerical Jacobian, fdjac1.h, asserts on line 25 in Eigen 3.0.3 when the input and values are not of the same length. This seems odd since it is normally permissible to have more values than inputs for a least-squares solver. I have attached toy code below that causes this assert.

Any hints on what I am doing wrong? The regression test does not seem to cover the case with numerically computed Jacobian and unequal input and value vector lengths.

Code: Select all

#include <unsupported/Eigen/NonLinearOptimization>

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Mat;
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> Vec;

class Functor {
public:
    typedef double Scalar;
    enum {
        InputsAtCompileTime = 3,
        ValuesAtCompileTime = Eigen::Dynamic
    };
    typedef Eigen::Matrix<Scalar, InputsAtCompileTime, 1> InputType;
    typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, 1> ValueType;
    typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, InputsAtCompileTime> JacobianType;

    Functor()
        : m_inputs(InputsAtCompileTime)
        , m_values(10)
        , xTruth(m_inputs)
        , vTruth(m_values)
    {
        xTruth << 1, 2, -3; // quadradic.
        op(xTruth, vTruth);
    }

    int inputs() { return m_inputs; }
    int values() { return m_values; }
   
    // Evaluate a polynomial over t (range 0 to 1)
    // y = x(0) + x(1) * t + x(2) * t^2
    void op(Vec const & x, Vec & value) const {
        value.resize(m_values);
        for (int i = 0; i < m_values; ++i) {
            Scalar t = static_cast<Scalar>(i) / (m_values - 1);
            Scalar tPower = t;
            Scalar v = x(0); // offset
            for (int k = 1; k < m_inputs; ++k) {
                v += x(k) * tPower;
                tPower *= t;
            }
            value(i) = v;
        }
    }

    int operator()(Vec const & x, Vec & fvec) const {
        op(x, fvec);
        fvec -= vTruth;
        return 0;
    }

private:
    int m_inputs;
    int m_values;
    Vec xTruth;
    Vec vTruth;
};

int main(int argc, char* argv[]) {
    Vec x = Vec::Zero(3);
    Functor functor;
    Eigen::HybridNonLinearSolver<Functor> solver(functor);
    int info = solver.hybrd1(x);
}


HybridNonLinearSolver is not a least-squares solver. As its documentation on http://eigen.tuxfamily.org/dox/unsuppor ... olver.html says:

Finds a zero of a system of n nonlinear functions in n variables by a modification of the Powell hybrid method ("dogleg").

So you need as many functions as variables.

indeed, if you want to solve for non linear least-square you should rather use the Levenberg-Marquart solver of the same module.

Thanks for the very quick replies! I originally tried the Levenberg-Marquart solver and also had compilation issues. I think that I have traced this LM issue down to my JacobianType specification: the LM method assumes that both dimensions are Eigen::Dynamic, but my code had one dimension specified at compilation time. I will post my results in the next couple of days. Thanks again!