Source code for pycgsp.linop.diff

# #############################################################################
# diff.py
# =======
# Author : Matthieu Simeoni [matthieu.simeoni@gmail.com]
# #############################################################################

r"""
Graph differential operators.
"""

import numpy as np
import pygsp

from pycsou.core.linop import LinearOperator
from pycsou.linop.base import SparseLinearOperator, IdentityOperator, PolynomialLinearOperator

[docs]class GraphLaplacian(LinearOperator): r""" Graph Laplacian. Normalised graph Laplacian for signals defined on graphs. Examples -------- .. plot:: import numpy as np from pygsp.graphs import Ring from pycgsp.linop.diff import GraphLaplacian np.random.seed(1) G = Ring(N=32, k=4) G.compute_laplacian(lap_type='normalized') G.set_coordinates(kind='spring') x = np.arange(G.N) signal = np.piecewise(x, [x < G.N//3, (x >= G.N//3) * (x< 2 * G.N//3), x>=2 * G.N//3], [lambda x: -x, lambda x: 3 * x - 4 * G.N//3, lambda x: -0.5 * x + G.N]) Lap = GraphLaplacian(Graph=G) lap_sig = Lap * signal plt.figure() ax=plt.gca() G.plot_signal(signal, ax=ax, backend='matplotlib') plt.title('Signal') plt.axis('equal') plt.figure() plt.plot(signal) plt.title('Signal') plt.figure() ax=plt.gca() G.plot_signal(lap_sig, ax=ax, backend='matplotlib') plt.title('Laplacian of signal') plt.axis('equal') plt.figure() plt.plot(-lap_sig) plt.title('Laplacian of signal') Notes ----- For undirected graphs, the normalized graph Laplacian is defined as .. math:: \mathbf{L} = \mathbf{I} - \mathbf{D}^{-1/2} \mathbf{W} \mathbf{D}^{-1/2}, where :math:`\mathbf{I}` is the identity matrix, :math:`\mathbf{W}` is the weighted adjacency matrix and :math:`\mathbf{D}` the weighted degree matrix. For directed graphs, the Laplacians are built from a symmetrized version of the weighted adjacency matrix that is the average of the weighted adjacency matrix and its transpose. As the Laplacian is defined as the divergence of the gradient, it is not affected by the orientation of the edges. For both Laplacians, the diagonal entries corresponding to disconnected nodes (i.e., nodes with degree zero) are set to zero. The ``GraphLaplacian`` operator is self-adjoint. See Also -------- :py:class:`~pycgsp.linop.diff.GraphGradient`, :py:func:`~pycgsp.linop.diff.GeneralisedGraphLaplacian` :py:class:`~pycgsp.linop.conv.GraphConvolution` """
[docs] def __init__(self, Graph: pygsp.graphs.Graph, dtype: type = np.float): 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>`_. 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: self.L = self.Graph.L.tocsc() super(GraphLaplacian, self).__init__(shape=self.Graph.W.shape, dtype=dtype, is_explicit=False, is_dense=False, is_sparse=False, is_dask=False, is_symmetric=True)
def __call__(self, x: np.ndarray) -> np.ndarray: return self.L.dot(x)
[docs] def adjoint(self, y: np.ndarray) -> np.ndarray: return self(y)
[docs]class GraphGradient(LinearOperator): r""" Graph gradient. Gradient operator for signals defined on graphs. Examples -------- .. testsetup:: import numpy as np from pygsp.graphs import Ring from pycgsp.linop.diff import GraphLaplacian, GraphGradient np.random.seed(1) .. doctest:: >>> G = Ring(N=32, k=4) >>> G.compute_laplacian(lap_type='normalized') >>> G.compute_differential_operator() >>> G.set_coordinates(kind='spring') >>> x = np.arange(G.N) >>> signal = np.piecewise(x, [x < G.N//3, (x >= G.N//3) * (x< 2 * G.N//3), x>=2 * G.N//3], [lambda x: -x, lambda x: 3 * x - 4 * G.N//3, lambda x: -0.5 * x + G.N]) >>> Lap = GraphLaplacian(Graph=G) >>> Grad = GraphGradient(Graph=G) >>> lap_sig = Lap * signal >>> lap_sig2 = Grad.adjoint(Grad(signal)) >>> np.allclose(lap_sig, lap_sig2) True Notes ----- The adjoint of the ``GraphGradient`` operator is called the graph divergence operator. Warnings -------- In the newest version of PyGSP (> 0.5.1) the convention is changed: ``Graph.D`` is the divergence operator and ``Graph.D.transpose()`` the gradient (see routine `Graph.compute_differential_operator <https://pygsp.readthedocs.io/en/latest/reference/graphs.html#pygsp.graphs.Graph.compute_differential_operator>`_). The code should be adapted when this new version is released. See Also -------- :py:class:`~pycgsp.linop.diff.GraphLaplacian`, :py:func:`~pycgsp.linop.diff.GeneralisedGraphLaplacian` """
[docs] def __init__(self, Graph: pygsp.graphs.Graph, dtype: type = np.float): 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 differential operator ``Graph.D`` precomputed (see `pygsp.graphs.Graph.compute_differential_operator() <https://pygsp.readthedocs.io/en/stable/reference/graphs.html#pygsp.graphs.Graph.compute_differential_operator>`_. dtype: type Type of the entries of the graph filer. Raises ------ AttributeError If ``Graph.D`` does not exist. """ self.Graph = Graph if Graph.D is None: raise AttributeError( r'Please compute the differential operator of the graph with the routine https://pygsp.readthedocs.io/en/stable/reference/graphs.html#pygsp.graphs.Graph.compute_differential_operator') else: self.D = self.Graph.D.tocsc() super(GraphGradient, self).__init__(shape=self.Graph.W.shape, dtype=dtype, is_explicit=False, is_dense=False, is_sparse=False, is_dask=False, is_symmetric=True)
def __call__(self, x: np.ndarray) -> np.ndarray: return self.D.dot(x)
[docs] def adjoint(self, y: np.ndarray) -> np.ndarray: return self.D.conj().transpose().dot(y)
[docs]def GeneralisedGraphLaplacian(Graph: pygsp.graphs.Graph, kind: str = 'iterated', **kwargs): r""" Generalised graph Laplacian operator. Generalised Laplacian operator signals defined on graphs. 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>`_. dtype: type Type of the entries of the graph filer. kind: str Type of generalised differential operator (``'iterated'``, ``'sobolev'``, ``'polynomial'``). Depending on the cases, the ``GeneralisedLaplacian`` operator is defined as follows: * ``'iterated'``: :math:`\mathscr{D}=\mathbf{L}^N`, * ``'sobolev'``: :math:`\mathscr{D}=(\alpha^2 \mathrm{Id}-\mathbf{L})^N`, with :math:`\alpha\in\mathbb{R}`, * ``'polynomial'``: :math:`\mathscr{D}=\sum_{n=0}^N \alpha_n \mathbf{L}^n`, with :math:`\{\alpha_0,\ldots,\alpha_N\} \subset\mathbb{R}`, where :math:`\mathbf{L}` is the :py:class:`~pycgsp.linop.diff.GraphLaplacian` operator. kwargs: Any Additional arguments depending on the value of ``kind``: * ``'iterated'``: ``kwargs={order: int}`` where ``order`` defines the exponent :math:`N`. * ``'sobolev'``: ``kwargs={order: int, constant: float}`` where ``order`` defines the exponent :math:`N` and ``constant`` the scalar :math:`\alpha\in\mathbb{R}`. * ``'polynomial'``: ``kwargs={coeffs: Union[np.ndarray, list, tuple]}`` where ``coeffs`` is an array containing the coefficients :math:`\{\alpha_0,\ldots,\alpha_N\} \subset\mathbb{R}`. Raises ------ AttributeError If ``Graph.L`` does not exist. NotImplementedError If ``Graph.lap_type`` is 'combinatorial'. NotImplementedError If ``kind`` is not 'iterated', 'sobolev' or 'polynomial'. Examples -------- .. plot:: import numpy as np from pygsp.graphs import Ring from pycgsp.linop.diff import GeneralisedGraphLaplacian np.random.seed(1) G = Ring(N=32, k=4) G.compute_laplacian(lap_type='normalized') G.set_coordinates(kind='spring') x = np.arange(G.N) signal = np.piecewise(x, [x < G.N//3, (x >= G.N//3) * (x< 2 * G.N//3), x>=2 * G.N//3], [lambda x: -x, lambda x: 3 * x - 4 * G.N//3, lambda x: -0.5 * x + G.N]) Dop = GeneralisedGraphLaplacian(Graph=G, kind='polynomial', coeffs=[1,-1,2]) gen_lap = Dop * 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(gen_lap, ax=ax, backend='matplotlib') plt.title('Generalized Laplacian of signal') plt.axis('equal') Notes ----- The ``GeneralisedGraphLaplacian`` operator is self-adjoint. See Also -------- :py:class:`~pycgsp.linop.diff.GraphLaplacian` """ 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 = Graph.L.tocsc() LapOp = SparseLinearOperator(L, is_symmetric=True) if kind == 'iterated': N = kwargs['order'] Dgen = LapOp ** N elif kind == 'sobolev': I = IdentityOperator(size=LapOp.shape[0] * LapOp.shape[1]) alpha = kwargs['constant'] N = kwargs['order'] Dgen = ((alpha ** 2) * I - LapOp) ** N elif kind == 'polynomial': coeffs = kwargs['coeffs'] Dgen = PolynomialLinearOperator(LinOp=LapOp, coeffs=coeffs) else: raise NotImplementedError( 'Supported generalised Laplacian types are: iterated, sobolev, polynomial.') return Dgen