|
| 1 | +///! Implements truncated singular value decomposition |
| 2 | +/// |
| 3 | +
|
| 4 | +use std::ops::DivAssign; |
| 5 | +use ndarray::prelude::*; |
| 6 | +use ndarray::stack; |
| 7 | +use ndarray_rand::rand_distr::Uniform; |
| 8 | +use ndarray_rand::RandomExt; |
| 9 | +use num_traits::{Float, NumCast}; |
| 10 | +use crate::{Scalar, Lapack}; |
| 11 | +use super::lobpcg::{lobpcg, EigResult, Order}; |
| 12 | +use crate::error::Result; |
| 13 | + |
| 14 | +#[derive(Debug)] |
| 15 | +pub struct TruncatedSvdResult<A> { |
| 16 | + eigvals: Array1<A>, |
| 17 | + eigvecs: Array2<A>, |
| 18 | + problem: Array2<A>, |
| 19 | + ngm: bool |
| 20 | +} |
| 21 | + |
| 22 | +impl<A: Float + PartialOrd + DivAssign<A> + 'static> TruncatedSvdResult<A> { |
| 23 | + fn singular_values_with_indices(&self) -> (Vec<A>, Vec<usize>) { |
| 24 | + let mut a = self.eigvals.iter() |
| 25 | + .map(|x| if *x < NumCast::from(1e-5).unwrap() { NumCast::from(0.0).unwrap() } else { *x }) |
| 26 | + .map(|x| x.sqrt()) |
| 27 | + .enumerate() |
| 28 | + .collect::<Vec<_>>(); |
| 29 | + |
| 30 | + a.sort_by(|(_,x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); |
| 31 | + |
| 32 | + a.into_iter().map(|(a,b)| (b,a)).unzip() |
| 33 | + } |
| 34 | + |
| 35 | + pub fn values(&self) -> Vec<A> { |
| 36 | + let (values, indices) = self.singular_values_with_indices(); |
| 37 | + |
| 38 | + values |
| 39 | + } |
| 40 | + |
| 41 | + pub fn values_vecs(&self) -> (Array2<A>, Vec<A>, Array2<A>) { |
| 42 | + let (values, indices) = self.singular_values_with_indices(); |
| 43 | + let n_values = values.iter().filter(|x| **x > NumCast::from(0.0).unwrap()).count(); |
| 44 | + |
| 45 | + if self.ngm { |
| 46 | + let vlarge = self.eigvecs.select(Axis(1), &indices); |
| 47 | + let mut ularge = self.problem.dot(&vlarge); |
| 48 | + |
| 49 | + ularge.gencolumns_mut().into_iter() |
| 50 | + .zip(values.iter()) |
| 51 | + .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); |
| 52 | + |
| 53 | + let vhlarge = vlarge.reversed_axes(); |
| 54 | + |
| 55 | + (vhlarge, values, ularge) |
| 56 | + } else { |
| 57 | + let ularge = self.eigvecs.select(Axis(1), &indices); |
| 58 | + |
| 59 | + let mut vlarge = ularge.dot(&self.problem); |
| 60 | + vlarge.gencolumns_mut().into_iter() |
| 61 | + .zip(values.iter()) |
| 62 | + .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); |
| 63 | + let vhlarge = vlarge.reversed_axes(); |
| 64 | + |
| 65 | + (vhlarge, values, ularge) |
| 66 | + } |
| 67 | + } |
| 68 | +} |
| 69 | + |
| 70 | +/// Truncated singular value decomposition |
| 71 | +/// |
| 72 | +/// This struct wraps the LOBPCG algorithm and provides convenient builder-pattern access to |
| 73 | +/// parameter like maximal iteration, precision and constraint matrix. Furthermore it allows |
| 74 | +/// conversion into a iterative solver where each iteration step yields a new eigenvalue/vector |
| 75 | +/// pair. |
| 76 | +pub struct TruncatedSvd<A: Scalar> { |
| 77 | + order: Order, |
| 78 | + problem: Array2<A>, |
| 79 | + precision: A::Real, |
| 80 | + maxiter: usize |
| 81 | +} |
| 82 | + |
| 83 | +impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedSvd<A> { |
| 84 | + pub fn new(problem: Array2<A>, order: Order) -> TruncatedSvd<A> { |
| 85 | + TruncatedSvd { |
| 86 | + precision: NumCast::from(1e-5).unwrap(), |
| 87 | + maxiter: problem.len_of(Axis(0)) * 2, |
| 88 | + order, |
| 89 | + problem |
| 90 | + } |
| 91 | + } |
| 92 | + |
| 93 | + pub fn precision(mut self, precision: A::Real) -> Self { |
| 94 | + self.precision = precision; |
| 95 | + |
| 96 | + self |
| 97 | + } |
| 98 | + |
| 99 | + pub fn maxiter(mut self, maxiter: usize) -> Self { |
| 100 | + self.maxiter = maxiter; |
| 101 | + |
| 102 | + self |
| 103 | + |
| 104 | + } |
| 105 | + |
| 106 | + // calculate the eigenvalues once |
| 107 | + pub fn once(&self, num: usize) -> Result<TruncatedSvdResult<A>> { |
| 108 | + let (n,m) = (self.problem.rows(), self.problem.ncols()); |
| 109 | + |
| 110 | + let x = Array2::random((usize::min(n,m), num), Uniform::new(0.0, 1.0)) |
| 111 | + .mapv(|x| NumCast::from(x).unwrap()); |
| 112 | + |
| 113 | + let res = if n > m { |
| 114 | + lobpcg(|y| self.problem.t().dot(&self.problem.dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) |
| 115 | + } else { |
| 116 | + lobpcg(|y| self.problem.dot(&self.problem.t().dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) |
| 117 | + }; |
| 118 | + |
| 119 | + match res { |
| 120 | + EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => { |
| 121 | + Ok(TruncatedSvdResult { |
| 122 | + problem: self.problem.clone(), |
| 123 | + eigvals: vals, |
| 124 | + eigvecs: vecs, |
| 125 | + ngm: n > m |
| 126 | + }) |
| 127 | + }, |
| 128 | + _ => panic!("") |
| 129 | + } |
| 130 | + } |
| 131 | +} |
| 132 | + |
| 133 | +#[cfg(test)] |
| 134 | +mod tests { |
| 135 | + use super::TruncatedSvd; |
| 136 | + use super::Order; |
| 137 | + use ndarray::{arr1, Array2}; |
| 138 | + |
| 139 | + #[test] |
| 140 | + fn test_truncated_svd() { |
| 141 | + let diag = arr1(&[ |
| 142 | + 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., |
| 143 | + ]); |
| 144 | + let a = Array2::from_diag(&diag); |
| 145 | + |
| 146 | + let res = TruncatedSvd::new(a, Order::Largest) |
| 147 | + .precision(1e-5) |
| 148 | + .maxiter(500) |
| 149 | + .once(3) |
| 150 | + .unwrap(); |
| 151 | + |
| 152 | + dbg!(&res.values()); |
| 153 | + } |
| 154 | +} |
0 commit comments