rusolver/
qr.rs

1// SPDX-License-Identifier: Apache-2.0
2use crate::{
3    Matrix, MatrixView, SolverError, Tolerance
4};
5use crate::numerics::{
6    checked, norm_iter, sum_iter, zeros
7};
8/// Economy QR with recomputed column norms and column pivoting.
9/// Convention: `A[:, permutation] = Q * R`, where Q has m rows and n columns.
10/// Householder reflectors are stored, not a full m*m Q. Only m >= n is supported.
11#[derive(Clone, Debug)]
12pub struct Qr {
13    packed: Matrix, tau: Vec<f64>, permutation: Vec<usize>, rank: usize
14}
15#[derive(Clone, Debug)]
16pub struct LeastSquares {
17    pub solution: Matrix, pub residual_norms: Vec<f64>, pub rank: usize
18}
19impl Qr {
20    pub fn factor(a: MatrixView<'_>, tolerance: Tolerance)->Result<Self, SolverError> {
21        let (m, n)=(a.rows(), a.columns());
22        if m<n {
23            return Err(SolverError::Shape("QR requires rows >= columns"));
24        }
25        tolerance.validate()?;
26        let mut packed=a.to_owned()?;
27        let mut tau=zeros(n)?;
28        let mut permutation: Vec<usize>=(0..n).collect();
29        let mut scale=0.0f64;
30        for j in 0..n {
31            scale=scale.max(norm_iter((0..m).map(|i|a.at(i, j)))?);
32        }
33        let cutoff=tolerance.threshold(scale)?;
34        for k in 0..n {
35            let mut best=k;
36            let mut best_norm=-1.0;
37            for j in k..n {
38                let v=norm_iter((k..m).map(|i|packed.data[i*n+j]))?;
39                if v>best_norm {
40                    best_norm=v;
41                    best=j;
42                }
43            }
44            if best!=k {
45                for i in 0..m {
46                    packed.data.swap(i*n+k, i*n+best);
47                }
48                permutation.swap(k, best);
49            }
50            if best_norm==0.0 {
51                continue;
52            }
53            let x0=packed.data[k*n+k];
54            let sign=if x0>=0.0 {
55                1.0
56            } else {
57                -1.0
58            };
59            let alpha=-sign*best_norm;
60            // Form x/(||x||) before subtracting sign to avoid x0-alpha overflow.
61            let divisor=x0/best_norm+sign;
62            tau[k]=1.0+x0.abs()/best_norm;
63            for i in k+1..m {
64                packed.data[i*n+k]=checked((packed.data[i*n+k]/best_norm)/divisor, "Householder vector")?;
65            }
66            packed.data[k*n+k]=alpha;
67            for j in k+1..n {
68                let w=checked(tau[k]*sum_iter(std::iter::once(packed.data[k*n+j])
69                .chain((k+1..m).map(|i|packed.data[i*n+k]*packed.data[i*n+j])))?, "Householder update")?;
70                packed.data[k*n+j]=checked(packed.data[k*n+j]-w, "QR leading update")?;
71                for i in k+1..m {
72                    packed.data[i*n+j]=checked((-packed.data[i*n+k]).mul_add(w, packed.data[i*n+j]), "QR trailing update")?;
73                }
74            }
75        }
76        // Rank is a tolerance-based diagnostic, not a singular-value estimate.
77        let rank=(0..n).take_while(|&k|packed.data[k*n+k].abs()>cutoff).count();
78        Ok(Self{
79            packed, tau, permutation, rank
80        })
81    }
82    pub fn rank(&self)->usize {
83        self.rank
84    }
85    pub fn permutation(&self)->&[usize] {
86        &self.permutation
87    }
88    pub fn r(&self)->Result<Matrix, SolverError> {
89        let n=self.packed.cols;
90        let mut r=Matrix::zeros(n, n)?;
91        for i in 0..n {
92            for j in i..n {
93                r.data[i*n+j]=self.packed.data[i*n+j];
94            }
95        }
96        Ok(r)
97    }
98    /// Materialize thin Q only on explicit request.
99    pub fn q(&self)->Result<Matrix, SolverError> {
100        let (m, n)=(self.packed.rows, self.packed.cols);
101        let mut q=Matrix::zeros(m, n)?;
102        for i in 0..n {
103            q.data[i*n+i]=1.0;
104        }
105        for k in (0..n).rev() {
106            self.apply_reflector(k, &mut q)?;
107        }
108        Ok(q)
109    }
110    fn apply_reflector(&self, k: usize, b: &mut Matrix)->Result<(), SolverError> {
111        let (m, n, p)=(self.packed.rows, self.packed.cols, b.cols);
112        if self.tau[k]==0.0 {
113            return Ok(());
114        }
115        for j in 0..p {
116            let w=checked(self.tau[k]*sum_iter(std::iter::once(b.data[k*p+j])
117            .chain((k+1..m).map(|i|self.packed.data[i*n+k]*b.data[i*p+j])))?, "Q application")?;
118            b.data[k*p+j]=checked(b.data[k*p+j]-w, "Q application")?;
119            for i in k+1..m {
120                b.data[i*p+j]=checked((-self.packed.data[i*n+k]).mul_add(w, b.data[i*p+j]), "Q application")?;
121            }
122        }
123        Ok(())
124    }
125    /// Full-column-rank least squares without forming A^T A or an inverse.
126    /// Rank-deficient and underdetermined minimum-norm problems are NOT approximated.
127    pub fn least_squares(&self, b: MatrixView<'_>)->Result<LeastSquares, SolverError> {
128        let (m, n, p)=(self.packed.rows, self.packed.cols, b.columns());
129        if b.rows()!=m {
130            return Err(SolverError::Shape("QR RHS row count"));
131        }
132        if self.rank<n {
133            return Err(SolverError::RankDeficient{
134                rank: self.rank, columns: n
135            });
136        }
137        let mut y=b.to_owned()?;
138        for k in 0..n {
139            self.apply_reflector(k, &mut y)?;
140        }
141        let mut residual_norms=zeros(p)?;
142        for j in 0..p {
143            residual_norms[j]=norm_iter((n..m).map(|i|y.data[i*p+j]))?;
144        }
145        for i in (0..n).rev() {
146            for c in 0..p {
147                let sum=sum_iter((i+1..n).map(|j|self.packed.data[i*n+j]*y.data[j*p+c]))?;
148                y.data[i*p+c]=checked((y.data[i*p+c]-sum)/self.packed.data[i*n+i], "QR triangular solve")?;
149            }
150        }
151        let mut x=Matrix::zeros(n, p)?;
152        for j in 0..n {
153            for c in 0..p {
154                x.data[self.permutation[j]*p+c]=y.data[j*p+c];
155            }
156        }
157        Ok(LeastSquares{
158            solution: x, residual_norms, rank: self.rank
159        })
160    }
161}
rusolver/qr.rs

rusolver/
qr.rs
