GSP Graph TV Demo - Reconstruction of missing sample on a graph using TV

Description

Reconstruction of missing sample on a graph using TV

In this demo, we try to reconstruct missing sample of a piece-wise smooth signal on a graph. To do so, we will minimize the well-known TV norm defined on the graph.

For this example, you need the pyunlocbox. You can download it from https://github.com/epfl-lts2/pyunlocbox and installing it.

We express the recovery problem as a convex optimization problem of the following form:

\[arg \min_x \|\nabla(x)\|_1 \text{ s. t. } \|Mx-b\|_2 \leq \epsilon\]

Where \(b\) represents the known measurements, \(M\) is an operator representing the mask and \(\epsilon\) is the radius of the l2 ball.

We set:

  • \(f_1(x)=||\nabla x ||_1\)

We define the prox of \(f_1\) as:

\[prox_{f1,\gamma} (z) = arg \min_{x} \frac{1}{2} \|x-z\|_2^2 + \gamma \| \nabla z \|_1\]
  • \(f_2\) is the indicator function of the set S define by :math:||Mx-b||_2 < epsilon

We define the prox of \(f_2\) as

\[prox_{f2,\gamma} (z) = arg \min_{x} \frac{1}{2} \|x-z\|_2^2 + i_S(x)\]

with \(i_S(x)\) is zero if \(x\) is in the set \(S\) and infinity otherwise. This previous problem has an identical solution as:

\[arg \min_{z} \|x - z\|_2^2 \hspace{1cm} such \hspace{0.25cm} that \hspace{1cm} \|Mz-b\|_2 \leq \epsilon\]

It is simply a projection on the B2-ball.

Results and code

>>> from pygsp import graphs, plotting
>>> import numpy as np
>>>
>>> # Create a random sensor graph
>>> G = graphs.Sensor(N=256, distribute=True)
>>> G.compute_fourier_basis()
>>>
>>> # Create signal
>>> graph_value = np.copysign(np.ones(np.shape(G.U[:, 3])[0]), G.U[:, 3])
>>>
>>> plotting.plt_plot_signal(G, graph_value, savefig=True, plot_name='doc/tutorials/img/original_signal')
tutorials/img/original_signal.*

This figure shows the original signal on graph.

>>> # Create the mask
>>> M = np.random.rand(G.U.shape[0], 1)
>>> M = M > 0.6  # Probability of having no label on a vertex.
>>>
>>> # Applying the mask to the data
>>> sigma = 0.0
>>> depleted_graph_value = M * (graph_value.reshape(graph_value.size, 1) + sigma * np.random.standard_normal((G.N, 1)))
>>>
>>> plotting.plt_plot_signal(G, depleted_graph_value, show_edges=True, savefig=True, plot_name='doc/tutorials/img/depleted_signal')
tutorials/img/depleted_signal.*

This figure shows the signal on graph after the application of the mask and addition of noise. More than half of the vertices are set to 0.

tutorials/img/tv_recons_signal.*

This figure shows the reconstructed signal thanks to the algorithm.

Comparison with Tikhonov regularization

We can also use the Tikhonov regularizer that will promote smoothness. In this case, we solve:

\[arg \min_x \tau \|\nabla(x)\|_2^2 \text{ s. t. } \|Mx-b\|_2 \leq \epsilon\]

The result is presented as following:

tutorials/img/tik_recons_signal.*

This figure shows the reconstructed signal thanks to the algorithm.