Skip to content

Commit 707c5a0

Browse files
Move QR helper functions to qr_helper.py
1 parent 706b3e2 commit 707c5a0

File tree

3 files changed

+87
-148
lines changed

3 files changed

+87
-148
lines changed

dpnp/tests/qr_helper.py

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
import numpy
2+
3+
from .helper import has_support_aspect64
4+
5+
6+
def gram(x, xp):
7+
# Return Gram matrix: X^H @ X
8+
return xp.conjugate(x).swapaxes(-1, -2) @ x
9+
10+
11+
def get_R_from_raw(h, m, n, xp):
12+
# Get reduced R from NumPy-style raw QR:
13+
# R = triu((tril(h))^T), shape (..., k, n)
14+
k = min(m, n)
15+
rt = xp.tril(h)
16+
r = xp.swapaxes(rt, -1, -2)
17+
r = xp.triu(r[..., :m, :n])
18+
return r[..., :k, :]
19+
20+
21+
def check_qr(a_np, a_xp, mode, xp):
22+
# QR is not unique:
23+
# element-wise comparison with NumPy may differ by sign/phase.
24+
# To verify correctness use mode-dependent functional checks:
25+
# complete/reduced: check decomposition Q @ R = A
26+
# raw/r: check invariant R^H @ R = A^H @ A
27+
if mode in ("complete", "reduced"):
28+
res = xp.linalg.qr(a_xp, mode)
29+
assert xp.allclose(res.Q @ res.R, a_xp, atol=1e-5)
30+
31+
# Since QR satisfies A = Q @ R with orthonormal Q (Q^H @ Q = I),
32+
# validate correctness via the invariant R^H @ R == A^H @ A
33+
# for raw/r modes
34+
elif mode == "raw":
35+
_, tau_np = numpy.linalg.qr(a_np, mode=mode)
36+
h_xp, tau_xp = xp.linalg.qr(a_xp, mode=mode)
37+
38+
m, n = a_np.shape[-2], a_np.shape[-1]
39+
Rraw_xp = get_R_from_raw(h_xp, m, n, xp)
40+
41+
# Use reduced QR as a reference:
42+
# reduced is validated via Q @ R == A
43+
exp_res = xp.linalg.qr(a_xp, mode="reduced")
44+
exp_r = exp_res.R
45+
assert xp.allclose(Rraw_xp, exp_r, atol=1e-4, rtol=1e-4)
46+
47+
exp_xp = gram(a_xp, xp)
48+
49+
# Compare R^H @ R == A^H @ A
50+
assert xp.allclose(gram(Rraw_xp, xp), exp_xp, atol=1e-4, rtol=1e-4)
51+
52+
assert tau_xp.shape == tau_np.shape
53+
if not has_support_aspect64(tau_xp.sycl_device):
54+
assert tau_xp.dtype.kind == tau_np.dtype.kind
55+
else:
56+
assert tau_xp.dtype == tau_np.dtype
57+
58+
else: # mode == "r"
59+
r_xp = xp.linalg.qr(a_xp, mode="r")
60+
61+
# Use reduced QR as a reference:
62+
# reduced is validated via Q @ R == A
63+
exp_res = xp.linalg.qr(a_xp, mode="reduced")
64+
exp_r = exp_res.R
65+
assert xp.allclose(r_xp, exp_r, atol=1e-4, rtol=1e-4)
66+
67+
exp_xp = gram(a_xp, xp)
68+
69+
# Compare R^H @ R == A^H @ A
70+
assert xp.allclose(gram(r_xp, xp), exp_xp, atol=1e-4, rtol=1e-4)

dpnp/tests/test_linalg.py

Lines changed: 6 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
has_support_aspect64,
2626
numpy_version,
2727
)
28+
from .qr_helper import check_qr
2829
from .third_party.cupy import testing
2930

3031

@@ -3135,85 +3136,6 @@ def test_error(self):
31353136

31363137

31373138
class TestQr:
3138-
def gram(self, X, xp):
3139-
# Return Gram matrix: X^H @ X
3140-
return xp.conjugate(X).swapaxes(-1, -2) @ X
3141-
3142-
def get_R_from_raw(self, h, m, n, xp):
3143-
# Get reduced R from NumPy-style raw QR:
3144-
# R = triu((tril(h))^T), shape (..., k, n)
3145-
k = min(m, n)
3146-
Rt = xp.tril(h)
3147-
R = xp.swapaxes(Rt, -1, -2)
3148-
R = xp.triu(R[..., :m, :n])
3149-
3150-
return R[..., :k, :]
3151-
3152-
# QR is not unique:
3153-
# element-wise comparison with NumPy may differ by sign/phase.
3154-
# To verify correctness use mode-dependent functional checks:
3155-
# complete/reduced: check decomposition Q @ R = A
3156-
# raw/r: check invariant R^H @ R = A^H @ A
3157-
def check_qr(self, a_np, a_dp, mode):
3158-
if mode in ("complete", "reduced"):
3159-
res = dpnp.linalg.qr(a_dp, mode)
3160-
assert dpnp.allclose(res.Q @ res.R, a_dp, atol=1e-5)
3161-
3162-
# Since QR satisfies A = Q @ R with orthonormal Q (Q^H @ Q = I),
3163-
# validate correctness via the invariant R^H @ R == A^H @ A
3164-
# for raw/r modes
3165-
elif mode == "raw":
3166-
h_np, tau_np = numpy.linalg.qr(a_np, mode=mode)
3167-
h_dp, tau_dp = dpnp.linalg.qr(a_dp, mode=mode)
3168-
3169-
m, n = a_np.shape[-2], a_np.shape[-1]
3170-
Rraw_np = self.get_R_from_raw(h_np, m, n, numpy)
3171-
Rraw_dp = self.get_R_from_raw(h_dp, m, n, dpnp)
3172-
3173-
# Use reduced QR as a reference:
3174-
# reduced is validated via Q @ R == A
3175-
exp_res = dpnp.linalg.qr(a_dp, mode="reduced")
3176-
exp_R = exp_res.R
3177-
assert_allclose(Rraw_dp, exp_R, atol=1e-4, rtol=1e-4)
3178-
3179-
exp_dp = self.gram(a_dp, dpnp).astype(Rraw_dp.dtype)
3180-
exp_np = self.gram(a_np, numpy).astype(Rraw_np.dtype)
3181-
3182-
# compare R^H @ R == A^H @ A
3183-
assert_allclose(
3184-
self.gram(Rraw_dp, dpnp), exp_dp, atol=1e-4, rtol=1e-4
3185-
)
3186-
assert_allclose(
3187-
self.gram(Rraw_np, numpy), exp_np, atol=1e-4, rtol=1e-4
3188-
)
3189-
3190-
assert tau_dp.shape == tau_np.shape
3191-
if not has_support_aspect64(tau_dp.sycl_device):
3192-
if tau_np.dtype == numpy.float64:
3193-
tau_np = tau_np.astype("float32")
3194-
elif tau_np.dtype == numpy.complex128:
3195-
tau_np = tau_np.astype("complex64")
3196-
assert tau_dp.dtype == tau_np.dtype
3197-
3198-
else: # mode == "r"
3199-
R_np = numpy.linalg.qr(a_np, mode="r")
3200-
R_dp = dpnp.linalg.qr(a_dp, mode="r")
3201-
3202-
# Use reduced QR as a reference:
3203-
# reduced is validated via Q @ R == A
3204-
exp_res = dpnp.linalg.qr(a_dp, mode="reduced")
3205-
exp_R = exp_res.R
3206-
assert_allclose(R_dp, exp_R, atol=1e-4, rtol=1e-4)
3207-
3208-
exp_dp = self.gram(a_dp, dpnp)
3209-
exp_np = self.gram(a_np, numpy)
3210-
3211-
# compare R^H @ R == A^H @ A
3212-
assert_allclose(self.gram(R_dp, dpnp), exp_dp, atol=1e-4, rtol=1e-4)
3213-
assert_allclose(
3214-
self.gram(R_np, numpy), exp_np, atol=1e-4, rtol=1e-4
3215-
)
3216-
32173139
@pytest.mark.parametrize("dtype", get_float_complex_dtypes())
32183140
@pytest.mark.parametrize(
32193141
"shape",
@@ -3245,7 +3167,7 @@ def test_qr(self, dtype, shape, mode):
32453167
a = generate_random_numpy_array(shape, dtype, seed_value=None)
32463168
ia = dpnp.array(a, dtype=dtype)
32473169

3248-
self.check_qr(a, ia, mode)
3170+
check_qr(a, ia, mode, dpnp)
32493171

32503172
@pytest.mark.parametrize("dtype", get_float_complex_dtypes())
32513173
@pytest.mark.parametrize(
@@ -3258,7 +3180,7 @@ def test_qr_large(self, dtype, shape, mode):
32583180
a = generate_random_numpy_array(shape, dtype, seed_value=81)
32593181
ia = dpnp.array(a)
32603182

3261-
self.check_qr(a, ia, mode)
3183+
check_qr(a, ia, mode, dpnp)
32623184

32633185
@pytest.mark.parametrize("dtype", get_float_complex_dtypes())
32643186
@pytest.mark.parametrize(
@@ -3278,17 +3200,17 @@ def test_qr_empty(self, dtype, shape, mode):
32783200
a = numpy.empty(shape, dtype=dtype)
32793201
ia = dpnp.array(a)
32803202

3281-
self.check_qr(a, ia, mode)
3203+
check_qr(a, ia, mode, dpnp)
32823204

32833205
@pytest.mark.parametrize("mode", ["complete", "reduced", "r", "raw"])
32843206
def test_qr_strides(self, mode):
32853207
a = generate_random_numpy_array((5, 5))
32863208
ia = dpnp.array(a)
32873209

32883210
# positive strides
3289-
self.check_qr(a[::2, ::2], ia[::2, ::2], mode)
3211+
check_qr(a[::2, ::2], ia[::2, ::2], mode, dpnp)
32903212
# negative strides
3291-
self.check_qr(a[::-2, ::-2], ia[::-2, ::-2], mode)
3213+
check_qr(a[::-2, ::-2], ia[::-2, ::-2], mode, dpnp)
32923214

32933215
def test_qr_errors(self):
32943216
a_dp = dpnp.array([[1, 2], [3, 5]], dtype="float32")

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

Lines changed: 11 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
from dpnp.tests.helper import (
1515
has_support_aspect64,
1616
)
17+
from dpnp.tests.qr_helper import check_qr
1718
from dpnp.tests.third_party.cupy import testing
1819
from dpnp.tests.third_party.cupy.testing import _condition
1920

@@ -167,19 +168,6 @@ def test_decomposition(self, dtype):
167168
)
168169
)
169170
class TestQRDecomposition(unittest.TestCase):
170-
def _gram(self, x, xp):
171-
# Gram matrix: X^H @ X
172-
return xp.conjugate(x).swapaxes(-1, -2) @ x
173-
174-
def _get_R_from_raw(self, h, m, n, xp):
175-
# Get reduced R from NumPy-style raw QR:
176-
# R = triu((tril(h))^T), shape (..., k, n)
177-
k = min(m, n)
178-
Rt = xp.tril(h)
179-
R = xp.swapaxes(Rt, -1, -2)
180-
R = xp.triu(R[..., :m, :n])
181-
return R[..., :k, :]
182-
183171
@testing.for_dtypes("fdFD")
184172
def check_mode(self, array, mode, dtype):
185173
# if runtime.is_hip and driver.get_build_version() < 307:
@@ -188,13 +176,20 @@ def check_mode(self, array, mode, dtype):
188176

189177
a_cpu = numpy.asarray(array, dtype=dtype)
190178
a_gpu = cupy.asarray(array, dtype=dtype)
191-
result_gpu = cupy.linalg.qr(a_gpu, mode=mode)
179+
# QR is not unique:
180+
# element-wise comparison with NumPy may differ by sign/phase.
181+
# To verify correctness use mode-dependent functional checks:
182+
# complete/reduced: check decomposition Q @ R = A
183+
# raw/r: check invariant R^H @ R = A^H @ A
184+
185+
# result_gpu = cupy.linalg.qr(a_gpu, mode=mode)
192186
if (
193187
mode != "raw"
194188
or numpy.lib.NumpyVersion(numpy.__version__) >= "1.22.0rc1"
195189
):
196-
result_cpu = numpy.linalg.qr(a_cpu, mode=mode)
197-
self._check_result(result_cpu, result_gpu, a_cpu, a_gpu, mode)
190+
# result_cpu = numpy.linalg.qr(a_cpu, mode=mode)
191+
# self._check_result(result_cpu, result_gpu, a_gpu, mode)
192+
check_qr(a_cpu, a_gpu, mode, cupy)
198193

199194
# def _check_result(self, result_cpu, result_gpu):
200195
# if isinstance(result_cpu, tuple):
@@ -205,54 +200,6 @@ def check_mode(self, array, mode, dtype):
205200
# assert result_cpu.dtype == result_gpu.dtype
206201
# testing.assert_allclose(result_cpu, result_gpu, atol=1e-4)
207202

208-
# QR is not unique:
209-
# element-wise comparison with NumPy may differ by sign/phase.
210-
# To verify correctness use mode-dependent functional checks:
211-
# complete/reduced: check decomposition Q @ R = A
212-
# raw/r: check invariant R^H @ R = A^H @ A
213-
def _check_result(self, result_cpu, result_gpu, a_cpu, a_gpu, mode):
214-
215-
if mode in ("complete", "reduced"):
216-
q_gpu, r_gpu = result_gpu
217-
testing.assert_allclose(q_gpu @ r_gpu, a_gpu, atol=1e-4)
218-
219-
elif mode == "raw":
220-
h_gpu, tau_gpu = result_gpu
221-
h_cpu, tau_cpu = result_cpu
222-
m, n = a_gpu.shape[-2], a_gpu.shape[-1]
223-
r_gpu = self._get_R_from_raw(h_gpu, m, n, cupy)
224-
r_cpu = self._get_R_from_raw(h_cpu, m, n, numpy)
225-
226-
exp_gpu = self._gram(a_gpu, cupy)
227-
exp_cpu = self._gram(a_cpu, numpy)
228-
229-
testing.assert_allclose(
230-
self._gram(r_gpu, cupy), exp_gpu, atol=1e-4, rtol=1e-4
231-
)
232-
testing.assert_allclose(
233-
self._gram(r_cpu, numpy), exp_cpu, atol=1e-4, rtol=1e-4
234-
)
235-
236-
assert tau_gpu.shape == tau_cpu.shape
237-
if not has_support_aspect64(tau_gpu.sycl_device):
238-
if tau_cpu.dtype == numpy.float64:
239-
tau_cpu = tau_cpu.astype("float32")
240-
elif tau_cpu.dtype == numpy.complex128:
241-
tau_cpu = tau_cpu.astype("complex64")
242-
assert tau_gpu.dtype == tau_cpu.dtype
243-
244-
else: # mode == "r"
245-
r_gpu = result_gpu
246-
r_cpu = result_cpu
247-
exp_gpu = self._gram(a_gpu, cupy)
248-
exp_cpu = self._gram(a_cpu, numpy)
249-
testing.assert_allclose(
250-
self._gram(r_gpu, cupy), exp_gpu, atol=1e-4, rtol=1e-4
251-
)
252-
testing.assert_allclose(
253-
self._gram(r_cpu, numpy), exp_cpu, atol=1e-4, rtol=1e-4
254-
)
255-
256203
@testing.fix_random()
257204
@_condition.repeat(3, 10)
258205
def test_mode(self):

0 commit comments

Comments
 (0)