Cubic Adjoint (1D)

The CubicAdjoint operator computes $\bar{\mathbf{f}} = W^\top \bar{\mathbf{y}}$ for 1D cubic spline interpolation — mapping query-space sensitivities back to grid-space.

It is query-baked and data-free: constructed from grid x and query points xq only. The same operator can be applied to any $\bar{\mathbf{y}}$ vector.

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 the adjoint operator (grid + queries only, no data)
adj = 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(ȳ)

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

Verify correctness via the dot-product identity $\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(ȳ)                   # 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

cubic_adjoint(x, x_query; bc=CubicFit()) -> CubicAdjoint

x: Grid points (sorted AbstractVector)
x_query: Query points (baked into the operator)
bc: Any boundary condition — all BC types are supported including PeriodicBC

See Adjoint API Reference for full docstrings.

Calling the Operator

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(f̄_buf, ȳ)
f̄_buf[1:5]

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

Size

size(adj)          # (n_grid, n_query)

(50, 5)

Matrix Materialization

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

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

(50, 5)

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

Matrix × ȳ ≈ adj(ȳ): true

The transpose gives the forward interpolation matrix:

W = Matrix(adj)'
@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

Construction: $O(n + m \log n)$ — Thomas factorization + anchor binary search (amortized via autocache)
Apply (per call): $O(n + m)$ — one transpose solve + one scatter pass
In-place: Zero allocation after construction
Cache sharing: Reuses the same CubicSplineCache as forward interpolation