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

Specific Domains
1

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

julia> answer == I2
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
2 Likes
2

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.

3

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!

3 Likes
4

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?

@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)
1 Like