Revision

**aaae58cc70f03ac357af64aa1300ab00eaf9bb6d**authored by JeanKossaifi on**23 October 2016, 18:50:41 UTC**, committed by JeanKossaifi on**23 October 2016, 18:50:41 UTC****1 parent**1e94ea4

_partial_svd.py

```
from scipy.linalg import svd
from scipy.sparse.linalg import eigsh
import numpy as np
# Author: Jean Kossaifi
def partial_svd(matrix, n_eigenvecs=None):
"""Computes a fast partial SVD on `matrix`
if `n_eigenvecs` is specified, sparse eigendecomposition
is used on either matrix.dot(matrix.T) or matrix.T.dot(matrix)
Parameters
----------
matrix : 2D-array
n_eigenvecs : int, optional, default is None
if specified, number of eigen[vectors-values] to return
Returns
-------
U : 2D-array
of shape (matrix.shape[0], n_eigenvecs)
contains the right singular vectors
S : 1D-array
of shape (n_eigenvecs, )
contains the singular values of `matrix`
V : 2D-array
of shape (n_eigenvecs, matrix.shape[1])
contains the left singular vectors
"""
# Check that matrix is... a matrix!
if matrix.ndim != 2:
raise ValueError('matrix be a matrix. matrix.ndim is {} != 2'.format(
matrix.ndim))
# Choose what to do depending on the params
dim_1, dim_2 = matrix.shape
min_dim = min(dim_1, dim_2)
if n_eigenvecs is None or n_eigenvecs >= min_dim:
# Default on standard SVD
U, S, V = svd(matrix)
U, S, V = U[:, :n_eigenvecs], S[:n_eigenvecs], V[:n_eigenvecs, :]
return U, S, V
else:
# We can perform a partial SVD
# First choose whether to use X * X.T or X.T *X
if dim_1 < dim_2:
S, U = eigsh(matrix.dot(matrix.T), k=n_eigenvecs, which='LM')
S = np.sqrt(S)
V = np.dot(matrix.T, U * 1/S[None, :])
else:
S, V = eigsh(matrix.T.dot(matrix), k=n_eigenvecs, which='LM')
S = np.sqrt(S)
U = np.dot(matrix, V) * 1/S[None, :]
# WARNING: here, V is still the transpose of what it should be
U, S, V = U[:, ::-1], S[::-1], V[:, ::-1]
return U, S, V.T
```

Computing file changes ...