From 2889bfb2f6bf041d29a53c7ddae6a016f63dfd10 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Fri, 24 Nov 2023 12:41:51 +0100
Subject: [PATCH] `LinearAlgebra.det` improvements
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Performance improvements, support for more types.

Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```julia-repl
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 11 on 8 virtual cores
Environment:
  JULIA_NUM_PRECOMPILE_TASKS = 3
  JULIA_PKG_PRECOMPILE_AUTO = 0

julia> using LinearAlgebra, DynamicPolynomials

julia> function f(n)
         @polyvar a b c d e
         diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
         )
       end
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  2.489404 seconds (46.04 M allocations: 2.244 GiB, 18.91% gc time, 13.04% compilation time)

julia> @time det(m15);
  2.231880 seconds (45.94 M allocations: 2.238 GiB, 21.50% gc time)

julia> @time det(m16);
  5.362580 seconds (107.70 M allocations: 5.243 GiB, 23.50% gc time)

julia> @time det(m16);
  5.405048 seconds (107.70 M allocations: 5.243 GiB, 23.65% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes #281
---
 src/det.jl  | 80 ++++++++++++++++++++++++++++++++++++++++++++++-------
 test/det.jl | 12 ++++++--
 2 files changed, 79 insertions(+), 13 deletions(-)

diff --git a/src/det.jl b/src/det.jl
index 66cd8dce..91f92faf 100644
--- a/src/det.jl
+++ b/src/det.jl
@@ -1,13 +1,73 @@
-function LinearAlgebra.det(M::Matrix{<:AbstractPolynomialLike})
-    m = size(M)[1]
-    if m > 2
-        return sum(
-            (-1)^(i - 1) * M[i, 1] * LinearAlgebra.det(M[1:end.!=i, 2:end]) for
-            i in 1:m
-        )
-    elseif m == 2
-        return M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2]
+# Scalar determinant, used for recursive computation of the determinant
+LinearAlgebra.det(p::AbstractPolynomialLike{<:Number}) = p
+
+# Matrix determinant by cofactor expansion, adapted from
+# `LinearAlgebraX.cofactor_det`.
+function det_impl(A::AbstractMatrix{T}) where {T}
+    r = first(size(A))
+    if 1 < r
+        total = LinearAlgebra.det(zero(T))
+        for i ∈ Base.OneTo(r)
+            a = LinearAlgebra.det(A[i, 1])
+            if !iszero(a)
+                ii = Base.OneTo(r) .!= i
+                jj = 2:r
+                B = A[ii, jj]
+                x = det_impl(B)
+                # XXX: should be faster but doesn't work:
+                #  a = MA.operate!!(*, a, x)
+                a *= x
+                if iseven(i)
+                    # XXX: not implemented yet:
+                    #  a = MA.operate!!(-, a)
+                    a = -a
+                end
+                total = MA.operate!!(+, total, a)
+            end
+        end
+        total
+    elseif isone(r)
+        LinearAlgebra.det(A[1, 1])
     else
-        return M[1, 1]
+        error("unexpected")
     end
 end
+
+collect_if_not_already_matrix(m::Matrix) = m
+collect_if_not_already_matrix(m::AbstractMatrix) = collect(m)
+
+function det_impl_outer(m::AbstractMatrix{T}) where {T}
+    if 0 < LinearAlgebra.checksquare(m)
+        det_impl(collect_if_not_already_matrix(m))
+    else
+        LinearAlgebra.det(one(T))
+    end
+end
+
+# Determinants of narrow integer type: `LinearAlgebra` seems to
+# promote these to `Float64` to prevent them from overflowing. We
+# instead promote to `BigInt` to keep things exact. In the case of
+# `Bool` we also need to promote for type stability.
+
+const NarrowIntegerTypes = Union{
+    Bool, UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64,
+    UInt128, Int128,
+}
+
+const NarrowIntegerPolynomialLike =
+  AbstractPolynomialLike{T} where {T<:NarrowIntegerTypes}
+
+promote_if_narrow(m::AbstractMatrix{<:AbstractPolynomialLike}) = m
+
+promote_if_narrow(m::AbstractMatrix{<:NarrowIntegerPolynomialLike}) =
+  map((p -> polynomial(p, BigInt)), m)
+
+# For type stability, we want to promote termlikes to polynomiallikes
+# before attempting to calculate the determinant.
+promote_if_termlike(m::AbstractMatrix{<:AbstractPolynomialLike}) = m
+promote_if_termlike(m::AbstractMatrix{<:AbstractTermLike}) = map(polynomial, m)
+
+promote_if_necessary(m) = promote_if_termlike(promote_if_narrow(m))
+
+LinearAlgebra.det(m::AbstractMatrix{<:AbstractPolynomialLike}) =
+    det_impl_outer(promote_if_necessary(m))
diff --git a/test/det.jl b/test/det.jl
index dad8f2fc..b946899a 100644
--- a/test/det.jl
+++ b/test/det.jl
@@ -1,7 +1,13 @@
 @testset "Det" begin
     Mod.@polyvar x y
 
-    @test det([x 1; y 1]) == x - y
-    @test det([x 1; x 1]) == 0
-    @test det([x 1 1 1; x y 1 2; 0 0 0 1; x 0 y 0]) == -x * y^2 + 2 * x * y - x
+    @testset "zero-one $T" for T ∈ (Bool, Int, Float64, BigInt, BigFloat)
+        z = zero(T)
+        o = one(T)
+        @test (@inferred det([x o; y o])) == x - y
+        @test (@inferred det([x o; x o])) == 0
+        @test (@inferred det([x+y y; z y])) == (x+y)*y
+    end
+
+    @test (@inferred det([x 1 1 1; x y 1 2; 0 0 0 1; x 0 y 0])) == -x * y^2 + 2 * x * y - x
 end