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	fn check_tiny_split<T: ComplexField>(a: MatRef<'_, T>) {
902		let (m, n) = a.shape();
903		let size = Ord::min(m, n);
904		let compute = |threshold, want_u, want_v| {
905			let params = Spec::new(SvdParams {
906				recursion_threshold: threshold,
907				..auto!(T)
908			});
909			let mut s = Col::zeros(size);
910			let mut u = Mat::zeros(m, if want_u { size } else { 0 });
911			let mut v = Mat::zeros(n, if want_v { size } else { 0 });
912			let vectors = |want| {
913				if want {
914					ComputeSvdVectors::Thin
915				} else {
916					ComputeSvdVectors::No
917				}
918			};
919			svd(
920				a,
921				s.as_mut().as_diagonal_mut(),
922				if want_u { Some(u.as_mut()) } else { None },
923				if want_v { Some(v.as_mut()) } else { None },
924				Par::Seq,
925				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
926					m,
927					n,
928					vectors(want_u),
929					vectors(want_v),
930					Par::Seq,
931					params,
932				))),
933				params,
934			)
935			.unwrap();
936			(s, u, v)
937		};
938		let (reference, _, _) = compute(usize::MAX, true, true);
939		for (want_u, want_v) in
940			[(true, true), (false, false), (true, false), (false, true)]
941		{
942			let (s, u, v) = compute(128, want_u, want_v);
943			assert!((&s - &reference).norm_max() < from_f64::<T::Real>(1e-13));
944			if want_u {
945				assert!(
946					(u.adjoint() * &u - Mat::<T>::identity(size, size))
947						.norm_l2() < from_f64::<T::Real>(1e-13)
948				);
949			}
950			if want_v {
951				assert!(
952					(v.adjoint() * &v - Mat::<T>::identity(size, size))
953						.norm_l2() < from_f64::<T::Real>(1e-13)
954				);
955			}
956			if want_u && want_v {
957				let back = &u * s.as_diagonal() * v.adjoint();
958				assert!(
959					(&back - a).norm_l2() / a.norm_l2()
960						< from_f64::<T::Real>(1e-13)
961				);
962			}
963		}
964	}
965	#[test]
966	fn test_tiny_diagonal_at_split() {
967		for n in [128, 256] {
968			for tiny in [0.0, 1e-15, 1e-13, 1e-10, 1e-6] {
969				let a = Mat::<f64>::from_fn(n, n, |i, j| {
970					if (i == j && i < 64) || (j == i + 1 && i < 64) {
971						1.0
972					} else if i == 64 && j == 64 {
973						tiny
974					} else {
975						0.0
976					}
977				});
978				check_tiny_split(a.as_ref());
979				let complex =
980					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
981				check_tiny_split(complex.as_ref());
982			}
983		}
984	}
985	#[test]
986	fn test_real() {
987		let rng = &mut StdRng::seed_from_u64(1);
988		for (m, n) in [
989			(3, 2),
990			(2, 2),
991			(4, 4),
992			(15, 10),
993			(10, 10),
994			(15, 15),
995			(50, 50),
996			(100, 100),
997			(150, 150),
998			(150, 20),
999			(20, 150),
1000		] {
1001			let mat = CwiseMatDistribution {
1002				nrows: m,
1003				ncols: n,
1004				dist: StandardNormal,
1005			}
1006			.rand::<Mat<f64>>(rng);
1007			test_svd(mat.as_ref());
1008		}
1009	}
1010	#[test]
1011	fn test_cplx() {
1012		let rng = &mut StdRng::seed_from_u64(1);
1013		for (m, n) in [
1014			(1, 1),
1015			(2, 2),
1016			(3, 2),
1017			(2, 2),
1018			(3, 3),
1019			(4, 4),
1020			(15, 10),
1021			(10, 10),
1022			(15, 15),
1023			(16, 16),
1024			(17, 17),
1025			(18, 18),
1026			(19, 19),
1027			(20, 20),
1028			(30, 30),
1029			(50, 50),
1030			(100, 100),
1031			(150, 150),
1032			(150, 20),
1033			(20, 150),
1034		] {
1035			let mat = CwiseMatDistribution {
1036				nrows: m,
1037				ncols: n,
1038				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
1039			}
1040			.rand::<Mat<c64>>(rng);
1041			test_svd(mat.as_ref());
1042		}
1043	}
1044	#[test]
1045	fn test_special() {
1046		for (m, n) in [
1047			(3, 2),
1048			(2, 2),
1049			(4, 4),
1050			(15, 10),
1051			(10, 10),
1052			(15, 15),
1053			(50, 50),
1054			(100, 100),
1055			(150, 150),
1056			(150, 20),
1057			(20, 150),
1058		] {
1059			test_svd(Mat::<f64>::zeros(m, n).as_ref());
1060			test_svd(Mat::<c64>::zeros(m, n).as_ref());
1061			test_svd(Mat::<f64>::full(m, n, 1.0).as_ref());
1062			test_svd(Mat::<c64>::full(m, n, c64::ONE).as_ref());
1063			test_svd(Mat::<f64>::identity(m, n).as_ref());
1064			test_svd(Mat::<c64>::identity(m, n).as_ref());
1065		}
1066	}
1067	#[test]
1068	fn test_zink() {
1069		let diag = [
1070			-9.931833701529301,
1071			-10.920807536026027,
1072			-52.33647796311243,
1073			2.3685025127736967,
1074			2.421701994236093,
1075			-0.5051763005624579,
1076			-0.04808263896606017,
1077			-0.003875251886338955,
1078			-0.0006413264967716465,
1079			-0.003381944152463707,
1080			2.981152313236375e-5,
1081			5.4290648208388795e-6,
1082			-6.329275972084404e-7,
1083			-6.879142344209158e-7,
1084			-5.265228263479126e-9,
1085			-2.941999902335516e-9,
1086			-1.3060984997930294e-10,
1087			7.07516117218088e-12,
1088			1.8657003929029376e-12,
1089			-6.216080089659131e-14,
1090		];
1091		let subdiag = [
1092			-57.8029649868477,
1093			17.67263066467847,
1094			8.884153814270894,
1095			-9.01998231080713,
1096			-1.028638150814966,
1097			0.22247719217200435,
1098			0.016389886745811315,
1099			-0.004090989452162578,
1100			0.00036818904090536926,
1101			-0.0031394146217732367,
1102			-7.571300829706796e-6,
1103			3.0045718718618155e-6,
1104			2.1329796886727743e-6,
1105			9.259701025627789e-8,
1106			2.2291214755992877e-9,
1107			-2.3017207713252894e-9,
1108			6.807967994979358e-11,
1109			2.1677299575405587e-12,
1110			-3.07282771050034e-13,
1111			0.0,
1112		];
1113		let n = diag.len();
1114		let params = SvdParams {
1115			recursion_threshold: 8,
1116			qr_ratio_threshold: 1.0,
1117			..auto!(f64)
1118		};
1119		let mut d = ColRef::from_slice(&diag).to_owned();
1120		let mut s = ColRef::from_slice(&subdiag).to_owned();
1121		compute_bidiag_real_svd(
1122			d.as_mut().try_as_col_major_mut().unwrap(),
1123			s.as_mut().try_as_col_major_mut().unwrap(),
1124			None,
1125			None,
1126			params,
1127			Par::Seq,
1128			MemStack::new(&mut MemBuffer::new(bidiag_real_svd_scratch::<f64>(
1129				n,
1130				false,
1131				false,
1132				Par::Seq,
1133				params,
1134			))),
1135		)
1136		.unwrap();
1137		assert!(d[n - 1] != 0.0);
1138	}
1139	#[test]
1140	fn test_pinv() {
1141		let rng = &mut StdRng::seed_from_u64(0);
1142		let m = 6;
1143		let n = 36;
1144		let mat = CwiseMatDistribution {
1145			nrows: m,
1146			ncols: n,
1147			dist: StandardNormal,
1148		}
1149		.rand::<Mat<f64>>(rng);
1150		let pinv = mat.svd().unwrap().pseudoinverse();
1151		let err = &mat * &pinv - Mat::<f64>::identity(m, m);
1152		assert!(err.norm_max() < 1e-10);
1153	}
1154	#[test]
1155	fn test_pinv2() {
1156		let rng = &mut StdRng::seed_from_u64(0);
1157		let m = 36;
1158		let n = 6;
1159		let mat = CwiseMatDistribution {
1160			nrows: m,
1161			ncols: n,
1162			dist: StandardNormal,
1163		}
1164		.rand::<Mat<f64>>(rng);
1165		let pinv = mat.svd().unwrap().pseudoinverse();
1166		let err = &pinv * &mat - Mat::<f64>::identity(n, n);
1167		assert!(err.norm_max() < 1e-10);
1168	}
1169}
faer/linalg/svd/mod.rs

faer/linalg/svd/
mod.rs
