C++ Собственное значение / векторное разложение, нужны только первые n векторов быстро

Question

У меня есть ковариационная матрица размером ~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() удостоверится в этом, но это не проблема).

Answer 1

В этой статье Саймон Функ показывает простой и эффективный способ оценки разложения по сингулярным числам (SVD) очень большой матрицы. В его случае матрица является разреженной, с размерами: 17 000 x 500 000.

Теперь, глядя здесь, описывается, как разложение по собственным значениям тесно связано с СВД. Таким образом, вы могли бы извлечь выгоду из рассмотрения модифицированной версии подхода Саймона Фанка, особенно если ваша матрица редкая. Кроме того, ваша матрица не только квадратная, но и симметричная (если вы это подразумеваете под ковариационной), что, вероятно, приводит к дополнительному упрощению.

... просто идея:)

Answer 2

Кажется, Spectra справится со своей задачей с хорошими характеристиками.

Вот пример из их документации для вычисления трех первых собственных значений плотной симметричной матрицы M (аналогично вашей ковариационной матрице):

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