diff --git a/dpnp/tests/test_product.py b/dpnp/tests/test_product.py
index 9d43b6453e0..88627b7d2cc 100644
--- a/dpnp/tests/test_product.py
+++ b/dpnp/tests/test_product.py
@@ -887,9 +887,14 @@ def test_order(self, dtype, order1, order2, order, shape1, shape2):
         ids=["-2", "2", "(-2, 2)", "(2, -2)"],
     )
     def test_strided1(self, dtype, stride):
-        # TODO: enable back when the root cause is identified
-        # for dim in [1, 2, 3, 4]:
-        for dim in [1, 2, 3]:
+        # dpnp copies the strided input into a c-contiguous array and runs the
+        # batched product through oneMKL gemm_batch, while NumPy uses OpenBLAS.
+        # The two backends accumulate in a different order, so the float32
+        # result deviates at the noise floor (~2e-5 for a length-20 contraction
+        # with dim=4), which exceeds the default tolerance. Integer input is
+        # exact, so the factor only affects the float32 case.
+        factor = 24 if dtype == numpy.float32 else 16
+        for dim in [1, 2, 3, 4]:
             shape = tuple(20 for _ in range(dim))
             A = generate_random_numpy_array(shape, dtype)
             iA = dpnp.array(A)
@@ -900,7 +905,7 @@ def test_strided1(self, dtype, stride):
             # the 2D base is not c-contiguous nor f-contigous
             result = dpnp.matmul(ia, ia)
             expected = numpy.matmul(a, a)
-            assert_dtype_allclose(result, expected, factor=16)
+            assert_dtype_allclose(result, expected, factor=factor)
 
             OUT = numpy.empty(shape, dtype=result.dtype)
             out = OUT[slices]
@@ -909,7 +914,7 @@ def test_strided1(self, dtype, stride):
             result = dpnp.matmul(ia, ia, out=iout)
             assert result is iout
             expected = numpy.matmul(a, a, out=out)
-            assert_dtype_allclose(result, expected, factor=16)
+            assert_dtype_allclose(result, expected, factor=factor)
 
     @pytest.mark.parametrize("dtype", _selected_dtypes)
     @pytest.mark.parametrize(
@@ -1545,7 +1550,9 @@ def test_axes(self, axes):
         result = dpnp.matvec(ia, ib, axes=axes)
         expected = numpy.matvec(a, b, axes=axes)
 
-        # TODO: check if failing with newer NumPy
+        # dpnp uses oneMKL gemm_batch while NumPy uses OpenBLAS gemv, so the
+        # summation order differs and the float64 result deviates at the noise
+        # floor (~1e-14), which exceeds the default tolerance
         assert_dtype_allclose(result, expected, factor=40)
 
     @pytest.mark.parametrize("xp", [numpy, dpnp])