⚠️ CI failed on Linux.ubuntu_atlas ⚠️

**CI is still failing on [Linux.ubuntu_atlas](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=46059&view=logs&j=0a287ed6-22f4-5cb4-88b1-d5fcdc4d8b7e)** (Aug 30, 2022)
- Test Collection Failure

COMMENT:
Seems like a real issue for global_random_seed=20. On my Ubuntu 20.04 machine (making sure atlas is used by chosing atlas with `sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu`) I get an error rather than a crash:

```
SKLEARN_RUN_FLOAT32_TESTS=1 SKLEARN_TESTS_GLOBAL_RANDOM_SEED=20 pytest sklearn/decomposition/tests/test_fastica.py -k test_fastica_simple -k False 
```

With an error like this which is due to some linalg operation being called on a array with NaNs:
```
ValueError: On entry to SLASCL parameter number 4 had an illegal value
```

<details>
<summary>full pytest failure info</summary>

```
_______________________________________ test_fastica_simple[20-float32-False] _______________________________________

add_noise = False, global_random_seed = 20, global_dtype = <class 'numpy.float32'>

    @pytest.mark.filterwarnings(
        "ignore:Starting in v1.3, whiten='unit-variance' will be used by default."
    )
    @pytest.mark.parametrize("add_noise", [True, False])
    def test_fastica_simple(add_noise, global_random_seed, global_dtype):
        # Test the FastICA algorithm on very simple data.
        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)
    
        if add_noise:
            m += 0.1 * rng.randn(2, 1000)
    
        center_and_norm(m)
    
        # function as fun arg
        def g_test(x):
            return x**3, (3 * x**2).mean(axis=-1)
    
        algos = ["parallel", "deflation"]
        nls = ["logcosh", "exp", "cube", g_test]
        whitening = ["arbitrary-variance", "unit-variance", False]
        for algo, nl, whiten in itertools.product(algos, nls, whitening):
            if whiten:
                kwargs = dict(X=m.T, fun=nl, whiten=whiten, algorithm=algo, random_state=rng)
                import copy
                kwargs_copy = copy.deepcopy(kwargs)
                try:
>                   k_, mixing_, s_ = fastica(
                       **kwargs
                    )

sklearn/decomposition/tests/test_fastica.py:113: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
sklearn/decomposition/_fastica.py:324: in fastica
    S = est._fit_transform(X, compute_sources=compute_sources)
sklearn/decomposition/_fastica.py:685: in _fit_transform
    self.mixing_ = linalg.pinv(self.components_, check_finite=False)
/usr/lib/python3/dist-packages/scipy/linalg/basic.py:1308: in pinv
    x, resids, rank, s = lstsq(a, b, cond=cond, check_finite=False)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

a = array([[-0.0178364 ,  0.02612615],
       [        nan,         nan]], dtype=float32)
b = array([[1., 0.],
       [0., 1.]], dtype=float32), cond = 2.384185791015625e-07, overwrite_a = False
overwrite_b = False, check_finite = False, lapack_driver = None

    def lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False,
              check_finite=True, lapack_driver=None):
        """
        Compute least-squares solution to equation Ax = b.
    
        Compute a vector x such that the 2-norm ``|b - A x|`` is minimized.
    
        Parameters
        ----------
        a : (M, N) array_like
            Left hand side array
        b : (M,) or (M, K) array_like
            Right hand side array
        cond : float, optional
            Cutoff for 'small' singular values; used to determine effective
            rank of a. Singular values smaller than
            ``rcond * largest_singular_value`` are considered zero.
        overwrite_a : bool, optional
            Discard data in `a` (may enhance performance). Default is False.
        overwrite_b : bool, optional
            Discard data in `b` (may enhance performance). Default is False.
        check_finite : bool, optional
            Whether to check that the input matrices contain only finite numbers.
            Disabling may give a performance gain, but may result in problems
            (crashes, non-termination) if the inputs do contain infinities or NaNs.
        lapack_driver : str, optional
            Which LAPACK driver is used to solve the least-squares problem.
            Options are ``'gelsd'``, ``'gelsy'``, ``'gelss'``. Default
            (``'gelsd'``) is a good choice.  However, ``'gelsy'`` can be slightly
            faster on many problems.  ``'gelss'`` was used historically.  It is
            generally slow but uses less memory.
    
            .. versionadded:: 0.17.0
    
        Returns
        -------
        x : (N,) or (N, K) ndarray
            Least-squares solution.  Return shape matches shape of `b`.
        residues : (K,) ndarray or float
            Square of the 2-norm for each column in ``b - a x``, if ``M > N`` and
            ``ndim(A) == n`` (returns a scalar if b is 1-D). Otherwise a
            (0,)-shaped array is returned.
        rank : int
            Effective rank of `a`.
        s : (min(M, N),) ndarray or None
            Singular values of `a`. The condition number of a is
            ``abs(s[0] / s[-1])``.
    
        Raises
        ------
        LinAlgError
            If computation does not converge.
    
        ValueError
            When parameters are not compatible.
    
        See Also
        --------
        scipy.optimize.nnls : linear least squares with non-negativity constraint
    
        Notes
        -----
        When ``'gelsy'`` is used as a driver, `residues` is set to a (0,)-shaped
        array and `s` is always ``None``.
    
        Examples
        --------
        >>> from scipy.linalg import lstsq
        >>> import matplotlib.pyplot as plt
    
        Suppose we have the following data:
    
        >>> x = np.array([1, 2.5, 3.5, 4, 5, 7, 8.5])
        >>> y = np.array([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6])
    
        We want to fit a quadratic polynomial of the form ``y = a + b*x**2``
        to this data.  We first form the "design matrix" M, with a constant
        column of 1s and a column containing ``x**2``:
    
        >>> M = x[:, np.newaxis]**[0, 2]
        >>> M
        array([[  1.  ,   1.  ],
               [  1.  ,   6.25],
               [  1.  ,  12.25],
               [  1.  ,  16.  ],
               [  1.  ,  25.  ],
               [  1.  ,  49.  ],
               [  1.  ,  72.25]])
    
        We want to find the least-squares solution to ``M.dot(p) = y``,
        where ``p`` is a vector with length 2 that holds the parameters
        ``a`` and ``b``.
    
        >>> p, res, rnk, s = lstsq(M, y)
        >>> p
        array([ 0.20925829,  0.12013861])
    
        Plot the data and the fitted curve.
    
        >>> plt.plot(x, y, 'o', label='data')
        >>> xx = np.linspace(0, 9, 101)
        >>> yy = p[0] + p[1]*xx**2
        >>> plt.plot(xx, yy, label='least squares fit, $y = a + bx^2$')
        >>> plt.xlabel('x')
        >>> plt.ylabel('y')
        >>> plt.legend(framealpha=1, shadow=True)
        >>> plt.grid(alpha=0.25)
        >>> plt.show()
    
        """
        a1 = _asarray_validated(a, check_finite=check_finite)
        b1 = _asarray_validated(b, check_finite=check_finite)
        if len(a1.shape) != 2:
            raise ValueError('Input array a should be 2-D')
        m, n = a1.shape
        if len(b1.shape) == 2:
            nrhs = b1.shape[1]
        else:
            nrhs = 1
        if m != b1.shape[0]:
            raise ValueError('Shape mismatch: a and b should have the same number'
                             ' of rows ({} != {}).'.format(m, b1.shape[0]))
        if m == 0 or n == 0:  # Zero-sized problem, confuses LAPACK
            x = np.zeros((n,) + b1.shape[1:], dtype=np.common_type(a1, b1))
            if n == 0:
                residues = np.linalg.norm(b1, axis=0)**2
            else:
                residues = np.empty((0,))
            return x, residues, 0, np.empty((0,))
    
        driver = lapack_driver
        if driver is None:
            driver = lstsq.default_lapack_driver
        if driver not in ('gelsd', 'gelsy', 'gelss'):
            raise ValueError('LAPACK driver "%s" is not found' % driver)
    
        lapack_func, lapack_lwork = get_lapack_funcs((driver,
                                                     '%s_lwork' % driver),
                                                     (a1, b1))
        real_data = True if (lapack_func.dtype.kind == 'f') else False
    
        if m < n:
            # need to extend b matrix as it will be filled with
            # a larger solution matrix
            if len(b1.shape) == 2:
                b2 = np.zeros((n, nrhs), dtype=lapack_func.dtype)
                b2[:m, :] = b1
            else:
                b2 = np.zeros(n, dtype=lapack_func.dtype)
                b2[:m] = b1
            b1 = b2
    
        overwrite_a = overwrite_a or _datacopied(a1, a)
        overwrite_b = overwrite_b or _datacopied(b1, b)
    
        if cond is None:
            cond = np.finfo(lapack_func.dtype).eps
    
        if driver in ('gelss', 'gelsd'):
            if driver == 'gelss':
                lwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond)
                v, x, s, rank, work, info = lapack_func(a1, b1, cond, lwork,
                                                        overwrite_a=overwrite_a,
                                                        overwrite_b=overwrite_b)
    
            elif driver == 'gelsd':
                if real_data:
                    lwork, iwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond)
>                   x, s, rank, info = lapack_func(a1, b1, lwork,
                                                   iwork, cond, False, False)
E                   ValueError: On entry to SLASCL parameter number 4 had an illegal value

/usr/lib/python3/dist-packages/scipy/linalg/basic.py:1210: ValueError
```

</details>

Trying to debug a bit further the code in fastica https://github.com/scikit-learn/scikit-learn/blob/b571d64066dd13b6f9b5eae399c6a20db39e5a83/sklearn/decomposition/_fastica.py#L85-L89

**atlas**

at one point w1 is exactly 0 and you divide by 0 so you get NaNs from there on:
```
component 0 iteration 0 w1 after decorrelation: [ 0.23501348 -0.19262002]
component 0 iteration 1 w1 before decorrelation: [0.04253763 0.40387636]
component 0 iteration 1 w1 after decorrelation: [0.04253763 0.40387636]
component 0 iteration 2 w1 before decorrelation: [ 0.13986371 -0.6734375 ]
component 0 iteration 2 w1 after decorrelation: [ 0.13986371 -0.6734375 ]
component 0 iteration 3 w1 before decorrelation: [-0.14534193  0.7702849 ]
component 0 iteration 3 w1 after decorrelation: [-0.14534193  0.7702849 ]
component 0 iteration 4 w1 before decorrelation: [ 0.14539295 -0.77054155]
component 0 iteration 4 w1 after decorrelation: [ 0.14539295 -0.77054155]
component 1 iteration 0 w1 before decorrelation: [-0.14537632  0.7704534 ]
component 1 iteration 0 w1 after decorrelation: [0. 0.]
component 1 iteration 1 w1 before decorrelation: [nan nan]
.
.
.
keeps going with nans until max_iter
```

**openblas**

instead of exactly zero you get a small value of the order of `1e-8` so smaller than the float32 precision)

```
component 0 iteration 0 w1 after decorrelation: [ 0.23501337 -0.19262013]
component 0 iteration 1 w1 before decorrelation: [0.04253733 0.4038767 ]
component 0 iteration 1 w1 after decorrelation: [0.04253733 0.4038767 ]
component 0 iteration 2 w1 before decorrelation: [ 0.13986367 -0.6734381 ]
component 0 iteration 2 w1 after decorrelation: [ 0.13986367 -0.6734381 ]
component 0 iteration 3 w1 before decorrelation: [-0.1453419   0.77028495]
component 0 iteration 3 w1 after decorrelation: [-0.1453419   0.77028495]
component 0 iteration 4 w1 before decorrelation: [ 0.14539292 -0.77054155]
component 0 iteration 4 w1 after decorrelation: [ 0.14539292 -0.77054155]
component 1 iteration 0 w1 before decorrelation: [-0.1453763  0.7704534]
component 1 iteration 0 w1 after decorrelation: [-1.4901161e-08  0.0000000e+00]
component 1 iteration 1 w1 before decorrelation: [-0.3199173  -0.06967697]
.
.
.
```

Suggestions more than welcome about how to fix or work around this issue!

Labels: Build / CI