Filter.filter: generalize analysis and synthesis

Standardize signal shapes as N_NODES x N_SIGNALS x N_FEATURES

Author: mdeff

pygsp/filters/__init__.py

@@ -16,6 +16,7 @@
 .. autosummary::
 
     Filter.evaluate
+    Filter.filter
     Filter.analysis
     Filter.synthesis
     Filter.compute_frame

pygsp/filters/filter.py

@@ -203,6 +203,177 @@ def evaluate(self, x):
             y[i] = g(x)
         return y
 
+    def filter(self, s, method='chebyshev', order=30):
+        r"""
+        Filter signals with the filter bank (analysis or synthesis).
+
+        A signal is defined as a rank-3 tensor of shape ``(N_NODES, N_SIGNALS,
+        N_FEATURES)``, where ``N_NODES`` is the number of nodes in the graph,
+        ``N_SIGNALS`` is the number of independent signals, and ``N_FEATURES``
+        is the number of features which compose a graph signal, or the
+        dimensionality of a graph signal. For example if you filter a signal
+        with a filter bank of 8 filters, you're extracting 8 features and
+        decomposing your signal into 8 parts. That is called analysis. Your are
+        thus transforming your signal tensor from ``(G.N, 1, 1)`` to ``(G.N, 1,
+        8)``. Now you may want to combine back the features to form an unique
+        signal. For this you apply again 8 filters, one filter per feature, and
+        sum the result up. As such you're transforming your ``(G.N, 1, 8)``
+        tensor signal back to ``(G.N, 1, 1)``. That is known as synthesis. More
+        generally, you may want to map a set of features to another, though
+        that is not implemented yet.
+
+        The method computes the transform coefficients of a signal :math:`s`,
+        where the atoms of the transform dictionary are generalized
+        translations of each graph spectral filter to each vertex on the graph:
+
+        .. math:: c = D^* s,
+
+        where the columns of :math:`D` are :math:`g_{i,m} = T_i g_m` and
+        :math:`T_i` is a generalized translation operator applied to each
+        filter :math:`\hat{g}_m(\cdot)`. Each column of :math:`c` is the
+        response of the signal to one filter.
+
+        In other words, this function is applying the analysis operator
+        :math:`D^*`, respectively the synthesis operator :math:`D`, associated
+        with the frame defined by the filter bank to the signals.
+
+        Parameters
+        ----------
+        s : ndarray
+            Graph signals, a tensor of shape ``(N_NODES, N_SIGNALS,
+            N_FEATURES)``, where ``N_NODES`` is the number of nodes in the
+            graph, ``N_SIGNALS`` the number of independent signals you want to
+            filter, and ``N_FEATURES`` is either 1 (analysis) or the number of
+            filters in the filter bank (synthesis).
+        method : {'exact', 'chebyshev'}
+            Whether to use the exact method (via the graph Fourier transform)
+            or the Chebyshev polynomial approximation. A Lanczos
+            approximation is coming.
+        order : int
+            Degree of the Chebyshev polynomials.
+
+        Returns
+        -------
+        s : ndarray
+            Graph signals, a tensor of shape ``(N_NODES, N_SIGNALS,
+            N_FEATURES)``, where ``N_NODES`` and ``N_SIGNALS`` are the number
+            of nodes and signals of the signal tensor that pas passed in, and
+            ``N_FEATURES`` is either 1 (synthesis) or the number of filters in
+            the filter bank (analysis).
+
+        References
+        ----------
+        See :cite:`hammond2011wavelets` for details on filtering graph signals.
+
+        Examples
+        --------
+
+        Create a bunch of smooth signals by low-pass filtering white noise:
+
+        >>> import matplotlib.pyplot as plt
+        >>> G = graphs.Ring(N=60)
+        >>> G.estimate_lmax()
+        >>> s = np.random.RandomState(42).uniform(size=(G.N, 10))
+        >>> taus = [1, 10, 100]
+        >>> s = filters.Heat(G, taus).filter(s)
+        >>> s.shape
+        (60, 10, 3)
+
+        Plot the 3 smoothed versions of the 10th signal:
+
+        >>> fig, ax = plt.subplots()
+        >>> G.set_coordinates('line1D')  # To visualize multiple signals in 1D.
+        >>> G.plot_signal(s[:, 9, :], ax=ax)
+        >>> legend = [r'$\tau={}$'.format(t) for t in taus]
+        >>> ax.legend(legend)  # doctest: +ELLIPSIS
+        <matplotlib.legend.Legend object at ...>
+
+        Low-pass filter a delta to create a localized smooth signal:
+
+        >>> G = graphs.Sensor(30, seed=42)
+        >>> G.compute_fourier_basis()  # Reproducible computation of lmax.
+        >>> s1 = np.zeros(G.N)
+        >>> s1[13] = 1
+        >>> s1 = filters.Heat(G, 3).filter(s1)
+        >>> s1.shape
+        (30, 1, 1)
+
+        Filter and reconstruct our signal:
+
+        >>> g = filters.MexicanHat(G, Nf=4)
+        >>> s2 = g.filter(s1)
+        >>> s2.shape
+        (30, 1, 4)
+        >>> s2 = g.filter(s2)
+        >>> s2.shape
+        (30, 1, 1)
+
+        Look how well we were able to reconstruct:
+
+        >>> fig, axes = plt.subplots(1, 2)
+        >>> G.plot_signal(s1, ax=axes[0])
+        >>> G.plot_signal(s2, ax=axes[1])
+        >>> print('{:.5f}'.format(np.linalg.norm(s1 - s2)))
+        0.29620
+
+        Perfect reconstruction with Itersine, a tight frame:
+
+        >>> g = filters.Itersine(G)
+        >>> s2 = g.filter(s1, method='exact')
+        >>> s2 = g.filter(s2, method='exact')
+        >>> np.linalg.norm(s1 - s2) < 1e-10
+        True
+
+        """
+        s = self.G.sanitize_signal(s)
+        N_NODES, N_SIGNALS, N_FEATURES_IN = s.shape
+
+        # TODO: generalize to 2D (m --> n) filter banks.
+        # Only 1 --> Nf (analysis) and Nf --> 1 (synthesis) for now.
+        if N_FEATURES_IN not in [1, self.Nf]:
+            raise ValueError('Last dimension (N_FEATURES) should either be '
+                             '1 or the number of filters (Nf), '
+                             'not {}.'.format(s.shape))
+        N_FEATURES_OUT = self.Nf if N_FEATURES_IN == 1 else 1
+
+        if method == 'exact':
+
+            axis = 1 if N_FEATURES_IN == 1 else 2
+            f = self.evaluate(self.G.e)
+            f = np.expand_dims(f.T, axis)
+            assert f.shape == (N_NODES, N_FEATURES_IN, N_FEATURES_OUT)
+
+            s = self.G.gft2(s)
+            s = np.matmul(s, f)
+            s = self.G.igft2(s)
+
+        elif method == 'chebyshev':
+
+            # TODO: update Chebyshev implementation (after 2D filter banks).
+            c = approximations.compute_cheby_coeff(self, m=order)
+
+            if N_FEATURES_IN == 1:  # Analysis.
+                s = s.squeeze(axis=2)
+                s = approximations.cheby_op(self.G, c, s)
+                s = s.reshape((N_NODES, N_FEATURES_OUT, N_SIGNALS), order='F')
+                s = s.swapaxes(1, 2)
+
+            elif N_FEATURES_IN == self.Nf:  # Synthesis.
+                s = s.swapaxes(1, 2)
+                s_in = s.reshape((N_NODES*N_FEATURES_IN, N_SIGNALS), order='F')
+                s = np.zeros((N_NODES, N_SIGNALS))
+                tmpN = np.arange(N_NODES, dtype=int)
+                for i in range(N_FEATURES_IN):
+                    s += approximations.cheby_op(self.G,
+                                                 c[i],
+                                                 s_in[i * N_NODES + tmpN])
+                s = np.expand_dims(s, 2)
+
+        else:
+            raise ValueError('Unknown method {}.'.format(method))
+
+        return s
+
     def inverse(self, c):
         r"""
         Not implemented yet.

pygsp/graphs/__init__.py

@@ -86,6 +86,7 @@
     Graph.set_coordinates
     Graph.subgraph
     Graph.extract_components
+    Graph.sanitize_signal
 
 Graph models
 ============

pygsp/graphs/fourier.py

@@ -139,6 +139,11 @@ def gft(self, s):
         """
         return np.dot(np.conjugate(self.U.T), s)  # True Hermitian here.
 
+    def gft2(self, s):
+        s = self.sanitize_signal(s)
+        U = np.conjugate(self.U)  # True Hermitian. (Although U is often real.)
+        return np.tensordot(U, s, ([0], [0]))
+
     def igft(self, s_hat):
         r"""Compute the inverse graph Fourier transform.
 
@@ -171,6 +176,10 @@ def igft(self, s_hat):
         """
         return np.dot(self.U, s_hat)
 
+    def igft2(self, s_hat):
+        s_hat = self.sanitize_signal(s_hat)
+        return np.tensordot(self.U, s_hat, ([1], [0]))
+
     def translate(self, f, i):
         r"""Translate the signal *f* to the node *i*.
 

pygsp/graphs/graph.py

@@ -725,6 +725,57 @@ def get_edge_list(self):
 
             return v_in, v_out, weights
 
+    def sanitize_signal(self, s):
+        r"""Standardize signal shape.
+
+        Add singleton dimensions at the end and check the resulting shape.
+
+        Parameters
+        ----------
+        s : ndarray
+            Signal tensor of shape ``(N_NODES)``, ``(N_NODES, N_SIGNALS)``, or
+            ``(N_NODES, N_SIGNALS, N_FEATURES)``.
+
+        Returns
+        -------
+        s : ndarray
+            Signal tensor of shape ``(N_NODES, N_SIGNALS, N_FEATURES)``.
+
+        Raises
+        ------
+        ValueError
+            If the passed signal tensor is more than 3 dimensions or if the
+            first dimension's size is not the number of nodes.
+
+        Examples
+        --------
+        >>> G = graphs.Logo()
+        >>> s = np.ones(G.N)  # One signal, one feature.
+        >>> G.sanitize_signal(s).shape
+        (1130, 1, 1)
+        >>> s = np.ones((G.N, 10))  # Ten signals of one feature.
+        >>> G.sanitize_signal(s).shape
+        (1130, 10, 1)
+        >>> s = np.ones((G.N, 10, 5))  # Ten signals of 5 features.
+        >>> G.sanitize_signal(s).shape
+        (1130, 10, 5)
+
+        """
+        if s.ndim == 1:
+            # Single signal, single feature.
+            s = np.expand_dims(s, axis=1)
+
+        if s.ndim == 2:
+            # Multiple signals, single feature.
+            s = np.expand_dims(s, axis=2)
+
+        if s.ndim != 3 or s.shape[0] != self.N:
+            raise ValueError('Signal must have shape N_NODES x N_SIGNALS x '
+                             'N_FEATURES, not {}. Last singleton dimensions '
+                             'may be omitted.'.format(s.shape))
+
+        return s
+
     def modulate(self, f, k):
         r"""Modulate the signal *f* to the frequency *k*.
 

pygsp/tests/test_filters.py

@@ -40,9 +40,30 @@ def _test_synthesis(self, f):
             self.assertRaises(NotImplementedError, f.synthesis, S,
                               method='lanczos')
 
+    def _test_filter(self, f, tight):
+        # Analysis.
+        s2 = f.filter(self._signal, method='exact')
+        s3 = f.filter(self._signal, method='chebyshev', order=100)
+
+        # Synthesis.
+        s4 = f.filter(s2, method='exact')
+        s5 = f.filter(s3, method='chebyshev', order=100)
+
+        if f.Nf < 100:
+            # TODO: does not pass for Gabor.
+            np.testing.assert_allclose(s2, s3, rtol=0.1, atol=0.01)
+            np.testing.assert_allclose(s4, s5, rtol=0.1, atol=0.01)
+
+        if tight:
+            A, _ = f.estimate_frame_bounds(use_eigenvalues=True)
+            np.testing.assert_allclose(s4.squeeze(), A * self._signal)
+            assert np.linalg.norm(s5.squeeze() - A * self._signal) < 0.1
+
     def _test_methods(self, f, tight):
         self.assertIs(f.G, self._G)
 
+        self._test_filter(f, tight)
+
         c_exact = f.analysis(self._signal, method='exact')
         c_cheby = f.analysis(self._signal, method='chebyshev')
         self.assertEqual(c_exact.shape, c_cheby.shape)

pygsp/tests/test_graphs.py

@@ -119,6 +119,23 @@ def test_set_coordinates(self):
         G.set_coordinates('community2D')
         self.assertRaises(ValueError, G.set_coordinates, 'invalid')
 
+    def test_sanitize_signal(self):
+        s1 = np.arange(self._G.N)
+        s2 = np.reshape(s1, (self._G.N, 1))
+        s3 = np.reshape(s1, (self._G.N, 1, 1))
+        s4 = np.arange(self._G.N*10).reshape((self._G.N, 10))
+        s5 = np.reshape(s4, (self._G.N, 10, 1))
+        s1 = self._G.sanitize_signal(s1)
+        s2 = self._G.sanitize_signal(s2)
+        s3 = self._G.sanitize_signal(s3)
+        s4 = self._G.sanitize_signal(s4)
+        s5 = self._G.sanitize_signal(s5)
+        np.testing.assert_equal(s2, s1)
+        np.testing.assert_equal(s3, s1)
+        np.testing.assert_equal(s5, s4)
+        self.assertRaises(ValueError, self._G.sanitize_signal,
+                          np.ones((2, 2, 2, 2)))
+
     def test_nngraph(self):
         Xin = np.arange(90).reshape(30, 3)
         dist_types = ['euclidean', 'manhattan', 'max_dist', 'minkowski']