犰狳 eigs_sym：分解失败

Posted 2023-02-17

【中文标题】犰狳 eigs_sym：分解失败

【英文标题】：Armadillo eigs_sym: decomposition failed

【发布时间】：2017-06-08 23:25:48

【问题描述】：

我正在 Visual Studio 中编写 C++ 代码并使用 Armadillo 7.900.1

我无法让 eigs_sym 函数在 Armadillo 下工作（我使用的是 Armadillo 附带的 lapack 和 blas 版本）。典型的错误信息如下：

error: lapack::stedc(): failed to compute all eigenvalues
warning: eigs_sym(): decomposition failed

产生这个输出的代码是：

sp_mat T(locations, values);
arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

在哪里

T is a 480,000x480,000 sparse matrix. 

如果 T 很小（即 2000x2000）并且 tol 很高（tol = 5），我可以让代码工作。但是一旦 T 约为 20000x20000，那么无论使用什么 tol 或 num_eigs_wanted 值都会出现上述错误。 （我显然改变了“值”和位置”来改变 T 的大小）。

矩阵 T 是对称的、实数的、正定的。

当我在Matlab中调用eigs函数时，完全相同的矩阵没有问题。

有什么想法吗？在我看来，犰狳中的 eigs_sym 似乎被破坏了......人们使用过任何替代品吗？

干杯

【问题讨论】：

你用的是哪个矩阵？

您确定要最小幅度"sm" 特征值吗？它们通常是垃圾并且难以可靠地计算。将"sm" 更改为"lm" 以获得最大幅度（如documentation 中所述）。如果你真的想要最小的，增加想要的特征值的数量：***.com/questions/36059755/…

不——我要求最小的是有原因的。无论如何它不应该工作吗？如果您必须知道，它也不适用于“lm”......相同的错误响应

【参考方案1】：

不是答案，但评论太长了。我可以用下面的代码重现。该矩阵只是一个对角矩阵，对角线上的值从 1 开始递增。特征值显然是 1,2,3,...，并且矩阵是正定且对称的。

#include <iostream>
#include <armadillo>

int main()

  constexpr int N = 20000;
  arma::umat locations(2,N);
  arma::vec values(N);

  for (int i = 0; i < N; ++i)
  
    locations(0,i) = i;
    locations(1,i) = i;
    values(i) = i+1;
  

  arma::sp_mat T(locations, values);

  int num_eigs_wanted = 1;
  arma::vec eigval(num_eigs_wanted);
  arma::mat eigvec(N,num_eigs_wanted);
  double tol = 1e-6;
  arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

  std::cout << eigval << '\n';

编译器标志是

clang++-5.0 -std=c++11 test.cpp -larmadillo

我正在使用带有 MKL 后端的犰狳，我收到了：

warning: eigs_sym(): decomposition failed
[matrix size: 0x1]

---

这似乎是犰狳的问题，因为当我使用 Eigen 的不受支持的 Arpack 支持时，它工作得很好。

#include <iostream>
#include <Eigen/Sparse>
#include <unsupported/Eigen/ArpackSupport>

int main()

  constexpr int N = 20000;

  std::vector< Eigen::Triplet<double> > triplets;
  triplets.reserve(N);
  for (int i = 0; i < N; ++i)
    triplets.push_back( i,i,i+1 );

  Eigen::SparseMatrix<double> T(N,N);
  T.setFromTriplets(triplets.begin(), triplets.end());

  int num_eigs_wanted = 1;
  double tol = 1e-6;
  Eigen::ArpackGeneralizedSelfAdjointEigenSolver< Eigen::SparseMatrix<double> > eigs_sym;
  eigs_sym.compute(T, num_eigs_wanted, "SM", Eigen::ComputeEigenvectors, tol);
  std::cout << eigs_sym.eigenvalues() << '\n';

您可能还想看看https://github.com/yixuan/spectra/

【讨论】：

干杯！ - 快速玩过 Spectra。看起来这可能是最好的解决方案。谢谢。将更新我的进度。

【参考方案2】：

感谢@Henri 的所有帮助！ Spectra (https://spectralib.org) 是必经之路！

在大小为 480,000 x 480,000 的稀疏矩阵上使用 SymEigsSolver 进行特征值/向量分解需要 3 分 50 秒（i7/16GB RAM）。

大小为 120,000x120,000 的稀疏矩阵需要 20 秒。

没有稳定性问题，程序员可以直接访问 Arnoldi 算法、Ritz 值和公差阈值。

由于 Spectra 完全基于标头，因此也易于使用。只需#include 头文件 - 很简单。您还需要 Eigen，它也可以在上述地址获得（并且它也是基于标头的）。