Is there a bug in large matrix inverse?

General Usage
1

I try to find the inverse of a symmetric matrix. For example,

a = randn(100)
b = a * transpose(a)
inv(b) * b

Here, matrix b is symmetric. However, inv(b) is not the correct inverse of matrix b. It is not symmetric and the product in the last line is far away from identity matrix. Is there anything wrong in my code?

2

Doesn’t your matrix have a rank of 1?

3

@zxjroger there are 2 problems here. The first is that although you know b is symmetric, Julia doesn’t. The second is that computing matrix inverses is an inherently numerically unstable operation.
However, if you look at how far inv(b) is from symmetric or from being an inverse, you will see that it is actually pretty close.

julia> norm(inv(b)*b-I)
4.7645273464508936e-8
julia> norm(inv(b) - inv(b)')
4.1066007247761773e-7

If we instead tell Julia that b is symmetric, we see the following

b_sym = Symmetric(b)
julia> norm(inv(b_sym)*b_sym-I)
5.868028292184572e-8

julia> norm(inv(b_sym) - inv(b_sym)')
0.0
1 Like
4 
julia> a = randn(100);
julia> b = a * transpose(a);
julia> using LinearAlgebra

julia> rank(b)
1

This matrix is not invertible.

5

Computing matrix inverses is numerically stable, in the sense that it is backwards stable. This doesn’t mean that it is always accurate, however — in particular, it can give a very inaccurate result if the matrix is badly conditioned (close to singular).

And a rank-1 matrix a * a' is very badly conditioned indeed! Because of roundoff errors, it is not exactly singular, but the condition number is still huge:

julia> a = randn(100);

julia> cond(a * a')
5.413008718471948e21

This means that the result of inverting it is basically garbage from roundoff errors.

4 Likes
6

I wonder if the OP was expecting randn(100) to produce a 100 X 100 matrix, as it does in matlab, as opposed to a length 100 vector, as it does in julia. As @PetrKryslUCSD and @stevengj point out, if a is a length 100 vector, then a*a’ is a 100 X 100 rank 1 matrix. But if we do

A = randn(100,100)
B = A*transpose(A)

then B will likely be a full rank 100 X 100 matrix. Of course @stevengj’s response is entirely correct. Just trying to clear up one more potential source of confusion.

5 Likes
8

Thanks to all. I tried to use pseudo inverse instead. Is there a function in Julia which computes the inverse of a symmetric and singular matrix? Currently, if I tell Julia that my matrix is symmetric, it will not be able to compute the inverse by pinv(). For example,

a = randn(100)
b = a * transpose(a)
b = Symmetric(b)

pinv(b) * b

I got error

ERROR: MethodError: no method matching svd(::Symmetric{Float64,Array{Float64,2}}; full=false)