Finding roots of polynomials in Haskell?

Posted on February 13, 2019 by Brent

tl;dr: There are several Haskell packages one can use to find an individual root of a function on a certain interval. But I’ve had less luck finding a suitable package for finding all the roots of a given polynomial. This blog post consists of my notes and a plea for help.

The diagrams-solve package contains miscellaneous numeric solving routines used in diagrams, including tridiagonal and cyclic tridiagonal linear system solvers (used for generating cubic splines and emulating Metafont paths) as well as functions for finding roots of low-dimensional polynomials (quadratic, cubic, and quartic). Solving quadratics is used in a bunch of places; solving cubics is needed specifically for doing interior/exterior testing on closed loops built from cubic Bezier curves; thankfully we have never needed to solve a quartic or higher.

Unfortunately, the polynomial root solvers are pretty naive: I simply transcribed some formulas (which is why it only goes up to quartics). This works OK a lot of the time, but the formulas are very numerically unstable, and it’s not hard to come up with example inputs where the returned roots are off by quite a bit. In fact, people regularly run into this when running the test suite. I am not specifically aware of any diagrams bugs that have arisen in actual practice due to the cubic solutions being off, but it’s probably just a matter of time.

So I decided it was finally time to look into better root-finding methods. This blog post is both a plea for help and a place to write down some of the things I’ve learned so far.

There’s root finding, and then there’s root finding

The first thing to get out of the way is that when you talk about “root finding” there are (at least!) two pretty different things you could mean:

  1. Given a function f and some particular interval, or an initial guess, find a value x in the interval/close to the initial guess for which f(x) = 0.
  2. Given a polynomial with {real, complex} coefficients, find all its {real, complex} roots.

If you want to do (1), there are several nice Haskell packages you could use. The Numeric.RootFinding module from the math-functions package is probably your best bet; it implements both Ridders’ method and the Newton-Raphson method, which both attempt to find a single root of a function on a given interval. They both work on any continuous Double -> Double function, not just polynomials (Newton-Raphson also needs to know the first derivative). But they don’t work well if you don’t already have an interval in mind to search; and if you want to find all the roots you have to call these multiple times (somehow coming up with an appropriate interval each time).

As for (2), I haven’t been able to find anything that would work well for diagrams-solve. Here are my notes:

  • The dsp package has a module Polynomial.Roots containing an implementation of Laguerre’s method, which finds all the (complex) roots of a polynomial. However, the dsp package is a rather heavyweight dependency to pull in just for the root-finding code; it’s also licensed under the GPL and I’d like to avoid having to “infect” the entire diagrams ecosystem with the GPL.

  • Laguerre’s method seems like it should be fairly easy to implement myself—but writing my own solver from scratch is how I got here in the first place; I’d really like to avoid it if possible. I am far from being an expert on numerical analysis, floating-point computation, etc.

  • The hmatrix-gsl package has the Numeric.GSL.Polynomials module, which has an interface to a root finding algorithm from GSL (apparently using something called “balanced-QR reduction”), but I’d like to avoid pulling in a C library as as dependency, and also, again, GPL.

  • From the Wikipedia page for Laguerre’s method I learned that the Jenkins-Traub algorithm is another widely used method for polynomial root-finding, and often preferred over Laguerre’s method. However, it seems rather complex, and the only Haskell implementation of Jenkins-Traub I have been able to fnid is this one which seems to be just a toy implementation; I don’t even know if it works correctly.

If you know of a good place where I can find polynomial solving code in Haskell, can you point me to it? Or if you know more about numerics than me, could you maybe whip up a quick implementation of Laguerre’s method and put it up on Hackage?

About Brent

Associate Professor of Computer Science at Hendrix College. Functional programmer, mathematician, teacher, pianist, follower of Jesus.
This entry was posted in haskell, math and tagged , , , , , . Bookmark the permalink.

14 Responses to Finding roots of polynomials in Haskell?

  1. Heinrich Apfelmus says:
    February 14, 2019 at 12:12 pm

    It probably doesn’t apply to your case, but just In case you want to find all roots of an arbitrary function, there is a MATLAB package called `chebfun` that offers interpolation by Chebychev polynomials, which I found highly useful. The roots of Chebychev polynomials can be found by means of a Fourier transform, but that is also hidden away in the package.

    Reply
  2. Jake says:
    February 14, 2019 at 7:18 pm

    I don’t know nearly enough Haskell to do it myself, but https://en.wikipedia.org/wiki/Durand–Kerner_method looks promising for finding the entire root set and easier to implement than Laguerre.

    But honestly what I think I’d do in that situation is bite the bullet about foreign dependencies and bring in hmatrix (NOT hmatrix-gsl). Then to find the roots of a polynomial, build the companion matrix : https://en.wikipedia.org/wiki/Companion_matrix and use hmatrix to compute the eigenvalues, which are exactly the roots.

    Reply
    • Brent says:
      February 14, 2019 at 10:42 pm

      Yeah, I tried implementing Durand-Kerner once but never really got it to work. I might try again.

      Thanks for the advice re: hmatrix. I may end up doing that!

      Reply
    • Alexander Adler says:
      February 17, 2019 at 11:12 am

      I also implemented Durand-Kerner-Weierstraß when I was in highschool. The implementation was in C and had to implement complex arithmetics and polynomials including input and output of them. Altogether about 300 lines. I couldn’t find any polynomial that wouldn’t very rapidly yield all roots (maybe I was just lucky). The code is lost, of course, long ago; it didn’t seem saving-worthy since it was no more than two hours of work. I would give DKW a shoot.

      Reply
  3. Jan Van lent says:
    February 17, 2019 at 4:47 pm

    Converting to an eigenvalue problem is now standard practice. A monomial basis representation leads to companion matrices. If you have an orthogonal basis representation then it is better not convert to monomials and use colleague/comrade matrices (this is the approach used in the Chebfun package already mentioned). Some names to look up: IJ Good, JP Boyd, V Noferini, DS Watkins, R Vandebril, J Aurentz.
    I found an article in the Monad.Reader that discusses a Haskell implementation of the Jacobi method for finding the eigenvalues of a matrix (using the Repa package) [1]. It seems an implementation is available [2].
    The other standard algorithms for finding all eigenvalues are QR and divide-and-conquer, but getting all the details right can be quite involved.

    [1] https://themonadreader.files.wordpress.com/2013/03/issue214.pdf
    [2] https://github.com/felipeZ/Haskell-abinitio/tree/master/Science/QuantumChemistry/NumericalTools

    Reply
  4. Dominic Steinitz (@idontgetoutmuch) says:
    February 19, 2019 at 3:30 am

    Just reading through the comments, you might be able to do some quick experiments using hmatrix. It can find eigenvalues (probably via LAPACK), Chebyshev polynomials (⚠️ – I never tried this) and roots. The latter two are in hmatrix-gsl which is a binding to GSL (GNU Scientific Library) which comes with its GPL!

    Reply
  5. Jan Van lent says:
    February 19, 2019 at 12:56 pm

    The implementation of Laguerre’s method in the DSP package is quite compact [1,2].
    Maybe the original author would be open to having it in a separate package (possibly using a more permissive license).

    [1] https://hackage.haskell.org/package/dsp-0.2.1/docs/src/Polynomial-Roots.html
    [2] http://hackage.haskell.org/package/dsp-0.2.4.1/docs/src/Polynomial.Basic.html

    Reply
  6. Jan Van lent says:
    February 21, 2019 at 9:33 am

    I thought the Alberth method looked quite elegant, so I wrote a basic implementation. This is not meant to be robust and it is not difficult to come up with situations where this basic implementation fails.


    — Alberth method for finding all roots of a (complex) polynomial
    — https://en.wikipedia.org/wiki/Aberth_method
    — The calculations have to use complex numbers.
    — Some optimizations are possible if the coefficients are real,
    — This is not implemented here.
    — Implementation was prompted by
    — https://byorgey.wordpress.com/2019/02/13/finding-roots-of-polynomials-in-haskell/
    — This is not meant to be robust implementation, just a proof of concept.
    — For analysis and implementation of this specific method,
    — see the work by D Bini and others.
    —
    — Jan Van lent
    — Started: 2019-02-20 Wed
    — Last updated: 2019-02-21 Thu
    import Data.Complex
    — |Given [p0, p1, …, pn] x, evaluate polynomial p0 + p1*x + … + pn*x^n
    polyval :: Num a => [a] -> a -> a
    polyval p x = rec p 1 0
    where rec [] z y = y
    rec (p0:p) z y = rec p (x*z) (y + p0*z)
    — alternative: Horner scheme, sparse polynomials, FFT, etc.
    diff :: Num a => [a] -> [a]
    diff p = zipWith (*) (map fromIntegral [1..]) (tail p)
    scales :: Fractional a => [a] -> [a]
    scales zs = [ sum [ 1/(zk-zj) | (zj, j) <- zip zs [1..], j /= k ] |
    (zk, k) <- zip zs [1..] ]
    offset :: Fractional a => a -> a -> a
    offset q s = -q/(1-q*s)
    offsets :: Fractional a => [a] -> [a] -> [a] -> [a]
    offsets fs gs ss = ws
    where
    qs = zipWith (/) fs gs
    ws = zipWith offset qs ss
    solve' :: RealFloat a => [Complex a] -> [Complex a] -> [Complex a] -> Int -> a -> [Complex a]
    solve' p p' zs maxit tol
    | maxit <= 0 = zs
    — | maxf <= tol = zs
    | maxw <= tol = newzs
    | otherwise = solve' p p' newzs (maxit-1) tol
    where
    fs = map (polyval p) zs
    maxf = maximum (map magnitude fs)
    gs = map (polyval p') zs
    ss = scales zs
    ws = offsets fs gs ss
    maxw = maximum (map magnitude ws)
    newzs = zipWith (+) zs ws
    solve :: RealFloat a => [Complex a] -> [Complex a] -> Int -> a -> [Complex a]
    solve p zs maxit tol = solve' p (diff p) zs maxit tol
    rootBound :: RealFloat a => [Complex a] -> a
    rootBound p = r
    where
    revp = reverse p
    num = maximum (map magnitude (tail revp))
    den = magnitude (head revp)
    r = 1 + num / den
    — choose initial points equally spaced on a circle with radius based on bounds
    — for the absolute values of the roots
    — for a polynomial with real coefficients it would make sense to start
    — with conjugate pairs, since this is preserved by the iteration
    — I suspect the literature on this method will have much more to say
    — about how to choose good initial values
    circleInit :: RealFloat a => [Complex a] -> [Complex a]
    circleInit p = [ mkPolar (r/2) (((fromIntegral i)+0.5)*step) | i <- [0..n-1] ]
    where
    step = 2*pi/(fromIntegral n)
    r = rootBound p
    n = (length p) – 1
    phi = ((sqrt 5) + 1)/2
    goldenAngle = 2 * pi / phi**2
    — similar to circleInit
    — the golden angle is used to generate a quasi-uniform distribution
    — this is still deterministic, but a bit more like random
    goldenCircleInit :: [Complex Double] -> [Complex Double]
    goldenCircleInit p = [ mkPolar (r/2) ((fromIntegral i)*goldenAngle) | i <- [0..n-1] ]
    where
    r = rootBound p
    n = (length p) – 1
    roots :: [Complex Double] -> Int -> Double -> [Complex Double]
    –roots p tol = solve p (circleInit p) tol
    roots p maxit tol = solve p (goldenCircleInit p) maxit tol
    rootpoly' :: Num a => [a] -> [a] -> [a]
    rootpoly' [] p = p
    rootpoly' (z:zs) p = rootpoly' zs (zipWith (-) (0:p) ((map (z*) p)++[0]))
    — |find the coefficients of a polynomial with given roots
    rootpoly :: Num a => [a] -> [a]
    rootpoly zs = rootpoly' zs [1]
    test1 = (z0s, ss, ws)
    where
    –p = [3:+0, 2:+0, 1:+0]
    p = [3, 2, 1]
    z0s = goldenCircleInit p
    p' = diff p
    ss = scales z0s
    fs = map (polyval p) z0s
    gs = map (polyval p') z0s
    ws = offsets fs gs ss
    test2 = zs
    where
    p = rootpoly [1:+0, 2:+0, 3:+0]
    zs = roots p 30 1e-3
    test3 n = zs
    where
    p = rootpoly [ i:+i^2 | i <- [1..n] ]
    zs = roots p 30 1e-3

    view raw

    alberth_method.hs

    hosted with ❤ by GitHub

    Reply

  7. Pingback: What’s the right way to QuickCheck floating-point routines? | blog :: Brent -> [String]

Leave a comment

This site uses Akismet to reduce spam. Learn how your comment data is processed.