Instability in fastica for float32

Originally seen in https://github.com/scikit-learn/scikit-learn/issues/24131#issuecomment-1208091119 for more details on a single random seed and only for Atlas but I have managed to reproduce with OpenBLAS see below.

For now we have a work-around for the CI in #24198 but it is still an issue and in an ideal world, someone would investigate and try to see whether it is fixable.

Other possibly related issues and PRs are mentioned in https://github.com/scikit-learn/scikit-learn/pull/24198#issuecomment-1252616127: https://github.com/scikit-learn/scikit-learn/issues/2735 and https://github.com/scikit-learn/scikit-learn/pull/2738.

Here is a snippet that reproduces the same problem on OpenBLAS:

```py
import numpy as np
from scipy import stats
from sklearn.decomposition import fastica
import sys

global_random_seed = 20
global_dtype = np.float32


def center_and_norm(x, axis=-1):
    x = np.rollaxis(x, axis)
    x -= x.mean(axis=0)
    x /= x.std(axis=0)


rng = np.random.RandomState(global_random_seed)
n_samples = 1000
# Generate two sources:
s1 = (2 * np.sin(np.linspace(0, 100, n_samples)) > 0) - 1
s2 = stats.t.rvs(1, size=n_samples, random_state=global_random_seed)
s = np.c_[s1, s2].T
center_and_norm(s)
s = s.astype(global_dtype)
s1, s2 = s

# Mixing angle
phi = 0.6
mixing = np.array([[np.cos(phi), np.sin(phi)], [np.sin(phi), -np.cos(phi)]])
mixing = mixing.astype(global_dtype)
m = np.dot(mixing, s)
center_and_norm(m)

algo = 'deflation'
nl = 'logcosh'
whiten = 'arbitrary-variance'

problematic_random_state = 13441

k_, mixing_, s_ = fastica(
    m.T, fun=nl, whiten=whiten, algorithm=algo, random_state=problematic_random_state
)
```

You get a warning with a division by zero:
```
/home/lesteve/dev/scikit-learn/sklearn/decomposition/_fastica.py:89: RuntimeWarning: invalid value encountered in divide
  w1 /= np.sqrt((w1**2).sum())
```

And then the following traceback because NaNs are passed into a BLAS routine eventually:

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 39
     35 whiten = 'arbitrary-variance'
     37 problematic_random_state = 13441 # both for openblas and atlas (also numpy < and >= 1.23)
---> 39 k_, mixing_, s_ = fastica(
     40     m.T, fun=nl, whiten=whiten, algorithm=algo, random_state=problematic_random_state
     41 )

File ~/dev/scikit-learn/sklearn/decomposition/_fastica.py:322, in fastica(X, n_components, algorithm, whiten, fun, fun_args, max_iter, tol, w_init, whiten_solver, random_state, return_X_mean, compute_sources, return_n_iter)
    174 """Perform Fast Independent Component Analysis.
    175 
    176 The implementation is based on [1]_.
   (...)
    308        pp. 411-430.
    309 """
    310 est = FastICA(
    311     n_components=n_components,
    312     algorithm=algorithm,
   (...)
    320     random_state=random_state,
    321 )
--> 322 S = est._fit_transform(X, compute_sources=compute_sources)
    324 if est._whiten in ["unit-variance", "arbitrary-variance"]:
    325     K = est.whitening_

File ~/dev/scikit-learn/sklearn/decomposition/_fastica.py:683, in FastICA._fit_transform(self, X, compute_sources)
    680 else:
    681     self.components_ = W
--> 683 self.mixing_ = linalg.pinv(self.components_, check_finite=False)
    684 self._unmixing = W
    686 return S

File ~/miniconda3/lib/python3.9/site-packages/scipy/linalg/_basic.py:1304, in pinv(a, atol, rtol, return_rank, check_finite, cond, rcond)
   1231 """
   1232 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
   1233 
   (...)
   1301 
   1302 """
   1303 a = _asarray_validated(a, check_finite=check_finite)
-> 1304 u, s, vh = _decomp_svd.svd(a, full_matrices=False, check_finite=False)
   1305 t = u.dtype.char.lower()
   1306 maxS = np.max(s)

File ~/miniconda3/lib/python3.9/site-packages/scipy/linalg/_decomp_svd.py:133, in svd(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)
    131     raise LinAlgError("SVD did not converge")
    132 if info < 0:
--> 133     raise ValueError('illegal value in %dth argument of internal gesdd'
    134                      % -info)
    135 if compute_uv:
    136     return u, s, v

ValueError: illegal value in 4th argument of internal gesdd
```

COMMENT:
since `float32` collapses under some seeds, maybe FastICA could normalize with an epsilon threshold or restart if the weight norm falls below tolerance. Casting internal steps to float64 might also stabilize whitening/SVD. Curious whether maintainers would accept such changes

Labels: Bug, module:decomposition