Null space using nalgebra

null_space.rs

use nalgebra::{DMatrix, DVector};

fn compute_nullspace_qr(a: &DMatrix<f64>, tolerance: f64) -> DMatrix<f64> {
    // Get dimensions
    let (m, n) = a.shape();
    
    // Compute QR decomposition with column pivoting
    let qr = a.clone().qr();
    let r = qr.r();
    
    // Find rank by counting non-zero diagonal elements of R
    let mut rank = 0;
    for i in 0..r.nrows().min(r.ncols()) {
        if r[(i, i)].abs() > tolerance {
            rank += 1;
        }
    }
    
    // If rank equals number of columns, nullspace is empty
    if rank == n {
        return DMatrix::zeros(n, 0);
    }
    
    // For each linearly dependent column, construct a nullspace vector
    let mut nullspace = Vec::new();
    
    for j in rank..n {
        // Create vector that will become a nullspace basis vector
        let mut x = DVector::zeros(n);
        
        // Set the free variable to 1
        x[j] = 1.0;
        
        // Back-substitute to find other components
        for i in (0..rank).rev() {
            let mut sum = 0.0;
            for k in (i + 1)..n {
                sum += r[(i, k)] * x[k];
            }
            x[i] = -sum / r[(i, i)];
        }
        
        nullspace.push(x);
    }
    
    // Combine vectors into matrix
    if nullspace.is_empty() {
        DMatrix::zeros(n, 0)
    } else {
        let first = &nullspace[0];
        let mut result = DMatrix::zeros(n, nullspace.len());
        for (j, vec) in nullspace.iter().enumerate() {
            result.set_column(j, vec);
        }
        result
    }
}

// Example usage and tests
fn main() {
    // Example matrix with known nullspace:
    // This matrix represents the linear dependency: row2 = 2*row1, row3 = 3*row1
    let a = DMatrix::from_row_slice(3, 4, &[
        1.0, 2.0, 3.0, 4.0,
        2.0, 4.0, 6.0, 8.0,
        3.0, 6.0, 9.0, 12.0,
    ]);
    
    let tolerance = 1e-10;
    let nullspace = compute_nullspace_qr(&a, tolerance);
    
    println!("Input matrix:\n{}", a);
    println!("\nNullspace basis vectors:\n{}", nullspace);
    
    // Verify result: A * nullspace should be approximately zero
    let verification = &a * &nullspace;
    println!("\nVerification (should be close to zero):\n{}", verification);
    
    // Verify orthogonality of nullspace vectors
    if nullspace.ncols() > 1 {
        let orthogonality = nullspace.transpose() * &nullspace;
        println!("\nOrthogonality check (should be diagonal-ish):\n{}", orthogonality);
    }
}