PyGSP: Graph Signal Processing in Python¶

The PyGSP is a Python package to ease Signal Processing on Graphs. It is a free software, distributed under the BSD license, and available on PyPI. The documentation is available on Read the Docs and development takes place on GitHub. (A Matlab counterpart exists.)

The PyGSP facilitates a wide variety of operations on graphs, like computing their Fourier basis, filtering or interpolating signals, plotting graphs, signals, and filters. Its core is spectral graph theory, and many of the provided operations scale to very large graphs. The package includes a wide range of graphs, from point clouds like the Stanford bunny and the Swiss roll; to networks like the Minnesota road network; to models for generating random graphs like stochastic block models, sensor networks, Erdős–Rényi model, Barabási-Albert model; to simple graphs like the path, the ring, and the grid. Many filter banks are also provided, e.g. various wavelets like the Mexican hat, Meyer, Half Cosine; some low-pass filters like the heat kernel and the exponential window; and Gabor filters. Despite all the pre-defined models, you can easily use a custom graph by defining its adjacency matrix, and a custom filter bank by defining a set of functions in the spectral domain.

The following demonstrates how to instantiate a graph and a filter, the two main objects of the package.

>>> from pygsp import graphs, filters
>>> G = graphs.Logo()
>>> G.estimate_lmax()
>>> g = filters.Heat(G, tau=100)

Let’s now create a graph signal: a set of three Kronecker deltas for that example. We can now look at one step of heat diffusion by filtering the deltas with the above defined filter. Note how the diffusion follows the local structure!

>>> import numpy as np
>>> DELTAS = [20, 30, 1090]
>>> s = np.zeros(G.N)
>>> s[DELTAS] = 1
>>> s = g.filter(s)
>>> G.plot_signal(s, highlight=DELTAS, backend='matplotlib')

You can try it online, look at the tutorials to learn how to use it, or look at the reference guide for an exhaustive documentation of the API. Enjoy the package!

Installation¶

The PyGSP is available on PyPI:

$ pip install pygsp

Note that you will need a recent version of pip and setuptools. Please run pip install --upgrade pip setuptools if you get any installation error.

Contributing¶

See the guidelines for contributing in CONTRIBUTING.rst.

Acknowledgments¶

The PyGSP was started in 2014 as an academic open-source project for research purpose at the EPFL LTS2 laboratory. This project has been partly funded by the Swiss National Science Foundation under grant 200021_154350 “Towards Signal Processing on Graphs”.

If you are using the library for your research, for the sake of reproducibility, please cite the version you used as indexed by Zenodo.

Tutorials¶

The following are some tutorials which explain how to use the toolbox and show how to use it to solve some problems.

Introduction to the PyGSP¶

This tutorial will show you the basic operations of the toolbox. After installing the package with pip, start by opening a python shell, e.g. a Jupyter notebook, and import the PyGSP. We will also need NumPy to create matrices and arrays.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pygsp import graphs, filters, plotting

We then set default plotting parameters. We’re using the matplotlib backend to embed plots in this tutorial. The pyqtgraph backend is best suited for interactive visualization.

>>> plotting.BACKEND = 'matplotlib'
>>> plt.rcParams['figure.figsize'] = (10, 5)

Graphs¶

Most likely, the first thing you would like to do is to create a graph. In this toolbox, a graph is encoded as an adjacency, or weight, matrix. That is because it’s the most efficient representation to deal with when using spectral methods. As such, you can construct a graph from any adjacency matrix as follows.

>>> rs = np.random.RandomState(42)  # Reproducible results.
>>> W = rs.uniform(size=(30, 30))  # Full graph.
>>> W[W < 0.93] = 0  # Sparse graph.
>>> W = W + W.T  # Symmetric graph.
>>> np.fill_diagonal(W, 0)  # No self-loops.
>>> G = graphs.Graph(W)
>>> print('{} nodes, {} edges'.format(G.N, G.Ne))
30 nodes, 60 edges

The pygsp.graphs.Graph class we just instantiated is the base class for all graph objects, which offers many methods and attributes.

Given, a graph object, we can test some properties.

>>> G.is_connected()
True
>>> G.is_directed()
False

We can retrieve our weight matrix, which is stored in a sparse format.

>>> (G.W == W).all()
True
>>> type(G.W)
<class 'scipy.sparse.lil.lil_matrix'>

We can access the graph Laplacian

>>> # The graph Laplacian (combinatorial by default).
>>> G.L.shape
(30, 30)

We can also compute and get the graph Fourier basis (see below).

>>> G.compute_fourier_basis()
>>> G.U.shape
(30, 30)

Or the graph differential operator, useful to e.g. compute the gradient or smoothness of a signal.

>>> G.compute_differential_operator()
>>> G.D.shape
(60, 30)

Note

Note that we called pygsp.graphs.Graph.compute_fourier_basis() and pygsp.graphs.Graph.compute_differential_operator() before accessing the Fourier basis pygsp.graphs.Graph.U and the differential operator pygsp.graphs.Graph.D. Doing so is however not mandatory as those matrices would have been computed when requested (lazy evaluation). Omitting to call the compute functions does print a warning to tell you that a potentially heavy computation is taking place under the hood (that’s also the reason those matrices are not computed when the graph object is instantiated). It is thus encouraged to call them so that you are aware of the involved computations.

To be able to plot a graph, we need to embed its nodes in a 2D or 3D space. While most included graph models define these coordinates, the graph we just created do not. We only passed a weight matrix after all. Let’s set some coordinates with pygsp.graphs.Graph.set_coordinates() and plot our graph.

>>> G.set_coordinates('ring2D')
>>> G.plot()

While we created our first graph ourselves, many standard models of graphs are implemented as subclasses of the Graph class and can be easily instantiated. Check the pygsp.graphs module to get a list of them and learn more about the Graph object.

Fourier basis¶

As in classical signal processing, the Fourier transform plays a central role in graph signal processing. Getting the Fourier basis is however computationally intensive as it needs to fully diagonalize the Laplacian. While it can be used to filter signals on graphs, a better alternative is to use one of the fast approximations (see pygsp.filters.Filter.filter()). Let’s compute it nonetheless to visualize the eigenvectors of the Laplacian. Analogous to classical Fourier analysis, they look like sinuses on the graph. Let’s plot the second and third eigenvectors (the first is constant). Those are graph signals, i.e. functions \(s: \mathcal{V} \rightarrow \mathbb{R}^d\) which assign a set of values (a vector in \(\mathbb{R}^d\)) at every node \(v \in \mathcal{V}\) of the graph.

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>>
>>> fig, axes = plt.subplots(1, 2, figsize=(10, 3))
>>> for i, ax in enumerate(axes):
...     G.plot_signal(G.U[:, i+1], vertex_size=30, ax=ax)
...     _ = ax.set_title('Eigenvector {}'.format(i+2))
...     ax.set_axis_off()
>>> fig.tight_layout()

The parallel with classical signal processing is best seen on a ring graph, where the graph Fourier basis is equivalent to the classical Fourier basis. The following plot shows some eigenvectors drawn on a 1D and 2D embedding of the ring graph. While the signals are easier to interpret on a 1D plot, the 2D plot best represents the graph.

>>> G2 = graphs.Ring(N=50)
>>> G2.compute_fourier_basis()
>>> fig, axes = plt.subplots(1, 2, figsize=(10, 4))
>>> G2.plot_signal(G2.U[:, 4], ax=axes[0])
>>> G2.set_coordinates('line1D')
>>> G2.plot_signal(G2.U[:, 1:4], ax=axes[1])
>>> fig.tight_layout()

Filters¶

To filter signals on graphs, we need to define filters. They are represented in the toolbox by the pygsp.filters.Filter class. Filters are usually defined in the spectral domain. Given the transfer function

\[g(x) = \frac{1}{1 + \tau x},\]

let’s define and plot that low-pass filter:

>>> tau = 1
>>> def g(x):
...     return 1. / (1. + tau * x)
>>> g = filters.Filter(G, g)
>>>
>>> fig, ax = plt.subplots()
>>> g.plot(plot_eigenvalues=True, ax=ax)
>>> _ = ax.set_title('Filter frequency response')

The filter is plotted along all the spectrum of the graph. The black crosses are the eigenvalues of the Laplacian. They are the points where the continuous filter will be evaluated to create a discrete filter.

Note

You can put multiple functions in a list to define a filter bank!

Note

The pygsp.filters module implements various standard filters.

Let’s create a graph signal and add some random noise.

>>> # Graph signal: each letter gets a different value + additive noise.
>>> s = np.zeros(G.N)
>>> s[G.info['idx_g']-1] = -1
>>> s[G.info['idx_s']-1] = 0
>>> s[G.info['idx_p']-1] = 1
>>> s += rs.uniform(-0.5, 0.5, size=G.N)

We can now try to denoise that signal by filtering it with the above defined low-pass filter.

>>> s2 = g.filter(s)
>>>
>>> fig, axes = plt.subplots(1, 2, figsize=(10, 3))
>>> G.plot_signal(s, vertex_size=30, ax=axes[0])
>>> _ = axes[0].set_title('Noisy signal')
>>> axes[0].set_axis_off()
>>> G.plot_signal(s2, vertex_size=30, ax=axes[1])
>>> _ = axes[1].set_title('Cleaned signal')
>>> axes[1].set_axis_off()
>>> fig.tight_layout()

While the noise is largely removed thanks to the filter, some energy is diffused between the letters. This is the typical behavior of a low-pass filter.

So here are the basics for the PyGSP. Please check the other tutorials and the reference guide for more. Enjoy!

Note

Please see the review article The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains for an overview of the methods this package leverages.

Introduction to spectral graph wavelets¶

This tutorial will show you how to easily construct a wavelet frame, a kind of filter bank, and apply it to a signal. This tutorial will walk you into computing the wavelet coefficients of a graph, visualizing filters in the vertex domain, and using the wavelets to estimate the curvature of a 3D shape.

As usual, we first have to import some packages.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pygsp import graphs, filters, plotting, utils

Then we can load a graph. The graph we’ll use is a nearest-neighbor graph of a point cloud of the Stanford bunny. It will allow us to get interesting visual results using wavelets.

>>> G = graphs.Bunny()

Note

At this stage we could compute the Fourier basis using pygsp.graphs.Graph.compute_fourier_basis(), but this would take some time, and can be avoided with a Chebychev polynomials approximation to graph filtering. See the documentation of the pygsp.filters.Filter.filter() filtering function and [hammond2011wavelets] for details on how it is down.

Simple filtering: heat diffusion¶

Before tackling wavelets, let’s observe the effect of a single filter on a graph signal. We first design a few heat kernel filters, each with a different scale.

>>> taus = [10, 25, 50]
>>> g = filters.Heat(G, taus)

Let’s create a signal as a Kronecker delta located on one vertex, e.g. the vertex 20. That signal is our heat source.

>>> s = np.zeros(G.N)
>>> DELTA = 20
>>> s[DELTA] = 1

We can now simulate heat diffusion by filtering our signal s with each of our heat kernels.

>>> s = g.filter(s, method='chebyshev')

And finally plot the filtered signal showing heat diffusion at different scales.

>>> fig = plt.figure(figsize=(10, 3))
>>> for i in range(g.Nf):
...     ax = fig.add_subplot(1, g.Nf, i+1, projection='3d')
...     G.plot_signal(s[:, i], colorbar=False, ax=ax)
...     title = r'Heat diffusion, $\tau={}$'.format(taus[i])
...     _ = ax.set_title(title)
...     ax.set_axis_off()
>>> fig.tight_layout()

Note

The pygsp.filters.Filter.localize() method can be used to visualize a filter in the vertex domain instead of doing it manually.

Visualizing wavelets atoms¶

Let’s now replace the Heat filter by a filter bank of wavelets. We can create a filter bank using one of the predefined filters, such as pygsp.filters.MexicanHat to design a set of Mexican hat wavelets.

>>> g = filters.MexicanHat(G, Nf=6)  # Nf = 6 filters in the filter bank.

Then plot the frequency response of those filters.

>>> fig, ax = plt.subplots(figsize=(10, 5))
>>> g.plot(ax=ax)
>>> _ = ax.set_title('Filter bank of mexican hat wavelets')

Note

We can see that the wavelet atoms are stacked on the low frequency part of the spectrum. A better coverage could be obtained by adapting the filter bank with pygsp.filters.WarpedTranslates or by using another filter bank like pygsp.filters.Itersine.

We can visualize the atoms as we did with the heat kernel, by filtering a Kronecker delta placed at one specific vertex.

>>> s = g.localize(DELTA)
>>>
>>> fig = plt.figure(figsize=(10, 2.5))
>>> for i in range(3):
...     ax = fig.add_subplot(1, 3, i+1, projection='3d')
...     G.plot_signal(s[:, i], ax=ax)
...     _ = ax.set_title('Wavelet {}'.format(i+1))
...     ax.set_axis_off()
>>> fig.tight_layout()

Curvature estimation¶

As a last and more applied example, let us try to estimate the curvature of the underlying 3D model by only using spectral filtering on the nearest-neighbor graph formed by its point cloud.

A simple way to accomplish that is to use the coordinates map \([x, y, z]\) and filter it using the above defined wavelets. Doing so gives us a 3-dimensional signal \([g_i(L)x, g_i(L)y, g_i(L)z], \ i \in [0, \ldots, N_f]\) which describes variation along the 3 coordinates.

>>> s = G.coords
>>> s = g.filter(s)

The curvature is then estimated by taking the \(\ell_1\) or \(\ell_2\) norm across the 3D position.

>>> s = np.linalg.norm(s, ord=2, axis=1)

Let’s finally plot the result to observe that we indeed have a measure of the curvature at different scales.

>>> fig = plt.figure(figsize=(10, 7))
>>> for i in range(4):
...     ax = fig.add_subplot(2, 2, i+1, projection='3d')
...     G.plot_signal(s[:, i], ax=ax)
...     title = 'Curvature estimation (scale {})'.format(i+1)
...     _ = ax.set_title(title)
...     ax.set_axis_off()
>>> fig.tight_layout()

Optimization problems: graph TV vs. Tikhonov regularization¶

Description¶

Modern signal processing often involves solving an optimization problem. Graph signal processing (GSP) consists roughly of working with linear operators defined by a graph (e.g., the graph Laplacian). The setting up of an optimization problem in the graph context is often then simply a matter of identifying which graph-defined linear operator is relevant to be used in a regularization and/or fidelity term.

This tutorial focuses on the problem of recovering a label signal on a graph from subsampled and noisy data, but most information should be fairly generally applied to other problems as well.

>>> import numpy as np
>>> from pygsp import graphs, plotting
>>>
>>> # Create a random sensor graph
>>> G = graphs.Sensor(N=256, distribute=True, seed=42)
>>> G.compute_fourier_basis()
>>>
>>> # Create label signal
>>> label_signal = np.copysign(np.ones(G.N), G.U[:, 3])
>>>
>>> G.plot_signal(label_signal)

The first figure shows a plot of the original label signal, that we wish to recover, on the graph.

>>> rs = np.random.RandomState(42)
>>>
>>> # Create the mask
>>> M = rs.rand(G.N)
>>> M = (M > 0.6).astype(float)  # Probability of having no label on a vertex.
>>>
>>> # Applying the mask to the data
>>> sigma = 0.1
>>> subsampled_noisy_label_signal = M * (label_signal + sigma * rs.standard_normal(G.N))
>>>
>>> G.plot_signal(subsampled_noisy_label_signal)

This figure shows the label signal on the graph after the application of the subsampling mask and the addition of noise. The label of more than half of the vertices has been set to \(0\).

Since the problem is ill-posed, we will use some regularization to reach a solution that is more in tune with what we expect a label signal to look like. We will compare two approaches, but they are both based on measuring local differences on the label signal. Those differences are essentially an edge signal: to each edge we can associate the difference between the label signals of its associated nodes. The linear operator that does such a mapping is called the graph gradient \(\nabla_G\), and, fortunately for us, it is available under the pygsp.graphs.Graph.D (for differential) attribute of any pygsp.graphs graph.

The reason for measuring local differences comes from prior knowledge: we assume label signals don’t vary too much locally. The precise measure of such variation is what distinguishes the two regularization approaches we’ll use.

The first one, shown below, is called graph total variation (TV) regularization. The quadratic fidelity term is multiplied by a regularization constant \(\gamma\) and its goal is to force the solution to stay close to the observed labels \(b\). The \(\ell_1\) norm of the action of the graph gradient is what’s called the graph TV. We will see that it promotes piecewise-constant solutions.

\[\operatorname*{arg\,min}_x \|\nabla_G x\|_1 + \gamma \|Mx-b\|_2^2\]

The second approach, called graph Tikhonov regularization, is to use a smooth (differentiable) quadratic regularizer. A consequence of this choice is that the solution will tend to have smoother transitions. The quadratic fidelity term is still the same.

\[\operatorname*{arg\,min}_x \|\nabla_G x\|_2^2 + \gamma \|Mx-b\|_2^2\]

Results and code¶

For solving the optimization problems we’ve assembled, you will need a numerical solver package. This part is implemented in this tutorial with the pyunlocbox, which is based on proximal splitting algorithms. Check also its documentation for more information about the parameters used here.

We start with the graph TV regularization. We will use the pyunlocbox.solvers.mlfbf solver from pyunlocbox. It is a primal-dual solver, which means for our problem that the regularization term will be written in terms of the dual variable \(u = \nabla_G x\), and the graph gradient \(\nabla_G\) will be passed to the solver as the primal-dual map. The value of \(3.0\) for the regularization parameter \(\gamma\) was chosen on the basis of the visual appeal of the returned solution.

>>> import pyunlocbox
>>>
>>> # Set the functions in the problem
>>> gamma = 3.0
>>> d = pyunlocbox.functions.dummy()
>>> r = pyunlocbox.functions.norm_l1()
>>> f = pyunlocbox.functions.norm_l2(w=M, y=subsampled_noisy_label_signal,
...                                  lambda_=gamma)
>>>
>>> # Define the solver
>>> G.compute_differential_operator()
>>> L = G.D.toarray()
>>> step = 0.999 / (1 + np.linalg.norm(L))
>>> solver = pyunlocbox.solvers.mlfbf(L=L, step=step)
>>>
>>> # Solve the problem
>>> x0 = subsampled_noisy_label_signal.copy()
>>> prob1 = pyunlocbox.solvers.solve([d, r, f], solver=solver,
...                                  x0=x0, rtol=0, maxit=1000)
Solution found after 1000 iterations:
    objective function f(sol) = 2.024594e+02
    stopping criterion: MAXIT
>>>
>>> G.plot_signal(prob1['sol'])

This figure shows the label signal recovered by graph total variation regularization. We can confirm here that this sort of regularization does indeed promote piecewise-constant solutions.

>>> # Set the functions in the problem
>>> r = pyunlocbox.functions.norm_l2(A=L, tight=False)
>>>
>>> # Define the solver
>>> step = 0.999 / np.linalg.norm(np.dot(L.T, L) + gamma * np.diag(M), 2)
>>> solver = pyunlocbox.solvers.gradient_descent(step=step)
>>>
>>> # Solve the problem
>>> x0 = subsampled_noisy_label_signal.copy()
>>> prob2 = pyunlocbox.solvers.solve([r, f], solver=solver,
...                                  x0=x0, rtol=0, maxit=1000)
Solution found after 1000 iterations:
    objective function f(sol) = 9.555135e+01
    stopping criterion: MAXIT
>>>
>>> G.plot_signal(prob2['sol'])

This last figure shows the label signal recovered by Tikhonov regularization. As expected, the recovered label signal has smoother transitions than the one obtained by graph TV regularization.

Graph multiresolution: Kron pyramid¶

In this demonstration file, we show how to reduce a graph using the PyGSP. Then we apply the pyramid to simple signal. To start open a python shell (IPython is recommended here) and import the required packages. You would probably also import numpy as you will need it to create matrices and arrays.

>>> import numpy as np
>>> from pygsp import graphs, reduction

For this demo we will be using a sensor graph with 512 nodes.

>>> G = graphs.Sensor(512, distribute=True)
>>> G.compute_fourier_basis()

The function graph_multiresolution computes the graph pyramid for you:

>>> levels = 5
>>> Gs = reduction.graph_multiresolution(G, levels, sparsify=False)

Next, we will compute the fourier basis of our different graph layers:

>>> for gr in Gs:
...     gr.compute_fourier_basis()

Those that were already computed are returning with an error, meaning that nothing happened. Let’s now create two signals and a filter, resp f, f2 and g:

>>> f = np.ones((G.N))
>>> f[np.arange(G.N//2)] = -1
>>> f = f + 10*Gs[0].U[:, 7]

>>> f2 = np.ones((G.N, 2))
>>> f2[np.arange(G.N//2)] = -1

>>> g = [lambda x: 5./(5 + x)]

We will run the analysis of the two signals on the pyramid and obtain a coarse approximation for each layer, with decreasing number of nodes. Additionally, we will also get prediction errors at each node at every layer.

>>> ca, pe = reduction.pyramid_analysis(Gs, f, h_filters=g, method='exact')
>>> ca2, pe2 = reduction.pyramid_analysis(Gs, f2, h_filters=g, method='exact')

Given the pyramid, the coarsest approximation and the prediction errors, we will now reconstruct the original signal on the full graph.

>>> f_pred, _ = reduction.pyramid_synthesis(Gs, ca[levels], pe, method='exact')
>>> f_pred2, _ = reduction.pyramid_synthesis(Gs, ca2[levels], pe2, method='exact')

Here are the final errors for each signal after reconstruction.

>>> err = np.linalg.norm(f_pred-f)/np.linalg.norm(f)
>>> err2 = np.linalg.norm(f_pred2-f2)/np.linalg.norm(f2)
>>> assert (err < 1e-10) & (err2 < 1e-10)

Reference guide¶

The pygsp package is mainly organized around the following two modules:

graphs to create and manipulate various kinds of graphs,
filters to create and manipulate various graph filters.

Moreover, the following modules provide additional functionality:

plotting to plot,
reduction to reduce a graph while keeping its structure,
features to compute features on graphs,
optimization to help solving convex optimization problems,
utils for various utilities.

Graphs¶

The pygsp.graphs module implements the graph class hierarchy. A graph object is either constructed from an adjacency matrix, or by instantiating one of the built-in graph models.

Interface¶

The Graph base class allows to construct a graph object from any adjacency matrix and provides a common interface to that object. Derived classes then allows to instantiate various standard graph models.

Matrix operators¶

`Graph.W`
`Graph.L`
`Graph.U`	Fourier basis (eigenvectors of the Laplacian).
`Graph.D`	Differential operator (for gradient and divergence).

Checks¶

`Graph.check_weights`(self)	Check the characteristics of the weights matrix.
`Graph.is_connected`(self[, recompute])	Check the strong connectivity of the graph (cached).
`Graph.is_directed`(self[, recompute])	Check if the graph has directed edges (cached).

Attributes computation¶

`Graph.compute_laplacian`(self[, lap_type])	Compute a graph Laplacian.
`Graph.estimate_lmax`(self[, recompute])	Estimate the Laplacian’s largest eigenvalue (cached).
`Graph.compute_fourier_basis`(self[, recompute])	Compute the Fourier basis of the graph (cached).
`Graph.compute_differential_operator`(self)	Compute the graph differential operator (cached).

Differential operators¶

`Graph.grad`(self, s)	Compute the gradient of a graph signal.
`Graph.div`(self, s)	Compute the divergence of a graph signal.

Localization¶

`Graph.modulate`(self, f, k)	Modulate the signal f to the frequency k.
`Graph.translate`(self, f, i)	Translate the signal f to the node i.

Transforms (frequency and vertex-frequency)¶

`Graph.gft`(self, s)	Compute the graph Fourier transform.
`Graph.igft`(self, s_hat)	Compute the inverse graph Fourier transform.
`Graph.gft_windowed`(self, g, f[, lowmemory])	Windowed graph Fourier transform.
`Graph.gft_windowed_gabor`(self, s, k)	Gabor windowed graph Fourier transform.
`Graph.gft_windowed_normalized`(self, g, f[, …])	Normalized windowed graph Fourier transform.

Plotting¶

`Graph.plot`(self, \\kwargs)	Plot the graph.
`Graph.plot_signal`(self, signal, \\kwargs)	Plot a signal on that graph.
`Graph.plot_spectrogram`(self, \\kwargs)	Plot the graph’s spectrogram.

Others¶

`Graph.get_edge_list`(self)	Return an edge list, an alternative representation of the graph.
`Graph.set_coordinates`(self[, kind])	Set node’s coordinates (their position when plotting).
`Graph.subgraph`(self, ind)	Create a subgraph given indices.
`Graph.extract_components`(self)	Split the graph into connected components.

Graph models¶

`Airfoil`(**kwargs)	Airfoil graph.
`BarabasiAlbert`([N, m0, m, seed])	Barabasi-Albert preferential attachment.
`Comet`([N, k])	Comet graph.
`Community`([N, Nc, min_comm, min_deg, …])	Community graph.
`DavidSensorNet`([N, seed])	Sensor network.
`ErdosRenyi`([N, p, directed, self_loops, …])	Erdos Renyi graph.
`FullConnected`([N])	Fully connected graph.
`Grid2d`([N1, N2])	2-dimensional grid graph.
`Logo`(**kwargs)	GSP logo.
`LowStretchTree`([k])	Low stretch tree.
`Minnesota`([connect])	Minnesota road network (from MatlabBGL).
`Path`([N])	Path graph.
`RandomRegular`([N, k, maxIter, seed])	Random k-regular graph.
`RandomRing`([N, seed])	Ring graph with randomly sampled nodes.
`Ring`([N, k])	K-regular ring graph.
`Sensor`([N, Nc, regular, n_try, distribute, …])	Random sensor graph.
`StochasticBlockModel`([N, k, z, M, p, q, …])	Stochastic Block Model (SBM).
`SwissRoll`([N, a, b, dim, thresh, s, noise, …])	Sampled Swiss roll manifold.
`Torus`([Nv, Mv])	Sampled torus manifold.

Nearest-neighbors graphs constructed from point clouds¶

`NNGraph`(Xin[, NNtype, use_flann, center, …])	Nearest-neighbor graph from given point cloud.
`Bunny`(**kwargs)	Stanford bunny (NN-graph).
`Cube`([radius, nb_pts, nb_dim, sampling, seed])	Hyper-cube (NN-graph).
`ImgPatches`(img[, patch_shape])	NN-graph between patches of an image.
`Grid2dImgPatches`(img[, aggregate])	Union of a patch graph with a 2D grid graph.
`Sphere`([radius, nb_pts, nb_dim, sampling, seed])	Spherical-shaped graph (NN-graph).
`TwoMoons`([moontype, dim, sigmag, N, sigmad, …])	Two Moons (NN-graph).

class pygsp.graphs.Graph(W, gtype='unknown', lap_type='combinatorial', coords=None, plotting={})[source]¶

Base graph class.

Provide a common interface (and implementation) to graph objects.
Can be instantiated to construct custom graphs from a weight matrix.
Initialize attributes for derived classes.

Parameters

Wsparse matrix or ndarray: The weight matrix which encodes the graph.
gtypestring: Graph type, a free-form string to help us recognize the kind of graph we are dealing with (default is ‘unknown’).
lap_type‘combinatorial’, ‘normalized’: The type of Laplacian to be computed by compute_laplacian() (default is ‘combinatorial’).
coordsndarray: Vertices coordinates (default is None).
plottingdict: Plotting parameters.

Examples

>>> W = np.arange(4).reshape(2, 2)
>>> G = graphs.Graph(W)

Attributes

Nint: the number of nodes / vertices in the graph.
Neint: the number of edges / links in the graph, i.e. connections between nodes.
Wsparse matrix: the weight matrix which contains the weights of the connections. It is represented as an N-by-N matrix of floats. \(W_{i,j} = 0\) means that there is no direct connection from i to j.
gtypestring: the graph type is a short description of the graph object designed to help sorting the graphs.
Lsparse matrix: the graph Laplacian, an N-by-N matrix computed from W.
lap_type‘normalized’, ‘combinatorial’: the kind of Laplacian that was computed by compute_laplacian().
coordsndarray: vertices coordinates in 2D or 3D space. Used for plotting only. Default is None.
plottingdict: plotting parameters.

check_weights(self)[source]¶

Check the characteristics of the weights matrix.

Returns

A dict of bools containing informations about the matrix
has_inf_valbool: True if the matrix has infinite values else false
has_nan_valuebool: True if the matrix has a “not a number” value else false
is_not_squarebool: True if the matrix is not square else false
diag_is_not_zerobool: True if the matrix diagonal has not only zeros else false

Examples

>>> W = np.arange(4).reshape(2, 2)
>>> G = graphs.Graph(W)
>>> cw = G.check_weights()
>>> cw == {'has_inf_val': False, 'has_nan_value': False,
...        'is_not_square': False, 'diag_is_not_zero': True}
True

compute_differential_operator(self)¶

Compute the graph differential operator (cached).

The differential operator is a matrix such that

\[L = D^T D,\]

where \(D\) is the differential operator and \(L\) is the graph Laplacian. It is used to compute the gradient and the divergence of a graph signal, see grad() and div().

The result is cached and accessible by the D property.

Filters¶

The pygsp.filters module implements methods used for filtering and defines commonly used filters that can be applied to pygsp.graphs. A filter is associated to a graph and is defined with one or several functions. We define by filter bank a list of filters, usually centered around different frequencies, applied to a single graph.

Interface¶

The Filter base class implements a common interface to all filters:

`Filter.evaluate`(self, x)	Evaluate the kernels at given frequencies.
`Filter.filter`(self, s[, method, order])	Filter signals (analysis or synthesis).
`Filter.analyze`(self, s[, method, order])	Convenience alias to `filter()`.
`Filter.synthesize`(self, s[, method, order])	Convenience wrapper around `filter()`.
`Filter.compute_frame`(self, \\kwargs)	Compute the associated frame.
`Filter.estimate_frame_bounds`(self[, min, …])	Estimate lower and upper frame bounds.
`Filter.plot`(self, \\kwargs)	Plot the filter bank’s frequency response.
`Filter.localize`(self, i, \\kwargs)	Localize the kernels at a node (to visualize them).

Filters¶

Then, derived classes implement various common graph filters.

Filter banks of N filters

`Abspline`(G[, Nf, lpfactor, scales])	Design an A B cubic spline wavelet filter bank.
`Gabor`(G, k)	Design a Gabor filter bank.
`HalfCosine`(G[, Nf])	Design an half cosine filter bank (tight frame).
`Itersine`(G[, Nf, overlap])	Design an itersine filter bank (tight frame).
`MexicanHat`(G[, Nf, lpfactor, scales, normalize])	Design a filter bank of Mexican hat wavelets.
`Meyer`(G[, Nf, scales])	Design a filter bank of Meyer wavelets (tight frame).
`SimpleTight`(G[, Nf, scales])	Design a simple tight frame filter bank (tight frame).

Filter banks of 2 filters: a low pass and a high pass

`Regular`(G[, d])	Design 2 filters with the regular construction (tight frame).
`Held`(G[, a])	Design 2 filters with the Held construction (tight frame).
`Simoncelli`(G[, a])	Design 2 filters with the Simoncelli construction (tight frame).
`Papadakis`(G[, a])	Design 2 filters with the Papadakis construction (tight frame).

Low pass filters

`Heat`(G[, tau, normalize])	Design a filter bank of heat kernels.
`Expwin`(G[, bmax, a])	Design an exponential window filter.

Approximations¶

Moreover, two approximation methods are provided for fast filtering. The computational complexity of filtering with those approximations is linear with the number of edges. The complexity of the exact solution, which is to use the Fourier basis, is quadratic with the number of nodes (without taking into account the cost of the necessary eigendecomposition of the graph Laplacian).

Chebyshev polynomials

`compute_cheby_coeff`(f[, m, N])	Compute Chebyshev coefficients for a Filterbank.
`compute_jackson_cheby_coeff`(filter_bounds, …)	To compute the m+1 coefficients of the polynomial approximation of an ideal band-pass between a and b, between a range of values defined by lambda_min and lambda_max.
`cheby_op`(G, c, signal, \\kwargs)	Chebyshev polynomial of graph Laplacian applied to vector.
`cheby_rect`(G, bounds, signal, \\kwargs)	Fast filtering using Chebyshev polynomial for a perfect rectangle filter.

Lanczos algorithm

`lanczos`(A, order, x)	TODO short description
`lanczos_op`(f, s[, order])	Perform the lanczos approximation of the signal s.

Plotting¶

The pygsp.plotting module implements functionality to plot PyGSP objects with a pyqtgraph or matplotlib drawing backend (which can be controlled by the BACKEND constant or individually for each plotting call):

graphs from pygsp.graphs with plot_graph(), plot_spectrogram(), and plot_signal(),
filters from pygsp.filters with plot_filter().

pygsp.plotting.BACKEND¶: Indicates which drawing backend to use if none are provided to the plotting functions. Should be either ‘matplotlib’ or ‘pyqtgraph’. In general pyqtgraph is better for interactive exploration while matplotlib is better at generating figures to be included in papers or elsewhere.

Reduction¶

The pygsp.reduction module implements functionalities for the reduction of graphs’ vertex set while keeping the graph structure.

`tree_multiresolution`(G, Nlevel[, …])	Compute a multiresolution of trees
`graph_multiresolution`(G, levels[, sparsify, …])	Compute a pyramid of graphs (by Kron reduction).
`kron_reduction`(G, ind)	Compute the Kron reduction.
`pyramid_analysis`(Gs, f, \\kwargs)	Compute the graph pyramid transform coefficients.
`pyramid_synthesis`(Gs, cap, pe[, order])	Synthesize a signal from its pyramid coefficients.
`interpolate`(G, f_subsampled, keep_inds[, …])	Interpolate a graph signal.
`graph_sparsify`(M, epsilon[, maxiter])	Sparsify a graph (with Spielman-Srivastava).

Features¶

The pygsp.features module implements different feature extraction techniques based on pygsp.graphs and pygsp.filters.

Optimization¶

The pygsp.optimization module provides tools for convex optimization on graphs.

Utilities¶

The pygsp.utils module implements some utility functions used throughout the package.

Contributing¶

Contributions are welcome, and they are greatly appreciated! The development of this package takes place on GitHub. Issues, bugs, and feature requests should be reported there. Code and documentation can be improved by submitting a pull request. Please add documentation and tests for any new code.

The package can be set up (ideally in a virtual environment) for local development with the following:

$ git clone https://github.com/epfl-lts2/pygsp.git
$ pip install -U -e pygsp[alldeps,test,doc,pkg]

You can improve or add functionality in the pygsp folder, along with corresponding unit tests in pygsp/tests/test_*.py (with reasonable coverage) and documentation in doc/reference/*.rst. If you have a nice example to demonstrate the use of the introduced functionality, please consider adding a tutorial in doc/tutorials.

Do not forget to update README.rst and doc/history.rst with e.g. new features. The version number needs to be updated in setup.py and pyunlocbox/__init__.py.

After making any change, please check the style, run the tests, and build the documentation with the following (enforced by Travis CI):

$ make lint
$ make test
$ make doc

Check the generated coverage report at htmlcov/index.html to make sure the tests reasonably cover the changes you’ve introduced.

Making a release¶

Update the version number and release date in setup.py, pygsp/__init__.py and doc/history.rst.
Create a git tag with git tag -a v0.5.0 -m "PyGSP v0.5.0".
Push the tag to GitHub with git push github v0.5.0. The tag should now appear in the releases and tags tab.
Create a release on GitHub and select the created tag. A DOI should then be issued by Zenodo.
Go on Zenodo and fix the metadata if necessary.
Build the distribution with make dist and check that the dist/PyGSP-0.5.0.tar.gz source archive contains all required files. The binary wheel should be found as dist/PyGSP-0.5.0-py2.py3-none-any.whl.

Test the upload and installation process:

$ twine upload --repository-url https://test.pypi.org/legacy/ dist/*
$ pip install --index-url https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple pygsp

Log in as the LTS2 user.

Build and upload the distribution to the real PyPI with make release.

Repository organization¶

LICENSE.txt         Project license
*.rst               Important documentation
Makefile            Targets for make
setup.py            Meta information about package (published on PyPI)
.gitignore          Files ignored by the git revision control system
.travis.yml         Defines testing on Travis continuous integration

pygsp/              Contains the modules (the actual toolbox implementation)
 __init.py__        Load modules at package import
 *.py               One file per module

pygsp/tests/        Contains the test suites (will be distributed to end user)
 __init.py__        Load modules at package import
 test_*.py          One test suite per module
 test_docstrings.py Test the examples in the docstrings (reference doc)
 test_tutorials.py  Test the tutorials in doc/tutorials
 test_all.py        Launch all the tests (docstrings, tutorials, modules)

doc/                Package documentation
 conf.py            Sphinx configuration
 index.rst          Documentation entry page
 *.rst              Include doc files from root directory

doc/reference/      Reference documentation
 index.rst          Reference entry page
 *.rst              Only directives, the actual doc is alongside the code

doc/tutorials/
 index.rst          Tutorials entry page
 *.rst              One file per tutorial

History¶

0.5.1 (2017-12-15)¶

The focus of this release was to ease installation by not requiring non-standard scientific Python packages to be installed. It was mostly a maintenance release. A conda package is now available in conda-forge. Moreover, the package can now be tried online thanks to binder.

The core functionality of this package only depends on numpy and scipy. Dependencies which are only required for particular usages are included in the alldeps extra dependency list. The alldeps list allows users to install dependencies to enable all the features. Finally, those optional packages are only loaded when needed, not when the PyGSP is imported. A nice side-effect is that importing the PyGSP is now much faster!

The following packages were made optional dependencies: * scikit-image, as it is only used to build patch graphs from images. The

problem was that scikit-image does not provide a wheel for Windows and its build is painful and error-prone. Moreover, scikit-image has a lot of dependencies.

pyqtgrpah, PyQt5 / PySide and PyOpenGl, as they are only used for interactive visualization, which not many users need. The problem was that pyqtgraph requires (via PyQt5, PySide, PyOpenGL) OpenGL (libGL.so) to be installed.
matplotlib: while it is a standard package for any scientific or data science workflow, it’s not necessary for users who only want to process data without plotting graphs, signals and filters.
pyflann, as it is only used for approximate kNN. The problem was that the source distribution would not build for Windows. On conda-forge, (py)flann is not built for Windows either.

Moreover, matplotlib is now the default drawing backend. It’s well integrated with the Jupyter environment for scientific and data science workflows, and most use cases do not require an interactive visualization. The pyqtgraph is still available for interactivity.

0.5.0 (2017-10-06)¶

Generalized the analysis and synthesis methods into the filter method.
Signals are now rank-3 tensors of N_NODES x N_SIGNALS x N_FEATURES.
Filter.evaluate returns a 2D array instead of a list of vectors.
The differential operator was integrated in the Graph class, as the Fourier basis and the Laplacian were already.
Removed the operators package. Transforms and differential operators went to the Graph class, the localization operator to the Filter class. These are now easier to use. Reduction became its own module.
Graph object uses properties to delay the computation (lazy evaluation) of the Fourier basis (G.U, G.e, G.mu), the estimation of the largest eigenvalue (G.lmax) and the computation of the differential operator (G.D). A warning is issued if client code don’t trigger computations with G.compute_*.
Approximations for filtering have been moved in the filters package.
PointCloud object removed. Functionality integrated in Graph object.
data_handling module merged into utils.
Fourier basis computed with eigh instead of svd (faster).
estimate_lmax uses Lanczos instead of Arnoldi (symmetric sparse).
Add a seed parameter to all non-deterministic graphs and filters.
Filter.Nf indicates the number of filters in the filter bank.
Don’t check connectedness on graph creation (can take a lot of time).
Erdos-Renyi now implemented as SBM with 1 block.
Many bug fixes (e.g. Minnesota graph, Meyer filter bank, Heat filter, Mexican hat filter bank, Gabor filter bank).
All GitHub issues fixed.

Plotting:

Much better handling of plotting parameters.
With matplotlib backend, plots are shown by default .
Allows to set a default plotting backend as plotting.BACKEND = ‘pyqtgraph’.
qtg_default=False becomes backend=’matplotlib’
Added coordinates for path, ring, and randomring graphs.
Set good default plotting parameters for most graphs.
Allows to plot multiple filters in 1D with set_coordinates(‘line1D’).
Allows to pass existing matplotlib axes to the plotting functions.
Show colorbar with matplotlib.
Allows to set a 3D view point.
Eigenvalues shown as vertical lines instead of crosses.
Vertices can be highlighted, e.g. to show where filters where localized.

Documentation:

More comprehensive documentation. Notably math definitions for operators.
Most filters and graphs are plotted in the API documentation.
List all methods and models at the top with autosummary.
Useful package and module-level documentation.
Doctests don’t need to import numpy and the pygsp every time.
Figures are automatically generated when building the documentation.
Build on RTD with conda and matplotlib 2 (prettier plots).
Intro and wavelets tutorials were updated.
Reference guide is completely auto-generated from automodule.
Added contribution guidelines.
Documentation reorganization.
Check that hyperlinks are valid.

Tests and infrastructure:

Start test coverage analysis.
Much more comprehensive tests. Coverage increased from 40% to 80%. Many bugs were uncovered.
Always test with virtual X framebuffer to avoid the opening of lots of windows.
Tested on Python 2.7, 3.4, 3.5, 3.6.
Clean configuration files.
Not using tox anymore (too painful to maintain multiple Pythons locally).
Sort out installation of dependencies. Plotting should now work right away.
Completely migrate development on GitHub.

0.4.2 (2017-04-27)¶

Improve documentation.
Various fixes.

0.4.1 (2016-09-06)¶

Added routines to compute coordinates for the graphs.
Added fast filtering of ideal band-pass.
Implemented graph spectrograms.
Added the Barabási-Albert model for graphs.
Renamed PointClouds features.
Various fixes.

0.4.0 (2016-06-17)¶

0.3.3 (2016-01-27)¶

Refactoring graphs using object programming and fail safe checks.
Refactoring filters to use only the Graph object used at the construction of the filter for all operations.
Refactoring Graph pyramid to match MATLAB implementation.
Removal of default coordinates (all vertices on the origin) for graphs that do not possess spatial meaning.
Correction of minor issues on Python3+ imports.
Various fixes.
Finalizing demos for the documentation.

0.3.2 (2016-01-14)¶

0.3.1 (2016-01-12)¶

0.3.0 (2015-12-01)¶

0.2.1 (2015-10-20)¶

Fix bug on pip installation.
Update full documentation.

0.2.0 (2015-10-12)¶

Adding functionalities to match the content of the Matlab GSP Box.
First release of the PyGSP.

0.1.0 (2015-07-02)¶

Main features of the box are present most of the graphs and filters can be used.
The utils and operators modules also have most of their features implemented.

References¶

Chu97: F. R. K. Chung. Spectral Graph Theory. Vol. 92 of the CBMS Regional Conference Series in Mathematics, American Mathematical Society, 1997.