# Why same calculations on matrices yield different results?

**URL:** <https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573>\
**Category:** General Usage\
**Tags:** question\
**Created:** [March 27, 2022, 9:58pm UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573 "2022-03-27T21:58:07Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![csubhodeep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/csubhodeep/32/34967_2.png) [@csubhodeep](https://discourse.julialang.org/u/csubhodeep)\
**Post date:** [March 27, 2022, 9:58pm UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573/1 "2022-03-27T21:58:07Z")

</div>

Hello all,

I am fairly new to Julia but I come from the Python world. While I am still learning the basic machinery of the `Matrix`. I noticed a small difference in the way a “self-declared identity” matrix is being used in the following matrix operations:

```julia
julia> self_defined_I = [1. 0.; 0. 1.]
2×2 Matrix{Float64}:
 1.0 0.0
 0.0 1.0

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1 2
 3 4

julia> inv(x)
2×2 Matrix{Float64}:
 -2.0 1.0
  1.5 -0.5

julia> x^(-1)
2×2 Matrix{Float64}:
 -2.0 1.0
  1.5 -0.5

julia> self_defined_I/x
2×2 Matrix{Float64}:
 -2.0 1.0
  1.5 -0.5

julia> x\self_defined_I
2×2 Matrix{Float64}:
 -2.0 1.0
  1.5 -0.5

julia> inv(x) == x^(-1) == self_defined_I/x == x\self_defined_I
false

julia> inv(x) == x^(-1) == x\self_defined_I
true

julia> x\self_defined_I == self_defined_I/x
false

```

However, when using the `isapprox` operator we get:

```julia
julia> x\self_defined_I ≈ self_defined_I/x
true

```

Also, when using the `LinearAlgebra.I` identity matrix such kind of inconsistencies don’t happen even with the exact `==` operator:

```julia
using LinearAlgebra

julia> LinearAlgebra.I == self_defined_I
true

julia> inv(x) == x^(-1) == x\self_defined_I == LinearAlgebra.I/x == x\LinearAlgebra.I
true

```

I can see that all the results are of type `Matrix{Float64}` then why specifically the `self_defined_I/x` operation results in a matrix that has a different precision?

```julia
julia> self_defined_I/x - inv(x)
2×2 Matrix{Float64}:
 -4.44089e-16 2.22045e-16
  2.22045e-16 -1.11022e-16

julia> x\self_defined_I - inv(x)
2×2 Matrix{Float64}:
 0.0 0.0
 0.0 0.0

```

I would also like to ask if it is something one should take care or be worried about and if this could have the potential to propagate in other kinds of matrix operations and cause significant impact in accuracy?

Another off category question - from a design perspective does it make sense to force the output of `x*inv(x)` to `LinearAlgebra.I`?

P.S. Following are the runtime details:

```julia
julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-10210U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 27, 2022, 11:37pm UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573/2 "2022-03-27T23:37:24Z")

</div>

To start with, it doesn’t really make sense to compare floats with `==`. The difference you’re seeing is on the order of `eps`, so the only way it changes the result of a calculation is if that calculation is chaotic.

And you don’t want to be inverting matrices anyway, use `\` directly or save an LU decomposition if you need to repeat it many times.

`LinearAlgebra.I` is of a specific type, so the compiler knows that `anything\I = inv(anything)`. It can’t know that about your home rolled `I` unless you make that information available via the type system.

It’s a similar story for `x*inv(x)`, the compiler cannot know that `inv(x)` is related to `x` unless you somehow made `inv(x)` a lazy type that upon multiplication with `y` checked if `x ≈ y`, which seems overkill since most of the time it will be false.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [March 28, 2022, 6:39am UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573/3 "2022-03-28T06:39:31Z")

</div>

> [@csubhodeep](#):
>
> I am fairly new to Julia but I come from the Python world.

I’m left wondering how and why [floating point behavior](https://www.itu.dk/~sestoft/bachelor/IEEE754_article.pdf) in Python is different from Julia? AFAIU `I` in Julia is a proxy for an identity matrix and dynamic dispatch can resolve `I * A = A` directly without using matrix multiplication.

Edit: granted, [`numpy.eye`](https://numpy.org/devdocs/reference/generated/numpy.eye.html) seems to generate the full matrix…

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [March 28, 2022, 7:39am UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573/4 "2022-03-28T07:39:14Z")

</div>

I doubt there is a real difference in the floating point behavior. It’s probably rather that Python doesn’t have `UniformScaling`, the type of `LinearAlgebra.I`. When OP rolls their own `I`, they remove this information from the type system, so dispatch _does_ need to perform the multiplications.

---

<div class="post-metadata">

**Author:** ![csubhodeep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/csubhodeep/32/34967_2.png) [@csubhodeep](https://discourse.julialang.org/u/csubhodeep)\
**Post date:** [March 28, 2022, 8:58pm UTC](https://discourse.julialang.org/t/why-same-calculations-on-matrices-yield-different-results/78573/5 "2022-03-28T20:58:22Z")

</div>

Thanks a lot for the answers !