犰狳:eigs_gen 最小特征值
Armadillo: eigs_gen for smallest eigenvalue
我正在使用犰狳的 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
,不完全是这里构造的矩阵,但我想问题应该是一样的。
问题在于重复的特征值;如果我将前两个矩阵元素更改为
A(1,2) = -1.00000001;
A(2,1) = -1.00000001;
达到预期效果
我认为这里的问题是您在对称矩阵上 运行 eigs_gen()
(调用 DNAUPD
)。 ARPACK 指出 DNAUPD
不适用于对称矩阵,但没有指定如果您仍然使用对称矩阵会发生什么:
NOTE: If the linear operator "OP" is real and symmetric with respect to the real positive semi-definite symmetric matrix B, i.e. B*OP = (OP')*B, then subroutine ssaupd should be used instead.
(来自 http://www.mathkeisan.com/usersguide/man/dnaupd.html )
我修改了 Armadillo 内部代码以将 "sa"(最小代数)传递给 eigs_sym() (sp_auxlib_meat.hpp) 中的 ARPACK 调用,并且我能够获得正确的特征值。我已经在上游提交了一个补丁,使 "sa" 和 "la" 支持可用于 eigs_sym(),我认为一旦新版本发布(或在某个时候发布),这应该可以解决您的问题未来)。
我正在使用犰狳的 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
,不完全是这里构造的矩阵,但我想问题应该是一样的。
问题在于重复的特征值;如果我将前两个矩阵元素更改为
A(1,2) = -1.00000001;
A(2,1) = -1.00000001;
达到预期效果
我认为这里的问题是您在对称矩阵上 运行 eigs_gen()
(调用 DNAUPD
)。 ARPACK 指出 DNAUPD
不适用于对称矩阵,但没有指定如果您仍然使用对称矩阵会发生什么:
NOTE: If the linear operator "OP" is real and symmetric with respect to the real positive semi-definite symmetric matrix B, i.e. B*OP = (OP')*B, then subroutine ssaupd should be used instead.
(来自 http://www.mathkeisan.com/usersguide/man/dnaupd.html )
我修改了 Armadillo 内部代码以将 "sa"(最小代数)传递给 eigs_sym() (sp_auxlib_meat.hpp) 中的 ARPACK 调用,并且我能够获得正确的特征值。我已经在上游提交了一个补丁,使 "sa" 和 "la" 支持可用于 eigs_sym(),我认为一旦新版本发布(或在某个时候发布),这应该可以解决您的问题未来)。