I want to evaluate and differentiate functions on a grid of points

that have an underlying tensor structure

g(x,y) = Lᵢ(xₖ) Aᵢⱼ Lⱼ(yₗ)

Here x = [x₁, …, xₙ] and y = [x₁, …, xₙ]

Here A is a n x n Matrix of coefficients and L(x) = [L₁(x) … Lₙ(x)]

and L(y) = [L₁(y) … Lₙ(y)] are sets of functions (e.g. Lagrange polynomials).

First consider evaluatation!

- evaluation of g(x,y) using matrix matrix products

number of flops ∝ n³ - evaluation of g(xₖ,xₗ) for all points in x and y one by one

number of flops ∝ n⁴

For algorithmic differentiation the flops counts are identical. It seems to me that the efficient version (1) is currently not possible in ForwardDiff

What would it take to add this kind of support?