faer/linalg/svd/
mod.rs

1//! low level implementation of the svd of a matrix
2//!
3//! the svd of a matrix $A$ of shape $(m, n)$ is a decomposition into three
4//! components $U$, $S$, and $V$, such that:
5//!
6//! - $U$ has shape $(m, m)$ and is a unitary matrix
7//! - $V$ has shape $(n, n)$ and is a unitary matrix
8//! - $S$ has shape $(m, n)$ and is zero everywhere except the main diagonal
9//! - and finally:
10//!
11//! $$A = U S V^H$$
12use crate::assert;
13use crate::internal_prelude::*;
14use bidiag::BidiagParams;
15use linalg::qr::no_pivoting::factor::QrParams;
16/// bidiagonalization
17pub mod bidiag;
18pub(crate) mod bidiag_svd;
19/// whether the singular vectors should be computed
20#[derive(Copy, Clone, Debug, PartialEq, Eq)]
21pub enum ComputeSvdVectors {
22	/// do not compute singular vectors
23	No,
24	/// compute the first $\min(\text{nrows}, \text{ncols})$ singular vectors
25	Thin,
26	/// compute singular vectors
27	Full,
28}
29/// svd error
30#[derive(Copy, Clone, Debug, PartialEq, Eq)]
31pub enum SvdError {
32	/// reached max iterations
33	NoConvergence,
34}
35/// svd tuning parameters
36#[derive(Debug, Copy, Clone)]
37pub struct SvdParams {
38	/// bidiagonalization parameters
39	pub bidiag: BidiagParams,
40	/// $QR$ parameters
41	pub qr: QrParams,
42	/// threshold at which the implementation should stop recursing
43	pub recursion_threshold: usize,
44	/// threshold at which parallelism should be disabled
45	pub qr_ratio_threshold: f64,
46	#[doc(hidden)]
47	pub non_exhaustive: NonExhaustive,
48}
49impl<T: ComplexField> Auto<T> for SvdParams {
50	fn auto() -> Self {
51		Self {
52			recursion_threshold: 128,
53			qr_ratio_threshold: 11.0 / 6.0,
54			bidiag: auto!(T),
55			qr: auto!(T),
56			non_exhaustive: NonExhaustive(()),
57		}
58	}
59}
60fn svd_imp_scratch<T: ComplexField>(
61	m: usize,
62	n: usize,
63	compute_u: ComputeSvdVectors,
64	compute_v: ComputeSvdVectors,
65	bidiag_svd_scratch: fn(
66		n: usize,
67		compute_u: bool,
68		compute_v: bool,
69		par: Par,
70		params: SvdParams,
71	) -> StackReq,
72	params: SvdParams,
73	par: Par,
74) -> StackReq {
75	assert!(m >= n);
76	let householder_block_size =
77		linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
78	let bid = temp_mat_scratch::<T>(m, n);
79	let householder_left = temp_mat_scratch::<T>(householder_block_size, n);
80	let householder_right = temp_mat_scratch::<T>(householder_block_size, n);
81	let compute_bidiag =
82		bidiag::bidiag_in_place_scratch::<T>(m, n, par, params.bidiag.into());
83	let diag = temp_mat_scratch::<T>(n, 1);
84	let subdiag = diag;
85	let compute_ub = compute_v != ComputeSvdVectors::No;
86	let compute_vb = compute_u != ComputeSvdVectors::No;
87	let u_b = temp_mat_scratch::<T>(if compute_ub { n + 1 } else { 2 }, n + 1);
88	let v_b = temp_mat_scratch::<T>(n, if compute_vb { n } else { 0 });
89	let compute_bidiag_svd =
90		bidiag_svd_scratch(n, compute_ub, compute_vb, par, params);
91	let apply_householder_u = linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
92		m,
93		householder_block_size,
94		match compute_u {
95			ComputeSvdVectors::No => 0,
96			ComputeSvdVectors::Thin => n,
97			ComputeSvdVectors::Full => m,
98		},
99	);
100	let apply_householder_v = linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
101		n - 1,
102		householder_block_size,
103		match compute_v {
104			ComputeSvdVectors::No => 0,
105			_ => n,
106		},
107	);
108	StackReq::all_of(&[
109		bid,
110		householder_left,
111		householder_right,
112		StackReq::any_of(&[
113			compute_bidiag,
114			StackReq::all_of(&[
115				diag,
116				subdiag,
117				u_b,
118				v_b,
119				StackReq::any_of(&[
120					compute_bidiag_svd,
121					StackReq::all_of(&[
122						apply_householder_u,
123						apply_householder_v,
124					]),
125				]),
126			]),
127		]),
128	])
129}
130fn bidiag_cplx_svd_scratch<T: ComplexField>(
131	n: usize,
132	compute_u: bool,
133	compute_v: bool,
134	par: Par,
135	params: SvdParams,
136) -> StackReq {
137	StackReq::all_of(&[
138		temp_mat_scratch::<T>(n, 1).array(4),
139		temp_mat_scratch::<T::Real>(n + 1, if compute_u { n + 1 } else { 0 }),
140		temp_mat_scratch::<T::Real>(n, if compute_v { n } else { 0 }),
141		bidiag_real_svd_scratch::<T::Real>(
142			n, compute_u, compute_v, par, params,
143		),
144	])
145}
146fn bidiag_real_svd_scratch<T: RealField>(
147	n: usize,
148	compute_u: bool,
149	compute_v: bool,
150	par: Par,
151	params: SvdParams,
152) -> StackReq {
153	if n < params.recursion_threshold {
154		StackReq::EMPTY
155	} else {
156		StackReq::all_of(&[
157			temp_mat_scratch::<T>(2, if compute_u { 0 } else { n + 1 }),
158			bidiag_svd::divide_and_conquer_scratch::<T>(
159				n,
160				params.recursion_threshold,
161				compute_u,
162				compute_v,
163				par,
164			),
165		])
166	}
167}
168fn compute_bidiag_cplx_svd<T: ComplexField>(
169	mut diag: ColMut<'_, T, usize, ContiguousFwd>,
170	subdiag: ColMut<'_, T, usize, ContiguousFwd>,
171	mut u: Option<MatMut<'_, T>>,
172	mut v: Option<MatMut<'_, T>>,
173	params: SvdParams,
174	par: Par,
175	stack: &mut MemStack,
176) -> Result<(), SvdError> {
177	let n = diag.nrows();
178	let (mut diag_real, stack) =
179		unsafe { temp_mat_uninit::<T::Real, _, _>(n, 1, stack) };
180	let (mut subdiag_real, stack) =
181		unsafe { temp_mat_uninit::<T::Real, _, _>(n, 1, stack) };
182	let (mut u_real, stack) = unsafe {
183		temp_mat_uninit::<T::Real, _, _>(
184			n + 1,
185			if u.is_some() { n + 1 } else { 0 },
186			stack,
187		)
188	};
189	let (mut v_real, stack) = unsafe {
190		temp_mat_uninit::<T::Real, _, _>(
191			n,
192			if v.is_some() { n } else { 0 },
193			stack,
194		)
195	};
196	let (mut col_mul, stack) =
197		unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
198	let (mut row_mul, stack) =
199		unsafe { temp_mat_uninit::<T, _, _>(n - 1, 1, stack) };
200	let mut u_real = u.rb().map(|_| u_real.as_mat_mut());
201	let mut v_real = v.rb().map(|_| v_real.as_mat_mut());
202	let mut diag_real = diag_real
203		.as_mat_mut()
204		.col_mut(0)
205		.try_as_col_major_mut()
206		.unwrap();
207	let mut subdiag_real = subdiag_real
208		.as_mat_mut()
209		.col_mut(0)
210		.try_as_col_major_mut()
211		.unwrap();
212	let mut col_mul = col_mul.as_mat_mut().col_mut(0);
213	let mut row_mul = row_mul.as_mat_mut().col_mut(0);
214	let normalized = |x: T| {
215		if x == zero() {
216			one()
217		} else {
218			let norm1 = x.real().abs().fmax(x.imag().abs());
219			let y = x.mul_real(norm1.recip());
220			y.mul_real(y.abs().recip())
221		}
222	};
223	let mut col_normalized = normalized(diag[0].conj());
224	col_mul[0] = col_normalized.copy();
225	diag_real[0] = diag[0].abs();
226	subdiag_real[n - 1] = zero();
227	for i in 1..n {
228		let row_normalized =
229			normalized((&subdiag[i - 1] * col_normalized).conj());
230		subdiag_real[i - 1] = subdiag[i - 1].abs();
231		row_mul[i - 1] = row_normalized.conj();
232		col_normalized = normalized((&diag[i] * row_normalized).conj());
233		diag_real[i] = diag[i].abs();
234		col_mul[i] = col_normalized.copy();
235	}
236	compute_bidiag_real_svd(
237		diag_real.rb_mut(),
238		subdiag_real.rb_mut(),
239		u_real.rb_mut(),
240		v_real.rb_mut(),
241		params,
242		par,
243		stack,
244	)?;
245	for i in 0..n {
246		diag[i] = diag_real[i].to_cplx();
247	}
248	let u_real = u_real.rb();
249	let v_real = v_real.rb();
250	if let (Some(mut u), Some(u_real)) = (u.rb_mut(), u_real) {
251		z!(u.rb_mut().row_mut(0), u_real.row(0))
252			.for_each(|uz!(u, r)| *u = r.to_cplx());
253		z!(u.rb_mut().row_mut(n), u_real.row(n))
254			.for_each(|uz!(u, r)| *u = r.to_cplx());
255		for j in 0..u.ncols() {
256			let mut u = u.rb_mut().col_mut(j).subrows_mut(1, n - 1);
257			let u_real = u_real.rb().col(j).subrows(1, n - 1);
258			z!(u.rb_mut(), u_real, row_mul.rb())
259				.for_each(|uz!(u, re, f)| *u = f.mul_real(re));
260		}
261	}
262	if let (Some(mut v), Some(v_real)) = (v.rb_mut(), v_real) {
263		for j in 0..v.ncols() {
264			let mut v = v.rb_mut().col_mut(j);
265			let v_real = v_real.rb().col(j);
266			z!(v.rb_mut(), v_real, col_mul.rb())
267				.for_each(|uz!(v, re, f)| *v = f.mul_real(re));
268		}
269	}
270	Ok(())
271}
272fn compute_bidiag_real_svd<T: RealField>(
273	mut diag: ColMut<'_, T, usize, ContiguousFwd>,
274	mut subdiag: ColMut<'_, T, usize, ContiguousFwd>,
275	mut u: Option<MatMut<'_, T, usize, usize>>,
276	mut v: Option<MatMut<'_, T, usize, usize>>,
277	params: SvdParams,
278	par: Par,
279	stack: &mut MemStack,
280) -> Result<(), SvdError> {
281	let n = diag.nrows();
282	for i in 0..n {
283		if !(diag[i].is_finite() && subdiag[i].is_finite()) {
284			return Err(SvdError::NoConvergence);
285		}
286	}
287	if n < params.recursion_threshold {
288		if let Some(mut u) = u.rb_mut() {
289			u.fill(zero());
290			u.diagonal_mut().fill(one());
291		}
292		if let Some(mut v) = v.rb_mut() {
293			v.fill(zero());
294			v.diagonal_mut().fill(one());
295		}
296		bidiag_svd::qr_algorithm(
297			diag.rb_mut(),
298			subdiag.rb_mut(),
299			u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
300			v.rb_mut(),
301		)?;
302		return Ok(());
303	} else {
304		let (mut u2, stack) = unsafe {
305			temp_mat_uninit::<T::Real, _, _>(
306				2,
307				if u.is_some() { 0 } else { n + 1 },
308				stack,
309			)
310		};
311		bidiag_svd::divide_and_conquer(
312			diag.as_row_shape_mut(n),
313			subdiag.as_row_shape_mut(n),
314			match u {
315				Some(u) => bidiag_svd::MatU::Full(u),
316				None => bidiag_svd::MatU::TwoRowsStorage(u2.as_mat_mut()),
317			},
318			v.map(|m| m.as_shape_mut(n, n)),
319			par,
320			stack,
321			params.recursion_threshold,
322		)
323	}
324}
325/// bidiag -> divide conquer svd / qr algo
326fn svd_imp<T: ComplexField>(
327	matrix: MatRef<'_, T>,
328	s: ColMut<'_, T>,
329	u: Option<MatMut<'_, T>>,
330	v: Option<MatMut<'_, T>>,
331	bidiag_svd: fn(
332		diag: ColMut<'_, T, usize, ContiguousFwd>,
333		subdiag: ColMut<'_, T, usize, ContiguousFwd>,
334		u: Option<MatMut<'_, T, usize, usize>>,
335		v: Option<MatMut<'_, T, usize, usize>>,
336		params: SvdParams,
337		par: Par,
338		stack: &mut MemStack,
339	) -> Result<(), SvdError>,
340	par: Par,
341	stack: &mut MemStack,
342	params: SvdParams,
343) -> Result<(), SvdError> {
344	assert!(matrix.nrows() >= matrix.ncols());
345	let m = matrix.nrows();
346	let n = matrix.ncols();
347	let bs = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
348	let (mut bid, stack) = unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
349	let mut bid = bid.as_mat_mut();
350	let (mut Hl, stack) = unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
351	let (mut Hr, stack) =
352		unsafe { temp_mat_uninit::<T, _, _>(bs, n - 1, stack) };
353	let mut Hl = Hl.as_mat_mut();
354	let mut Hr = Hr.as_mat_mut();
355	bid.copy_from(matrix);
356	bidiag::bidiag_in_place(
357		bid.rb_mut(),
358		Hl.rb_mut(),
359		Hr.rb_mut(),
360		par,
361		stack,
362		params.bidiag.into(),
363	);
364	let (mut diag, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
365	let (mut subdiag, stack) =
366		unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
367	let mut diag = diag.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
368	let mut subdiag = subdiag
369		.as_mat_mut()
370		.col_mut(0)
371		.try_as_col_major_mut()
372		.unwrap();
373	let (mut ub, stack) = unsafe {
374		temp_mat_uninit::<T, _, _>(
375			n + 1,
376			if v.is_some() { n + 1 } else { 0 },
377			stack,
378		)
379	};
380	let (mut vb, stack) = unsafe {
381		temp_mat_uninit::<T, _, _>(n, if u.is_some() { n } else { 0 }, stack)
382	};
383	let mut ub = ub.as_mat_mut();
384	let mut vb = vb.as_mat_mut();
385	for i in 0..n {
386		diag[i] = bid[(i, i)].conj();
387		if i + 1 < n {
388			subdiag[i] = bid[(i, i + 1)].conj();
389		} else {
390			subdiag[i] = zero();
391		}
392	}
393	bidiag_svd(
394		diag.rb_mut(),
395		subdiag.rb_mut(),
396		v.rb().map(|_| ub.rb_mut()),
397		u.rb().map(|_| vb.rb_mut()),
398		params,
399		par,
400		stack,
401	)?;
402	{ s }.copy_from(diag);
403	if let Some(mut u) = u {
404		let ncols = u.ncols();
405		u.rb_mut().submatrix_mut(0, 0, n, n).copy_from(vb.rb());
406		u.rb_mut().submatrix_mut(n, 0, m - n, ncols).fill(zero());
407		u.rb_mut().submatrix_mut(0, n, n, ncols - n).fill(zero());
408		u.rb_mut()
409			.submatrix_mut(n, n, ncols - n, ncols - n)
410			.diagonal_mut()
411			.fill(one());
412		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(bid.rb(), Hl.rb(), Conj::No, u, par, stack);
413	}
414	if let Some(mut v) = v {
415		v.copy_from(ub.rb().submatrix(0, 0, n, n));
416		for j in 1..n {
417			for i in 0..j {
418				bid[(j, i)] = bid[(i, j)].copy();
419			}
420		}
421		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
422			bid.rb().submatrix(1, 0, n - 1, n - 1),
423			Hr.rb(),
424			Conj::Yes,
425			v.subrows_mut(1, n - 1),
426			par,
427			stack,
428		);
429	}
430	Ok(())
431}
432fn compute_squareish_svd<T: ComplexField>(
433	matrix: MatRef<'_, T>,
434	s: ColMut<'_, T>,
435	u: Option<MatMut<'_, T>>,
436	v: Option<MatMut<'_, T>>,
437	par: Par,
438	stack: &mut MemStack,
439	params: SvdParams,
440) -> Result<(), SvdError> {
441	if const { T::IS_REAL } {
442		svd_imp::<T::Real>(
443			unsafe { core::mem::transmute(matrix) },
444			unsafe { core::mem::transmute(s) },
445			unsafe { core::mem::transmute(u) },
446			unsafe { core::mem::transmute(v) },
447			compute_bidiag_real_svd::<T::Real>,
448			par,
449			stack,
450			params,
451		)
452	} else {
453		svd_imp::<T>(
454			matrix,
455			s,
456			u,
457			v,
458			compute_bidiag_cplx_svd::<T>,
459			par,
460			stack,
461			params,
462		)
463	}
464}
465/// computes the layout of the workspace required to compute a matrix's svd
466pub fn svd_scratch<T: ComplexField>(
467	nrows: usize,
468	ncols: usize,
469	compute_u: ComputeSvdVectors,
470	compute_v: ComputeSvdVectors,
471	par: Par,
472	params: Spec<SvdParams, T>,
473) -> StackReq {
474	let params = params.config;
475	let mut m = nrows;
476	let mut n = ncols;
477	let mut compute_u = compute_u;
478	let mut compute_v = compute_v;
479	if n > m {
480		core::mem::swap(&mut m, &mut n);
481		core::mem::swap(&mut compute_u, &mut compute_v);
482	}
483	if n == 0 {
484		return StackReq::EMPTY;
485	}
486	let bidiag_svd_scratch = if const { T::IS_REAL } {
487		bidiag_real_svd_scratch::<T::Real>
488	} else {
489		bidiag_cplx_svd_scratch::<T>
490	};
491	if m as f64 / n as f64 <= params.qr_ratio_threshold {
492		svd_imp_scratch::<T>(
493			m,
494			n,
495			compute_u,
496			compute_v,
497			bidiag_svd_scratch,
498			params,
499			par,
500		)
501	} else {
502		let bs =
503			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
504		StackReq::all_of(&[
505			temp_mat_scratch::<T>(m, n),
506			temp_mat_scratch::<T>(bs, n),
507			StackReq::any_of(&[
508				StackReq::all_of(&[
509					temp_mat_scratch::<T>(n, n),
510					svd_imp_scratch::<T>(n, n, compute_u, compute_v, bidiag_svd_scratch, params, par),
511				]),
512				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
513					m,
514					bs,
515					match compute_u {
516						ComputeSvdVectors::No => 0,
517						ComputeSvdVectors::Thin => n,
518						ComputeSvdVectors::Full => m,
519					},
520				),
521			]),
522		])
523	}
524}
525/// computes the svd of $A$, with the singular vectors being omitted, thin or
526/// full
527///
528/// the singular are stored in $S$, and the singular vectors in $U$ and $V$ such
529/// that the singular values are sorted in nonincreasing order
530pub fn svd<T: ComplexField>(
531	A: MatRef<'_, T>,
532	s: DiagMut<'_, T>,
533	u: Option<MatMut<'_, T>>,
534	v: Option<MatMut<'_, T>>,
535	par: Par,
536	stack: &mut MemStack,
537	params: Spec<SvdParams, T>,
538) -> Result<(), SvdError> {
539	let params = params.config;
540	let (m, n) = A.shape();
541	let size = Ord::min(m, n);
542	assert!(s.dim() == size);
543	let s = s.column_vector_mut();
544	if let Some(u) = u.rb() {
545		assert!(all(
546			u.nrows() == A.nrows(),
547			any(u.ncols() == A.nrows(), u.ncols() == size),
548		));
549	}
550	if let Some(v) = v.rb() {
551		assert!(all(
552			v.nrows() == A.ncols(),
553			any(v.ncols() == A.ncols(), v.ncols() == size),
554		));
555	}
556	#[cfg(feature = "perf-warn")]
557	match (u.rb(), v.rb()) {
558		(Some(matrix), _) | (_, Some(matrix)) => {
559			if matrix.row_stride().unsigned_abs() != 1
560				&& crate::__perf_warn!(QR_WARN)
561			{
562				if matrix.col_stride().unsigned_abs() == 1 {
563					log::warn!(
564						target : "faer_perf",
565						"SVD prefers column-major singular vector matrices. Found row-major matrix."
566					);
567				} else {
568					log::warn!(
569						target : "faer_perf",
570						"SVD prefers column-major singular vector matrices. Found matrix with generic strides."
571					);
572				}
573			}
574		},
575		_ => {},
576	}
577	let mut u = u;
578	let mut v = v;
579	let mut matrix = A;
580	let do_transpose = n > m;
581	if do_transpose {
582		matrix = matrix.transpose();
583		core::mem::swap(&mut u, &mut v)
584	}
585	let (m, n) = matrix.shape();
586	if n == 0 {
587		if let Some(mut u) = u {
588			u.fill(zero());
589			u.rb_mut().diagonal_mut().fill(one());
590		}
591		return Ok(());
592	}
593	if m as f64 / n as f64 <= params.qr_ratio_threshold {
594		compute_squareish_svd(
595			matrix,
596			s,
597			u.rb_mut(),
598			v.rb_mut(),
599			par,
600			stack,
601			params,
602		)?;
603	} else {
604		let bs =
605			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
606		let (mut qr, stack) =
607			unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
608		let mut qr = qr.as_mat_mut();
609		let (mut householder, stack) =
610			unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
611		let mut householder = householder.as_mat_mut();
612		{
613			qr.copy_from(matrix.rb());
614			linalg::qr::no_pivoting::factor::qr_in_place(
615				qr.rb_mut(),
616				householder.rb_mut(),
617				par,
618				stack,
619				params.qr.into(),
620			);
621		}
622		{
623			let (mut r, stack) =
624				unsafe { temp_mat_uninit::<T, _, _>(n, n, stack) };
625			let mut r = r.as_mat_mut();
626			z!(r.rb_mut()).for_each_triangular_lower(
627				linalg::zip::Diag::Skip,
628				|uz!(dst)| *dst = zero(),
629			);
630			z!(r.rb_mut(), qr.rb().submatrix(0, 0, n, n))
631				.for_each_triangular_upper(
632					linalg::zip::Diag::Include,
633					|uz!(dst, src)| *dst = src.copy(),
634				);
635			compute_squareish_svd(
636				r.rb(),
637				s,
638				u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
639				v.rb_mut(),
640				par,
641				stack,
642				params,
643			)?;
644		}
645		if let Some(mut u) = u.rb_mut() {
646			u.rb_mut().subrows_mut(n, m - n).fill(zero());
647			if u.ncols() == m {
648				u.rb_mut()
649					.submatrix_mut(n, n, m - n, m - n)
650					.diagonal_mut()
651					.fill(one());
652			}
653			linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
654				qr.rb(),
655				householder.rb(),
656				Conj::No,
657				u.rb_mut(),
658				par,
659				stack,
660			);
661		}
662	}
663	if do_transpose {
664		if let Some(u) = u.rb_mut() {
665			z!(u).for_each(|uz!(u)| *u = u.conj())
666		}
667		if let Some(v) = v.rb_mut() {
668			z!(v).for_each(|uz!(v)| *v = v.conj())
669		}
670	}
671	Ok(())
672}
673/// computes the layout of the workspace required to compute a matrix's
674/// pseudoinverse, given the svd
675pub fn pseudoinverse_from_svd_scratch<T: ComplexField>(
676	nrows: usize,
677	ncols: usize,
678	par: Par,
679) -> StackReq {
680	_ = par;
681	let size = Ord::min(nrows, ncols);
682	StackReq::all_of(&[
683		temp_mat_scratch::<T>(nrows, size),
684		temp_mat_scratch::<T>(ncols, size),
685	])
686}
687/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
688/// $U$ and $V$
689pub fn pseudoinverse_from_svd<T: ComplexField>(
690	pinv: MatMut<'_, T>,
691	s: DiagRef<'_, T>,
692	u: MatRef<'_, T>,
693	v: MatRef<'_, T>,
694	par: Par,
695	stack: &mut MemStack,
696) {
697	pseudoinverse_from_svd_with_tolerance(
698		pinv,
699		s,
700		u,
701		v,
702		zero(),
703		eps::<T::Real>()
704			* from_f64::<T::Real>(Ord::max(u.nrows(), v.nrows()) as f64),
705		par,
706		stack,
707	);
708}
709/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
710/// $U$ and $V$, and tolerance parameters for determining zero singular values
711pub fn pseudoinverse_from_svd_with_tolerance<T: ComplexField>(
712	pinv: MatMut<'_, T>,
713	s: DiagRef<'_, T>,
714	u: MatRef<'_, T>,
715	v: MatRef<'_, T>,
716	abs_tol: T::Real,
717	rel_tol: T::Real,
718	par: Par,
719	stack: &mut MemStack,
720) {
721	let mut pinv = pinv;
722	let m = u.nrows();
723	let n = v.nrows();
724	let size = Ord::min(m, n);
725	assert!(all(
726		u.nrows() == m,
727		v.nrows() == n,
728		u.ncols() >= size,
729		v.ncols() >= size,
730		s.dim() >= size
731	));
732	let s = s.column_vector();
733	let u = u.get(.., ..size);
734	let v = v.get(.., ..size);
735	let smax = s.norm_max();
736	let tol = abs_tol.fmax(rel_tol * smax);
737	let (mut u_trunc, stack) =
738		unsafe { temp_mat_uninit::<T, _, _>(m, size, stack) };
739	let (mut vp_trunc, _) =
740		unsafe { temp_mat_uninit::<T, _, _>(n, size, stack) };
741	let mut u_trunc = u_trunc.as_mat_mut();
742	let mut vp_trunc = vp_trunc.as_mat_mut();
743	let mut len = 0;
744	for j in 0..size {
745		let x = s[j].absmax();
746		if x > tol {
747			let ref p = s[j].real().recip();
748			u_trunc.rb_mut().col_mut(len).copy_from(u.col(j));
749			z!(vp_trunc.rb_mut().col_mut(len), v.col(j))
750				.for_each(|uz!(dst, src)| *dst = src.mul_real(p));
751			len += 1;
752		}
753	}
754	let u_trunc = u_trunc.get(.., ..len);
755	let vp_trunc = vp_trunc.get(.., ..len);
756	linalg::matmul::matmul(
757		pinv.rb_mut(),
758		Accum::Replace,
759		vp_trunc,
760		u_trunc.adjoint(),
761		one(),
762		par,
763	);
764}
765#[cfg(test)]
766mod tests {
767	use super::*;
768	use crate::assert;
769	use crate::stats::prelude::*;
770	use crate::utils::approx::*;
771	use dyn_stack::MemBuffer;
772	#[track_caller]
773	fn test_svd<T: ComplexField>(mat: MatRef<'_, T>) {
774		let (m, n) = mat.shape();
775		let params = Spec::new(SvdParams {
776			recursion_threshold: 8,
777			qr_ratio_threshold: 1.0,
778			..auto!(T)
779		});
780		let approx_eq = CwiseMat(
781			ApproxEq::<T::Real>::eps()
782				* from_f64::<T::Real>(8.0 * Ord::max(m, n) as f64).sqrt(),
783		);
784		{
785			let mut s = Mat::zeros(m, n);
786			let mut u = Mat::zeros(m, m);
787			let mut v = Mat::zeros(n, n);
788			svd(
789				mat.as_ref(),
790				s.as_mut().diagonal_mut(),
791				Some(u.as_mut()),
792				Some(v.as_mut()),
793				Par::Seq,
794				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
795					m,
796					n,
797					ComputeSvdVectors::Full,
798					ComputeSvdVectors::Full,
799					Par::Seq,
800					params,
801				))),
802				params,
803			)
804			.unwrap();
805			let reconstructed = &u * &s * v.adjoint();
806			assert!(reconstructed ~ mat);
807		}
808		let size = Ord::min(m, n);
809		let mut s = Mat::zeros(size, size);
810		let mut u = Mat::zeros(m, size);
811		let mut v = Mat::zeros(n, size);
812		{
813			svd(
814				mat.as_ref(),
815				s.as_mut().diagonal_mut(),
816				Some(u.as_mut()),
817				Some(v.as_mut()),
818				Par::Seq,
819				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
820					m,
821					n,
822					ComputeSvdVectors::Thin,
823					ComputeSvdVectors::Thin,
824					Par::Seq,
825					params,
826				))),
827				params,
828			)
829			.unwrap();
830			let reconstructed = &u * &s * v.adjoint();
831			assert!(reconstructed ~ mat);
832		}
833		{
834			let mut s2 = Mat::zeros(size, size);
835			let mut u2 = Mat::zeros(m, size);
836			svd(
837				mat.as_ref(),
838				s2.as_mut().diagonal_mut(),
839				Some(u2.as_mut()),
840				None,
841				Par::Seq,
842				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
843					m,
844					n,
845					ComputeSvdVectors::Thin,
846					ComputeSvdVectors::No,
847					Par::Seq,
848					params,
849				))),
850				params,
851			)
852			.unwrap();
853			assert!(s2 ~ s);
854			assert!(u2 ~ u);
855		}
856		{
857			let mut s2 = Mat::zeros(size, size);
858			let mut v2 = Mat::zeros(n, size);
859			svd(
860				mat.as_ref(),
861				s2.as_mut().diagonal_mut(),
862				None,
863				Some(v2.as_mut()),
864				Par::Seq,
865				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
866					m,
867					n,
868					ComputeSvdVectors::No,
869					ComputeSvdVectors::Thin,
870					Par::Seq,
871					params,
872				))),
873				params,
874			)
875			.unwrap();
876			assert!(s2 ~ s);
877			assert!(v2 ~ v);
878		}
879		{
880			let mut s2 = Mat::zeros(size, size);
881			svd(
882				mat.as_ref(),
883				s2.as_mut().diagonal_mut(),
884				None,
885				None,
886				Par::Seq,
887				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
888					m,
889					n,
890					ComputeSvdVectors::No,
891					ComputeSvdVectors::No,
892					Par::Seq,
893					params,
894				))),
895				params,
896			)
897			.unwrap();
898			assert!(s2 ~ s);
899		}
900	}
901	#[test]
902	fn test_real() {
903		let rng = &mut StdRng::seed_from_u64(1);
904		for (m, n) in [
905			(3, 2),
906			(2, 2),
907			(4, 4),
908			(15, 10),
909			(10, 10),
910			(15, 15),
911			(50, 50),
912			(100, 100),
913			(150, 150),
914			(150, 20),
915			(20, 150),
916		] {
917			let mat = CwiseMatDistribution {
918				nrows: m,
919				ncols: n,
920				dist: StandardNormal,
921			}
922			.rand::<Mat<f64>>(rng);
923			test_svd(mat.as_ref());
924		}
925	}
926	#[test]
927	fn test_cplx() {
928		let rng = &mut StdRng::seed_from_u64(1);
929		for (m, n) in [
930			(1, 1),
931			(2, 2),
932			(3, 2),
933			(2, 2),
934			(3, 3),
935			(4, 4),
936			(15, 10),
937			(10, 10),
938			(15, 15),
939			(16, 16),
940			(17, 17),
941			(18, 18),
942			(19, 19),
943			(20, 20),
944			(30, 30),
945			(50, 50),
946			(100, 100),
947			(150, 150),
948			(150, 20),
949			(20, 150),
950		] {
951			let mat = CwiseMatDistribution {
952				nrows: m,
953				ncols: n,
954				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
955			}
956			.rand::<Mat<c64>>(rng);
957			test_svd(mat.as_ref());
958		}
959	}
960	#[test]
961	fn test_special() {
962		for (m, n) in [
963			(3, 2),
964			(2, 2),
965			(4, 4),
966			(15, 10),
967			(10, 10),
968			(15, 15),
969			(50, 50),
970			(100, 100),
971			(150, 150),
972			(150, 20),
973			(20, 150),
974		] {
975			test_svd(Mat::<f64>::zeros(m, n).as_ref());
976			test_svd(Mat::<c64>::zeros(m, n).as_ref());
977			test_svd(Mat::<f64>::full(m, n, 1.0).as_ref());
978			test_svd(Mat::<c64>::full(m, n, c64::ONE).as_ref());
979			test_svd(Mat::<f64>::identity(m, n).as_ref());
980			test_svd(Mat::<c64>::identity(m, n).as_ref());
981		}
982	}
983	#[test]
984	fn test_zink() {
985		let diag = [
986			-9.931833701529301,
987			-10.920807536026027,
988			-52.33647796311243,
989			2.3685025127736967,
990			2.421701994236093,
991			-0.5051763005624579,
992			-0.04808263896606017,
993			-0.003875251886338955,
994			-0.0006413264967716465,
995			-0.003381944152463707,
996			2.981152313236375e-5,
997			5.4290648208388795e-6,
998			-6.329275972084404e-7,
999			-6.879142344209158e-7,
1000			-5.265228263479126e-9,
1001			-2.941999902335516e-9,
1002			-1.3060984997930294e-10,
1003			7.07516117218088e-12,
1004			1.8657003929029376e-12,
1005			-6.216080089659131e-14,
1006		];
1007		let subdiag = [
1008			-57.8029649868477,
1009			17.67263066467847,
1010			8.884153814270894,
1011			-9.01998231080713,
1012			-1.028638150814966,
1013			0.22247719217200435,
1014			0.016389886745811315,
1015			-0.004090989452162578,
1016			0.00036818904090536926,
1017			-0.0031394146217732367,
1018			-7.571300829706796e-6,
1019			3.0045718718618155e-6,
1020			2.1329796886727743e-6,
1021			9.259701025627789e-8,
1022			2.2291214755992877e-9,
1023			-2.3017207713252894e-9,
1024			6.807967994979358e-11,
1025			2.1677299575405587e-12,
1026			-3.07282771050034e-13,
1027			0.0,
1028		];
1029		let n = diag.len();
1030		let params = SvdParams {
1031			recursion_threshold: 8,
1032			qr_ratio_threshold: 1.0,
1033			..auto!(f64)
1034		};
1035		let mut d = ColRef::from_slice(&diag).to_owned();
1036		let mut s = ColRef::from_slice(&subdiag).to_owned();
1037		compute_bidiag_real_svd(
1038			d.as_mut().try_as_col_major_mut().unwrap(),
1039			s.as_mut().try_as_col_major_mut().unwrap(),
1040			None,
1041			None,
1042			params,
1043			Par::Seq,
1044			MemStack::new(&mut MemBuffer::new(bidiag_real_svd_scratch::<f64>(
1045				n,
1046				false,
1047				false,
1048				Par::Seq,
1049				params,
1050			))),
1051		)
1052		.unwrap();
1053		assert!(d[n - 1] != 0.0);
1054	}
1055	#[test]
1056	fn test_pinv() {
1057		let rng = &mut StdRng::seed_from_u64(0);
1058		let m = 6;
1059		let n = 36;
1060		let mat = CwiseMatDistribution {
1061			nrows: m,
1062			ncols: n,
1063			dist: StandardNormal,
1064		}
1065		.rand::<Mat<f64>>(rng);
1066		let pinv = mat.svd().unwrap().pseudoinverse();
1067		let err = &mat * &pinv - Mat::<f64>::identity(m, m);
1068		assert!(err.norm_max() < 1e-10);
1069	}
1070	#[test]
1071	fn test_pinv2() {
1072		let rng = &mut StdRng::seed_from_u64(0);
1073		let m = 36;
1074		let n = 6;
1075		let mat = CwiseMatDistribution {
1076			nrows: m,
1077			ncols: n,
1078			dist: StandardNormal,
1079		}
1080		.rand::<Mat<f64>>(rng);
1081		let pinv = mat.svd().unwrap().pseudoinverse();
1082		let err = &pinv * &mat - Mat::<f64>::identity(n, n);
1083		assert!(err.norm_max() < 1e-10);
1084	}
1085}
faer/linalg/svd/mod.rs

faer/linalg/svd/
mod.rs
