tensor4all_core/
matrix_luci.rs

1//! Matrix LU-based Cross Interpolation (MatrixLUCI) implementation.
2//!
3//! [`MatrixLUCI`] provides a higher-level [`Matrix`] 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::matrixlu::{rrlu, rrlu_mut, RrLU, RrLUOptions};
9use crate::matrixluci::block_rook::LazyBlockRookKernel;
10use crate::matrixluci::factors::CrossFactors;
11use crate::matrixluci::source::LazyMatrixSource;
12use crate::matrixluci::types::{PivotKernelOptions, PivotSelectionCore};
13use crate::matrixluci::PivotKernel;
14use crate::scalar::Scalar;
15use crate::traits::AbstractMatrixCI;
16use tensor4all_tensorbackend::{mat_mul_owned, submatrix, triangular_solve_matrix_owned, Matrix};
17
18/// Matrix LU-based Cross Interpolation.
19///
20/// This is a higher-level [`Matrix`] wrapper around the lower-level `matrixluci`
21/// substrate.
22///
23/// # Examples
24///
25/// ```
26/// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
27/// use tensor4all_tensorbackend::from_vec2d;
28///
29/// let m = from_vec2d(vec![
30///     vec![1.0_f64, 2.0, 3.0],
31///     vec![4.0, 5.0, 6.0],
32///     vec![7.0, 8.0, 9.0],
33/// ]);
34///
35/// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
36/// // The approximation must have at most rank min(nrows, ncols)
37/// assert!(ci.rank() <= 3);
38/// // Reconstructed matrix should match original at pivot positions
39/// let row_indices = ci.row_indices().to_vec();
40/// let col_indices = ci.col_indices().to_vec();
41/// for (&i, &j) in row_indices.iter().zip(col_indices.iter()) {
42///     let approx = ci.evaluate(i, j);
43///     let exact = m[[i, j]];
44///     assert!((approx - exact).abs() < 1e-10);
45/// }
46/// ```
47#[derive(Debug, Clone)]
48pub struct MatrixLUCI<T: Scalar + crate::MatrixLuciScalar> {
49    nrows: usize,
50    ncols: usize,
51    row_indices: Vec<usize>,
52    col_indices: Vec<usize>,
53    left: Matrix<T>,
54    right: Matrix<T>,
55    pivot_errors: Vec<f64>,
56}
57
58/// High-level factors produced by MatrixLUCI.
59///
60/// This is the public result of the LUCI factorization facade. It exposes
61/// the selected pivot metadata and the left and right factors
62/// needed by higher-level tensor-network code without exposing the low-level
63/// pivot kernel substrate.
64///
65/// Related types: [`MatrixLUCI`] is the owning CI wrapper, while
66/// [`MatrixACA`](crate::MatrixACA) and [`RrLU`](crate::RrLU) are related
67/// matrix factorization entry points.
68///
69/// # Examples
70///
71/// ```
72/// use tensor4all_core::matrix_luci_factors_from_matrix;
73/// use tensor4all_tensorbackend::from_vec2d;
74///
75/// let m = from_vec2d(vec![
76///     vec![1.0_f64, 2.0],
77///     vec![3.0, 4.0],
78/// ]);
79/// let factors = matrix_luci_factors_from_matrix(&m, None).unwrap();
80/// assert!(factors.rank >= 1);
81/// assert_eq!(factors.row_indices.len(), factors.rank);
82/// assert_eq!(factors.left.nrows(), m.nrows());
83/// assert_eq!(factors.right.ncols(), m.ncols());
84/// ```
85#[derive(Debug, Clone)]
86pub struct MatrixLuciFactors<T> {
87    /// Selected row indices.
88    pub row_indices: Vec<usize>,
89    /// Selected column indices.
90    pub col_indices: Vec<usize>,
91    /// Pivot error history.
92    pub pivot_errors: Vec<f64>,
93    /// Selected rank.
94    pub rank: usize,
95    /// Left factor.
96    pub left: Matrix<T>,
97    /// Right factor.
98    pub right: Matrix<T>,
99}
100
101pub(crate) fn map_backend_error(err: crate::matrixluci::MatrixLuciError) -> MatrixCIError {
102    match err {
103        crate::matrixluci::MatrixLuciError::InvalidArgument { message } => {
104            MatrixCIError::InvalidArgument { message }
105        }
106    }
107}
108
109fn factors_to_public<T>(
110    selection: PivotSelectionCore,
111    factors: CrossFactors<T>,
112    left_orthogonal: bool,
113) -> Result<MatrixLuciFactors<T>>
114where
115    T: Scalar + crate::MatrixLuciScalar,
116{
117    let left = if left_orthogonal {
118        factors.cols_solve_pivot().map_err(map_backend_error)?
119    } else {
120        factors.pivot_cols.clone()
121    };
122    let right = if left_orthogonal {
123        factors.pivot_rows.clone()
124    } else {
125        factors.solve_pivot_rows().map_err(map_backend_error)?
126    };
127
128    Ok(MatrixLuciFactors {
129        row_indices: selection.row_indices,
130        col_indices: selection.col_indices,
131        pivot_errors: selection.pivot_errors,
132        rank: selection.rank,
133        left,
134        right,
135    })
136}
137
138fn backend_linalg_error(message: impl Into<String>) -> MatrixCIError {
139    MatrixCIError::InvalidArgument {
140        message: message.into(),
141    }
142}
143
144fn index_range(start: usize, end: usize) -> Vec<usize> {
145    (start..end).collect()
146}
147
148fn identity_rect<T: Scalar>(nrows: usize, ncols: usize) -> Matrix<T> {
149    let mut result = Matrix::zeros(nrows, ncols);
150    for i in 0..nrows.min(ncols) {
151        result[[i, i]] = T::one();
152    }
153    result
154}
155
156fn apply_row_permutation<T: Scalar>(matrix: &Matrix<T>, permutation: &[usize]) -> Matrix<T> {
157    let mut result = Matrix::zeros(matrix.nrows(), matrix.ncols());
158    for col in 0..matrix.ncols() {
159        for (new_row, &old_row) in permutation.iter().enumerate().take(matrix.nrows()) {
160            result[[old_row, col]] = matrix[[new_row, col]];
161        }
162    }
163    result
164}
165
166fn apply_col_permutation<T: Scalar>(matrix: &Matrix<T>, permutation: &[usize]) -> Matrix<T> {
167    let mut result = Matrix::zeros(matrix.nrows(), matrix.ncols());
168    for (new_col, &old_col) in permutation.iter().enumerate().take(matrix.ncols()) {
169        for row in 0..matrix.nrows() {
170            result[[row, old_col]] = matrix[[row, new_col]];
171        }
172    }
173    result
174}
175
176pub(crate) fn rrlu_colmatrix<T>(lu: &RrLU<T>) -> Result<Matrix<T>>
177where
178    T: Scalar + crate::MatrixLuciScalar,
179{
180    let rank = lu.npivots();
181    if rank == 0 {
182        return Ok(Matrix::zeros(lu.nrows(), 0));
183    }
184    let rows = index_range(0, rank);
185    let cols = index_range(0, rank);
186    let right_pivot_cols = submatrix(lu.right_unpermuted(), &rows, &cols);
187    mat_mul_owned(lu.left(true), right_pivot_cols)
188        .map_err(|err| backend_linalg_error(format!("MatrixLUCI colmatrix multiply failed: {err}")))
189}
190
191pub(crate) fn rrlu_rowmatrix<T>(lu: &RrLU<T>) -> Result<Matrix<T>>
192where
193    T: Scalar + crate::MatrixLuciScalar,
194{
195    let rank = lu.npivots();
196    if rank == 0 {
197        return Ok(Matrix::zeros(0, lu.ncols()));
198    }
199    let rows = index_range(0, rank);
200    let cols = index_range(0, rank);
201    let left_pivot_rows = submatrix(lu.left_unpermuted(), &rows, &cols);
202    mat_mul_owned(left_pivot_rows, lu.right(true))
203        .map_err(|err| backend_linalg_error(format!("MatrixLUCI rowmatrix multiply failed: {err}")))
204}
205
206pub(crate) fn rrlu_cols_times_pivot_solve<T>(lu: &RrLU<T>) -> Result<Matrix<T>>
207where
208    T: Scalar + crate::MatrixLuciScalar,
209{
210    let rank = lu.npivots();
211    let mut result = identity_rect(lu.nrows(), rank);
212    if rank > 0 && rank < lu.nrows() {
213        let pivot_rows = index_range(0, rank);
214        let pivot_cols = index_range(0, rank);
215        let rest_rows = index_range(rank, lu.nrows());
216        let pivot = submatrix(lu.left_unpermuted(), &pivot_rows, &pivot_cols);
217        let rest = submatrix(lu.left_unpermuted(), &rest_rows, &pivot_cols);
218        let solved = triangular_solve_matrix_owned(pivot, rest, false, true, false, false)
219            .map_err(|err| {
220                backend_linalg_error(format!("MatrixLUCI lower triangular solve failed: {err}"))
221            })?;
222        for row in 0..solved.nrows() {
223            for col in 0..solved.ncols() {
224                result[[rank + row, col]] = solved[[row, col]];
225            }
226        }
227    }
228    Ok(apply_row_permutation(&result, lu.row_permutation()))
229}
230
231pub(crate) fn rrlu_pivot_solve_times_rows<T>(lu: &RrLU<T>) -> Result<Matrix<T>>
232where
233    T: Scalar + crate::MatrixLuciScalar,
234{
235    let rank = lu.npivots();
236    let mut result = identity_rect(rank, lu.ncols());
237    if rank > 0 && rank < lu.ncols() {
238        let pivot_rows = index_range(0, rank);
239        let pivot_cols = index_range(0, rank);
240        let rest_cols = index_range(rank, lu.ncols());
241        let pivot = submatrix(lu.right_unpermuted(), &pivot_rows, &pivot_cols);
242        let rest = submatrix(lu.right_unpermuted(), &pivot_rows, &rest_cols);
243        let solved = triangular_solve_matrix_owned(pivot, rest, true, false, false, false)
244            .map_err(|err| {
245                backend_linalg_error(format!("MatrixLUCI upper triangular solve failed: {err}"))
246            })?;
247        for row in 0..solved.nrows() {
248            for col in 0..solved.ncols() {
249                result[[row, rank + col]] = solved[[row, col]];
250            }
251        }
252    }
253    Ok(apply_col_permutation(&result, lu.col_permutation()))
254}
255
256fn factors_from_rrlu<T>(lu: &RrLU<T>) -> Result<MatrixLuciFactors<T>>
257where
258    T: Scalar + crate::MatrixLuciScalar,
259{
260    let left = if lu.is_left_orthogonal() {
261        rrlu_cols_times_pivot_solve(lu)?
262    } else {
263        rrlu_colmatrix(lu)?
264    };
265    let right = if lu.is_left_orthogonal() {
266        rrlu_rowmatrix(lu)?
267    } else {
268        rrlu_pivot_solve_times_rows(lu)?
269    };
270
271    Ok(MatrixLuciFactors {
272        row_indices: lu.row_indices(),
273        col_indices: lu.col_indices(),
274        pivot_errors: lu.pivot_errors(),
275        rank: lu.npivots(),
276        left,
277        right,
278    })
279}
280
281pub(crate) fn dense_matrix_luci_factors_from_matrix<T>(
282    a: &Matrix<T>,
283    options: RrLUOptions,
284) -> Result<MatrixLuciFactors<T>>
285where
286    T: Scalar + crate::MatrixLuciScalar,
287{
288    let lu = rrlu(a, Some(options))?;
289    factors_from_rrlu(&lu)
290}
291
292pub(crate) fn dense_matrix_luci_factors_from_matrix_owned<T>(
293    mut a: Matrix<T>,
294    options: RrLUOptions,
295) -> Result<MatrixLuciFactors<T>>
296where
297    T: Scalar + crate::MatrixLuciScalar,
298{
299    let lu = rrlu_mut(&mut a, Some(options))?;
300    factors_from_rrlu(&lu)
301}
302
303pub(crate) fn lazy_matrix_luci_factors_from_blocks<T, F>(
304    nrows: usize,
305    ncols: usize,
306    fill_block: F,
307    options: RrLUOptions,
308) -> Result<MatrixLuciFactors<T>>
309where
310    T: Scalar + crate::MatrixLuciScalar,
311    F: Fn(&[usize], &[usize], &mut [T]),
312    LazyBlockRookKernel: PivotKernel<T>,
313{
314    let source = LazyMatrixSource::new(nrows, ncols, fill_block);
315    let kernel_options = PivotKernelOptions {
316        max_bond_dim: options.max_bond_dim,
317        rel_tol: options.rel_tol,
318        abs_tol: options.abs_tol,
319        left_orthogonal: options.left_orthogonal,
320    };
321
322    let selection = LazyBlockRookKernel
323        .factorize(&source, &kernel_options)
324        .map_err(map_backend_error)?;
325    let factors = CrossFactors::from_source(&source, &selection).map_err(map_backend_error)?;
326    factors_to_public(selection, factors, options.left_orthogonal)
327}
328
329/// Factorize a dense matrix with MatrixLUCI.
330///
331/// Returns the selected pivot metadata together with left and
332/// right factors. This is the public facade used by higher-level crates.
333///
334/// # Arguments
335///
336/// * `a` - Dense matrix to factorize.
337/// * `options` - Optional rank and tolerance controls. `None` uses the
338///
339///   default LUCI settings.
340///
341/// # Returns
342///
343/// A [`MatrixLuciFactors`] value containing the selected pivot indices,
344/// error history, rank, and factors.
345///
346/// # Errors
347///
348/// Returns a [`MatrixCIError`] if the factorization fails, for example if the
349/// pivot block is singular or the backend rejects the input.
350///
351/// # Examples
352///
353/// ```
354/// use tensor4all_core::matrix_luci_factors_from_matrix;
355/// use tensor4all_tensorbackend::from_vec2d;
356///
357/// let m = from_vec2d(vec![
358///     vec![1.0_f64, 0.0],
359///     vec![0.0, 2.0],
360/// ]);
361/// let factors = matrix_luci_factors_from_matrix(&m, None).unwrap();
362/// assert_eq!(factors.left.nrows(), 2);
363/// assert_eq!(factors.right.ncols(), 2);
364/// assert_eq!(factors.row_indices.len(), factors.rank);
365/// ```
366pub fn matrix_luci_factors_from_matrix<T>(
367    a: &Matrix<T>,
368    options: Option<RrLUOptions>,
369) -> Result<MatrixLuciFactors<T>>
370where
371    T: Scalar + crate::MatrixLuciScalar,
372{
373    <T as crate::MatrixLuciScalar>::matrix_luci_factors_from_matrix(a, options.unwrap_or_default())
374}
375
376/// Factorize a dense matrix with MatrixLUCI while consuming the input matrix.
377///
378/// This is the owned-buffer counterpart of [`matrix_luci_factors_from_matrix`].
379/// It reuses the local matrix storage for the rrLU pivot selection step.
380///
381/// # Errors
382///
383/// Returns a [`MatrixCIError`] if the factorization fails, for example if the
384/// pivot block is singular or the backend rejects a factor solve.
385pub fn matrix_luci_factors_from_matrix_owned<T>(
386    a: Matrix<T>,
387    options: Option<RrLUOptions>,
388) -> Result<MatrixLuciFactors<T>>
389where
390    T: Scalar + crate::MatrixLuciScalar,
391{
392    dense_matrix_luci_factors_from_matrix_owned(a, options.unwrap_or_default())
393}
394
395/// Factorize a lazily supplied matrix with MatrixLUCI block-rook search.
396///
397/// The caller provides a block-fill closure that receives row and column
398/// index lists and writes the corresponding matrix block in column-major
399/// order.
400///
401/// # Arguments
402///
403/// * `nrows` - Number of matrix rows.
404/// * `ncols` - Number of matrix columns.
405/// * `fill_block` - Closure that fills `out` with `A[rows, cols]` in
406///
407///   column-major order.
408/// * `options` - Rank and tolerance controls.
409///
410/// # Returns
411///
412/// A [`MatrixLuciFactors`] value containing the pivot metadata and
413/// factors.
414///
415/// # Errors
416///
417/// Returns a [`MatrixCIError`] if the lazy factorization fails or if the
418/// block callback is inconsistent.
419///
420/// # Examples
421///
422/// ```
423/// use tensor4all_core::matrix_luci_factors_from_blocks;
424/// use tensor4all_core::RrLUOptions;
425///
426/// let factors = matrix_luci_factors_from_blocks(
427///     2,
428///     2,
429///     |rows, cols, out| {
430///         for (j, &col) in cols.iter().enumerate() {
431///             for (i, &row) in rows.iter().enumerate() {
432///                 out[i + rows.len() * j] = if row == col { 1.0 } else { 0.0 };
433///             }
434///         }
435///     },
436///     RrLUOptions::default(),
437/// ).unwrap();
438/// assert_eq!(factors.rank, 2);
439/// ```
440pub fn matrix_luci_factors_from_blocks<T, F>(
441    nrows: usize,
442    ncols: usize,
443    fill_block: F,
444    options: RrLUOptions,
445) -> Result<MatrixLuciFactors<T>>
446where
447    T: Scalar + crate::MatrixLuciScalar,
448    F: Fn(&[usize], &[usize], &mut [T]),
449{
450    <T as crate::MatrixLuciScalar>::matrix_luci_factors_from_blocks(
451        nrows, ncols, fill_block, options,
452    )
453}
454
455impl<T> MatrixLUCI<T>
456where
457    T: Scalar + crate::MatrixLuciScalar,
458{
459    /// Create a MatrixLUCI from a dense matrix.
460    ///
461    /// # Errors
462    ///
463    /// Returns an error when the construction or conversion fails (a shape or
464    /// /// index mismatch, or a backend failure).
465    ///
466    /// # Examples
467    ///
468    /// ```
469    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
470    /// use tensor4all_tensorbackend::from_vec2d;
471    ///
472    /// let m = from_vec2d(vec![
473    ///     vec![2.0_f64, 0.0],
474    ///     vec![0.0, 3.0],
475    /// ]);
476    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
477    /// assert!(ci.rank() >= 1);
478    /// ```
479    pub fn from_matrix(a: &Matrix<T>, options: Option<RrLUOptions>) -> Result<Self> {
480        let factors = matrix_luci_factors_from_matrix(a, options)?;
481
482        Ok(Self {
483            nrows: a.nrows(),
484            ncols: a.ncols(),
485            row_indices: factors.row_indices,
486            col_indices: factors.col_indices,
487            left: factors.left,
488            right: factors.right,
489            pivot_errors: factors.pivot_errors,
490        })
491    }
492
493    /// Left CI factor (shape: `nrows x rank`).
494    ///
495    /// The approximation is `left * right`.
496    ///
497    /// # Examples
498    ///
499    /// ```
500    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
501    /// use tensor4all_tensorbackend::{from_vec2d, mat_mul};
502    ///
503    /// let m = from_vec2d(vec![
504    ///     vec![1.0_f64, 2.0],
505    ///     vec![3.0, 4.0],
506    /// ]);
507    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
508    /// let reconstructed = mat_mul(&ci.left(), &ci.right()).unwrap();
509    /// for i in 0..2 {
510    ///     for j in 0..2 {
511    ///         assert!((reconstructed[[i, j]] - m[[i, j]]).abs() < 1e-10);
512    ///     }
513    /// }
514    /// ```
515    pub fn left(&self) -> Matrix<T> {
516        self.left.clone()
517    }
518
519    /// Right CI factor (shape: `rank x ncols`).
520    ///
521    /// The approximation is `left * right`.
522    ///
523    /// # Examples
524    ///
525    /// ```
526    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
527    /// use tensor4all_tensorbackend::{from_vec2d, mat_mul};
528    ///
529    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
530    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
531    /// let r = ci.right();
532    /// assert_eq!(r.nrows(), ci.rank());
533    /// assert_eq!(r.ncols(), ci.ncols());
534    /// // left * right reconstructs the matrix
535    /// let recon = mat_mul(&ci.left(), &r).unwrap();
536    /// for i in 0..2 {
537    ///     for j in 0..2 {
538    ///         assert!((recon[[i, j]] - m[[i, j]]).abs() < 1e-10);
539    ///     }
540    /// }
541    /// ```
542    pub fn right(&self) -> Matrix<T> {
543        self.right.clone()
544    }
545
546    /// Pivot error history (one entry per pivot, plus a final residual estimate).
547    ///
548    /// # Examples
549    ///
550    /// ```
551    /// use tensor4all_core::MatrixLUCI;
552    /// use tensor4all_tensorbackend::from_vec2d;
553    ///
554    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
555    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
556    /// let errs = ci.pivot_errors();
557    /// assert!(!errs.is_empty());
558    /// // All errors are non-negative
559    /// for &e in &errs {
560    ///     assert!(e >= 0.0);
561    /// }
562    /// ```
563    pub fn pivot_errors(&self) -> Vec<f64> {
564        self.pivot_errors.clone()
565    }
566
567    /// Last pivot error (the residual estimate after all pivots).
568    ///
569    /// # Examples
570    ///
571    /// ```
572    /// use tensor4all_core::MatrixLUCI;
573    /// use tensor4all_tensorbackend::from_vec2d;
574    ///
575    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
576    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
577    /// let err = ci.last_pivot_error();
578    /// assert!(err >= 0.0);
579    /// ```
580    pub fn last_pivot_error(&self) -> f64 {
581        self.pivot_errors.last().copied().unwrap_or(0.0)
582    }
583}
584
585impl<T> AbstractMatrixCI<T> for MatrixLUCI<T>
586where
587    T: Scalar + crate::MatrixLuciScalar,
588{
589    fn nrows(&self) -> usize {
590        self.nrows
591    }
592
593    fn ncols(&self) -> usize {
594        self.ncols
595    }
596
597    fn rank(&self) -> usize {
598        self.row_indices.len()
599    }
600
601    fn row_indices(&self) -> &[usize] {
602        &self.row_indices
603    }
604
605    fn col_indices(&self) -> &[usize] {
606        &self.col_indices
607    }
608
609    fn evaluate(&self, i: usize, j: usize) -> T {
610        let mut sum = T::zero();
611        for k in 0..self.rank() {
612            sum = sum + self.left[[i, k]] * self.right[[k, j]];
613        }
614        sum
615    }
616
617    fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
618        let r = self.rank();
619        let mut result = Matrix::zeros(rows.len(), cols.len());
620        for (j_out, &col) in cols.iter().enumerate() {
621            for (i_out, &row) in rows.iter().enumerate() {
622                let mut sum = T::zero();
623                for k in 0..r {
624                    sum = sum + self.left[[row, k]] * self.right[[k, col]];
625                }
626                result[[i_out, j_out]] = sum;
627            }
628        }
629        result
630    }
631}
632
633#[cfg(test)]
634mod tests;
tensor4all_core/matrix_luci.rs

tensor4all_core/
matrix_luci.rs
