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?
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:
(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.
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 .
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:
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
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.
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.
