I want to determine the spectral radius of a large non-symmetric matrix $A$ whose dominant eigenvalues are a pair of complex conjugates. My first instinct was to use a power iteration with a starting vector $x$ in the complex plane. I do recall that this method tends to have trouble when the dominant eigenvalue is complex and indeed my numerical experiments seem agree with this assessment. Is there a way I can modify the power iteration to converge to one of the two dominant eigenvalues in the complex plane? Is it more advisable to use a arnoldi algorithm to determine the spectral radius?

  • $\begingroup$ $A$ itself is real but nonsymmetric? $\endgroup$
    – hardmath
    Jun 11 '13 at 10:53
  • $\begingroup$ @hardmath: Yes. $\endgroup$
    – Paul
    Jun 11 '13 at 13:20

The power method indeed does not converge in the presence of multiple dominant eigenvalues of the same magnitude. (If you follow the proof, you can see that the iterates will get closer and closer to the subspace spanned by the dominant eigenvectors, but not to any particular vector in this eigenspace).

A restarted Arnoldi method is in fact the way to go, but it can be quite a bit simplified since you are only interested in the dominant eigenvalue. Recall that the Arnoldi method consists in computing a unitary matrix $Q\in\mathrm{C}^{n\times k}$ (consisting of basis vectors of the Krylov space spanned by $v,Av,A^2v,\dots A^{k-1}v$ for some given vector $v$) and then solving the smaller eigenvalue problem for $Q^*AQ$ (using the fact that it is in upper Hessenberg form).

The idea is now to restart every second step, i.e., given an iterate $x^k$, compute an orthonormal basis of the span of $\{x^k, Ax^k\}$ (since $x^k$ is assumed to be normalized, you only need to orthogonalize $Ax^k$ against $x^k$). Since projection methods best approximate the extremal eigenvalues, the eigenvalues of the projected matrix $Q^*AQ$ will converge to the complex conjugate pair of dominant eigenvalues. In fact, since you are only interested in one of the two eigenvalues, you can just take the Rayleigh quotient for $x^k$; the full iteration (started from a normalized complex vector x) is

w = A*x;
l = w'*x;
w = w-l*x;
x = w/norm(w);

where l converges to one of the dominant eigenvalues.

(This is Problem P-4.2 in Saad's book Numerical Methods for Large Eigenvalue Problems).

Edit: Paul explicitly asked about power methods, but for the sake of later readers I should caution that the above works only if $A$ is real and in fact has a pair of complex conjugate dominant eigenvalues. If this is not known a priori, it's better to do several power iterations (with a real starting vector), orthonormalize the last two iterates, compute the eigenvalues of the projected $2\times 2$ matrix, and repeat if necessary:

for k = 1:kmax
    x = w/norm(w);
    w = A*x;
w = w-(x'*w)*x;
w = w/norm(w);
Q = [x,w];
l = max(eig(Q'*A*Q));

(Depending on the available routines for dense linear algebra, that might even be faster for the original question as well due to the better convergence of Krylov methods. In fact, if you have access to ARPACK or something equivalent, eigs(A,1) (or its equivalent) will be pretty hard to beat, especially if you need high accuracy.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.