tensor4all_tcicore/
matrix_luci.rs

1//! Matrix LU-based Cross Interpolation (MatrixLUCI) implementation.
2//!
3//! [`MatrixLUCI`] provides a higher-level row-major API over the lower-level
4//! `matrixluci` substrate. It decomposes a matrix into left and right factors
5//! via LU cross interpolation and implements [`AbstractMatrixCI`].
6
7use crate::error::{MatrixCIError, Result};
8use crate::matrix::{submatrix, zeros, Matrix};
9use crate::matrixlu::RrLUOptions;
10use crate::matrixluci::{
11    CrossFactors, DenseFaerLuKernel, DenseMatrixSource, PivotKernel, PivotKernelOptions,
12    PivotSelectionCore,
13};
14use crate::scalar::Scalar;
15use crate::traits::AbstractMatrixCI;
16
17/// Matrix LU-based Cross Interpolation.
18///
19/// This is a higher-level row-major wrapper around the lower-level `matrixluci`
20/// substrate.
21///
22/// # Examples
23///
24/// ```
25/// use tensor4all_tcicore::{AbstractMatrixCI, MatrixLUCI, from_vec2d};
26///
27/// let m = from_vec2d(vec![
28///     vec![1.0_f64, 2.0, 3.0],
29///     vec![4.0, 5.0, 6.0],
30///     vec![7.0, 8.0, 9.0],
31/// ]);
32///
33/// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
34/// // The approximation must have at most rank min(nrows, ncols)
35/// assert!(ci.rank() <= 3);
36/// // Reconstructed matrix should match original at pivot positions
37/// let row_indices = ci.row_indices().to_vec();
38/// let col_indices = ci.col_indices().to_vec();
39/// for (&i, &j) in row_indices.iter().zip(col_indices.iter()) {
40///     let approx = ci.evaluate(i, j);
41///     let exact = m[[i, j]];
42///     assert!((approx - exact).abs() < 1e-10);
43/// }
44/// ```
45#[derive(Debug, Clone)]
46pub struct MatrixLUCI<T: Scalar + crate::matrixluci::Scalar> {
47    nrows: usize,
48    ncols: usize,
49    row_indices: Vec<usize>,
50    col_indices: Vec<usize>,
51    left: Matrix<T>,
52    right: Matrix<T>,
53    pivot_errors: Vec<f64>,
54}
55
56pub(crate) fn map_backend_error(err: crate::matrixluci::MatrixLuciError) -> MatrixCIError {
57    match err {
58        crate::matrixluci::MatrixLuciError::InvalidArgument { message } => {
59            MatrixCIError::InvalidArgument { message }
60        }
61        crate::matrixluci::MatrixLuciError::SingularPivotBlock => MatrixCIError::SingularMatrix,
62    }
63}
64
65pub(crate) fn to_column_major<T: Scalar>(matrix: &Matrix<T>) -> Vec<T> {
66    let mut out = Vec::with_capacity(matrix.nrows() * matrix.ncols());
67    for col in 0..matrix.ncols() {
68        for row in 0..matrix.nrows() {
69            out.push(matrix[[row, col]]);
70        }
71    }
72    out
73}
74
75pub(crate) fn to_row_major<T: Scalar + crate::matrixluci::Scalar>(
76    matrix: &crate::matrixluci::DenseOwnedMatrix<T>,
77) -> Matrix<T> {
78    let mut out = zeros(matrix.nrows(), matrix.ncols());
79    for col in 0..matrix.ncols() {
80        for row in 0..matrix.nrows() {
81            out[[row, col]] = matrix[[row, col]];
82        }
83    }
84    out
85}
86
87pub(crate) fn dense_selection_from_matrix<T>(
88    a: &Matrix<T>,
89    options: RrLUOptions,
90) -> Result<(PivotSelectionCore, CrossFactors<T>)>
91where
92    T: Scalar + crate::matrixluci::Scalar,
93    DenseFaerLuKernel: PivotKernel<T>,
94{
95    let data = to_column_major(a);
96    let source = DenseMatrixSource::from_column_major(&data, a.nrows(), a.ncols());
97    let kernel_options = PivotKernelOptions {
98        max_rank: options.max_rank,
99        rel_tol: options.rel_tol,
100        abs_tol: options.abs_tol,
101        left_orthogonal: options.left_orthogonal,
102    };
103
104    let selection = DenseFaerLuKernel
105        .factorize(&source, &kernel_options)
106        .map_err(map_backend_error)?;
107    let factors = CrossFactors::from_source(&source, &selection).map_err(map_backend_error)?;
108    Ok((selection, factors))
109}
110
111impl<T> MatrixLUCI<T>
112where
113    T: Scalar + crate::matrixluci::Scalar,
114    DenseFaerLuKernel: PivotKernel<T>,
115{
116    /// Create a MatrixLUCI from a dense row-major matrix.
117    ///
118    /// # Examples
119    ///
120    /// ```
121    /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixLUCI, from_vec2d};
122    ///
123    /// let m = from_vec2d(vec![
124    ///     vec![2.0_f64, 0.0],
125    ///     vec![0.0, 3.0],
126    /// ]);
127    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
128    /// assert!(ci.rank() >= 1);
129    /// ```
130    pub fn from_matrix(a: &Matrix<T>, options: Option<RrLUOptions>) -> Result<Self> {
131        let options = options.unwrap_or_default();
132        let left_orthogonal = options.left_orthogonal;
133        let (selection, factors) = dense_selection_from_matrix(a, options)?;
134
135        let left = if left_orthogonal {
136            to_row_major(&factors.cols_times_pivot_inv().map_err(map_backend_error)?)
137        } else {
138            to_row_major(&factors.pivot_cols)
139        };
140        let right = if left_orthogonal {
141            to_row_major(&factors.pivot_rows)
142        } else {
143            to_row_major(&factors.pivot_inv_times_rows().map_err(map_backend_error)?)
144        };
145
146        Ok(Self {
147            nrows: a.nrows(),
148            ncols: a.ncols(),
149            row_indices: selection.row_indices,
150            col_indices: selection.col_indices,
151            left,
152            right,
153            pivot_errors: selection.pivot_errors,
154        })
155    }
156
157    /// Left CI factor (shape: `nrows x rank`).
158    ///
159    /// The approximation is `left * right`.
160    ///
161    /// # Examples
162    ///
163    /// ```
164    /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixLUCI, from_vec2d, matrix::mat_mul};
165    ///
166    /// let m = from_vec2d(vec![
167    ///     vec![1.0_f64, 2.0],
168    ///     vec![3.0, 4.0],
169    /// ]);
170    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
171    /// let reconstructed = mat_mul(&ci.left(), &ci.right());
172    /// for i in 0..2 {
173    ///     for j in 0..2 {
174    ///         assert!((reconstructed[[i, j]] - m[[i, j]]).abs() < 1e-10);
175    ///     }
176    /// }
177    /// ```
178    pub fn left(&self) -> Matrix<T> {
179        self.left.clone()
180    }
181
182    /// Right CI factor (shape: `rank x ncols`).
183    ///
184    /// The approximation is `left * right`.
185    ///
186    /// # Examples
187    ///
188    /// ```
189    /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixLUCI, from_vec2d, matrix::mat_mul};
190    ///
191    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
192    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
193    /// let r = ci.right();
194    /// assert_eq!(r.nrows(), ci.rank());
195    /// assert_eq!(r.ncols(), ci.ncols());
196    /// // left * right reconstructs the matrix
197    /// let recon = mat_mul(&ci.left(), &r);
198    /// for i in 0..2 {
199    ///     for j in 0..2 {
200    ///         assert!((recon[[i, j]] - m[[i, j]]).abs() < 1e-10);
201    ///     }
202    /// }
203    /// ```
204    pub fn right(&self) -> Matrix<T> {
205        self.right.clone()
206    }
207
208    /// Pivot error history (one entry per pivot, plus a final residual estimate).
209    ///
210    /// # Examples
211    ///
212    /// ```
213    /// use tensor4all_tcicore::{MatrixLUCI, from_vec2d};
214    ///
215    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
216    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
217    /// let errs = ci.pivot_errors();
218    /// assert!(!errs.is_empty());
219    /// // All errors are non-negative
220    /// for &e in &errs {
221    ///     assert!(e >= 0.0);
222    /// }
223    /// ```
224    pub fn pivot_errors(&self) -> Vec<f64> {
225        self.pivot_errors.clone()
226    }
227
228    /// Last pivot error (the residual estimate after all pivots).
229    ///
230    /// # Examples
231    ///
232    /// ```
233    /// use tensor4all_tcicore::{MatrixLUCI, from_vec2d};
234    ///
235    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
236    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
237    /// let err = ci.last_pivot_error();
238    /// assert!(err >= 0.0);
239    /// ```
240    pub fn last_pivot_error(&self) -> f64 {
241        self.pivot_errors.last().copied().unwrap_or(0.0)
242    }
243}
244
245impl<T> AbstractMatrixCI<T> for MatrixLUCI<T>
246where
247    T: Scalar + crate::matrixluci::Scalar,
248{
249    fn nrows(&self) -> usize {
250        self.nrows
251    }
252
253    fn ncols(&self) -> usize {
254        self.ncols
255    }
256
257    fn rank(&self) -> usize {
258        self.row_indices.len()
259    }
260
261    fn row_indices(&self) -> &[usize] {
262        &self.row_indices
263    }
264
265    fn col_indices(&self) -> &[usize] {
266        &self.col_indices
267    }
268
269    fn evaluate(&self, i: usize, j: usize) -> T {
270        let mut sum = T::zero();
271        for k in 0..self.rank() {
272            sum = sum + self.left[[i, k]] * self.right[[k, j]];
273        }
274        sum
275    }
276
277    fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
278        let r = self.rank();
279        let left_sub = submatrix(&self.left, rows, &(0..r).collect::<Vec<_>>());
280        let right_sub = submatrix(&self.right, &(0..r).collect::<Vec<_>>(), cols);
281
282        crate::matrix::mat_mul(&left_sub, &right_sub)
283    }
284}
285
286#[cfg(test)]
287mod tests;
tensor4all_tcicore/matrix_luci.rs

tensor4all_tcicore/
matrix_luci.rs
