# Is there a bug in large matrix inverse?

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?

Doesn’t your matrix have a rank of 1?

@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
``````julia> a = randn(100);
julia> b = a * transpose(a);
julia> using LinearAlgebra

julia> rank(b)
1
``````

This matrix is not invertible.

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

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

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)
``````