ruvector-solver

Sublinear-time sparse solvers -- O(log n) PageRank, spectral methods, and linear systems in Rust and WASM.

Most numerical libraries use dense solvers that slow down dramatically as data grows. ruvector-solver provides seven specialized sparse algorithms that run in O(log n) to O(sqrt(n)) time, automatically picks the best one for your problem, and works in native Rust or in the browser via WebAssembly. It powers the graph analytics and AI coherence layers inside RuVector.

Quick Start

[dependencies]
ruvector-solver = "0.1"

use ruvector_solver::types::CsrMatrix;
use ruvector_solver::neumann::NeumannSolver;

// Build a sparse 3x3 system and solve it
let a = CsrMatrix::<f32>::from_coo(3, 3, vec![
    (0, 0, 2.0_f32), (0, 1, -0.5),
    (1, 0, -0.5),    (1, 1, 2.0),  (1, 2, -0.5),
    (2, 1, -0.5),    (2, 2, 2.0),
]);
let b = vec![1.0_f32, 0.0, 1.0];

let solver = NeumannSolver::new(1e-6, 500);
let result = solver.solve(&a, &b).unwrap();

println!("solution: {:?}", result.solution);
println!("iterations: {}", result.iterations);
println!("residual:   {:.2e}", result.residual_norm);

Or let the router pick the best algorithm automatically:

use ruvector_solver::router::{SolverRouter, QueryType};
use ruvector_solver::types::{CsrMatrix, ComputeBudget};

let router = SolverRouter::new();
let (algo, result) = router.solve(&matrix, &rhs, &ComputeBudget::default(), QueryType::LinearSystem).unwrap();
println!("Router selected: {:?}", algo);

Key Features

Algorithms

Feature Flags

Enable all algorithms:

[dependencies]
ruvector-solver = { version = "0.1", features = ["all-algorithms"] }

Performance Optimisations

Bounds-check-free SpMV (spmv_unchecked)

The inner SpMV loop is the single hottest path in every iterative solver. The spmv_unchecked method on CsrMatrix<f32> and CsrMatrix<f64> uses raw pointers to eliminate per-element bounds checks, relying on a one-time CSR structural validation (validation::validate_csr_matrix) performed before entering the solve loop.

Fused residual + norm computation (fused_residual_norm_sq)

Standard implementations compute the residual r = b - Ax and its squared norm ||r||^2 in separate passes (SpMV, vector subtraction, dot product -- three full memory traversals). fused_residual_norm_sq computes both in a single pass, reducing memory traffic by roughly 3x per iteration.

AVX2 8-wide SIMD SpMV

When the simd feature is enabled on x86_64, spmv_simd dispatches to an AVX2 kernel that processes 8 f32 values per instruction using _mm256 intrinsics with a horizontal sum reduction at the end of each row. Falls back to a portable scalar loop on other architectures.

4-wide unrolled Jacobi update

The Neumann iteration's update step x[j] += d_inv[j] * r[j] is manually unrolled 4-wide for instruction-level parallelism, with a scalar remainder loop for dimensions not divisible by 4.

Arena allocator

SolverArena provides bump allocation for per-solve scratch buffers. All temporary vectors are allocated from a single contiguous backing buffer and reclaimed in O(1) via arena.reset(), eliminating heap allocation overhead inside the iteration loop.

Architecture

                          +-------------------+
                          |   SolverRouter    |
                          | (algorithm select)|
                          +--------+----------+
                                   |
            +----------+-----------+-----------+----------+
            |          |           |           |          |
       +----v---+ +----v---+ +----v------+ +--v----+ +---v----+
       |Neumann | |   CG   | |ForwardPush| | TRUE  | | BMSSP  |
       |Solver  | | Solver | |  Solver   | |Solver | |Solver  |
       +----+---+ +----+---+ +-----+-----+ +--+----+ +---+----+
            |          |            |          |          |
            +-----+----+-----+-----+-----+----+-----+---+
                  |          |           |           |
             +----v---+ +---v----+ +----v----+ +----v-----+
             |types.rs| |simd.rs | |arena.rs | |budget.rs |
             |CsrMatrix| |AVX2   | |Bump     | |ComputeBudget|
             +--------+ |SpMV   | |Alloc    | |enforcement|
                         +-------+ +--------+ +----------+
                  |          |           |
             +----v---+ +---v------+ +--v---------+
             |traits.rs| |validate.rs| |error.rs    |
             |SolverEngine| |CSR check| |SolverError |
             +--------+ +---------+ +-----------+

The SolverRouter analyses the matrix SparsityProfile and QueryType to select the optimal algorithm. When the selected algorithm fails, SolverOrchestrator::solve_with_fallback tries a deterministic fallback chain: selected -> CG -> Dense.

API Overview

Core types (types.rs)

Traits (traits.rs)

Error hierarchy (error.rs)

Infrastructure modules

Testing

The crate includes 177 tests (138 unit tests + 39 integration/doctests):

# Run all tests
cargo test -p ruvector-solver

# Run tests with all algorithms enabled
cargo test -p ruvector-solver --features all-algorithms

# Run a specific test module
cargo test -p ruvector-solver -- neumann::tests

Benchmarks

Five Criterion benchmark groups are provided:

# Run all benchmarks
cargo bench -p ruvector-solver

# Run a specific benchmark
cargo bench -p ruvector-solver --bench solver_neumann

use ruvector_solver::types::CsrMatrix;

// 4x4 tridiagonal system (diagonally dominant)
let a = CsrMatrix::<f32>::from_coo(4, 4, vec![
    (0, 0, 3.0), (0, 1, -1.0),
    (1, 0, -1.0), (1, 1, 3.0), (1, 2, -1.0),
    (2, 1, -1.0), (2, 2, 3.0), (2, 3, -1.0),
    (3, 2, -1.0), (3, 3, 3.0),
]);
let b = vec![2.0f32, 1.0, 1.0, 2.0];

use ruvector_solver::neumann::NeumannSolver;

let solver = NeumannSolver::new(1e-6, 500);
let result = solver.solve(&a, &b).unwrap();

println!("Solution:   {:?}", result.solution);
println!("Iterations: {}", result.iterations);
println!("Residual:   {:.2e}", result.residual_norm);

use ruvector_solver::router::{SolverRouter, QueryType};
use ruvector_solver::types::{CsrMatrix, ComputeBudget};

let a64 = CsrMatrix::<f64>::from_coo(4, 4, vec![/* ... */]);
let b64 = vec![2.0, 1.0, 1.0, 2.0];
let budget = ComputeBudget::default();

let router = SolverRouter::new();
let (algo, result) = router.solve(&a64, &b64, &budget, QueryType::LinearSystem).unwrap();
println!("Router selected: {:?}", algo);

use ruvector_solver::validation::validate_csr_matrix;

let errors = validate_csr_matrix(&a);
assert!(errors.is_empty(), "CSR validation failed: {:?}", errors);

cargo bench -p ruvector-solver --bench solver_neumann
cargo bench -p ruvector-solver --bench solver_e2e

use ruvector_solver::forward_push::ForwardPushSolver;
use ruvector_solver::types::CsrMatrix;

// Build adjacency matrix for a small graph
let adj = CsrMatrix::<f32>::from_coo(4, 4, vec![
    (0, 1, 1.0), (1, 0, 1.0),
    (1, 2, 1.0), (2, 1, 1.0),
    (2, 3, 1.0), (3, 2, 1.0),
    (0, 3, 1.0), (3, 0, 1.0),
]);

let solver = ForwardPushSolver::new(0.85, 1e-6);  // alpha=0.85
let ppr = solver.ppr(&adj, 0);  // PPR from node 0

println!("PPR scores: {:?}", ppr);

Related Crates

Minimum Supported Rust Version

Rust 1.77 or later.

License

Licensed under the MIT License.

Part of RuVector -- the self-learning vector database.