From c6ef5b2ce7cfa9cb9c4598f85da795c56a4c0cc6 Mon Sep 17 00:00:00 2001 From: Masaki Saito Date: Wed, 30 Aug 2017 12:17:23 +0900 Subject: [PATCH 1/4] Implement cupy.linalg.pinv() --- cupy/linalg/__init__.py | 1 + cupy/linalg/solve.py | 28 ++++++++++++++- docs/source/reference/linalg.rst | 1 + tests/cupy_tests/linalg_tests/test_solve.py | 39 +++++++++++++++++++++ 4 files changed, 68 insertions(+), 1 deletion(-) diff --git a/cupy/linalg/__init__.py b/cupy/linalg/__init__.py index 32368a77b0b..cc4a7212742 100644 --- a/cupy/linalg/__init__.py +++ b/cupy/linalg/__init__.py @@ -20,5 +20,6 @@ from cupy.linalg.eigenvalue import eigh # NOQA from cupy.linalg.eigenvalue import eigvalsh # NOQA +from cupy.linalg.solve import pinv # NOQA from cupy.linalg.solve import solve # NOQA from cupy.linalg.solve import tensorsolve # NOQA diff --git a/cupy/linalg/solve.py b/cupy/linalg/solve.py index 351b58079ed..6ca02ee58ad 100644 --- a/cupy/linalg/solve.py +++ b/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: @@ -134,7 +136,31 @@ def tensorsolve(a, b, axes=None): # TODO(okuta): Implement inv -# 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` + ''' + util._assert_rank2(a) + + u, s, vt = decomposition.svd(a, full_matrices=False) + cutoff = rcond * s.max() + s1 = 1 / s + s1[s <= cutoff] = 0 + return core.dot(vt.transpose(), s1[:, None] * u.transpose()) # TODO(okuta): Implement tensorinv diff --git a/docs/source/reference/linalg.rst b/docs/source/reference/linalg.rst index d94dc6d59d4..6a2d68f47fb 100644 --- a/docs/source/reference/linalg.rst +++ b/docs/source/reference/linalg.rst @@ -57,3 +57,4 @@ Solving linear equations cupy.linalg.solve cupy.linalg.tensorsolve + cupy.linalg.pinv diff --git a/tests/cupy_tests/linalg_tests/test_solve.py b/tests/cupy_tests/linalg_tests/test_solve.py index 4a1bcc197cb..c7e43bc97c4 100644 --- a/tests/cupy_tests/linalg_tests/test_solve.py +++ b/tests/cupy_tests/linalg_tests/test_solve.py @@ -65,3 +65,42 @@ def check_x(self, a_shape, b_shape, dtype): def test_solve(self): self.check_x((2, 3, 6), (2, 3)) self.check_x((3, 4, 4, 3), (3, 4)) + + +@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) From ca30189d7e4fe22b3a78649042c6181e65f80f7a Mon Sep 17 00:00:00 2001 From: Masaki Saito Date: Thu, 31 Aug 2017 10:51:28 +0900 Subject: [PATCH 2/4] Add missing components --- tests/cupy_tests/linalg_tests/test_solve.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/cupy_tests/linalg_tests/test_solve.py b/tests/cupy_tests/linalg_tests/test_solve.py index cddc15f843c..b50775335b2 100644 --- a/tests/cupy_tests/linalg_tests/test_solve.py +++ b/tests/cupy_tests/linalg_tests/test_solve.py @@ -69,10 +69,12 @@ def test_solve(self): @unittest.skipUnless( cuda.cusolver_enabled, 'Only cusolver in CUDA 8.0 is supported') +@testing.gpu class TestInv(unittest.TestCase): _multiprocess_can_split_ = True + @testing.for_float_dtypes(no_float16=True) def check_x(self, a_shape, dtype): a_cpu = numpy.random.randint(0, 10, size=a_shape).astype(dtype) a_gpu = cupy.asarray(a_cpu) @@ -99,10 +101,12 @@ def test_invalid_shape(self): @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) From 3e76f3b112fcec18b3e55841c55f639daf96eca7 Mon Sep 17 00:00:00 2001 From: Masaki Saito Date: Thu, 31 Aug 2017 11:04:03 +0900 Subject: [PATCH 3/4] Remove unnecessary assert --- cupy/linalg/solve.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/cupy/linalg/solve.py b/cupy/linalg/solve.py index 81097276106..d4cf81a2e66 100644 --- a/cupy/linalg/solve.py +++ b/cupy/linalg/solve.py @@ -185,8 +185,6 @@ def pinv(a, rcond=1e-15): .. seealso:: :func:`numpy.linalg.pinv` ''' - util._assert_rank2(a) - u, s, vt = decomposition.svd(a, full_matrices=False) cutoff = rcond * s.max() s1 = 1 / s From a7adcd8222f462ce3d3a055e7197120a0d1615f8 Mon Sep 17 00:00:00 2001 From: Masaki Saito Date: Tue, 5 Sep 2017 21:14:03 +0900 Subject: [PATCH 4/4] Use .T instead of .transpose() --- cupy/linalg/solve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cupy/linalg/solve.py b/cupy/linalg/solve.py index d4cf81a2e66..b9a7bdede72 100644 --- a/cupy/linalg/solve.py +++ b/cupy/linalg/solve.py @@ -189,7 +189,7 @@ def pinv(a, rcond=1e-15): cutoff = rcond * s.max() s1 = 1 / s s1[s <= cutoff] = 0 - return core.dot(vt.transpose(), s1[:, None] * u.transpose()) + return core.dot(vt.T, s1[:, None] * u.T) # TODO(okuta): Implement tensorinv