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 a matrix-vector product with the weight matrix:

Because every row of has exactly k nonzeros (one per stencil neighbor), the weights are stored stencil-wise in ELL format (StencilWeights) and the product is evaluated as a row-parallel gather kernel,    — multithreaded (with SIMD) on CPU and a KernelAbstractions kernel on GPU backends. This is roughly 7.6× faster than the sparse matrix-vector product it replaced (measured at  ,  , 13 threads).

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
│   ├── VirtualPartial   ∂f/∂xᵢ via offset-point finite differences
│   └── 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 = rand(SVector{2,Float64}, 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 StencilWeights (or tuple 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!.

Getting a sparse matrix

weights(op) returns the operator's StencilWeights — a tuple of them for gradient-family operators. When you need a standard sparse matrix, e.g. for global system assembly or an implicit solve, convert with sparse(op) (or SparseMatrixCSC(op)); this is the supported path, and the result can go straight into sparse(op) \ rhs or any sparse linear-algebra library:

using SparseArrays
A = sparse(op)  # SparseMatrixCSC, ready for global assembly or `A \ rhs`
typeof(A)

SparseArrays.SparseMatrixCSC{Float64, Int64}

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)

3.7107950630558952

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.

This is also why hand-written operators that compose functors with arithmetic need two methods — see Dual dispatch for composed functors for the user-facing consequence.

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.

Virtual Operators

Virtual operators (∂virtual) take a different route to a derivative: instead of solving for weights that approximate directly, they interpolate the field at points offset by along the axis and apply a standard finite difference formula to those interpolated values. This is useful for schemes that need the derivative to be consistent with a particular FD stencil.

# Virtual partial derivative in x-direction with spacing Δ=0.01
virtual_dx = ∂virtual(x, 1, 0.01)
result = virtual_dx(u)
size(result)

(100,)

The offset is a genuine tuning parameter: too large and the finite difference truncation error dominates; too small and interpolation error is amplified by the factor.

Next Steps

• Building PDE Operators — compose operators with @operator

• Custom Operators — the closure escape hatch for anything the macro can't express

• Boundary Conditions — Hermite interpolation at boundary nodes