LinearAlgebra: Add bareiss det for BigInt Matrices (#40128) :: Take 2 #40868
Conversation
stdlib/LinearAlgebra/src/generic.jl
Outdated
| # Resolve Issue #40128 | ||
| function det(A::AbstractMatrix{BigInt}) | ||
| m, n = size(A) | ||
| if m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))")) |
There was a problem hiding this comment.
Seems like this test should be in det_bareiss!, and then just have:
det(A::AbstractMatrix{BigInt}) = det_bareiss(A)Note also that you can simply call checksquare(A) in det_bareiss!
There was a problem hiding this comment.
Is there any reasoning behind this design choice? I'm trying to learn as I go along, so It would be helpful if you could shed some light on this.
I'll add a commit correcting this, thanks!
There was a problem hiding this comment.
It shouldn't matter whether someone calls det_bareiss! directly, or through det_bareiss, or through det — in all cases, it should check whether the matrix is square, because otherwise that function doesn't make sense.
stdlib/LinearAlgebra/test/generic.jl
Outdated
| @@ -355,6 +355,10 @@ end | |||
| @test [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] ≈ [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] | |||
| end | |||
|
|
|||
| @testset "Issue 40128" begin | |||
| @test LinearAlgebra.det_bareiss(BigInt[1 2; 3 4]) == -2 | |||
There was a problem hiding this comment.
Fix the indenting here. (4 spaces)
There was a problem hiding this comment.
Also check that the result is a BigInt
There was a problem hiding this comment.
Note that det([1 2; 3 4]) == -2 returns true as well, so this is not a great test.
There was a problem hiding this comment.
Specifically, you should probably test something like BigInt[1 big(2)^65+1; 3 4]) which will ensure that no Float based method would produce the right answer.
There was a problem hiding this comment.
You could try det(BigInt[9 1 8 0; 0 0 8 7; 7 6 8 3; 2 9 7 7])::BigInt == -1, since for Int this is off by about 1e-12 (and even for BigInt it currently gives a non-integer result).
There was a problem hiding this comment.
Specifically, you should probably test something like BigInt[1 big(2)^65+1; 3 4]) which will ensure that no Float based method would produce the right answer.
Actually the current LU algorithm gives the exact answer here:
julia> det(BigInt[1 big(2)^65+1; 3 4]) - (4 - 3*(big(2)^65+1))
0.0(Remember that it uses BigFloat.)
| return sign * M[end,end] | ||
| end | ||
| """ | ||
| LinearAlgebra.det_bareiss(M) |
There was a problem hiding this comment.
Don't put the LinearAlgebra. in the docstring, I think.
There was a problem hiding this comment.
I'll change it, butLinearAlgebra.jl contains docstrings like
This is the one example I could find though, should I change this too?
There was a problem hiding this comment.
I would leave the peakflops case out of this PR.
|
@stevengj thanks for the test case. Turns out the Thanks! |
|
Needs a NEWS.md entry, otherwise looking good |
We have several implementations of the Bareiss fraction-free row-reduction algorithm in the Julia ecosystem. The one in Base was added in #40868 to compute exact determinants. We also have implementations in MTK [1] and Modia [2]. The MTK and Modia versions additionally have support for custom pivot selection, open-code a sparse matrix data structure adapted to their domains and support rank-deficient matrices. I would like to separate out the algorithmic and data-structures concerns so that they may tested independently. Of course this function isn't particularly large, but implementing it correctly and performantly is still not trivial, so it seems better to have one implementation rather than three. [1] https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/alias_elimination.jl#L236 [2] https://github.com/ModiaSim/ModiaBase.jl/blob/6c341eed72d9867553cb9565330d8ae85221b343/src/LinearIntegerEquations.jl#L204
We have several implementations of the Bareiss fraction-free row-reduction algorithm in the Julia ecosystem. The one in Base was added in #40868 to compute exact determinants. We also have implementations in MTK [1] and Modia [2]. The MTK and Modia versions additionally have support for custom pivot selection, open-code a sparse matrix data structure adapted to their domains and support rank-deficient matrices. I would like to separate out the algorithmic and data-structures concerns so that they may tested independently. Of course this function isn't particularly large, but implementing it correctly and performantly is still not trivial, so it seems better to have one implementation rather than three. [1] https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/alias_elimination.jl#L236 [2] https://github.com/ModiaSim/ModiaBase.jl/blob/6c341eed72d9867553cb9565330d8ae85221b343/src/LinearIntegerEquations.jl#L204
We have several implementations of the Bareiss fraction-free row-reduction algorithm in the Julia ecosystem. The one in Base was added in #40868 to compute exact determinants. We also have implementations in MTK [1] and Modia [2]. The MTK and Modia versions additionally have support for custom pivot selection, open-code a sparse matrix data structure adapted to their domains and support rank-deficient matrices. I would like to separate out the algorithmic and data-structures concerns so that they may tested independently. Of course this function isn't particularly large, but implementing it correctly and performantly is still not trivial, so it seems better to have one implementation rather than three. [1] https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/alias_elimination.jl#L236 [2] https://github.com/ModiaSim/ModiaBase.jl/blob/6c341eed72d9867553cb9565330d8ae85221b343/src/LinearIntegerEquations.jl#L204
New users are commonly surprised that determinants of integer matrices give an approximate floating-point answer (#40128), and are unaware that an exact algorithm is implemented for `BigInt` matrices (#40868). This PR comments on both of these facts in the `det` documentation. (At some point, we may want to mark `LinearAlgebra.det_bareiss` as `public`, and document it, but that can be done in a future PR.)
New users are commonly surprised that determinants of integer matrices give an approximate floating-point answer (JuliaLang#40128), and are unaware that an exact algorithm is implemented for `BigInt` matrices (JuliaLang#40868). This PR comments on both of these facts in the `det` documentation. (At some point, we may want to mark `LinearAlgebra.det_bareiss` as `public`, and document it, but that can be done in a future PR.)
As said in the previous PR with the same name, I messed up my julia fork trying to rebase and push clean code.
This fixes #40128
Do suggest ways to improve this PR.
cc: @simeonschaub