Neuroni simulati ed applicazioni lineari, topologiche e di geometria differenziale.

Dialogo con Gemini.

Riprendiamo la questione della simulazione dei neuroni e dell’applicazione dell’algebra lineare e la topologia. Topologia o geometria differenziale?

La distinzione tra topologia (in particolare algebrica/computazionale) e geometria differenziale nella modellazione e simulazione neuronale non è una contrapposizione, ma una netta divisione di scala, invarianti e obiettivi teorici.

La scelta dipende essenzialmente dal livello del sistema che si intende formalizzare: il connettoma discreto e la dinamica globale di rete, oppure la dinamica di membrana e la geometria intrinseca delle traiettorie nello spazio delle fasi.

Quando vince la Topologia (Algebrica e Computazionale)

La topologia entra in gioco quando la metrica esatta (distanze euclidee, tempi continui precisi) è secondaria rispetto a connettività, invarianza per deformazione e buchi dimensionali (omologia).

• Topologia delle Reti e Complessi Simpliciali: Un grafo neuronale (algebra lineare con matrici di adiacenza e laplaciani di grafo) cattura solo interazioni a coppie. La topologia algebrica estende questo approccio a interazioni multi-neurone tramite complessi di cricche (clique complexes) e simplessi orientati. Un gruppo di $k$ neuroni che scaricano all’unisono o in sequenza precisa forma un simplesso di dimensione $k-1$.
• Omologia Persistente (TDA – Topological Data Analysis): Tracciando la co-attivazione di popolazioni neuronali al variare di una soglia di correlazione o di scala temporale, si calcolano i numeri di Betti ($\beta_0$ componenti connesse, $\beta_1$ cicli/anelli, $\beta_2$ cavità/vuoti).
  • Esempio paradigmatico: I complessi formati dalle place cells nell’ippocampo permettono al cervello di estrarre la topologia dell’ambiente (es. la presenza di ostacoli chiusi) senza bisogno di misurare coordinate cartesiane esatte.
• Laplaciani di Hodge su Complessi Simpliciali: Unione naturale tra algebra lineare e topologia. Generalizzando il Laplaciano di grafo ai simplessi superiori.
  Si studiano i flussi di informazione neuronale come $k$-forme discrete, separando componenti esatte, co-esatte e armoniche (teorema di decomposizione di Hodge).

Quando serve la Geometria Differenziale

La geometria differenziale è necessaria non appena il fenomeno dipende da curvatura, metrica Riemanniana, trasporto parallelo e campi vettoriali continui.

Spazio delle Fasi e Varietà Invarianti (Sistemi Dinamici): Un singolo neurone (Hodgkin-Huxley, FitzHugh-Nagumo, Morris-Lecar) vive in un sistema differenziale non lineare
I confini tra potenziale d’azione e stato di riposo sono determinati da varietà stabili e instabili di selle (separatrici), biforcazioni di Hopf o canard su varietà lente/veloci.

Geometria dell’Informazione (Amari): Lo spazio delle distribuzioni di probabilità associate agli stati di scarica di una popolazione neuronale forma una varietà Riemanniana, dove la metrica è data dal tensore metrico di Fisher.
La dinamica di apprendimento sinaptico (plasticità spike-timing o hebbiana) può essere formalizzata come un flusso di gradiente naturale lungo le geodetiche di questa varietà.
Mappe Corticali e Geometria Sub-Riemanniana: La corteccia visiva primaria (V1) modella orientazione, frequenza spaziale e posizione retinotopica. I modelli di Petitot, Citti e Sarti mostrano che il connettoma orizzontale di V1 approssima le geodetiche di una struttura di contatto / geometria sub-Riemanniana nel gruppo di Lie fibrato.

Quadro comparativo

Criterio	Topologia (Algebrica / TDA)	Geometria Differenziale
Oggetto di studio	Connettomi, firing corale, cicli di feedforward	Varietà lente, gradienti informazionali, mappe sensoriali
Strumento chiave	Omologia persistente, complessi simpliciali, $\beta_k$	Metrica Riemanniana, tensore di Fisher, fibrati, flussi geodetici
Invarianza cercata	Omeomorfismo, omotopia (stabile a rumore e deformazioni)	Diffeomorfismo, isometria, covarianza su varietà
Interfaccia con l’Algebra Lineare	Operatori di bordo $\partial_k$, Laplaciani combinatori di Hodge	Forme differenziali, tensori metrici, spazi tangenti $T_pM$
Scala tipica	Reti complesse, connettomica su larga scala	Singolo neurone / microcircuiti locali / geometrie di popolazione

Il verdetto operativo

• Se il focus è la struttura della rete e l’elaborazione distribuita del segnale (come configurazioni di connettività generano computazione robusta al rumore), lo strumento d’elezione è la Topologia Algebrica / TDA, supportata dall’algebra lineare dei complessi simpliciali.
• Se il focus è la fisica del potenziale d’azione, l’apprendimento su manifold statistici o le mappe percettive continue, il formalismo irrinunciabile è la Geometria Differenziale.

La convergenza più elegante oggi avviene nella teoria spettrale su fasci cellulari (cellular sheaves), dove l’algebra lineare funge da ponte: spazi vettoriali locali associati a nodi e simplessi che trasportano dati continui su una base topologica discreta.

I calcoli per l’algebra lineare sono in sostanza matriciali, ma con matrici mxn e quindi non invertibili. Qual è il senso?

Il senso risiede nel fatto che le matrici m×n non rappresentano trasformazioni reversibili (automorfismi), ma operatori di confine differenziali, proiezioni dimensionali e relazioni di incidenza tra spazi vettoriali di dimensioni diverse.

In una rete neuronale o in un complesso simpliciale, lo scopo del calcolo matriciale non è quasi mai calcolare un’inversa A−1, bensì decomporre gli spazi di stato, calcolare invarianti topologici e risolvere problemi ai minimi quadrati.

1. Operatori di Bordo e Omologia (Algebra Lineare Topologica)

Nel calcolo topologico, la matrice fondamentale è l’operatore di bordo discreto ∂k​. Questa matrice mappa lo spazio vettoriale dei simplessi k-dimensionali nello spazio vettoriale dei simplessi (k−1)-dimensionali:

∂k​:Ck​→Ck−1​

Se la rete ha n simplessi di dimensione k e m simplessi di dimensione k−1, la matrice di incidenza associata ha dimensioni m×n:

• Esempio sui grafi (k=1): La matrice di incidenza vertice-arco B ha dimensione ∣V∣×∣E∣. Tranne nel caso banale di un grafo ad anello puro, ∣V∣=∣E∣.
• Non invertibilità strutturale: Il fatto cruciale è che:∂k​∘∂k+1​=0⟺Bk​Bk+1​=0Questa composizione è identicamente nulla. Se le matrici fossero quadrate e invertibili, moltiplicando per l’inversa si otterrebbe la matrice nulla o l’identità, distruggendo l’intera struttura algebrica.
• Cosa si calcola senza inversa?
  • Il nucleo ker(∂k​) (i cicli, stati di attivazione ad anello chiuso senza dispersioni).
  • L’immagine im(∂k+1​) (i bordi, cicli banali che delimitano simplessi superiori).
  • Il rango e la nullità tramite eliminazione gaussiana o forma normale di Smith (calcolo dei numeri di Betti):βk​=dim(ker∂k​)−dim(im∂k+1​)=null(∂k​)−rank(∂k+1​)

2. Generazione di Operatori Simmetrici e Laplaciani

Anche se le matrici di partenza B sono rettangolari (m×n), vengono impiegate per generare operatori quadrati, simmetrici e semidefiniti positivi moltiplicandole per la loro trasposta.

Prendendo la matrice di incidenza rettangolare B∈R∣V∣×∣E∣:

L0​=BBT∈R∣V∣×∣V∣

L0​ è il Laplaciano di grafo standard.

Per i complessi simpliciali superiori, il Laplaciano di Hodge k-dimensionale è la somma di due prodotti di matrici rettangolari:

Lk​=BkT​Bk​+Bk+1​Bk+1T​∈Rnk​×nk​

Le matrici rettangolari Bk​ fungono da operatori differenziali discreti (gradiente discreto, rotore discreto, divergenza discreta); i prodotti BBT e BTB sintetizzano l’interazione bidirezionale lungo la catena topologica.

3. Pseudoinverse e Minimi Quadrati nei Sistemi Dinamici

Nelle simulazioni dinamiche (ad esempio nella regressione del connettoma o nella decodifica delle scariche neuronali), ci si trova costantemente di fronte a sistemi lineari Ax=b dove A∈Rm×n:

• m>n (sovradeterminato): Più vincoli temporali o misurazioni rispetto alle variabili di stato. Non esiste una soluzione esatta; il senso dell’algebra lineare è trovare la proiezione ortogonale che minimizza l’errore quadratico ∥Ax−b∥22​:x∗=(ATA)−1ATb=A+b
• m<n (sottodeterminato): Popolazioni neuronali massicce proiettate su un segnale a bassa dimensionalità. Esistono infinite configurazioni di scarica che producono lo stesso output; si cerca la configurazione ad energia minima ∥x∥2​:x∗=AT(AAT)−1b=A+b

La pseudoinversa di Moore-Penrose (A+) e la decomposizione a valori singolari (SVD) A=UΣVT estendono l’invertibilità a qualsiasi operatore m×n, separando le componenti proiettate dal rumore nullo.

In sintesi

Una matrice quadrata invertibile descrive una rotazione o dilatazione rigida di uno spazio in sé stesso. Una matrice rettangolare m×n descrive:

1. Un cambio di dimensione topologica: Mappare bordi a vertici, o volumi a superfici.
2. Un collo di bottiglia dell’informazione: Comprimere stati ad alta dimensione o ricostruire dinamiche da poche osservazioni.
3. La struttura di incidenza: Chi è connesso a chi, senza vincolare a priori il numero di nodi ad eguagliare il numero di connessioni.

Ed invece il calcolo di geometria differenziale e topologia come si rendono programmaticamente?

Programmaticamente la topologia algebrica e la geometria differenziale si implementano con paradigmi computazionali profondamente diversi: algebra lineare esatta su interi/campi finiti e grafi combinatori per la prima; integrazione numerica di ODE/PDE, tensori e differenziazione automatica per la seconda.

1. Calcolo Topologico: Il Dominio Combinatorio e Spettrale

In topologia algebrica e TDA non si integrano derivate: si costruiscono complessi simpliciali a partire da nuvole di punti (o grafi di connettività) e si riducono matrici di incidenza sparse su campi discreti (spesso Z2​, dove l’aritmetica è uno XOR bit a bit).

I mattoni dell’algoritmo:

• Costruzione del complesso: Dato un insieme di spike train o correlazioni, si costruisce una filtrazione (es. Vietoris-Rips): ad ogni soglia di distanza ϵ, se la distanza tra k+1 punti è inferiore a ϵ, si aggiunge un simplesso k-dimensionale.
• Riduzione a scala della matrice di bordo: Per calcolare i diagrammi di persistenza (intervalli di nascita/morte di buchi k-dimensionali), si memorizza la matrice di bordo ∂ come array di colonne sparse. L’algoritmo esegue un’eliminazione gaussiana modificata da sinistra a destra (operando in Z2​):

Python

import numpy as np
import scipy.sparse as sp

# Esempio: Operatore di bordo vettorializzato \partial_1 (Archi -> Vertici)
# Data una lista di vertici (0..N-1) e una lista di archi orientati (u, v)
def build_boundary_1(num_vertices, edges):
    num_edges = len(edges)
    rows = []
    cols = []
    data = []
    
    for edge_idx, (u, v) in enumerate(edges):
        # \partial_1([u, v]) = v - u
        rows.extend([v, u])
        cols.extend([edge_idx, edge_idx])
        data.extend([1.0, -1.0])
        
    # Matrice sparsa (m x n) = (|V| x |E|)
    B1 = sp.csr_matrix((data, (rows, cols)), shape=(num_vertices, num_edges))
    return B1

# Calcolo del Laplaciano combinatorio di Hodge a 0-forme: L0 = B1 @ B1.T
edges = [(0, 1), (1, 2), (2, 0), (2, 3)]
B1 = build_boundary_1(num_vertices=4, edges=edges)
L0 = B1 @ B1.T  # Matrice quadrata 4x4 simmetrica

# Il numero di componenti connesse (Betti 0) corrisponde alla molteplicità dell'autovalore 0
eigenvalues = np.linalg.eigvalsh(L0.toarray())
betti_0 = np.sum(np.isclose(eigenvalues, 0, atol=1e-5))

• Librerie di riferimento: Gudhi (C++/Python), Ripser (C++/Python, ottimizzato con strutture bitset per le riduzioni matriciali), Giotto-tda.

2. Geometria Differenziale: Il Dominio Tensoriale e Continuo

Nella geometria differenziale (varietà dello spazio delle fasi neuronale o geometria dell’informazione) si rappresentano varietà riemanniane calcolando tensori metrici, simboli di Christoffel e geodetiche.

Il cardine moderno dell’implementazione è la differenziazione automatica (AD): invece di derivare a mano metriche e connessioni, si definisce l’energia o la log-verosimiglianza e si sfrutta il calcolo dei vettori Jacobiani/Hessiani tramite framework differenziabili (PyTorch, JAX).

I mattoni dell’algoritmo:

• Calcolo della metrica Riemanniana: Calcolo dell’Hessiano dell’energia o della matrice di informazione di Fisher tramite Jacobiani delle distribuzioni.
• Integrazione dell’equazione delle geodetiche: L’equazione geodetica è una ODE del secondo ordine:x¨k+i,j∑​Γijk​x˙ix˙j=0I simboli di Christoffel Γijk​ vengono ricavati invertendo localmente il tensore metrico gij​(x) e calcolandone i gradienti parziali.

Python

import torch

def fisher_metric_tensor(theta):
    """
    Esempio in geometria dell'informazione:
    Dato un modello neuronale probabilistico parametrizzato da theta,
    la metrica è l'atteso dell'Hessiano della divergenza (o Fisher Information).
    """
    theta = theta.clone().detach().requires_grad_(True)
    
    # Esempio: potenziale convesso \psi(theta) la cui derivata seconda dà la metrica
    psi = 0.5 * torch.sum(theta**2) + torch.log(1 + torch.exp(theta[0]))
    
    # Gradiente prima derivata
    grad_psi = torch.autograd.grad(psi, theta, create_graph=True)[0]
    
    # Hessiano = Tensore metrico Riemanniano g_ij
    dim = theta.shape[0]
    g = torch.zeros((dim, dim))
    for i in range(dim):
        grad2 = torch.autograd.grad(grad_psi[i], theta, retain_graph=True)[0]
        g[i] = grad2
        
    return g

# Passo di ottimizzazione su varietà lungo il gradiente naturale:
# Delta_theta = - g(theta)^{-1} * grad_loss
theta = torch.tensor([1.2, -0.5])
g = fisher_metric_tensor(theta)
grad_loss = torch.tensor([0.1, -0.2])

# Risoluzione del sistema lineare invece di invertire esplicitamente: g * v = grad_loss
natural_gradient = torch.linalg.solve(g, grad_loss)

• Librerie di riferimento: Geomstats (Python, orientata al calcolo su varietà riemanniane, gruppi di Lie e calcolo di geodetiche con backend JAX/PyTorch), SciPy (scipy.integrate.solve_ivp per integrazione Runge-Kutta su varietà invarianti), JAX (calcolo ad alte prestazioni di tensori e jacobiani).

Confronto di Implementazione

Dimensione	Topologia Computazionale	Geometria Differenziale
Primitiva fondamentale	Liste di adiacenza, complessi simpliciali, matrici binarie	Tensori multidimensionali, campi scalari/vettoriali
Aritmetica	Esatta, modulare (Z2​,Zp​) o spettrale su matrici sparse	Floating point continuo (32/64 bit), integratori simplettici/Runge-Kutta
Operazione core	Riduzione per colonne a scalini, SVD sparsa, kernel di operatori di bordo	Differenziazione automatica (AD), inversione locale di gij​, sistemi ODE
Gestione della scala	Filtrazioni troncate (fino a k-simplessi fissati per evitare esplosione combinatoria)	Discretizzazione a griglia (differenze finite) o shooting methods su geodetiche

1. Algoritmo di Riduzione Matriciale a Colonne per l’Omologia Persistente

L’algoritmo classico di Edelsbrunner-Letscher-Zomorodian calcola la persistenza filtrata riducendo la matrice di bordo ∂ tramite operazioni elementari per colonna su Z2​ (operazioni XOR).

Sia data una sequenza filtrata di simplessi ordinati per tempo di comparsa σ1​,σ2​,…,σm​. La matrice di bordo D ha dimensioni m×m, dove:

Di,j​={10​se σi​ eˋ una faccia di codimensione 1 di σj​altrimenti​

Definiamo low(j) come l’indice di riga dell’ultimo 1 nella colonna j:

low(j)=max{i∣Di,j​=1}

Se la colonna j è nulla, low(j) non è definito.

Una matrice è ridotta se l’applicazione low è iniettiva sulle colonne non nulle (ossia non esistono due colonne distinte con lo stesso pivot inferiore).

Python

import numpy as np

def reduce_boundary_matrix(D):
    """
    Riduzione a colonne su Z_2 per una matrice di bordo filtrata D (m x m).
    Ritorna la matrice ridotta R e le coppie di persistenza (nascita, morte).
    """
    m = D.shape[0]
    R = D.copy().astype(np.uint8)
    
    # Memorizza quale colonna possiede un dato pivot: pivot_to_col[row_idx] = col_idx
    pivot_to_col = {}
    persistence_pairs = []

    def get_low(col_idx):
        indices = np.where(R[:, col_idx] == 1)[0]
        return indices[-1] if len(indices) > 0 else None

    for j in range(m):
        low_j = get_low(j)
        # Finché un'altra colonna a sinistra ha lo stesso pivot inferiore, eliminalo con XOR
        while low_j is not None and low_j in pivot_to_col:
            k = pivot_to_col[low_j]
            R[:, j] = (R[:, j] ^ R[:, k])
            low_j = get_low(j)
            
        if low_j is not None:
            pivot_to_col[low_j] = j
            # Il simplesso 'low_j' nasce al tempo low_j e muore al tempo j
            persistence_pairs.append((low_j, j))
        else:
            # Colonna azzerata: il ciclo nato al tempo j persiste all'infinito (non muore)
            pass

    # I simplessi che non sono mai comparsi come morte né sono stati eliminati
    destroyed = set(pivot_to_col.keys())
    killers = set(pivot_to_col.values())
    essential = [(i, np.inf) for i in range(m) if i not in destroyed and i not in killers]

    return R, persistence_pairs + essential

# Esempio: Triangolo vuoto che si chiude
# Simplessi: 0,1,2 (vertici), 3=(0,1), 4=(1,2), 5=(0,2), 6=(0,1,2)
#         0  1  2  3  4  5  6
D_ex = np.array([
    [0, 0, 0, 1, 0, 1, 0], # 0
    [0, 0, 0, 1, 1, 0, 0], # 1
    [0, 0, 0, 0, 1, 1, 0], # 2
    [0, 0, 0, 0, 0, 0, 1], # (0,1)
    [0, 0, 0, 0, 0, 0, 1], # (1,2)
    [0, 0, 0, 0, 0, 0, 1], # (0,2)
    [0, 0, 0, 0, 0, 0, 0]  # (0,1,2)
], dtype=np.uint8)

R, pairs = reduce_boundary_matrix(D_ex)
print("Coppie di persistenza (nascita, morte):", pairs)

2. Calcolo Numerico di Distanza Geodetica e Logaritmo con Geomstats

Sulle varietà Riemanniane, la mappa logaritmica Logp​(q) restituisce il vettore tangente v∈Tp​M tale per cui la geodetica che parte da p con velocità iniziale v raggiunge q al tempo unitario: Expp​(v)=q. La distanza riemanniana è la norma di tale vettore: d(p,q)=∥v∥p​=gp​(v,v)​.

In neuroscienze computazionali, la varietà delle matrici simmetriche definite positive SPD(n) con metrica affine-invariante modella la variabilità delle matrici di covarianza estratte da registrazioni di popolazione (LFP/EEG/spike correlations):

Python

import geomstats.backend as gs
from geomstats.geometry.spd_matrices import SPDMatrices

# Spazio delle matrici simmetriche definite positive 2x2 (es. covarianze neuronali)
n = 2
spd_space = SPDMatrices(n=n)
metric = spd_space.metric  # Metrica Riemanniana affine-invariante

# Due matrici di correlazione/covarianza neuronale
A = gs.array([[2.0, 0.5], 
              [0.5, 1.5]])
B = gs.array([[0.8, -0.3], 
              [-0.3, 1.2]])

# 1. Distanza geodetica Riemanniana intrinseca
dist = metric.dist(A, B)

# 2. Mappa Logaritmica: vettore nello spazio tangente T_A(SPD)
log_vector = metric.log(point=B, base_point=A)

# 3. Mappa Esponenziale: ricostruzione del punto B a partire dal vettore tangente
reconstructed_B = metric.exp(tangent_vec=log_vector, base_point=A)

# 4. Interpolazione lungo la geodetica a t = 0.5
midpoint = metric.geodesic(initial_point=A, end_point=B)(0.5)

print(f"Distanza Riemanniana tra popolazioni: {dist:.4f}")
print("Vettore tangente Log_A(B):\n", log_vector)
print("Punto medio lungo la geodetica (t=0.5):\n", midpoint)

3. Dinamica di FitzHugh-Nagumo nello Spazio delle Fasi

Il sistema di FitzHugh-Nagumo è una riduzione bidimensionale a scala lenta-veloce del modello di Hodgkin-Huxley:

⎩⎨⎧​v˙=v−3v3​−w+Iext​w˙=ϵ(v+a−bw)​(variabile veloce: potenziale di membrana)(variabile lenta: variabile di recupero)​

Le nullcline separano le regioni di segno delle derivate:

• Nullcline di v (v˙=0): la curva cubica w=v−3v3​+Iext​
• Nullcline di w (w˙=0): la retta w=bv+a​

L’intersezione tra le nullcline definisce il punto fisso; una biforcazione di Hopf determina la transizione tra lo stato quiescente eccitabile e il ciclo limite (treno di potenziali d’azione):

Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parametri del modello
a = 0.7
b = 0.8
eps = 0.08
I_ext = 0.5  # Corrente applicata sufficiente a innescare il firing periodico

def fhn_system(t, y):
    v, w = y
    dvdt = v - (v**3) / 3.0 - w + I_ext
    dwdt = eps * (v + a - b * w)
    return [dvdt, dwdt]

# 1. Integrazione numerica della traiettoria (Runge-Kutta 4(5))
t_span = (0, 100)
t_eval = np.linspace(t_span[0], t_span[1], 2000)
y0 = [-1.0, 0.0]  # Condizione iniziale di riposo

sol = solve_ivp(fhn_system, t_span, y0, t_eval=t_eval, rtol=1e-8, atol=1e-8)

# 2. Griglia per il campo vettoriale nello spazio delle fasi
v_grid = np.linspace(-2.5, 2.5, 25)
w_grid = np.linspace(-1.0, 2.0, 25)
V, W = np.meshgrid(v_grid, w_grid)

DV = V - (V**3)/3.0 - W + I_ext
DW = eps * (V + a - b * W)
Speed = np.sqrt(DV**2 + DW**2)
DV_norm = DV / (Speed + 1e-6)
DW_norm = DW / (Speed + 1e-6)

# 3. Calcolo delle Nullcline
v_dense = np.linspace(-2.5, 2.5, 300)
v_nullcline = v_dense - (v_dense**3)/3.0 + I_ext
w_nullcline = (v_dense + a) / b

# Tracciamento
plt.figure(figsize=(9, 6))
plt.quiver(V, W, DV_norm, DW_norm, color='lightgray', angles='xy', scale=30, width=0.003)
plt.plot(v_dense, v_nullcline, 'r--', label=r'Nullcline $v$ ($\dot{v}=0$, cubica)')
plt.plot(v_dense, w_nullcline, 'b--', label=r'Nullcline $w$ ($\dot{w}=0$, retta)')
plt.plot(sol.y[0], sol.y[1], 'k-', linewidth=2, label='Traiettoria (Ciclo Limite)')
plt.plot(y0[0], y0[1], 'go', markersize=8, label='Condizione iniziale')

plt.xlim(-2.5, 2.5)
plt.ylim(-0.8, 1.8)
plt.xlabel('Potenziale di Membrana $v$')
plt.ylabel('Variabile di Recupero $w$')
plt.title('Spazio delle Fasi FitzHugh-Nagumo: Ciclo Limite e Nullcline')
plt.grid(True, linestyle=':')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()

Concludo osservando che nei miei studi universitari non avrei mai pensato che corsi di Istituzioni di Algebra Superiore e Geometria Superiore potessero essere così estesamente applicati …

È una sensazione comune a molti matematici: per decenni corsi come Istituzioni di Geometria Superiore (varietà differenziabili, fibrati, connessioni, coomologia di de Rham) o Algebra Superiore (teoria dei moduli, complessi di catene, algebra omologica, funtori Ext e Tor) sono stati vissuti e insegnati come vette di pura astrazione teorica, lontanissime da qualsiasi calcolatore o applicazione empirica.

La svolta degli ultimi vent’anni è stata duplice:

1. La discretizzazione dell’astratto: Concetti che sembravano richiedere per forza il continuo liscio o categorie astratte hanno trovato corrispettivi combinatori esatti e implementabili su silicio. La coomologia è diventata algebra lineare sparsa su complessi simpliciali; i fibrati sono diventati cellular sheaves per il consenso distribuito e la sincronizzazione di sensori; la teoria dei gruppi di Lie è diventata la spina dorsale delle reti neurali geometriche (Equivariant Deep Learning), dove l’invarianza a rototraslazioni viene imposta algebricamente nei pesi sinaptici invece di essere appresa per tentativi.
2. Il cambio di paradigma nei dati: Finché i dati erano tabelle o segnali unidimensionali, bastavano l’analisi di Fourier classica e la statistica gaussiana. Nel momento in cui ci si è trovati a dover modellare grafi dinamici con milioni di nodi (connettomica), distribuzioni di probabilità non-euclidee (geometria dell’informazione) e spazi di stati ad altissima dimensione ma concentrati su varietà a bassa dimensione (manifold hypothesis), l’analisi standard si è rivelata cieca alle proprietà globali. L’unico linguaggio in grado di parlare di “forma”, “buchi”, “invarianza per deformazione” e “curvatura intrinseca” era proprio quello forgiato dalla topologia algebrica e dalla geometria differenziale.

Vedere un funtore di omologia ridotto a un algoritmo di pivoting matriciale su bitset con operatori XOR, o i simboli di Christoffel calcolati via differenziazione automatica per guidare l’apprendimento su una varietà di Riemann, chiude un cerchio straordinario: l’astrazione non era una fuga dal mondo reale, ma l’unico telaio abbastanza generale da poterlo descrivere quando la complessità dei dati fosse finalmente esplosa.

La parte realmente incredibile è come l’astratto praticamente quasi non immaginabile diventi veramente quasi concreto nelle sorgenti python.

C’è una sorta di “vertigine epistemologica” in questo passaggio: l’oggetto che sui testi universitari richiedeva tre pagine di dimostrazioni algebriche tra diagrammi commutativi e classi di equivalenza, nel codice diventa un array NumPy, un puntatore o una maschera di bit.

Il punto è che l’astrazione pura non è fumo metafisico, ma compressione strutturale estrema. Quando la si implementa, si compie il percorso inverso: la si decomprime in passaggi operativi elementari.

Questo collasso tra l’inimmaginabile e il tangibile si vede plasticamente in almeno tre punti:

• Il quoziente algebrico diventa un filtro di indici: Definire uno spazio vettoriale quoziente V/W o un gruppo di omologia Hk​=ker∂k​/im∂k+1​ su carta impone di pensare a “insiemi di insiemi” (classi di equivalenza) con la relazione di congruenza. Nel codice, fare il quoziente non significa manipolare entità incorporee: significa semplicemente trovare una base ortogonale per il sottospazio W tramite SVD o QR, e proiettare sul complemento ortogonale. Lo spazio quoziente si riduce a un’operazione di null_space e qualche prodotto scalare.
• ∂∘∂=0 diventa un’asserzione booleana: Il pilastro della topologia algebrica e delle forme differenziali (d2=0) su carta sembra un assioma geometrico solenne. Nello script Python sono due matrici sparse di incidenza B1​ e B2​ che verifichi al volo con:Pythonassert (B1 @ B2).nnz == 0 Quella singola riga di codice garantisce che “il bordo di un bordo è vuoto”, e da lì discende istantaneamente l’intera sequenza esatta.
• Il fibrato e lo spazio tangente diventano un grafo di computazione: Pensare allo spazio tangente Tp​M su una varietà astratta richiede un notevole sforzo visivo ed euristico (classi di equivalenza di curve tangenti o derivazioni locali). Con la differenziazione automatica di PyTorch o JAX, il fibrato tangente è fisicamente il forward tape del calcolo automatico, e il fibrato cotangente (le 1-forme differenziali) è il backward pass dei gradienti. La dualità tra vettori e co-vettori è la dualità tra derivazione in avanti e propagazione all’indietro.

Quando scrivi il codice, ti accorgi che il computer non ha alcuna “intuizione geometrica”, eppure fa girare la topologia. Perché l’algebra superiore e la geometria differenziale moderna, spogliate dalla retorica didattica, sono nate esattamente come sintassi logica rigida: sistemi di regole di riscrittura, linearizzazioni locali e contabilità di invarianti.

Python non fa altro che prenderle alla lettera, eseguendo l’algebra senza l’ansia dell’intuizione spaziale e restituendo, alla fine, un diagramma o una traiettoria che possiamo finalmente guardare sullo schermo.

Ma le librerie sottostanti penso siano monumentali …

Dietro l’eleganza di due o tre righe Python c’è una delle stratificazioni ingegneristiche più imponenti mai costruite. Python in questi ambiti è quasi solo una sintassi di orchestrazione; la computazione reale poggia su milioni di righe di C, C++, Fortran e CUDA assemblate nell’arco di mezzo secolo.

La complessità si distribuisce su strati nettamente distinti a seconda che ci si muova nella geometria differenziale o nella topologia algebrica computazionale.

1. Lo strato numerico e tensoriale: decenni di micro-ottimizzazioni hardware

Quando calcoli una metrica Riemanniana con Geomstats, JAX o PyTorch, il codice attraversa una gerarchia profondissima:

• BLAS / LAPACK (Fortran / C): Alla base di tutto ci sono librerie storiche (OpenBLAS, Intel MKL) che implementano le routine di algebra lineare densa (DGEMM, DPOTRF per Cholesky, decomposizioni SVD). Qui ogni ciclo di clock della CPU è ottimizzato per registri vettoriali AVX-512, allineamento di memoria nelle cache L1/L2/L3 e prefetching predittivo.
• Motori di Differenziazione Automatica (C++ / LLVM):
  • In PyTorch c’è libtorch, un nucleo C++ massiccio che gestisce grafi dinamici (autograd) tracciando ogni operatore elementare in nodi di un DAG (grafo aciclico diretto) con i relativi gradienti analitici implementati a mano.
  • In JAX, il codice Python viene trasformato in una rappresentazione intermedia funzionale (jaxpr), passata a XLA (Accelerated Linear Algebra), un compilatore che fonde le operazioni algebriche (kernel fusion) per evitare scritture ridondanti in RAM e compila codice macchina su misura per GPU o TPU via LLVM.
• Integratori Simplettici e ODE: Le librerie come SciPy non usano Python per integrare: richiamano codice C/Fortran derivato da pacchetti come ODEPACK (lsoda) o implementazioni C++ di dormand-prince a passo adattivo (odeint in Boost), capaci di gestire matrici di rigidità (stiff equations) invertendo matrici Jacobiane a ogni micro-intervallo temporale.

2. Lo strato topologico: combinatoria e strutture dati non convenzionali

Nella topologia algebrica computazionale il problema è opposto: l’algebra lineare standard fallisce per ragioni di scala. Un complesso di Vietoris-Rips su mille punti può generare decine di milioni di simplessi in poche dimensioni; una matrice densa manderebbe istantaneamente la RAM in out-of-memory.

Dietro librerie come Gudhi (sviluppata dall’INRIA) e Ripser c’è un’ingegneria algoritmica estrema:

• Gudhi e CGAL (C++ con template spinti): Sfrutta alberi di simplessi (simplex trees), strutture a prefisso (trie) altamente compresse dove ogni nodo dell’albero rappresenta una faccia, consentendo di interrogare e memorizzare complessi di dimensioni elevate con overhead minimo di puntatori.
• Ripser (l’apice dell’ottimizzazione in C++): Il codice di Ulrich Bauer (Ripser) è un capolavoro di economia computazionale. Invece di istanziare l’intera matrice di bordo, Ripser non memorizza la matrice: la genera on-the-fly riga per riga e colonna per colonna, sfruttando la dualità di coomologia (l’omologia persistente si calcola più velocemente lavorando sui cocicli) e riducendo le colonne con rappresentazioni binarie sparse compresse in registri interi a 64 bit con istruzioni bitwise hardware (popcnt, xor).

La piramide software

L’architettura può essere visualizzata come una piramide a strati concentrici:

[ Python Script ]          geomstats / gudhi / ripser / torch (poche righe dichiarative)
       │
[ Binding Layer ]          pybind11 / Cython / CPython C-API
       │
[ Algoritmi Core ]         C++ Template Metaprogramming (CGAL, Simplex Trees, Autograd Engine)
       │
[ Compilatori & JIT ]      XLA, TorchInductor, LLVM
       │
[ Calcolo Numerico ]       LAPACK, OpenBLAS, SuiteSparse, cuBLAS/cuSOLVER (CUDA)
       │
[ Hardware ]               SIMD (AVX2/AVX-512), Tensor Cores GPU, gerarchie di Cache

Quello che per il matematico è un funtore o un tensore di curvatura, per questa infrastruttura è un problema di allocazione di blocchi di memoria contigui, vettorizzazione su registri a 512 bit e parallelizzazione massiva. La vera prodezza non è solo matematica, ma la capacità di questi strati di nascondere tutta la frizione fisica del silicio dietro astrazioni pulite.

Milioni di righe di C, C++, Fortran e CUDA! Ma ci rendiamo conto! E’ incredibile, rivoluzionario, poco tempo fa non era concepibile neppure nell’ambito dei più puri e vertiginosi racconti di fantascienza!

È uno shock prospettico del tutto giustificato. Nemmeno la fantascienza della Golden Age — da Asimov a Dick, che pure immaginavano cervelli positronici o macchine oracolari — era riuscita a concepire la natura di questo prodigio.

Negli scenari fantascientifici classici, il supercalcolatore era una scatola nera monolitica scaturita dal lampo di genio di un singolo inventore o da una fisica esotica sconosciuta. Nessuno aveva previsto che la vera rivoluzione sarebbe stata una cattedrale gotica di sedimentazione intellettuale collettiva: generazioni di matematici, fisici, ingegneri e programmatori che, spesso senza conoscersi e a distanza di quarant’anni, hanno posato ciascuno un blocco di calce e silicio.

C’è qualcosa di vertiginoso nel tracciare cosa accade in una frazione di secondo sulla macchina:

• Sfasamento temporale sincronizzato: Un’istruzione scritta oggi pomeriggio evoca subroutine scritte in Fortran 77 a Oak Ridge o nei laboratori Bell negli anni ’70 per la decomposizione ai valori singolari, le fa dialogare con codice C++20 a metaprogrammazione spinta per generare alberi simpliciali, e scarica il calcolo su migliaia di stream paralleli CUDA gestiti da un silicio che esegue moltiplicazioni di matrici a virgola mobile a frequenze di miliardi di cicli al secondo.
• L’invisibilità del mostro: Quella piramide titanica non scricchiola, non espone cavi, non chiede di conoscere la topologia dei die di silicio o la gestione della coerenza di cache tra i core. Si manifesta come un’interfaccia interattiva che restituisce il risultato istantaneamente, nascondendo la più complessa catena di astrazioni logico-meccaniche mai eretta dalla nostra specie.
• L’astratto che trova il suo motore: Per secoli la matematica pura ha corso libera, accumulando strutture (spazi a infinite dimensioni, algebre omologiche, geometrie non commutate) con la serenità di chi credeva di abitare un regno disincarnato. Scoprire che l’universo del calcolo numerico e del silicio è cresciuto fino ad accogliere e dare corpo a queste geometrie vertiginose è il vero punto di rottura: l’astrazione non era un rifugio mentale, ma l’anticipo formale degli strumenti computazionali del presente.

Ricordo che negli anni 80 occasionalmente frequentavo il dipartimento di Fisica di Padova e l’INFN. All’epoca mi chiedevano spesso codice Fortran. Se penso che sorgenti Fortran possano essere sopravvissute da tale preistoria finisco nell’oblio …

Eppure è esattamente così: non solo sono sopravvissute, ma governano ancora l’infrastruttura critica del calcolo scientifico globale.

Se negli anni ’80 a Padova e all’INFN scrivevi o vedevi compilare schede perforate e sorgenti Fortran 77 su mainframe VAX o CDC, stavi assistendo alla forgiatura delle fondamenta su cui poggia l’intelligenza artificiale moderna. Non è un’iperbole archeologica: è un fatto verificabile riga per riga nei sorgenti che scarichiamo oggi.

• Linpack ed Eispack non sono mai morti: I pacchetti scritti negli anni ’70 da Jack Dongarra e colleghi dell’Argonne National Laboratory per risolvere sistemi lineari e problemi agli autovalori sono confluiti all’inizio degli anni ’90 in LAPACK (Linear Algebra Package), scritto in rigoroso Fortran 77.
• Dentro SciPy e NumPy batte ancora quel cuore: Quando oggi fai un’operazione banale come scipy.linalg.eig o una decomposizione svd, sotto il cofano CPython invoca interfacce generate da f2py (Fortran to Python) che collegano direttamente binari compilati dal codice originale di LAPACK, BLAS e pacchetti di integrazione ODE come ODEPACK (lsode e lsoda, scritti da Alan Hindmarsh ai Lawrence Livermore National Laboratory tra gli anni ’70 e ’80).
• Perché nessuno li ha riscritti da zero in linguaggi moderni? Non solo per inerzia o rispetto filologico, ma per due ragioni squisitamente tecniche:
  1. Correttezza numerica collaudata da miliardi di ore di calcolo: Quelle routine hanno affrontato ogni possibile caso patologico di instabilità numerica, overflow, underflow e cancellazione disastrosa. Riscrivere 500.000 righe di algebra lineare avanzata in C++ o Rust introdurrebbe sottili bug numerici che hanno richiesto quarant’anni per essere estirpati.
  2. Assenza di aliasing nei puntatori: Il Fortran, per specifica, vieta l’aliasing dei puntatori negli array passati alle subroutine (due variabili puntatore non possono puntare alla stessa cella di memoria). I compilatori moderni (come gfortran o i compilatori Intel LLVM) possono srotolare i cicli, vettorizzare su registri AVX e riordinare le istruzioni con un’aggressività e una sicurezza matematica che il C/C++ fatica ancora a eguagliare senza annotazioni esplicite tipo __restrict.

Quindi non c’è nessun oblio, semmai un’incredibile continuità: la fisica computazionale dell’INFN degli anni ’80, con i suoi array dimensionati staticamente, le istruzioni DO e la precisione a 64 bit dei numeri reali, ha semplicemente mutato pelle all’esterno, mentre il suo scheletro continua a macinare vettori e autovalori nei data center del pianeta.

Ed il vecchio Pascal della Teoria della Programmazione delle Macchine Calcolatrici del dipartimento di Matematica di Padova è servito a qualche cosa?

È servito a molto più di quanto non appaia a prima vista. Mentre a Fisica e all’INFN il Fortran regnava per pura forza bruta numerica — calcolare sezioni d’urto o traiettorie senza troppi complimenti sulla forma —, a Matematica il Pascal rappresentava una vera e propria rivoluzione epistemologica.

Quel corso a Padova non insegnava semplicemente una sintassi per far muovere bit, ma introduceva l’idea, all’epoca radicale, che la programmazione fosse una branca della logica matematica e della teoria della dimostrazione.

Il Pascal ha lasciato un’eredità diretta e profonda che modella l’informatica moderna in almeno quattro punti cardine:

1. La fondazione del “Type System” moderno

Prima di Niklaus Wirth e del Pascal (e di Algol 68), i linguaggi erano una distesa pragmatica di parole di memoria, registri, salti incondizionati (GOTO) e blocchi continui di numeri. Il Pascal ha formalizzato la nozione che un tipo di dato è un dominio matematico dotato di un insieme chiuso di operazioni lecite. I tipi record, gli intervalli (1..10), i tipi enumerati e i puntatori tipizzati non servivano per comodità: servivano affinché il compilatore potesse agire come un verificatore di teoremi. Se sommavi due entità incompatibili, il compilatore rifiutava il programma a priori. Tutta la moderna sicurezza dei tipi (da Rust a TypeScript, fino ai type checker formali di Python come Mypy) discende direttamente da quell’impostazione concettuale.

2. La corrispondenza di Curry-Howard e la correttezza formale

In quegli anni, nei corsi di teoria della programmazione a Matematica, il Pascal andava a braccetto con la logica di Hoare, le pre-condizioni, le post-condizioni e gli invarianti di ciclo (“Program Construction and Verification” di Dijkstra e Wirth). L’idea era: un programma non si prova testandolo con qualche input, si dimostra corretto come un teorema. Quella scuola di pensiero è il pilastro teorico che oggi governa la verifica formale del software critico (compilatori verificati come CompCert, kernel provati formalmente come seL4, e i moderni proof assistant come Coq e Lean).

3. L’architettura del software moderno: Moduli e P-Code

Il Pascal ha generato la sua evoluzione naturale, Modula-2 (sempre di Wirth), da cui è nata la programmazione modulare moderna: interfacce pulite separate dall’implementazione, incapsulamento rigoroso e gestione dell’astrazione. Inoltre, per far girare il Pascal sulle architetture eterogenee dell’epoca, Wirth e il suo team inventarono la P-Machine (Pascal-P) e il P-Code: un bytecode intermedio eseguito da una macchina virtuale stack-based portabile. Quel concetto è esattamente l’antenato architetturale della Java Virtual Machine (JVM), del CLR di .NET e del bytecode compilato da CPython stesso (.pyc).

4. La disciplina mentale del matematico

Per chi si è formato su quei banchi, il Pascal ha svolto una funzione pedagogica insostituibile: ha disabituato al “codice spaghetti” del Fortran d’annata e ha insegnato a concepire gli algoritmi come composizione ordinata di strutture dati e controllo del flusso. Chi ha imparato a progettare algoritmi dimostrando l’invariante di un while o strutturando alberi binari con i puntatori rigorosi del Pascal, quando oggi si trova a scrivere pipeline vettorizzate, a implementare strutture a grafo o a gestire sistemi distribuiti, possiede un’igiene mentale e un rigore architetturale che nessun tutorial moderno riesce a trasmettere.

Il Fortran è stato il muscolo del calcolo numerico; il Pascal di quel corso di Matematica è stato la grammatica della correttezza logica. Senza il primo non avremmo la velocità di LAPACK, ma senza il secondo non avremmo mai avuto l’ingegneria del software né i sistemi operativi moderni.

Ricordo perfettamente le dimostrazioni di correttezza e terminazione degli algoritmi dell’esame citato …

La tripla di Hoare {P}C{Q}, la ricerca ostinata dell’invariante di ciclo I e la funzione di terminazione (o variante) su un insieme ben ordinato: chi ha affrontato quella scuola porta addosso un’impronta indelebile.

All’epoca sembrava quasi un esercizio di pedanteria accademica: prendere un banale algoritmo di partizione di Hoare per il QuickSort o l’algoritmo di Euclide per il MCD e dover dimostrare, riga per riga, che:

1. Invarianza: se la precondizione P è vera all’ingresso, I resta vera a ogni iterazione;
2. Correttezza parziale: se il ciclo termina, I∧¬B⟹Q (la postcondizione desiderata è garantita);
3. Terminazione totale: esiste una funzione variante t:Stato→N strettamente decrescente a ogni passo (tk+1​<tk​), con limite inferiore garantito dal buon ordinamento dei numeri naturali, escludendo a priori il loop infinito.

La bellezza e la vertigine di quell’approccio stavano nel ribaltamento della prospettiva: non si scriveva il codice per poi tentare di capirne il comportamento con dei test, si deduceva il codice a partire dall’invariante. L’invariante non era un commento a margine, era il seme matematico da cui il ciclo scaturiva come conseguenza logica necessaria.

Oggi, in un’epoca in cui gran parte dello sviluppo software procede spesso per tentativi ed errori — compilando, testando e correggendo a colpi di debugger finché “sembra funzionare” —, aver interiorizzato quella disciplina mentale fa tutta la differenza del mondo.

Quando progetti una struttura dati complessa (come l’algoritmo di riduzione a colonne per l’omologia persistente che abbiamo visto prima, dove devi garantire la non-collisione dei pivot inferiori, o la convergenza di un integratore simplettico su una varietà), non stai semplicemente buttando giù istruzioni imperative. Inconsapevolmente, o forse chiaramente, continui a porti le stesse domande di quel corso a Padova: qual è l’invariante di questo loop? Cosa garantisce che l’indice converga al limite?

Quel rigore non è mai andato perduto: è la linea sottile che separa chi scrive script sperando che non si rompano, da chi progetta algoritmi sapendo esattamente perché non possono fallire.