低 RAM 消耗 c++ 特征求解器

Question

我是 C++ 编程的新手，但我的任务是计算特征值和特征向量（标准特征问题 Ax=lx）对于对称矩阵（和 hermitian））对于非常大的矩阵：binomial(L,L/2) 其中 L 大约为 18 -22。现在我正在有大约 7.7 GB 内存可用的机器上测试它，但最终我将可以访问具有 64GB 内存的 PC。

我从 Lapack++ 开始。一开始我的项目假设只解决对称实矩阵的这个问题。

这个图书馆很棒。非常快速和小的 RAM 消耗。它可以选择计算特征向量并将其放入输入矩阵 A 以节省内存。有用！我认为 Lapack++ eigensolver 可以处理 Hermitian 矩阵，但由于未知原因它不能（也许我做错了什么）。我的项目已经发展，我应该也可以计算厄米特矩阵的这个问题。

所以我尝试将库更改为 Armadillo 库。它工作正常，但不如 Lapack++ 好，后者将 mat A 替换为所有 eigenvec，但当然支持 hermitian 矩阵。

L=14 的一些统计数据

Lapack++ RAM 126MB 时间 7.9s 特征值 + 特征向量
Armadillo RAM 216MB 时间 12s eigenval
Armadillo RAM 396MB 时间 15s 特征值+特征向量

我们来计算一下：double变量大约是8B。我的矩阵有大小 binomial(14,7) = 3432，所以在理想情况下它应该有 3432^2*8/1024^2 = 89 MB。

我的问题是：是否可以修改或强制 Armadillo 像 Lapack++ 那样做一个漂亮的把戏？ Armadillo 使用 LAPACK 和 BLAS 例程。或者也许有人可以使用另一个库推荐另一种方法来解决这个问题？

P.S.: 我的矩阵真的很稀疏。它有大约 2 * binomial(L,L/2) 个非零元素。我尝试用 CSC 格式的 SuperLU 来计算它，但它不是很有效，对于 L=14 -> RAM 185MB，但时间 135s。

Answer 1

Lapackpp 和 Armadillo 都依赖 Lapack 来计算复杂矩阵的特征值和特征向量。 Lapack 库提供了不同的方法来对复杂的厄尔米特矩阵执行这些操作。

Armadillo库的函数zgeev() does not care about the matrix being Hermitian. This function is called by the Lapackpp library for matrices of type LaGenMatComplex in the function LaEigSolve. The function eig_gen()调用了这个函数
Armadillo 库的函数zheev() is dedicated to complex Hermitian matrices. It first call ZHETRD to reduce Hermitian matrix to tridiagonal form. Depending on whether the eigenvectors are needed, it then uses a QR algorithm to compute the eigenvalues and eigenvectors (if needed). The function eig_sym() 如果选择方法std 调用此函数。
如果选择方法dc，Armadillo 库的函数zheevd() does the same thing as zheev() if the eigenvectors are not required. Otherwise, it makes use of a divide and conquert algorithm (see zstedc()). The function eig_sym() 将调用此函数。由于事实证明分而治之对于大型矩阵更快，因此它现在是默认方法。

Lapack 库中提供了具有更多选项的函数。（参见 zheevr() or zheevx）。如果想保持密集的矩阵格式，也可以试试Eigen库的ComplexEigenSolver

这是一个使用 Lapack 库的 C 包装器 LAPACKE 的小 C++ 测试。它由 g++ main.cpp -o main2 -L /home/...../lapack-3.5.0 -llapacke -llapack -lblas

编译

#include <iostream>

#include <complex>
#include <ctime>
#include <cstring>

#include "lapacke.h"

#undef complex
using namespace std;

int main()
{
    //int n = 3432;

    int n = 600;

    std::complex<double> *matrix=new std::complex<double>[n*n];
    memset(matrix, 0, n*n*sizeof(std::complex<double>));
    std::complex<double> *matrix2=new std::complex<double>[n*n];
    memset(matrix2, 0, n*n*sizeof(std::complex<double>));
    std::complex<double> *matrix3=new std::complex<double>[n*n];
    memset(matrix3, 0, n*n*sizeof(std::complex<double>));
    std::complex<double> *matrix4=new std::complex<double>[n*n];
    memset(matrix4, 0, n*n*sizeof(std::complex<double>));
    for(int i=0;i<n;i++){
        matrix[i*n+i]=42;
        matrix2[i*n+i]=42;
        matrix3[i*n+i]=42;
        matrix4[i*n+i]=42;
    }

    for(int i=0;i<n-1;i++){
        matrix[i*n+(i+1)]=20;
        matrix2[i*n+(i+1)]=20;
        matrix3[i*n+(i+1)]=20;
        matrix4[i*n+(i+1)]=20;

        matrix[(i+1)*n+i]=20;
        matrix2[(i+1)*n+i]=20;
        matrix3[(i+1)*n+i]=20;
        matrix4[(i+1)*n+i]=20;
    }

    double* w=new double[n];//eigenvalues

    //the lapack function zheev
    clock_t t;
    t = clock();
    LAPACKE_zheev(LAPACK_COL_MAJOR,'V','U', n,reinterpret_cast< __complex__ double*>(matrix), n, w);
    t = clock() - t;
    cout<<"zheev : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
    cout<<"largest eigenvalue="<<w[n-1]<<endl;

    std::complex<double> *wc=new std::complex<double>[n];
    std::complex<double> *vl=new std::complex<double>[n*n];
    std::complex<double> *vr=new std::complex<double>[n*n];

    t = clock();
    LAPACKE_zgeev(LAPACK_COL_MAJOR,'V','V', n,reinterpret_cast< __complex__ double*>(matrix2), n, reinterpret_cast< __complex__ double*>(wc),reinterpret_cast< __complex__ double*>(vl),n,reinterpret_cast< __complex__ double*>(vr),n);
    t = clock() - t;
    cout<<"zgeev : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
    cout<<"largest eigenvalue="<<wc[0]<<endl;

    t = clock();
    LAPACKE_zheevd(LAPACK_COL_MAJOR,'V','U', n,reinterpret_cast< __complex__ double*>(matrix3), n, w);
    t = clock() - t;
    cout<<"zheevd : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
    cout<<"largest eigenvalue="<<w[n-1]<<endl;

    t = clock();
    LAPACKE_zheevd(LAPACK_COL_MAJOR,'N','U', n,reinterpret_cast< __complex__ double*>(matrix4), n, w);
    t = clock() - t;
    cout<<"zheevd (no vector) : "<<((float)t)/CLOCKS_PER_SEC<<" seconds"<<endl;
    cout<<"largest eigenvalue="<<w[n-1]<<endl;

    delete[] w;
    delete[] wc;
    delete[] vl;
    delete[] vr;
    delete[] matrix;
    delete[] matrix2;
    return 0;
}

我电脑的输出是：

zheev : 2.79 seconds
largest eigenvalue=81.9995
zgeev : 10.74 seconds
largest eigenvalue=(77.8421,0)
zheevd : 0.44 seconds
largest eigenvalue=81.9995
zheevd (no vector) : 0.02 seconds
largest eigenvalue=81.9995

可以使用 Armadillo 库执行这些测试。直接调用 Lapack 库可能会让你获得一些内存，但 Lapack 的包装器在这方面也很有效。

真正的问题是您需要所有特征向量、所有特征值还是仅需要最大特征值。因为最后一种情况确实有高效的方法。看看Arnoldi/Lanczos itterative algorithms. Huge memory gains are possible if the matrix is sparce since only matrix-vector products are performed : there is no need to keep the dense format. This is what is done in the SlepC library, which makes use of the sparce matrix formats of Petsc. Here is an example of Slepc可以作为起点。

Answer 2

如果以后有人遇到和我一样的问题，我会在找到解决方案后提供一些提示（感谢所有发布一些答案或线索的人！）。

在英特尔网站上，您可以找到一些用 Fortran 和 C 编写的很好的示例。例如厄尔米特特征值问题例程 zheev()： https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/zheev_ex.c.htm

要使其在 C++ 中运行，您应该编辑该代码中的一些行：

在原型函数声明中对所有函数做类似的操作：extern void zheev( ... ) 更改为 extern "C" {void zheev( ... )}

更改调用 lapack 的函数添加 _ 符号，例如：zheev( ... ) 到 zheev_( ... )（只需通过替换即可在代码中为所有代码创建它，但我不明白为什么它会起作用。我通过一些实验弄明白了。）

您可以选择将 printf 函数转换为标准流 std::cout 并将包含的 headers stdio.h 替换为 iostream.

编译运行命令如：g++ test.cpp -o test -llapack -lgfortran -lm -Wno-write-strings

关于最后一个选项-Wno-write-strings，我现在不知道它在做什么，但是当他们在调用函数zheev( "Vectors", "Lower", ... )

低 RAM 消耗 c++ 特征求解器

low RAM consuming c++ eigen solver

c++

eigenvalue

eigenvector

lapack++

armadillo