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: ParShape,
70		params: SvdParams,
71	) -> StackReq,
72	params: SvdParams,
73	par: ParShape,
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: ParShape,
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: ParShape,
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			// a subnormal component scale has an infinite reciprocal; rescale it the
220			// way reference `dlarfg` does instead of dividing by a zero product
221			let ref shift = if norm1 >= min_positive::<T::Real>() {
222				from_f64::<T::Real>(1.0)
223			} else {
224				from_f64::<T::Real>(4503599627370496.0)
225			};
226			let y = x.mul_real(shift).mul_real((&norm1 * shift).recip());
227			y.mul_real(y.abs().recip())
228		}
229	};
230	let magnitude = |x: &T| {
231		// `abs` squares components scaled by `sqrt(min_positive)` or
232		// `sqrt(max_positive)`, which flushes a magnitude below `min_positive`
233		// to zero and overflows above about `max_positive / 2`. Outside
234		// `[min_positive, 1 / min_positive]`, rescale by that exact power of
235		// two first; other entries keep the historical `abs` bit for bit.
236		let sml = min_positive::<T::Real>();
237		let big = sml.recip();
238		let m = x.real().abs().fmax(x.imag().abs());
239		if m < sml && m != zero() {
240			x.mul_real(&big).abs() * sml
241		} else if m > big && m.is_finite() {
242			x.mul_real(&sml).abs() * big
243		} else {
244			x.abs()
245		}
246	};
247	let mut col_normalized = normalized(diag[0].conj());
248	col_mul[0] = col_normalized.copy();
249	diag_real[0] = magnitude(&diag[0]);
250	subdiag_real[n - 1] = zero();
251	for i in 1..n {
252		let row_normalized =
253			normalized((&subdiag[i - 1] * col_normalized).conj());
254		subdiag_real[i - 1] = magnitude(&subdiag[i - 1]);
255		row_mul[i - 1] = row_normalized.conj();
256		col_normalized = normalized((&diag[i] * row_normalized).conj());
257		diag_real[i] = magnitude(&diag[i]);
258		col_mul[i] = col_normalized.copy();
259	}
260	compute_bidiag_real_svd(
261		diag_real.rb_mut(),
262		subdiag_real.rb_mut(),
263		u_real.rb_mut(),
264		v_real.rb_mut(),
265		params,
266		par,
267		stack,
268	)?;
269	for i in 0..n {
270		diag[i] = diag_real[i].to_cplx();
271	}
272	let u_real = u_real.rb();
273	let v_real = v_real.rb();
274	if let (Some(mut u), Some(u_real)) = (u.rb_mut(), u_real) {
275		z!(u.rb_mut().row_mut(0), u_real.row(0))
276			.for_each(|uz!(u, r)| *u = r.to_cplx());
277		z!(u.rb_mut().row_mut(n), u_real.row(n))
278			.for_each(|uz!(u, r)| *u = r.to_cplx());
279		for j in 0..u.ncols() {
280			let mut u = u.rb_mut().col_mut(j).subrows_mut(1, n - 1);
281			let u_real = u_real.rb().col(j).subrows(1, n - 1);
282			z!(u.rb_mut(), u_real, row_mul.rb())
283				.for_each(|uz!(u, re, f)| *u = f.mul_real(re));
284		}
285	}
286	if let (Some(mut v), Some(v_real)) = (v.rb_mut(), v_real) {
287		for j in 0..v.ncols() {
288			let mut v = v.rb_mut().col_mut(j);
289			let v_real = v_real.rb().col(j);
290			z!(v.rb_mut(), v_real, col_mul.rb())
291				.for_each(|uz!(v, re, f)| *v = f.mul_real(re));
292		}
293	}
294	Ok(())
295}
296fn compute_bidiag_real_svd<T: RealField>(
297	mut diag: ColMut<'_, T, usize, ContiguousFwd>,
298	mut subdiag: ColMut<'_, T, usize, ContiguousFwd>,
299	mut u: Option<MatMut<'_, T, usize, usize>>,
300	mut v: Option<MatMut<'_, T, usize, usize>>,
301	params: SvdParams,
302	par: Par,
303	stack: &mut MemStack,
304) -> Result<(), SvdError> {
305	let n = diag.nrows();
306	for i in 0..n {
307		if !(diag[i].is_finite() && subdiag[i].is_finite()) {
308			return Err(SvdError::NoConvergence);
309		}
310	}
311	if n < params.recursion_threshold {
312		if let Some(mut u) = u.rb_mut() {
313			u.fill(zero());
314			u.diagonal_mut().fill(one());
315		}
316		if let Some(mut v) = v.rb_mut() {
317			v.fill(zero());
318			v.diagonal_mut().fill(one());
319		}
320		bidiag_svd::qr_algorithm(
321			diag.rb_mut(),
322			subdiag.rb_mut(),
323			u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
324			v.rb_mut(),
325		)?;
326		return Ok(());
327	} else {
328		let (mut u2, stack) = unsafe {
329			temp_mat_uninit::<T::Real, _, _>(
330				2,
331				if u.is_some() { 0 } else { n + 1 },
332				stack,
333			)
334		};
335		bidiag_svd::divide_and_conquer(
336			diag.as_row_shape_mut(n),
337			subdiag.as_row_shape_mut(n),
338			match u {
339				Some(u) => bidiag_svd::MatU::Full(u),
340				None => bidiag_svd::MatU::TwoRowsStorage(u2.as_mat_mut()),
341			},
342			v.map(|m| m.as_shape_mut(n, n)),
343			par,
344			stack,
345			params.recursion_threshold,
346		)
347	}
348}
349/// Workspace for [`bidiagonal_svd`], including the padded
350/// lower-bidiagonal storage used by the existing divide-and-conquer
351/// implementation.
352///
353/// # Examples
354/// ```
355/// use faer::Par;
356/// use faer::linalg::svd::bidiagonal_svd_scratch;
357/// let _buffer = dyn_stack::MemBuffer::new(bidiagonal_svd_scratch::<f64>(
358/// 	3,
359/// 	faer::ParShape::Seq,
360/// 	Default::default(),
361/// ));
362/// ```
363pub fn bidiagonal_svd_scratch<T: ComplexField>(
364	dim: usize,
365	par: ParShape,
366	params: Spec<SvdParams, T>,
367) -> StackReq {
368	let n = dim;
369	if n == 0 {
370		return StackReq::EMPTY;
371	}
372	StackReq::all_of(&[
373		temp_mat_scratch::<T>(n, 1).array(2),
374		temp_mat_scratch::<T>(n + 1, n + 1),
375		temp_mat_scratch::<T>(n, n),
376		bidiag_cplx_svd_scratch::<T>(n, true, true, par, *params),
377	])
378}
379/// Computes a square bidiagonal SVD `B = U * diag(s) * V.adjoint()`.
380/// `diag` is overwritten with nonnegative, descending singular values (with
381/// zero imaginary part). `offdiag` has length `n.saturating_sub(1)` and is
382/// unchanged. `upper` selects the superdiagonal rather than the subdiagonal.
383/// Both `u` and `v` are `n` by `n`; the output is **V**, not its adjoint.
384/// Real and complex entries are supported via the existing phase reduction.
385///
386/// # Panics
387/// Panics on incompatible dimensions or insufficient scratch.
388///
389/// # Errors
390/// Returns [`SvdError::NoConvergence`] for nonfinite input, iteration
391/// failure, or a result that would not be finite (for example when a complex
392/// entry's modulus overflows). Finite input never yields `Ok` with nonfinite
393/// output. Outputs are unspecified on error; `offdiag` is always unchanged.
394///
395/// # Examples
396/// ```
397/// use dyn_stack::{MemBuffer, MemStack};
398/// use faer::linalg::svd::*;
399/// use faer::{Mat, Par};
400/// let mut d = [3.0f64, -2.0];
401/// let (mut u, mut v) = (Mat::zeros(2, 2), Mat::zeros(2, 2));
402/// let mut mem = MemBuffer::new(bidiagonal_svd_scratch::<f64>(
403/// 	2,
404/// 	faer::ParShape::Seq,
405/// 	Default::default(),
406/// ));
407/// bidiagonal_svd(
408/// 	&mut d,
409/// 	&[0.0],
410/// 	true,
411/// 	u.as_mut(),
412/// 	v.as_mut(),
413/// 	Par::Seq,
414/// 	MemStack::new(&mut mem),
415/// 	Default::default(),
416/// )
417/// .unwrap();
418/// assert_eq!(d, [3.0, 2.0]);
419/// ```
420pub fn bidiagonal_svd<T: ComplexField>(
421	diag: &mut [T],
422	offdiag: &[T],
423	upper: bool,
424	mut u: MatMut<'_, T>,
425	mut v: MatMut<'_, T>,
426	par: Par,
427	stack: &mut MemStack,
428	params: Spec<SvdParams, T>,
429) -> Result<(), SvdError> {
430	let n = diag.len();
431	assert!(offdiag.len() == n.saturating_sub(1));
432	assert!(
433		u.nrows() == n && u.ncols() == n && v.nrows() == n && v.ncols() == n
434	);
435	if !diag.iter().chain(offdiag).all(|x| x.is_finite()) {
436		return Err(SvdError::NoConvergence);
437	}
438	if n == 0 {
439		return Ok(());
440	}
441	// SAFETY: Every entry is written below before the kernel reads the
442	// diagonal.
443	let (mut d, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
444	let (mut e, stack) = temp_mat_zeroed::<T, _, _>(n, 1, stack);
445	// SAFETY: The phase-reduction kernel fully writes each factor on success;
446	// its error path returns before either factor is read here.
447	let (mut left, stack) =
448		unsafe { temp_mat_uninit::<T, _, _>(n + 1, n + 1, stack) };
449	// SAFETY: Same full-overwrite contract as the left factor above.
450	let (mut right, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, n, stack) };
451	let mut d = d.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
452	let mut e = e.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
453	for i in 0..n {
454		d[i] = if upper {
455			diag[i].conj()
456		} else {
457			diag[i].copy()
458		};
459		if i + 1 < n {
460			e[i] = if upper {
461				offdiag[i].conj()
462			} else {
463				offdiag[i].copy()
464			};
465		}
466	}
467	compute_bidiag_cplx_svd(
468		d.rb_mut(),
469		e.rb_mut(),
470		Some(left.as_mat_mut()),
471		Some(right.as_mat_mut()),
472		*params,
473		par,
474		stack,
475	)?;
476	// Backstop: finite input never reports success with nonfinite output.
477	if !(d.rb().is_all_finite()
478		&& left.as_mat_ref().submatrix(0, 0, n, n).is_all_finite()
479		&& right.as_mat_ref().is_all_finite())
480	{
481		return Err(SvdError::NoConvergence);
482	}
483	for i in 0..n {
484		diag[i] = d[i].copy();
485	}
486	if upper {
487		u.copy_from(right.as_mat_ref());
488		v.copy_from(left.as_mat_ref().submatrix(0, 0, n, n));
489	} else {
490		u.copy_from(left.as_mat_ref().submatrix(0, 0, n, n));
491		v.copy_from(right.as_mat_ref());
492	}
493	Ok(())
494}
495/// bidiag -> divide conquer svd / qr algo
496fn svd_imp<T: ComplexField>(
497	matrix: MatRef<'_, T>,
498	s: ColMut<'_, T>,
499	u: Option<MatMut<'_, T>>,
500	v: Option<MatMut<'_, T>>,
501	bidiag_svd: fn(
502		diag: ColMut<'_, T, usize, ContiguousFwd>,
503		subdiag: ColMut<'_, T, usize, ContiguousFwd>,
504		u: Option<MatMut<'_, T, usize, usize>>,
505		v: Option<MatMut<'_, T, usize, usize>>,
506		params: SvdParams,
507		par: Par,
508		stack: &mut MemStack,
509	) -> Result<(), SvdError>,
510	mut par: Par,
511	stack: &mut MemStack,
512	params: SvdParams,
513) -> Result<(), SvdError> {
514	assert!(matrix.nrows() >= matrix.ncols());
515	let m = matrix.nrows();
516	let n = matrix.ncols();
517	let bs = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
518	let (mut bid, stack) = unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
519	let mut bid = bid.as_mat_mut();
520	let (mut Hl, stack) = unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
521	let (mut Hr, stack) =
522		unsafe { temp_mat_uninit::<T, _, _>(bs, n - 1, stack) };
523	let mut Hl = Hl.as_mat_mut();
524	let mut Hr = Hr.as_mat_mut();
525	bid.copy_from(matrix);
526	bidiag::bidiag_in_place(
527		bid.rb_mut(),
528		Hl.rb_mut(),
529		Hr.rb_mut(),
530		par.reborrow(),
531		stack,
532		params.bidiag.into(),
533	);
534	let (mut diag, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
535	let (mut subdiag, stack) =
536		unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
537	let mut diag = diag.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
538	let mut subdiag = subdiag
539		.as_mat_mut()
540		.col_mut(0)
541		.try_as_col_major_mut()
542		.unwrap();
543	let (mut ub, stack) = unsafe {
544		temp_mat_uninit::<T, _, _>(
545			n + 1,
546			if v.is_some() { n + 1 } else { 0 },
547			stack,
548		)
549	};
550	let (mut vb, stack) = unsafe {
551		temp_mat_uninit::<T, _, _>(n, if u.is_some() { n } else { 0 }, stack)
552	};
553	let mut ub = ub.as_mat_mut();
554	let mut vb = vb.as_mat_mut();
555	for i in 0..n {
556		diag[i] = bid[(i, i)].conj();
557		if i + 1 < n {
558			subdiag[i] = bid[(i, i + 1)].conj();
559		} else {
560			subdiag[i] = zero();
561		}
562	}
563	bidiag_svd(
564		diag.rb_mut(),
565		subdiag.rb_mut(),
566		v.rb().map(|_| ub.rb_mut()),
567		u.rb().map(|_| vb.rb_mut()),
568		params,
569		par.reborrow(),
570		stack,
571	)?;
572	{ s }.copy_from(diag);
573	if let Some(mut u) = u {
574		let ncols = u.ncols();
575		u.rb_mut().submatrix_mut(0, 0, n, n).copy_from(vb.rb());
576		u.rb_mut().submatrix_mut(n, 0, m - n, ncols).fill(zero());
577		u.rb_mut().submatrix_mut(0, n, n, ncols - n).fill(zero());
578		u.rb_mut()
579			.submatrix_mut(n, n, ncols - n, ncols - n)
580			.diagonal_mut()
581			.fill(one());
582		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(bid.rb(), Hl.rb(), Conj::No, u, par.reborrow(), stack);
583	}
584	if let Some(mut v) = v {
585		v.copy_from(ub.rb().submatrix(0, 0, n, n));
586		for j in 1..n {
587			for i in 0..j {
588				bid[(j, i)] = bid[(i, j)].copy();
589			}
590		}
591		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
592			bid.rb().submatrix(1, 0, n - 1, n - 1),
593			Hr.rb(),
594			Conj::Yes,
595			v.subrows_mut(1, n - 1),
596			par,
597			stack,
598		);
599	}
600	Ok(())
601}
602fn compute_squareish_svd<T: ComplexField>(
603	matrix: MatRef<'_, T>,
604	s: ColMut<'_, T>,
605	u: Option<MatMut<'_, T>>,
606	v: Option<MatMut<'_, T>>,
607	par: Par,
608	stack: &mut MemStack,
609	params: SvdParams,
610) -> Result<(), SvdError> {
611	if const { T::IS_REAL } {
612		svd_imp::<T::Real>(
613			unsafe { core::mem::transmute(matrix) },
614			unsafe { core::mem::transmute(s) },
615			unsafe { core::mem::transmute(u) },
616			unsafe { core::mem::transmute(v) },
617			compute_bidiag_real_svd::<T::Real>,
618			par,
619			stack,
620			params,
621		)
622	} else {
623		svd_imp::<T>(
624			matrix,
625			s,
626			u,
627			v,
628			compute_bidiag_cplx_svd::<T>,
629			par,
630			stack,
631			params,
632		)
633	}
634}
635/// computes the layout of the workspace required to compute a matrix's svd
636pub fn svd_scratch<T: ComplexField>(
637	nrows: usize,
638	ncols: usize,
639	compute_u: ComputeSvdVectors,
640	compute_v: ComputeSvdVectors,
641	par: ParShape,
642	params: Spec<SvdParams, T>,
643) -> StackReq {
644	let params = params.config;
645	let mut m = nrows;
646	let mut n = ncols;
647	let mut compute_u = compute_u;
648	let mut compute_v = compute_v;
649	if n > m {
650		core::mem::swap(&mut m, &mut n);
651		core::mem::swap(&mut compute_u, &mut compute_v);
652	}
653	if n == 0 {
654		return StackReq::EMPTY;
655	}
656	let bidiag_svd_scratch = if const { T::IS_REAL } {
657		bidiag_real_svd_scratch::<T::Real>
658	} else {
659		bidiag_cplx_svd_scratch::<T>
660	};
661	if m as f64 / n as f64 <= params.qr_ratio_threshold {
662		svd_imp_scratch::<T>(
663			m,
664			n,
665			compute_u,
666			compute_v,
667			bidiag_svd_scratch,
668			params,
669			par,
670		)
671	} else {
672		let bs =
673			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
674		StackReq::all_of(&[
675			temp_mat_scratch::<T>(m, n),
676			temp_mat_scratch::<T>(bs, n),
677			StackReq::any_of(&[
678				StackReq::all_of(&[
679					temp_mat_scratch::<T>(n, n),
680					svd_imp_scratch::<T>(n, n, compute_u, compute_v, bidiag_svd_scratch, params, par),
681				]),
682				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
683					m,
684					bs,
685					match compute_u {
686						ComputeSvdVectors::No => 0,
687						ComputeSvdVectors::Thin => n,
688						ComputeSvdVectors::Full => m,
689					},
690				),
691			]),
692		])
693	}
694}
695/// computes the svd of $A$, with the singular vectors being omitted, thin or
696/// full
697///
698/// the singular are stored in $S$, and the singular vectors in $U$ and $V$ such
699/// that the singular values are sorted in nonincreasing order
700pub fn svd<T: ComplexField>(
701	A: MatRef<'_, T>,
702	s: DiagMut<'_, T>,
703	u: Option<MatMut<'_, T>>,
704	v: Option<MatMut<'_, T>>,
705	mut par: Par,
706	stack: &mut MemStack,
707	params: Spec<SvdParams, T>,
708) -> Result<(), SvdError> {
709	let params = params.config;
710	let (m, n) = A.shape();
711	let size = Ord::min(m, n);
712	assert!(s.dim() == size);
713	let s = s.column_vector_mut();
714	if let Some(u) = u.rb() {
715		assert!(all(
716			u.nrows() == A.nrows(),
717			any(u.ncols() == A.nrows(), u.ncols() == size),
718		));
719	}
720	if let Some(v) = v.rb() {
721		assert!(all(
722			v.nrows() == A.ncols(),
723			any(v.ncols() == A.ncols(), v.ncols() == size),
724		));
725	}
726	#[cfg(feature = "perf-warn")]
727	match (u.rb(), v.rb()) {
728		(Some(matrix), _) | (_, Some(matrix)) => {
729			if matrix.row_stride().unsigned_abs() != 1
730				&& crate::__perf_warn!(QR_WARN)
731			{
732				if matrix.col_stride().unsigned_abs() == 1 {
733					log::warn!(
734						target : "faer_perf",
735						"SVD prefers column-major singular vector matrices. Found row-major matrix."
736					);
737				} else {
738					log::warn!(
739						target : "faer_perf",
740						"SVD prefers column-major singular vector matrices. Found matrix with generic strides."
741					);
742				}
743			}
744		},
745		_ => {},
746	}
747	let mut u = u;
748	let mut v = v;
749	let mut matrix = A;
750	let do_transpose = n > m;
751	if do_transpose {
752		matrix = matrix.transpose();
753		core::mem::swap(&mut u, &mut v)
754	}
755	let (m, n) = matrix.shape();
756	if n == 0 {
757		if let Some(mut u) = u {
758			u.fill(zero());
759			u.rb_mut().diagonal_mut().fill(one());
760		}
761		return Ok(());
762	}
763	if m as f64 / n as f64 <= params.qr_ratio_threshold {
764		compute_squareish_svd(
765			matrix,
766			s,
767			u.rb_mut(),
768			v.rb_mut(),
769			par,
770			stack,
771			params,
772		)?;
773	} else {
774		let bs =
775			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
776		let (mut qr, stack) =
777			unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
778		let mut qr = qr.as_mat_mut();
779		let (mut householder, stack) =
780			unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
781		let mut householder = householder.as_mat_mut();
782		{
783			qr.copy_from(matrix.rb());
784			linalg::qr::no_pivoting::factor::qr_in_place(
785				qr.rb_mut(),
786				householder.rb_mut(),
787				par.reborrow(),
788				stack,
789				params.qr.into(),
790			);
791		}
792		{
793			let (mut r, stack) =
794				unsafe { temp_mat_uninit::<T, _, _>(n, n, stack) };
795			let mut r = r.as_mat_mut();
796			z!(r.rb_mut()).for_each_triangular_lower(
797				linalg::zip::Diag::Skip,
798				|uz!(dst)| *dst = zero(),
799			);
800			z!(r.rb_mut(), qr.rb().submatrix(0, 0, n, n))
801				.for_each_triangular_upper(
802					linalg::zip::Diag::Include,
803					|uz!(dst, src)| *dst = src.copy(),
804				);
805			compute_squareish_svd(
806				r.rb(),
807				s,
808				u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
809				v.rb_mut(),
810				par.reborrow(),
811				stack,
812				params,
813			)?;
814		}
815		if let Some(mut u) = u.rb_mut() {
816			u.rb_mut().subrows_mut(n, m - n).fill(zero());
817			if u.ncols() == m {
818				// the extra columns are `Q * [0; I]`, so their top block has
819				// to start at zero instead of at the caller's old entries
820				u.rb_mut().submatrix_mut(0, n, n, m - n).fill(zero());
821				u.rb_mut()
822					.submatrix_mut(n, n, m - n, m - n)
823					.diagonal_mut()
824					.fill(one());
825			}
826			linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
827				qr.rb(),
828				householder.rb(),
829				Conj::No,
830				u.rb_mut(),
831				par,
832				stack,
833			);
834		}
835	}
836	if do_transpose {
837		if let Some(u) = u.rb_mut() {
838			z!(u).for_each(|uz!(u)| *u = u.conj())
839		}
840		if let Some(v) = v.rb_mut() {
841			z!(v).for_each(|uz!(v)| *v = v.conj())
842		}
843	}
844	Ok(())
845}
846/// computes the layout of the workspace required to compute a matrix's
847/// pseudoinverse, given the svd
848pub fn pseudoinverse_from_svd_scratch<T: ComplexField>(
849	nrows: usize,
850	ncols: usize,
851	par: ParShape,
852) -> StackReq {
853	_ = par;
854	let size = Ord::min(nrows, ncols);
855	StackReq::all_of(&[
856		temp_mat_scratch::<T>(nrows, size),
857		temp_mat_scratch::<T>(ncols, size),
858	])
859}
860/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
861/// $U$ and $V$
862pub fn pseudoinverse_from_svd<T: ComplexField>(
863	pinv: MatMut<'_, T>,
864	s: DiagRef<'_, T>,
865	u: MatRef<'_, T>,
866	v: MatRef<'_, T>,
867	par: Par,
868	stack: &mut MemStack,
869) {
870	pseudoinverse_from_svd_with_tolerance(
871		pinv,
872		s,
873		u,
874		v,
875		zero(),
876		eps::<T::Real>()
877			* from_f64::<T::Real>(Ord::max(u.nrows(), v.nrows()) as f64),
878		par,
879		stack,
880	);
881}
882/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
883/// $U$ and $V$, and tolerance parameters for determining zero singular values
884pub fn pseudoinverse_from_svd_with_tolerance<T: ComplexField>(
885	pinv: MatMut<'_, T>,
886	s: DiagRef<'_, T>,
887	u: MatRef<'_, T>,
888	v: MatRef<'_, T>,
889	abs_tol: T::Real,
890	rel_tol: T::Real,
891	par: Par,
892	stack: &mut MemStack,
893) {
894	let mut pinv = pinv;
895	let m = u.nrows();
896	let n = v.nrows();
897	let size = Ord::min(m, n);
898	assert!(all(
899		u.nrows() == m,
900		v.nrows() == n,
901		u.ncols() >= size,
902		v.ncols() >= size,
903		s.dim() >= size
904	));
905	let s = s.column_vector();
906	let u = u.get(.., ..size);
907	let v = v.get(.., ..size);
908	let smax = s.norm_max();
909	let tol = abs_tol.fmax(rel_tol * smax);
910	let (mut u_trunc, stack) =
911		unsafe { temp_mat_uninit::<T, _, _>(m, size, stack) };
912	let (mut vp_trunc, _) =
913		unsafe { temp_mat_uninit::<T, _, _>(n, size, stack) };
914	let mut u_trunc = u_trunc.as_mat_mut();
915	let mut vp_trunc = vp_trunc.as_mat_mut();
916	let mut len = 0;
917	for j in 0..size {
918		let x = s[j].absmax();
919		if x > tol {
920			let ref p = s[j].real().recip();
921			u_trunc.rb_mut().col_mut(len).copy_from(u.col(j));
922			z!(vp_trunc.rb_mut().col_mut(len), v.col(j))
923				.for_each(|uz!(dst, src)| *dst = src.mul_real(p));
924			len += 1;
925		}
926	}
927	let u_trunc = u_trunc.get(.., ..len);
928	let vp_trunc = vp_trunc.get(.., ..len);
929	linalg::matmul::matmul(
930		pinv.rb_mut(),
931		Accum::Replace,
932		vp_trunc,
933		u_trunc.adjoint(),
934		one(),
935		par,
936	);
937}
938#[cfg(test)]
939mod tests {
940	// Instantiate the same cases and checks at each native precision.
941	mod lapack_svd_f64 {
942		type Real = f64;
943		type Scalar = Real;
944		const GUARD: Scalar = -12345.0;
945		const ROTATION_PHASE: Scalar = 1.0;
946		fn diagonal_phase(i: usize) -> Scalar {
947			if i % 2 == 0 { 1.0 } else { -1.0 }
948		}
949		use super::super::{
950			bidiag_real_svd_scratch as bidiagonal_scratch,
951			compute_bidiag_real_svd as compute_bidiagonal,
952		};
953		mod lapack {
954			include!("tests/lapack.rs");
955		}
956		mod svd_views {
957			include!("tests/svd_views.rs");
958		}
959	}
960	mod lapack_svd_f32 {
961		type Real = f32;
962		type Scalar = Real;
963		const GUARD: Scalar = -12345.0;
964		const ROTATION_PHASE: Scalar = 1.0;
965		fn diagonal_phase(i: usize) -> Scalar {
966			if i % 2 == 0 { 1.0 } else { -1.0 }
967		}
968		use super::super::{
969			bidiag_real_svd_scratch as bidiagonal_scratch,
970			compute_bidiag_real_svd as compute_bidiagonal,
971		};
972		mod lapack {
973			include!("tests/lapack.rs");
974		}
975		mod svd_views {
976			include!("tests/svd_views.rs");
977		}
978		#[test]
979		fn f32_fixtures_use_single_precision_limits() {
980			let (_, large) = lapack::diagonal_case(2, 2, 6);
981			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
982			let (_, small) = lapack::diagonal_case(2, 2, 7);
983			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
984			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
985			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
986			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
987			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
988		}
989	}
990	mod lapack_svd_c64 {
991		type Real = f64;
992		type Scalar = crate::c64;
993		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
994		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
995		fn diagonal_phase(i: usize) -> Scalar {
996			match i % 4 {
997				0 => Scalar::new(0.0, 1.0),
998				1 => Scalar::new(0.6, 0.8),
999				2 => Scalar::new(-0.8, 0.6),
1000				_ => Scalar::new(-1.0, 0.0),
1001			}
1002		}
1003		use super::super::{
1004			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1005			compute_bidiag_cplx_svd as compute_bidiagonal,
1006		};
1007		mod lapack {
1008			include!("tests/lapack.rs");
1009		}
1010		mod svd_views {
1011			include!("tests/svd_views.rs");
1012		}
1013		mod complex_tests {
1014			include!("tests/complex.rs");
1015		}
1016	}
1017	mod lapack_svd_c32 {
1018		type Real = f32;
1019		type Scalar = crate::c32;
1020		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
1021		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
1022		fn diagonal_phase(i: usize) -> Scalar {
1023			match i % 4 {
1024				0 => Scalar::new(0.0, 1.0),
1025				1 => Scalar::new(0.6, 0.8),
1026				2 => Scalar::new(-0.8, 0.6),
1027				_ => Scalar::new(-1.0, 0.0),
1028			}
1029		}
1030		use super::super::{
1031			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1032			compute_bidiag_cplx_svd as compute_bidiagonal,
1033		};
1034		mod lapack {
1035			include!("tests/lapack.rs");
1036		}
1037		mod svd_views {
1038			include!("tests/svd_views.rs");
1039		}
1040		mod complex_tests {
1041			include!("tests/complex.rs");
1042		}
1043		#[test]
1044		fn c32_fixtures_use_single_precision_limits() {
1045			let (_, large) = lapack::diagonal_case(2, 2, 6);
1046			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
1047			let (_, small) = lapack::diagonal_case(2, 2, 7);
1048			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
1049			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
1050			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
1051			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
1052			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
1053		}
1054	}
1055	use super::*;
1056	use crate::assert;
1057	use crate::stats::prelude::*;
1058	use crate::utils::approx::*;
1059	use dyn_stack::MemBuffer;
1060	#[track_caller]
1061	fn test_svd<T: ComplexField>(mat: MatRef<'_, T>) {
1062		let (m, n) = mat.shape();
1063		let params = Spec::new(SvdParams {
1064			recursion_threshold: 8,
1065			qr_ratio_threshold: 1.0,
1066			..auto!(T)
1067		});
1068		let approx_eq = CwiseMat(
1069			ApproxEq::<T::Real>::eps()
1070				* from_f64::<T::Real>(8.0 * Ord::max(m, n) as f64).sqrt(),
1071		);
1072		{
1073			let mut s = Mat::zeros(m, n);
1074			let mut u = Mat::zeros(m, m);
1075			let mut v = Mat::zeros(n, n);
1076			svd(
1077				mat.as_ref(),
1078				s.as_mut().diagonal_mut(),
1079				Some(u.as_mut()),
1080				Some(v.as_mut()),
1081				Par::Seq,
1082				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1083					m,
1084					n,
1085					ComputeSvdVectors::Full,
1086					ComputeSvdVectors::Full,
1087					ParShape::Seq,
1088					params,
1089				))),
1090				params,
1091			)
1092			.unwrap();
1093			let reconstructed = &u * &s * v.adjoint();
1094			assert!(reconstructed ~ mat);
1095		}
1096		let size = Ord::min(m, n);
1097		let mut s = Mat::zeros(size, size);
1098		let mut u = Mat::zeros(m, size);
1099		let mut v = Mat::zeros(n, size);
1100		{
1101			svd(
1102				mat.as_ref(),
1103				s.as_mut().diagonal_mut(),
1104				Some(u.as_mut()),
1105				Some(v.as_mut()),
1106				Par::Seq,
1107				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1108					m,
1109					n,
1110					ComputeSvdVectors::Thin,
1111					ComputeSvdVectors::Thin,
1112					ParShape::Seq,
1113					params,
1114				))),
1115				params,
1116			)
1117			.unwrap();
1118			let reconstructed = &u * &s * v.adjoint();
1119			assert!(reconstructed ~ mat);
1120		}
1121		{
1122			let mut s2 = Mat::zeros(size, size);
1123			let mut u2 = Mat::zeros(m, size);
1124			svd(
1125				mat.as_ref(),
1126				s2.as_mut().diagonal_mut(),
1127				Some(u2.as_mut()),
1128				None,
1129				Par::Seq,
1130				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1131					m,
1132					n,
1133					ComputeSvdVectors::Thin,
1134					ComputeSvdVectors::No,
1135					ParShape::Seq,
1136					params,
1137				))),
1138				params,
1139			)
1140			.unwrap();
1141			assert!(s2 ~ s);
1142			assert!(u2 ~ u);
1143		}
1144		{
1145			let mut s2 = Mat::zeros(size, size);
1146			let mut v2 = Mat::zeros(n, size);
1147			svd(
1148				mat.as_ref(),
1149				s2.as_mut().diagonal_mut(),
1150				None,
1151				Some(v2.as_mut()),
1152				Par::Seq,
1153				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1154					m,
1155					n,
1156					ComputeSvdVectors::No,
1157					ComputeSvdVectors::Thin,
1158					ParShape::Seq,
1159					params,
1160				))),
1161				params,
1162			)
1163			.unwrap();
1164			assert!(s2 ~ s);
1165			assert!(v2 ~ v);
1166		}
1167		{
1168			let mut s2 = Mat::zeros(size, size);
1169			svd(
1170				mat.as_ref(),
1171				s2.as_mut().diagonal_mut(),
1172				None,
1173				None,
1174				Par::Seq,
1175				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1176					m,
1177					n,
1178					ComputeSvdVectors::No,
1179					ComputeSvdVectors::No,
1180					ParShape::Seq,
1181					params,
1182				))),
1183				params,
1184			)
1185			.unwrap();
1186			assert!(s2 ~ s);
1187		}
1188	}
1189	fn check_tiny_split<T: ComplexField>(a: MatRef<'_, T>) {
1190		let (m, n) = a.shape();
1191		let size = Ord::min(m, n);
1192		let compute = |threshold, want_u, want_v| {
1193			let params = Spec::new(SvdParams {
1194				recursion_threshold: threshold,
1195				..auto!(T)
1196			});
1197			let mut s = Col::zeros(size);
1198			let mut u = Mat::zeros(m, if want_u { size } else { 0 });
1199			let mut v = Mat::zeros(n, if want_v { size } else { 0 });
1200			let vectors = |want| {
1201				if want {
1202					ComputeSvdVectors::Thin
1203				} else {
1204					ComputeSvdVectors::No
1205				}
1206			};
1207			svd(
1208				a,
1209				s.as_mut().as_diagonal_mut(),
1210				if want_u { Some(u.as_mut()) } else { None },
1211				if want_v { Some(v.as_mut()) } else { None },
1212				Par::Seq,
1213				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1214					m,
1215					n,
1216					vectors(want_u),
1217					vectors(want_v),
1218					ParShape::Seq,
1219					params,
1220				))),
1221				params,
1222			)
1223			.unwrap();
1224			(s, u, v)
1225		};
1226		let (reference, _, _) = compute(usize::MAX, true, true);
1227		for (want_u, want_v) in
1228			[(true, true), (false, false), (true, false), (false, true)]
1229		{
1230			let (s, u, v) = compute(128, want_u, want_v);
1231			assert!((&s - &reference).norm_max() < from_f64::<T::Real>(1e-13));
1232			if want_u {
1233				assert!(
1234					(u.adjoint() * &u - Mat::<T>::identity(size, size))
1235						.norm_l2() < from_f64::<T::Real>(1e-13)
1236				);
1237			}
1238			if want_v {
1239				assert!(
1240					(v.adjoint() * &v - Mat::<T>::identity(size, size))
1241						.norm_l2() < from_f64::<T::Real>(1e-13)
1242				);
1243			}
1244			if want_u && want_v {
1245				let back = &u * s.as_diagonal() * v.adjoint();
1246				assert!(
1247					(&back - a).norm_l2() / a.norm_l2()
1248						< from_f64::<T::Real>(1e-13)
1249				);
1250			}
1251		}
1252	}
1253	#[test]
1254	fn test_tiny_diagonal_at_split() {
1255		for n in [128, 256] {
1256			for tiny in [0.0, 1e-15, 1e-13, 1e-10, 1e-6] {
1257				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1258					if (i == j && i < 64) || (j == i + 1 && i < 64) {
1259						1.0
1260					} else if i == 64 && j == 64 {
1261						tiny
1262					} else {
1263						0.0
1264					}
1265				});
1266				check_tiny_split(a.as_ref());
1267				let complex =
1268					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1269				check_tiny_split(complex.as_ref());
1270			}
1271		}
1272	}
1273	/// A zero block at a divide-and-conquer split must still hand its parent merge the
1274	/// identity boundary entry that carries the nonzero part of the spectrum, whatever
1275	/// singular vectors the caller asked for (#33).
1276	#[test]
1277	fn test_zero_block_at_split() {
1278		fn check<T: ComplexField>(
1279			n: usize,
1280			threshold: usize,
1281			want: (bool, bool),
1282			tol: f64,
1283		) {
1284			let (want_u, want_v) = want;
1285			let a = Mat::<T>::from_fn(n, n, |i, j| {
1286				if i == n / 2 && j == n / 2 {
1287					one()
1288				} else {
1289					zero()
1290				}
1291			});
1292			let params = Spec::new(SvdParams {
1293				recursion_threshold: threshold,
1294				..auto!(T)
1295			});
1296			let vectors = |want| {
1297				if want {
1298					ComputeSvdVectors::Full
1299				} else {
1300					ComputeSvdVectors::No
1301				}
1302			};
1303			let mut s = Col::<T>::zeros(n);
1304			let mut u = Mat::<T>::zeros(n, if want_u { n } else { 0 });
1305			let mut v = Mat::<T>::zeros(n, if want_v { n } else { 0 });
1306			svd(
1307				a.as_ref(),
1308				s.as_mut().as_diagonal_mut(),
1309				if want_u { Some(u.as_mut()) } else { None },
1310				if want_v { Some(v.as_mut()) } else { None },
1311				Par::Seq,
1312				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1313					n,
1314					n,
1315					vectors(want_u),
1316					vectors(want_v),
1317					ParShape::Seq,
1318					params,
1319				))),
1320				params,
1321			)
1322			.unwrap();
1323			// `norm_max` skips NaN entries, so finiteness is stated separately
1324			assert!(s.is_all_finite());
1325			if want_u {
1326				assert!(u.is_all_finite());
1327			}
1328			if want_v {
1329				assert!(v.is_all_finite());
1330			}
1331			let expected =
1332				Col::<T>::from_fn(n, |i| if i == 0 { one() } else { zero() });
1333			let tol = from_f64::<T::Real>(tol);
1334			assert!((&s - &expected).norm_max() <= tol.copy());
1335			let identity = Mat::<T>::identity(n, n);
1336			if want_u && want_v {
1337				let back = &u * s.as_diagonal() * v.adjoint();
1338				assert!((&back - &a).norm_max() <= tol.copy());
1339			}
1340			if want_u {
1341				assert!(
1342					(u.adjoint() * &u - &identity).norm_max() <= tol.copy()
1343				);
1344			}
1345			if want_v {
1346				assert!(
1347					(v.adjoint() * &v - &identity).norm_max() <= tol.copy()
1348				);
1349			}
1350		}
1351		for n in [8usize, 16, 32] {
1352			for threshold in [4, 8] {
1353				for want in
1354					[(true, true), (false, false), (true, false), (false, true)]
1355				{
1356					check::<f32>(n, threshold, want, 1e-5);
1357					check::<c32>(n, threshold, want, 1e-5);
1358					check::<f64>(n, threshold, want, 1e-13);
1359					check::<c64>(n, threshold, want, 1e-13);
1360				}
1361			}
1362		}
1363		// the reported `singular_values()` call runs at the default threshold
1364		let threshold = <SvdParams as Auto<f64>>::auto().recursion_threshold;
1365		check::<f64>(256, threshold, (false, false), 1e-13);
1366	}
1367	/// An extreme but finite scale must not overflow or underflow the reduction that
1368	/// bidiagonalisation and the SVD share (#36): the spectrum and the reconstruction
1369	/// stay finite and accurate.
1370	#[test]
1371	fn test_extreme_scale() {
1372		fn check<T: ComplexField>(s: f64, tol: f64) {
1373			let scale = from_f64::<T>(s);
1374			let a = Mat::<T>::from_fn(2, 2, |i, j| match (i, j) {
1375				(0, 0) => scale.copy(),
1376				(0, 1) => scale.mul_real(from_f64::<T::Real>(0.5)),
1377				(1, 1) => scale.copy(),
1378				_ => zero(),
1379			});
1380			let params = Spec::new(SvdParams {
1381				recursion_threshold: 4,
1382				..auto!(T)
1383			});
1384			let mut sigma = Mat::<T>::zeros(2, 2);
1385			let mut u = Mat::<T>::zeros(2, 2);
1386			let mut v = Mat::<T>::zeros(2, 2);
1387			svd(
1388				a.as_ref(),
1389				sigma.as_mut().diagonal_mut(),
1390				Some(u.as_mut()),
1391				Some(v.as_mut()),
1392				Par::Seq,
1393				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1394					2,
1395					2,
1396					ComputeSvdVectors::Thin,
1397					ComputeSvdVectors::Thin,
1398					ParShape::Seq,
1399					params,
1400				))),
1401				params,
1402			)
1403			.unwrap();
1404			let rec = &u * &sigma * v.adjoint();
1405			let err = (&rec - &a).norm_max();
1406			assert!(rec.is_all_finite());
1407			assert!(err <= a.norm_max() * from_f64::<T::Real>(tol));
1408		}
1409		for s in [1.0f64, 1e160, 1e-300] {
1410			check::<f64>(s, 1e-14);
1411			check::<c64>(s, 1e-14);
1412		}
1413		for s in [1.0f64, 1e30, 1e-35] {
1414			check::<f32>(s, 1e-6);
1415			check::<c32>(s, 1e-6);
1416		}
1417	}
1418	/// A rank-one matrix's bidiagonalisation cascades down to subnormal entries, where the
1419	/// row norm's reciprocal (and, for complex input, the phase normalisation's) overflows.
1420	#[test]
1421	fn test_rank_one_subnormal_cascade() {
1422		for n in [16usize, 32, 150] {
1423			test_svd(Mat::<f32>::full(n, n, 1.0).as_ref());
1424			test_svd(Mat::<c32>::full(n, n, c32::ONE).as_ref());
1425		}
1426		for n in [32usize, 150] {
1427			test_svd(Mat::<f64>::full(n, n, 1.0).as_ref());
1428			test_svd(Mat::<c64>::full(n, n, c64::ONE).as_ref());
1429		}
1430	}
1431	#[test]
1432	fn test_zero_child_at_split() {
1433		for n in [128, 256] {
1434			for pos in [0, n / 2, n - 1] {
1435				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1436					if i == pos && j == pos { 1.0 } else { 0.0 }
1437				});
1438				check_tiny_split(a.as_ref());
1439				let complex =
1440					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1441				check_tiny_split(complex.as_ref());
1442			}
1443		}
1444	}
1445	#[test]
1446	fn test_real() {
1447		let rng = &mut StdRng::seed_from_u64(1);
1448		for (m, n) in [
1449			(3, 2),
1450			(2, 2),
1451			(4, 4),
1452			(15, 10),
1453			(10, 10),
1454			(15, 15),
1455			(50, 50),
1456			(100, 100),
1457			(150, 150),
1458			(150, 20),
1459			(20, 150),
1460		] {
1461			let mat = CwiseMatDistribution {
1462				nrows: m,
1463				ncols: n,
1464				dist: StandardNormal,
1465			}
1466			.rand::<Mat<f64>>(rng);
1467			test_svd(mat.as_ref());
1468		}
1469	}
1470	#[test]
1471	fn test_cplx() {
1472		let rng = &mut StdRng::seed_from_u64(1);
1473		for (m, n) in [
1474			(1, 1),
1475			(2, 2),
1476			(3, 2),
1477			(2, 2),
1478			(3, 3),
1479			(4, 4),
1480			(15, 10),
1481			(10, 10),
1482			(15, 15),
1483			(16, 16),
1484			(17, 17),
1485			(18, 18),
1486			(19, 19),
1487			(20, 20),
1488			(30, 30),
1489			(50, 50),
1490			(100, 100),
1491			(150, 150),
1492			(150, 20),
1493			(20, 150),
1494		] {
1495			let mat = CwiseMatDistribution {
1496				nrows: m,
1497				ncols: n,
1498				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
1499			}
1500			.rand::<Mat<c64>>(rng);
1501			test_svd(mat.as_ref());
1502		}
1503	}
1504	#[test]
1505	fn test_special() {
1506		for (m, n) in [
1507			(3, 2),
1508			(2, 2),
1509			(4, 4),
1510			(15, 10),
1511			(10, 10),
1512			(15, 15),
1513			(50, 50),
1514			(100, 100),
1515			(150, 150),
1516			(150, 20),
1517			(20, 150),
1518		] {
1519			test_svd(Mat::<f64>::zeros(m, n).as_ref());
1520			test_svd(Mat::<c64>::zeros(m, n).as_ref());
1521			test_svd(Mat::<f64>::full(m, n, 1.0).as_ref());
1522			test_svd(Mat::<c64>::full(m, n, c64::ONE).as_ref());
1523			test_svd(Mat::<f64>::identity(m, n).as_ref());
1524			test_svd(Mat::<c64>::identity(m, n).as_ref());
1525		}
1526	}
1527	#[test]
1528	fn test_zink() {
1529		let diag = [
1530			-9.931833701529301,
1531			-10.920807536026027,
1532			-52.33647796311243,
1533			2.3685025127736967,
1534			2.421701994236093,
1535			-0.5051763005624579,
1536			-0.04808263896606017,
1537			-0.003875251886338955,
1538			-0.0006413264967716465,
1539			-0.003381944152463707,
1540			2.981152313236375e-5,
1541			5.4290648208388795e-6,
1542			-6.329275972084404e-7,
1543			-6.879142344209158e-7,
1544			-5.265228263479126e-9,
1545			-2.941999902335516e-9,
1546			-1.3060984997930294e-10,
1547			7.07516117218088e-12,
1548			1.8657003929029376e-12,
1549			-6.216080089659131e-14,
1550		];
1551		let subdiag = [
1552			-57.8029649868477,
1553			17.67263066467847,
1554			8.884153814270894,
1555			-9.01998231080713,
1556			-1.028638150814966,
1557			0.22247719217200435,
1558			0.016389886745811315,
1559			-0.004090989452162578,
1560			0.00036818904090536926,
1561			-0.0031394146217732367,
1562			-7.571300829706796e-6,
1563			3.0045718718618155e-6,
1564			2.1329796886727743e-6,
1565			9.259701025627789e-8,
1566			2.2291214755992877e-9,
1567			-2.3017207713252894e-9,
1568			6.807967994979358e-11,
1569			2.1677299575405587e-12,
1570			-3.07282771050034e-13,
1571			0.0,
1572		];
1573		let n = diag.len();
1574		let params = SvdParams {
1575			recursion_threshold: 8,
1576			qr_ratio_threshold: 1.0,
1577			..auto!(f64)
1578		};
1579		let mut d = ColRef::from_slice(&diag).to_owned();
1580		let mut s = ColRef::from_slice(&subdiag).to_owned();
1581		compute_bidiag_real_svd(
1582			d.as_mut().try_as_col_major_mut().unwrap(),
1583			s.as_mut().try_as_col_major_mut().unwrap(),
1584			None,
1585			None,
1586			params,
1587			Par::Seq,
1588			MemStack::new(&mut MemBuffer::new(bidiag_real_svd_scratch::<f64>(
1589				n,
1590				false,
1591				false,
1592				ParShape::Seq,
1593				params,
1594			))),
1595		)
1596		.unwrap();
1597		assert!(d[n - 1] != 0.0);
1598	}
1599	#[test]
1600	fn test_pinv() {
1601		let rng = &mut StdRng::seed_from_u64(0);
1602		let m = 6;
1603		let n = 36;
1604		let mat = CwiseMatDistribution {
1605			nrows: m,
1606			ncols: n,
1607			dist: StandardNormal,
1608		}
1609		.rand::<Mat<f64>>(rng);
1610		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1611		let err = &mat * &pinv - Mat::<f64>::identity(m, m);
1612		assert!(err.norm_max() < 1e-10);
1613	}
1614	#[test]
1615	fn test_pinv2() {
1616		let rng = &mut StdRng::seed_from_u64(0);
1617		let m = 36;
1618		let n = 6;
1619		let mat = CwiseMatDistribution {
1620			nrows: m,
1621			ncols: n,
1622			dist: StandardNormal,
1623		}
1624		.rand::<Mat<f64>>(rng);
1625		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1626		let err = &pinv * &mat - Mat::<f64>::identity(n, n);
1627		assert!(err.norm_max() < 1e-10);
1628	}
1629	/// The phase reduction's magnitude and phase helpers outside the normal
1630	/// range: for an exact power of two `s`, the singular values of `s B` are
1631	/// `s` times those of `B`. Entries below `min_positive` were flushed to
1632	/// zero or made `inf`/NaN, entries above `1 / min_positive` overflowed.
1633	fn cplx_phase_reduction_is_scale_exact<T: ComplexField>(
1634		make: fn(f64, f64) -> T,
1635	) {
1636		let n = 3;
1637		let base_d = [make(1.8, -0.4), make(-0.3, 0.7), make(0.5, 0.5)];
1638		let base_e = [make(0.2, 0.6), make(-0.8, 0.1)];
1639		let params: SvdParams = auto!(T);
1640		let run = |s: &T::Real| {
1641			let mut d = Col::<T>::from_fn(n, |i| base_d[i].mul_real(s));
1642			let mut e = Col::<T>::from_fn(n, |i| {
1643				if i + 1 < n {
1644					base_e[i].mul_real(s)
1645				} else {
1646					zero()
1647				}
1648			});
1649			let mut u = Mat::<T>::zeros(n + 1, n + 1);
1650			let mut v = Mat::<T>::zeros(n, n);
1651			compute_bidiag_cplx_svd(
1652				d.as_mut().try_as_col_major_mut().unwrap(),
1653				e.as_mut().try_as_col_major_mut().unwrap(),
1654				Some(u.as_mut()),
1655				Some(v.as_mut()),
1656				params,
1657				Par::Seq,
1658				MemStack::new(&mut MemBuffer::new(
1659					bidiag_cplx_svd_scratch::<T>(
1660						n,
1661						true,
1662						true,
1663						ParShape::Seq,
1664						params,
1665					),
1666				)),
1667			)
1668			.unwrap();
1669			assert!(u.is_all_finite() && v.is_all_finite());
1670			d
1671		};
1672		let reference = run(&one());
1673		let sml = min_positive::<T::Real>();
1674		for s in [sml.recip(), &sml * from_f64::<T::Real>(0.0625)] {
1675			let d = run(&s);
1676			for i in 0..n {
1677				let scaled_back = d[i].real() / &s;
1678				assert!(
1679					(scaled_back - reference[i].real()).abs()
1680						<= from_f64::<T::Real>(1e3)
1681							* eps::<T::Real>() * reference[0].real()
1682				);
1683			}
1684		}
1685	}
1686	#[test]
1687	fn cplx_phase_reduction_outside_normal_range() {
1688		cplx_phase_reduction_is_scale_exact::<c64>(c64::new);
1689		cplx_phase_reduction_is_scale_exact::<c32>(|r, i| {
1690			c32::new(r as f32, i as f32)
1691		});
1692		cplx_phase_reduction_is_scale_exact::<f64>(|r, _| r);
1693		cplx_phase_reduction_is_scale_exact::<f32>(|r, _| r as f32);
1694	}
1695	/// full factors must be orthogonal even when the caller's output buffers
1696	/// hold old entries; the tall/wide QR path used to leave the top-right block
1697	/// of `U` (and, after transposition, of `V`) as it found it
1698	#[test]
1699	fn full_factors_ignore_old_output_entries() {
1700		for (m, n) in [(4, 2), (2, 4)] {
1701			let a = Mat::<f64>::identity(m, n);
1702			let size = Ord::min(m, n);
1703			let mut s = Diag::<f64>::zeros(size);
1704			let mut u = Mat::<f64>::full(m, m, 1.0);
1705			let mut v = Mat::<f64>::full(n, n, 1.0);
1706			svd(
1707				a.as_ref(),
1708				s.as_mut(),
1709				Some(u.as_mut()),
1710				Some(v.as_mut()),
1711				Par::Seq,
1712				MemStack::new(&mut MemBuffer::new(svd_scratch::<f64>(
1713					m,
1714					n,
1715					ComputeSvdVectors::Full,
1716					ComputeSvdVectors::Full,
1717					ParShape::Seq,
1718					default(),
1719				))),
1720				default(),
1721			)
1722			.unwrap();
1723			let orth_u =
1724				(u.transpose() * &u - Mat::<f64>::identity(m, m)).norm_max();
1725			let orth_v =
1726				(v.transpose() * &v - Mat::<f64>::identity(n, n)).norm_max();
1727			assert!(orth_u < 1e-12, "U orthogonality error {orth_u}");
1728			assert!(orth_v < 1e-12, "V orthogonality error {orth_v}");
1729			let rec = u.as_ref().get(.., ..size) * s.as_ref()
1730				* v.as_ref().get(.., ..size).transpose();
1731			assert!((&rec - &a).norm_max() < 1e-12);
1732		}
1733	}
1734}
faer/linalg/svd/mod.rs

faer/linalg/svd/
mod.rs
