Merge 9cda6f6 into 50a9a34

Author: vlad-perevezentsev

dpnp/linalg/dpnp_iface_linalg.py

@@ -51,6 +51,7 @@
     dpnp_det,
     dpnp_eigh,
     dpnp_inv,
+    dpnp_pinv,
     dpnp_qr,
     dpnp_slogdet,
     dpnp_solve,
@@ -69,6 +70,7 @@
     "matrix_rank",
     "multi_dot",
     "norm",
+    "pinv",
     "qr",
     "solve",
     "svd",
@@ -474,6 +476,56 @@ 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, 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. Broadcasts against the stack of matrices.
+        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)
+    dpnp.check_supported_arrays_type(rcond, scalar_type=True)
+    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.

dpnp/linalg/dpnp_utils_linalg.py

@@ -39,6 +39,7 @@
     "dpnp_det",
     "dpnp_eigh",
     "dpnp_inv",
+    "dpnp_pinv",
     "dpnp_qr",
     "dpnp_slogdet",
     "dpnp_solve",
@@ -998,6 +999,46 @@ 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).
+
+    """
+
+    if a.size == 0:
+        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)
+
+    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
+    cutoff = rcond * dpnp.max(s, axis=-1)
+    leq = s <= cutoff[..., None]
+    dpnp.reciprocal(s, out=s)
+    s[leq] = 0
+
+    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"):
     """
     dpnp_qr_batch(a, mode="reduced")

tests/test_linalg.py

@@ -1167,3 +1167,148 @@ 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-04
+        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_dp, 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)
+
+        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(
+        "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`
+        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)

tests/test_sycl_queue.py

@@ -1665,3 +1665,63 @@ 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, rcond_as_array",
+    [
+        ((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)",
+    ],
+)
+@pytest.mark.parametrize(
+    "device",
+    valid_devices,
+    ids=[device.filter_string for device in valid_devices],
+)
+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
+    else:
+        a_np = numpy.random.randn(*shape)
+
+    a_dp = dpnp.array(a_np, device=device)
+
+    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
+
+    assert_sycl_queue_equal(B_queue, a_dp.sycl_queue)

tests/test_usm_type.py

@@ -829,6 +829,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",