MATLAB: Am I unable to use non-Hermitian matrices with EIGS for a Generalized Complex Eigenvalue problem

MATLAB

When I try to use a non-Hermitian matrix in the EIGS command for Generalized Complex Eigenvalue Problems, I get an error saying that B should be a symmetric positive (semi-) definite matrix.
But EIGS uses ARPACK and ARPACK can take in non-Hermitian matrices and yet I am getting this error.

Best Answer

ARPACK does NOT accept a non-Hermitian matrix B when attempting to solve the generalized eigenvalue problem. This is due to the choice of algorithm, namely, Implicitly restarted Arnoldi with respect to the B-inner product. When B is non-Hermitian (or Hermitian indefinite), it does not define an inner product. So EIGS will need a Hermitian matrix as the second input for a Generalized Complex Eigenvalue Problem.
If you want to use a non-Hermitian matrix B, one of the following workarounds can be used:
1. If A is Hermitian positive definite, the user could solve Bx = \mu
Ax, and then transform \mu into \lambda.
2. If you want to compute the eigenvalues of largest magnitude,
he must apply B^{-1} every time he needs to apply A. That is, call
eigs(@(x) B\(A*x), size(A,1), ...);
Or instead,
[L,U,P,Q] = lu(B);
myMatVec = @(x) Q*(U\(L\(P*(A*x))));  %  Compute B\(A*x)
eigs(myMatVec, size(A,1), ...);
3. If you want to compute the eigenvalues nearest to some target
c, you should work with the operator C := (A - c*B)^{-1}*B. That is, you
could perform the following:
% Suppose user wants 6 eigenvalues nearest to c = 2.
[L,U,P,Q] = lu(A - c*B);  %  Provided A - c*B is sparse.  Otherwise [L,
U, P] = lu(A - c*B);
myMatVec = @(x) Q*(U\(L\(P*(B*x))));
sigma = eigs(myMatVec, size(A,1), 6, 'lm');
lambda = (1./sigma) + c;  % These are the actual eigenvalues.
To see how this works, you can try it with small enough (50 x 50, say)
dense examples:
A = rand(50);
B = rand(50);
c = 2;
nev = 6;  %  Number of eigenvalues.
actual = eig(A,B);
[v, ind] = sort(abs(actual - c));
actLam = actual(ind(1:nev));
[L,U,P] = lu(A - c*B);
myMatVec = @(x) U\(L\(P*(B*x)));
sigma = eigs(myMatVec, 50, nev, 'lm');  %  Compute transformed
eigenvalues.
lambda = 1./sigma + c;  %  Transform to actual eigenvalues.
%  Compare lambda to actLam
Related Question