如何用Ceres求解大规模非线性优化问题

我需要优化由点的2D网格表示的曲面，以生成与所提供的目标法线向量对齐的曲面的法线向量。网格大小可能介于201x201和1001x1001之间。这意味着变量的数量将是40000到1000000，因为我只是在修改网格点的z坐标。

我使用的是Ceres框架，因为它应该擅长于大规模非线性优化问题。我已经尝试过MATLAB的fmincon，但它使用了大量的内存。我写了一个适用于小网格的目标函数（成功地实现了3x3和31x31）。然而，当我试图用大网格大小（157x200）编译代码时，我看到了下面的错误。我读到这是Eigen的一个限制。然而，当我告诉Ceres使用LAPACK而不是Eigen时，对于大矩阵，我会得到同样的错误。我试过这些线路：

options.dense_linear_algebra_library_type = ceres::LAPACK;
options.linear_solver_type = ceres::DENSE_QR;

这些命令告诉求解器使用LAPACK和DENSE_QR，如使用3x3网格的输出所示：

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

然而，当我使用大参数时，我仍然会得到Eigen的误差。

无论如何，我真的需要一些帮助。如何让Ceres优化大量变量（>30000）？提前感谢

Ceres链接：http://ceres-solver.org

链接到Eigen：http://eigen.tuxfamily.org/dox/

错误：

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<T, Size, MatrixOrArrayOptions, Alignment>::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<T, Size, _Rows, _Cols, _Options>::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<Derived>::PlainObjectBase() [with Derived = Eigen::Matrix<double, 31400, 1, 2, 31400, 1>]’
/usr/include/eigen3/Eigen/src/Core/Matrix.h:203:41:   required from ‘Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::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<T, N>::Jet() [with T = double; int N = 31400]’
/usr/local/include/ceres/internal/fixed_array.h:138:10:   required from ‘ceres::internal::FixedArray<T, inline_elements>::FixedArray(ceres::internal::FixedArray<T, inline_elements>::size_type) [with T = ceres::Jet<double, 31400>; long int inline_elements = 0l; ceres::internal::FixedArray<T, inline_elements>::size_type = long unsigned int]’
/usr/local/include/ceres/internal/autodiff.h:233:70:   required from ‘static bool ceres::internal::AutoDiff<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Differentiate(const Functor&, 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<CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::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<false>’
 EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 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<T, Size, MatrixOrArrayOptions, 16>::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<T, Size, _Rows, _Cols, _Options>::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<Derived>::PlainObjectBase(Eigen::PlainObjectBase<Derived>::Index, Eigen::PlainObjectBase<Derived>::Index, Eigen::PlainObjectBase<Derived>::Index) [with Derived = Eigen::Matrix<double, 1, 31400, 1, 1, 31400>; Eigen::PlainObjectBase<Derived>::Index = long int]’
/usr/include/eigen3/Eigen/src/Core/Matrix.h:273:76:   required from ‘Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Matrix(const Eigen::MatrixBase<OtherDerived>&) [with OtherDerived = Eigen::Transpose<const Eigen::Matrix<double, 31400, 1, 2, 31400, 1> >; _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<Derived>::EvalReturnType Eigen::DenseBase<Derived>::eval() const [with Derived = Eigen::Transpose<const Eigen::Matrix<double, 31400, 1, 2, 31400, 1> >; Eigen::DenseBase<Derived>::EvalReturnType = const Eigen::Matrix<double, 1, 31400, 1, 1, 31400>]’
/usr/include/eigen3/Eigen/src/Core/IO.h:244:69:   required from ‘std::ostream& Eigen::operator<<(std::ostream&, const Eigen::DenseBase<Derived>&) [with Derived = Eigen::Transpose<const Eigen::Matrix<double, 31400, 1, 2, 31400, 1> >; std::ostream = std::basic_ostream<char>]’
/usr/local/include/ceres/jet.h:632:35:   required from ‘std::ostream& ceres::operator<<(std::ostream&, const ceres::Jet<T, N>&) [with T = double; int N = 31400; std::ostream = std::basic_ostream<char>]’
/home/ubuntu/code/surfaceopt/surfaceopt.cc:103:50:   required from ‘bool ComputeEint::operator()(const T*, T*) const [with T = ceres::Jet<double, 31400>]’
/usr/local/include/ceres/internal/variadic_evaluate.h:175:26:   required from ‘static bool ceres::internal::VariadicEvaluate<Functor, T, N0, 0, 0, 0, 0, 0, 0, 0, 0, 0>::Call(const Functor&, const T* const*, T*) [with Functor = ComputeEint; T = ceres::Jet<double, 31400>; int N0 = 31400]’
/usr/local/include/ceres/internal/autodiff.h:283:45:   required from ‘static bool ceres::internal::AutoDiff<Functor, T, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::Differentiate(const Functor&, 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<CostFunctor, kNumResiduals, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>::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<false>’
 EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 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

我的代码是这样的（缩写为去掉不相关的材料）：

#define TEXT true
#define VERBOSE false 
#define NV 31400
#define NF 62088
#define NX 157
#define NY 200
#define MAXNB 6
#include <math.h>
#include <ceres/ceres.h>
#include <ceres/rotation.h>
#include "glog/logging.h"
#include <iostream>
#include <fstream>
#include <iterator>
#include <algorithm>
#include <string>
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(&XY[0][0], &XY[0][0]+NV*2,&XY_[0][0]);
...
template <typename T>
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<ComputeEint, 1, NV>(new ComputeEint(XY, Tc, AF, Ntgt, num_AF));
  std::cout << "Created cost function" << "n";
  problem.AddResidualBlock(cost_function, NULL, &z[0]);
  std::cout << "Added residual block" << "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, &problem, &summary);
  std::cout << summary.FullReport() << "n";
  //Write output of optimization to file
  ...
  return 0;
}