Struct ResidualProjectionSolver 

pub struct ResidualProjectionSolver { /* private fields */ }

Residual-driven probabilistic solver for dense SPD systems.

At iteration k, the normalized residual of the current posterior mean, s_k = r_k / ||r_k||, defines the next exact projection observation s_k^T A x = s_k^T b.

Implementations

impl ResidualProjectionSolver

pub fn new( residual_tolerance: f64, covariance_trace_tolerance: f64, max_iterations: usize, ) -> Result<Self, LinearSolverError>

Construct the solver with explicit stopping tolerances.

Errors

Returns LinearSolverError when either tolerance is non-finite or negative.

Examples found in repository

examples/probabilistic_linear_solver.rs (line 22)

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");
    println!("  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(),
        );
    }

    // A-conjugate policy: the same information subspace in an A-orthogonal basis.
    let conjugate = AConjugateProjectionSolver::new(1.0e-10, 0.0, 3)?;
    let result = conjugate.solve(&system, &prior)?;
    println!("A-conjugate policy");
    println!("  mean            = {:?}", result.belief().mean());
    println!("  stopped because = {:?}", result.termination());

    // 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 policy (budget of two projections)");
    println!("  mean            = {:?}", result.belief().mean());
    println!("  stopped because = {:?}", result.termination());
    for (index, step) in result.steps().iter().enumerate() {
        println!(
            "  step {index}: direction {:?}, predicted trace reduction {:.3e}, posterior trace {:.3e}",
            step.direction(),
            step.predicted_trace_reduction(),
            step.posterior_trace(),
        );
    }
    Ok(())
}

pub fn solve( &self, system: &SpdLinearSystem, initial_belief: &GaussianLinearBelief, ) -> Result<ProbabilisticLinearSolveResult, LinearSolverError>

Solve from an initial Gaussian belief.

Errors

Returns LinearSolverError when the belief dimension does not match the system or when a conditioning step fails for a reason other than a degenerate residual projection.

Examples found in repository

examples/probabilistic_linear_solver.rs (line 23)

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");
    println!("  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(),
        );
    }

    // A-conjugate policy: the same information subspace in an A-orthogonal basis.
    let conjugate = AConjugateProjectionSolver::new(1.0e-10, 0.0, 3)?;
    let result = conjugate.solve(&system, &prior)?;
    println!("A-conjugate policy");
    println!("  mean            = {:?}", result.belief().mean());
    println!("  stopped because = {:?}", result.termination());

    // 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 policy (budget of two projections)");
    println!("  mean            = {:?}", result.belief().mean());
    println!("  stopped because = {:?}", result.termination());
    for (index, step) in result.steps().iter().enumerate() {
        println!(
            "  step {index}: direction {:?}, predicted trace reduction {:.3e}, posterior trace {:.3e}",
            step.direction(),
            step.predicted_trace_reduction(),
            step.posterior_trace(),
        );
    }
    Ok(())
}

Trait Implementations

impl Clone for ResidualProjectionSolver

fn clone(&self) -> ResidualProjectionSolver

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Copy for ResidualProjectionSolver

impl Debug for ResidualProjectionSolver

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl PartialEq for ResidualProjectionSolver

fn eq(&self, other: &ResidualProjectionSolver) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl StructuralPartialEq for ResidualProjectionSolver