c++ - 如何用Ceres解决大规模非线性优化问题？

Question

我需要优化由点的二维网格表示的表面，以生成与提供的目标法 vector 对齐的表面法 vector 。网格大小可能在 201x201 和 1001x1001 之间。这意味着变量的数量将为 40,000 到 1,000,000，因为我只修改网格点的 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 优化大量变量(> 30,000)？提前致谢

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

我的代码看起来像这样(为了去掉不相关的 Material 而简化):

#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;
}

Answer 1

两件事

1. 您正在使用 DENSE_QR 作为线性求解器，这会产生密集的雅可比矩阵。这是一个坏主意。将您的线性求解器更改为 SPARSE_NORMAL_CHOLESKY，您应该能够很容易地解决这种规模的问题。

如果您使用的是 1.9 或更早版本，您将需要 SuiteSpare/CXSparse。如果您使用最新的候选版本或 git 版本，您应该也可以使用 Eigen 来执行稀疏线性代数。

1. 您正在为整个问题创建一个成本函数。这意味着您没有暴露任何问题的稀疏性。这是导致堆栈分配问题的原因，因为自动微分涉及堆栈上的数据。

查看 ceres 附带的示例代码，例如 denoising.cc，它对整个图像进行去噪，并具有类似的网格状结构。

更一般地说，为问题中的每个顶点创建一个残差 block 。