Matlab: complessità e stabilità degli algoritmi. Alcuni esempi.

Ángeles Martínez Calomardo e Alvise Sommariva

Università degli Studi di Padova

6 novembre 2013

Esempio di algoritmo instabile

Formula risolutiva dell’equazione di secondo grado

Dato x^{2} + 2 px − q, con p^{2} + q ≥ 0 eseguiamo un primo algoritmo Matlab che valuta la radice via:

y = −p + \sqrt{p^{2} + q}.  (1)

p^{2} + q ≥ 0 implica radici reali.

Potenzialmente instabile per p ≫ q a causa della sottrazione tra p e \sqrt{p^{2} + q} (cancellazione).

Valutiamo la radice con un secondo algoritmo stabile via razionalizzazione di (1):

y = −p + \sqrt{p^{2} + q} = \frac{(−p + \sqrt{p^{2} + q})(p + \sqrt{p^{2} + q})}{(p + \sqrt{p^{2} + q})}
= \frac{q}{(p + \sqrt{p^{2} + q})}  (2)

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  2/ 37

Codice stabilità : algoritmo 1

Salviamo il seguente codice in radicesecgrado.m.

p =1000; q =0.018000000081; sol =0.9∗10ˆ ( −5) ;
% ALGORITMO 1
s=p ˆ 2 ;

t=s+q ;

i f t >=0
u=s q r t ( t ) ;

    s1=−p+u ;

e l s e

f p r i n t f ( ’ \n \ t [ RADICI COMPLESSE] ’ ) ;
s1 = [ ] ;

end

Codice stabilità : algoritmo 2

% ALGORITMO 2
s=p ˆ 2 ;

t=s+q ;

i f t >=0
u=s q r t ( t ) ;

v=p+u ;
t1=q/v ;

e l s e

f p r i n t f ( ’ \n \ t [ RADICI COMPLESSE] ’ ) ;
t1 = [ ] ;

end

Codice stabilità : stampa risultati

f p r i n t f ( ’ \n \ t [ ALG . 1 ] : %10.19 f ’ , s1 ) ;
f p r i n t f ( ’ \n \ t [ ALG . 2 ] : %10.19 f ’ , t1 ) ;
i f l e n g t h ( sol ) > 0 & ( sol ˜ = 0) & l e n g t h ( s1 ) > 0 & …
l e n g t h ( t1 ) > 0
rerr1 =abs ( s1−sol ) / abs ( sol ) ;
rerr2=abs ( t1−sol ) / abs ( sol ) ;
f p r i n t f ( ’ \n \ t [ REL . ERR . ALG . 1 ] : %2.2 e ’ , rerr1 ) ;
f p r i n t f ( ’ \n \ t [ REL . ERR . ALG . 2 ] : %2.2 e ’ , rerr2 ) ;

end

Test.

Come previsto, il secondo algoritmo si comporta notevolmente meglio del primo, che compie un errore relativo dell’ordine di circa 10^{−9}. Infatti:

>> radicesecgrado

[ ALG . 1 ] [ 1 ] : 0.0000089999999772772
[ ALG . 2 ] [ 1 ] : 0.0000090000000000000
[ REL . ERR . ] [ ALG . 1 ] : 2 . 5 2 e−009
[ REL . ERR . ] [ ALG . 2 ] : 0 . 0 0 e+000

>>

Calcolo π

Eseguiamo un codice Matlab che valuti le successioni {u_{n}}, {z_{n}}, definite rispettivamente come

s_{1} = 1, s_{2} = 1 + \frac{1}{4}

u_{1} = 1, u_{2} = 1 + \frac{1}{4}
s_{n+1} = s_{n} + \frac{1}{(n+1)^{2}}
u_{n+1} = \sqrt{6 s_{n+1}}

e

\{ \begin{matrix}z_{1} = 1, z_{2} = 2 \\ z_{n+1} = 2^{n− \frac{1}{2}} \sqrt{1 − \sqrt{1 − 4^{1−n} · z_{n}^{2}}} \end{matrix}  (3)

che teoricamente convergono a π.

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  7/ 37

Calcolo π

Implementiamo poi la successione, diciamo {y_{n}}, che si ottiene razionalizzando (3), cioè moltiplicando numeratore e denominatore di

z_{n+1} = 2^{n− \frac{1}{2}} \sqrt{1 − \sqrt{1 − 4^{1−n} · z_{n}^{2}}}

per

\sqrt{1 + \sqrt{1 − 4^{1−n} · z_{n}^{2}}}

e calcoliamo u_{m}, z_{m} e y_{m} per m = 2, 3, …, 40 (che teoricamente dovrebbero approssimare π). Infine disegniamo in un unico grafico l’andamento dell’errore relativo di u_{n}, z_{n} e y_{n} rispetto a π aiutandoci con l’help di Matlab relativo al comando semilogy.

Calcolo π: metodo 1

In seguito scriviamo un’implementazione di quanto richiesto commentando i risultati. Si salvi in un file pigreco.m il codice

% SEQUENZE CONVERGENTI ” PI GRECO” .
% METODO 1 .
s ( 1 ) =1; u ( 1 ) =1;
s ( 2 ) =1.25; u ( 2 )=s ( 2 ) ;
f o r n =2:40
s ( n+1)=s ( n )+(n+1)ˆ ( −2) ;
u ( n+1)=s q r t (6∗s ( n+1) ) ;

end

rel_err_u=abs ( u−p i ) / p i ;
f p r i n t f ( ’ \n ’ ) ;

Calcolo π: metodo 2

% METODO 2 .
format long
z ( 1 ) =1;
z ( 2 ) =2;
f o r n =2:40
c=(4ˆ (1−n ) ) ∗ ( z ( n ) ) ˆ 2 ; inner_sqrt=s q r t (1−c ) ;
z ( n+1)=(2ˆ (n −0.5) ) ∗ s q r t ( 1−inner_sqrt ) ;

end

rel_err_z=abs ( z−p i ) / p i ;
f p r i n t f ( ’ \n ’ ) ;

Calcolo π: metodo 3

% METODO 3 .
y ( 1 ) =1;
y ( 2 ) =2;
f o r n =2:40
num =(2ˆ (1/2) ) ∗ abs ( y ( n ) ) ;
c=(4ˆ (1−n ) ) ∗ ( z ( n ) ) ˆ 2 ;

inner_sqrt=s q r t (1−c ) ;
den=s q r t ( 1+inner_sqrt ) ;
y ( n+1)=num / den ;

end

rel_err_y=abs ( y−p i ) / p i ;

Calcolo π: plots

% SEMILOGY PLOT.
h o l d on ;
s e m i l o g y ( 1 : l e n g t h ( u ) , rel_err_u , ’ k . ’ ) ;
s e m i l o g y ( 1 : l e n g t h ( z ) , rel_err_z , ’m+ ’ ) ;
s e m i l o g y ( 1 : l e n g t h ( y ) , rel_err_y , ’ ro ’ ) ;
h o l d off ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> pigreco

Plot risultati

Figura: Errore relativo commesso con le 3 successioni, rappresentate rispettivamente da ., + e o.

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  13/ 37

Discussione risultati.

La prima successione converge molto lentamente a π, la seconda diverge mentre la terza converge velocemente a π.

Per alcuni valori {z_{n}} e {y_{n}} coincidono per alcune iterazioni per poi rispettivamente divergere e convergere a π. Tutto ciò è naturale poichè le due sequenze sono analiticamente (ma non numericamente) equivalenti.

Dal grafico dell’errore relativo, la terza successione, dopo aver raggiunto errori relativi prossimi alla precisione di macchina, si assesta ad un errore relativo di circa 10^{−15} (probabilmente per questioni di arrotondamento).

L’algoritmo 2 in dettaglio

Successione approssimante π

Nell’approssimare il valore di π con la seguente formula ricorsiva

z_{2}  = 2
z_{n+1} = 2^{n−0.5}\sqrt{1 − \sqrt{1 − 4^{1−n}z_{n}^{2}}}, n = 2, 3, …,

si ottiene la seguente successione di valori (dove si è posto c = 4^{1−n}z_{n}^{2}).

n + 1  c  1 − \sqrt{1 − c}  z_{n+1}  \frac{|zn+1−π|}{π}
⋯  ⋯  ⋯  ⋯  ⋯
10 1.505e-04 7.529e-05 3.14157294036 6.27e-06
11 3.764e-05 1.882e-05 3.14158772527 1.57e-06
12 9.412e-06 4.706e-06 3.14159142150 3.92e-07
13 2.353e-06 1.176e-06 3.14159234561 9.80e-08
14 5.882e-07 2.941e-07 3.14159257654 2.45e-08
15 1.470e-07 7.353e-08 3.14159263346 6.41e-09
16 3.676e-08 1.838e-08 3.14159265480 3.88e-10
17 9.191e-09 4.595e-09 3.14159264532 2.63e-09
18 2.297e-09 1.148e-09 3.14159260737 1.47e-08
19 5.744e-10 2.872e-10 3.14159291093 8.19e-08
⋯ ⋯ ⋯ ⋯ ⋯
28 2.220e-15 1.110e-15 3.16227766016 6.58e-03
29 5.551e-16 3.330e-16 3.46410161513 1.03e-01
30 1.665e-16 1.110e-16 4.00000000000 2.73e-01
31 5.551e-17 0.000e+00 0.00000000000 1.00e+00
32 0.000e+00 0.000e+00 0.00000000000 1.00e+00

Una successione ricorrente.

Consideriamo la successione {I_{n}} definita da

I_{n} = e^{−1} ∫_{0}^{1} x^{n} e^{x} dx  (4)

n = 0: I_{0} = e^{−1} ∫_{0}^{1} e^{x} dx = e^{−1}(e^{1} − 1). integrando per parti

I_{n+1} = e^{−1} (x^{n+1} e^{x} |_{0}^{1} −(n + 1) ∫_{0}^{1} x^{n} e^{x} dx)
= 1 − (n + 1) I_{n}.

I_{n} > 0, decrescente (x^{n+1} < x^{n} per x ∈ (0, 1)) e si prova che I_{n} → 0 come 1/n.

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  16/ 37

Problema.

Calcoliamo I_{n} per n = 1, …, 99:

mediante la successione in avanti

\{ \begin{matrix}I_{0} = e^{−1}(e^{1} − 1) \\ I_{n+1} = 1 − (n + 1) I_{n}. \end{matrix}  (5)

mediante la successione all’indietro

\{ \begin{matrix}t_{1000} = 0 \\ t_{n−1} = (1 − t_{n})/n. \end{matrix}

Si noti che se I_{n+1} = 1 − (n + 1) I_{n} allora I_{n} = (1 − I_{n+1})/(n + 1) e quindi I_{n−1} = (1 − I_{n})/n.

Successione ricorrente in Matlab

Scriviamo il codice in un file succricorrente.m.

% SUCCESSIONE RICORRENTE .
c l e a r a l l ;
% SUCCESSIONE ” s n ” .
s ( 1 )=exp (−1) ;
f o r n =1:99
s ( n+1)=1−(n+1)∗s ( n ) ;

end

% SUCCESSIONE ” t n ” .
M =1000;
t=z e r o s ( M , 1 ) ; % INIZIALIZZAZIONE ” t ” .
f o r n=M : −1:2
j=n −1;
t ( j )=(1−t ( n ) ) /n ;

end

Successione ricorrente in Matlab

% PLOT SEMI−LOGARITMICO .
c l f ;
h o l d on ;
s e m i l o g y ( 1 : l e n g t h ( s ) , abs ( s ) , ’ k− ’ ) ;
s e m i l o g y ( 1 : l e n g t h ( s ) , abs ( t ( 1 : l e n g t h ( s ) ) ) , ’m− ’ ) ;
h o l d off ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> succricorrente

Plot risultati

Figura: Grafico che illustra i valori assoluti assunti dalla successione in avanti (in nero) e all’indietro (in rosa magenta).

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  20/ 37

Instabilità della formula ricorsiva

La formula I_{n} = 1 − n I_{n−1} è instabile, quindi amplifica l’errore ad ogni passo.

Infatti, nel calcolatore

(I_{n} + ε_{n}) = 1 − n(I_{n−1} + ε_{n−1}).

Sottraendo dalla precedente equazione la relazione I_{n} = 1 − nI_{n−1} si può quantificare l’errore:

ε_{n} = −n ε_{n−1}, e per induzione |ε_{n}| = n! |ε_{0.}|

Il fattore n! amplifica l’errore di rappresentazione iniziale (su I_{0}), ε_{0}.

Esempio. Nel calcolo di I_{20} l’errore è ε_{20} = 20! ε_{0} ≈ 2.7 · 10^{2}.

Formula alternativa stabile

La formula all’indietro smorza l’errore

Per l’errore al passo n − 1 si trova

ε_{n−1} = \frac{−1}{n} ε_{n}.

Partendo da m

|ε_{m−1}| = \frac{|ε_{m}|}{m}, |ε_{m−2}| = \frac{|ε_{m}|}{m(m − 1)}, …, |ε_{m−k}| = \frac{|ε_{m}|}{m(m − 1) ⋯ (m − k + 1)}.

La produttoria al denominatore abbatte rapidamente l’errore iniziale! Per esempio, per calcolare I_{25} partendo da I_{40} = 0.5, l’errore iniziale |ε_{40}| < 0.5 verrebbe abbattuto di un fattore

40 · 39 ⋯ 27 · 26 = 5.2602 · 10^{22}

.

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  22/ 37

Complessità : potenza di matrice.

Problema: calcolare la potenza p-esima di una matrice quadrata A di ordine n cioè

A^{p} := \underbrace{A ∗ … ∗ A}_{p volte}

senza usare l’operatore di elevamento a potenza ∧.

Primo algoritmo. Si può implementare il seguente pseudocodice

B=I ;

f o r i=1:p

     B=B∗A ;

end

in cui I è la matrice identica di ordine n e ∗ è il classico prodotto tra matrici.

Secondo algoritmo. Alternativamente (in maniera più stabile ed efficiente) si può decomporre p come

p = ∑_{i=0}^{M}c_{i}2^{i}

ove M = ⌊log_{2} p⌋ e c_{i} = 0 oppure c_{i} = 1. Si osserva facilmente che questa è la classica rappresentazione di p in base 2. Usando la proprietà della potenze

B = A^{p} = A^{∑_{i=0}^{M} c_{i}2^{i}} = A^{∑_{i=0}^{M} 2^{i}c_{i}}
= (A^{2^{0}})^{c_{0}} ∗ … ∗ (A^{2^{M}})^{c_{M}} = ∏_{i=0}^{M}(A^{2^{i}})^{c_{i}}  (6)

ove ogni termine A^{2^{i}} può essere calcolato come A^{2^{i−1}}A^{2^{i−1}}.

Per capire il funzionamento dei metodi, introduciamo un ulteriore esempio esplicativo. Sia A una matrice di ordine 7 e p = 22. Nel primo algoritmo si calcola A^{p} come

A^{p} = A ∗ ⋯ ∗ A

e quindi sono necessari 21 prodotti tra matrici. Nel secondo algoritmo, essendo p = 22, Il valore assunto da c è quindi [0 1 1 0 1] poiché, come si legge dalla rappresentazione binaria (da sinistra verso destra),

22 = 0 · 2^{0} + 1 · 2^{1} + 1 · 2^{2} + 0 · 2^{3} + 1 · 2^{4}.

Notiamo che il vettore c ha lunghezza M + 1 = floor(log2(22)) + 1 = 5. Vogliamo quindi calcolare

B = (A)^{0} ∗ (A^{2})^{1} ∗ (A^{4})^{1} ∗ (A^{8})^{0} ∗ (A^{16})^{1}.  (7)

Calcoliamo B = (A)^{0} ∗ (A^{2})^{1} ∗ (A^{4})^{1} ∗ (A^{8})^{0} ∗ (A^{16})^{1}. Al primo passo, dobbiamo calcolare B = (A)^{0}. Poniamo C = A. Come suggerito, se inizialmente B = I_{7}, a questo punto non c’è nulla da fare poichè A ha ordine 7 e A^{0} = I_{7}. A questa altezza abbiamo immagazzinato B, C ed A. Al secondo passo, ricordando (7), dobbiamo moltiplicare il valore di B ottenuto al passo precedente per (A^{2})^{1}. Essendo C = A osserviamo che A^{2} = A ∗ A = C ∗ C e quindi posto il nuovo C uguale a C ∗ C (per cui C = A^{2}) si ha B = B ∗ C. A questa altezza abbiamo immagazzinato B, C ed A. Al terzo passo, ricordando (7), dobbiamo moltiplicare il valore di B ottenuto al passo precedente per (A^{4})^{1}. Essendo C = A^{2} osserviamo che A^{4} = A^{2} ∗ A^{2} = C ∗ C e quindi posto il nuovo C uguale a C ∗ C (per cui C = A^{4}) si ha B = B ∗ C. A questa altezza abbiamo immagazzinato B, C ed A.

Al quarto passo, ricordando (7), dobbiamo moltiplicare il valore di B ottenuto al passo precedente per (A^{8})^{0}. Essendo C = A^{4} osserviamo che A^{8} = A^{4} ∗ A^{4} = C ∗ C. Poniamo il nuovo C uguale a C ∗ C (e quindi C = A^{8}). A questo punto non c’è nulla da fare poichè (A^{8})^{0} = I_{7} e quindi B = B ∗ I_{7}. A questa altezza abbiamo immagazzinato B, C ed A. Al passo M = 5, ricordando (7), dobbiamo moltiplicare il valore di B ottenuto al passo precedente per (A^{16})^{1}. Essendo C = A^{8} osserviamo che A^{16} = A^{8} ∗ A^{8} = C ∗ C e quindi posto il nuovo C uguale a C ∗ C (per cui C = A^{16}) si ha B = B ∗ C. A questa altezza abbiamo immagazzinato B, C ed A.

Notiamo, che una volta finito il processo, A^{16} è immagazzinata in B.

Complessità : potenza di matrice.

Un codice che produce la decomposizione in potenze di 2 di un numero p è il seguente:

q=p ;

M=f l o o r ( l o g 2 ( p ) ) +1;
c = [ ] ;
f o r i =1:1: M
c ( i )=mod ( q , 2 ) ;
q=f l o o r ( q /2) ;

end

Complessità : potenza di matrice.

Uno pseudocodice che implementa il secondo algoritmo è

p =100; n =200;

c=trasforma_in_binario ( p ) ;

A=rand ( n ) ;
B=eye ( n ) ;

C=A ;

M=f l o o r ( l o g 2 ( p ) ) ;
% B c o n t i e n e l a potenza d i A f i n o r a c a l c o l a t a .
% C c o n t i e n e l a potenza Aˆ (2ˆ i n d e x ) ) f i n o r a c a l c o l a t a .

f o r index =0:M
    j=index +1;

i f c ( j ) = 1

       B=B∗C ;

end

    C=C∗C ;

end

Esercizio 1

Si implementino i due algoritmi proposti per il calcolo della potenza di matrice tramite due functions (senza usare l’operatore ∧) e si calcoli l’errore relativo in norma infinito rispetto all’elevamento a potenza di Matlab/Octave per diverse matrici e potenze (n = 25, 50 e p = 20, 40).

Esercizio 2

Si confrontino poi i tempi di esecuzione delle due functions per il calcolo di A^{100}, con A matrice di numeri casuali di dimensione 200 × 200.

Norme matriciali

Nell’introduzione a Matlab/Octave abbiamo visto che il comando norm(v) calcola la norma del vettore v.

Analogamente, per calcolare in Matlab/Octave la norma di una matrice A si usa il comando norm(A).

Se nessun ulteriore parametro viene specificato tale comando restituisce la norma 2 della matrice ovvero:

∥A∥_{2} = \sqrt{ρ(A^{T}A)}

dove ρ(A) è il raggio spettrale della matrice A. Altre possibilità sono:

norma 1

∥A∥_{1} = max_{j} ∑_{i=1}^{n} |a_{ij}|, in Matlab/Octave norm(A,1); norma infinito ∥A∥_{∞} = max_{i} ∑_{j=1}^{n} |a_{ij}|, in Matlab/Octave norm(A,inf); norma di Frobenius ∥A∥_{F} = \sqrt{∑_{i=1}^{n} ∑_{j=1}^{n} |a_{ij}|^{2}}, in Matlab/Octave

    norm(A,'fro');

Ángeles Martínez Calomardo e Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  31/ 37

Norme matriciali

Esempi

>>> A = [ 5 −4 2 ; 1 7 −6; 1 1 9 ]

A =

5 −4 2
1  7 −6
1  1  9
>>> norm ( A , 1 )
ans = 17
>>> norm ( A , inf )
ans = 14
>>> norm ( A , ’ f r o ’ )
ans = 14.6287388383278
>>> norm ( A )
ans = 12.0560586095913

Come misurare la durata di un programma

Per confrontare due programmi che risolvono lo stesso problema è utile misurare il tempo di CPU (wallclock time) impiegato per eseguirli. In Matlab/Octave questo tempo si misura in secondi con il comando: cputime.

Esempio:

>> A = rand (5000) ;
t = cputime ;
det ( A ) ;
tfin=cputime ;

cpu=tfin−t

>> cpu = 32.500

Facoltativo. Complessità : algoritmo di Horner.

Ci poniamo il problema di valutare il polinomio

p(x) = a_{0} + a_{1} · x + … + a_{n} · x^{n}  (8)

in un punto x. Osserviamo che

p(x) = a_{0} + x · (a_{1} + x · (a_{2} + … + x · (a_{n−1} + x · a_{n}))) (9)

Supponiamo sia a = (a_{0}, …, a_{n}) il vettore di dimensione n + 1 delle componenti del polinomio. Possiamo valutare il polinomio tramite i seguenti due algoritmi, il primo che valuta direttamente il polinomio secondo quanto descritto in (8), il secondo che effettua la stessa operazione come descritto in (9) calcolando dapprima s_{1} = a_{n−1} + x · a_{n}, poi s_{2} = a_{n−2} + x · s_{1} e così via.

Facoltativo: Complessità : algoritmo di Horner.

In Matlab avremo allora

f u n c t i o n s=algoritmo1 ( a , x )

xk =1; s=a ( 1 ) ;
f o r i=2: l e n g t h ( a )

   xk=xk∗x ;
   s=s+a ( i ) ∗xk ;

end

e

f u n c t i o n s=algoritmo2 ( a , x )

L=l e n g t h ( a ) ;
s=a ( L ) ; % COMPONENTE a n IMMAGAZZINATA IN a ( n+1) .
f o r i=L −1: −1:1
s=a ( i )+x∗s ;

end

Facoltativo: Complessità : algoritmo di Horner.

Se lanciamo il codice demo horner per la valutazione di p(x) = 1 + 2 · x + 3 · x^{2} + 4 · x^{3} in x = π

c l e a r a l l ;
a=[1 2 3 4 ] ;
x=p i ;

y1=algoritmo1 ( a , x ) ;
y2=algoritmo2 ( a , x ) ;

format long ;

y1
y2

otteniamo

>> demo_horner

ans = 1.609171052316469 e+02
y1 = 1.609171052316469 e+02
y2 =1.609171052316469 e+02

>>

Complessità : algoritmo di Horner.

La differenza sta nella complessità computazionale e non nel risultato numerico. Il primo codice richiede 2n moltiplicazioni e n somme, mentre il secondo algoritmo n moltiplicazioni e n somme.