Try demmel's implementation.

report/main.tex

@@ -8,21 +8,17 @@
 \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,
@@ -46,14 +42,6 @@
 \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}
@@ -66,8 +54,8 @@
 \newtheorem{remark}[theorem]{Remark}
 
 \title{%
-    Accelerated Filtering on Graphs using Lanczos method
+    Accelerated filtering on graphs using Lanczos method
-    \\ \large Scientific Computing project report}
+    \\ \large Relazione del progetto di Calcolo Scientifico}
 \author{Alberto Defendi}
 
 \date{}
@@ -85,170 +73,24 @@
 {\setlength{\parskip}{0em}
 \tableofcontents}
 
-\section{Introduction}
+\section{Introduzione}
 
-We introduce basic graph theory concepts and briefly overview the results used in the project experiments.
+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 
 
-\subsection{Signal processing on graphs}
+Nell'analisi consideriamo i grafi di Erdo''s-Reiny (Figura)
 
-We consider a weighted undirected graph $ \G = (\V, \E, \mathcal{W}) $,
+Consideriamo un grafo non diretto e pesato $ G = (V, E, 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}
-
-
-
-\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(\L)s$.
+attraverso la valutazione $g(\mathcal{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}

src/afgl/main.py

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

src/afgl/test_1.py

@@ -14,7 +14,6 @@ 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
@@ -53,7 +52,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(1 / 2 * np.pi * (np.cos(np.pi * t) ** 2))
+    return np.sin(0.5 * np.pi * np.cos(np.pi * t) ** 2)
 
 
 """
@@ -71,21 +70,19 @@ 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
 
 
@@ -142,14 +139,14 @@ def run_comparison_1_for_graph(
     j = 3
     V, alp, beta = lanczos(L, s, M_MAX + j)
 
-    lanczos_err = np.zeros(M_MAX)
+    lanczos_err = np.zeros(M_MAX + j)
-    true_err = np.zeros(M_MAX)
+    true_err = np.zeros(M_MAX + j)
 
     GLs = filter_signal_with_fourier(G, s)
 
-    for M in range(1, M_MAX + 1):
+    for M in range(2, M_MAX + j):
-        g_M = compute_g_M(V[:, :M], alp[:M], beta[: M - 1], s)
+        g_M = compute_g_M(V[:, 0:M], alp[0:M], beta[0 : M - 1], s)
-        g_Mj = compute_g_M(V[:, : M + j], alp[: M + j], beta[: M + j - 1], s)
+        g_Mj = compute_g_M(V[:, 0 : M + j], alp[0 : M + j], beta[0 : M + j - 1], s)
 
         lanczos_err[M - 1] = LA.norm(g_Mj - g_M)
         true_err[M - 1] = LA.norm(GLs - g_M)
@@ -160,10 +157,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 = 200
+    M_MAX = 200
     p = 0.04
 
     s = np.random.randint(1, 10000, N).astype(float)
@@ -173,8 +170,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)
+    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)
+    l_err_S, t_err_S = run_comparison_1_for_graph(G_S, s, M_MAX)
 
     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)

src/afgl/util/build_T_matrix.py

@@ -2,14 +2,4 @@ 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)

src/afgl/util/lanczos.py

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

tests/test_main.py

@@ -1,31 +1,17 @@
 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():
+def test_lanczos_return_correct_solution():
-    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 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
@@ -42,35 +28,3 @@ def test_lanczos_return_correct_solution_with_dense():
     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