∑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 implement gauss_turán_weights_nodes (which could simply tabulate wᵢ, xᵢ from the paper) and gauss_turán_integrand which would have to compute gradients of f with an AD backend selected from sensealg (e.g. ForwardDiff.jl) and values of aᵢₖ (which again could be tabulated). Here, Integrals.QuadratureRule is documented and in the public API and Integrals.ChangeOfVariables is 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 the Integrals.jl API 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!

Indeed, it is great that you demonstrate that your nodes and weights get the same accuracy for the first 2(s+1)*n polynomials, which is something I do not get for the n=5,s=1 rule in [1], for which I observe a degradation in accuracy after the first 2n polynomials. 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 AbstractFloat type. I learned that BigFloat makes 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 for ForwardDiff.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 @eval is very suspect though. Is there a reason to not use TaylorDiff for this part? ForwardDiff stacking has bad complexity.