What is matlab doing under the hood when I solve this generalized eigenvalue problem?

I have the following code for a non-hermitian generalized eigenvalue problem that runs extremely fast on matlab. This seems like it should be a hard problem, given its size and lack of nice properties such as symmetry and hermiticity. Can someone illucidate to me what's going on here to make this work so nicely?

dx = 0.0025;
xMin = -100;
xMax = 20;

a = 0.9; 
mu = 0.01;

n   = 1;

rp = 1 + sqrt(1 - a^2);
rm = 1 - sqrt(1 - a^2);

omegac = a / (2 * rp);
omegan = mu * (1 - mu^2 / (2 * (n + 2)^2));

x = (xMin : dx : xMax)';
xLeft  = x + dx;
xRight = x - dx;

r = rp * (exp(x) + 1);
rLeft  = rp * (exp(xLeft ) + 1);
rRight = rp * (exp(xRight) + 1);

N = length(x);

e = ones(N,1);

HElements        =  ((r      - rp)./(r      - rm));
HElementsLeft    =  ((rLeft  - rp)./(rLeft  - rm));
HElementsRight   =  ((rRight - rp)./(rRight - rm));
Vt0Elements      =  (a^2)./((r - rm).^2) - HElements .*(mu^2 * r.^2 + 2);
Vt1Elements      = -(4 * a * r)./((r - rm).^2);
Vt2Elements      =  ((r.^2 +  a^2).^2)./((r - rm).^2) - HElements * a^2;
B0               = 1 + dx * omegac * 2 * 1i * (rp/(rp - rm));
B1               =            - dx * 2 * 1i * (rp/(rp - rm));

D2Left   = (1/(dx^2) - HElementsLeft /(2 * dx));
D2Middle = (Vt0Elements - 2 /(dx^2));
D2Right  = (1/(dx^2) + HElementsRight/(2 * dx));

DLeft   = D2Left;
DMiddle = D2Middle;
DRight  = D2Right;
D2 = spdiags([DLeft DMiddle DRight] ,-1:1,N,N);
D2(1,1) = 0;
D2(1,2) = 0;
Vt1 = spdiags(Vt1Elements,0,N,N);
Vt1(1,1) = 0;
Vt2 = spdiags(Vt2Elements,0,N,N);
Vt2(1,1) = 0;
Id = speye(N,N);
Z = sparse(N,N);
S11 = Z;
S12 = Z;
S11(1,1) = 1;
S12(1,2) = 1;


A11 = S11 * B0 - S12 + D2;
A12 = S11 * B1 + Vt1 + omegan * Vt2;
A21 = Id * omegan;
A22 = Id * (-1);

A = [A11 A12;A21 A22];

B11 = Z;
B12 = Vt2 * (-1);
B21 = Id  * (-1);
B22 = Z;

B = [B11 B12;B21 B22];

[V,D] = eigs(A,B,1,'smallestabs');

Warning: First input matrix is close to singular or badly scaled (RCOND = 2.060539e-16) and results may be inaccurate. Consider specifying a small nonzero numeric sigma value instead of 'smallestabs' to improve the condition of the matrix.

disp(D)

  -5.1680e-12 + 5.8500e-21i

Answer 1

[Christine Tobler]

You can use the Display option to get some more information on what's going on inside of eigs.

load matrices.mat
[V,D] = eigs(A,B,1,"smallestabs", Display=true);

=== Generalized eigenvalue problem A*x = lambda*B*x ===

The eigenvalue problem is complex non-Hermitian.
Matrix B is not Hermitian positive (semi-)definite.

Computing 1 eigenvalues of type 'smallestabs'.


Parameters passed to Krylov-Schur method:
  Maximum number of iterations: 300
  Tolerance: 1e-14
  Subspace Dimension: 20

Find eigenvalues of A\(B*x) = mu*x.
Compute decomposition of A...

Warning: First input matrix is close to singular or badly scaled (RCOND = 2.060539e-16) and results may be inaccurate. Consider specifying a small nonzero numeric sigma value instead of 'smallestabs' to improve the condition of the matrix.

--- Start of Krylov-Schur method ---
Iteration   1: 1 of 1 eigenvalues converged.

D

D = -5.1680e-12 + 5.8500e-21i

norm(A*V - B*V*D)

ans = 1.5547e-06

So for the case of a non-hermitian generalized eigenvalue problem A*x = lambda*B*x, eigs iteratively solves the simple eigenvalue problems A\(B*x) = mu*x.

This works well unless A is close to singular. Unfortunately, that's the case for your matrix A, and it looks like the shifted matrix A - sigma*B is also singular for any value of sigma, so there's not much to do about this unless you have some way to project out the singularity of this system.

However, based on the residual the found eigenvalue and eigenvector are quite accurate.

For cases where eigs gives a warning about the first input matrix being singular, I would recommended

a) trying to choose a shift sigma for which A is not singular (not possible here)

b) making sure to double-check the output to see the effect of this ill-conditioning on the computed result.

[David Cyncynates]

This is really helpful, thank you so much!