CODEWITHSUNDEEP
juliadifferential-equationspolynomial-math
Yifei
Yifei

Reputation: 103

How to solve polynomial equations fast in Julia?

I'm working on evolution simulations in biological systems. I have to solve polynomial equations, finding the root(u*X^3 - N*p*r*X^2 - N*p*X^2 + K^2*u*X - N*K^2*p), where u and K are constants, N is an constant array, and p, r are evolving parameters. Basically, for each individual in the population at each generation, I need to do the below calculations (length(N)>>length(p)):

for i = 1:length(p)
    for j = 1:length(N)
         X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
    end
end

I'm aware that I can parallel the code by solving the equations for different individuals using different cores, and even within each individual I can parallel solve each X[j,i]. I'm wondering what is the best practice/fast way to deal with this situation? Is it possible to solve a single equation in a much faster way?

Upvotes: 3

Views: 1018

Answers (1)

Dan Getz
Dan Getz

Reputation: 18227

This answer compares Roots, GSL and PolynomialRoots packages for solving a version of the problem in the question. The timings can vary on machines and it is best to run it locally. But essentially, GSL is fastest and about 1000x faster than Roots, with PolynomialRoots somewhere in between (but closer to GSL). See comments for more information.

The code:

# using Julia 0.5.1-pre+4
# Pkg.add("Roots")
# Pkg.add("GSL")
# Pkg.add("PolynomialRoots")

using GSL
using Roots
using PolynomialRoots

function test1(p,r,N,u,K)
  X = Matrix{Float64}(length(N),length(p))
  for i = 1:length(p)
    for j = 1:length(N)
      X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
    end
  end
  return X
end

function test2(p,r,N,u,K)
  X = Matrix{Float64}(length(N),length(p))
  uinv = inv(u)
  for i = 1:length(p)
    for j = 1:length(N)
      X[j,i] = mean(poly_solve_cubic(-uinv*(p[i]*N[j]*r[i]+p[i]*N[j]),K^2,-uinv*p[i]*K^2*N[j]) )
    end
  end
  return X
end

function test3(p,r,N,u,K)
  X = Matrix{Float64}(length(N),length(p))
  for i = 1:length(p)
    for j = 1:length(N)
      X[j,i] = mean(filter(x->abs(imag(x))<1.0e-10,
        PolynomialRoots.solve_cubic_eq(Complex{Float64}[- p[i]*K^2*N[j], K^2*u, - p[i]*N[j]*r[i] - p[i]*N[j],u])))
    end
  end
  return X
end

K = 1.0
u = 1.0
N = rand(1000)+1.0
p = rand(10)+1.0
r = rand(10)+1.0

res1 = test1(p,r,N,u,K);
res2 = test2(p,r,N,u,K);
res3 = test3(p,r,N,u,K);

using Base.Test
@test_approx_eq res1 res2
@test_approx_eq res1 res3

@time test1(p,r,N,u,K); # Roots
@time test2(p,r,N,u,K); # GSL
@time test3(p,r,N,u,K); # PolynomialRoots

nothing

My timings:

 20.664363 seconds (225.67 M allocations: 13.650 GB, 18.81% gc time) # Roots
  0.010303 seconds (80.01 k allocations: 4.044 MB, 75.30% gc time) # GSL
  0.215453 seconds (394.90 k allocations: 9.917 MB) # PolynomialRoots

It would be nice to hear a report on the real problem timings if you compare the methods. Thanks.

Upvotes: 6

Related Questions