From f0e9163d700030a816666f48d29f3958c13ac43d Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Thu, 7 Nov 2024 13:57:48 +0000 Subject: [PATCH 01/29] Interpolative decomposition and null space for column and row spaces implemented --- examples/rlst_interpolative_decomposition.rs | 75 ++++++ examples/rlst_null_space.rs | 15 ++ src/dense/linalg.rs | 9 +- .../linalg/interpolative_decomposition.rs | 234 ++++++++++++++++++ src/dense/linalg/null_space.rs | 189 ++++++++++++++ src/prelude.rs | 2 + 6 files changed, 521 insertions(+), 3 deletions(-) create mode 100644 examples/rlst_interpolative_decomposition.rs create mode 100644 examples/rlst_null_space.rs create mode 100644 src/dense/linalg/interpolative_decomposition.rs create mode 100644 src/dense/linalg/null_space.rs diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs new file mode 100644 index 00000000..8756e3ef --- /dev/null +++ b/examples/rlst_interpolative_decomposition.rs @@ -0,0 +1,75 @@ +//! Demo the inverse of a matrix + +pub use rlst::prelude::*; + +//Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere. +fn low_rank_matrix(n: usize, arr: &mut Array, 2>, 2>){ + //Obtain n equally distributed angles 0 LinAlg +impl LinAlg for T { } diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs new file mode 100644 index 00000000..5850b5f9 --- /dev/null +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -0,0 +1,234 @@ +use crate::dense::array::Array; +use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, UnsafeRandomAccessByRef, MultIntoResize, DefaultIterator}; +use crate::dense::types::{RlstResult, RlstScalar}; +use crate::dense::traits::accessors::RandomAccessMut; +use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, empty_array, rlst_dynamic_array2}; +use crate::dense::types::{c32, c64}; +use num::One; + +/// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. +/// +/// The matrix interpolative decomposition is defined for a two dimensional 'long' array `arr` of +/// shape `[m, n]`, where `n>m`. +/// +/// # Example +/// +/// The following command computes the interpolative decomposition of an array `a` for a given tolerance, tol. +/// ``` +/// # use rlst::rlst_dynamic_array2; +/// # let tol: f64 = 1e-5; +/// # let mut a = rlst_dynamic_array2!(f64, [50, 100]); +/// # a.fill_from_seed_equally_distributed(0); +/// # let res = a.view_mut().into_id_alloc(tol, None).unwrap(); +/// ``` + +/// This method allocates memory for the interpolative decomposition. +pub trait MatrixId: RlstScalar { + ///This method allocates space for ID + fn into_id_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + >( + arr: Array, tol: ::Real, k: Option + ) -> RlstResult>; +} + +macro_rules! implement_into_id { + ($scalar:ty) => { + impl MatrixId for $scalar { + fn into_id_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + >( + arr: Array, tol: ::Real, k: Option + ) -> RlstResult> { + IdDecomposition::<$scalar, ArrayImpl>::new(arr, tol, k) + } + } + }; +} + +implement_into_id!(f32); +implement_into_id!(f64); +implement_into_id!(c32); +implement_into_id!(c64); + +impl< + Item: RlstScalar + MatrixId, + ArrayImplId: UnsafeRandomAccessByValue<2, Item = Item> + + Stride<2> + + RawAccessMut + + Shape<2>, + > Array +{ + /// Compute the interpolative decomposition of a given 2-dimensional array. + pub fn into_id_alloc(self, tol: ::Real, k: Option) -> RlstResult> { + ::into_id_alloc(self, tol, k) + } +} + +/// Compute the matrix interpolative decomposition +pub trait MatrixIdDecomposition: Sized { + /// Item type + type Item: RlstScalar; + /// Array implementaion + type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + RawAccessMut + + Shape<2>; + + /// Create a new Interpolative Decomposition from a given array. + fn new(arr: Array, tol: ::Real, k: Option) -> RlstResult; + + ///Compute the permutation matrix associated to the Interpolative Decomposition + fn get_p< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + UnsafeRandomAccessByRef<2, Item = Self::Item>, + >( + arr: Array, + perm: Vec + ); + + +} + +///Stores the relevant features regarding interpolative decomposition. +pub struct IdDecomposition< + Item: RlstScalar, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, +> { + /// arr: permuted array + pub arr: Array, + /// perm_mat: permutation matrix associated to the interpolative decomposition + pub perm_mat: Array, 2>, 2>, + /// rank: rank of the matrix associated to the interpolative decomposition for a given tolerance + pub rank: usize, + ///id_mat: interpolative matrix calculated for a given tolerance + pub id_mat: Array, 2>, 2>, +} + + +macro_rules! impl_id { + ($scalar:ty) => { + impl< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + > MatrixIdDecomposition for IdDecomposition<$scalar, ArrayImpl> + { + type Item = $scalar; + type ArrayImpl = ArrayImpl; + + fn new(mut arr: Array<$scalar, ArrayImpl, 2>, tol: <$scalar as RlstScalar>::Real, k: Option) -> RlstResult{ + //We assume that for a matrix of m rows and n columns, n>m, so we apply ID the transpose + let mut arr_trans: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); + arr_trans.fill_from(arr.view().conj().transpose()); + + //We compute the QR decomposition using rlst QR decomposition + let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); + arr_qr.fill_from(arr.view().conj().transpose()); + let arr_qr_shape = arr_qr.shape(); + let qr = arr_qr.view_mut().into_qr_alloc().unwrap(); + + //We obtain relevant parameters of the decomposition: the permutation induced by the pivoting and the R matrix + let perm = qr.get_perm(); + let mut r_mat = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); + qr.get_r(r_mat.view_mut()); + + + //The maximum rank is given by the number of columns of the transposed matrix + let dim: usize = arr_qr_shape[1]; + let rank: usize; + + //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. + match k { + Some(k) => rank = k, + None => { + //We extract the diagonal to calculate the rank of the matrix + let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.view_mut()); + let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); + + //We compute the rank of the matrix + if max.re() > 0.0{ + let alpha: $scalar = (1.0/max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag.iter().filter(|el| el.abs() > tol.into() ).collect::>(); + rank = aux_vec.len(); + } + else{ + rank = dim; + } + }, + } + + let mut perm_mat = rlst_dynamic_array2!($scalar, [dim, dim]); + Self::get_p(perm_mat.view_mut(), perm); + + + let mut perm_arr = empty_array::<$scalar, 2>() + .simple_mult_into_resize(perm_mat.view_mut(), arr.view()); + + + for col in 0..arr.shape()[1]{ + for row in 0..arr.shape()[0]{ + arr.data_mut()[col*arr_qr_shape[1] + row] = *perm_arr.get_mut([row, col]).unwrap() + } + } + + + //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. + //If not, compute the Interpolative Decomposition matrix. + if rank == 0 || rank >= dim{ + let mut id_mat = rlst_dynamic_array2!($scalar, [dim, dim]); + id_mat.set_identity(); + Ok(Self{arr, perm_mat, rank, id_mat}) + } + else{ + + let shape: [usize; 2] = r_mat.shape(); + let mut id_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>= rlst_dynamic_array2!($scalar, [dim-rank, rank]); + + let mut k11: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [rank, rank]); + k11.fill_from(r_mat.view_mut().into_subview([0, 0], [rank, rank])); + k11.view_mut().into_inverse_alloc().unwrap(); + let mut k12: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [rank, dim-rank]); + k12.fill_from(r_mat.view_mut().into_subview([0, rank], [rank, shape[1]-rank])); + + let res = empty_array().simple_mult_into_resize(k11.view(), k12.view()); + id_mat.fill_from(res.view().conj().transpose().view()); + Ok(Self{arr, perm_mat, rank, id_mat}) + } + } + + fn get_p< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = $scalar> + + UnsafeRandomAccessByRef<2, Item = $scalar>, + >( + mut arr: Array<$scalar, ArrayImplMut, 2>, + perm: Vec + ) { + let m = arr.shape()[0]; + + arr.set_zero(); + for col in 0..m { + arr[[col, perm[col]]] = <$scalar as One>::one(); + } + } + } + }; +} + +impl_id!(f64); +impl_id!(f32); +impl_id!(c32); +impl_id!(c64); diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs new file mode 100644 index 00000000..8a0b0d22 --- /dev/null +++ b/src/dense/linalg/null_space.rs @@ -0,0 +1,189 @@ +use crate::dense::array::Array; +use crate::dense::traits::{ + RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, +}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; +use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, rlst_dynamic_array2}; +use crate::dense::traits::DefaultIterator; +use crate::empty_array; +use itertools::min; + +/// Compute the matrix nullspace. +/// +/// The matrix nullspace is defined for a two dimensional array `arr` of +/// shape `[m, n]`. +/// +/// # Example +/// +/// The following command computes the nullspace of an array `a`. +/// The nullspace is found in +/// ``` +/// # use rlst::rlst_dynamic_array2; +/// # use rlst::dense::linalg::null_space::{NullSpaceType, MatrixNull}; +/// # let mut a = rlst_dynamic_array2!(f64, [3, 4]); +/// # a.fill_from_seed_equally_distributed(0); +/// # let null_res = a.view_mut().into_null_alloc(NullSpaceType::Row).unwrap(); +/// ``` + +/// This method allocates memory for the nullspace computation. + +// pub trait MatrixNull { +// /// Compute the matrix null space +// fn into_null_alloc(arr, null_space_type) -> RlstResult>; +// } + +pub trait MatrixNull: RlstScalar { + /// Compute the matrix null space + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + null_space_type: NullSpaceType + ) -> RlstResult>; +} + +macro_rules! implement_into_null { + ($scalar:ty) => { + impl MatrixNull for $scalar { + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + null_space_type: NullSpaceType + ) -> RlstResult> { + NullSpace::<$scalar>::new(arr, null_space_type) + } + } + }; +} + +implement_into_null!(f32); +implement_into_null!(f64); +implement_into_null!(c32); +implement_into_null!(c64); + +impl< + Item: RlstScalar + MatrixNull, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + + Stride<2> + + RawAccessMut + + Shape<2>, + > Array +{ + /// Compute the Column or Row nullspace of a given 2-dimensional array. + pub fn into_null_alloc(self, null_space_type: NullSpaceType) -> RlstResult> { + ::into_null_alloc(self, null_space_type) + } +} + +///Null space +#[derive(Clone, Copy)] +#[repr(u8)] +pub enum NullSpaceType { + /// Row Nullspace + Row = b'R', + /// Column Nullspace + Column = b'C', +} + +/// QR decomposition +pub struct NullSpace< + Item: RlstScalar +> { + ///Row or column nullspace + pub null_space_type: NullSpaceType, + ///Computed null space + pub null_space_arr: Array, 2>, 2>, +} + +macro_rules! implement_null_space { + ($scalar:ty) => { + impl NullSpace<$scalar> + { + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + >(arr: Array<$scalar, ArrayImpl, 2>, null_space_type: NullSpaceType) -> RlstResult { + + let shape = arr.shape(); + + match null_space_type { + NullSpaceType::Row => { + let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); + arr_qr.fill_from(arr.view().conj().transpose()); + let mut null_space_arr = empty_array(); + Self::find_null_space(arr_qr, &mut null_space_arr); + Ok(Self {null_space_type, null_space_arr}) + + }, + NullSpaceType::Column => { + let mut null_space_arr = empty_array(); + Self::find_null_space(arr, &mut null_space_arr); + Ok(Self {null_space_type, null_space_arr}) + }, + } + + } + + fn find_null_space< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + >(arr: Array<$scalar, ArrayImpl, 2>, null_space_arr: &mut Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>) + { + let shape = arr.shape(); + let dim: usize = min(shape).unwrap(); + let qr = arr.into_qr_alloc().unwrap(); + + //We compute the QR decomposition to find a linearly independent basis of the space. + let mut q: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); + let _ = qr.get_q_alloc(q.view_mut()); + + let mut r_mat = rlst_dynamic_array2!($scalar, [dim, dim]); + qr.get_r(r_mat.view_mut()); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = Self::find_matrix_rank(r_mat, dim); + + null_space_arr.fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0]-rank])); + } + + fn find_matrix_rank(r_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>, dim: usize)->usize{ + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.view_mut()); + + let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); + let rank: usize; + + if max.re() > 0.0{ + let alpha: $scalar = (1.0/max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag.iter().filter(|el| el.abs() > 1e-15 ).collect::>(); + rank = aux_vec.len(); + } + else{ + rank = dim; + } + rank + + } + } + }; +} + +implement_null_space!(f64); +implement_null_space!(f32); +implement_null_space!(c64); +implement_null_space!(c32); + diff --git a/src/prelude.rs b/src/prelude.rs index 529dad7e..f1da0c8f 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -64,6 +64,8 @@ pub use crate::dense::linalg::lu::{LuDecomposition, MatrixLuDecomposition}; pub use crate::dense::linalg::pseudo_inverse::MatrixPseudoInverse; pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; +pub use crate::dense::linalg::interpolative_decomposition::{MatrixId, IdDecomposition}; +pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; pub use crate::dense::array::rank1_array::Rank1Array; From c042a1228625bc510747552d07e208744396b056 Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Thu, 14 Nov 2024 20:01:33 +0000 Subject: [PATCH 02/29] Elementary matrices implemented --- examples/rlst_elementary_matrices.rs | 69 +++++++ src/dense/linalg.rs | 7 +- src/dense/linalg/elementary_matrix.rs | 255 ++++++++++++++++++++++++++ src/doc/dense_linear_algebra.rs | 35 ++++ src/prelude.rs | 1 + 5 files changed, 364 insertions(+), 3 deletions(-) create mode 100644 examples/rlst_elementary_matrices.rs create mode 100644 src/dense/linalg/elementary_matrix.rs diff --git a/examples/rlst_elementary_matrices.rs b/examples/rlst_elementary_matrices.rs new file mode 100644 index 00000000..1d7610ba --- /dev/null +++ b/examples/rlst_elementary_matrices.rs @@ -0,0 +1,69 @@ +//! Demo the inverse of a matrix +use rlst::dense::linalg::elementary_matrix::ElementaryOperations; +pub use rlst::prelude::*; +use rlst::dense::linalg::elementary_matrix::{OpType, RowOpType}; + +pub fn main() { + //Example 1: use elementary matrices to perform one step of an LU block decomposition + //Here we use a 9x9 matrix and we partition it in 3x3 blocks + let mut arr = rlst_dynamic_array2!(f64, [9, 9]); + arr.fill_from_seed_equally_distributed(0); + + //We use the block 11 to eliminate block 21 using the LU with elementary matrices + let mut block_11 = rlst_dynamic_array2!(f64, [3, 3]); + let mut block_21 = rlst_dynamic_array2!(f64, [3, 3]); + block_11.fill_from(arr.view().into_subview([0, 0], [3, 3])); + block_21.fill_from(arr.view().into_subview([3, 0], [3, 3])); + block_11.view_mut().into_inverse_alloc().unwrap(); + + //We define the elementary matrix using the coordinates of the square blocks block_11 and block_21 and the dimension of the elementary matrix + //(in this case it's a 9x9 matrix and we assume that elementary matrices are squared) + let mut el_mat_fact = empty_array().simple_mult_into_resize(block_21.view(), block_11.view()); + let row_indices: Vec = (3..6).collect(); + let col_indices: Vec = (0..3).collect(); + let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Row, false).unwrap(); + el_mat.el_mul(arr.view_mut(), Some(RowOpType::Sub), None).unwrap(); + + //As a result we se that the block 21 has been eliminated + let res = arr.view().into_subview([3, 0], [3, 3]); + println!("L2 of block 21 after LU elimination through elementary matrices {}", res.view_flat().norm_2()); + + + ////////////////////////////////////////////////////////////////////////////////// + //Example 2: Use elementary matrices for row scaling + //Using the matrix of the previous example, we can re-scale the rows corresponding to blocks 11, 12, 13: + let row_indices: Vec = (0..3).collect(); //The row and col indices must be the same, since the relevant dimension is the row_indices + let col_indices: Vec = (0..3).collect(); + + //We simply use a 1x1 matrix for this case, since the only relevant information is contained in the scaling factor + let mut el_mat_fact = rlst_dynamic_array2!(f64, [1, 1]); + let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Mul, false).unwrap(); + + //We extract the elements of the matrix before scaling: + let mut bef_scaling = rlst_dynamic_array2!(f64, [3, 9]); + bef_scaling.fill_from(arr.view().into_subview([0, 0], [3, 9])); + el_mat.el_mul(arr.view_mut(), None, Some(5.0)).unwrap(); + let aft_scaling = arr.view().into_subview([0, 0], [3, 9]); + let norm_comp = aft_scaling.view_flat().norm_2()/bef_scaling.view_flat().norm_2(); + println!("The norm of the rows corresponding to blocks 11, 12, 13 after scaling is {} times larger", norm_comp); + + + ////////////////////////////////////////////////////////////////////////////////// + //Example 3: Use elementary matrices for row permutation + // The permutation is as follows col_indices->row_indices. In this case: + // row 0 moves into row 2, row 1 moves into row 0 and row 2 moves into row 1. + let col_indices: Vec = [0, 1, 2].to_vec(); + let row_indices: Vec = [2, 0, 1].to_vec(); + + let mut el_mat_fact = rlst_dynamic_array2!(f64, [1, 1]); + let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Perm, false).unwrap(); + let mut bef_perm = rlst_dynamic_array2!(f64, [3, 9]); + bef_perm.fill_from(arr.view().into_subview([0, 0], [3, 9])); + el_mat.el_mul(arr.view_mut(), None, None).unwrap(); + let aft_perm = arr.view().into_subview([0, 0], [3, 9]); + + //Here we check if row 0 moved into row 2: + let res = bef_perm.into_subview([0, 0], [1, 9]) - aft_perm.into_subview([2, 0], [1, 9]); + + println!("The difference between the row 0 before permutation and the row 2 after the permutation is {}",res.view_flat().norm_2()); +} diff --git a/src/dense/linalg.rs b/src/dense/linalg.rs index b8327b1f..11f3979b 100644 --- a/src/dense/linalg.rs +++ b/src/dense/linalg.rs @@ -1,7 +1,7 @@ //! Linear algebra routines -use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar}; +use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar, ElementaryMatrixData}; use self::lu::MatrixLu; @@ -11,6 +11,7 @@ pub mod pseudo_inverse; pub mod qr; pub mod svd; pub mod interpolative_decomposition; +pub mod elementary_matrix; pub mod null_space; /// Return true if stride is column major as required by Lapack. @@ -23,9 +24,9 @@ pub fn assert_lapack_stride(stride: [usize; 2]) { } /// Marker trait for objects that support Matrix decompositions. -pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull{} +pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull +ElementaryMatrixData{} -impl LinAlg +impl LinAlg for T { } diff --git a/src/dense/linalg/elementary_matrix.rs b/src/dense/linalg/elementary_matrix.rs new file mode 100644 index 00000000..6d3dd1a3 --- /dev/null +++ b/src/dense/linalg/elementary_matrix.rs @@ -0,0 +1,255 @@ +//! Elementary matrices (row swapping, row multiplication and row addition) +use crate::dense::array::Array; +use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, MultIntoResize}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar, RlstError}; +use crate::{empty_array, rlst_dynamic_array2}; +use crate::dense::traits::accessors::RandomAccessMut; + +///Allocate space for the elementary matrix +pub trait ElementaryMatrixData: RlstScalar { + ///This is the method that allocates space for the elements of the elementary matrix + fn into_el_mat_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool + ) -> RlstResult>; +} + +macro_rules! implement_into_el_mat { + ($scalar:ty) => { + impl ElementaryMatrixData for $scalar { + fn into_el_mat_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool + ) -> RlstResult> { + ElementaryMatrix::<$scalar, ArrayImpl>::new(arr, dim, row_indices, col_indices, op_type, trans) + } + } + }; +} + +implement_into_el_mat!(f32); +implement_into_el_mat!(f64); +implement_into_el_mat!(c32); +implement_into_el_mat!(c64); + +impl< + Item: RlstScalar + ElementaryMatrixData, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + + Stride<2> + + RawAccessMut + + Shape<2>, + > Array +{ + ///This computes an elementary (row addition/substraction, row scaling and row permutation) matrix given a series of parameters + pub fn into_el_mat_alloc(self, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult> { + ::into_el_mat_alloc(self, dim, row_indices, col_indices, op_type, trans) + } +} + +///These are the methods associated to elementary matrices +pub trait ElementaryOperations: Sized { + /// Item type + type Item: RlstScalar; + /// Array implementaion + type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + RawAccessMut + + Shape<2>; + + /// We create the Elementary matrix (E), so it can be applied to a matrix A. In other words, we implement E(A). + /// arr: is needed when row additions/substractions on A are performed. + /// dim: is the dimension of the elementary matrix, given by the row dimension of A. + /// row_indices and col_indices are respectively the domain and range of this application; i.e. we take the rows of A + /// corresponding to col_indices of A and we place the result of this operation in row_indices of A. + /// op_type: indicates if we are adding/substracting rows of A, scaling rows of A by a scalar or if we are permuting rows of A. + /// trans: indicates if we are applying E^T + fn new(arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult; + + ///This method performs E(A). Here: + /// right_arr: matrix A. + /// row_op_type: indicates substraction or addition of rows + /// alpha: is the scaling parameter of a scaling is applied + fn el_mul(self, right_arr: Array, row_op_type: Option, alpha: Option)-> RlstResult<()> ; + + ///This method implements the row addition/substraction + fn row_ops(self, right_arr: Array, op_type: RowOpType); + + ///This method implements the row scaling + fn row_mul(self, right_arr: Array, alpha: Self::Item); + + ///This method implements the row permutation + fn row_perm(self, right_arr: Array); + +} + +/// Container for the LU Decomposition of a matrix. +pub struct ElementaryMatrix< + Item: RlstScalar, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, +> { + dim: usize, + arr: Array, + row_indices: Vec, + col_indices: Vec, + op_type: OpType, + trans: bool + +} + +/// Type of row operation +#[derive(Clone, Copy)] +#[repr(u8)] +pub enum RowOpType { + /// Row addition + Add = b'A', + /// Row substraction + Sub = b'S', +} + +/// Type of Elementary matrix +#[derive(Clone, Copy)] +#[repr(u8)] +pub enum OpType { + /// Row operation (addition or substraction) + Row = b'R', + /// Row scaling + Mul = b'M', + /// Row permutation + Perm = b'P', +} + + +macro_rules! impl_el_ops { + ($scalar:ty) => { + impl< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + > ElementaryOperations for ElementaryMatrix<$scalar, ArrayImpl> + { + type Item = $scalar; + type ArrayImpl = ArrayImpl; + + fn new(arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult { + + match op_type{ + OpType::Row => {}, + OpType::Mul => {assert_eq!(row_indices.len(), col_indices.len())}, + OpType::Perm => {assert_eq!(row_indices.len(), col_indices.len())} + } + + if trans{ + Ok(Self{arr, dim, row_indices: col_indices, col_indices: row_indices, op_type, trans}) + } + else{ + Ok(Self{arr, dim, row_indices, col_indices, op_type, trans}) + } + } + + fn el_mul(self, right_arr: Array, row_op_type: Option, alpha: Option)-> RlstResult<()> + { + match self.op_type{ + OpType::Row => { + match(row_op_type){ + Some(row_op_type)=> Ok(self.row_ops(right_arr, row_op_type)), + None =>Err(RlstError::IoError("Missing parameter".to_string())) + }}, + OpType::Mul => { + match(alpha){ + Some(alpha) => Ok(self.row_mul(right_arr, alpha)), + None =>Err(RlstError::IoError("Missing parameter".to_string())) + }}, + OpType::Perm => {Ok(self.row_perm(right_arr))} + } + } + + fn row_ops(self, mut right_arr: Array, row_op_type: RowOpType){ + let mut subarr_cols = rlst_dynamic_array2!($scalar, [self.col_indices.len(), self.dim]); + let mut subarr_rows = rlst_dynamic_array2!($scalar, [self.row_indices.len(), self.dim]); + let mut aux_arr = empty_array(); + + if self.trans{ + aux_arr.fill_from_resize(self.arr.view().conj().transpose()); + } + else{ + aux_arr.fill_from_resize(self.arr.view()); + } + + let right_arr_shape = right_arr.view().shape(); + let dim = self.dim; + let col_indices = self.col_indices; + let row_indices = self.row_indices; + + for col in 0..dim{ + for row in 0..col_indices.len(){ + *subarr_cols.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + col_indices[row]]; + } + for row in 0..row_indices.len(){ + *subarr_rows.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]]; + } + } + + let res_mul = empty_array::<$scalar, 2>().simple_mult_into_resize(aux_arr, subarr_cols.view_mut()); + let mut add_res = rlst_dynamic_array2!($scalar, subarr_rows.shape()); + + match row_op_type{ + RowOpType::Add => {add_res.fill_from(subarr_rows.view() + res_mul.view());}, + RowOpType::Sub => {add_res.fill_from(subarr_rows.view() - res_mul.view());} + } + + for col in 0..dim{ + for row in 0..row_indices.len(){ + right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = *add_res.get_mut([row, col]).unwrap(); + } + } + } + + fn row_mul(self, mut right_arr: Array, alpha: Self::Item){ + let right_arr_shape = right_arr.view().shape(); + let dim = self.dim; + let row_indices = self.row_indices; + + for col in 0..dim{ + for row in 0..row_indices.len(){ + right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = alpha*right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] + } + } + + } + + fn row_perm(self, mut right_arr: Array){ + let dim = self.dim; + let col_indices = self.col_indices; + let row_indices = self.row_indices; + let right_arr_shape = right_arr.view().shape(); + let mut subarr_cols = rlst_dynamic_array2!($scalar, [col_indices.len(), dim]); + + for col in 0..dim{ + for row in 0..col_indices.len(){ + *subarr_cols.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + col_indices[row]]; + } + } + + for col in 0..dim{ + for row in 0..row_indices.len(){ + right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = *subarr_cols.get_mut([row, col]).unwrap(); + } + } + } + } + }; +} + +impl_el_ops!(f64); +impl_el_ops!(f32); +impl_el_ops!(c64); +impl_el_ops!(c32); diff --git a/src/doc/dense_linear_algebra.rs b/src/doc/dense_linear_algebra.rs index cb1eb41c..827caa71 100644 --- a/src/doc/dense_linear_algebra.rs +++ b/src/doc/dense_linear_algebra.rs @@ -514,6 +514,41 @@ //! ``` //! To compute the full SVD use the parameter [SvdMode::Full](crate::SvdMode::Full). //! +//! # Interpolative decomposition +//! +//! To compute the Interpolative Decomposition of a two-dimensional array (long or square) `arr`, for a given tolerance, use +//! +//! ``` +//! # use rlst::prelude::*; +//! let mut rand = rand::thread_rng(); +//! let mut arr = rlst_dynamic_array2!(f64, [5, 8]); +//! let tol: f64 = 1e-5; +//! let mut res = arr.view_mut().into_id_alloc(tol, None).unwrap(); +//! ``` +//! +//! You can take a look at the interpolation matrix using +//! +//! ``` +//! res.id_mat.pretty_print(); +//! ``` +//! +//! With this, the approximation of the low-rank elements of `arr` is given by +//! +//! ``` +//! let arr_app = empty_array().simple_mult_into_resize(res.id_mat.view(), res.arr.view_mut().into_subview([0, 0], [res.rank, n])); +//! ``` +//! +//! You can also obtain the Interpolative Decomposition by giving the algorithm a fixed rank, and the tolerance parameter is ignored, like so +//! +//! ``` +//! # use rlst::prelude::*; +//! let mut rand = rand::thread_rng(); +//! let mut arr = rlst_dynamic_array2!(f64, [5, 8]); +//! let tol: f64 = 1e-5; +//! k = 2; +//! let mut res = arr.view_mut().into_id_alloc(tol, k).unwrap(); +//! ``` +//! //! # Other vector functions //! //! The following other functions for arrays with a single dimension (i.e. vectors) are provided. diff --git a/src/prelude.rs b/src/prelude.rs index f1da0c8f..b759a4cc 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -66,6 +66,7 @@ pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; pub use crate::dense::linalg::interpolative_decomposition::{MatrixId, IdDecomposition}; pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; +pub use crate::dense::linalg::elementary_matrix::{ElementaryMatrixData, ElementaryMatrix}; pub use crate::dense::array::rank1_array::Rank1Array; From a8045100947746e2055ba147f35fe82ee64ea474 Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Fri, 15 Nov 2024 14:34:57 +0000 Subject: [PATCH 03/29] updates to transposition of elementary matrix --- src/dense/linalg/elementary_matrix.rs | 34 ++++++++++++++++++--------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/src/dense/linalg/elementary_matrix.rs b/src/dense/linalg/elementary_matrix.rs index 6d3dd1a3..93b446c7 100644 --- a/src/dense/linalg/elementary_matrix.rs +++ b/src/dense/linalg/elementary_matrix.rs @@ -100,7 +100,9 @@ pub struct ElementaryMatrix< row_indices: Vec, col_indices: Vec, op_type: OpType, - trans: bool + + ///Indicates if the operation is to be transposed or not + pub trans: bool } @@ -147,12 +149,7 @@ macro_rules! impl_el_ops { OpType::Perm => {assert_eq!(row_indices.len(), col_indices.len())} } - if trans{ - Ok(Self{arr, dim, row_indices: col_indices, col_indices: row_indices, op_type, trans}) - } - else{ - Ok(Self{arr, dim, row_indices, col_indices, op_type, trans}) - } + Ok(Self{arr, dim, row_indices, col_indices, op_type, trans}) } fn el_mul(self, right_arr: Array, row_op_type: Option, alpha: Option)-> RlstResult<()> @@ -176,18 +173,23 @@ macro_rules! impl_el_ops { let mut subarr_cols = rlst_dynamic_array2!($scalar, [self.col_indices.len(), self.dim]); let mut subarr_rows = rlst_dynamic_array2!($scalar, [self.row_indices.len(), self.dim]); let mut aux_arr = empty_array(); + let col_indices :Vec; + let row_indices :Vec; if self.trans{ aux_arr.fill_from_resize(self.arr.view().conj().transpose()); + col_indices = self.row_indices; + row_indices = self.col_indices; } else{ aux_arr.fill_from_resize(self.arr.view()); + col_indices = self.col_indices; + row_indices = self.row_indices; } let right_arr_shape = right_arr.view().shape(); let dim = self.dim; - let col_indices = self.col_indices; - let row_indices = self.row_indices; + for col in 0..dim{ for row in 0..col_indices.len(){ @@ -228,8 +230,18 @@ macro_rules! impl_el_ops { fn row_perm(self, mut right_arr: Array){ let dim = self.dim; - let col_indices = self.col_indices; - let row_indices = self.row_indices; + let col_indices :Vec; + let row_indices :Vec; + + if self.trans{ + col_indices = self.row_indices; + row_indices = self.col_indices; + } + else{ + col_indices = self.col_indices; + row_indices = self.row_indices; + } + let right_arr_shape = right_arr.view().shape(); let mut subarr_cols = rlst_dynamic_array2!($scalar, [col_indices.len(), dim]); From ae4dfe7daa5d1fa7ddca8c73a0bf83bd53600e43 Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Fri, 22 Nov 2024 12:24:48 +0000 Subject: [PATCH 04/29] fixes to interpolative decomposition and null space --- examples/rlst_elementary_matrices.rs | 69 ----- examples/rlst_interpolative_decomposition.rs | 24 +- examples/rlst_null_space.rs | 7 +- src/dense/linalg.rs | 7 +- src/dense/linalg/elementary_matrix.rs | 267 ------------------ .../linalg/interpolative_decomposition.rs | 182 ++++++------ src/dense/linalg/null_space.rs | 101 ++----- src/prelude.rs | 1 - 8 files changed, 152 insertions(+), 506 deletions(-) delete mode 100644 examples/rlst_elementary_matrices.rs delete mode 100644 src/dense/linalg/elementary_matrix.rs diff --git a/examples/rlst_elementary_matrices.rs b/examples/rlst_elementary_matrices.rs deleted file mode 100644 index 1d7610ba..00000000 --- a/examples/rlst_elementary_matrices.rs +++ /dev/null @@ -1,69 +0,0 @@ -//! Demo the inverse of a matrix -use rlst::dense::linalg::elementary_matrix::ElementaryOperations; -pub use rlst::prelude::*; -use rlst::dense::linalg::elementary_matrix::{OpType, RowOpType}; - -pub fn main() { - //Example 1: use elementary matrices to perform one step of an LU block decomposition - //Here we use a 9x9 matrix and we partition it in 3x3 blocks - let mut arr = rlst_dynamic_array2!(f64, [9, 9]); - arr.fill_from_seed_equally_distributed(0); - - //We use the block 11 to eliminate block 21 using the LU with elementary matrices - let mut block_11 = rlst_dynamic_array2!(f64, [3, 3]); - let mut block_21 = rlst_dynamic_array2!(f64, [3, 3]); - block_11.fill_from(arr.view().into_subview([0, 0], [3, 3])); - block_21.fill_from(arr.view().into_subview([3, 0], [3, 3])); - block_11.view_mut().into_inverse_alloc().unwrap(); - - //We define the elementary matrix using the coordinates of the square blocks block_11 and block_21 and the dimension of the elementary matrix - //(in this case it's a 9x9 matrix and we assume that elementary matrices are squared) - let mut el_mat_fact = empty_array().simple_mult_into_resize(block_21.view(), block_11.view()); - let row_indices: Vec = (3..6).collect(); - let col_indices: Vec = (0..3).collect(); - let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Row, false).unwrap(); - el_mat.el_mul(arr.view_mut(), Some(RowOpType::Sub), None).unwrap(); - - //As a result we se that the block 21 has been eliminated - let res = arr.view().into_subview([3, 0], [3, 3]); - println!("L2 of block 21 after LU elimination through elementary matrices {}", res.view_flat().norm_2()); - - - ////////////////////////////////////////////////////////////////////////////////// - //Example 2: Use elementary matrices for row scaling - //Using the matrix of the previous example, we can re-scale the rows corresponding to blocks 11, 12, 13: - let row_indices: Vec = (0..3).collect(); //The row and col indices must be the same, since the relevant dimension is the row_indices - let col_indices: Vec = (0..3).collect(); - - //We simply use a 1x1 matrix for this case, since the only relevant information is contained in the scaling factor - let mut el_mat_fact = rlst_dynamic_array2!(f64, [1, 1]); - let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Mul, false).unwrap(); - - //We extract the elements of the matrix before scaling: - let mut bef_scaling = rlst_dynamic_array2!(f64, [3, 9]); - bef_scaling.fill_from(arr.view().into_subview([0, 0], [3, 9])); - el_mat.el_mul(arr.view_mut(), None, Some(5.0)).unwrap(); - let aft_scaling = arr.view().into_subview([0, 0], [3, 9]); - let norm_comp = aft_scaling.view_flat().norm_2()/bef_scaling.view_flat().norm_2(); - println!("The norm of the rows corresponding to blocks 11, 12, 13 after scaling is {} times larger", norm_comp); - - - ////////////////////////////////////////////////////////////////////////////////// - //Example 3: Use elementary matrices for row permutation - // The permutation is as follows col_indices->row_indices. In this case: - // row 0 moves into row 2, row 1 moves into row 0 and row 2 moves into row 1. - let col_indices: Vec = [0, 1, 2].to_vec(); - let row_indices: Vec = [2, 0, 1].to_vec(); - - let mut el_mat_fact = rlst_dynamic_array2!(f64, [1, 1]); - let el_mat = el_mat_fact.view_mut().into_el_mat_alloc(9, row_indices, col_indices, OpType::Perm, false).unwrap(); - let mut bef_perm = rlst_dynamic_array2!(f64, [3, 9]); - bef_perm.fill_from(arr.view().into_subview([0, 0], [3, 9])); - el_mat.el_mul(arr.view_mut(), None, None).unwrap(); - let aft_perm = arr.view().into_subview([0, 0], [3, 9]); - - //Here we check if row 0 moved into row 2: - let res = bef_perm.into_subview([0, 0], [1, 9]) - aft_perm.into_subview([2, 0], [1, 9]); - - println!("The difference between the row 0 before permutation and the row 2 after the permutation is {}",res.view_flat().norm_2()); -} diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 8756e3ef..44a106fb 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -1,6 +1,7 @@ //! Demo the inverse of a matrix pub use rlst::prelude::*; +use rlst::dense::linalg::interpolative_decomposition::{MatrixIdDecomposition, Accuracy}; //Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere. fn low_rank_matrix(n: usize, arr: &mut Array, 2>, 2>){ @@ -54,22 +55,31 @@ pub fn main() { let tol: f64 = 1e-5; //Create a low rank matrix - let mut arr = rlst_dynamic_array2!(f64, [slice, n]); - low_rank_matrix(n,&mut arr); + let mut arr: DynamicArray = rlst_dynamic_array2!(f64, [slice, n]); + low_rank_matrix(n, &mut arr); - let mut res = arr.view_mut().into_id_alloc(tol, None).unwrap(); + let res: IdDecomposition = arr.view_mut().into_id_alloc(Accuracy::Tol(tol)).unwrap(); - println!("The permuted matrix is:"); - res.arr.pretty_print(); + println!("The skeleton of the matrix is given by"); + res.skel.pretty_print(); println!("The interpolative decomposition matrix is:"); res.id_mat.pretty_print(); println!("The rank of this matrix is {}\n", res.rank); - let a_rs_app = empty_array().simple_mult_into_resize(res.id_mat.view(), res.arr.view_mut().into_subview([0, 0], [res.rank, n])); - let a_rs = res.arr.view_mut().into_subview([res.rank, 0], [slice-res.rank, n]); + //We extract the residuals of the matrix + let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); + res.get_p(perm_mat.view_mut()); + let perm_arr: DynamicArray = empty_array::() + .simple_mult_into_resize(perm_mat.view_mut(), arr.view()); + let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice-res.rank, n]); + a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice-res.rank, n])); + + //We compute an approximation of the residual columns of the matrix + let a_rs_app: DynamicArray = empty_array().simple_mult_into_resize(res.id_mat.view(), res.skel); + let error: f64 = a_rs.view().sub(a_rs_app.view()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); } diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index 12e29b1b..5c0afd7c 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -1,14 +1,13 @@ //! Demo the inverse of a matrix -use rlst::dense::linalg::null_space::NullSpaceType; pub use rlst::prelude::*; //Here we compute the null space (B) of the rowspace of a short matrix (A). pub fn main() { let mut arr = rlst_dynamic_array2!(f64, [3, 4]); arr.fill_from_seed_equally_distributed(0); - - let null_res = arr.view_mut().into_null_alloc(NullSpaceType::Row).unwrap(); - let res = empty_array().simple_mult_into_resize(arr.view_mut(), null_res.null_space_arr.view()); + let tol = 1e-15; + let null_res = arr.view_mut().into_null_alloc(tol).unwrap(); + let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.view_mut(), null_res.null_space_arr.view()); println!("Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2()); diff --git a/src/dense/linalg.rs b/src/dense/linalg.rs index 11f3979b..b8327b1f 100644 --- a/src/dense/linalg.rs +++ b/src/dense/linalg.rs @@ -1,7 +1,7 @@ //! Linear algebra routines -use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar, ElementaryMatrixData}; +use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar}; use self::lu::MatrixLu; @@ -11,7 +11,6 @@ pub mod pseudo_inverse; pub mod qr; pub mod svd; pub mod interpolative_decomposition; -pub mod elementary_matrix; pub mod null_space; /// Return true if stride is column major as required by Lapack. @@ -24,9 +23,9 @@ pub fn assert_lapack_stride(stride: [usize; 2]) { } /// Marker trait for objects that support Matrix decompositions. -pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull +ElementaryMatrixData{} +pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull{} -impl LinAlg +impl LinAlg for T { } diff --git a/src/dense/linalg/elementary_matrix.rs b/src/dense/linalg/elementary_matrix.rs deleted file mode 100644 index 93b446c7..00000000 --- a/src/dense/linalg/elementary_matrix.rs +++ /dev/null @@ -1,267 +0,0 @@ -//! Elementary matrices (row swapping, row multiplication and row addition) -use crate::dense::array::Array; -use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, MultIntoResize}; -use crate::dense::types::{c32, c64, RlstResult, RlstScalar, RlstError}; -use crate::{empty_array, rlst_dynamic_array2}; -use crate::dense::traits::accessors::RandomAccessMut; - -///Allocate space for the elementary matrix -pub trait ElementaryMatrixData: RlstScalar { - ///This is the method that allocates space for the elements of the elementary matrix - fn into_el_mat_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool - ) -> RlstResult>; -} - -macro_rules! implement_into_el_mat { - ($scalar:ty) => { - impl ElementaryMatrixData for $scalar { - fn into_el_mat_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool - ) -> RlstResult> { - ElementaryMatrix::<$scalar, ArrayImpl>::new(arr, dim, row_indices, col_indices, op_type, trans) - } - } - }; -} - -implement_into_el_mat!(f32); -implement_into_el_mat!(f64); -implement_into_el_mat!(c32); -implement_into_el_mat!(c64); - -impl< - Item: RlstScalar + ElementaryMatrixData, - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> - + Stride<2> - + RawAccessMut - + Shape<2>, - > Array -{ - ///This computes an elementary (row addition/substraction, row scaling and row permutation) matrix given a series of parameters - pub fn into_el_mat_alloc(self, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult> { - ::into_el_mat_alloc(self, dim, row_indices, col_indices, op_type, trans) - } -} - -///These are the methods associated to elementary matrices -pub trait ElementaryOperations: Sized { - /// Item type - type Item: RlstScalar; - /// Array implementaion - type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + RawAccessMut - + Shape<2>; - - /// We create the Elementary matrix (E), so it can be applied to a matrix A. In other words, we implement E(A). - /// arr: is needed when row additions/substractions on A are performed. - /// dim: is the dimension of the elementary matrix, given by the row dimension of A. - /// row_indices and col_indices are respectively the domain and range of this application; i.e. we take the rows of A - /// corresponding to col_indices of A and we place the result of this operation in row_indices of A. - /// op_type: indicates if we are adding/substracting rows of A, scaling rows of A by a scalar or if we are permuting rows of A. - /// trans: indicates if we are applying E^T - fn new(arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult; - - ///This method performs E(A). Here: - /// right_arr: matrix A. - /// row_op_type: indicates substraction or addition of rows - /// alpha: is the scaling parameter of a scaling is applied - fn el_mul(self, right_arr: Array, row_op_type: Option, alpha: Option)-> RlstResult<()> ; - - ///This method implements the row addition/substraction - fn row_ops(self, right_arr: Array, op_type: RowOpType); - - ///This method implements the row scaling - fn row_mul(self, right_arr: Array, alpha: Self::Item); - - ///This method implements the row permutation - fn row_perm(self, right_arr: Array); - -} - -/// Container for the LU Decomposition of a matrix. -pub struct ElementaryMatrix< - Item: RlstScalar, - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, -> { - dim: usize, - arr: Array, - row_indices: Vec, - col_indices: Vec, - op_type: OpType, - - ///Indicates if the operation is to be transposed or not - pub trans: bool - -} - -/// Type of row operation -#[derive(Clone, Copy)] -#[repr(u8)] -pub enum RowOpType { - /// Row addition - Add = b'A', - /// Row substraction - Sub = b'S', -} - -/// Type of Elementary matrix -#[derive(Clone, Copy)] -#[repr(u8)] -pub enum OpType { - /// Row operation (addition or substraction) - Row = b'R', - /// Row scaling - Mul = b'M', - /// Row permutation - Perm = b'P', -} - - -macro_rules! impl_el_ops { - ($scalar:ty) => { - impl< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - > ElementaryOperations for ElementaryMatrix<$scalar, ArrayImpl> - { - type Item = $scalar; - type ArrayImpl = ArrayImpl; - - fn new(arr: Array, dim:usize, row_indices: Vec, col_indices: Vec, op_type: OpType, trans: bool) -> RlstResult { - - match op_type{ - OpType::Row => {}, - OpType::Mul => {assert_eq!(row_indices.len(), col_indices.len())}, - OpType::Perm => {assert_eq!(row_indices.len(), col_indices.len())} - } - - Ok(Self{arr, dim, row_indices, col_indices, op_type, trans}) - } - - fn el_mul(self, right_arr: Array, row_op_type: Option, alpha: Option)-> RlstResult<()> - { - match self.op_type{ - OpType::Row => { - match(row_op_type){ - Some(row_op_type)=> Ok(self.row_ops(right_arr, row_op_type)), - None =>Err(RlstError::IoError("Missing parameter".to_string())) - }}, - OpType::Mul => { - match(alpha){ - Some(alpha) => Ok(self.row_mul(right_arr, alpha)), - None =>Err(RlstError::IoError("Missing parameter".to_string())) - }}, - OpType::Perm => {Ok(self.row_perm(right_arr))} - } - } - - fn row_ops(self, mut right_arr: Array, row_op_type: RowOpType){ - let mut subarr_cols = rlst_dynamic_array2!($scalar, [self.col_indices.len(), self.dim]); - let mut subarr_rows = rlst_dynamic_array2!($scalar, [self.row_indices.len(), self.dim]); - let mut aux_arr = empty_array(); - let col_indices :Vec; - let row_indices :Vec; - - if self.trans{ - aux_arr.fill_from_resize(self.arr.view().conj().transpose()); - col_indices = self.row_indices; - row_indices = self.col_indices; - } - else{ - aux_arr.fill_from_resize(self.arr.view()); - col_indices = self.col_indices; - row_indices = self.row_indices; - } - - let right_arr_shape = right_arr.view().shape(); - let dim = self.dim; - - - for col in 0..dim{ - for row in 0..col_indices.len(){ - *subarr_cols.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + col_indices[row]]; - } - for row in 0..row_indices.len(){ - *subarr_rows.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]]; - } - } - - let res_mul = empty_array::<$scalar, 2>().simple_mult_into_resize(aux_arr, subarr_cols.view_mut()); - let mut add_res = rlst_dynamic_array2!($scalar, subarr_rows.shape()); - - match row_op_type{ - RowOpType::Add => {add_res.fill_from(subarr_rows.view() + res_mul.view());}, - RowOpType::Sub => {add_res.fill_from(subarr_rows.view() - res_mul.view());} - } - - for col in 0..dim{ - for row in 0..row_indices.len(){ - right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = *add_res.get_mut([row, col]).unwrap(); - } - } - } - - fn row_mul(self, mut right_arr: Array, alpha: Self::Item){ - let right_arr_shape = right_arr.view().shape(); - let dim = self.dim; - let row_indices = self.row_indices; - - for col in 0..dim{ - for row in 0..row_indices.len(){ - right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = alpha*right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] - } - } - - } - - fn row_perm(self, mut right_arr: Array){ - let dim = self.dim; - let col_indices :Vec; - let row_indices :Vec; - - if self.trans{ - col_indices = self.row_indices; - row_indices = self.col_indices; - } - else{ - col_indices = self.col_indices; - row_indices = self.row_indices; - } - - let right_arr_shape = right_arr.view().shape(); - let mut subarr_cols = rlst_dynamic_array2!($scalar, [col_indices.len(), dim]); - - for col in 0..dim{ - for row in 0..col_indices.len(){ - *subarr_cols.get_mut([row, col]).unwrap() = right_arr.data_mut()[col*right_arr_shape[0] + col_indices[row]]; - } - } - - for col in 0..dim{ - for row in 0..row_indices.len(){ - right_arr.data_mut()[col*right_arr_shape[0] + row_indices[row]] = *subarr_cols.get_mut([row, col]).unwrap(); - } - } - } - } - }; -} - -impl_el_ops!(f64); -impl_el_ops!(f32); -impl_el_ops!(c64); -impl_el_ops!(c32); diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 5850b5f9..7740eef8 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -1,10 +1,9 @@ use crate::dense::array::Array; -use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, UnsafeRandomAccessByRef, MultIntoResize, DefaultIterator}; -use crate::dense::types::{RlstResult, RlstScalar}; -use crate::dense::traits::accessors::RandomAccessMut; +use crate::dense::traits::{RandomAccessMut, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, UnsafeRandomAccessByRef, MultIntoResize, DefaultIterator}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, empty_array, rlst_dynamic_array2}; -use crate::dense::types::{c32, c64}; -use num::One; +use crate::DynamicArray; + /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. /// @@ -31,8 +30,8 @@ pub trait MatrixId: RlstScalar { + Shape<2> + RawAccessMut >( - arr: Array, tol: ::Real, k: Option - ) -> RlstResult>; + arr: Array, rank_param: Accuracy<::Real> + ) -> RlstResult>; } macro_rules! implement_into_id { @@ -44,9 +43,9 @@ macro_rules! implement_into_id { + Shape<2> + RawAccessMut >( - arr: Array, tol: ::Real, k: Option - ) -> RlstResult> { - IdDecomposition::<$scalar, ArrayImpl>::new(arr, tol, k) + arr: Array, rank_param: Accuracy<::Real> + ) -> RlstResult> { + IdDecomposition::<$scalar>::new(arr, rank_param) } } }; @@ -66,8 +65,8 @@ impl< > Array { /// Compute the interpolative decomposition of a given 2-dimensional array. - pub fn into_id_alloc(self, tol: ::Real, k: Option) -> RlstResult> { - ::into_id_alloc(self, tol, k) + pub fn into_id_alloc(self, rank_param: Accuracy<::Real>) -> RlstResult> { + ::into_id_alloc(self, rank_param) } } @@ -75,15 +74,21 @@ impl< pub trait MatrixIdDecomposition: Sized { /// Item type type Item: RlstScalar; - /// Array implementaion - type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + RawAccessMut - + Shape<2>; - /// Create a new Interpolative Decomposition from a given array. - fn new(arr: Array, tol: ::Real, k: Option) -> RlstResult; + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item=Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut>(arr: Array, rank_param: Accuracy<::Real>) -> RlstResult; + /// Compute the rank of the decomposition from a given tolerance + fn rank_from_tolerance< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + UnsafeRandomAccessByRef<2, Item = Self::Item>, + >(r_diag: Array, tol: ::Real)->usize; + ///Compute the permutation matrix associated to the Interpolative Decomposition fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> @@ -91,120 +96,136 @@ pub trait MatrixIdDecomposition: Sized { + UnsafeRandomAccessMut<2, Item = Self::Item> + UnsafeRandomAccessByRef<2, Item = Self::Item>, >( - arr: Array, - perm: Vec + &self, arr: Array ); + } ///Stores the relevant features regarding interpolative decomposition. pub struct IdDecomposition< - Item: RlstScalar, - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, + Item: RlstScalar > { - /// arr: permuted array - pub arr: Array, - /// perm_mat: permutation matrix associated to the interpolative decomposition - pub perm_mat: Array, 2>, 2>, + /// skel: skeleton of the interpolative decomposition + pub skel: DynamicArray, + /// perm: permutation associated to the pivoting indiced interpolative decomposition + pub perm: Vec, /// rank: rank of the matrix associated to the interpolative decomposition for a given tolerance pub rank: usize, ///id_mat: interpolative matrix calculated for a given tolerance pub id_mat: Array, 2>, 2>, } +///Options to decide the matrix rank +pub enum Accuracy{ + /// Indicates that the rank of the decomposition will be computed from a given tolerance + Tol(T), + /// Indicates that the rank of the decomposition is given beforehand by the user + FixedRank(usize), + /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range + MaxRank(T, usize) +} + macro_rules! impl_id { ($scalar:ty) => { - impl< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - > MatrixIdDecomposition for IdDecomposition<$scalar, ArrayImpl> + impl MatrixIdDecomposition for IdDecomposition<$scalar> { type Item = $scalar; - type ArrayImpl = ArrayImpl; - fn new(mut arr: Array<$scalar, ArrayImpl, 2>, tol: <$scalar as RlstScalar>::Real, k: Option) -> RlstResult{ - //We assume that for a matrix of m rows and n columns, n>m, so we apply ID the transpose - let mut arr_trans: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); - arr_trans.fill_from(arr.view().conj().transpose()); + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + >(arr: Array<$scalar, ArrayImpl, 2>, rank_param: Accuracy<<$scalar as RlstScalar>::Real>) -> RlstResult{ //We compute the QR decomposition using rlst QR decomposition - let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); + let mut arr_qr: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); arr_qr.fill_from(arr.view().conj().transpose()); let arr_qr_shape = arr_qr.shape(); let qr = arr_qr.view_mut().into_qr_alloc().unwrap(); - + //We obtain relevant parameters of the decomposition: the permutation induced by the pivoting and the R matrix let perm = qr.get_perm(); let mut r_mat = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); qr.get_r(r_mat.view_mut()); - //The maximum rank is given by the number of columns of the transposed matrix let dim: usize = arr_qr_shape[1]; let rank: usize; //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. - match k { - Some(k) => rank = k, - None => { - //We extract the diagonal to calculate the rank of the matrix - let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.view_mut()); - let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); - - //We compute the rank of the matrix - if max.re() > 0.0{ - let alpha: $scalar = (1.0/max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag.iter().filter(|el| el.abs() > tol.into() ).collect::>(); - rank = aux_vec.len(); - } - else{ - rank = dim; - } + + match rank_param{ + Accuracy::Tol(tol) =>{ + rank = Self::rank_from_tolerance(r_mat.view_mut(), tol); + }, + Accuracy::FixedRank(k) => rank = k, + Accuracy::MaxRank(tol, k) =>{ + rank = std::cmp::max(k, Self::rank_from_tolerance(r_mat.view_mut(), tol)); }, } - - let mut perm_mat = rlst_dynamic_array2!($scalar, [dim, dim]); - Self::get_p(perm_mat.view_mut(), perm); - - - let mut perm_arr = empty_array::<$scalar, 2>() - .simple_mult_into_resize(perm_mat.view_mut(), arr.view()); + //We permute arr to extract the columns belonging to the skeleton + let mut permutation = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); + permutation.set_zero(); - for col in 0..arr.shape()[1]{ - for row in 0..arr.shape()[0]{ - arr.data_mut()[col*arr_qr_shape[1] + row] = *perm_arr.get_mut([row, col]).unwrap() - } + for (index, &elem) in perm.iter().enumerate() { + *permutation.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); } + let perm_arr = empty_array::<$scalar, 2>() + .simple_mult_into_resize(permutation.view_mut(), arr.view()); + + let mut skel = empty_array(); + skel.fill_from_resize(perm_arr.into_subview([0,0],[rank, arr_qr_shape[0]])); //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. //If not, compute the Interpolative Decomposition matrix. if rank == 0 || rank >= dim{ let mut id_mat = rlst_dynamic_array2!($scalar, [dim, dim]); id_mat.set_identity(); - Ok(Self{arr, perm_mat, rank, id_mat}) + Ok(Self{skel, perm, rank, id_mat}) } else{ - let shape: [usize; 2] = r_mat.shape(); - let mut id_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>= rlst_dynamic_array2!($scalar, [dim-rank, rank]); + let shape: [usize; 2] = [arr_qr_shape[1], arr_qr_shape[1]]; + let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim-rank, rank]); - let mut k11: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [rank, rank]); + let mut k11: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [rank, rank]); k11.fill_from(r_mat.view_mut().into_subview([0, 0], [rank, rank])); k11.view_mut().into_inverse_alloc().unwrap(); - let mut k12: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [rank, dim-rank]); + let mut k12: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [rank, dim-rank]); k12.fill_from(r_mat.view_mut().into_subview([0, rank], [rank, shape[1]-rank])); let res = empty_array().simple_mult_into_resize(k11.view(), k12.view()); id_mat.fill_from(res.view().conj().transpose().view()); - Ok(Self{arr, perm_mat, rank, id_mat}) + Ok(Self{skel, perm, rank, id_mat}) + } + } + + fn rank_from_tolerance< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = $scalar> + + UnsafeRandomAccessByRef<2, Item = $scalar>, + >(r_mat: Array<$scalar, ArrayImplMut, 2>, tol: <$scalar as RlstScalar>::Real)->usize{ + let dim = r_mat.shape()[0]; + let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.view_mut()); + let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); + + //We compute the rank of the matrix + if max.re() > 0.0{ + let alpha: $scalar = (1.0/max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag.iter().filter(|el| el.abs() > tol.into() ).collect::>(); + aux_vec.len() + } + else{ + dim } } @@ -214,14 +235,11 @@ macro_rules! impl_id { + UnsafeRandomAccessMut<2, Item = $scalar> + UnsafeRandomAccessByRef<2, Item = $scalar>, >( - mut arr: Array<$scalar, ArrayImplMut, 2>, - perm: Vec + &self, mut arr: Array<$scalar, ArrayImplMut, 2> ) { - let m = arr.shape()[0]; - arr.set_zero(); - for col in 0..m { - arr[[col, perm[col]]] = <$scalar as One>::one(); + for (index, &elem) in self.perm.iter().enumerate() { + *arr.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); } } } diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index 8a0b0d22..7d5d626c 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,11 +1,7 @@ use crate::dense::array::Array; -use crate::dense::traits::{ - RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, -}; +use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, rlst_dynamic_array2}; -use crate::dense::traits::DefaultIterator; -use crate::empty_array; +use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, rlst_dynamic_array2, empty_array}; use itertools::min; /// Compute the matrix nullspace. @@ -26,12 +22,6 @@ use itertools::min; /// ``` /// This method allocates memory for the nullspace computation. - -// pub trait MatrixNull { -// /// Compute the matrix null space -// fn into_null_alloc(arr, null_space_type) -> RlstResult>; -// } - pub trait MatrixNull: RlstScalar { /// Compute the matrix null space fn into_null_alloc< @@ -41,7 +31,7 @@ pub trait MatrixNull: RlstScalar { + RawAccessMut, >( arr: Array, - null_space_type: NullSpaceType + tol: ::Real ) -> RlstResult>; } @@ -55,9 +45,9 @@ macro_rules! implement_into_null { + RawAccessMut, >( arr: Array, - null_space_type: NullSpaceType + tol: ::Real ) -> RlstResult> { - NullSpace::<$scalar>::new(arr, null_space_type) + NullSpace::<$scalar>::new(arr, tol) } } }; @@ -77,27 +67,17 @@ impl< > Array { /// Compute the Column or Row nullspace of a given 2-dimensional array. - pub fn into_null_alloc(self, null_space_type: NullSpaceType) -> RlstResult> { - ::into_null_alloc(self, null_space_type) + pub fn into_null_alloc(self, tol: ::Real) -> RlstResult> { + ::into_null_alloc(self, tol) } } -///Null space -#[derive(Clone, Copy)] -#[repr(u8)] -pub enum NullSpaceType { - /// Row Nullspace - Row = b'R', - /// Column Nullspace - Column = b'C', -} +type RealScalar = ::Real; /// QR decomposition pub struct NullSpace< Item: RlstScalar > { - ///Row or column nullspace - pub null_space_type: NullSpaceType, ///Computed null space pub null_space_arr: Array, 2>, 2>, } @@ -106,70 +86,47 @@ macro_rules! implement_null_space { ($scalar:ty) => { impl NullSpace<$scalar> { + fn new< ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + Stride<2> + Shape<2> + RawAccessMut, - >(arr: Array<$scalar, ArrayImpl, 2>, null_space_type: NullSpaceType) -> RlstResult { + >(arr: Array<$scalar, ArrayImpl, 2>, tol: RealScalar<$scalar>) -> RlstResult { let shape = arr.shape(); - - match null_space_type { - NullSpaceType::Row => { - let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); - arr_qr.fill_from(arr.view().conj().transpose()); - let mut null_space_arr = empty_array(); - Self::find_null_space(arr_qr, &mut null_space_arr); - Ok(Self {null_space_type, null_space_arr}) - - }, - NullSpaceType::Column => { - let mut null_space_arr = empty_array(); - Self::find_null_space(arr, &mut null_space_arr); - Ok(Self {null_space_type, null_space_arr}) - }, - } - - } - - fn find_null_space< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >(arr: Array<$scalar, ArrayImpl, 2>, null_space_arr: &mut Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>) - { - let shape = arr.shape(); let dim: usize = min(shape).unwrap(); - let qr = arr.into_qr_alloc().unwrap(); - - //We compute the QR decomposition to find a linearly independent basis of the space. - let mut q: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); - let _ = qr.get_q_alloc(q.view_mut()); + let mut singular_values = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); + let mode = crate::dense::linalg::svd::SvdMode::Full; + let mut u = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); + let mut vt = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); - let mut r_mat = rlst_dynamic_array2!($scalar, [dim, dim]); - qr.get_r(r_mat.view_mut()); + arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); //For a full rank rectangular matrix, then rank = dim. //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(r_mat, dim); + let rank: usize = Self::find_matrix_rank(singular_values, dim, tol); - null_space_arr.fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0]-rank])); + //The null space is given by the last shape[1]-rank columns of V + let mut null_space_arr = empty_array(); + null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); + + Ok(Self {null_space_arr}) + } - fn find_matrix_rank(r_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>, dim: usize)->usize{ + fn find_matrix_rank(singular_values: Array, BaseArray, VectorContainer>, 1>, 1>, dim: usize, tol:<$scalar as RlstScalar>::Real)->usize{ //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.view_mut()); + let mut singular_values_copy: Array, BaseArray, VectorContainer>, 1>, 1> = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); + singular_values_copy.fill_from(singular_values.view()); - let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); + let max: RealScalar<$scalar> = singular_values.view().iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); let rank: usize; if max.re() > 0.0{ - let alpha: $scalar = (1.0/max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag.iter().filter(|el| el.abs() > 1e-15 ).collect::>(); + let alpha: RealScalar<$scalar> = (1.0/max); + singular_values_copy.scale_inplace(alpha); + let aux_vec = singular_values_copy.iter().filter(|el| el.abs() > tol ).collect::>(); rank = aux_vec.len(); } else{ diff --git a/src/prelude.rs b/src/prelude.rs index b759a4cc..f1da0c8f 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -66,7 +66,6 @@ pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; pub use crate::dense::linalg::interpolative_decomposition::{MatrixId, IdDecomposition}; pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; -pub use crate::dense::linalg::elementary_matrix::{ElementaryMatrixData, ElementaryMatrix}; pub use crate::dense::array::rank1_array::Rank1Array; From 9d8a9659604d1f8a5a2fc9af52116b98b52cb188 Mon Sep 17 00:00:00 2001 From: ignacia-fp Date: Fri, 20 Dec 2024 18:47:54 +0000 Subject: [PATCH 05/29] null space implemented for RlstScalar --- src/dense/linalg/null_space.rs | 158 ++++++++++++++++----------------- 1 file changed, 76 insertions(+), 82 deletions(-) diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index 7d5d626c..94fb5465 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,8 +1,10 @@ use crate::dense::array::Array; use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; -use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, rlst_dynamic_array2, empty_array}; +use crate::dense::types::{RlstResult, RlstScalar}; +use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, SvdMode, VectorContainer}; use itertools::min; +use num::{One, Zero}; + /// Compute the matrix nullspace. /// @@ -35,28 +37,20 @@ pub trait MatrixNull: RlstScalar { ) -> RlstResult>; } -macro_rules! implement_into_null { - ($scalar:ty) => { - impl MatrixNull for $scalar { - fn into_null_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, - tol: ::Real - ) -> RlstResult> { - NullSpace::<$scalar>::new(arr, tol) - } - } - }; -} -implement_into_null!(f32); -implement_into_null!(f64); -implement_into_null!(c32); -implement_into_null!(c64); +impl MatrixNull for T { + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + tol: ::Real + ) -> RlstResult> { + NullSpace::::new(arr, tol) + } +} impl< Item: RlstScalar + MatrixNull, @@ -74,7 +68,7 @@ impl< type RealScalar = ::Real; -/// QR decomposition +/// Null Space pub struct NullSpace< Item: RlstScalar > { @@ -82,65 +76,65 @@ pub struct NullSpace< pub null_space_arr: Array, 2>, 2>, } -macro_rules! implement_null_space { - ($scalar:ty) => { - impl NullSpace<$scalar> - { - - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >(arr: Array<$scalar, ArrayImpl, 2>, tol: RealScalar<$scalar>) -> RlstResult { - - let shape = arr.shape(); - let dim: usize = min(shape).unwrap(); - let mut singular_values = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); - let mode = crate::dense::linalg::svd::SvdMode::Full; - let mut u = rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); - let mut vt = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); - - arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); - - //For a full rank rectangular matrix, then rank = dim. - //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(singular_values, dim, tol); - - //The null space is given by the last shape[1]-rank columns of V - let mut null_space_arr = empty_array(); - null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); - - Ok(Self {null_space_arr}) - - } - - fn find_matrix_rank(singular_values: Array, BaseArray, VectorContainer>, 1>, 1>, dim: usize, tol:<$scalar as RlstScalar>::Real)->usize{ - //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let mut singular_values_copy: Array, BaseArray, VectorContainer>, 1>, 1> = rlst_dynamic_array1!(RealScalar<$scalar>, [dim]); - singular_values_copy.fill_from(singular_values.view()); - - let max: RealScalar<$scalar> = singular_values.view().iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); - let rank: usize; - - if max.re() > 0.0{ - let alpha: RealScalar<$scalar> = (1.0/max); - singular_values_copy.scale_inplace(alpha); - let aux_vec = singular_values_copy.iter().filter(|el| el.abs() > tol ).collect::>(); - rank = aux_vec.len(); - } - else{ - rank = dim; - } - rank - - } +///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct +pub trait NullSpaceComputation{ + ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) + type Item: RlstScalar; + + ///We create the null space decomposition + fn new + + Stride<2> + RawAccessMut + Shape<2>> + (arr: Array, tol: RealScalar) -> RlstResult where Self: Sized; + + ///This function helps us to find the matrix rank + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:RealScalar)->usize; +} + + +impl NullSpaceComputation for NullSpace { + type Item = T; + + fn new + + Stride<2> + + RawAccessMut + + Shape<2>>(arr: Array, tol: RealScalar) -> RlstResult { + + let shape: [usize; 2] = arr.shape(); + let dim: usize = min(shape).unwrap(); + let mut singular_values: DynamicArray, 1> = rlst_dynamic_array1!(RealScalar, [dim]); + let mode: SvdMode = SvdMode::Full; + let mut u: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); + let mut vt: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); + + arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); + + //The null space is given by the last shape[1]-rank columns of V + let mut null_space_arr: DynamicArray = empty_array(); + null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); + + Ok(Self{null_space_arr}) + + } + + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:::Real)->usize{ + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let max: RealScalar = singular_values.view().iter().max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()).unwrap().abs(); + let mut rank: usize = dim; + + if max.re() > as Zero>::zero(){ + let alpha: RealScalar = as One>::one()/max; + singular_values.scale_inplace(alpha); + let aux_vec: Vec> = singular_values.iter().filter(|el| el.abs() > tol ).collect::>>(); + rank = aux_vec.len(); } - }; + + rank + + } } -implement_null_space!(f64); -implement_null_space!(f32); -implement_null_space!(c64); -implement_null_space!(c32); From 0c7d81286fad528c44afa596467f43c726b4b856 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 12:35:28 +0100 Subject: [PATCH 06/29] fixed implementation of id and null space --- examples/rlst_interpolative_decomposition.rs | 70 ++--- examples/rlst_null_space.rs | 13 +- src/dense/linalg.rs | 26 +- .../linalg/interpolative_decomposition.rs | 283 ++++++++++-------- src/dense/linalg/null_space.rs | 243 +++++++++------ src/doc/dense_linear_algebra.rs | 18 +- src/prelude.rs | 4 +- 7 files changed, 382 insertions(+), 275 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 44a106fb..47a19f9b 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -1,10 +1,9 @@ -//! Demo the inverse of a matrix +//! Interpolative decomposition. pub use rlst::prelude::*; -use rlst::dense::linalg::interpolative_decomposition::{MatrixIdDecomposition, Accuracy}; //Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere. -fn low_rank_matrix(n: usize, arr: &mut Array, 2>, 2>){ +fn low_rank_matrix(n: usize, arr: &mut Array, 2>, 2>) { //Obtain n equally distributed angles 0 = rlst_dynamic_array2!(f64, [slice, n]); + let mut arr = rlst_dynamic_array2!(f64, [slice, n]); low_rank_matrix(n, &mut arr); - let res: IdDecomposition = arr.view_mut().into_id_alloc(Accuracy::Tol(tol)).unwrap(); + let mut res = arr.r_mut().into_id_alloc(tol, None).unwrap(); - println!("The skeleton of the matrix is given by"); - res.skel.pretty_print(); + println!("The permuted matrix is:"); + res.arr.pretty_print(); println!("The interpolative decomposition matrix is:"); res.id_mat.pretty_print(); println!("The rank of this matrix is {}\n", res.rank); - //We extract the residuals of the matrix - let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); - res.get_p(perm_mat.view_mut()); - let perm_arr: DynamicArray = empty_array::() - .simple_mult_into_resize(perm_mat.view_mut(), arr.view()); - - let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice-res.rank, n]); - a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice-res.rank, n])); + let a_rs_app = empty_array().simple_mult_into_resize( + res.id_mat.r(), + res.arr.r_mut().into_subview([0, 0], [res.rank, n]), + ); + let a_rs = res + .arr + .r_mut() + .into_subview([res.rank, 0], [slice - res.rank, n]); - //We compute an approximation of the residual columns of the matrix - let a_rs_app: DynamicArray = empty_array().simple_mult_into_resize(res.id_mat.view(), res.skel); - - let error: f64 = a_rs.view().sub(a_rs_app.view()).view_flat().norm_2(); + let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); } diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index 5c0afd7c..955ce568 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -1,14 +1,17 @@ -//! Demo the inverse of a matrix +//! Demo the null space of a matrix. +use rlst::dense::linalg::null_space::NullSpaceType; pub use rlst::prelude::*; //Here we compute the null space (B) of the rowspace of a short matrix (A). pub fn main() { let mut arr = rlst_dynamic_array2!(f64, [3, 4]); arr.fill_from_seed_equally_distributed(0); - let tol = 1e-15; - let null_res = arr.view_mut().into_null_alloc(tol).unwrap(); - let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.view_mut(), null_res.null_space_arr.view()); - println!("Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2()); + let null_res = arr.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); + let res = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); + println!( + "Value of |A*B|_2, where B=null(A): {}", + res.view_flat().norm_2() + ); } diff --git a/src/dense/linalg.rs b/src/dense/linalg.rs index b8327b1f..7e2b41f5 100644 --- a/src/dense/linalg.rs +++ b/src/dense/linalg.rs @@ -1,17 +1,18 @@ //! Linear algebra routines - -use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixId, MatrixSvd, MatrixNull, RlstScalar}; +use crate::{ + MatrixId, MatrixInverse, MatrixNull, MatrixPseudoInverse, MatrixQr, MatrixSvd, RlstScalar, +}; use self::lu::MatrixLu; +pub mod interpolative_decomposition; pub mod inverse; pub mod lu; +pub mod null_space; pub mod pseudo_inverse; pub mod qr; pub mod svd; -pub mod interpolative_decomposition; -pub mod null_space; /// Return true if stride is column major as required by Lapack. pub fn assert_lapack_stride(stride: [usize; 2]) { @@ -23,10 +24,21 @@ pub fn assert_lapack_stride(stride: [usize; 2]) { } /// Marker trait for objects that support Matrix decompositions. -pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull{} +pub trait LinAlg: + MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull +{ +} -impl LinAlg - for T +impl< + T: RlstScalar + + MatrixInverse + + MatrixQr + + MatrixSvd + + MatrixLu + + MatrixPseudoInverse + + MatrixId + + MatrixNull, + > LinAlg for T { } diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 7740eef8..ac54c664 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -1,9 +1,14 @@ +//! Demo the null space of a matrix. use crate::dense::array::Array; -use crate::dense::traits::{RandomAccessMut, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, UnsafeRandomAccessByRef, MultIntoResize, DefaultIterator}; -use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::{BaseArray, VectorContainer, rlst_dynamic_array1, empty_array, rlst_dynamic_array2}; -use crate::DynamicArray; - +use crate::dense::traits::accessors::RandomAccessMut; +use crate::dense::traits::{ + DefaultIterator, MultIntoResize, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, + UnsafeRandomAccessByValue, UnsafeRandomAccessMut, +}; +use crate::dense::types::{c32, c64}; +use crate::dense::types::{RlstResult, RlstScalar}; +use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, VectorContainer}; +use num::One; /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. /// @@ -12,15 +17,14 @@ use crate::DynamicArray; /// /// # Example /// -/// The following command computes the interpolative decomposition of an array `a` for a given tolerance, tol. +/// The following command computes the interpolative decomposition of an array `a` for a given tolerance, tol. /// ``` /// # use rlst::rlst_dynamic_array2; /// # let tol: f64 = 1e-5; /// # let mut a = rlst_dynamic_array2!(f64, [50, 100]); /// # a.fill_from_seed_equally_distributed(0); -/// # let res = a.view_mut().into_id_alloc(tol, None).unwrap(); +/// # let res = a.r_mut().into_id_alloc(tol, None).unwrap(); /// ``` - /// This method allocates memory for the interpolative decomposition. pub trait MatrixId: RlstScalar { ///This method allocates space for ID @@ -28,10 +32,12 @@ pub trait MatrixId: RlstScalar { ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + Stride<2> + Shape<2> - + RawAccessMut + + RawAccessMut, >( - arr: Array, rank_param: Accuracy<::Real> - ) -> RlstResult>; + arr: Array, + tol: ::Real, + k: Option, + ) -> RlstResult>; } macro_rules! implement_into_id { @@ -41,11 +47,13 @@ macro_rules! implement_into_id { ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + Stride<2> + Shape<2> - + RawAccessMut + + RawAccessMut, >( - arr: Array, rank_param: Accuracy<::Real> - ) -> RlstResult> { - IdDecomposition::<$scalar>::new(arr, rank_param) + arr: Array, + tol: ::Real, + k: Option, + ) -> RlstResult> { + IdDecomposition::<$scalar, ArrayImpl>::new(arr, tol, k) } } }; @@ -65,8 +73,12 @@ impl< > Array { /// Compute the interpolative decomposition of a given 2-dimensional array. - pub fn into_id_alloc(self, rank_param: Accuracy<::Real>) -> RlstResult> { - ::into_id_alloc(self, rank_param) + pub fn into_id_alloc( + self, + tol: ::Real, + k: Option, + ) -> RlstResult> { + ::into_id_alloc(self, tol, k) } } @@ -74,21 +86,19 @@ impl< pub trait MatrixIdDecomposition: Sized { /// Item type type Item: RlstScalar; - /// Create a new Interpolative Decomposition from a given array. - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item=Self::Item> + /// Array implementaion + type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + Stride<2> - + Shape<2> - + RawAccessMut>(arr: Array, rank_param: Accuracy<::Real>) -> RlstResult; + + RawAccessMut + + Shape<2>; + + /// Create a new Interpolative Decomposition from a given array. + fn new( + arr: Array, + tol: ::Real, + k: Option, + ) -> RlstResult; - /// Compute the rank of the decomposition from a given tolerance - fn rank_from_tolerance< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = Self::Item> - + UnsafeRandomAccessByRef<2, Item = Self::Item>, - >(r_diag: Array, tol: ::Real)->usize; - ///Compute the permutation matrix associated to the Interpolative Decomposition fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> @@ -96,136 +106,158 @@ pub trait MatrixIdDecomposition: Sized { + UnsafeRandomAccessMut<2, Item = Self::Item> + UnsafeRandomAccessByRef<2, Item = Self::Item>, >( - &self, arr: Array + arr: Array, + perm: Vec, ); - - - } -///Stores the relevant features regarding interpolative decomposition. +///Stores the relevant features regarding interpolative decomposition. pub struct IdDecomposition< - Item: RlstScalar + Item: RlstScalar, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, > { - /// skel: skeleton of the interpolative decomposition - pub skel: DynamicArray, - /// perm: permutation associated to the pivoting indiced interpolative decomposition - pub perm: Vec, + /// arr: permuted array + pub arr: Array, + /// perm_mat: permutation matrix associated to the interpolative decomposition + pub perm_mat: Array, 2>, 2>, /// rank: rank of the matrix associated to the interpolative decomposition for a given tolerance pub rank: usize, ///id_mat: interpolative matrix calculated for a given tolerance pub id_mat: Array, 2>, 2>, } -///Options to decide the matrix rank -pub enum Accuracy{ - /// Indicates that the rank of the decomposition will be computed from a given tolerance - Tol(T), - /// Indicates that the rank of the decomposition is given beforehand by the user - FixedRank(usize), - /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range - MaxRank(T, usize) -} - - macro_rules! impl_id { ($scalar:ty) => { - impl MatrixIdDecomposition for IdDecomposition<$scalar> + impl< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + > MatrixIdDecomposition for IdDecomposition<$scalar, ArrayImpl> { type Item = $scalar; + type ArrayImpl = ArrayImpl; - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >(arr: Array<$scalar, ArrayImpl, 2>, rank_param: Accuracy<<$scalar as RlstScalar>::Real>) -> RlstResult{ + fn new( + mut arr: Array<$scalar, ArrayImpl, 2>, + tol: <$scalar as RlstScalar>::Real, + k: Option, + ) -> RlstResult { + //We assume that for a matrix of m rows and n columns, n>m, so we apply ID the transpose + let mut arr_trans: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + > = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); + arr_trans.fill_from(arr.r().conj().transpose()); //We compute the QR decomposition using rlst QR decomposition - let mut arr_qr: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); - arr_qr.fill_from(arr.view().conj().transpose()); + let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = + rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); + arr_qr.fill_from(arr.r().conj().transpose()); let arr_qr_shape = arr_qr.shape(); - let qr = arr_qr.view_mut().into_qr_alloc().unwrap(); - + let qr = arr_qr.r_mut().into_qr_alloc().unwrap(); + //We obtain relevant parameters of the decomposition: the permutation induced by the pivoting and the R matrix let perm = qr.get_perm(); let mut r_mat = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); - qr.get_r(r_mat.view_mut()); + qr.get_r(r_mat.r_mut()); //The maximum rank is given by the number of columns of the transposed matrix let dim: usize = arr_qr_shape[1]; let rank: usize; //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. + match k { + Some(k) => rank = k, + None => { + //We extract the diagonal to calculate the rank of the matrix + let mut r_diag: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 1>, + 1, + > = rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.r_mut()); + let max: $scalar = r_diag + .iter() + .max_by(|a, b| a.abs().total_cmp(&b.abs())) + .unwrap() + .abs() + .into(); - match rank_param{ - Accuracy::Tol(tol) =>{ - rank = Self::rank_from_tolerance(r_mat.view_mut(), tol); - }, - Accuracy::FixedRank(k) => rank = k, - Accuracy::MaxRank(tol, k) =>{ - rank = std::cmp::max(k, Self::rank_from_tolerance(r_mat.view_mut(), tol)); - }, + //We compute the rank of the matrix + if max.re() > 0.0 { + let alpha: $scalar = (1.0 / max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag + .iter() + .filter(|el| el.abs() > tol.into()) + .collect::>(); + rank = aux_vec.len(); + } else { + rank = dim; + } + } } - //We permute arr to extract the columns belonging to the skeleton - let mut permutation = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); - permutation.set_zero(); - - for (index, &elem) in perm.iter().enumerate() { - *permutation.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); - } + let mut perm_mat = rlst_dynamic_array2!($scalar, [dim, dim]); + Self::get_p(perm_mat.r_mut(), perm); - let perm_arr = empty_array::<$scalar, 2>() - .simple_mult_into_resize(permutation.view_mut(), arr.view()); + let mut perm_arr = + empty_array::<$scalar, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); - let mut skel = empty_array(); - skel.fill_from_resize(perm_arr.into_subview([0,0],[rank, arr_qr_shape[0]])); + for col in 0..arr.shape()[1] { + for row in 0..arr.shape()[0] { + arr.data_mut()[col * arr_qr_shape[1] + row] = + *perm_arr.get_mut([row, col]).unwrap() + } + } //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. //If not, compute the Interpolative Decomposition matrix. - if rank == 0 || rank >= dim{ + if rank == 0 || rank >= dim { let mut id_mat = rlst_dynamic_array2!($scalar, [dim, dim]); id_mat.set_identity(); - Ok(Self{skel, perm, rank, id_mat}) - } - else{ - - let shape: [usize; 2] = [arr_qr_shape[1], arr_qr_shape[1]]; - let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim-rank, rank]); - - let mut k11: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [rank, rank]); - k11.fill_from(r_mat.view_mut().into_subview([0, 0], [rank, rank])); - k11.view_mut().into_inverse_alloc().unwrap(); - let mut k12: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [rank, dim-rank]); - k12.fill_from(r_mat.view_mut().into_subview([0, rank], [rank, shape[1]-rank])); - - let res = empty_array().simple_mult_into_resize(k11.view(), k12.view()); - id_mat.fill_from(res.view().conj().transpose().view()); - Ok(Self{skel, perm, rank, id_mat}) - } - } + Ok(Self { + arr, + perm_mat, + rank, + id_mat, + }) + } else { + let shape: [usize; 2] = r_mat.shape(); + let mut id_mat: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + > = rlst_dynamic_array2!($scalar, [dim - rank, rank]); - fn rank_from_tolerance< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = $scalar> - + UnsafeRandomAccessByRef<2, Item = $scalar>, - >(r_mat: Array<$scalar, ArrayImplMut, 2>, tol: <$scalar as RlstScalar>::Real)->usize{ - let dim = r_mat.shape()[0]; - let mut r_diag:Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.view_mut()); - let max: $scalar = r_diag.iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap().abs().into(); - - //We compute the rank of the matrix - if max.re() > 0.0{ - let alpha: $scalar = (1.0/max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag.iter().filter(|el| el.abs() > tol.into() ).collect::>(); - aux_vec.len() - } - else{ - dim + let mut k11: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + > = rlst_dynamic_array2!($scalar, [rank, rank]); + k11.fill_from(r_mat.r_mut().into_subview([0, 0], [rank, rank])); + k11.r_mut().into_inverse_alloc().unwrap(); + let mut k12: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + > = rlst_dynamic_array2!($scalar, [rank, dim - rank]); + k12.fill_from( + r_mat + .r_mut() + .into_subview([0, rank], [rank, shape[1] - rank]), + ); + + let res = empty_array().simple_mult_into_resize(k11.r(), k12.r()); + id_mat.fill_from(res.r().conj().transpose().r()); + Ok(Self { + arr, + perm_mat, + rank, + id_mat, + }) } } @@ -235,11 +267,14 @@ macro_rules! impl_id { + UnsafeRandomAccessMut<2, Item = $scalar> + UnsafeRandomAccessByRef<2, Item = $scalar>, >( - &self, mut arr: Array<$scalar, ArrayImplMut, 2> + mut arr: Array<$scalar, ArrayImplMut, 2>, + perm: Vec, ) { + let m = arr.shape()[0]; + arr.set_zero(); - for (index, &elem) in self.perm.iter().enumerate() { - *arr.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); + for col in 0..m { + arr[[col, perm[col]]] = <$scalar as One>::one(); } } } diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index 94fb5465..b9131e1c 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,10 +1,11 @@ +//! Null space. use crate::dense::array::Array; -use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; -use crate::dense::types::{RlstResult, RlstScalar}; -use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, SvdMode, VectorContainer}; +use crate::dense::traits::DefaultIterator; +use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; +use crate::empty_array; +use crate::{rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, VectorContainer}; use itertools::min; -use num::{One, Zero}; - /// Compute the matrix nullspace. /// @@ -13,16 +14,15 @@ use num::{One, Zero}; /// /// # Example /// -/// The following command computes the nullspace of an array `a`. -/// The nullspace is found in +/// The following command computes the nullspace of an array `a`. +/// The nullspace is found in /// ``` /// # use rlst::rlst_dynamic_array2; /// # use rlst::dense::linalg::null_space::{NullSpaceType, MatrixNull}; /// # let mut a = rlst_dynamic_array2!(f64, [3, 4]); /// # a.fill_from_seed_equally_distributed(0); -/// # let null_res = a.view_mut().into_null_alloc(NullSpaceType::Row).unwrap(); +/// # let null_res = a.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); /// ``` - /// This method allocates memory for the nullspace computation. pub trait MatrixNull: RlstScalar { /// Compute the matrix null space @@ -33,25 +33,33 @@ pub trait MatrixNull: RlstScalar { + RawAccessMut, >( arr: Array, - tol: ::Real + null_space_type: NullSpaceType, ) -> RlstResult>; } - -impl MatrixNull for T { - fn into_null_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, - tol: ::Real - ) -> RlstResult> { - NullSpace::::new(arr, tol) - } +macro_rules! implement_into_null { + ($scalar:ty) => { + impl MatrixNull for $scalar { + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + null_space_type: NullSpaceType, + ) -> RlstResult> { + NullSpace::<$scalar>::new(arr, null_space_type) + } + } + }; } +implement_into_null!(f32); +implement_into_null!(f64); +implement_into_null!(c32); +implement_into_null!(c64); + impl< Item: RlstScalar + MatrixNull, ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> @@ -61,80 +69,137 @@ impl< > Array { /// Compute the Column or Row nullspace of a given 2-dimensional array. - pub fn into_null_alloc(self, tol: ::Real) -> RlstResult> { - ::into_null_alloc(self, tol) + pub fn into_null_alloc(self, null_space_type: NullSpaceType) -> RlstResult> { + ::into_null_alloc(self, null_space_type) } } +///Null space +#[derive(Clone, Copy)] +#[repr(u8)] +pub enum NullSpaceType { + /// Row Nullspace + Row = b'R', + /// Column Nullspace + Column = b'C', +} -type RealScalar = ::Real; -/// Null Space -pub struct NullSpace< - Item: RlstScalar -> { +/// QR decomposition +pub struct NullSpace { + ///Row or column nullspace + pub null_space_type: NullSpaceType, ///Computed null space pub null_space_arr: Array, 2>, 2>, } -///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct -pub trait NullSpaceComputation{ - ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) - type Item: RlstScalar; - - ///We create the null space decomposition - fn new - + Stride<2> + RawAccessMut + Shape<2>> - (arr: Array, tol: RealScalar) -> RlstResult where Self: Sized; - - ///This function helps us to find the matrix rank - fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:RealScalar)->usize; -} - - -impl NullSpaceComputation for NullSpace { - type Item = T; - - fn new - + Stride<2> - + RawAccessMut - + Shape<2>>(arr: Array, tol: RealScalar) -> RlstResult { - - let shape: [usize; 2] = arr.shape(); - let dim: usize = min(shape).unwrap(); - let mut singular_values: DynamicArray, 1> = rlst_dynamic_array1!(RealScalar, [dim]); - let mode: SvdMode = SvdMode::Full; - let mut u: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); - let mut vt: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); - - arr.into_svd_alloc(u.view_mut(), vt.view_mut(), singular_values.data_mut(), mode).unwrap(); - - //For a full rank rectangular matrix, then rank = dim. - //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); - - //The null space is given by the last shape[1]-rank columns of V - let mut null_space_arr: DynamicArray = empty_array(); - null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); - - Ok(Self{null_space_arr}) - - } - - fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:::Real)->usize{ - //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let max: RealScalar = singular_values.view().iter().max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()).unwrap().abs(); - let mut rank: usize = dim; - - if max.re() > as Zero>::zero(){ - let alpha: RealScalar = as One>::one()/max; - singular_values.scale_inplace(alpha); - let aux_vec: Vec> = singular_values.iter().filter(|el| el.abs() > tol ).collect::>>(); - rank = aux_vec.len(); +macro_rules! implement_null_space { + ($scalar:ty) => { + impl NullSpace<$scalar> { + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array<$scalar, ArrayImpl, 2>, + null_space_type: NullSpaceType, + ) -> RlstResult { + let shape = arr.shape(); + + match null_space_type { + NullSpaceType::Row => { + let mut arr_qr: Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + > = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); + arr_qr.fill_from(arr.r().conj().transpose()); + let mut null_space_arr = empty_array(); + Self::find_null_space(arr_qr, &mut null_space_arr); + Ok(Self { + null_space_type, + null_space_arr, + }) + } + NullSpaceType::Column => { + let mut null_space_arr = empty_array(); + Self::find_null_space(arr, &mut null_space_arr); + Ok(Self { + null_space_type, + null_space_arr, + }) + } + } + } + + fn find_null_space< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array<$scalar, ArrayImpl, 2>, + null_space_arr: &mut Array< + $scalar, + BaseArray<$scalar, VectorContainer<$scalar>, 2>, + 2, + >, + ) { + let shape = arr.shape(); + let dim: usize = min(shape).unwrap(); + let qr = arr.into_qr_alloc().unwrap(); + + //We compute the QR decomposition to find a linearly independent basis of the space. + let mut q: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = + rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); + let _ = qr.get_q_alloc(q.r_mut()); + + let mut r_mat = rlst_dynamic_array2!($scalar, [dim, dim]); + qr.get_r(r_mat.r_mut()); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = Self::find_matrix_rank(r_mat, dim); + + null_space_arr + .fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0] - rank])); + } + + fn find_matrix_rank( + r_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>, + dim: usize, + ) -> usize { + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let mut r_diag: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = + rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.r_mut()); + + let max: $scalar = r_diag + .iter() + .max_by(|a, b| a.abs().total_cmp(&b.abs())) + .unwrap() + .abs() + .into(); + let rank: usize; + + if max.re() > 0.0 { + let alpha: $scalar = (1.0 / max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag + .iter() + .filter(|el| el.abs() > 1e-15) + .collect::>(); + rank = aux_vec.len(); + } else { + rank = dim; + } + rank + } } - - rank - - } + }; } - +implement_null_space!(f64); +implement_null_space!(f32); +implement_null_space!(c64); +implement_null_space!(c32); diff --git a/src/doc/dense_linear_algebra.rs b/src/doc/dense_linear_algebra.rs index 2965c0b9..39ffc587 100644 --- a/src/doc/dense_linear_algebra.rs +++ b/src/doc/dense_linear_algebra.rs @@ -515,9 +515,9 @@ //! To compute the full SVD use the parameter [SvdMode::Full](crate::SvdMode::Full). //! //! # Interpolative decomposition -//! +//! //! To compute the Interpolative Decomposition of a two-dimensional array (long or square) `arr`, for a given tolerance, use -//! +//! //! ``` //! # use rlst::prelude::*; //! let mut rand = rand::thread_rng(); @@ -525,21 +525,21 @@ //! let tol: f64 = 1e-5; //! let mut res = arr.view_mut().into_id_alloc(tol, None).unwrap(); //! ``` -//! +//! //! You can take a look at the interpolation matrix using -//! +//! //! ``` //! res.id_mat.pretty_print(); //! ``` -//! +//! //! With this, the approximation of the low-rank elements of `arr` is given by -//! +//! //! ``` //! let arr_app = empty_array().simple_mult_into_resize(res.id_mat.view(), res.arr.view_mut().into_subview([0, 0], [res.rank, n])); //! ``` -//! +//! //! You can also obtain the Interpolative Decomposition by giving the algorithm a fixed rank, and the tolerance parameter is ignored, like so -//! +//! //! ``` //! # use rlst::prelude::*; //! let mut rand = rand::thread_rng(); @@ -548,7 +548,7 @@ //! k = 2; //! let mut res = arr.view_mut().into_id_alloc(tol, k).unwrap(); //! ``` -//! +//! //! # Other vector functions //! //! The following other functions for arrays with a single dimension (i.e. vectors) are provided. diff --git a/src/prelude.rs b/src/prelude.rs index 7e4228ae..8c6dbc5a 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -61,13 +61,13 @@ pub use crate::dense::data_container::{ pub use crate::dense::array::empty_axis::AxisPosition; +pub use crate::dense::linalg::interpolative_decomposition::{IdDecomposition, MatrixId}; pub use crate::dense::linalg::inverse::MatrixInverse; pub use crate::dense::linalg::lu::{LuDecomposition, MatrixLuDecomposition}; +pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; pub use crate::dense::linalg::pseudo_inverse::MatrixPseudoInverse; pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; -pub use crate::dense::linalg::interpolative_decomposition::{MatrixId, IdDecomposition}; -pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; pub use crate::dense::array::rank1_array::Rank1Array; From 841eaf55c3076524d58e657f315f55958bddcfe0 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 14:26:50 +0100 Subject: [PATCH 07/29] syn updated --- proc-macro/Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/proc-macro/Cargo.toml b/proc-macro/Cargo.toml index 73a4adab..39ec5d94 100644 --- a/proc-macro/Cargo.toml +++ b/proc-macro/Cargo.toml @@ -19,7 +19,7 @@ categories = ["mathematics", "science"] proc-macro = true [dependencies] -syn = { version = "1.0", features = ["full"]} +syn = { version = "2.0.100", features = ["full"]} quote = "1.0" darling = "0.20" From 183b91dad6e1698713b8787acde8ad9292cb3c02 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 14:47:28 +0100 Subject: [PATCH 08/29] fixes to interpolative decomposition --- examples/rlst_interpolative_decomposition.rs | 31 ++- .../linalg/interpolative_decomposition.rs | 242 +++++++++--------- 2 files changed, 141 insertions(+), 132 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 47a19f9b..e779de7a 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -1,5 +1,6 @@ -//! Interpolative decomposition. +//! Demo interpolative decomposition of a matrix. +use rlst::dense::linalg::interpolative_decomposition::{Accuracy, MatrixIdDecomposition}; pub use rlst::prelude::*; //Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere. @@ -50,27 +51,31 @@ pub fn main() { let tol: f64 = 1e-5; //Create a low rank matrix - let mut arr = rlst_dynamic_array2!(f64, [slice, n]); + let mut arr: DynamicArray = rlst_dynamic_array2!(f64, [slice, n]); low_rank_matrix(n, &mut arr); - let mut res = arr.r_mut().into_id_alloc(tol, None).unwrap(); + let res: IdDecomposition = arr.r_mut().into_id_alloc(Accuracy::Tol(tol)).unwrap(); - println!("The permuted matrix is:"); - res.arr.pretty_print(); + println!("The skeleton of the matrix is given by"); + res.skel.pretty_print(); println!("The interpolative decomposition matrix is:"); res.id_mat.pretty_print(); println!("The rank of this matrix is {}\n", res.rank); - let a_rs_app = empty_array().simple_mult_into_resize( - res.id_mat.r(), - res.arr.r_mut().into_subview([0, 0], [res.rank, n]), - ); - let a_rs = res - .arr - .r_mut() - .into_subview([res.rank, 0], [slice - res.rank, n]); + //We extract the residuals of the matrix + let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); + res.get_p(perm_mat.r_mut()); + let perm_arr: DynamicArray = + empty_array::().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); + + let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice - res.rank, n]); + a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n])); + + //We compute an approximation of the residual columns of the matrix + let a_rs_app: DynamicArray = + empty_array().simple_mult_into_resize(res.id_mat.r(), res.skel); let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index ac54c664..f3b3920b 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -1,14 +1,12 @@ -//! Demo the null space of a matrix. +//! Interpolative decomposition of a matrix. use crate::dense::array::Array; -use crate::dense::traits::accessors::RandomAccessMut; use crate::dense::traits::{ - DefaultIterator, MultIntoResize, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, - UnsafeRandomAccessByValue, UnsafeRandomAccessMut, + DefaultIterator, MultIntoResize, RandomAccessMut, RawAccessMut, Shape, Stride, + UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; -use crate::dense::types::{c32, c64}; -use crate::dense::types::{RlstResult, RlstScalar}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; +use crate::DynamicArray; use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, VectorContainer}; -use num::One; /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. /// @@ -35,9 +33,8 @@ pub trait MatrixId: RlstScalar { + RawAccessMut, >( arr: Array, - tol: ::Real, - k: Option, - ) -> RlstResult>; + rank_param: Accuracy<::Real>, + ) -> RlstResult>; } macro_rules! implement_into_id { @@ -50,10 +47,9 @@ macro_rules! implement_into_id { + RawAccessMut, >( arr: Array, - tol: ::Real, - k: Option, - ) -> RlstResult> { - IdDecomposition::<$scalar, ArrayImpl>::new(arr, tol, k) + rank_param: Accuracy<::Real>, + ) -> RlstResult> { + IdDecomposition::<$scalar>::new(arr, rank_param) } } }; @@ -75,10 +71,9 @@ impl< /// Compute the interpolative decomposition of a given 2-dimensional array. pub fn into_id_alloc( self, - tol: ::Real, - k: Option, - ) -> RlstResult> { - ::into_id_alloc(self, tol, k) + rank_param: Accuracy<::Real>, + ) -> RlstResult> { + ::into_id_alloc(self, rank_param) } } @@ -86,19 +81,28 @@ impl< pub trait MatrixIdDecomposition: Sized { /// Item type type Item: RlstScalar; - /// Array implementaion - type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + RawAccessMut - + Shape<2>; - /// Create a new Interpolative Decomposition from a given array. - fn new( - arr: Array, - tol: ::Real, - k: Option, + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + rank_param: Accuracy<::Real>, ) -> RlstResult; + /// Compute the rank of the decomposition from a given tolerance + fn rank_from_tolerance< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + UnsafeRandomAccessByRef<2, Item = Self::Item>, + >( + r_diag: Array, + tol: ::Real, + ) -> usize; + ///Compute the permutation matrix associated to the Interpolative Decomposition fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> @@ -106,53 +110,49 @@ pub trait MatrixIdDecomposition: Sized { + UnsafeRandomAccessMut<2, Item = Self::Item> + UnsafeRandomAccessByRef<2, Item = Self::Item>, >( + &self, arr: Array, - perm: Vec, ); } ///Stores the relevant features regarding interpolative decomposition. -pub struct IdDecomposition< - Item: RlstScalar, - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, -> { - /// arr: permuted array - pub arr: Array, - /// perm_mat: permutation matrix associated to the interpolative decomposition - pub perm_mat: Array, 2>, 2>, +pub struct IdDecomposition { + /// skel: skeleton of the interpolative decomposition + pub skel: DynamicArray, + /// perm: permutation associated to the pivoting indiced interpolative decomposition + pub perm: Vec, /// rank: rank of the matrix associated to the interpolative decomposition for a given tolerance pub rank: usize, ///id_mat: interpolative matrix calculated for a given tolerance pub id_mat: Array, 2>, 2>, } +///Options to decide the matrix rank +pub enum Accuracy { + /// Indicates that the rank of the decomposition will be computed from a given tolerance + Tol(T), + /// Indicates that the rank of the decomposition is given beforehand by the user + FixedRank(usize), + /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range + MaxRank(T, usize), +} + macro_rules! impl_id { ($scalar:ty) => { - impl< + impl MatrixIdDecomposition for IdDecomposition<$scalar> { + type Item = $scalar; + + fn new< ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + Stride<2> + Shape<2> + RawAccessMut, - > MatrixIdDecomposition for IdDecomposition<$scalar, ArrayImpl> - { - type Item = $scalar; - type ArrayImpl = ArrayImpl; - - fn new( - mut arr: Array<$scalar, ArrayImpl, 2>, - tol: <$scalar as RlstScalar>::Real, - k: Option, + >( + arr: Array<$scalar, ArrayImpl, 2>, + rank_param: Accuracy<<$scalar as RlstScalar>::Real>, ) -> RlstResult { - //We assume that for a matrix of m rows and n columns, n>m, so we apply ID the transpose - let mut arr_trans: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - > = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); - arr_trans.fill_from(arr.r().conj().transpose()); - //We compute the QR decomposition using rlst QR decomposition - let mut arr_qr: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = + let mut arr_qr: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); arr_qr.fill_from(arr.r().conj().transpose()); let arr_qr_shape = arr_qr.shape(); @@ -168,82 +168,54 @@ macro_rules! impl_id { let rank: usize; //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. - match k { - Some(k) => rank = k, - None => { - //We extract the diagonal to calculate the rank of the matrix - let mut r_diag: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 1>, - 1, - > = rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.r_mut()); - let max: $scalar = r_diag - .iter() - .max_by(|a, b| a.abs().total_cmp(&b.abs())) - .unwrap() - .abs() - .into(); - //We compute the rank of the matrix - if max.re() > 0.0 { - let alpha: $scalar = (1.0 / max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag - .iter() - .filter(|el| el.abs() > tol.into()) - .collect::>(); - rank = aux_vec.len(); - } else { - rank = dim; - } + match rank_param { + Accuracy::Tol(tol) => { + rank = Self::rank_from_tolerance(r_mat.r_mut(), tol); + } + Accuracy::FixedRank(k) => rank = k, + Accuracy::MaxRank(tol, k) => { + rank = std::cmp::max(k, Self::rank_from_tolerance(r_mat.r_mut(), tol)); } } - let mut perm_mat = rlst_dynamic_array2!($scalar, [dim, dim]); - Self::get_p(perm_mat.r_mut(), perm); + //We permute arr to extract the columns belonging to the skeleton + let mut permutation = + rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); + permutation.set_zero(); - let mut perm_arr = - empty_array::<$scalar, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); - - for col in 0..arr.shape()[1] { - for row in 0..arr.shape()[0] { - arr.data_mut()[col * arr_qr_shape[1] + row] = - *perm_arr.get_mut([row, col]).unwrap() - } + for (index, &elem) in perm.iter().enumerate() { + *permutation.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); } + let perm_arr = empty_array::<$scalar, 2>() + .simple_mult_into_resize(permutation.r_mut(), arr.r()); + + let mut skel = empty_array(); + skel.fill_from_resize(perm_arr.into_subview([0, 0], [rank, arr_qr_shape[0]])); + //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. //If not, compute the Interpolative Decomposition matrix. if rank == 0 || rank >= dim { let mut id_mat = rlst_dynamic_array2!($scalar, [dim, dim]); id_mat.set_identity(); Ok(Self { - arr, - perm_mat, + skel, + perm, rank, id_mat, }) } else { - let shape: [usize; 2] = r_mat.shape(); - let mut id_mat: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - > = rlst_dynamic_array2!($scalar, [dim - rank, rank]); + let shape: [usize; 2] = [arr_qr_shape[1], arr_qr_shape[1]]; + let mut id_mat: DynamicArray<$scalar, 2> = + rlst_dynamic_array2!($scalar, [dim - rank, rank]); - let mut k11: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - > = rlst_dynamic_array2!($scalar, [rank, rank]); + let mut k11: DynamicArray<$scalar, 2> = + rlst_dynamic_array2!($scalar, [rank, rank]); k11.fill_from(r_mat.r_mut().into_subview([0, 0], [rank, rank])); k11.r_mut().into_inverse_alloc().unwrap(); - let mut k12: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - > = rlst_dynamic_array2!($scalar, [rank, dim - rank]); + let mut k12: DynamicArray<$scalar, 2> = + rlst_dynamic_array2!($scalar, [rank, dim - rank]); k12.fill_from( r_mat .r_mut() @@ -253,28 +225,60 @@ macro_rules! impl_id { let res = empty_array().simple_mult_into_resize(k11.r(), k12.r()); id_mat.fill_from(res.r().conj().transpose().r()); Ok(Self { - arr, - perm_mat, + skel, + perm, rank, id_mat, }) } } + fn rank_from_tolerance< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = $scalar> + + UnsafeRandomAccessByRef<2, Item = $scalar>, + >( + r_mat: Array<$scalar, ArrayImplMut, 2>, + tol: <$scalar as RlstScalar>::Real, + ) -> usize { + let dim = r_mat.shape()[0]; + let mut r_diag: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = + rlst_dynamic_array1!($scalar, [dim]); + r_mat.get_diag(r_diag.r_mut()); + let max: $scalar = r_diag + .iter() + .max_by(|a, b| a.abs().total_cmp(&b.abs())) + .unwrap() + .abs() + .into(); + + //We compute the rank of the matrix + if max.re() > 0.0 { + let alpha: $scalar = (1.0 / max) as $scalar; + r_diag.scale_inplace(alpha); + let aux_vec = r_diag + .iter() + .filter(|el| el.abs() > tol.into()) + .collect::>(); + aux_vec.len() + } else { + dim + } + } + fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + Shape<2> + UnsafeRandomAccessMut<2, Item = $scalar> + UnsafeRandomAccessByRef<2, Item = $scalar>, >( + &self, mut arr: Array<$scalar, ArrayImplMut, 2>, - perm: Vec, ) { - let m = arr.shape()[0]; - arr.set_zero(); - for col in 0..m { - arr[[col, perm[col]]] = <$scalar as One>::one(); + for (index, &elem) in self.perm.iter().enumerate() { + *arr.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); } } } From 2d186be7edefd8d0dcd82a354d1ce9b7bcc31b27 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 15:03:20 +0100 Subject: [PATCH 09/29] fixes to null space --- examples/rlst_null_space.rs | 12 +- src/dense/linalg/null_space.rs | 240 ++++++++++++--------------------- 2 files changed, 93 insertions(+), 159 deletions(-) diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index 955ce568..b6bc3f50 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -1,17 +1,15 @@ //! Demo the null space of a matrix. -use rlst::dense::linalg::null_space::NullSpaceType; pub use rlst::prelude::*; //Here we compute the null space (B) of the rowspace of a short matrix (A). pub fn main() { let mut arr = rlst_dynamic_array2!(f64, [3, 4]); arr.fill_from_seed_equally_distributed(0); + let tol = 1e-15; + let null_res = arr.r_mut().into_null_alloc(tol).unwrap(); + let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); - let null_res = arr.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); - let res = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); + println!("Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2()); - println!( - "Value of |A*B|_2, where B=null(A): {}", - res.view_flat().norm_2() - ); } + diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index b9131e1c..ebdd8c1c 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,11 +1,11 @@ //! Null space. use crate::dense::array::Array; -use crate::dense::traits::DefaultIterator; -use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue}; -use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::empty_array; -use crate::{rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, VectorContainer}; +use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; +use crate::dense::types::{RlstResult, RlstScalar}; +use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, SvdMode, VectorContainer}; use itertools::min; +use num::{One, Zero}; + /// Compute the matrix nullspace. /// @@ -14,8 +14,8 @@ use itertools::min; /// /// # Example /// -/// The following command computes the nullspace of an array `a`. -/// The nullspace is found in +/// The following command computes the nullspace of an array `a`. +/// The nullspace is found in /// ``` /// # use rlst::rlst_dynamic_array2; /// # use rlst::dense::linalg::null_space::{NullSpaceType, MatrixNull}; @@ -23,6 +23,7 @@ use itertools::min; /// # a.fill_from_seed_equally_distributed(0); /// # let null_res = a.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); /// ``` + /// This method allocates memory for the nullspace computation. pub trait MatrixNull: RlstScalar { /// Compute the matrix null space @@ -33,32 +34,24 @@ pub trait MatrixNull: RlstScalar { + RawAccessMut, >( arr: Array, - null_space_type: NullSpaceType, + tol: ::Real ) -> RlstResult>; } -macro_rules! implement_into_null { - ($scalar:ty) => { - impl MatrixNull for $scalar { - fn into_null_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array, - null_space_type: NullSpaceType, - ) -> RlstResult> { - NullSpace::<$scalar>::new(arr, null_space_type) - } - } - }; -} -implement_into_null!(f32); -implement_into_null!(f64); -implement_into_null!(c32); -implement_into_null!(c64); +impl MatrixNull for T { + fn into_null_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + arr: Array, + tol: ::Real + ) -> RlstResult> { + NullSpace::::new(arr, tol) + } +} impl< Item: RlstScalar + MatrixNull, @@ -69,137 +62,80 @@ impl< > Array { /// Compute the Column or Row nullspace of a given 2-dimensional array. - pub fn into_null_alloc(self, null_space_type: NullSpaceType) -> RlstResult> { - ::into_null_alloc(self, null_space_type) + pub fn into_null_alloc(self, tol: ::Real) -> RlstResult> { + ::into_null_alloc(self, tol) } } -///Null space -#[derive(Clone, Copy)] -#[repr(u8)] -pub enum NullSpaceType { - /// Row Nullspace - Row = b'R', - /// Column Nullspace - Column = b'C', -} -/// QR decomposition -pub struct NullSpace { - ///Row or column nullspace - pub null_space_type: NullSpaceType, +type RealScalar = ::Real; +/// Null Space +pub struct NullSpace< + Item: RlstScalar +> { ///Computed null space pub null_space_arr: Array, 2>, 2>, } -macro_rules! implement_null_space { - ($scalar:ty) => { - impl NullSpace<$scalar> { - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array<$scalar, ArrayImpl, 2>, - null_space_type: NullSpaceType, - ) -> RlstResult { - let shape = arr.shape(); - - match null_space_type { - NullSpaceType::Row => { - let mut arr_qr: Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - > = rlst_dynamic_array2!($scalar, [shape[1], shape[0]]); - arr_qr.fill_from(arr.r().conj().transpose()); - let mut null_space_arr = empty_array(); - Self::find_null_space(arr_qr, &mut null_space_arr); - Ok(Self { - null_space_type, - null_space_arr, - }) - } - NullSpaceType::Column => { - let mut null_space_arr = empty_array(); - Self::find_null_space(arr, &mut null_space_arr); - Ok(Self { - null_space_type, - null_space_arr, - }) - } - } - } - - fn find_null_space< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - arr: Array<$scalar, ArrayImpl, 2>, - null_space_arr: &mut Array< - $scalar, - BaseArray<$scalar, VectorContainer<$scalar>, 2>, - 2, - >, - ) { - let shape = arr.shape(); - let dim: usize = min(shape).unwrap(); - let qr = arr.into_qr_alloc().unwrap(); - - //We compute the QR decomposition to find a linearly independent basis of the space. - let mut q: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2> = - rlst_dynamic_array2!($scalar, [shape[0], shape[0]]); - let _ = qr.get_q_alloc(q.r_mut()); - - let mut r_mat = rlst_dynamic_array2!($scalar, [dim, dim]); - qr.get_r(r_mat.r_mut()); - - //For a full rank rectangular matrix, then rank = dim. - //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(r_mat, dim); - - null_space_arr - .fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0] - rank])); - } - - fn find_matrix_rank( - r_mat: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 2>, 2>, - dim: usize, - ) -> usize { - //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let mut r_diag: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = - rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.r_mut()); - - let max: $scalar = r_diag - .iter() - .max_by(|a, b| a.abs().total_cmp(&b.abs())) - .unwrap() - .abs() - .into(); - let rank: usize; - - if max.re() > 0.0 { - let alpha: $scalar = (1.0 / max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag - .iter() - .filter(|el| el.abs() > 1e-15) - .collect::>(); - rank = aux_vec.len(); - } else { - rank = dim; - } - rank - } +///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct +pub trait NullSpaceComputation{ + ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) + type Item: RlstScalar; + + ///We create the null space decomposition + fn new + + Stride<2> + RawAccessMut + Shape<2>> + (arr: Array, tol: RealScalar) -> RlstResult where Self: Sized; + + ///This function helps us to find the matrix rank + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:RealScalar)->usize; +} + + +impl NullSpaceComputation for NullSpace { + type Item = T; + + fn new + + Stride<2> + + RawAccessMut + + Shape<2>>(arr: Array, tol: RealScalar) -> RlstResult { + + let shape: [usize; 2] = arr.shape(); + let dim: usize = min(shape).unwrap(); + let mut singular_values: DynamicArray, 1> = rlst_dynamic_array1!(RealScalar, [dim]); + let mode: SvdMode = SvdMode::Full; + let mut u: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); + let mut vt: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); + + arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode).unwrap(); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); + + //The null space is given by the last shape[1]-rank columns of V + let mut null_space_arr: DynamicArray = empty_array(); + null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); + + Ok(Self{null_space_arr}) + + } + + fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:::Real)->usize{ + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let max: RealScalar = singular_values.r().iter().max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()).unwrap().abs(); + let mut rank: usize = dim; + + if max.re() > as Zero>::zero(){ + let alpha: RealScalar = as One>::one()/max; + singular_values.scale_inplace(alpha); + let aux_vec: Vec> = singular_values.iter().filter(|el| el.abs() > tol ).collect::>>(); + rank = aux_vec.len(); } - }; + + rank + + } } -implement_null_space!(f64); -implement_null_space!(f32); -implement_null_space!(c64); -implement_null_space!(c32); + From 02a455e2b97a27a0c2658e59bd7f5b0fd2bb874b Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 15:05:16 +0100 Subject: [PATCH 10/29] format fixes --- examples/rlst_null_space.rs | 10 +-- src/dense/linalg/null_space.rs | 116 +++++++++++++++++++++------------ 2 files changed, 79 insertions(+), 47 deletions(-) diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index b6bc3f50..a79c9bbb 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -7,9 +7,11 @@ pub fn main() { arr.fill_from_seed_equally_distributed(0); let tol = 1e-15; let null_res = arr.r_mut().into_null_alloc(tol).unwrap(); - let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); - - println!("Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2()); + let res: Array, 2>, 2> = + empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); + println!( + "Value of |A*B|_2, where B=null(A): {}", + res.view_flat().norm_2() + ); } - diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index ebdd8c1c..ba9986d1 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -1,12 +1,16 @@ //! Null space. use crate::dense::array::Array; -use crate::dense::traits::{RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, DefaultIterator}; +use crate::dense::traits::{ + DefaultIterator, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, +}; use crate::dense::types::{RlstResult, RlstScalar}; -use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, SvdMode, VectorContainer}; +use crate::{ + empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, + SvdMode, VectorContainer, +}; use itertools::min; use num::{One, Zero}; - /// Compute the matrix nullspace. /// /// The matrix nullspace is defined for a two dimensional array `arr` of @@ -14,8 +18,8 @@ use num::{One, Zero}; /// /// # Example /// -/// The following command computes the nullspace of an array `a`. -/// The nullspace is found in +/// The following command computes the nullspace of an array `a`. +/// The nullspace is found in /// ``` /// # use rlst::rlst_dynamic_array2; /// # use rlst::dense::linalg::null_space::{NullSpaceType, MatrixNull}; @@ -23,7 +27,6 @@ use num::{One, Zero}; /// # a.fill_from_seed_equally_distributed(0); /// # let null_res = a.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); /// ``` - /// This method allocates memory for the nullspace computation. pub trait MatrixNull: RlstScalar { /// Compute the matrix null space @@ -34,12 +37,11 @@ pub trait MatrixNull: RlstScalar { + RawAccessMut, >( arr: Array, - tol: ::Real + tol: ::Real, ) -> RlstResult>; } - -impl MatrixNull for T { +impl MatrixNull for T { fn into_null_alloc< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + Stride<2> @@ -47,7 +49,7 @@ impl MatrixNull for T { + RawAccessMut, >( arr: Array, - tol: ::Real + tol: ::Real, ) -> RlstResult> { NullSpace::::new(arr, tol) } @@ -67,75 +69,103 @@ impl< } } - type RealScalar = ::Real; /// Null Space -pub struct NullSpace< - Item: RlstScalar -> { +pub struct NullSpace { ///Computed null space pub null_space_arr: Array, 2>, 2>, } ///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct -pub trait NullSpaceComputation{ +pub trait NullSpaceComputation { ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) type Item: RlstScalar; ///We create the null space decomposition - fn new - + Stride<2> + RawAccessMut + Shape<2>> - (arr: Array, tol: RealScalar) -> RlstResult where Self: Sized; - + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + RawAccessMut + + Shape<2>, + >( + arr: Array, + tol: RealScalar, + ) -> RlstResult + where + Self: Sized; + ///This function helps us to find the matrix rank - fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:RealScalar)->usize; + fn find_matrix_rank( + singular_values: &mut DynamicArray, 1>, + dim: usize, + tol: RealScalar, + ) -> usize; } - -impl NullSpaceComputation for NullSpace { +impl NullSpaceComputation for NullSpace { type Item = T; - fn new - + Stride<2> - + RawAccessMut - + Shape<2>>(arr: Array, tol: RealScalar) -> RlstResult { - + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + RawAccessMut + + Shape<2>, + >( + arr: Array, + tol: RealScalar, + ) -> RlstResult { let shape: [usize; 2] = arr.shape(); let dim: usize = min(shape).unwrap(); - let mut singular_values: DynamicArray, 1> = rlst_dynamic_array1!(RealScalar, [dim]); + let mut singular_values: DynamicArray, 1> = + rlst_dynamic_array1!(RealScalar, [dim]); let mode: SvdMode = SvdMode::Full; - let mut u: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); - let mut vt: DynamicArray = rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); + let mut u: DynamicArray = + rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); + let mut vt: DynamicArray = + rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); - arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode).unwrap(); + arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode) + .unwrap(); - //For a full rank rectangular matrix, then rank = dim. + //For a full rank rectangular matrix, then rank = dim. //find_matrix_rank checks if the matrix is full rank and recomputes the rank. let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); //The null space is given by the last shape[1]-rank columns of V let mut null_space_arr: DynamicArray = empty_array(); - null_space_arr.fill_from_resize(vt.conj().transpose().into_subview([0, rank], [shape[1], shape[1]-rank])); - - Ok(Self{null_space_arr}) + null_space_arr.fill_from_resize( + vt.conj() + .transpose() + .into_subview([0, rank], [shape[1], shape[1] - rank]), + ); + Ok(Self { null_space_arr }) } - fn find_matrix_rank(singular_values: &mut DynamicArray, 1>, dim: usize, tol:::Real)->usize{ + fn find_matrix_rank( + singular_values: &mut DynamicArray, 1>, + dim: usize, + tol: ::Real, + ) -> usize { //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let max: RealScalar = singular_values.r().iter().max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()).unwrap().abs(); + let max: RealScalar = singular_values + .r() + .iter() + .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()) + .unwrap() + .abs(); let mut rank: usize = dim; - if max.re() > as Zero>::zero(){ - let alpha: RealScalar = as One>::one()/max; + if max.re() > as Zero>::zero() { + let alpha: RealScalar = as One>::one() / max; singular_values.scale_inplace(alpha); - let aux_vec: Vec> = singular_values.iter().filter(|el| el.abs() > tol ).collect::>>(); + let aux_vec: Vec> = singular_values + .iter() + .filter(|el| el.abs() > tol) + .collect::>>(); rank = aux_vec.len(); } rank - } } - - From fcaf4cb46d523a8f43135ca187fcbe0fe38425e7 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 30 Mar 2025 15:40:23 +0100 Subject: [PATCH 11/29] test nullspace update --- examples/rlst_null_space.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index a79c9bbb..7e2b4ef4 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -9,7 +9,7 @@ pub fn main() { let null_res = arr.r_mut().into_null_alloc(tol).unwrap(); let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); - + //Result println!( "Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2() From 002a4259b35bfcdb2fc36602b8b60a2be8bc0c58 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Sun, 6 Apr 2025 15:42:12 +0100 Subject: [PATCH 12/29] qr implementation for item and nullification via qr --- examples/rlst_null_space.rs | 19 +- .../linalg/interpolative_decomposition.rs | 1 + src/dense/linalg/null_space.rs | 227 ++++++++++++------ src/dense/linalg/qr.rs | 143 ++++++++--- src/prelude.rs | 2 +- 5 files changed, 284 insertions(+), 108 deletions(-) diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index 7e2b4ef4..640b36e0 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -1,4 +1,5 @@ //! Demo the null space of a matrix. +use rlst::dense::linalg::null_space::Method; pub use rlst::prelude::*; //Here we compute the null space (B) of the rowspace of a short matrix (A). @@ -6,12 +7,26 @@ pub fn main() { let mut arr = rlst_dynamic_array2!(f64, [3, 4]); arr.fill_from_seed_equally_distributed(0); let tol = 1e-15; - let null_res = arr.r_mut().into_null_alloc(tol).unwrap(); + let null_res = arr.r_mut().into_null_alloc(tol, Method::Svd).unwrap(); let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); //Result println!( - "Value of |A*B|_2, where B=null(A): {}", + "Nullification via SVD. Value of |A*B|_2, where B=null(A): {}", + res.view_flat().norm_2() + ); + + let mut arr = rlst_dynamic_array2!(f64, [3, 4]); + arr.fill_from_seed_equally_distributed(0); + let shape = arr.shape(); + let mut arr_qr = rlst_dynamic_array2!(f64, [shape[1], shape[0]]); + arr_qr.fill_from(arr.r().conj().transpose()); + let null_res = arr_qr.r_mut().into_null_alloc(tol, Method::Qr).unwrap(); + let res: Array, 2>, 2> = + empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); + //Result + println!( + "Nullification via QR. Value of |A*B|_2, where B=null(A): {}", res.view_flat().norm_2() ); } diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index f3b3920b..6e5aa9d2 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -1,5 +1,6 @@ //! Interpolative decomposition of a matrix. use crate::dense::array::Array; +use crate::dense::linalg::qr::MatrixQrDecomposition; use crate::dense::traits::{ DefaultIterator, MultIntoResize, RandomAccessMut, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index ba9986d1..8b6734eb 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -5,11 +5,11 @@ use crate::dense::traits::{ }; use crate::dense::types::{RlstResult, RlstScalar}; use crate::{ - empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixSvd, - SvdMode, VectorContainer, + empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, DynamicArray, MatrixQr, + MatrixQrDecomposition, MatrixSvd, QrDecomposition, SvdMode, VectorContainer, }; use itertools::min; -use num::{One, Zero}; +use num::Zero; /// Compute the matrix nullspace. /// @@ -38,10 +38,13 @@ pub trait MatrixNull: RlstScalar { >( arr: Array, tol: ::Real, - ) -> RlstResult>; + method: Method, + ) -> RlstResult> + where + QrDecomposition: MatrixQrDecomposition; } -impl MatrixNull for T { +impl MatrixNull for T { fn into_null_alloc< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + Stride<2> @@ -50,8 +53,12 @@ impl MatrixNull for T { >( arr: Array, tol: ::Real, - ) -> RlstResult> { - NullSpace::::new(arr, tol) + method: Method, + ) -> RlstResult> + where + QrDecomposition: MatrixQrDecomposition, + { + NullSpace::::new(arr, tol, method) } } @@ -64,8 +71,15 @@ impl< > Array { /// Compute the Column or Row nullspace of a given 2-dimensional array. - pub fn into_null_alloc(self, tol: ::Real) -> RlstResult> { - ::into_null_alloc(self, tol) + pub fn into_null_alloc( + self, + tol: ::Real, + method: Method, + ) -> RlstResult> + where + QrDecomposition: MatrixQrDecomposition, + { + ::into_null_alloc(self, tol, method) } } @@ -76,6 +90,14 @@ pub struct NullSpace { pub null_space_arr: Array, 2>, 2>, } +///Method defines the way to compute the null space +pub enum Method { + ///SVD method + Svd, + ///QR method + Qr, +} + ///NullSpaceComputation creates the null space decomposition and saves it in the NullSpace Struct pub trait NullSpaceComputation { ///This trait is implemented for RlstScalar (ie. f32, f64, c32, c64) @@ -90,19 +112,14 @@ pub trait NullSpaceComputation { >( arr: Array, tol: RealScalar, + method: Method, ) -> RlstResult where - Self: Sized; - - ///This function helps us to find the matrix rank - fn find_matrix_rank( - singular_values: &mut DynamicArray, 1>, - dim: usize, - tol: RealScalar, - ) -> usize; + Self: Sized, + QrDecomposition: MatrixQrDecomposition; } -impl NullSpaceComputation for NullSpace { +impl NullSpaceComputation for NullSpace { type Item = T; fn new< @@ -113,59 +130,131 @@ impl NullSpaceComputation for NullSpace { >( arr: Array, tol: RealScalar, - ) -> RlstResult { - let shape: [usize; 2] = arr.shape(); - let dim: usize = min(shape).unwrap(); - let mut singular_values: DynamicArray, 1> = - rlst_dynamic_array1!(RealScalar, [dim]); - let mode: SvdMode = SvdMode::Full; - let mut u: DynamicArray = - rlst_dynamic_array2!(Self::Item, [shape[0], shape[0]]); - let mut vt: DynamicArray = - rlst_dynamic_array2!(Self::Item, [shape[1], shape[1]]); - - arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode) - .unwrap(); - - //For a full rank rectangular matrix, then rank = dim. - //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = Self::find_matrix_rank(&mut singular_values, dim, tol); - - //The null space is given by the last shape[1]-rank columns of V - let mut null_space_arr: DynamicArray = empty_array(); - null_space_arr.fill_from_resize( - vt.conj() - .transpose() - .into_subview([0, rank], [shape[1], shape[1] - rank]), - ); - - Ok(Self { null_space_arr }) + method: Method, + ) -> RlstResult + where + QrDecomposition: MatrixQrDecomposition, + { + match method { + Method::Svd => { + let null_space_arr = svd_nullification(arr, tol); + Ok(Self { null_space_arr }) + } + Method::Qr => { + let null_space_arr = qr_nullification(arr, tol); + Ok(Self { null_space_arr }) + } + } } +} + +fn svd_nullification< + Item: RlstScalar + MatrixSvd, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + RawAccessMut + Shape<2>, +>( + arr: Array, + tol: RealScalar, +) -> DynamicArray { + let shape: [usize; 2] = arr.shape(); + let dim: usize = min(shape).unwrap(); + let mut singular_values: DynamicArray, 1> = + rlst_dynamic_array1!(RealScalar, [dim]); + let mode: SvdMode = SvdMode::Full; + let mut u: DynamicArray = rlst_dynamic_array2!(Item, [shape[0], shape[0]]); + let mut vt: DynamicArray = rlst_dynamic_array2!(Item, [shape[1], shape[1]]); + + arr.into_svd_alloc(u.r_mut(), vt.r_mut(), singular_values.data_mut(), mode) + .unwrap(); + + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = find_matrix_rank_svd::(&mut singular_values, dim, tol); - fn find_matrix_rank( - singular_values: &mut DynamicArray, 1>, - dim: usize, - tol: ::Real, - ) -> usize { - //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. - let max: RealScalar = singular_values - .r() + //The null space is given by the last shape[1]-rank columns of V + let mut null_space_arr: DynamicArray = empty_array(); + null_space_arr.fill_from_resize( + vt.conj() + .transpose() + .into_subview([0, rank], [shape[1], shape[1] - rank]), + ); + + null_space_arr +} + +fn qr_nullification< + Item: RlstScalar + MatrixQr, + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> + Shape<2> + RawAccessMut, +>( + arr: Array, + tol: RealScalar, +) -> DynamicArray +where + QrDecomposition: MatrixQrDecomposition, +{ + let shape = arr.shape(); + let dim: usize = min(shape).unwrap(); + let qr = arr.into_qr_alloc().unwrap(); + //We compute the QR decomposition to find a linearly independent basis of the space. + let mut q = rlst_dynamic_array2!(Item, [shape[0], shape[0]]); + let _ = qr.get_q_alloc(q.r_mut()); + let mut r_mat: DynamicArray = rlst_dynamic_array2!(Item, [dim, dim]); + qr.get_r(r_mat.r_mut()); + //For a full rank rectangular matrix, then rank = dim. + //find_matrix_rank checks if the matrix is full rank and recomputes the rank. + let rank: usize = find_matrix_rank_qr::(&r_mat, dim, tol); + let mut null_space_arr = empty_array(); + null_space_arr.fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0] - rank])); + null_space_arr +} + +fn find_matrix_rank_svd( + singular_values: &mut DynamicArray, 1>, + dim: usize, + tol: ::Real, +) -> usize { + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let max: RealScalar = singular_values + .r() + .iter() + .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()) + .unwrap() + .abs(); + + let rank: usize = if max.re() > as Zero>::zero() { + let aux_vec: Vec> = singular_values .iter() - .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()) - .unwrap() - .abs(); - let mut rank: usize = dim; - - if max.re() > as Zero>::zero() { - let alpha: RealScalar = as One>::one() / max; - singular_values.scale_inplace(alpha); - let aux_vec: Vec> = singular_values - .iter() - .filter(|el| el.abs() > tol) - .collect::>>(); - rank = aux_vec.len(); - } + .filter(|el| el.abs() > tol * max) + .collect(); + aux_vec.len() + } else { + dim + }; - rank - } + rank +} + +fn find_matrix_rank_qr( + r_mat: &DynamicArray, + dim: usize, + tol: ::Real, +) -> usize { + //We compute the rank of the matrix by expecting the values of the elements in the diagonal of R. + let mut r_diag = rlst_dynamic_array1!(Item, [dim]); + r_mat.get_diag(r_diag.r_mut()); + + let max: RealScalar = r_diag + .r() + .iter() + .max_by(|a, b| (a.abs().partial_cmp(&b.abs())).unwrap()) + .unwrap() + .abs(); + + let rank: usize = if max.re() > as Zero>::zero() { + let aux_vec: Vec = r_diag.iter().filter(|el| el.abs() > tol * max).collect(); + aux_vec.len() + } else { + dim + }; + + rank } diff --git a/src/dense/linalg/qr.rs b/src/dense/linalg/qr.rs index f932b864..3c6e200f 100644 --- a/src/dense/linalg/qr.rs +++ b/src/dense/linalg/qr.rs @@ -90,6 +90,94 @@ pub struct QrDecomposition< jpvt: Vec, } +/// Compute the QR decomposition of a matrix. +/// +/// The QR Decomposition of an `(m,n)` matrix `A` is defined +/// by `AP = QR`, where `P` is an `(n, n)` permutation matrix, +/// `Q` is a `(m, m)` orthogonal matrix, and `R` is +/// an `(m, n)` upper triangular matrix. +pub trait MatrixQrDecomposition: Sized { + /// Item type + type Item: RlstScalar; + /// Array implementaion + type ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + RawAccessMut + + Shape<2>; + + /// Create a new QR Decomposition. + fn new(arr: Array) -> RlstResult; + + /// Return the permuation vector for the QR decomposition. + /// + /// If `perm[i] = j` then the ith column of QR corresponds + /// to the jth column of the original array. + fn get_perm(&self) -> Vec; + + /// Return the R matrix of the QR decomposition. + /// + /// If `A` has dimension `(m, n)` then R has + /// dimension `(k, n)` with `k=min(m, n)`. + fn get_r< + ArrayImplR: UnsafeRandomAccessByValue<2, Item = Self::Item> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + RawAccessMut + + Shape<2> + + Stride<2>, + >( + &self, + arr: Array, + ); + + /// Return the permuation matrix `P`. + /// + /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`. + fn get_p< + ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + RawAccessMut + + Shape<2> + + Stride<2>, + >( + &self, + arr: Array, + ); + + /// Return the Q matrix of the QR decomposition. + /// + /// If `A` has dimension `(m, n)` then `arr` needs + /// to be of dimension `(m, r)`, where `r<= m`` + /// is the desired number of columns of `Q`. + /// + /// This method allocates temporary memory during execution. + fn get_q_alloc< + ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + RawAccessMut + + Shape<2> + + Stride<2>, + >( + &self, + arr: Array, + ) -> RlstResult<()>; + + /// Apply Q to a given matrix. + /// + /// This method allocates temporary memory during execution. + fn apply_q_alloc< + ArrayImplQ: UnsafeRandomAccessByValue<2, Item = Self::Item> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + RawAccessMut + + Shape<2> + + Stride<2>, + >( + &self, + arr: Array, + side: ApplyQSide, + trans: ApplyQTrans, + ) -> RlstResult<()>; +} + macro_rules! implement_qr_real { ($scalar:ty, $geqp3:expr, $ormqr:expr) => { impl< @@ -97,10 +185,12 @@ macro_rules! implement_qr_real { + Stride<2> + Shape<2> + RawAccessMut, - > QrDecomposition<$scalar, ArrayImpl> + > MatrixQrDecomposition for QrDecomposition<$scalar, ArrayImpl> { - /// Create a new QR Decomposition. - pub fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult { + type Item = $scalar; + type ArrayImpl = ArrayImpl; + + fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult { let stride = arr.stride(); let shape = arr.shape(); @@ -165,22 +255,14 @@ macro_rules! implement_qr_real { } } - /// Return the permuation vector for the QR decomposition. - /// - /// If `perm[i] = j` then the ith column of QR corresponds - /// to the jth column of the original array. - pub fn get_perm(&self) -> Vec { + fn get_perm(&self) -> Vec { self.jpvt .iter() .map(|&elem| elem as usize - 1) .collect_vec() } - /// Return the R matrix of the QR decomposition. - /// - /// If `A` has dimension `(m, n)` then R has - /// dimension `(k, n)` with `k=min(m, n)`. - pub fn get_r< + fn get_r< ArrayImplR: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -205,10 +287,7 @@ macro_rules! implement_qr_real { } } - /// Return the permuation matrix `P`. - /// - /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`. - pub fn get_p< + fn get_p< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -226,14 +305,7 @@ macro_rules! implement_qr_real { } } - /// Return the Q matrix of the QR decomposition. - /// - /// If `A` has dimension `(m, n)` then `arr` needs - /// to be of dimension `(m, r)`, where `r<= m`` - /// is the desired number of columns of `Q`. - /// - /// This method allocates temporary memory during execution. - pub fn get_q_alloc< + fn get_q_alloc< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -249,10 +321,7 @@ macro_rules! implement_qr_real { self.apply_q_alloc(arr, ApplyQSide::Left, ApplyQTrans::NoTrans) } - /// Apply Q to a given matrix. - /// - /// This method allocates temporary memory during execution. - pub fn apply_q_alloc< + fn apply_q_alloc< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -360,10 +429,12 @@ macro_rules! implement_qr_complex { + Stride<2> + Shape<2> + RawAccessMut, - > QrDecomposition<$scalar, ArrayImpl> + > MatrixQrDecomposition for QrDecomposition<$scalar, ArrayImpl> { + type Item = $scalar; + type ArrayImpl = ArrayImpl; /// Create a new QR Decomposition. - pub fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult { + fn new(mut arr: Array<$scalar, ArrayImpl, 2>) -> RlstResult { let stride = arr.stride(); let shape = arr.shape(); @@ -437,7 +508,7 @@ macro_rules! implement_qr_complex { /// /// If `perm[i] = j` then the ith column of QR corresponds /// to the jth column of the original array. - pub fn get_perm(&self) -> Vec { + fn get_perm(&self) -> Vec { self.jpvt .iter() .map(|&elem| elem as usize - 1) @@ -448,7 +519,7 @@ macro_rules! implement_qr_complex { /// /// If `A` has dimension `(m, n)` then R has /// dimension `(k, n)` with `k=min(m, n)`. - pub fn get_r< + fn get_r< ArrayImplR: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -476,7 +547,7 @@ macro_rules! implement_qr_complex { /// Return the permuation matrix `P`. /// /// For `A` an `(m,n)` matrix `P` has dimension `(n, n)`. - pub fn get_p< + fn get_p< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -501,7 +572,7 @@ macro_rules! implement_qr_complex { /// is the desired number of columns of `Q`. /// /// This method allocates temporary memory during execution. - pub fn get_q_alloc< + fn get_q_alloc< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut @@ -520,7 +591,7 @@ macro_rules! implement_qr_complex { /// Apply Q to a given matrix. /// /// This method allocates temporary memory during execution. - pub fn apply_q_alloc< + fn apply_q_alloc< ArrayImplQ: UnsafeRandomAccessByValue<2, Item = $scalar> + UnsafeRandomAccessMut<2, Item = $scalar> + RawAccessMut diff --git a/src/prelude.rs b/src/prelude.rs index 8c6dbc5a..7067534d 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -66,7 +66,7 @@ pub use crate::dense::linalg::inverse::MatrixInverse; pub use crate::dense::linalg::lu::{LuDecomposition, MatrixLuDecomposition}; pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; pub use crate::dense::linalg::pseudo_inverse::MatrixPseudoInverse; -pub use crate::dense::linalg::qr::{MatrixQr, QrDecomposition}; +pub use crate::dense::linalg::qr::{MatrixQr, MatrixQrDecomposition, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; pub use crate::dense::array::rank1_array::Rank1Array; From e0aa57b5b73a8af93f27ddc5003de2c5a9d9dcc1 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Thu, 17 Apr 2025 15:50:51 +0100 Subject: [PATCH 13/29] use of upper triangular matrices in id --- examples/rlst_interpolative_decomposition.rs | 5 +- .../linalg/interpolative_decomposition.rs | 229 ++++++++++++++---- 2 files changed, 191 insertions(+), 43 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index e779de7a..91cf0583 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -48,7 +48,7 @@ pub fn main() { let n: usize = 100; let slice: usize = 50; //Tolerance given for the - let tol: f64 = 1e-5; + let tol: f64 = 1e-4; //Create a low rank matrix let mut arr: DynamicArray = rlst_dynamic_array2!(f64, [slice, n]); @@ -64,12 +64,13 @@ pub fn main() { println!("The rank of this matrix is {}\n", res.rank); + println!("12"); //We extract the residuals of the matrix let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); res.get_p(perm_mat.r_mut()); let perm_arr: DynamicArray = empty_array::().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); - + println!("13"); let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice - res.rank, n]); a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n])); diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 6e5aa9d2..e2e497b1 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -1,14 +1,17 @@ //! Interpolative decomposition of a matrix. use crate::dense::array::Array; use crate::dense::linalg::qr::MatrixQrDecomposition; +use crate::dense::traits::ResizeInPlace; use crate::dense::traits::{ - DefaultIterator, MultIntoResize, RandomAccessMut, RawAccessMut, Shape, Stride, - UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, + DefaultIterator, MultIntoResize, RandomAccessByRef, RandomAccessMut, RawAccessMut, Shape, + Stride, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::DynamicArray; -use crate::{empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, VectorContainer}; - +use crate::{ + empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, RlstError, VectorContainer, +}; +use crate::{DynamicArray, RawAccess}; +use blas::{ctrsm, dtrsm, strsm, ztrsm}; /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. /// /// The matrix interpolative decomposition is defined for a two dimensional 'long' array `arr` of @@ -29,6 +32,7 @@ pub trait MatrixId: RlstScalar { ///This method allocates space for ID fn into_id_alloc< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + UnsafeRandomAccessMut<2, Item = Self> + Stride<2> + Shape<2> + RawAccessMut, @@ -43,6 +47,7 @@ macro_rules! implement_into_id { impl MatrixId for $scalar { fn into_id_alloc< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + UnsafeRandomAccessMut<2, Item = Self> + Stride<2> + Shape<2> + RawAccessMut, @@ -64,6 +69,7 @@ implement_into_id!(c64); impl< Item: RlstScalar + MatrixId, ArrayImplId: UnsafeRandomAccessByValue<2, Item = Item> + + UnsafeRandomAccessMut<2, Item = Item> + Stride<2> + RawAccessMut + Shape<2>, @@ -85,6 +91,7 @@ pub trait MatrixIdDecomposition: Sized { /// Create a new Interpolative Decomposition from a given array. fn new< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + UnsafeRandomAccessMut<2, Item = Self::Item> + Stride<2> + Shape<2> + RawAccessMut, @@ -100,7 +107,7 @@ pub trait MatrixIdDecomposition: Sized { + UnsafeRandomAccessMut<2, Item = Self::Item> + UnsafeRandomAccessByRef<2, Item = Self::Item>, >( - r_diag: Array, + ut_mat: Array, tol: ::Real, ) -> usize; @@ -138,6 +145,149 @@ pub enum Accuracy { MaxRank(T, usize), } +///Interface to obtain the upper-triangular part of the matrix +pub trait UpperTriangular: RlstScalar { + ///Compute the upper triangular + fn into_upper_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self>, + >( + arr: Array, + ) -> RlstResult>; +} + +macro_rules! implement_into_upper_triangular { + ($scalar:ty) => { + impl UpperTriangular for $scalar { + fn into_upper_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self>, + >( + arr: Array, + ) -> RlstResult> { + UpperTriangularMatrix::<$scalar>::new(arr) + } + } + }; +} + +implement_into_upper_triangular!(f32); +implement_into_upper_triangular!(f64); +implement_into_upper_triangular!(c32); +implement_into_upper_triangular!(c64); + +/// A struct representing an upper triangular matrix. +pub struct UpperTriangularMatrix { + /// The upper triangular part of the matrix. + pub upper: DynamicArray, +} + +///Solver for upper-triangular matrix. +pub trait UppperTriangularOperations: Sized { + ///Defines an abstract type for the matrix to solve + type Item: RlstScalar; + ///Extract the upper triangular and save it + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self::Item>, + >( + arr: Array, + ) -> RlstResult; + ///Solves the upper-triangular system + fn solve_upper_triangular< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array, + ); +} + +macro_rules! implement_solve_upper_triangular { + ($scalar:ty, $trsm:expr) => { + impl UppperTriangularOperations for UpperTriangularMatrix<$scalar> { + type Item = $scalar; + + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self::Item>, + >( + arr: Array, + ) -> RlstResult { + let shape = arr.shape(); + let mut upper = rlst_dynamic_array2!(Self::Item, shape); + for i in 0..shape[0] { + for j in i..shape[1] { + *upper.get_mut([i, j]).unwrap() = *arr.get([i, j]).unwrap(); + } + } + + println!("Diff: {}", (arr.r() - upper.r()).norm_1()); + + if (arr.r() - upper.r()).norm_1() < 1e-15 { + Ok(Self { upper }) + } else { + Err(RlstError::OperationFailed( + "Not upper triangular".to_string(), + )) + } + } + + fn solve_upper_triangular< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array<$scalar, ArrayImpl, 2>, + ) { + let lda = self.upper.stride()[1] as i32; + let m = arr_b.shape()[0] as i32; + let n = arr_b.shape()[1] as i32; + let ldb = arr_b.stride()[1] as i32; + let alpha = 1.0; + + unsafe { + // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb) + $trsm( + b'L', // Side: A is on the left + b'U', // Uplo: A is upper triangular + b'N', // TransA: no transpose + b'N', // Diag: A is not unit triangular + m, + n, + alpha.into(), + self.upper.data(), + lda, + arr_b.data_mut(), + ldb, + ); + } + } + } + }; +} + +implement_solve_upper_triangular!(f32, strsm); +implement_solve_upper_triangular!(f64, dtrsm); +implement_solve_upper_triangular!(c32, ctrsm); +implement_solve_upper_triangular!(c64, ztrsm); + macro_rules! impl_id { ($scalar:ty) => { impl MatrixIdDecomposition for IdDecomposition<$scalar> { @@ -145,6 +295,7 @@ macro_rules! impl_id { fn new< ArrayImpl: UnsafeRandomAccessByValue<2, Item = $scalar> + + UnsafeRandomAccessMut<2, Item = Self::Item> + Stride<2> + Shape<2> + RawAccessMut, @@ -153,48 +304,47 @@ macro_rules! impl_id { rank_param: Accuracy<<$scalar as RlstScalar>::Real>, ) -> RlstResult { //We compute the QR decomposition using rlst QR decomposition - let mut arr_qr: DynamicArray<$scalar, 2> = - rlst_dynamic_array2!($scalar, [arr.shape()[1], arr.shape()[0]]); - arr_qr.fill_from(arr.r().conj().transpose()); - let arr_qr_shape = arr_qr.shape(); - let qr = arr_qr.r_mut().into_qr_alloc().unwrap(); + let mut arr_work = empty_array(); + let mut u_tri = empty_array(); + arr_work.fill_from_resize(arr.r().conj().transpose()); + let shape = arr_work.shape(); + u_tri.resize_in_place([shape[1], shape[1]]); + let dim = shape[1]; + let qr = arr_work.r_mut().into_qr_alloc().unwrap(); //We obtain relevant parameters of the decomposition: the permutation induced by the pivoting and the R matrix let perm = qr.get_perm(); - let mut r_mat = rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); - qr.get_r(r_mat.r_mut()); + qr.get_r(u_tri.r_mut()); //The maximum rank is given by the number of columns of the transposed matrix - let dim: usize = arr_qr_shape[1]; let rank: usize; //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. - match rank_param { Accuracy::Tol(tol) => { - rank = Self::rank_from_tolerance(r_mat.r_mut(), tol); + rank = Self::rank_from_tolerance(u_tri.r_mut(), tol); } Accuracy::FixedRank(k) => rank = k, Accuracy::MaxRank(tol, k) => { - rank = std::cmp::max(k, Self::rank_from_tolerance(r_mat.r_mut(), tol)); + rank = std::cmp::max(k, Self::rank_from_tolerance(u_tri.r_mut(), tol)); } } - //We permute arr to extract the columns belonging to the skeleton - let mut permutation = - rlst_dynamic_array2!($scalar, [arr_qr_shape[1], arr_qr_shape[1]]); + let mut permutation = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); permutation.set_zero(); for (index, &elem) in perm.iter().enumerate() { *permutation.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); } - - let perm_arr = empty_array::<$scalar, 2>() + let mut perm_arr = empty_array::<$scalar, 2>(); + perm_arr + .r_mut() .simple_mult_into_resize(permutation.r_mut(), arr.r()); - let mut skel = empty_array(); - skel.fill_from_resize(perm_arr.into_subview([0, 0], [rank, arr_qr_shape[0]])); + //We permute arr to extract the columns belonging to the skeleton + let mut skel = empty_array(); + skel.fill_from_resize(perm_arr.into_subview([0, 0], [rank, shape[0]])); //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. //If not, compute the Interpolative Decomposition matrix. if rank == 0 || rank >= dim { @@ -207,24 +357,21 @@ macro_rules! impl_id { id_mat, }) } else { - let shape: [usize; 2] = [arr_qr_shape[1], arr_qr_shape[1]]; + let shape: [usize; 2] = [shape[1], shape[1]]; let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim - rank, rank]); + println!("9"); + let r11 = UpperTriangularMatrix::<$scalar>::new( + u_tri.r_mut().into_subview([0, 0], [rank, rank]), + ) + .unwrap(); - let mut k11: DynamicArray<$scalar, 2> = - rlst_dynamic_array2!($scalar, [rank, rank]); - k11.fill_from(r_mat.r_mut().into_subview([0, 0], [rank, rank])); - k11.r_mut().into_inverse_alloc().unwrap(); - let mut k12: DynamicArray<$scalar, 2> = - rlst_dynamic_array2!($scalar, [rank, dim - rank]); - k12.fill_from( - r_mat - .r_mut() - .into_subview([0, rank], [rank, shape[1] - rank]), - ); + let mut r12 = u_tri + .r_mut() + .into_subview([0, rank], [rank, shape[1] - rank]); + r11.solve_upper_triangular(&mut r12); - let res = empty_array().simple_mult_into_resize(k11.r(), k12.r()); - id_mat.fill_from(res.r().conj().transpose().r()); + id_mat.fill_from(r12.r().conj().transpose().r()); Ok(Self { skel, perm, @@ -240,13 +387,13 @@ macro_rules! impl_id { + UnsafeRandomAccessMut<2, Item = $scalar> + UnsafeRandomAccessByRef<2, Item = $scalar>, >( - r_mat: Array<$scalar, ArrayImplMut, 2>, + ut_mat: Array<$scalar, ArrayImplMut, 2>, tol: <$scalar as RlstScalar>::Real, ) -> usize { - let dim = r_mat.shape()[0]; + let dim = ut_mat.shape()[0]; let mut r_diag: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = rlst_dynamic_array1!($scalar, [dim]); - r_mat.get_diag(r_diag.r_mut()); + ut_mat.get_diag(r_diag.r_mut()); let max: $scalar = r_diag .iter() .max_by(|a, b| a.abs().total_cmp(&b.abs())) From f618500f8cdc34ab108c4b5becfeec6558afe955 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Thu, 17 Apr 2025 15:53:47 +0100 Subject: [PATCH 14/29] remove message --- src/dense/linalg/interpolative_decomposition.rs | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index e2e497b1..04c0dc45 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -236,8 +236,6 @@ macro_rules! implement_solve_upper_triangular { } } - println!("Diff: {}", (arr.r() - upper.r()).norm_1()); - if (arr.r() - upper.r()).norm_1() < 1e-15 { Ok(Self { upper }) } else { @@ -360,7 +358,6 @@ macro_rules! impl_id { let shape: [usize; 2] = [shape[1], shape[1]]; let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim - rank, rank]); - println!("9"); let r11 = UpperTriangularMatrix::<$scalar>::new( u_tri.r_mut().into_subview([0, 0], [rank, rank]), ) From f09dd2fd8d5b675514441f9151cd2abc62f48479 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 18 Apr 2025 11:21:13 +0100 Subject: [PATCH 15/29] rank computation simplified --- examples/rlst_interpolative_decomposition.rs | 2 - .../linalg/interpolative_decomposition.rs | 59 +++++++++---------- 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 91cf0583..590db5d6 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -64,13 +64,11 @@ pub fn main() { println!("The rank of this matrix is {}\n", res.rank); - println!("12"); //We extract the residuals of the matrix let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); res.get_p(perm_mat.r_mut()); let perm_arr: DynamicArray = empty_array::().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); - println!("13"); let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice - res.rank, n]); a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n])); diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 04c0dc45..6d2a2890 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -3,13 +3,11 @@ use crate::dense::array::Array; use crate::dense::linalg::qr::MatrixQrDecomposition; use crate::dense::traits::ResizeInPlace; use crate::dense::traits::{ - DefaultIterator, MultIntoResize, RandomAccessByRef, RandomAccessMut, RawAccessMut, Shape, - Stride, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, + MultIntoResize, RandomAccessByRef, RandomAccessMut, RawAccessMut, Shape, Stride, + UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; -use crate::{ - empty_array, rlst_dynamic_array1, rlst_dynamic_array2, BaseArray, RlstError, VectorContainer, -}; +use crate::{empty_array, rlst_dynamic_array2, BaseArray, VectorContainer}; use crate::{DynamicArray, RawAccess}; use blas::{ctrsm, dtrsm, strsm, ztrsm}; /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. @@ -156,6 +154,7 @@ pub trait UpperTriangular: RlstScalar { + RandomAccessByRef<2, Item = Self>, >( arr: Array, + extract: bool, ) -> RlstResult>; } @@ -170,8 +169,9 @@ macro_rules! implement_into_upper_triangular { + RandomAccessByRef<2, Item = Self>, >( arr: Array, + extract: bool, ) -> RlstResult> { - UpperTriangularMatrix::<$scalar>::new(arr) + UpperTriangularMatrix::<$scalar>::new(arr, extract) } } }; @@ -201,6 +201,7 @@ pub trait UppperTriangularOperations: Sized { + RandomAccessByRef<2, Item = Self::Item>, >( arr: Array, + extract: bool, ) -> RlstResult; ///Solves the upper-triangular system fn solve_upper_triangular< @@ -227,21 +228,22 @@ macro_rules! implement_solve_upper_triangular { + RandomAccessByRef<2, Item = Self::Item>, >( arr: Array, + extract: bool, ) -> RlstResult { let shape = arr.shape(); - let mut upper = rlst_dynamic_array2!(Self::Item, shape); - for i in 0..shape[0] { - for j in i..shape[1] { - *upper.get_mut([i, j]).unwrap() = *arr.get([i, j]).unwrap(); - } - } - if (arr.r() - upper.r()).norm_1() < 1e-15 { + if extract { + let mut upper = rlst_dynamic_array2!(Self::Item, shape); + for i in 0..shape[0] { + for j in i..shape[1] { + *upper.get_mut([i, j]).unwrap() = *arr.get([i, j]).unwrap(); + } + } Ok(Self { upper }) } else { - Err(RlstError::OperationFailed( - "Not upper triangular".to_string(), - )) + let mut upper = empty_array::<$scalar, 2>(); + upper.fill_from_resize(arr.r()); + Ok(Self { upper }) } } @@ -360,6 +362,7 @@ macro_rules! impl_id { rlst_dynamic_array2!($scalar, [dim - rank, rank]); let r11 = UpperTriangularMatrix::<$scalar>::new( u_tri.r_mut().into_subview([0, 0], [rank, rank]), + false, ) .unwrap(); @@ -388,25 +391,17 @@ macro_rules! impl_id { tol: <$scalar as RlstScalar>::Real, ) -> usize { let dim = ut_mat.shape()[0]; - let mut r_diag: Array<$scalar, BaseArray<$scalar, VectorContainer<$scalar>, 1>, 1> = - rlst_dynamic_array1!($scalar, [dim]); - ut_mat.get_diag(r_diag.r_mut()); - let max: $scalar = r_diag - .iter() - .max_by(|a, b| a.abs().total_cmp(&b.abs())) - .unwrap() - .abs() - .into(); + let max = ut_mat.get([0, 0]).unwrap().abs(); //We compute the rank of the matrix if max.re() > 0.0 { - let alpha: $scalar = (1.0 / max) as $scalar; - r_diag.scale_inplace(alpha); - let aux_vec = r_diag - .iter() - .filter(|el| el.abs() > tol.into()) - .collect::>(); - aux_vec.len() + let mut rank = 0; + for i in 0..dim { + if ut_mat.get([i, i]).unwrap().abs() > tol * max { + rank += 1; + } + } + rank } else { dim } From 3cd677f4fe3cf750fa977377a189e16043518aa8 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 18 Apr 2025 11:32:38 +0100 Subject: [PATCH 16/29] improved memory management --- .../linalg/interpolative_decomposition.rs | 37 ++++++++----------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 6d2a2890..af0e24be 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -3,8 +3,8 @@ use crate::dense::array::Array; use crate::dense::linalg::qr::MatrixQrDecomposition; use crate::dense::traits::ResizeInPlace; use crate::dense::traits::{ - MultIntoResize, RandomAccessByRef, RandomAccessMut, RawAccessMut, Shape, Stride, - UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, + MultIntoResize, RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, + UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; use crate::{empty_array, rlst_dynamic_array2, BaseArray, VectorContainer}; @@ -154,7 +154,6 @@ pub trait UpperTriangular: RlstScalar { + RandomAccessByRef<2, Item = Self>, >( arr: Array, - extract: bool, ) -> RlstResult>; } @@ -169,9 +168,8 @@ macro_rules! implement_into_upper_triangular { + RandomAccessByRef<2, Item = Self>, >( arr: Array, - extract: bool, ) -> RlstResult> { - UpperTriangularMatrix::<$scalar>::new(arr, extract) + UpperTriangularMatrix::<$scalar>::new(arr) } } }; @@ -201,7 +199,6 @@ pub trait UppperTriangularOperations: Sized { + RandomAccessByRef<2, Item = Self::Item>, >( arr: Array, - extract: bool, ) -> RlstResult; ///Solves the upper-triangular system fn solve_upper_triangular< @@ -228,23 +225,16 @@ macro_rules! implement_solve_upper_triangular { + RandomAccessByRef<2, Item = Self::Item>, >( arr: Array, - extract: bool, ) -> RlstResult { let shape = arr.shape(); - - if extract { - let mut upper = rlst_dynamic_array2!(Self::Item, shape); - for i in 0..shape[0] { - for j in i..shape[1] { - *upper.get_mut([i, j]).unwrap() = *arr.get([i, j]).unwrap(); - } + let mut upper = rlst_dynamic_array2!(Self::Item, shape); + let mut view = upper.r_mut(); + for i in 0..shape[0] { + for j in i..shape[1] { + view[[i, j]] = arr[[i, j]]; } - Ok(Self { upper }) - } else { - let mut upper = empty_array::<$scalar, 2>(); - upper.fill_from_resize(arr.r()); - Ok(Self { upper }) } + Ok(Self { upper }) } fn solve_upper_triangular< @@ -333,8 +323,10 @@ macro_rules! impl_id { let mut permutation = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); permutation.set_zero(); + let mut view = permutation.r_mut(); + for (index, &elem) in perm.iter().enumerate() { - *permutation.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); + view[[index, elem]] = <$scalar as num::One>::one(); } let mut perm_arr = empty_array::<$scalar, 2>(); perm_arr @@ -362,7 +354,6 @@ macro_rules! impl_id { rlst_dynamic_array2!($scalar, [dim - rank, rank]); let r11 = UpperTriangularMatrix::<$scalar>::new( u_tri.r_mut().into_subview([0, 0], [rank, rank]), - false, ) .unwrap(); @@ -417,8 +408,10 @@ macro_rules! impl_id { mut arr: Array<$scalar, ArrayImplMut, 2>, ) { arr.set_zero(); + let mut view = arr.r_mut(); + for (index, &elem) in self.perm.iter().enumerate() { - *arr.get_mut([index, elem]).unwrap() = <$scalar as num::One>::one(); + view[[index, elem]] = <$scalar as num::One>::one(); } } } From 981834b9385ce7ca466d884ed86700d4e5935152 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Mon, 28 Apr 2025 11:25:20 +0100 Subject: [PATCH 17/29] transpose mode for interpolative decomposition --- examples/rlst_interpolative_decomposition.rs | 55 ++++++++++++++++--- .../linalg/interpolative_decomposition.rs | 33 ++++++++--- 2 files changed, 73 insertions(+), 15 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 590db5d6..24ee4470 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -44,17 +44,45 @@ fn low_rank_matrix(n: usize, arr: &mut Array = rlst_dynamic_array2!(f64, [slice, n]); + low_rank_matrix(n, &mut arr); - //Create a low rank matrix + let res: IdDecomposition = arr.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap(); + + println!("The skeleton of the matrix is given by"); + res.skel.pretty_print(); + + println!("The interpolative decomposition matrix is:"); + res.id_mat.pretty_print(); + + println!("The rank of this matrix is {}\n", res.rank); + + //We extract the residuals of the matrix + let mut perm_mat: DynamicArray = rlst_dynamic_array2!(f64, [slice, slice]); + res.get_p(perm_mat.r_mut()); + let perm_arr: DynamicArray = + empty_array::().simple_mult_into_resize(perm_mat.r_mut(), arr.r()); + let mut a_rs: DynamicArray = rlst_dynamic_array2!(f64, [slice - res.rank, n]); + a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n])); + + //We compute an approximation of the residual columns of the matrix + let a_rs_app: DynamicArray = + empty_array().simple_mult_into_resize(res.id_mat.r(), res.skel); + + let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); + println!("Interpolative Decomposition L2 absolute error: {}", error); + +} + +fn test_skinny_matrix(slice: usize, n: usize, tol: f64){ let mut arr: DynamicArray = rlst_dynamic_array2!(f64, [slice, n]); low_rank_matrix(n, &mut arr); - let res: IdDecomposition = arr.r_mut().into_id_alloc(Accuracy::Tol(tol)).unwrap(); + let mut arr_trans = empty_array(); + arr_trans.fill_from_resize(arr.r().transpose()); + + let res: IdDecomposition = arr_trans.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::Trans).unwrap(); println!("The skeleton of the matrix is given by"); res.skel.pretty_print(); @@ -78,4 +106,17 @@ pub fn main() { let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); + +} + +pub fn main() { + let n: usize = 100; + let slice: usize = 50; + //Tolerance given for the + let tol: f64 = 1e-4; + + //Create a low rank matrix + test_fat_matrix(slice, n, tol); + test_skinny_matrix(slice, n, tol); + } diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index af0e24be..8710bfee 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -6,6 +6,7 @@ use crate::dense::traits::{ MultIntoResize, RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; +use crate::dense::types::TransMode; // Import TransMode from the appropriate module use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; use crate::{empty_array, rlst_dynamic_array2, BaseArray, VectorContainer}; use crate::{DynamicArray, RawAccess}; @@ -37,6 +38,7 @@ pub trait MatrixId: RlstScalar { >( arr: Array, rank_param: Accuracy<::Real>, + trans_mode: TransMode, ) -> RlstResult>; } @@ -52,8 +54,9 @@ macro_rules! implement_into_id { >( arr: Array, rank_param: Accuracy<::Real>, + trans_mode: TransMode, ) -> RlstResult> { - IdDecomposition::<$scalar>::new(arr, rank_param) + IdDecomposition::<$scalar>::new(arr, rank_param, trans_mode) } } }; @@ -77,8 +80,9 @@ impl< pub fn into_id_alloc( self, rank_param: Accuracy<::Real>, + trans_mode: TransMode, ) -> RlstResult> { - ::into_id_alloc(self, rank_param) + ::into_id_alloc(self, rank_param, trans_mode) } } @@ -96,6 +100,7 @@ pub trait MatrixIdDecomposition: Sized { >( arr: Array, rank_param: Accuracy<::Real>, + trans_mode: TransMode, ) -> RlstResult; /// Compute the rank of the decomposition from a given tolerance @@ -292,11 +297,19 @@ macro_rules! impl_id { >( arr: Array<$scalar, ArrayImpl, 2>, rank_param: Accuracy<<$scalar as RlstScalar>::Real>, + trans_mode: TransMode, ) -> RlstResult { //We compute the QR decomposition using rlst QR decomposition let mut arr_work = empty_array(); let mut u_tri = empty_array(); - arr_work.fill_from_resize(arr.r().conj().transpose()); + + match trans_mode { + TransMode::Trans => arr_work.fill_from_resize(arr.r().conj()), + TransMode::NoTrans => arr_work.fill_from_resize(arr.r().conj().transpose()), + TransMode::ConjNoTrans => arr_work.fill_from_resize(arr.r()), + TransMode::ConjTrans => arr_work.fill_from_resize(arr.r().transpose()), + }; + let shape = arr_work.shape(); u_tri.resize_in_place([shape[1], shape[1]]); let dim = shape[1]; @@ -324,17 +337,21 @@ macro_rules! impl_id { permutation.set_zero(); let mut view = permutation.r_mut(); - for (index, &elem) in perm.iter().enumerate() { view[[index, elem]] = <$scalar as num::One>::one(); } + let mut perm_arr = empty_array::<$scalar, 2>(); - perm_arr - .r_mut() - .simple_mult_into_resize(permutation.r_mut(), arr.r()); + perm_arr.r_mut().mult_into_resize( + TransMode::NoTrans, + trans_mode, + num::One::one(), + permutation.r_mut(), + arr.r(), + num::Zero::zero(), + ); //We permute arr to extract the columns belonging to the skeleton - let mut skel = empty_array(); skel.fill_from_resize(perm_arr.into_subview([0, 0], [rank, shape[0]])); //In the case the matrix is full rank or we get a matrix of rank 0, then return the identity matrix. From 8e85986adbba01ad63b0adcaacec6e76bff59e01 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Mon, 28 Apr 2025 11:27:51 +0100 Subject: [PATCH 18/29] format --- examples/rlst_interpolative_decomposition.rs | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/rlst_interpolative_decomposition.rs b/examples/rlst_interpolative_decomposition.rs index 24ee4470..d76fd1f8 100644 --- a/examples/rlst_interpolative_decomposition.rs +++ b/examples/rlst_interpolative_decomposition.rs @@ -44,11 +44,14 @@ fn low_rank_matrix(n: usize, arr: &mut Array = rlst_dynamic_array2!(f64, [slice, n]); low_rank_matrix(n, &mut arr); - let res: IdDecomposition = arr.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap(); + let res: IdDecomposition = arr + .r_mut() + .into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans) + .unwrap(); println!("The skeleton of the matrix is given by"); res.skel.pretty_print(); @@ -72,17 +75,19 @@ fn test_fat_matrix(slice: usize, n: usize, tol: f64){ let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); - } -fn test_skinny_matrix(slice: usize, n: usize, tol: f64){ +fn test_skinny_matrix(slice: usize, n: usize, tol: f64) { let mut arr: DynamicArray = rlst_dynamic_array2!(f64, [slice, n]); low_rank_matrix(n, &mut arr); let mut arr_trans = empty_array(); arr_trans.fill_from_resize(arr.r().transpose()); - let res: IdDecomposition = arr_trans.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::Trans).unwrap(); + let res: IdDecomposition = arr_trans + .r_mut() + .into_id_alloc(Accuracy::Tol(tol), TransMode::Trans) + .unwrap(); println!("The skeleton of the matrix is given by"); res.skel.pretty_print(); @@ -106,7 +111,6 @@ fn test_skinny_matrix(slice: usize, n: usize, tol: f64){ let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2(); println!("Interpolative Decomposition L2 absolute error: {}", error); - } pub fn main() { @@ -118,5 +122,4 @@ pub fn main() { //Create a low rank matrix test_fat_matrix(slice, n, tol); test_skinny_matrix(slice, n, tol); - } From a7daebca70e6e44d5545d2ff296d9f53229f73e0 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Mon, 28 Apr 2025 13:00:34 +0100 Subject: [PATCH 19/29] remove conj --- src/dense/linalg/interpolative_decomposition.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 8710bfee..5dbbcf24 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -304,10 +304,10 @@ macro_rules! impl_id { let mut u_tri = empty_array(); match trans_mode { - TransMode::Trans => arr_work.fill_from_resize(arr.r().conj()), - TransMode::NoTrans => arr_work.fill_from_resize(arr.r().conj().transpose()), - TransMode::ConjNoTrans => arr_work.fill_from_resize(arr.r()), - TransMode::ConjTrans => arr_work.fill_from_resize(arr.r().transpose()), + TransMode::Trans => arr_work.fill_from_resize(arr.r()), + TransMode::NoTrans => arr_work.fill_from_resize(arr.r().transpose()), + TransMode::ConjNoTrans => arr_work.fill_from_resize(arr.r().conj()), + TransMode::ConjTrans => arr_work.fill_from_resize(arr.r().transpose().conj()), }; let shape = arr_work.shape(); From 1909047a96955ac5825283d2ec262965f7efe290 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Tue, 29 Apr 2025 13:26:07 +0100 Subject: [PATCH 20/29] relaxed tolerance implemented --- .../linalg/interpolative_decomposition.rs | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 5dbbcf24..c899c142 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -114,6 +114,17 @@ pub trait MatrixIdDecomposition: Sized { tol: ::Real, ) -> usize; + /// Relax tolerance when we have a plateau of eigenvalues + fn rank_from_tolerance_relaxed< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Self::Item> + + UnsafeRandomAccessByRef<2, Item = Self::Item>, + >( + ut_mat: Array, + tol: ::Real, + ) -> usize; + ///Compute the permutation matrix associated to the Interpolative Decomposition fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> @@ -146,6 +157,8 @@ pub enum Accuracy { FixedRank(usize), /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range MaxRank(T, usize), + /// Uses the relaxed tolerance method + RelaxedTol(T), } ///Interface to obtain the upper-triangular part of the matrix @@ -331,6 +344,9 @@ macro_rules! impl_id { Accuracy::MaxRank(tol, k) => { rank = std::cmp::max(k, Self::rank_from_tolerance(u_tri.r_mut(), tol)); } + Accuracy::RelaxedTol(tol) => { + rank = Self::rank_from_tolerance_relaxed(u_tri.r_mut(), tol); + } } let mut permutation = rlst_dynamic_array2!($scalar, [shape[1], shape[1]]); @@ -415,6 +431,42 @@ macro_rules! impl_id { } } + fn rank_from_tolerance_relaxed< + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = $scalar> + + UnsafeRandomAccessByRef<2, Item = $scalar>, + >( + ut_mat: Array<$scalar, ArrayImplMut, 2>, + tol: <$scalar as RlstScalar>::Real, + ) -> usize { + let dim = ut_mat.shape()[0]; + let max = ut_mat.get([0, 0]).unwrap().abs(); + + //We compute the rank of the matrix + let mut rank = 0; + if max.re() > 0.0 { + for i in 0..dim { + if ut_mat.get([i, i]).unwrap().abs() > tol * max { + rank += 1; + } + } + } else { + rank = dim + } + + for i in (0..rank - 1).rev() { + let diff = (ut_mat.get([i, i]).unwrap().abs() + - ut_mat.get([rank, rank]).unwrap().abs()) + .abs(); + if diff > 1e-1 { + return i + 1; + } + } + + rank + } + fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + Shape<2> From e075210f0f420af312ce085aef87bf43547bf211 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Tue, 29 Apr 2025 13:45:58 +0100 Subject: [PATCH 21/29] relaxation implemented --- src/dense/linalg/interpolative_decomposition.rs | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index c899c142..f5453bb3 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -455,12 +455,11 @@ macro_rules! impl_id { rank = dim } - for i in (0..rank - 1).rev() { - let diff = (ut_mat.get([i, i]).unwrap().abs() - - ut_mat.get([rank, rank]).unwrap().abs()) - .abs(); - if diff > 1e-1 { - return i + 1; + for i in (0..rank + 1).rev() { + if ut_mat.get([i, i]).unwrap().abs() + > 1e-1 * ut_mat.get([i - 1, i - 1]).unwrap().abs() + { + return i - 1; } } From 7bc60cc225ecbe304c687bb5d8a4c631e49e8199 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Tue, 29 Apr 2025 17:00:56 +0100 Subject: [PATCH 22/29] clone id accuracy --- src/dense/linalg/interpolative_decomposition.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index f5453bb3..0355a9cb 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -149,6 +149,7 @@ pub struct IdDecomposition { pub id_mat: Array, 2>, 2>, } +#[derive(Debug, Clone)] ///Options to decide the matrix rank pub enum Accuracy { /// Indicates that the rank of the decomposition will be computed from a given tolerance From 80fa3c3955b8f90ca2ca06b186fd7853bb1cb97c Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Wed, 30 Apr 2025 12:15:49 +0100 Subject: [PATCH 23/29] triangular matrices implemented --- src/dense/linalg.rs | 1 + .../linalg/interpolative_decomposition.rs | 15 +- src/dense/linalg/triangular_arrays.rs | 243 ++++++++++++++++++ src/dense/types.rs | 18 ++ src/prelude.rs | 6 +- 5 files changed, 275 insertions(+), 8 deletions(-) create mode 100644 src/dense/linalg/triangular_arrays.rs diff --git a/src/dense/linalg.rs b/src/dense/linalg.rs index 7e2b41f5..0dc24457 100644 --- a/src/dense/linalg.rs +++ b/src/dense/linalg.rs @@ -13,6 +13,7 @@ pub mod null_space; pub mod pseudo_inverse; pub mod qr; pub mod svd; +pub mod triangular_arrays; /// Return true if stride is column major as required by Lapack. pub fn assert_lapack_stride(stride: [usize; 2]) { diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 0355a9cb..e6a796c1 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -6,11 +6,11 @@ use crate::dense::traits::{ MultIntoResize, RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByRef, UnsafeRandomAccessByValue, UnsafeRandomAccessMut, }; -use crate::dense::types::TransMode; // Import TransMode from the appropriate module use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; +use crate::dense::types::{Side, TransMode, TriangularType}; // Import TransMode from the appropriate module +use crate::DynamicArray; use crate::{empty_array, rlst_dynamic_array2, BaseArray, VectorContainer}; -use crate::{DynamicArray, RawAccess}; -use blas::{ctrsm, dtrsm, strsm, ztrsm}; +use crate::{TriangularMatrix, TriangularOperations}; // Import TriangularType from the appropriate module /// Compute the matrix interpolative decomposition, by providing a rank and an interpolation matrix. /// /// The matrix interpolative decomposition is defined for a two dimensional 'long' array `arr` of @@ -161,7 +161,7 @@ pub enum Accuracy { /// Uses the relaxed tolerance method RelaxedTol(T), } - +/* ///Interface to obtain the upper-triangular part of the matrix pub trait UpperTriangular: RlstScalar { ///Compute the upper triangular @@ -295,7 +295,7 @@ macro_rules! implement_solve_upper_triangular { implement_solve_upper_triangular!(f32, strsm); implement_solve_upper_triangular!(f64, dtrsm); implement_solve_upper_triangular!(c32, ctrsm); -implement_solve_upper_triangular!(c64, ztrsm); +implement_solve_upper_triangular!(c64, ztrsm);*/ macro_rules! impl_id { ($scalar:ty) => { @@ -386,15 +386,16 @@ macro_rules! impl_id { let shape: [usize; 2] = [shape[1], shape[1]]; let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim - rank, rank]); - let r11 = UpperTriangularMatrix::<$scalar>::new( + let r11 = TriangularMatrix::<$scalar>::new( u_tri.r_mut().into_subview([0, 0], [rank, rank]), + TriangularType::Upper, ) .unwrap(); let mut r12 = u_tri .r_mut() .into_subview([0, rank], [rank, shape[1] - rank]); - r11.solve_upper_triangular(&mut r12); + r11.solve(&mut r12, Side::Left, TransMode::NoTrans); id_mat.fill_from(r12.r().conj().transpose().r()); Ok(Self { diff --git a/src/dense/linalg/triangular_arrays.rs b/src/dense/linalg/triangular_arrays.rs new file mode 100644 index 00000000..21fbc970 --- /dev/null +++ b/src/dense/linalg/triangular_arrays.rs @@ -0,0 +1,243 @@ +//! Interpolative decomposition of a matrix. +use crate::dense::array::Array; +use crate::dense::traits::{ + RandomAccessByRef, RawAccessMut, Shape, Stride, UnsafeRandomAccessByValue, +}; +use crate::dense::types::{c32, c64, RlstResult, RlstScalar}; +use crate::{rlst_dynamic_array2, Side, TransMode, TriangularType}; +use crate::{DynamicArray, RawAccess}; +use blas::{ctrmm, ctrsm, dtrmm, dtrsm, strmm, strsm, ztrmm, ztrsm}; + +///Interface to obtain the upper-triangular part of the matrix +pub trait Triangular: RlstScalar { + ///Compute the upper triangular + fn into_triangular_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self>, + >( + arr: Array, + triangular_type: TriangularType, + ) -> RlstResult>; +} + +macro_rules! implement_into_triangular { + ($scalar:ty) => { + impl Triangular for $scalar { + fn into_triangular_alloc< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self>, + >( + arr: Array, + triangular_type: TriangularType, + ) -> RlstResult> { + TriangularMatrix::<$scalar>::new(arr, triangular_type) + } + } + }; +} + +implement_into_triangular!(f32); +implement_into_triangular!(f64); +implement_into_triangular!(c32); +implement_into_triangular!(c64); + +/// A struct representing an upper triangular matrix. +pub struct TriangularMatrix { + /// The upper triangular part of the matrix. + pub tri: DynamicArray, + /// Defines if the matrix is lower or upper triangular + pub triangular_type: TriangularType, +} + +///Solver for upper-triangular matrix. +pub trait TriangularOperations: Sized { + ///Defines an abstract type for the matrix to solve + type Item: RlstScalar; + ///Extract the upper triangular and save it + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self::Item>, + >( + arr: Array, + triangular_type: TriangularType, + ) -> RlstResult; + ///Solves the upper-triangular system + fn solve< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array, + side: Side, + trans: TransMode, + ); + + /// Multiplies the triangular matrix by a regular matrix + fn mul< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array, + side: Side, + trans: TransMode, + ); +} + +macro_rules! implement_solve_upper_triangular { + ($scalar:ty, $trsm:expr, $trmm:expr) => { + impl TriangularOperations for TriangularMatrix<$scalar> { + type Item = $scalar; + + fn new< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut + + RandomAccessByRef<2, Item = Self::Item>, + >( + arr: Array, + triangular_type: TriangularType, + ) -> RlstResult { + let shape = arr.shape(); + let mut tri = rlst_dynamic_array2!(Self::Item, shape); + let mut view = tri.r_mut(); + for i in 0..shape[0] { + for j in i..shape[1] { + view[[i, j]] = arr[[i, j]]; + } + } + Ok(Self { + tri, + triangular_type, + }) + } + + fn solve< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array<$scalar, ArrayImpl, 2>, + side: Side, + trans: TransMode, + ) { + let lda = self.tri.stride()[1] as i32; + let m = arr_b.shape()[0] as i32; + let n = arr_b.shape()[1] as i32; + let ldb = arr_b.stride()[1] as i32; + let alpha = 1.0; + + let t_type_char = match self.triangular_type { + TriangularType::Upper => b'U', // A is on the left + TriangularType::Lower => b'L', // A is on the right + }; + + let side_char = match side { + Side::Left => b'L', // A is on the left + Side::Right => b'R', // A is on the right + }; + + let trans_char = match trans { + TransMode::NoTrans => b'N', + TransMode::ConjNoTrans => { + panic!("TransMode::ConjNoTrans not supported for trsm implementation.") + } + TransMode::Trans => b'T', + TransMode::ConjTrans => b'C', + }; + + unsafe { + // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb) + $trsm( + side_char, // Side: A is on the left + t_type_char, // Uplo: A is upper triangular + trans_char, // TransA: no transpose + b'N', // Diag: A is not unit triangular + m, + n, + alpha.into(), + self.tri.data(), + lda, + arr_b.data_mut(), + ldb, + ); + } + } + + fn mul< + ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + + Stride<2> + + Shape<2> + + RawAccessMut, + >( + &self, + arr_b: &mut Array<$scalar, ArrayImpl, 2>, + side: Side, + trans: TransMode, + ) { + let lda = self.tri.stride()[1] as i32; + let m = arr_b.shape()[0] as i32; + let n = arr_b.shape()[1] as i32; + let ldb = arr_b.stride()[1] as i32; + let alpha = 1.0; + + let t_type_char = match self.triangular_type { + TriangularType::Upper => b'U', // A is on the left + TriangularType::Lower => b'L', // A is on the right + }; + + let side_char = match side { + Side::Left => b'L', // A is on the left + Side::Right => b'R', // A is on the right + }; + + let trans_char = match trans { + TransMode::NoTrans => b'N', + TransMode::ConjNoTrans => { + panic!("TransMode::ConjNoTrans not supported for trsm implementation.") + } + TransMode::Trans => b'T', + TransMode::ConjTrans => b'C', + }; + + unsafe { + // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb) + $trmm( + side_char, // Side: A is on the left + t_type_char, // Uplo: A is upper triangular + trans_char, // TransA: no transpose + b'N', // Diag: A is not unit triangular + m, + n, + alpha.into(), + self.tri.data(), + lda, + arr_b.data_mut(), + ldb, + ); + } + } + } + }; +} + +implement_solve_upper_triangular!(f32, strsm, strmm); +implement_solve_upper_triangular!(f64, dtrsm, dtrmm); +implement_solve_upper_triangular!(c32, ctrsm, ctrmm); +implement_solve_upper_triangular!(c64, ztrsm, ztrmm); diff --git a/src/dense/types.rs b/src/dense/types.rs index 815ff0b4..9cafa02d 100644 --- a/src/dense/types.rs +++ b/src/dense/types.rs @@ -431,3 +431,21 @@ pub enum TransMode { /// Conjugate transpose of matrix. ConjTrans, } + +/// Defines if the matrix is lower or upper triangular +#[derive(Clone, Copy, PartialEq)] +pub enum TriangularType { + /// The matrix is upper-triangular + Upper, + /// The matrix is lower-triangular + Lower, +} + +/// Defines the side through which the upper or lower triangular matrix acts +#[derive(Clone, Copy, PartialEq)] +pub enum Side { + /// Left side + Left, + /// Right side + Right, +} diff --git a/src/prelude.rs b/src/prelude.rs index 7067534d..9238c7dd 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -51,7 +51,8 @@ pub use crate::dense::traits::{ }; pub use crate::dense::types::{ - c32, c64, DataChunk, RlstBase, RlstError, RlstNum, RlstResult, RlstScalar, TransMode, + c32, c64, DataChunk, RlstBase, RlstError, RlstNum, RlstResult, RlstScalar, Side, TransMode, + TriangularType, }; pub use crate::dense::base_array::BaseArray; @@ -68,6 +69,9 @@ pub use crate::dense::linalg::null_space::{MatrixNull, NullSpace}; pub use crate::dense::linalg::pseudo_inverse::MatrixPseudoInverse; pub use crate::dense::linalg::qr::{MatrixQr, MatrixQrDecomposition, QrDecomposition}; pub use crate::dense::linalg::svd::{MatrixSvd, SvdMode}; +pub use crate::dense::linalg::triangular_arrays::{ + Triangular, TriangularMatrix, TriangularOperations, +}; pub use crate::dense::array::rank1_array::Rank1Array; From 848be4b08e62fa5ac8b41a48d826bf7a691c9324 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Wed, 30 Apr 2025 17:09:24 +0100 Subject: [PATCH 24/29] remove unnecessary constraints --- .../linalg/interpolative_decomposition.rs | 137 +----------------- src/dense/linalg/triangular_arrays.rs | 10 +- 2 files changed, 5 insertions(+), 142 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index e6a796c1..439efa3f 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -161,141 +161,6 @@ pub enum Accuracy { /// Uses the relaxed tolerance method RelaxedTol(T), } -/* -///Interface to obtain the upper-triangular part of the matrix -pub trait UpperTriangular: RlstScalar { - ///Compute the upper triangular - fn into_upper_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut - + RandomAccessByRef<2, Item = Self>, - >( - arr: Array, - ) -> RlstResult>; -} - -macro_rules! implement_into_upper_triangular { - ($scalar:ty) => { - impl UpperTriangular for $scalar { - fn into_upper_alloc< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> - + Stride<2> - + Shape<2> - + RawAccessMut - + RandomAccessByRef<2, Item = Self>, - >( - arr: Array, - ) -> RlstResult> { - UpperTriangularMatrix::<$scalar>::new(arr) - } - } - }; -} - -implement_into_upper_triangular!(f32); -implement_into_upper_triangular!(f64); -implement_into_upper_triangular!(c32); -implement_into_upper_triangular!(c64); - -/// A struct representing an upper triangular matrix. -pub struct UpperTriangularMatrix { - /// The upper triangular part of the matrix. - pub upper: DynamicArray, -} - -///Solver for upper-triangular matrix. -pub trait UppperTriangularOperations: Sized { - ///Defines an abstract type for the matrix to solve - type Item: RlstScalar; - ///Extract the upper triangular and save it - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + Shape<2> - + RawAccessMut - + RandomAccessByRef<2, Item = Self::Item>, - >( - arr: Array, - ) -> RlstResult; - ///Solves the upper-triangular system - fn solve_upper_triangular< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - &self, - arr_b: &mut Array, - ); -} - -macro_rules! implement_solve_upper_triangular { - ($scalar:ty, $trsm:expr) => { - impl UppperTriangularOperations for UpperTriangularMatrix<$scalar> { - type Item = $scalar; - - fn new< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + Shape<2> - + RawAccessMut - + RandomAccessByRef<2, Item = Self::Item>, - >( - arr: Array, - ) -> RlstResult { - let shape = arr.shape(); - let mut upper = rlst_dynamic_array2!(Self::Item, shape); - let mut view = upper.r_mut(); - for i in 0..shape[0] { - for j in i..shape[1] { - view[[i, j]] = arr[[i, j]]; - } - } - Ok(Self { upper }) - } - - fn solve_upper_triangular< - ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Stride<2> - + Shape<2> - + RawAccessMut, - >( - &self, - arr_b: &mut Array<$scalar, ArrayImpl, 2>, - ) { - let lda = self.upper.stride()[1] as i32; - let m = arr_b.shape()[0] as i32; - let n = arr_b.shape()[1] as i32; - let ldb = arr_b.stride()[1] as i32; - let alpha = 1.0; - - unsafe { - // dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb) - $trsm( - b'L', // Side: A is on the left - b'U', // Uplo: A is upper triangular - b'N', // TransA: no transpose - b'N', // Diag: A is not unit triangular - m, - n, - alpha.into(), - self.upper.data(), - lda, - arr_b.data_mut(), - ldb, - ); - } - } - } - }; -} - -implement_solve_upper_triangular!(f32, strsm); -implement_solve_upper_triangular!(f64, dtrsm); -implement_solve_upper_triangular!(c32, ctrsm); -implement_solve_upper_triangular!(c64, ztrsm);*/ macro_rules! impl_id { ($scalar:ty) => { @@ -387,7 +252,7 @@ macro_rules! impl_id { let mut id_mat: DynamicArray<$scalar, 2> = rlst_dynamic_array2!($scalar, [dim - rank, rank]); let r11 = TriangularMatrix::<$scalar>::new( - u_tri.r_mut().into_subview([0, 0], [rank, rank]), + &u_tri.r_mut().into_subview([0, 0], [rank, rank]), TriangularType::Upper, ) .unwrap(); diff --git a/src/dense/linalg/triangular_arrays.rs b/src/dense/linalg/triangular_arrays.rs index 21fbc970..9a0dbe7b 100644 --- a/src/dense/linalg/triangular_arrays.rs +++ b/src/dense/linalg/triangular_arrays.rs @@ -18,7 +18,7 @@ pub trait Triangular: RlstScalar { + RawAccessMut + RandomAccessByRef<2, Item = Self>, >( - arr: Array, + arr: &Array, triangular_type: TriangularType, ) -> RlstResult>; } @@ -33,7 +33,7 @@ macro_rules! implement_into_triangular { + RawAccessMut + RandomAccessByRef<2, Item = Self>, >( - arr: Array, + arr: &Array, triangular_type: TriangularType, ) -> RlstResult> { TriangularMatrix::<$scalar>::new(arr, triangular_type) @@ -64,10 +64,9 @@ pub trait TriangularOperations: Sized { ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + Stride<2> + Shape<2> - + RawAccessMut + RandomAccessByRef<2, Item = Self::Item>, >( - arr: Array, + arr: &Array, triangular_type: TriangularType, ) -> RlstResult; ///Solves the upper-triangular system @@ -106,10 +105,9 @@ macro_rules! implement_solve_upper_triangular { ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self::Item> + Stride<2> + Shape<2> - + RawAccessMut + RandomAccessByRef<2, Item = Self::Item>, >( - arr: Array, + arr: &Array, triangular_type: TriangularType, ) -> RlstResult { let shape = arr.shape(); From 9bddc53bf1e9d979b7c588861bcd4202d4f0117b Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Thu, 1 May 2025 14:48:37 +0100 Subject: [PATCH 25/29] fix lower triangular --- src/dense/linalg/triangular_arrays.rs | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/dense/linalg/triangular_arrays.rs b/src/dense/linalg/triangular_arrays.rs index 9a0dbe7b..766360e3 100644 --- a/src/dense/linalg/triangular_arrays.rs +++ b/src/dense/linalg/triangular_arrays.rs @@ -113,9 +113,21 @@ macro_rules! implement_solve_upper_triangular { let shape = arr.shape(); let mut tri = rlst_dynamic_array2!(Self::Item, shape); let mut view = tri.r_mut(); - for i in 0..shape[0] { - for j in i..shape[1] { - view[[i, j]] = arr[[i, j]]; + + match triangular_type { + TriangularType::Upper => { + for i in 0..shape[0] { + for j in i..shape[1] { + view[[i, j]] = arr[[i, j]]; + } + } + } + TriangularType::Lower => { + for i in 0..shape[0] { + for j in 0..=i { + view[[i, j]] = arr[[i, j]]; + } + } } } Ok(Self { From 398c87793befed8aac19aaa7b4c91f179204d406 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 9 May 2025 14:12:47 +0100 Subject: [PATCH 26/29] format changes --- .../linalg/interpolative_decomposition.rs | 120 +++++------------- src/dense/linalg/null_space.rs | 8 +- 2 files changed, 34 insertions(+), 94 deletions(-) diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 439efa3f..5256bd39 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -103,28 +103,6 @@ pub trait MatrixIdDecomposition: Sized { trans_mode: TransMode, ) -> RlstResult; - /// Compute the rank of the decomposition from a given tolerance - fn rank_from_tolerance< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = Self::Item> - + UnsafeRandomAccessByRef<2, Item = Self::Item>, - >( - ut_mat: Array, - tol: ::Real, - ) -> usize; - - /// Relax tolerance when we have a plateau of eigenvalues - fn rank_from_tolerance_relaxed< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = Self::Item> - + UnsafeRandomAccessByRef<2, Item = Self::Item>, - >( - ut_mat: Array, - tol: ::Real, - ) -> usize; - ///Compute the permutation matrix associated to the Interpolative Decomposition fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Self::Item> @@ -158,8 +136,6 @@ pub enum Accuracy { FixedRank(usize), /// Computes the rank from the tolerance, and if this one is smaller than a user set range, then we stick to the user set range MaxRank(T, usize), - /// Uses the relaxed tolerance method - RelaxedTol(T), } macro_rules! impl_id { @@ -204,14 +180,11 @@ macro_rules! impl_id { //The rank can be given a priori, in which case, we do not need to compute the rank using the tolerance parameter. match rank_param { Accuracy::Tol(tol) => { - rank = Self::rank_from_tolerance(u_tri.r_mut(), tol); + rank = rank_from_tolerance(u_tri.r_mut(), tol); } Accuracy::FixedRank(k) => rank = k, Accuracy::MaxRank(tol, k) => { - rank = std::cmp::max(k, Self::rank_from_tolerance(u_tri.r_mut(), tol)); - } - Accuracy::RelaxedTol(tol) => { - rank = Self::rank_from_tolerance_relaxed(u_tri.r_mut(), tol); + rank = std::cmp::max(k, rank_from_tolerance(u_tri.r_mut(), tol)); } } @@ -272,67 +245,6 @@ macro_rules! impl_id { } } - fn rank_from_tolerance< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = $scalar> - + UnsafeRandomAccessByRef<2, Item = $scalar>, - >( - ut_mat: Array<$scalar, ArrayImplMut, 2>, - tol: <$scalar as RlstScalar>::Real, - ) -> usize { - let dim = ut_mat.shape()[0]; - let max = ut_mat.get([0, 0]).unwrap().abs(); - - //We compute the rank of the matrix - if max.re() > 0.0 { - let mut rank = 0; - for i in 0..dim { - if ut_mat.get([i, i]).unwrap().abs() > tol * max { - rank += 1; - } - } - rank - } else { - dim - } - } - - fn rank_from_tolerance_relaxed< - ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> - + Shape<2> - + UnsafeRandomAccessMut<2, Item = $scalar> - + UnsafeRandomAccessByRef<2, Item = $scalar>, - >( - ut_mat: Array<$scalar, ArrayImplMut, 2>, - tol: <$scalar as RlstScalar>::Real, - ) -> usize { - let dim = ut_mat.shape()[0]; - let max = ut_mat.get([0, 0]).unwrap().abs(); - - //We compute the rank of the matrix - let mut rank = 0; - if max.re() > 0.0 { - for i in 0..dim { - if ut_mat.get([i, i]).unwrap().abs() > tol * max { - rank += 1; - } - } - } else { - rank = dim - } - - for i in (0..rank + 1).rev() { - if ut_mat.get([i, i]).unwrap().abs() - > 1e-1 * ut_mat.get([i - 1, i - 1]).unwrap().abs() - { - return i - 1; - } - } - - rank - } - fn get_p< ArrayImplMut: UnsafeRandomAccessByValue<2, Item = $scalar> + Shape<2> @@ -357,3 +269,31 @@ impl_id!(f64); impl_id!(f32); impl_id!(c32); impl_id!(c64); + +/// Compute the rank of the decomposition from a given tolerance +fn rank_from_tolerance< + Item: RlstScalar, + ArrayImplMut: UnsafeRandomAccessByValue<2, Item = Item> + + Shape<2> + + UnsafeRandomAccessMut<2, Item = Item> + + UnsafeRandomAccessByRef<2, Item = Item>, +>( + ut_mat: Array, + tol: ::Real, +) -> usize { + let dim = ut_mat.shape()[0]; + let max = ut_mat.get([0, 0]).unwrap().abs(); + + //We compute the rank of the matrix + if max.re() > num::Zero::zero() { + let mut rank = 0; + for i in 0..dim { + if ut_mat.get([i, i]).unwrap().abs() > tol * max { + rank += 1; + } + } + rank + } else { + dim + } +} diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index 8b6734eb..fb7dbc8b 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -168,7 +168,7 @@ fn svd_nullification< //For a full rank rectangular matrix, then rank = dim. //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = find_matrix_rank_svd::(&mut singular_values, dim, tol); + let rank: usize = find_svd_rank::(&mut singular_values, dim, tol); //The null space is given by the last shape[1]-rank columns of V let mut null_space_arr: DynamicArray = empty_array(); @@ -201,13 +201,13 @@ where qr.get_r(r_mat.r_mut()); //For a full rank rectangular matrix, then rank = dim. //find_matrix_rank checks if the matrix is full rank and recomputes the rank. - let rank: usize = find_matrix_rank_qr::(&r_mat, dim, tol); + let rank: usize = find_qr_rank::(&r_mat, dim, tol); let mut null_space_arr = empty_array(); null_space_arr.fill_from_resize(q.into_subview([0, shape[1]], [shape[0], shape[0] - rank])); null_space_arr } -fn find_matrix_rank_svd( +fn find_svd_rank( singular_values: &mut DynamicArray, 1>, dim: usize, tol: ::Real, @@ -233,7 +233,7 @@ fn find_matrix_rank_svd( rank } -fn find_matrix_rank_qr( +fn find_qr_rank( r_mat: &DynamicArray, dim: usize, tol: ::Real, From c5ab27ff93a080ad4c83f560ea610a78165c4154 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 9 May 2025 14:37:35 +0100 Subject: [PATCH 27/29] docs fixed --- examples/rlst_null_space.rs | 4 ++-- .../linalg/interpolative_decomposition.rs | 4 +++- src/dense/linalg/null_space.rs | 9 +++---- src/doc/dense_linear_algebra.rs | 24 ++++++------------- 4 files changed, 17 insertions(+), 24 deletions(-) diff --git a/examples/rlst_null_space.rs b/examples/rlst_null_space.rs index 640b36e0..91e2eb9e 100644 --- a/examples/rlst_null_space.rs +++ b/examples/rlst_null_space.rs @@ -16,11 +16,11 @@ pub fn main() { res.view_flat().norm_2() ); - let mut arr = rlst_dynamic_array2!(f64, [3, 4]); + let mut arr = rlst_dynamic_array2!(f64, [4, 3]); arr.fill_from_seed_equally_distributed(0); let shape = arr.shape(); let mut arr_qr = rlst_dynamic_array2!(f64, [shape[1], shape[0]]); - arr_qr.fill_from(arr.r().conj().transpose()); + arr_qr.fill_from(arr.r().transpose()); let null_res = arr_qr.r_mut().into_null_alloc(tol, Method::Qr).unwrap(); let res: Array, 2>, 2> = empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r()); diff --git a/src/dense/linalg/interpolative_decomposition.rs b/src/dense/linalg/interpolative_decomposition.rs index 5256bd39..93693ff2 100644 --- a/src/dense/linalg/interpolative_decomposition.rs +++ b/src/dense/linalg/interpolative_decomposition.rs @@ -21,10 +21,12 @@ use crate::{TriangularMatrix, TriangularOperations}; // Import TriangularType fr /// The following command computes the interpolative decomposition of an array `a` for a given tolerance, tol. /// ``` /// # use rlst::rlst_dynamic_array2; +/// # use rlst::dense::types::TransMode; +/// # use rlst::dense::linalg::interpolative_decomposition::Accuracy; /// # let tol: f64 = 1e-5; /// # let mut a = rlst_dynamic_array2!(f64, [50, 100]); /// # a.fill_from_seed_equally_distributed(0); -/// # let res = a.r_mut().into_id_alloc(tol, None).unwrap(); +/// # let res = a.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap(); /// ``` /// This method allocates memory for the interpolative decomposition. pub trait MatrixId: RlstScalar { diff --git a/src/dense/linalg/null_space.rs b/src/dense/linalg/null_space.rs index fb7dbc8b..27263c95 100644 --- a/src/dense/linalg/null_space.rs +++ b/src/dense/linalg/null_space.rs @@ -21,11 +21,12 @@ use num::Zero; /// The following command computes the nullspace of an array `a`. /// The nullspace is found in /// ``` -/// # use rlst::rlst_dynamic_array2; -/// # use rlst::dense::linalg::null_space::{NullSpaceType, MatrixNull}; -/// # let mut a = rlst_dynamic_array2!(f64, [3, 4]); +/// # use rlst::prelude::*; +/// # use rlst::dense::linalg::null_space::{Method}; +/// # let mut a = rlst_dynamic_array2!(f64, [4, 3]); +/// # let tol = 1e-10; /// # a.fill_from_seed_equally_distributed(0); -/// # let null_res = a.r_mut().into_null_alloc(NullSpaceType::Row).unwrap(); +/// # let null_res = a.r_mut().into_null_alloc(tol, Method::Qr).unwrap(); /// ``` /// This method allocates memory for the nullspace computation. pub trait MatrixNull: RlstScalar { diff --git a/src/doc/dense_linear_algebra.rs b/src/doc/dense_linear_algebra.rs index 39ffc587..c92bb685 100644 --- a/src/doc/dense_linear_algebra.rs +++ b/src/doc/dense_linear_algebra.rs @@ -487,7 +487,7 @@ //! qr.get_q_alloc(q_mat.r_mut()); //! qr.get_p(p_mat.r_mut()); //! ```` -//! The content of `arr` is overwritten with the QR decomposition. The method [get_q_alloc](crate::dense::linalg::qr::QrDecomposition::get_q_alloc) +//! The content of `arr` is overwritten with the QR decomposition. The method [get_q_alloc](crate::dense::linalg::qr::MatrixQrDecomposition::get_q_alloc) //! needs to allocate additional temporary memory on the heap. This is why it is annoted with `_alloc`. //! //! # Singular value decomposition @@ -520,33 +520,23 @@ //! //! ``` //! # use rlst::prelude::*; +//! # use rlst::dense::linalg::interpolative_decomposition::Accuracy; //! let mut rand = rand::thread_rng(); //! let mut arr = rlst_dynamic_array2!(f64, [5, 8]); //! let tol: f64 = 1e-5; -//! let mut res = arr.view_mut().into_id_alloc(tol, None).unwrap(); -//! ``` -//! -//! You can take a look at the interpolation matrix using -//! -//! ``` -//! res.id_mat.pretty_print(); -//! ``` -//! -//! With this, the approximation of the low-rank elements of `arr` is given by -//! -//! ``` -//! let arr_app = empty_array().simple_mult_into_resize(res.id_mat.view(), res.arr.view_mut().into_subview([0, 0], [res.rank, n])); +//! let res = arr.r_mut().into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans).unwrap(); //! ``` //! //! You can also obtain the Interpolative Decomposition by giving the algorithm a fixed rank, and the tolerance parameter is ignored, like so //! //! ``` //! # use rlst::prelude::*; +//! # use rlst::dense::linalg::interpolative_decomposition::Accuracy; +//! //! let mut rand = rand::thread_rng(); //! let mut arr = rlst_dynamic_array2!(f64, [5, 8]); -//! let tol: f64 = 1e-5; -//! k = 2; -//! let mut res = arr.view_mut().into_id_alloc(tol, k).unwrap(); +//! let k = 2; +//! let res = arr.r_mut().into_id_alloc(Accuracy::FixedRank(k), TransMode::NoTrans).unwrap(); //! ``` //! //! # Other vector functions From 891d4a026dde6666b325c02969eb9865a6a3b134 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 9 May 2025 14:45:36 +0100 Subject: [PATCH 28/29] style check fixed --- src/doc/dense_linear_algebra.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/doc/dense_linear_algebra.rs b/src/doc/dense_linear_algebra.rs index c92bb685..366b4684 100644 --- a/src/doc/dense_linear_algebra.rs +++ b/src/doc/dense_linear_algebra.rs @@ -532,7 +532,6 @@ //! ``` //! # use rlst::prelude::*; //! # use rlst::dense::linalg::interpolative_decomposition::Accuracy; -//! //! let mut rand = rand::thread_rng(); //! let mut arr = rlst_dynamic_array2!(f64, [5, 8]); //! let k = 2; From 29ed540379572f4d0f89d3725eda576235116082 Mon Sep 17 00:00:00 2001 From: Ignacia Piccardo Date: Fri, 9 May 2025 14:54:58 +0100 Subject: [PATCH 29/29] over indentation fixed --- src/dense/linalg/pseudo_inverse.rs | 8 ++++---- src/dense/linalg/svd.rs | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/dense/linalg/pseudo_inverse.rs b/src/dense/linalg/pseudo_inverse.rs index eec186b6..5b0cc40e 100644 --- a/src/dense/linalg/pseudo_inverse.rs +++ b/src/dense/linalg/pseudo_inverse.rs @@ -45,9 +45,9 @@ pub trait MatrixPseudoInverse: RlstScalar + MatrixSvd { /// /// # Parameters /// - `pinv`: Array to store the pseudo-inverse in. If `self` has shape `[m, n]` then - /// `pinv` must have shape `[n, m]`. + /// `pinv` must have shape `[n, m]`. /// - `tol`: The relative tolerance. Singular values smaller or equal `tol * s\[0\]` will be discarded, - /// where s\[0\] is the largest singular value. + /// where s\[0\] is the largest singular value. fn into_pseudo_inverse_alloc< ArrayImpl: UnsafeRandomAccessByValue<2, Item = Self> + Stride<2> @@ -103,9 +103,9 @@ impl< /// /// # Parameters /// - `pinv`: Array to store the pseudo-inverse in. If `self` has shape `[m, n]` then - /// `pinv` must have shape `[n, m]`. + /// `pinv` must have shape `[n, m]`. /// - `tol`: The relative tolerance. Singular values smaller or equal `tol * s\[0\]` will be discarded, - /// where s\[0\] is the largest singular value. + /// where s\[0\] is the largest singular value. pub fn into_pseudo_inverse_alloc< ArrayImplPInv: UnsafeRandomAccessByValue<2, Item = Item> + Stride<2> diff --git a/src/dense/linalg/svd.rs b/src/dense/linalg/svd.rs index 7d4b9ea8..a29ab1a5 100644 --- a/src/dense/linalg/svd.rs +++ b/src/dense/linalg/svd.rs @@ -34,11 +34,11 @@ pub trait MatrixSvd: RlstScalar { /// # Parameters /// /// - `u` - Stores the matrix `U`. For the full SVD the shape - /// needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`. + /// needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`. /// - `vt` - Stores the matrix `Vt`. For the full SVD the shape needs to be `(n, n)`. - /// For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores - /// the complex conjugate transpose of the matrix of right singular vectors. - /// Hence, the columns of `vt.transpose().conj()` will be the right singular vectors. + /// For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores + /// the complex conjugate transpose of the matrix of right singular vectors. + /// Hence, the columns of `vt.transpose().conj()` will be the right singular vectors. /// - `singular_values` - Stores the `k` singular values of `A`. /// - `mode` - Choose between full SVD [SvdMode::Full] or reduced SVD [SvdMode::Reduced]. /// @@ -504,11 +504,11 @@ impl< /// # Parameters /// /// - `u` - Stores the matrix `U`. For the full SVD the shape - /// needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`. + /// needs to be `(m, m)`. For the reduced SVD it needs to be `(m, k)`. /// - `vt` - Stores the matrix `Vt`. For the full SVD the shape needs to be `(n, n)`. - /// For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores - /// the complex conjugate transpose of the matrix of right singular vectors. - /// Hence, the columns of `vt.transpose().conj()` will be the right singular vectors. + /// For the reduced SVD it needs to be `(k, n)`. Note that `vt` stores + /// the complex conjugate transpose of the matrix of right singular vectors. + /// Hence, the columns of `vt.transpose().conj()` will be the right singular vectors. /// - `singular_values` - Stores the `k` singular values of `A`. /// - `mode` - Choose between full SVD [SvdMode::Full] or reduced SVD [SvdMode::Reduced]. ///