Operators & Type Hierarchy

Operators are the core abstraction in RadialBasisFunctions.jl for RBF-FD differentiation on scattered data. This page explains the operator type system, rank semantics, and how everything fits together.

For basic usage, see Getting Started. For the underlying math, see Radial Basis Functions Theory.

using RadialBasisFunctions
using StaticArrays

Math Refresher

An RBF-FD operator approximates a differential operator at a point using a weighted sum over its local stencil:

The weights are precomputed by solving a local collocation system (see Radial Basis Functions Theory for the full derivation). Once computed, applying the operator is just a sparse matrix-vector multiply:

This is what RadialBasisOperator stores and evaluates.

The Type Hierarchy

All operators inherit from AbstractOperator, where the parameter N is the tensor rank added to the output:

AbstractOperator{N}
├── N=0 (rank-preserving)
│   ├── Partial          ∂ⁿf/∂xᵢⁿ
│   ├── MixedPartial     ∂²f/∂xᵢ∂xⱼ
│   ├── Laplacian        ∇²f
│   ├── Directional      ∇f⋅v
│   ├── Divergence       ∇⋅u (vector field → scalar)
│   ├── Curl             ∇×u (vector field → scalar/vector)
│   ├── Identity         f (function itself)
│   ├── ScaledOperator   α * op (algebra result, any rank N)
│   ├── SumOperator      op₁ + op₂ (algebra result, any rank N)
│   ├── StrainRate       ½(∇u + (∇u)ᵀ)
│   ├── RotationRate     ½(∇u − (∇u)ᵀ)
│   ├── Regrid           interpolation to new points
│   └── Custom{0}        user-defined
├── N=1 (rank-adding)
│   ├── Jacobian          [∂fᵢ/∂xⱼ]
│   └── Custom{1}         user-defined
└── N=2 (rank+2)
    └── Hessian           [∂²f/∂xᵢ∂xⱼ]

Understanding Rank (N)

The parameter N captures whether differentiation adds a tensor index to the output.

N=0 (rank-preserving): The output has the same shape as the input. The operator stores a single weight matrix W, and evaluation is W * u.

x = [SVector{2}(rand(2)) for _ in 1:100]
u = sin.(getindex.(x, 1))

lap = laplacian(x)
result = lap(u)
size(result)  # (100,) — same shape as input

(100,)

N=1 (rank-adding): The output gains a trailing dimension of size D (spatial dimension). The operator stores a tuple of D weight matrices (W₁, W₂, …, W_D), one per spatial dimension.

jac = jacobian(x)
result = jac(u)
size(result)  # (100, 2) — trailing dimension added

(100, 2)

When a rank-1 operator is applied to a vector field (matrix input), the output gains yet another dimension:

v = hcat(sin.(getindex.(x, 1)), cos.(getindex.(x, 2)))  # (100, 2) vector field
result = jac(v)
size(result)  # (100, 2, 2) — full Jacobian tensor

(100, 2, 2)

Input/Output Shape Summary

Operator rank	Input shape	Output shape	Example
N=0	(N,)	(N,)	laplacian, partial
N=0	(N, D)	(N, D)	laplacian on vector field
N=0	(N, D)	(N,)	divergence (vector field → scalar)
N=0	(N, 2)	(N,)	curl in 2D (vector field → scalar)
N=0	(N, 3)	(N, 3)	curl in 3D (vector field → vector)
N=1	(N,)	(N, D)	jacobian on scalar field
N=1	(N, D)	(N, D, D)	jacobian on vector field
N=2	(N,)	(N, D, D)	hessian on scalar field

RadialBasisOperator: The Wrapper

RadialBasisOperator wraps an operator with everything needed to compute and apply it:

op = laplacian(x)

RadialBasisOperator
├─Operator: Laplacian (∇²f)
├─Properties: self-adjoint
├─Data type: StaticArraysCore.SVector{2, Float64}
├─Number of points: 100
├─Stencil size: 12
└─Basis: Polyharmonic spline (r³) with degree 2 polynomial augmentation

Key fields:

Field	Description
ℒ	The operator type (e.g., Laplacian())
weights	Precomputed weight matrix (or tuple of matrices for multi-component operators)
data	Source points used to build stencils
eval_points	Points where the operator is evaluated
adjl	Adjacency list (neighbor indices per stencil)
basis	RBF basis function used

Weight Caching and Invalidation

Weights are computed eagerly during construction and cached. If you mutate the underlying data (e.g., move points), invalidate the cache to trigger recomputation on next evaluation:

# Manually invalidate if data changes
RadialBasisFunctions.invalidate_cache!(op)

# Next call recomputes weights automatically
result = op(u)
typeof(result)

Vector{Float64} (alias for Array{Float64, 1})

You can also force an immediate recomputation with update_weights!.

Basis Derivative Functors

When you call an operator type on a basis, it returns a functor — a callable struct that evaluates the differentiated basis function at two points (x, xᵢ). These functors are the building blocks for both built-in and custom operators.

basis = PHS(3; poly_deg=2)

# Laplacian() applied to a basis returns a ∇² functor
lap_functor = Laplacian()(basis)
typeof(lap_functor)

RadialBasisFunctions.∇²{PHS3{Int64}}

# Partial(1, 1) applied to a basis returns a ∂ functor
partial_functor = Partial(1, 1)(basis)
typeof(partial_functor)

RadialBasisFunctions.∂{PHS3{Int64}}

These functors are callable as (x, xᵢ) -> scalar:

x1 = SVector(0.5, 0.3)
x2 = SVector(0.1, 0.2)

# Evaluate ∇²ϕ(‖x₁ - x₂‖)
lap_functor(x1, x2)

4.947726750741193

The Jacobian operator returns a tuple of functors (one per spatial dimension):

jac_functors = Jacobian{2}()(basis)
typeof(jac_functors)

Tuple{RadialBasisFunctions.∂{PHS3{Int64}}, RadialBasisFunctions.∂{PHS3{Int64}}}

Available functor types (accessed via RadialBasisFunctions.∂ etc.):

Functor	Constructor	Signature
∂	∂(basis, dim)	(x, xᵢ) -> scalar
∂²	∂²(basis, dim)	(x, xᵢ) -> scalar
∂mixed	∂mixed(basis, dim1, dim2)	(x, xᵢ) -> scalar
∇	∇(basis)	(x, xᵢ) -> vector
∇²	∇²(basis)	(x, xᵢ) -> scalar
H	H(basis)	(x, xᵢ) -> matrix

These functors are the interface between operators and Custom Operators. See that page for how to use them.

Two differentiation protocols

The names ∂, ∂², ∇, ∇², H, and ∂mixed deliberately serve two protocols, selected by the basis type:

Basis	What ∂(basis, dim) returns	How it evaluates
AbstractRadialBasis (PHS, IMQ, Gaussian)	A functor struct (∂, ∇², …, defined in src/basis/basis.jl)	(x, xᵢ) -> scalar
MonomialBasis	A factory function result: ℒMonomialBasis (defined in src/operators/monomial/monomial.jl)	in-place (b, x) fill of the differentiated monomial vector

This dual identity is load-bearing: an operator action like (op::Partial)(basis) = ∂(basis, op.order, op.dim) is written once, and weight building applies it to the RBF basis for the collocation rows and to the MonomialBasis for the polynomial-augmentation rows — each returning the form the assembly kernels expect. When adding a basis or operator, implement both halves; do not try to unify them.

Operator Algebra

Built RadialBasisOperators can be combined with + and -. This operates on precomputed weights and returns a new operator:

∂x = partial(x, 1, 1)
∂y = partial(x, 1, 2)

combined = ∂x + ∂y  # ∂f/∂x + ∂f/∂y
result = combined(u)
typeof(result)

Vector{Float64} (alias for Array{Float64, 1})

Both operands must share the same data, stencils, and rank N.

See Custom Operators for more on building your own operators.