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

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


person MShekow    schedule 16.06.2012    source источник
comment
Под «первыми 30 собственными значениями / векторами» вы подразумеваете собственные значения с наибольшими модулями, наибольшими действительными частями или чем-то еще? После поиска в Google оказалось, что SLEPc может иметь то, что вы ищете.   -  person James    schedule 16.06.2012
comment
Я ищу 30 собственных векторов, которые соответствуют 30 наибольшим действительным собственным значениям, полученным в результате собственного разложения реальной симметричной матрицы.   -  person MShekow    schedule 16.06.2012
comment
Я бы использовал ARPACK для этого. Вы мгновенно получите свои 30 собственных векторов.   -  person David Heffernan    schedule 16.06.2012
comment
Какую конкретную реализацию вы рекомендуете? Мне нужно, чтобы эта штука работала в Windows (конечно, я готов скомпилировать ее сам, если это необходимо).   -  person MShekow    schedule 16.06.2012
comment
Используйте веб-поиск, чтобы узнать об ARPACK.   -  person David Heffernan    schedule 16.06.2012
comment
Да, ARPACK кажется очень модным и свежим (включая сарказм xD). Мне скорее интересно, есть ли какая-нибудь современная библиотека, которая может делать то, что мне нужно, например, eigen (я посмотрел на их API, похоже, они не поддерживают ограничение количества определяемых собственных векторов) .   -  person MShekow    schedule 16.06.2012
comment
А, я думал, тебе нужно выступление. Если вы хотите что-то блестящее, то ARPACK не для вас. Имейте в виду, этого достаточно для MATLAB. Я не знаю, как вы пришли к выводу, что ARPACK не будет извлекать подмножества собственных векторов. Вы ошиблись.   -  person David Heffernan    schedule 16.06.2012
comment
Теоретически это должно позволить значительно ускорить вычисления, верно? Это зависит от используемого алгоритма, по крайней мере, это не совсем так. Но да, есть алгоритмы, которые позволяют быстро вычислять только собственные векторы в некотором диапазоне собственных значений. Арнольди/Ланцош — самые выдающиеся, так что ARPACK в некотором роде каноничен. То, что он такой старый, не означает, что он плохой, он определенно хорош, если вам действительно нужна производительность; но я согласен, что эти библиотеки Fortran довольно ужасны в использовании.   -  person leftaroundabout    schedule 16.06.2012
comment
Ах, я думаю, вы имеете в виду, что eigen не делает то, что вы хотите, а не ARPACK. Извиняюсь. В любом случае, все, что я могу сказать, это то, что раньше я использовал прямой алгоритм, который мог занять много минут и часто потреблял всю системную память для получения каких-либо результатов. С помощью ARPACK я могу получить первые несколько сотен собственных пар за секунду или две, не требуя больших системных ресурсов.   -  person David Heffernan    schedule 16.06.2012
comment
Да, у него ужасный интерфейс, и мне пришлось сильно его настроить, чтобы сделать его потокобезопасным. Но это, безусловно, превосходно с точки зрения алгоритма. Это фантастически высокое качество. Просто немного сложно использовать. Если вы не являетесь сильным программистом, вам лучше использовать что-то более простое.   -  person David Heffernan    schedule 16.06.2012
comment
@DavidHeffernan да, я действительно имел в виду, что eigen (библиотека) не может указать нет. векторов. Учитывая, что я не эксперт по C++, я совершенно уверен, что время, необходимое мне для понимания и использования ARPACK, НАМНОГО больше, чем простое ожидание нескольких сотен декомпозиций матриц, которые мне нужно выполнить. для моих экспериментов. Тем не менее, если кто-нибудь знает более высокоуровневую (C/C++) библиотеку с моими указанными требованиями, пожалуйста, стреляйте :)   -  person MShekow    schedule 16.06.2012
comment
Я мог бы исследовать, почему OpenCV 2.3.1 НЕ говорит (в своей документации), что параметры низкого и высокого индекса игнорируются (вместо этого они объясняют их значение), а 2.4.1 делает это.   -  person MShekow    schedule 16.06.2012
comment
Оказывается, OpenCV 2.3 также не использовал эти два параметра. В любом случае, я создал новую тему, в которой пытаюсь заставить ARPACK++ работать в Windows, см. stackoverflow.com/questions/11071998/   -  person MShekow    schedule 17.06.2012
comment
Обновил исходный пост, наконец-то он работает... :)   -  person MShekow    schedule 17.06.2012
comment
@NameZero912 NameZero912, если у вас есть ответ, но никто его не предоставил, ответьте сами и примите его. Это уберет этот вопрос из списка вопросов без ответа. :)   -  person Yakk - Adam Nevraumont    schedule 24.12.2012
comment
@NameZero912: Пингую, чтобы вы выполнили предложение Якка. Пожалуйста, отредактируйте ответ в своем вопросе и опубликуйте его отдельно.   -  person GManNickG    schedule 16.01.2013
comment
Привет! Вы сравнивали ARPACK и Intel MKL ssyevx?   -  person Glenn Yu    schedule 28.01.2013


Ответы (2)


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

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

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

person arr_sea    schedule 23.01.2013

Похоже, что Spectra справится со своей задачей с хорошей производительностью.

Вот пример из их документации для вычисления 3 первых собственных значений плотной симметричной матрицы 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;
}
person Cryckx    schedule 11.08.2020