Preprocessing when solving generalized eigenvalue problems using the eigs function

I am solving a generalized eigenvalue problem using the eigs function. Both matrices A and B are symmetric and B is real positive. They are also sparse matrices.In this case, there is a difference in convergence between

[V, D] = eigs(A,B,50, "smallestabs")

and

[V, D] = eigs(A,B,50, "smallestreal")

According to the documentation, this depends on whether or not an inverse matrix is used, and is displayed as follows when opts.disp=2.

smallestabs:

Find eigenvalues of R*A\(R'*y) = mu*y, with y = R*x, B = R'*R.

No need to compute R, as it is only used implicitly.

Compute decomposition of A...

smallestreal:

...

R'\(A*(R\x)) = lambda*x, with B = R'*R.

...

It seems that the standard eigenvalue problem for the new matrix Anew is transformed and then its inverse is taken.

How is Anew calculated in this case?

Is it possible without calculating R?

I would like to use Anew in other solvers of my own making.

[Christine Tobler]

The matrix R here is computed using Cholesky decomposition, R = chol(B);

(eigs also detects the case of a singular, positive semi-definite B, for the 'smallestabs' case, in a different way - but this is only useful for some rare cases).

Keep in mind the matrices R*inv(A)*R' or inv(R)*A*inv(R)' are likely to be not very sparse at all, so you wouldn't want to store them explicitly. Usually people will uses these implicitly, by passing a function handle that computes R*A\(R'*x) or R'\(A*(R\x)) and treat this like a function handle that applies Anew*x or Anew\x.

If you solve a linear system with the same matrix multiple times, the decomposition object can make this cheaper by computing the factorization once and then reusing it for each solve operation.

To avoid computing R at all, usually you have to reformulate your algorithm. Let me show you with the example of the inverse iteration:

• Start off using Anew = inv(R)*A*inv(R)', with B = R'*R
• The algorithm looks like this:

for i=1:maxit
    x = Anew \ x;
    x = x / norm(x);
end

Now we insert the definition of Anew:

for i=1:maxit
    x = R*(A\(R'*x));
    x = x / norm(x);
end

The trick now is to substitute y = inv(R)*x (meaning x = R*y) in the algorithm above (using inv here since it's easier for rearranging equations - will use backslash for computation since inv would be too slow):

xnew = R*inv(A)*R'*xold
   => inv(R)*xnew = inv(A)*R'*(R*inv(R)*xold)
   => ynew = inv(A) * (R'*R) * xold
   => ynew = inv(A)*B*xold

xnew = xold / norm(xold) 
   => xnew = xold / sqrt(xold'*xold)
   => R*ynew = R*yold / sqrt(xold'*R'*R*xold)
   => ynew = yold / sqrt(xold'*B*xold)   % inner product in B

This means we can write the whole for-loop in terms of y instead of x, and avoid needing an explicit R anywhere:

for i=1:maxit
    y = A\(B*y);
    y = y / sqrt(y'*B*y);
end

You may need R in the end to get back to x from y - but in many cases (for example in eigs) y is actually what you want, since Anew was only constructed as a helper matrix originally.