Accuracy of eig in support of complex step differentation: Derivative estimate accuracy degrades when step gets too small

17 views (last 30 days)
Mark Stone
Mark Stone on 15 May 2015
Edited: Mark Stone on 21 May 2015
Consider the uses of complex step differentiation to estimate the derivative of an eigenvalue of a real non-symmetric matrix, using eig in MATLAB double precision, for the calculation. The traditional recommendation is to use step = 1e-20 * 1i. But perhaps extended precision calculation of eig is needed to support this?
I noticed that for steps smaller than perhaps about 1e-14 *1i (or maybe even less), the accuracy of the complex step differentiation estimate seems to degrade, and become unstable. Is this due to accuracy limitation in eig in MATLAB double precision? The matrix in question has various norms in the ballpark of 1. Furthermore, the same calculation showed greater degradation in MATLAB R2013A WIN32 than in R2014A WIN 64. Is there any reason to think the accuracy of eig should be different between these configurations?
The calculation is carried out as d = step length (say 1e-14 or 1e-20 or whatever). For simplicity, I'll show this for the maximum eigenvalue. Consider a real matrix A. Then to compute the (i,j) component of the gradient of the maximum eigenvalue with respect to A,
B = A;
B(i,j) = A(i,j) + d*1i;
derivative_estimate = imag(max(eig(B))/d)
Complex step differentiation is not "supposed to" have a problem with small step size, and in fact, that's supposed to be its advantage over forward or central differences, whose accuracy is limited by a step size below which accuracy degrades.
Thanks.
  3 Comments
Mark Stone
Mark Stone on 19 May 2015
Matt,
I stated that the matrix A is real. Therefore, the machine precision limitation is not encountered on adding a small imaginary number. So imag(B(i,j)) = d*1i, and there is no loss of precision in this addition. That's the only reason why anyone would ever recommend a step magnitude of 1e-20 when calculations are being performed in double precision. Therefore, I believe the key matter here is the accuracy of eig, and what that means about smallest step size.

Sign in to comment.

Accepted Answer

Cleve Moler
Cleve Moler on 19 May 2015
I have to admit that I am unsure about what is going on here. Complex step differentiation relies upon the unperturbed function values for real arguments being real. Then the complex perturbation causes the software to produce a complex result, so numerical differencing is not a concern. But eig of a real nonsymmetric matrix can be complex and taking max just makes things worse. Consider A = gallery(5). The exact eigenvalues are all zero. But eig(A) doesn't reveal that. And a complex step isn't any help. The eigenvalue condition numbers are crucial. For multiple eigenvalues without a full set of eigenvectors, such as gallery(5), the condition numbers are infinite. Taking a complex step, or any step, in the vicinity of a badly conditioned eigenvalue will have numerical difficulties. -- Cleve

More Answers (3)

Cleve Moler
Cleve Moler on 19 May 2015
But you should not have to resort to complex step differentiation because an analytic expression for the derivative of an eigenvalue is known. If y and x are the left and right eigenvectors of a matrix A corresponding to a particular eigenvalue lambda, then d(lambda) = (y'*d(A)*x)/(y'*x)
  3 Comments
Mark Stone
Mark Stone on 21 May 2015
Well, after additional consideration, i guess that
{V,D] = eig(A);
U = inv(V)';
U will work as a perfectly matched left eigenvector matrix, so long as V is not defective (i.e., full rank, therefore non-singular), so this is not really a general purpose approach.
Prior to my accidental discovery that left eigenvectors are not supportd in eig in R2013A, I had been thinking that they were supported since time immemorial, for instance back to when I first used the Fortran-based MATLAB as a student in early 1982. I guess this passed through the cracks of my knowledge because i hadn't really started dealing with left eigenvectors until fairly recently.

Sign in to comment.


Cleve Moler
Cleve Moler on 20 May 2015
I had forgotten about Nick Higham's comment after my blog post referencing his paper, "The complex step approximation to the Fréchet derivative of a matrix function", http://link.springer.com/article/10.1007/s11075-009-9323-y. Unfortunately, it's behind a Springer pay wall. But I think it says that Mark is on the right track.
  1 Comment
Mark Stone
Mark Stone on 20 May 2015
The as submitted, but not as revised Tech Report version of that paper is at http://eprints.ma.man.ac.uk/1260/01/covered/MIMS_ep2009_31.pdf .
It contains this interesting quote in the concluding remarks: "The main weakness of the CS approximation is that it is prone to damaging cancellation when the underlying method for evaluating f employs complex arithmetic. But for many algorithms, such as those based on real polynomial and rational approximations or matrix iterations, this is not a concern." So I guess this was what you were referring to about eig, the use of complex arithmetic, which therefore makes it subject to non-trivial cancelation error at small step sizes, effectively limiting minimum step size. I now have a deeper understanding of complex step differentiation and its limitations. Actually, in the scalar case anyway, as neat as it seems, I'm not really sure what the big attraction is, since, in the best case, it seems to amount essentially to automatic differentiation, which does not suffer from these little quirks. Of course, back before automatic differentiation (such as in 1967), it would have been a different story.

Sign in to comment.


Matt J
Matt J on 15 May 2015
I don't think eigenvalues are in general, differentiable as a function of the matrix elements. However, even where they are, shouldn't your derivative estimate involve the difference eig(B)-eig(A) ?
  3 Comments
Mark Stone
Mark Stone on 20 May 2015
In the case I was looking at, in fact the eigenvalues in question were real and unique in magnitude. There was no tie or ambiguity in order, even though there were ties for other eigenvalues. Such an eigenvalue is supposed to be infinitely differentiable.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by