# #############################################################################
# conv.py
# =======
# Author : Matthieu Simeoni [matthieu.simeoni@gmail.com]
# #############################################################################
r"""
Graph convolution operator.
"""
from typing import Union
import numpy as np
import pygsp
from pycsou.linop.base import SparseLinearOperator, PolynomialLinearOperator
[docs]class GraphConvolution(PolynomialLinearOperator):
r"""
Graph convolution.
Convolve a signal :math:`\mathbf{u}\in\mathbb{C}^N` defined on a graph with a polynomial filter :math:`\mathbf{D}:\mathbb{C}^N\rightarrow \mathbb{C}^N`
of the form:
.. math::
\mathbf{D}=\sum_{k=0}^K \theta_k \mathbf{L}^k,
where :math:`\mathbf{L}:\mathbb{C}^N\rightarrow \mathbb{C}^N` is the *normalised graph Laplacian* (see [FuncSphere]_ Section 2.3 of Chapter 6).
Examples
--------
.. testsetup::
import numpy as np
from pygsp.graphs import RandomRegular
from pycgsp.linop.conv import GraphConvolution
np.random.seed(0)
.. doctest::
>>> G = RandomRegular(seed=0)
>>> G.compute_laplacian(lap_type='normalized')
>>> signal = np.random.binomial(n=1,p=0.2,size=G.N)
>>> coefficients = np.ones(shape=(3,))
>>> ConvOp = GraphConvolution(Graph=G, coefficients=coefficients)
>>> filtered = ConvOp * signal
.. plot::
import numpy as np
from pygsp.graphs import Ring
from pycgsp.linop.conv import GraphConvolution
np.random.seed(0)
G = Ring(N=32, k=2)
G.compute_laplacian(lap_type='normalized')
G.set_coordinates(kind='spring')
signal = np.random.binomial(n=1,p=0.2,size=G.N)
coefficients = np.ones(3)
ConvOp = GraphConvolution(Graph=G, coefficients=coefficients)
e1 = np.zeros(shape=G.N)
e1[0] = 1
filter = ConvOp * e1
filtered = ConvOp * signal
plt.figure()
ax=plt.gca()
G.plot_signal(signal, ax=ax, backend='matplotlib')
plt.title('Signal')
plt.axis('equal')
plt.figure()
ax=plt.gca()
G.plot_signal(filter, ax=ax, backend='matplotlib')
plt.title('Filter')
plt.axis('equal')
plt.figure()
ax=plt.gca()
G.plot_signal(filtered, ax=ax, backend='matplotlib')
plt.title('Filtered Signal')
plt.axis('equal')
Notes
-----
The ``GraphConvolution`` operator is self-adjoint and operates in a matrix-free fashion, as described in Section 4.3, Chapter 7 of [FuncSphere]_.
See Also
--------
:py:class:`~pycgsp.linop.diff.GraphLaplacian`
"""
[docs] def __init__(self, Graph: pygsp.graphs.Graph, coefficients: Union[np.ndarray, list, tuple]):
r"""
Parameters
----------
Graph: `pygsp.graphs.Graph <https://pygsp.readthedocs.io/en/stable/reference/graphs.html#pygsp.graphs.Graph>`_
Graph on which the signal is defined, with normalised Laplacian ``Graph.L`` precomputed (see `pygsp.graphs.Graph.compute_laplacian(lap_type='normalized') <https://pygsp.readthedocs.io/en/stable/reference/graphs.html#pygsp.graphs.Graph.compute_laplacian>`_.
coefficients: Union[np.ndarray, list, tuple]
Coefficients :math:`\{\theta_k, \,k=0,\ldots,K\}\subset \mathbb{C}` of the polynomial filter.
dtype: type
Type of the entries of the graph filer.
Raises
------
AttributeError
If ``Graph.L`` does not exist.
NotImplementedError
If ``Graph.lap_type`` is 'combinatorial'.
"""
self.Graph = Graph
if Graph.L is None:
raise AttributeError(
r'Please compute the normalised Laplacian of the graph with the routine https://pygsp.readthedocs.io/en/stable/reference/graphs.html#pygsp.graphs.Graph.compute_laplacian')
elif Graph.lap_type != 'normalized':
raise NotImplementedError(r'Combinatorial graph Laplacians are not supported.')
else:
L = self.Graph.L.tocsc()
Lop = SparseLinearOperator(L, is_symmetric=True)
self.coefficients = coefficients
super(GraphConvolution, self).__init__(LinOp=Lop, coeffs=self.coefficients)