【问题标题】:low RAM consuming c++ eigen solver低 RAM 消耗 C++ 特征求解器
【发布时间】:2015-11-22 23:49:30
【问题描述】:

我是 C++ 编程的新手,但我的任务是计算对称矩阵(和厄米特矩阵)的特征值和特征向量(标准特征问题 Ax=lx))对于尺寸非常大的矩阵:binomial(L,L/2),其中 L 约为 18-22。现在我在有大约 7.7 GB 内存的机器上测试它,但最终我可以使用 64GB 内存的电脑。

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

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

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

L=14 的一些统计数据

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

  • 犰狳 RAM 216MB 时间 12s 特征值

  • 犰狳 RAM 396MB 时间 15s eigenval+eigenvectors

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

我的问题是:是否可以修改或强制 ArmadilloLapack++ 那样做一个很好的把戏? 犰狳 使用LAPACKBLAS 例程。或者也许有人可以推荐使用另一个库来解决这个问题的另一种方法?

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

【问题讨论】:

  • 您是否考虑过将 SLEPcEigenSolver 与 PETSc 一起使用?
  • @symphonic 我听说过,但我认为这个库对我来说太难了。如果我有时间,我会尝试测试它。
  • 在 LAPACK 中,您有例程 zheev (apfel.mathematik.uni-ulm.de/~lehn/FLENS/cxxlapack/netlib/lapack/…) 用于计算厄米矩阵的特征值和 - 向量。但是,您在 Lapack++ 中似乎是对的,它们似乎不支持复值矩阵(或至少不支持厄米矩阵)。所以也许直接从 C++ 调用 LAPACK 例程是个好主意。
  • @MichaelLehn 感谢您的回复。我知道LAPACK有这样的例程,但是我还没有学习如何在C++中直接调用它,所以我会先尝试学习如何调用它! :)
  • 在 FLENS 中有用于 LAPACK 的 C++ 包装器。包括对于zheev。所以基本上你可以撕掉FLENS。您也可以先检查 FLENS 是否在为您完成这项工作。这里 (apfel.mathematik.uni-ulm.de/~lehn/FLENS/flens/examples/…) 是计算对称矩阵的特征值/向量的示例。对于 Hermitian 矩阵,github 中只有一个示例 (github.com/michael-lehn/FLENS/blob/public/flens/examples/…),但在线教程中还没有。本教程中的这些示例仅设置小型演示案例。

标签: c++ armadillo eigenvector eigenvalue lapack++


【解决方案1】:

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

  • zgeev() 函数不关心 Hermitian 矩阵。 Lapackpp 库在函数LaEigSolve 中为LaGenMatComplex 类型的矩阵调用此函数。 Armadillo 库的函数eig_gen() 调用了这个函数。

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

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

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 的包装器在这方面也很有效。

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

【讨论】:

    【解决方案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 并将包含的标头stdio.h 替换为iostream

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

    关于最后一个选项-Wno-write-strings,我现在不知道它在做什么,但是当英特尔的示例将字符串而不是字符放入调用函数zheev( "Vectors", "Lower", ... )时,可能存在问题

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2012-05-08
      • 2018-11-29
      • 1970-01-01
      • 1970-01-01
      • 2019-12-30
      相关资源
      最近更新 更多