Symbolic linear system solving not supported by ModelKit

[luanborelli]

Suppose we want to solve a system $Ax = b$, where $A$ is a polynomial matrix and $b$ is a vector of polynomials. ModelKit seems to be unable to solve such a system. Here's an example:

using HomotopyContinuation

@polyvar x y

A1 = x^2 + y^2
A2 = y + x^2 
A3 = x*y 
A4 = y^2 

b1 = 5*y - 4*x
b2 = -3*x - y 

A = [A1 A2; A3 A4]
b = [b1; b2]

A\b

Output:

ERROR: MethodError: no method matching MultivariatePolynomials.RationalPoly{DynamicPolynomials.Polynomial{true, Int64}, DynamicPolynomials.Polynomial{true, Int64}}(::MultivariatePolynomials.RationalPoly{DynamicPolynomials.Polynomial{true, Int64}, DynamicPolynomials.Polynomial{true, Int64}})

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Closest candidates are:
  MultivariatePolynomials.RationalPoly{NT, DT}(::Any, ::Any) where {NT<:(MultivariatePolynomials.AbstractPolynomialLike), DT<:(MultivariatePolynomials.AbstractPolynomialLike)} at D:\Users\b46525\.julia\packages\MultivariatePolynomials\sWAOE\src\rational.jl:5
  MultivariatePolynomials.RationalPoly{NT, DT}(::Bool) where {NT, DT} at D:\Users\b46525\.julia\packages\MultivariatePolynomials\sWAOE\src\rational.jl:10
Stacktrace:
 [1] oneunit(#unused#::Type{MultivariatePolynomials.RationalPoly{DynamicPolynomials.Polynomial{true, Int64}, DynamicPolynomials.Polynomial{true, Int64}}})
   @ Base .\number.jl:358
 [2] lutype(T::Type)
   @ LinearAlgebra D:\Users\b46525\AppData\Local\Programs\Julia-1.8.5\share\julia\stdlib\v1.8\LinearAlgebra\src\lu.jl:205
 [3] lu(A::Matrix{DynamicPolynomials.Polynomial{true, Int64}}, pivot::RowMaximum; check::Bool)
   @ LinearAlgebra D:\Users\b46525\AppData\Local\Programs\Julia-1.8.5\share\julia\stdlib\v1.8\LinearAlgebra\src\lu.jl:279
 [4] lu(A::Matrix{DynamicPolynomials.Polynomial{true, Int64}}, pivot::RowMaximum) (repeats 2 times)
   @ LinearAlgebra D:\Users\b46525\AppData\Local\Programs\Julia-1.8.5\share\julia\stdlib\v1.8\LinearAlgebra\src\lu.jl:278
 [5] \(A::Matrix{DynamicPolynomials.Polynomial{true, Int64}}, B::Vector{DynamicPolynomials.Polynomial{true, Int64}})
   @ LinearAlgebra D:\Users\b46525\AppData\Local\Programs\Julia-1.8.5\share\julia\stdlib\v1.8\LinearAlgebra\src\generic.jl:1110
 [6] top-level scope
   @ Untitled-2:16

This is a type of operation that other symbolic systems seem to handle with ease. I have already performed such an operation in Matlab, for example. In Julia, other systems like Symbolics (suggested in #456 for integration --- I endorse the suggestion!) are also capable of solving this kind of problem. Here's an example:

using Symbolics

@variables x y

A1 = x^2 + y^2
A2 = y + x^2 
A3 = x*y 
A4 = y^2 

b1 = 5*y - 4*x
b2 = -3*x - y 

A = [A1 A2; A3 A4]
b = [b1; b2]

A\b

Output:

2-element Vector{Num}:
 (5y + (-(y + x^2)*((-x*y*(5y - 4x)) / (x^2 + y^2) - 3x - y)) / (y^2 + (-x*y*(y + x^2)) / (x^2 + y^2)) - 4x) / (x^2 + y^2)
                                        ((-x*y*(5y - 4x)) / (x^2 + y^2) - 3x - y) / (y^2 + (-x*y*(y + x^2)) / (x^2 + y^2))

Obs.: I have a rather complicated system that depends on this type of operation to be built. I need to solve it using homotopy and I'm trying to rebuild it so that it's ModelKit compatible, but I'm stuck on this operation. Any suggestions on how to get around this problem (at least temporarily) would also be greatly appreciated. I have been trying to parse Symbolics expressions to ModelKit expressions, but without success so far. Thanks.

[saschatimme]

ModelKit is not intended as a full blown symbolic algebra package and does therefore not support many symbolic algorithms.

I have been trying to parse Symbolics expressions to ModelKit expressions, but without success so far. Thanks.

Maybe you can use a similar approach as used here #520

Personally, I would try to avoid using a direct symbolic solution to a linear system in any problem formulation since the expression will be numerically very unstable. You can introduce additionally variable and instead if f - A^{-1}b=0 you can solve the equivalent problem [f - y; A * y] = [0; b].

[luanborelli]

Thanks for your reply, @saschatimme. I understand that ModelKit is not intended to be a complete symbolic algebra package and that my problem is quite specific. I found a way to port my system from Mathematica to ModelKit and at least temporarily it solves my problem.