update dpnp.cond (#1773)

* update dpnp.cond

* address comments

* update pinv function

* Update dpnp/linalg/dpnp_utils_linalg.py

---------

Co-authored-by: Anton <100830759+antonwolfy@users.noreply.github.com>

Author: vtavana

dpnp/linalg/dpnp_algo_linalg.pyx

@@ -45,7 +45,6 @@ cimport numpy
 cimport dpnp.dpnp_utils as utils
 
 __all__ = [
-    "dpnp_cond",
     "dpnp_eig",
     "dpnp_eigvals",
 ]
@@ -60,30 +59,6 @@ ctypedef c_dpctl.DPCTLSyclEventRef(*custom_linalg_2in_1out_func_ptr_t)(c_dpctl.D
                                                                        const c_dpctl.DPCTLEventVectorRef)
 
 
-cpdef object dpnp_cond(object input, object p):
-    if p in ('f', 'fro'):
-        # TODO: change order='K' when support is implemented
-        input = dpnp.ravel(input, order='C')
-        sqnorm = dpnp.dot(input, input)
-        res = dpnp.sqrt(sqnorm)
-        ret = dpnp.array([res])
-    elif p == dpnp.inf:
-        dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=1)
-        ret = dpnp.max(dpnp_sum_val)
-    elif p == -dpnp.inf:
-        dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=1)
-        ret = dpnp.min(dpnp_sum_val)
-    elif p == 1:
-        dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=0)
-        ret = dpnp.max(dpnp_sum_val)
-    elif p == -1:
-        dpnp_sum_val = dpnp.sum(dpnp.abs(input), axis=0)
-        ret = dpnp.min(dpnp_sum_val)
-    else:
-        ret = dpnp.array([input.item(0)])
-    return ret
-
-
 cpdef tuple dpnp_eig(utils.dpnp_descriptor x1):
     cdef shape_type_c x1_shape = x1.shape
 

dpnp/linalg/dpnp_iface_linalg.py

@@ -48,6 +48,7 @@
     check_stacked_2d,
     check_stacked_square,
     dpnp_cholesky,
+    dpnp_cond,
     dpnp_det,
     dpnp_eigh,
     dpnp_inv,
@@ -145,32 +146,61 @@ def cholesky(a, upper=False):
     return dpnp_cholesky(a, upper=upper)
 
 
-def cond(input, p=None):
+def cond(x, p=None):
     """
     Compute the condition number of a matrix.
 
     For full documentation refer to :obj:`numpy.linalg.cond`.
 
-    Limitations
-    -----------
-    Input array is supported as :obj:`dpnp.ndarray`.
-    Parameter p=[None, 1, -1, 2, -2, dpnp.inf, -dpnp.inf, 'fro'] is supported.
+    Parameters
+    ----------
+    x : {dpnp.ndarray, usm_ndarray}
+        The matrix whose condition number is sought.
+    p : {None, 1, -1, 2, -2, inf, -inf, "fro"}, optional
+        Order of the norm used in the condition number computation.
+        inf means the `dpnp.inf` object, and the Frobenius norm is
+        the root-of-sum-of-squares norm. The default is ``None``.
+
+    Returns
+    -------
+    out : dpnp.ndarray
+        The condition number of the matrix. May be infinite.
 
     See Also
     --------
-    :obj:`dpnp.norm` : Matrix or vector norm.
-    """
+    :obj:`dpnp.linalg.norm` : Matrix or vector norm.
 
-    if not use_origin_backend(input):
-        if p in [None, 1, -1, 2, -2, dpnp.inf, -dpnp.inf, "fro"]:
-            result_obj = dpnp_cond(input, p)
-            result = dpnp.convert_single_elem_array_to_scalar(result_obj)
+    Examples
+    --------
+    >>> import dpnp as np
+    >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
+    >>> a
+    array([[ 1,  0, -1],
+           [ 0,  1,  0],
+           [ 1,  0,  1]])
+    >>> np.linalg.cond(a)
+    array(1.41421356)
+    >>> np.linalg.cond(a, 'fro')
+    array(3.16227766)
+    >>> np.linalg.cond(a, np.inf)
+    array(2.)
+    >>> np.linalg.cond(a, -np.inf)
+    array(1.)
+    >>> np.linalg.cond(a, 1)
+    array(2.)
+    >>> np.linalg.cond(a, -1)
+    array(1.)
+    >>> np.linalg.cond(a, 2)
+    array(1.41421356)
+    >>> np.linalg.cond(a, -2)
+    array(0.70710678) # may vary
+    >>> min(np.linalg.svd(a, compute_uv=False))*min(np.linalg.svd(np.linalg.inv(a), compute_uv=False))
+    array(0.70710678) # may vary
 
-            return result
-        else:
-            pass
+    """
 
-    return call_origin(numpy.linalg.cond, input, p)
+    dpnp.check_supported_arrays_type(x)
+    return dpnp_cond(x, p)
 
 
 def det(a):
@@ -409,6 +439,11 @@ def inv(a):
     out : (..., M, M) dpnp.ndarray
         (Multiplicative) inverse of the matrix a.
 
+    See Also
+    --------
+    :obj:`dpnp.linalg.cond` : Compute the condition number of a matrix.
+    :obj:`dpnp.linalg.svd` : Compute the singular value decomposition.
+
     Examples
     --------
     >>> import dpnp as np
@@ -710,11 +745,11 @@ def norm(x, ord=None, axis=None, keepdims=False):
            [ 2,  3,  4]])
 
     >>> np.linalg.norm(a)
-    array(7.745966692414834)
+    array(7.74596669)
     >>> np.linalg.norm(b)
-    array(7.745966692414834)
+    array(7.74596669)
     >>> np.linalg.norm(b, 'fro')
-    array(7.745966692414834)
+    array(7.74596669)
     >>> np.linalg.norm(a, np.inf)
     array(4.)
     >>> np.linalg.norm(b, np.inf)
@@ -733,16 +768,16 @@ def norm(x, ord=None, axis=None, keepdims=False):
     >>> np.linalg.norm(b, -1)
     array(6.)
     >>> np.linalg.norm(a, 2)
-    array(7.745966692414834)
+    array(7.74596669)
     >>> np.linalg.norm(b, 2)
-    array(7.3484692283495345)
+    array(7.34846923)
 
     >>> np.linalg.norm(a, -2)
     array(0.)
     >>> np.linalg.norm(b, -2)
-    array(1.8570331885190563e-016) # may vary
+    array(4.35106603e-18) # may vary
     >>> np.linalg.norm(a, 3)
-    array(5.8480354764257312) # may vary
+    array(5.84803548) # may vary
     >>> np.linalg.norm(a, -3)
     array(0.)
 
@@ -763,7 +798,7 @@ def norm(x, ord=None, axis=None, keepdims=False):
     >>> np.linalg.norm(m, axis=(1,2))
     array([  3.74165739,  11.22497216])
     >>> np.linalg.norm(m[0, :, :]), np.linalg.norm(m[1, :, :])
-    (array(3.7416573867739413), array(11.224972160321824))
+    (array(3.74165739), array(11.22497216))
 
     """
 

dpnp/linalg/dpnp_utils_linalg.py

@@ -38,6 +38,7 @@
     "check_stacked_2d",
     "check_stacked_square",
     "dpnp_cholesky",
+    "dpnp_cond",
     "dpnp_det",
     "dpnp_eigh",
     "dpnp_inv",
@@ -199,6 +200,11 @@ def _common_inexact_type(default_dtype, *dtypes):
     return dpnp.result_type(*inexact_dtypes)
 
 
+def _is_empty_2d(arr):
+    # check size first for efficiency
+    return arr.size == 0 and numpy.prod(arr.shape[-2:]) == 0
+
+
 def _lu_factor(a, res_type):
     """
     Compute pivoted LU decomposition.
@@ -841,6 +847,36 @@ def dpnp_cholesky(a, upper):
     return a_h
 
 
+def dpnp_cond(x, p=None):
+    """Compute the condition number of a matrix."""
+
+    if _is_empty_2d(x):
+        raise dpnp.linalg.LinAlgError("cond is not defined on empty arrays")
+    if p is None or p == 2 or p == -2:
+        s = dpnp.linalg.svd(x, compute_uv=False)
+        if p == -2:
+            r = s[..., -1] / s[..., 0]
+        else:
+            r = s[..., 0] / s[..., -1]
+    else:
+        result_t = _common_type(x)
+        # The result array will contain nans in the entries
+        # where inversion failed
+        invx = dpnp.linalg.inv(x)
+        r = dpnp.linalg.norm(x, p, axis=(-2, -1)) * dpnp.linalg.norm(
+            invx, p, axis=(-2, -1)
+        )
+        r = r.astype(result_t, copy=False)
+
+    # Convert nans to infs unless the original array had nan entries
+    nan_mask = dpnp.isnan(r)
+    if nan_mask.any():
+        nan_mask &= ~dpnp.isnan(x).any(axis=(-2, -1))
+        r[nan_mask] = dpnp.inf
+
+    return r
+
+
 def dpnp_det(a):
     """
     dpnp_det(a)
@@ -1255,18 +1291,18 @@ def dpnp_multi_dot(n, arrays, out=None):
     """Compute the dot product of two or more arrays in a single function call."""
 
     if not arrays[0].ndim in [1, 2]:
-        raise numpy.linalg.LinAlgError(
+        raise dpnp.linalg.LinAlgError(
             f"{arrays[0].ndim}-dimensional array given. First array must be 1-D or 2-D."
         )
 
     if not arrays[-1].ndim in [1, 2]:
-        raise numpy.linalg.LinAlgError(
+        raise dpnp.linalg.LinAlgError(
             f"{arrays[-1].ndim}-dimensional array given. Last array must be 1-D or 2-D."
         )
 
     for arr in arrays[1:-1]:
         if arr.ndim != 2:
-            raise numpy.linalg.LinAlgError(
+            raise dpnp.linalg.LinAlgError(
                 f"{arr.ndim}-dimensional array given. Inner arrays must be 2-D."
             )
 
@@ -1401,13 +1437,9 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False):
 
     """
 
-    if a.size == 0:
+    if _is_empty_2d(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)
+        return dpnp.empty_like(a, shape=(a.shape[:-2] + (n, m)))
 
     if dpnp.is_supported_array_type(rcond):
         # Check that `a` and `rcond` are allocated on the same device
@@ -1417,7 +1449,7 @@ def dpnp_pinv(a, rcond=1e-15, hermitian=False):
         # 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)
+    u, s, vt = dpnp_svd(dpnp.conj(a), full_matrices=False, hermitian=hermitian)
 
     # discard small singular values
     cutoff = rcond * dpnp.max(s, axis=-1)

tests/skipped_tests.tbl

@@ -37,19 +37,6 @@ tests/third_party/cupy/fft_tests/test_fft.py::TestFftn_param_23_{axes=None, norm
 
 tests/third_party/intel/test_zero_copy_test1.py::test_dpnp_interaction_with_dpctl_memory
 
-tests/test_linalg.py::test_cond[-1-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[1-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[-2-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[2-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[-2-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[2-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond["fro"-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond["fro"-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[None-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[None-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[-numpy.inf-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[numpy.inf-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-float64]
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-float32]
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-int64]

tests/skipped_tests_gpu.tbl

@@ -158,19 +158,6 @@ tests/third_party/cupy/random_tests/test_distributions.py::TestDistributionsMult
 
 tests/third_party/intel/test_zero_copy_test1.py::test_dpnp_interaction_with_dpctl_memory
 
-tests/test_linalg.py::test_cond[-1-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[1-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[-2-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[2-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[-2-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[2-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond["fro"-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond["fro"-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[None-[[1, 0, -1], [0, 1, 0], [1, 0, 1]]]
-tests/test_linalg.py::test_cond[None-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[-numpy.inf-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-tests/test_linalg.py::test_cond[numpy.inf-[[1, 2, 3], [4, 5, 6], [7, 8, 9]]]
-
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-float64]
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-float32]
 tests/test_linalg.py::test_matrix_rank[None-[0, 1]-int64]