PyGSP: Graph Signal Processing in Python¶

The PyGSP is a Python package to ease Signal Processing on Graphs. The documentation is available on Read the Docs and development takes place on GitHub. A (mostly unmaintained) Matlab version 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.

While NetworkX and graph-tool are tools to analyze the topology of graphs, the aim of the PyGSP is to analyze graph signals, also known as features or properties (i.e., not the graph itself). Those three tools are complementary and work well together with the provided import / export facility.

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.compute_fourier_basis()  # Fourier to plot the eigenvalues.
>>> # G.estimate_lmax() is otherwise sufficient.
>>> g = filters.Heat(G, scale=50)
>>> fig, ax = g.plot()

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)
>>> fig, ax = G.plot(s, highlight=DELTAS)

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!

Installation¶

The PyGSP is available on PyPI:

$ pip install pygsp

The PyGSP is available on conda-forge:

$ conda install -c conda-forge pygsp

The PyGSP is available in the Arch User Repository:

$ git clone https://aur.archlinux.org/python-pygsp.git
$ cd python-pygsp
$ makepkg -csi

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

It is released under the terms of the BSD 3-Clause license.

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

@misc{pygsp,
  title = {PyGSP: Graph Signal Processing in Python},
  author = {Defferrard, Micha\"el and Martin, Lionel and Pena, Rodrigo and Perraudin, Nathana\"el},
  doi = {10.5281/zenodo.1003157},
  url = {https://github.com/epfl-lts2/pygsp/},
}

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.

>>> rng = np.random.default_rng(1)  # Reproducible results.
>>> W = rng.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, 64 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.csr.csr_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
(30, 64)

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')
>>> fig, ax = 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(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(G2.U[:, 4], ax=axes[0])
>>> G2.set_coordinates('line1D')
>>> _ = G2.plot(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 = g.plot(eigenvalues=True, 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
>>> s[G.info['idx_s']] = 0
>>> s[G.info['idx_p']] = 1
>>> s += rng.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(s, vertex_size=30, title='noisy', ax=axes[0])
>>> axes[0].set_axis_off()
>>> _ = G.plot(s2, vertex_size=30, title='cleaned', ax=axes[1])
>>> 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')
...     title = r'Heat diffusion, $\tau={}$'.format(taus[i])
...     _ = G.plot(s[:, i], colorbar=False, title=title, ax=ax)
...     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(title='Filter bank of mexican hat wavelets', ax=ax)

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(s[:, i], title='Wavelet {}'.format(i+1), ax=ax)
...     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')
...     title = 'Curvature estimation (scale {})'.format(i+1)
...     _ = G.plot(s[:, i], title=title, ax=ax)
...     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, distributed=True, seed=42)
>>> G.compute_fourier_basis()
>>>
>>> # Create label signal
>>> label_signal = np.copysign(np.ones(G.N), G.U[:, 3])
>>>
>>> fig, ax = G.plot(label_signal)

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

>>> rng = np.random.default_rng(42)
>>>
>>> # Create the mask
>>> M = rng.uniform(size=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 * rng.standard_normal(G.N))
>>>
>>> fig, ax = G.plot(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.T.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.213139e+02
    stopping criterion: MAXIT
>>>
>>> fig, ax = G.plot(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) = 6.422673e+01
    stopping criterion: MAXIT
>>>
>>> fig, ax = G.plot(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 400 nodes.

>>> G = graphs.Sensor(400, distributed=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)

Examples¶

Kernel localization¶

Note

Click here to download the full example code

Kernel localization¶

In classical signal processing, a filter can be translated in the vertex domain. We cannot do that on graphs. Instead, we can localize() a filter kernel. Note how on classic structures (like the ring), the localized kernel is the same everywhere, while it changes when localized on irregular graphs.

$heat kernel, $g(L) \delta_{0}$, $g(L) \delta_{10}$, $g(L) \delta_{20}$, heat kernel, $g(L) \delta_{0}$, $g(L) \delta_{10}$, $g(L) \delta_{20}$$

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

fig, axes = plt.subplots(2, 4, figsize=(10, 4))

graphs = [
    pg.graphs.Ring(40),
    pg.graphs.Sensor(64, seed=42),
]

locations = [0, 10, 20]

for graph, axs in zip(graphs, axes):
    graph.compute_fourier_basis()
    g = pg.filters.Heat(graph)
    g.plot(ax=axs[0], title='heat kernel')
    axs[0].set_xlabel(r'eigenvalues $\lambda$')
    axs[0].set_ylabel(r'$g(\lambda) = \exp \left( \frac{{-{}\lambda}}{{\lambda_{{max}}}} \right)$'.format(g.scale[0]))
    maximum = 0
    for loc in locations:
        x = g.localize(loc)
        maximum = np.maximum(maximum, x.max())
    for loc, ax in zip(locations, axs[1:]):
        graph.plot(g.localize(loc), limits=[0, maximum], highlight=loc, ax=ax,
                   title=r'$g(L) \delta_{{{}}}$'.format(loc))
        ax.set_axis_off()

fig.tight_layout()

Total running time of the script: ( 0 minutes 1.012 seconds)

Estimated memory usage: 25 MB

Gallery generated by Sphinx-Gallery

Fourier basis¶

Note

Click here to download the full example code

Fourier basis¶

The eigenvectors of the graph Laplacian form the Fourier basis. The eigenvalues are a measure of variation of their corresponding eigenvector. The lower the eigenvalue, the smoother the eigenvector. They are hence a measure of “frequency”.

In classical signal processing, Fourier modes are completely delocalized, like on the grid graph. For general graphs however, Fourier modes might be localized. See pygsp.graphs.Graph.coherence.

$u_1^\top L u_1 = 0.00$, $u_2^\top L u_2 = 0.10$, $u_3^\top L u_3 = 0.10$, $u_4^\top L u_4 = 0.20$, $u_5^\top L u_5 = 0.38$, $u_6^\top L u_6 = 0.38$, $u_7^\top L u_7 = 0.48$, $u_1^\top L u_1 = 0.00$, $u_2^\top L u_2 = 0.10$, $u_3^\top L u_3 = 0.23$, $u_4^\top L u_4 = 0.29$, $u_5^\top L u_5 = 0.52$, $u_6^\top L u_6 = 1.01$, $u_7^\top L u_7 = 1.20$

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

n_eigenvectors = 7

fig, axes = plt.subplots(2, 7, figsize=(15, 4))

def plot_eigenvectors(G, axes):
    G.compute_fourier_basis(n_eigenvectors)
    limits = [f(G.U) for f in (np.min, np.max)]
    for i, ax in enumerate(axes):
        G.plot(G.U[:, i], limits=limits, colorbar=False, vertex_size=50, ax=ax)
        energy = abs(G.dirichlet_energy(G.U[:, i]))
        ax.set_title(r'$u_{0}^\top L u_{0} = {1:.2f}$'.format(i+1, energy))
        ax.set_axis_off()

G = pg.graphs.Grid2d(10, 10)
plot_eigenvectors(G, axes[0])
fig.subplots_adjust(hspace=0.5, right=0.8)
cax = fig.add_axes([0.82, 0.60, 0.01, 0.26])
fig.colorbar(axes[0, -1].collections[1], cax=cax, ticks=[-0.2, 0, 0.2])

G = pg.graphs.Sensor(seed=42)
plot_eigenvectors(G, axes[1])
fig.subplots_adjust(hspace=0.5, right=0.8)
cax = fig.add_axes([0.82, 0.16, 0.01, 0.26])
_ = fig.colorbar(axes[1, -1].collections[1], cax=cax, ticks=[-0.4, 0, 0.4])

Total running time of the script: ( 0 minutes 0.566 seconds)

Estimated memory usage: 16 MB

Gallery generated by Sphinx-Gallery

Concentration of the eigenvalues¶

Note

Click here to download the full example code

Concentration of the eigenvalues¶

The eigenvalues of the graph Laplacian concentrates to the same value as the graph becomes full.

Ring(n_vertices=17, n_edges=17, k=1), Ring(n_vertices=17, n_edges=34, k=2), Ring(n_vertices=17, n_edges=85, k=5), Ring(n_vertices=17, n_edges=136, k=8), k=1, k=2, k=5, k=8

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

n_neighbors = [1, 2, 5, 8]
fig, axes = plt.subplots(3, len(n_neighbors), figsize=(15, 8))

for k, ax in zip(n_neighbors, axes.T):

    graph = pg.graphs.Ring(17, k=k)
    graph.compute_fourier_basis()
    graph.plot(graph.U[:, 1], ax=ax[0])
    ax[0].axis('equal')
    ax[1].spy(graph.W)
    ax[2].plot(graph.e, '.')
    ax[2].set_title('k={}'.format(k))
    #graph.set_coordinates('line1D')
    #graph.plot(graph.U[:, :4], ax=ax[3], title='')

    # Check that the DFT matrix is an eigenbasis of the Laplacian.
    U = np.fft.fft(np.identity(graph.n_vertices))
    LambdaM = (graph.L.todense().dot(U)) / (U + 1e-15)
    # Eigenvalues should be real.
    assert np.all(np.abs(np.imag(LambdaM)) < 1e-10)
    LambdaM = np.real(LambdaM)
    # Check that the eigenvectors are really eigenvectors of the laplacian.
    Lambda = np.mean(LambdaM, axis=0)
    assert np.all(np.abs(LambdaM - Lambda) < 1e-10)

fig.tight_layout()

Total running time of the script: ( 0 minutes 1.414 seconds)

Estimated memory usage: 14 MB

Gallery generated by Sphinx-Gallery

Fourier transform¶

Note

Click here to download the full example code

Fourier transform¶

The graph Fourier transform pygsp.graphs.Graph.gft() transforms a signal from the vertex domain to the spectral domain. The smoother the signal (see pygsp.graphs.Graph.dirichlet_energy()), the lower in the frequencies its energy is concentrated.

$x^T L x = 0.30$, $x^T L x = 1.55$, $x^T L x = 5.43$

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

G = pg.graphs.Sensor(seed=42)
G.compute_fourier_basis()

scales = [10, 3, 0]
limit = 0.44

fig, axes = plt.subplots(2, len(scales), figsize=(12, 4))
fig.subplots_adjust(hspace=0.5)

x0 = np.random.default_rng(1).normal(size=G.N)
for i, scale in enumerate(scales):
    g = pg.filters.Heat(G, scale)
    x = g.filter(x0).squeeze()
    x /= np.linalg.norm(x)
    x_hat = G.gft(x).squeeze()

    assert np.all((-limit < x) & (x < limit))
    G.plot(x, limits=[-limit, limit], ax=axes[0, i])
    axes[0, i].set_axis_off()
    axes[0, i].set_title('$x^T L x = {:.2f}$'.format(G.dirichlet_energy(x)))

    axes[1, i].plot(G.e, np.abs(x_hat), '.-')
    axes[1, i].set_xticks(range(0, 16, 4))
    axes[1, i].set_xlabel(r'graph frequency $\lambda$')
    axes[1, i].set_ylim(-0.05, 0.95)

axes[1, 0].set_ylabel(r'frequency content $\hat{x}(\lambda)$')

# axes[0, 0].set_title(r'$x$: signal in the vertex domain')
# axes[1, 0].set_title(r'$\hat{x}$: signal in the spectral domain')

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.596 seconds)

Estimated memory usage: 10 MB

Gallery generated by Sphinx-Gallery

Heat diffusion¶

Note

Click here to download the full example code

Heat diffusion¶

Solve the heat equation by filtering the initial conditions with the heat kernel pygsp.filters.Heat.

$\hat{f}(0) = g_{1,0} \odot \hat{f}(0)$, $\hat{f}(5) = g_{1,5} \odot \hat{f}(0)$, $\hat{f}(10) = g_{1,10} \odot \hat{f}(0)$, $\hat{f}(20) = g_{1,20} \odot \hat{f}(0)$, $f(0)$, $f(5)$, $f(10)$, $f(20)$

from os import path

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

n_side = 13
G = pg.graphs.Grid2d(n_side)
G.compute_fourier_basis()

sources = [
    (n_side//4 * n_side) + (n_side//4),
    (n_side*3//4 * n_side) + (n_side*3//4),
]
x = np.zeros(G.n_vertices)
x[sources] = 5

times = [0, 5, 10, 20]

fig, axes = plt.subplots(2, len(times), figsize=(12, 5))
for i, t in enumerate(times):
    g = pg.filters.Heat(G, scale=t)
    title = r'$\hat{{f}}({0}) = g_{{1,{0}}} \odot \hat{{f}}(0)$'.format(t)
    g.plot(alpha=1, ax=axes[0, i], title=title)
    axes[0, i].set_xlabel(r'$\lambda$')
#    axes[0, i].set_ylabel(r'$g(\lambda)$')
    if i > 0:
        axes[0, i].set_ylabel('')
    y = g.filter(x)
    line, = axes[0, i].plot(G.e, G.gft(y))
    labels = [r'$\hat{{f}}({})$'.format(t), r'$g_{{1,{}}}$'.format(t)]
    axes[0, i].legend([line, axes[0, i].lines[-3]], labels, loc='lower right')
    G.plot(y, edges=False, highlight=sources, ax=axes[1, i], title=r'$f({})$'.format(t))
    axes[1, i].set_aspect('equal', 'box')
    axes[1, i].set_axis_off()

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.950 seconds)

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery

Wave propagation¶

Note

Click here to download the full example code

Wave propagation¶

Solve the wave equation by filtering the initial conditions with the wave kernel pygsp.filters.Wave.

$\hat{f}(0) = g_{1,0} \odot \hat{f}(0)$, $\hat{f}(5) = g_{1,5} \odot \hat{f}(0)$, $\hat{f}(10) = g_{1,10} \odot \hat{f}(0)$, $\hat{f}(20) = g_{1,20} \odot \hat{f}(0)$, $f(0)$, $f(5)$, $f(10)$, $f(20)$

from os import path

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

n_side = 13
G = pg.graphs.Grid2d(n_side)
G.compute_fourier_basis()

sources = [
    (n_side//4 * n_side) + (n_side//4),
    (n_side*3//4 * n_side) + (n_side*3//4),
]
x = np.zeros(G.n_vertices)
x[sources] = 5

times = [0, 5, 10, 20]

fig, axes = plt.subplots(2, len(times), figsize=(12, 5))
for i, t in enumerate(times):
    g = pg.filters.Wave(G, time=t, speed=1)
    title = r'$\hat{{f}}({0}) = g_{{1,{0}}} \odot \hat{{f}}(0)$'.format(t)
    g.plot(alpha=1, ax=axes[0, i], title=title)
    axes[0, i].set_xlabel(r'$\lambda$')
#    axes[0, i].set_ylabel(r'$g(\lambda)$')
    if i > 0:
        axes[0, i].set_ylabel('')
    y = g.filter(x)
    line, = axes[0, i].plot(G.e, G.gft(y))
    labels = [r'$\hat{{f}}({})$'.format(t), r'$g_{{1,{}}}$'.format(t)]
    axes[0, i].legend([line, axes[0, i].lines[-3]], labels, loc='lower right')
    G.plot(y, edges=False, highlight=sources, ax=axes[1, i], title=r'$f({})$'.format(t))
    axes[1, i].set_aspect('equal', 'box')
    axes[1, i].set_axis_off()

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.960 seconds)

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery

Localization of Fourier modes¶

Note

Click here to download the full example code

Localization of Fourier modes¶

The Fourier modes (the eigenvectors of the graph Laplacian) can be localized in the spacial domain. As a consequence, graph signals can be localized in both space and frequency (which is impossible for Euclidean domains or manifolds, by the Heisenberg’s uncertainty principle).

This example demonstrates that the more isolated a node is, the more a Fourier mode will be localized on it.

The mutual coherence between the basis of Kronecker deltas and the basis formed by the eigenvectors of the Laplacian, pygsp.graphs.Graph.coherence, is a measure of the localization of the Fourier modes. The larger the value, the more localized the eigenvectors can be.

See Global and Local Uncertainty Principles for Signals on Graphs for details.

import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
import pygsp as pg

fig, axes = plt.subplots(2, 2, figsize=(8, 8))

for w, ax in zip([10, 1, 0.1, 0.01], axes.flatten()):

    adjacency = [
        [0, w, 0, 0],
        [w, 0, 1, 0],
        [0, 1, 0, 1],
        [0, 0, 1, 0],
    ]
    graph = pg.graphs.Graph(adjacency)
    graph.compute_fourier_basis()

    # Plot eigenvectors.
    ax.plot(graph.U)
    ax.set_ylim(-1, 1)
    ax.set_yticks([-1, 0, 1])
    ax.legend([f'$u_{i}(v)$, $\lambda_{i}={graph.e[i]:.1f}$' for i in
              range(graph.n_vertices)], loc='upper right')

    ax.text(0, -0.9, f'coherence = {graph.coherence:.2f}'
            f'$\in [{1/np.sqrt(graph.n_vertices)}, 1]$')

    # Plot vertices.
    ax.set_xticks(range(graph.n_vertices))
    ax.set_xticklabels([f'$v_{i}$' for i in range(graph.n_vertices)])

    # Plot graph.
    x, y = np.arange(0, graph.n_vertices), -1.20*np.ones(graph.n_vertices)
    line = mpl.lines.Line2D(x, y, lw=3, color='k', marker='.', markersize=20)
    line.set_clip_on(False)
    ax.add_line(line)

    # Plot edge weights.
    for i in range(graph.n_vertices - 1):
        j = i+1
        ax.text(i+0.5, -1.15, f'$w_{{{i}{j}}} = {adjacency[i][j]}$',
                horizontalalignment='center')

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.559 seconds)

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery

Filtering a signal¶

Note

Click here to download the full example code

Filtering a signal¶

A graph signal is filtered by transforming it to the spectral domain (via the Fourier transform), performing a point-wise multiplication (motivated by the convolution theorem), and transforming it back to the vertex domain (via the inverse graph Fourier transform).

Note

In practice, filtering is implemented in the vertex domain to avoid the computationally expensive graph Fourier transform. To do so, filters are implemented as polynomials of the eigenvalues / Laplacian. Hence, filtering a signal reduces to its multiplications with sparse matrices (the graph Laplacian).

input signal $x$ in the vertex domain, signals in the spectral domain, filtered signal $y$ in the vertex domain

import numpy as np
from matplotlib import pyplot as plt
import pygsp as pg

G = pg.graphs.Sensor(seed=42)
G.compute_fourier_basis()

#g = pg.filters.Rectangular(G, band_max=0.2)
g = pg.filters.Expwin(G, band_max=0.5)

fig, axes = plt.subplots(1, 3, figsize=(12, 4))
fig.subplots_adjust(hspace=0.5)

x = np.random.default_rng(1).normal(size=G.N)
#x = np.random.default_rng(42).uniform(-1, 1, size=G.N)
x = 3 * x / np.linalg.norm(x)
y = g.filter(x)
x_hat = G.gft(x).squeeze()
y_hat = G.gft(y).squeeze()

limits = [x.min(), x.max()]

G.plot(x, limits=limits, ax=axes[0], title='input signal $x$ in the vertex domain')
axes[0].text(0, -0.1, '$x^T L x = {:.2f}$'.format(G.dirichlet_energy(x)))
axes[0].set_axis_off()

g.plot(ax=axes[1], alpha=1)
line_filt = axes[1].lines[-2]
line_in, = axes[1].plot(G.e, np.abs(x_hat), '.-')
line_out, = axes[1].plot(G.e, np.abs(y_hat), '.-')
#axes[1].set_xticks(range(0, 16, 4))
axes[1].set_xlabel(r'graph frequency $\lambda$')
axes[1].set_ylabel(r'frequency content $\hat{x}(\lambda)$')
axes[1].set_title(r'signals in the spectral domain')
axes[1].legend(['input signal $\hat{x}$'])
labels = [
    r'input signal $\hat{x}$',
    'kernel $g$',
    r'filtered signal $\hat{y}$',
]
axes[1].legend([line_in, line_filt, line_out], labels, loc='upper right')

G.plot(y, limits=limits, ax=axes[2], title='filtered signal $y$ in the vertex domain')
axes[2].text(0, -0.1, '$y^T L y = {:.2f}$'.format(G.dirichlet_energy(y)))
axes[2].set_axis_off()

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.620 seconds)

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery

Random walks¶

Note

Click here to download the full example code

Random walks¶

Probability of a random walker to be on any given vertex after a given number of steps starting from a given distribution.

# sphinx_gallery_thumbnail_number = 2

import numpy as np
from scipy import sparse
from matplotlib import pyplot as plt
import pygsp as pg

N = 7
steps = [0, 1, 2, 3]

graph = pg.graphs.Grid2d(N)
delta = np.zeros(graph.N)
delta[N//2*N + N//2] = 1

probability = sparse.diags(graph.dw**(-1)).dot(graph.W)

fig, axes = plt.subplots(1, len(steps), figsize=(12, 3))
for step, ax in zip(steps, axes):
    state = (probability**step).__rmatmul__(delta)  ## = delta @ probability**step
    graph.plot(state, ax=ax, title=r'$\delta P^{}$'.format(step))
    ax.set_axis_off()

fig.tight_layout()

$\delta P^0$, $\delta P^1$, $\delta P^2$, $\delta P^3$

Stationary distribution.

graphs = [
    pg.graphs.Ring(10),
    pg.graphs.Grid2d(5),
    pg.graphs.Comet(8, 4),
    pg.graphs.BarabasiAlbert(20, seed=42),
]

fig, axes = plt.subplots(1, len(graphs), figsize=(12, 3))

for graph, ax in zip(graphs, axes):

    if not hasattr(graph, 'coords'):
        graph.set_coordinates(seed=10)

    P = sparse.diags(graph.dw**(-1)).dot(graph.W)

#    e, u = np.linalg.eig(P.T.toarray())
#    np.testing.assert_allclose(np.linalg.inv(u.T) @ np.diag(e) @ u.T,
#                               P.toarray(), atol=1e-10)
#    np.testing.assert_allclose(np.abs(e[0]), 1)
#    stationary = np.abs(u.T[0])

    e, u = sparse.linalg.eigs(P.T, k=1, which='LR')
    np.testing.assert_allclose(e, 1)
    stationary = np.abs(u).squeeze()
    assert np.all(stationary < 0.71)

    colorbar = False if type(graph) is pg.graphs.Ring else True
    graph.plot(stationary, colorbar=colorbar, ax=ax, title='$xP = x$')
    ax.set_axis_off()

fig.tight_layout()

Total running time of the script: ( 0 minutes 0.965 seconds)

Estimated memory usage: 9 MB

Gallery generated by Sphinx-Gallery

API reference¶

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,
learning to solve learning problems,
optimization to help solving convex optimization problems,
utils for various utilities.

pygsp.test()[source]¶: Run the test suite.

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.

Attributes¶

Matrix operators

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

Vectors

`Graph.d`	The degree (number of neighbors) of vertices.
`Graph.dw`	The weighted degree of vertices.
`Graph.e`	Eigenvalues of the Laplacian (square of graph frequencies).

Scalars

`Graph.lmax`	Largest eigenvalue of the graph Laplacian.
`Graph.coherence`	Coherence of the Fourier basis.

Attributes computation¶

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

Differential operators¶

`Graph.grad`(x)	Compute the gradient of a signal defined on the vertices.
`Graph.div`(y)	Compute the divergence of a signal defined on the edges.
`Graph.dirichlet_energy`(x)	Compute the Dirichlet energy of a signal defined on the vertices.

Transforms¶

`Graph.gft`(s)	Compute the graph Fourier transform.
`Graph.igft`(s_hat)	Compute the inverse graph Fourier transform.

Vertex-frequency transforms are implemented as filter banks and are found in pygsp.filters (such as Gabor and Modulation).

Checks¶

`Graph.is_weighted`()	Check if the graph is weighted.
`Graph.is_connected`()	Check if the graph is connected (cached).
`Graph.is_directed`()	Check if the graph has directed edges (cached).
`Graph.has_loops`()	Check if any vertex is connected to itself.

Plotting¶

`Graph.plot`([vertex_color, vertex_size, …])	Plot a graph with signals as color or vertex size.
`Graph.plot_spectrogram`([node_idx])	Plot the graph’s spectrogram.

Import and export (I/O)¶

We provide import and export facility to two well-known Python packages for network analysis: NetworkX and graph-tool. Those packages and the PyGSP are fundamentally different in their goals (graph analysis versus graph signal analysis) and graph representations (if in the PyGSP everything is an ndarray, in NetworkX everything is a dictionary). Those tools are complementary and good interoperability is necessary to exploit the strengths of each tool. We ourselves leverage NetworkX and graph-tool to save and load graphs.

Note: to tie a signal with the graph, such that they are exported together, attach it first with Graph.set_signal().

`Graph.load`(path[, fmt, backend])	Load a graph from a file.
`Graph.save`(path[, fmt, backend])	Save the graph to a file.
`Graph.from_networkx`(graph[, weight])	Import a graph from NetworkX.
`Graph.to_networkx`()	Export the graph to NetworkX.
`Graph.from_graphtool`(graph[, weight])	Import a graph from graph-tool.
`Graph.to_graphtool`()	Export the graph to graph-tool.

Others¶

`Graph.get_edge_list`()	Return an edge list, an alternative representation of the graph.
`Graph.set_signal`(signal, name)	Attach a signal to the graph.
`Graph.set_coordinates`([kind, seed])	Set node’s coordinates (their position when plotting).
`Graph.subgraph`(vertices)	Create a subgraph from a list of vertices.
`Graph.extract_components`()	Split the graph into connected components.

Graph models¶

In addition to the below graphs, useful resources are the random graph generators from NetworkX (see NetworkX’s documentation) and graph-tool (see graph_tool.generation), as well as graph-tool’s assortment of standard networks (see graph_tool.collection). Any graph created by NetworkX or graph-tool can be imported in the PyGSP with Graph.from_networkx() and Graph.from_graphtool().

Graphs built from other graphs¶

LineGraph(graph, **kwargs)

Build the line graph of a graph.

Generated graphs¶

`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, diagonal])	2-dimensional grid graph.
`Logo`(**kwargs)	GSP logo.
`LowStretchTree`([k])	Low stretch tree.
`Minnesota`([connected])	Minnesota road network (from MatlabBGL).
`Path`([N, directed])	Path graph.
`RandomRegular`([N, k, max_iter, seed])	Random k-regular graph.
`RandomRing`([N, angles, seed])	Ring graph with randomly sampled vertices.
`Ring`([N, k])	K-regular ring graph.
`Star`([N])	Star 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.
`Sensor`([N, k, distributed, seed])	Random sensor 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(adjacency, lap_type='combinatorial', coords=None, plotting={})[source]¶

Base graph class.

Instantiate it to construct a graph from a (weighted) adjacency matrix.
Provide a common interface (and implementation) for graph objects.
Initialize attributes for derived classes.

Parameters

adjacencysparse matrix or array_like: The (weighted) adjacency matrix of size n_vertices by n_vertices that encodes the graph. The data is copied except if it is a sparse matrix in CSR format.
lap_type{‘combinatorial’, ‘normalized’}: The kind of Laplacian to be computed by compute_laplacian().
coordsarray_like: A matrix of size n_vertices by d that represents the coordinates of the vertices in a d-dimensional embedding space.
plottingdict: Plotting parameters.

Examples

Define a simple graph.

>>> graph = graphs.Graph([
...     [0., 2., 0.],
...     [2., 0., 5.],
...     [0., 5., 0.],
... ])
>>> graph
Graph(n_vertices=3, n_edges=2)
>>> graph.n_vertices, graph.n_edges
(3, 2)
>>> graph.W.toarray()
array([[0., 2., 0.],
       [2., 0., 5.],
       [0., 5., 0.]])
>>> graph.d
array([1, 2, 1], dtype=int32)
>>> graph.dw
array([2., 7., 5.])
>>> graph.L.toarray()
array([[ 2., -2.,  0.],
       [-2.,  7., -5.],
       [ 0., -5.,  5.]])

Add some coordinates to plot it.

>>> import matplotlib.pyplot as plt
>>> graph.set_coordinates([
...     [0, 0],
...     [0, 1],
...     [1, 0],
... ])
>>> fig, ax = graph.plot()

Attributes

n_vertices or Nint: The number of vertices (nodes) in the graph.
n_edges or Neint: The number of edges (links) in the graph.
Wscipy.sparse.csr_matrix: Weighted adjacency matrix of the graph.
Lscipy.sparse.csr_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().
signalsdict (string -> numpy.ndarray): Signals attached to the graph.
coordsnumpy.ndarray: Vertices coordinates in 2D or 3D space. Used for plotting only.
plottingdict: Plotting parameters.

compute_differential_operator()¶

Compute the graph differential operator (cached).

The differential operator is the matrix $D$ such that

\[L = D D^\top,\]

where $L$ is the graph Laplacian (combinatorial or normalized). It is used to compute the gradient and the divergence of graph signals (see grad() and div()).

The differential operator computes the gradient and divergence of signals, and the Laplacian computes the divergence of the gradient, as follows:

\[z = L x = D y = D D^\top x,\]

where $y = D^\top x = \nabla_\mathcal{G} x$ is the gradient of $x$ and $z = D y = \operatorname{div}_\mathcal{G} y$ is the divergence of $y$. See grad() and div() for details.

The difference operator is actually an incidence matrix of the graph, defined as

\[\begin{split}D[i, k] = \begin{cases} -\sqrt{W[i, j] / 2} & \text{if } e_k = (v_i, v_j) \text{ for some } j, \\ +\sqrt{W[i, j] / 2} & \text{if } e_k = (v_j, v_i) \text{ for some } j, \\ 0 & \text{otherwise} \end{cases}\end{split}\]

for the combinatorial Laplacian, and

\[\begin{split}D[i, k] = \begin{cases} -\sqrt{W[i, j] / 2 / d[i]} & \text{if } e_k = (v_i, v_j) \text{ for some } j, \\ +\sqrt{W[i, j] / 2 / d[i]} & \text{if } e_k = (v_j, v_i) \text{ for some } j, \\ 0 & \text{otherwise} \end{cases}\end{split}\]

for the normalized Laplacian, where $v_i \in \mathcal{V}$ is a vertex, $e_k = (v_i, v_j) \in \mathcal{E}$ is an edge from $v_i$ to $v_j$, $W[i, j]$ is the weight W of the edge $(v_i, v_j)$, $d[i]$ is the degree dw of vertex $v_i$.

For undirected graphs, only half the edges are kept (the upper triangular part of the adjacency matrix) in the interest of space and time. In that case, the $1/\sqrt{2}$ factor disappears from the above equations for $L = D D^\top$ to stand at all times.

The result is cached and accessible by the D property.

See also

grad: compute the gradient
div: compute the divergence

Examples

The difference operator is an incidence matrix. Example with a undirected graph.

>>> graph = graphs.Graph([
...     [0, 2, 0],
...     [2, 0, 1],
...     [0, 1, 0],
... ])
>>> graph.compute_laplacian('combinatorial')
>>> graph.compute_differential_operator()
>>> graph.D.toarray()
array([[-1.41421356,  0.        ],
       [ 1.41421356, -1.        ],
       [ 0.        ,  1.        ]])
>>> graph.compute_laplacian('normalized')
>>> graph.compute_differential_operator()
>>> graph.D.toarray()
array([[-1.        ,  0.        ],
       [ 0.81649658, -0.57735027],
       [ 0.        ,  1.        ]])

Example with a directed graph.

>>> graph = graphs.Graph([
...     [0, 2, 0],
...     [2, 0, 1],
...     [0, 0, 0],
... ])
>>> graph.compute_laplacian('combinatorial')
>>> graph.compute_differential_operator()
>>> graph.D.toarray()
array([[-1.        ,  1.        ,  0.        ],
       [ 1.        , -1.        , -0.70710678],
       [ 0.        ,  0.        ,  0.70710678]])
>>> graph.compute_laplacian('normalized')
>>> graph.compute_differential_operator()
>>> graph.D.toarray()
array([[-0.70710678,  0.70710678,  0.        ],
       [ 0.63245553, -0.63245553, -0.4472136 ],
       [ 0.        ,  0.        ,  1.        ]])

The graph Laplacian acts on a signal as the divergence of the gradient.

>>> G = graphs.Logo()
>>> G.compute_differential_operator()
>>> s = np.random.default_rng().normal(size=G.N)
>>> s_grad = G.D.T.dot(s)
>>> s_lap = G.D.dot(s_grad)
>>> np.linalg.norm(s_lap - G.L.dot(s)) < 1e-10
True

compute_fourier_basis(n_eigenvectors=None)¶

Compute the (partial) Fourier basis of the graph (cached).

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

Parameters

n_eigenvectorsint or None: Number of eigenvectors to compute. If None, all eigenvectors are computed. (default: None)

Notes

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

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

or a partial eigendecomposition of the graph Laplacian $L$ such that:

\[L \approx 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 n_eigenvectors $\le$ 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.coherence.

References

See [Chu97].

Examples

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

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

Compute a graph Laplacian.

For undirected graphs, the combinatorial Laplacian is defined as

\[L = D - W,\]

where $W$ is the weighted adjacency matrix and $D$ the weighted degree matrix. The normalized Laplacian is defined as

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

where $I$ is the identity matrix.

For directed graphs, the Laplacians are built from a symmetrized version of the weighted adjacency matrix that is the average of the weighted adjacency matrix and its transpose. As the Laplacian is defined as the divergence of the gradient, it is not affected by the orientation of the edges.

For both Laplacians, the diagonal entries corresponding to disconnected nodes (i.e., nodes with degree zero) are set to zero.

Once computed, the Laplacian is accessible by the attribute L.

Parameters

lap_type{‘combinatorial’, ‘normalized’}: The kind of Laplacian to compute. Default is combinatorial.

Examples

Combinatorial and normalized Laplacians of an undirected graph.

>>> graph = graphs.Graph([
...     [0, 2, 0],
...     [2, 0, 1],
...     [0, 1, 0],
... ])
>>> graph.compute_laplacian('combinatorial')
>>> graph.L.toarray()
array([[ 2., -2.,  0.],
       [-2.,  3., -1.],
       [ 0., -1.,  1.]])
>>> graph.compute_laplacian('normalized')
>>> graph.L.toarray()
array([[ 1.        , -0.81649658,  0.        ],
       [-0.81649658,  1.        , -0.57735027],
       [ 0.        , -0.57735027,  1.        ]])

Combinatorial and normalized Laplacians of a directed graph.

>>> graph = graphs.Graph([
...     [0, 2, 0],
...     [2, 0, 1],
...     [0, 0, 0],
... ])
>>> graph.compute_laplacian('combinatorial')
>>> graph.L.toarray()
array([[ 2. , -2. ,  0. ],
       [-2. ,  2.5, -0.5],
       [ 0. , -0.5,  0.5]])
>>> graph.compute_laplacian('normalized')
>>> graph.L.toarray()
array([[ 1.        , -0.89442719,  0.        ],
       [-0.89442719,  1.        , -0.4472136 ],
       [ 0.        , -0.4472136 ,  1.        ]])

The Laplacian is defined as the divergence of the gradient. See compute_differential_operator() for details.

>>> graph = graphs.Path(20)
>>> graph.compute_differential_operator()
>>> L = graph.D.dot(graph.D.T)
>>> np.all(L.toarray() == graph.L.toarray())
True

The Laplacians have a bounded spectrum.

>>> G = graphs.Sensor(50)
>>> G.compute_laplacian('combinatorial')
>>> G.compute_fourier_basis()
>>> -1e-10 < G.e[0] < 1e-10 < G.e[-1] < 2*np.max(G.dw)
True
>>> G.compute_laplacian('normalized')
>>> G.compute_fourier_basis()
>>> -1e-10 < G.e[0] < 1e-10 < G.e[-1] < 2
True

dirichlet_energy(x)[source]¶

Compute the Dirichlet energy of a signal defined on the vertices.

The Dirichlet energy of a signal $x$ is defined as

\[x^\top L x = \| \nabla_\mathcal{G} x \|_2^2 = \frac12 \sum_{i,j} W[i, j] (x[j] - x[i])^2\]

for the combinatorial Laplacian, and

\[x^\top L x = \| \nabla_\mathcal{G} x \|_2^2 = \frac12 \sum_{i,j} W[i, j] \left( \frac{x[j]}{d[j]} - \frac{x[i]}{d[i]} \right)^2\]

for the normalized Laplacian, where $d$ is the weighted degree dw, $\nabla_\mathcal{G} x = D^\top x$ and $D$ is the differential operator D. See grad() for the definition of the gradient $\nabla_\mathcal{G}$.

Parameters

xarray_like: Signal of length n_vertices living on the vertices.

Returns

energyfloat: The Dirichlet energy of the graph signal.

See also

grad: compute the gradient of a vertex signal

Examples

Non-directed graph:

>>> graph = graphs.Path(5, directed=False)
>>> signal = [0, 2, 2, 4, 4]
>>> graph.dirichlet_energy(signal)
8.0
>>> # The Dirichlet energy is indeed the squared norm of the gradient.
>>> graph.compute_differential_operator()
>>> graph.grad(signal)
array([2., 0., 2., 0.])

Directed graph:

>>> graph = graphs.Path(5, directed=True)
>>> signal = [0, 2, 2, 4, 4]
>>> graph.dirichlet_energy(signal)
4.0
>>> # The Dirichlet energy is indeed the squared norm of the gradient.
>>> graph.compute_differential_operator()
>>> graph.grad(signal)
array([1.41421356, 0.        , 1.41421356, 0.        ])

div(y)¶

Compute the divergence of a signal defined on the edges.

The divergence $z$ of a signal $y$ is defined as

\[z = \operatorname{div}_\mathcal{G} y = D y,\]

where $D$ is the differential operator D.

The value of the divergence on the vertex $v_i$ is

\[z[i] = \sum_k D[i, k] y[k] = \sum_{\{k,j | e_k=(v_j, v_i) \in \mathcal{E}\}} \sqrt{\frac{W[j, i]}{2}} y[k] - \sum_{\{k,j | e_k=(v_i, v_j) \in \mathcal{E}\}} \sqrt{\frac{W[i, j]}{2}} y[k]\]

for the combinatorial Laplacian, and

\[z[i] = \sum_k D[i, k] y[k] = \sum_{\{k,j | e_k=(v_j, v_i) \in \mathcal{E}\}} \sqrt{\frac{W[j, i]}{2 d[i]}} y[k] - \sum_{\{k,j | e_k=(v_i, v_j) \in \mathcal{E}\}} \sqrt{\frac{W[i, j]}{2 d[i]}} y[k]\]

for the normalized Laplacian.

For undirected graphs, only half the edges are kept and the $1/\sqrt{2}$ factor disappears from the above equations. See compute_differential_operator() for details.

Parameters

yarray_like: Signal of length n_edges living on the edges.

Returns

zndarray: Divergence signal of length n_vertices living on the vertices.

See also

compute_differential_operator
grad: compute the gradient of a vertex signal

Examples

Non-directed graph and combinatorial Laplacian:

>>> graph = graphs.Path(4, directed=False, lap_type='combinatorial')
>>> graph.compute_differential_operator()
>>> graph.div([2, -2, 0])
array([-2.,  4., -2.,  0.])

Directed graph and combinatorial Laplacian:

>>> graph = graphs.Path(4, directed=True, lap_type='combinatorial')
>>> graph.compute_differential_operator()
>>> graph.div([2, -2, 0])
array([-1.41421356,  2.82842712, -1.41421356,  0.        ])

Non-directed graph and normalized Laplacian:

>>> graph = graphs.Path(4, directed=False, lap_type='normalized')
>>> graph.compute_differential_operator()
>>> graph.div([2, -2, 0])
array([-2.        ,  2.82842712, -1.41421356,  0.        ])

Directed graph and normalized Laplacian:

>>> graph = graphs.Path(4, directed=True, lap_type='normalized')
>>> graph.compute_differential_operator()
>>> graph.div([2, -2, 0])
array([-2.        ,  2.82842712, -1.41421356,  0.        ])

estimate_lmax(method='lanczos')[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

method{‘lanczos’, ‘bounds’}: Whether to estimate the largest eigenvalue with the implicitly restarted Lanczos method, or to return an upper bound on the spectrum of the Laplacian.

Notes

Runs the implicitly restarted Lanczos method (as implemented in scipy.sparse.linalg.eigsh()) with a large tolerance, then increases the calculated largest eigenvalue by 1 percent. For much of the PyGSP machinery, we need to approximate filter 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.

A faster but less tight alternative is to use known algebraic bounds on the graph Laplacian.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()  # True value.
>>> print('{:.2f}'.format(G.lmax))
13.78
>>> G.estimate_lmax(method='lanczos')  # Estimate.
>>> print('{:.2f}'.format(G.lmax))
13.92
>>> G.estimate_lmax(method='bounds')  # Upper bound.
>>> print('{:.2f}'.format(G.lmax))
18.58

extract_components()[source]¶

Split the graph into connected components.

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

Returns

graphslist: 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)
>>> components = G.extract_components()
>>> has_sinks = 'sink' in components[0].info
>>> sinks_0 = components[0].info['sink'] if has_sinks else []

classmethod from_graphtool(graph, weight='weight')¶

Import a graph from graph-tool.

Edge weights are retrieved as an edge property, under the name specified by the weight parameter.

Signals are retrieved from node properties, and stored in the signals dictionary under the property name. N-dimensional signals that were broken during export are joined.

Parameters

graphgraph_tool.Graph: A graph-tool graph object.
weightstring: The edge property that holds the numerical values used as the edge weights. All edge weights are set to 1 if None, or not found.

Returns

graphGraph: A PyGSP graph object.

See also

from_networkx: import from NetworkX
load: load from a file

Notes

If the graph has multiple edge connecting the same two nodes, a sum over the edges is taken to merge them.

Examples

>>> import graph_tool as gt
>>> graph = gt.Graph(directed=False)
>>> e1 = graph.add_edge(0, 1)
>>> e2 = graph.add_edge(1, 2)
>>> v = graph.add_vertex()
>>> eprop = graph.new_edge_property("double")
>>> eprop[e1] = 0.2
>>> eprop[graph.edge(1, 2)] = 0.9
>>> graph.edge_properties["weight"] = eprop
>>> vprop = graph.new_vertex_property("double", val=np.nan)
>>> vprop[3] = 3.1416
>>> graph.vertex_properties["sig"] = vprop
>>> graph = graphs.Graph.from_graphtool(graph)
>>> graph.W.toarray()
array([[0. , 0.2, 0. , 0. ],
       [0.2, 0. , 0.9, 0. ],
       [0. , 0.9, 0. , 0. ],
       [0. , 0. , 0. , 0. ]])
>>> graph.signals
{'sig': PropertyArray([   nan,    nan,    nan, 3.1416])}

classmethod from_networkx(graph, weight='weight')¶

Import a graph from NetworkX.

Edge weights are retrieved as an edge attribute, under the name specified by the weight parameter.

Signals are retrieved from node attributes, and stored in the signals dictionary under the attribute name. N-dimensional signals that were broken during export are joined.

Parameters

graphnetworkx.Graph: A NetworkX graph object.
weightstring or None, optional: The edge attribute that holds the numerical values used as the edge weights. All edge weights are set to 1 if None, or not found.

Returns

graphGraph: A PyGSP graph object.

See also

from_graphtool: import from graph-tool
load: load from a file

Notes

The nodes are ordered according to networkx.Graph.nodes().

In NetworkX, node attributes need not be set for every node. If a node attribute is not set for a node, a NaN is assigned to the corresponding signal for that node.

If the graph is a networkx.MultiGraph, multiedges are aggregated by summation.

Examples

>>> import networkx as nx
>>> graph = nx.Graph()
>>> graph.add_edge(1, 2, weight=0.2)
>>> graph.add_edge(2, 3, weight=0.9)
>>> graph.add_node(4, sig=3.1416)
>>> graph.nodes()
NodeView((1, 2, 3, 4))
>>> graph = graphs.Graph.from_networkx(graph)
>>> graph.W.toarray()
array([[0. , 0.2, 0. , 0. ],
       [0.2, 0. , 0.9, 0. ],
       [0. , 0.9, 0. , 0. ],
       [0. , 0. , 0. , 0. ]])
>>> graph.signals
{'sig': array([   nan,    nan,    nan, 3.1416])}

get_edge_list()[source]¶

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

Each edge $e_k = (v_i, v_j) \in \mathcal{E}$ from $v_i$ to $v_j$ is associated with the weight $W[i, j]$. For each edge $e_k$, the method returns $(i, j, W[i, j])$ as (sources[k], targets[k], weights[k]), with $i \in [0, |\mathcal{V}|-1], j \in [0, |\mathcal{V}|-1], k \in [0, |\mathcal{E}|-1]$.

Returns

sourcesvector of int: Source node indices.
targetsvector of int: Target node indices.
weightsvector of float: Edge weights.

Notes

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.

Edge orientation (i.e., which node is the source or the target) is arbitrary for undirected graphs. The implementation uses the upper triangular part of the adjacency matrix, hence $i \leq j \ \forall k$.

Examples

Edge list of a directed graph.

>>> graph = graphs.Graph([
...     [0, 3, 0],
...     [3, 0, 4],
...     [0, 0, 0],
... ])
>>> sources, targets, weights = graph.get_edge_list()
>>> list(sources), list(targets), list(weights)
([0, 1, 1], [1, 0, 2], [3, 3, 4])

Edge list of an undirected graph.

>>> graph = graphs.Graph([
...     [0, 3, 0],
...     [3, 0, 4],
...     [0, 4, 0],
... ])
>>> sources, targets, weights = graph.get_edge_list()
>>> list(sources), list(targets), list(weights)
([0, 1], [1, 2], [3, 4])

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

sarray_like: Graph signal in the vertex domain.

Returns

s_hatndarray: Representation of s in the Fourier domain.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> s = np.random.default_rng().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

grad(x)¶

Compute the gradient of a signal defined on the vertices.

The gradient $y$ of a signal $x$ is defined as

\[y = \nabla_\mathcal{G} x = D^\top x,\]

where $D$ is the differential operator D.

The value of the gradient on the edge $e_k = (v_i, v_j)$ from $v_i$ to $v_j$ with weight $W[i, j]$ is

\[y[k] = D[i, k] x[i] + D[j, k] x[j] = \sqrt{\frac{W[i, j]}{2}} (x[j] - x[i])\]

for the combinatorial Laplacian, and

\[y[k] = \sqrt{\frac{W[i, j]}{2}} \left( \frac{x[j]}{\sqrt{d[j]}} - \frac{x[i]}{\sqrt{d[i]}} \right)\]

for the normalized Laplacian.

For undirected graphs, only half the edges are kept and the $1/\sqrt{2}$ factor disappears from the above equations. See compute_differential_operator() for details.

Parameters

xarray_like: Signal of length n_vertices living on the vertices.

Returns

yndarray: Gradient signal of length n_edges living on the edges.

See also

compute_differential_operator
div: compute the divergence of an edge signal
dirichlet_energy: compute the norm of the gradient

Examples

Non-directed graph and combinatorial Laplacian:

>>> graph = graphs.Path(4, directed=False, lap_type='combinatorial')
>>> graph.compute_differential_operator()
>>> graph.grad([0, 2, 4, 2])
array([ 2.,  2., -2.])

Directed graph and combinatorial Laplacian:

>>> graph = graphs.Path(4, directed=True, lap_type='combinatorial')
>>> graph.compute_differential_operator()
>>> graph.grad([0, 2, 4, 2])
array([ 1.41421356,  1.41421356, -1.41421356])

Non-directed graph and normalized Laplacian:

>>> graph = graphs.Path(4, directed=False, lap_type='normalized')
>>> graph.compute_differential_operator()
>>> graph.grad([0, 2, 4, 2])
array([ 1.41421356,  1.41421356, -0.82842712])

Directed graph and normalized Laplacian:

>>> graph = graphs.Path(4, directed=True, lap_type='normalized')
>>> graph.compute_differential_operator()
>>> graph.grad([0, 2, 4, 2])
array([ 1.41421356,  1.41421356, -0.82842712])

has_loops()[source]¶

Check if any vertex is connected to itself.

A graph has self-loops if and only if the diagonal entries of its adjacency matrix are not all zero.

Returns

loopsbool: True if the graph has self-loops, False otherwise.

Examples

Without self-loops:

>>> graph = graphs.Graph([
...     [0, 3, 0],
...     [3, 0, 4],
...     [0, 0, 0],
... ])
>>> graph.has_loops()
False

With a self-loop:

>>> graph = graphs.Graph([
...     [1, 3, 0],
...     [3, 0, 4],
...     [0, 0, 0],
... ])
>>> graph.has_loops()
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_hatarray_like: Graph signal in the Fourier domain.

Returns

sndarray: Representation of s_hat in the vertex domain.

Examples

>>> G = graphs.Logo()
>>> G.compute_fourier_basis()
>>> s_hat = np.random.default_rng().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()[source]¶

Check if the graph is connected (cached).

A graph is connected if and only if there exists a (directed) path between any two vertices.

Returns

connectedbool: True if the graph is connected, False otherwise.

Notes

For undirected graphs, starting at a vertex and trying to visit all the others is enough. For directed graphs, one needs to check that a vertex can both be visited by all the others and visit all the others.

Examples

Connected graph:

>>> graph = graphs.Graph([
...     [0, 3, 0, 0],
...     [3, 0, 4, 0],
...     [0, 4, 0, 2],
...     [0, 0, 2, 0],
... ])
>>> graph.is_connected()
True

Disconnected graph:

>>> graph = graphs.Graph([
...     [0, 3, 0, 0],
...     [3, 0, 4, 0],
...     [0, 0, 0, 2],
...     [0, 0, 2, 0],
... ])
>>> graph.is_connected()
False

is_directed()[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 not symmetric.

Returns

directedbool: True if the graph is directed, False otherwise.

Examples

Directed graph:

>>> graph = graphs.Graph([
...     [0, 3, 0],
...     [3, 0, 4],
...     [0, 0, 0],
... ])
>>> graph.is_directed()
True

Undirected graph:

>>> graph = graphs.Graph([
...     [0, 3, 0],
...     [3, 0, 4],
...     [0, 4, 0],
... ])
>>> graph.is_directed()
False

is_weighted()[source]¶

Check if the graph is weighted.

A graph is unweighted (binary) if and only if all the entries in the adjacency matrix are either zero or one.

Returns

weightedbool: True if the graph is weighted, False otherwise.

Examples

Unweighted (binary) graph:

>>> graph = graphs.Graph([
...     [0, 1, 0],
...     [1, 0, 1],
...     [0, 1, 0],
... ])
>>> graph.is_weighted()
False

Weighted graph:

>>> graph = graphs.Graph([
...     [0, 2, 0],
...     [2, 0, 1],
...     [0, 1, 0],
... ])
>>> graph.is_weighted()
True

classmethod load(path, fmt=None, backend=None)¶

Load a graph from a file.

Edge weights are retrieved as an edge attribute named “weight”.

Signals are retrieved from node attributes, and stored in the signals dictionary under the attribute name. N-dimensional signals that were broken during export are joined.

Parameters

pathstring: Path to the file from which to load the graph.
fmt{‘graphml’, ‘gml’, ‘gexf’, None}, optional: Format in which the graph is saved. Guessed from the filename extension if None.
backend{‘networkx’, ‘graph-tool’, None}, optional: Library used to load the graph. Automatically chosen if None.

Returns

graphGraph: The loaded graph.

See also

save: save a graph to a file
from_networkx: load with NetworkX then import in the PyGSP
from_graphtool: load with graph-tool then import in the PyGSP

Notes

A lossless round-trip is only guaranteed if the graph (and its signals) is saved and loaded with the same backend.

Loading from other formats is possible by loading in NetworkX or graph-tool, and importing to the PyGSP. The proposed formats are however tested for faithful round-trips.

Examples

>>> graph = graphs.Logo()
>>> graph.save('logo.graphml')
>>> graph = graphs.Graph.load('logo.graphml')
>>> import os
>>> os.remove('logo.graphml')

plot(vertex_color=None, vertex_size=None, highlight=[], edges=None, edge_color=None, edge_width=None, indices=False, colorbar=True, limits=None, ax=None, title=None, backend=None)[source]¶

Plot a graph with signals as color or vertex size.

Parameters

vertex_colorarray_like or color: Signal to plot as vertex color (length is the number of vertices). If None, vertex color is set to graph.plotting[‘vertex_color’]. Alternatively, a color can be set in any format accepted by matplotlib. Each vertex color can by specified by an RGB(A) array of dimension n_vertices x 3 (or 4).
vertex_sizearray_like or int: Signal to plot as vertex size (length is the number of vertices). Vertex size ranges from 0.5 to 2 times graph.plotting[‘vertex_size’]. If None, vertex size is set to graph.plotting[‘vertex_size’]. Alternatively, a size can be passed as an integer. The pyqtgraph backend only accepts an integer size.
highlightiterable: List of indices of vertices to be highlighted. Useful for example to show where a filter was localized. Only available with the matplotlib backend.
edgesbool: Whether to draw edges in addition to vertices. Default to True if less than 10,000 edges to draw. Note that drawing many edges can be slow.
edge_colorarray_like or color: Signal to plot as edge color (length is the number of edges). Edge color is given by graph.plotting[‘edge_color’] and transparency ranges from 0.2 to 0.9. If None, edge color is set to graph.plotting[‘edge_color’]. Alternatively, a color can be set in any format accepted by matplotlib. Each edge color can by specified by an RGB(A) array of dimension n_edges x 3 (or 4). Only available with the matplotlib backend.
edge_widtharray_like or int: Signal to plot as edge width (length is the number of edges). Edge width ranges from 0.5 to 2 times graph.plotting[‘edge_width’]. If None, edge width is set to graph.plotting[‘edge_width’]. Alternatively, a width can be passed as an integer. Only available with the matplotlib backend.
indicesbool: Whether to print the node indices (in the adjacency / Laplacian matrix and signal vectors) on top of each node. Useful to locate a node of interest. Only available with the matplotlib backend.
colorbarbool: Whether to plot a colorbar indicating the signal’s amplitude. Only available with the matplotlib backend.
limits[vmin, vmax]: Map colors from vmin to vmax. Defaults to signal minimum and maximum value. Only available with the matplotlib backend.
axmatplotlib.axes.Axes: Axes where to draw the graph. Optional, created if not passed. Only available with the matplotlib backend.
titlestr: Title of the figure.
backend: {‘matplotlib’, ‘pyqtgraph’, None}: Defines the drawing backend to use. Defaults to pygsp.plotting.BACKEND.

Returns

figmatplotlib.figure.Figure: The figure the plot belongs to. Only with the matplotlib backend.
axmatplotlib.axes.Axes: The axes the plot belongs to. Only with the matplotlib backend.

Notes

The orientation of directed edges is not shown. If edges exist in both directions, they will be drawn on top of each other.

Examples

>>> import matplotlib
>>> graph = graphs.Sensor(20, seed=42)
>>> graph.compute_fourier_basis(n_eigenvectors=4)
>>> _, _, weights = graph.get_edge_list()
>>> fig, ax = graph.plot(graph.U[:, 1], vertex_size=graph.dw,
...                      edge_color=weights)
>>> graph.plotting['vertex_size'] = 300
>>> graph.plotting['edge_width'] = 5
>>> graph.plotting['edge_style'] = '--'
>>> fig, ax = graph.plot(edge_width=weights, edge_color=(0, .8, .8, .5),
...                      vertex_color='black')
>>> fig, ax = graph.plot(vertex_size=graph.dw, indices=True,
...                      highlight=[17, 3, 16], edges=False)

plot_signal(*args, **kwargs)[source]¶: Deprecated, use plot() instead.

plot_spectrogram(node_idx=None)[source]¶

Plot the graph’s spectrogram.

Parameters

node_idxndarray: Order to sort the nodes in the spectrogram. By default, does not reorder the nodes.

Notes

This function is only implemented for the pyqtgraph backend at the moment.

Examples

>>> G = graphs.Ring(15)
>>> G.plot_spectrogram()

save(path, fmt=None, backend=None)¶

Save the graph to a file.

Edge weights are stored as an edge attribute, under the name “weight”.

Signals are stored as node attributes, under their name in the signals dictionary. N-dimensional signals are broken into N 1-dimensional signals. They will eventually be joined back together on import.

Supported formats are:

GraphML, a comprehensive XML format. Supported by NetworkX, graph-tool, NetworKit, igraph, Gephi, Cytoscape, SocNetV.
GML (Graph Modelling Language), a simple non-XML format. Supported by NetworkX, graph-tool, NetworKit, igraph, Gephi, Cytoscape, SocNetV, Tulip.
GEXF (Graph Exchange XML Format), Gephi’s XML format. Supported by NetworkX, NetworKit, Gephi, Tulip, ngraph.

If unsure, we recommend GraphML.

Parameters

pathstring: Path to the file where the graph is to be saved.
fmt{‘graphml’, ‘gml’, ‘gexf’, None}, optional: Format in which to save the graph. Guessed from the filename extension if None.
backend{‘networkx’, ‘graph-tool’, None}, optional: Library used to load the graph. Automatically chosen if None.

See also

load: load a graph from a file
to_networkx: export as a NetworkX graph, and save with NetworkX
to_graphtool: export as a graph-tool graph, and save with graph-tool

Notes

A lossless round-trip is only guaranteed if the graph (and its signals) is saved and loaded with the same backend.

Saving in other formats is possible by exporting to NetworkX or graph-tool, and using their respective saving functionality. The proposed formats are however tested for faithful round-trips.

Edge weights and signal values are rounded at the sixth decimal when saving in fmt='gml' with backend='graph-tool'.

Examples

>>> graph = graphs.Logo()
>>> graph.save('logo.graphml')
>>> graph = graphs.Graph.load('logo.graphml')
>>> import os
>>> os.remove('logo.graphml')

set_coordinates(kind='spring', seed=None, **kwargs)¶

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

Parameters

kindstring 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, laplacian_eigenmap2D, laplacian_eigenmap3D. Default is ‘spring’.
seedint: Seed for the random number generator when kind is ‘random’, ‘community’, or ‘spring’.
kwargsdict: Additional parameters to be passed to the Fruchterman-Reingold force-directed algorithm when kind is spring.

Examples

>>> G = graphs.ErdosRenyi()
>>> G.set_coordinates()
>>> fig, ax = G.plot()

set_signal(signal, name)[source]¶

Attach a signal to the graph.

Attached signals can be accessed (and modified or deleted) through the signals dictionary.

Parameters

signalarray_like: A sequence that assigns a value to each vertex. The value of the signal at vertex i is signal[i].
nameString: Name of the signal used as a key in the signals dictionary.

Examples

>>> graph = graphs.Sensor(10)
>>> signal = np.arange(graph.n_vertices)
>>> graph.set_signal(signal, 'mysignal')
>>> graph.signals
{'mysignal': array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])}

subgraph(vertices)[source]¶

Create a subgraph from a list of vertices.

Parameters

verticeslist: Vertices to keep. Either a list of indices or an indicator function.

Returns

subgraphGraph: Subgraph.

Examples

>>> graph = graphs.Graph([
...     [0., 3., 0., 0.],
...     [3., 0., 4., 0.],
...     [0., 4., 0., 2.],
...     [0., 0., 2., 0.],
... ])
>>> graph = graph.subgraph([0, 2, 1])
>>> graph.W.toarray()
array([[0., 0., 3.],
       [0., 0., 4.],
       [3., 4., 0.]])

to_graphtool()¶

Export the graph to graph-tool.

Edge weights are stored as an edge property map, under the name “weight”.

Signals are stored as vertex property maps, under their name in the signals dictionary. N-dimensional signals are broken into N 1-dimensional signals. They will eventually be joined back together on import.

Returns

graphgraph_tool.Graph: A graph-tool graph object.

See also

to_networkx: export to NetworkX
save: save to a file

Examples

>>> import graph_tool as gt
>>> import graph_tool.draw
>>> from matplotlib import pyplot as plt
>>> graph = graphs.Path(4, directed=True)
>>> graph.set_signal(np.full(4, 2.3), 'signal')
>>> graph = graph.to_graphtool()
>>> graph.is_directed()
True
>>> graph.vertex_properties['signal'][2]
2.3
>>> graph.edge_properties['weight'][graph.edge(0, 1)]
1.0
>>> # gt.draw.graph_draw(graph, vertex_text=graph.vertex_index)

Another common goal is to use graph-tool to compute some properties to be imported back in the PyGSP as signals.

>>> import graph_tool as gt
>>> import graph_tool.centrality
>>> from matplotlib import pyplot as plt
>>> graph = graphs.Sensor(100, seed=42)
>>> graph.set_signal(graph.coords, 'coords')
>>> graph = graph.to_graphtool()
>>> vprop, eprop = gt.centrality.betweenness(
...     graph, weight=graph.edge_properties['weight'])
>>> graph.vertex_properties['betweenness'] = vprop
>>> graph = graphs.Graph.from_graphtool(graph)
>>> graph.compute_fourier_basis()
>>> graph.set_coordinates(graph.signals['coords'])
>>> fig, axes = plt.subplots(1, 2)
>>> _ = graph.plot(graph.signals['betweenness'], ax=axes[0])
>>> _ = axes[1].plot(graph.e, graph.gft(graph.signals['betweenness']))

to_networkx()¶

Export the graph to NetworkX.

Edge weights are stored as an edge attribute, under the name “weight”.

Signals are stored as node attributes, under their name in the signals dictionary. N-dimensional signals are broken into N 1-dimensional signals. They will eventually be joined back together on import.

Returns

graphnetworkx.Graph: A NetworkX graph object.

See also

to_graphtool: export to graph-tool
save: save to a file

Examples

>>> import networkx as nx
>>> from matplotlib import pyplot as plt
>>> graph = graphs.Path(4, directed=True)
>>> graph.set_signal(np.full(4, 2.3), 'signal')
>>> graph = graph.to_networkx()
>>> print(nx.info(graph))
DiGraph named 'Path' with 4 nodes and 3 edges
>>> nx.is_directed(graph)
True
>>> graph.nodes()
NodeView((0, 1, 2, 3))
>>> graph.edges()
OutEdgeView([(0, 1), (1, 2), (2, 3)])
>>> graph.nodes()[2]
{'signal': 2.3}
>>> graph.edges()[(0, 1)]
{'weight': 1.0}
>>> # nx.draw(graph, with_labels=True)

Another common goal is to use NetworkX to compute some properties to be be imported back in the PyGSP as signals.

>>> import networkx as nx
>>> from matplotlib import pyplot as plt
>>> graph = graphs.Sensor(100, seed=42)
>>> graph.set_signal(graph.coords, 'coords')
>>> graph = graph.to_networkx()
>>> betweenness = nx.betweenness_centrality(graph, weight='weight')
>>> nx.set_node_attributes(graph, betweenness, 'betweenness')
>>> graph = graphs.Graph.from_networkx(graph)
>>> graph.compute_fourier_basis()
>>> graph.set_coordinates(graph.signals['coords'])
>>> fig, axes = plt.subplots(1, 2)
>>> _ = graph.plot(graph.signals['betweenness'], ax=axes[0])
>>> _ = axes[1].plot(graph.e, graph.gft(graph.signals['betweenness']))

property 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$.

property D¶

Differential operator (for gradient and divergence).

Is computed by compute_differential_operator().

property U¶

Fourier basis (eigenvectors of the Laplacian).

Is computed by compute_fourier_basis().

property W¶: Weighted adjacency matrix of the graph.

property coherence¶

Coherence of the Fourier basis.

The mutual coherence between the basis of Kronecker deltas and the basis formed by the eigenvectors of the graph Laplacian is defined as

\[\mu = \max_{\ell,i} \langle U_\ell, \delta_i \rangle = \max_{\ell,i} | U_{\ell, i} | \in \left[ \frac{1}{\sqrt{N}}, 1 \right],\]

where $N$ is the number of vertices, $\delta_i \in \mathbb{R}^N$ denotes the Kronecker delta that is non-zero on vertex $v_i$, and $U_\ell \in \mathbb{R}^N$ denotes the $\ell^\text{th}$ eigenvector of the graph Laplacian (i.e., the $\ell^\text{th}$ Fourier mode).

The coherence is a measure of the localization of the Fourier modes (Laplacian eigenvectors). The larger the value, the more localized the eigenvectors can be. The extreme is a node that is disconnected from the rest of the graph: an eigenvector will be localized as a Kronecker delta there (i.e., $\mu = 1$). In the classical setting, Fourier modes (which are complex exponentials) are completely delocalized, and the coherence is minimal.

The value is computed by compute_fourier_basis().

Examples

Delocalized eigenvectors.

>>> graph = graphs.Path(100)
>>> graph.compute_fourier_basis()
>>> minimum = 1 / np.sqrt(graph.n_vertices)
>>> print('{:.2f} in [{:.2f}, 1]'.format(graph.coherence, minimum))
0.14 in [0.10, 1]
>>>
>>> # Plot some delocalized eigenvectors.
>>> import matplotlib.pyplot as plt
>>> graph.set_coordinates('line1D')
>>> _ = graph.plot(graph.U[:, :5])

Localized eigenvectors.

>>> graph = graphs.Sensor(64, seed=10)
>>> graph.compute_fourier_basis()
>>> minimum = 1 / np.sqrt(graph.n_vertices)
>>> print('{:.2f} in [{:.2f}, 1]'.format(graph.coherence, minimum))
0.84 in [0.12, 1]
>>>
>>> # Plot the most localized eigenvector.
>>> import matplotlib.pyplot as plt
>>> idx = np.argmax(np.abs(graph.U))
>>> idx_vertex, idx_fourier = np.unravel_index(idx, graph.U.shape)
>>> _ = graph.plot(graph.U[:, idx_fourier], highlight=idx_vertex)

property d¶

The degree (number of neighbors) of vertices.

For undirected graphs, the degree of a vertex is the number of vertices it is connected to. For directed graphs, the degree is the average of the in and out degrees, where the in degree is the number of incoming edges, and the out degree the number of outgoing edges.

In both cases, the degree of the vertex $v_i$ is the average between the number of non-zero values in the $i$-th column (the in degree) and the $i$-th row (the out degree) of the weighted adjacency matrix W.

Examples

Undirected graph:

>>> graph = graphs.Graph([
...     [0, 1, 0],
...     [1, 0, 2],
...     [0, 2, 0],
... ])
>>> print(graph.d)  # Number of neighbors.
[1 2 1]
>>> print(graph.dw)  # Weighted degree.
[1 3 2]

Directed graph:

>>> graph = graphs.Graph([
...     [0, 1, 0],
...     [0, 0, 2],
...     [0, 2, 0],
... ])
>>> print(graph.d)  # Number of neighbors.
[0.5 1.5 1. ]
>>> print(graph.dw)  # Weighted degree.
[0.5 2.5 2. ]

property dw¶

The weighted degree of vertices.

For undirected graphs, the weighted degree of the vertex $v_i$ is defined as

\[d[i] = \sum_j W[j, i] = \sum_j W[i, j],\]

where $W$ is the weighted adjacency matrix W.

For directed graphs, the weighted degree of the vertex $v_i$ is defined as

\[d[i] = \frac12 (d^\text{in}[i] + d^\text{out}[i]) = \frac12 (\sum_j W[j, i] + \sum_j W[i, j]),\]

i.e., as the average of the in and out degrees.

Examples

Undirected graph:

>>> graph = graphs.Graph([
...     [0, 1, 0],
...     [1, 0, 2],
...     [0, 2, 0],
... ])
>>> print(graph.d)  # Number of neighbors.
[1 2 1]
>>> print(graph.dw)  # Weighted degree.
[1 3 2]

Directed graph:

>>> graph = graphs.Graph([
...     [0, 1, 0],
...     [0, 0, 2],
...     [0, 2, 0],
... ])
>>> print(graph.d)  # Number of neighbors.
[0.5 1.5 1. ]
>>> print(graph.dw)  # Weighted degree.
[0.5 2.5 2. ]

property e¶

Eigenvalues of the Laplacian (square of graph frequencies).

Is computed by compute_fourier_basis().

property lmax¶

Largest eigenvalue of the graph Laplacian.

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

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.complement`([frame_bound])	Return the filter that makes the frame tight.
`Filter.inverse`()	Return the pseudo-inverse filter bank.
`Filter.compute_frame`(**kwargs)	Compute the associated frame.
`Filter.estimate_frame_bounds`([x])	Estimate lower and upper frame bounds.
`Filter.plot`([n, eigenvalues, sum, labels, …])	Plot the spectral response of a filter bank.
`Filter.localize`(i, **kwargs)	Localize the kernels at a node (to visualize them).

Filters¶

Then, derived classes implement various common graph filters.

Filters that solve differential equations

The following filters solve partial differential equations (PDEs) on graphs, which model processes such as heat diffusion or wave propagation.

`Heat`(G[, scale, normalize])	Design a filter bank of heat kernels.
`Wave`(G[, time, speed])	Design a filter bank of wave kernels.

Low-pass filters

Heat(G[, scale, normalize])

Design a filter bank of heat kernels.

Band-pass filters

These filters can be configured to be low-pass, high-pass, or band-pass.

`Expwin`(G[, band_min, band_max, slope])	Design an exponential window filter.
`Rectangular`(G[, band_min, band_max])	Design a rectangular filter.

Filter banks of two filters: a low-pass and a high-pass

`Regular`(G[, degree])	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).

Filter banks composed of dilated or translated filters

`Abspline`(G[, Nf, lpfactor, scales])	Design an A B cubic spline wavelet 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 for vertex-frequency analyzes

Those filter banks are composed of shifted versions of a mother filter, one per graph frequency (Laplacian eigenvalue). They can analyze frequency content locally, as a windowed graph Fourier transform.

`Gabor`(graph, kernel)	Design a filter bank with a kernel centered at each frequency.
`Modulation`(graph, kernel[, modulation_first])	Design a filter bank with a kernel centered at each frequency.

Approximations¶

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

Chebyshev polynomials

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

Lanczos algorithm

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

Plotting¶

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

Most users won’t use this module directly. Graphs (from pygsp.graphs) are to be plotted with pygsp.graphs.Graph.plot() and pygsp.graphs.Graph.plot_spectrogram(). Filters (from pygsp.filters) are to be plotted with pygsp.filters.Filter.plot().

pygsp.plotting.BACKEND¶: The default 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 last created figure, alias to plt.close().

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

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

Show created figures, alias to plt.show().

By default, showing plots does not block the prompt. Calling this function will block execution.

Learning¶

The pygsp.learning module provides functions to solve learning problems.

Semi-supervized learning¶

Those functions help to solve a semi-supervized learning problem, i.e., a problem where only some values of a graph signal are known and the others shall be inferred.

`regression_tikhonov`(G, y, M[, tau])	Solve a regression problem on graph via Tikhonov minimization.
`classification_tikhonov`(G, y, M[, tau])	Solve a classification problem on graph via Tikhonov minimization.
`classification_tikhonov_simplex`(G, y, M[, tau])	Solve a classification problem on graph via Tikhonov minimization with simple constraints.

pygsp.learning.classification_tikhonov(G, y, M, tau=0)[source]¶

Solve a classification problem on graph via Tikhonov minimization.

The function first transforms $y$ in logits $Y$, then solves

\[\operatorname*{arg min}_X \| M X - Y \|_2^2 + \tau \ tr(X^T L X)\]

if $\tau > 0$, and

\[\operatorname*{arg min}_X tr(X^T L X) \ \text{ s. t. } \ Y = M X\]

otherwise, where $X$ and $Y$ are logits. The function returns the maximum of the logits.

Parameters

Gpygsp.graphs.Graph
yarray, length G.n_vertices: Measurements.
Marray of boolean, length G.n_vertices: Masking vector.
taufloat: Regularization parameter.

Returns

logitsarray, length G.n_vertices: The logits $X$.

Examples

>>> from pygsp import graphs, learning
>>> import matplotlib.pyplot as plt
>>>
>>> G = graphs.Logo()

Create a ground truth signal:

>>> signal = np.zeros(G.n_vertices)
>>> signal[G.info['idx_s']] = 1
>>> signal[G.info['idx_p']] = 2

Construct a measurement signal from a binary mask:

>>> rng = np.random.default_rng(42)
>>> mask = rng.uniform(0, 1, G.n_vertices) > 0.5
>>> measures = signal.copy()
>>> measures[~mask] = np.nan

Solve the classification problem by reconstructing the signal:

>>> recovery = learning.classification_tikhonov(G, measures, mask, tau=0)

Plot the results. Note that we recover the class with np.argmax(recovery, axis=1).

>>> prediction = np.argmax(recovery, axis=1)
>>> fig, ax = plt.subplots(2, 3, sharey=True, figsize=(10, 6))
>>> _ = G.plot(signal, ax=ax[0, 0], title='Ground truth')
>>> _ = G.plot(measures, ax=ax[0, 1], title='Measurements')
>>> _ = G.plot(prediction, ax=ax[0, 2], title='Recovered class')
>>> _ = G.plot(recovery[:, 0], ax=ax[1, 0], title='Logit 0')
>>> _ = G.plot(recovery[:, 1], ax=ax[1, 1], title='Logit 1')
>>> _ = G.plot(recovery[:, 2], ax=ax[1, 2], title='Logit 2')
>>> _ = fig.tight_layout()

pygsp.learning.classification_tikhonov_simplex(G, y, M, tau=0.1, **kwargs)[source]¶

Solve a classification problem on graph via Tikhonov minimization with simple constraints.

The function first transforms $y$ in logits $Y$, then solves

\[\operatorname*{arg min}_X \| M X - Y \|_2^2 + \tau \ tr(X^T L X) \text{ s.t. } sum(X) = 1 \text{ and } X >= 0,\]

where $X$ and $Y$ are logits.

Parameters

Gpygsp.graphs.Graph
yarray, length G.n_vertices: Measurements.
Marray of boolean, length G.n_vertices: Masking vector.
taufloat: Regularization parameter.
kwargsdict: Parameters for pyunlocbox.solvers.solve().

Returns

logitsarray, length G.n_vertices: The logits $X$.

Examples

>>> from pygsp import graphs, learning
>>> import matplotlib.pyplot as plt
>>>
>>> G = graphs.Logo()
>>> G.estimate_lmax()

Create a ground truth signal:

>>> signal = np.zeros(G.n_vertices)
>>> signal[G.info['idx_s']] = 1
>>> signal[G.info['idx_p']] = 2

Construct a measurement signal from a binary mask:

>>> rng = np.random.default_rng(42)
>>> mask = rng.uniform(0, 1, G.n_vertices) > 0.5
>>> measures = signal.copy()
>>> measures[~mask] = np.nan

Solve the classification problem by reconstructing the signal:

>>> recovery = learning.classification_tikhonov_simplex(
...     G, measures, mask, tau=0.1, verbosity='NONE')

Plot the results. Note that we recover the class with np.argmax(recovery, axis=1).

>>> prediction = np.argmax(recovery, axis=1)
>>> fig, ax = plt.subplots(2, 3, sharey=True, figsize=(10, 6))
>>> _ = G.plot(signal, ax=ax[0, 0], title='Ground truth')
>>> _ = G.plot(measures, ax=ax[0, 1], title='Measurements')
>>> _ = G.plot(prediction, ax=ax[0, 2], title='Recovered class')
>>> _ = G.plot(recovery[:, 0], ax=ax[1, 0], title='Logit 0')
>>> _ = G.plot(recovery[:, 1], ax=ax[1, 1], title='Logit 1')
>>> _ = G.plot(recovery[:, 2], ax=ax[1, 2], title='Logit 2')
>>> _ = fig.tight_layout()

pygsp.learning.regression_tikhonov(G, y, M, tau=0)[source]¶

Solve a regression problem on graph via Tikhonov minimization.

The function solves

\[\operatorname*{arg min}_x \| M x - y \|_2^2 + \tau \ x^T L x\]

if $\tau > 0$, and

\[\operatorname*{arg min}_x x^T L x \ \text{ s. t. } \ y = M x\]

otherwise.

Parameters

Gpygsp.graphs.Graph
yarray, length G.n_vertices: Measurements.
Marray of boolean, length G.n_vertices: Masking vector.
taufloat: Regularization parameter.

Returns

xarray, length G.n_vertices: Recovered values $x$.

Examples

>>> from pygsp import graphs, filters, learning
>>> import matplotlib.pyplot as plt
>>>
>>> G = graphs.Sensor(N=100, seed=42)
>>> G.estimate_lmax()

Create a smooth ground truth signal:

>>> filt = lambda x: 1 / (1 + 10*x)
>>> filt = filters.Filter(G, filt)
>>> rng = np.random.default_rng(42)
>>> signal = filt.analyze(rng.normal(size=G.n_vertices))

Construct a measurement signal from a binary mask:

>>> mask = rng.uniform(0, 1, G.n_vertices) > 0.5
>>> measures = signal.copy()
>>> measures[~mask] = np.nan

Solve the regression problem by reconstructing the signal:

>>> recovery = learning.regression_tikhonov(G, measures, mask, tau=0)

Plot the results:

>>> fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(10, 3))
>>> limits = [signal.min(), signal.max()]
>>> _ = G.plot(signal, ax=ax1, limits=limits, title='Ground truth')
>>> _ = G.plot(measures, ax=ax2, limits=limits, title='Measures')
>>> _ = G.plot(recovery, ax=ax3, limits=limits, title='Recovery')
>>> _ = fig.tight_layout()

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, seed])	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

GGraph structure: The graph to reduce.
levelsint: Number of level of decomposition
lambdfloat: Stability parameter. It adds self loop to the graph to give the algorithm some stability (default = 0.025). [UNUSED?!]
sparsifybool: To perform a spectral sparsification step immediately after the graph reduction (default is True).
sparsify_epsfloat: 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_methodstring: The graph reduction method (default is ‘kron’)
compute_full_eigenbool: To also compute the graph Laplacian eigenvalues and eigenvectors for every graph in the multiresolution sequence (default is False).
reg_epsfloat: 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

Gslist: 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):
...     fig, ax = Gs[idx].plot(title='Reduction level: {}'.format(idx))

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

Sparsify a graph (with Spielman-Srivastava).

Parameters

MGraph or sparse matrix: Graph structure or a Laplacian matrix
epsilonfloat: Sparsification parameter, which must be between 1/sqrt(N) and 1.
maxiterint, optional: Maximum number of iterations.
seed{None, int, RandomState, Generator}, optional: Seed for the random number generator (for reproducible sparsification).

Returns

MnewGraph or sparse matrix: New graph structure or sparse matrix

References

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

Examples

>>> from pygsp import reduction
>>> from matplotlib import pyplot as plt
>>> G = graphs.Sensor(100, k=20, distributed=True, seed=1)
>>> Gs = reduction.graph_sparsify(G, epsilon=0.4, seed=1)
>>> fig, axes = plt.subplots(1, 2)
>>> _ = G.plot(ax=axes[0], title='original')
>>> Gs.coords = G.coords
>>> _ = Gs.plot(ax=axes[1], title='sparsified')

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

Interpolate a graph signal.

Parameters

GGraph
f_subsampledndarray: A graph signal on the graph G.
keep_indsndarray: List of indices on which the signal is sampled.
orderint: Degree of the Chebyshev approximation (default = 100).
reg_epsfloat: 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_interpolatedndarray: 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

GGraph or sparse matrix: Graph structure or weight matrix
indlist: indices of the nodes to keep

Returns

GnewGraph 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

Gslist of graphs: A multiresolution sequence of graph structures.
fndarray: Graph signal to analyze.
h_filterslist: 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

candarray: Coarse approximation at each level
pendarray: Prediction error at each level
h_filterslist: 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

GsArray of Graphs: A multiresolution sequence of graph structures.
capndarray: Coarsest approximation of the original signal.
pendarray: Prediction error at each level.
use_exactbool: To use exact graph spectral filtering instead of the Chebyshev approximation.
orderint: Degree of the Chebyshev approximation (default=30).
least_squaresbool: To use the least squares synthesis (default=False).
h_filtersndarray: The filters used in the analysis operator. These are required for least squares synthesis, but not for the direct synthesis method.
use_landweberbool: To use the Landweber iteration approximation in the least squares synthesis.
reg_epsfloat: Interpolation parameter.
landweber_itsint: Number of iterations in the Landweber approximation for least squares synthesis.
landweber_taufloat: Parameter for the Landweber iteration.

Returns

reconstructionndarray: The reconstructed signal.
candarray: 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

GGraph: Graph structure of a tree.
NlevelNumber of times to downsample and coarsen the tree
rootint: The index of the root of the tree. (default = 1)
reduction_methodstr: The graph reduction method (default = ‘resistance_distance’)
compute_full_eigenbool: To also compute the graph Laplacian eigenvalues for every tree in the sequence

Returns

Gsndarray: Ndarray, with each element containing a graph structure represent a reduced tree.
subsampled_vertex_indicesndarray: 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

GGraph: Graph on which to compute the spectrogram.
atomfunc: Kernel to use in the spectrogram (default = exp(-M*(x/lmax)²)).
Mint (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 to solve convex optimization problems 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

lminfloat: Smallest non-zero eigenvalue.
lmaxfloat: Largest eigenvalue, i.e. pygsp.graphs.Graph.lmax.
Nscalesint: Number of scales.

Returns

scalesndarray: 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

xndarray: First colon vector
yndarray: Second colon vector

Returns

dndarray: 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.loadmat(path)[source]¶

Load a matlab data file.

Parameters

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

Returns

datadict: 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.rescale_center(x)[source]¶

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

Parameters

xndarray: Data to be rescaled.

Returns

rndarray: 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

GGraph or sparse matrix: Graph structure or Laplacian matrix (L)

Returns

rdsparse matrix: distance matrix

References

[KRandic93]

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

Symmetrize a square matrix.

Parameters

Warray_like

Square matrix to be symmetrized

methodstring

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

Notes

You can have the sum by multiplying the average by two. It is however not a good candidate for this function as it modifies an already symmetric matrix.

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.]])
>>> 2 * utils.symmetrize(W, method='average')
array([[0., 6., 4.],
       [6., 2., 8.],
       [4., 8., 6.]])
>>> 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.]])

Changelog¶

All notable changes to this project will be documented in this file. The format is based on Keep a Changelog and this project adheres to Semantic Versioning.

Unreleased¶

print(graph) and print(filters) now show valuable information.
Building a graph object is much faster.
New rectangular filter (low-pass and band-pass).
The exponential window has been updated from low-pass only to band-pass.
Much better documentation for the coherence of the Fourier basis.
Removed translate and modulate (they were not working and have no real use).
Fixed and documented vertex-frequency transforms. They are now implemented as filter banks.
Directed graphs are now completely supported.
The differential operator (D, grad, div) is better tested and documented.
G.dirichlet_energy(sig) computes the Dirichlet energy of a signal.
Better documentation of the frame and its bounds.
g.inverse() returns the pseudo-inverse of the filter bank.
g.complement() returns the filter that makes the frame tight.
Wave filter bank which application simulates the propagation of a wave.
Continuous integration with Python 3.7, 3.8, 3.9. Dropped 2.7, 3.4, 3.5, 3.6.
New implementation of the Sensor graph that is simpler and scales better.
A new learning module with three functions to solve standard semi-supervised classification and regression problems.
Import and export graphs and their signals to NetworkX and graph-tool.
Save and load graphs and theirs signals to / from GraphML, GML, and GEXF.
Documentation: path graph linked to DCT, ring graph linked to DFT.
We now have a gallery of examples! That is convenient for users to get a taste of what the library can do, and to start working from a code snippet.
Merged all the extra requirements in a single dev requirement.
New star graph (implemented as a comet without a tail).

Experimental filter API (to be tested and validated):

evaluate a filter bank with g(values)
filter with g(graph) @ signal
get an array representation (the frame) with g.toarray()
index the len(g) filters of a filter bank with g[idx]
concatenate filter banks with g + h

Plotting:

The plotting interface was updated to be more user-friendly. First, the documentation is now shown for filter.plot(), graph.plot(), and co. Second, the API in the plotting library has been deprecated. That module is now mostly for implementation only. Third, graph.plot() and graph.plot_signal() have been merged. As such, plot_signal() is deprecated. Finally, the following parameter names were changed:

plot_name => title
plot_eigenvalues => eigenvalues
show_sum => sum
show_edges => edges
vertex_size => size
npoints => n
save_as was removed

Other changes regarding plotting:

Plotting functions return matplotlib figures and axes.
Nodes, edges, and filters are plotted in transparency to avoid occlusion.
The node index can be printed on top of nodes to identify them easily.
Two vertex signals can now be plotted together as vertex color and size.
Two edges signals can be plotted as edge color and width.
Much faster (10 to 100 times faster) edge plotting with matplotlib.

There are many other small changes, look at the git history for the details.

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.

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 fresh virtual environment) for local development with the following:

$ git clone https://github.com/epfl-lts2/pygsp.git
$ pip install --upgrade --editable pygsp[dev]

The dev “extras requirement” ensures that dependencies required for development (to run the test suite and build the documentation) are installed. Only graph-tool will be missing: install it manually as it cannot be installed by pip.

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

Update README.rst and CHANGELOG.rst if applicable.

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.

To iterate faster, you can partially run the test suite, at various degrees of granularity, as follows:

$ python -m unittest pygsp.tests.test_docstrings.suite_reference
$ python -m unittest pygsp.tests.test_graphs.TestImportExport
$ python -m unittest pygsp.tests.test_graphs.TestImportExport.test_save_load

Making a release¶

Update the version number and release date in setup.py, pygsp/__init__.py and CHANGELOG.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.
Update the conda feedstock (at least the version number and sha256 in recipe/meta.yaml) by sending a PR to conda-forge.

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

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.
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.
KRandic93: Douglas J Klein and M Randić. Resistance distance. Journal of Mathematical Chemistry, 12(1):81–95, 1993.
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.