LinearAlgebra.det and symbolic variables

General Usage Performance
, ,
1

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[2]
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.

2

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[2]
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
3

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.

4

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.

5

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
6

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
7

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