scanpy.pp.pca fails on small datasets

[scottgigante]

scanpy.pp.pca fails if n_samples < 50 < n_features

import numpy as np
import scanpy as sc
import anndata

adata = anndata.AnnData(np.random.normal(0, 1, (40, 100)))
sc.pp.pca(adata)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/scottgigante/.local/lib/python3.8/site-packages/scanpy/preprocessing/_simple.py", line 531, in pca
    X_pca = pca_.fit_transform(X)
  File "/usr/lib/python3.8/site-packages/sklearn/decomposition/_pca.py", line 369, in fit_transform
    U, S, V = self._fit(X)
  File "/usr/lib/python3.8/site-packages/sklearn/decomposition/_pca.py", line 418, in _fit
    return self._fit_truncated(X, n_components, self._fit_svd_solver)
  File "/usr/lib/python3.8/site-packages/sklearn/decomposition/_pca.py", line 497, in _fit_truncated
    raise ValueError("n_components=%r must be between 1 and "
ValueError: n_components=50 must be between 1 and min(n_samples, n_features)=40 with svd_solver='arpack'

Versions:

scanpy==1.2.3.dev1409+g7ca201d.d20200112 anndata==0.6.22.post1 umap==0.3.10 numpy==1.18.0 scipy==1.4.1 pandas==0.25.3 scikit-learn==0.22 statsmodels==0.11.0rc1 python-igraph==0.7.1 louvain==0.6.1

[LuckyMD]

That's not a bug, that's a feature ;). You can only compute as many PCs as the minimum number of dimensions in n_samples and n_features. Do you feel as though the error message is unclear on this? I feel as though changing the default to match the setting can be dangerous as may not recall how many PCs were used then.

[scottgigante]

Three things:

1. If that is a feature, then this is a bug (runs without error):

   import numpy as np
   import scanpy as sc
   import anndata
   
   adata = anndata.AnnData(np.random.normal(0, 1, (40, 10)))
   sc.pp.pca(adata)
   

2. Defaults should work without tuning.

3. I feel as though changing the default to match the setting can be dangerous as may not recall how many PCs were used then.

   Given that I'm running sc.pp.pca without setting n_comps, I contend that the average user does not remember what the default value is. It would make more sense in both cases (n_features <= n_comps and n_samples <= n_comps) to throw a warning and set n_comps to the maximum allowable value.

[LuckyMD]

Hmm... you are right, that should also create an error by the above logic.

It does look like this check is done for the case that adata.n_vars < n_comps here:
https://github.com/theislab/scanpy/blob/be1a0555252cfd97b9d00f51dc5fbab462588da0/scanpy/preprocessing/_simple.py#L472-L477

I'm not sure why that wasn't also done for n_obs. @Koncopd you made this fix at the time... any reason for not also checking adata.n_obs in the same way? Could quickly add a check for adata.n_obs unless there is a reason not to @ivirshup, @flying-sheep?

[ivirshup]

I'm thinking n_comps should default to None, and we'd define behaviour like:

if n_comps is None:
    min_dim = min(adata.n_vars, adata.n_obs)
    if 50 >= min_dim:
        n_comps = min_dim - 1
    else:
        n_comps = 50

and we let sklearn throw an error if the user specified a number of components that doesn't work.

[gokceneraslan]

It's not only about adata.n_vars, don't forget about use_highly_variable option e. g. sc.pp.highly_variable_genes(adata, n_top_genes=49).

[ivirshup]

Thanks for the catch! It should be min_dim = min(*X.shape) where X is the selected representation.

[ivirshup]

Ha, actually we implement filtering for highly_variable_genes as taking a subset of the anndata object, so it is min(adata.n_vars, adata.n_obs)

[scottgigante]

Thanks all!