JuliaAlgebra / MultivariatePolynomials.jl

Multivariate polynomials interface
https://juliaalgebra.github.io/MultivariatePolynomials.jl/stable/
Other
134 stars 27 forks source link

use the LinearAlgebraX package for linear algebra? #281

Open nsajko opened 8 months ago

nsajko commented 8 months ago

The LinearAlgebra.det method defined in this package is not very good:

  1. It's restricted to Matrix, so LinearAlgebra.det fails for, e.g., LinearAlgebra.Symmetric polynomial matrices
  2. It doesn't check that the matrix is square
  3. It uses recursion with depth that's not statically known
  4. It uses sum without an init argument, causing a dynamic dispatch warning with JET

The LinearAlgebraX package, however, provides two different methods that both work with polynomial matrices (as far as I tried): detx and cofactor_det.

Perhaps LinearAlgebraX could be used to improve the linear algebra functionality in this package.

nsajko commented 8 months ago

Slight correction: only cofactor_det works for polynomial matrices, while detx uses it as a fallback in this case, after a try-catch. So I think that LinearAlgebra.det should be implemented using LinearAlgebraX.cofactor_det here.

nsajko commented 8 months ago

Even for Matrix, which is supported by both det and cofactor_det, the latter is much faster. I also tried improving the current code for det, but it was still much slower than cofactor_det for the matrix I tried:

function det_mapreduce_map(M)
    f = i -> Int8(iseven(i) ? -1 : 1)
    let M = M, f = f
        i -> f(i) * M[i, 1] * LinearAlgebra.det(M[1:end.!=i, 2:end])
    end
end

function det_impl(M::Matrix{<:AbstractPolynomialLike})
    m = size(M)[1]
    if m > 2
        mapreduce(
            det_mapreduce_map(M), +, Base.OneTo(m), init = zero(eltype(M)),
        )
    elseif m == 2
        return M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2]
    else
        return M[1, 1]
    end
end

function LinearAlgebra.det(M::Matrix{<:AbstractPolynomialLike})
    LinearAlgebra.checksquare(M)
    det_impl(M)
end

I guess there are two options now:

  1. Depend on LinearAlgebraX
  2. Copy the (short) cofactor_det implementation into this package

What do you think?

blegat commented 8 months ago

I'm a bit hesitant to add this dependency but I would accept a PR improving the current implementation