rusolver/
svd.rs

1// SPDX-License-Identifier: Apache-2.0
2//! Thin real SVD by scaled one-sided cyclic Jacobi rotations. No A^T A is formed.
3use crate::{Matrix, MatrixView, SolverError, Tolerance};
4use crate::numerics::{checked, norm_iter, sum_iter, zeros};
5
6#[derive(Clone, Copy, Debug)]
7pub struct SvdOptions {
8    /// Maximum absolute cosine between non-negligible columns.
9    pub orthogonality_tolerance: f64,
10    pub max_sweeps: usize,
11    /// Singular values at/below this threshold are explicitly truncated to zero.
12    /// The returned factorization is then a rank-truncated approximation.
13    pub rank_tolerance: Tolerance,
14}
15impl Default for SvdOptions {
16    fn default() -> Self { Self { orthogonality_tolerance: 1e-12, max_sweeps: 100,
17        rank_tolerance: Tolerance { absolute: 0.0, relative: 1e-14 } } }
18}
19#[derive(Clone, Debug)]
20pub struct Svd {
21    /// m by min(m,n), including an orthonormal completion of the null columns.
22    pub u: Matrix,
23    /// Nonnegative, descending, explicitly thresholded singular values.
24    pub singular_values: Vec<f64>,
25    /// min(m,n) by n; transpose, not conjugate transpose (real inputs only).
26    pub vt: Matrix,
27    pub rank: usize,
28    pub sweeps: usize,
29    pub max_column_correlation: f64,
30}
31impl Svd {
32    pub fn factor(a: MatrixView<'_>, options: SvdOptions) -> Result<Self, SolverError> {
33        options.rank_tolerance.validate()?;
34        if !options.orthogonality_tolerance.is_finite() || options.orthogonality_tolerance <= 0.0
35            || options.orthogonality_tolerance >= 1.0 || options.max_sweeps == 0 {
36            return Err(SolverError::InvalidOption("SVD sweep budget/tolerance"));
37        }
38        if a.rows() < a.columns() {
39            let tall = Self::factor(a.transpose(), options)?;
40            return Ok(Self { u: tall.vt.transpose()?, vt: tall.u.transpose()?,
41                singular_values: tall.singular_values, rank: tall.rank, sweeps: tall.sweeps,
42                max_column_correlation: tall.max_column_correlation });
43        }
44        let (m,n) = (a.rows(), a.columns());
45        let scale = a.max_abs();
46        let mut b = a.to_owned()?;
47        if scale > 0.0 { for x in &mut b.data { *x /= scale; } }
48        let mut v = Matrix::identity(n)?;
49        let norm = norm_iter(b.data.iter().copied())?;
50        let cutoff = if scale == 0.0 { 0.0 } else {
51            (options.rank_tolerance.absolute/scale).max(options.rank_tolerance.relative*norm)
52        };
53        let mut sweeps = 0;
54        let correlation;
55        loop {
56            let mut largest = 0.0f64;
57            for p in 0..n { for q in p+1..n {
58                let np = norm_iter((0..m).map(|i| b.data[i*n+p]))?;
59                let nq = norm_iter((0..m).map(|i| b.data[i*n+q]))?;
60                if np <= cutoff || nq <= cutoff { continue; }
61                let rho = sum_iter((0..m).map(|i| (b.data[i*n+p]/np)*(b.data[i*n+q]/nq)))?;
62                largest = largest.max(rho.abs());
63                if rho.abs() <= options.orthogonality_tolerance || sweeps == options.max_sweeps { continue; }
64                // Scale the 2x2 Gram block locally. Unlike A^T A, this is never
65                // stored as a global squared-condition-number eigenproblem.
66                let s = np.max(nq);
67                let alpha = (np/s)*(np/s);
68                let beta = (nq/s)*(nq/s);
69                let gamma = rho*(np/s)*(nq/s);
70                let delta = 0.5*(beta-alpha);
71                let t = if delta == 0.0 { gamma.signum() } else {
72                    gamma/(delta+delta.hypot(gamma).copysign(delta))
73                };
74                let c = 1.0/(1.0+t*t).sqrt();
75                let sn = t*c;
76                for i in 0..m {
77                    let (x,y)=(b.data[i*n+p],b.data[i*n+q]);
78                    b.data[i*n+p]=c*x-sn*y; b.data[i*n+q]=sn*x+c*y;
79                }
80                for i in 0..n {
81                    let (x,y)=(v.data[i*n+p],v.data[i*n+q]);
82                    v.data[i*n+p]=c*x-sn*y; v.data[i*n+q]=sn*x+c*y;
83                }
84            }}
85            if largest <= options.orthogonality_tolerance { correlation=largest; break; }
86            if sweeps >= options.max_sweeps {
87                return Err(SolverError::NonConvergence { iterations: sweeps, residual: largest });
88            }
89            sweeps += 1;
90        }
91        let mut norms=zeros(n)?;
92        for j in 0..n { norms[j]=norm_iter((0..m).map(|i| b.data[i*n+j]))?; }
93        let mut order: Vec<usize>=(0..n).collect();
94        order.sort_by(|&p,&q| norms[q].total_cmp(&norms[p]));
95        let mut u=Matrix::zeros(m,n)?;
96        let mut vt=Matrix::zeros(n,n)?;
97        let mut singular_values=zeros(n)?;
98        let mut rank=0;
99        for (j,&old) in order.iter().enumerate() {
100            for i in 0..n { vt.data[j*n+i]=v.data[i*n+old]; }
101            if norms[old] > cutoff && scale > 0.0 {
102                singular_values[j]=checked(norms[old]*scale,"SVD rescale")?;
103                for i in 0..m { u.data[i*n+j]=b.data[i*n+old]/norms[old]; }
104                rank+=1;
105            } else { complete_column(&mut u,j)?; }
106        }
107        Ok(Self {u,vt,singular_values,rank,sweeps,max_column_correlation:correlation})
108    }
109    /// Moore-Penrose minimum-norm solve using the explicitly truncated SVD.
110    /// Supports over/underdetermined and rank-deficient problems, multiple RHS.
111    pub fn solve(&self, b: MatrixView<'_>) -> Result<Matrix,SolverError> {
112        if b.rows()!=self.u.rows { return Err(SolverError::Shape("SVD RHS rows")); }
113        let (m,n,k,r)=(self.u.rows,self.vt.cols,self.singular_values.len(),b.columns());
114        let mut temp=Matrix::zeros(k,r)?;
115        for i in 0..k { if self.singular_values[i]>0.0 { for c in 0..r {
116            temp.data[i*r+c]=checked(sum_iter((0..m).map(|j|self.u.data[j*k+i]*b.at(j,c)))?
117                /self.singular_values[i],"SVD solve projection")?;
118        }}}
119        let mut x=Matrix::zeros(n,r)?;
120        for i in 0..n { for c in 0..r {
121            x.data[i*r+c]=sum_iter((0..k).map(|j|self.vt.data[j*n+i]*temp.data[j*r+c]))?;
122        }}
123        Ok(x)
124    }
125    pub fn pseudo_inverse(&self)->Result<Matrix,SolverError> {
126        self.solve(Matrix::identity(self.u.rows)?.view())
127    }
128}
129fn complete_column(q: &mut Matrix, column: usize)->Result<(),SolverError> {
130    let (m,n)=(q.rows,q.cols);
131    let mut best=zeros(m)?; let mut best_norm=0.0f64;
132    for axis in 0..m {
133        let mut z=zeros(m)?; z[axis]=1.0;
134        for _ in 0..2 { for j in 0..column {
135            let d=sum_iter((0..m).map(|i|q.data[i*n+j]*z[i]))?;
136            for i in 0..m { z[i]-=d*q.data[i*n+j]; }
137        }}
138        let norm=norm_iter(z.iter().copied())?;
139        if norm>best_norm {best_norm=norm; best=z;}
140    }
141    if best_norm<=64.0*f64::EPSILON {return Err(SolverError::Breakdown("SVD null-space completion"));}
142    for i in 0..m {q.data[i*n+column]=best[i]/best_norm;}
143    Ok(())
144}
rusolver/svd.rs

rusolver/
svd.rs
