feat: Implement analytical polynomial coefficient formulas (Skala 2026)

[kylebeggs]

Summary

This issue tracks potential improvements to RadialBasisFunctions.jl based on analytical formulas for polynomial coefficients in RBF interpolation. The key insight is that polynomial coefficients in RBF interpolation can be computed analytically without solving the full linear system, enabling new diagnostics and optimizations.

Reference

Paper: Skala, V. (2026). "Analytical formulas for polynomial coefficients in radial basis function interpolation"

Key Result: For RBF interpolation with polynomial augmentation, the polynomial coefficients can be expressed analytically:

a₀ = (1/N) Σᵢ fᵢ                    # constant term = mean of data
a₁ = (XᵀX)⁻¹ Xᵀf                    # linear terms = least squares fit

where X is the Vandermonde-like matrix of evaluation points.

Proposed Improvements

1. Stencil Quality Diagnostics

What: Add optional diagnostic that compares numerically computed a₀ against analytical mean(f).

Why: Large deviation indicates numerical issues (ill-conditioning, near-singular stencils).

Implementation:

function stencil_quality(operator::RadialBasisOperator, f::AbstractVector)
    # Compare numerical a₀ vs analytical mean(f) per stencil
    # Return vector of relative errors
end

Current state: We have no stencil quality diagnostics currently.

2. Optional Polynomial Coefficient Retention

What: Allow operators to optionally return polynomial coefficients alongside RBF weights.

Why:

• The a₁ coefficients approximate the gradient (useful for visualization/debugging)
• Enables hybrid schemes that use polynomial part for extrapolation

Implementation:

struct RadialBasisOperator{...}
    # existing fields...
    poly_coeffs::Union{Nothing, Vector}  # optional storage
end

function (op::RadialBasisOperator)(f; return_poly_coeffs=false)
    # existing evaluation...
    if return_poly_coeffs
        return result, poly_coeffs
    end
    return result
end

Current state: We discard polynomial coefficients after solving.

3. Preconditioning Infrastructure

What: Use analytical formulas to construct preconditioners for iterative solvers.

Why: The analytical structure provides natural scaling that could improve conditioning.

Implementation:

function build_preconditioner(points, basis)
    # Use analytical coefficient formulas to construct
    # block-diagonal or other structured preconditioner
end

Current state: We use direct solves only (\ operator). No iterative solver support.

4. Boundary Condition Improvements

What: Leverage analytical polynomial structure for better boundary treatment.

Why: Near boundaries, polynomial part dominates; analytical formulas could improve extrapolation.

Current state: We have Hermite interpolation support in src/solve/ but polynomial coefficients aren't exposed.

Priority Assessment

Improvement	Complexity	Value	Priority
Stencil diagnostics	Low	Medium	P1
Poly coefficient retention	Medium	Medium	P2
Preconditioning	High	High (for large problems)	P3
Boundary improvements	High	Medium	P4

Related Code

• Weight computation: src/solve/assembly.jl
• Operator structure: src/operators/operators.jl:10-31
• Hermite interpolation: src/solve/api.jl

Next Steps

1. Implement stencil quality diagnostic as proof-of-concept
2. Benchmark to determine if overhead is acceptable
3. Consider exposing polynomial coefficients in operator API