uncertain-numerics

Probabilistic numerical methods in Rust, with explicit uncertainty over computational quantities.

uncertain-numerics treats numerical computation as an inference problem. Instead of returning only a point estimate of an integral or of the solution of a linear system, every method returns a validated Gaussian posterior: its mean is the estimate and its variance states how much the computation still does not know. Each statistical and numerical assumption behind that posterior is explicit, tested against analytic references, and checked for calibration.

Highlights

• Bayesian quadrature in one dimension with a Gaussian-process prior (RBF kernel, Gaussian integration measure) and closed-form kernel integrals.
• Active Bayesian quadrature: a posterior-variance-reduction acquisition rule, deterministic candidate selection, a sequential evaluation loop, and explicit stopping rules.
• Probabilistic linear solvers for dense symmetric positive-definite systems: exact projection conditioning of a Gaussian solution belief with residual-driven, A-conjugate, and covariance-greedy search policies.
• Auditable numerics: Cholesky solves instead of explicit inverses, jitter as a visible parameter that is never escalated silently, and typed errors instead of NaN for every invalid input.
• Scientific validation shipped as tests: analytic fixtures, independent deterministic quadrature, coverage and calibration studies, misspecification and numerical-stability regressions.
• Small footprint: one dependency (nalgebra), unsafe forbidden, MSRV 1.85.

Installation

cargo add uncertain-numerics

or, in Cargo.toml:

[dependencies]
uncertain-numerics = "0.1"

Quick start

All three snippets below are compiled and run as doctests, and each has a runnable counterpart in examples/.

Bayesian quadrature

Infer $I = \int f(x),p(x),dx$ from a handful of function evaluations and get the posterior uncertainty alongside the estimate.

use uncertain_numerics::{BayesianQuadrature, GaussianMeasure, RbfKernel};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Prior over the integrand: zero-mean Gaussian process with an RBF kernel.
    let kernel = RbfKernel::new(1.0, 1.0)?; // signal variance, length scale
    // Integration measure p(x) = N(0, 1).
    let measure = GaussianMeasure::new(0.0, 1.0)?;
    // Jitter is an explicit, fixed diagonal regularizer. It is never escalated silently.
    let quadrature = BayesianQuadrature::new(kernel, measure, 1.0e-10);

    // Observe f(x) = cos(x) at seven nodes. Exactly, E[cos X] = exp(-1/2) for X ~ N(0, 1).
    let nodes = [-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0];
    let values: Vec<f64> = nodes.iter().copied().map(f64::cos).collect();
    let exact = (-0.5_f64).exp();

    let posterior = quadrature.posterior(&nodes, &values)?;
    println!("E[I | y]  = {:.6}   (exact {exact:.6})", posterior.mean());
    println!("sd[I | y] = {:.3e}", posterior.standard_deviation());

    // The posterior is honest about its own error here: the exact value lies well
    // inside the reported uncertainty.
    assert!((posterior.mean() - exact).abs() < 3.0 * posterior.standard_deviation());
    Ok(())
}

Active Bayesian quadrature

Let the acquisition rule choose where to evaluate next. Every step records the selected point, the predicted variance reduction, and the resulting posterior variance, so stopping decisions are auditable.

use uncertain_numerics::{
    ActiveBayesianQuadrature, BayesianQuadrature, GaussianMeasure, RbfKernel,
};

// A Gaussian bump centred at 0.5. Against p(x) = N(0, 1) its integral has a closed form.
fn integrand(x: f64) -> f64 {
    let centered = x - 0.5;
    (-centered * centered).exp()
}

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let kernel = RbfKernel::new(1.0, 1.0)?;
    let measure = GaussianMeasure::new(0.0, 1.0)?;
    let quadrature = BayesianQuadrature::new(kernel, measure, 1.0e-10);
    let active = ActiveBayesianQuadrature::new(quadrature);

    // Start from three evaluations; allow up to six more from a fixed candidate grid,
    // stopping early once the posterior variance of the integral drops below 1e-6.
    let initial_nodes = [-1.0, 0.0, 1.0];
    let initial_values: Vec<f64> = initial_nodes.iter().copied().map(integrand).collect();
    let candidates: Vec<f64> = (0..=23).map(|i| -2.875 + 0.25 * f64::from(i)).collect();

    let result = active.run(
        &initial_nodes,
        &initial_values,
        &candidates,
        6,
        1.0e-6,
        integrand,
    )?;

    for step in result.steps() {
        println!(
            "x = {:+.3}  predicted reduction = {:.3e}  posterior variance = {:.3e}",
            step.point(),
            step.predicted_variance_reduction(),
            step.posterior_variance(),
        );
    }

    let exact = (1.0_f64 / 3.0).sqrt() * (-0.25_f64 / 3.0).exp();
    let posterior = result.posterior();
    println!("stopped because: {:?}", result.termination());
    println!(
        "E[I | y] = {:.6} ± {:.3e}   (exact {exact:.6})",
        posterior.mean(),
        posterior.standard_deviation(),
    );
    Ok(())
}

Probabilistic linear solvers

Represent the unknown solution of $A x = b$ as a Gaussian belief and condition it on exact projections $s^{\top A x = s}\top b$. The posterior covariance trace is a diagnostic of the uncertainty that remains, distinct from the residual norm.

use uncertain_numerics::{
    CovarianceGreedyProjectionSolver, GaussianLinearBelief, ResidualProjectionSolver,
    SpdLinearSystem,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // A x = b with A symmetric positive definite, stored row-major.
    let system = SpdLinearSystem::new(
        &[4.0, 1.0, 0.0, 1.0, 3.0, 1.0, 0.0, 1.0, 2.0],
        &[1.0, 2.0, 3.0],
        3,
    )?;
    // Prior belief x ~ N(0, I).
    let identity = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
    let prior = GaussianLinearBelief::new(&[0.0; 3], &identity, 3)?;

    // Residual-driven policy: each step observes the projection along the normalized residual.
    let solver = ResidualProjectionSolver::new(1.0e-10, 0.0, 3)?;
    let result = solver.solve(&system, &prior)?;
    println!("residual policy: mean = {:?}", result.belief().mean());
    println!("stopped because: {:?}", result.termination());
    for (index, step) in result.steps().iter().enumerate() {
        println!(
            "  step {index}: residual {:.3e} -> {:.3e}, covariance trace {:.3e}",
            step.residual_norm_before(),
            step.residual_norm_after(),
            step.covariance_trace_after(),
        );
    }

    // Covariance-greedy policy: directions are chosen without looking at b, so the
    // posterior covariance keeps its calibrated interpretation under the assumed prior.
    let candidates: Vec<Vec<f64>> = (0..3)
        .map(|axis| {
            let mut direction = vec![0.0; 3];
            direction[axis] = 1.0;
            direction
        })
        .collect();
    let greedy = CovarianceGreedyProjectionSolver::new(0.0, 2)?;
    let result = greedy.solve(&system, &prior, &candidates)?;
    println!("covariance-greedy: mean = {:?}", result.belief().mean());
    println!("stopped because: {:?}", result.termination());
    Ok(())
}

Mathematical background

Bayesian quadrature posterior

Place a zero-mean Gaussian-process prior with covariance kernel $k$ on the integrand $f$. For observation nodes $X$ with values $y = f(X)$, the integral $I = \int f(x),p(x),dx$ is itself Gaussian a posteriori, with

$$ \mathbb{E}[I \mid y] = z^{\top (K + \lambda I)}{-1} y, \qquad \operatorname{Var}(I \mid y) = \kappa - z^{\top (K + \lambda I)}{-1} z, $$

where $K$ is the Gram matrix, $\lambda$ is the explicit jitter, $z_i = \int k(x_i, x),p(x),dx$ is the kernel mean, and $\kappa = \iint k(x, x'),p(x),p(x'),dx,dx'$ is the prior integral variance. For the RBF kernel and a Gaussian measure both $z$ and $\kappa$ are available in closed form and are tested against independent deterministic quadrature. The inverse is never formed; both solves reuse one Cholesky factorization.

Active-design criterion

For a candidate node $x_*$ the expected reduction in posterior integral variance is

$$ \Delta(x_) = \frac{\left(z_ - k_^{\top A}{-1} z\right)^2}{k(x_, x_) + \lambda - k_^{\top A}{-1} k_*}, \qquad A = K + \lambda I, $$

which does not depend on observed function values. The sequential design selects the first maximum over a finite candidate set, reusing a single factorization across all candidates, and stops on a variance tolerance, an evaluation budget, or candidate exhaustion.

Projection conditioning

For an SPD system $A x = b$ with belief $x \sim \mathcal{N}(m, \Sigma)$, an exact observation along direction $s$ with $h = A^\top s$ gives the rank-one update

$$ m^{+ = m + \frac{\Sigma h}{h}\top \Sigma h}\left(s^{\top b - h}\top m\right), \qquad \Sigma^{+ = \Sigma - \frac{(\Sigma h)(\Sigma h)}\top}{h^\top \Sigma h}. $$

Three policies choose $s$: the normalized residual, the residual orthogonalized in the $A$-inner product against previous directions, and a covariance-greedy rule that maximizes the exact trace reduction $h^{\top \Sigma}2 h / h^\top \Sigma h$ over a candidate set without using $b$. The first two produce equivalent posteriors at equal informative budgets; only the third keeps a calibrated interpretation under repeated prior draws, because its selections do not depend on the unknown solution. See the documents linked below.

Validation

The integration tests are scientific studies rather than smoke tests. Each documents a claim the library makes about its own behavior.

Run everything with:

cargo test --all-features

Documentation

• API reference: docs.rs/uncertain-numerics (or cargo doc --open).
• Bayesian quadrature under misspecification
• Calibration of probabilistic linear-solver uncertainty
• Linear-policy equivalence under exact Gaussian conditioning
• Covariance-greedy probabilistic linear solver
• Roadmap and changelog

Design principles

• Expose uncertainty over the computational quantity, not only a numerical estimate.
• Keep statistical and numerical assumptions explicit and named.
• Validate algorithms against analytic reference problems wherever possible.
• Treat calibration and numerical stability as part of correctness.
• Keep the core independent of dataframe and machine-learning ecosystems.
• Add dependencies only after the numerical abstractions are stable.

Status and scope

The crate is pre-1.0. The public API is deliberately small while the mathematical contracts are established, and minor versions may contain breaking changes that are called out in the changelog.

Current limitations worth knowing before you depend on it:

• Bayesian quadrature supports one dimension, the RBF kernel, a Gaussian measure, and a zero prior mean.
• Linear solvers operate on dense SPD systems and use dense covariance matrices, so they are meant for moderate dimensions.
• Posterior uncertainty from the adaptive residual policies is conservative under repeated prior draws; use the covariance-greedy policy when calibration matters.

The project will not become a general machine-learning framework, Gaussian-process package, or replacement for classical numerical-integration crates.

Minimum supported Rust version

The MSRV is Rust 1.85 (edition 2024). It is checked in CI, and raising it is treated as at least a minor version bump.

Contributing

Contributions are welcome. Please read CONTRIBUTING.md for the development workflow and the testing standards expected for numerical claims. Security issues should follow SECURITY.md.

License

Licensed under either of

• Apache License, Version 2.0 (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or http://opensource.org/licenses/MIT)

at your option.

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.