def compute_fourier_basis(self, n_eigenvectors=None):
        r"""Compute the (partial) Fourier basis of the graph (cached).

        The result is cached and accessible by the :attr:`U`, :attr:`e`,
        :attr:`lmax`, and :attr:`coherence` properties.

        Parameters
        ----------
        n_eigenvectors : int or `None`
            Number of eigenvectors to compute. If `None`, all eigenvectors
            are computed. (default: None)

        Notes
        -----
        'G.compute_fourier_basis()' computes a full eigendecomposition of
        the graph Laplacian :math:`L` such that:

        .. math:: L = U \Lambda U^*,

        or a partial eigendecomposition of the graph Laplacian :math:`L`
        such that:

        .. math:: L \approx U \Lambda U^*,

        where :math:`\Lambda` is a diagonal matrix of eigenvalues and the
        columns of :math:`U` are the eigenvectors.

        *G.e* is a vector of length `n_eigenvectors` :math:`\le` *G.N*
        containing the Laplacian eigenvalues. The largest eigenvalue is stored
        in *G.lmax*. The eigenvectors are stored as column vectors of *G.U* in
        the same order that the eigenvalues. Finally, the coherence of the
        Fourier basis is found in *G.coherence*.

        References
        ----------
        See :cite:`chung1997spectral`.

        Examples
        --------
        >>> G = graphs.Torus()
        >>> G.compute_fourier_basis(n_eigenvectors=64)
        >>> G.U.shape
        (256, 64)
        >>> G.e.shape
        (64,)
        >>> G.compute_fourier_basis()
        >>> G.U.shape
        (256, 256)
        >>> G.e.shape
        (256,)
        >>> G.lmax == G.e[-1]
        True
        >>> G.coherence < 1
        True

        """
        if n_eigenvectors is None:
            n_eigenvectors = self.n_vertices

        if (self._U is not None and n_eigenvectors <= len(self._e)):
            return

        assert self.L.shape == (self.n_vertices, self.n_vertices)
        if self.n_vertices**2 * n_eigenvectors > 3000**3:
            self.logger.warning(
                'Computing the {0} eigendecomposition of a large matrix ({1} x'
                ' {1}) is expensive. Consider decreasing n_eigenvectors '
                'or, if using the Fourier basis to filter, using a '
                'polynomial filter instead.'.format(
                    'full' if n_eigenvectors == self.N else 'partial',
                    self.N))

        # TODO: handle non-symmetric Laplacians. Test lap_type?
        if n_eigenvectors == self.n_vertices:
            self._e, self._U = np.linalg.eigh(self.L.toarray())
        else:
            # fast partial eigendecomposition of hermitian matrices
            self._e, self._U = sparse.linalg.eigsh(self.L,
                                                   n_eigenvectors,
                                                   which='SM')
        # Columns are eigenvectors. Sorted in ascending eigenvalue order.

        # Smallest eigenvalue should be zero: correct numerical errors.
        # Eigensolver might sometimes return small negative values, which
        # filter's implementations may not anticipate. Better for plotting too.
        assert -1e-5 < self._e[0] < 1e-5
        self._e[0] = 0

        # Bounded spectrum.
        assert self._e[-1] <= self._get_upper_bound() + 1e-5

        assert np.max(self._e) == self._e[-1]
        if n_eigenvectors == self.N:
            self._lmax = self._e[-1]
            self._lmax_method = 'fourier'
            self._coherence = np.max(np.abs(self._U))