From 195082858377374329a2757608c99be2230bc2df Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 8 Feb 2024 15:59:03 +0100 Subject: [PATCH 01/16] Impl dpnp.linalg.pinv() --- dpnp/linalg/dpnp_iface_linalg.py | 51 ++++++++++++++++++++++++++++++++ dpnp/linalg/dpnp_utils_linalg.py | 31 +++++++++++++++++++ 2 files changed, 82 insertions(+) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 2b8506130ad8..1ff7935fbd0f 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -51,6 +51,7 @@ dpnp_det, dpnp_eigh, dpnp_inv, + dpnp_pinv, dpnp_slogdet, dpnp_solve, dpnp_svd, @@ -68,6 +69,7 @@ "matrix_rank", "multi_dot", "norm", + "pinv", "qr", "solve", "svd", @@ -473,6 +475,55 @@ def multi_dot(arrays, out=None): return result +def pinv(a, rcond=1e-15, hermitian=False): + """ + Compute the (Moore-Penrose) pseudo-inverse of a matrix. + + Calculate the generalized inverse of a matrix using its + singular-value decomposition (SVD) and including all large singular values. + + For full documentation refer to :obj:`numpy.linalg.inv`. + + Parameters + ---------- + a : (..., M, N) {dpnp.ndarray, usm_ndarray} + Matrix or stack of matrices to be pseudo-inverted. + rcond : float or dpnp.ndarray of float, optional + Cutoff for small singular values. + Singular values less than or equal to ``rcond * largest_singular_value`` + are set to zero. + Default: ``1e-15``. + hermitian : bool, optional + If ``True``, a is assumed to be Hermitian (symmetric if real-valued), + enabling a more efficient method for finding singular values. + Default: ``False``. + + Returns + ------- + out : (..., N, M) dpnp.ndarray + The pseudo-inverse of a. + + Examples + -------- + The following example checks that ``a * a+ * a == a`` and + ``a+ * a * a+ == a+``: + + >>> import dpnp as np + >>> a = np.random.randn(9, 6) + >>> B = np.linalg.pinv(a) + >>> np.allclose(a, np.dot(a, np.dot(B, a))) + array([ True]) + >>> np.allclose(B, np.dot(B, np.dot(a, B))) + array([ True]) + + """ + + dpnp.check_supported_arrays_type(a) + check_stacked_2d(a) + + return dpnp_pinv(a, rcond=rcond, hermitian=hermitian) + + def norm(x1, ord=None, axis=None, keepdims=False): """ Matrix or vector norm. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 93f418831331..0b4644a5028d 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -39,6 +39,7 @@ "dpnp_det", "dpnp_eigh", "dpnp_inv", + "dpnp_pinv", "dpnp_slogdet", "dpnp_solve", "dpnp_svd", @@ -955,6 +956,36 @@ def dpnp_inv(a): return b_f +def dpnp_pinv(a, rcond=1e-15, hermitian=False): + """ + dpnp_pinv(a, rcond=1e-15, hermitian=False): + + 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). + + """ + + rcond = dpnp.array(rcond, device=a.sycl_device, sycl_queue=a.sycl_queue) + if a.size == 0: + res_type = _common_type(a) + m, n = a.shape[-2:] + if m == 0 or n == 0: + res_type = a.dtype + return dpnp.empty_like(a, shape=(a.shape[:-2] + (n, m)), dtype=res_type) + + u, s, vt = dpnp_svd(a.conj(), full_matrices=False, hermitian=hermitian) + + # discard small singular values + cutoff = rcond * dpnp.amax(s, axis=-1) + leq = s <= cutoff[..., None] + dpnp.reciprocal(s, out=s) + s[leq] = 0 + + return dpnp.matmul(vt.swapaxes(-2, -1), s[..., None] * u.swapaxes(-2, -1)) + + def dpnp_solve(a, b): """ dpnp_solve(a, b) From 095da9af686edf73971214d8c0b648dffc88705c Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 8 Feb 2024 16:00:35 +0100 Subject: [PATCH 02/16] Add cupy tests for dpnp.linalg.pinv() --- .../cupy/linalg_tests/test_solve.py | 42 +++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index cd397f6c9e1b..d0d7cc295a76 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -166,3 +166,45 @@ def test_batched_inv(self, dtype): assert a.ndim >= 3 # CuPy internally uses a batched function. with pytest.raises(xp.linalg.LinAlgError): xp.linalg.inv(a) + + +class TestPinv(unittest.TestCase): + @testing.for_dtypes("ifdFD") + @_condition.retry(10) + def check_x(self, a_shape, rcond, dtype): + a_gpu = testing.shaped_random(a_shape, dtype=dtype) + a_cpu = cupy.asnumpy(a_gpu) + a_gpu_copy = a_gpu.copy() + if not isinstance(rcond, float): + rcond = numpy.asarray(rcond) + result_cpu = numpy.linalg.pinv(a_cpu, rcond=rcond) + if not isinstance(rcond, float): + rcond = cupy.asarray(rcond) + result_gpu = cupy.linalg.pinv(a_gpu, rcond=rcond) + + assert_dtype_allclose(result_gpu, result_cpu) + testing.assert_array_equal(a_gpu_copy, a_gpu) + + 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_pinv_batched(self): + self.check_x((2, 3, 4), rcond=1e-15) + self.check_x((2, 3, 4, 5), rcond=1e-15) + + def test_pinv_batched_vector_rcond(self): + self.check_x((2, 3, 4), rcond=[0.2, 0.8]) + self.check_x((2, 3, 4, 5), rcond=[[0.2, 0.9, 0.1], [0.7, 0.2, 0.5]]) + + def test_pinv_size_0(self): + self.check_x((3, 0), rcond=1e-15) + self.check_x((0, 3), rcond=1e-15) + self.check_x((0, 0), rcond=1e-15) + self.check_x((0, 2, 3), rcond=1e-15) + self.check_x((2, 0, 3), rcond=1e-15) From b7167040ffa9a37c84c5001197f7299bbf3d5e26 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 9 Feb 2024 16:32:12 +0100 Subject: [PATCH 03/16] Add tests to test_sycl_queue and test_usm_type --- tests/test_sycl_queue.py | 42 ++++++++++++++++++++++++++++++++++++++++ tests/test_usm_type.py | 35 ++++++++++++++++++++++++++++++++- 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index de2437444034..6fca75ae3fa3 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1632,3 +1632,45 @@ def test_slogdet(shape, is_empty, device): assert_sycl_queue_equal(sign_queue, dpnp_x.sycl_queue) assert_sycl_queue_equal(logdet_queue, dpnp_x.sycl_queue) + + +@pytest.mark.parametrize( + "shape, hermitian", + [ + ((4, 4), False), + ((2, 0), False), + ((4, 4), True), + ((2, 2, 3), False), + ((0, 2, 3), False), + ((1, 0, 3), False), + ], + ids=[ + "(4, 4)", + "(2, 0)", + "(2, 2), hermitian)", + "(2, 2, 3)", + "(0, 2, 3)", + "(1, 0, 3)", + ], +) +@pytest.mark.parametrize( + "device", + valid_devices, + ids=[device.filter_string for device in valid_devices], +) +def test_pinv(shape, hermitian, device): + if hermitian: + a_np = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) + a_np = numpy.conj(a_np.T) @ a_np + else: + a_np = numpy.random.randn(*shape) + + a_dp = dpnp.array(a_np, device=device) + + B_result = dpnp.linalg.pinv(a_dp, hermitian=hermitian) + B_expected = numpy.linalg.pinv(a_np, hermitian=hermitian) + assert_allclose(B_expected, B_result, rtol=1e-3, atol=1e-4) + + B_queue = B_result.sycl_queue + + assert_sycl_queue_equal(B_queue, a_dp.sycl_queue) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 56e2a68756a3..4c920c8d0dc7 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -798,6 +798,39 @@ def test_svd(usm_type, shape, full_matrices_param, compute_uv_param): assert x.usm_type == s.usm_type +@pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) +@pytest.mark.parametrize( + "shape, hermitian", + [ + ((4, 4), False), + ((2, 0), False), + ((4, 4), True), + ((2, 2, 3), False), + ((0, 2, 3), False), + ((1, 0, 3), False), + ], + ids=[ + "(4, 4)", + "(2, 0)", + "(2, 2), hermitian)", + "(2, 2, 3)", + "(0, 2, 3)", + "(1, 0, 3)", + ], +) +def test_pinv(shape, hermitian, usm_type): + if hermitian: + a = dp.random.randn(*shape) + 1j * dp.random.randn(*shape) + a = dp.conj(a.T) @ a + else: + a = dp.random.randn(*shape) + + a = dp.array(a, usm_type=usm_type) + B = dp.lialg.pinv(a, hermitian=hermitian) + + assert a.usm_type == B.usm_type + + @pytest.mark.parametrize("usm_type", list_of_usm_types, ids=list_of_usm_types) @pytest.mark.parametrize( "shape", @@ -821,7 +854,7 @@ def test_svd(usm_type, shape, full_matrices_param, compute_uv_param): ["r", "raw", "complete", "reduced"], ids=["r", "raw", "complete", "reduced"], ) -def test_qr(shape, mode, usm_type): +def test_pinv(shape, mode, usm_type): count_elems = numpy.prod(shape) a = dp.arange(count_elems, usm_type=usm_type).reshape(shape) From d0bcbbdb51b9c2bd2a430b4801cf6586bf42e01b Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 9 Feb 2024 17:13:46 +0100 Subject: [PATCH 04/16] Add TestPinv to test_linalg.py --- tests/test_linalg.py | 129 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 129 insertions(+) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 8e32b867b851..79eed0154a4c 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1167,3 +1167,132 @@ def test_svd_errors(self): # a.ndim < 2 a_dp_ndim_1 = a_dp.flatten() assert_raises(inp.linalg.LinAlgError, inp.linalg.svd, a_dp_ndim_1) + + +class TestPinv: + def get_tol(self, dtype): + tol = 1e-06 + if dtype in (inp.float32, inp.complex64): + tol = 1e-04 + elif not has_support_aspect64() and dtype in ( + inp.int32, + inp.int64, + None, + ): + tol = 1e-05 + self._tol = tol + + def check_types_shapes(self, dp_B, np_B): + if has_support_aspect64(): + assert dp_B.dtype == np_B.dtype + else: + assert dp_B.dtype.kind == np_B.dtype.kind + + assert dp_B.shape == np_B.shape + + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + @pytest.mark.parametrize( + "shape", + [(2, 2), (3, 4), (5, 3), (16, 16), (2, 2, 2), (2, 4, 2), (2, 2, 4)], + ids=[ + "(2, 2)", + "(3, 4)", + "(5, 3)", + "(16, 16)", + "(2, 2, 2)", + "(2, 4, 2)", + "(2, 2, 4)", + ], + ) + def test_pinv(self, dtype, shape): + a = numpy.random.rand(*shape).astype(dtype) + a_dp = inp.array(a) + + B = numpy.linalg.pinv(a) + B_dp = inp.linalg.pinv(a_dp) + + self.check_types_shapes(B_dp, B) + self.get_tol(dtype) + tol = self._tol + assert_allclose(B_dp, B, rtol=tol, atol=tol) + + if a.ndim == 2: + reconstructed = inp.dot(a_dp, inp.dot(B_dp, a_dp)) + else: # a.ndim > 2 + reconstructed = inp.matmul(a_dp, inp.matmul(B_dp, a_dp)) + + assert_allclose(reconstructed, a, rtol=tol, atol=tol) + + @pytest.mark.parametrize("dtype", get_complex_dtypes()) + @pytest.mark.parametrize( + "shape", + [(2, 2), (16, 16)], + ids=["(2,2)", "(16, 16)"], + ) + def test_pinv_hermitian(self, dtype, shape): + a = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) + a = numpy.conj(a.T) @ a + + a = a.astype(dtype) + a_dp = inp.array(a) + + B = numpy.linalg.pinv(a) + B_dp = inp.linalg.pinv(a_dp) + + self.check_types_shapes(B_dp, B) + self.get_tol(dtype) + tol = self._tol + assert_allclose(B_dp, B, rtol=tol, atol=tol) + + reconstructed = inp.dot(a_dp, inp.dot(B_dp, a_dp)) + assert_allclose(reconstructed, a, rtol=tol, atol=1e-03) + + @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) + @pytest.mark.parametrize( + "shape", + [(0, 0), (0, 2), (2, 0), (2, 0, 3), (2, 3, 0), (0, 2, 3)], + ids=[ + "(0, 0)", + "(0, 2)", + "(2 ,0)", + "(2, 0, 3)", + "(2, 3, 0)", + "(0, 2, 3)", + ], + ) + def test_pinv_empty(self, dtype, shape): + a = numpy.empty(shape, dtype=dtype) + a_dp = inp.array(a) + + B = numpy.linalg.pinv(a) + B_dp = inp.linalg.pinv(a_dp) + + assert_dtype_allclose(B_dp, B) + + def test_pinv_strides(self): + a = numpy.random.rand(5, 5) + a_dp = inp.array(a) + + self.get_tol(a_dp.dtype) + tol = self._tol + + # positive strides + B = numpy.linalg.pinv(a[::2, ::2]) + B_dp = inp.linalg.pinv(a_dp[::2, ::2]) + assert_allclose(B_dp, B, rtol=tol, atol=tol) + + # negative strides + B = numpy.linalg.pinv(a[::-2, ::-2]) + B_dp = inp.linalg.pinv(a_dp[::-2, ::-2]) + assert_allclose(B_dp, B, rtol=tol, atol=tol) + + def test_pinv_errors(self): + a_dp = inp.array([[1, 2], [3, 4]], dtype="float32") + + # unsupported type + a_np = inp.asnumpy(a_dp) + assert_raises(TypeError, inp.linalg.pinv, a_np) + + # a.ndim < 2 + a_dp_ndim_1 = a_dp.flatten() + assert_raises(inp.linalg.LinAlgError, inp.linalg.pinv, a_dp_ndim_1) From b7ee11154df32ef1f82919353b13af0de7d98ae6 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 12 Feb 2024 12:09:18 +0100 Subject: [PATCH 05/16] Address remarks --- dpnp/linalg/dpnp_iface_linalg.py | 2 +- dpnp/linalg/dpnp_utils_linalg.py | 7 ++++--- tests/test_usm_type.py | 2 +- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index cbed5c4e0a1c..7fa4cd50b897 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -489,7 +489,7 @@ def pinv(a, rcond=1e-15, hermitian=False): ---------- a : (..., M, N) {dpnp.ndarray, usm_ndarray} Matrix or stack of matrices to be pseudo-inverted. - rcond : float or dpnp.ndarray of float, optional + rcond : {float, array_like}, optional Cutoff for small singular values. Singular values less than or equal to ``rcond * largest_singular_value`` are set to zero. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index a5c2b027ec69..1aab242aa51d 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -1010,18 +1010,19 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False): """ - rcond = dpnp.array(rcond, device=a.sycl_device, sycl_queue=a.sycl_queue) + rcond = dpnp.array(rcond, usm_type=a.usm_type, sycl_queue=a.sycl_queue) if a.size == 0: - res_type = _common_type(a) m, n = a.shape[-2:] if m == 0 or n == 0: res_type = a.dtype + else: + res_type = _common_type(a) return dpnp.empty_like(a, shape=(a.shape[:-2] + (n, m)), dtype=res_type) u, s, vt = dpnp_svd(a.conj(), full_matrices=False, hermitian=hermitian) # discard small singular values - cutoff = rcond * dpnp.amax(s, axis=-1) + cutoff = rcond * dpnp.max(s, axis=-1) leq = s <= cutoff[..., None] dpnp.reciprocal(s, out=s) s[leq] = 0 diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 32ae4836338b..251855e81158 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -880,7 +880,7 @@ def test_pinv(shape, hermitian, usm_type): ["r", "raw", "complete", "reduced"], ids=["r", "raw", "complete", "reduced"], ) -def test_pinv(shape, mode, usm_type): +def test_qr(shape, mode, usm_type): count_elems = numpy.prod(shape) a = dp.arange(count_elems, usm_type=usm_type).reshape(shape) From 7f03075883668ef0eb25084ec7aba3b83686f644 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 12 Feb 2024 13:32:56 +0100 Subject: [PATCH 06/16] Add additional checks for rcond parameter --- dpnp/linalg/dpnp_iface_linalg.py | 5 +++-- dpnp/linalg/dpnp_utils_linalg.py | 9 ++++++++- tests/test_linalg.py | 18 +++++++++++++++++- 3 files changed, 28 insertions(+), 4 deletions(-) diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 7fa4cd50b897..14853cd991e8 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -489,10 +489,10 @@ def pinv(a, rcond=1e-15, hermitian=False): ---------- a : (..., M, N) {dpnp.ndarray, usm_ndarray} Matrix or stack of matrices to be pseudo-inverted. - rcond : {float, array_like}, optional + rcond : {float, dpnp.ndarray, usm_ndarray}, optional Cutoff for small singular values. Singular values less than or equal to ``rcond * largest_singular_value`` - are set to zero. + are set to zero. Broadcasts against the stack of matrices. Default: ``1e-15``. hermitian : bool, optional If ``True``, a is assumed to be Hermitian (symmetric if real-valued), @@ -520,6 +520,7 @@ def pinv(a, rcond=1e-15, hermitian=False): """ dpnp.check_supported_arrays_type(a) + dpnp.check_supported_arrays_type(rcond, scalar_type=True) check_stacked_2d(a) return dpnp_pinv(a, rcond=rcond, hermitian=hermitian) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 1aab242aa51d..006e7fef011e 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -1010,7 +1010,6 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False): """ - rcond = dpnp.array(rcond, usm_type=a.usm_type, sycl_queue=a.sycl_queue) if a.size == 0: m, n = a.shape[-2:] if m == 0 or n == 0: @@ -1019,6 +1018,14 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False): res_type = _common_type(a) return dpnp.empty_like(a, shape=(a.shape[:-2] + (n, m)), dtype=res_type) + if dpnp.is_supported_array_type(rcond): + # Check that `a` and `rcond` are allocated on the same device + # and have the same queue. Otherwise, `ValueError`` will be raised. + get_usm_allocations([a, rcond]) + else: + # Allocate dpnp.ndarray if rcond is a scalar + rcond = dpnp.array(rcond, usm_type=a.usm_type, sycl_queue=a.sycl_queue) + u, s, vt = dpnp_svd(a.conj(), full_matrices=False, hermitian=hermitian) # discard small singular values diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 79eed0154a4c..9a7119b9bfce 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1289,10 +1289,26 @@ def test_pinv_strides(self): def test_pinv_errors(self): a_dp = inp.array([[1, 2], [3, 4]], dtype="float32") - # unsupported type + # unsupported type `a` a_np = inp.asnumpy(a_dp) assert_raises(TypeError, inp.linalg.pinv, a_np) + # unsupported type `rcond` + rcond = numpy.array(0.5, dtype="float32") + assert_raises(TypeError, inp.linalg.pinv, a_dp, rcond) + assert_raises(TypeError, inp.linalg.pinv, a_dp, [0.5]) + + # non-broadcastable `rcond` + rcond_dp = inp.array([0.5], dtype="float32") + assert_raises(ValueError, inp.linalg.pinv, a_dp, rcond_dp) + # a.ndim < 2 a_dp_ndim_1 = a_dp.flatten() assert_raises(inp.linalg.LinAlgError, inp.linalg.pinv, a_dp_ndim_1) + + # diffetent queue + a_queue = dpctl.SyclQueue() + rcond_queue = dpctl.SyclQueue() + a_dp_q = inp.array(a_dp, sycl_queue=a_queue) + rcond_dp_q = inp.array([0.5], dtype="float32", sycl_queue=rcond_queue) + assert_raises(ValueError, inp.linalg.pinv, a_dp_q, rcond_dp_q) From 34aa37c2117c0fab26a213fb771ce74b9dde38df Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 12 Feb 2024 13:37:07 +0100 Subject: [PATCH 07/16] Add a more efficient implementation --- dpnp/linalg/dpnp_utils_linalg.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index 006e7fef011e..b92dcae0f47f 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -1034,7 +1034,9 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False): dpnp.reciprocal(s, out=s) s[leq] = 0 - return dpnp.matmul(vt.swapaxes(-2, -1), s[..., None] * u.swapaxes(-2, -1)) + u = u.swapaxes(-2, -1) + dpnp.multiply(s[..., None], u, out=u) + return dpnp.matmul(vt.swapaxes(-2, -1), u) def dpnp_qr_batch(a, mode="reduced"): From 0766347912c5332f198ba1e0d708c6c928798330 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 12 Feb 2024 14:00:55 +0100 Subject: [PATCH 08/16] Update test_sycl_queue --- tests/test_sycl_queue.py | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 56fcfd2763af..52c3a603dbb0 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1663,20 +1663,26 @@ def test_slogdet(shape, is_empty, device): @pytest.mark.parametrize( - "shape, hermitian", + "shape, hermitian, rcond_as_array", [ - ((4, 4), False), - ((2, 0), False), - ((4, 4), True), - ((2, 2, 3), False), - ((0, 2, 3), False), - ((1, 0, 3), False), + ((4, 4), False, False), + ((4, 4), False, True), + ((2, 0), False, False), + ((4, 4), True, False), + ((4, 4), True, True), + ((2, 2, 3), False, False), + ((2, 2, 3), False, True), + ((0, 2, 3), False, False), + ((1, 0, 3), False, False), ], ids=[ "(4, 4)", + "(4, 4), rcond_as_array", "(2, 0)", "(2, 2), hermitian)", + "(2, 2), hermitian, rcond_as_array)", "(2, 2, 3)", + "(2, 2, 3), rcond_as_array", "(0, 2, 3)", "(1, 0, 3)", ], @@ -1686,7 +1692,7 @@ def test_slogdet(shape, is_empty, device): valid_devices, ids=[device.filter_string for device in valid_devices], ) -def test_pinv(shape, hermitian, device): +def test_pinv(shape, hermitian, rcond_as_array, device): if hermitian: a_np = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) a_np = numpy.conj(a_np.T) @ a_np @@ -1695,8 +1701,20 @@ def test_pinv(shape, hermitian, device): a_dp = dpnp.array(a_np, device=device) - B_result = dpnp.linalg.pinv(a_dp, hermitian=hermitian) - B_expected = numpy.linalg.pinv(a_np, hermitian=hermitian) + if rcond_as_array: + rcond_np = numpy.array(1e-15) + rcond_dp = dpnp.array(1e-15, device=device) + + B_result = dpnp.linalg.pinv(a_dp, rcond=rcond_dp, hermitian=hermitian) + B_expected = numpy.linalg.pinv( + a_np, rcond=rcond_np, hermitian=hermitian + ) + + else: + # rcond == 1e-15 by default + B_result = dpnp.linalg.pinv(a_dp, hermitian=hermitian) + B_expected = numpy.linalg.pinv(a_np, hermitian=hermitian) + assert_allclose(B_expected, B_result, rtol=1e-3, atol=1e-4) B_queue = B_result.sycl_queue From bad4ad119779a510d2c5bcda34ce187a603dc749 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Mon, 12 Feb 2024 16:29:40 +0100 Subject: [PATCH 09/16] Update TestPinv in test_linalg.py --- tests/test_linalg.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 9a7119b9bfce..eaed5d59bfdb 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1179,7 +1179,7 @@ def get_tol(self, dtype): inp.int64, None, ): - tol = 1e-05 + tol = 1e-04 self._tol = tol def check_types_shapes(self, dp_B, np_B): @@ -1221,7 +1221,7 @@ def test_pinv(self, dtype, shape): else: # a.ndim > 2 reconstructed = inp.matmul(a_dp, inp.matmul(B_dp, a_dp)) - assert_allclose(reconstructed, a, rtol=tol, atol=tol) + assert_allclose(reconstructed, a_dp, rtol=tol, atol=tol) @pytest.mark.parametrize("dtype", get_complex_dtypes()) @pytest.mark.parametrize( @@ -1241,11 +1241,11 @@ def test_pinv_hermitian(self, dtype, shape): self.check_types_shapes(B_dp, B) self.get_tol(dtype) - tol = self._tol - assert_allclose(B_dp, B, rtol=tol, atol=tol) - reconstructed = inp.dot(a_dp, inp.dot(B_dp, a_dp)) - assert_allclose(reconstructed, a, rtol=tol, atol=1e-03) + reconstructed = inp.dot(inp.dot(a_dp, B_dp), a_dp) + # TODO : calculation accuracy decreases for matrix shape (16,16) + # Find out why + assert_allclose(reconstructed, a_dp, rtol=1e-02, atol=1e-02) @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) @pytest.mark.parametrize( From 35dede1d487ec3a4238716678be87db80b13a8ee Mon Sep 17 00:00:00 2001 From: Anton <100830759+antonwolfy@users.noreply.github.com> Date: Mon, 12 Feb 2024 20:57:47 +0100 Subject: [PATCH 10/16] Update tests/test_usm_type.py --- tests/test_usm_type.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 5dad4bc15de6..5cbf1cf2f0c7 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -857,7 +857,7 @@ def test_pinv(shape, hermitian, usm_type): a = dp.random.randn(*shape) a = dp.array(a, usm_type=usm_type) - B = dp.lialg.pinv(a, hermitian=hermitian) + B = dp.linalg.pinv(a, hermitian=hermitian) assert a.usm_type == B.usm_type From e7b899a73d89f48bf153c8c14e0e9bfe437a76a7 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 13 Feb 2024 19:04:20 +0100 Subject: [PATCH 11/16] Update test_pinv_hermitian --- tests/test_linalg.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index eaed5d59bfdb..e895f6206347 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -15,6 +15,7 @@ assert_dtype_allclose, get_all_dtypes, get_complex_dtypes, + get_float_complex_dtypes, has_support_aspect64, is_cpu_device, ) @@ -1223,29 +1224,29 @@ def test_pinv(self, dtype, shape): assert_allclose(reconstructed, a_dp, rtol=tol, atol=tol) - @pytest.mark.parametrize("dtype", get_complex_dtypes()) + @pytest.mark.parametrize("dtype", get_float_complex_dtypes()) @pytest.mark.parametrize( "shape", [(2, 2), (16, 16)], - ids=["(2,2)", "(16, 16)"], + ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): - a = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) - a = numpy.conj(a.T) @ a + a = numpy.random.randn(*shape).astype(dtype) + if numpy.issubdtype(dtype, numpy.complexfloating): + a += 1j * numpy.random.randn(*shape) + a = (a + a.conj().T) / 2 - a = a.astype(dtype) a_dp = inp.array(a) - B = numpy.linalg.pinv(a) - B_dp = inp.linalg.pinv(a_dp) + B = numpy.linalg.pinv(a, hermitian=True) + B_dp = inp.linalg.pinv(a_dp, hermitian=True) self.check_types_shapes(B_dp, B) self.get_tol(dtype) + tol = self._tol reconstructed = inp.dot(inp.dot(a_dp, B_dp), a_dp) - # TODO : calculation accuracy decreases for matrix shape (16,16) - # Find out why - assert_allclose(reconstructed, a_dp, rtol=1e-02, atol=1e-02) + assert_allclose(reconstructed, a_dp, rtol=tol, atol=tol) @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True)) @pytest.mark.parametrize( From a7bb8e05a01d0b9b415750a7ec7c65be39ed1083 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Tue, 13 Feb 2024 19:10:15 +0100 Subject: [PATCH 12/16] Update test_svd_hermitian --- tests/test_linalg.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index e895f6206347..a28d1b459ac5 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1122,18 +1122,19 @@ def test_svd(self, dtype, shape): dp_a, dp_u, dp_s, dp_vt, np_u, np_s, np_vt, True ) - @pytest.mark.parametrize("dtype", get_complex_dtypes()) + @pytest.mark.parametrize("dtype", get_float_complex_dtypes()) @pytest.mark.parametrize("compute_vt", [True, False], ids=["True", "False"]) @pytest.mark.parametrize( "shape", [(2, 2), (16, 16)], - ids=["(2,2)", "(16, 16)"], + ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): - a = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) - a = numpy.conj(a.T) @ a + a = numpy.random.randn(*shape).astype(dtype) + if numpy.issubdtype(dtype, numpy.complexfloating): + a += 1j * numpy.random.randn(*shape) + a = (a + a.conj().T) / 2 - a = a.astype(dtype) dp_a = inp.array(a) if compute_vt: From 16c292c19d71e30ec66edae3eb84e2f87f147c91 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 14 Feb 2024 15:54:49 +0100 Subject: [PATCH 13/16] Use numpy.random.seed in test_linalg --- tests/test_linalg.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index a28d1b459ac5..d80ac790c14f 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -703,7 +703,10 @@ class TestQr: ids=["r", "raw", "complete", "reduced"], ) def test_qr(self, dtype, shape, mode): - a = numpy.random.rand(*shape).astype(dtype) + numpy.random.seed(76) + a = numpy.random.randn(*shape).astype(dtype) + if numpy.issubdtype(dtype, numpy.complexfloating): + a += 1j * numpy.random.randn(*shape) ia = inp.array(a) if mode == "r": @@ -773,7 +776,8 @@ def test_qr_empty(self, dtype, shape, mode): ids=["r", "raw", "complete", "reduced"], ) def test_qr_strides(self, mode): - a = numpy.random.rand(5, 5) + numpy.random.seed(76) + a = numpy.random.randn(5, 5) ia = inp.array(a) # positive strides @@ -1130,6 +1134,7 @@ def test_svd(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): + numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -1207,7 +1212,10 @@ def check_types_shapes(self, dp_B, np_B): ], ) def test_pinv(self, dtype, shape): - a = numpy.random.rand(*shape).astype(dtype) + numpy.random.seed(76) + a = numpy.random.randn(*shape).astype(dtype) + if numpy.issubdtype(dtype, numpy.complexfloating): + a += 1j * numpy.random.randn(*shape) a_dp = inp.array(a) B = numpy.linalg.pinv(a) @@ -1232,6 +1240,7 @@ def test_pinv(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): + numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -1272,7 +1281,8 @@ def test_pinv_empty(self, dtype, shape): assert_dtype_allclose(B_dp, B) def test_pinv_strides(self): - a = numpy.random.rand(5, 5) + numpy.random.seed(76) + a = numpy.random.randn(5, 5) a_dp = inp.array(a) self.get_tol(a_dp.dtype) From c08ea5e1caadefaef8f909f4a2e74a247e07e55f Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Wed, 14 Feb 2024 16:39:30 +0100 Subject: [PATCH 14/16] Use setup_method to rudece code duplication in test_linalg --- tests/test_linalg.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index d80ac790c14f..167d78a8ee60 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -14,7 +14,6 @@ from .helper import ( assert_dtype_allclose, get_all_dtypes, - get_complex_dtypes, get_float_complex_dtypes, has_support_aspect64, is_cpu_device, @@ -679,6 +678,11 @@ def test_norm3(array, ord, axis): class TestQr: + # Set numpy.random.seed for test methods to prevent + # random generation of the input singular matrix + def setup_method(self): + numpy.random.seed(76) + # TODO: New packages that fix issue CMPLRLLVM-53771 are only available in internal CI. # Skip the tests on cpu until these packages are available for the external CI. # Specifically dpcpp_linux-64>=2024.1.0 @@ -703,7 +707,6 @@ class TestQr: ids=["r", "raw", "complete", "reduced"], ) def test_qr(self, dtype, shape, mode): - numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -776,7 +779,6 @@ def test_qr_empty(self, dtype, shape, mode): ids=["r", "raw", "complete", "reduced"], ) def test_qr_strides(self, mode): - numpy.random.seed(76) a = numpy.random.randn(5, 5) ia = inp.array(a) @@ -1037,6 +1039,11 @@ def test_slogdet_errors(self): class TestSvd: + # Set numpy.random.seed for test methods to prevent + # random generation of the input singular matrix + def setup_method(self): + numpy.random.seed(76) + def get_tol(self, dtype): tol = 1e-06 if dtype in (inp.float32, inp.complex64): @@ -1134,7 +1141,6 @@ def test_svd(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_svd_hermitian(self, dtype, compute_vt, shape): - numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -1177,6 +1183,11 @@ def test_svd_errors(self): class TestPinv: + # Set numpy.random.seed for test methods to prevent + # random generation of the input singular matrix + def setup_method(self): + numpy.random.seed(76) + def get_tol(self, dtype): tol = 1e-06 if dtype in (inp.float32, inp.complex64): @@ -1212,7 +1223,6 @@ def check_types_shapes(self, dp_B, np_B): ], ) def test_pinv(self, dtype, shape): - numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -1240,7 +1250,6 @@ def test_pinv(self, dtype, shape): ids=["(2, 2)", "(16, 16)"], ) def test_pinv_hermitian(self, dtype, shape): - numpy.random.seed(76) a = numpy.random.randn(*shape).astype(dtype) if numpy.issubdtype(dtype, numpy.complexfloating): a += 1j * numpy.random.randn(*shape) @@ -1281,7 +1290,6 @@ def test_pinv_empty(self, dtype, shape): assert_dtype_allclose(B_dp, B) def test_pinv_strides(self): - numpy.random.seed(76) a = numpy.random.randn(5, 5) a_dp = inp.array(a) From a975fc4be1aa84060207d759c82a83187a5d1aa8 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Thu, 15 Feb 2024 11:36:32 +0100 Subject: [PATCH 15/16] Change the seed number --- tests/test_linalg.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 167d78a8ee60..5cf226762af0 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -681,7 +681,7 @@ class TestQr: # Set numpy.random.seed for test methods to prevent # random generation of the input singular matrix def setup_method(self): - numpy.random.seed(76) + numpy.random.seed(81) # TODO: New packages that fix issue CMPLRLLVM-53771 are only available in internal CI. # Skip the tests on cpu until these packages are available for the external CI. @@ -1042,7 +1042,7 @@ class TestSvd: # Set numpy.random.seed for test methods to prevent # random generation of the input singular matrix def setup_method(self): - numpy.random.seed(76) + numpy.random.seed(81) def get_tol(self, dtype): tol = 1e-06 @@ -1186,7 +1186,7 @@ class TestPinv: # Set numpy.random.seed for test methods to prevent # random generation of the input singular matrix def setup_method(self): - numpy.random.seed(76) + numpy.random.seed(81) def get_tol(self, dtype): tol = 1e-06 From 32f8f68ac2f43cdd09bf19294c3b555629ecc624 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Feb 2024 10:59:29 +0100 Subject: [PATCH 16/16] Use numpy.random.seed in test_usm_type and test_sycl_queue --- tests/test_sycl_queue.py | 1 + tests/test_usm_type.py | 1 + 2 files changed, 2 insertions(+) diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index 2954a9eab8ee..888891d80f6e 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -1698,6 +1698,7 @@ def test_slogdet(shape, is_empty, device): ids=[device.filter_string for device in valid_devices], ) def test_pinv(shape, hermitian, rcond_as_array, device): + numpy.random.seed(81) if hermitian: a_np = numpy.random.randn(*shape) + 1j * numpy.random.randn(*shape) a_np = numpy.conj(a_np.T) @ a_np diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index 5cbf1cf2f0c7..43f526ebcb46 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -850,6 +850,7 @@ def test_svd(usm_type, shape, full_matrices_param, compute_uv_param): ], ) def test_pinv(shape, hermitian, usm_type): + numpy.random.seed(81) if hermitian: a = dp.random.randn(*shape) + 1j * dp.random.randn(*shape) a = dp.conj(a.T) @ a