title: How to solve large-scale nonlinear optimization problems with Ceres?
body: <p>I need to optimize a surface represented by a 2D grid of points to produce normal vectors of the surface that align with provided target normal vectors. The grid size is likely to be between 201x201 and 1001x1001. That means that the number of variables will be 40,000 to 1,000,000, as I am only modifying the z-coordinates of the mesh points. </p>
  
  <p>I am using the Ceres framework, as it is supposed to excel at large-scale nonlinear optimization problems. I already tried MATLAB's fmincon, but it uses an incredible amount of memory. I wrote an objective function that works for small meshes (successful at 3x3 and 31x31). However, when I try to compile the code with a large mesh size (157x200), I see the error below. I have read that this is a limitation of Eigen. However, when I tell Ceres to use LAPACK instead of Eigen, I get the same error for large matrices. I tried these lines:</p>
  
  <pre><code>options.dense_linear_algebra_library_type = ceres::LAPACK;
  
  options.linear_solver_type = ceres::DENSE_QR;
  </code></pre>
  
  <p>These tell the solver to use LAPACK and DENSE_QR, as the output using a 3x3 mesh shows:</p>
  
  <pre><code>Minimizer                        TRUST_REGION
  
  Dense linear algebra library           LAPACK
  Trust region strategy     LEVENBERG_MARQUARDT
  
                                      Given                     Used
  Linear solver                        DENSE_QR                 DENSE_QR
  Threads                                     1                        1
  Linear solver threads                       1                        1
  </code></pre>
  
  <p>However, when I use large parameters, I still get the errors for Eigen. </p>
  
  <p>Anyways, I could really use some help with this. How can I get Ceres to optimize a large number of variables (> 30,000)? Thanks in advance</p>
  
  <p>Link to Ceres: <a href="http://ceres-solver.org" rel="nofollow">http://ceres-solver.org</a></p>
  
  <p>Link to Eigen: <a href="http://eigen.tuxfamily.org/dox/" rel="nofollow">http://eigen.tuxfamily.org/dox/</a></p>
  
  <p>Error:</p>
  
  <pre><code>In file included from /usr/include/eigen3/Eigen/Core:254:0,
               from /usr/local/include/ceres/jet.h:165,
               from /usr/local/include/ceres/internal/autodiff.h:145,
               from /usr/local/include/ceres/autodiff_cost_function.h:132,
               from /usr/local/include/ceres/ceres.h:37,
               from /home/ubuntu/code/surfaceopt/surfaceopt.cc:10:
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h: In instantiation of ‘Eigen::internal::plain_array&lt;T, Size, MatrixOrArrayOptions, Alignment&gt;::plain_array()     [with T = double; int Size = 31400; int MatrixOrArrayOptions = 2; int Alignment = 0]’:
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:117:27:   required from     ‘Eigen::DenseStorage&lt;T, Size, _Rows, _Cols, _Options&gt;::DenseStorage() [with T = double; int Size = 31400; int _Rows = 31400; int _Cols = 1; int _Options = 2]’
  /usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:421:55:   required from ‘Eigen::PlainObjectBase&lt;Derived&gt;::PlainObjectBase() [with Derived = Eigen::Matrix&lt;double, 31400, 1, 2, 31400, 1&gt;]’
  /usr/include/eigen3/Eigen/src/Core/Matrix.h:203:41:   required from ‘Eigen::Matrix&lt;_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols&gt;::Matrix() [with _Scalar = double; int _Rows = 31400; int _Cols = 1; int _Options = 2; int _MaxRows = 31400; int _MaxCols = 1]’
  /usr/local/include/ceres/jet.h:179:13:   required from ‘ceres::Jet&lt;T, N&gt;::Jet() [with T = double; int N = 31400]’
  /usr/local/include/ceres/internal/fixed_array.h:138:10:   required from ‘ceres::internal::FixedArray&lt;T, inline_elements&gt;::FixedArray(ceres::internal::FixedArray&lt;T, inline_elements&gt;::size_type) [with T = ceres::Jet&lt;double, 31400&gt;; long int inline_elements = 0l; ceres::internal::FixedArray&lt;T, inline_elements&gt;::size_type = long unsigned int]’
  /usr/local/include/ceres/internal/autodiff.h:233:70:   required from ‘static bool ceres::internal::AutoDiff&lt;Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9&gt;::Differentiate(const Functor&amp;, const T* const*, int, T*, T**) [with Functor = ComputeEint; T = double; int N0 = 31400; int N1 = 0; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
  /usr/local/include/ceres/autodiff_cost_function.h:218:25:   required from ‘bool ceres::AutoDiffCostFunction&lt;CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9&gt;::Evaluate(const double* const*, double*, double**) const [with CostFunctor = ComputeEint; int kNumResiduals = 1; int N0 = 31400; int N1 = 0; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
  /home/ubuntu/code/surfaceopt/surfaceopt.cc:367:1:   required from here
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:41:5: error: ‘OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG’ is not a member of ‘Eigen::internal::static_assertion&lt;false&gt;’
   EIGEN_STATIC_ASSERT(Size * sizeof(T) &lt;= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG);
   ^
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h: In instantiation of ‘Eigen::internal::plain_array&lt;T, Size, MatrixOrArrayOptions, 16&gt;::plain_array() [with T = double; int Size = 31400; int MatrixOrArrayOptions = 1]’:
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:120:59:   required from ‘Eigen::DenseStorage&lt;T, Size, _Rows, _Cols, _Options&gt;::DenseStorage(Eigen::DenseIndex, Eigen::DenseIndex, Eigen::DenseIndex) [with T = double; int Size = 31400; int _Rows = 1; int _Cols = 31400; int _Options = 1; Eigen::DenseIndex = long int]’
  /usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:438:41:   required from ‘Eigen::PlainObjectBase&lt;Derived&gt;::PlainObjectBase(Eigen::PlainObjectBase&lt;Derived&gt;::Index, Eigen::PlainObjectBase&lt;Derived&gt;::Index, Eigen::PlainObjectBase&lt;Derived&gt;::Index) [with Derived = Eigen::Matrix&lt;double, 1, 31400, 1, 1, 31400&gt;; Eigen::PlainObjectBase&lt;Derived&gt;::Index = long int]’
  /usr/include/eigen3/Eigen/src/Core/Matrix.h:273:76:   required from ‘Eigen::Matrix&lt;_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols&gt;::Matrix(const Eigen::MatrixBase&lt;OtherDerived&gt;&amp;) [with OtherDerived = Eigen::Transpose&lt;const Eigen::Matrix&lt;double, 31400, 1, 2, 31400, 1&gt; &gt;; _Scalar = double; int _Rows = 1; int _Cols = 31400; int _Options = 1; int _MaxRows = 1; int _MaxCols = 31400]’
  /usr/include/eigen3/Eigen/src/Core/DenseBase.h:367:62:   required from ‘Eigen::DenseBase&lt;Derived&gt;::EvalReturnType Eigen::DenseBase&lt;Derived&gt;::eval() const [with Derived = Eigen::Transpose&lt;const Eigen::Matrix&lt;double, 31400, 1, 2, 31400, 1&gt; &gt;; Eigen::DenseBase&lt;Derived&gt;::EvalReturnType = const Eigen::Matrix&lt;double, 1, 31400, 1, 1, 31400&gt;]’
  /usr/include/eigen3/Eigen/src/Core/IO.h:244:69:   required from ‘std::ostream&amp; Eigen::operator&lt;&lt;(std::ostream&amp;, const Eigen::DenseBase&lt;Derived&gt;&amp;) [with Derived = Eigen::Transpose&lt;const Eigen::Matrix&lt;double, 31400, 1, 2, 31400, 1&gt; &gt;; std::ostream = std::basic_ostream&lt;char&gt;]’
  /usr/local/include/ceres/jet.h:632:35:   required from ‘std::ostream&amp; ceres::operator&lt;&lt;(std::ostream&amp;, const ceres::Jet&lt;T, N&gt;&amp;) [with T = double; int N = 31400; std::ostream = std::basic_ostream&lt;char&gt;]’
  /home/ubuntu/code/surfaceopt/surfaceopt.cc:103:50:   required from ‘bool ComputeEint::operator()(const T*, T*) const [with T = ceres::Jet&lt;double, 31400&gt;]’
  /usr/local/include/ceres/internal/variadic_evaluate.h:175:26:   required from ‘static bool ceres::internal::VariadicEvaluate&lt;Functor, T, N0, 0, 0, 0, 0, 0, 0, 0, 0, 0&gt;::Call(const Functor&amp;, const T* const*, T*) [with Functor = ComputeEint; T = ceres::Jet&lt;double, 31400&gt;; int N0 = 31400]’
  /usr/local/include/ceres/internal/autodiff.h:283:45:   required from ‘static bool ceres::internal::AutoDiff&lt;Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9&gt;::Differentiate(const Functor&amp;, const T* const*, int, T*, T**) [with Functor = ComputeEint; T = double; int N0 = 31400; int N1 = 0; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
  /usr/local/include/ceres/autodiff_cost_function.h:218:25:   required from ‘bool ceres::AutoDiffCostFunction&lt;CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9&gt;::Evaluate(const double* const*, double*, double**) const [with CostFunctor = ComputeEint; int kNumResiduals = 1; int N0 = 31400; int N1 = 0; int N2 = 0; int N3 = 0; int N4 = 0; int N5 = 0; int N6 = 0; int N7 = 0; int N8 = 0; int N9 = 0]’
  /home/ubuntu/code/surfaceopt/surfaceopt.cc:367:1:   required from here
  /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:79:5: error: ‘OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG’ is not a member of ‘Eigen::internal::static_assertion&lt;false&gt;’
   EIGEN_STATIC_ASSERT(Size * sizeof(T) &lt;= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG);
   ^
  make[2]: *** [CMakeFiles/surfaceopt.dir/surfaceopt.cc.o] Error 1
  make[1]: *** [CMakeFiles/surfaceopt.dir/all] Error 2
  make: *** [all] Error 2
  </code></pre>
  
  <p>My code looks like this (abbreviated to take out irrelevant material):</p>
  
  <pre><code>#define TEXT true
  #define VERBOSE false 
  #define NV 31400
  #define NF 62088
  #define NX 157
  #define NY 200
  #define MAXNB 6
  
  #include &lt;math.h&gt;
  #include &lt;ceres/ceres.h&gt;
  #include &lt;ceres/rotation.h&gt;
  #include "glog/logging.h"
  #include &lt;iostream&gt;
  #include &lt;fstream&gt;
  #include &lt;iterator&gt;
  #include &lt;algorithm&gt;
  #include &lt;string&gt;
  
  using ceres::AutoDiffCostFunction;
  using ceres::CostFunction;
  using ceres::Problem;
  using ceres::Solver;
  using ceres::Solve;
  using ceres::CrossProduct;
  ...
  
  class ComputeEint {
  
  private:
    double XY_ [NV][2];         // X and Y coords
    int C_ [NF][3];             // Connectivity list
    int AF_ [NV][MAXNB];        // List of adjacent faces to each vertex
    double Ntgt_ [NV][3];       // Target normal vectors
    int num_AF_ [NV];           // Number of adjacent faces to each vertex
  public:
  
  //Constructor
  ComputeEint(double XY[][2], int C[][3], int AF[][MAXNB], double Ntgt[][3], int num_AF[NV]) {
  
  std::copy(&amp;XY[0][0], &amp;XY[0][0]+NV*2,&amp;XY_[0][0]);
  ...
  
  template &lt;typename T&gt;
  bool operator()(const T* const z, T* e) const {
    e[0] = T(0);
    ... 
    //Computes vertex normals for triangulated surface by averaging adjacent face normals
    ...
    e[0] = e[0]/T(NV);
    return true;
    }
  };
  
  int main(int argc, char** argv) {
    google::InitGoogleLogging(argv[0]);
  
    double Tp [NV][3];            //Points in mesh
    int    Tc [NF][3];            //Mesh connectivity list
    double Ntgt [NV][3];          //Target normals
    int    AF [NV][MAXNB];        //List of adjacent faces of each vertex
    int    num_AF [NV];           //Number of adjacent faces for each vertex
  
    int nx = NX;
    int ny = NY;
  
    //Read Tp, Tc, Ntgt, AF, num_AF from file
    ...
    // Set up XY for cost functor
    double XY [NV][2];
    double z [NV];
    //Copy first two columns of Tp into XY
    Problem problem;
  
    // Set up the only cost function (also known as residual). This uses
    // auto-differentiation to obtain the derivative (jacobian).
    CostFunction* cost_function =
    new AutoDiffCostFunction&lt;ComputeEint, 1, NV&gt;(new ComputeEint(XY, Tc, AF, Ntgt, num_AF));
  
    std::cout &lt;&lt; "Created cost function" &lt;&lt; "\n";
    problem.AddResidualBlock(cost_function, NULL, &amp;z[0]);
  
    std::cout &lt;&lt; "Added residual block" &lt;&lt; "\n";
  
    // Run the solver!
    Solver::Options options;
    options.minimizer_progress_to_stdout = true;
    options.max_num_iterations = 50;
    options.function_tolerance = 1e-4;
    options.dense_linear_algebra_library_type = ceres::LAPACK;
    Solver::Summary summary;
    Solve(options, &amp;problem, &amp;summary);
  
    std::cout &lt;&lt; summary.FullReport() &lt;&lt; "\n";
  
    //Write output of optimization to file
    ...
    return 0;
  }
  </code></pre>