Thin SVD

[lmmx]

I was studying PICARD (a faster version of ICA than "FastICA") and saw some mention of various types of SVD, which I'm already familiar with. One of the variants was "thin" SVD. I can see that there are other libs that save computation by only computing this "thin" variant and I was wondering if maybe ndarray-linalg could too.

I did a search of the code earlier and saw no mention of it, I'm curious if it'd be worth developing.

My understanding is that there are separate LAPACK functions/instructions you'd use for 'economy'/'thin' SVD, gesvdx being the relevant one

• Intel docs on gesvdkx https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-0/gesvdx.html

  • note the x part is what differs from the others already implemented, i.e. what makes it 'thin'

For reference I am interested in its use in whitening in this crate

• svd function:

  ndarray-linalg/lax/src/lib.rs

  Lines 422 to 437 in 49e964b

  	fn svd(
  	l: MatrixLayout,
  	calc_u: bool,
  	calc_vt: bool,
  	a: &mut [Self],
  	) -> Result<SvdOwned<Self>> {
  	use svd::*;
  	let work = SvdWork::<$s>::new(l, calc_u, calc_vt)?;
  	work.eval(a)
  	}
  	fn svddc(layout: MatrixLayout, jobz: JobSvd, a: &mut [Self]) -> Result<SvdOwned<Self>> {
  	use svddc::*;
  	let work = SvdDcWork::<$s>::new(layout, jobz)?;
  	work.eval(a)
  	}

• SvdWork struct:

  ndarray-linalg/lax/src/svd.rs

  Line 15 in 49e964b

  pub struct SvdWork<T: Scalar> {

• SVD trait:

  ndarray-linalg/ndarray-linalg/src/svd.rs

  Lines 8 to 19 in 49e964b

  	/// singular-value decomposition of matrix reference
  	pub trait SVD {
  	type U;
  	type VT;
  	type Sigma;
  	fn svd(
  	&self,
  	calc_u: bool,
  	calc_vt: bool,
  	) -> Result<(Option<Self::U>, Self::Sigma, Option<Self::VT>)>;
  	}

So two things I see there are that

1. there's this line saying it is unimplemented

ndarray-linalg/lax/src/svd.rs

Lines 63 to 67 in 49e964b

	let mut u = match ju {
	JobSvd::All => Some(vec_uninit((m * m) as usize)),
	JobSvd::None => None,
	_ => unimplemented!("SVD with partial vector output is not supported yet"),
	};

2. there's this comment making v clear which LAPACK instructions are supported

ndarray-linalg/lax/src/svd.rs

Lines 6 to 8 in 49e964b

	//! | f32 | f64 | c32 | c64 |
	//! |:-------|:-------|:-------|:-------|
	//! | sgesvd | dgesvd | cgesvd | zgesvd |

and (from right to left in that table)

ndarray-linalg/lax/src/svd.rs

Lines 180 to 181 in 49e964b

	impl_svd_work_c!(c64, lapack_sys::zgesvd_);
	impl_svd_work_c!(c32, lapack_sys::cgesvd_);

ndarray-linalg/lax/src/svd.rs

Lines 313 to 314 in 49e964b

	impl_svd_work_r!(f64, lapack_sys::dgesvd_);
	impl_svd_work_r!(f32, lapack_sys::sgesvd_);

I'd be interested in this but for now just pointing out where the relevant parts of code are, have not dug deeper! 👍