Try demmel's implementation.

Errors are distant

pyproject.toml

@@ -7,6 +7,7 @@ dependencies = [
     "matplotlib>=3.10.8",
     "numpy>=2.0.0",
     "pygsp>=0.6.1",
+    "scienceplots>=2.2.1",
     "scipy>=1.10.0",
 ]
 requires-python = ">=3.10"

report/main.tex

@@ -1,45 +1,102 @@
-% !TeX program = pdflatex
-\documentclass{mathreport} % Uses our custom mathreport.cls
+\documentclass[11pt]{article}
+\usepackage[margin=1.2in]{geometry}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage[english]{babel}
+\usepackage{fourier}
+\usepackage{amsthm}
+\usepackage{amssymb}
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{latexsym}
+\usepackage{graphicx}
+\usepackage{float}
+\usepackage{etoolbox}
+\usepackage{hyperref}
+\usepackage{tikz}
+\usepackage{lipsum}
+\usepackage{algorithm}
+\usepackage{algpseudocode}
+\usepackage{mathtools}
+\usepackage{nccmath}
+\usepackage[most]{tcolorbox}
+\newtcolorbox[auto counter]{problem}[1][]{%
+    enhanced,
+    breakable,
+    colback=white,
+    colbacktitle=white,
+    coltitle=black,
+    fonttitle=\bfseries,
+    boxrule=.6pt,
+    titlerule=.2pt,
+    toptitle=3pt,
+    bottomtitle=3pt,
+    title=GitHub repository of this project}
 
-% Load bibliography
-\addbibresource{references.bib}
-\usepackage{lipsum} % Just for the demo text
+\RequirePackage[activate={true,nocompatibility},final,tracking=true,kerning=true,spacing=true]{microtype}
+\SetTracking{encoding={*}, shape=sc}{40} % Spacing for Small Caps
 
-% --- Document Metadata ---
-\title{\normalfont\scshape\Large The Geometry of Complex Systems}
-\author{\normalfont\itshape A. N. Other}
-\date{\small\today}
+
+\newcommand{\R}{\mathbb{R}}
+\newcommand{\N}{\mathbb{N}}
+\newcommand{\Z}{\mathbb{Z}}
+\newcommand{\Q}{\mathbb{Q}}
+\newcommand{\C}{\mathbb{C}}
+
+\newtheorem{theorem}{Theorem}[section]
+\newtheorem{lemma}[theorem]{Lemma}
+\newtheorem{proposition}[theorem]{Proposition}
+\newtheorem{corollary}[theorem]{Corollary}
+\theoremstyle{definition}
+\newtheorem{definition}[theorem]{Definition}
+\newtheorem{example}[theorem]{Example}
+\theoremstyle{remark}
+\newtheorem{remark}[theorem]{Remark}
+
+\title{%
+    Accelerated filtering on graphs using Lanczos method
+    \\ \large Relazione del progetto di Calcolo Scientifico}
+\author{Alberto Defendi}
+
+\date{}
+
+\setlength{\parskip}{1em}
+\setlength{\parindent}{0em}
 
 \begin{document}
-
 \maketitle
 
 \begin{abstract}
-    \noindent \small \textbf{\textit{Abstract.}} \lipsum[1][1-4]
+    \noindent The Lanczos algorithm \ldots 
 \end{abstract}
 
-\section{Introduction}
-This document uses a custom class file (\texttt{mathreport.cls}). This keeps the main file clean. The typography is set to Bringhurst's standards: wide margins, Palatino font, and old-style figures (e.g., 12345).
+{\setlength{\parskip}{0em}
+\tableofcontents}
 
-\section{Mathematical Theory}
-We define our primary operator in the Hilbert space $\mathcal{H}$.
+\section{Introduzione}
 
-\begin{definition}[Compact Operator]
-    An operator $T: X \to Y$ is compact if $\overline{T(B_X)}$ is compact in $Y$.
-\end{definition}
+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 
 
-\begin{theorem}[Spectral Theorem]
-    There exists an orthonormal basis of eigenvectors.
-\end{theorem}
+Nell'analisi consideriamo i grafi di Erdo''s-Reiny (Figura)
 
-Consider the harmonic series shown in \cref{eq:harmonic}. The styling is handled entirely by the external class file.
-\begin{equation} \label{eq:harmonic}
-    H_n = \sum_{k=1}^n \frac{1}{k} \approx \ln n + \gamma
-\end{equation}
+Consideriamo un grafo non diretto e pesato $ G = (V, E, W)$.
+
+
+\section{Esperimento 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$.
+test
+
+
+\clearpage
+\bibliographystyle{unsrt}
+\bibliography{ref}
+\nocite{*}
 
-\section{Results and Discussion}
-\lipsum[2-4]
 
-\printbibliography
 
 \end{document}

report/ref.bib

@@ -0,0 +1,9 @@
+@article{susnjara2015,
+      title={Accelerated filtering on graphs using Lanczos method}, 
+      author={Ana Susnjara and Nathanael Perraudin and Daniel Kressner and Pierre Vandergheynst},
+      year={2015},
+      eprint={1509.04537},
+      archivePrefix={arXiv},
+      primaryClass={math.NA},
+      url={https://arxiv.org/abs/1509.04537}, 
+}

report/references.bib

@@ -1,7 +0,0 @@
-@book{bringhurst2004,
-  author    = {Robert Bringhurst},
-  title     = {The Elements of Typographic Style},
-  year      = {2004},
-  publisher = {Hartley \& Marks},
-  address   = {Vancouver}
-}

src/afgl/lanczos.py

@@ -1,40 +0,0 @@
-import numpy as np
-import numpy.linalg as LA
-
-"""
-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.
-"""
-
-
-def lanczos(L, s, M):
-    N = len(s)
-    alp = np.zeros(M)
-    beta = np.zeros(M - 1)
-    V = np.zeros((N, M))
-    V[:, 0] = s / LA.norm(s)
-
-    for j in range(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]
-
-        if j < M - 1:
-            beta[j] = LA.norm(v_tilde)
-            V[:, j + 1] = v_tilde / beta[j]
-
-    return [V, alp, beta]

src/afgl/main.py

@@ -1,10 +1,13 @@
-# src/afgl/main.py
 import sys
 
+import afgl.test_2 as t_2
+from afgl.util.plot import plot_setup
 
-def run():
-    """Main execution function."""
-    print("Starting the Accelerated Filtering Graphs Lanczos (AFGL) process...")
+
+def run() -> None:
+    plot_setup()
+    # t_1.run()
+    t_2.run()
 
 
 if __name__ == "__main__":

src/afgl/test_1.py

@@ -0,0 +1,177 @@
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+import numpy as np
+import numpy.linalg as LA
+
+# Requires latex installed
+import scienceplots  # noqa: F401
+import scipy
+from pygsp import graphs
+
+from afgl.util.build_T_matrix import build_T_matrix
+from afgl.util.lanczos import lanczos
+from afgl.util.plot import latex_log_formatter
+
+
+def plot_graphs(G_ER, G_Sensor, s: np.ndarray, N: int, p: float) -> None:
+    fig, axs = plt.subplots(2, 2, figsize=(6.6, 5))
+
+    # Set coordinates
+    G_ER.set_coordinates()
+    G_Sensor.set_coordinates()
+
+    signal_ER = filter_signal_with_fourier(G_ER, s)
+    signal_S = filter_signal_with_fourier(G_Sensor, s)
+
+    # TOP LEFT
+    G_ER.plot(s, ax=axs[0, 0], vertex_size=15, edge_width=0.5, edge_color="gray")
+    axs[0, 0].set_title(rf"Erdős-Rényi Graph $(N = {N}, p = {p})$", pad=20)
+    axs[0, 0].set_axis_off()
+
+    # BOTTOM LEFT
+    G_ER.plot(
+        signal_ER, ax=axs[1, 0], vertex_size=15, edge_width=0.5, edge_color="gray"
+    )
+    axs[1, 0].set_title("", pad=20)
+    axs[1, 0].set_axis_off()
+
+    # TOP RIGHT
+    G_Sensor.plot(s, ax=axs[0, 1], vertex_size=15, edge_width=0.5, edge_color="gray")
+    axs[0, 1].set_title(rf"Sensor Network $(N = {N})$", pad=20)
+    axs[0, 1].set_axis_off()
+
+    # BOTTOM RIGHT
+    G_Sensor.plot(
+        signal_S, ax=axs[1, 1], vertex_size=15, edge_width=0.5, edge_color="gray"
+    )
+    axs[1, 1].set_title("", pad=20)
+    axs[1, 1].set_axis_off()
+
+    # Prevent label/title overlap
+    plt.savefig("./out/printed_graphs.pdf", bbox_inches="tight")
+
+
+def g_extended(t: np.ndarray) -> np.ndarray:
+    return np.sin(0.5 * np.pi * np.cos(np.pi * t) ** 2)
+
+
+"""
+Evaluates the function sin(0.5*pi*cos(pi*t)^2)chi_[-1/2,1/2] where chi_I is the
+characteristic function of I, as defined in example 1 (see [1]).
+"""
+
+
+def g(T: np.ndarray) -> np.ndarray:
+    if scipy.sparse.issparse(T):
+        # Operator & not supporting sparse matrix
+        T = T.toarray()
+    Chi = ((T >= -1 / 2) & (T <= 1 / 2)).astype(int)
+    # Apply g_extended where Chi is True, else output 0
+    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:
+    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)
+    return V @ y
+
+
+def filter_signal_with_fourier(G, s: np.ndarray) -> np.ndarray:
+    G.compute_fourier_basis()
+    U = G.U
+    return (U @ np.diag(g(G.e)) @ U.T) @ s
+
+
+def plot_error_comparison(
+    l_err_ER: np.ndarray, t_err_ER: np.ndarray, l_err_S: np.ndarray, t_err_S: np.ndarray
+) -> None:
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6.6, 2.5))
+
+    # Left plot (Erdos-Renyi)
+    ax1.plot(l_err_ER, label=r"$\left\lVert g_{M+3} - g_M \right\rVert_2$")
+    ax1.plot(t_err_ER, label=r"$\left\lVert e_M \right\rVert_2$")
+    ax1.set_title("Erdős-Rényi graph")
+
+    # Right plot (Sensor)
+    ax2.plot(l_err_S, label=r"$\left\lVert g_{M+3} - g_M \right\rVert_2$")
+    ax2.plot(t_err_S, label=r"$\left\lVert e_M \right\rVert_2$")
+    ax2.set_title("Sensor graph")
+
+    # Apply identical formatting to both subplots
+    for ax in (ax1, ax2):
+        ax.xaxis.set_major_locator(ticker.MultipleLocator(50))
+        ax.set_yscale("log")
+        ax.yaxis.set_major_formatter(ticker.FuncFormatter(latex_log_formatter))
+        ax.legend()
+
+    # Prevents overlapping of labels between the subplots
+    plt.tight_layout()
+
+    plt.savefig("./out/ex1_estimate.pdf", bbox_inches="tight")
+
+
+def run_comparison_1_for_graph(
+    G, s: np.ndarray, M_MAX: int
+) -> tuple[np.ndarray, np.ndarray]:
+    """Compares the error with error generated by Lanczos.
+
+    Args:
+        G: Graph
+        s: Signal vector
+
+    Returns:
+        [lanczos_err, true_err]: Vectors of errors norm(g_{M+3} - g_M) and norm(e_M)
+        as defined in [1]
+    """
+    G.compute_laplacian("combinatorial")
+    L = G.L
+
+    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)
+
+    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)
+
+        lanczos_err[M - 1] = LA.norm(g_Mj - g_M)
+        true_err[M - 1] = LA.norm(GLs - g_M)
+
+    return lanczos_err, true_err
+
+
+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]}.
+    """
+    N = 500
+    M_MAX = 200
+    p = 0.04
+
+    s = np.random.randint(1, 10000, N).astype(float)
+    # Normalize s as in request
+    s /= LA.norm(s)
+
+    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)
+
+    plot_error_comparison(l_err_ER, t_err_ER, l_err_S, t_err_S)
+    plot_graphs(G_ER, G_S, s, N, p)

src/afgl/test_2.py

@@ -0,0 +1,36 @@
+import numpy as np
+import numpy.linalg as LA
+from pygsp import graphs
+
+from afgl.util.lanczos import lanczos
+
+
+def run() -> None:
+    """Genera i grafi di di Erdos-Reny di grandezza crescente (ad esempio 250,
+    500,1000, 2000, 4000) e parametro p = 0.04 e misura il tempo
+    computazionale del metodo di Lanczos utilizzando come soglia per il criterio
+    d'arresto epsilon = 10^-2 (o una soglia a scelta).
+
+     Args:
+         None
+
+     Returns:
+         None
+    """
+    n = 6
+    p = 0.04
+
+    M = 200
+
+    N_VALUES = 250 * (2 ** np.arange(n))
+
+    for N in N_VALUES:
+        s = np.random.randint(1, 10000, N).astype(float)
+        # Normalize s as in request
+        s /= LA.norm(s)
+
+        G = graphs.ErdosRenyi(N, p)
+        G.compute_laplacian("combinatorial")
+        L = G.L
+
+        lanczos(L, s, M)

src/afgl/util/build_T_matrix.py

@@ -0,0 +1,5 @@
+import numpy as np
+
+
+def build_T_matrix(alp, beta):
+    return np.diag(alp) + np.diag(beta, -1) + np.diag(beta, 1)

src/afgl/util/lanczos.py

@@ -0,0 +1,43 @@
+import numpy as np
+import numpy.linalg as LA
+
+"""
+Classic Lanczos method (without re-orthogonalization)
+Using Demmel's book version.
+
+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.
+"""
+
+
+def lanczos(L, s, M):
+    N = len(s)
+    alp = np.zeros(M)
+    beta = np.zeros(M)
+    V = np.zeros((N, M + 1))
+    V[:, 1] = s / LA.norm(s)
+
+    for j in range(1, M):
+        w = L @ V[:, j]
+        alp[j] = np.dot(V[:, j], w)
+
+        w = w - V[:, j] * alp[j] - V[:, j - 1] * beta[j - 1]
+
+        beta[j] = LA.norm(w)
+        if beta[j] == 0:
+            print("Breakdown")
+            break
+        V[:, j + 1] = w / beta[j]
+
+    return [V[:, 1:], alp, beta[1:]]

src/afgl/util/plot.py

@@ -0,0 +1,28 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import scienceplots  # noqa: F401
+from pygsp import plotting
+
+
+def latex_sci(val: float, decimals: int = 2) -> str:
+    """Converts a value to LaTeX scientific notation A x 10^{B}."""
+    if val == 0:
+        return "0"
+    exponent = int(np.floor(np.log10(abs(val))))
+    mantissa = val / 10**exponent
+    return rf"{mantissa:.{decimals}f} \times 10^{{{exponent}}}"
+
+
+def latex_log_formatter(y: float, pos: int) -> str:
+    """Custom formatter to render tick labels as LaTeX 10^{n}."""
+    if y <= 0:
+        return ""
+    # Extract the exponent using log10
+    n = int(np.round(np.log10(y)))
+    return f"$10^{{{n}}}$"
+
+
+def plot_setup() -> None:
+    plotting.BACKEND = "matplotlib"
+    plt.style.use(["science"])
+    # TODO match font with document

tests/test_main.py

@@ -1,6 +1,7 @@
 import numpy as np
 import numpy.linalg as LA
-from afgl.lanczos import lanczos
+from afgl.util.build_T_matrix import build_T_matrix
+from afgl.util.lanczos import lanczos
 
 """
 Todo: better test case
@@ -18,7 +19,7 @@ def test_lanczos_return_correct_solution():
     s = np.random.randint(1, 10, N)
     [V, alp, beta] = lanczos(L, s, M)
 
-    T = np.diag(alp) + np.diag(beta, -1) + np.diag(beta, 1)
+    T = build_T_matrix(alp, beta)
 
     x = LA.solve(L, s)
     e_1 = np.zeros(M)

uv.lock

@@ -15,6 +15,7 @@ dependencies = [
     { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },
     { name = "numpy", version = "2.4.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" },
     { name = "pygsp" },
+    { name = "scienceplots" },
     { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },
     { name = "scipy", version = "1.17.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" },
 ]
@@ -24,6 +25,7 @@ requires-dist = [
     { name = "matplotlib", specifier = ">=3.10.8" },
     { name = "numpy", specifier = ">=2.0.0" },
     { name = "pygsp", specifier = ">=0.6.1" },
+    { name = "scienceplots", specifier = ">=2.2.1" },
     { name = "scipy", specifier = ">=1.10.0" },
 ]
 
@@ -737,6 +739,18 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/ec/57/56b9bcc3c9c6a792fcbaf139543cee77261f3651ca9da0c93f5c1221264b/python_dateutil-2.9.0.post0-py2.py3-none-any.whl", hash = "sha256:a8b2bc7bffae282281c8140a97d3aa9c14da0b136dfe83f850eea9a5f7470427", size = 229892, upload-time = "2024-03-01T18:36:18.57Z" },
 ]
 
+[[package]]
+name = "scienceplots"
+version = "2.2.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "matplotlib" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/ae/a5/5f858668ca1a513033a7f0d55cd12b0940a12e822f9f61f317ce344e07c6/scienceplots-2.2.1.tar.gz", hash = "sha256:51ad98c420e499d3284d07b6447a4b3aedd8ec122aca39ba91c58205226c408a", size = 17966, upload-time = "2026-02-25T01:26:54.603Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/cc/22/14e7b20f7d11f3e9ea9b5ae6acf9ea695c423eee785906cd5c5914a841dc/scienceplots-2.2.1-py3-none-any.whl", hash = "sha256:a1a9f670cbf5b59d92cdd3250be85b079ee4cf08bcaf2ace5ef57995be0b6c42", size = 30237, upload-time = "2026-02-25T01:26:53.298Z" },
+]
+
 [[package]]
 name = "scipy"
 version = "1.15.3"