ORDINARE IL MONDO

Ordinare il mondo

La matematica intesa come una razionalizzazione dell’esperienza, secondo la concezione del filosofo e matematico italiano Federigo Enriques (1871-1946), ha sempre cercato di descrivere i più svariati fenomeni della natura con calcoli ed equazioni, così come ha sempre elaborato, più in generale, una scienza certa e positiva del sistema di idee e concetti di cui sono intessuti il nostro sapere e il nostro comportamento. In un certo senso si tratta sempre di modelli, cioè di schemi atti a rappresentare, riprodurre o simulare varie specie di fenomeni in modo astratto e universale, con l’intento di includervi il maggior numero di cose possibili e diversificate tra loro. La rappresentazione consiste non solo in sistemi di equazioni o in altri strumenti (come grafi, matrici o sistemi dinamici) applicabili alla fisica, alla biologia o all’economia, ma anche in complesse costruzioni teoriche con cui si tenta di definire e limitare l’ambito di oggetti che rientrano in una certa categoria. I modelli (equazioni) della fisica matematica sono definiti di solito sul continuo, cioè con variabili che assumono valori nel campo dei numeri reali (razionali e irrazionali), ma esistono diverse teorie plausibili del continuo e del concetto di numero reale: J.W. Richard Dedekind (1831-1916) lo concepì come una sezione del corpo razionale, Georg F.L.Ph. Cantor (1845-1918) come una classe di equivalenza di insiemi di numeri razionali, Karl T.W. Weierstrass (1815-1897) come una classe di equivalenza di insiemi di razionali con somma limitata. Analogamente vi sono diverse definizioni possibili, equivalenti tra loro, del concetto di algoritmo: le funzioni ricorsive, il λ-calcolo, il formalismo di Andrej A. Markov (1903-1979) e la macchina di Alan M. Turing (1912-1954).

Occorre distinguere tra i modelli in diverse discipline. Nella matematica applicata il modello consiste di solito in un sistema di equazioni o in altri formalismi con cui si tenta di riprodurre, per astrazione, un fenomeno reale. Nella logica matematica i modelli consistono, invece, in sistemi di oggetti atti a rappresentare un sistema formale astratto: per es., rappresentazioni di uno stesso sistema formale dell’aritmetica sono i numeri naturali come cardinali di insiemi finiti, oppure i soli numeri pari (assumendo l’addizione +2 come rappresentativa dell’operazione di successore). In questa seconda accezione i modelli dovrebbero differire soltanto per una relazione di isomorfismo, risultando infine indistinguibili rispetto alle proprietà che contano e che definiscono astrattamente un sistema formale o una struttura algebrica. A questo proposito, il matematico tedesco Hermann Weyl (1885-1955) dichiarava che la scienza rimane del tutto indifferente all’‘essenza’ dei suoi oggetti e che sarebbe una pazzia aspettarsi che l’atto conoscitivo riveli all’intuizione una qualche essenza segreta delle cose; in questo senso l’isomorfismo, egli sosteneva, traccia i confini insormontabili della nostra conoscenza per mezzo di modelli. Un esempio è il fatto che un qualsiasi campo completo (rispetto alle operazioni di somma e prodotto) e con una struttura d’ordine (≤) è isomorfo al campo dei numeri reali, motivo che autorizza a concludere, in base all’importanza centrale della proprietà di completezza nel calcolo differenziale e integrale, che il campo reale si presenta, in un certo senso, ‘una sola volta’, anche se non richiede necessariamente l’apparato denotazionale più familiare e intuitivo. Il calcolo differenziale e integrale fornisce a sua volta lo strumento per una prima elaborazione dei modelli della natura. Tuttavia, il calcolo algebrico e numerico da cui è nato il calcolo differenziale ha rivelato, nel 20° sec., uno statuto autonomo e perfino prioritario: nel calcolo scientifico i modelli differenziali definiti nel continuo dei numeri reali si approssimano, nel discreto, con modelli algebrici e calcoli aritmetici, e la natura di questa approssimazione, la trasmissione di contenuti e di strutture matematiche dal continuo al discreto deve essere ancora analizzata in tutti i suoi svariati aspetti. Nella matematica computazionale e applicata la funzionalità di un modello si misura soprattutto in base alla capacità predittiva sull’andamento di un fenomeno. Per questo non ha importanza soltanto la maggiore o minore adeguatezza di un’equazione a descrivere un fenomeno, ma anche la possibilità di una sua effettiva risoluzione numerica.

Modelli e algoritmi nel calcolo scientifico

Il calcolo scientifico si fonda su due pilastri principali: le equazioni della fisica matematica, su cui si basa la consueta accezione di modello, e gli algoritmi. Questi si devono quindi considerare come un prolungamento necessario dell’idea classica di modello fondata sulla costruzione di equazioni differenziali o integrali. Informalmente un algoritmo consiste in un processo caratterizzato da almeno tre requisiti principali: la precisione prescrittiva che non lasci spazio all’arbitrio, ovvero un requisito di definitezza; la possibilità di avviare il calcolo da un insieme di dati iniziali, variabili entro certi limiti, ovvero un requisito di generalità; la necessità di ottenere un risultato, calcolabile in base ai dati iniziali, ossia di orientare il processo affinché sia conclusivo.

Il terzo punto può anche essere tradotto in un’esigenza di effettività. Il matematico statunitense Alonzo Church (1903-1995) usò, infatti, l’espressione effettivamente calcolabile per denotare una proprietà che deve rientrare in ogni possibile caratterizzazione formale. In diverse definizioni del concetto di algoritmo si evidenzia il suo carattere ripetitivo: ogni algoritmo è un processo essenzialmente iterativo, in cui uno stesso operatore è applicato molte volte. Nel calcolatore si devono considerare ‘materializzabili’ gli algoritmi e, in generale, le formule della matematica costruttiva. Tutta la logica e tutta la matematica di Luitzen E.J. Brouwer (1881-1966, il principale esponente della scuola intuizionista del primo Novecento) – osservava John von Neumann in una lettera a Norbert Wiener del 1946 – possono essere date da un appropriato meccanismo, possono cioè tradursi in reti elettroniche, in sistemi nervosi idealizzati che danno contenuto reale alle formule con cui la matematica costruisce effettivamente le soluzioni di un’equazione. Un presupposto di questa meccanizzazione delle formule doveva essere un procedimento di approssimazione: i modelli, precisavano sia Herman H. Goldstine (1913-2004) sia J. von Neumann (1913-1957), consistono di solito in problemi analitici in cui le variabili assumono valori nel continuo, cioè nel campo dei numeri reali, e che hanno un carattere totalmente o parzialmente implicito. Ma nel calcolo digitale questi si devono approssimare con procedure puramente aritmetiche, finite (finitistic) ed esplicite, di carattere essenzialmente iterativo. Una precisazione che ha trovato conferma soprattutto nei metodi iterativi per la risoluzione numerica di sistemi di equazioni lineari e non lineari, nella teoria dei sistemi dinamici e nel calcolo del minimo di una funzione. Prototipo di questi metodi è il criterio generale del punto fisso: si esprime un problema nella forma x=g(x), in cui il punto x rimane invariato (fisso) sotto l’azione dell’operatore g, e si calcola passo per passo x_k=g(x_k−1), per k=1, 2, 3,…, iniziando da un valore assegnato x0, in modo che x_k converga alla soluzione esatta. Il grande vantaggio di uno schema iterativo consiste nella semplicità e uniformità delle operazioni richieste a ogni passo che consentono una completa meccanizzazione del processo di risoluzione numerica di problemi di elevate dimensioni. Questa meccanizzazione, soprattutto nella seconda metà del Novecento, ha contribuito a un processo di trasformazione dell’arte del calcolo, imperniata su artifici e singoli espedienti, in una scienza computazionale, basata sullo studio sistematico della stabilità e della complessità di classi di algoritmi, e sulla definizione di precisi modelli di calcolo sequenziale e parallelo.

La scienza del calcolo non si è sempre sviluppata in modo unitario. Se negli anni Sessanta e Settanta l’informatica mirava a una fondazione teorica del concetto di algoritmo, il calcolo scientifico era più orientato a fini pratici e applicativi. Ora si è maggiormente propensi a includere nel concetto teorico di effettività, elaborato in ambito logico e informatico, l’efficienza tipica delle procedure della matematica applicata. Un algoritmo numerico risulta efficiente se la propagazione dell’errore nei calcoli non compromette il senso del risultato e se la complessità computazionale risulta accettabile. Semplici calcoli ricorsivi hanno spesso una complessità troppo elevata. Per risolvere un sistema di 50 equazioni in 50 incognite con il metodo di Cramer un moderno calcolatore digitale impiegherebbe un tempo – ancorché finito – paragonabile all’età dell’Universo. Ma la matematica applicata richiede la risoluzione di decine di migliaia di equazioni: il calcolo scientifico è su grande scala, e per l’enorme dimensione dei calcoli serve un’attenta analisi dell’efficienza delle procedure. Specialmente nella seconda metà del 20° sec. le richieste pratiche di effettività ed efficienza hanno acquisito un significato teorico di presupposti necessari alla definizione di un algoritmo.

Aritmetizzazione

Gli scopi applicativi s’intrecciano regolarmente con istanze teoriche e fondazionali. L’aritmetizzazione messa in atto dal calcolo scientifico insegue uno scopo analogo, se pur di segno contrario, al progetto di fine Novecento di fondare l’edificio matematico su proprietà di insiemi di interi. Inoltre, lo studio dei limiti di complessità computazionale si serve spesso delle tecniche per lo studio di classici problemi indecidibili, come la risolubilità delle equazioni diofantee. Il problema P-NP e la questione se esistono procedure efficienti per venire a capo di problemi di natura combinatoria risolubili in tempo finito ma di enorme complessità, dimostrano una sostanziale affinità con ricerche di natura fondazionale, in cui s’individuano problemi che non hanno soluzione algoritmica neppure in linea di principio. Il progetto di rappresentare l’infinito con un sistema coerente di simboli si è in parte spostato sulla questione di come si possa venire a capo di una quantità finita, ma ingente, di operazioni elementari finite.

La matematica computazionale della seconda metà del Novecento ha reso più esplicito un fatto di cui erano già consapevoli i matematici francesi e gli esponenti della scuola intuizionista di inizio secolo: la tesi di fondo di David Hilbert (1862-1943), che ogni concetto, teorema e dimostrazione devono potersi ricondurre al finito, era insufficiente, o almeno si è chiarito come il finito debba consistere essenzialmente in un processo che, sulla base di dati iniziali, abbia termine in un numero limitato di passi. La riduzione al finito deve cioè prendere corpo in un processo algoritmico, che possa poi diventare l’effettivo lavoro di un computer. La medesima insistenza, da parte di Goldstine e di von Neumann, sulla necessità di approssimare i modelli nel finito mediante calcoli aritmetici di tipo iterativo (nel termine finitistic è implicito il richiamo a Hilbert) ha spostato il programma hilbertiano – la fondazione dell’infinito matematico sul finito – nell’ambito di una scienza pratica degli algoritmi.

La riduzione a procedure aritmetiche rientra pure nel classico obiettivo di ricondursi a modelli semplificati di calcolo. Lo sviluppo di funzioni in serie di potenze e la risoluzione di un’equazione algebrica per radicali fanno parte di questo disegno. Il calcolo su grande scala si basa, inoltre, sull’esecuzione ripetuta di operazioni aritmetiche e sulla risoluzione di sottoproblemi che rimangono completamente nascosti, nel senso che il soggetto umano non è in grado di vedere né i dati in ingresso né quelli in uscita, spesso generati da un programma e immediatamente utilizzati da un altro. Si possono usare in questo modo, per es., operazioni elementari di aritmetica reale o complessa, trasformate veloci di Fourier per l’elaborazione di immagini, prodotti matrice×vettore in procedure per il calcolo degli autovalori di una matrice o per il calcolo del minimo di una funzione. È evidente che il calcolo relativamente nascosto di queste funzioni elementari, in cui si articolano le più complesse strategie computazionali, richiede un grado di affidabilità maggiore delle operazioni di cui si possono controllare in modo più diretto i risultati. Soltanto in tempi recenti si è avviata una ricerca sistematica sul grado di efficienza di questo calcolo ‘interno’ e relativamente inaccessibile: uno studio necessario per giustificare l’atto di delega, da parte del soggetto umano, a un processo completamente automatico.

Strutture e complessità

Il concetto di struttura, in matematica, ha per lo più una forte connotazione astratta e teorica, di derivazione bourbakista. Tuttavia, uno dei principali motivi che inducono a concepire nuove strutture deriva da esigenze computazionali. Gli aspetti computazionali sono già impliciti nel concetto generale di campo di numeri. Informalmente, un insieme F è un campo di numeri se, operando sugli elementi di F con le ordinarie operazioni aritmetiche, si ottengono ancora elementi di F. Il campo dei numeri razionali gode di questa proprietà di chiusura: operando su numeri razionali con somme, moltiplicazioni e divisioni si ottengono sempre e solo numeri razionali; inoltre, ogni campo di numeri è un’estensione del campo dei razionali. Un’analoga proprietà di chiusura (rispetto a operazioni di sostituzione e diagonalizzazione) è necessaria – come dimostrò Dedekind nel suo celebre teorema di ricorsione (Was sind und was sollen die Zahlen?, 1888) – per definire in modo rigoroso le operazioni fondamentali dell’aritmetica e il concetto teorico di calcolabilità. Le strutture di gruppo, anello o campo, sono servite a fondare l’algebra moderna fin dalla seconda metà del 19° sec., ma anche, in tempi più recenti, a dare un assetto rigoroso allo studio della complessità computazionale di importanti classi di funzioni. Lo studio della complessità algebrica e numerica si è sviluppato intorno alla metà del 20° sec., con i primi progressi dell’informatica e del calcolo scientifico. Lo scopo era misurare la difficoltà di calcolo di elementi di un anello o di un campo più o meno estesi, per es. di un generico polinomio oppure di un insieme di funzioni razionali di grado assegnato, le funzioni del calcolo ‘interno’ che vanno a formare la trama molecolare del calcolo su grande scala. La struttura specifica del campo in cui variano i coefficienti di una funzione razionale f può intervenire, in modo essenziale, nella valutazione del numero di operazioni aritmetiche necessarie o sufficienti al calcolo di f.

Non sempre le classiche procedure dell’algebra e della matematica numerica – si pensi al metodo di Horner-Ruffini per il calcolo di un polinomio – raggiungono la complessità minima. Un celebre risultato ottenuto dal matematico tedesco Volker Strassen (n. 1936) nel 1969 ha stabilito che il metodo di Gauss per risolvere un sistema di n equazioni lineari in n incognite, che usa un numero di operazioni dell’ordine di n3, non è ottimo, perché esistono algoritmi di complessità O(n^α) con α<3. Strassen ha dimostrato che si possono moltiplicare due matrici quadrate 2×2 con sole 7 moltiplicazioni (invece di 8) su un anello non commutativo e che si possono, quindi, moltiplicare due matrici n×n, dimezzando in modo ricorsivo la loro dimensione, in O(n^α) operazioni con α=log27≈2,81. Questo implica a sua volta che si può invertire una matrice non singolare, o risolvere un sistema di n equazioni lineari, in O(n^α) operazioni aritmetiche. Dal risultato di Strassen si deduce che la complessità asintotica (al crescere della dimensione n) della risoluzione di un sistema di n equazioni lineari dipende dalla complessità del prodotto di matrici di piccola dimensione, definito da semplici funzioni razionali bilineari. Tuttavia, la complessità di queste funzioni, ancorché semplici, è in genere di difficile valutazione, a causa di un’esplosione combinatoria simile a quella di problemi di ottimizzazione classificati come NP-completi, che costituiscono ora una delle sfide più importanti per la scienza del calcolo. Valutare la minima complessità asintotica (al crescere di n) dell’inversione di una matrice per mezzo di un algoritmo diretto sequenziale, è una delle principali questioni aperte della matematica computazionale e applicata, dal momento che quasi tutti i problemi, per via di successive approssimazioni, si riconducono alla risoluzione numerica di sistemi di equazioni lineari. La struttura di determinate matrici risulta, a questo proposito, di centrale importanza. Si dimostra, innanzitutto, che la complessità moltiplicativa di un insieme di forme bilineari, a^TM_kb – dove M_k (k=1, 2,…, p) sono matrici a elementi in un campo F e a, b sono vettori (verticali) di indeterminate su F (T indica la trasposizione) –, è data precisamente, nel caso generale (non commutativo), dal minimo numero q per cui ciascuna M_k si può riscrivere come prodotto di tre matrici, cioè M_k=UD_kV, ove D_k è diagonale q×q e U, V sono matrici che non dipendono da k. Il calcolo di questo numero q, che si chiama rango tensoriale delle matrici M_k, è in genere difficile e, in semplici casi come, per es., il prodotto di matrici 3×3, è un problema NP-completo. Per alcuni spazi di matrici con struttura speciale, tuttavia, come le matrici di Hankel e di Toeplitz (v. oltre), esso è perfettamente noto per determinati campi F. Una conseguenza è la possibile semplificazione di numerose operazioni fondamentali dell’algebra dei numeri e dei polinomi.

Matrici con struttura e trasformate veloci

Il termine struttura, in ambito computazionale, è usato ormai abitualmente soprattutto nella letteratura dell’algebra lineare numerica, anche se manca una sua precisa definizione formale. Non è questa una restrizione troppo grave: nell’ingegneria, nell’informatica e nella fisica l’algebra lineare numerica occupa ora uno spazio rilevante, dovuto in parte, come ha notato il matematico statunitense Gilbert Strang (n. 1934), allo spostamento di attenzione dall’analogico al digitale. Alle funzioni si sostituiscono vettori, e si combina lo studio degli spazi n-dimensionali con l’applicazione di matrici. L’algebra delle matrici diventa, quindi, uno dei principali strumenti di analisi dei modelli della natura e della loro traduzione in pura informazione numerica.

La struttura può riferirsi innanzitutto all’eventualità che tra gli n2 elementi di una matrice A n×n esista una relazione funzionale che permetta di individuare A con un numero di parametri k minore di n2. Inizialmente la struttura aveva il carattere di una semplice rappresentazione, un criterio per memorizzare una matrice nel modo più semplice ed economico possibile. George E. Forsythe (1917-1972) usava per questo, in un articolo del 1967 (Today’s computation­al methods of linear algebra, «SIAM review», 1967, 9, 3, pp. 489-515), il concetto di contenuto informativo, ripreso successivamente da altri autori, con cui intendeva il numero di celle necessario per memorizzare, nel calcolatore, i dati e il programma sufficienti a generare tutti gli elementi di una matrice. In seguito la struttura ha assunto un significato più spiccatamente computazionale, legato al numero di operazioni aritmetiche richieste nelle ordinarie operazioni matriciali: matrice×vettore, matrice×matrice e inversione.

Un esempio di matrice n×n con struttura sono le matrici circolanti C, in cui ogni riga i si ottiene dalla riga precedente i−1 con una permutazione ciclica dei suoi elementi: se la prima riga è [c1 c2… c_n], la seconda riga è [c_n c1… c_n−1], la terza è [c_n−1 c_n… c_n−2] e l’ultima è [c2 c3… c_n c1]. Quindi n parametri sono sufficienti a definire C. Queste matrici godono di proprietà speciali, che semplificano tutte le operazioni: la somma e il prodotto di matrici circolanti sono circolanti e l’inversa di una matrice circolante è circolante. Le matrici circolanti formano un’algebra commutativa i cui elementi sono riconducibili a una forma diagonale mediante una trasformata veloce nota come fast Fourier transform (FFT). Si ha cioè C=FDF⁻1, dove D rappresenta la matrice diagonale n×n degli autovalori, e F ha elementi F_rs=ω(r−1)(s−1) (r, s=1, 2,…, n), con ω=eⁱ2^π/n e i=√−−−1. Il contenuto informativo di C è, quindi, racchiuso negli elementi (autovalori) non nulli di D, mentre la struttura dell’algebra è definita dalla matrice F di Fourier. Il prodotto di C per un vettore b è dato da FDF⁻1b, e si riconduce, quindi, al calcolo di due trasformate veloci (una diretta e una inversa) più il prodotto di una matrice diagonale per un vettore. Studi recenti hanno fatto scoprire diverse trasformate analoghe alla FFT che definiscono altrettante algebre di matrici commutative simultaneamente diagonalizzabili. La risoluzione effettiva di molti problemi numerici dipende dall’uso di queste algebre, che permette di ricondurre il calcolo a poche trasformate veloci. L’individuazione di una matrice A n×n con un numero di parametri minore di n2 può anche dipendere dall’eventualità che A abbia diversi elementi nulli. Tra le prime matrici di cui si è analizzata la struttura figurano infatti le matrici tridiagonali, i cui elementi non nulli sono disposti sulla diagonale principale e sulle due diagonali adiacenti (cioè a_ij≠0 se e solo se |i−j|≤1). In questo caso A dipende da soli 3n−2 parametri (2n−1 se A è simmetrica). Matrici tridiagonali – o che differiscono poco dalla forma tridiagonale –, a elementi o a blocchi, intervengono, per es., nella discretizzazione, con il metodo delle differenze finite di problemi al contorno per equazioni alle derivate parziali di tipo ellittico, come pure nell’approssimazione con funzioni spline. Le funzioni spline in una dimensione sono funzioni coincidenti con un polinomio in ciascuno dei sottointervalli che definiscono una partizione di un dato intervallo sulla retta reale e che possiedono un certo grado di regolarità nei punti di congiunzione dei sottointervalli. Queste funzioni erano usate inizialmente in problemi di architettura navale per disegnare la linea di uno scafo unendo tra loro asticelle o stecche di legno (da cui il termine spline). La loro utilità si è poi estesa ad altri ambiti, come l’industria automobilistica, ove sia richiesto il disegno automatico di curve di profilo.

Se è ovvio che 3n−2 parametri individuano una matrice tridiagonale T n×n, meno evidente è che lo stesso numero di parametri basta a definire la sua inversa T⁻1, che è una matrice piena (priva di zeri). Si ha precisamente (T⁻1)_ij=a_ib_j per j≥i e (T⁻1)_ij=c_id_j per j≤i, con i vincoli a_ib_i=c_id_i e a1=c1=1 che riducono i 4n parametri a, b, c, d a soli 3n−2. Questa invarianza di contenuto informativo rispetto a operazioni d’inversione vale più in generale per le matrici A con struttura a banda, con gli elementi non nulli di A collocati sulle p diagonali sopra e sotto la diagonale principale, per un certo intero positivo p<n; ovvero A_ij=0 per |i−j|>p (per p=1, A è tridiagonale). L’inversa di una matrice a banda A si può descrivere con lo stesso numero di parametri sufficienti a definire A, anche se è una matrice piena. Matrici che hanno questa struttura si chiamano semiseparabili. Oltre al fatto che la loro inversa ha una struttura a banda, le matrici A semiseparabili sono caratterizzate da una proprietà di particolare interesse computazionale: ogni sottomatrice situata al di sopra della p-esima diagonale superiore di A ha rango ≤p, e lo stesso accade per ogni sottomatrice situata al di sotto della p-esima diagonale inferiore. Per p=1 (l’inversa ha forma tridiagonale) i blocchi B superiori e inferiori della matrice semiseparabile hanno rango ≤1, ovvero B ha la forma uv^T, dove u e v sono due vettori di dimensione, rispettivamente, r e s se B è una matrice r×s. Questo implica che per calcolare un prodotto matrice×vettore Bx, riscritto nella forma u(v^Tx), bastano r+s moltiplicazioni, invece di rs. Partizionando ripetutamente (in modo ricorsivo) una matrice A semiseparabile la cui inversa abbia forma tridiagonale, si trova che il prodotto matrice×vettore Ay si calcola in 2nlog2n, invece di n2, moltiplicazioni (se n=1000, l’ordine è delle decine di migliaia invece del milione di operazioni). La struttura di una matrice A è spesso nascosta, non immediatamente riconoscibile nella disposizione o nel valore degli elementi A_ij, come per le inverse di matrici a banda e di Toeplitz. Le matrici di Toeplitz intervengono nell’aritmetica polinomiale, nella teoria delle serie temporali, nelle operazioni di filtraggio di segnali e nella discretizzazione di operatori integrali. Una matrice di Toeplitz n×n ha elementi uguali su ciascuna diagonale (sulla diagonale principale e sulle diagonali parallele alla principale), ovvero A_ij=A_i−j, ed è quindi definita da 2n−1 parametri indipendenti. Le matrici circolanti sono matrici di Toeplitz a cui si aggiunge una condizione di periodicità. Questo si traduce nel fatto che il prodotto di una matrice A di Toeplitz per un vettore b si risolve in trasformate veloci: la matrice A, di ordine n, si ‘immerge’ in una matrice C di ordine 2n e il vettore b si completa con n zeri, in modo da ritrovare Ab nelle prime n componenti del prodotto di C per il nuovo vettore. L’inversa di una matrice di Toeplitz non singolare non è di Toeplitz, ma è descrivibile con un numero lineare di parametri. Da questo si deduce che la struttura consiste non tanto nell’uguaglianza degli elementi sulle diagonali, quanto in una proprietà più riposta e non immediatamente visibile. Il concetto chiave è dato dalla nozione di rango di spostamento. Un operatore di spostamento è una matrice Z che ha elementi uguali a 1 sulla diagonale sottostante la diagonale principale e 0 altrove. La matrice ZAZ^T si ottiene allora da A spostando di una posizione verso il basso a destra i suoi elementi e collocando alcuni zeri sulla prima riga e sulla prima colonna. Il rango di spostamento viene definito dal rango r della differenza A−ZAZ^T, tra la matrice A e la matrice ‘spostata’ ZAZ^T. La proprietà di cui godono sia una matrice di Toeplitz simmetrica A sia la sua inversa A⁻1 è di avere un rango r di spostamento piccolo, precisamente r=2. Questo fatto permette di scrivere sia A sia A⁻1 come somma di soli 2 prodotti di elementi di speciali algebre di matrici riducibili a forma diagonale mediante trasformate veloci. Il numero di operazioni aritmetiche per operazioni matrice×vettore si riduce allora drasticamente, tipicamente da O(n2) a O(nlogn).

La struttura di una matrice è spesso legata alla possibilità di ottenere vantaggi computazionali per mezzo di trasformate veloci. La piena consapevolezza della potenza di semplificazione della FFT risale al celebre contributo An algorithm for the machine calculation of complex Fourier series («Mathematics of computation», 1965, 19, 90, pp. 297-301) di James W. Cooley (n. 1926) e John W. Tukey (1915-2000) in cui si è dimostrato che per la discrete Fourier transform (DFT) su n punti è sufficiente un numero di operazioni proporzionale a nlogn anziché a n2. L’algoritmo Cooley-Tukey ha comunque diverse varianti, tra cui una, dovuta a Shmuel Winograd (n. 1936), dipende dal fatto che la parte ‘essenziale’ di una matrice di Fourier si può riscrivere come una matrice circolante, circostanza che permette di risolvere la DFT su n punti, per n altamente composito (n=prodotto di potenze di numeri primi piccoli), con una serie di convoluzioni (prodotti polinomiali) di piccola dimensione. Un’osservazione di Beresford N. Parlett (n. 1932) suggerisce di riformulare l’algoritmo di Winograd come una decomposizione della parte essenziale della matrice di Fourier in un prodotto UDV, ove U e V sono matrici ‘semplici’, che non contribuiscono alla complessità di calcolo, e D è una matrice diagonale a blocchi 2×2 che contiene l’informazione essenziale sul numero di moltiplicazioni della DFT. Per es., per n=5, la struttura di D consente di ridurre il numero di moltiplicazioni da 25 a 4.

La DFT si è rivelata uno strumento universale, utile in ogni circostanza e in ogni applicazione. L’aritmetica degli interi e dei polinomi si semplifica grazie alla FFT, con una riduzione di complessità delle operazioni più importanti da O(n2) a O(nlogn). La più tipica applicazione della FFT è tuttavia il filtraggio di segnali. Per es., se f(x) è una funzione periodica rappresentabile come serie trigonometrica, cioè come somma di infiniti termini della forma α_kcoskx+β_ksenkx, e si ha un segnale della forma g(x)=f(x)+r(x), ove r(x) rappresenta un rumore ad alte frequenze, si può recuperare buona parte dell’informazione contenuta in f. A questo scopo con una DFT si calcola il polinomio trigonometrico di g(x) su una serie di punti e si arresta la somma a un numero m di termini in maniera tale da eliminare le frequenze più alte. Ammesso che le ampiezze massime delle componenti di f, ovverosia (α_k2+β_k2)1^/2, siano trascurabili per k sufficientemente elevato, ci si aspetta che la somma così troncata assomigli alla funzione f(x), cioè al segnale da recuperare. Diverse affinità si riscontrano tra la DFT e la trasformata di Radon: le moderne tecniche diagnostiche della tomografia computerizzata si basano su un’importante teoria del matematico austriaco Johann Radon (1887-1956) del 1917. La teoria delle ondine (wavelets) offre ora uno strumento particolarmente potente per le operazioni di filtraggio, consentendo la separazione di un segnale in bande di frequenza in tempi ottimali. Le operazioni aritmetiche fondamentali coinvolgono regolarmente matrici a banda di Toeplitz. Matematicamente, un filtro digitale consiste, infatti, al variare di k, in una somma di prodotti h(k)x(n−k) rappresentabile nella forma di un prodotto matrice×vettore, ove la matrice è di Toeplitz triangolare inferiore, con elementi uguali sulla diagonale principale e su alcune diagonali sottostanti.

Le algebre di matrici simultaneamente diagonizzabili (come le circolanti) sono pure applicabili nelle tecniche di precondizionamento. L’efficienza di diversi algoritmi per la risoluzione numerica di sistemi lineari Ax=b, con A matrice di Toeplitz simmetrica, dipende dalla distribuzione degli autovalori di A, e si può tentare di modificare questa distribuzione per ottenere vantaggi computazionali. In uno dei metodi più efficienti per la risoluzione del sistema Ax=b, il metodo del gradiente coniugato, si calcola una matrice P di precondizionamento in modo che nel sistema equivalente P⁻1Ax=P⁻1b la distribuzione degli autovalori di P⁻1A comporti un’accelerazione della convergenza e, quindi, una riduzione della complessità. Si dimostra che P può essere scelta all’interno di uno spazio di matrici strutturate, per es. come la matrice circolante C più vicina ad A (con un’opportuna definizione della distanza tra due matrici, l’effettiva esistenza di C è una conseguenza del teorema della proiezione di Hilbert).

Vengono qui di seguito riassunti tre distinti filoni di ricerca, in cui l’algebra matriciale ha un ruolo decisivo nell’elaborazione di metodi di risoluzione numerica.

Problemi differenziali

Fu Joseph-Louis Lagrange (1736-1813) a introdurre il concetto di energia potenziale, dimostrando la sua utilità per la dinamica e per la teoria della gravitazione. Divenne poi chiaro che idee concernenti diversi ambiti della fisica, come temperatura, pressione e voltaggio, hanno un ruolo paragonabile, rispettivamente, nella teoria del calore, nella fluidodinamica e nell’elettricità. Questo fatto condusse a una teoria generale del potenziale con Pierre-Simon Laplace, Siméon-Denis Poisson, Carl F. Gauss e George Green nel 19° secolo. In questa teoria si colloca anche il celebre problema di Dirichlet, prototipo di una serie di problemi al contorno per equazioni alle derivate parziali di tipo ellittico. Nello spazio fisico a tre dimensioni il problema di Dirichlet riguarda la distribuzione della temperatura, in condizioni di equilibrio, in un solido omogeneo che occupi una regione Ω la cui superficie Γ sia mantenuta a una data temperatura g(s), ove s è un punto variabile su Γ. Si tratta cioè di calcolare, dopo un tempo sufficiente affinché il flusso di calore si stabilizzi, il valore della temperatura in Ω. Se si denota con u(r) questo valore (r è un vettore di componenti x, y, z), u deve soddisfare l’equazione di Laplace Δ2u=0 in ogni punto r di Ω, dove Δ2u indica la somma delle derivate seconde di u rispetto a x, y, z. Il problema di Dirichlet è allora il seguente: data g(s) trovare u che soddisfi l’equazione di Laplace al variare delle coordinate x, y, z e tale che, per ogni punto s sulla superficie Γ, u(r) tenda a g(s) per r che tende a s. Una simile funzione u(r) si dice armonica in Ω.

Uno dei fondamentali mutamenti di prospettiva del 20° sec. è stato l’allargamento del concetto di soluzione, che ha implicato a sua volta l’allargamento del concetto di modello a tutto il procedimento di algebrizzazione dell’equazione differenziale, fino alla sua approssimazione con un sistema di equazioni lineari e allo svolgimento di procedure puramente aritmetiche. La soluzione analitica di un problema differenziale è una funzione che soddisfa l’equazione in ogni punto di una regione Ω e in ogni punto del suo contorno Γ. La soluzione numerica è, invece, la soluzione del problema aritmetico con cui si approssima l’equazione, e non può quindi considerarsi come una semplice stima diretta della soluzione analitica in un punto di Ω o di Γ. Tra le principali strategie per ottenere la soluzione analitica figurano la separazione delle variabili, la funzione di Green e la formulazione variazionale. Quest’ultima, nel caso del problema di Dirichlet, permette di interpretare la soluzione u come la funzione che minimizza un integrale di energia, secondo il celebre principio di Dirichlet. Il metodo di Green è, invece, il modo più generale per invertire un’equazione alle derivate parziali Lu=f, dove L è un operatore differenziale (L=Δ2 nel caso dell’equazione di Laplace). Si esprime, infatti, la soluzione u nella forma L⁻1f, ottenendo una rappresentazione di L⁻1 come un operatore integrale il cui nucleo è la funzione di Green relativa al problema Lu=f: una strategia che ricorda la possibile risoluzione di un sistema di equazioni lineari Ax=f nella forma A⁻1f.

Diverso è il criterio per calcolare una soluzione numerica. Il matematico e fisico tedesco Carl D.T. Runge (1856-1927), uno dei principali promotori dello sviluppo dell’analisi numerica e del calcolo scientifico nel Novecento, propose, nel 1908, di risolvere il problema di Laplace in un dominio regolare Ω del piano (x, y), con condizioni al contorno, con un metodo alle differenze finite. Runge introdusse una griglia Ωh di passo h, di punti x_i, y_j, e sostituì alle derivate di u dei quozienti alle differenze. Runge pensò di trascurare l’errore di discretizzazione, ottenendo un sistema di equazioni lineari al posto del problema differenziale. Ricerche successive perfezionarono questo risultato, utile anche per dimostrare un teorema costruttivo di esistenza di una soluzione del problema di Dirichlet: un passo importante per rafforzare la tesi che la matematica può essere aritmetizzata almeno in linea di principio. Con il metodo alle differenze finite, assumendo n=2 e ricoprendo con una griglia discreta di punti x_i, y_j una regione rettangolare Ω sul piano x, y assieme al contorno Γ, facendo variare i e j tra 1 e n, il problema Δ2u=f(x, y), con la condizione u=0 su Γ, si traduce in un sistema di equazioni lineari Az=b, in cui la matrice A, di n2 righe e n2 colonne, ha una struttura tridiagonale in termini di blocchi (matrici n×n) che sono a loro volta diagonali o tridiagonali. La matrice A è simmetrica, irriducibile, definita positiva (ha tutti gli autovalori positivi) e possiede una predominanza diagonale (in valore assoluto ogni termine sulla diagonale è maggiore o uguale, e in qualche caso maggiore, della somma degli altri termini sulla stessa riga): proprietà che consentono di risolvere il sistema con i classici metodi iterativi di Jacobi o di Gauss-Seidel, oppure con un metodo diretto basato sulla decomposizione di A nel prodotto LL^T, dove L è triangolare inferiore (decomposizione di Cholesky). Un’approssimativa ma sufficiente valutazione dello spettro di A, da cui si deduce l’applicabilità di tali metodi, si deve in questo, come in altri casi, ai celebri teoremi di Gershgorin sulla collocazione degli autovalori di una matrice sul piano complesso.

Nel metodo alle differenze finite la soluzione numerica dipende dalla struttura della matrice del sistema di equazioni algebriche che approssima nel discreto il problema differenziale. Questa struttura determina l’efficienza degli algoritmi numerici: la loro convergenza, la complessità aritmetica e la limitazione dell’errore nei calcoli. Una circostanza analoga si verifica con altri metodi numerici, basati ancora sull’indebolimento del concetto di soluzione, come le tecniche variazionali applicate a problemi al contorno per equazioni differenziali lineari di tipo ellittico. Si tratta di problemi della forma Lu=f in Ω, con condizioni al contorno di Dirichlet u=g (oppure di Neumann ∂u/∂n=h) su Γ, ove Lu denota un’espressione con derivate prime e seconde di u e con espressioni lineari in u (il termine ellittico si riferisce a condizioni sui coefficienti delle derivate seconde). Questi modelli descrivono diversi processi, come gli spostamenti di una membrana elastica vincolata al bordo Γ dovuti a una forza ortogonale, l’evoluzione nel tempo della temperatura di un corpo conduttore di calore e fenomeni di convezione-diffusione. In questi ultimi, con due sole variabili x e y, si ha un fluido incomprimibile che si muove sul piano in modo stazionario (la velocità dipende dallo spazio e non varia nel tempo). Se una sostanza inquinante viene sciolta nel liquido, la sua concentrazione in ogni punto dipende da due fenomeni distinti: la diffusione, causata dai moti browniani delle molecole, e la convezione, dovuta al fatto che le medesime particelle sono trasportate per inerzia dalla corrente del fluido. Una sorta di bilancio locale tra i due effetti separati di diffusione e convezione si rappresenta come un problema Lu=0 non omogeneo (con derivate prime e seconde di u) ove si assegni sul bordo la concentrazione di inquinante.

Per tali problemi può essere restrittivo cercare una soluzione analitica che sia regolare in ogni punto di Ω e di Γ e che verifichi le condizioni richieste dal modello. In molti casi non si richiede neppure questa regolarità; l’esistenza di tale soluzione non è nemmeno, in generale, facilmente dimostrabile. Conviene allora indebolire il concetto di soluzione mediante una tecnica variazionale, una strategia che risale a Lord Rayleigh (John W. Strutt, 1842-1919), Walter Ritz (1878-1909) e Boris G. Galerkin (1871-1945), consistente nel trovare un’espressione integrale J[φ], definibile per tutte le funzioni φ di una certa classe e che assuma un minimo esattamente per la funzione y che è soluzione del problema al contorno, e nel tradurre il problema di minimo in un sistema di equazioni lineari. In particolare, il metodo di Ritz (1909) consiste nel trovare una soluzione approssimata del problema di minimo sostituendo al posto di φ una combinazione lineare di funzioni coordinate φ1, φ2,…, φ_n, in modo che J[φ] sia interpretabile come una funzione dei coefficienti della combinazione. Si eguagliano quindi a zero le derivate di J[φ] rispetto ai coefficienti, che diventano le incognite di un sistema di equazioni lineari.

Il problema al contorno Lu=f in Ω, u=g su Γ può ricondursi, più in generale, a un problema variazionale astratto della forma seguente: trovare u in uno spazio di Hilbert V in modo che sia a(u, v)=F(v) per ogni v in V, dove a indica una forma bilineare che soddisfa determinate condizioni. Il metodo di Galerkin (1915) consiste nel ricondurre questo problema a una famiglia di problemi di dimensione finita per ottenere, infine, risolvendo un sistema di equazioni lineari, un’approssimazione numerica della soluzione u. A questo scopo si sostituisce a V una famiglia di sottospazi {V_h} di V, di dimensione finita, dipendenti da un parametro di discretizzazione h positivo che viene fatto tendere a 0, e si formula il problema variazionale in termini degli spazi V_h, cioè nella forma a(u_h, v_h)=F(v_h), dove u_h e v_h appartengono a V_h. Fissando un numero finito di funzioni coordinate φ_i nella cui combinazione lineare si possa esprimere ogni elemento di V_h, e riscrivendo le nuove relazioni variazionali con queste combinazioni al posto di u_h e con φ_i al posto di v_h, si ottiene infine un sistema di equazioni lineari Au_h=F_h, dove la matrice A ha elementi a(φ_j, φ_i) e F_h è il vettore che ha per elementi i valori assunti da F nelle φ_j. La costruzione effettiva degli spazi V_h e della base di funzioni φ_i si realizza con uno dei metodi numerici più diffusi per la risoluzione di problemi al contorno per equazioni alle derivate parziali, il cosiddetto metodo degli elementi finiti.

Una questione importante è come risolvere il sistema Au_h=F_h. Se la matrice A è simmetrica e definita positiva si può usare la decomposizione di Cholesky, un metodo diretto di cui si sa controllare la stabilità e la complessità. Ma la simmetria di A dipende dalla simmetria della forma bilineare a del problema variazionale continuo, non sempre verificata. Se A non è simmetrica si può usare il metodo iterativo di Richardson, prototipo di altri metodi più complessi e avanzati, come il metodo del gradiente coniugato e il generalized minimal residual method (GMRES). Il metodo di Richardson, applicato a un sistema Ax=b, calcola ripetutamente un’approssimazione numerica x_k della soluzione secondo lo schema x_k=x_k−1+ω(b−Ax_k−1), in base a un valore iniziale x0, per k=1, 2, 3,…, fino alla precisione richiesta. Al parametro ω viene assegnato un valore da cui dipende la velocità di convergenza della successione x_k alla soluzione, e quindi la complessità in termini di numero di operazioni aritmetiche, dominata essenzialmente dal costo di un prodotto ma­trice×vettore per ogni passo. Il metodo di Richardson, che risale al primo Novecento, ha il vantaggio di non richiedere speciali proprietà della matrice A; tuttavia, si dimostra che la velocità di convergenza dipende dal numero di condizionamento di A, un indice della sensibilità del problema Ax=b rispetto a errori sui dati, cioè sugli elementi di A e di b. L’indice di condizionamento dipende a sua volta dalla distribuzione degli autovalori di A e dalla scelta della base {φ_i}. Per realizzare il processo di aritmetizzazione in linea di fatto, cioè per ricondursi in modo effettivo a pure informazioni numeriche, occorre quindi poter analizzare la struttura (in questo caso le proprietà spettrali) di operatori lineari che intervengono nell’ultimo stadio della discretizzazione del modello differenziale.

Motori di ricerca su rete

Matrici di grandi dimensioni si usano regolarmente nei modelli per il funzionamento dei motori di ricerca su rete. I collegamenti che vengono realizzati da tali motori sono rappresentabili mediante un grafo orientato assieme alla sua matrice L di adiacenza. Più precisamente, se G={1, 2,…, n} è l’insieme dei nodi del grafo, si considera la matrice L i cui elementi L_ij valgono 1 se (i, j) è un arco orientato del grafo, e 0 altrimenti. La matrice L ci informa quindi circa l’esistenza di un collegamento dal nodo i al nodo j.

Sia ora deg(i) il numero degli archi che escono da i. Il numero deg(i) coincide con la somma degli elementi sulla riga i di L. Si consideri quindi la matrice P i cui elementi P_ij sono uguali a 1/deg(i) se L_ij=1, e a 0 altrimenti. Il numero (<1) 1/deg(i) si può interpretare come la probabilità che i si colleghi effettivamente a uno dei nodi j per cui risulti L_ij=1; sicché si presume che tutte queste probabilità siano uguali (indipendenti da j). Se non ha righe nulle, la matrice P, detta matrice di transizione, è stocastica per righe, cioè la somma degli elementi su ciascuna riga è 1. Si può allora definire intuitivamente il grado d’importanza di un generico nodo j, che deve essere proporzionale all’importanza dei nodi i da cui escono archi che arrivano a j, e inversamente proporzionale al numero complessivo deg(i) di archi che escono da ciascun nodo i. Infatti, se questo numero è grande, diminuisce la probabilità che, a un certo istante, si vada dal nodo i al nodo j. L’importanza p_j di j è definibile quindi come la somma dei quozienti p_i/deg(i) per tutti gli indici i per cui vale la relazione L_ij=1, tali cioè che ci sia collegamento da i a j. Nei termini della matrice di transizione l’importanza p_j risulta uguale alla somma dei prodotti P_ij p_i per tutti gli indici i per cui si abbia L_ij=1. In termini di matrici e vettori, se p è il vettore di elementi p_i, si ottiene quindi la formula p=P^Tp, da cui risulta che il vettore p è un punto fisso dell’operatore P^T.

Un teorema del matematico tedesco Oskar Perron (1880-1975) enunciato nel 1907 stabilisce che una matrice con elementi reali positivi ha un unico autovalore reale λ di modulo massimo. Il risultato si può estendere, come ha dimostrato nel 1912 il matematico tedesco Ferdinand G. Frobenius (1849-1917), a matrici irriducibili non negative, cioè con elementi reali non negativi e tali da non potersi ridurre, con permutazioni di righe e colonne, a una matrice divisibile in 4 blocchi di cui uno, sotto la diagonale principale, di elementi nulli (irriducibilità). Il caso della matrice stocastica P^T rientra nella teoria di Perron-Frobenius, in quanto P si può modificare in modo da ottenere una matrice irriducibile non negativa. L’autovalore di modulo massimo è uguale esattamente a 1 (mentre gli altri n−1 autovalori hanno modulo minore di 1). Se si chiama ancora P questa matrice modificata, la formula p=P^Tp dice che il vettore p è l’autovettore che corrisponde all’autovalore di modulo massimo di P^T.

La formula ottenuta si può interpretare come un modello del processo di ricerca su rete. Questo processo avviene lungo il grafo G. Se al passo k si sta alla pagina i, al passo k+1 ci si sposta in modo casuale a una delle pagine j per cui esiste un arco che congiunge i a j, cioè L_ij=1. Le pagine (nodi) possono considerarsi come gli n stati di una catena di Markov con matrice di transizione P: p_ij=1/deg(i) è la probabilità di collegarsi da i a j in un istante qualsiasi. Quindi, se p_i^(k) è la probabilità di trovarsi, al passo k, alla pagina i, si ha che p_i^(k+1⁾ è la somma dei quozienti p_i^(k)/deg(i) per tutti gli indici i per cui esiste un collegamento da i a j. Infatti, la probabilità che ci si trovi in un istante qualsiasi del processo alla pagina i è la somma delle probabilità che si attivi il collegamento (link) da i a j; mentre ciascuna di queste probabilità è data a sua volta dal prodotto della probabilità p_i^(k) di trovarsi alla pagina i per la probabilità 1/deg(i) di passare da i a j.

Si ha quindi, in termini matriciali, p^(k+1⁾=P^Tp^(k). Questa formula definisce a sua volta un algoritmo iterativo, noto come metodo delle potenze, che calcola, a ogni passo k, un’approssimazione dell’autovettore corrispondente all’autovalore di modulo massimo della matrice P^T. Il fatto che occorre appurare è la convergenza della successione di approssimazioni p^(k) alla soluzione, che è definita dal punto (vettore) fisso p. Ora è ben noto che questa proprietà di convergenza vale precisamente qualora vi sia un unico autovalore (reale) di modulo massimo, una condizione garantita, in questo caso, dalla teoria di Perron-Frobenius sulle matrici non negative.

La teoria di Perron-Frobenius rappresenta quindi una base teorica per ogni futura ricerca nel campo. Un fatto indicativo della natura dell’informatica, definita dall’informatico statunitense Donald Knuth come una scienza degli algoritmi. Per ironia della storia, Perron era alquanto scettico sulla consistenza della nuova scienza dei calcolatori e si chiedeva ancora, in una lettera del 1972, se l’informatica non fosse per caso un semplice sinonimo o un ramo secondario della ‘spionistica’ (Spionatik), riferendosi forse al lavoro di decrittazione dei messaggi militari tedeschi in cui si era impegnato Turing durante la Seconda guerra mondiale. Oltre che nella modellizzazione dei motori di ricerca, la teoria di Perron-Frobenius si è rivelata di fondamentale importanza nella risoluzione numerica di problemi differenziali, come i problemi al contorno per le equazioni di Laplace.

Metodi newtoniani, reti neuronali, potenza del calcolo

Nel calcolo su grande scala si eseguono operazioni aritmetiche approssimate tra numeri razionali. Tuttavia, i numeri non sono rapporti o frazioni p/q, con p e q interi, bensì numeri di macchina, ovvero numeri della forma .d1d2...d_t×B^p, dove B è la base della rappresentazione. Le frazioni p/q si prestano, infatti, all’inconveniente di una crescita eccessiva di p e di q, con conseguente instabilità del calcolo. I numeri di macchina sono finiti, con un numero finito di cifre, e non godono della maggior parte delle proprietà dell’aritmetica dei numeri reali.

Il fatto che il calcolo scientifico si basa su operazioni approssimate tra numeri di macchina suggerisce che l’errore non compromette, in linea di principio, l’efficienza delle procedure. Al contrario, poiché l’errore è inevitabile e un calcolo fondato sulle classiche proprietà del campo reale e su formule analitiche esatte non ha senso, si possono introdurre intenzionalmente errori di approssimazione allo scopo di semplificare il problema fino al punto da non compromettere il senso del risultato. Questa semplificazione consiste spesso nella seguente strategia: approssimare un problema privo di struttura, oppure con una struttura debole, con un altro di struttura più forte, da cui possa conseguire una semplificazione dei calcoli. È necessario però dimostrare che l’errore di approssimazione che ne segue non risulti fatale alla precisione richiesta.

Un esempio di questa strategia si ha nei metodi newtoniani per il calcolo del minimo non vincolato di una funzione f(x), dove x è un vettore di n componenti reali. Per calcolare un punto di minimo locale α si approssima la f con un suo modello lineare, definito dalla serie di Taylor troncata al termine di primo grado. Questa prima approssimazione consente di cercare il punto di minimo mediante il metodo di Newton definito dalla formula iterativa x_k+1=x_k−λH(x_k)⁻1g_k con un valore iniziale x0 prefissato, per k=0, 1, 2,…., dove H(x_k) è l’hessiano di f (la matrice il cui elemento i,j è la derivata parziale di f rispetto alla i-esima e alla j-esima componente di x), g_k è il gradiente di f calcolato in x_k (il vettore delle derivate parziali di f in x_k) e λ è un parametro opportuno che definisce la lunghezza del passo. Questo metodo converge velocemente ad α quando la distanza di x0 da α è sufficientemente piccola, tuttavia ha un grave inconveniente: l’assenza di struttura dell’hessiano, per cui a ogni passo si richiede un numero di operazioni proporzionale a n3 e la discesa della funzione non è garantita. Si evita l’inconveniente approssimando a ogni passo k l’hessiano con una matrice A_k, calcolata dalla precedente A_k−1 con una formula iterativa A_k=φ(A_k−1, s_k−1, y_k−1), scegliendo in modo opportuno una matrice iniziale A0, dove s_k=x_k−x_k−1 e y_k=g_k−g_k−1. Si ottengono così i metodi newtoniani, studiati fin dagli anni Settanta del 20° sec., che abbassano la complessità a O(n2) e fanno in modo che la matrice A_k abbia le proprietà necessarie per garantire la discesa della funzione. Ma un’ulteriore approssimazione è possibile: A_k si può approssimare a sua volta (grazie al teorema della proiezione in spazi di Hilbert) mediante un operatore lineare B_k appartenente a un’algebra di matrici simultaneamente diagonalizzabili con una trasformata veloce (come l’algebra delle matrici circolanti). Quest’ulteriore approssimazione non compromette, se si definisce in modo opportuno la direzione di discesa, la convergenza del metodo, e la minore velocità di avvicinamento ad α viene compensata dalla riduzione di complessità a O(nlogn) operazioni aritmetiche per passo. Inoltre, l’iterazione matriciale che permetteva di calcolare A_k diventa un’iterazione vettoriale sui soli autovalori delle matrici coinvolte, con una riduzione degli spazi di memoria da O(n2) a O(n). Questa riduzione della complessità consente di risolvere problemi, quali l’apprendimento automatico mediante reti neurali, altrimenti irrisolvibili a causa dell’elevata dimensione. La teoria delle reti neurali, nata da un lavoro pionieristico di Warren McCulloch e Walter Pitts (A lo­gical calculus of the ideas immanent in nervous activity, «Bulletin of mathematical biophysics», 1943, 5, pp. 115-33), si basa su un’idea di calcolo che dipende dall’interazione tra neuroni configurabili come nodi di un grafo, in cui a ogni cammino (sinapsi) tra il nodo i e il nodo j corrisponde un peso w_ij. Il modello è orientato a risolvere non tanto classici problemi computazionali (come addizionare o moltiplicare due numeri) quanto problemi di carattere cognitivo; ma può ricondursi a sua volta a problemi numerici, come il calcolo del minimo di una funzione. L’apprendimento automatico si basa, infatti, su una riduzione progressiva della discrepanza tra la risposta esatta e quella calcolata dalla rete (come nel caso del riconoscimento di una forma), e a questo scopo occorre definire metodi di discesa della funzione che esprime questa discrepanza, aggiornando a ogni passo i pesi che regolano le comunicazioni tra i nodi della rete. Le matrici che intervengono nel modello newtoniano sono spesso di dimensioni enormi, tali da rendere una procedura iterativa di complessità O(n2) per passo, pressoché inutilizzabile. La strategia migliore, applicabile anche in altri contesti, consiste allora in una combinazione di due criteri: possibile indebolimento delle condizioni che assicurano la convergenza del metodo e approssimazione di operatori privi di struttura con matrici diagonalizzabili con trasformate veloci. Gli errori che ne conseguono non sono fatali e sono compensati da strutture matematiche senza le quali il calcolo su grande scala non sarebbe nemmeno pensabile.

Bibliografia

Th. Kailath, A.H. Sayed, Displacement structure. Theory and applications, «SIAM Review», 1995, 37, 3, pp. 297-386.

R.H. Chan, M.K. Ng, Conjugate gradient methods for Toeplitz systems, «SIAM Review», 1996, 38, 3, pp. 427-82.

G. Strang, T. Nguyen, Wavelets and filter banks, Wellesley (Mass.) 1996.

R.B. Bapat, T.E.S. Raghavan, Nonnegative matrices and applications, Cambridge-New York 1997.

M. Bertero, P. Boccacci, Introduction to inverse problems in imaging, Bristol-Philadelphia 1998.

L. Blum, F. Cucker, M. Shub, S. Smale, Complexity and real computation, New York 1998.

J. Stoer, R. Bulirsch, Introduction to numerical analysis, New York 20023.

C. Di Fiore, S. Fanelli, F. Lepore, P. Zellini, Matrix algebras in quasi-Newton methods for unconstrained minimization, «Numerische mathematik», 2003, 94, 3, pp. 479-500.

A.N. Langville, C.D. Meyer, A survey of eigenvector methods for web information retrieval, «SIAM Review», 2005, 47, 1, pp. 135-61.