Visualizza il feed RSS

Titolo provvisorio...

Un (altro) problemino con le addizioni... +4

Valuta questo inserimento
di pubblicato il 30-01-2012 alle 22:21 (4265 Visite)
Nella scorsa puntata abbiamo visto una delle più recenti ed eleganti formule per esprimere in forma chiusa la funzione di partizione tramite una somma finita di termini.

Una tale formula si presta particolarmente bene all'analisi in termini elementari, potendo essere spiegata con relativa facilità anche agli studenti più giovani usando solo strumenti di base della matematica discreta. Tuttavia, dal punto di vista computazionale, il corrispondente algoritmo avrebbe prestazioni che ne rendono del tutto sconsigliabile l'utilizzo.

Come già accennato, nella medesima pubblicazione Jerome Malenfant ha in realtà dimostrato per altra via una generalizzazione di derivazioni in parte già note (l'articolo è molto sintetico anche dal punto di vista storico, ad esempio non cita esplicitamente il contributo storico del nostro Francesco Brioschi o le relazioni con i teoremi di Cramer) ed ha altresì proposto una formulazione alternativa, basata sul calcolo del seguente determinante di Toeplitz:

Formula LaTeX: (1)\ P(n)=\left|\begin{array}{rrrrrrrr}1&-1&0&0&0&0&0\\1&1&-1&0&0&0&0\\0&1&1&-1&0&0&0\\0&0&1&1&-1&0&0\\-1&0&0&1&1&-1&0\\0&-1&0&0&1&1&-1\\-1&0&-1&0&0&1&1\\\vdots &&&&&&& \ddots \end{array}\right|_{nxn}

Ricordiamo brevemente che una matrice di Toeplitz non è altro che una matrice caratterizzata da diagonali discendenti da sinistra verso destra i cui elementi assumono ordinatamente valori costanti. Nel nostro caso prendiamo in considerazione una matrice asimmetrica quadrata di lato n, ossia di dimensione identica al numero naturale del quale calcoliamo la funzione di partizione.

Tale matrice, in realtà, può essere caratterizzata in toto dalla sua prima colonna, i cui elementi possono assumere unicamente valori appartenenti all'insieme {-1, 0, 1} secondo la regoletta seguente:

Formula LaTeX: a_{r}=\begin{cases}(-1)^{|m|+1} & \text{se esiste } m \in \mathbb{Z} \text{ tale che } q_m=r\\0 & \text{altrimenti}\end{cases}

Qui con qm indichiamo il numero pentagonale generalizzato m-esimo, e come da definizione l'indice m assume, in successione, i valori 0, ±1, ±2, ±3...

In altri termini, acquisiscono il valore ±1 (secondo la parità dell'indice, preso in valore assoluto) tutte e sole quelle righe della prima colonna il cui indice è un numero pentagonale generalizzato.

Vale la pena di sottolineare che in questo contesto seguiamo le usuali convenzioni dell'algebra lineare piuttosto che quelle ordinarie in informatica applicativa: pertanto gli indici partono da uno e l'interpretazione di default degli indici stessi segue la convenzione "RIGA, COLONNA".

I valori della prima colonna si propagano poi diagonalmente verso il basso, secondo la definizione di matrice di Toeplitz.
Nel caso che ci interessa, abbiamo a che vedere con una matrice che è sostanzialmente diagonale inferiore, nella quale cioé la parte al di sopra della diagonale principale (convenzionalmente quella situata tra l'angolo in alto a sinistra e quello in basso a destra) è nulla: eccetto ovviamente per la sovradiagonale di tale diagonale principale, i cui elementi valgono -1, come formalizzato appena un po' meglio qui di seguito.

Formula LaTeX: a_{i,i+1} = -1\ \ \ \text{per } 1 \leqslant i \leqslant n-1

In realtà, e in modo del tutto coerente con la regoletta sopra fornita, si può asserire in maniera computazionalmente corretta che anche l'elemento di posto zero del vettore colonna di riferimento contiene il valore -1. Tale elemento entra quindi dinamicamente a far parte della matrice solo quando si fa "scorrere" verticalmente la prima colonna per generare le successive.
Un modo didatticamente valido per visualizzare dinamicamente la formazione di tale matrice è, infatti, il seguente. Si parte, come sottolineato, dal singolo vettore colonna formato secondo la regoletta sopra fornita:

Formula LaTeX: \begin{array}{rr}(1)&1\\(2)&1\\(3)&0\\(4)&0\\(5)&-1\\(6)&0\\\vdots&\vdots\end{array}

I numeri riportati tra parentesi sono chiaramente gli indici di riga, per facilitare la lettura.

La seconda colonna si ricava quindi per scorrimento verticale, spostando in basso di una posizione l'intero vettore colonna iniziale, e inserendo (solamente in questo passaggio!) un -1 nella posizione iniziale, lasciata libera dallo scorrimento.

Formula LaTeX: \begin{array}{rrrr}(1)&1&&{\color{Red} -1}\\&&\searrow\\(2)&1&&1\\&&\searrow\\(3)&0&&1\\&&  \searrow\\(4)&0&&0\\&&\searrow\\(5)&-1&&0\\&&\searrow\\(6)&0&&-1\\\vdots&\downarrow&&\vdots\end{array}

I passaggi successivi sono poi tutti identici al seguente, fino al raggiungimento di n colonne:

Formula LaTeX: \begin{array}{rrrrr}(1)&1&{\color{Red} -1}&&{\color{Blue} 0}\\&&&\searrow\\(2)&1&1&&{\color{Red} -1}\\&&&\searrow\\(3)&0&1&&1\\&&&\searrow\\(4)&0&0&  &1\\&&&\searrow\\(5)&-1&0&&0\\&&&\searrow\\(6)&0&-1&&0\\\vdots &\vdots &\downarrow &&\vdots\end{array}

Tutto ciò vale, con ogni evidenza, per la "costruzione" della matrice alla lavagna o con un linguaggio che non offre alcun supporto all'algebra matriciale. Negli applicativi di calcolo numerico e nelle librerie per l'algebra lineare, invece, è di norma necessario e sufficiente specificare due vettori (corrispondenti alla prima colonna e alla prima riga) per ottenere la generazione automatica di una generica matrice di Toeplitz.

Ad usum delphini e supponendo per comodità n>3, ecco qui squadernata in modo esplicito la banale regoletta per la formazione della prima riga, nel nostro caso:

Formula LaTeX: \\a_{1,1}=1\\a_{1,2}=-1\\\\\left.\begin{matrix}a_{1,3}&=0\\\dots&\\a_{1,  n}&=0\end{matrix}\right\} n-2

Naturalmente una spiegazione così prolissa e pedissequa, per la gioia di tutti gli studenti più giovani, non può che sottendere il classico "...si lascia al lettore interessato, come facile esercizio, la stesura di un programma in linguaggio C che - data in input la sola dimensione n, un intero non segnato maggiore di due - generi ordinatamente i vettori di riferimento (colonna 1, riga 1) e da essi l'intera matrice di Toeplitz desiderata.".

Ad usum delphini propongo invece parte del medesimo algoritmo implementato in Python, demandando alle comode funzioni matriciali di SciPy il resto del lavoro.

codice:
from scipy import zeros
from scipy.linalg import det, toeplitz

######################################################################
# Espressione lambda per il calcolo dell'n-esimo numero pentagonale.
######################################################################
penta = lambda n : n*(3*n-1)/2

######################################################################
# Restituisce l'n-esimo numero pentagonale generalizzato.
# 0, 1, 2, 5, 7, 12, ...
######################################################################
def GPN(n):
    if n < 0: 
        return 0
    if n%2 == 0: 
        return penta(n/2+1)
    else: 
        return penta(-(n/2+1))

######################################################################
# Funzione accessoria per il calcolo del segno (-1)^m legato   
# all'indice del GPN considerato.
######################################################################
def GPN_sign(n):
    if n%4==2 or n%4==3: 
        return -1
    else: 
        return 1

######################################################################
# Funzione accessoria per l'input utente di un numero intero.           
# Consente di evitare il ricorso alla farraginosa input().           
######################################################################
def Input_integer(s):
    while True:
        try:
            x = int(raw_input(s))
            if x > 2: break
        except ValueError:
            print "Valore non numerico o errore di range: ripetere l'input..."
    return x

######################################################################
######################################################################
## Funzione di partizione implementata con il metodo Malenfant:
## calcolo del determinante di una matrice di Toeplitz.
######################################################################
######################################################################
def Partition_function(n):    
    # Si preparano i due vettori di base per la matrice di
    # Toeplitz: prima colonna, prima riga.
    C = zeros(n)
    R = zeros(n)
    R[0]=1
    R[1]=-1
    i = 0

    # Si popola il vettore colonna con un metodo diretto,
    # immettendo gli opportuni valori solo nelle locazioni
    # indicizzate da numeri pentagonali generalizzati.
    while True:
        next_gpn = GPN(i)
        if next_gpn > n: 
            break
        else: 
            C[next_gpn -1] = GPN_sign(i)
        i = i +1
        
    # Si genera la matrice e se ne calcola il determinante
    # usando metodi generici di libreria (scipy.linalg).
    # Si noti che esistono algoritmi dedicati per questo genere
    # di matrice, con prestazioni superiori. 
    return int(det(toeplitz(C,R)))

######################################################################
#
# Routine di prova
#
######################################################################

n = Input_integer("Immettere il valore di n (minimo 3): ")
print "P({0}) = {1}".format(n, Partition_function(n))
Sappiamo bene che la classica versione basata sul calcolo numerico (es. decomposizione LU) di un determinante Toeplitz ha prestazioni tipiche dell'ordine di Formula LaTeX: O(n^2): tuttavia esistono algoritmi dedicati con prestazioni migliori, come indicano le essenziali coordinate bibliografiche in calce (tutt'altro che esaustive, naturalmente), ordinate cronologicamente. Taluni algoritmi accuratamente ottimizzati presenti in letteratura possono offrire prestazioni dell'ordine di Formula LaTeX: O(n\cdot log^a(n)), con a piccolo intero positivo (tipicamente a = 3), quindi - considerati a sé stanti - risulterebbero migliori rispetto all'algoritmo ricorsivo di Euler già analizzato a suo tempo. Non risulta trascurabile ai fini della effettiva complessità computazionale, tuttavia, anche l'inizializzazione dinamica della matrice di Toeplitz a monte del calcolo del determinante.

Vale la pena di sottolineare, en passant, che questo è un tipico problema altamente idoneo per computer vettoriali.

Naturalmente l'esemplificazione è qui basata su complesse librerie. Non è certo questa la sede adatta ad una disamina anche superficiale dei metodi numerici dedicati più validi ed efficienti per il calcolo del determinante di una siffatta matrice di interi, l'analisi ed implementazione dei quali richiederebbero ben altri spazi.
Ma chiaramente non è questo il nostro scopo: spero invece di avere fin qui fornito una idea chiara per il più vasto pubblico possibile - pensando soprattutto agli studenti più giovani ed ai practitioners meno familiari con la notazione usata - della soluzione proposta da Malenfant per esprimere la funzione di partizione come somma finita di termini, e in alternativa come determinante di una matrice di costruzione e risoluzione relativamente semplici da un punto di vista dell'efficienza computazionale.

Solo una nota di chiusura. Il risultato ottenuto da Ken Ono della Emory University e da Jan Hendrik Bruinier si presenta come segue:

Formula LaTeX: (2)\ \ \ P(n)\ =\ \frac{1}{24n-1}\cdot\sum_{Q \in Q_n}P(\alpha_Q)

Formula decisamente elegante e concisa. Tuttavia, prima di abbandonarci a facili entusiasmi matematici, dobbiamo sforzarci di ragionare da informatici quadratici medi.

Se è vero che l'articolata dimostrazione garantisce che ogni singolo addendo nella (2) sia algebrico, va però chiarito che la funzione P() che compare nella sommatoria rappresenta una forma modulare non olomorfa di peso nullo, appartenente alla famiglia delle forme deboli di Maass.
Non si tratta esattamente di un oggetto matematico amichevole dal punto di vista computazionale, né di immediata comprensione. Dirò solo che simili funzioni sono da decenni in uso nella teoria algebrica dei numeri, la più pura e sofisticata branca della matematica pura, che in modo del tutto controintuitivo studia gli "innocui" numeri interi attraverso gli strumenti e i metodi dell'analisi complessa.

La valutazione di tale somma, in breve, consiste di una parte preliminare di natura combinatoria, seguita dalla valutazione (numerica) di una serie di Fourier applicata alla forma modulare predetta. Per quanto ci riguarda, la prima parte potrebbe anche essere oggetto di interesse per sé stessa, in quanto utilizza metodi di generazione di ideali algebrici e/o di corrispondenti forme quadratiche binarie definite positive con determinante della forma -24n+1 e ulteriori vincoli sui coefficienti, applicando poi la corrispondenza di Gross-Kohnen-Zagier.

Tuttavia, da quel punto in poi, temo che sarebbe un lavoro improbo proseguire sia nella esemplificazione con poche righe di codice in un HLL idoneo, sia nella spiegazione elementare del risultato, mantenendola comprensibile per la maggioranza dei lettori.

Al di là degli aspetti computazionali (sono tuttora in corso ulteriori studi anche dal punto di vista algoritmico), l'importanza matematica di questa formula è veramente notevole, perché conferma delle ipotesi "visionarie" risalenti a Ramanujan e svela un inatteso comportamento macroscopico di tipo frattale nella funzione di partizione. Invito tutti gli interessati a seguire la lezione di Ken Ono - in english, of course - che richiede un background matematico minimale.

A causa delle considerazioni sopra riassunte, seppur a malincuore, ho dovuto pertanto optare per limitare drasticamente lo spazio dedicato al lavoro di Ono e soci, che invece meriterebbe ben altra attenzione per il suo portato filosofico ed epistemico...

Stewart M. (2003), "A superfast Toeplitz solver with improved numerical stability", SIAM Journal on Matrix Analysis and Applications, 25: 669-693.

Kravanja P., Van Barel M. (2006) "Coupled Vandermonde matrices and the superfast computation of Toeplitz determinants", NUMERICAL ALGORITHMS, Vol. 24, No. 1-2, 99-116.

Chandrasekeran S., Gu M., Sun X., Xia J., Zhu J. (2007), "A superfast algorithm for Toeplitz systems of linear equations", SIAM Journal on Matrix Analysis and Applications, 29: 1247-1266.

Hsuan-Chu Li (2011), "On Calculating the Determinants of Toeplitz Matrices", Journal of Applied Mathematics & Bioinformatics, vol. 1, no. 1, 55-64.

Commenti