Adjoint Operators (1D)

FastInterpolations provides native adjoint operators for all 1D interpolant types. Each operator is query-baked and data-free: constructed from grid and query points only, then applied to any sensitivity vector $\bar{\mathbf{y}}$.

Constructor	Struct	Notes
`constant_adjoint(x, xq; ...)`	`ConstantAdjoint`	Pure scatter
`linear_adjoint(x, xq; ...)`	`LinearAdjoint`	Pure scatter
`quadratic_adjoint(x, xq; ...)`	`QuadraticAdjoint`	Recurrence adjoint + scatter
`cubic_adjoint(x, xq; ...)`	`CubicAdjoint`	Transpose Thomas solve + scatter

Quick Start

using FastInterpolations
using LinearAlgebra: dot

x = collect(range(0, 1, 50))
xq = [0.1, 0.25, 0.5, 0.75, 0.9]
f = sin.(2π .* x)

# Build adjoint operators (grid + queries only, no data)
adj_const  = constant_adjoint(x, xq)
adj_lin    = linear_adjoint(x, xq)
adj_quad   = quadratic_adjoint(x, xq)
adj_cubic  = cubic_adjoint(x, xq)

CubicAdjoint{Float64}
├─ Grid:  Vector, 50 points ∈ [0, 1]
├─ Query:  5 points
└─ BC:     CubicFit

# Apply to a sensitivity vector
ȳ = collect(1.0:length(xq))
f̄ = adj_cubic(ȳ)

50-element Vector{Float64}:
  0.00018469921830076248
 -0.0014248225411773103
  0.0064117014352978955
 -0.024644152841844597
  0.09313852432946604
  0.9791428267292088
 -0.06670983124630073
  0.0166964982559946
 -7.616177767766135e-5
 -0.016391851145283968
  ⋮
 -0.014594242814458095
  0.0872909694569015
 -0.3345696350131474
  4.895987570595675
  0.46561935263044585
 -0.12320116705522254
  0.03205340855933533
 -0.007122979679852251
  0.0009233492177586249

Dot-Product Identity

Verify $\langle W f, \bar{y} \rangle = \langle f, W^\top \bar{y} \rangle$:

itp = cubic_interp(x, f)
Wf = itp(xq)                  # forward: W·f
WTy = adj_cubic(ȳ)            # adjoint: Wᵀ·ȳ

@assert dot(Wf, ȳ) ≈ dot(f, WTy)
println("⟨Wf, ȳ⟩ = ", dot(Wf, ȳ))
println("⟨f, Wᵀȳ⟩ = ", dot(f, WTy))

⟨Wf, ȳ⟩ = -4.351140062672254
⟨f, Wᵀȳ⟩ = -4.351140062672254

Constructor Signatures

constant_adjoint(x, xq; side=NearestSide(), extrap=NoExtrap())
linear_adjoint(x, xq; extrap=NoExtrap())
quadratic_adjoint(x, xq; bc=Left(QuadraticFit()), extrap=NoExtrap())
cubic_adjoint(x, xq; bc=CubicFit(), extrap=NoExtrap())

Each constructor accepts the same keyword arguments as its forward interpolation counterpart. See Adjoint API Reference for full docstrings.

Calling the Operator

All adjoint operators share a unified calling convention:

Allocating

f̄ = adj(ȳ)                          # value adjoint
f̄ = adj(ȳ; deriv=DerivOp(1))        # 1st derivative adjoint
f̄ = adj(ȳ; deriv=DerivOp(2))        # 2nd derivative adjoint

In-Place (Zero Allocation)

adj(f̄, ȳ)                           # zeros f̄, then accumulates
adj(f̄, ȳ; deriv=DerivOp(1))         # derivative adjoint, in-place

f̄_buf = zeros(length(x))
adj_cubic(f̄_buf, ȳ)
f̄_buf[1:5]

5-element Vector{Float64}:
  0.00018469921830076248
 -0.0014248225411773103
  0.0064117014352978955
 -0.024644152841844597
  0.09313852432946604

Size

size(adj_cubic)          # (n_grid, n_query)

(50, 5)

Matrix Materialization

For verification and debugging, materialize the full weight matrix. This is $O(n \times m)$ and not intended for production use.

Wᵀ = Matrix(adj_cubic)
size(Wᵀ)          # (n_grid, n_query)

(50, 5)

@assert Wᵀ * ȳ ≈ adj_cubic(ȳ)
println("Matrix × ȳ ≈ adj(ȳ): ", isapprox(Wᵀ * ȳ, adj_cubic(ȳ)))

Matrix × ȳ ≈ adj(ȳ): true

The transpose gives the forward interpolation matrix:

W = Matrix(adj_cubic)'
@assert W * f ≈ itp(xq)
println("W * f ≈ itp(xq): ", isapprox(W * f, itp(xq)))

W * f ≈ itp(xq): true

A convenience method builds the forward matrix directly from an interpolant:

W_direct = Matrix(itp, xq)
@assert W_direct * f ≈ itp(xq)
println("Matrix(itp, xq) * f ≈ itp(xq): ", isapprox(W_direct * f, itp(xq)))

Matrix(itp, xq) * f ≈ itp(xq): true

Affine BCs

For BCs with non-zero values (e.g., Deriv1(0.5)), the forward operator is affine: $y = W f + c$. The adjoint computes $W^\top \bar{y}$ — the constant $c$ does not appear.

Performance

Method	Construction	Apply (per call)
Constant / Linear	$O(m \log n)$ — anchor search	$O(m)$ — pure scatter
Quadratic	$O(n + m \log n)$ — anchor + slope setup	$O(n + m)$ — recurrence adjoint + scatter
Cubic	$O(n + m \log n)$ — Thomas factorization + anchor search	$O(n + m)$ — transpose solve + scatter

In-place: Zero allocation after construction (all types)
Cache sharing: Cubic and quadratic reuse the same spline caches as forward interpolation