tensor4all_tcicore/matrixaca.rs
1//! Adaptive Cross Approximation (ACA) implementation.
2//!
3//! [`MatrixACA`] provides the ACA algorithm for incrementally building a
4//! low-rank approximation of a matrix. It stores the approximation in
5//! factored form `A ~ U * diag(alpha) * V` and supports adding pivots
6//! one at a time.
7
8use crate::error::{MatrixCIError, Result};
9use crate::matrix::{
10 append_col, append_row, from_vec2d, get_col, get_row, mat_mul, ncols, nrows, submatrix,
11 submatrix_argmax, zeros, Matrix,
12};
13use crate::scalar::Scalar;
14use crate::traits::AbstractMatrixCI;
15
16/// Adaptive Cross Approximation representation.
17///
18/// Stores a low-rank approximation `A ~ U * diag(alpha) * V`, where
19/// `alpha[k] = 1/delta[k]`. Implements [`AbstractMatrixCI`] for
20/// element-wise evaluation and submatrix extraction.
21///
22/// # Examples
23///
24/// ```
25/// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
26///
27/// let m = from_vec2d(vec![
28/// vec![1.0_f64, 2.0, 3.0],
29/// vec![4.0, 5.0, 6.0],
30/// vec![7.0, 8.0, 9.0],
31/// ]);
32///
33/// // Build ACA from a starting pivot
34/// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (1, 1)).unwrap();
35/// assert_eq!(aca.rank(), 1);
36///
37/// // Add more pivots
38/// aca.add_pivot(&m, (0, 0)).unwrap();
39/// assert_eq!(aca.rank(), 2);
40///
41/// // Evaluate the approximation
42/// let approx_11 = aca.evaluate(1, 1);
43/// assert!((approx_11 - 5.0).abs() < 1e-10);
44/// ```
45#[derive(Debug, Clone)]
46pub struct MatrixACA<T: Scalar> {
47 /// Row indices (I set)
48 row_indices: Vec<usize>,
49 /// Column indices (J set)
50 col_indices: Vec<usize>,
51 /// U matrix: u_k(x) for all x, shape (nrows, npivots)
52 u: Matrix<T>,
53 /// V matrix: v_k(y) for all y, shape (npivots, ncols)
54 v: Matrix<T>,
55 /// Alpha values: 1/delta for each pivot
56 alpha: Vec<T>,
57}
58
59impl<T: Scalar> MatrixACA<T> {
60 /// Create an empty MatrixACA for a matrix of given size
61 ///
62 /// # Examples
63 ///
64 /// ```
65 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA};
66 ///
67 /// let aca = MatrixACA::<f64>::new(3, 4);
68 /// assert_eq!(aca.nrows(), 3);
69 /// assert_eq!(aca.ncols(), 4);
70 /// assert_eq!(aca.rank(), 0);
71 /// assert!(aca.is_empty());
72 /// ```
73 pub fn new(nr: usize, nc: usize) -> Self {
74 Self {
75 row_indices: Vec::new(),
76 col_indices: Vec::new(),
77 u: zeros(nr, 0),
78 v: zeros(0, nc),
79 alpha: Vec::new(),
80 }
81 }
82
83 /// Create a MatrixACA from a matrix with an initial pivot
84 ///
85 /// Returns an error if the pivot value is zero or near-zero.
86 ///
87 /// # Examples
88 ///
89 /// ```
90 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
91 ///
92 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
93 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
94 /// assert_eq!(aca.rank(), 1);
95 /// assert_eq!(aca.row_indices(), &[0]);
96 /// assert_eq!(aca.col_indices(), &[0]);
97 /// ```
98 pub fn from_matrix_with_pivot(a: &Matrix<T>, first_pivot: (usize, usize)) -> Result<Self> {
99 let (i, j) = first_pivot;
100 let pivot_val = a[[i, j]];
101
102 if pivot_val.abs_sq() < f64::EPSILON * f64::EPSILON {
103 return Err(MatrixCIError::SingularMatrix);
104 }
105
106 // u = A[:, j]
107 let u_col = get_col(a, j);
108 let u = from_vec2d((0..nrows(a)).map(|r| vec![u_col[r]]).collect());
109
110 // v = A[i, :]
111 let v_row = get_row(a, i);
112 let v = from_vec2d(vec![v_row]);
113
114 let alpha = vec![T::one() / pivot_val];
115
116 Ok(Self {
117 row_indices: vec![i],
118 col_indices: vec![j],
119 u,
120 v,
121 alpha,
122 })
123 }
124
125 /// Number of pivots
126 ///
127 /// # Examples
128 ///
129 /// ```
130 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
131 ///
132 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
133 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
134 /// assert_eq!(aca.npivots(), 1);
135 /// aca.add_pivot(&m, (1, 1)).unwrap();
136 /// assert_eq!(aca.npivots(), 2);
137 /// ```
138 pub fn npivots(&self) -> usize {
139 self.alpha.len()
140 }
141
142 /// Get reference to U matrix
143 ///
144 /// # Examples
145 ///
146 /// ```
147 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
148 ///
149 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
150 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
151 /// let u = aca.u();
152 /// assert_eq!(u.nrows(), 2); // nrows of original matrix
153 /// assert_eq!(u.ncols(), 1); // one pivot
154 /// ```
155 pub fn u(&self) -> &Matrix<T> {
156 &self.u
157 }
158
159 /// Get reference to V matrix
160 ///
161 /// # Examples
162 ///
163 /// ```
164 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
165 ///
166 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
167 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
168 /// let v = aca.v();
169 /// assert_eq!(v.nrows(), 1); // one pivot
170 /// assert_eq!(v.ncols(), 2); // ncols of original matrix
171 /// ```
172 pub fn v(&self) -> &Matrix<T> {
173 &self.v
174 }
175
176 /// Get reference to alpha values
177 ///
178 /// # Examples
179 ///
180 /// ```
181 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
182 ///
183 /// let m = from_vec2d(vec![vec![2.0_f64, 1.0], vec![1.0, 3.0]]);
184 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
185 /// let alpha = aca.alpha();
186 /// assert_eq!(alpha.len(), 1);
187 /// // alpha = 1 / pivot_value = 1 / m[0,0] = 0.5
188 /// assert!((alpha[0] - 0.5).abs() < 1e-10);
189 /// ```
190 pub fn alpha(&self) -> &[T] {
191 &self.alpha
192 }
193
194 /// Compute u_k(x) for all x (the k-th column of U after adding a new pivot column)
195 fn compute_uk(&self, a: &Matrix<T>) -> Result<Vec<T>> {
196 let k = self.col_indices.len();
197 let yk = self.col_indices[k - 1];
198
199 // result = A[:, yk]
200 let mut result: Vec<T> = get_col(a, yk);
201
202 // Subtract contribution from previous pivots
203 for l in 0..(k - 1) {
204 let xl = self.row_indices[l];
205 let v_l_yk = self.v[[l, yk]];
206 let u_xl_l = self.u[[xl, l]];
207 if u_xl_l.abs_sq() < f64::EPSILON * f64::EPSILON {
208 return Err(MatrixCIError::SingularMatrix);
209 }
210 let factor = v_l_yk / u_xl_l;
211
212 for (i, res) in result.iter_mut().enumerate() {
213 *res = *res - factor * self.u[[i, l]];
214 }
215 }
216
217 Ok(result)
218 }
219
220 /// Compute v_k(y) for all y (the k-th row of V after adding a new pivot row)
221 fn compute_vk(&self, a: &Matrix<T>) -> Result<Vec<T>> {
222 let k = self.row_indices.len();
223 let xk = self.row_indices[k - 1];
224
225 // result = A[xk, :]
226 let mut result: Vec<T> = get_row(a, xk);
227
228 // Subtract contribution from previous pivots
229 for l in 0..(k - 1) {
230 let xl = self.row_indices[l];
231 let u_xk_l = self.u[[xk, l]];
232 let u_xl_l = self.u[[xl, l]];
233 if u_xl_l.abs_sq() < f64::EPSILON * f64::EPSILON {
234 return Err(MatrixCIError::SingularMatrix);
235 }
236 let factor = u_xk_l / u_xl_l;
237
238 for (j, res) in result.iter_mut().enumerate() {
239 *res = *res - factor * self.v[[l, j]];
240 }
241 }
242
243 Ok(result)
244 }
245
246 /// Add a pivot column
247 ///
248 /// # Examples
249 ///
250 /// ```
251 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
252 ///
253 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
254 /// let mut aca = MatrixACA::<f64>::new(2, 2);
255 /// aca.add_pivot_col(&m, 0).unwrap();
256 /// assert_eq!(aca.u().ncols(), 1);
257 /// ```
258 pub fn add_pivot_col(&mut self, a: &Matrix<T>, col_index: usize) -> Result<()> {
259 if col_index >= self.ncols() {
260 return Err(MatrixCIError::IndexOutOfBounds {
261 row: 0,
262 col: col_index,
263 nrows: self.nrows(),
264 ncols: self.ncols(),
265 });
266 }
267
268 self.col_indices.push(col_index);
269 let uk = self.compute_uk(a)?;
270 self.u = append_col(&self.u, &uk);
271
272 Ok(())
273 }
274
275 /// Add a pivot row
276 ///
277 /// # Examples
278 ///
279 /// ```
280 /// use tensor4all_tcicore::{MatrixACA, from_vec2d};
281 ///
282 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
283 /// let mut aca = MatrixACA::<f64>::new(2, 2);
284 /// aca.add_pivot_col(&m, 0).unwrap();
285 /// aca.add_pivot_row(&m, 0).unwrap();
286 /// assert_eq!(aca.npivots(), 1);
287 /// assert_eq!(aca.v().nrows(), 1);
288 /// ```
289 pub fn add_pivot_row(&mut self, a: &Matrix<T>, row_index: usize) -> Result<()> {
290 if row_index >= self.nrows() {
291 return Err(MatrixCIError::IndexOutOfBounds {
292 row: row_index,
293 col: 0,
294 nrows: self.nrows(),
295 ncols: self.ncols(),
296 });
297 }
298
299 self.row_indices.push(row_index);
300 let vk = self.compute_vk(a)?;
301 self.v = append_row(&self.v, &vk);
302
303 // Update alpha
304 let xk = row_index;
305 let u_xk_last = self.u[[xk, ncols(&self.u) - 1]];
306 if u_xk_last.abs_sq() < f64::EPSILON * f64::EPSILON {
307 return Err(MatrixCIError::SingularMatrix);
308 }
309 self.alpha.push(T::one() / u_xk_last);
310
311 Ok(())
312 }
313
314 /// Add a pivot at the given position
315 ///
316 /// # Examples
317 ///
318 /// ```
319 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
320 ///
321 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
322 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
323 /// aca.add_pivot(&m, (1, 1)).unwrap();
324 /// assert_eq!(aca.rank(), 2);
325 /// // Exact reconstruction at all positions
326 /// for i in 0..2 {
327 /// for j in 0..2 {
328 /// assert!((aca.evaluate(i, j) - m[[i, j]]).abs() < 1e-10);
329 /// }
330 /// }
331 /// ```
332 pub fn add_pivot(&mut self, a: &Matrix<T>, pivot: (usize, usize)) -> Result<()> {
333 self.add_pivot_col(a, pivot.1)?;
334 self.add_pivot_row(a, pivot.0)?;
335 Ok(())
336 }
337
338 /// Add a pivot that maximizes the error using ACA heuristic
339 ///
340 /// # Examples
341 ///
342 /// ```
343 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
344 ///
345 /// let m = from_vec2d(vec![
346 /// vec![1.0_f64, 2.0, 3.0],
347 /// vec![4.0, 5.0, 6.0],
348 /// vec![7.0, 8.0, 10.0],
349 /// ]);
350 /// let mut aca = MatrixACA::<f64>::new(3, 3);
351 /// let (r, c) = aca.add_best_pivot(&m).unwrap();
352 /// assert!(r < 3);
353 /// assert!(c < 3);
354 /// assert_eq!(aca.rank(), 1);
355 /// ```
356 pub fn add_best_pivot(&mut self, a: &Matrix<T>) -> Result<(usize, usize)> {
357 if self.is_empty() {
358 // Find global maximum
359 let (i, j, _) = submatrix_argmax(a, 0..self.nrows(), 0..self.ncols());
360 self.add_pivot(a, (i, j))?;
361 return Ok((i, j));
362 }
363
364 // ACA heuristic: find max in last row of V, then max in corresponding column of U
365 let avail_cols = self.available_cols();
366 if avail_cols.is_empty() {
367 return Err(MatrixCIError::FullRank);
368 }
369
370 // Find column with max |v[last, j]| among available columns
371 let last_row = nrows(&self.v) - 1;
372 let mut max_val = self.v[[last_row, avail_cols[0]]].abs_sq();
373 let mut best_col = avail_cols[0];
374 for &c in &avail_cols {
375 let val = self.v[[last_row, c]].abs_sq();
376 if val > max_val {
377 max_val = val;
378 best_col = c;
379 }
380 }
381
382 self.add_pivot_col(a, best_col)?;
383
384 // Find row with max |u[i, last]| among available rows
385 let avail_rows = self.available_rows();
386 if avail_rows.is_empty() {
387 return Err(MatrixCIError::FullRank);
388 }
389
390 let last_col = ncols(&self.u) - 1;
391 let mut max_val = self.u[[avail_rows[0], last_col]].abs_sq();
392 let mut best_row = avail_rows[0];
393 for &r in &avail_rows {
394 let val = self.u[[r, last_col]].abs_sq();
395 if val > max_val {
396 max_val = val;
397 best_row = r;
398 }
399 }
400
401 self.add_pivot_row(a, best_row)?;
402
403 Ok((best_row, best_col))
404 }
405
406 /// Set columns with new pivot rows and permutation
407 ///
408 /// Updates the column indices and V matrix according to the given
409 /// permutation and new pivot row data.
410 ///
411 /// # Examples
412 ///
413 /// ```
414 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
415 ///
416 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
417 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
418 /// // Identity permutation keeps columns the same
419 /// let pivot_rows = from_vec2d(vec![vec![1.0_f64, 2.0]]);
420 /// aca.set_cols(&pivot_rows, &[0, 1]);
421 /// assert_eq!(aca.col_indices(), &[0]);
422 /// ```
423 pub fn set_cols(&mut self, new_pivot_rows: &Matrix<T>, permutation: &[usize]) {
424 // Permute column indices
425 self.col_indices = self.col_indices.iter().map(|&c| permutation[c]).collect();
426
427 // Permute V matrix columns
428 let mut temp_v = zeros(nrows(&self.v), ncols(new_pivot_rows));
429 for i in 0..nrows(&self.v) {
430 for (new_j, &old_j) in permutation.iter().enumerate() {
431 if old_j < ncols(&self.v) {
432 temp_v[[i, new_j]] = self.v[[i, old_j]];
433 }
434 }
435 }
436 self.v = temp_v;
437
438 // Insert new elements
439 let new_indices: Vec<usize> = (0..ncols(new_pivot_rows))
440 .filter(|j| !permutation.contains(j))
441 .collect();
442
443 for k in 0..nrows(new_pivot_rows) {
444 for &j in &new_indices {
445 let mut val = new_pivot_rows[[k, j]];
446 for l in 0..k {
447 let factor = self.u[[self.row_indices[k], l]] * self.alpha[l];
448 val = val - self.v[[l, j]] * factor;
449 }
450 self.v[[k, j]] = val;
451 }
452 }
453 }
454
455 /// Set rows with new pivot columns and permutation
456 ///
457 /// Updates the row indices and U matrix according to the given
458 /// permutation and new pivot column data.
459 ///
460 /// # Examples
461 ///
462 /// ```
463 /// use tensor4all_tcicore::{AbstractMatrixCI, MatrixACA, from_vec2d};
464 ///
465 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
466 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
467 /// // Identity permutation keeps rows the same
468 /// let pivot_cols = from_vec2d(vec![vec![1.0_f64], vec![3.0]]);
469 /// aca.set_rows(&pivot_cols, &[0, 1]);
470 /// assert_eq!(aca.row_indices(), &[0]);
471 /// ```
472 pub fn set_rows(&mut self, new_pivot_cols: &Matrix<T>, permutation: &[usize]) {
473 // Permute row indices
474 self.row_indices = self.row_indices.iter().map(|&r| permutation[r]).collect();
475
476 // Permute U matrix rows
477 let mut temp_u = zeros(nrows(new_pivot_cols), ncols(&self.u));
478 for (new_i, &old_i) in permutation.iter().enumerate() {
479 if old_i < nrows(&self.u) {
480 for j in 0..ncols(&self.u) {
481 temp_u[[new_i, j]] = self.u[[old_i, j]];
482 }
483 }
484 }
485 self.u = temp_u;
486
487 // Insert new elements
488 let new_indices: Vec<usize> = (0..nrows(new_pivot_cols))
489 .filter(|i| !permutation.contains(i))
490 .collect();
491
492 for k in 0..ncols(new_pivot_cols) {
493 for &i in &new_indices {
494 let mut val = new_pivot_cols[[i, k]];
495 for l in 0..k {
496 let factor = self.v[[l, self.col_indices[k]]] * self.alpha[l];
497 val = val - self.u[[i, l]] * factor;
498 }
499 self.u[[i, k]] = val;
500 }
501 }
502 }
503}
504
505impl<T: Scalar> AbstractMatrixCI<T> for MatrixACA<T> {
506 fn nrows(&self) -> usize {
507 nrows(&self.u)
508 }
509
510 fn ncols(&self) -> usize {
511 ncols(&self.v)
512 }
513
514 fn rank(&self) -> usize {
515 self.row_indices.len()
516 }
517
518 fn row_indices(&self) -> &[usize] {
519 &self.row_indices
520 }
521
522 fn col_indices(&self) -> &[usize] {
523 &self.col_indices
524 }
525
526 fn evaluate(&self, i: usize, j: usize) -> T {
527 if self.is_empty() {
528 return T::zero();
529 }
530
531 let mut sum = T::zero();
532 for k in 0..self.rank() {
533 sum = sum + self.u[[i, k]] * self.alpha[k] * self.v[[k, j]];
534 }
535 sum
536 }
537
538 fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
539 if self.is_empty() {
540 return zeros(rows.len(), cols.len());
541 }
542
543 let r = self.rank();
544 let pivot_indices: Vec<usize> = (0..r).collect();
545 let mut left = submatrix(&self.u, rows, &pivot_indices);
546 for i in 0..rows.len() {
547 for k in 0..r {
548 left[[i, k]] = left[[i, k]] * self.alpha[k];
549 }
550 }
551 let right = submatrix(&self.v, &pivot_indices, cols);
552 mat_mul(&left, &right)
553 }
554}
555
556#[cfg(test)]
557mod tests;