Matrix multiplied by inverse matrix does not reproduce identity matrix

General Usage
,
1

Hello,
I found that multiplying a matrix by its inverse does not yield an identity matrix. Please see below a MWE:

matrix = [5.93419 4.8067; 4.8067 5.93419]
invmatrix = inv(matrix)
matrix*invmatrix

The last line generates the following output:

2×2 Matrix{Float64}:
1.0 4.44089e-16
4.44089e-16 1.0

Is there any way to make the off-diagonal elements zero instead of some value that is very close to zero? Thanks!

2

One solution is posted here, for display purposes.

1 Like
3

Thanks! That solves the previous MWE, but does not solve the issue I am having. Please see below a new MWE with the method you suggested:

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.4f", f)
Σzz = [5.9342  5.3408  4.8067  4.3260; 5.3408  5.9342  5.3408  4.8067; 4.8067  5.3408  5.9342  5.3408; 4.3260  4.8067  5.3408  5.9342]
Σzz⁺ = [5.3408  4.8067  4.3260; 5.9342  5.3408  4.8067; 5.3408  5.9342  5.3408; 4.8067  5.3408  5.9342]
Σz⁺z⁺ = [5.9342  5.3408  4.8067; 5.3408  5.9342  5.3408; 4.8067  5.3408  5.9342]
iΣz⁺z⁺ = inv(Σz⁺z⁺)
tΣzz⁺ = transpose(Σzz⁺)
varc .= Σzz - Σzz⁺ * iΣz⁺z⁺ * tΣzz⁺
issymmetric(varc)

Σzz⁺ and Σz⁺z⁺ are just partitions of Σzz.
varc is supposed to be symmetric, but the last line of code shows that it is not.
Is there any way to make varc symmetric?
Thank you!

4

It’s nearly symmetric. Because of the imprecision inherent in floating point numbers, it’s best to use approximate comparisons like this:

julia> varc ≈ varc'
true

The approximate equality symbol can be typed in the REPL with \approx<TAB>, or you can use the equivalent function call: isapprox(varc, varc'). If you know your matrix is structurally symmetric, you can use the Symmetric wrapper type:

julia> Symmetric(varc)
4×4 Symmetric{Float64, Matrix{Float64}}:
  1.1275  -0.0000  -0.0000  -0.0000
 -0.0000  -0.0000  -0.0000  -0.0000
 -0.0000  -0.0000  -0.0000  -0.0000
 -0.0000  -0.0000  -0.0000  -0.0000

julia> issymmetric(ans)
true
5 Likes
5

You can zero them out manually, eg as in

M = matrix*invmatrix;
M[abs.(M) .≤ 1e-10] .= 0;
M

but this is hard to do in a disciplined manner in general.

1 Like