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::{
9    rrlu, rrlu_left_mut, rrlu_mut, validate_luci_options, validate_owned_luci_input, RrLU,
10    RrLUOptions,
11};
12use crate::matrixluci::block_rook::LazyBlockRookKernel;
13use crate::matrixluci::factors::CrossFactors;
14use crate::matrixluci::source::LazyMatrixSource;
15use crate::matrixluci::types::{PivotKernelOptions, PivotSelectionCore};
16use crate::matrixluci::PivotKernel;
17use crate::scalar::Scalar;
18use crate::traits::AbstractMatrixCI;
19use tensor4all_tensorbackend::{
20    explicit::Session, mat_mul_owned, submatrix, triangular_solve_matrix_owned,
21    CpuExecutionContext, Matrix,
22};
23
24/// Where MatrixLUCI factor assembly runs its matrix products and solves.
25///
26/// `Global` keeps the compatibility routes of the legacy entry points;
27/// `Context` runs every operation in one explicit concrete session of the
28/// caller's context and never consults the process-global default.
29pub(crate) enum Linalg<'a, 'session> {
30    Global,
31    Context(&'a CpuExecutionContext),
32    /// A caller-supplied session, entered by the caller rather than opened here.
33    ///
34    /// This is the form a worker lane uses: the session was entered by the phase that owns
35    /// the worker, so re-entering a context here would either be refused or pick another
36    /// backend.
37    Session(&'a mut Session<'session>),
38}
39
40impl Linalg<'_, '_> {
41    fn mat_mul<T>(&mut self, a: Matrix<T>, b: Matrix<T>) -> std::result::Result<Matrix<T>, String>
42    where
43        T: Scalar + crate::MatrixLuciScalar,
44    {
45        match self {
46            Self::Global => mat_mul_owned(a, b).map_err(|err| err.to_string()),
47            Self::Context(context) => context
48                .with_concrete_session(|session| session.mat_mul(&a, &b))
49                .map_err(|err| err.to_string())?
50                .map_err(|err| err.to_string()),
51            Self::Session(session) => session.mat_mul(&a, &b).map_err(|err| err.to_string()),
52        }
53    }
54
55    fn triangular_solve<T>(
56        &mut self,
57        a: Matrix<T>,
58        b: Matrix<T>,
59        left_side: bool,
60        lower: bool,
61    ) -> std::result::Result<Matrix<T>, String>
62    where
63        T: Scalar + crate::MatrixLuciScalar,
64    {
65        match self {
66            Self::Global => triangular_solve_matrix_owned(a, b, left_side, lower, false, false)
67                .map_err(|err| err.to_string()),
68            Self::Context(context) => context
69                .with_concrete_session(|session| {
70                    T::matrix_luci_triangular_solve_in(session, &a, &b, left_side, lower)
71                })
72                .map_err(|err| err.to_string())?
73                .map_err(|err| err.to_string()),
74            Self::Session(session) => {
75                T::matrix_luci_triangular_solve_in(session, &a, &b, left_side, lower)
76                    .map_err(|err| err.to_string())
77            }
78        }
79    }
80}
81
82/// Matrix LU-based Cross Interpolation.
83///
84/// This is a higher-level [`Matrix`] wrapper around the lower-level `matrixluci`
85/// substrate.
86///
87/// # Examples
88///
89/// ```
90/// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
91/// use tensor4all_tensorbackend::from_vec2d;
92///
93/// let m = from_vec2d(vec![
94///     vec![1.0_f64, 2.0, 3.0],
95///     vec![4.0, 5.0, 6.0],
96///     vec![7.0, 8.0, 9.0],
97/// ]);
98///
99/// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
100/// // The approximation must have at most rank min(nrows, ncols)
101/// assert!(ci.rank() <= 3);
102/// // Reconstructed matrix should match original at pivot positions
103/// let row_indices = ci.row_indices().to_vec();
104/// let col_indices = ci.col_indices().to_vec();
105/// for (&i, &j) in row_indices.iter().zip(col_indices.iter()) {
106///     let approx = ci.evaluate(i, j);
107///     let exact = m[[i, j]];
108///     assert!((approx - exact).abs() < 1e-10);
109/// }
110/// ```
111#[derive(Debug, Clone)]
112pub struct MatrixLUCI<T: Scalar + crate::MatrixLuciScalar> {
113    nrows: usize,
114    ncols: usize,
115    row_indices: Vec<usize>,
116    col_indices: Vec<usize>,
117    left: Matrix<T>,
118    right: Matrix<T>,
119    pivot_errors: Vec<f64>,
120}
121
122/// High-level factors produced by MatrixLUCI.
123///
124/// This is the public result of the LUCI factorization facade. It exposes
125/// the selected pivot metadata and the left and right factors
126/// needed by higher-level tensor-network code without exposing the low-level
127/// pivot kernel substrate.
128///
129/// Related types: [`MatrixLUCI`] is the owning CI wrapper, while
130/// [`MatrixACA`](crate::MatrixACA) and [`RrLU`](crate::RrLU) are related
131/// matrix factorization entry points.
132///
133/// # Examples
134///
135/// ```
136/// use tensor4all_core::matrix_luci_factors_from_matrix;
137/// use tensor4all_tensorbackend::from_vec2d;
138///
139/// let m = from_vec2d(vec![
140///     vec![1.0_f64, 2.0],
141///     vec![3.0, 4.0],
142/// ]);
143/// let factors = matrix_luci_factors_from_matrix(&m, None).unwrap();
144/// assert!(factors.rank >= 1);
145/// assert_eq!(factors.row_indices.len(), factors.rank);
146/// assert_eq!(factors.left.nrows(), m.nrows());
147/// assert_eq!(factors.right.ncols(), m.ncols());
148/// ```
149#[derive(Debug, Clone)]
150pub struct MatrixLuciFactors<T> {
151    /// Selected row indices.
152    pub row_indices: Vec<usize>,
153    /// Selected column indices.
154    pub col_indices: Vec<usize>,
155    /// Pivot error history, or one measured maximum residual for fixed pivots.
156    pub pivot_errors: Vec<f64>,
157    /// Selected rank.
158    pub rank: usize,
159    /// Left factor.
160    pub left: Matrix<T>,
161    /// Right factor.
162    pub right: Matrix<T>,
163}
164
165/// Row interpolation `A ~= interpolation * A[rows, :]` selected by MatrixLUCI.
166/// Unlike [`MatrixLuciFactors`], this does not build an unused right factor.
167///
168/// # Examples
169/// ```
170/// use tensor4all_core::matrix_luci_row_interpolation_owned_in;
171/// use tensor4all_tensorbackend::{default_cpu_execution_context, Matrix};
172/// let a = Matrix::from_col_major_vec(2, 1, vec![2.0_f64, 4.0]);
173/// let id = matrix_luci_row_interpolation_owned_in(a, None, &default_cpu_execution_context())?;
174/// assert_eq!(id.rows, vec![1]);
175/// assert_eq!(id.interpolation.as_col_major_slice(), &[0.5, 1.0]);
176/// # Ok::<(), Box<dyn std::error::Error>>(())
177/// ```
178#[derive(Clone, Debug)]
179pub struct MatrixLuciRowInterpolation<T> {
180    /// Original matrix row indices, in pivot order.
181    pub rows: Vec<usize>,
182    /// Row weights, with an identity submatrix at the selected rows.
183    pub interpolation: Matrix<T>,
184    /// Accepted LU pivot magnitudes followed by the maximum absolute entry of
185    /// the remaining Schur complement. This is a local matrix diagnostic.
186    pub pivot_magnitudes: Vec<f64>,
187}
188
189/// Selects a row interpolant without constructing the complementary right factor.
190///
191/// Consumes column-major `a`, uses `options` for pivot rank/tolerances, and
192/// performs the factor solve in the supplied CPU `context`. `left_orthogonal`
193/// is always set to true because this operation returns row interpolation.
194/// A zero matrix yields zero columns; the caller chooses its zero-rank convention.
195///
196/// # Errors
197/// Returns [`MatrixCIError::InvalidArgument`] for negative or nonfinite
198/// tolerances, infinite input magnitudes, or backend solve failures; returns
199/// [`MatrixCIError::NaNEncountered`] for NaN input or nonfinite factors or
200/// residual magnitudes.
201///
202/// # Examples
203/// ```
204/// use tensor4all_core::{matrix_luci_row_interpolation_owned_in, RrLUOptions};
205/// use tensor4all_tensorbackend::{default_cpu_execution_context, Matrix};
206/// let a = Matrix::from_col_major_vec(2, 2, vec![1.0_f64, 2.0, 3.0, 6.0]);
207/// let id = matrix_luci_row_interpolation_owned_in(a, Some(RrLUOptions::default()),
208///     &default_cpu_execution_context())?;
209/// assert_eq!(id.rows, vec![1]);
210/// assert_eq!(id.interpolation.as_col_major_slice(), &[0.5, 1.0]);
211/// # Ok::<(), Box<dyn std::error::Error>>(())
212/// ```
213pub fn matrix_luci_row_interpolation_owned_in<T>(
214    mut a: Matrix<T>,
215    options: Option<RrLUOptions>,
216    context: &CpuExecutionContext,
217) -> Result<MatrixLuciRowInterpolation<T>>
218where
219    T: Scalar + crate::MatrixLuciScalar,
220{
221    let mut options = options.unwrap_or_default();
222    options.left_orthogonal = true;
223    validate_owned_luci_input(&a, &options)?;
224    let lu = rrlu_left_mut(&mut a, Some(options))?;
225    let rank = lu.n_pivot;
226    Ok(MatrixLuciRowInterpolation {
227        rows: lu.row_permutation[..rank].to_vec(),
228        interpolation: cols_times_pivot_solve(
229            &lu.l,
230            &lu.row_permutation,
231            rank,
232            &mut Linalg::Context(context),
233        )?,
234        pivot_magnitudes: lu.pivot_magnitudes,
235    })
236}
237
238// Validation for the new owned/context facades lives with the rrLU entry
239// points (`matrixlu::validate_owned_luci_input`) so the selection-only path
240// checks inputs and tolerances identically.
241
242pub(crate) fn map_backend_error(err: crate::matrixluci::MatrixLuciError) -> MatrixCIError {
243    match err {
244        crate::matrixluci::MatrixLuciError::NaNEncountered { matrix } => {
245            MatrixCIError::NaNEncountered {
246                matrix: matrix.to_string(),
247            }
248        }
249        crate::matrixluci::MatrixLuciError::InvalidArgument { message } => {
250            MatrixCIError::InvalidArgument { message }
251        }
252    }
253}
254
255fn factors_to_public<T>(
256    selection: PivotSelectionCore,
257    factors: CrossFactors<T>,
258    left_orthogonal: bool,
259) -> Result<MatrixLuciFactors<T>>
260where
261    T: Scalar + crate::MatrixLuciScalar,
262{
263    let left = if left_orthogonal {
264        factors.cols_solve_pivot().map_err(map_backend_error)?
265    } else {
266        factors.pivot_cols.clone()
267    };
268    let right = if left_orthogonal {
269        factors.pivot_rows.clone()
270    } else {
271        factors.solve_pivot_rows().map_err(map_backend_error)?
272    };
273
274    Ok(MatrixLuciFactors {
275        row_indices: selection.row_indices,
276        col_indices: selection.col_indices,
277        pivot_errors: selection.pivot_errors,
278        rank: selection.rank,
279        left,
280        right,
281    })
282}
283
284fn backend_linalg_error(message: impl Into<String>) -> MatrixCIError {
285    MatrixCIError::InvalidArgument {
286        message: message.into(),
287    }
288}
289
290fn index_range(start: usize, end: usize) -> Vec<usize> {
291    (start..end).collect()
292}
293
294fn identity_rect<T: Scalar>(nrows: usize, ncols: usize) -> Matrix<T> {
295    let mut result = Matrix::zeros(nrows, ncols);
296    for i in 0..nrows.min(ncols) {
297        result[[i, i]] = T::one();
298    }
299    result
300}
301
302fn apply_row_permutation<T: Scalar>(matrix: &Matrix<T>, permutation: &[usize]) -> Matrix<T> {
303    let mut result = Matrix::zeros(matrix.nrows(), matrix.ncols());
304    for col in 0..matrix.ncols() {
305        for (new_row, &old_row) in permutation.iter().enumerate().take(matrix.nrows()) {
306            result[[old_row, col]] = matrix[[new_row, col]];
307        }
308    }
309    result
310}
311
312fn apply_col_permutation<T: Scalar>(matrix: &Matrix<T>, permutation: &[usize]) -> Matrix<T> {
313    let mut result = Matrix::zeros(matrix.nrows(), matrix.ncols());
314    for (new_col, &old_col) in permutation.iter().enumerate().take(matrix.ncols()) {
315        for row in 0..matrix.nrows() {
316            result[[row, old_col]] = matrix[[row, new_col]];
317        }
318    }
319    result
320}
321
322pub(crate) fn rrlu_colmatrix<T>(lu: &RrLU<T>, linalg: &mut Linalg<'_, '_>) -> Result<Matrix<T>>
323where
324    T: Scalar + crate::MatrixLuciScalar,
325{
326    let rank = lu.npivots();
327    if rank == 0 {
328        return Ok(Matrix::zeros(lu.nrows(), 0));
329    }
330    let rows = index_range(0, rank);
331    let cols = index_range(0, rank);
332    let right_pivot_cols = submatrix(lu.right_unpermuted(), &rows, &cols);
333    linalg
334        .mat_mul(lu.left(true), right_pivot_cols)
335        .map_err(|err| backend_linalg_error(format!("MatrixLUCI colmatrix multiply failed: {err}")))
336}
337
338pub(crate) fn rrlu_rowmatrix<T>(lu: &RrLU<T>, linalg: &mut Linalg<'_, '_>) -> Result<Matrix<T>>
339where
340    T: Scalar + crate::MatrixLuciScalar,
341{
342    let rank = lu.npivots();
343    if rank == 0 {
344        return Ok(Matrix::zeros(0, lu.ncols()));
345    }
346    let rows = index_range(0, rank);
347    let cols = index_range(0, rank);
348    let left_pivot_rows = submatrix(lu.left_unpermuted(), &rows, &cols);
349    linalg
350        .mat_mul(left_pivot_rows, lu.right(true))
351        .map_err(|err| backend_linalg_error(format!("MatrixLUCI rowmatrix multiply failed: {err}")))
352}
353
354pub(crate) fn rrlu_cols_times_pivot_solve<T>(
355    lu: &RrLU<T>,
356    linalg: &mut Linalg<'_, '_>,
357) -> Result<Matrix<T>>
358where
359    T: Scalar + crate::MatrixLuciScalar,
360{
361    cols_times_pivot_solve(
362        lu.left_unpermuted(),
363        lu.row_permutation(),
364        lu.npivots(),
365        linalg,
366    )
367}
368
369fn cols_times_pivot_solve<T>(
370    left: &Matrix<T>,
371    row_permutation: &[usize],
372    rank: usize,
373    linalg: &mut Linalg<'_, '_>,
374) -> Result<Matrix<T>>
375where
376    T: Scalar + crate::MatrixLuciScalar,
377{
378    let mut result = identity_rect(left.nrows(), rank);
379    if rank > 0 && rank < left.nrows() {
380        let pivot_rows = index_range(0, rank);
381        let pivot_cols = index_range(0, rank);
382        let rest_rows = index_range(rank, left.nrows());
383        let pivot = submatrix(left, &pivot_rows, &pivot_cols);
384        let rest = submatrix(left, &rest_rows, &pivot_cols);
385        let solved = linalg
386            .triangular_solve(pivot, rest, false, true)
387            .map_err(|err| {
388                backend_linalg_error(format!("MatrixLUCI lower triangular solve failed: {err}"))
389            })?;
390        for row in 0..solved.nrows() {
391            for col in 0..solved.ncols() {
392                result[[rank + row, col]] = solved[[row, col]];
393            }
394        }
395    }
396    Ok(apply_row_permutation(&result, row_permutation))
397}
398
399pub(crate) fn rrlu_pivot_solve_times_rows<T>(
400    lu: &RrLU<T>,
401    linalg: &mut Linalg<'_, '_>,
402) -> Result<Matrix<T>>
403where
404    T: Scalar + crate::MatrixLuciScalar,
405{
406    let rank = lu.npivots();
407    let mut result = identity_rect(rank, lu.ncols());
408    if rank > 0 && rank < lu.ncols() {
409        let pivot_rows = index_range(0, rank);
410        let pivot_cols = index_range(0, rank);
411        let rest_cols = index_range(rank, lu.ncols());
412        let pivot = submatrix(lu.right_unpermuted(), &pivot_rows, &pivot_cols);
413        let rest = submatrix(lu.right_unpermuted(), &pivot_rows, &rest_cols);
414        let solved = linalg
415            .triangular_solve(pivot, rest, true, false)
416            .map_err(|err| {
417                backend_linalg_error(format!("MatrixLUCI upper triangular solve failed: {err}"))
418            })?;
419        for row in 0..solved.nrows() {
420            for col in 0..solved.ncols() {
421                result[[row, rank + col]] = solved[[row, col]];
422            }
423        }
424    }
425    Ok(apply_col_permutation(&result, lu.col_permutation()))
426}
427
428fn factors_from_rrlu<T>(lu: &RrLU<T>, linalg: &mut Linalg<'_, '_>) -> Result<MatrixLuciFactors<T>>
429where
430    T: Scalar + crate::MatrixLuciScalar,
431{
432    let left = if lu.is_left_orthogonal() {
433        rrlu_cols_times_pivot_solve(lu, linalg)?
434    } else {
435        rrlu_colmatrix(lu, linalg)?
436    };
437    let right = if lu.is_left_orthogonal() {
438        rrlu_rowmatrix(lu, linalg)?
439    } else {
440        rrlu_pivot_solve_times_rows(lu, linalg)?
441    };
442
443    Ok(MatrixLuciFactors {
444        row_indices: lu.row_indices(),
445        col_indices: lu.col_indices(),
446        pivot_errors: lu.pivot_errors(),
447        rank: lu.npivots(),
448        left,
449        right,
450    })
451}
452
453pub(crate) fn dense_matrix_luci_factors_from_matrix<T>(
454    a: &Matrix<T>,
455    options: RrLUOptions,
456) -> Result<MatrixLuciFactors<T>>
457where
458    T: Scalar + crate::MatrixLuciScalar,
459{
460    validate_owned_luci_input(a, &options)?;
461    let lu = rrlu(a, Some(options))?;
462    factors_from_rrlu(&lu, &mut Linalg::Global)
463}
464
465pub(crate) fn dense_matrix_luci_factors_from_matrix_owned<T>(
466    mut a: Matrix<T>,
467    options: RrLUOptions,
468    mut linalg: Linalg<'_, '_>,
469) -> Result<MatrixLuciFactors<T>>
470where
471    T: Scalar + crate::MatrixLuciScalar,
472{
473    let lu = rrlu_mut(&mut a, Some(options))?;
474    factors_from_rrlu(&lu, &mut linalg)
475}
476
477pub(crate) fn lazy_matrix_luci_factors_from_blocks<T, F>(
478    nrows: usize,
479    ncols: usize,
480    fill_block: F,
481    options: RrLUOptions,
482) -> Result<MatrixLuciFactors<T>>
483where
484    T: Scalar + crate::MatrixLuciScalar,
485    F: Fn(&[usize], &[usize], &mut [T]),
486    LazyBlockRookKernel: PivotKernel<T>,
487{
488    validate_luci_options(&options)?;
489    let source = LazyMatrixSource::new(nrows, ncols, fill_block);
490    let kernel_options = PivotKernelOptions {
491        max_bond_dim: options.max_bond_dim,
492        rel_tol: options.rel_tol,
493        abs_tol: options.abs_tol,
494        left_orthogonal: options.left_orthogonal,
495    };
496
497    let selection = LazyBlockRookKernel
498        .factorize(&source, &kernel_options)
499        .map_err(map_backend_error)?;
500    let factors = CrossFactors::from_source(&source, &selection).map_err(map_backend_error)?;
501    factors_to_public(selection, factors, options.left_orthogonal)
502}
503
504/// Factorize a dense matrix with MatrixLUCI.
505///
506/// Returns the selected pivot metadata together with left and
507/// right factors. This is the public facade used by higher-level crates.
508///
509/// # Arguments
510///
511/// * `a` - Dense matrix to factorize.
512/// * `options` - Optional rank and tolerance controls. `None` uses the
513///
514///   default LUCI settings.
515///
516/// # Returns
517///
518/// A [`MatrixLuciFactors`] value containing the selected pivot indices,
519/// error history, rank, and factors.
520///
521/// # Errors
522///
523/// Returns [`MatrixCIError::InvalidArgument`] for negative or non-finite
524/// tolerances, infinite input magnitudes, or backend failures; returns
525/// [`MatrixCIError::NaNEncountered`] for NaN input or non-finite factors.
526///
527/// # Examples
528///
529/// ```
530/// use tensor4all_core::matrix_luci_factors_from_matrix;
531/// use tensor4all_tensorbackend::from_vec2d;
532///
533/// let m = from_vec2d(vec![
534///     vec![1.0_f64, 0.0],
535///     vec![0.0, 2.0],
536/// ]);
537/// let factors = matrix_luci_factors_from_matrix(&m, None).unwrap();
538/// assert_eq!(factors.left.nrows(), 2);
539/// assert_eq!(factors.right.ncols(), 2);
540/// assert_eq!(factors.row_indices.len(), factors.rank);
541/// ```
542pub fn matrix_luci_factors_from_matrix<T>(
543    a: &Matrix<T>,
544    options: Option<RrLUOptions>,
545) -> Result<MatrixLuciFactors<T>>
546where
547    T: Scalar + crate::MatrixLuciScalar,
548{
549    <T as crate::MatrixLuciScalar>::matrix_luci_factors_from_matrix(a, options.unwrap_or_default())
550}
551
552/// Factorize a dense matrix with MatrixLUCI while consuming the input matrix.
553///
554/// This is the owned-buffer counterpart of [`matrix_luci_factors_from_matrix`].
555/// It reuses the local matrix storage for the rrLU pivot selection step.
556///
557/// # Errors
558///
559/// Returns a [`MatrixCIError`] if the factorization fails, for example if the
560/// pivot block is singular or the backend rejects a factor solve.
561pub fn matrix_luci_factors_from_matrix_owned<T>(
562    a: Matrix<T>,
563    options: Option<RrLUOptions>,
564) -> Result<MatrixLuciFactors<T>>
565where
566    T: Scalar + crate::MatrixLuciScalar,
567{
568    let options = options.unwrap_or_default();
569    validate_owned_luci_input(&a, &options)?;
570    dense_matrix_luci_factors_from_matrix_owned(a, options, Linalg::Global)
571}
572
573/// Factorize a dense matrix with MatrixLUCI in a caller-owned CPU context.
574///
575/// Context-scoped counterpart of [`matrix_luci_factors_from_matrix_owned`].
576/// The rank-revealing pivot search is identical; the left and right factors
577/// are assembled with `context`'s configured backend (matrix products and
578/// triangular solves). No process-global context is consulted.
579///
580/// # Arguments
581///
582/// * `a` - Dense column-major matrix to factorize; its storage is reused.
583/// * `options` - Optional rank and tolerance controls. `None` uses the
584///   default LUCI settings.
585/// * `context` - Caller-owned CPU execution context running the factor
586///   assembly.
587///
588/// # Returns
589///
590/// A [`MatrixLuciFactors`] value with the selected pivots, pivot error
591/// history, rank, and factors satisfying `a ~= left * right`.
592///
593/// # Errors
594///
595/// Returns [`MatrixCIError::InvalidArgument`] for negative or nonfinite
596/// tolerances, infinite input magnitudes, or backend solve failures; returns
597/// [`MatrixCIError::NaNEncountered`] for NaN input or factors.
598///
599/// # Examples
600///
601/// ```
602/// use tensor4all_core::{matrix_luci_factors_from_matrix_owned_in, RrLUOptions};
603/// use tensor4all_tensorbackend::{from_vec2d, mat_mul, CpuExecutionContext};
604/// use tensor4all_tensorbackend::CpuBackend;
605///
606/// let context = CpuExecutionContext::from_backend(CpuBackend::new());
607/// // Rank-1 matrix [1 2; 2 4].
608/// let a = from_vec2d(vec![vec![1.0_f64, 2.0], vec![2.0, 4.0]]);
609/// let options = RrLUOptions { left_orthogonal: true, ..RrLUOptions::default() };
610/// let factors = matrix_luci_factors_from_matrix_owned_in(a.clone(), Some(options), &context)?;
611/// assert_eq!(factors.rank, 1);
612/// let rebuilt = mat_mul(&factors.left, &factors.right)?;
613/// for (x, y) in rebuilt.as_col_major_slice().iter().zip(a.as_col_major_slice()) {
614///     assert!((x - y).abs() < 1e-12);
615/// }
616/// # Ok::<(), Box<dyn std::error::Error>>(())
617/// ```
618pub fn matrix_luci_factors_from_matrix_owned_in<T>(
619    a: Matrix<T>,
620    options: Option<RrLUOptions>,
621    context: &CpuExecutionContext,
622) -> Result<MatrixLuciFactors<T>>
623where
624    T: Scalar + crate::MatrixLuciScalar,
625{
626    let options = options.unwrap_or_default();
627    validate_owned_luci_input(&a, &options)?;
628    dense_matrix_luci_factors_from_matrix_owned(a, options, Linalg::Context(context))
629}
630
631/// Factorize a dense matrix with MatrixLUCI in a session the caller already entered.
632///
633/// Session-taking counterpart of [`matrix_luci_factors_from_matrix_owned_in`], for callers that own their
634/// execution rather than a context: a worker lane entered by the phase that owns it, for example, whose session
635/// must be used as it is. The rank-revealing pivot search is identical; the factor assembly's matrix products
636/// and triangular solves run on `session`, and neither a process-global context nor a fresh session entry is
637/// consulted.
638///
639/// # Arguments
640///
641/// * `a` - Dense column-major matrix to factorize; its storage is reused.
642/// * `options` - Optional rank and tolerance controls. `None` uses the default LUCI settings.
643/// * `session` - An already-entered explicit session that runs the factor assembly.
644///
645/// # Returns
646///
647/// A [`MatrixLuciFactors`] value with the selected pivots, pivot error history, rank, and factors satisfying
648/// `a ~= left * right`.
649///
650/// # Errors
651///
652/// Returns [`MatrixCIError::InvalidArgument`] for negative or nonfinite tolerances, infinite input magnitudes,
653/// or backend solve failures; returns [`MatrixCIError::NaNEncountered`] for NaN input or factors.
654///
655/// # Examples
656///
657/// ```
658/// use tensor4all_core::{matrix_luci_factors_from_matrix_owned_on_session, RrLUOptions};
659/// use tensor4all_tensorbackend::{from_vec2d, mat_mul, CpuExecutionContext};
660/// use tensor4all_tensorbackend::CpuBackend;
661///
662/// let context = CpuExecutionContext::from_backend(CpuBackend::new());
663/// let a = from_vec2d(vec![vec![1.0_f64, 2.0], vec![2.0, 4.0]]);
664/// let factors = context.with_concrete_session(|session| {
665///     matrix_luci_factors_from_matrix_owned_on_session(
666///         a.clone(),
667///         Some(RrLUOptions { left_orthogonal: true, ..RrLUOptions::default() }),
668///         session,
669///     )
670/// })??;
671/// assert_eq!(factors.rank, 1);
672/// # Ok::<(), Box<dyn std::error::Error>>(())
673/// ```
674pub fn matrix_luci_factors_from_matrix_owned_on_session<T>(
675    a: Matrix<T>,
676    options: Option<RrLUOptions>,
677    session: &mut Session<'_>,
678) -> Result<MatrixLuciFactors<T>>
679where
680    T: Scalar + crate::MatrixLuciScalar,
681{
682    let options = options.unwrap_or_default();
683    validate_owned_luci_input(&a, &options)?;
684    dense_matrix_luci_factors_from_matrix_owned(a, options, Linalg::Session(session))
685}
686
687/// Conservative logical working-byte estimate for owned dense LUCI factors.
688///
689/// Includes the consumed matrix, rrLU L/U, factor construction, permutation
690/// and triangular-solve operands/results, and pivot/index metadata. The
691/// factor rank is bounded by `min(nrows, ncols, max_rank)`. The estimate sums
692/// potentially overlapping factor buffers conservatively; it excludes
693/// allocator overhead and provider-private workspace. Callers must add their
694/// own buffers retained across the factorization.
695///
696/// Returns `None` when checked shape/byte arithmetic overflows `usize`.
697///
698/// # Examples
699///
700/// ```
701/// use tensor4all_core::matrix_luci_factors_working_bytes;
702/// let bytes = matrix_luci_factors_working_bytes::<f64>(32, 32, 32).unwrap();
703/// // Even just the matrix and full-rank L/U require three matrix payloads.
704/// assert!(bytes >= 3 * 32 * 32 * std::mem::size_of::<f64>());
705/// assert!(matrix_luci_factors_working_bytes::<f64>(usize::MAX, 2, 1).is_none());
706/// ```
707pub fn matrix_luci_factors_working_bytes<T: Scalar>(
708    nrows: usize,
709    ncols: usize,
710    max_rank: usize,
711) -> Option<usize> {
712    let rank = max_rank.min(nrows).min(ncols);
713    let axes = nrows.checked_add(ncols)?;
714    let matrix = nrows.checked_mul(ncols)?;
715    // L/U, both output factors, permutation copies, and triangular solve
716    // operands/results. Charge both orientations rather than assuming which
717    // side has the larger solve, including Matrix/backend operand copies.
718    let factors = axes.checked_mul(rank)?.checked_mul(6)?;
719    let pivots = rank.checked_mul(rank)?.checked_mul(2)?;
720    let scalar_bytes = matrix
721        .checked_add(factors)?
722        .checked_add(pivots)?
723        .checked_mul(std::mem::size_of::<T>())?;
724    let index_bytes = axes
725        .checked_add(rank)?
726        .checked_mul(4)?
727        .checked_mul(std::mem::size_of::<usize>())?;
728    let error_bytes = rank
729        .checked_add(1)?
730        .checked_mul(std::mem::size_of::<f64>())?;
731    scalar_bytes
732        .checked_add(index_bytes)?
733        .checked_add(error_bytes)
734}
735
736/// Factorize a lazily supplied matrix with MatrixLUCI block-rook search.
737///
738/// The caller provides a block-fill closure that receives row and column
739/// index lists and writes the corresponding matrix block in column-major
740/// order.
741/// Exactly zero starting residual columns are skipped. An all-zero residual
742/// can require inspecting all remaining entries, without materializing the
743/// full candidate matrix.
744///
745/// # Arguments
746///
747/// * `nrows` - Number of matrix rows.
748/// * `ncols` - Number of matrix columns.
749/// * `fill_block` - Closure that fills `out` with `A[rows, cols]` in
750///
751///   column-major order.
752/// * `options` - Rank and tolerance controls.
753///
754/// # Returns
755///
756/// A [`MatrixLuciFactors`] value containing the pivot metadata and
757/// factors.
758///
759/// # Errors
760///
761/// Returns [`MatrixCIError::InvalidArgument`] for negative or non-finite
762/// tolerances; returns a [`MatrixCIError`] if the lazy factorization fails or
763/// if the block callback is inconsistent.
764///
765/// # Examples
766///
767/// ```
768/// use tensor4all_core::matrix_luci_factors_from_blocks;
769/// use tensor4all_core::RrLUOptions;
770///
771/// let factors = matrix_luci_factors_from_blocks(
772///     2,
773///     2,
774///     |rows, cols, out| {
775///         for (j, &col) in cols.iter().enumerate() {
776///             for (i, &row) in rows.iter().enumerate() {
777///                 out[i + rows.len() * j] = if row == col { 1.0 } else { 0.0 };
778///             }
779///         }
780///     },
781///     RrLUOptions::default(),
782/// ).unwrap();
783/// assert_eq!(factors.rank, 2);
784/// ```
785pub fn matrix_luci_factors_from_blocks<T, F>(
786    nrows: usize,
787    ncols: usize,
788    fill_block: F,
789    options: RrLUOptions,
790) -> Result<MatrixLuciFactors<T>>
791where
792    T: Scalar + crate::MatrixLuciScalar,
793    F: Fn(&[usize], &[usize], &mut [T]),
794{
795    <T as crate::MatrixLuciScalar>::matrix_luci_factors_from_blocks(
796        nrows, ncols, fill_block, options,
797    )
798}
799
800impl<T> MatrixLUCI<T>
801where
802    T: Scalar + crate::MatrixLuciScalar,
803{
804    /// Create a MatrixLUCI from a dense matrix.
805    ///
806    /// # Errors
807    ///
808    /// Returns [`MatrixCIError`] when `a` is not a valid matrix for LU cross
809    /// interpolation or when the factor construction fails.
810    ///
811    /// # Examples
812    ///
813    /// ```
814    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
815    /// use tensor4all_tensorbackend::from_vec2d;
816    ///
817    /// let m = from_vec2d(vec![
818    ///     vec![2.0_f64, 0.0],
819    ///     vec![0.0, 3.0],
820    /// ]);
821    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
822    /// assert!(ci.rank() >= 1);
823    /// ```
824    pub fn from_matrix(a: &Matrix<T>, options: Option<RrLUOptions>) -> Result<Self> {
825        let factors = matrix_luci_factors_from_matrix(a, options)?;
826
827        Ok(Self {
828            nrows: a.nrows(),
829            ncols: a.ncols(),
830            row_indices: factors.row_indices,
831            col_indices: factors.col_indices,
832            left: factors.left,
833            right: factors.right,
834            pivot_errors: factors.pivot_errors,
835        })
836    }
837
838    /// Left CI factor (shape: `nrows x rank`).
839    ///
840    /// The approximation is `left * right`.
841    ///
842    /// # Examples
843    ///
844    /// ```
845    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
846    /// use tensor4all_tensorbackend::{from_vec2d, mat_mul};
847    ///
848    /// let m = from_vec2d(vec![
849    ///     vec![1.0_f64, 2.0],
850    ///     vec![3.0, 4.0],
851    /// ]);
852    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
853    /// let reconstructed = mat_mul(&ci.left(), &ci.right()).unwrap();
854    /// for i in 0..2 {
855    ///     for j in 0..2 {
856    ///         assert!((reconstructed[[i, j]] - m[[i, j]]).abs() < 1e-10);
857    ///     }
858    /// }
859    /// ```
860    pub fn left(&self) -> Matrix<T> {
861        self.left.clone()
862    }
863
864    /// Right CI factor (shape: `rank x ncols`).
865    ///
866    /// The approximation is `left * right`.
867    ///
868    /// # Examples
869    ///
870    /// ```
871    /// use tensor4all_core::{AbstractMatrixCI, MatrixLUCI};
872    /// use tensor4all_tensorbackend::{from_vec2d, mat_mul};
873    ///
874    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
875    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
876    /// let r = ci.right();
877    /// assert_eq!(r.nrows(), ci.rank());
878    /// assert_eq!(r.ncols(), ci.ncols());
879    /// // left * right reconstructs the matrix
880    /// let recon = mat_mul(&ci.left(), &r).unwrap();
881    /// for i in 0..2 {
882    ///     for j in 0..2 {
883    ///         assert!((recon[[i, j]] - m[[i, j]]).abs() < 1e-10);
884    ///     }
885    /// }
886    /// ```
887    pub fn right(&self) -> Matrix<T> {
888        self.right.clone()
889    }
890
891    /// Pivot error history (one entry per pivot, plus a final residual estimate).
892    ///
893    /// # Examples
894    ///
895    /// ```
896    /// use tensor4all_core::MatrixLUCI;
897    /// use tensor4all_tensorbackend::from_vec2d;
898    ///
899    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
900    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
901    /// let errs = ci.pivot_errors();
902    /// assert!(!errs.is_empty());
903    /// // All errors are non-negative
904    /// for &e in &errs {
905    ///     assert!(e >= 0.0);
906    /// }
907    /// ```
908    pub fn pivot_errors(&self) -> Vec<f64> {
909        self.pivot_errors.clone()
910    }
911
912    /// Last pivot error (the residual estimate after all pivots).
913    ///
914    /// # Examples
915    ///
916    /// ```
917    /// use tensor4all_core::MatrixLUCI;
918    /// use tensor4all_tensorbackend::from_vec2d;
919    ///
920    /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
921    /// let ci = MatrixLUCI::from_matrix(&m, None).unwrap();
922    /// let err = ci.last_pivot_error();
923    /// assert!(err >= 0.0);
924    /// ```
925    pub fn last_pivot_error(&self) -> f64 {
926        self.pivot_errors.last().copied().unwrap_or(0.0)
927    }
928}
929
930impl<T> AbstractMatrixCI<T> for MatrixLUCI<T>
931where
932    T: Scalar + crate::MatrixLuciScalar,
933{
934    fn nrows(&self) -> usize {
935        self.nrows
936    }
937
938    fn ncols(&self) -> usize {
939        self.ncols
940    }
941
942    fn rank(&self) -> usize {
943        self.row_indices.len()
944    }
945
946    fn row_indices(&self) -> &[usize] {
947        &self.row_indices
948    }
949
950    fn col_indices(&self) -> &[usize] {
951        &self.col_indices
952    }
953
954    fn evaluate(&self, i: usize, j: usize) -> T {
955        let mut sum = T::zero();
956        for k in 0..self.rank() {
957            sum = sum + self.left[[i, k]] * self.right[[k, j]];
958        }
959        sum
960    }
961
962    fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
963        let r = self.rank();
964        let mut result = Matrix::zeros(rows.len(), cols.len());
965        for (j_out, &col) in cols.iter().enumerate() {
966            for (i_out, &row) in rows.iter().enumerate() {
967                let mut sum = T::zero();
968                for k in 0..r {
969                    sum = sum + self.left[[row, k]] * self.right[[k, col]];
970                }
971                result[[i_out, j_out]] = sum;
972            }
973        }
974        result
975    }
976}
977
978#[cfg(test)]
979mod tests;
tensor4all_core/matrix_luci.rs

tensor4all_core/
matrix_luci.rs
