企业官网建设流程全解析

• 计算机视觉
• 科研

Spectra是 openMVG 仓库内src/third_party/spectra/目录下捆绑的 C++ 大规模特征值求解库，全称SparseEigenvalueComputationToolkit as aRedesignedARPACK。它以头文件（header-only）形式随 openMVG 分发，并在 LiGT 全局优化 等模块中被实际用于求解特征值问题。阅读本文后，你将掌握 Spectra 的设计原理、8 类求解器的选用方法、三种典型调用范式（稠密对称 / 稀疏一般 / 自定义矩阵运算）、特征值选择规则与 shift-and-invert 模式，并理解它是如何被集成进 openMVG 的。

Spectra 是什么：构建在 Eigen 之上的 C++ 特征值工具包

根据仓库内 Spectra 官方总览 与 README 的说明：

• 定位：面向大规模特征值问题的 C++ 库，构建于开源线性代数库 Eigen 之上；
• 形态：纯 header-only 实现，唯一依赖 Eigen 同样是 header-only 库，因此可以极其轻量地嵌入任何需要计算大型矩阵特征值的 C++ 项目；
• 适用范围：当需要从大型方阵中求出少量特征值时，Spectra 通常比计算完整的谱分解（spectral decomposition）高效得多。

与 ARPACK 的关系：重设计而非克隆

ARPACK 是 FORTRAN 编写的大规模特征值求解软件。Spectra 的开发深受其启发——从全称即可看出它是 ARPACK 的 C++ 重设计（redesign）：

• Spectra 基于 ARPACK Users' Guide 描述的隐式重启 Arnoldi/Lanczos 方法（implicitly restarted Arnoldi/Lanczos method）；
• 但 Spectra不使用 ARPACK 的代码，也不是 ARPACK 的 C++ 克隆；
• 它实现了 ARPACK 的主要算法，却提供了完全不同的接口，并且不依赖 ARPACK（文档原文强调 "NOT a clone of ARPACK for C++"）。

这一设计带来的直接好处是：用户不需要链接任何 FORTRAN 运行时或 ARPACK 库，仅凭 Eigen 即可完成大规模特征值计算。

核心设计思想：只算 k 个特征值，只暴露矩阵运算

Spectra 被设计为计算大型方阵 $A$ 中指定数量（$k$）个特征值。通常 $k$ 远小于矩阵规模（$n$），因此只计算少数特征值和特征向量，一般比计算整个谱分解更高效。

其最关键的抽象是：用户不需要直接提供整个矩阵，算法只要求定义在 $A$ 上的某些运算。在基本设定下，这个运算就是矩阵-向量乘法：

$$y = Ax$$

因此，只要矩阵-向量积 $Ax$ 能被高效计算——例如 $A$ 是稀疏矩阵——Spectra 就能在大规模特征值问题上发挥威力。这也是它被命名为 "Sparse Eigenvalue Computation Toolkit" 的原因：矩阵运算被当作黑盒，稀疏性带来的加速由用户端的运算实现自然继承。

使用 Spectra 的两大步：矩阵运算类 + 求解器对象

官方文档给出了明确的两步使用流程：

1. 定义一个实现特定矩阵运算的类，例如矩阵-向量乘法 $y=Ax$，或 shift-solve 运算 $y=(A-\sigma I)^{-1}x$。Spectra 提供了大量 helper 类来快速从矩阵对象构造这类运算，例如Spectra::DenseGenMatProd、Spectra::DenseSymShiftSolve等；
2. 创建某个特征值求解器类的对象，例如面向对称矩阵的Spectra::SymEigsSolver、面向一般矩阵的Spectra::GenEigsSolver，然后调用其成员函数完成计算并取回特征值与特征向量。

8 类求解器全景

仓库src/third_party/spectra/include/Spectra/下实际存在的求解器头文件与官方文档列出的求解器一一对应：

其中SymEigsSolver与GenEigsSolver的默认模板参数分别是DenseSymMatProd<double>与DenseGenMatProd<double>（见 SymEigsSolver.h），这也印证了文档“两步走”中的默认用法。

示例一：稠密对称矩阵——SymEigsSolver 入门

以下示例取自官方文档（Overview.md），演示对称矩阵特征值求解：

#include <Eigen/Core> #include <Spectra/SymEigsSolver.h> // <Spectra/MatOp/DenseSymMatProd.h> is implicitly included #include <iostream> using namespace Spectra; int main() { // We are going to calculate the eigenvalues of M Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10); Eigen::MatrixXd M = A + A.transpose(); // Construct matrix operation object using the wrapper class DenseSymMatProd DenseSymMatProd<double> op(M); // Construct eigen solver object, requesting the largest three eigenvalues SymEigsSolver<DenseSymMatProd<double>> eigs(op, 3, 6); // Initialize and compute eigs.init(); int nconv = eigs.compute(SortRule::LargestAlge); // Retrieve results Eigen::VectorXd evalues; if(eigs.info() == CompInfo::Successful) evalues = eigs.eigenvalues(); std::cout << "Eigenvalues found:\n" << evalues << std::endl; return 0; }

参数含义与取值约束（源码级说明）

从 SymEigsSolver.h 的构造函数文档可以确认三个构造参数：

• op：矩阵运算对象，实现 $Av$ 运算。可用DenseSymMatProd/SparseSymMatProd等包装类，也可自定义（需定义Scalar类型并实现与DenseSymMatProd相同的公有成员）；
• nev：请求的特征值个数，必须满足 $1 \le nev \le n-1$；
• ncv：控制算法收敛速度的参数（Krylov 子空间维数）。通常ncv越大收敛越快，但内存占用与每轮迭代的矩阵运算量也更大。必须满足 $nev < ncv \le n$，建议取 $ncv \ge 2 \cdot nev$。

以上约束在 SymEigsBase.h 的构造函数中通过std::invalid_argument强制校验，不满足会直接抛异常。

DenseSymMatProd的底层实现（DenseSymMatProd.h）利用 Eigen 的selfadjointView<Uplo>()只读取对称矩阵的下三角（默认Eigen::Lower）完成 $y = A x$，因此即使矩阵只填了一半也能正确处理对称结构。

示例二：稀疏一般矩阵——SparseGenMatProd 与复数特征值

一般实矩阵（非对称）的特征值可能为复数，因此eigenvalues()返回Eigen::VectorXcd。稀疏矩阵通过SparseGenMatProd、SparseSymMatProd等类支持：

#include <Eigen/Core> #include <Eigen/SparseCore> #include <Spectra/GenEigsSolver.h> #include <Spectra/MatOp/SparseGenMatProd.h> #include <iostream> using namespace Spectra; int main() { // A band matrix with 1 on the main diagonal, 2 on the below-main subdiagonal, // and 3 on the above-main subdiagonal const int n = 10; Eigen::SparseMatrix<double> M(n, n); M.reserve(Eigen::VectorXi::Constant(n, 3)); for(int i = 0; i < n; i++) { M.insert(i, i) = 1.0; if(i > 0) M.insert(i - 1, i) = 3.0; if(i < n - 1) M.insert(i + 1, i) = 2.0; } // Construct matrix operation object using the wrapper class SparseGenMatProd SparseGenMatProd<double> op(M); // Construct eigen solver object, requesting the largest three eigenvalues GenEigsSolver<SparseGenMatProd<double>> eigs(op, 3, 6); // Initialize and compute eigs.init(); int nconv = eigs.compute(SortRule::LargestMagn); // Retrieve results Eigen::VectorXcd evalues; if(eigs.info() == CompInfo::Successful) evalues = eigs.eigenvalues(); std::cout << "Eigenvalues found:\n" << evalues << std::endl; return 0; }

GenEigsSolver 与 SymEigsSolver 的参数约束差异

一般矩阵的nev/ncv约束与对称情形不同（GenEigsSolver.h）：

• nev需满足 $1 \le nev \le n-2$；
• ncv需满足 $nev+2 \le ncv \le n$，建议取 $ncv \ge 2 \cdot nev + 1$。

这一差异源于一般矩阵特征值为复数时算法内部需要保留共轭对，Krylov 子空间需要更大的余量。

示例三：自定义矩阵运算类——不持有矩阵也能求解

Spectra 最灵活的特性是：只要实现矩阵运算接口，甚至不需要真正构造矩阵。下面的例子中矩阵以“对角线元素为 1..10”的隐式形式存在（Overview.md）：

#include <Eigen/Core> #include <Spectra/SymEigsSolver.h> #include <iostream> using namespace Spectra; // M = diag(1, 2, ..., 10) class MyDiagonalTen { public: using Scalar = double; // A typedef named "Scalar" is required int rows() const { return 10; } int cols() const { return 10; } // y_out = M * x_in void perform_op(const double *x_in, double *y_out) const { for(int i = 0; i < rows(); i++) { y_out[i] = x_in[i] * (i + 1); } } }; int main() { MyDiagonalTen op; SymEigsSolver<MyDiagonalTen> eigs(op, 3, 6); eigs.init(); eigs.compute(SortRule::LargestAlge); if(eigs.info() == CompInfo::Successful) { Eigen::VectorXd evalues = eigs.eigenvalues(); std::cout << "Eigenvalues found:\n" << evalues << std::endl; } return 0; }

该程序将得到(10, 9, 8)三个最大特征值（注释同样出现在 SymEigsSolver.h 的类文档中）。自定义类只需满足三个要求：

1. 提供using Scalar = ...;类型定义（元素类型）；
2. 提供rows()、cols()返回矩阵维度；
3. 提供perform_op(const Scalar* x_in, Scalar* y_out)实现核心矩阵运算。

正是这种“运算即矩阵”的抽象，使得 Spectra 可以轻松接入任何能够高效计算 $Ax$ 的领域代码——例如 openMVG 中由 LiGT_algorithm.cpp 构造的矩阵运算类。

特征值选择规则 SortRule：9 种规则与适用边界

compute()的第一个参数selection决定要提取哪部分特征值。仓库 SelectionRule.h 完整定义了 9 种规则：

底层实现机制

从源码看，排序通过 SortingTarget 的特化模板将每个特征值映射为一个“目标值”后升序排序（std::sort），例如：

• LargestMagn：目标为-abs(val)（负号是因为升序排序，最小的目标值对应最大的模）；
• LargestAlge：目标为-val；
• BothEnds：先按LargestAlge排序，再通过 argsort 交错重排为“最大、最小、次大、次小……”的顺序，保证无论nev取何值，前k个元素都是期望的集合。

若使用不兼容的规则（例如对一般矩阵使用LargestAlge），SelectionRule.h 会抛出std::invalid_argument("incompatible selection rule")异常。

求解器核心 API 与计算流程

对称系求解器的全部公有接口在基类 SymEigsBase.h 中定义（SymEigsSolver、SymEigsShiftSolver、SymGEigsSolver均继承自它），一般矩阵系对应 GenEigsBase.h。核心成员函数如下：

compute() 的四个参数

compute(SortRule selection, Index maxit, Scalar tol, SortRule sorting)（SymEigsBase.h）中：

• selection：选择规则，决定在全谱中选取哪些特征值（如最大的 k 个）；
• maxit：允许的最大迭代次数（默认 1000）；
• tol：特征值的精度参数（默认 1e-10），收敛判定阈值为tol * max(eps^(2/3), |θ|)，其中 θ 为 Ritz 值（见 SymEigsBase.h）；
• sorting：对最终结果的排序规则，仅支持LargestAlge/LargestMagn/SmallestAlge/SmallestMagn四种（SymEigsBase.h）。

计算状态 CompInfo

CompInfo.h 定义了四种状态：

典型判读模式：compute()返回值等于请求的nev时全部收敛；info()为Successful时方可安全读取结果。注意eigenvalues()只返回已收敛的特征值，未收敛部分不会混入结果。

底层算法骨架：隐式重启 Lanczos

对称系求解器内部执行“m 步 Lanczos 分解 → 计算 Ritz 对 → 重启”的循环（SymEigsBase.h）：

1. factorize_from(1, ncv, nmatop)建立 Lanczos 分解；
2. retrieve_ritzpair(selection)计算并按选择规则排序 Ritz 值/向量；
3. 检查收敛数nconv，若未达到nev则restart(nev_adj, selection)重启（隐式重启核心见 SymEigsBase.h：对H - μI做 QR 分解、压缩 H 与 V、再扩展分解）；
4. 达到收敛或maxit上限后按sorting规则排序并返回。

配套的线性代数基础设施位于 LinAlg/（Lanczos、TridiagEigen、UpperHessenbergQR 等），矩阵运算抽象位于 MatOp/。

Shift-and-invert 模式：寻找靠近 σ 的特征值

当需要找最接近某个数 $\sigma$ 的特征值时——例如求正定矩阵的最小特征值（此时 $\sigma=0$）——官方文档明确建议使用 shift-and-invert 模式。

数学原理

如果 $(\lambda, x)$ 是 $A$ 的特征对，即 $Ax = \lambda x$，则对任意 $\sigma$ 有：

$$(A-\sigma I)^{-1}x = \nu x, \quad \nu = \frac{1}{\lambda - \sigma}$$

也就是说 $(\nu, x)$ 是 $(A-\sigma I)^{-1}$ 的特征对。把矩阵运算 $Ay$ 替换为 $(A-\sigma I)^{-1}y$ 传给求解器，就能得到 $\nu$，再通过 $\lambda = \sigma + \nu^{-1}$ 还原原问题特征值。

为什么需要它

Spectra（以及 ARPACK）的算法擅长找大模特征值，但在寻找接近零的特征值时可能失效。设 $\sigma=0$：此时找 $A^{-1}$ 的最大特征值 $\nu$，对应 $A$ 的最小特征值 $\lambda$（因为 $\nu$ 最大意味着 $\lambda$ 最小）。

模式要点（源码确认）

• 在 shift-and-invert 模式下，选择规则作用于 $\nu = 1/(\lambda-\sigma)$ 而非 $\lambda$。因此LargestMagn+ 位移 $\sigma$ 找到的是 $A$ 中最接近 $\sigma$的特征值；
• 但eigenvalues()始终返回原问题的特征值 $\lambda$（而非 $\nu$）；特征向量在两种问题下相同；
• 还原逻辑在 SymEigsShiftSolver.h 的sort_ritzpair()重写中实现：m_ritz_val = 1 / m_ritz_val + m_sigma。

实际使用：SymEigsShiftSolver

#include <Eigen/Core> #include <Spectra/SymEigsShiftSolver.h> // <Spectra/MatOp/DenseSymShiftSolve.h> is implicitly included #include <iostream> using namespace Spectra; int main() { // A size-10 diagonal matrix with elements 1, 2, ..., 10 Eigen::MatrixXd M = Eigen::MatrixXd::Zero(10, 10); for (int i = 0; i < M.rows(); i++) M(i, i) = i + 1; // Construct matrix operation object using the wrapper class DenseSymShiftSolve<double> op(M); // Construct eigen solver object with shift 0 // This will find eigenvalues that are closest to 0 SymEigsShiftSolver<DenseSymShiftSolve<double>> eigs(op, 3, 6, 0.0); eigs.init(); eigs.compute(SortRule::LargestMagn); if (eigs.info() == CompInfo::Successful) { Eigen::VectorXd evalues = eigs.eigenvalues(); // Will get (3.0, 2.0, 1.0) std::cout << "Eigenvalues found:\n" << evalues << std::endl; } return 0; }

SymEigsShiftSolver的构造参数在 SymEigsShiftSolver.h 中为(op, nev, ncv, sigma)，构造函数内部会调用op.set_shift(m_sigma)把位移写入运算对象。

Shift-solve 运算类的底层实现

DenseSymShiftSolve（DenseSymShiftSolve.h）通过set_shift(sigma)对 $A - \sigma I$ 做BKLDLT 分解（带改进的 LDLT，见 LinAlg/BKLDLT.h），perform_op则调用m_solver.solve(x)完成 $(A-\sigma I)^{-1}x$。若分解失败（例如位移使矩阵奇异），set_shift会抛出std::invalid_argument异常（DenseSymShiftSolve.h）。

自定义 shift-solve 运算类

与自定义perform_op类似，shift-solve 运算类还需额外实现set_shift(Scalar sigma)方法。官方文档给出了MyDiagonalTenShiftSolve示例（Overview.md）：

// M = diag(1, 2, ..., 10) class MyDiagonalTenShiftSolve { private: double sigma_; public: using Scalar = double; // A typedef named "Scalar" is required int rows() const { return 10; } int cols() const { return 10; } void set_shift(double sigma) { sigma_ = sigma; } // y_out = inv(A - sigma * I) * x_in // inv(A - sigma * I) = diag(1/(1-sigma), 1/(2-sigma), ...) void perform_op(double *x_in, double *y_out) const { for (int i = 0; i < rows(); i++) { y_out[i] = x_in[i] / (i + 1 - sigma_); } } }; // 使用：找最接近 3.14 的三个特征值（得到 4.0, 3.0, 2.0） SymEigsShiftSolver<MyDiagonalTenShiftSolve> eigs(op, 3, 6, 3.14);

广义特征值问题：SymGEigsSolver 的两种模式

SymGEigsSolver解决 $Ax = \lambda Bx$（$A$ 对称、$B$ 正定对称）的广义特征值问题。由 SymGEigsSolver.h 的文档可知，它由模板参数Mode决定两种工作模式（枚举定义见 GEigsMode.h）：

• Cholesky 模式（GEigsMode::Cholesky）：假设 $B$ 可用 Cholesky 分解，是优先推荐模式；第二个运算对象用DenseCholesky/SparseCholesky创建；
• RegularInverse 模式（GEigsMode::RegularInverse）：要求 $Bv$ 与 $B^{-1}v$ 两种运算，仅在 Cholesky 分解难以实现、或 $B^{-1}v$ 计算远快于 Cholesky 分解时使用；第二个运算对象用SparseRegularInverse创建。

GEigsMode枚举还包含ShiftInvert、Buckling、Cayley三种模式（GEigsMode.h），供对应的广义 shift-and-invert 系列求解器（如SymGEigsShiftSolver）使用。

openMVG 中的实际集成：LiGT 的全局优化

Spectra 并非孤立捆绑的第三方库——它已被 openMVG 的核心算法实际调用。在 LiGT 全局优化实现 中：

• 第 23 行包含third_party/spectra/include/Spectra/SymEigsShiftSolver.h，第 27 行using namespace Spectra;；
• 第 262 行注释// ========================= Solve Problem by Spectra's Eigs =======================标记了特征值求解入口。

这证明了 openMVG 在 LiGT（一种用于相机全局位姿优化的方法）中，正是利用 Spectra 的SymEigsShiftSolver完成大规模特征值求解，是"以矩阵运算抽象替代整矩阵存储"设计思想的典型生产级应用。若你需要在 openMVG 其他模块中做特征值分解，可直接复用这一集成路径：包含src/third_party/spectra/include/Spectra/下对应头文件即可，无需额外安装外部依赖。

在 openMVG 中的构建与安装方式

Spectra 位于src/third_party/spectra/，其自身的 CMakeLists.txt 记录了版本与集成细节：

• 项目版本1.0.1（project (Spectra VERSION 1.0.1 LANGUAGES CXX)）；
• 作为INTERFACE 库导出（纯头文件，无编译产物），target_link_libraries(Spectra INTERFACE Eigen3::Eigen)；
• 可选构建开关：BUILD_TESTS（测试，见 test/ 下的 SymEigs.cpp、GenEigs.cpp、SymEigsShift.cpp、SparseSymMatProd.cpp 等）与BUILD_EXAMPLES（示例，见 examples/ 的 DavidsonSymEigs_example.cpp）；
• 安装后通过find_package生成Spectra::SpectraCMake target 供其他项目链接；
• 需要 Eigen 3.x 且 C++11 及以上（set(CMAKE_CXX_STANDARD 11)）。

由于是 header-only，在 openMVG 内最直接的用法就是直接包含头文件路径：#include <Spectra/SymEigsSolver.h>，并保证 Eigen 头文件在 include 路径中（openMVG 已内置 Eigen，开箱即用）。

许可证

Spectra采用MPL2（Mozilla Public License 2.0）开源协议，与 Eigen 相同。许可证文件见 LICENSE，版本变更历史见 CHANGELOG.md（1.0.0 起存在 API 破坏性变更，迁移说明见 MIGRATION.md）。

总结

Spectra 以"隐式重启 Arnoldi/Lanczos 方法"为核心算法，用 header-only 的轻量形态和"矩阵运算对象 + 求解器对象"的两段式接口，把大规模特征值计算的门槛降到了仅依赖 Eigen 的程度。在 openMVG 中，它不仅是捆绑依赖，更是 LiGT 全局优化等模块的运行时引擎。掌握本文的 8 类求解器选型、9 种SortRule选择规则、nev/ncv参数约束与 shift-and-invert 变换，即可在 openMVG 及你自己的 C++ 项目中高效复用这套能力。

• 计算机视觉
• 科研