add support for solving arbitrary linear systems

[Suavesito-Olimpiada]

This adds the ability to solve arbitrary consistent linear systems. It throws if an inconsistent system is passed (e.j. $a+b=2$ and $a+b=1$ on $(a,b)$ ), otherwise it returns a vector with the solutions.

It uses lu factorization for square rank-complete systems; for under-determined systems it uses lu factorization first, with a rref based solver later, and finally rref directly for over-determined systems.

Usage example

julia> vars = @variables x y z
3-element Vector{Num}:
 x
 y
 z

julia> eqs = [x + y + z ~ 3,
              x + y ~ 2]
2-element Vector{Equation}:
 x + y + z ~ 3
 x + y ~ 2

julia> Symbolics.solve_for(eqs, vars) # `ℝ` express that `y` is a free variable, maybe just put `y`
3-element Vector{Any}:
  2.0 - y
  ℝ
 1.0

There is documentation to add and amend. I would like to add efficient sparse algorithms [1] for getting to rref, but I need to read more about this (maybe here?). There is also the problem that it is slow, mainly because of the need for _sym_urref! and dynamic dispatch for every *(::Num, ::Num) and -(::Num, ::Num) call.

(1) In Julia v1.9, SparseMatrix{Num} just works, finally. 🥳

This is what a ProfileView of Symbolics.solve_for looks like for a system of $91 \times 190$.

julia> length(eqs)
91

julia> length(vars)
190

julia> @time Symbolics.solve_for(eqs, vars);
  1.069665 seconds (5.90 M allocations: 148.025 MiB, 2.76% gc time)

You can notice that 8-11 and 14-17 are *(::Num, ::Num) or -(::Num, ::Num), marked red for being dynamic dispatch.

 1. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:271, MethodInstance for Symbolics.solve_for(::Vector{Equation}, ::Vector{Num})
 2. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:274, MethodInstance for Symbolics.var"#solve_for#328"(::Bool, ::Bool, ::typeof(Symbolics.solve_for), ::Vector{Equation}, ::Vector{Num})
 3. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:430, MethodInstance for Symbolics.linear_expansion(::Vector{Equation}, ::Vector{Num})
 4. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:279, MethodInstance for Symbolics.var"#solve_for#328"(::Bool, ::Bool, ::typeof(Symbolics.solve_for), ::Vector{Equation}, ::Vector{Num})
 5. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:300, MethodInstance for Symbolics._solve(::Matrix{Num}, ::Vector{Num}, ::Bool, ::Vector{Num})
 6. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:220, _factorize [inlined]
 7. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:74, MethodInstance for Symbolics._sym_lu(::Matrix{Num})
 8. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:53, MethodInstance for *(::Num, ::Num)
 9. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:67, MethodInstance for *(::Num, ::Num)
10. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:53, MethodInstance for -(::Num, ::Num)
11. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:67, MethodInstance for -(::Num, ::Num)
12. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:323, MethodInstance for Symbolics._solve(::Matrix{Num}, ::Vector{Num}, ::Bool, ::Vector{Num})
13. /media/data/Programming/Langs/Julia/Symbolics.jl/src/linear_algebra.jl:200, MethodInstance for Symbolics._sym_urref!(::Matrix{Num}, ::Vector{Num}, ::Vector{Int64})
14. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:53, MethodInstance for *(::Num, ::Num)
15. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:67, MethodInstance for *(::Num, ::Num)
16. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:53, MethodInstance for -(::Num, ::Num)
17. /home/jose/.julia/packages/SymbolicUtils/qulQp/src/methods.jl:67, MethodInstance for -(::Num, ::Num)

[codecov-commenter]

Codecov Report

Merging #757 (21673e5) into master (08ebc48) will decrease coverage by 2.13%. The diff coverage is 35.03%.

:mega: This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##           master     #757      +/-   ##
==========================================
- Coverage   76.84%   74.71%   -2.13%     
==========================================
  Files          26       26              
  Lines        3239     3370     +131     
==========================================
+ Hits         2489     2518      +29     
- Misses        750      852     +102     

Impacted Files	Coverage Δ
src/linear_algebra.jl	64.24% <35.03%> (-25.90%)	:arrow_down:

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[Suavesito-Olimpiada]

I don't think the problem in CI is related to this PR. I think code-wise, it is ready (modulus review). It just needs documentation and testing (specially the second).

[Suavesito-Olimpiada]

Bump

[ChrisRackauckas]

@shashi

[shashi]

@Suavesito-Olimpiada did you ever try to profile it after the Unityper update? I'm running CI and will merge after that.

[shashi]

The tests don't seem to hit a lot of cases... Could you add some more, for example to chatch that it errors for inconsistent systems?

[Suavesito-Olimpiada]

I haven't tried to profile it again, but I'm doing it. I'll add more test.

[Suavesito-Olimpiada]

I have profiled it again, here are the results @shashi.

julia> length(eqs)
91

julia> length(vars)
190

julia> @time Symbolics.solve_for(eqs, vars);
  1.221608 seconds (6.97 M allocations: 186.474 MiB, 5.23% gc time)

1.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:268, MethodInstance for Symbolics.solve_for(::Vector{Equation}, ::Vector{Num})
2.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:271, MethodInstance for Symbolics.var"#solve_for#341"(::Bool, ::Bool, ::typeof(Symbolics.solve_for), ::Vector{Equation}, ::Vector{Num})
3.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:427, MethodInstance for Symbolics.linear_expansion(::Vector{Equation}, ::Vector{Num})
4.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:276, MethodInstance for Symbolics.var"#solve_for#341"(::Bool, ::Bool, ::typeof(Symbolics.solve_for), ::Vector{Equation}, ::Vector{Num})
5.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:297, MethodInstance for Symbolics._solve(::Matrix{Num}, ::Vector{Num}, ::Vector{Num}, ::Bool)
6.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:217, _factorize [inlined]
7.  ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:77, MethodInstance for Symbolics._sym_lu(::Matrix{Num})
8.  ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:54, MethodInstance for *(::Num, ::Num)
9.  ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:68, MethodInstance for *(::Num, ::Num)
10. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:54, MethodInstance for -(::Num, ::Num)
11. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:68, MethodInstance for -(::Num, ::Num)
12. ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:320, MethodInstance for Symbolics._solve(::Matrix{Num}, ::Vector{Num}, ::Vector{Num}, ::Bool)
13. ~/.julia/packages/Symbolics/ThKGL/src/linear_algebra.jl:197, MethodInstance for Symbolics._sym_urref!(::Matrix{Num}, ::Vector{Num}, ::Vector{Int64})
14. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:54, MethodInstance for *(::Num, ::Num)
15. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:68, MethodInstance for *(::Num, ::Num)
16. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:54, MethodInstance for -(::Num, ::Num)
17. ~/.julia/packages/SymbolicUtils/1JRDc/src/methods.jl:68, MethodInstance for -(::Num, ::Num)

Where I installed Symbolics.jl doing (test) pkg> add https://github.com/Suavesito-Olimpiada/Symbolics.jl#suave/linalg.

[dom-linkevicius]

Hi @Suavesito-Olimpiada I was wondering if this is currently usable in any way? I'm looking for specifically this feature, tried pulling the branch but errors for the example case in the initial comment as

ERROR: DimensionMismatch: matrix is not square: dimensions are (10, 9)

[shashi]

Hi @Suavesito-Olimpiada I'm sorry we didn't finish this branch up. I think this looks good to me now. Maybe could use a bit more test coverage.

cc @ChrisRackauckas

[dom-linkevicius]

Hi,

I tried again using this PR, found my mistake in trying to install it the previous time, but the PR itself seems to throw incorrectly for over-determined systems. For example,

using Symbolics

vars = @variables a b c

eqs = [a + 2*b ~ 0,
c + 2*a ~ 0,
a+b+c-3 ~ 0]

Symbolics.solve_for(eqs, vars)

correctly returns

3-element Vector{Float64}:
 -2.0
  1.0
  4.0

but then if I do

eqs = [a + 2*b ~ 0,
c + 2*a ~ 0,
4*b - c ~ 0,
a+b+c-3 ~ 0]

Symbolics.solve_for(eqs, vars)

which ought to have the same solution instead throws

ERROR: ArgumentError: Inconsistent linear system
Stacktrace:
 [1] _solve(A::Matrix{Num}, b::Vector{Num}, vars::Vector{Num}, do_simplify::Bool)
   @ Symbolics C:\Users\domli\.julia\packages\Symbolics\ThKGL\src\linear_algebra.jl:305
 [2] solve_for(eq::Vector{Equation}, var::Vector{Num}; simplify::Bool, check::Bool)
   @ Symbolics C:\Users\domli\.julia\packages\Symbolics\ThKGL\src\linear_algebra.jl:276
 [3] solve_for(eq::Vector{Equation}, var::Vector{Num})
   @ Symbolics C:\Users\domli\.julia\packages\Symbolics\ThKGL\src\linear_algebra.jl:268
 [4] top-level scope
   @ REPL[33]:1

@Suavesito-Olimpiada

[dom-linkevicius]

There also seem to be some numerical issues. Taking as an example:

var = @variables a, b, c
eqn = [a + b + c ~ 2,
       6a - 4b + 5c ~ 31,
       5a + 2b + 2c ~ 13]
Symbolics.solve_for(eqn, var)

gives

3-element Vector{Float64}:
  3.0
 -2.0
  1.0

whereas

var = @variables a, b, c
eqn = [a + b + c ~ 2,
       6a - 4b + 5c ~ 31,
       5a + 2b + 2c ~ 13,
       a ~ 3c]
Symbolics.solve_for(eqn, var)

throws

ERROR: ArgumentError: Inconsistent linear system
Stacktrace:
 [1] _solve(A::Matrix{Num}, b::Vector{Num}, vars::Vector{Num}, do_simplify::Bool)
   @ Symbolics C:\Users\domli\.julia\dev\Symbolics.jl\src\linear_algebra.jl:315
 [2] solve_for(eq::Vector{Equation}, var::Vector{Num}; simplify::Bool, check::Bool)
   @ Symbolics C:\Users\domli\.julia\dev\Symbolics.jl\src\linear_algebra.jl:284
 [3] solve_for(eq::Vector{Equation}, var::Vector{Num})
   @ Symbolics C:\Users\domli\.julia\dev\Symbolics.jl\src\linear_algebra.jl:276
 [4] top-level scope
   @ REPL[64]:1

inspecting the b array returned from _factorize in _solve it is

b: Num[3.0, -2.0, 0.9999999999999997, 1.3322676295501878e-15]

I think there's one bug in addition to this which I will try to find a test case for, but linear algebra is not my strongest suit.