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

accuracycomplex step differentiationeig

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.

Best Answer

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

Related Solutions

MATLAB: Eigenvectors of a nonsymmetric matrix.

Below W are the left eigenvectors, use the one associate with the largest eigenvalue.

If A is your matrix

   [W,D] = eig(A.'); 
   W = conj(W)

MATLAB: Is the rank function reliable

Many things to discuss here.

1. Yes, rank is quite reliable. Is it perfect? Well, no numerical method can be perfect. One can trivially generate a matrix that is not singular, but in double precision determining the rank using any simple scheme will fail to determine the singular nature properly.

A = diag([1 1e-20])
A =
            1            0
            0        1e-20
rank(A)
ans =
     1

So A is clearly diagonal, with non-zero elements, but in double precision, it is numerically singular as far as rank can tell.

However, rank IS essentially as good a scheme as you can get in double precision. If rank tells you the matrix is singular in double precision, then effectively it is so. Anyway, floating point arithmetic limits what you can know about any such matrix.

2. Can you use cond instead of rank? Of course you can, but if rank tells you the matrix has full rank, then cond would do so too. The extra thing that cond tells you is how CLOSE the matrix is to singular, but you still need to make a judgment about the singularity. Similarly, looking at the singular values of the matrix would tell you something, IF you understand what you are looking at.

Regardless, just use rank. For the other tools that require judgement about how close a number is to big, if you don't know what that number should be, then you are FAR better off to just use rank! (My brother-in-law once asked me what a certain tool did, and if he needed it in his shop. My response was essentially that if he had no idea what the tool was, he did not need it.)

3. Evaluating det(A'*A), and wondering how close it need to be to zero is INSANE. The determinant is NEVER a good test of singularity. This is a myth perpetrated on students by authors of textbooks and papers, authors who sadly don't know any better. They keep reading papers by other authors of the same ilk, and repeating what they see. Of course, this is an infinite loop.

The determinant is NEVER a good test of singularity. PERIOD. PERIOD squared.

Once can trivially create a matrix that has any determinant you wish, simply by scaling that matrix by any arbitrary constant. Clearly scaling a matrix by a constant has NO effect on the rank of that matrix, however it DOES strongly affect the determinant. In fact, for an NxN matrix A

det(A) == 1

then

det(k*A) == k^N

So if k is small, then det(k*A) will be tiny! Likewise, if k is large, then det(k*A) will be immense. Since the value of k here has no impact on the singularity of k*A, clearly the det is meaningless to know if the matrix is singular. PERIOD. I don't care how many papers you have read by authors who have no understanding of the issues, how many books. There are a lot of fools (sorry, but true) out there who will endlessly repeat what they have read, not understanding the mathematics.

4. Of course, computing det(A'*A) is worse yet. Here we are squaring the condition of the matrix, so if we could not determine if A was singular using det, why in the name of god and little green apples would forming A'*A help?

5. Similarly, comparing rank(A) to rank(A'*A) is meaningless.

A = diag([1 1e-10]);
rank(A)
ans =
     2
rank(A'*A)
ans =
     1

Was A singular? Of course not. Rank told you so, but you can screw things up by forming A'*A when there was no need to do so. In fact, it is almost NEVER a good idea to form A'*A. Almost any problem that "needs" A'*A will almost always be better served if you read some of your old linear algebra 101 notes, in terms of a factorization of A.

DON'T compute A'*A. Instead, you might compute the QR factorization of A, then use that to form what is essentially a Cholesky factorization of A'*A, BUT keep it in factored form! I.e., if

A = Q*R

then

A'*A = R'*Q'*Q*R = R'*R

Again, DON'T compute R'*R. That would just get you in the same hot water as before. Keep R'*R in factored form. Anything you need to do with that can be gotten from R, and more efficiently AND more accurately. (READ YOUR NOTES! If you never took linear algebra, it is high time to do so, or at least read a book on the subject if you are asking questions like this.

Sorry about the rant. It is not directed at you personally. But if nobody ever tells you the issues and the reasons to avoid things like det, etc., you will never learn, and you will likely perpetuate the myth. There are GOOD numerical schemes to be found for many of the things we do. At the same time, there are many schemes out there that are numerically useless or worse, because they are dangerous if people use them without understanding the issues.

set('soapbox','off','rant','off')

Best Answer

Related Solutions

MATLAB: Eigenvectors of a nonsymmetric matrix.

MATLAB: Is the rank function reliable

Related Question