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

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

``````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> B = inv(A) * A;

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

true

``````

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

Using the Caley-Hamilton method (Wikipedia):

``````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
``````
1 Like

A side comment: last time I checked, quite a few trigonometric identities where missing in Symbolics.jl. I tried to make a PR but it got stuck. Feel free to revive it if that’s helpful.

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!

2 Likes