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	mod qr_stall;
941	// Instantiate the same cases and checks at each native precision.
942	mod lapack_svd_f64 {
943		type Real = f64;
944		type Scalar = Real;
945		const GUARD: Scalar = -12345.0;
946		const ROTATION_PHASE: Scalar = 1.0;
947		fn diagonal_phase(i: usize) -> Scalar {
948			if i % 2 == 0 { 1.0 } else { -1.0 }
949		}
950		use super::super::{
951			bidiag_real_svd_scratch as bidiagonal_scratch,
952			compute_bidiag_real_svd as compute_bidiagonal,
953		};
954		mod lapack {
955			include!("tests/lapack.rs");
956		}
957		mod svd_views {
958			include!("tests/svd_views.rs");
959		}
960	}
961	mod lapack_svd_f32 {
962		type Real = f32;
963		type Scalar = Real;
964		const GUARD: Scalar = -12345.0;
965		const ROTATION_PHASE: Scalar = 1.0;
966		fn diagonal_phase(i: usize) -> Scalar {
967			if i % 2 == 0 { 1.0 } else { -1.0 }
968		}
969		use super::super::{
970			bidiag_real_svd_scratch as bidiagonal_scratch,
971			compute_bidiag_real_svd as compute_bidiagonal,
972		};
973		mod lapack {
974			include!("tests/lapack.rs");
975		}
976		mod svd_views {
977			include!("tests/svd_views.rs");
978		}
979		#[test]
980		fn f32_fixtures_use_single_precision_limits() {
981			let (_, large) = lapack::diagonal_case(2, 2, 6);
982			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
983			let (_, small) = lapack::diagonal_case(2, 2, 7);
984			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
985			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
986			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
987			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
988			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
989		}
990	}
991	mod lapack_svd_c64 {
992		type Real = f64;
993		type Scalar = crate::c64;
994		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
995		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
996		fn diagonal_phase(i: usize) -> Scalar {
997			match i % 4 {
998				0 => Scalar::new(0.0, 1.0),
999				1 => Scalar::new(0.6, 0.8),
1000				2 => Scalar::new(-0.8, 0.6),
1001				_ => Scalar::new(-1.0, 0.0),
1002			}
1003		}
1004		use super::super::{
1005			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1006			compute_bidiag_cplx_svd as compute_bidiagonal,
1007		};
1008		mod lapack {
1009			include!("tests/lapack.rs");
1010		}
1011		mod svd_views {
1012			include!("tests/svd_views.rs");
1013		}
1014		mod complex_tests {
1015			include!("tests/complex.rs");
1016		}
1017	}
1018	mod lapack_svd_c32 {
1019		type Real = f32;
1020		type Scalar = crate::c32;
1021		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
1022		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
1023		fn diagonal_phase(i: usize) -> Scalar {
1024			match i % 4 {
1025				0 => Scalar::new(0.0, 1.0),
1026				1 => Scalar::new(0.6, 0.8),
1027				2 => Scalar::new(-0.8, 0.6),
1028				_ => Scalar::new(-1.0, 0.0),
1029			}
1030		}
1031		use super::super::{
1032			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1033			compute_bidiag_cplx_svd as compute_bidiagonal,
1034		};
1035		mod lapack {
1036			include!("tests/lapack.rs");
1037		}
1038		mod svd_views {
1039			include!("tests/svd_views.rs");
1040		}
1041		mod complex_tests {
1042			include!("tests/complex.rs");
1043		}
1044		#[test]
1045		fn c32_fixtures_use_single_precision_limits() {
1046			let (_, large) = lapack::diagonal_case(2, 2, 6);
1047			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
1048			let (_, small) = lapack::diagonal_case(2, 2, 7);
1049			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
1050			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
1051			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
1052			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
1053			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
1054		}
1055	}
1056	use super::*;
1057	use crate::assert;
1058	use crate::stats::prelude::*;
1059	use crate::utils::approx::*;
1060	use dyn_stack::MemBuffer;
1061	#[track_caller]
1062	fn test_svd<T: ComplexField>(mat: MatRef<'_, T>) {
1063		let (m, n) = mat.shape();
1064		let params = Spec::new(SvdParams {
1065			recursion_threshold: 8,
1066			qr_ratio_threshold: 1.0,
1067			..auto!(T)
1068		});
1069		let approx_eq = CwiseMat(
1070			ApproxEq::<T::Real>::eps()
1071				* from_f64::<T::Real>(8.0 * Ord::max(m, n) as f64).sqrt(),
1072		);
1073		{
1074			let mut s = Mat::zeros(m, n);
1075			let mut u = Mat::zeros(m, m);
1076			let mut v = Mat::zeros(n, n);
1077			svd(
1078				mat.as_ref(),
1079				s.as_mut().diagonal_mut(),
1080				Some(u.as_mut()),
1081				Some(v.as_mut()),
1082				Par::Seq,
1083				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1084					m,
1085					n,
1086					ComputeSvdVectors::Full,
1087					ComputeSvdVectors::Full,
1088					ParShape::Seq,
1089					params,
1090				))),
1091				params,
1092			)
1093			.unwrap();
1094			let reconstructed = &u * &s * v.adjoint();
1095			assert!(reconstructed ~ mat);
1096		}
1097		let size = Ord::min(m, n);
1098		let mut s = Mat::zeros(size, size);
1099		let mut u = Mat::zeros(m, size);
1100		let mut v = Mat::zeros(n, size);
1101		{
1102			svd(
1103				mat.as_ref(),
1104				s.as_mut().diagonal_mut(),
1105				Some(u.as_mut()),
1106				Some(v.as_mut()),
1107				Par::Seq,
1108				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1109					m,
1110					n,
1111					ComputeSvdVectors::Thin,
1112					ComputeSvdVectors::Thin,
1113					ParShape::Seq,
1114					params,
1115				))),
1116				params,
1117			)
1118			.unwrap();
1119			let reconstructed = &u * &s * v.adjoint();
1120			assert!(reconstructed ~ mat);
1121		}
1122		{
1123			let mut s2 = Mat::zeros(size, size);
1124			let mut u2 = Mat::zeros(m, size);
1125			svd(
1126				mat.as_ref(),
1127				s2.as_mut().diagonal_mut(),
1128				Some(u2.as_mut()),
1129				None,
1130				Par::Seq,
1131				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1132					m,
1133					n,
1134					ComputeSvdVectors::Thin,
1135					ComputeSvdVectors::No,
1136					ParShape::Seq,
1137					params,
1138				))),
1139				params,
1140			)
1141			.unwrap();
1142			assert!(s2 ~ s);
1143			assert!(u2 ~ u);
1144		}
1145		{
1146			let mut s2 = Mat::zeros(size, size);
1147			let mut v2 = Mat::zeros(n, size);
1148			svd(
1149				mat.as_ref(),
1150				s2.as_mut().diagonal_mut(),
1151				None,
1152				Some(v2.as_mut()),
1153				Par::Seq,
1154				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1155					m,
1156					n,
1157					ComputeSvdVectors::No,
1158					ComputeSvdVectors::Thin,
1159					ParShape::Seq,
1160					params,
1161				))),
1162				params,
1163			)
1164			.unwrap();
1165			assert!(s2 ~ s);
1166			assert!(v2 ~ v);
1167		}
1168		{
1169			let mut s2 = Mat::zeros(size, size);
1170			svd(
1171				mat.as_ref(),
1172				s2.as_mut().diagonal_mut(),
1173				None,
1174				None,
1175				Par::Seq,
1176				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1177					m,
1178					n,
1179					ComputeSvdVectors::No,
1180					ComputeSvdVectors::No,
1181					ParShape::Seq,
1182					params,
1183				))),
1184				params,
1185			)
1186			.unwrap();
1187			assert!(s2 ~ s);
1188		}
1189	}
1190	fn check_tiny_split<T: ComplexField>(a: MatRef<'_, T>) {
1191		let (m, n) = a.shape();
1192		let size = Ord::min(m, n);
1193		let compute = |threshold, want_u, want_v| {
1194			let params = Spec::new(SvdParams {
1195				recursion_threshold: threshold,
1196				..auto!(T)
1197			});
1198			let mut s = Col::zeros(size);
1199			let mut u = Mat::zeros(m, if want_u { size } else { 0 });
1200			let mut v = Mat::zeros(n, if want_v { size } else { 0 });
1201			let vectors = |want| {
1202				if want {
1203					ComputeSvdVectors::Thin
1204				} else {
1205					ComputeSvdVectors::No
1206				}
1207			};
1208			svd(
1209				a,
1210				s.as_mut().as_diagonal_mut(),
1211				if want_u { Some(u.as_mut()) } else { None },
1212				if want_v { Some(v.as_mut()) } else { None },
1213				Par::Seq,
1214				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1215					m,
1216					n,
1217					vectors(want_u),
1218					vectors(want_v),
1219					ParShape::Seq,
1220					params,
1221				))),
1222				params,
1223			)
1224			.unwrap();
1225			(s, u, v)
1226		};
1227		let (reference, _, _) = compute(usize::MAX, true, true);
1228		for (want_u, want_v) in
1229			[(true, true), (false, false), (true, false), (false, true)]
1230		{
1231			let (s, u, v) = compute(128, want_u, want_v);
1232			assert!((&s - &reference).norm_max() < from_f64::<T::Real>(1e-13));
1233			if want_u {
1234				assert!(
1235					(u.adjoint() * &u - Mat::<T>::identity(size, size))
1236						.norm_l2() < from_f64::<T::Real>(1e-13)
1237				);
1238			}
1239			if want_v {
1240				assert!(
1241					(v.adjoint() * &v - Mat::<T>::identity(size, size))
1242						.norm_l2() < from_f64::<T::Real>(1e-13)
1243				);
1244			}
1245			if want_u && want_v {
1246				let back = &u * s.as_diagonal() * v.adjoint();
1247				assert!(
1248					(&back - a).norm_l2() / a.norm_l2()
1249						< from_f64::<T::Real>(1e-13)
1250				);
1251			}
1252		}
1253	}
1254	#[test]
1255	fn test_v_only_rank_one_exact_scratch() {
1256		fn check<T: ComplexField>(entry: T, tol: f64) {
1257			let n = 8;
1258			let a = Mat::<T>::from_fn(n, n, |i, j| {
1259				if i == 0 && j == 0 {
1260					entry.copy()
1261				} else {
1262					zero()
1263				}
1264			});
1265			let mut params = Spec::<SvdParams, T>::default();
1266			params.recursion_threshold = 4;
1267			let mut s = Col::<T>::zeros(n);
1268			let mut v = Mat::<T>::zeros(n, n);
1269			let mut mem = MemBuffer::new(svd_scratch::<T>(
1270				n,
1271				n,
1272				ComputeSvdVectors::No,
1273				ComputeSvdVectors::Thin,
1274				ParShape::Seq,
1275				params,
1276			));
1277			svd(
1278				a.as_ref(),
1279				s.as_mut().as_diagonal_mut(),
1280				None,
1281				Some(v.as_mut()),
1282				Par::Seq,
1283				MemStack::new(&mut mem),
1284				params,
1285			)
1286			.unwrap();
1287			assert!(s.is_all_finite() && v.is_all_finite());
1288			let expected =
1289				Col::<T>::from_fn(n, |i| if i == 0 { one() } else { zero() });
1290			let tol = from_f64::<T::Real>(tol);
1291			assert!((&s - &expected).norm_max() <= tol.copy());
1292			assert!(
1293				(v.adjoint() * &v - Mat::<T>::identity(n, n)).norm_max() <= tol
1294			);
1295		}
1296		check::<f32>(1.0, 1e-5);
1297		check::<c32>(c32::new(0.0, 1.0), 1e-5);
1298		check::<f64>(1.0, 1e-13);
1299		check::<c64>(c64::new(0.0, 1.0), 1e-13);
1300	}
1301	#[test]
1302	fn test_tiny_diagonal_at_split() {
1303		for n in [128, 256] {
1304			for tiny in [0.0, 1e-15, 1e-13, 1e-10, 1e-6] {
1305				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1306					if (i == j && i < 64) || (j == i + 1 && i < 64) {
1307						1.0
1308					} else if i == 64 && j == 64 {
1309						tiny
1310					} else {
1311						0.0
1312					}
1313				});
1314				check_tiny_split(a.as_ref());
1315				let complex =
1316					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1317				check_tiny_split(complex.as_ref());
1318			}
1319		}
1320	}
1321	/// A zero block at a divide-and-conquer split must still hand its parent merge the
1322	/// identity boundary entry that carries the nonzero part of the spectrum, whatever
1323	/// singular vectors the caller asked for (#33).
1324	#[test]
1325	fn test_zero_block_at_split() {
1326		fn check<T: ComplexField>(
1327			n: usize,
1328			threshold: usize,
1329			want: (bool, bool),
1330			tol: f64,
1331		) {
1332			let (want_u, want_v) = want;
1333			let a = Mat::<T>::from_fn(n, n, |i, j| {
1334				if i == n / 2 && j == n / 2 {
1335					one()
1336				} else {
1337					zero()
1338				}
1339			});
1340			let params = Spec::new(SvdParams {
1341				recursion_threshold: threshold,
1342				..auto!(T)
1343			});
1344			let vectors = |want| {
1345				if want {
1346					ComputeSvdVectors::Full
1347				} else {
1348					ComputeSvdVectors::No
1349				}
1350			};
1351			let mut s = Col::<T>::zeros(n);
1352			let mut u = Mat::<T>::zeros(n, if want_u { n } else { 0 });
1353			let mut v = Mat::<T>::zeros(n, if want_v { n } else { 0 });
1354			svd(
1355				a.as_ref(),
1356				s.as_mut().as_diagonal_mut(),
1357				if want_u { Some(u.as_mut()) } else { None },
1358				if want_v { Some(v.as_mut()) } else { None },
1359				Par::Seq,
1360				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1361					n,
1362					n,
1363					vectors(want_u),
1364					vectors(want_v),
1365					ParShape::Seq,
1366					params,
1367				))),
1368				params,
1369			)
1370			.unwrap();
1371			// `norm_max` skips NaN entries, so finiteness is stated separately
1372			assert!(s.is_all_finite());
1373			if want_u {
1374				assert!(u.is_all_finite());
1375			}
1376			if want_v {
1377				assert!(v.is_all_finite());
1378			}
1379			let expected =
1380				Col::<T>::from_fn(n, |i| if i == 0 { one() } else { zero() });
1381			let tol = from_f64::<T::Real>(tol);
1382			assert!((&s - &expected).norm_max() <= tol.copy());
1383			let identity = Mat::<T>::identity(n, n);
1384			if want_u && want_v {
1385				let back = &u * s.as_diagonal() * v.adjoint();
1386				assert!((&back - &a).norm_max() <= tol.copy());
1387			}
1388			if want_u {
1389				assert!(
1390					(u.adjoint() * &u - &identity).norm_max() <= tol.copy()
1391				);
1392			}
1393			if want_v {
1394				assert!(
1395					(v.adjoint() * &v - &identity).norm_max() <= tol.copy()
1396				);
1397			}
1398		}
1399		for n in [8usize, 16, 32] {
1400			for threshold in [4, 8] {
1401				for want in
1402					[(true, true), (false, false), (true, false), (false, true)]
1403				{
1404					check::<f32>(n, threshold, want, 1e-5);
1405					check::<c32>(n, threshold, want, 1e-5);
1406					check::<f64>(n, threshold, want, 1e-13);
1407					check::<c64>(n, threshold, want, 1e-13);
1408				}
1409			}
1410		}
1411		// the reported `singular_values()` call runs at the default threshold
1412		let threshold = <SvdParams as Auto<f64>>::auto().recursion_threshold;
1413		check::<f64>(256, threshold, (false, false), 1e-13);
1414	}
1415	/// An extreme but finite scale must not overflow or underflow the reduction that
1416	/// bidiagonalisation and the SVD share (#36): the spectrum and the reconstruction
1417	/// stay finite and accurate.
1418	#[test]
1419	fn test_extreme_scale() {
1420		fn check<T: ComplexField>(s: f64, tol: f64) {
1421			let scale = from_f64::<T>(s);
1422			let a = Mat::<T>::from_fn(2, 2, |i, j| match (i, j) {
1423				(0, 0) => scale.copy(),
1424				(0, 1) => scale.mul_real(from_f64::<T::Real>(0.5)),
1425				(1, 1) => scale.copy(),
1426				_ => zero(),
1427			});
1428			let params = Spec::new(SvdParams {
1429				recursion_threshold: 4,
1430				..auto!(T)
1431			});
1432			let mut sigma = Mat::<T>::zeros(2, 2);
1433			let mut u = Mat::<T>::zeros(2, 2);
1434			let mut v = Mat::<T>::zeros(2, 2);
1435			svd(
1436				a.as_ref(),
1437				sigma.as_mut().diagonal_mut(),
1438				Some(u.as_mut()),
1439				Some(v.as_mut()),
1440				Par::Seq,
1441				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1442					2,
1443					2,
1444					ComputeSvdVectors::Thin,
1445					ComputeSvdVectors::Thin,
1446					ParShape::Seq,
1447					params,
1448				))),
1449				params,
1450			)
1451			.unwrap();
1452			let rec = &u * &sigma * v.adjoint();
1453			let err = (&rec - &a).norm_max();
1454			assert!(rec.is_all_finite());
1455			assert!(err <= a.norm_max() * from_f64::<T::Real>(tol));
1456		}
1457		for s in [1.0f64, 1e160, 1e-300] {
1458			check::<f64>(s, 1e-14);
1459			check::<c64>(s, 1e-14);
1460		}
1461		for s in [1.0f64, 1e30, 1e-35] {
1462			check::<f32>(s, 1e-6);
1463			check::<c32>(s, 1e-6);
1464		}
1465	}
1466	/// A rank-one matrix's bidiagonalisation cascades down to subnormal entries, where the
1467	/// row norm's reciprocal (and, for complex input, the phase normalisation's) overflows.
1468	#[test]
1469	fn test_rank_one_subnormal_cascade() {
1470		for n in [16usize, 32, 150] {
1471			test_svd(Mat::<f32>::full(n, n, 1.0).as_ref());
1472			test_svd(Mat::<c32>::full(n, n, c32::ONE).as_ref());
1473		}
1474		for n in [32usize, 150] {
1475			test_svd(Mat::<f64>::full(n, n, 1.0).as_ref());
1476			test_svd(Mat::<c64>::full(n, n, c64::ONE).as_ref());
1477		}
1478	}
1479	#[test]
1480	fn test_zero_child_at_split() {
1481		for n in [128, 256] {
1482			for pos in [0, n / 2, n - 1] {
1483				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1484					if i == pos && j == pos { 1.0 } else { 0.0 }
1485				});
1486				check_tiny_split(a.as_ref());
1487				let complex =
1488					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1489				check_tiny_split(complex.as_ref());
1490			}
1491		}
1492	}
1493	#[test]
1494	fn test_real() {
1495		let rng = &mut StdRng::seed_from_u64(1);
1496		for (m, n) in [
1497			(3, 2),
1498			(2, 2),
1499			(4, 4),
1500			(15, 10),
1501			(10, 10),
1502			(15, 15),
1503			(50, 50),
1504			(100, 100),
1505			(150, 150),
1506			(150, 20),
1507			(20, 150),
1508		] {
1509			let mat = CwiseMatDistribution {
1510				nrows: m,
1511				ncols: n,
1512				dist: StandardNormal,
1513			}
1514			.rand::<Mat<f64>>(rng);
1515			test_svd(mat.as_ref());
1516		}
1517	}
1518	#[test]
1519	fn test_cplx() {
1520		let rng = &mut StdRng::seed_from_u64(1);
1521		for (m, n) in [
1522			(1, 1),
1523			(2, 2),
1524			(3, 2),
1525			(2, 2),
1526			(3, 3),
1527			(4, 4),
1528			(15, 10),
1529			(10, 10),
1530			(15, 15),
1531			(16, 16),
1532			(17, 17),
1533			(18, 18),
1534			(19, 19),
1535			(20, 20),
1536			(30, 30),
1537			(50, 50),
1538			(100, 100),
1539			(150, 150),
1540			(150, 20),
1541			(20, 150),
1542		] {
1543			let mat = CwiseMatDistribution {
1544				nrows: m,
1545				ncols: n,
1546				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
1547			}
1548			.rand::<Mat<c64>>(rng);
1549			test_svd(mat.as_ref());
1550		}
1551	}
1552	#[test]
1553	fn test_special() {
1554		for (m, n) in [
1555			(3, 2),
1556			(2, 2),
1557			(4, 4),
1558			(15, 10),
1559			(10, 10),
1560			(15, 15),
1561			(50, 50),
1562			(100, 100),
1563			(150, 150),
1564			(150, 20),
1565			(20, 150),
1566		] {
1567			test_svd(Mat::<f64>::zeros(m, n).as_ref());
1568			test_svd(Mat::<c64>::zeros(m, n).as_ref());
1569			test_svd(Mat::<f64>::full(m, n, 1.0).as_ref());
1570			test_svd(Mat::<c64>::full(m, n, c64::ONE).as_ref());
1571			test_svd(Mat::<f64>::identity(m, n).as_ref());
1572			test_svd(Mat::<c64>::identity(m, n).as_ref());
1573		}
1574	}
1575	#[test]
1576	fn test_zink() {
1577		let diag = [
1578			-9.931833701529301,
1579			-10.920807536026027,
1580			-52.33647796311243,
1581			2.3685025127736967,
1582			2.421701994236093,
1583			-0.5051763005624579,
1584			-0.04808263896606017,
1585			-0.003875251886338955,
1586			-0.0006413264967716465,
1587			-0.003381944152463707,
1588			2.981152313236375e-5,
1589			5.4290648208388795e-6,
1590			-6.329275972084404e-7,
1591			-6.879142344209158e-7,
1592			-5.265228263479126e-9,
1593			-2.941999902335516e-9,
1594			-1.3060984997930294e-10,
1595			7.07516117218088e-12,
1596			1.8657003929029376e-12,
1597			-6.216080089659131e-14,
1598		];
1599		let subdiag = [
1600			-57.8029649868477,
1601			17.67263066467847,
1602			8.884153814270894,
1603			-9.01998231080713,
1604			-1.028638150814966,
1605			0.22247719217200435,
1606			0.016389886745811315,
1607			-0.004090989452162578,
1608			0.00036818904090536926,
1609			-0.0031394146217732367,
1610			-7.571300829706796e-6,
1611			3.0045718718618155e-6,
1612			2.1329796886727743e-6,
1613			9.259701025627789e-8,
1614			2.2291214755992877e-9,
1615			-2.3017207713252894e-9,
1616			6.807967994979358e-11,
1617			2.1677299575405587e-12,
1618			-3.07282771050034e-13,
1619			0.0,
1620		];
1621		let n = diag.len();
1622		let params = SvdParams {
1623			recursion_threshold: 8,
1624			qr_ratio_threshold: 1.0,
1625			..auto!(f64)
1626		};
1627		let mut d = ColRef::from_slice(&diag).to_owned();
1628		let mut s = ColRef::from_slice(&subdiag).to_owned();
1629		compute_bidiag_real_svd(
1630			d.as_mut().try_as_col_major_mut().unwrap(),
1631			s.as_mut().try_as_col_major_mut().unwrap(),
1632			None,
1633			None,
1634			params,
1635			Par::Seq,
1636			MemStack::new(&mut MemBuffer::new(bidiag_real_svd_scratch::<f64>(
1637				n,
1638				false,
1639				false,
1640				ParShape::Seq,
1641				params,
1642			))),
1643		)
1644		.unwrap();
1645		assert!(d[n - 1] != 0.0);
1646	}
1647	#[test]
1648	fn test_pinv() {
1649		let rng = &mut StdRng::seed_from_u64(0);
1650		let m = 6;
1651		let n = 36;
1652		let mat = CwiseMatDistribution {
1653			nrows: m,
1654			ncols: n,
1655			dist: StandardNormal,
1656		}
1657		.rand::<Mat<f64>>(rng);
1658		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1659		let err = &mat * &pinv - Mat::<f64>::identity(m, m);
1660		assert!(err.norm_max() < 1e-10);
1661	}
1662	#[test]
1663	fn test_pinv2() {
1664		let rng = &mut StdRng::seed_from_u64(0);
1665		let m = 36;
1666		let n = 6;
1667		let mat = CwiseMatDistribution {
1668			nrows: m,
1669			ncols: n,
1670			dist: StandardNormal,
1671		}
1672		.rand::<Mat<f64>>(rng);
1673		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1674		let err = &pinv * &mat - Mat::<f64>::identity(n, n);
1675		assert!(err.norm_max() < 1e-10);
1676	}
1677	/// The phase reduction's magnitude and phase helpers outside the normal
1678	/// range: for an exact power of two `s`, the singular values of `s B` are
1679	/// `s` times those of `B`. Entries below `min_positive` were flushed to
1680	/// zero or made `inf`/NaN, entries above `1 / min_positive` overflowed.
1681	fn cplx_phase_reduction_is_scale_exact<T: ComplexField>(
1682		make: fn(f64, f64) -> T,
1683	) {
1684		let n = 3;
1685		let base_d = [make(1.8, -0.4), make(-0.3, 0.7), make(0.5, 0.5)];
1686		let base_e = [make(0.2, 0.6), make(-0.8, 0.1)];
1687		let params: SvdParams = auto!(T);
1688		let run = |s: &T::Real| {
1689			let mut d = Col::<T>::from_fn(n, |i| base_d[i].mul_real(s));
1690			let mut e = Col::<T>::from_fn(n, |i| {
1691				if i + 1 < n {
1692					base_e[i].mul_real(s)
1693				} else {
1694					zero()
1695				}
1696			});
1697			let mut u = Mat::<T>::zeros(n + 1, n + 1);
1698			let mut v = Mat::<T>::zeros(n, n);
1699			compute_bidiag_cplx_svd(
1700				d.as_mut().try_as_col_major_mut().unwrap(),
1701				e.as_mut().try_as_col_major_mut().unwrap(),
1702				Some(u.as_mut()),
1703				Some(v.as_mut()),
1704				params,
1705				Par::Seq,
1706				MemStack::new(&mut MemBuffer::new(
1707					bidiag_cplx_svd_scratch::<T>(
1708						n,
1709						true,
1710						true,
1711						ParShape::Seq,
1712						params,
1713					),
1714				)),
1715			)
1716			.unwrap();
1717			assert!(u.is_all_finite() && v.is_all_finite());
1718			d
1719		};
1720		let reference = run(&one());
1721		let sml = min_positive::<T::Real>();
1722		for s in [sml.recip(), &sml * from_f64::<T::Real>(0.0625)] {
1723			let d = run(&s);
1724			for i in 0..n {
1725				let scaled_back = d[i].real() / &s;
1726				assert!(
1727					(scaled_back - reference[i].real()).abs()
1728						<= from_f64::<T::Real>(1e3)
1729							* eps::<T::Real>() * reference[0].real()
1730				);
1731			}
1732		}
1733	}
1734	#[test]
1735	fn cplx_phase_reduction_outside_normal_range() {
1736		cplx_phase_reduction_is_scale_exact::<c64>(c64::new);
1737		cplx_phase_reduction_is_scale_exact::<c32>(|r, i| {
1738			c32::new(r as f32, i as f32)
1739		});
1740		cplx_phase_reduction_is_scale_exact::<f64>(|r, _| r);
1741		cplx_phase_reduction_is_scale_exact::<f32>(|r, _| r as f32);
1742	}
1743	/// `svd` promises nonnegative singular values in nonincreasing order, for any
1744	/// finite input. This 128 x 128 matrix with three nonzero entries
1745	/// bidiagonalises to a cascade whose spectrum reaches the subnormal range; the
1746	/// divide-and-conquer merge used to classify a nonzero subnormal value with the
1747	/// deflated zeros and return `[0.596..., 1.9e-310, 0, ..., 6.3e-4, ...]` (#82).
1748	#[test]
1749	fn test_sorted_with_subnormal_at_split() {
1750		let n = 128;
1751		// (row, column, real part, imaginary part) of each nonzero entry, as
1752		// reported in #82.
1753		let entries = [
1754			(41, 80, -1.3368711811580607e-64, 3.48452997681476e-64),
1755			(67, 80, -5.871970588424443e-1, -1.065196368760547e-1),
1756			(5, 112, -2.8153454475366503e-4, 5.652970997533713e-4),
1757		];
1758		let mut a = Mat::<c64>::zeros(n, n);
1759		for (row, col, re, im) in entries {
1760			a[(row, col)] = c64::new(re, im);
1761		}
1762		let svd = a.as_ref().svd(Par::Seq).unwrap();
1763		let s = svd.S().column_vector();
1764		let values: std::vec::Vec<f64> =
1765			(0..s.nrows()).map(|i| s[i].re).collect();
1766		for i in 0..n {
1767			assert!(values[i] >= 0.0);
1768			assert!(s[i].im == 0.0);
1769			if i + 1 < n {
1770				assert!(
1771					values[i] >= values[i + 1],
1772					"s[{i}] = {:e} < s[{}] = {:e}",
1773					values[i],
1774					i + 1,
1775					values[i + 1],
1776				);
1777			}
1778		}
1779		// the two nonzero values of the input, as `zgesdd` reports them
1780		assert!((values[0] - 0.596780377486914).abs() <= 1e-14);
1781		assert!((values[1] - 6.315239590698253e-4).abs() <= 1e-17);
1782		let u = svd.U();
1783		let v = svd.V();
1784		let rec = &u * svd.S() * v.adjoint();
1785		assert!((&rec - &a).norm_max() <= 1e-14);
1786		assert!(
1787			(u.adjoint() * &u - Mat::<c64>::identity(n, n)).norm_max() <= 1e-14
1788		);
1789		assert!(
1790			(v.adjoint() * &v - Mat::<c64>::identity(n, n)).norm_max() <= 1e-14
1791		);
1792	}
1793	/// full factors must be orthogonal even when the caller's output buffers
1794	/// hold old entries; the tall/wide QR path used to leave the top-right block
1795	/// of `U` (and, after transposition, of `V`) as it found it
1796	#[test]
1797	fn full_factors_ignore_old_output_entries() {
1798		for (m, n) in [(4, 2), (2, 4)] {
1799			let a = Mat::<f64>::identity(m, n);
1800			let size = Ord::min(m, n);
1801			let mut s = Diag::<f64>::zeros(size);
1802			let mut u = Mat::<f64>::full(m, m, 1.0);
1803			let mut v = Mat::<f64>::full(n, n, 1.0);
1804			svd(
1805				a.as_ref(),
1806				s.as_mut(),
1807				Some(u.as_mut()),
1808				Some(v.as_mut()),
1809				Par::Seq,
1810				MemStack::new(&mut MemBuffer::new(svd_scratch::<f64>(
1811					m,
1812					n,
1813					ComputeSvdVectors::Full,
1814					ComputeSvdVectors::Full,
1815					ParShape::Seq,
1816					default(),
1817				))),
1818				default(),
1819			)
1820			.unwrap();
1821			let orth_u =
1822				(u.transpose() * &u - Mat::<f64>::identity(m, m)).norm_max();
1823			let orth_v =
1824				(v.transpose() * &v - Mat::<f64>::identity(n, n)).norm_max();
1825			assert!(orth_u < 1e-12, "U orthogonality error {orth_u}");
1826			assert!(orth_v < 1e-12, "V orthogonality error {orth_v}");
1827			let rec = u.as_ref().get(.., ..size) * s.as_ref()
1828				* v.as_ref().get(.., ..size).transpose();
1829			assert!((&rec - &a).norm_max() < 1e-12);
1830		}
1831	}
1832}
faer/linalg/svd/mod.rs

faer/linalg/svd/
mod.rs
