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

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

Start with pca_variance_ratio in scanpy/plotting/tools/init.py and the ranking logic in scanpy/plotting/anndata.py, then reproduce the low-component case using svd_solver='arpack'. Check the variance-ratio values and their ordering before plotting; done means the warning and non-finite axis limits are resolved and the plotted PCs appear in the expected order.