ValueError in PLSRegression.fit() with zero-variance predictor

### Describe the bug

Related: https://github.com/scipy/scipy/commit/5bc3d8814d566ef328f41cfa69ccd797c68b0d02

When fitting a PLSRegression model, if the input array X contains a feature with zero variance (i.e., a constant column), the fit method raises a ValueError: illegal value in 4th argument of internal gesdd.

This results in a division by zero when a predictor has no variance, creating NaN values likely in the intermediate matrices. These NaN values are then passed to the SciPy function, which in turn calls the LAPACK gesdd routine for SVD, causing it to crash.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.cross_decomposition import PLSRegression

n_samples = 20
y = np.arange(n_samples, dtype=float)

# This feature has zero variance.
X = np.ones((n_samples, 1))

# This will raise the error.
pls = PLSRegression(n_components=1)
pls.fit(X, y)
```

### Expected Results

The model should either fit successfully (e.g., perhaps assigning a zero weight to the zero-variance feature) or raise a more informative ValueError indicating that a predictor has zero variance.

### Actual Results

We get "ValueError: illegal value in 4th argument of internal gesdd"

```python
/home/user/.local/lib/python3.13/site-packages/sklearn/cross_decomposition/_pls.py:99: RuntimeWarning: invalid value encountered in divide
  y_weights = np.dot(Y.T, x_score) / np.dot(x_score.T, x_score)
/home/user/.local/lib/python3.13/site-packages/sklearn/cross_decomposition/_pls.py:368: RuntimeWarning: invalid value encountered in divide
  x_loadings = np.dot(x_scores, Xk) / np.dot(x_scores, x_scores)
/home/user/.local/lib/python3.13/site-packages/sklearn/cross_decomposition/_pls.py:377: RuntimeWarning: invalid value encountered in divide
  y_loadings = np.dot(x_scores, yk) / np.dot(x_scores, x_scores)
Traceback (most recent call last):
  File "/home/user/agents/test/f.py", line 14, in <module>
    pls.fit(X, y)
    ~~~~~~~^^^^^^
  File "/home/user/.local/lib/python3.13/site-packages/sklearn/cross_decomposition/_pls.py", line 719, in fit
    super().fit(X, y)
    ~~~~~~~~~~~^^^^^^
  File "/home/user/.local/lib/python3.13/site-packages/sklearn/base.py", line 1389, in wrapper
    return fit_method(estimator, *args, **kwargs)
  File "/home/user/.local/lib/python3.13/site-packages/sklearn/cross_decomposition/_pls.py", line 396, in fit
    pinv2(np.dot(self.x_loadings_.T, self.x_weights_), check_finite=False),
    ~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/.local/lib/python3.13/site-packages/scipy/linalg/_basic.py", line 1607, in pinv
    u, s, vh = _decomp_svd.svd(a, full_matrices=False, check_finite=False)                                                     ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/.local/lib/python3.13/site-packages/scipy/linalg/_decomp_svd.py", line 168, in svd               raise ValueError('illegal value in %dth argument of internal gesdd'                                                              % -info)
ValueError: illegal value in 4th argument of internal gesdd
```

### Versions

```shell
System:                                                     python: 3.13.5 (main, Jun 12 2025, 00:00:00) [GCC 14.3.1 20250523 (Red Hat 14.3.1-1)]
executable: /usr/bin/python3                               machine: Linux-6.14.8-400.asahi.fc41.aarch64+16k-aarch64-with-glibc2.40                                                                                              Python dependencies:
      sklearn: 1.6.1
          pip: 25.1.1                                      setuptools: 69.2.0
        numpy: 1.26.4                                           scipy: 1.15.3
       Cython: None                                            pandas: 2.2.3
   matplotlib: 3.10.3
       joblib: 1.5.0                                    threadpoolctl: 3.6.0
                                                        Built with OpenMP: True
threadpoolctl info:                                            user_api: blas
   internal_api: openblas                                   num_threads: 1
         prefix: libopenblas                                   filepath: /usr/lib64/libopenblas-r0.3.26.so
        version: 0.3.26                                 threading_layer: disabled
   architecture: neoversen1
       user_api: openmp                                    internal_api: openmp
    num_threads: 8                                               prefix: libgomp
       filepath: /usr/lib64/libgomp.so.1.0.0                    version: None
                                                               user_api: blas
   internal_api: openblas
    num_threads: 8                                               prefix: libscipy_openblas
       filepath: /home/user/.local/lib/python3.13/site-packages/scipy.libs/libscipy_openblas-9778f98e.so
        version: 0.3.28                                 threading_layer: pthreads
   architecture: neoversen1

       user_api: openmp                                    internal_api: openmp
    num_threads: 8                                               prefix: libgomp
       filepath: /home/user/.local/lib/python3.13/site-packages/scikit_learn.libs/libgomp-d22c30c5.so.1.0.0             version: None
```

COMMENT:
PLS (NIPALS) extracts one component at a time and deflates the data after each extraction, so the original data could have good variance and rank, but the residual matrix might be numerically zero after the extraction of the informative directions, like if we request many components (more than the numerical rank). (When the residual matrix is near-zero, the the next score is the zero vector, so we end up dividing by zero.)

It seems like if the data are trivial/rank-0 case at the start, then ValueError is good and expected, but because PLS is iterative, we can get the problem from non-degenerate data.

I originally encountered this on non-trivial genomic data (specifically, here: https://github.com/SauersML/ferromic/blob/main/stats/linked.py). I was assessing how different component numbers perform, and suspected at first that I was requesting component numbers higher than the numerical rank, leading to the observed problem.

Trivial replication:
Code:
```
import numpy as np
from sklearn.cross_decomposition import PLSRegression

n_samples = 20
y = np.arange(n_samples, dtype=float)

# This feature has zero variance.
X = np.ones((n_samples, 1))

# This will raise the error.
pls = PLSRegression(n_components=1)
pls.fit(X, y)
```

Result (default Google Colab):
```
/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py:99: RuntimeWarning: invalid value encountered in divide
  y_weights = np.dot(Y.T, x_score) / np.dot(x_score.T, x_score)
/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py:368: RuntimeWarning: invalid value encountered in divide
  x_loadings = np.dot(x_scores, Xk) / np.dot(x_scores, x_scores)
/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py:377: RuntimeWarning: invalid value encountered in divide
  y_loadings = np.dot(x_scores, yk) / np.dot(x_scores, x_scores)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[/tmp/ipython-input-1137519010.py](https://localhost:8080/#) in <cell line: 0>()
     10 # This will raise the error.
     11 pls = PLSRegression(n_components=1)
---> 12 pls.fit(X, y)

6 frames
[/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py](https://localhost:8080/#) in fit(self, X, y, Y)
    717         y = _deprecate_Y_when_required(y, Y)
    718 
--> 719         super().fit(X, y)
    720         # expose the fitted attributes `x_scores_` and `y_scores_`
    721         self.x_scores_ = self._x_scores

[/usr/local/lib/python3.12/dist-packages/sklearn/base.py](https://localhost:8080/#) in wrapper(estimator, *args, **kwargs)
   1387                 )
   1388             ):
-> 1389                 return fit_method(estimator, *args, **kwargs)
   1390 
   1391         return wrapper

[/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py](https://localhost:8080/#) in fit(self, X, y, Y)
    394         self.x_rotations_ = np.dot(
    395             self.x_weights_,
--> 396             pinv2(np.dot(self.x_loadings_.T, self.x_weights_), check_finite=False),
    397         )
    398         self.y_rotations_ = np.dot(

[/usr/local/lib/python3.12/dist-packages/scipy/_lib/_util.py](https://localhost:8080/#) in wrapper(*args, **kwargs)
   1231             # Early exit if call is not batched
   1232             if not any(batch_shapes):
-> 1233                 return f(*arrays, *other_args, **kwargs)
   1234 
   1235             # Determine broadcasted batch shape

[/usr/local/lib/python3.12/dist-packages/scipy/linalg/_basic.py](https://localhost:8080/#) in pinv(a, atol, rtol, return_rank, check_finite)
   1631     """
   1632     a = _asarray_validated(a, check_finite=check_finite)
-> 1633     u, s, vh = _decomp_svd.svd(a, full_matrices=False, check_finite=False)
   1634     t = u.dtype.char.lower()
   1635     maxS = np.max(s, initial=0.)

[/usr/local/lib/python3.12/dist-packages/scipy/_lib/_util.py](https://localhost:8080/#) in wrapper(*args, **kwargs)
   1231             # Early exit if call is not batched
   1232             if not any(batch_shapes):
-> 1233                 return f(*arrays, *other_args, **kwargs)
   1234 
   1235             # Determine broadcasted batch shape

[/usr/local/lib/python3.12/dist-packages/scipy/linalg/_decomp_svd.py](https://localhost:8080/#) in svd(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)
    172         if lapack_driver == "gesdd" and info == -4:
    173             msg = "A has a NaN entry"
--> 174             raise ValueError(msg)
    175         raise ValueError(f'illegal value in {-info}th argument of internal gesdd')
    176     if compute_uv:

ValueError: A has a NaN entry
```

I've tried a few ways and cannot replicate this in a more non-trivial fashion. I might go back to my genomic data and analyze what makes the cases that fail with this error different from the cases that succeed.

Example code:
```
import numpy as np
from sklearn.cross_decomposition import PLSRegression

rng = np.random.default_rng(0)
n = 300
Z = np.stack([rng.normal(size=n) for _ in range(2)], axis=1)
W = rng.normal(size=(2, 8))
X = Z @ W
y = Z @ rng.normal(size=2)

print("X shape:", X.shape, " y shape:", y.shape)
print("matrix_rank(X):", np.linalg.matrix_rank(X))

pls = PLSRegression(n_components=3)
pls.fit(X, y)
y_pred = pls.predict(X)
print("First 10 predictions:", y_pred[:10].ravel())
```

Results:
```
X shape: (300, 8)  y shape: (300,)
matrix_rank(X): 2
First 10 predictions: [-1.74180058 -0.44105344 -2.21793852  4.25196395  0.98597734 -0.8453282
 -2.98242326 -1.54044121  1.6381953   2.57236178]
/usr/local/lib/python3.12/dist-packages/sklearn/cross_decomposition/_pls.py:308: UserWarning: y residual is constant at iteration 2
  warnings.warn(f"y residual is constant at iteration {k}")
```

I also saw this warning in real data; however, this occured at different times and did not cause a crash.

Labels: Bug