probabilistic_linear_solver/
probabilistic_linear_solver.rs

1// Probabilistic linear solvers: a Gaussian belief over x conditioned on exact projections.
2//
3// Run with: cargo run --example probabilistic_linear_solver
4
5use uncertain_numerics::{
6    AConjugateProjectionSolver, CovarianceGreedyProjectionSolver, GaussianLinearBelief,
7    ResidualProjectionSolver, SpdLinearSystem,
8};
9
10fn main() -> Result<(), Box<dyn std::error::Error>> {
11    // A x = b with A symmetric positive definite, stored row-major.
12    let system = SpdLinearSystem::new(
13        &[4.0, 1.0, 0.0, 1.0, 3.0, 1.0, 0.0, 1.0, 2.0],
14        &[1.0, 2.0, 3.0],
15        3,
16    )?;
17    // Prior belief x ~ N(0, I).
18    let identity = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
19    let prior = GaussianLinearBelief::new(&[0.0; 3], &identity, 3)?;
20
21    // Residual-driven policy: each step observes the projection along the normalized residual.
22    let solver = ResidualProjectionSolver::new(1.0e-10, 0.0, 3)?;
23    let result = solver.solve(&system, &prior)?;
24    println!("residual policy");
25    println!("  mean            = {:?}", result.belief().mean());
26    println!("  stopped because = {:?}", result.termination());
27    for (index, step) in result.steps().iter().enumerate() {
28        println!(
29            "  step {index}: residual {:.3e} -> {:.3e}, covariance trace {:.3e}",
30            step.residual_norm_before(),
31            step.residual_norm_after(),
32            step.covariance_trace_after(),
33        );
34    }
35
36    // A-conjugate policy: the same information subspace in an A-orthogonal basis.
37    let conjugate = AConjugateProjectionSolver::new(1.0e-10, 0.0, 3)?;
38    let result = conjugate.solve(&system, &prior)?;
39    println!("A-conjugate policy");
40    println!("  mean            = {:?}", result.belief().mean());
41    println!("  stopped because = {:?}", result.termination());
42
43    // Covariance-greedy policy: directions are chosen without looking at b, so the
44    // posterior covariance keeps its calibrated interpretation under the assumed prior.
45    let candidates: Vec<Vec<f64>> = (0..3)
46        .map(|axis| {
47            let mut direction = vec![0.0; 3];
48            direction[axis] = 1.0;
49            direction
50        })
51        .collect();
52    let greedy = CovarianceGreedyProjectionSolver::new(0.0, 2)?;
53    let result = greedy.solve(&system, &prior, &candidates)?;
54    println!("covariance-greedy policy (budget of two projections)");
55    println!("  mean            = {:?}", result.belief().mean());
56    println!("  stopped because = {:?}", result.termination());
57    for (index, step) in result.steps().iter().enumerate() {
58        println!(
59            "  step {index}: direction {:?}, predicted trace reduction {:.3e}, posterior trace {:.3e}",
60            step.direction(),
61            step.predicted_trace_reduction(),
62            step.posterior_trace(),
63        );
64    }
65    Ok(())
66}
probabilistic_linear_solver/probabilistic_linear_solver.rs

probabilistic_linear_solver/
probabilistic_linear_solver.rs
