Why do I see a singular matrix warning here?

I've posted this question as a "response" to a question posed by someone else, but as an answer to another question. Since I cannot move that answer into a question by that person, I'll post the question myself, since this seems to be a question I have seen too often.

A = [2 1 -5; 0 1 -2; 0 0 2];
[V, D] = eig(A);

% Matrix of eigenvectors returned
V

V = 3×3

1.0000 -0.7071 1.0000 0 0.7071 -0.0000 0 0 0.0000

% Eigenvalues for those eigenvectors
diag(D)

ans = 3×1

2 1 2

% But now, when I try to recover A from the eigenvalues and eigenvectors,
% this fails, and I get a singular matrix warning? Why is that?
V\D*V

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.344132e-17.

ans = 3×3

2.0000 -0.7071 0.0000 0 1.0000 0.0000 0 0 2.0000

Even worse, the result does not reproduce A? What is wrong with eig?

A

A = 3×3

2 1 -5 0 1 -2 0 0 2

[Clara Casado-Coterillo]

when a jacobian in nlparci gives back this warning in nlparci, what should we do?

Answer 1

[John D'Errico]

The issue is not a problem with eig, except in our expectations of what eig can do. And those expecations are fueled by our experience that eig always seems to produce a set of linearly independent, as well as orthogonal, eigenvectors for our matrices. Great, when it works. But when does eig fail to work as we expect?

Again, we have expectations about eig, because most of the matrices we see as statisticians, as engineers, etc., are symmetric, and are even positive definite. Even the example problems we are given as students seem to be of that sort. But what happpens when A is NOT symmetric, or even triangular? What happens when the matrix has eigenvalues of multiplicity higher than 1?

https://en.wikipedia.org/wiki/Defective_matrix

The classic example has a matrix like the matrix A in my question:

A = [2 1 -5; 0 1 -2; 0 0 2]

A = 3×3

2 1 -5 0 1 -2 0 0 2

So is upper triangular, so clearly not the normal symmetric matrix we come to expect. Worse, it has an eigenvalue of multiplicity 2.

[V,D] = eig(A)

V = 3×3

1.0000 -0.7071 1.0000 0 0.7071 -0.0000 0 0 0.0000

D = 3×3

2 0 0 0 1 0 0 0 2

We see that 2 appears twice in the set of eigenvalues. This is our first clue that A falls in the class of defective matrices. We don't yet know that A is defective. Is it?

The second eigenvector listed is indeed an eigenvector, corresponding to eigenvalue 1.

[A*V(:,2) , V(:,2)*1]

ans = 3×2

-0.7071 -0.7071 0.7071 0.7071 0 0

As you can see, the second column of V is indeed an eigenvector. How about the first column of V?

[A*V(:,1) , V(:,1)*2]

ans = 3×2

2 2 0 0 0 0

And that clearly works too. Th problem arises when we look at that third column of V. If fact, down to floating point trash, it is the same thing as V(:,1).

norm(V(:,1) - V(:,3))

ans = 1.4186e-16

Here we see evidence of a defective matrix. While we see an eigenvalue of multiplicity 2, it lacks two independent eigenvectors for that eigenvalue.

Why then did the expression V \ D * V fail, generating a singular matrix warning? It failed, because V is a singular matrix. The columns of V do not span the R^3 space we normallty expect them to span. It is still true that this next result is zero, (or effectively zero, but for floating point trash):

norm(A*V - V*D)

ans = 4.4409e-16

But the point is, we do NOT have the relationship that V'*V is the identity matrix. In fact, the matrix V has rank 2. So we do not have the ability to write V\D*V and get something that makes sense.

rank(V)

ans = 2

And that is we we see a singular matrix warning when we try it.

V\D*V

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.344132e-17.

ans = 3×3

2.0000 -0.7071 0.0000 0 1.0000 0.0000 0 0 2.0000

Defective matrices are not the normal, run of the mill matrix we expect to see. Certainly not so based on things like covariance matrices from statistics, or stiffness or stress & strain matrices from engineeering. What is broken here is not eig, since there does not exist a complete set of eigenvectors for that matrix. (I suppose it would be a nice feature if eig would throw out a warning message when that happens, but it does not. Such is life.)

In that Wikipedia link I gave above, we see the statement

"Any nontrivial Jordan block of size 2 × 2 or larger (that is, not completely diagonal) is defective. (A diagonal matrix is a special case of the Jordan normal form and is not defective.)"

So what does Jordan tell us here?

jordan(A)

ans = 3×3

1 0 0 0 2 1 0 0 2

The Jordan form for A indeed has a 2x2 non-diagonal block in the lower right corner.

Finally, is the problem that A is upper triangular? We can change A simply enough, like this:

B = A;
B(3,3) = 3

B = 3×3

2 1 -5 0 1 -2 0 0 3

What would Jordan predict here?

jordan(B)

ans = 3×3

3 0 0 0 1 0 0 0 2

Ah, so now the Jordan form is diagonal. It would predict that everything will work now. Well, I hope it will.

[Vb,Db] = eig(B)

Vb = 3×3

1.0000 -0.7071 -0.9733 0 0.7071 -0.1622 0 0 0.1622

Db = 3×3

2 0 0 0 1 0 0 0 3

Again, V clearly has rank 3.

rank(Vb)

ans = 3

But does the form Vb\Db*Vb recover B?

Vb \ Db * Vb

ans = 3×3

2.0000 -0.7071 1.2978 0 1.0000 0.4588 0 0 3.0000

Ugh. This still fails to work. What happened now? Is there a problem with eig? Again, the problem lies in our expectations. Well, maybe in YOUR expectations. My expectations are low here. In the help for eig, we see the statement:

[V,D] = eig(A) produces a diagonal matrix D of eigenvalues and

a full matrix V whose columns are the corresponding eigenvectors

so that A*V = V*D.

So at best, we can expect this result will be effectively zero:

norm(B*Vb - Vb*Db)

ans = 0

And that is all we can hope to see.

[Paul]

"We see that 2 appears twice in the set of eigenvalues. And so A is a member of that class of defective matrices."

This statement as written suggest that A must be defective because it has a an eigenvalue of multiplicity two. But as stated later, and in the wikipedia link, defectiveness is determined by whether or not the eigenvectors are linearly independent.

So perhaps it would be better to say "And so A might be a member of that class of defective matrices." And then show that it is defective based on the eigenvectors.

[John D'Errico]

I knew there would be something I said that was not quite correct.

Good point by @Paul. Yes, just because a matrix has eigenvalues of higher multiplicity does not make the matrix defective. It may be, but that need not always be the case. For example, we could consider the diagonal matrix:

A = diag([2 2 1])

A = 3×3

2 0 0 0 2 0 0 0 1

A will not be defective, and eig will be exquisitely happy.

[Bruno Luong]

Paul "defectiveness is determined by whether or not the eigenvectors are linearly independent."

This statement is not correct interpretation, the Wikipedia states the equivalent condition is there is no n linear independent eigenvectors, where n is the size of matrix.

The simplest matrix one can think of is

J = [0 1; 
     0 0]

[Paul]

Agreed. When I wrote that I should have said " ... whether or not there exists a full set of linearly independent eigenvectors."

[Bruno Luong]

But all eigen vector of J is {[1; 0]} (at constant multiplication), and it is independent. (EDIT typo)

MATLAB EIG returns the same eigen vector twice, one cannot consider columns of P as two eigen vectors (in a math sense).

[Paul]

Then we agree that J is defective. I also agree with your point about eig(), but I didn't say anything about what eig() returns, so not sure if that statement was in response to any of my comments or if that was just statement of fact.

[Christine Tobler]

Great post, John! I just wanted to point out two MATLAB functions that can be used to check if a matrix is defective.

If you're using the symbolic toolbox, calling eig on a symbolic matrix will tell you this exactly: V only contains linearly independent eigenvectors, while D gives the multiplicities of each eigenvalue. The third output P can be used to map from the eigenvalues in D to the eigenvectors in P.

A = [2 1 -5; 0 1 -2; 0 0 2];
[V, D, P] = eig(sym(A));
V

V = 

diag(D)

ans = 

P

P = 1×2

1 2

If we're doing numeric computations, "is this matrix defective" can't be answered easily, because any tiny change to the matrix would change the answer. This is similar to the question of "is this matrix low-rank" - adding some tiny random noise will make the matrix be full-rank in exact arithmetic, but when doing numerical computations, we would want to know if the matrix is close up to round-off to a low-rank matrix.

So here also, we can answer the question "is this matrix close to being defective". In a way, the warning about U being ill-conditioned was already doing that, telling us that the eigenvectors found are close to being linearly dependent.

Some more detailed information is given by condeig:

[U, D, s] = condeig(A);
diag(D)

ans = 3×1

2 1 2

log10(s)

ans = 3×1

16.1976 0.5000 16.1976

Each value in s gives the sensitivity of the corresponding eigenvalue towards a small change in A. If s(i) is very large it means that the matrix A is close to a matrix for which this i'th eigenvalue "is defective" (that is, that it has larger algebraic than geometric multiplicity).

[Bruno Luong]

"But does the form Vb\Db*Vb recover B?"

Actually Vb*Db/Vb recovers B

B=[2     1    -5
    0     1    -2
    0     0     3]

B = 3×3

2 1 -5 0 1 -2 0 0 3

[Vb,Db] = eig(B)

Vb = 3×3

1.0000 -0.7071 -0.9733 0 0.7071 -0.1622 0 0 0.1622

Db = 3×3

2 0 0 0 1 0 0 0 3

Br = Vb*Db/Vb

Br = 3×3

2.0000 1.0000 -5.0000 0 1.0000 -2.0000 0 0 3.0000

norm(Br - B)

ans = 6.2804e-16