Checking invertibility of matrices and solvability of linear equations

[hwjsnc]

The functions to invert a matrix or solve a system of linear equations such as solve, solve_left, or inv currently throw an unspecific ErrorException. For example:

julia> M = matrix(ZZ, 2, 2, [1, 1, 1, 1])
[1  1]
[1  1]

julia> inv(M)
ERROR: Singular matrix in solve_fflu
Stacktrace:
 [1] solve_fflu at /home/user/.julia/packages/AbstractAlgebra/7m2YN/src/generic/Matrix.jl:1902 [inlined]
 [2] pseudo_inv at /home/user/.julia/packages/AbstractAlgebra/7m2YN/src/generic/Matrix.jl:2390 [inlined]
 [3] inv(::AbstractAlgebra.Generic.MatSpaceElem{BigInt}) at /home/user/.julia/packages/AbstractAlgebra/7m2YN/src/generic/Matrix.jl:2409
 [4] top-level scope at REPL[3]:1
 [5] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.0/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288

julia> x = matrix(ZZ, 1, 2, [1, 0])
[1  0]

julia> solve_left(M, x)
ERROR: Unable to solve linear system
Stacktrace:
 [1] solve_left(::AbstractAlgebra.Generic.MatSpaceElem{BigInt}, ::AbstractAlgebra.Generic.MatSpaceElem{BigInt}) at /home/user/.julia/packages/AbstractAlgebra/7m2YN/src/generic/Matrix.jl:2216
 [2] top-level scope at REPL[5]:1
 [3] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.0/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288

There doesn't appear to be a convenient way to check beforehand if a solution is exists.

In a user program, it is hard to catch the exception and figure out the source of the problem: Is it an unsolvable problem, a bug, a memory exception, ...?

For this reason, I believe it would be useful to have functions that can be called to tell if the problem is solvable (for example isinvertible or hasinv); and/or functions that return a union type (for example try_solve).

[thofma]

Yes, the linear algebra solving part needs some love. There will be a can_solve (something in that direction), which gives you both the bool and a solution. We have this already in some other parts of the eco system (e.g. Nemo), but not here yet.

[wbhart]

Agreed. I think isinvertible would be the canonical name. And as @thofma has suggested, can_solve seems to be the name of choice elsewhere for solve.

If you want to implement something like this, feel free. Extra points if it is relatively fast with early termination in the case that it is not invertible/soluble.

[thofma]

The methods are in principle at https://github.com/thofma/Hecke.jl/blob/master/src/Misc/Matrix.jl#L1759-L1787, but as you/we noticed they are probably transposing too much :)

It is a bit messy to get all the combinations of left/right and w/o kernel right.

In view of Nemocas/Nemo.jl#862, it should probably be can_solve_with_solution and can_solve_with_solution_and_kernel and the left/right a keyword argument.

[hwjsnc]

Thanks for the quick replies, I will try to implement this in the next few days. (No promises though, because I don't know how much effort it may be.)