计算机视觉科研【免费下载链接】openMVGopen Multiple View Geometry library. Basis for 3D computer vision and Structure from Motion.项目地址https://gitcode.com/gh_mirrors/op/openMVG点击查看免费下载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 重设计redesignSpectra 基于 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 的两大步矩阵运算类 求解器对象官方文档给出了明确的两步使用流程定义一个实现特定矩阵运算的类例如矩阵-向量乘法 $yAx$或 shift-solve 运算 $y(A-\sigma I)^{-1}x$。Spectra 提供了大量 helper 类来快速从矩阵对象构造这类运算例如Spectra::DenseGenMatProd、Spectra::DenseSymShiftSolve等创建某个特征值求解器类的对象例如面向对称矩阵的Spectra::SymEigsSolver、面向一般矩阵的Spectra::GenEigsSolver然后调用其成员函数完成计算并取回特征值与特征向量。8 类求解器全景仓库src/third_party/spectra/include/Spectra/下实际存在的求解器头文件与官方文档列出的求解器一一对应求解器类适用问题说明SymEigsSolver实对称矩阵 $Ax\lambda x$基础模式见 SymEigsSolver.hGenEigsSolver一般实矩阵 $Ax\lambda x$特征值/特征向量可为复数见 GenEigsSolver.hSymEigsShiftSolver实对称矩阵shift-and-invert 模式找最接近 $\sigma$ 的特征值见 SymEigsShiftSolver.hGenEigsRealShiftSolver一般实矩阵shift-and-invert 模式实数位移见 GenEigsRealShiftSolver.hGenEigsComplexShiftSolver一般实矩阵shift-and-invert 模式复数位移见 GenEigsComplexShiftSolver.hSymGEigsSolver广义特征值问题 $Ax\lambda Bx$实对称支持 Cholesky / RegularInverse 两种模式见 SymGEigsSolver.hSymGEigsShiftSolver广义特征值问题实对称shift-and-invert 模式见 SymGEigsShiftSolver.hDavidsonSymEigsSolver实对称矩阵Jacobi-Davidson 算法DPR 校正见 DavidsonSymEigsSolver.h其中SymEigsSolver与GenEigsSolver的默认模板参数分别是DenseSymMatProddouble与DenseGenMatProddouble见 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 DenseSymMatProddouble op(M); // Construct eigen solver object, requesting the largest three eigenvalues SymEigsSolverDenseSymMatProddouble 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 的selfadjointViewUplo()只读取对称矩阵的下三角默认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::SparseMatrixdouble 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 SparseGenMatProddouble op(M); // Construct eigen solver object, requesting the largest three eigenvalues GenEigsSolverSparseGenMatProddouble 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.hnev需满足 $1 \le nev \le n-2$ncv需满足 $nev2 \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; SymEigsSolverMyDiagonalTen 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 的类文档中。自定义类只需满足三个要求提供using Scalar ...;类型定义元素类型提供rows()、cols()返回矩阵维度提供perform_op(const Scalar* x_in, Scalar* y_out)实现核心矩阵运算。正是这种“运算即矩阵”的抽象使得 Spectra 可以轻松接入任何能够高效计算 $Ax$ 的领域代码——例如 openMVG 中由 LiGT_algorithm.cpp 构造的矩阵运算类。特征值选择规则 SortRule9 种规则与适用边界compute()的第一个参数selection决定要提取哪部分特征值。仓库 SelectionRule.h 完整定义了 9 种规则SortRule 枚举值含义适用求解器LargestMagn模绝对值/复数范数最大的特征值对称与一般求解器LargestReal实部最大的特征值仅一般求解器LargestImag虚部按模最大的特征值仅一般求解器LargestAlge代数值最大的特征值考虑负号仅对称求解器SmallestMagn模最小的特征值对称与一般求解器SmallestReal实部最小的特征值仅一般求解器SmallestImag虚部按模最小的特征值仅一般求解器SmallestAlge代数值最小的特征值仅对称求解器BothEnds谱的两端各取一半nev为奇数时高端多取一个仅对称求解器底层实现机制从源码看排序通过 SortingTarget 的特化模板将每个特征值映射为一个“目标值”后升序排序std::sort例如LargestMagn目标为-abs(val)负号是因为升序排序最小的目标值对应最大的模LargestAlge目标为-valBothEnds先按LargestAlge排序再通过 argsort 交错重排为“最大、最小、次大、次小……”的顺序保证无论nev取何值前k个元素都是期望的集合。若使用不兼容的规则例如对一般矩阵使用LargestAlgeSelectionRule.h 会抛出std::invalid_argument(incompatible selection rule)异常。求解器核心 API 与计算流程对称系求解器的全部公有接口在基类 SymEigsBase.h 中定义SymEigsSolver、SymEigsShiftSolver、SymGEigsSolver均继承自它一般矩阵系对应 GenEigsBase.h。核心成员函数如下成员函数作用init(const Scalar* init_resid)用用户提供的初始残差向量初始化init()用随机初始残差向量初始化元素服从独立的 Uniform(-0.5, 0.5) 分布固定随机种子见 SymEigsBase.hcompute(selection, maxit, tol, sorting)执行主要计算返回收敛的特征值个数默认参数为maxit1000、tol1e-10、sortingSortRule::LargestAlgeinfo()返回计算状态CompInfonum_iterations()返回迭代次数num_operations()返回调用的矩阵运算次数eigenvalues()返回已收敛的特征值向量eigenvectors(nvec)/eigenvectors()返回已收敛的特征向量矩阵按列排列compute() 的四个参数compute(SortRule selection, Index maxit, Scalar tol, SortRule sorting)SymEigsBase.h中selection选择规则决定在全谱中选取哪些特征值如最大的 k 个maxit允许的最大迭代次数默认 1000tol特征值的精度参数默认 1e-10收敛判定阈值为tol * max(eps^(2/3), |θ|)其中 θ 为 Ritz 值见 SymEigsBase.hsorting对最终结果的排序规则仅支持LargestAlge/LargestMagn/SmallestAlge/SmallestMagn四种SymEigsBase.h。计算状态 CompInfoCompInfo.h 定义了四种状态枚举值含义Successful计算成功NotComputed尚未调用compute()NotConverging部分特征值未收敛compute()会返回已收敛个数NumericalIssue数值问题如 Cholesky 分解遇到非正定矩阵典型判读模式compute()返回值等于请求的nev时全部收敛info()为Successful时方可安全读取结果。注意eigenvalues()只返回已收敛的特征值未收敛部分不会混入结果。底层算法骨架隐式重启 Lanczos对称系求解器内部执行“m 步 Lanczos 分解 → 计算 Ritz 对 → 重启”的循环SymEigsBase.hfactorize_from(1, ncv, nmatop)建立 Lanczos 分解retrieve_ritzpair(selection)计算并按选择规则排序 Ritz 值/向量检查收敛数nconv若未达到nev则restart(nev_adj, selection)重启隐式重启核心见 SymEigsBase.h对H - μI做 QR 分解、压缩 H 与 V、再扩展分解达到收敛或maxit上限后按sorting规则排序并返回。配套的线性代数基础设施位于 LinAlg/Lanczos、TridiagEigen、UpperHessenbergQR 等矩阵运算抽象位于 MatOp/。Shift-and-invert 模式寻找靠近 σ 的特征值当需要找最接近某个数 $\sigma$ 的特征值时——例如求正定矩阵的最小特征值此时 $\sigma0$——官方文档明确建议使用 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的算法擅长找大模特征值但在寻找接近零的特征值时可能失效。设 $\sigma0$此时找 $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 DenseSymShiftSolvedouble op(M); // Construct eigen solver object with shift 0 // This will find eigenvalues that are closest to 0 SymEigsShiftSolverDenseSymShiftSolvedouble 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 运算类的底层实现DenseSymShiftSolveDenseSymShiftSolve.h通过set_shift(sigma)对 $A - \sigma I$ 做BKLDLT 分解带改进的 LDLT见 LinAlg/BKLDLT.hperform_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 SymEigsShiftSolverMyDiagonalTenShiftSolve eigs(op, 3, 6, 3.14);广义特征值问题SymGEigsSolver 的两种模式SymGEigsSolver解决 $Ax \lambda Bx$$A$ 对称、$B$ 正定对称的广义特征值问题。由 SymGEigsSolver.h 的文档可知它由模板参数Mode决定两种工作模式枚举定义见 GEigsMode.hCholesky 模式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 Spectras Eigs 标记了特征值求解入口。这证明了 openMVG 在 LiGT一种用于相机全局位姿优化的方法中正是利用 Spectra 的SymEigsShiftSolver完成大规模特征值求解是以矩阵运算抽象替代整矩阵存储设计思想的典型生产级应用。若你需要在 openMVG 其他模块中做特征值分解可直接复用这一集成路径包含src/third_party/spectra/include/Spectra/下对应头文件即可无需额外安装外部依赖。在 openMVG 中的构建与安装方式Spectra 位于src/third_party/spectra/其自身的 CMakeLists.txt 记录了版本与集成细节项目版本1.0.1project (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 且 C11 及以上set(CMAKE_CXX_STANDARD 11)。由于是 header-only在 openMVG 内最直接的用法就是直接包含头文件路径#include Spectra/SymEigsSolver.h并保证 Eigen 头文件在 include 路径中openMVG 已内置 Eigen开箱即用。许可证Spectra采用MPL2Mozilla Public License 2.0开源协议与 Eigen 相同。许可证文件见 LICENSE版本变更历史见 CHANGELOG.md1.0.0 起存在 API 破坏性变更迁移说明见 MIGRATION.md。总结Spectra 以隐式重启 Arnoldi/Lanczos 方法为核心算法用 header-only 的轻量形态和矩阵运算对象 求解器对象的两段式接口把大规模特征值计算的门槛降到了仅依赖 Eigen 的程度。在 openMVG 中它不仅是捆绑依赖更是 LiGT 全局优化等模块的运行时引擎。掌握本文的 8 类求解器选型、9 种SortRule选择规则、nev/ncv参数约束与 shift-and-invert 变换即可在 openMVG 及你自己的 C 项目中高效复用这套能力。赞分享计算机视觉科研【免费下载链接】openMVGopen Multiple View Geometry library. Basis for 3D computer vision and Structure from Motion.项目地址https://gitcode.com/gh_mirrors/op/openMVG点击查看免费下载相关推荐OneUptime Host Monitor 完全指南用 OpenTelemetry 主机指标构建服务器监控与告警OneUptime Host Monitor 完全指南用 OpenTelemetry 主机指标构建服务器监控与告警 本篇技术指南围绕 OneUptime 的计算机视觉科研SciPy 稀疏特征值问题教程用 ARPACK 的 eigs/eigsh 高效求解大规模特征值SciPy 稀疏特征值问题教程用 ARPACK 的 eigs/eigsh 高效求解大规模特征值 导读 本文深入讲解 SciPy 中基于 ARPACK 的大规模科学计算数据科学高性能计算【亲测免费】 探索Spectra大规模稀疏矩阵的高效特征值计算库探索Spectra大规模稀疏矩阵的高效特征值计算库 如果你在寻找一个可以处理大型稀疏矩阵并计算其特征值的C库那么Spectra绝对值得你关注。这个基于E上一篇5大核心功能3种使用场景开源IPTV播放器IPTVnator完整指南下一篇如何给 KernelSU 装上 meta-overlayfs 元模块让模块真的改得动 /system创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考