Better determinant? det looks horrible

using SymEngine
using LinearAlgebra

[a*z2  2*z2 + a*z1    z4    z3;
2*z1*z3    -6*b*z2^2  z1^2     0;
2*z1 - z2          -z1  2*z3     0;
    -z3         2*z2   -z1  2*z4]|>det

The output is horrible

(-(-(2*a*z2^2 + z3*(2*z2 + a*z1))*(a*z2*z1^2 - 2*z1*z4*z3) + (z4*z3 - a*z2*z1)*(-6*a*b*z2^3 - 2*z1*z3*(2*z2 + a*z1)))*(2*z1*z3^2*(-(2*z1 - z2)*(2*z2 + a*z1) 
- a*z2*z1) - (2*z1 - z2)*z3*(-6*a*b*z2^3 - 2*z1*z3*(2*z2 + a*z1)))/(a^2*z2^2)
 + ((2*a*z2*z4 + z3^2)*(-6*a*b*z2^3 - 2*z1*z3*(2*z2 + a*z1)) 
+ 2*z1*z3^2*(2*a*z2^2 + z3*(2*z2 + a*z1)))*((-(2*z1 - z2)*z4 + 2*a*z2*z3)*(-6*a*b*z2^3 - 2*z1*z3*(2*z2 + a*z1)) 
- (a*z2*z1^2 - 2*z1*z4*z3)*(-(2*z1 - z2)*(2*z2 + a*z1) - a*z2*z1))/(a^2*z2^2))/(-6*a*b*z2^3 - 2*z1*z3*(2*z2 + a*z1))

I am pretty sure the determinate should have more compact form. Why the det in LinearAlgebra generator so horrible form? Is there any other function can generate simple form?

What happens if you simplify it?

It will simplify to

-z1^3 (z3^2 + 4 z2 (z3 - 2 z4)) + 4 a z1^4 z4 - 
 12 b z2^2 (z3^3 + 2 a z2 z3 z4 + z2 z4^2) + 
 2 z1^2 (z2^2 (z3 + 6 b z3 - 2 z4) - 2 z3 z4 (2 a z3 + z4)) + 
 2 z1 z2 (4 z3^2 (z3 - 2 z4) - 3 b z2 (z2 z3 - 4 z4^2))

by Mathematica.

The denominator from det in LinearAlgebra looks wired to me.

Probably because LinearAlgebra.det is not using a division-free algorithm (since division-free determinant algorithms don’t scale), whereas you want a final expression with no divisions, and SymEngine is not able to simplify well enough to realize that the divisions exactly cancel (you might need some option to tell it not to worry about dividing by zero?).

Symbolics.jl does a better job, however, and is able to eliminate all the divisions:

julia> @variables a b z1 z2 z3 z4;

julia> [a*z2  2*z2 + a*z1    z4    z3;
              2*z1*z3    -6*b*z2^2  z1^2     0;
              2*z1 - z2          -z1  2*z3     0;
                  -z3         2*z2   -z1  2*z4] |> det
z4*(12b*z4*(z2^2)*(2z1 - z2) - 4z3*z4*(z1^2)) + (-2z2 - a*z1)*(8z1*z4*(z3^2) - 2z4*(z1^2)*(2z1 - z2)) + a*z2*(2z4*(z1^3) - 24b*z3*z4*(z2^2)) - z3*((z1^2)*(2z2*(2z1 - z2) - z1*z3) + 6b*(z2^2)*(2(z3^2) - z1*(2z1 - z2)) + 2z1*z3*(z1^2 - 4z2*z3))

(Division-free determinant algorithms are something it might be nice to have in a package, but they are fairly specialized in their practical utility.)

Realize that the number of terms in this kind of expanded division-free form will generally grow faster than exponentially with the matrix size.

2 Likes

Seems det in LinearAlgebra is only suitable for numerical problems. For Symbolic, it is not quite match now.

Matlab has some options in det function. Maybe the developers can extend it in LinearAlgebra similarly.
B = det(A,‘Algorithm’,‘minor-expansion’) uses the minor expansion algorithm to evaluate the determinant of A .

It seems to me the Bareiss determinant has been added to Julia last year

Bareiss is exact for integer matrices, but is not division-free.

The minor-expansion algorithm has \Theta(n!) complexity, no? So it’s not really practical except for pedagogy (or tiny matrices).

There are O(n^4) (or slightly better) algorithms for division-free determinants, such as this one by Bird (2011) or this one by Kaltoven (1994). The Bird algorithm in particular is really beautifully simple:

6 Likes

Bareiss will work for any ring where exactdiv is implemented.

Yes, but that’s not the same thing as division-free.

For example, here is an implementation of the Bird algorithm, which is truly division-free (only uses ± and ×), at the cost of worse complexity than Bareiss:

using LinearAlgebra
function birddet(A::AbstractMatrix)
    n = LinearAlgebra.checksquare(A)
    
    function μ!(M, X)
        for j = 2:n, i = 1:j-1
            M[i,j] = X[i,j]
        end
        for j = 1:n-1, i = j+1:n
            M[i,j] = zero(eltype(X))
        end
        M[n,n] = zero(eltype(X))
        for k = n-1:-1:1
            M[k,k] = M[k+1,k+1] - X[k+1,k+1]
        end
        return M
    end
    
    M = similar(A)
    X = copy(A)
    for i = 1:n-1
        mul!(X, μ!(M, X), A)
    end
    return isodd(n) ? X[1,1] : -X[1,1]
end
5 Likes

Thanks for the birddet and Bareiss determinant. Now I understand, the LinearAlgebra does have different determinant. It just renamed, it does not use det with options.

Another question: is there any different inverse function suitable for large symbolic matrix?
I tried inv(AA), it is extremely slow.

AA=Basic[0 0 0 0 0 0 0 0 0 0 0 0 4*a 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 4*a 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*a 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4*a 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2*a 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2*a 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2*a 4 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4*a 4 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4*a 12 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2*a 4 0 0 0 0 0 0 0 0 12 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*a 12 0 0 0 0 0 24 0; 0 0 0 0 -48*b 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 -12*b 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 -48*b 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -24*b 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -12*b 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -48*b 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -12*b 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -24*b 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 120 -24 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 24 -12 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 24 -6 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 6 -4 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 12 -4 0 0 0 0 0 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 4 -4 0 0 0 0 0 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 -6 0 0 0 120 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 24 -6 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 -4 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 -2 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 -2 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -2 0 0 0 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 -4 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -4 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -2 0 0 0 12 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 -6 0 0 0 12 0 0 0 0; 0 0 12 0 0 0 -24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 12 0 0 0 -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 24 0 0 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 120 0 0 0 -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 4 0 0 -12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 6 0 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 24 0 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 -12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 -24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 24 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 -4 0 0 0 0 0 0 0 0 0 0 0 24 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 -2 0 0 0 0 0 0 0 0 0 0 0 24 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -4 0 0 0 0 0 0 0 0 0 0 24 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 -6 0 0 0 0 0 120; 0 0 0 0 0 0 0 -24 8 + 48*b -36*b 0 -12 0 0 0 0 48 -144*b 0 0 0 96*a 48 -16 0 0 0 0 0 -144*a*b -32*a -32 0 0 0 0 0 0 96*b -144*b -16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

There is the comatrix algorithm (see Cohen, “Computational algebraic number theory 2.2.7”) which
computes simultaneously the comatrix (the matrix of minors) and the coefficients of the characteristic
polynomial. It is reasonably fast but again is not completely division free (it uses exact division which is easy to implement for multivariate polynomials).

function charpolyandcomatrix(m)
  C=one(m)
  res=C
  n=size(m,1)
  a=ones(eltype(m),n+1)
  for i in 1:n
    if i==n res=(-1)^(n-1)*C end
    C*=m
    a[n+1-i]=exactdiv(-sum(C[i,i] for i in axes(C,1)),i)
    if i!=n C+=a[n+1-i]*one(C) end
  end
  a,res
end

In the result a is the vector of ceofficients of the characteristic polynomial and res the comatrix.

1 Like