# Is there a bug in large matrix inverse?

**URL:** <https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787>\
**Category:** General Usage\
**Tags:** question\
**Created:** [May 19, 2020, 11:35pm UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787 "2020-05-19T23:35:41Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![zxjroger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zxjroger/32/11826_2.png) [@zxjroger](https://discourse.julialang.org/u/zxjroger)\
**Post date:** [May 19, 2020, 11:35pm UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/1 "2020-05-19T23:35:41Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [May 19, 2020, 11:41pm UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/2 "2020-05-19T23:41:26Z")

</div>

Doesn’t your matrix have a rank of 1?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 19, 2020, 11:57pm UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/3 "2020-05-19T23:57:32Z")

</div>

@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
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

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

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [May 20, 2020, 12:11am UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/4 "2020-05-20T00:11:07Z")

</div>

```julia
julia> a = randn(100);
julia> b = a * transpose(a);
julia> using LinearAlgebra

julia> rank(b)
1

```

This matrix is not invertible.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 20, 2020, 1:17am UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/5 "2020-05-20T01:17:07Z")

</div>

> [@Oscar\_Smith](#):
>
> The second is that computing matrix inverses is an inherently numerically unstable operation.

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
julia> a = randn(100);

julia> cond(a * a')
5.413008718471948e21

```

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

---

<div class="post-metadata">

**Author:** ![Pbellive](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pbellive/32/3604_2.png) [@Pbellive](https://discourse.julialang.org/u/Pbellive)\
**Post date:** [May 20, 2020, 2:19am UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/6 "2020-05-20T02:19:51Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![zxjroger](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zxjroger/32/11826_2.png) [@zxjroger](https://discourse.julialang.org/u/zxjroger)\
**Post date:** [May 20, 2020, 10:33pm UTC](https://discourse.julialang.org/t/is-there-a-bug-in-large-matrix-inverse/39787/8 "2020-05-20T22:33:18Z")

</div>

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,

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

pinv(b) * b

```

I got error

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

```