C++特征值/向量分解，只需要快速的前n个向量

Posted 2023-02-16

【中文标题】C++特征值/向量分解，只需要快速的前n个向量

【英文标题】：C++ eigenvalue/vector decomposition, only need first n vectors fast

【发布时间】：2012-06-16 15:17:02

【问题描述】：

我有一个约 3000x3000 的协方差矩阵，我在该矩阵上计算特征值-特征向量分解（它是一个 OpenCV 矩阵，我使用 cv::eigen() 来完成工作）。

但是，我实际上只需要前 30 个特征值/向量，其余的我不关心。从理论上讲，这应该可以显着加快计算速度，对吧？我的意思是，这意味着它需要计算的特征向量少了 2970 个。

哪个 C++ 库允许我这样做？请注意，OpenCV 的 eigen() 方法确实有参数，但是文档说它们被忽略了，我自己测试过，它们确实被忽略了：D

更新： 我设法用 ARPACK 做到了。我设法为 Windows 编译它，甚至使用它。结果看起来很有希望，可以在这个玩具示例中看到一个插图：

#include "ardsmat.h"
#include "ardssym.h"
int     n = 3;           // Dimension of the problem.
    double* EigVal = NULL;  // Eigenvalues.
    double* EigVec = NULL; // Eigenvectors stored sequentially.


    int lowerHalfElementCount = (n*n+n) / 2;
    //whole matrix:
    /*
    2  3  8
    3  9  -7
    8  -7 19
    */
    double* lower = new double[lowerHalfElementCount]; //lower half of the matrix
    //to be filled with COLUMN major (i.e. one column after the other, always starting from the diagonal element)
    lower[0] = 2; lower[1] = 3; lower[2] = 8; lower[3] = 9; lower[4] = -7; lower[5] = 19;
    //params: dimensions (i.e. width/height), array with values of the lower or upper half (sequentially, row major), 'L' or 'U' for upper or lower
    ARdsSymMatrix<double> mat(n, lower, 'L');

    // Defining the eigenvalue problem.
    int noOfEigVecValues = 2;
    //int maxIterations = 50000000;
    //ARluSymStdEig<double> dprob(noOfEigVecValues, mat, "LM", 0, 0.5, maxIterations);
    ARluSymStdEig<double> dprob(noOfEigVecValues, mat);

    // Finding eigenvalues and eigenvectors.

    int converged = dprob.EigenValVectors(EigVec, EigVal);
    for (int eigValIdx = 0; eigValIdx < noOfEigVecValues; eigValIdx++) 
        std::cout << "Eigenvalue: " << EigVal[eigValIdx] << "\nEigenvector: ";

        for (int i = 0; i < n; i++) 
            int idx = n*eigValIdx+i;
            std::cout << EigVec[idx] << " ";
        
        std::cout << std::endl;
    

结果是：

9.4298, 24.24059

对于特征值，并且

-0.523207, -0.83446237, -0.17299346
0.273269, -0.356554, 0.893416

分别为 2 个特征向量（每行一个特征向量） 代码找不到 3 个特征向量（在这种情况下它只能找到 1-2 个，assert() 可以确保这一点，但这不是问题）。

【问题讨论】：

“前 30 个特征值/向量”是指具有最大模数、最大实部还是其他的特征值？谷歌搜索后，看起来SLEPc 可能有你要找的东西。

我会为此使用 ARPACK。您将立即获得 30 个特征向量。

“理论上，这应该可以显着加快计算速度，对吧？”这取决于所使用的算法，至少并非如此。但是，是的，有些算法只允许快速计算某些特征值范围内的特征向量。 Arnoldi/Lanczos 是最突出的，所以 ARPACK 是一种规范。它这么老并不意味着它很糟糕，如果你真的想要性能，它肯定很棒；但我同意这些 Fortran 库使用起来相当糟糕。

啊，我想你的意思是 eigen 不会做你想做的事，而不是 ARPACK。对不起。无论如何，我只能说我以前使用的直接算法可能需要几分钟，并且经常消耗所有系统内存来产生任何结果。使用 ARPACK，我可以在一两秒内获得前几百个特征对，而无需大量系统资源。

@NameZero912，如果您有答案但没有人提供，请自己回答并接受。它会将这个问题从未回答的问题列表中删除。 :)

【参考方案1】：

在this 文章中，Simon Funk 展示了一种简单有效的方法来估计一个非常大的矩阵的奇异值分解 (SVD)。在他的例子中，矩阵是稀疏的，尺寸为：17,000 x 500,000。

现在，看here，描述了特征值分解如何与 SVD 密切相关。因此，考虑 Simon Funk 方法的修改版本，您可能会受益，特别是如果您的矩阵是稀疏的。此外，您的矩阵不仅是方形的，而且是对称的（如果这就是您所说的协方差），这可能会导致额外的简化。

...只是一个想法:)

【参考方案2】：

看来Spectra 会表现出色。

这是他们文档中的一个示例，用于计算密集对称矩阵 M（同样是您的协方差矩阵）的 3 个第一个特征值：

#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< double, LARGEST_ALGE, DenseSymMatProd<double> > eigs(&op, 3, 6);
    // Initialize and compute
    eigs.init();
    int nconv = eigs.compute();
    // Retrieve results
    Eigen::VectorXd evalues;
    if(eigs.info() == SUCCESSFUL)
        evalues = eigs.eigenvalues();
    std::cout << "Eigenvalues found:\n" << evalues << std::endl;
    return 0;