Use FastGaussQuadrature.jl for computation of Gauss quadrature nodes and weights?

[BenjaminBorn]

It seems that there is a considerable performance difference between, e.g., qnwlege() from this package and gaussianlegendre() from the FastGaussQuadrature package.

using QuantEcon, FastGaussQuadrature

@time qnwlege(10000,-1,1)
@time gausslegendre(10000)

1.424605 seconds (300.10 k allocations: 5.601 GiB, 30.51% gc time)
0.000436 seconds (15 allocations: 234.875 KiB)

The timings are from the second run of the code.

A caveat is that FastGaussianQuadrature doesn't allow for integrals over [a, b] but that could be dealt with.

[sglyon]

Thanks for letting us know about this!

I think you are right that we could wrap the FastGaussQuadrature routines in a function that shifts/scales the nodes and weights they supply to be on [a, b]

I personally don't have the time to attack this right now, but if you are up for it I could try to offer some suggestions on where to start.

Thanks again

[tbeason]

So I did a little bit of this in my own work just because I was trying to cut out dependencies. Here is what I did, if that helps get things rolling. Benchmark shows it is faster than QE for small n, but does slow down as n grows. I suspect somebody more experienced than I am could fix that.

# qnwnorm (like from QuantEcon/CompEcon) replacement
function qnwnorm_new(n::Int,mu::Real,variance::Real)
	x,w = gausshermite(n)
	x = mu + sqrt(2*variance)*x
	w = w/sqrt(pi)
	return x::Vector{Float64},w::Vector{Float64}
end

[sglyon]

Thanks @tbeason. Did you collect benchmark data? If so, would you mind sharing it.

When I try to run it on n=10 nodes I see almost no difference in performance:

julia> @benchmark qnwnorm_new(10, 5.0, 4.0)
BenchmarkTools.Trial:
  memory estimate:  8.88 KiB
  allocs estimate:  87
  --------------
  minimum time:     19.567 μs (0.00% GC)
  median time:      21.022 μs (0.00% GC)
  mean time:        24.029 μs (5.45% GC)
  maximum time:     3.957 ms (97.16% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark qnwnorm(10, 5.0, 4.0)
BenchmarkTools.Trial:
  memory estimate:  5.41 KiB
  allocs estimate:  134
  --------------
  minimum time:     20.339 μs (0.00% GC)
  median time:      23.227 μs (0.00% GC)
  mean time:        26.142 μs (3.77% GC)
  maximum time:     5.151 ms (98.52% GC)
  --------------
  samples:          10000
  evals/sample:     1

[BenjaminBorn]

I can confirm that they are similar in performance (with a slight edge for qnwnorm)

@benchmark qnwnorm(10, 5.0, 4.0)

BenchmarkTools.Trial: 
  memory estimate:  5.41 KiB
  allocs estimate:  134
  --------------
  minimum time:     20.247 μs (0.00% GC)
  median time:      23.601 μs (0.00% GC)
  mean time:        28.616 μs (3.78% GC)
  maximum time:     5.987 ms (97.05% GC)
  --------------
  samples:          10000
  evals/sample:     1

@benchmark qnwnorm_new(10, 5.0, 4.0)

BenchmarkTools.Trial: 
  memory estimate:  9.28 KiB
  allocs estimate:  94
  --------------
  minimum time:     21.954 μs (0.00% GC)
  median time:      25.949 μs (0.00% GC)
  mean time:        35.466 μs (3.75% GC)
  maximum time:     4.601 ms (95.28% GC)
  --------------
  samples:          10000
  evals/sample:     1

  @benchmark qnwnorm(50, 5.0, 4.0)

  BenchmarkTools.Trial: 
  memory estimate:  200.89 KiB
  allocs estimate:  998
  --------------
  minimum time:     181.330 μs (0.00% GC)
  median time:      201.811 μs (0.00% GC)
  mean time:        271.024 μs (8.00% GC)
  maximum time:     4.435 ms (76.89% GC)
  --------------
  samples:          10000
  evals/sample:     1

  @benchmark qnwnorm_new(50, 5.0, 4.0)

  BenchmarkTools.Trial: 
  memory estimate:  200.89 KiB
  allocs estimate:  998
  --------------
  minimum time:     183.272 μs (0.00% GC)
  median time:      219.670 μs (0.00% GC)
  mean time:        298.312 μs (7.69% GC)
  maximum time:     7.164 ms (89.91% GC)
  --------------
  samples:          10000
  evals/sample:     1

And it's not just the shifting/rescaling. The pure function is also not faster.

  @benchmark gausshermite(50)
  
  BenchmarkTools.Trial: 
  memory estimate:  199.22 KiB
  allocs estimate:  991
  --------------
  minimum time:     182.641 μs (0.00% GC)
  median time:      209.245 μs (0.00% GC)
  mean time:        293.728 μs (8.21% GC)
  maximum time:     7.643 ms (86.71% GC)
  --------------
  samples:          10000
  evals/sample:     1

For Gauss-Legendre on the other hand, the difference is considerable

@benchmark qnwlege(10,-1,1)

  BenchmarkTools.Trial: 
    memory estimate:  32.22 KiB
    allocs estimate:  258
    --------------
    minimum time:     9.961 μs (0.00% GC)
    median time:      11.401 μs (0.00% GC)
    mean time:        15.917 μs (15.83% GC)
    maximum time:     2.749 ms (96.53% GC)
    --------------
    samples:          10000
    evals/sample:     1

  @benchmark gausslegendre(10)

  BenchmarkTools.Trial: 
    memory estimate:  1.66 KiB
    allocs estimate:  17
    --------------
    minimum time:     1.561 μs (0.00% GC)
    median time:      1.785 μs (0.00% GC)
    mean time:        2.365 μs (12.38% GC)
    maximum time:     502.211 μs (98.03% GC)
    --------------
    samples:          10000
    evals/sample:     10

  @benchmark qnwlege(10000,-1,1)

  BenchmarkTools.Trial: 
  memory estimate:  5.60 GiB
  allocs estimate:  300093
  --------------
  minimum time:     1.537 s (37.76% GC)
  median time:      1.748 s (38.27% GC)
  mean time:        1.849 s (38.79% GC)
  maximum time:     2.261 s (39.89% GC)
  --------------
  samples:          3
  evals/sample:     1

  @benchmark gausslegendre(10000)

  BenchmarkTools.Trial: 
    memory estimate:  234.72 KiB
    allocs estimate:  11
    --------------
    minimum time:     346.701 μs (0.00% GC)
    median time:      361.636 μs (0.00% GC)
    mean time:        434.014 μs (4.73% GC)
    maximum time:     4.264 ms (86.25% GC)
    --------------
    samples:          10000
    evals/sample:     1

It might make sense to do a comprehensive benchmarking exercise of the QuantEcon vs. FaustGaussQuadrature routines.

[sglyon]

Glad someone confirmed that qnwnorm isn't significantly slower than the competition, thank you @BenjaminBorn for doing that

As for qnwlege, I honestly think that even though the difference in time is large, it probably isn't a huge deal in practice. When I use these routines I typically have relatively few quadrature weights and nodes (between 5 and 30), so the difference in time reported by your first example when n=10 becomes the relevant consideration for me. I'm usually ok with an 8 micro-second performance change.

That being said, I would be open to seeing a new version of qnwlege that leverages FaustGaussQuadrature under the hood as long as it behaves the same as qnwlege for all arguments.