Why $X^TX$ not always positive semidefinite using numpy.linalg.eig?

Question

We can use $z^T (X^TX)z = (Xz)^T(Xz)=||Xz||^2_2 \geq 0$ to prove $X^TX$ is positive semidefinite. However, when I use numpy.linalg.eig to compute the eigenvalues of dot product matrix, I cannot get all positive eigenvalues. How to explain for it?

import numpy as np
import math
np.random.seed(33)
X = np.random.randint(10, size=(2,6))
print(X)
Kernel_Matrix = np.dot(X.T,X)
eigenvalue, eigenvector = np.linalg.eig(Kernel_Matrix)
print(eigenvalue > 0)
print(eigenvalue)

Output:

[[4 7 8 2 2 9]
 [9 3 6 3 3 1]]
[[ 97  55  86  35  35  45]
 [ 55  58  74  23  23  66]
 [ 86  74 100  34  34  78]
 [ 35  23  34  13  13  21]
 [ 35  23  34  13  13  21]
 [ 45  66  78  21  21  82]]
[ True  True False False  True False]
[ 3.12680220e+02  5.03197805e+01 -9.39575362e-15 -1.44635182e-14
  1.43755791e-16 -7.87904824e-32]

Did you check whether the fourth eigenvalue was actually zero (which wouldn't spoil positive semidefiniteness)? — Ian, Oct 01 '20 at 14:53
@N.S. I have print all eigenvalues in new edit, the fourth eigenvalue is negative — Z Mario, Oct 01 '20 at 14:55
First, a quick note that you have been computing $XX^T$ rather than $X^TX$. Note that because $XX^T$ has rank $2$ and size $6$, it must have the eigenvalue $0$ with multiplicity $4$. It is likely that one of the $0$ eigenvalues is computed to be negative due to numerical error (i.e. rounding in floating point arithmetic) — Ben Grossmann, Oct 01 '20 at 14:56
@Ian, I have print all eigenvalues in new edit, the fourth eigenvalue is negative — Z Mario, Oct 01 '20 at 14:56
Indeed, if you use the statement print(eigenvalue > -1e-12), you should find that all eigenvalues are "positive up to a small margin of error". — Ben Grossmann, Oct 01 '20 at 15:00
@BenGrossmann , I use X = np.random.randint(10, size=(6,2)) and Kernel_Matrix = np.dot(X.T,X) to build $X^TX$ as a size 6 matrix. I still get same result. — Z Mario, Oct 01 '20 at 15:08
Also, try numpy.linalg.eigh. This is a specialized routine for complex-hermitian or real-symmetric matrices. This may provide a bit more accuracy, and will certainly suppress those imaginary numbers, to convince you that your matrix is semidefinite. — Charlie S, Oct 01 '20 at 15:09
@ZMario That doesn't make sense. $X^TX$ is a $2 \times 2$ matrix, not $6 \times 6$ — Ben Grossmann, Oct 01 '20 at 15:10
You can check by hand that the rank of $X^TX$ is 2. This means that 0 is an eigenvalue (of multiplicity 4) and it doesn't appear in your list. — N. S., Oct 01 '20 at 15:15

score 9 · Answer 1 · answered Oct 01 '20 at 15:00

You've printed out the result of computing 6 numbers using IEEE floating point arithmetic; you have not printed out the 6 eigenvalues of $X^t X$ for any real matrix $X$, because those are all nonnegative real numbers.

Notice that your "negative" fourth eigenvalue is smaller in absolute value than $10^{-14}$; that's pretty much the IEEE float world whispering in your ear "this is zero, but disguised by roundoff errors."

There's a difference between computers and mathematics, and it's good to know (and embrace) that as early as possible.

I recommend a read-though of Trefethen and Bao's section on floating-point arithmetic in Numerical Linear Algebra; it's wonderfully compact and informative.

score 2 · Accepted Answer · edited Oct 02 '20 at 13:51

As Ben Grossmann pointed out, for this matrix, it is actually easy to calculate by hand the eigenvalues and convince yourself that the error lies within the floating arithmetic.

First, you should calculate the rank of the matrix, and you'll see that it is 2 (I doublechecked it with some online software and the rank is indeed 2). Note that this can be calculated by hand, and it does not involve any approximation, while the numerical computation of eigenvalues of large matrices involves errors.

The above means that $\lambda=0$ is an eigenvalue of multiplicity 4 [you can even calculate the eigenspace]. Anyhow, at this point you can already see that there is a mistake with the output: 4 of your eigenvalues are 0

Then $$\det(xI-A)= (x^2+ax+b)x^4$$

You can easily figure out what $a,b$ are by simply plugging in two values for $x$ and solving. Or even better, you can recall that $a=-\operatorname{tr}(A)$ and then you only need to figure out what $b$ is.

Why $X^TX$ not always positive semidefinite using numpy.linalg.eig?

2 Answers2