Struct SpdLinearSystem

Source

pub struct SpdLinearSystem { /* private fields */ }

Expand description

Dense symmetric positive-definite linear system A x = b.

Implementations§

Source §

impl SpdLinearSystem

Source

pub fn new( matrix: &[f64], rhs: &[f64], dimension: usize, ) -> Result<Self, LinearSolverError>

Construct and validate an SPD linear system from a row-major matrix.

§Errors

Returns LinearSolverError for invalid dimensions, non-finite values, non-symmetry, or failure of Cholesky positive-definiteness validation.

Examples found in repository ?

examples/probabilistic_linear_solver.rs (lines 12-16)

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}