Inversion of anti-symmetric matrix fails

[michakraus]

Inversion of (even-dimensional) anti-symmetric matrices fails:

julia> inv([ Num(0)  -Num(1)
            +Num(1)   Num(0)])
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] checknonsingular at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:20 [inlined]
 [2] generic_lufact!(::Array{Num,2}, ::Val{false}; check::Bool) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:178
 [3] #lu!#134 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:130 [inlined]
 [4] #lu#136 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:273 [inlined]
 [5] lu at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:272 [inlined]
 [6] #lu#5 at ~/.julia/packages/ModelingToolkit/QPQyE/src/ModelingToolkit.jl:126 [inlined]
 [7] lu at ~/.julia/packages/ModelingToolkit/QPQyE/src/ModelingToolkit.jl:126 [inlined]
 [8] inv(::Array{Num,2}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/dense.jl:781
 [9] top-level scope at REPL[4]:1

while

julia> inv([ 0  -1
            +1   0])
2×2 Array{Float64,2}:
  0.0  1.0
 -1.0  0.0

[michakraus]

I just found that solve_for (which would be all I really need for now) almost works:

(α,) = @variables α[1:2]
(β,) = @variables β[1:2]

Ω = [ Num(0)  -Num(1)
     +Num(1)   Num(0)]

lhs = simplify.(Ω * α)
equs = [Equation(lhs[i], β[i]) for i in eachindex(lhs,β)]

ModelingToolkit.solve_for(equs, α)

returns

2-element Array{SymbolicUtils.Symbolic{Real},1}:
 β₂
 -1 * β₁

which is fine. However, if the matrix Ω is non-constant (which it is in my application), I get an error, namely

@parameters t
@variables x₁(t), x₂(t)

Ω = [ Num(0)  -x₁*x₂
     +x₁*x₂   Num(0)]

lhs = simplify.(Ω * α)
equs = [Equation(lhs[i], β[i]) for i in eachindex(lhs,β)]

ModelingToolkit.solve_for(equs, α)

throws

TypeError: non-boolean (Term{Bool}) used in boolean context

Stacktrace:
 [1] (::ModelingToolkit.var"#40#41"{Array{Any,2},Int64})(::Int64) at ~/.julia/packages/ModelingToolkit/QPQyE/src/solve.jl:22
 [2] iterate at ./generator.jl:47 [inlined]
 [3] collect_to!(::Array{Float64,1}, ::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#40#41"{Array{Any,2},Int64}}, ::Int64, ::Int64) at ./array.jl:732
 [4] collect_to_with_first!(::Array{Float64,1}, ::Float64, ::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#40#41"{Array{Any,2},Int64}}, ::Int64) at ./array.jl:710
 [5] _collect(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#40#41"{Array{Any,2},Int64}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:704
 [6] collect_similar at ./array.jl:628 [inlined]
 [7] map at ./abstractarray.jl:2162 [inlined]
 [8] sym_lu(::Array{Any,2}) at ~/.julia/packages/ModelingToolkit/QPQyE/src/solve.jl:22
 [9] _solve(::Array{Num,2}, ::Array{Num,1}) at ~/.julia/packages/ModelingToolkit/QPQyE/src/solve.jl:75
 [10] solve_for(::Array{Equation,1}, ::Array{Num,1}) at ~/.julia/packages/ModelingToolkit/QPQyE/src/solve.jl:69
 [11] top-level scope at In[343]:1

The corresponding line in solve.jl is

_, i = findmin(map(ii->iszero(U[ii, k]) ? Inf : nterms(U[ii,k]), k:n))

[YingboMa]

This is fixed by #683

[michakraus]

Thanks, solve_for does indeed work now for a non-constant anti-symmetric matrix. However, with MTK v4.5.0 I still get the ZeroPivotException for

> inv(Num.([ 0  -1
             1   0 ]))

and the same for the non-constant matrix from above.

[YingboMa]

I cannot reproduce that. Did you restart Julia after upgrading?

See

julia> (α,) = @variables α[1:2]
(Num[α₁, α₂],)

julia> (β,) = @variables β[1:2]
(Num[β₁, β₂],)

julia> @parameters t
(t,)

julia> @variables x₁(t), x₂(t)
(x₁(t), x₂(t))

julia> Ω = [ Num(0)  -x₁*x₂
            +x₁*x₂   Num(0)]
2×2 Matrix{Num}:
           0  -1x₁(t)*x₂(t)
 x₁(t)*x₂(t)              0

julia> lhs = simplify.(Ω * α)
2-element Vector{Num}:
 -1x₁(t)*x₂(t)*α₂
   x₁(t)*x₂(t)*α₁

julia> equs = [Equation(lhs[i], β[i]) for i in eachindex(lhs,β)]
2-element Vector{Equation}:
 Equation(-1x₁(t)*x₂(t)*α₂, β₁)
 Equation(x₁(t)*x₂(t)*α₁, β₂)

julia> ModelingToolkit.solve_for(equs, α)
2-element Vector{SymbolicUtils.Mul{Number, Int64, Dict{Any, Number}}}:
 (x₁(t)^-1)*(x₂(t)^-1)*β₂
 -1(x₁(t)^-1)*(x₂(t)^-1)*β₁

and

julia> inv(Num.([ 0  -1
                    1   0 ]))
2×2 Matrix{Num}:
  0.0   1.0
 -1.0  -0.0

[michakraus]

Ha, there is something strange going on. Indeed I restarted, and I ran in a Jupyter notebook. If I run in the REPL, everything works, but in the Jupyter Notebook I get a ZeroPivotException. Both running Julia v1.5.3, both freshly started... solve_for works in both, but not inv.

Addendum I: I guess the notebook is not using the anticipated environment, but I don't see why... this needs some more investigation.

Addendum II: I found the issue. Works in both REPL and notebook now. Sorry for the noise.

[michakraus]

I found another issue with the 4d version of the above example:

julia> @parameters t
(t,)

julia> @variables x₁(t), x₂(t), x₃(t), x₄(t)
(x₁(t), x₂(t), x₃(t), x₄(t))

julia> Ω = [ 0        -1/x₁/x₂  +1/x₁/x₃   -1/x₁/x₄
            +1/x₂/x₁   0        -1/x₂/x₃   +1/x₂/x₄
            -1/x₃/x₁  +1/x₃/x₂   0         -1/x₃/x₄
            +1/x₄/x₁  -1/x₄/x₂  +1/x₄/x₃    0      ]
4×4 Matrix{Num}:
                                                                                                                                                                0  …  (-1 / x₁(t)) / x₄(t)
  (1 / x₂(t)) / x₁(t)      (1 / x₂(t)) / x₄(t)
 (-1 / x₃(t)) / x₁(t)     (-1 / x₃(t)) / x₄(t)
  (1 / x₄(t)) / x₁(t)                                                                                                                                                                    0

julia> P = simplify.(inv(Ω))
4×4 Array{Num,2}:
(0 ^ -1) * (x₁(t) ^ 2) * (x₃(t) ^ -2)                                                                                                                                                                                                                                                                                                                                                                                                                                                  …  x₁(t) * x₄(t) * (x₃(t) ^ -2)
                                                                                                                                                                                                                                                                                                   -1.0 * x₁(t) * x₂(t) * (1 + (-1 * (0 ^ -1) * (x₃(t) ^ -2)))                                                                                                                                                                                                                              0.0
                                                         x₁(t) * (0 ^ -1) * (x₃(t) ^ -1)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               x₄(t) * (x₃(t) ^ -1)
                                                                                                                                                                                                                                                                                                        -0.0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            0.0

The problem being the (0 ^ -1)'s, which should be 0. If I didn't introduce some stupid typo, that thing should be invertible.

Addendum: I verified with SymEngine which produces the correct result:

julia> (x₁, x₂, x₃, x₄) = symbols("x₁, x₂, x₃, x₄")
(x₁, x₂, x₃, x₄)

julia> Ω = [ 0        -1/x₁/x₂  +1/x₁/x₃   -1/x₁/x₄
            +1/x₂/x₁   0        -1/x₂/x₃   +1/x₂/x₄
            -1/x₃/x₁  +1/x₃/x₂   0         -1/x₃/x₄
            +1/x₄/x₁  -1/x₄/x₂  +1/x₄/x₃    0      ]
4×4 Array{Basic,2}:
          0  -1/(x₂*x₁)   1/(x₃*x₁)  -1/(x₄*x₁)
  1/(x₂*x₁)           0  -1/(x₂*x₃)   1/(x₄*x₂)
 -1/(x₃*x₁)   1/(x₂*x₃)           0  -1/(x₄*x₃)
  1/(x₄*x₁)  -1/(x₄*x₂)   1/(x₄*x₃)           0

julia> expand.(inv(Ω))
4×4 Array{Basic,2}:
      0   x₂*x₁   x₃*x₁  x₄*x₁
 -x₂*x₁       0   x₂*x₃  x₄*x₂
 -x₃*x₁  -x₂*x₃       0  x₄*x₃
 -x₄*x₁  -x₄*x₂  -x₄*x₃      0

Tested with ModelingToolkit v4.5.0 on Julia v1.6-beta1.

[YingboMa]

The master branch of MTK gives

julia> P = simplify.(inv(Ω))
4×4 Matrix{Num}:
                   0          x₁(t)*x₂(t)          x₁(t)*x₃(t)  x₁(t)*x₄(t)
 (-1//1)*x₁(t)*x₂(t)                    0          x₂(t)*x₃(t)  x₂(t)*x₄(t)
 (-1//1)*x₁(t)*x₃(t)  (-1//1)*x₂(t)*x₃(t)                    0  x₃(t)*x₄(t)
 (-1//1)*x₁(t)*x₄(t)  (-1//1)*x₂(t)*x₄(t)  (-1//1)*x₃(t)*x₄(t)            0

[ChrisRackauckas]

(-1//1) intentional?

[YingboMa]

Yes, it's deliberate. SymbolicUtils.jl's canonicalization doesn't have operator -. There's only -1*x

Also, this problem cannot be solved in general, since it amounts to test whether a set of arbitrary expressions is 0 under +, -, /, *.

[ChrisRackauckas]

But I would assume // can get simplified away when it has a 1 or a zero? That looks like a missing rule.

[YingboMa]

Yeah, that can be done.