diff --git a/dpnp/linalg/dpnp_iface_linalg.py b/dpnp/linalg/dpnp_iface_linalg.py index 2582df59a471..8d285217e22c 100644 --- a/dpnp/linalg/dpnp_iface_linalg.py +++ b/dpnp/linalg/dpnp_iface_linalg.py @@ -42,6 +42,7 @@ import dpnp from .dpnp_utils_linalg import ( + check_2d, check_stacked_2d, check_stacked_square, dpnp_cholesky, @@ -49,6 +50,7 @@ dpnp_det, dpnp_eigh, dpnp_inv, + dpnp_lstsq, dpnp_matrix_power, dpnp_matrix_rank, dpnp_multi_dot, @@ -69,6 +71,7 @@ "eigvals", "eigvalsh", "inv", + "lstsq", "matrix_power", "matrix_rank", "multi_dot", @@ -594,6 +597,79 @@ def inv(a): return dpnp_inv(a) +def lstsq(a, b, rcond=None): + """ + Return the least-squares solution to a linear matrix equation. + + For full documentation refer to :obj:`numpy.linalg.lstsq`. + + Parameters + ---------- + a : (M, N) {dpnp.ndarray, usm_ndarray} + "Coefficient" matrix. + b : {(M,), (M, K)} {dpnp.ndarray, usm_ndarray} + Ordinate or "dependent variable" values. + If `b` is two-dimensional, the least-squares solution + is calculated for each of the `K` columns of `b`. + rcond : {int, float, None}, optional + Cut-off ratio for small singular values of `a`. + For the purposes of rank determination, singular values are treated + as zero if they are smaller than `rcond` times the largest singular + value of `a`. + The default uses the machine precision times ``max(M, N)``. Passing + ``-1`` will use machine precision. + + Returns + ------- + x : {(N,), (N, K)} dpnp.ndarray + Least-squares solution. If `b` is two-dimensional, + the solutions are in the `K` columns of `x`. + residuals : {(1,), (K,), (0,)} dpnp.ndarray + Sums of squared residuals: Squared Euclidean 2-norm for each column in + ``b - a @ x``. + If the rank of `a` is < N or M <= N, this is an empty array. + If `b` is 1-dimensional, this is a (1,) shape array. + Otherwise the shape is (K,). + rank : int + Rank of matrix `a`. + s : (min(M, N),) dpnp.ndarray + Singular values of `a`. + + Examples + -------- + Fit a line, ``y = mx + c``, through some noisy data-points: + + >>> import dpnp as np + >>> x = np.array([0, 1, 2, 3]) + >>> y = np.array([-1, 0.2, 0.9, 2.1]) + + By examining the coefficients, we see that the line should have a + gradient of roughly 1 and cut the y-axis at, more or less, -1. + + We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]`` + and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`: + + >>> A = np.vstack([x, np.ones(len(x))]).T + >>> A + array([[0., 1.], + [1., 1.], + [2., 1.], + [3., 1.]]) + + >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0] + >>> m, c + (array(1.), array(-0.95)) # may vary + + """ + + dpnp.check_supported_arrays_type(a, b) + check_2d(a) + if rcond is not None and not isinstance(rcond, (int, float)): + raise TypeError("rcond must be integer, floating type, or None") + + return dpnp_lstsq(a, b, rcond=rcond) + + def matrix_power(a, n): """ Raise a square matrix to the (integer) power `n`. diff --git a/dpnp/linalg/dpnp_utils_linalg.py b/dpnp/linalg/dpnp_utils_linalg.py index accb18981a59..c1c6e7686881 100644 --- a/dpnp/linalg/dpnp_utils_linalg.py +++ b/dpnp/linalg/dpnp_utils_linalg.py @@ -42,6 +42,7 @@ "dpnp_det", "dpnp_eigh", "dpnp_inv", + "dpnp_lstsq", "dpnp_matrix_power", "dpnp_matrix_rank", "dpnp_multi_dot", @@ -503,6 +504,35 @@ def _multi_svd_norm(x, row_axis, col_axis, op): return result +def _nrm2_last_axis(x): + """ + Calculate the sum of squares along the last axis of an array. + + This function handles arrays containing real or complex numbers. + For complex data types, it computes the sum of squared magnitudes; + For real adata types, it sums the squares of the elements. + + Parameters + ---------- + x : {dpnp.ndarray, usm_ndarray} + + Returns + ------- + out : dpnp.ndarray + Sum of squares calculated along the last axis. + + """ + + real_dtype = _real_type(x.dtype) + # TODO: use dpnp.sum(dpnp.square(dpnp.view(x)), axis=-1, dtype=real_dtype) + # w/a since dpnp.view() in not implemented yet + # Сalculate and sum the squares of both real and imaginary parts for compelex array. + if dpnp.issubdtype(x.dtype, dpnp.complexfloating): + return dpnp.sum(dpnp.abs(x) ** 2, axis=-1, dtype=real_dtype) + else: + return dpnp.sum(dpnp.square(x), axis=-1, dtype=real_dtype) + + def _real_type(dtype, device=None): """ Returns the real data type corresponding to a given dpnp data type. @@ -653,6 +683,37 @@ def _triu_inplace(a, host_tasks, depends=None): return out +def check_2d(*arrays): + """ + Return ``True`` if each array in `arrays` is exactly two dimensions. + + If any array is not two-dimensional, `dpnp.linalg.LinAlgError` will be raised. + + Parameters + ---------- + arrays : {dpnp.ndarray, usm_ndarray} + A sequence of input arrays to check for dimensionality. + + Returns + ------- + out : bool + ``True`` if each array in `arrays` is exactly two-dimensional. + + Raises + ------ + dpnp.linalg.LinAlgError + If any array in `arrays` is not exactly two-dimensional. + + """ + + for a in arrays: + if a.ndim != 2: + raise dpnp.linalg.LinAlgError( + f"{a.ndim}-dimensional array given. The input " + "array must be exactly two-dimensional" + ) + + def check_stacked_2d(*arrays): """ Return ``True`` if each array in `arrays` has at least two dimensions. @@ -1224,6 +1285,58 @@ def dpnp_inv(a): return b_f +def dpnp_lstsq(a, b, rcond=None): + """ + dpnp_lstsq(a, b, rcond=None) + + Return the least-squares solution to a linear matrix equation. + + """ + + if b.ndim > 2: + raise dpnp.linalg.LinAlgError( + f"{b.ndim}-dimensional array given. The input " + "array must be exactly two-dimensional" + ) + + m, n = a.shape[-2:] + m2 = b.shape[0] + if m != m2: + raise dpnp.linalg.LinAlgError("Incompatible dimensions") + + u, s, vh = dpnp_svd(a, full_matrices=False, related_arrays=[b]) + + if rcond is None: + rcond = dpnp.finfo(s.dtype).eps * max(m, n) + elif rcond <= 0 or rcond >= 1: + # some doc of gelss/gelsd says "rcond < 0", but it's not true! + rcond = dpnp.finfo(s.dtype).eps + + # number of singular values and matrix rank + s1 = 1 / s + rank = dpnp.array( + s.size, dtype="int32", sycl_queue=s.sycl_queue, usm_type=s.usm_type + ) + if s.size > 0: + cutoff = rcond * s.max() + sing_vals = s <= cutoff + s1[sing_vals] = 0 + rank -= sing_vals.sum(dtype="int32") + + # Solve the least-squares solution + # x = vh.T.conj() @ diag(s1) @ u.T.conj() @ b + z = (dpnp.dot(b.T, u.conj()) * s1).T + x = dpnp.dot(vh.T.conj(), z) + # Calculate squared Euclidean 2-norm for each column in b - a*x + if m <= n or rank != n: + resids = dpnp.empty_like(s, shape=(0,)) + else: + e = b - a.dot(x) + resids = dpnp.atleast_1d(_nrm2_last_axis(e.T)) + + return x, resids, rank, s + + def dpnp_matrix_power(a, n): """ dpnp_matrix_power(a, n) @@ -2031,16 +2144,25 @@ def dpnp_slogdet(a): ) -def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): +def dpnp_svd_batch( + a, uv_type, s_type, full_matrices=True, compute_uv=True, related_arrays=None +): """ - dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True) + dpnp_svd_batch( + a, uv_type, s_type, full_matrices=True, compute_uv=True, related_arrays=None + ) Return the batched singular value decomposition (SVD) of a stack of matrices. """ - a_usm_type = a.usm_type - a_sycl_queue = a.sycl_queue + # Set USM type and SYCL queue to be used based on `a` + # and optionally provided `related_arrays`. + # If `related_arrays` is not provided, default to USM type and SYCL queue of `a`. + # Otherwise, determine USM type and SYCL queue using + # compute-follows-data execution model for `a` and `related arrays`. + usm_type, exec_q = get_usm_allocations([a] + (related_arrays or [])) + reshape = False batch_shape_orig = a.shape[:-2] @@ -2057,8 +2179,8 @@ def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): s = dpnp.empty( batch_shape_orig + (k,), dtype=s_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) if compute_uv: if full_matrices: @@ -2071,14 +2193,14 @@ def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): u = dpnp.empty( u_shape, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) vt = dpnp.empty( vt_shape, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) return u, s, vt else: @@ -2087,8 +2209,8 @@ def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): s = dpnp.empty( batch_shape_orig + (0,), dtype=s_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) if compute_uv: if full_matrices: @@ -2096,28 +2218,28 @@ def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): batch_shape_orig, m, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) vt = _stacked_identity( batch_shape_orig, n, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) else: u = dpnp.empty( batch_shape_orig + (m, 0), dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) vt = dpnp.empty( batch_shape_orig + (0, n), dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) return u, s, vt else: @@ -2165,10 +2287,22 @@ def dpnp_svd_batch(a, uv_type, s_type, full_matrices=True, compute_uv=True): def dpnp_svd( - a, full_matrices=True, compute_uv=True, hermitian=False, batch_call=False + a, + full_matrices=True, + compute_uv=True, + hermitian=False, + batch_call=False, + related_arrays=None, ): """ - dpnp_svd(a, full_matrices=True, compute_uv=True, hermitian=False, batch_call=False) + dpnp_svd( + a, + full_matrices=True, + compute_uv=True, + hermitian=False, + batch_call=False, + related_arrays=None, + ) Return the singular value decomposition (SVD). @@ -2201,22 +2335,38 @@ def dpnp_svd( s = dpnp.abs(s) return dpnp.sort(s)[..., ::-1] - uv_type = _common_type(a) + uv_type = ( + _common_type(a) + if not related_arrays + else _common_type(a, *related_arrays) + ) s_type = _real_type(uv_type) if a.ndim > 2: - return dpnp_svd_batch(a, uv_type, s_type, full_matrices, compute_uv) + return dpnp_svd_batch( + a, + uv_type, + s_type, + full_matrices, + compute_uv, + related_arrays=related_arrays, + ) + + # Set USM type and SYCL queue to be used based on `a` + # and optionally provided `related_arrays`. + # If `related_arrays` is not provided, default to USM type and SYCL queue of `a`. + # Otherwise, determine USM type and SYCL queue using + # compute-follows-data execution model for `a` and `related arrays`. + usm_type, exec_q = get_usm_allocations([a] + (related_arrays or [])) - a_usm_type = a.usm_type - a_sycl_queue = a.sycl_queue m, n = a.shape if m == 0 or n == 0: s = dpnp.empty( (0,), dtype=s_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) if compute_uv: if full_matrices: @@ -2229,14 +2379,14 @@ def dpnp_svd( u = dpnp.eye( *u_shape, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) vt = dpnp.eye( *vt_shape, dtype=uv_type, - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, + usm_type=usm_type, + sycl_queue=exec_q, ) return u, s, vt else: @@ -2244,14 +2394,16 @@ def dpnp_svd( # oneMKL LAPACK gesvd destroys `a` and assumes fortran-like array as input. # Allocate 'F' order memory for dpnp arrays to comply with these requirements. - a_h = dpnp.empty_like(a, order="F", dtype=uv_type) + a_h = dpnp.empty_like( + a, order="F", dtype=uv_type, usm_type=usm_type, sycl_queue=exec_q + ) a_usm_arr = dpnp.get_usm_ndarray(a) # use DPCTL tensor function to fill the сopy of the input array # from the input array a_ht_copy_ev, a_copy_ev = ti._copy_usm_ndarray_into_usm_ndarray( - src=a_usm_arr, dst=a_h.get_array(), sycl_queue=a_sycl_queue + src=a_usm_arr, dst=a_h.get_array(), sycl_queue=exec_q ) k = min(m, n) @@ -2273,26 +2425,20 @@ def dpnp_svd( # oneMKL LAPACK assumes fortran-like array as input. # Allocate 'F' order memory for dpnp output arrays to comply with these requirements. - u_h = dpnp.empty( - u_shape, - dtype=uv_type, + u_h = dpnp.empty_like( + a_h, + shape=u_shape, order="F", - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, ) - vt_h = dpnp.empty( - vt_shape, - dtype=uv_type, + vt_h = dpnp.empty_like( + a_h, + shape=vt_shape, order="F", - usm_type=a_usm_type, - sycl_queue=a_sycl_queue, - ) - s_h = dpnp.empty( - k, dtype=s_type, usm_type=a_usm_type, sycl_queue=a_sycl_queue ) + s_h = dpnp.empty_like(a_h, shape=(k,), dtype=s_type) ht_lapack_ev, _ = li._gesvd( - a_sycl_queue, + exec_q, jobu, jobvt, a_h.get_array(), diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 690e1666be07..7e0b99443d61 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -706,6 +706,89 @@ def test_inv_errors(self): assert_raises(inp.linalg.LinAlgError, inp.linalg.inv, a_dp) +class TestLstsq: + @pytest.mark.parametrize("dtype", get_all_dtypes()) + @pytest.mark.parametrize( + "a_shape, b_shape", + [ + ((2, 2), (2, 2)), + ((2, 2), (2, 3)), + ((2, 2), (2,)), + ((4, 2), (4, 2)), + ((4, 2), (4, 6)), + ((4, 2), (4,)), + ((2, 4), (2, 4)), + ((2, 4), (2, 6)), + ((2, 4), (2,)), + ], + ) + def test_lstsq(self, a_shape, b_shape, dtype): + a_np = numpy.random.rand(*a_shape).astype(dtype) + b_np = numpy.random.rand(*b_shape).astype(dtype) + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + result = inp.linalg.lstsq(a_dp, b_dp) + expected = numpy.linalg.lstsq(a_np, b_np) + + for param_dp, param_np in zip(result, expected): + assert_dtype_allclose(param_dp, param_np) + + @pytest.mark.parametrize("a_dtype", get_all_dtypes()) + @pytest.mark.parametrize("b_dtype", get_all_dtypes()) + def test_lstsq_diff_type(self, a_dtype, b_dtype): + a_np = numpy.array([[1, 2], [3, -5]], dtype=a_dtype) + b_np = numpy.array([4, 1], dtype=b_dtype) + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + expected = numpy.linalg.lstsq(a_np, b_np) + result = inp.linalg.lstsq(a_dp, b_dp) + + for param_dp, param_np in zip(result, expected): + assert_dtype_allclose(param_dp, param_np) + + @pytest.mark.parametrize("dtype", get_all_dtypes()) + @pytest.mark.parametrize( + ["m", "n", "nrhs"], + [(0, 4, 1), (0, 4, 2), (4, 0, 1), (4, 0, 2), (4, 2, 0), (0, 0, 0)], + ) + def test_lstsq_empty(self, m, n, nrhs, dtype): + a_np = numpy.arange(m * n).reshape(m, n).astype(dtype) + b_np = numpy.ones((m, nrhs)).astype(dtype) + + a_dp = inp.array(a_np) + b_dp = inp.array(b_np) + + result = inp.linalg.lstsq(a_dp, b_dp) + expected = numpy.linalg.lstsq(a_np, b_np) + + for param_dp, param_np in zip(result, expected): + assert_dtype_allclose(param_dp, param_np) + + def test_lstsq_errors(self): + a_dp = inp.array([[1, 0.5], [0.5, 1]], dtype="float32") + b_dp = inp.array(a_dp, dtype="float32") + + # diffetent queue + a_queue = dpctl.SyclQueue() + b_queue = dpctl.SyclQueue() + a_dp_q = inp.array(a_dp, sycl_queue=a_queue) + b_dp_q = inp.array(b_dp, sycl_queue=b_queue) + assert_raises(ValueError, inp.linalg.lstsq, a_dp_q, b_dp_q) + + # unsupported type `a` and `b` + a_np = inp.asnumpy(a_dp) + b_np = inp.asnumpy(b_dp) + assert_raises(TypeError, inp.linalg.lstsq, a_np, b_dp) + assert_raises(TypeError, inp.linalg.lstsq, a_dp, b_np) + + # unsupported type `rcond` + assert_raises(TypeError, inp.linalg.lstsq, a_dp, b_dp, [-1]) + + class TestMatrixPower: @pytest.mark.parametrize("dtype", get_all_dtypes()) @pytest.mark.parametrize( diff --git a/tests/test_sycl_queue.py b/tests/test_sycl_queue.py index cdd3f9f23db7..9b0f136d0e89 100644 --- a/tests/test_sycl_queue.py +++ b/tests/test_sycl_queue.py @@ -2022,3 +2022,35 @@ def test_tensorsolve(device): result_queue = result.sycl_queue assert_sycl_queue_equal(result_queue, a_dp.sycl_queue) + + +@pytest.mark.parametrize( + "device", + valid_devices, + ids=[device.filter_string for device in valid_devices], +) +@pytest.mark.parametrize( + ["m", "n", "nrhs"], + [ + (4, 2, 2), + (4, 0, 1), + (4, 2, 0), + (0, 0, 0), + ], +) +def test_lstsq(m, n, nrhs, device): + a_np = numpy.arange(m * n).reshape(m, n).astype(dpnp.default_float_type()) + b_np = numpy.ones((m, nrhs)).astype(dpnp.default_float_type()) + + a_dp = dpnp.array(a_np, device=device) + b_dp = dpnp.array(b_np, device=device) + + result_dp = dpnp.linalg.lstsq(a_dp, b_dp) + result = numpy.linalg.lstsq(a_np, b_np) + + for param_dp, param_np in zip(result_dp, result): + assert_dtype_allclose(param_dp, param_np) + + for param_dp in result_dp: + assert_sycl_queue_equal(param_dp.sycl_queue, a_dp.sycl_queue) + assert_sycl_queue_equal(param_dp.sycl_queue, b_dp.sycl_queue) diff --git a/tests/test_usm_type.py b/tests/test_usm_type.py index b74e927f66dd..0bdd969df6c3 100644 --- a/tests/test_usm_type.py +++ b/tests/test_usm_type.py @@ -1170,3 +1170,28 @@ def test_tensorsolve(usm_type_a, usm_type_b): assert a.usm_type == usm_type_a assert b.usm_type == usm_type_b assert result.usm_type == du.get_coerced_usm_type([usm_type_a, usm_type_b]) + + +@pytest.mark.parametrize("usm_type_a", list_of_usm_types, ids=list_of_usm_types) +@pytest.mark.parametrize("usm_type_b", list_of_usm_types, ids=list_of_usm_types) +@pytest.mark.parametrize( + ["m", "n", "nrhs"], + [ + (4, 2, 2), + (4, 0, 1), + (4, 2, 0), + (0, 0, 0), + ], +) +def test_lstsq(m, n, nrhs, usm_type_a, usm_type_b): + a = dp.arange(m * n, usm_type=usm_type_a).reshape(m, n) + b = dp.ones((m, nrhs), usm_type=usm_type_b) + + result = dp.linalg.lstsq(a, b) + + assert a.usm_type == usm_type_a + assert b.usm_type == usm_type_b + for param in result: + assert param.usm_type == du.get_coerced_usm_type( + [usm_type_a, usm_type_b] + ) diff --git a/tests/third_party/cupy/linalg_tests/test_solve.py b/tests/third_party/cupy/linalg_tests/test_solve.py index 38b2ae31f1c2..dae2dd1c9eef 100644 --- a/tests/third_party/cupy/linalg_tests/test_solve.py +++ b/tests/third_party/cupy/linalg_tests/test_solve.py @@ -19,6 +19,7 @@ } ) ) +@testing.fix_random() class TestSolve(unittest.TestCase): # TODO: add get_batched_gesv_limit # def setUp(self): @@ -137,9 +138,7 @@ def check_x(self, a_shape, dtype): def check_shape(self, a_shape): a = cupy.random.rand(*a_shape) - with self.assertRaises( - (numpy.linalg.LinAlgError, cupy.linalg.LinAlgError) - ): + with self.assertRaises(cupy.linalg.LinAlgError): cupy.linalg.inv(a) def test_inv(self): @@ -226,6 +225,98 @@ def test_pinv_size_0(self): self.check_x((2, 0, 3), rcond=1e-15) +class TestLstsq: + @testing.for_dtypes("ifdFD") + @testing.numpy_cupy_allclose(atol=1e-3, type_check=has_support_aspect64()) + def check_lstsq_solution( + self, a_shape, b_shape, seed, rcond, xp, dtype, singular=False + ): + if singular: + m, n = a_shape + rank = min(m, n) - 1 + a = xp.matmul( + testing.shaped_random( + (m, rank), xp, dtype=dtype, scale=3, seed=seed + ), + testing.shaped_random( + (rank, n), xp, dtype=dtype, scale=3, seed=seed + 42 + ), + ) + else: + a = testing.shaped_random(a_shape, xp, dtype=dtype, seed=seed) + b = testing.shaped_random(b_shape, xp, dtype=dtype, seed=seed + 37) + a_copy = a.copy() + b_copy = b.copy() + results = xp.linalg.lstsq(a, b, rcond) + if xp is cupy: + testing.assert_array_equal(a_copy, a) + testing.assert_array_equal(b_copy, b) + return results + + def check_invalid_shapes(self, a_shape, b_shape): + a = cupy.random.rand(*a_shape) + b = cupy.random.rand(*b_shape) + with pytest.raises(cupy.linalg.LinAlgError): + cupy.linalg.lstsq(a, b, rcond=None) + + def test_lstsq_solutions(self): + # Compares numpy.linalg.lstsq and cupy.linalg.lstsq solutions for: + # a shapes range from (3, 3) to (5, 3) and (3, 5) + # b shapes range from (i, 3) to (i, ) + # sets a random seed for deterministic testing + for i in range(3, 6): + for j in range(3, 6): + for k in range(2, 4): + seed = i + j + k + # check when b has shape (i, k) + self.check_lstsq_solution((i, j), (i, k), seed, rcond=-1) + self.check_lstsq_solution((i, j), (i, k), seed, rcond=None) + self.check_lstsq_solution((i, j), (i, k), seed, rcond=0.5) + self.check_lstsq_solution( + (i, j), (i, k), seed, rcond=1e-6, singular=True + ) + # check when b has shape (i, ) + self.check_lstsq_solution((i, j), (i,), seed + 1, rcond=-1) + self.check_lstsq_solution((i, j), (i,), seed + 1, rcond=None) + self.check_lstsq_solution((i, j), (i,), seed + 1, rcond=0.5) + self.check_lstsq_solution( + (i, j), (i,), seed + 1, rcond=1e-6, singular=True + ) + + @pytest.mark.parametrize("rcond", [-1, None, 0.5]) + @pytest.mark.parametrize("k", [None, 0, 1, 4]) + @pytest.mark.parametrize(("i", "j"), [(0, 0), (3, 0), (0, 7)]) + @testing.for_dtypes("ifdFD") + @testing.numpy_cupy_allclose(rtol=1e-6, type_check=has_support_aspect64()) + def test_lstsq_empty_matrix(self, xp, dtype, i, j, k, rcond): + a = xp.empty((i, j), dtype=dtype) + assert a.size == 0 + b_shape = (i,) if k is None else (i, k) + b = testing.shaped_random(b_shape, xp, dtype) + return xp.linalg.lstsq(a, b, rcond) + + def test_invalid_shapes(self): + self.check_invalid_shapes((4, 3), (3,)) + self.check_invalid_shapes((3, 3, 3), (2, 2)) + self.check_invalid_shapes((3, 3, 3), (3, 3)) + self.check_invalid_shapes((3, 3), (3, 3, 3)) + self.check_invalid_shapes((2, 2), (10,)) + self.check_invalid_shapes((3, 3), (2, 2)) + self.check_invalid_shapes((4, 3), (10, 3, 3)) + + # dpnp.linalg.lstsq() does not raise a FutureWarning + # because dpnp did not have a previous implementation of dpnp.linalg.lstsq() + # and there is no need to get rid of old deprecated behavior as numpy did. + @pytest.mark.skip("No support of deprecated behavior") + @testing.for_float_dtypes(no_float16=True) + @testing.numpy_cupy_allclose(atol=1e-3) + def test_warn_rcond(self, xp, dtype): + a = testing.shaped_random((3, 3), xp, dtype) + b = testing.shaped_random((3,), xp, dtype) + with testing.assert_warns(FutureWarning): + return xp.linalg.lstsq(a, b) + + class TestTensorInv(unittest.TestCase): @testing.for_dtypes("ifdFD") @_condition.retry(10) @@ -240,9 +331,7 @@ def check_x(self, a_shape, ind, dtype): def check_shape(self, a_shape, ind): a = cupy.random.rand(*a_shape) - with self.assertRaises( - (numpy.linalg.LinAlgError, cupy.linalg.LinAlgError) - ): + with self.assertRaises(cupy.linalg.LinAlgError): cupy.linalg.tensorinv(a, ind=ind) def check_ind(self, a_shape, ind): diff --git a/tests/third_party/cupy/random_tests/test_sample.py b/tests/third_party/cupy/random_tests/test_sample.py index 532f4f122f4b..2ef8f2605e5d 100644 --- a/tests/third_party/cupy/random_tests/test_sample.py +++ b/tests/third_party/cupy/random_tests/test_sample.py @@ -38,7 +38,7 @@ def test_zero_sizes(self): numpy.testing.assert_array_equal(a, cupy.array(())) -# @testing.fix_random() +@testing.fix_random() class TestRandint2(unittest.TestCase): @pytest.mark.usefixtures("allow_fall_back_on_numpy") @_condition.repeat(3, 10)