add reconstruction tests for lu_factor
@@ -11,7 +11,7 @@ if cupyx.scipy._scipy_available:
@testing.gpu
@testing.parameterize(*testing.product({
- 'shape': [(1, 1), (2, 2), (3, 3), (5, 5)],
+ 'shape': [(1, 1), (2, 2), (3, 3), (5, 5), (1, 5), (5, 1), (2, 5), (5, 2)],
}))
@testing.fix_random()
@unittest.skipUnless(
@@ -21,6 +21,9 @@ class TestLUFactor(unittest.TestCase):
@testing.for_float_dtypes(no_float16=True)
def test_lu_factor(self, dtype):
+ if self.shape[0] != self.shape[1]:
+ # skip non-square tests since scipy.lu_factor requires square
+ return
array = numpy.random.randn(*self.shape)
a_cpu = numpy.asarray(array, dtype=dtype)
a_gpu = cupy.asarray(array, dtype=dtype)
@@ -32,6 +35,34 @@ class TestLUFactor(unittest.TestCase):
cupy.testing.assert_allclose(result_cpu[0], result_gpu[0], atol=1e-5)
cupy.testing.assert_array_equal(result_cpu[1], result_gpu[1])
+ @testing.for_float_dtypes(no_float16=True)
+ def test_lu_factor_reconstruction(self, dtype):
+ m, n = self.shape
+ array = cupy.random.randn(m, n, dtype=dtype)
+ lu, piv = cupyx.scipy.linalg.lu_factor(array)
+ # extract ``L`` and ``U`` from ``lu``
+ L = cupy.tril(lu, k=-1)
+ cupy.fill_diagonal(L, 1.)
+ if m < n:
+ L = L[:, :m]
+ U = cupy.triu(lu)
+ if m > n:
+ U = U[:n, :]
+ # apply pivot (on CPU since slaswp is not available in cupy)
+ piv = cupy.asnumpy(piv)
+ rows = numpy.arange(array.shape[0])
+ for i, row in enumerate(piv):
+ if i != row:
+ rows[i], rows[row] = rows[row], rows[i]
+ # revert pivot
+ reversed_piv = numpy.empty_like(rows)
+ reversed_piv[rows] = numpy.arange(rows.size)
+ # swap L
+ L = L[reversed_piv]
+ # check that reconstruction is close to original
+ reconstructed = L.dot(U)
+ cupy.testing.assert_allclose(reconstructed, array, atol=1e-5)
+
@testing.gpu
@testing.parameterize(*testing.product({
