PatWie / CppNumericalSolvers

a lightweight header-only C++17 library of numerical optimization methods for nonlinear functions based on Eigen
MIT License
867 stars 201 forks source link

L-BFGS-B - cannot instantiate a solver with Dim=Eigen::Dynamic #53

Closed nikitaved closed 8 years ago

nikitaved commented 8 years ago

Hi!

Eigen/src/Core/CwiseNullaryOp.h failed at YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR assertion. It is a constant matrix (Eigen::DenseBase::ConstantReturnType) which is expected somewhere (getters of bounds, I presume) and its size must be known beforehand.

My problem is defined like this:

class OptWrtTProblem : public cppoptlib::BoundedProblem<Double> {
    public:
        using Superclass    = cppoptlib::BoundedProblem<Double>;
        using DVector       = TVector;
        using EigenMD       = Map<Matrix>;
        using EigenMDconst  = Map<const Matrix>;
        using DVectorM      = Map<DVector>;
        using DVectorMconst = Map<const DVector>;

    private:
        const Matrix& mD;
        const Matrix& mTfixed;
        const Matrix& mA;
        double lambda;

        const size_t m;
        const size_t n;
        const size_t r;

        /* some time savers */
        const Matrix& onesmr;
        Matrix mAAt;
        Matrix mDAt;
        Double regTfixed;
        Matrix gradRegTfixed;
        Double dotGradTfixed;

    public:
        OptWrtTProblem(const Matrix& mDIn,
                       const Matrix& mTfixedIn,
                       const Matrix& mAIn,
                       double lambdaIn,
                       const Matrix& onesmrIn,
                       const DVector& lb,
                       const DVector& ub
        )
        : Superclass(lb, ub),
          mD{mDIn},
          mTfixed{mTfixedIn},
          mA{mAIn},
          lambda{lambdaIn},
          onesmr{onesmrIn},
          m(mD.rows()),
          n(mD.cols()),
          r(mA.rows())
        {
            mAAt = mA * mA.transpose();
            mDAt = mD * mA.transpose();

            regTfixed = mTfixed.sum() - mTfixed.squaredNorm();
            gradRegTfixed = onesmr - 2 * mTfixed;
            dotGradTfixed = gradRegTfixed.cwiseProduct(mTfixed).sum();
        }

        Double value(const DVector& mTVec) {
            EigenMDconst mT(mTVec.data(), m, r);
            Double loss = (mD - mT * mA).squaredNorm();
            return 0.5 * (loss + lambda * (regTfixed
                        + gradRegTfixed.cwiseProduct(mT).sum() - dotGradTfixed));
        }

        void gradient(const DVector& mTVec, DVector& gradVec) {
            EigenMDconst mT(mTVec.data(), m, r);
            EigenMD grad(gradVec.data(), m, r);
            grad = (mT * mAAt - mDAt) + 0.5 * lambda * gradRegTfixed;
        }
};

And this is how I create an instance:

        OptWrtTProblem problem(mD, mTfixed, mA, lambda, onesmr,
                OptWrtTProblem::DVector::Zero(DIM),
                OptWrtTProblem::DVector::Ones(DIM));
        cppoptlib::LbfgsbSolver<OptWrtTProblem> solver;

Am I doing something wrong? My c++ is quite poor...

Thank you very much for your time!

spinicist commented 8 years ago

Which line of your code fails the assertion? If you follow the compiler messages all the way back up the stack you should be able to locate it.

When I introduced the BoundedProblem class recently I was concerned that dynamic-sized problems wouldn’t work. Your cost function looks a bit odd (compared to how I use the library). Do you not know the number of parameters beforehand?

On 8 Jul 2016, at 15:45, nikitaved notifications@github.com wrote:

Hi!

Eigen/src/Core/CwiseNullaryOp.h failed at YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR assertion. It is a constant matrix (Eigen::DenseBase::ConstantReturnType) which is expected somewhere (getters of bounds, I presume) and its size must be known beforehand.

My problem is defined like this:

class OptWrtTProblem : public cppoptlib::BoundedProblem { public: using Superclass = cppoptlib::BoundedProblem; using DVector = TVector; using EigenMD = Map; using EigenMDconst = Map; using DVectorM = Map; using DVectorMconst = Map;

private:
    const Matrix& mD;
    const Matrix& mTfixed;
    const Matrix& mA;
    double lambda;

    const size_t m;
    const size_t n;
    const size_t r;

    /* some time savers */
    const Matrix& onesmr;
    Matrix mAAt;
    Matrix mDAt;
    Double regTfixed;
    Matrix gradRegTfixed;
    Double dotGradTfixed;

public:
    OptWrtTProblem(const Matrix& mDIn,
                   const Matrix& mTfixedIn,
                   const Matrix& mAIn,
                   double lambdaIn,
                   const Matrix& onesmrIn,
                   const DVector& lb,
                   const DVector& ub
    )
    : Superclass(lb, ub),
      mD{mDIn},
      mTfixed{mTfixedIn},
      mA{mAIn},
      lambda{lambdaIn},
      onesmr{onesmrIn},
      m(mD.rows()),
      n(mD.cols()),
      r(mA.rows())
    {
        mAAt = mA * mA.transpose();
        mDAt = mD * mA.transpose();

        regTfixed = mTfixed.sum() - mTfixed.squaredNorm();
        gradRegTfixed = onesmr - 2 * mTfixed;
        dotGradTfixed = gradRegTfixed.cwiseProduct(mTfixed).sum();
    }

    Double value(const DVector& mTVec) {
        EigenMDconst mT(mTVec.data(), m, r);
        Double loss = (mD - mT * mA).squaredNorm();
        return 0.5 * (loss + lambda * (regTfixed
                    + gradRegTfixed.cwiseProduct(mT).sum() - dotGradTfixed));
    }

    void gradient(const DVector& mTVec, DVector& gradVec) {
        EigenMDconst mT(mTVec.data(), m, r);
        EigenMD grad(gradVec.data(), m, r);
        grad = (mT * mAAt - mDAt) + 0.5 * lambda * gradRegTfixed;
    }

}; And this is how I create an instance:

    OptWrtTProblem problem(mD, mTfixed, mA, lambda, onesmr,
            OptWrtTProblem::DVector::Zero(DIM),
            OptWrtTProblem::DVector::Ones(DIM));
    cppoptlib::LbfgsbSolver<OptWrtTProblem> solver;

Am I doing something wrong? My c++ is quite poor...

Thank you very much for your time!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/PatWie/CppNumericalSolvers/issues/53, or mute the thread https://github.com/notifications/unsubscribe/AKie-Ku_Ws22gzMZ4JOkvO1HTt7k093Tks5qTmJ_gaJpZM4JIGR5.

nikitaved commented 8 years ago

I use mex. The output is quite cryptic, but the problem is within LBfgsbSolver type initiation/ constructor for sure...

Error using mex
In file included from ../eigen/Eigen/Core:254:0,
                 from ../eigen/Eigen/Dense:1,
                 from /home/nikita/tafact/tafact2_L_BFGS_B.cpp:7:
../eigen/Eigen/src/Core/CwiseNullaryOp.h: In instantiation of ‘static const ConstantReturnType Eigen::DenseBase<Derived>::Constant(const Scalar&) [with Derived = Eigen::Matrix<double,
-1, 1>; Eigen::DenseBase<Derived>::ConstantReturnType = Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1> >; typename
Eigen::internal::traits<T>::Scalar = double; Eigen::DenseBase<Derived>::Scalar = double]’:
../eigen/Eigen/src/Core/CwiseNullaryOp.h:470:28:   required from ‘static const ConstantReturnType Eigen::DenseBase<Derived>::Zero() [with Derived = Eigen::Matrix<double, -1, 1>;
Eigen::DenseBase<Derived>::ConstantReturnType = Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1> >; typename
Eigen::internal::traits<T>::Scalar = double]’
/home/nikita/tafact/../CppNumericalSolvers/include/cppoptlib/solver/lbfgsbsolver.h:228:43:   required from ‘void cppoptlib::LbfgsbSolver<TProblem>::minimize(TProblem&, typename
cppoptlib::LbfgsbSolver<TProblem>::Superclass::TVector&) [with TProblem = OptWrtTProblem; typename cppoptlib::LbfgsbSolver<TProblem>::Superclass::TVector = Eigen::Matrix<double, -1,
1>]’
/home/nikita/tafact/tafact2_L_BFGS_B.cpp:396:1:   required from here
../eigen/Eigen/src/Core/util/StaticAssert.h:32:40: error: static assertion failed: YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR
     #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
                                        ^
../eigen/Eigen/src/Core/util/StaticAssert.h:134:3: note: in expansion of macro ‘EIGEN_STATIC_ASSERT’
   EIGEN_STATIC_ASSERT(TYPE::SizeAtCompileTime!=Eigen::Dynamic, \
   ^
../eigen/Eigen/src/Core/CwiseNullaryOp.h:219:3: note: in expansion of macro ‘EIGEN_STATIC_ASSERT_FIXED_SIZE’
   EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived)
   ^

I do matrix factorizations, therefore I have to switch between vector-matrix representations to use solely matrix notation. No, I do not know the number of parameters beforehand. I plan to solve thousands of these kinds of problems.

nikitaved commented 8 years ago

Another problem with making DIM to be known in advance is that you cannot solve big problems, because they would not fit into stack memory, but this then destroys the whole purpose of L-BFGS, which was designed for huge problems where Hessian inversion is expensive and probably infeasible. It is better to stick to the all-mighty Newton's method or the Nesterov's accelerated gradient descent then...

PatWie commented 8 years ago

DIM can be Eigen::Dynamic. There should be no issue with having large vectors.

spinicist commented 8 years ago

Actually, I think I found the same bug as @nikitaved by experimenting on the nonnegls example. I have pushed an update to master.

The constructor for BoundedProblem now takes a run-time dimension argument, which can be omitted for fixed-size problems. This should allow fixed- and dynamic-sized problems to exist.

@nikitaved please take a look at the nonnegls example to see how this works and test with your code.

nikitaved commented 8 years ago

Yep, now it compiles. Thank you!