Repository: scikit-learn/scikit-learn
Pull request: FIX Add input array check to `randomized_range_finder`
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #30736. See also #30737.

#### What does this implement/fix? Explain your changes.

Adds a check of the input array to `randomized_range_finder` so that complex valued inputs are rejected. As a result, complex inputs are also rejected from `randomized_svd`, which calls `randomized_range_finder`.


<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
https://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
`randomized_svd` incorrect for complex valued matrices
### Describe the bug

The `randomized_svd` utility function accepts complex valued inputs without error, but the result is inconsistent with `scipy.linalg.svd`.

### Steps/Code to Reproduce

```python
import numpy as np
from scipy import linalg
from sklearn.utils.extmath import randomized_svd

rng = np.random.RandomState(42)
X = rng.randn(100, 20) + 1j * rng.randn(100, 20)

_, s, _ = linalg.svd(X)
_, s2, _ = randomized_svd(X, n_components=5)

print("s:", s[:5])
print("s2:", s2[:5])
```

### Expected Results

I expected the singular values to be numerically close.

### Actual Results

```
s: [19.81481515 18.69019042 17.62107998 17.23689681 16.3148512 ]
s2: [11.25690754  9.97157079  9.01542947  8.06160863  7.54068744]
```

### Versions

```shell
System:
    python: 3.11.4 (main, Jul  5 2023, 08:40:20) [Clang 14.0.6 ]
executable: /Users/clane/miniconda3/bin/python
   machine: macOS-13.7-arm64-arm-64bit

Python dependencies:
      sklearn: 1.7.dev0
          pip: 25.0
   setuptools: 65.5.0
        numpy: 2.2.2
        scipy: 1.15.1
       Cython: 3.0.11
       pandas: 2.2.3
   matplotlib: 3.10.0
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 8
         prefix: libscipy_openblas
       filepath: /Users/clane/Projects/misc/scikit-learn/.venv/lib/python3.11/site-packages/numpy/.dylibs/libscipy_openblas64_.dylib
        version: 0.3.28
threading_layer: pthreads
   architecture: neoversen1

       user_api: blas
   internal_api: openblas
    num_threads: 8
         prefix: libscipy_openblas
       filepath: /Users/clane/Projects/misc/scikit-learn/.venv/lib/python3.11/site-packages/scipy/.dylibs/libscipy_openblas.dylib
        version: 0.3.28
threading_layer: pthreads
   architecture: neoversen1

       user_api: openmp
   internal_api: openmp
    num_threads: 8
         prefix: libomp
       filepath: /opt/homebrew/Cellar/libomp/19.1.3/lib/libomp.dylib
        version: None
```

Hints from issue discussion:

Commit message: iter

Patch to review:
diff --git a/doc/whats_new/upcoming_changes/array-api/30819.feature.rst b/doc/whats_new/upcoming_changes/array-api/30819.feature.rst
new file mode 100644
index 0000000000000..40c4fc306a90c
--- /dev/null
+++ b/doc/whats_new/upcoming_changes/array-api/30819.feature.rst
@@ -0,0 +1,2 @@
+- :func:`sklearn.utils.extmath.randomized_svd` now support Array API compatible inputs.
+  By :user:`Connor Lane <clane9>`.
\ No newline at end of file
diff --git a/doc/whats_new/upcoming_changes/sklearn.utils/30819.fix.rst b/doc/whats_new/upcoming_changes/sklearn.utils/30819.fix.rst
new file mode 100644
index 0000000000000..81c7564023ac1
--- /dev/null
+++ b/doc/whats_new/upcoming_changes/sklearn.utils/30819.fix.rst
@@ -0,0 +1,4 @@
+- :func:`utils.extmath.randomized_svd` and :func:`utils.extmath.randomized_range_finder`
+  now validate their input array to fail early with an informative error message on
+  invalid input.
+  By :user:`Connor Lane <clane9>`.
diff --git a/sklearn/utils/extmath.py b/sklearn/utils/extmath.py
index b4af090344d74..9a189f69bb03c 100644
--- a/sklearn/utils/extmath.py
+++ b/sklearn/utils/extmath.py
@@ -219,7 +219,7 @@ def randomized_range_finder(
 
     Parameters
     ----------
-    A : 2D array
+    A : {array-like, sparse matrix} of shape (n_samples, n_features)
         The input data matrix.
 
     size : int
@@ -246,9 +246,9 @@ def randomized_range_finder(
 
     Returns
     -------
-    Q : ndarray
-        A (size x size) projection matrix, the range of which
-        approximates well the range of the input matrix A.
+    Q : ndarray of shape (size, size)
+        A projection matrix, the range of which approximates well the range of the
+        input matrix A.
 
     Notes
     -----
@@ -274,6 +274,7 @@ def randomized_range_finder(
            [-0.82..., -0.38...]])
     """
     xp, is_array_api_compliant = get_namespace(A)
+    A = check_array(A, accept_sparse=True)
     random_state = check_random_state(random_state)
 
     # Generating normal random vectors with shape: (A.shape[1], size)
@@ -344,7 +345,7 @@ def randomized_range_finder(
 
 @validate_params(
     {
-        "M": [np.ndarray, "sparse matrix"],
+        "M": ["array-like", "sparse matrix"],
         "n_components": [Interval(Integral, 1, None, closed="left")],
         "n_oversamples": [Interval(Integral, 0, None, closed="left")],
         "n_iter": [Interval(Integral, 0, None, closed="left"), StrOptions({"auto"})],
@@ -381,7 +382,7 @@ def randomized_svd(
 
     Parameters
     ----------
-    M : {ndarray, sparse matrix}
+    M : {array-like, sparse matrix} of shape (n_samples, n_features)
         Matrix to decompose.
 
     n_components : int
@@ -499,6 +500,9 @@ def randomized_svd(
     >>> U.shape, s.shape, Vh.shape
     ((3, 2), (2,), (2, 4))
     """
+    xp, is_array_api_compliant = get_namespace(M)
+    M = check_array(M, accept_sparse=True)
+
     if sparse.issparse(M) and M.format in ("lil", "dok"):
         warnings.warn(
             "Calculating SVD of a {} is expensive. "
@@ -533,7 +537,6 @@ def randomized_svd(
     B = Q.T @ M
 
     # compute the SVD on the thin matrix: (k + p) wide
-    xp, is_array_api_compliant = get_namespace(B)
     if is_array_api_compliant:
         Uhat, s, Vt = xp.linalg.svd(B, full_matrices=False)
     else:
diff --git a/sklearn/utils/tests/test_extmath.py b/sklearn/utils/tests/test_extmath.py
index 74cb47388692f..907de11702af2 100644
--- a/sklearn/utils/tests/test_extmath.py
+++ b/sklearn/utils/tests/test_extmath.py
@@ -9,10 +9,18 @@
 from scipy.linalg import eigh
 from scipy.sparse.linalg import eigsh
 
+from sklearn import config_context
 from sklearn.datasets import make_low_rank_matrix, make_sparse_spd_matrix
 from sklearn.utils import gen_batches
 from sklearn.utils._arpack import _init_arpack_v0
+from sklearn.utils._array_api import (
+    _convert_to_numpy,
+    _get_namespace_device_dtype_ids,
+    get_namespace,
+    yield_namespace_device_dtype_combinations,
+)
 from sklearn.utils._testing import (
+    _array_api_for_tests,
     assert_allclose,
     assert_allclose_dense_sparse,
     assert_almost_equal,
@@ -28,6 +36,7 @@
     _safe_accumulator_op,
     cartesian,
     density,
+    randomized_range_finder,
     randomized_svd,
     row_norms,
     safe_sparse_dot,
@@ -1060,3 +1069,53 @@ def test_approximate_mode():
     # 25% * 99.000 = 24.750
     # 25% *  1.000 =    250
     assert_array_equal(ret, [24750, 250])
+
+
+@pytest.mark.parametrize(
+    "array_namespace, device, dtype",
+    yield_namespace_device_dtype_combinations(),
+    ids=_get_namespace_device_dtype_ids,
+)
+def test_randomized_svd_array_api_compliance(array_namespace, device, dtype):
+    xp = _array_api_for_tests(array_namespace, device)
+
+    rng = np.random.RandomState(0)
+    X = rng.normal(size=(30, 10)).astype(dtype)
+    X_xp = xp.asarray(X, device=device)
+    n_components = 5
+    atol = 1e-5 if dtype == "float32" else 0
+
+    with config_context(array_api_dispatch=True):
+        u_np, s_np, vt_np = randomized_svd(X, n_components, random_state=0)
+        u_xp, s_xp, vt_xp = randomized_svd(X_xp, n_components, random_state=0)
+
+        assert get_namespace(u_xp)[0].__name__ == xp.__name__
+        assert get_namespace(s_xp)[0].__name__ == xp.__name__
+        assert get_namespace(vt_xp)[0].__name__ == xp.__name__
+
+        assert_allclose(_convert_to_numpy(u_xp, xp), u_np, atol=atol)
+        assert_allclose(_convert_to_numpy(s_xp, xp), s_np, atol=atol)
+        assert_allclose(_convert_to_numpy(vt_xp, xp), vt_np, atol=atol)
+
+
+@pytest.mark.parametrize(
+    "array_namespace, device, dtype",
+    yield_namespace_device_dtype_combinations(),
+    ids=_get_namespace_device_dtype_ids,
+)
+def test_randomized_range_finder_array_api_compliance(array_namespace, device, dtype):
+    xp = _array_api_for_tests(array_namespace, device)
+
+    rng = np.random.RandomState(0)
+    X = rng.normal(size=(30, 10)).astype(dtype)
+    X_xp = xp.asarray(X, device=device)
+    size = 5
+    n_iter = 10
+    atol = 1e-5 if dtype == "float32" else 0
+
+    with config_context(array_api_dispatch=True):
+        Q_np = randomized_range_finder(X, size=size, n_iter=n_iter, random_state=0)
+        Q_xp = randomized_range_finder(X_xp, size=size, n_iter=n_iter, random_state=0)
+
+        assert get_namespace(Q_xp)[0].__name__ == xp.__name__
+        assert_allclose(_convert_to_numpy(Q_xp, xp), Q_np, atol=atol)


Code review comments:
diff_hunk: @@ -274,6 +274,7 @@ def randomized_range_finder(
            [-0.82..., -0.38...]])
     """
     xp, is_array_api_compliant = get_namespace(A)
+    A = check_array(A, accept_sparse=True)
path: sklearn/utils/extmath.py
text: @user1:
Similar comment as for `_randomized_svd`: we need to introduce a private `_randomized_range_finder` to remove the overhead of multiple input checks.