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 [HVG11] 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`()	Check the characteristics of the weights matrix.
`Graph.is_connected`([recompute])	Check the strong connectivity of the graph (cached).
`Graph.is_directed`([recompute])	Check if the graph has directed edges (cached).

Attributes computation¶

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

Differential operators¶

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

Localization¶

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

Transforms (frequency and vertex-frequency)¶

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

Plotting¶

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

Others¶

`Graph.get_edge_list`()	Return an edge list, an alternative representation of the graph.
`Graph.set_coordinates`([kind])	Set node’s coordinates (their position when plotting).
`Graph.subgraph`(ind)	Create a subgraph given indices.
`Graph.extract_components`()	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:

W : sparse matrix or ndarray

The weight matrix which encodes the graph.

gtype : string

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’).

coords : ndarray

Vertices coordinates (default is None).

plotting : dict

Plotting parameters.

Examples

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

Attributes

N	(int) the number of nodes / vertices in the graph.
Ne	(int) the number of edges / links in the graph, i.e. connections between nodes.
W	(sparse 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.
gtype	(string) the graph type is a short description of the graph object designed to help sorting the graphs.
L	(sparse 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()`.
coords	(ndarray) vertices coordinates in 2D or 3D space. Used for plotting only. Default is None.
plotting	(dict) plotting parameters.

check_weights()[source]¶

Check the characteristics of the weights matrix.

Returns:

A dict of bools containing informations about the matrix

has_inf_val : bool

True if the matrix has infinite values else false

has_nan_value : bool

True if the matrix has a “not a number” value else false

is_not_square : bool

True if the matrix is not square else false

diag_is_not_zero : bool

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()¶

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.

See also

grad: compute the gradient
div: compute the divergence

Examples

>>> G = graphs.Logo()
>>> G.N, G.Ne
(1130, 3131)
>>> G.compute_differential_operator()
>>> G.D.shape == (G.Ne, G.N)
True

compute_fourier_basis(recompute=False)¶

Compute the Fourier basis of the graph (cached).

The result is cached and accessible by the U, e, lmax, and mu properties.

Parameters:

recompute: bool

Force to recompute the Fourier basis if already existing.

Notes

‘G.compute_fourier_basis()’ computes a full eigendecomposition of the graph Laplacian $L$ such that:

\[L = U \Lambda U^*,\]

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

G.e is a vector of length 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.mu.

References

See [Chu97].

Examples

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

compute_laplacian(lap_type='combinatorial')[source]¶

Compute a graph Laplacian.

The result is accessible by the L attribute.

Parameters:

lap_type : ‘combinatorial’, ‘normalized’

The type of Laplacian to compute. Default is combinatorial.

Notes

For undirected graphs, the combinatorial Laplacian is defined as

\[L = D - W,\]

where $W$ is the weight matrix and $D$ the degree matrix, and the normalized Laplacian is defined as

\[L = I - D^{-1/2} W D^{-1/2},\]

where $I$ is the identity matrix.

Examples

>>> G = graphs.Sensor(50)
>>> G.L.shape
(50, 50)
>>>
>>> G.compute_laplacian('combinatorial')
>>> G.compute_fourier_basis()
>>> -1e-10 < G.e[0] < 1e-10  # Smallest eigenvalue close to 0.
True
>>>
>>> G.compute_laplacian('normalized')
>>> G.compute_fourier_basis(recompute=True)
>>> -1e-10 < G.e[0] < 1e-10 < G.e[-1] < 2  # Spectrum in [0, 2].
True

div(s)¶

Compute the divergence of a graph signal.

The divergence of a signal $s$ is defined as

\[y = D^T s,\]

where $D$ is the differential operator D.

Parameters:

s : ndarray

Signal of length G.Ne/2 living on the edges (non-directed graph).

Returns:

s_div : ndarray

Divergence signal of length G.N living on the nodes.

See also

compute_differential_operator

grad: compute the gradient

Examples

>>> G = graphs.Logo()
>>> G.N, G.Ne
(1130, 3131)
>>> s = np.random.normal(size=G.Ne)
>>> s_div = G.div(s)
>>> s_grad = G.grad(s_div)

estimate_lmax(recompute=False)[source]¶

Estimate the Laplacian’s largest eigenvalue (cached).

The result is cached and accessible by the lmax property.

Exact value given by the eigendecomposition of the Laplacian, see compute_fourier_basis(). That estimation is much faster than the eigendecomposition.

Parameters:

recompute : boolean

Force to recompute the largest eigenvalue. Default is false.

Notes

Runs the implicitly restarted Lanczos method with a large tolerance, then increases the calculated largest eigenvalue by 1 percent. For much of the PyGSP machinery, we need to approximate wavelet kernels on an interval that contains the spectrum of L. The only cost of using a larger interval is that the polynomial approximation over the larger interval may be a slightly worse approximation on the actual spectrum. As this is a very mild effect, it is not necessary to obtain very tight bounds on the spectrum of L.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> print('{:.2f}'.format(G.lmax))
13.78
>>> G = graphs.Logo()
>>> G.estimate_lmax(recompute=True)
>>> print('{:.2f}'.format(G.lmax))
13.92

extract_components()[source]¶

Split the graph into connected components.

See is_connected() for the method used to determine connectedness.

Returns:

graphs : list

A list of graph structures. Each having its own node list and weight matrix. If the graph is directed, add into the info parameter the information about the source nodes and the sink nodes.

Examples

>>> from scipy import sparse
>>> W = sparse.rand(10, 10, 0.2)
>>> W = utils.symmetrize(W)
>>> G = graphs.Graph(W=W)
>>> components = G.extract_components()
>>> has_sinks = 'sink' in components[0].info
>>> sinks_0 = components[0].info['sink'] if has_sinks else []

get_edge_list()[source]¶

Return an edge list, an alternative representation of the graph.

The weighted adjacency matrix is the canonical form used in this package to represent a graph as it is the easiest to work with when considering spectral methods.

Returns:

v_in : vector of int

v_out : vector of int

weights : vector of float

Examples

>>> G = graphs.Logo()
>>> v_in, v_out, weights = G.get_edge_list()
>>> v_in.shape, v_out.shape, weights.shape
((3131,), (3131,), (3131,))

gft(s)¶

Compute the graph Fourier transform.

The graph Fourier transform of a signal $s$ is defined as

\[\hat{s} = U^* s,\]

where $U$ is the Fourier basis attr:U and $U^*$ denotes the conjugate transpose or Hermitian transpose of $U$.

Parameters:

s : ndarray

Graph signal in the vertex domain.

Returns:

s_hat : ndarray

Representation of s in the Fourier domain.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> s = np.random.normal(size=(G.N, 5, 1))
>>> s_hat = G.gft(s)
>>> s_star = G.igft(s_hat)
>>> np.all((s - s_star) < 1e-10)
True

gft_windowed(g, f, lowmemory=True)¶

Windowed graph Fourier transform.

Parameters:

g : ndarray or Filter

Window (graph signal or kernel).

f : ndarray

Graph signal in the vertex domain.

lowmemory : bool

Use less memory (default=True).

Returns:

C : ndarray

Coefficients.

gft_windowed_gabor(s, k)¶

Gabor windowed graph Fourier transform.

Parameters:

s : ndarray

Graph signal in the vertex domain.

k : function

Gabor kernel. See pygsp.filters.Gabor.

Returns:

s : ndarray

Vertex-frequency representation of the signals.

Examples

>>> G = graphs.Logo()
>>> s = np.random.normal(size=(G.N, 2))
>>> s = G.gft_windowed_gabor(s, lambda x: x/(1.-x))
>>> s.shape
(1130, 2, 1130)

gft_windowed_normalized(g, f, lowmemory=True)¶

Normalized windowed graph Fourier transform.

Parameters:

g : ndarray

Window.

f : ndarray

Graph signal in the vertex domain.

lowmemory : bool

Use less memory. (default = True)

Returns:

C : ndarray

Coefficients.

grad(s)¶

Compute the gradient of a graph signal.

The gradient of a signal $s$ is defined as

\[y = D s,\]

where $D$ is the differential operator D.

Parameters:

s : ndarray

Signal of length G.N living on the nodes.

Returns:

s_grad : ndarray

Gradient signal of length G.Ne/2 living on the edges (non-directed graph).

See also

compute_differential_operator

div: compute the divergence

Examples

>>> G = graphs.Logo()
>>> G.N, G.Ne
(1130, 3131)
>>> s = np.random.normal(size=G.N)
>>> s_grad = G.grad(s)
>>> s_div = G.div(s_grad)
>>> np.linalg.norm(s_div - G.L.dot(s)) < 1e-10
True

igft(s_hat)¶

Compute the inverse graph Fourier transform.

The inverse graph Fourier transform of a Fourier domain signal $\hat{s}$ is defined as

\[s = U \hat{s},\]

where $U$ is the Fourier basis U.

Parameters:

s_hat : ndarray

Graph signal in the Fourier domain.

Returns:

s : ndarray

Representation of s_hat in the vertex domain.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> s_hat = np.random.normal(size=(G.N, 5, 1))
>>> s = G.igft(s_hat)
>>> s_hat_star = G.gft(s)
>>> np.all((s_hat - s_hat_star) < 1e-10)
True

is_connected(recompute=False)[source]¶

Check the strong connectivity of the graph (cached).

It uses DFS travelling on graph to ensure that each node is visited. For undirected graphs, starting at any vertex and trying to access all others is enough. For directed graphs, one needs to check that a random vertex is accessible by all others and can access all others. Thus, we can transpose the adjacency matrix and compute again with the same starting point in both phases.

Parameters:

recompute: bool

Force to recompute the connectivity if already known.

Returns:

connected : bool

True if the graph is connected.

Examples

>>> from scipy import sparse
>>> W = sparse.rand(10, 10, 0.2)
>>> G = graphs.Graph(W=W)
>>> connected = G.is_connected()

is_directed(recompute=False)[source]¶

Check if the graph has directed edges (cached).

In this framework, we consider that a graph is directed if and only if its weight matrix is non symmetric.

Parameters:

recompute : bool

Force to recompute the directedness if already known.

Returns:

directed : bool

True if the graph is directed.

Notes

Can also be used to check if a matrix is symmetrical

Examples

>>> from scipy import sparse
>>> W = sparse.rand(10, 10, 0.2)
>>> G = graphs.Graph(W=W)
>>> directed = G.is_directed()

modulate(f, k)[source]¶

Modulate the signal f to the frequency k.

Parameters:

f : ndarray

Signal (column)

k : int

Index of frequencies

Returns:

fm : ndarray

Modulated signal

plot(**kwargs)[source]¶

Plot the graph.

See pygsp.plotting.plot_graph().

plot_signal(signal, **kwargs)[source]¶

Plot a signal on that graph.

See pygsp.plotting.plot_signal().

plot_spectrogram(**kwargs)[source]¶

Plot the graph’s spectrogram.

See pygsp.plotting.plot_spectrogram().

set_coordinates(kind='spring', **kwargs)[source]¶

Set node’s coordinates (their position when plotting).

Parameters:

kind : string or array-like

Kind of coordinates to generate. It controls the position of the nodes when plotting the graph. Can either pass an array of size Nx2 or Nx3 to set the coordinates manually or the name of a layout algorithm. Available algorithms: community2D, random2D, random3D, ring2D, line1D, spring. Default is ‘spring’.

kwargs : dict

Additional parameters to be passed to the Fruchterman-Reingold force-directed algorithm when kind is spring.

Examples

>>> G = graphs.ErdosRenyi()
>>> G.set_coordinates()
>>> G.plot()

subgraph(ind)[source]¶

Create a subgraph given indices.

Parameters:

ind : list

Nodes to keep

Returns:

sub_G : Graph

Subgraph

Examples

>>> W = np.arange(16).reshape(4, 4)
>>> G = graphs.Graph(W)
>>> ind = [1, 3]
>>> sub_G = G.subgraph(ind)

translate(f, i)¶

Translate the signal f to the node i.

Parameters:

f : ndarray

Signal

i : int

Indices of vertex

Returns:

ft : translate signal

A¶

Graph adjacency matrix (the binary version of W).

The adjacency matrix defines which edges exist on the graph. It is represented as an N-by-N matrix of booleans. $A_{i,j}$ is True if $W_{i,j} > 0$.

D¶

Differential operator (for gradient and divergence).

Is computed by compute_differential_operator().

U¶

Fourier basis (eigenvectors of the Laplacian).

Is computed by compute_fourier_basis().

d¶: The degree (the number of neighbors) of each node.

dw¶: The weighted degree (the sum of weighted edges) of each node.

e¶

Eigenvalues of the Laplacian (square of graph frequencies).

Is computed by compute_fourier_basis().

lmax¶

Largest eigenvalue of the graph Laplacian.

Can be exactly computed by compute_fourier_basis() or approximated by estimate_lmax().

mu¶

Coherence of the Fourier basis.

Is computed by compute_fourier_basis().

class pygsp.graphs.Airfoil(**kwargs)[source]¶

Airfoil graph.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Airfoil()
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.5)
>>> G.plot(show_edges=True, ax=axes[1])

class pygsp.graphs.BarabasiAlbert(N=1000, m0=1, m=1, seed=None, **kwargs)[source]¶

Barabasi-Albert preferential attachment.

The Barabasi-Albert graph is constructed by connecting nodes in two steps. First, m0 nodes are created. Then, nodes are added one by one.

By lack of clarity, we take the liberty to create it as follows:

the m0 initial nodes are disconnected,

each node is connected to m of the older nodes with a probability distribution depending of the node-degrees of the other nodes, $p_n(i) = \frac{1 + k_i}{\sum_j{1 + k_j}}$.

Parameters:

N : int

Number of nodes (default is 1000)

m0 : int

Number of initial nodes (default is 1)

m : int

Number of connections at each step (default is 1) m can never be larger than m0.

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.BarabasiAlbert(N=150, seed=42)
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Comet(N=32, k=12, **kwargs)[source]¶

Comet graph.

The comet graph is a path graph with a star of degree k at its end.

Parameters:

N : int

Number of nodes.

k : int

Degree of center vertex.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Comet(15, 10)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Community(N=256, Nc=None, min_comm=None, min_deg=0, comm_sizes=None, size_ratio=1, world_density=None, comm_density=None, k_neigh=None, epsilon=None, seed=None, **kwargs)[source]¶

Community graph.

Parameters:

N : int

Number of nodes (default = 256).

Nc : int (optional)

Number of communities (default = $\lfloor \sqrt{N}/2 \rceil$).

min_comm : int (optional)

Minimum size of the communities (default = $\lfloor N/Nc/3 \rceil$).

min_deg : int (optional)

NOT IMPLEMENTED. Minimum degree of each node (default = 0).

comm_sizes : int (optional)

Size of the communities (default = random).

size_ratio : float (optional)

Ratio between the radius of world and the radius of communities (default = 1).

world_density : float (optional)

Probability of a random edge between two different communities (default = 1/N).

comm_density : float (optional)

Probability of a random edge inside any community (default = None, which implies k_neigh or epsilon will be used to determine intra-edges).

k_neigh : int (optional)

Number of intra-community connections. Not used if comm_density is defined (default = None, which implies comm_density or epsilon will be used to determine intra-edges).

epsilon : float (optional)

Largest distance at which two nodes sharing a community are connected. Not used if k_neigh or comm_density is defined (default = $\sqrt{2\sqrt{N}}/2$).

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Community(N=250, Nc=3, comm_sizes=[50, 120, 80], seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.5)
>>> G.plot(ax=axes[1])

class pygsp.graphs.DavidSensorNet(N=64, seed=None, **kwargs)[source]¶

Sensor network.

Parameters:

N : int

Number of vertices (default = 64). Values of 64 and 500 yield pre-computed and saved graphs. Other values yield randomly generated graphs.

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.DavidSensorNet()
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.ErdosRenyi(N=100, p=0.1, directed=False, self_loops=False, connected=False, max_iter=10, seed=None, **kwargs)[source]¶

Erdos Renyi graph.

The Erdos Renyi graph is constructed by randomly connecting nodes. Each edge is included in the graph with probability p, independently from any other edge. All edge weights are equal to 1.

Parameters:

N : int

Number of nodes (default is 100).

p : float

Probability to connect a node with another one.

directed : bool

Allow directed edges if True (default is False).

self_loops : bool

Allow self loops if True (default is False).

connected : bool

Force the graph to be connected (default is False).

max_iter : int

Maximum number of trials to get a connected graph (default is 10).

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.ErdosRenyi(N=64, seed=42)
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.FullConnected(N=10, **kwargs)[source]¶

Fully connected graph.

All weights are set to 1. There is no self-connections.

Parameters:

N : int

Number of vertices (default = 10)

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.FullConnected(N=20)
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=5)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Grid2d(N1=16, N2=None, **kwargs)[source]¶

2-dimensional grid graph.

Parameters:

N1 : int

Number of vertices along the first dimension.

N2 : int

Number of vertices along the second dimension (default N1).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Grid2d(N1=5, N2=4)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Logo(**kwargs)[source]¶

GSP logo.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Logo()
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.5)
>>> G.plot(ax=axes[1])

class pygsp.graphs.LowStretchTree(k=6, **kwargs)[source]¶

Low stretch tree.

Build the root of a low stretch tree on a grid of points. There are $2k$ points on each side of the grid, and therefore $2^{2k}$ vertices in total. The edge weights are all equal to 1.

Parameters:

k : int

$2^k$ points on each side of the grid of vertices.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.LowStretchTree(k=2)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Minnesota(connect=True, **kwargs)[source]¶

Minnesota road network (from MatlabBGL).

Parameters:

connect : bool

If True, the adjacency matrix is adjusted so that all edge weights are equal to 1, and the graph is connected. Set to False to get the original disconnected graph.

References

See [Gle].

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Minnesota()
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.5)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Path(N=16, **kwargs)[source]¶

Path graph.

Parameters:

N : int

Number of vertices.

References

See [Str99] for more informations.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Path(N=10)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])

class pygsp.graphs.RandomRegular(N=64, k=6, maxIter=10, seed=None, **kwargs)[source]¶

Random k-regular graph.

The random regular graph has the property that every node is connected to k other nodes. That graph is simple (without loops or double edges), k-regular (each vertex is adjacent to k nodes), and undirected.

Parameters:

N : int

Number of nodes (default is 64)

k : int

Number of connections, or degree, of each node (default is 6)

maxIter : int

Maximum number of iterations (default is 10)

seed : int

Seed for the random number generator (for reproducible graphs).

Notes

The pairing model algorithm works as follows. First create n*d half edges. Then repeat as long as possible: pick a pair of half edges and if it’s legal (doesn’t create a loop nor a double edge) add it to the graph.

References

See [KV03]. This code has been adapted from matlab to python.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.RandomRegular(N=64, k=5, seed=42)
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

is_regular()[source]¶: Troubleshoot a given regular graph.

class pygsp.graphs.RandomRing(N=64, seed=None, **kwargs)[source]¶

Ring graph with randomly sampled nodes.

Parameters:

N : int

Number of vertices.

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.RandomRing(N=10, seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])
>>> _ = axes[1].set_xlim(-1.1, 1.1)
>>> _ = axes[1].set_ylim(-1.1, 1.1)

class pygsp.graphs.Ring(N=64, k=1, **kwargs)[source]¶

K-regular ring graph.

Parameters:

N : int

Number of vertices.

k : int

Number of neighbors in each direction.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=10)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Sensor(N=64, Nc=2, regular=False, n_try=50, distribute=False, connected=True, seed=None, **kwargs)[source]¶

Random sensor graph.

Parameters:

N : int

Number of nodes (default = 64)

Nc : int

Minimum number of connections (default = 2)

regular : bool

Flag to fix the number of connections to nc (default = False)

n_try : int

Number of attempt to create the graph (default = 50)

distribute : bool

To distribute the points more evenly (default = False)

connected : bool

To force the graph to be connected (default = True)

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Sensor(N=64, seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.StochasticBlockModel(N=1024, k=5, z=None, M=None, p=0.7, q=None, directed=False, self_loops=False, connected=False, max_iter=10, seed=None, **kwargs)[source]¶

Stochastic Block Model (SBM).

The Stochastic Block Model graph is constructed by connecting nodes with a probability which depends on the cluster of the two nodes. One can define the clustering association of each node, denoted by vector z, but also the probability matrix M. All edge weights are equal to 1. By default, Mii > Mjk and nodes are uniformly clusterized.

Parameters:

N : int

Number of nodes (default is 1024).

k : float

Number of classes (default is 5).

z : array_like

the vector of length N containing the association between nodes and classes (default is random uniform).

M : array_like

the k by k matrix containing the probability of connecting nodes based on their class belonging (default using p and q).

p : float or array_like

the diagonal value(s) for the matrix M. If scalar they all have the same value. Otherwise expect a length k vector (default is p = 0.7).

q : float or array_like

the off-diagonal value(s) for the matrix M. If scalar they all have the same value. Otherwise expect a k x k matrix, diagonal will be discarded (default is q = 0.3/k).

directed : bool

Allow directed edges if True (default is False).

self_loops : bool

Allow self loops if True (default is False).

connected : bool

Force the graph to be connected (default is False).

max_iter : int

Maximum number of trials to get a connected graph (default is 10).

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.StochasticBlockModel(
...     100, k=3, p=[0.4, 0.6, 0.3], q=0.02, seed=42)
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.8)
>>> G.plot(ax=axes[1])

class pygsp.graphs.SwissRoll(N=400, a=1, b=4, dim=3, thresh=1e-06, s=None, noise=False, srtype='uniform', seed=None, **kwargs)[source]¶

Sampled Swiss roll manifold.

Parameters:

N : int

Number of vertices (default = 400)

a : int

(default = 1)

b : int

(default = 4)

dim : int

(default = 3)

thresh : float

(default = 1e-6)

s : float

sigma (default = sqrt(2./N))

noise : bool

Wether to add noise or not (default = False)

srtype : str

Swiss roll Type, possible arguments are ‘uniform’ or ‘classic’ (default = ‘uniform’)

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.SwissRoll(N=200, seed=42)
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)
>>> ax2 = fig.add_subplot(122, projection='3d')
>>> _ = ax1.spy(G.W, markersize=1)
>>> G.plot(ax=ax2)

class pygsp.graphs.Torus(Nv=16, Mv=None, **kwargs)[source]¶

Sampled torus manifold.

Parameters:

Nv : int

Number of vertices along the first dimension (default is 16)

Mv : int

Number of vertices along the second dimension (default is Nv)

References

See [Str99] for more informations.

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Torus(10)
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)
>>> ax2 = fig.add_subplot(122, projection='3d')
>>> _ = ax1.spy(G.W, markersize=1.5)
>>> G.plot(ax=ax2)
>>> _ = ax2.set_zlim(-1.5, 1.5)

class pygsp.graphs.NNGraph(Xin, NNtype='knn', use_flann=False, center=True, rescale=True, k=10, sigma=0.1, epsilon=0.01, gtype=None, plotting={}, symmetrize_type='average', dist_type='euclidean', order=0, **kwargs)[source]¶

Nearest-neighbor graph from given point cloud.

Parameters:

Xin : ndarray

Input points, Should be an N-by-d matrix, where N is the number of nodes in the graph and d is the dimension of the feature space.

NNtype : string, optional

Type of nearest neighbor graph to create. The options are ‘knn’ for k-Nearest Neighbors or ‘radius’ for epsilon-Nearest Neighbors (default is ‘knn’).

use_flann : bool, optional

Use Fast Library for Approximate Nearest Neighbors (FLANN) or not. (default is False)

center : bool, optional

Center the data so that it has zero mean (default is True)

rescale : bool, optional

Rescale the data so that it lies in a l2-sphere (default is True)

k : int, optional

Number of neighbors for knn (default is 10)

sigma : float, optional

Width parameter of the similarity kernel (default is 0.1)

epsilon : float, optional

Radius for the epsilon-neighborhood search (default is 0.01)

gtype : string, optional

The type of graph (default is ‘nearest neighbors’)

plotting : dict, optional

Dictionary of plotting parameters. See pygsp.plotting. (default is {})

symmetrize_type : string, optional

Type of symmetrization to use for the adjacency matrix. See pygsp.utils.symmetrization() for the options. (default is ‘average’)

dist_type : string, optional

Type of distance to compute. See pyflann.index.set_distance_type() for possible options. (default is ‘euclidean’)

order : float, optional

Only used if dist_type is ‘minkowski’; represents the order of the Minkowski distance. (default is 0)

Examples

>>> import matplotlib.pyplot as plt
>>> X = np.random.RandomState(42).uniform(size=(30, 2))
>>> G = graphs.NNGraph(X)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=5)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Bunny(**kwargs)[source]¶

Stanford bunny (NN-graph).

References

See [TL94].

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Bunny()
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)
>>> ax2 = fig.add_subplot(122, projection='3d')
>>> _ = ax1.spy(G.W, markersize=0.1)
>>> G.plot(ax=ax2)

class pygsp.graphs.Cube(radius=1, nb_pts=300, nb_dim=3, sampling='random', seed=None, **kwargs)[source]¶

Hyper-cube (NN-graph).

Parameters:

radius : float

Edge lenght (default = 1)

nb_pts : int

Number of vertices (default = 300)

nb_dim : int

Dimension (default = 3)

sampling : string

Variance of the distance kernel (default = ‘random’) (Can now only be ‘random’)

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Cube(seed=42)
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)
>>> ax2 = fig.add_subplot(122, projection='3d')
>>> _ = ax1.spy(G.W, markersize=0.5)
>>> G.plot(ax=ax2)

class pygsp.graphs.ImgPatches(img, patch_shape=(3, 3), **kwargs)[source]¶

NN-graph between patches of an image.

Extract a feature vector in the form of a patch for every pixel of an image, then construct a nearest-neighbor graph between these feature vectors. The feature matrix, i.e. the patches, can be found in Xin.

Parameters:

img : array

Input image.

patch_shape : tuple, optional

Dimensions of the patch window. Syntax: (height, width), or (height,), in which case width = height.

kwargs : dict

Parameters passed to NNGraph.

Notes

The feature vector of a pixel i will consist of the stacking of the intensity values of all pixels in the patch centered at i, for all color channels. So, if the input image has d color channels, the dimension of the feature vector of each pixel is (patch_shape[0] * patch_shape[1] * d).

Examples

>>> import matplotlib.pyplot as plt
>>> from skimage import data, img_as_float
>>> img = img_as_float(data.camera()[::64, ::64])
>>> G = graphs.ImgPatches(img, patch_shape=(3, 3))
>>> print('{} nodes ({} x {} pixels)'.format(G.Xin.shape[0], *img.shape))
64 nodes (8 x 8 pixels)
>>> print('{} features per node'.format(G.Xin.shape[1]))
9 features per node
>>> G.set_coordinates(kind='spring', seed=42)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Grid2dImgPatches(img, aggregate=<function Grid2dImgPatches.<lambda>>, **kwargs)[source]¶

Union of a patch graph with a 2D grid graph.

Parameters:

img : array

Input image.

aggregate: callable, optional

Function to aggregate the weights Wp of the patch graph and the Wg of the grid graph. Default is lambda Wp, Wg: Wp + Wg.

kwargs : dict

Parameters passed to ImgPatches.

Examples

>>> import matplotlib.pyplot as plt
>>> from skimage import data, img_as_float
>>> img = img_as_float(data.camera()[::64, ::64])
>>> G = graphs.Grid2dImgPatches(img)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=2)
>>> G.plot(ax=axes[1])

class pygsp.graphs.Sphere(radius=1, nb_pts=300, nb_dim=3, sampling='random', seed=None, **kwargs)[source]¶

Spherical-shaped graph (NN-graph).

Parameters:

radius : flaot

Radius of the sphere (default = 1)

nb_pts : int

Number of vertices (default = 300)

nb_dim : int

Dimension (default = 3)

sampling : sting

Variance of the distance kernel (default = ‘random’) (Can now only be ‘random’)

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.Sphere(nb_pts=100, seed=42)
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)
>>> ax2 = fig.add_subplot(122, projection='3d')
>>> _ = ax1.spy(G.W, markersize=1.5)
>>> G.plot(ax=ax2)

class pygsp.graphs.TwoMoons(moontype='standard', dim=2, sigmag=0.05, N=400, sigmad=0.07, d=0.5, seed=None, **kwargs)[source]¶

Two Moons (NN-graph).

Parameters:

moontype : ‘standard’ or ‘synthesized’

You have the freedom to chose if you want to create a standard two_moons graph or a synthesized one (default is ‘standard’). ‘standard’ : Create a two_moons graph from a based graph. ‘synthesized’ : Create a synthesized two_moon

sigmag : float

Variance of the distance kernel (default = 0.05)

dim : int

The dimensionality of the points (default = 2). Only valid for moontype == ‘standard’.

N : int

Number of vertices (default = 2000) Only valid for moontype == ‘synthesized’.

sigmad : float

Variance of the data (do not set it too high or you won’t see anything) (default = 0.05) Only valid for moontype == ‘synthesized’.

d : float

Distance of the two moons (default = 0.5) Only valid for moontype == ‘synthesized’.

seed : int

Seed for the random number generator (for reproducible graphs).

Examples

>>> import matplotlib.pyplot as plt
>>> G = graphs.TwoMoons()
>>> fig, axes = plt.subplots(1, 2)
>>> _ = axes[0].spy(G.W, markersize=0.5)
>>> G.plot(show_edges=True, ax=axes[1])

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`(x)	Evaluate the kernels at given frequencies.
`Filter.filter`(s[, method, order])	Filter signals (analysis or synthesis).
`Filter.analyze`(s[, method, order])	Convenience alias to `filter()`.
`Filter.synthesize`(s[, method, order])	Convenience wrapper around `filter()`.
`Filter.compute_frame`(**kwargs)	Compute the associated frame.
`Filter.estimate_frame_bounds`([min, max, N, …])	Estimate lower and upper frame bounds.
`Filter.plot`(**kwargs)	Plot the filter bank’s frequency response.
`Filter.localize`(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.

class pygsp.filters.Filter(G, kernels)[source]¶

The base Filter class.

Provide a common interface (and implementation) to filter objects.
Can be instantiated to construct custom filters from functions.
Initialize attributes for derived classes.

Parameters:

G : graph

The graph to which the filter bank is tailored.

kernels : function or list of functions

A (list of) function(s) defining the filter bank in the Fourier domain. One function per filter.

Examples

>>>
>>> G = graphs.Logo()
>>> my_filter = filters.Filter(G, lambda x: x / (1. + x))
>>>
>>> # Signal: Kronecker delta.
>>> signal = np.zeros(G.N)
>>> signal[42] = 1
>>>
>>> filtered_signal = my_filter.filter(signal)

Attributes

G	(Graph) The graph to which the filter bank was tailored. It is a reference to the graph passed when instantiating the class.
kernels	(function or list of functions) A (list of) function defining the filter bank. One function per filter. Either passed by the user when instantiating the base class, either constructed by the derived classes.
Nf	(int) Number of filters in the filter bank.

analyze(s, method='chebyshev', order=30)[source]¶: Convenience alias to filter().

can_dual()[source]¶: Creates a dual graph form a given graph

compute_frame(**kwargs)[source]¶

Compute the associated frame.

The size of the returned matrix operator $D$ is N x MN, where M is the number of filters and N the number of nodes. Multiplying this matrix with a set of signals is equivalent to analyzing them with the associated filterbank. Though computing this matrix is a rather inefficient way of doing it.

The frame is defined as follows:

\[g_i(L) = U g_i(\Lambda) U^*,\]

where $g$ is the filter kernel, $L$ is the graph Laplacian, $\Lambda$ is a diagonal matrix of the Laplacian’s eigenvalues, and $U$ is the Fourier basis, i.e. its columns are the eigenvectors of the Laplacian.

Parameters:

kwargs: dict

Parameters to be passed to the analyze() method.

Returns:

frame : ndarray

Matrix of size N x MN.

See also

filter: more efficient way to filter signals

Examples

Filtering signals as a matrix multiplication.

>>> G = graphs.Sensor(N=1000, seed=42)
>>> G.estimate_lmax()
>>> f = filters.MexicanHat(G, Nf=6)
>>> s = np.random.uniform(size=G.N)
>>>
>>> frame = f.compute_frame()
>>> frame.shape
(1000, 1000, 6)
>>> frame = frame.reshape(G.N, -1).T
>>> s1 = np.dot(frame, s)
>>> s1 = s1.reshape(G.N, -1)
>>>
>>> s2 = f.filter(s)
>>> np.all((s1 - s2) < 1e-10)
True

estimate_frame_bounds(min=0, max=None, N=1000, use_eigenvalues=False)[source]¶

Estimate lower and upper frame bounds.

The frame bounds are estimated using the vector np.linspace(min, max, N) with min=0 and max=G.lmax by default. The eigenvalues G.e can be used instead if you set use_eigenvalues to True.

Parameters:

min : float

The lowest value the filter bank is evaluated at. By default filtering is bounded by the eigenvalues of G, i.e. min = 0.

max : float

The largest value the filter bank is evaluated at. By default filtering is bounded by the eigenvalues of G, i.e. max = G.lmax.

N : int

Number of points where the filter bank is evaluated. Default is 1000.

use_eigenvalues : bool

Set to True to use the Laplacian eigenvalues instead.

Returns:

A : float

Lower frame bound of the filter bank.

B : float

Upper frame bound of the filter bank.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> f = filters.MexicanHat(G)

Bad estimation:

>>> A, B = f.estimate_frame_bounds(min=1, max=20, N=100)
>>> print('A={:.3f}, B={:.3f}'.format(A, B))
A=0.125, B=0.270

Better estimation:

>>> A, B = f.estimate_frame_bounds()
>>> print('A={:.3f}, B={:.3f}'.format(A, B))
A=0.177, B=0.270

Best estimation:

>>> G.compute_fourier_basis()
>>> A, B = f.estimate_frame_bounds(use_eigenvalues=True)
>>> print('A={:.3f}, B={:.3f}'.format(A, B))
A=0.177, B=0.270

The Itersine filter bank defines a tight frame:

>>> f = filters.Itersine(G)
>>> A, B = f.estimate_frame_bounds(use_eigenvalues=True)
>>> print('A={:.3f}, B={:.3f}'.format(A, B))
A=1.000, B=1.000

evaluate(x)[source]¶

Evaluate the kernels at given frequencies.

Parameters:

x : ndarray

Graph frequencies at which to evaluate the filter.

Returns:

y : ndarray

Frequency response of the filters. Shape (G.Nf, len(x)).

Examples

Frequency response of a low-pass filter:

>>> import matplotlib.pyplot as plt
>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> f = filters.Expwin(G)
>>> G.compute_fourier_basis()
>>> y = f.evaluate(G.e)
>>> plt.plot(G.e, y[0])  
[<matplotlib.lines.Line2D object at ...>]

filter(s, method='chebyshev', order=30)[source]¶

Filter signals (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 $s$, where the atoms of the transform dictionary are generalized translations of each graph spectral filter to each vertex on the graph:

\[c = D^* s,\]

where the columns of $D$ are $g_{i,m} = T_i g_m$ and $T_i$ is a generalized translation operator applied to each filter $\hat{g}_m(\cdot)$. Each column of $c$ is the response of the signal to one filter.

In other words, this function is applying the analysis operator $D^*$, respectively the synthesis operator $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 [HVG11] 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)  
<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,)

Filter and reconstruct our signal:

>>> g = filters.MexicanHat(G, Nf=4)
>>> s2 = g.analyze(s1)
>>> s2.shape
(30, 4)
>>> s2 = g.synthesize(s2)
>>> s2.shape
(30,)

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.analyze(s1, method='exact')
>>> s2 = g.synthesize(s2, method='exact')
>>> np.linalg.norm(s1 - s2) < 1e-10
True

localize(i, **kwargs)[source]¶

Localize the kernels at a node (to visualize them).

That is particularly useful to visualize a filter in the vertex domain.

A kernel is localized by filtering a Kronecker delta, i.e.

\[\begin{split}g(L) s = g(L)_i, \text{ where } s_j = \delta_{ij} = \begin{cases} 0 \text{ if } i \neq j \\ 1 \text{ if } i = j \end{cases}\end{split}\]

Parameters:

i : int

Index of the node where to localize the kernel.

kwargs: dict

Parameters to be passed to the analyze() method.

Returns:

s : ndarray

Kernel localized at vertex i.

Examples

Visualize heat diffusion on a grid by localizing the heat kernel.

>>> import matplotlib
>>> N = 20
>>> DELTA = N//2 * (N+1)
>>> G = graphs.Grid2d(N)
>>> G.estimate_lmax()
>>> g = filters.Heat(G, 100)
>>> s = g.localize(DELTA)
>>> G.plot_signal(s, highlight=DELTA)

plot(**kwargs)[source]¶

Plot the filter bank’s frequency response.

See pygsp.plotting.plot_filter().

synthesize(s, method='chebyshev', order=30)[source]¶

Convenience wrapper around filter().

Will be an alias to adjoint().filter() in the future.

class pygsp.filters.Abspline(G, Nf=6, lpfactor=20, scales=None, **kwargs)[source]¶

Design an A B cubic spline wavelet filter bank.

Parameters:

G : graph

Nf : int

Number of filters from 0 to lmax (default = 6)

lpfactor : int

Low-pass factor lmin=lmax/lpfactor will be used to determine scales, the scaling function will be created to fill the lowpass gap. (default = 20)

scales : ndarray

Vector of scales to be used. By default, initialized with pygsp.utils.compute_log_scales().

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Abspline(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Expwin(G, bmax=0.2, a=1.0, **kwargs)[source]¶

Design an exponential window filter.

Parameters:

G : graph

bmax : float

Maximum relative band (default = 0.2)

a : int

Slope parameter (default = 1)

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Expwin(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Gabor(G, k)[source]¶

Design a Gabor filter bank.

Design a filter bank where the kernel k is placed at each graph frequency.

Parameters:

G : graph

k : lambda function

kernel

Examples

>>> G = graphs.Logo()
>>> k = lambda x: x / (1. - x)
>>> g = filters.Gabor(G, k);

class pygsp.filters.HalfCosine(G, Nf=6, **kwargs)[source]¶

Design an half cosine filter bank (tight frame).

Parameters:

G : graph

Nf : int

Number of filters from 0 to lmax (default = 6)

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.HalfCosine(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Heat(G, tau=10, normalize=False, **kwargs)[source]¶

Design a filter bank of heat kernels.

The filter is defined in the spectral domain as

\[\hat{g}(x) = \exp \left( \frac{-\tau x}{\lambda_{\text{max}}} \right),\]

and is as such a low-pass filter. An application of this filter to a signal simulates heat diffusion.

Parameters:

G : graph

tau : int or list of ints

Scaling parameter. If a list, creates a filter bank with one filter per value of tau.

normalize : bool

Normalizes the kernel. Needs the eigenvalues.

Examples

Regular heat kernel.

>>> G = graphs.Logo()
>>> g = filters.Heat(G, tau=[5, 10])
>>> print('{} filters'.format(g.Nf))
2 filters
>>> y = g.evaluate(G.e)
>>> print('{:.2f}'.format(np.linalg.norm(y[0])))
9.76

Normalized heat kernel.

>>> g = filters.Heat(G, tau=[5, 10], normalize=True)
>>> y = g.evaluate(G.e)
>>> print('{:.2f}'.format(np.linalg.norm(y[0])))
1.00

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Heat(G, tau=[5, 10, 100])
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Held(G, a=0.6666666666666666, **kwargs)[source]¶

Design 2 filters with the Held construction (tight frame).

This function create a parseval filterbank of $2$ filters. The low-pass filter is defined by the function

\[\begin{split}f_{l}=\begin{cases} 1 & \mbox{if }x\leq a\\ \sin\left(2\pi\mu\left(\frac{x}{8a}\right)\right) & \mbox{if }a<x\leq2a\\ 0 & \mbox{if }x>2a \end{cases}\end{split}\]

with

\[\mu(x) = -1+24x-144*x^2+256*x^3\]

The high pass filter is adaptated to obtain a tight frame.

Parameters:

G : graph

a : float

See equation above for this parameter The spectrum is scaled between 0 and 2 (default = 2/3)

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Held(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Itersine(G, Nf=6, overlap=2.0, **kwargs)[source]¶

Design an itersine filter bank (tight frame).

Create an itersine half overlap filter bank of Nf filters. Going from 0 to lambda_max.

Parameters:

G : graph

Nf : int (optional)

Number of filters from 0 to lmax. (default = 6)

overlap : int (optional)

(default = 2)

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.HalfCosine(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.MexicanHat(G, Nf=6, lpfactor=20, scales=None, normalize=False, **kwargs)[source]¶

Design a filter bank of Mexican hat wavelets.

The Mexican hat wavelet is the second oder derivative of a Gaussian. Since we express the filter in the Fourier domain, we find:

\[\hat{g}_b(x) = x * \exp(-x)\]

for the band-pass filter. Note that in our convention the eigenvalues of the Laplacian are equivalent to the square of graph frequencies, i.e. $x = \lambda^2$.

The low-pass filter is given by

Parameters:

G : graph

Nf : int

Number of filters to cover the interval [0, lmax].

lpfactor : int

Low-pass factor. lmin=lmax/lpfactor will be used to determine scales. The scaling function will be created to fill the low-pass gap.

scales : array-like

Scales to be used. By default, initialized with pygsp.utils.compute_log_scales().

normalize : bool

Whether to normalize the wavelet by the factor sqrt(scales).

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.MexicanHat(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Meyer(G, Nf=6, scales=None, **kwargs)[source]¶

Design a filter bank of Meyer wavelets (tight frame).

Parameters:

G : graph

Nf : int

Number of filters from 0 to lmax (default = 6).

scales : ndarray

Vector of scales to be used (default: log scale).

References

Use of this kernel for SGWT proposed by Nora Leonardi and Dimitri Van De Ville in [LVDV11].

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Meyer(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Papadakis(G, a=0.75, **kwargs)[source]¶

Design 2 filters with the Papadakis construction (tight frame).

This function creates a Parseval filter bank of 2 filters. The low-pass filter is defined by the function

\[\begin{split}f_{l}=\begin{cases} 1 & \mbox{if }x\leq a\\ \sqrt{1-\frac{\sin\left(\frac{3\pi}{2a}x\right)}{2}} & \mbox{if }a<x\leq \frac{5a}{3} \\ 0 & \mbox{if }x>\frac{5a}{3} \end{cases}\end{split}\]

The high pass filter is adapted to obtain a tight frame.

Parameters:

G : graph

a : float

See above equation for this parameter. The spectrum is scaled between 0 and 2 (default = 3/4).

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Papadakis(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Regular(G, d=3, **kwargs)[source]¶

Design 2 filters with the regular construction (tight frame).

This function creates a Parseval filter bank of 2 filters. The low-pass filter is defined by a function $f_l(x)$ between $0$ and $2$. For $d = 0$.

\[f_{l}= \sin\left( \frac{\pi}{4} x \right)\]

For $d = 1$

\[f_{l}= \sin\left( \frac{\pi}{4} \left( 1+ \sin\left(\frac{\pi}{2}(x-1)\right) \right) \right)\]

For $d = 2$

\[f_{l}= \sin\left( \frac{\pi}{4} \left( 1+ \sin\left(\frac{\pi}{2} \sin\left(\frac{\pi}{2}(x-1)\right)\right) \right) \right)\]

And so forth for other degrees $d$.

The high pass filter is adapted to obtain a tight frame.

Parameters:

G : graph

d : float

Degree (default = 3). See above equations.

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Regular(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.Simoncelli(G, a=0.6666666666666666, **kwargs)[source]¶

Design 2 filters with the Simoncelli construction (tight frame).

This function creates a Parseval filter bank of 2 filters. The low-pass filter is defined by the function

\[\begin{split}f_{l}=\begin{cases} 1 & \mbox{if }x\leq a\\ \cos\left(\frac{\pi}{2}\frac{\log\left(\frac{x}{2}\right)}{\log(2)}\right) & \mbox{if }a<x\leq2a\\ 0 & \mbox{if }x>2a \end{cases}\end{split}\]

The high pass filter is adapted to obtain a tight frame.

Parameters:

G : graph

a : float

See above equation for this parameter. The spectrum is scaled between 0 and 2 (default = 2/3).

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.Simoncelli(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

class pygsp.filters.SimpleTight(G, Nf=6, scales=None, **kwargs)[source]¶

Design a simple tight frame filter bank (tight frame).

These filters have been designed to be a simple tight frame wavelet filter bank. The kernel is similar to Meyer, but simpler. The function is essentially $\sin^2(x)$ in ascending part and $\cos^2$ in descending part.

Parameters:

G : graph

Nf : int

Number of filters to cover the interval [0, lmax].

scales : array-like

Scales to be used. Defaults to log scale.

Examples

Filter bank’s representation in Fourier and time (ring graph) domains.

>>> import matplotlib.pyplot as plt
>>> G = graphs.Ring(N=20)
>>> G.estimate_lmax()
>>> G.set_coordinates('line1D')
>>> g = filters.SimpleTight(G)
>>> s = g.localize(G.N // 2)
>>> fig, axes = plt.subplots(1, 2)
>>> g.plot(ax=axes[0])
>>> G.plot_signal(s, ax=axes[1])

pygsp.filters.compute_cheby_coeff(f, m=30, N=None, *args, **kwargs)[source]¶

Compute Chebyshev coefficients for a Filterbank.

Parameters:

f : Filter

Filterbank with at least 1 filter

m : int

Maximum order of Chebyshev coeff to compute (default = 30)

N : int

Grid order used to compute quadrature (default = m + 1)

i : int

Index of the Filterbank element to compute (default = 0)

Returns:

c : ndarray

Matrix of Chebyshev coefficients

pygsp.filters.compute_jackson_cheby_coeff(filter_bounds, delta_lambda, m)[source]¶

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.

Parameters:

filter_bounds : list

[a, b]

delta_lambda : list

[lambda_min, lambda_max]

m : int

Returns:

ch : ndarray

jch : ndarray

References

[TPGV16]

pygsp.filters.cheby_op(G, c, signal, **kwargs)[source]¶

Chebyshev polynomial of graph Laplacian applied to vector.

Parameters:

G : Graph

c : ndarray or list of ndarrays

Chebyshev coefficients for a Filter or a Filterbank

signal : ndarray

Signal to filter

Returns:

r : ndarray

Result of the filtering

pygsp.filters.cheby_rect(G, bounds, signal, **kwargs)[source]¶

Fast filtering using Chebyshev polynomial for a perfect rectangle filter.

Parameters:

G : Graph

bounds : array-like

The bounds of the pass-band filter

signal : array-like

Signal to filter

order : int (optional)

Order of the Chebyshev polynomial (default: 30)

Returns:

r : array-like

Result of the filtering

pygsp.filters.lanczos(A, order, x)[source]¶

TODO short description

Parameters:	A: ndarray

pygsp.filters.lanczos_op(f, s, order=30)[source]¶

Perform the lanczos approximation of the signal s.

Parameters:

f: Filter

s : ndarray

Signal to approximate.

order : int

Degree of the lanczos approximation. (default = 30)

Returns:

L : ndarray

lanczos approximation of 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.

pygsp.plotting.close(*args, **kwargs)[source]¶

Close created figures.

Alias to plt.close().

pygsp.plotting.close_all()[source]¶: Close all opened windows.

pygsp.plotting.plot(O, **kwargs)[source]¶

Main plotting function.

This convenience function either calls plot_graph() or plot_filter() given the type of the passed object. Parameters can be passed to those functions.

Parameters:

O : Graph, Filter

object to plot

Examples

>>> from pygsp import plotting
>>> G = graphs.Logo()
>>> plotting.plot(G)

pygsp.plotting.plot_filter(obj, *args, **kwargs)[source]¶

pygsp.plotting.plot_graph(G, backend=None, **kwargs)[source]¶

Plot a graph or a list of graphs.

Parameters:

G : Graph

Graph to plot.

show_edges : bool

True to draw edges, false to only draw vertices. Default True if less than 10,000 edges to draw. Note that drawing a large number of edges might be particularly slow.

backend: {‘matplotlib’, ‘pyqtgraph’}

Defines the drawing backend to use. Defaults to BACKEND.

vertex_size : float

Size of circle representing each node.

plot_name : str

name of the plot

save_as : str

Whether to save the plot as save_as.png and save_as.pdf. Shown in a window if None (default). Only available with the matplotlib backend.

ax : matplotlib.axes

Axes where to draw the graph. Optional, created if not passed. Only available with the matplotlib backend.

Examples

>>> from pygsp import plotting
>>> G = graphs.Logo()
>>> plotting.plot_graph(G)

pygsp.plotting.plot_signal(G, signal, backend=None, **kwargs)[source]¶

Plot a signal on top of a graph.

Parameters:

G : Graph

Graph to plot a signal on top.

signal : array of int

Signal to plot. Signal length should be equal to the number of nodes.

show_edges : bool

True to draw edges, false to only draw vertices. Default True if less than 10,000 edges to draw. Note that drawing a large number of edges might be particularly slow.

cp : list of int

NOT IMPLEMENTED. Camera position when plotting a 3D graph.

vertex_size : float

Size of circle representing each node.

highlight : iterable

List of indices of vertices to be highlighted. Useful to e.g. show where a filter was localized. Only available with the matplotlib backend.

colorbar : bool

Whether to plot a colorbar indicating the signal’s amplitude. Only available with the matplotlib backend.

limits : [vmin, vmax]

Maps colors from vmin to vmax. Defaults to signal minimum and maximum value. Only available with the matplotlib backend.

bar : boolean

NOT IMPLEMENTED. Signal values are displayed using colors when False, and bars when True (default False).

bar_width : int

NOT IMPLEMENTED. Width of the bar (default 1).

backend: {‘matplotlib’, ‘pyqtgraph’}

Defines the drawing backend to use. Defaults to BACKEND.

plot_name : string

Name of the plot.

save_as : str

Whether to save the plot as save_as.png and save_as.pdf. Shown in a window if None (default). Only available with the matplotlib backend.

ax : matplotlib.axes

Axes where to draw the graph. Optional, created if not passed. Only available with the matplotlib backend.

Examples

>>> from pygsp import plotting
>>> G = graphs.Grid2d(4)
>>> signal = np.sin((np.arange(16) * 2*np.pi/16))
>>> plotting.plot_signal(G, signal)

pygsp.plotting.plot_spectrogram(G, node_idx=None)[source]¶

Plot the spectrogram of the given graph.

Parameters:

G : Graph

Graph to analyse.

node_idx : ndarray

Order to sort the nodes in the spectrogram

Examples

>>> from pygsp import plotting
>>> G = graphs.Ring(15)
>>> plotting.plot_spectrogram(G)

pygsp.plotting.show(*args, **kwargs)[source]¶

Show created figures.

Alias to plt.show(). By default, showing plots does not block the prompt.

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).

pygsp.reduction.graph_multiresolution(G, levels, sparsify=True, sparsify_eps=None, downsampling_method='largest_eigenvector', reduction_method='kron', compute_full_eigen=False, reg_eps=0.005)[source]¶

Compute a pyramid of graphs (by Kron reduction).

‘graph_multiresolution(G,levels)’ computes a multiresolution of graph by repeatedly downsampling and performing graph reduction. The default downsampling method is the largest eigenvector method based on the polarity of the components of the eigenvector associated with the largest graph Laplacian eigenvalue. The default graph reduction method is Kron reduction followed by a graph sparsification step. param is a structure of optional parameters.

Parameters:

G : Graph structure

The graph to reduce.

levels : int

Number of level of decomposition

lambd : float

Stability parameter. It adds self loop to the graph to give the algorithm some stability (default = 0.025). [UNUSED?!]

sparsify : bool

To perform a spectral sparsification step immediately after the graph reduction (default is True).

sparsify_eps : float

Parameter epsilon used in the spectral sparsification (default is min(10/sqrt(G.N),.3)).

downsampling_method: string

The graph downsampling method (default is ‘largest_eigenvector’).

reduction_method : string

The graph reduction method (default is ‘kron’)

compute_full_eigen : bool

To also compute the graph Laplacian eigenvalues and eigenvectors for every graph in the multiresolution sequence (default is False).

reg_eps : float

The regularized graph Laplacian is $\bar{L}=L+\epsilon I$. A smaller epsilon may lead to better regularization, but will also require a higher order Chebyshev approximation. (default is 0.005)

Returns:

Gs : list

A list of graph layers.

Examples

>>> from pygsp import reduction
>>> levels = 5
>>> G = graphs.Sensor(N=512)
>>> G.compute_fourier_basis()
>>> Gs = reduction.graph_multiresolution(G, levels, sparsify=False)
>>> for idx in range(levels):
...     Gs[idx].plotting['plot_name'] = 'Reduction level: {}'.format(idx)
...     Gs[idx].plot()

pygsp.reduction.graph_sparsify(M, epsilon, maxiter=10)[source]¶

Sparsify a graph (with Spielman-Srivastava).

Parameters:

M : Graph or sparse matrix

Graph structure or a Laplacian matrix

epsilon : int

Sparsification parameter

Returns:

Mnew : Graph or sparse matrix

New graph structure or sparse matrix

Notes

Epsilon should be between 1/sqrt(N) and 1

References

See [SS11], [Rud99] and [RV07]. for more informations

Examples

>>> from pygsp import reduction
>>> G = graphs.Sensor(256, Nc=20, distribute=True)
>>> epsilon = 0.4
>>> G2 = reduction.graph_sparsify(G, epsilon)

pygsp.reduction.interpolate(G, f_subsampled, keep_inds, order=100, reg_eps=0.005, **kwargs)[source]¶

Interpolate a graph signal.

Parameters:

G : Graph

f_subsampled : ndarray

A graph signal on the graph G.

keep_inds : ndarray

List of indices on which the signal is sampled.

order : int

Degree of the Chebyshev approximation (default = 100).

reg_eps : float

The regularized graph Laplacian is $bar{L}=L+epsilon I$. A smaller epsilon may lead to better regularization, but will also require a higher order Chebyshev approximation.

Returns:

f_interpolated : ndarray

Interpolated graph signal on the full vertex set of G.

References

See [Pes09]

pygsp.reduction.kron_reduction(G, ind)[source]¶

Compute the Kron reduction.

This function perform the Kron reduction of the weight matrix in the graph G, with boundary nodes labeled by ind. This function will create a new graph with a weight matrix Wnew that contain only boundary nodes and is computed as the Schur complement of the original matrix with respect to the selected indices.

Parameters:

G : Graph or sparse matrix

Graph structure or weight matrix

ind : list

indices of the nodes to keep

Returns:

Gnew : Graph or sparse matrix

New graph structure or weight matrix

References

See [DB13]

pygsp.reduction.pyramid_analysis(Gs, f, **kwargs)[source]¶

Compute the graph pyramid transform coefficients.

Parameters:

Gs : list of graphs

A multiresolution sequence of graph structures.

f : ndarray

Graph signal to analyze.

h_filters : list

A list of filter that will be used for the analysis and sythesis operator. If only one filter is given, it will be used for all levels. Default is h(x) = 1 / (2x+1)

Returns:

ca : ndarray

Coarse approximation at each level

pe : ndarray

Prediction error at each level

h_filters : list

Graph spectral filters applied

References

See [SFV13] and [Pes09].

pygsp.reduction.pyramid_synthesis(Gs, cap, pe, order=30, **kwargs)[source]¶

Synthesize a signal from its pyramid coefficients.

Parameters:

Gs : Array of Graphs

A multiresolution sequence of graph structures.

cap : ndarray

Coarsest approximation of the original signal.

pe : ndarray

Prediction error at each level.

use_exact : bool

To use exact graph spectral filtering instead of the Chebyshev approximation.

order : int

Degree of the Chebyshev approximation (default=30).

least_squares : bool

To use the least squares synthesis (default=False).

h_filters : ndarray

The filters used in the analysis operator. These are required for least squares synthesis, but not for the direct synthesis method.

use_landweber : bool

To use the Landweber iteration approximation in the least squares synthesis.

reg_eps : float

Interpolation parameter.

landweber_its : int

Number of iterations in the Landweber approximation for least squares synthesis.

landweber_tau : float

Parameter for the Landweber iteration.

Returns:

reconstruction : ndarray

The reconstructed signal.

ca : ndarray

Coarse approximations at each level

pygsp.reduction.tree_multiresolution(G, Nlevel, reduction_method='resistance_distance', compute_full_eigen=False, root=None)[source]¶

Compute a multiresolution of trees

Parameters:

G : Graph

Graph structure of a tree.

Nlevel : Number of times to downsample and coarsen the tree

root : int

The index of the root of the tree. (default = 1)

reduction_method : str

The graph reduction method (default = ‘resistance_distance’)

compute_full_eigen : bool

To also compute the graph Laplacian eigenvalues for every tree in the sequence

Returns:

Gs : ndarray

Ndarray, with each element containing a graph structure represent a reduced tree.

subsampled_vertex_indices : ndarray

Indices of the vertices of the previous tree that are kept for the subsequent tree.

Features¶

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

pygsp.features.compute_avg_adj_deg(G)[source]¶

Compute the average adjacency degree for each node.

The average adjacency degree is the average of the degrees of a node and its neighbors.

Parameters:

G: Graph

Graph on which the statistic is extracted

pygsp.features.compute_norm_tig(g, **kwargs)[source]¶

Compute the $\ell_2$ norm of the Tig. See compute_tig().

Parameters:

g: Filter

The filter or filter bank.

kwargs: dict

Additional parameters to be passed to the pygsp.filters.Filter.filter() method.

pygsp.features.compute_spectrogram(G, atom=None, M=100, **kwargs)[source]¶

Compute the norm of the Tig for all nodes with a kernel shifted along the spectral axis.

Parameters:

G : Graph

Graph on which to compute the spectrogram.

atom : func

Kernel to use in the spectrogram (default = exp(-M*(x/lmax)²)).

M : int (optional)

Number of samples on the spectral scale. (default = 100)

kwargs: dict

Additional parameters to be passed to the pygsp.filters.Filter.filter() method.

pygsp.features.compute_tig(g, **kwargs)[source]¶

Compute the Tig for a given filter or filter bank.

\[T_ig(n) = g(L)_{i, n}\]

Parameters:

g: Filter

One of pygsp.filters.

kwargs: dict

Additional parameters to be passed to the pygsp.filters.Filter.filter() method.

Optimization¶

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

pygsp.optimization.prox_tv(x, gamma, G, A=None, At=None, nu=1, tol=0.001, maxit=200, use_matrix=True)[source]¶

Total Variation proximal operator for graphs.

This function computes the TV proximal operator for graphs. The TV norm is the one norm of the gradient. The gradient is defined in the function pygsp.graphs.Graph.grad(). This function requires the PyUNLocBoX to be executed.

This function solves:

$sol = \min_{z} \frac{1}{2} \|x - z\|_2^2 + \gamma \|x\|_{TV}$

Parameters:

x: int

Input signal

gamma: ndarray

Regularization parameter

G: graph object

Graphs structure

A: lambda function

Forward operator, this parameter allows to solve the following problem: $sol = \min_{z} \frac{1}{2} \|x - z\|_2^2 + \gamma \| A x\|_{TV}$ (default = Id)

At: lambda function

Adjoint operator. (default = Id)

nu: float

Bound on the norm of the operator (default = 1)

tol: float

Stops criterion for the loop. The algorithm will stop if : $\frac{n(t) - n(t - 1)} {n(t)} < tol$ where $n(t) = f(x) + 0.5 \|x-y\|_2^2$ is the objective function at iteration $t$ (default = $10e-4$)

maxit: int

Maximum iteration. (default = 200)

use_matrix: bool

If a matrix should be used. (default = True)

Returns:

sol: solution

Utilities¶

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

pygsp.utils.build_logger(name)[source]¶

pygsp.utils.compute_log_scales(lmin, lmax, Nscales, t1=1, t2=2)[source]¶

Compute logarithm scales for wavelets.

Parameters:

lmin : float

Smallest non-zero eigenvalue.

lmax : float

Largest eigenvalue, i.e. pygsp.graphs.Graph.lmax.

Nscales : int

Number of scales.

Returns:

scales : ndarray

List of scales of length Nscales.

Examples

>>> from pygsp import utils
>>> utils.compute_log_scales(1, 10, 3)
array([ 2.       ,  0.4472136,  0.1      ])

pygsp.utils.distanz(x, y=None)[source]¶

Calculate the distance between two colon vectors.

Parameters:

x : ndarray

First colon vector

y : ndarray

Second colon vector

Returns:

d : ndarray

Distance between x and y

Examples

>>> from pygsp import utils
>>> x = np.arange(3)
>>> utils.distanz(x, x)
array([[ 0.,  1.,  2.],
       [ 1.,  0.,  1.],
       [ 2.,  1.,  0.]])

pygsp.utils.filterbank_handler(func)[source]¶

pygsp.utils.graph_array_handler(func)[source]¶

pygsp.utils.import_classes(names, src, dst)[source]¶: Import classes in package from their implementation modules.

pygsp.utils.import_functions(names, src, dst)[source]¶: Import functions in package from their implementation modules.

pygsp.utils.import_modules(names, src, dst)[source]¶: Import modules in package.

pygsp.utils.loadmat(path)[source]¶

Load a matlab data file.

Parameters:

path : string

Path to the mat file from the data folder, without the .mat extension.

Returns:

data : dict

dictionary with variable names as keys, and loaded matrices as values.

Examples

>>> from pygsp import utils
>>> data = utils.loadmat('pointclouds/bunny')
>>> data['bunny'].shape
(2503, 3)

pygsp.utils.repmatline(A, ncol=1, nrow=1)[source]¶

Repeat the matrix A in a specific manner.

Parameters:

A : ndarray

ncol : int

default is 1

nrow : int

default is 1

Returns:

Ar : ndarray

Examples

>>> from pygsp import utils
>>> x = np.array([[1, 2], [3, 4]])
>>> x
array([[1, 2],
       [3, 4]])
>>> utils.repmatline(x, nrow=2, ncol=3)
array([[1, 1, 1, 2, 2, 2],
       [1, 1, 1, 2, 2, 2],
       [3, 3, 3, 4, 4, 4],
       [3, 3, 3, 4, 4, 4]])

pygsp.utils.rescale_center(x)[source]¶

Rescale and center data, e.g. embedding coordinates.

Parameters:

x : ndarray

Data to be rescaled.

Returns:

r : ndarray

Rescaled data.

Examples

>>> from pygsp import utils
>>> x = np.array([[1, 6], [2, 5], [3, 4]])
>>> utils.rescale_center(x)
array([[-1. ,  1. ],
       [-0.6,  0.6],
       [-0.2,  0.2]])

pygsp.utils.resistance_distance(G)[source]¶

Compute the resistance distances of a graph.

Parameters:

G : Graph or sparse matrix

Graph structure or Laplacian matrix (L)

Returns:

rd : sparse matrix

distance matrix

References

[KRandic93]

pygsp.utils.sparsifier(func)[source]¶

pygsp.utils.symmetrize(W, method='average')[source]¶

Symmetrize a square matrix.

Parameters:

W : array_like

Square matrix to be symmetrized

method : string

‘average’ : symmetrize by averaging with the transpose. Most useful when transforming a directed graph to an undirected one.

‘maximum’ : symmetrize by taking the maximum with the transpose. Similar to ‘fill’ except that ambiguous entries are resolved by taking the largest value.

‘fill’ : symmetrize by filling in the zeros in both the upper and lower triangular parts. Ambiguous entries are resolved by averaging the values.

‘tril’ : symmetrize by considering the lower triangular part only.

‘triu’ : symmetrize by considering the upper triangular part only.

Examples

>>> from pygsp import utils
>>> W = np.array([[0, 3, 0], [3, 1, 6], [4, 2, 3]], dtype=float)
>>> W
array([[ 0.,  3.,  0.],
       [ 3.,  1.,  6.],
       [ 4.,  2.,  3.]])
>>> utils.symmetrize(W, method='average')
array([[ 0.,  3.,  2.],
       [ 3.,  1.,  4.],
       [ 2.,  4.,  3.]])
>>> utils.symmetrize(W, method='maximum')
array([[ 0.,  3.,  4.],
       [ 3.,  1.,  6.],
       [ 4.,  6.,  3.]])
>>> utils.symmetrize(W, method='fill')
array([[ 0.,  3.,  4.],
       [ 3.,  1.,  4.],
       [ 4.,  4.,  3.]])
>>> utils.symmetrize(W, method='tril')
array([[ 0.,  3.,  4.],
       [ 3.,  1.,  2.],
       [ 4.,  2.,  3.]])
>>> utils.symmetrize(W, method='triu')
array([[ 0.,  3.,  0.],
       [ 3.,  1.,  6.],
       [ 0.,  6.,  3.]])

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.

[DB13]

Florian Dorfler and Francesco Bullo. Kron reduction of graphs with applications to electrical networks. Circuits and Systems I: Regular Papers, IEEE Transactions on, 60(1):150–163, 2013.

[Gle]

D. Gleich. The MatlabBGL Matlab library. http://www.cs.purdue.edu/homes/dgleich/packages/matlab_bgl/index.html.

[HVG11]

D. K. Hammond, P. Vandergheynst, and R. Gribonval. Wavelets on graphs via spectral graph theory. Appl. Comput. Harmon. Anal., 30(2):129–150, Mar. 2011.

[KV03]

Jeong Han Kim and Van H Vu. Generating random regular graphs. In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing. 2003.

[KRandic93]

Douglas J Klein and M Randić. Resistance distance. Journal of Mathematical Chemistry, 12(1):81–95, 1993.

[LVDV11]

Nora Leonardi and Dimitri Van De Ville. Wavelet frames on graphs defined by fmri functional connectivity. In Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on, 2136–2139. IEEE, 2011.

[Pes09]

Isaac Pesenson. Variational splines and paley–wiener spaces on combinatorial graphs. Constructive Approximation, 29(1):1–21, 2009.

[Rud99]

Mark Rudelson. Random vectors in the isotropic position. Journal of Functional Analysis, 164(1):60–72, 1999.

[RV07]

Mark Rudelson and Roman Vershynin. Sampling from large matrices: an approach through geometric functional analysis. Journal of the ACM (JACM), 54(4):21, 2007.

[SFV13]

David I Shuman, Mohammad Javad Faraji, and Pierre Vandergheynst. A framework for multiscale transforms on graphs. arXiv preprint arXiv:1308.4942, 2013.

[SS11]

Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. SIAM Journal on Computing, 40(6):1913–1926, 2011.

[Str99]

Gilbert Strang. The discrete cosine transform. SIAM review, 41(1):135–147, 1999.

[TPGV16]

Nicolas Tremblay, Gilles Puy, Rémi Gribonval, and Pierre Vandergheynst. Compressive spectral clustering. arXiv preprint arXiv:1602.02018, 2016.

[TL94]

Greg Turk and Marc Levoy. Zippered polygon meshes from range images. In Proceedings of the 21st annual conference on Computer graphics and interactive techniques, 311–318. ACM, 1994.