example code: arbitrary-precision Gaussian quadrature weights in Julia

[Steven G. Johnson]

I was playing around with this code this evening and thought I would post it as a potentially useful code snippet.  It computes Gaussian quadrature nodes and weights for arbitrary number types, including BigFloat, which allows you to perform arbitrary-precision integrals with exponential accuracy (for analytic integrands).

It would be pretty trivial except that Julia does not yet have a BigFloat eig() function for arbitrary precision. I didn't implement a general routine though; this problem has some special properties that permit a few shortcuts (and has overall Õ(N^2 * p * polylog(p)) compiexity in the number N of quadrature weights and the precision p, not including logarithmic factors).

e.g.
     (x, w) = gauss(BigFloat, 30); dot(cos(x), w)
will compute the integral of cos(x) from -1 to 1 to the default BigFloat precision (256 bits, about 1e-77 error).

(Thanks to Alessandro Andrioni for the greatly improved BigFloat support based on MPFR.  This code may or may not work with Julia 0.1, I haven't checked.)

Steven

PS. I also have working routines to get Kronrod points/weights as well (for Gauss-Kronrod adaptive quadrature in arbitrary precision), but since I based it on http://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html out of laziness I can't distribute the Kronrod routines as free/open-source software.  (OPQ was published in the semi-evil TOMS journal, which requires authors to assign ACM copyright to their code and then slaps a semi-free license on it.)  On the other hand, the computed quadrature weights & points are uncopyrightable mathematical facts, so I was thinking of precomputing a table at reasonably high precision and using it to distribute a Julia quadgk routine.

##################################################################################################################
# Copyright (c) 2013 Steven G. Johnson.
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# return (x, w) pair of N quadrature points x[i] and weights w[i] to           
# integrate functions on the interval (-1, 1).  i.e. dot(f(x), w)              
# approximates the integral.  Uses the method described in Trefethen &         
# Bau, Numerical Linear Algebra to find the N-point Gaussian quadrature 
function gauss(T, N::Integer)
    two = convert(T, 2)
    Ja = zeros(T, N)
    Jb = T[ 0.5 / sqrt(1 - (two*n)^-2) for n = 1:N-1 ]
    d = eignewt(Ja,Jb)
    v = Array(T,length(Ja))
    for i = 1:length(Ja)
        v[i] = eigvec1(Ja,Jb,d[i])[1]
    end
    return (d, 2*v.^2)
end

# Given a symmetric tridiagonal matrix H with H[i,i] = a[i] and                
# H[i-1,i] = H[i,i-1] = b[i-1], compute p(z) = det(H - z I) and its            
# derivative p'(z), returning (p,p').                                          
function eigpoly(a,b,z)
    m = length(a)
    d1 = a[1] - z
    d2 = one(z);
    d1deriv = convert(typeof(z), -1)
    d2deriv = zero(z);
    for i = 2:m
        b2 = b[i-1]^2
        az = (a[i] - z)
        d = az * d1 - b2 * d2;
        dderiv = az * d1deriv - b2 * d2deriv - d1;
        d2 = d1;
        d1 = d;
        d2deriv = d1deriv;
        d1deriv = dderiv;
    end
    return (d1, d1deriv)
end

# compute the eigenvalues of the symmetric tridiagonal matrix H                
# (defined from a,b as in eigpoly) using a Newton iteration                    
# on det(H - lambda I).  Unlike eig, handles BigFlaot.                         
function eignewt(a,b)
    m = length(a)

    # get initial guess from eig on Float64 matrix                             
    H = diagm(float64(a))
    for i in 2:m
        H[i-1,i] = H[i,i-1] = b[i-1]
    end
    lambda0 = sort(eigvals(H))

    lambda = similar(a)
    for i = 1:m
        lambda[i] = lambda0[i]
        for k = 1:1000
            (p,pderiv) = eigpoly(a,b,lambda[i])
            lambda[i] = (lamold = lambda[i]) - p / pderiv
            if abs(lambda[i] - lamold) < 10 * eps(lambda[i]) * abs(lambda[i])
                break
            end
        end
    end
    return lambda
end

# given an eigenvalue z and the matrix H(a,b) from above, return               
# the corresponding eigenvector, normalized to 1.                              
function eigvec1(a,b,z::Number)
    # "cheat" and use the fact that our eigenvector v must have a              
    # nonzero first entries (since it is a quadrature weight), so we           
    # can set v[1] = 1 to solve for the rest of the components:.               
    v = similar(a)
    v[1] = 1
    if length(a) > 1
        s = v[1]
        v[2] = -(a[1]-z)*v[1] / b[1]
        s += v[2]^2
        for i = 3:length(a)
            # b[i-2]*v[i-2] + (a[i-1]-z)*v[i-1] + b[i-1]*v[i] = 0              
            v[i] = - (b[i-2]*v[i-2] + (a[i-1]-z)*v[i-1]) / b[i-1]
            s += v[i]^2
        end
        scale!(v, 1 / sqrt(s))
    end
    return v
end

[Alessandro "Jake" Andrioni]

On 9 May 2013 00:43, Steven G. Johnson <steve...@gmail.com> wrote:

(Thanks to Alessandro Andrioni for the greatly improved BigFloat support based on MPFR.  This code may or may not work with Julia 0.1, I haven't checked.)

It won't be relevant on Julia 0.1, as the old BigFloat type was using the default precision (i.e. 53 bits). If there is demand, I can backport the new BigFloat in my unreleased MPFR.jl package to allow Julia 0.1 users to use them.

[Steven G. Johnson]

I ended up re-implementing the Kronrod-weight computation from scratch, and it's the basis of a new quadgk(f,a,b) adaptive Gauss-Kronrod integration scheme (including arbitrary-precision support) that I just submitted a pull request for:

    https://github.com/JuliaLang/julia/pull/3140

--SGJ