# svd_robust¶

Functions

 svd(a[, full_matrices, compute_uv, …]) Wrapper around scipy.linalg.svd() with gesvd backup plan. svd_gesvd(a[, full_matrices, compute_uv, …]) svd with LAPACK’s ‘#gesvd’ (with # = d/z for float/complex).

Module description

(More) robust version of singular value decomposition.

We often need to perform an SVD. In general, an SVD is a matrix factorization that is always well defined and should also work for ill-conditioned matrices. But sadly, both numpy.linalg.svd() and scipy.linalg.svd() fail from time to time, raising LinalgError("SVD did not converge"). The reason is that both of them call the LAPACK function #gesdd (where # depends on the data type), which takes an iterative approach that can fail. However, it is usually much faster than the alternative (and robust) #gesvd.

Our workaround is as follows: we provide a function svd() with call signature as scipy’s svd. This function is basically just a wrapper around scipy’s svd, i.e., we keep calling the faster dgesdd. But if that fails, we can still use dgesvd as a backup.

Sadly, dgesvd and zgesvd were not included into scipy until version ‘0.18.0’ (nor in numpy), which is as the time of this writing the latest stable scipy version. For scipy version newer than ‘0.18.0’, we make use of the new keyword ‘lapack_driver’ for svd, otherwise we (try to) load dgesvd and zgesvd from shared LAPACK libraries.

The tribute for the dgesvd wrapper code goes to ‘jgarcke’, originally posted at http://projects.scipy.org/numpy/ticket/990, which is now hosted at https://github.com/numpy/numpy/issues/1588 He explains a bit more in detail what fails.

The include of dgesvd to scipy was done in https://github.com/scipy/scipy/pull/5994.

Examples

The idea is that you just import the svd from this module and use it as replacement for np.linalg.svd or scipy.linalg.svd:

>>> from svd_robust import svd
>>> U, S, VT = svd([[1., 1.], [0., [1.]])