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

Question

我是C++编程的新手，但我有一个任务来计算对称矩阵(和埃尔米特矩阵)的特征值和特征向量(标准特征问题Ax=lx))对于尺寸非常大的矩阵:二项式(L,L/2)，其中L约为18-22。现在我正在具有大约 7.7 GB RAM 可用的机器上进行测试，但最终我将可以访问具有 64GB RAM 的 PC。

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

这个图书馆很棒。非常快且占用内存小。它可以选择计算特征向量并将其放入输入矩阵 A 以节省内存。有用!我认为 Lapack++ 特征求解器可以处理 Hermitian 矩阵，但由于未知原因它不能处理(也许我做错了什么)。我的项目已经发展，我应该也能够计算埃尔米特矩阵的这个问题。

所以我尝试将库更改为 Armadillo 库。它工作得很好，但它不如 Lapack++ 那么好，它用所有 eigenvec 替换 mat A，但当然支持埃尔米特矩阵。

L=14 的一些统计

• Lapack++ RAM 126MB 时间 7.9s 特征值 + 特征向量

• Armadillo RAM 216MB 时间 12 秒特征值

• Armadillo RAM 396MB 时间 15s 特征值+特征向量

我们来计算一下:double 变量约为 8B。我的矩阵有大小二项式(14,7) = 3432，因此在理想情况下，它应该具有3432^2*8/1024^2 = 89 MB。

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

附注:我的矩阵非常稀疏。它有大约 2 * 二项式(L,L/2) 个非零元素。我尝试使用 CSC 格式的 SuperLU 进行计算，但效果不是很好，L=14 -> RAM 185MB，但时间为 135 秒。

Answer 1

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

• 函数zgeev()不关心矩阵是 Hermitian 矩阵。此函数由 Lapackpp 库在函数 LaEigSolve 中针对 LaGenMatComplex 类型的矩阵调用。 。函数eig_gen() Armadillo 库的 调用此函数。

• 函数zheev()致力于复杂的埃尔米特矩阵。它首先调用 ZHETRD 将 Hermitian 矩阵简化为三对角形式。根据是否需要特征向量，它然后使用 QR algorithm计算特征值和特征向量(如果需要)。函数eig_sym()如果选择了方法 std，则 Armadillo 库的 调用此函数。

• 函数zheevd()如果不需要特征向量，则与 zheev() 执行相同的操作。否则，它会使用分治算法(请参阅 zstedc() )。函数eig_sym()如果选择了方法dc，则 Armadillo 库的 调用此函数。由于事实证明分而治之对于大型矩阵更快，因此它现在是默认方法。

Lapack 库中提供了具有更多选项的函数。 (参见 zheevr() 或 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 的包装器在这方面也可以高效。

真正的问题是您是否需要所有特征向量、所有特征值或仅需要最大特征值。因为最后一种情况确实有有效的方法。看看阿诺尔迪/Lanczos迭代算法。如果矩阵是稀疏的，则可能会获得巨大的内存增益，因为仅执行矩阵 vector 乘积:无需保持密集格式。这就是 SlepC 库中所做的事情，它利用了 Petsc 的稀疏矩阵格式。 Here is an example of Slepc可以作为起点。