It might make sense to decide how to solve depending on p / n, right?

Yes this is correct. But I think that the region boundary is hard to perfectly determine.

There are actually quite a lot of calls to linalg.svd in scikit-learn where you do not need the three outputs, but I don't know if it's the same story as for ICA where it is in 99% of the cases much cheaper to simply do eigh(X.dot(X.T)).

It's possible to at least roughly determine a good cut-off for the heuristic, while saying "people usually use this with n>>p" is a hard to support statement and we'd do a disservice to everybody using it in a different setting, even if that's a small group.

Ok, I have pushed a new commit where the heuristic is simply n < p. On my machine eigh is faster by a factor ~2 when n=p.

PS: Just figured out that what I called nin the above conversation is actually p in the code and vice versa.

After discussing with @agramfort , I have rather added an option to choose different solvers for svd. Indeed, linalg.svd is numerically more stable than the eighmethod, so it might be a problem when the input data are degenerate.

Here is a simple example to convince oneself that it accelerates fastica:

import numpy as np
from sklearn.decomposition import fastica


p, n = 100, 50000  # Typical EEG size problem
n_components = 10


S = np.random.laplace(size=(p, n))  # Generate sources
A = np.random.randn(p, p)  # Generate mixing matrix
X = np.dot(A, S).T  # Generate signals

%timeit fastica(X, n_components=n_components, svd_solver='svd')
1.02 s ± 55.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit fastica(X, n_components=n_components, svd_solver='eigh')
191 ms ± 68.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Continuing in #22527 if anyone is interested!