# LinearAlgebra.det and symbolic variables

I am interested in the computation of the determinant of a (very sparse) matrix involving a symbolic variable.

My code is

``````using Symbolics, LinearAlgebra, SparseArrays

N = 11

FloatType = Float64

c = zeros(FloatType, N)
for k in 2:N
c[k] = (k - 1)/(2.0*(N - 1))
end

A = zeros(FloatType, N, N)
A[2, 1] = c
for k in 3:N
A[k, 1]   = c[k] - 1
A[k, k-1] = 1
end

A = sparse(A)
display(A); println()

b = zeros(FloatType, N)
b[N] = 1.0
b = sparse(b)

@variables z
display(I - z * A + z * ones(FloatType, N) * b'); println()
println("Start computing determinant")
Determinant = det(I - z * A + z * ones(FloatType, N) * b')
println("Start expanding")
Determinant = expand(Determinant)

println(Determinant)

``````

So I use the `det` from Julia Linear Algebra and the symbolic variable from `Symbolics.jl`.

For small `N` the program runs somewhat quick, but already for the given `N=11` it takes literally ages. Is there any hope to speed this up?

I can run the same problem for instance in Mathematica in basically no time, even for `N=100` it returns instantly.

Here’s a solution using Nemo.jl, which finishes under 3 seconds for N=100. Nemo has its own matrix type, which is a little different from Julia’s own `Matrix`. I also use Nemo’s rational number type `QQ` instead of `Float64`. Unfortunately, Nemo does not support sparse matrices, so I’m computing with a dense matrix, which will cause memory usage problems if `N` is really large.

``````using Nemo, LinearAlgebra

N = 100

r, z = PolynomialRing(QQ, "z")

c = fill(r(0), N) # `r(0)` is the zero polynomial. `c` is a vector of polynomials
for k in 2:N
c[k] = r((k - 1)//(2*(N - 1)))
end

A = MatrixSpace(r, N, N)() # This is using Nemo's matrix type. Intialized to zero matrix
A[2, 1] = c
for k in 3:N
A[k, 1]   = c[k] - r(1)
A[k, k-1] = r(1)
end

b = MatrixSpace(r, 1, N)()
b[1, N] = r(1)

matrix = MatrixSpace(r, N, N)(1) - z * A + z * MatrixSpace(r, N, 1)(r(1)) * b

display(Matrix(matrix)) # convert back to Julia's native Matrix type before calling `display`, as Nemo's printing of matrices is not as pretty

@time Determinant = det(matrix) # Nemo expands polynomials by default. No need to expand further

``````
2 Likes

This doesn’t have a specialization. It’s just using a basic fallback. It’s worth an issue to document that a specialization should be made.

Thanks, do you also know a function how to convert the coefficients of the determinant of type `fmpq` back to something like `Float64` ? I couldn’t find something immediately in the docs.

Nemo provides conversion between `fmpq` and Julia’s `Rational{BigInt}`, so you can manually write a conversion function to `Float64`:

``````function Base.convert(::Type{Float64}, a::Nemo.fmpq)
convert(Float64, convert(Rational{BigInt}, a))
end
``````

Usage example:

``````julia> convert(Float64, fmpq(1//2))
0.5
``````
1 Like

I found a small bug. `MatrixSpace(r, N, 1)(r(1))` gives not a whole vector of ones, but initializes only the first entry to one. My solution to this is

``````ones = MatrixSpace(r, N, 1)()
for k = 1:N
ones[k, 1] = r(1)
end

matrix = MatrixSpace(r, N, N)(1) - z * AMat1 + z * ones * b
``````

You’re right. The constructor only fills the diagonal elements, which is just one element in this special case. You could also first construct a Julia matrix and then convert it.

``````ones = MatrixSpace(r, N, 1)(fill(r(1), N, 1))
``````
1 Like