Merge pull request #459 from rezoo/pinv

Implement cupy.linalg.pinv

cupy/linalg/__init__.py

@@ -21,5 +21,6 @@
 from cupy.linalg.eigenvalue import eigvalsh  # NOQA
 
 from cupy.linalg.solve import inv  # NOQA
+from cupy.linalg.solve import pinv  # NOQA
 from cupy.linalg.solve import solve  # NOQA
 from cupy.linalg.solve import tensorsolve  # NOQA

cupy/linalg/solve.py

@@ -3,9 +3,11 @@
 import six
 
 import cupy
+from cupy.core import core
 from cupy import cuda
 from cupy.cuda import cublas
 from cupy.cuda import device
+from cupy.linalg import decomposition
 from cupy.linalg import util
 
 if cuda.cusolver_enabled:
@@ -165,7 +167,29 @@ def inv(a):
     return solve(a, b)
 
 
-# TODO(okuta): Implement pinv
+def pinv(a, rcond=1e-15):
+    '''Compute the Moore-Penrose pseudoinverse of a matrix.
+
+    It computes a pseudoinverse of a matrix ``a``, which is a generalization
+    of the inverse matrix with Singular Value Decomposition (SVD).
+    Note that it automatically removes small singular values for stability.
+
+    Args:
+        a (cupy.ndarray): The matrix with dimension ``(M, N)``
+        rcond (float): Cutoff parameter for small singular values.
+            For stability it computes the largest singular value denoted by
+            ``s``, and sets all singular values smaller than ``s`` to zero.
+
+    Returns:
+        cupy.ndarray: The pseudoinverse of ``a`` with dimension ``(N, M)``.
+
+    .. seealso:: :func:`numpy.linalg.pinv`
+    '''
+    u, s, vt = decomposition.svd(a, full_matrices=False)
+    cutoff = rcond * s.max()
+    s1 = 1 / s
+    s1[s <= cutoff] = 0
+    return core.dot(vt.T, s1[:, None] * u.T)
 
 
 # TODO(okuta): Implement tensorinv

docs/source/reference/linalg.rst

@@ -58,3 +58,4 @@ Solving linear equations
    cupy.linalg.solve
    cupy.linalg.tensorsolve
    cupy.linalg.inv
+   cupy.linalg.pinv

tests/cupy_tests/linalg_tests/test_solve.py

@@ -91,11 +91,50 @@ def check_shape(self, a_shape):
             cupy.linalg.inv(a)
 
     @condition.retry(10)
-    def test_solve(self):
+    def test_inv(self):
         self.check_x((3, 3))
         self.check_x((4, 4))
         self.check_x((5, 5))
 
     def test_invalid_shape(self):
         self.check_shape((2, 3))
         self.check_shape((4, 1))
+
+
+@unittest.skipUnless(
+    cuda.cusolver_enabled, 'Only cusolver in CUDA 8.0 is supported')
+@testing.gpu
+class TestPinv(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
+    @testing.for_float_dtypes(no_float16=True)
+    def check_x(self, a_shape, rcond, dtype):
+        a_cpu = numpy.random.randint(0, 10, size=a_shape).astype(dtype)
+        a_gpu = cupy.asarray(a_cpu)
+        result_cpu = numpy.linalg.pinv(a_cpu, rcond=rcond)
+        result_gpu = cupy.linalg.pinv(a_gpu, rcond=rcond)
+
+        self.assertEqual(result_cpu.dtype, result_gpu.dtype)
+        cupy.testing.assert_allclose(result_cpu, result_gpu, atol=1e-3)
+
+    def check_shape(self, a_shape, rcond):
+        a = cupy.random.rand(*a_shape)
+        with self.assertRaises(numpy.linalg.LinAlgError):
+            cupy.linalg.pinv(a)
+
+    @condition.retry(10)
+    def test_pinv(self):
+        self.check_x((3, 3), rcond=1e-15)
+        self.check_x((2, 4), rcond=1e-15)
+        self.check_x((3, 2), rcond=1e-15)
+
+        self.check_x((4, 4), rcond=0.3)
+        self.check_x((2, 5), rcond=0.5)
+        self.check_x((5, 3), rcond=0.6)
+
+    def test_invalid_shape(self):
+        self.check_shape((2, 3, 4), rcond=1e-15)
+        self.check_shape((2, 3, 4), rcond=0.5)
+        self.check_shape((4, 3, 2, 1), rcond=1e-14)
+        self.check_shape((4, 3, 2, 1), rcond=0.1)