Repository navigation
Gauss-Turán Quadratures #245
Description
Activity
∑wᵢ*∑(aᵢₖ∂ᵏf(xᵢ))Note that at a high level
∑wᵢ*∑(aᵢₖ∂ᵏf(xᵢ)) = ∑wᵢ* g(xᵢ)so this can be implemented as a meta-algorithm, e.g.GaussTuránRule(sensealg, n, k) == Integrals.ChangeOfVariables((f, domain) -> (gauss_turán_integrand(f, n, k, sensealg), domain), Integrals.QuadratureRule(gauss_turán_weights_nodes, n=n))where some work would have to be done to implementgauss_turán_weights_nodes(which could simply tabulatewᵢ, xᵢfrom the paper) andgauss_turán_integrandwhich would have to compute gradients offwith an AD backend selected fromsensealg(e.g. ForwardDiff.jl) and values ofaᵢₖ(which again could be tabulated). Here,Integrals.QuadratureRuleis documented and in the public API andIntegrals.ChangeOfVariablesis internal and recently merged to the main branch, and not actually necessary for this case but a reference as another meta-algorithm.This would be an interesting contribution for anyone interested in:
- automating the calculation of higher-order derivatives of integrands
- tabulating known Gauss-Turán rules
- optionally implementing a numerical method for the calculation of Gauss-Turán rules
So I went down quite the rabbit hole to implement an algorithm for computing Gauss-Turán quadratures 😄 Find my code here.
Here is the example from [1]:
using DoubleFloats FT = Double64 n = 5 s = 1 a = FT(0.0) b = FT(1.0) derivs! = @eval @generate_derivs($(2s)) """ Example functions ϕⱼ from [1] """ function f(x::T, j::Int)::T where {T<:Number} pow = if j % 2 == 0 j / 2 - 1 else (j - 1) / 2 - 1 / 3 end x^pow end deriv! = @eval @generate_derivs($(2s)) I = GaussTuran( f, deriv!, a, b, n, s; optimization_options = Optim.Options(; x_abstol = FT(1e-250), g_tol = FT(1e-250), show_trace = true, show_every = 100, iterations = 10_000, ), ) # Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ max_int_error_upper = maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) for j = 1:I.cache.N) @show max_int_error_upper # max_int_error_upper = 3.0814879110195774e-32 # Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ max_int_error_lower = maximum( abs(I(x -> f(x, j)) - I.cache.rhs_lower[j-I.cache.N]) for j = (I.cache.N+1):(I.cache.N+I.cache.n) ) @show max_int_error_lower # max_int_error_lower = 2.8858134286698342e-30 # Example with eˣ exp_error = abs(I(Base.exp) - (Base.exp(1) - 1)) @show exp_error; # exp_error = -6.436621416953051e-18
The weights for the quadrature rule are not shown here but they are known (
I.A). My results for the nodes differ from those in [1], but the integrals are correct so I claim that this is because I have achieved higher accuracy :) I am not familiar with theIntegrals.jlAPI so I haven't incorporated that yet.I didn't implement the construction algorithm from [1], but derived a loss function which I think is simpler given that we can use automatic differentiation these days.
automating the calculation of higher-order derivatives of integrands
I struggled with that, I currently hardcode the derivatives of the integrands. I tried something like this, but that is a nightmare in terms of type stability.
That looks like a really good start!
Reacted by Bart de KoningIndeed, it is great that you demonstrate that your nodes and weights get the same accuracy for the first
2(s+1)*npolynomials, which is something I do not get for then=5,s=1rule in [1], for which I observe a degradation in accuracy after the first2npolynomials. I suppose the optimization problem is rather costly, so make sure to save your results! Packages like GeneralizedGauss.jl tabulate them as well.Regarding cost: I updated my code and the example to work with any
AbstractFloattype. I learned thatBigFloatmakes for very slow code because it leads to a lot of allocation even for scalar operations. Obtaining the results shown above took only half a minute, and I'm sure there's room for optimization.I updated the code and the example above again. I added a macro for computing the first
$n$ derivatives of a scalar function at once, maybe a nice feature forForwardDiff.jl@ChrisRackauckas?The only thing I don't understand is why initializing the optimization problem takes so long, more specifically computing the first Hessian.
The code needs to be shared to see, that
@evalis very suspect though. Is there a reason to not use TaylorDiff for this part? ForwardDiff stacking has bad complexity.
What kind of problems is it mostly used for? Please describe.
As most of those familiar with Gauss-Legendre quadratures, a one-dimensional integral can be aproximated to
∑(f(xᵢ)*wᵢ). wherexᵢandwᵢdepend on various factors. this equation notably uses just zero-order derivative information (just the function evaluated at the nodes). Gauss-Turán quadratures are an extension to this scheme, allowing the use of higher derivative information. This allows to higher precision with less nodes (given that you also have derivative information at the nodes).Describe the algorithm you’d like
The Gauss-Turán Quadrature is just expressed as a double sum:
where
∂ᵏf(xᵢ)is just the k-th derivative evaluated at the i-th node, andwᵢ,aᵢₖare just precomputed weights.The main algorithm needed is the one that calculates those weights and nodes, for one or more polynomial interpolants (like Legendre or Chebyshev).
Other implementations to know about
I haven't found other implementations online, but i did whip up an example with the weights and nodes in [1]. the code below is for a 5-point, k = 2 (value at the point, first and second derivative, the paper allows odd values only, so in the paper notation, it is 5-node, n = k/2 = 1 order quadrature)
References
[1] Kovačević, M. A. (2000). Construction of generalized Gauss-Turán quadrature rules for systems of arbitrary functions. Computers & Mathematics with Applications (Oxford, England: 1987), 40(8–9), 895–909. doi:10.1016/s0898-1221(00)85001-4
[2] Gori, L., & Micchelli, C. (1996). On weight functions which admit explicit Gauss-Turán quadrature formulas. Mathematics of computation, 65(216), 1567-1581. link