Compare commits

...

6 Commits

Author SHA1 Message Date
alberto fdbe1b2d55 Major fixes to pipeline.
Error now are close in plot. Still Lanczos orthogonalization selection
must be fixed.
1 month ago
alberto 23a95fe1f2 Minor fix. 1 month ago
alberto ec756aaf8d Minor fix to breakdown. 1 month ago
alberto a57b1a4496 Add theoretical introduction. 1 month ago
alberto 710f3da157 Match order with pseudocode. 1 month ago
alberto dcc0916063 Rename v_tilde to w. 1 month ago

@ -8,17 +8,21 @@
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{mathtools}
\usepackage{latexsym}
\usepackage{graphicx}
\usepackage{float}
\usepackage{etoolbox}
\usepackage{hyperref}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.18}
\usepackage{lipsum}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{mathtools}
\usepackage{nccmath}
\usepackage{eucal}
\usepackage[most]{tcolorbox}
\newtcolorbox[auto counter]{problem}[1][]{%
enhanced,
@ -42,6 +46,14 @@
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\V}{\mathcal{V}}
\newcommand{\E}{\mathcal{E}}
\newcommand{\G}{\mathcal{G}}
\renewcommand{\L}{\mathcal{L}}
\newcommand{\defeq}{\vcentcolon=}
\newcommand{\eqdef}{=\vcentcolon}
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
@ -54,8 +66,8 @@
\newtheorem{remark}[theorem]{Remark}
\title{%
Accelerated filtering on graphs using Lanczos method
\\ \large Relazione del progetto di Calcolo Scientifico}
Accelerated Filtering on Graphs using Lanczos method
\\ \large Scientific Computing project report}
\author{Alberto Defendi}
\date{}
@ -73,24 +85,170 @@
{\setlength{\parskip}{0em}
\tableofcontents}
\section{Introduzione}
\section{Introduction}
Introduciamo alcuni concetti di teoria dei grafi e alcuni risultati del corso che verranno usati nel corso della sperimentazione.
Scopo del progetto è verificare numericamente i risultati
We introduce basic graph theory concepts and briefly overview the results used in the project experiments.
Nell'analisi consideriamo i grafi di Erdo''s-Reiny (Figura)
\subsection{Signal processing on graphs}
Consideriamo un grafo non diretto e pesato $ G = (V, E, W)$.
We consider a weighted undirected graph $ \G = (\V, \E, \mathcal{W}) $,
where $ \V \subset \R^N $ is the set of vertices, $ \E \subset \R^M $ is the set of edges, and $ \mathcal{W}
: \V \times \mathcal{V} \to \R$ is a weight function. The weight function can be
represented with a $ N \times N $ matrix $ W $ such that
$W_{i,j} =
\begin{cases}
W(v_i, v_j), & \text{if } (v_i, v_j) \in \E \\
0, & \text{otherwise}
\end{cases}
\quad \text{for all } i,j = 1, \dots, |\V|
$,
and for our needs, we assume $ W $ to be symmetric, that is $ W_{i,j} = W_{j,i} $.
In the study of signal processing on graph, we model the idea of sending a signal to a node as
assigning a value to a vertex with a function $ s : \V \to \R$, whose values can be represented as a vector
$ s = s(\V) = [s_1, \dots, s_N] \in \R^N $ where each entry $ s_i \defeq s(v_i)$ represents
the signal sent over a node $ v_i \in \V$. We also keep track of the weight of each vertex
with a function $ d(i) \defeq \sum_{j=1}^N W_{i,j} $, that we represent with the matrix $D =
\text{diag}(d(1), \dots, d(N))$.
In this setting, we introduce the graph Laplacian $\L$ defined as
$ \L = D - W$, where $D $ is the diagonal degree matrix with entries $D_{ii} = d(i)$. By
construction, $\L^T = \mathcal{L}$, thus the Lanczos method can be safely applied to our
problem.
\begin{remark}
The Laplacian $\L$ is symmetric and semi-definite (it is diagonally dominant), hence by the
spectral theorem it admits the decomposition
\begin{equation*}
\L = U \Lambda U^{*},
\end{equation*}
where $ U = \left[ u_0,\dots,u_{N-1} \right] \in O(N)$ is called Fourier basis, and $ \Lambda =
\text{diag}\left(\lambda_0, \dots, \lambda_{N-1}\right)$ is the matrix of eigenvalues of $ \L $,
without loss of generality we can assume $ 0 =\lambda_0 \leq \lambda_1 \dots \leq \lambda_{N-1} $
and the vectors in $ U $ to be in the same order that the eigenvalues.
\footnote{this assumption is aligned with the PyGSP $ \text{compute\_fourier\_basis()} $ function
implementation, used in this work to compute the matrix $ \Lambda $. Clarifying this assumption is
outside our scope.}.
\end{remark}
\begin{definition}[Graph signal]
A graph signal is a continuous function $ g : \R^+ \to \R $.
\end{definition}
Diagonalising the Laplacian would give an easy way to compute the function $ g(\L) $ by
evaluating it on the eigenvalues of $ \L $, which means computing $ g(\L) = U g(\Lambda) U^{*} $.
In matrix notation, applying a filter to a graph $ \G $ corresponds to the operation
\begin{equation*}
s^\prime \defeq g(\L)s = U g(\L) U^{*}s.
\end{equation*}
\begin{remark}[Computational cost]
Computing the Fourier basis for $ \L $ is computationally expensive for large graphs, thus the
choice of using the Lanczos method for $ g(\L) $, with its computational cost of $ O(M \dot |\E|) $,
it offers an efficient alternative in computing $ g(\L) $. However, storing the basis $ V_M $ costs MN
additional memory, that could be avoided using a two-step implementation, that we leave for future
work.
\end{remark}
\subsection{The Lanczos method}
\begin{definition}
Given a matrix $ A \in \R^{N \times N} $ and a vector $ b \in R^N $, the Krylov subspace of order $j$ is defined as the set $ \mathcal{K}_j (A,b) = \{ b, Ab, A^2b, \dots,
A^{j-1}b\} $. We represent the basis of this subspace in a matrix $ V_M = \left[ v_1,\dots,v_M
\right] \in \R^{N \times M}$.
\end{definition}
We now consider the Arnoldi relation
\begin{equation*}
AV_j = V_jH_j + h_{j+1,j}v_{j+1}e_j^{*}, \\ \hspace{20pt}
H_j = V_j^{*}AV_j =
\begin{bmatrix}
h_{11} & \dots & \dots & h_{1,j} \\
h_{21} & h_{22} & & \vdots \\
& \ddots & \ddots & \vdots \\
& & h_{j,j-1} & h_{j,j}.
\end{bmatrix}
\end{equation*}
Letting
\begin{equation*}
\alpha_j \defeq h_{j,j},\hspace{20pt}\beta_j \defeq h_{j-1,j},
\end{equation*}
if $ A $ is symmetric and positive defined, then so is $ H_j $. Because $ H_j $ is both symmetric
and upper and lower Hessenberg matrix, then it is tridiagonal, and we refer to it as $ T_j $, hence
the Arnoldi relation takes the form:
\begin{equation*}
AV_j = V_j \alpha_j + \beta_j v_{j+1}e_j^{T}, \\ \hspace{20pt}
T_j = V_j^{\top} A V_j = \begin{bmatrix}
\alpha_1 & \beta_1 \\
\beta_1 & \ddots & \ddots \\
& \ddots & \ddots & \beta_{j-1} \\
& & \beta_{j-1} & \alpha_j
\end{bmatrix}.
\end{equation*}.
HERE PUT ALGORITHM
\section{Esperimento 1}
Studiamo i grafi di Erdos-Reiny e di tipo Sensors. Dal plot possiamo
\section{Experiments}
\subsection{}
Consider the filter $ g : [0, \lambda_{\text{max}}] \to \R$ and a signal vector $ s \in \R^N $, by a
a result of Gallopolus and Saad (see ?) it holds\footnote{This results holds outside the graph
signal processing context, that is for any function $f :\Omega \to \C$, where $ \Omega \subset \Lambda(A) $.} that
\begin{equation}
g(\L)s \approx \norm{s}_2 V_M g(T_M) e_1 \eqdef g_M \label{eq:1}
\end{equation}
\begin{definition}(Errors)
We define the Lanczos iteration error as $\norm{g_{M+j} - g_M}_2 $, where $ j $ is small, and the true
error as $ \norm{e_M} = \norm{g(\L)s - g_M} $.
\end{definition}
The scope of this experiment is verifying numerically equation \eqref{eq:1}
Studiamo i grafi di Erdos-Reiny e di tipo Sensors. Dal plot possiamo
Figura (dida: Grafi di ER e sensor colorati in base al segnale (non filtrato, sopra) e filtrato
attraverso la valutazione $g(\mathcal{L})s$.
attraverso la valutazione $g(\L)s$.
test
\begin{tikzpicture}
\begin{axis}[
% Size parameters suitable for standard 2-column IEEE/Springer papers
width=8cm,
height=5.5cm,
axis lines=middle,
xlabel={$t$},
ylabel={$f(t)$},
xmin=-1.2, xmax=1.2,
ymin=-0.2, ymax=1.2,
% Clean fractional ticks
xtick={-1, -0.5, 0, 0.5, 1},
xticklabels={$-1$, $-\frac{1}{2}$, $0$, $\frac{1}{2}$, $1$},
ytick={0, 1},
ticklabel style={font=\footnotesize},
xlabel style={anchor=west},
ylabel style={anchor=south},
enlargelimits=false,
% Define the function natively.
% Note: pgfplots evaluates trig functions in degrees by default.
% 1/2 pi rad = 90 deg; pi*t rad = 180*t deg.
declare function={
f(\t) = (abs(\t) <= 0.5) * sin(90 * (cos(180*\t))^2);
}
]
\addplot [
domain=-1.2:1.2,
samples=300, % High sample count for smooth paper-quality rendering
thick,
blue!80!black % Professional dark blue (prints better than pure blue)
] {f(x)};
\end{axis}
\end{tikzpicture}
\clearpage
\bibliographystyle{unsrt}

@ -14,6 +14,7 @@ from afgl.util.plot import latex_log_formatter
def plot_graphs(G_ER, G_Sensor, s: np.ndarray, N: int, p: float) -> None:
"""Visualization of signal being filtered of two different types of graphs."""
fig, axs = plt.subplots(2, 2, figsize=(6.6, 5))
# Set coordinates
@ -52,7 +53,7 @@ def plot_graphs(G_ER, G_Sensor, s: np.ndarray, N: int, p: float) -> None:
def g_extended(t: np.ndarray) -> np.ndarray:
return np.sin(0.5 * np.pi * np.cos(np.pi * t) ** 2)
return np.sin(1 / 2 * np.pi * (np.cos(np.pi * t) ** 2))
"""
@ -70,19 +71,21 @@ def g(T: np.ndarray) -> np.ndarray:
return np.where(Chi, g_extended(T), 0)
"""
Computes the approximation g_M (see [1]) using Lanczos
"""
def compute_g_M(
V: np.ndarray, alp: np.ndarray, beta: np.ndarray, s: np.ndarray
) -> np.ndarray:
"""
Computes the approximation g_M (see [1]) using Lanczos
"""
M = len(alp)
e_1 = np.zeros(M)
e_1[0] = 1
T = build_T_matrix(alp, beta)
y = LA.norm(s) * (g(T) @ e_1)
eigvals, eigvecs = LA.eigh(T)
g_T = eigvecs @ np.diag(g(eigvals)) @ eigvecs.T
y = LA.norm(s) * (g_T @ e_1)
return V @ y
@ -139,14 +142,14 @@ def run_comparison_1_for_graph(
j = 3
V, alp, beta = lanczos(L, s, M_MAX + j)
lanczos_err = np.zeros(M_MAX + j)
true_err = np.zeros(M_MAX + j)
lanczos_err = np.zeros(M_MAX)
true_err = np.zeros(M_MAX)
GLs = filter_signal_with_fourier(G, s)
for M in range(2, M_MAX + j):
g_M = compute_g_M(V[:, 0:M], alp[0:M], beta[0 : M - 1], s)
g_Mj = compute_g_M(V[:, 0 : M + j], alp[0 : M + j], beta[0 : M + j - 1], s)
for M in range(1, M_MAX + 1):
g_M = compute_g_M(V[:, :M], alp[:M], beta[: M - 1], s)
g_Mj = compute_g_M(V[:, : M + j], alp[: M + j], beta[: M + j - 1], s)
lanczos_err[M - 1] = LA.norm(g_Mj - g_M)
true_err[M - 1] = LA.norm(GLs - g_M)
@ -157,10 +160,10 @@ def run_comparison_1_for_graph(
def run() -> None:
"""Ripete il test corrispondente ad Example 1 dell'articolo limitandosi al
metodo di Lanczos (no Chebyshev) e utilizzando come funzione g(t) = sin(0.5π
cos(πt)2) * \chi_{[-0.5, 0.5]}.
cos(πt)2) * chi_{[-0.5, 0.5]}.
"""
N = 500
M_MAX = 200
M = 200
p = 0.04
s = np.random.randint(1, 10000, N).astype(float)
@ -170,8 +173,8 @@ def run() -> None:
G_ER = graphs.ErdosRenyi(N, p)
G_S = graphs.Sensor(N)
l_err_ER, t_err_ER = run_comparison_1_for_graph(G_ER, s, M_MAX)
l_err_S, t_err_S = run_comparison_1_for_graph(G_S, s, M_MAX)
l_err_ER, t_err_ER = run_comparison_1_for_graph(G_ER, s, M)
l_err_S, t_err_S = run_comparison_1_for_graph(G_S, s, M)
plot_error_comparison(l_err_ER, t_err_ER, l_err_S, t_err_S)
plot_graphs(G_ER, G_S, s, N, p)
# plot_graphs(G_ER, G_S, s, N, p)

@ -1,13 +1,13 @@
import sys
import afgl.test_2 as t_2
import afgl.ex_1 as ex_1
from afgl.util.plot import plot_setup
def run() -> None:
plot_setup()
# t_1.run()
t_2.run()
ex_1.run()
# t_2.run()
if __name__ == "__main__":

@ -2,4 +2,14 @@ import numpy as np
def build_T_matrix(alp, beta):
"""Constructs Lanczos T tridiagonal matrix.
Args:
alp: Vector of alphas of size N
beta: Vector of betas of size N-1
Returns:
T : Matrix T
"""
return np.diag(alp) + np.diag(beta, -1) + np.diag(beta, 1)

@ -1,26 +1,39 @@
import numpy as np
import numpy.linalg as LA
"""
Classic Lanczos method (without re-orthogonalization)
Arguments
L : Real valued NxN symmetric matrix
s : vector of size N
M : natural number indicating basis size
def double_orthogonalization(V, w, j):
# Why it works well until j+1? it should be until j-1 from Demmel
for _ in range(2):
w -= V[:, : j + 1] @ (V[:, : j + 1].T @ w)
return w
Returns
-------
V : ndarray
M-dimensional vector with orthonormal columns.
alp : ndarray
M-dimensional array of scalars.
beta : ndarray
M-dimensional array of scalars.
"""
def no_orthogonalization(V, w, alp, beta, j):
w = w - alp[j] * V[:, j]
if j > 0:
w = w - beta[j - 1] * V[:, j - 1]
return w
def lanczos(L, s, M):
"""
Classic Lanczos method (without re-orthogonalization)
Arguments
L : Real valued NxN symmetric matrix
s : vector of size N
M : natural number indicating basis size
Returns
-------
V : ndarray
M-dimensional vector with orthonormal columns.
alp : ndarray
M-dimensional array of scalars.
beta : ndarray
M-dimensional array of scalars.
"""
N = len(s)
alp = np.zeros(M)
beta = np.zeros(M - 1)
@ -31,14 +44,13 @@ def lanczos(L, s, M):
w = L @ V[:, j]
alp[j] = np.dot(V[:, j], w)
v_tilde = w - V[:, j] * alp[j]
if j > 0:
v_tilde = v_tilde - V[:, j - 1] * beta[j - 1]
w = no_orthogonalization(V, w, alp, beta, j)
if j < M - 1:
beta[j] = LA.norm(v_tilde)
if beta[j] == 0:
break
V[:, j + 1] = v_tilde / beta[j]
beta[j] = LA.norm(w)
if beta[j] < 1e-14:
print("BREAKDOWN")
return V[:, : j + 1], alp[: j + 1], beta[:j]
V[:, j + 1] = w / beta[j]
return [V, alp, beta]
return V, alp, beta

@ -1,17 +1,31 @@
import numpy as np
import numpy.linalg as LA
from afgl.ex_1 import compute_g_M, filter_signal_with_fourier, g
from afgl.util.build_T_matrix import build_T_matrix
from afgl.util.lanczos import lanczos
from pygsp import graphs
"""
Todo: better test case
"""
def g_evaluation_should_respect_chi():
A = np.array([[1, 0, 1], [0, 0, 0], [1, 0, 0]])
expected_gA = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 1]])
assert (expected_gA == g(A).astype(int)).all()
def test_lanczos_return_correct_solution():
def g_evaluation_should_return_matrix_of_zeros():
A = 1 / 2 * np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
assert (g(A).astype(int) == np.zeros((3, 3))).all()
def test_lanczos_return_correct_solution_with_dense():
"""Tests correctness of solution of Lx=s comparing Lanczos projected
solution with numpy solve function.
"""
N = 1000
M = 999
# Generate a good conditioned matrix
eigvals = np.random.uniform(10000, 100000, N)
Q, _ = LA.qr(np.random.randn(N, N))
L = Q @ np.diag(eigvals) @ Q.T
@ -28,3 +42,35 @@ def test_lanczos_return_correct_solution():
x_lanczos = V @ y
assert LA.norm(x - x_lanczos) < 1e-10
def test_function_g_with_graph_laplacian():
N = 1000
p = 0.04
j = 3
n = 5
M_VALS = 25 * (2 ** np.arange(n))
M_VALS = [200]
for M in M_VALS:
s = np.random.randint(1, 10000, N).astype(float)
# Normalize s as in request
s /= LA.norm(s)
GRAPHS = [graphs.ErdosRenyi(N, p), graphs.Sensor(N)]
for G in GRAPHS:
G.compute_laplacian()
L = G.L
V, alp, beta = lanczos(L, s, M + j)
GLs = filter_signal_with_fourier(G, s)
g_M = compute_g_M(V[:, 0:M], alp[0:M], beta[0 : M - 1], s)
g_Mj = compute_g_M(V[:, 0 : M + j], alp[0 : M + j], beta[0 : M + j - 1], s)
diff = LA.norm(g_Mj - g_M)
e_M = LA.norm(GLs - g_M)
assert abs(diff - e_M) < 1e-2

Loading…
Cancel
Save