犰狳：最小特征值的 eigs_gen答案

【问题标题】：Armadillo: eigs_gen for smallest eigenvalue犰狳：最小特征值的 eigs_gen
【发布时间】：2015-01-20 18:42:43
【问题描述】：

我正在使用犰狳的 eigs_gen 来查找稀疏矩阵的最小代数特征值。

如果我只为最小的特征值请求函数，结果是不正确的，但如果我为 2 个最小的特征值请求它，结果是正确的。代码是：

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;


int
main(int argc, char** argv)
  {
  cout << "Armadillo version: " << arma_version::as_string() << endl;

  sp_mat A(5,5);

  A(1,2) = -1;
  A(2,1) = -1;
  A(3,4) = -1;
  A(4,3) = -1;

  cx_vec eigval;
  cx_mat eigvec;

  eigs_gen(eigval, eigvec, A, 1, "sr");  // find smallest eigenvalue ---> INCORRECT RESULTS
  eigval.print("Smallest real eigval:");

  eigs_gen(eigval, eigvec, A, 2, "sr");  // find 2 smallest eigenvalues ---> ALMOST CORRECT RESULTS 
  eigval.print("Two smallest real eigvals:");

  return 0;
  }

我的编译命令是：

 g++ file.cpp -o file.exe -O2 -I/path-to-armadillo/armadillo-4.600.3/include -DARMA_DONT_USE_WRAPPER -lblas -llapack -larpack

输出是：

Armadillo version: 4.600.3 (Off The Reservation)
Smallest real eigval:
    (+1.000e+00,+0.000e+00)
Two smallest real eigvals:
    (-1.000e+00,+0.000e+00)
    (-1.164e-17,+0.000e+00)

感谢任何关于为什么会发生这种情况以及如何克服这种情况的想法。

注意：第二个结果几乎是正确的，因为我们期望 -1、-1 作为两个最低特征值，但可能会忽略重复的特征值。

更新：包括一个测试矩阵构造，在 ryan 将“sa”选项包含到库中之后，该构造似乎没有收敛：

    #define ARMA_64BIT_WORD
    #include <armadillo>
    #include <iostream>
    #include <vector>
    #include <stdio.h>

    using namespace arma;
    using namespace std;

    int main(){

    size_t l(3), ls(l*l*l);
    sp_mat A = sprandn<sp_mat>(ls, ls, 0.01);
    sp_mat B = A.t()*A;

    vec eigval;
    mat eigvec;
    eigs_sym(eigval, eigvec, B, 1, "sa");

    return 0;

    }

感兴趣的矩阵大小要大得多，例如ls = 8000 - 27000，并不是这里构造的矩阵，但我认为问题应该是一样的。

【问题讨论】：

标签： sparse-matrix armadillo eigenvalue

【解决方案1】：

我认为这里的问题是您在对称矩阵上运行 eigs_gen()（调用 DNAUPD）。 ARPACK 指出 DNAUPD 不适用于对称矩阵，但没有说明如果您使用对称矩阵会发生什么：

注意：如果线性算子“OP”是实数并且相对于实半正定对称矩阵 B 对称，即 B*OP = (OP')*B，则应使用子程序 ssaupd。

（来自http://www.mathkeisan.com/usersguide/man/dnaupd.html）

我修改了内部 Armadillo 代码以将“sa”（最小代数）传递给 eigs_sym() (sp_auxlib_meat.hpp) 中的 ARPACK 调用，并且我能够获得正确的特征值。我已经向上游提交了一个补丁，以使 eigs_sym() 可以使用“sa”和“la”支持，我认为一旦发布新版本（或在未来的某个时间），这应该可以解决您的问题。

【讨论】：

谢谢，这实际上是对图书馆的一个非常期望和实质性的改进。
实际上“sa”选项似乎存在收敛问题，而“lm”选项则没有。我已经包含了一个代码更新，显示了一个测试矩阵结构，这对我来说似乎没有收敛。
我认为这里的问题不是犰狳对 ARPACK 的包装，而是 ARPACK 使用的数值方法。更好地调查输入矩阵的属性以及隐式重新启动的 Arnoldi 方法是否适合您的特定矩阵可能是接下来要做的事情。

【解决方案2】：

问题在于重复的特征值；如果我将前两个矩阵元素更改为

A(1,2) = -1.00000001;
A(2,1) = -1.00000001;

得到了预期的结果。

【讨论】：