Update QR tests to avoid element-wise comparisons

dpnp/tests/test_linalg.py

@@ -3135,7 +3135,86 @@ def test_error(self):
 
 
 class TestQr:
-    @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True))
+    def gram(self, X, xp):
+        # Return Gram matrix: X^H @ X
+        return xp.conjugate(X).swapaxes(-1, -2) @ X
+
+    def get_R_from_raw(self, h, m, n, xp):
+        # Get reduced R from NumPy-style raw QR:
+        # R = triu((tril(h))^T), shape (..., k, n)
+        k = min(m, n)
+        Rt = xp.tril(h)
+        R = xp.swapaxes(Rt, -1, -2)
+        R = xp.triu(R[..., :m, :n])
+
+        return R[..., :k, :]
+
+    # QR is not unique:
+    # element-wise comparison with NumPy may differ by sign/phase.
+    # To verify correctness use mode-dependent functional checks:
+    # complete/reduced: check decomposition Q @ R = A
+    # raw/r: check invariant R^H @ R = A^H @ A
+    def check_qr(self, a_np, a_dp, mode):
+        if mode in ("complete", "reduced"):
+            res = dpnp.linalg.qr(a_dp, mode)
+            assert dpnp.allclose(res.Q @ res.R, a_dp, atol=1e-5)
+
+        # Since QR satisfies A = Q @ R with orthonormal Q (Q^H @ Q = I),
+        # validate correctness via the invariant R^H @ R == A^H @ A
+        # for raw/r modes
+        elif mode == "raw":
+            h_np, tau_np = numpy.linalg.qr(a_np, mode=mode)
+            h_dp, tau_dp = dpnp.linalg.qr(a_dp, mode=mode)
+
+            m, n = a_np.shape[-2], a_np.shape[-1]
+            Rraw_np = self.get_R_from_raw(h_np, m, n, numpy)
+            Rraw_dp = self.get_R_from_raw(h_dp, m, n, dpnp)
+
+            # Use reduced QR as a reference:
+            # reduced is validated via Q @ R == A
+            exp_res = dpnp.linalg.qr(a_dp, mode="reduced")
+            exp_R = exp_res.R
+            assert_allclose(Rraw_dp, exp_R, atol=1e-4, rtol=1e-4)
+
+            exp_dp = self.gram(a_dp, dpnp).astype(Rraw_dp.dtype)
+            exp_np = self.gram(a_np, numpy).astype(Rraw_np.dtype)
+
+            # compare R^H @ R == A^H @ A
+            assert_allclose(
+                self.gram(Rraw_dp, dpnp), exp_dp, atol=1e-4, rtol=1e-4
+            )
+            assert_allclose(
+                self.gram(Rraw_np, numpy), exp_np, atol=1e-4, rtol=1e-4
+            )
+
+            assert tau_dp.shape == tau_np.shape
+            if not has_support_aspect64(tau_dp.sycl_device):
+                if tau_np.dtype == numpy.float64:
+                    tau_np = tau_np.astype("float32")
+                elif tau_np.dtype == numpy.complex128:
+                    tau_np = tau_np.astype("complex64")
+            assert tau_dp.dtype == tau_np.dtype
+
+        else:  # mode == "r"
+            R_np = numpy.linalg.qr(a_np, mode="r")
+            R_dp = dpnp.linalg.qr(a_dp, mode="r")
+
+            # Use reduced QR as a reference:
+            # reduced is validated via Q @ R == A
+            exp_res = dpnp.linalg.qr(a_dp, mode="reduced")
+            exp_R = exp_res.R
+            assert_allclose(R_dp, exp_R, atol=1e-4, rtol=1e-4)
+
+            exp_dp = self.gram(a_dp, dpnp)
+            exp_np = self.gram(a_np, numpy)
+
+            # compare R^H @ R == A^H @ A
+            assert_allclose(self.gram(R_dp, dpnp), exp_dp, atol=1e-4, rtol=1e-4)
+            assert_allclose(
+                self.gram(R_np, numpy), exp_np, atol=1e-4, rtol=1e-4
+            )
+
+    @pytest.mark.parametrize("dtype", get_float_complex_dtypes())
     @pytest.mark.parametrize(
         "shape",
         [
@@ -3161,60 +3240,27 @@ class TestQr:
             "(2, 2, 4)",
         ],
     )
-    @pytest.mark.parametrize("mode", ["r", "raw", "complete", "reduced"])
+    @pytest.mark.parametrize("mode", ["complete", "reduced", "r", "raw"])
     def test_qr(self, dtype, shape, mode):
-        a = generate_random_numpy_array(shape, dtype, seed_value=81)
-        ia = dpnp.array(a)
-
-        if mode == "r":
-            np_r = numpy.linalg.qr(a, mode)
-            dpnp_r = dpnp.linalg.qr(ia, mode)
-        else:
-            np_q, np_r = numpy.linalg.qr(a, mode)
-
-            # check decomposition
-            if mode in ("complete", "reduced"):
-                result = dpnp.linalg.qr(ia, mode)
-                dpnp_q, dpnp_r = result.Q, result.R
-                assert dpnp.allclose(
-                    dpnp.matmul(dpnp_q, dpnp_r), ia, atol=1e-05
-                )
-            else:  # mode=="raw"
-                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
-                assert_dtype_allclose(dpnp_q, np_q, factor=24)
+        a = generate_random_numpy_array(shape, dtype, seed_value=None)
+        ia = dpnp.array(a, dtype=dtype)
 
-        if mode in ("raw", "r"):
-            assert_dtype_allclose(dpnp_r, np_r, factor=24)
+        self.check_qr(a, ia, mode)
 
-    @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True))
+    @pytest.mark.parametrize("dtype", get_float_complex_dtypes())
     @pytest.mark.parametrize(
         "shape",
         [(32, 32), (8, 16, 16)],
         ids=["(32, 32)", "(8, 16, 16)"],
     )
-    @pytest.mark.parametrize("mode", ["r", "raw", "complete", "reduced"])
+    @pytest.mark.parametrize("mode", ["complete", "reduced", "r", "raw"])
     def test_qr_large(self, dtype, shape, mode):
         a = generate_random_numpy_array(shape, dtype, seed_value=81)
         ia = dpnp.array(a)
 
-        if mode == "r":
-            np_r = numpy.linalg.qr(a, mode)
-            dpnp_r = dpnp.linalg.qr(ia, mode)
-        else:
-            np_q, np_r = numpy.linalg.qr(a, mode)
-
-            # check decomposition
-            if mode in ("complete", "reduced"):
-                result = dpnp.linalg.qr(ia, mode)
-                dpnp_q, dpnp_r = result.Q, result.R
-                assert dpnp.allclose(dpnp.matmul(dpnp_q, dpnp_r), ia, atol=1e-5)
-            else:  # mode=="raw"
-                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
-                assert_allclose(dpnp_q, np_q, atol=1e-4)
-        if mode in ("raw", "r"):
-            assert_allclose(dpnp_r, np_r, atol=1e-4)
+        self.check_qr(a, ia, mode)
 
-    @pytest.mark.parametrize("dtype", get_all_dtypes(no_bool=True))
+    @pytest.mark.parametrize("dtype", get_float_complex_dtypes())
     @pytest.mark.parametrize(
         "shape",
         [(0, 0), (0, 2), (2, 0), (2, 0, 3), (2, 3, 0), (0, 2, 3)],
@@ -3227,65 +3273,22 @@ def test_qr_large(self, dtype, shape, mode):
             "(0, 2, 3)",
         ],
     )
-    @pytest.mark.parametrize("mode", ["r", "raw", "complete", "reduced"])
+    @pytest.mark.parametrize("mode", ["complete", "reduced", "r", "raw"])
     def test_qr_empty(self, dtype, shape, mode):
         a = numpy.empty(shape, dtype=dtype)
         ia = dpnp.array(a)
 
-        if mode == "r":
-            np_r = numpy.linalg.qr(a, mode)
-            dpnp_r = dpnp.linalg.qr(ia, mode)
-        else:
-            np_q, np_r = numpy.linalg.qr(a, mode)
-
-            if mode in ("complete", "reduced"):
-                result = dpnp.linalg.qr(ia, mode)
-                dpnp_q, dpnp_r = result.Q, result.R
-            else:
-                dpnp_q, dpnp_r = dpnp.linalg.qr(ia, mode)
-
-            assert_dtype_allclose(dpnp_q, np_q)
+        self.check_qr(a, ia, mode)
 
-        assert_dtype_allclose(dpnp_r, np_r)
-
-    @pytest.mark.parametrize("mode", ["r", "raw", "complete", "reduced"])
+    @pytest.mark.parametrize("mode", ["complete", "reduced", "r", "raw"])
     def test_qr_strides(self, mode):
         a = generate_random_numpy_array((5, 5))
         ia = dpnp.array(a)
 
         # positive strides
-        if mode == "r":
-            np_r = numpy.linalg.qr(a[::2, ::2], mode)
-            dpnp_r = dpnp.linalg.qr(ia[::2, ::2], mode)
-        else:
-            np_q, np_r = numpy.linalg.qr(a[::2, ::2], mode)
-
-            if mode in ("complete", "reduced"):
-                result = dpnp.linalg.qr(ia[::2, ::2], mode)
-                dpnp_q, dpnp_r = result.Q, result.R
-            else:
-                dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::2, ::2], mode)
-
-            assert_dtype_allclose(dpnp_q, np_q)
-
-        assert_dtype_allclose(dpnp_r, np_r)
-
+        self.check_qr(a[::2, ::2], ia[::2, ::2], mode)
         # negative strides
-        if mode == "r":
-            np_r = numpy.linalg.qr(a[::-2, ::-2], mode)
-            dpnp_r = dpnp.linalg.qr(ia[::-2, ::-2], mode)
-        else:
-            np_q, np_r = numpy.linalg.qr(a[::-2, ::-2], mode)
-
-            if mode in ("complete", "reduced"):
-                result = dpnp.linalg.qr(ia[::-2, ::-2], mode)
-                dpnp_q, dpnp_r = result.Q, result.R
-            else:
-                dpnp_q, dpnp_r = dpnp.linalg.qr(ia[::-2, ::-2], mode)
-
-            assert_dtype_allclose(dpnp_q, np_q)
-
-        assert_dtype_allclose(dpnp_r, np_r)
+        self.check_qr(a[::-2, ::-2], ia[::-2, ::-2], mode)
 
     def test_qr_errors(self):
         a_dp = dpnp.array([[1, 2], [3, 5]], dtype="float32")

dpnp/tests/third_party/cupy/linalg_tests/test_decomposition.py

@@ -12,9 +12,7 @@
 # from cupy.cuda import runtime
 # from cupy.linalg import _util
 from dpnp.tests.helper import (
-    LTS_VERSION,
     has_support_aspect64,
-    is_lts_driver,
 )
 from dpnp.tests.third_party.cupy import testing
 from dpnp.tests.third_party.cupy.testing import _condition
@@ -169,6 +167,18 @@ def test_decomposition(self, dtype):
     )
 )
 class TestQRDecomposition(unittest.TestCase):
+    def _gram(self, x, xp):
+        # Gram matrix: X^H @ X
+        return xp.conjugate(x).swapaxes(-1, -2) @ x
+
+    def _get_R_from_raw(self, h, m, n, xp):
+        # Get reduced R from NumPy-style raw QR:
+        # R = triu((tril(h))^T), shape (..., k, n)
+        k = min(m, n)
+        Rt = xp.tril(h)
+        R = xp.swapaxes(Rt, -1, -2)
+        R = xp.triu(R[..., :m, :n])
+        return R[..., :k, :]
 
     @testing.for_dtypes("fdFD")
     def check_mode(self, array, mode, dtype):
@@ -184,16 +194,64 @@ def check_mode(self, array, mode, dtype):
             or numpy.lib.NumpyVersion(numpy.__version__) >= "1.22.0rc1"
         ):
             result_cpu = numpy.linalg.qr(a_cpu, mode=mode)
-            self._check_result(result_cpu, result_gpu)
+            self._check_result(result_cpu, result_gpu, a_cpu, a_gpu, mode)
+
+    # def _check_result(self, result_cpu, result_gpu):
+    #     if isinstance(result_cpu, tuple):
+    #         for b_cpu, b_gpu in zip(result_cpu, result_gpu):
+    #             assert b_cpu.dtype == b_gpu.dtype
+    #             testing.assert_allclose(b_cpu, b_gpu, atol=1e-4)
+    #     else:
+    #         assert result_cpu.dtype == result_gpu.dtype
+    #         testing.assert_allclose(result_cpu, result_gpu, atol=1e-4)
+
+    # QR is not unique:
+    # element-wise comparison with NumPy may differ by sign/phase.
+    # To verify correctness use mode-dependent functional checks:
+    # complete/reduced: check decomposition Q @ R = A
+    # raw/r: check invariant R^H @ R = A^H @ A
+    def _check_result(self, result_cpu, result_gpu, a_cpu, a_gpu, mode):
+
+        if mode in ("complete", "reduced"):
+            q_gpu, r_gpu = result_gpu
+            testing.assert_allclose(q_gpu @ r_gpu, a_gpu, atol=1e-4)
+
+        elif mode == "raw":
+            h_gpu, tau_gpu = result_gpu
+            h_cpu, tau_cpu = result_cpu
+            m, n = a_gpu.shape[-2], a_gpu.shape[-1]
+            r_gpu = self._get_R_from_raw(h_gpu, m, n, cupy)
+            r_cpu = self._get_R_from_raw(h_cpu, m, n, numpy)
+
+            exp_gpu = self._gram(a_gpu, cupy)
+            exp_cpu = self._gram(a_cpu, numpy)
+
+            testing.assert_allclose(
+                self._gram(r_gpu, cupy), exp_gpu, atol=1e-4, rtol=1e-4
+            )
+            testing.assert_allclose(
+                self._gram(r_cpu, numpy), exp_cpu, atol=1e-4, rtol=1e-4
+            )
 
-    def _check_result(self, result_cpu, result_gpu):
-        if isinstance(result_cpu, tuple):
-            for b_cpu, b_gpu in zip(result_cpu, result_gpu):
-                assert b_cpu.dtype == b_gpu.dtype
-                testing.assert_allclose(b_cpu, b_gpu, atol=1e-4)
-        else:
-            assert result_cpu.dtype == result_gpu.dtype
-            testing.assert_allclose(result_cpu, result_gpu, atol=1e-4)
+            assert tau_gpu.shape == tau_cpu.shape
+            if not has_support_aspect64(tau_gpu.sycl_device):
+                if tau_cpu.dtype == numpy.float64:
+                    tau_cpu = tau_cpu.astype("float32")
+                elif tau_cpu.dtype == numpy.complex128:
+                    tau_cpu = tau_cpu.astype("complex64")
+            assert tau_gpu.dtype == tau_cpu.dtype
+
+        else:  # mode == "r"
+            r_gpu = result_gpu
+            r_cpu = result_cpu
+            exp_gpu = self._gram(a_gpu, cupy)
+            exp_cpu = self._gram(a_cpu, numpy)
+            testing.assert_allclose(
+                self._gram(r_gpu, cupy), exp_gpu, atol=1e-4, rtol=1e-4
+            )
+            testing.assert_allclose(
+                self._gram(r_cpu, numpy), exp_cpu, atol=1e-4, rtol=1e-4
+            )
 
     @testing.fix_random()
     @_condition.repeat(3, 10)
@@ -202,19 +260,13 @@ def test_mode(self):
         self.check_mode(numpy.random.randn(3, 3), mode=self.mode)
         self.check_mode(numpy.random.randn(5, 4), mode=self.mode)
 
-    @pytest.mark.skipif(
-        is_lts_driver(version=LTS_VERSION.V1_6), reason="SAT-8375"
-    )
     @testing.with_requires("numpy>=1.22")
     @testing.fix_random()
     def test_mode_rank3(self):
         self.check_mode(numpy.random.randn(3, 2, 4), mode=self.mode)
         self.check_mode(numpy.random.randn(4, 3, 3), mode=self.mode)
         self.check_mode(numpy.random.randn(2, 5, 4), mode=self.mode)
 
-    @pytest.mark.skipif(
-        is_lts_driver(version=LTS_VERSION.V1_6), reason="SAT-8375"
-    )
     @testing.with_requires("numpy>=1.22")
     @testing.fix_random()
     def test_mode_rank4(self):