PCA, "RuntimeWarning: divide by zero encountered in log" and wrong PCs ordering with svd_solver arpack

[fbrundu]

Hi,
I am trying to evaluate a PCA on a dataset with a very low number of Principal Components.

I am getting the following error when trying to plot the PCA:

sc.pl.pca_variance_ratio(adata_h, log=True, save=True)

Result:

/Users/user/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/scanpy/plotting/anndata.py:530: RuntimeWarning: divide by zero encountered in log
  if log: scores = np.log(scores)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-7-48cc676a34cc> in <module>()
      1 # log is natural logarithm
----> 2 sc.pl.pca_variance_ratio(adata_h, log=True, save=True)

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/scanpy/plotting/tools/__init__.py in pca_variance_ratio(adata, log, show, save)
    157         default filename. Infer the filetype if ending on {{'.pdf', '.png', '.svg'}}.
    158     """
--> 159     ranking(adata, 'uns', 'variance_ratio', dictionary='pca', labels='PC', log=log)
    160     utils.savefig_or_show('pca_variance_ratio', show=show, save=save)
    161 

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/scanpy/plotting/anndata.py in ranking(adata, attr, keys, dictionary, indices, labels, color, n_points, log, show)
    555         score_min, score_max = np.min(score[indices]), np.max(score[indices])
    556         pl.ylim((0.95 if score_min > 0 else 1.05) * score_min,
--> 557                 (1.05 if score_max > 0 else 0.95) * score_max)
    558     if show == False: return gs
    559 

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/pyplot.py in ylim(*args, **kwargs)
   1588     if not args and not kwargs:
   1589         return ax.get_ylim()
-> 1590     ret = ax.set_ylim(*args, **kwargs)
   1591     return ret
   1592 

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axes/_base.py in set_ylim(self, bottom, top, emit, auto, **kw)
   3455             bottom, top = bottom
   3456 
-> 3457         bottom = self._validate_converted_limits(bottom, self.convert_yunits)
   3458         top = self._validate_converted_limits(top, self.convert_yunits)
   3459 

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axes/_base.py in _validate_converted_limits(self, limit, convert)
   3066                     (not np.isreal(converted_limit) or
   3067                         not np.isfinite(converted_limit))):
-> 3068                 raise ValueError("Axis limits cannot be NaN or Inf")
   3069             return converted_limit
   3070 

ValueError: Axis limits cannot be NaN or Inf

posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values
posx and posy should be finite values

/Users/user/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:83: RuntimeWarning: invalid value encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/pylabtools.py in <lambda>(fig)
    241         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    242     if 'retina' in formats or 'png2x' in formats:
--> 243         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))
    244     if 'jpg' in formats or 'jpeg' in formats:
    245         jpg_formatter.for_type(Figure, lambda fig: print_figure(fig, 'jpg', **kwargs))

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/pylabtools.py in retina_figure(fig, **kwargs)
    131 def retina_figure(fig, **kwargs):
    132     """format a figure as a pixel-doubled (retina) PNG"""
--> 133     pngdata = print_figure(fig, fmt='retina', **kwargs)
    134     # Make sure that retina_figure acts just like print_figure and returns
    135     # None when the figure is empty.

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    123 
    124     bytes_io = BytesIO()
--> 125     fig.canvas.print_figure(bytes_io, **kw)
    126     data = bytes_io.getvalue()
    127     if fmt == 'svg':

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2261                 orientation=orientation,
   2262                 bbox_inches_restore=_bbox_inches_restore,
-> 2263                 **kwargs)
   2264         finally:
   2265             if bbox_inches and restore_bbox:

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py in print_png(self, filename_or_obj, *args, **kwargs)
    515 
    516     def print_png(self, filename_or_obj, *args, **kwargs):
--> 517         FigureCanvasAgg.draw(self)
    518         renderer = self.get_renderer()
    519         original_dpi = renderer.dpi

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py in draw(self)
    435             # if toolbar:
    436             #     toolbar.set_cursor(cursors.WAIT)
--> 437             self.figure.draw(self.renderer)
    438             # A GUI class may be need to update a window using this draw, so
    439             # don't forget to call the superclass.

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/figure.py in draw(self, renderer)
   1491 
   1492             mimage._draw_list_compositing_images(
-> 1493                 renderer, self, artists, self.suppressComposite)
   1494 
   1495             renderer.close_group('figure')

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    139     if not_composite or not has_images:
    140         for a in artists:
--> 141             a.draw(renderer)
    142     else:
    143         # Composite any adjacent images together

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
   2633             renderer.stop_rasterizing()
   2634 
-> 2635         mimage._draw_list_compositing_images(renderer, self, artists)
   2636 
   2637         renderer.close_group('axes')

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    139     if not_composite or not has_images:
    140         for a in artists:
--> 141             a.draw(renderer)
    142     else:
    143         # Composite any adjacent images together

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axis.py in draw(self, renderer, *args, **kwargs)
   1188         renderer.open_group(__name__)
   1189 
-> 1190         ticks_to_draw = self._update_ticks(renderer)
   1191         ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw,
   1192                                                                 renderer)

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axis.py in _update_ticks(self, renderer)
   1026 
   1027         interval = self.get_view_interval()
-> 1028         tick_tups = list(self.iter_ticks())  # iter_ticks calls the locator
   1029         if self._smart_bounds and tick_tups:
   1030             # handle inverted limits

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axis.py in iter_ticks(self)
    969         Iterate through all of the major and minor ticks.
    970         """
--> 971         majorLocs = self.major.locator()
    972         majorTicks = self.get_major_ticks(len(majorLocs))
    973         self.major.formatter.set_locs(majorLocs)

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/ticker.py in __call__(self)
   1952     def __call__(self):
   1953         vmin, vmax = self.axis.get_view_interval()
-> 1954         return self.tick_values(vmin, vmax)
   1955 
   1956     def tick_values(self, vmin, vmax):

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/ticker.py in tick_values(self, vmin, vmax)
   1960         vmin, vmax = mtransforms.nonsingular(
   1961             vmin, vmax, expander=1e-13, tiny=1e-14)
-> 1962         locs = self._raw_ticks(vmin, vmax)
   1963 
   1964         prune = self._prune

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/ticker.py in _raw_ticks(self, vmin, vmax)
   1906         if self._nbins == 'auto':
   1907             if self.axis is not None:
-> 1908                 nbins = np.clip(self.axis.get_tick_space(),
   1909                                 max(1, self._min_n_ticks - 1), 9)
   1910             else:

~/.pyenv/versions/3.6.5/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axis.py in get_tick_space(self)
   2127         size = tick.label1.get_size() * 3
   2128         if size > 0:
-> 2129             return int(np.floor(length / size))
   2130         else:
   2131             return 2**31 - 1

ValueError: cannot convert float NaN to integer

<Figure size 800x800 with 1 Axes>

I think the problem is that the variance explained for some components is 0, so taking the log (with log=True will ask scanpy to plot the variance explained starting from minus infinity. Would it be possible to plot only the non-zero variance components as a default?

Thanks,
Francesco

[fbrundu]

Additionally, this only happens when PCA has been computed with svd_solver='arpack' parameter. The current default (svd_solver='auto') produces an always non-zero variance explained for the components. Additionally the order of variance explained is wrong with the svd_solver='arpack' parameter. I attach some examples:

>>> sc.tl.pca(adata_h, svd_solver='arpack')

>>> adata_h.uns['pca']

{'variance': array([   0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      ,    0.      ,    0.      ,    0.      ,
           0.      ,    0.      , 2633.797   ,  457.86526 ,  316.44687 ,
         237.71556 ,  143.87927 ,  119.6577  ,  105.01371 ,   91.51559 ,
          66.951355,   61.23979 ,   59.957714,   58.998177,   57.82413 ],
       dtype=float32),
 'variance_ratio': array([0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        , 0.5971161 , 0.10380401, 0.07174262,
        0.05389321, 0.0326193 , 0.02712796, 0.02380797, 0.02074778,
        0.01517874, 0.01388386, 0.01359319, 0.01337565, 0.01310948],
       dtype=float32)}

>>> sc.tl.pca(adata_h)
Note that scikit-learn's randomized PCA might not be exactly reproducible across different computational platforms. For exact reproducibility, choose `svd_solver='arpack'.` This will likely become the Scanpy default in the future.

>>> adata_h.uns['pca']

{'variance': array([2.63379761e+03, 4.57865112e+02, 3.16446930e+02, 2.37715851e+02,
        1.43879318e+02, 1.19657700e+02, 1.05013855e+02, 9.15156784e+01,
        6.69513855e+01, 6.12398453e+01, 5.99577942e+01, 5.89982376e+01,
        5.78241539e+01, 6.29622976e-09, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12, 9.35712879e-12, 9.35712879e-12,
        9.35712879e-12, 9.35712879e-12], dtype=float32),
 'variance_ratio': array([5.9711617e-01, 1.0380398e-01, 7.1742639e-02, 5.3893276e-02,
        3.2619312e-02, 2.7127957e-02, 2.3808008e-02, 2.0747796e-02,
        1.5178750e-02, 1.3883870e-02, 1.3593213e-02, 1.3375668e-02,
        1.3109488e-02, 1.4274373e-12, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15, 2.1213828e-15, 2.1213828e-15,
        2.1213828e-15, 2.1213828e-15], dtype=float32)}

[falexwolf]

Arrrgh, this prepending of the 0s is a bug in scikit-learn. Thank you very much for pointing it out. Let me briefly think about solving this elegantly...

[ivirshup]

Still occurring, though it's a really weird case. Here's some code to reproduce:

a = sc.AnnData(np.ones((100, 100)))
sc.pp.pca(a)
# RuntimeWarning: invalid value encountered in true_divide
#  self.explained_variance_ / total_var.sum()
sc.pl.pca(a)  # Throws the same error as above
a.uns["pca"]

{'params': {'zero_center': True, 'use_highly_variable': False},
 'variance': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       dtype=float32),
 'variance_ratio': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
        nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       dtype=float32)}

Though I do get the same behaviour with svd_solver="randomized"