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