Esempio di utilizzo del linguaggio C per ilcalcolo scientificor.vincelli@campus.unimib.it

Introduzione

Parte dell’esame e consistita nella stesura di un programma per risolvere sistemilineari ed autosistemi su matrici date. E’ stata lasciata liberta di scelta riguardol’implementazione ma e stato posto il vincolo di non utilizzare un computeralgebra system (es. Wolfram Mathematica, MathWorks MATLAB...). La miascelta e ricaduta sul linguaggio C. Il motivo principale e stato il desiderio diriprenderlo ed approfondirlo ma e anche vero che e il linguaggio di scelta perprogrammi realtime e performanti; inoltre il linguaggio C e da tempo utilizzabileper la stesura di codice parallelo attraverso differenti flavors (es. MPI) e sipresta quindi naturalmente all’implementazione di molti algoritmi numerici incontesti di calcolo ad alte prestazioni (HPC). Mi sono attenuto in particolareallo standard (ANSI C) eccezion fatta ovviamente per l’adozione di una libreriascientifica specializzata, GNU GSL. GNU Scientific Library e sicuramente lalibreria libera di riferimento. Qualita sono, oltre al fatto di essere opensource emantenuta da esperti, il grande numero di algoritmi implementati, la presenza distrutture dati accessorie per la gestione di matrici, il che salva dal dovere gestirea mano la memoria, una documentazione chiara ed esaustiva. GNU GSL si basasu di una core library, CBLAS, che e un’implementazione C dell’interfacciaapplicativa standard di fatto per le operazioni base su matrici, BLAS (BasicLinear Algebra Subprograms).La maggior parte degli algoritmi di algebra lineare implementati in GNU GSLsono presi direttamente da [1].

Sistemi lineari

Un sistema lineare e un sistema di equazioni lineari di forma:

$'''&'''%

a1,1x1 � a1,2x2 � � � � � a1,nxn � b1a2,1x1 � a2,2x2 � � � � � a2,nxn � b2

...am,1x1 � am,2x2 � � � � � am,nxn � bm

In forma matriciale si puo scrivere:

�����a11 a12 � � � a1na21 a22 � � � a2n...

.... . .

...am1 am2 � � � amn

�����

�����x1x2...xn

����� �

�����b1b2...bm

�����

e, in modo piu compatto:

Ax � b

1

Algoritmi

Se e m � n la matrice A e quadrata ed il sistema ha una ed una sola soluzionese e solo se A e anche nonsingolare/invertibile. A lo e se esiste una matrice A�1

tale che AA�1 � Id. Se A e quadrata ed invertibile allora esistono sempre trematrici P , di permutazione sulle righe, L triangolare inferiore (lower), U trian-golare superiore (upper) tali che vale PA � LU . Se vogliamo scambiare le righei e j della matrice A prendiamo Id e scambiamo in essa queste: chiamata questanuova matrice P fare PA determina lo scambio di queste righe in A (duale perle colonne facendo AP ). Una matrice triangolare inferiore (superiore) e unamatrice quadrata che ha tutti gli elementi al di sopra (al di sotto) della dia-gonale principale a zero. Questo procedimento si chiama LUP decomposizioneoppure LU decomposizione con pivoting parziale. Ottenuta la tripla di matricipP, U, Lq il sistema si puo riscrivere:

LUx � Pb

e puo essere risolto mediante questi due passi:

1. risolvere Ly � Pb

2. risolvere Ux � y

Sono principalmente coinvolti quindi due algoritmi:

1. eliminazione gaussiana con pivoting parziale, per trovare pP, U, Lq

2. sostituzione avanti-indietro (forward and backward substitution) per risol-vere i due sistemi, essendo che L ed U sono triangolari

L’algoritmo di eliminazione gaussiana con pivoting parziale opera sulla matriceaumentata B e consiste nei seguenti passi:

1. scambiare la prima riga con la riga che ha il primo elemento massimo inmodulo della prima colonna, e salvare lo scambio (in caso che ce ne sianodi piu uguali, prendere la riga piu alta); se tutte le righe hanno il primoelemento nullo, andare al punto 3

2. per ogni riga Bi, i ¡ 1 con primo elemento non nullo, moltiplicarla perun coefficiente ki scelto in maniera tale che la differenza tra Bi e la primariga scalata per ki abbia il primo elemento nullo (il coefficiente e percioki � Bi1{B11); sostituire Bi con la somma appena ricavata ma porre, nellaprima posizione, ki invece che 0

3. adesso sulla prima colonna tutte le cifre, eccetto la prima, hanno ki; a que-sto punto ritornare al punto 1 considerando la sottomatrice che s’ ottienecancellando la prima riga e la prima colonna

4. terminata questa procedura, si scarta l’ultima colonna e si fa una secondacopia di questa nuova matrice A1; in una poniamo la diagonale ugualead 1, sopra di essa 0 e teniamo invariata la parte inferiore, ottenendoL; nell’altra mettiamo a 0 la parte inferiore e lasciamo invariato il resto,ottenendo U

2

5. P e facilmente costruibile a partire dagli scambi salvati: per ogni coppiasi scambiano le rispettive righe in Id

Esempio. Consideriamo il sistema:

$&%

x1 � 2x2 � x3 � 24x1 � 3x2 � x3 � 32x1 � 2x2 � 3x3 � 5

La forma matriciale e: ��1 2 �1

4 3 12 2 3

����x1x2x3

�� �

��2

35

��

e la matrice aumentata: ��1 2 �1 2

4 3 1 32 2 3 5

��

il pivot sarebbe 1 ma il massimo in modulo della prima colonna e 4 quindiabbiamo lo scambio p1, 2q ed otteniamo:

��4 3 1 3

1 2 �1 22 2 3 5

��

il primo coefficiente e k1 � 1{4 e la matrice diventa:

�� 4 3 1 3k1 5{4 �5{4 5{42 2 3 5

��

e con il secondo k2 � 1{2 si ha:

�� 4 3 1 3k1 5{4 �5{4 5{4k2 1{2 5{2 7{2

��

il che porta a ridurre il problema alla sottomatrice

�5{4 �5{4 5{41{2 5{2 7{2

�

per la quale non sono necessari scambi ed il coefficiente risulta k3 � 2{5 con cuisi ha �

5{4 �5{4 5{4k3 3 3

�

e l’eliminazione termina; la matrice finale e

�� 4 3 1 3k1 5{4 �5{4 5{4k2 k3 3 3

��

3

mentre la triangolare L e �� 1 0 0

1{4 1 01{2 2{5 1

��

accompagnata dalla triangolare U��4 3 1

0 5{4 �5{40 0 3

��

ed infine P e ��0 1 0

1 0 00 0 1

��

Come ci si aspetta vale PA � LU , infatti LU da��4 3 1

1 2 �12 2 3

��

e calcolare PA significa scambiare la prima con la seconda riga di A.

L’algoritmo di sostituzione avanti-indietro ha nel nome quel che in pratica vafatto per risolvere i due sistemi. Nel risolvere quello su L l’incognita y1 e l’unicanon nulla e le successive si ottengono operando una sola sostituzione e svolgendoi calcoli; per quello su U e duale, e l’incognita xn ad essere isolata, e si sostitui-sce non in avanti ma all’indietro.Proseguendo l’esempio, il primo sistema da risolvere e in avanti, Ly � Pb�

� 1 0 01{4 1 01{2 2{5 1

����y1y2y3

�� �

��3

25

��

ed ha soluzione y � p3, 5{4, 3q, mentre il secondo Ux � y si concretizza in��4 3 1

0 5{4 �5{40 0 3

����x1x2x3

�� �

�� 3

5{43

��

che ha soluzione x � p�1, 2, 1q, che e la soluzione finale del sistema.

Autosistemi

Data una matrice quadrata A risolvere l’autosistema (eigensystem) significa tro-varne autovettori ed autovalori. v � 0 e un autovettore con associato autovaloreλ per la matrice A se vale la seguente:

Av � λv

La libreria usa due metodi diversi per il calcolo degli autovalori a seconda chesi tratti di una matrice reale simmetrica o nonsimmetrica.Se A e simmetrica la procedura e incentrata sui seguenti due passi:

4

• calcolo di T , la tridiagonalizzazione di A

• esecuzione dell’algoritmo QR simmetrico su T

Se A invece e nonsimmetrica:

• calcolo di H, la forma di Hessenberg superiore di A

• esecuzione dell’algoritmo QR con passo di Francis

Una matrice tridiagonale come suggerisce il nome e nulla in tutte le posizioninon appartenenti alla diagonale ed alle diagonali inferiore e superiore. Essapuo essere ottenuta tramite la trasformazione di Householder. Una matrice siHessenberg e una matrice tridiagonale superiore, nel senso che e nulla solamen-te al di sotto della diagonale inferiore. Il passaggio a queste forme ha comeobiettivo quello di minimizzare l’informazione portata avanti, infatti salviamosolo le diagonali, ottenendo pero una matrice dagli stessi autovalori, poiche so-no simili. Due matrici M ed N sono simili se esiste una P invertibile tale cheM � P�1NP ; informalmente due matrici simili sono la stessa trasformazione egli autovalori coincidono (in generale gli autovettori no).La QR decomposizione di una matrice A e il calcolo di una coppia di matricipQ,Rq con Q ortogonale ed R triangolare superiore tali che A � QR. Unamatrice ortogonale e una dove i vettori colonna e riga sono un insieme di vet-tori ortonormali, ossia a norma unitaria e perpendicolari tra loro (due vettorilo sono sse il loro prodotto scalare e nullo). Poniamo A0 � T ; l’iterazione tdell’algoritmo QR consiste nel porre At � Rt�1Qt�1 dove pQt�1, Rt�1q e laQR decomposizione di At�1. Ad ogni iterazione la diagonale di At contieneun’approssimazione degli autovalori di A. Questo e vero poiche la successionetAtu converge alla forma di Schur di A; essa e la matrice triangolare superioreU nell’omonima decomposizione assieme alla matrice unitaria S per cui valeA � SUS�1. Una matrice unitaria e una che se moltiplicata a destra per lapropria trasposta coniugata da Id; la trasposta coniugata si ottiene prendendocome valori i coniugati complessi della trasposta, cioe negando le parti immagi-narie ma non reali dei valori. Ora U e simile ad A e la relazione e simmetricaquindi vale anche U � SAS�1; possiamo allora affermare che A ed U hanno glistessi autovalori e siccome U e triangolare essi si trovano sulla sua diagonale.Le implementazioni della libreria non incarnano la versione naive dell’algoritmoQR presentata, ed anzitutto hanno una chiara condizione di terminazione chegarantisce stabilita (vedi dopo) rispetto ad una tolleranza data. Una differenzatra la versione simmetrica e nonsimmetrica e la complessita computazionale,piu chiaramente determinabile ed inferiore per la prima. Per i dettagli delleprocedure il riferimento sono i capp. 7 ed 8 di [1].Esempio. I seguenti comandi MATLAB permettono di provare la bonta dellaprocedura ingenua, in modo diretto e comodo.

• A=[1 2 -1; 4 3 1; 2 2 3]); creazione matrice (e nonsimmetrica)

• H=hess(A) trasformazione di Hessenberg

• [Q R]=qr(H); H=RQ; iterazione; basta eseguirla piu volte per proce-dere con l’algoritmo

Gia dopo una ventina di iterazioni si ha la soluzione esatta.

5

Considerazioni numeriche

Stabilita

Consideriamo l’algoritmo una funzione f : X Ñ Y che associa ad ogni istanzadi un problema x P X la sua soluzione y P Y ; questa definizione e di uso comunein informatica teorica. In analisi numerica identifichiamo due tipi di algoritmi,diretti ed iterativi. I primi calcolano fpxq in un numero finito di passi, mentreper i secondi non ci si puo aspettare la terminazione in un numero finito dipassi, convergono con una certa velocita alla soluzione esatta, che raggiungonoal limite. Gli algoritmi iterativi sono i piu diffusi. L’eliminazione gaussiana e unalgoritmo diretto mentre la procedura QR e iterativa. Per gli algoritmi iterativila terminazione si ha raggiunta una certa precisione ossia la garanzia che la so-luzione calcolata non si discosta per piu di un parametro ε dato dall’esatta. Lasoluzione calcolata dai metodi diretti e precisa in un contesto di calcolo ideale,ma non e cosı per via di come sono memorizzati i numeri al calcolatore. D’altrocanto per la stessa ragione la convergenza degli algoritmi iterativi puo arrestarsiad un valore che puo anche essere quello esatto: la differenza e inferiore a quantopuo essere precisa la macchina e quindi risultano uguali.L’aritmetica a virgola mobile e adottata nella stragrande maggioranza dei cal-colatori usati per applicazioni scientifiche, e per sua stessa natura e fonte diimprecisioni negli algoritmi. Le principali cause di queste imprecisioni sono l’er-rore inerente/di arrotondamento, i numeri periodici ed irrazionali sono forzati aduna rappresentazione finita nel calcolatore, e l’errore di troncamento, una som-ma infinita e arrestata. Queste imprecisioni sono a loro volta identificabili comeerrori commessi dall’algoritmo, determinati appunto da scenari necessari percheinerenti all’insieme numerico a virgola mobile. In algebra lineare numerica diinteresse sono gli errori:

• in avanti: la differenza tra il risultato dell’algoritmo e la soluzione vera;chiamando y� il primo e � y� � y (forward stability)

• indietro: la minima quantita da aggiungere all’input perche la soluzionecalcolata eseguendo l’algoritmo su di esso coincida con la soluzione ve-ra dell’algoritmo eseguito sull’input sommato alla suddetta quantita; e ilminimo ∆x tale che fpx� ∆xq � y� (backward stability)

Essendo ∆x � x� � x l’errore all’indietro ci dice quale istanza del problema x�

ha soluzione vera y�. Queste quantita sono riferite ad errori assoluti ma, speciese l’esecuzione tipica lavora su numeri molto grandi, si declinano sull’errorerelativo; la prima diventa |y��y|{|y| e nella seconda sostituiamo a ∆x |∆x|{|x|.Un algoritmo e stabile rispetto ad uno od all’altro errore se le rispettive quantitasono piccole. Una richiesta tipica e che coincida o almeno sia dello stesso ordine-esponente della precisione dell’aritmetica a virgola mobile sottostante. Questaquantita e comunemente chiamata epsilon della macchina; in ANSI C sono postea 10�5 e 10�9 almeno per il tipo float e per il tipo double rispettivamentee sono precalcolate ed accessibili dall’intestazione float.h. Nella pratica unalgoritmo si dice stabile se e stabile in un’accezione a meta strada tra le due, seesiste un ∆x piccolo tale che fpx�∆xq�y� sia piccolo pure ((mixed) stability),intuitivamente se risolve con sufficiente precisione da un’istanza sufficientementevicina alla desiderata. Questa stabilita e implicata dalla stabilita all’indietro.

6

Stabilita degli algoritmi adoperati

L’eliminazione gaussiana con pivoting parziale non e in generale stabile all’in-dietro. Se pero A e una matrice simmetrica e definita positiva (le matrici easy,bad e verybad nel test) lo e. D’altro canto nelle applicazioni non si incon-trano matrici problematiche; in [2] e presentato un esperimento con il quale siargomenta che, generato casualmente un miliardo di matrici, la probabilita cheanche solo una determini un’istanza d’esecuzione instabile e inferiore a .01.Per quanto riguarda la procedura di calcolo degli autovalori di una matrice sim-metrica, per la tridiagonalizzazione esistono algoritmi stabili (vedi per es. [3])e l’algoritmo QR e argomentato essere stabile in [4] poiche si basa sulla QRdecomposizione che e una trasformazione di similarita ortogonale. Piu in det-taglio la backward stability di entrambe le varianti dell’algoritmo QR risiede nelfatto che e per x � A e calcolata la soluzione esatta dell’istanza x� � A � Econ ||E||2 � u||A||2 dove u e l’epsilon della macchina e || � ||2 e la p-normaeuclidea per matrici. Essendo la matrice A quadrata e chiamata norma spet-trale e coincide con l’autovalore massimo di A�A con A� coniugata traspostadi A (gli autovalori di A�A sono detti valori singolari di A). Questo valore ecomunemente molto piccolo e giustifica quindi l’affermazione di stabilita.

7

Bibliografia

[1] Gene H. Golub and Charles F. Van Loan. Matrix Computations. The JohnsHopkins University Press, 3rd edition, 1996.

[2] Lloyd N. Trefethen and David Bau. Numerical Linear Algebra. SIAM:Society for Industrial and Applied Mathematics, 1997.

[3] Symmetric Matrices Chengshu, Chengshu Guo, and Sanzheng Qiao. A sta-ble lanczos tridiagonalization of complex. In Advanced Signal ProcessingAlgorithms, Architectures, and Implementations XV, Proc. of SPIE, page591010, 2005.

[4] Alston S. Householder and Friedrich L. Bauer. On certain methods forexpanding the characteristic polynomial. j-NUM-MATH, 1:29–37, 1959.

8