add support for solving arbitrary linear systems #757
Conversation
Codecov Report
📣 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
📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
|
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
changed the title
|
Bump |
|
@Suavesito-Olimpiada did you ever try to profile it after the Unityper update? I'm running CI and will merge after that. |
|
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
force-pushed
the
suave/linalg
branch
from
046cbc1 to
21673e5
Compare
|
I haven't tried to profile it again, but I'm doing it. I'll add more test. |
|
I have profiled it again, here are the results @shashi.
Where I installed Symbolics.jl doing |
|
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 |
|
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. |
|
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.0but 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 |
| # swap rows, only needed to swap lead:end | ||
| if k != kp | ||
| for j = lead:n | ||
| F[k, j], F[kp, j] = F[kp, j], F[k, j] |
There was a problem hiding this comment.
Rows of b are not being swapped, which results in incorrect throwing of inconsistent system for consistent over-determined systems. Adding
b[k], b[kp] = b[kp, b[k]fixes this
|
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.0whereas 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]:1inspecting the 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. |
| for i = 1:m | ||
| if i != k | ||
| # this line occupies most of the time, distributed in the | ||
| # following methods |
There was a problem hiding this comment.
This isn't quite correct as well, because it doesn't set the columns of the lead index (other than the row of the pivot) to zero (it does so below, but don't see why splitting is necessary/better), as well as can result in numerical issues due to Floating points. Something like
for i = 1:m
if i != k
b[i] = b[i] - F[i, lead] * b[k]
F[i, :] .= F[i, :] - F[i, lead] * F[k, :]
F[abs.(F) .< atol] .= 0
b[abs.(b) .< atol] .= 0
end
endwhere atol is some reasonable cut off value, e.g. 1e-10. Not the cleanest solution, but does the trick on the example cases in the comments that it failed previously on.
There was a problem hiding this comment.
One might also need to add something like
b[i] = simplify.(b[i] - F[i, lead] * b[k])
F[i, :] .= simplify.(F[i, :] - F[i, lead] * F[k, :])because for larger/more complex fully symbolic systems the unsimplified expressions are different enough so that they don't get set to 0 at appropriate iterations leading to incorrect solutions.
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
lufactorization for square rank-complete systems; for under-determined systems it useslufactorization first, with arrefbased solver later, and finallyrrefdirectly for over-determined systems.Usage example
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$91 \times 190$ .
ProfileViewofSymbolics.solve_forlooks like for a system ofYou can notice that 8-11 and 14-17 are
*(::Num, ::Num)or-(::Num, ::Num), marked red for being dynamic dispatch.