# Symbolics.jl: Why doesn't a matrix multiplied by its inverse simplify to the identity?

**URL:** <https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978>\
**Category:** Specific Domains\
**Tags:** symbolics\
**Created:** [August 7, 2021, 3:46am UTC](https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978 "2021-08-07T03:46:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![garrek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/garrek/32/27937_2.png) [@garrek](https://discourse.julialang.org/u/garrek)\
**Post date:** [August 7, 2021, 3:46am UTC](https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978/1 "2021-08-07T03:46:13Z")

</div>

**Expected behavior** : Using `simplify` in the Symbolics.jl package, return the identity for the operation `A^-1 * A` for a square matrix involving multiplication of trigonometric functions in the matrix elements (assuming `det(A) != 0`).

**Actual behavior** : Using `simplify` does not reduce the final expression to the identity, but substituting numerical values at the end of symbolic calculations does evaluate to the identity matrix.

I think the issue has to do with how `lu()` is computed, but I’m not sure what it’s doing or how to fix it.  
Is there a preferred way to do these kinds of operations? I’ve gone through the [Symbolics.jl](https://symbolics.juliasymbolics.org/dev/) documentation and can’t find much about matrices that include trigonometric functions as elements.

Here’s a simple case with the 2x2 rotation matrix (I’ve tried other variations, too):

```julia
julia> using LinearAlgebra, Symbolics, SymbolicUtils

julia> @variables x::Real;
julia> I2 = [1 0; 0 1]; #identity matrix

julia> A = [cos(x) -sin(x); sin(x) cos(x)]
2×2 Matrix{Num}:
 cos(x) -sin(x)
 sin(x) cos(x)

julia> simplify.(inv(A) * A)
2×2 Matrix{Num}:
 1 … sin(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) - (sin(x)*(true - ((cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)))*(cos(x)^-1))
 0 cos(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) + (cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)

```

By eye, the result should simplify to `I2` (except for those `true` values). Trying to simplify `inv(A)` first does not work. Directly accessing the LU factorization method and using left-division instead of `inv` also does not work (i.e. solving `A^-1 = A / I`. See [Base.inv](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.inv-Tuple%7BAbstractMatrix%7BT%7D%20where%20T%7D)):

```julia
julia> simplify.(lu(A) \ I2 * A)
2×2 Matrix{Num}:
 1 … sin(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) - (sin(x)*(1 - ((cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)))*(cos(x)^-1))
 0 cos(x)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1) + (cos(x)^-1)*(sin(x)^2)*((cos(x) + (cos(x)^-1)*(sin(x)^2))^-1)

```

Again, it evaluates the expression correctly, but can’t simplify further (even though it should be able to).  
When setting `x = π/4` (or any other value) using `substitute`, the expression correctly reduces to the identity matrix, but some of the elements are converted to floats:

```julia
julia> B = inv(A) * A;

julia> answer = simplify.(substitute.(B, (Dict(x => π/4),)))
2×2 Matrix{Num}:
 1 0.0
 0 1.0

julia> answer == I2
true

```

I found one workaround, bypassing `inv` altogether, but it only works for 2x2 matrices.

Using the Caley-Hamilton method ([Wikipedia](https://en.wikipedia.org/wiki/Invertible_matrix#Analytic_solution)):

```julia
julia> A_inv = det(A) * (tr(A) * I2 - A)
2×2 Matrix{Num}:
  cos(x)*(cos(x)^2 + sin(x)^2) sin(x)*(cos(x)^2 + sin(x)^2)
 -sin(x)*(cos(x)^2 + sin(x)^2) cos(x)*(cos(x)^2 + sin(x)^2)

julia> simplify.(A_inv * A)
2×2 Matrix{Num}:
 1 0
 0 1

julia> simplify.(A_inv * A) == I2
true

```

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [August 7, 2021, 8:58am UTC](https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978/2 "2021-08-07T08:58:33Z")

</div>

A side comment: last time I checked, quite a few trigonometric identities where missing in Symbolics.jl. I tried to make a [PR](https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/158) but it got stuck. Feel free to revive it if that’s helpful.

---

<div class="post-metadata">

**Author:** ![garrek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/garrek/32/27937_2.png) [@garrek](https://discourse.julialang.org/u/garrek)\
**Post date:** [August 8, 2021, 8:40am UTC](https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978/3 "2021-08-08T08:40:37Z")

</div>

I actually ran into your pull request when I was trying to figure out trig identities in Symbolics.jl. Today I cloned your branch onto my machine and was able to reproduce the inconsistencies from your March 6 post. I also ran the tests (copied and pasted into rulesets.jl) you mentioned in that post and the same tests passed/failed. I’ll revive your PR and do a little tinkering.

I just started learning Julia (I use Python in my research) so I’m not sure how much progress I can make on this, but I’ll give it a try!

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [February 20, 2023, 7:16pm UTC](https://discourse.julialang.org/t/symbolics-jl-why-doesnt-a-matrix-multiplied-by-its-inverse-simplify-to-the-identity/65978/4 "2023-02-20T19:16:36Z")

</div>

Hi all, this is still an issue even for examples without trigonometric  
functions. Is there an answer from the developers how to get the correct answer?

```julia
@variables x::Real y::Real z::Real
a = [x^2+y^2, x^2+z^2, z^2 + y^2, x^2+y^2+z^2]
b = reshape(a,2,2)
simplify(b*inv(b), expand=true)
2×2 Matrix{Num}:
     1 … 0
 ((true + ((y^2 + z^2)*((x^2 + z^2) / (x^2 + y^2))) / (x^2 + y^2 + z^2 + (-(x^2 + z^2)*(y^2 + z^2)) / (x^2 + y^2)))*(x^2 + z^2)) / (x^2 + y^2) + (-(x^2 + y^2 + z^2)*((x^2 + z^2) / (x^2 + y^2))) / (x^2 + y^2 + z^2 + (-(x^2 + z^2)*(y^2 + z^2)) / (x^2 + y^2)) (x^2 + y^2 + z^2) / (x^2 + y^2 + z^2 + (-(x^2 + z^2)*(y^2 + z^2)) / (x^2 + y^2)) + (-(x^2 + z^2)*((y^2 + z^2) / (x^2 + y^2 + z^2 + (-(x^2 + z^2)*(y^2 + z^2)) / (x^2 + y^2)))) / (x^2 + y^2)

```