Sistemi lineari sovradeterminati e SVD

27 febbraio 2009

1 Sistemi lineari sovradeterminati

Si consideri il problema

\begin{matrix}_{1} _{2} 1 \\ _{1} _{2} 0 \\ _{1} 3 _{2} 0 \end{matrix}

Risulta chiaro che il sistema non ha soluzione. Infatti l’unica soluzione delle prime due equazioni è _{1} _{2} 1/2, che però non verifica _{1} 3 _{2} 0.

Ritrascrivendo tutto in termini matriciali, possiamo dunque affermare che il problema , con

{ \atop }\begin{matrix}1 & 1 \\ 1 & 1 \\ 1 & 3 \end{matrix} { \atop }

e (1 0 0) non ha soluzione.

In alternativa, risulta quindi ragionevole calcolare ^{}( _{1}^{}_{2}^{}) tale che

: min_{^{}}_{2}  ^{}_{2}  (1)

dove

2 \overline{_{}_{}^{2}}

La prima questione riguarda l’esistenza e unicità di tale minimo ^{}. A tal proposito citiamo il segue teorema (cf. [1], p. 432)

Teorema 1.1 Sia l’insieme dei vettori di ^{}tali che ˆ se e solo se

min_{^{}}_{2}  _{  ˆ 2}  (2)
Supponiamo  ^{}con  e  ^{}. Valgono le seguenti proprietà

1. se e solo se

(3)

cioé risolve il sistema delle equazioni normali.

2. è un insieme non vuoto, chiuso e convesso.

3. Esiste ^{}tale che

2 min_{}_{2}

Tale ^{}è detto soluzione di minima norma.

4. L’insieme si riduce ad un solo elemento ^{}se e solo se la matrice ha rango massimo.

In altre parole, se ha rango allora ha un unico elemento, mentre se ha rango minore di allora ha un unico elemento di minima norma 2.

2 Alcune fattorizzazioni di matrici rettangolari.

2.1 Fattorizzazione QR

Data una matrice ^{}esistono ^{}unitaria (cioė ^{}_{}, ^{}_{}) ed ^{}triangolare superiore tali che . Si osservi che a seconda degli autori la fattorizzazione QR sopraindicata viene sostituita con la fattorizzazione

con  ^{}unitaria ed  ^{}triangolare superiore cioè

\begin{matrix}_{1} \\ 0 \end{matrix}

con _{1} ^{}triangolare superiore.

2.2 Fattorizzazione SVD

Sussiste il seguente teorema [7]:

Teorema 2.1 Sia  ^{}. Allora esistono due matrici unitarie  ^{},

e una matrice diagonale _{}^{}avente elementi _{}nulli per e uguali a _{}_{}per con

1 _{2}  _{}0  (4)

tali che

Il termine sta per singular value decomposition e sottolinea la presenza dei valori singolari _{}. Se _{}e _{}sono rispettivamente gli autovalori di ^{}in ordine decrescente e i valori singolari di in ordine decrescente, allora

\overline{_{}}1

3 Fattorizzazione QR ed SVD in Matlab

Vediamo di seguito alcuni dettagli di come tali fattorizzazioni sono implementate in Matlab.

3.1 QR in Matlab

In merito alla fattorizzazione QR, l’help di Matlab fornisce le seguenti indicazioni:

>> help qr
QR     Orthogonal-triangular decomposition.
[Q,R] = QR(A) produces an upper triangular matrix R of the same
dimension as A and a unitary matrix Q so that A = Q*R.
[Q,R,E] = QR(A) produces a permutation matrix E, an upper
triangular R and a unitary Q so that A*E = Q*R.  The column
permutation E is chosen so that abs(diag(R)) is decreasing.
[Q,R] = QR(A,0) produces the "economy size" decomposition.
If A is m-by-n with m > n, then only the first n columns of Q
are computed.
[Q,R,E] = QR(A,0) produces an "economy size" decomposition in
which E is a permutation vector, so that Q*R = A(:,E).  The column
permutation E is chosen so that abs(diag(R)) is decreasing.
By itself, QR(A) is the output of LAPACK'S DGEQRF or ZGEQRF routine.
TRIU(QR(A)) is R.
For sparse matrices, QR can compute a "Q-less QR decomposition",
which has the following slightly different behavior.
R = QR(A) returns only R.  Note that R = chol(A'*A).
[Q,R] = QR(A) returns both Q and R, but Q is often nearly full.
[C,R] = QR(A,B), where B has as many rows as A, returns C = Q'*B.
R = QR(A,0) and [C,R] = QR(A,B,0) produce economy size results.
The sparse version of QR does not do column permutations.
The full version of QR does not return C.
The least squares approximate solution to A*x = b can be found
with the Q-less QR decomposition and one step of iterative refinement:
    x = R\(R'\(A'*b))
    r = b - A*x
    e = R\(R'\(A'*r))
    x = x + e;
See also LU, NULL, ORTH, QRDELETE, QRINSERT, QRUPDATE.
>>

Nel caso della matrice

{ \atop }\begin{matrix}1 & 1 \\ 1 & 1 \\ 1 & 3 \end{matrix} { \atop }

vediamo alcune fattorizzazioni QR.

>> A=[1 1; 1 -1; 1 3]
A =
     1     1
     1    -1
     1     3
>> % PRIMA FATTORIZZAZIONE QR.
>> [Q,R]=qr(A)
Q =
   -0.5774    0.0000   -0.8165
   -0.5774   -0.7071    0.4082
   -0.5774    0.7071    0.4082
R =
   -1.7321   -1.7321
         0    2.8284
         0         0
>> Q'*Q
ans =
    1.0000         0   -0.0000
         0    1.0000   -0.0000
   -0.0000   -0.0000    1.0000
>> norm(A-Q*R)
ans =
  4.9938e-016
>> % SECONDA FATTORIZZAZIONE QR.
>> [Q,R]=qr(A,0)
Q =
   -0.5774    0.0000
   -0.5774   -0.7071
   -0.5774    0.7071
R =
   -1.7321   -1.7321
         0    2.8284
>> Q'*Q
ans =
    1.0000         0
         0    1.0000
>> norm(A-Q*R)
ans =
  4.9938e-016
>> % TERZA FATTORIZZAZIONE QR.

>> [Q,R,E]=qr(A)
Q =
   -0.3015   -0.4924   -0.8165
    0.3015   -0.8616    0.4082
   -0.9045   -0.1231    0.4082
R =
   -3.3166   -0.9045
         0   -1.4771
         0         0
E =
     0     1
     1     0
>> Q'*Q
ans =
    1.0000         0    0.0000
         0    1.0000    0.0000
    0.0000    0.0000    1.0000
>> norm(A*E-Q*R)
ans =
  5.6322e-016
>> A*E-Q*R
ans =
  1.0e-015 *
    0.2220    0.1110
    0.2220    0.2220
   -0.4441         0
>>

Nota facoltativa. Dal punto di vista implementativo è facile calcolare la fattorizzazione di anche quando il rango non è massimo, utilizzando [Q,R,E]=qr(A). Si stabilisce facilmente il rango di cercando il numero di elementi diagonali di non nulli.

3.2 SVD in Matlab

Osserviamo che Matlab dispone del comando svd. La descrizione dell’help è la seguente:

>> help svd
 SVD    Singular value decomposition.
    [U,S,V] = SVD(X) produces a diagonal matrix S, of the same
    dimension as X and with nonnegative diagonal elements in

decreasing order, and unitary matrices U and V so that
X = U*S*V'.
S = SVD(X) returns a vector containing the singular values.
[U,S,V] = SVD(X,0) produces the "economy size"
decomposition. If X is m-by-n with m > n, then only the
first n columns of U are computed and S is n-by-n.
See also SVDS, GSVD.
Overloaded methods
help sym/svd.m
>>

Per fare pratica con questo comando sia

{ \atop }\begin{matrix}1 & 1 \\ 1 & 1 \\ 1 & 3 \end{matrix} { \atop }

e utilizziamo [U,S,V] = svd(A) per avere la decomposizione SVD della matrice . Il risultato è

>> A=[1 1; 1 -1; 1 3]
A =
 1     1
 1    -1
 1     3
>> [U,S,V] = svd(A)
U =
0.3651    0.4472   -0.8165
-0.1826    0.8944    0.4082
0.9129   -0.0000    0.4082
S =
3.4641         0
     0    1.4142
     0         0

V =
    0.3162    0.9487
    0.9487   -0.3162
>> % U, V SONO UNITARIE.
>> U'*U
ans =
    1.0000   -0.0000   -0.0000
   -0.0000    1.0000    0.0000
   -0.0000    0.0000    1.0000
>> V'*V
ans =
     1     0
     0     1
>>
>> X = U*S*V'
X =
    1.0000    1.0000
    1.0000   -1.0000
    1.0000    3.0000
>>

Si puó dimostrare che il rango della matrice è esattamente il numero di _{}non nulli. In realtá la tentazione di determinarlo in virtú dei termini diagonali di puó risultare fallace in quanto il computer fornisce solo un’approssimazione dei _{}e qualora questi siano molto piccoli non si puó determinare con sicurezza il rango della matrice. Nota Facoltativa. Esiste un legame tra valori singolari e autovalori. Più precisamente se _{}sono gli autovalori di ^{}allora _{}\overline{_{}}. Vediamolo con un esempio:

>> A=[1 1; 1 -1; 1 3]
A =
     1     1
     1    -1
     1     3
>> [U,S,V]=svd(A);
>> S
S =
    3.4641         0
         0    1.4142
         0         0
>> % LEGAME VALORI SINGOLARI E AUTOVETTORI.

>> lambda=eig(A'*A)
lambda =
     2
    12
>> sqrt(12)
ans =
    3.4641
>> sqrt(2)
ans =
    1.4142
>>

4 Calcolo della soluzione ai minimi quadrati

Mostriamo ora tre metodi per risolvere il problema sopra descritto, cioè relativo al calcolo di ^{}che minimizza la norma 2 del residuo, cioè _{2}.

4.1 Risoluzione equazioni normali

Supponiamo �^{}abbia rango . Dal teorema 1.1sappiamo che la soluzione ai minimi quadrati esiste, è unica e risolve ^{}^{}. Nel caso abbia coefficienti in , ció si riduce a risolvere ^{}^{}. In questo caso ^{}è definita positiva e si pò utilizzare il metodo di Cholesky. Se ^{}^{}per triangolare inferiore (fattorizzazione di Cholesky) basta risolvere

Il costo computazionale è di ^{2} /2 operazioni moltiplicative per la costruzione di ^{}e di ^{3}/6 moltiplicative per la soluzione del sistema, e quindi il costo totale è di

2 /2 ^{3}/6

operazioni moltiplicative. Osserviamo che se non ha rango massimo non si puó applicare il metodo di Cholesky ma il metodo di Gauss con variante di massimo pivot. Nota facoltativa. Se non ha rango massimo non si può applicare il metodo di Cholesky, ma l’eliminazione gaussiana con la variante del massimo pivot.

4.2 Metodo QR

Supponiamo il rango di �^{}sia . Nota la fattorizzazione si deduce che .

Osserviamo infatti che essendo ^{}_{} e quindi prima calcoliamo ^{}e poi risolviamo ^{}. Il passaggio

è delicato, in quanto ^{}è unitaria (cioė ^{}_{}, ^{}_{}) ma rettan-

golare. Dice che posto  ^{}, l’equazione diventa  . Infatti, le corrispettive
equazioni normali (si ricordi che  ^{}con  è unitaria) sono

Se ha rango massimo allora è non singolare e quindi il problema ^{} risolto facilmente per sostituzione all’indietro. Il costo computazionale per la fattorizzazione è di ^{2}( /3), il costo computazionale della risoluzione del sistema triangolare è ^{2}/2 operazioni moltiplicative. Quindi il costo totale è

2( /3) ^{2}/2

Nota facoltativa. Se non ha rango massimo allora con matrice di permutazione,

\begin{matrix}_{1} & \\ 0 & 0 \end{matrix}

dove ^{}, _{1} ^{}sono matrici triangolari superiori. Si dimostra (non facile) che posto

\begin{matrix}_{1} \\ _{2} \end{matrix} _{1} ^{}_{2} ^{}

la soluzione ai minimi quadrati di minima norma risolve il problema

1 _{1}

4.3 Metodo SVD

Siano _{}, _{}le i-sime righe rispettivamente di e . Per capire come risolvere un sistema lineare sovradeterminato tramite questa fattorizzazione ci affidiamo al seguente teorema

Teorema 4.1 Sia  ^{}di rango con  e sia

la decomposizione ai valori singolari di . Allora la soluzione di minima norma è data da

_{1}^{}\frac{_{}^{}}{_{}}_{}

e

2 _{1}^{}_{}^{}^{2}

dove al solito

: min_{^{}}_{2}  ^{}_{2}  (5)

5 Un esempio: soluzione di sistemi sovradeterminati.

Proviamo a risolvere il problema iniziale

\begin{matrix}_{1} _{2} 1 \\ _{1} _{2} 0 \\ _{1} 3 _{2} 0 \end{matrix}

coi metodi sopracitati.

Lavoriamo direttamente nella shell di Matlab/Octave:

>> % EQUAZIONI NORMALI.
>> A=[1 1; 1 -1; 1 3];
>> b=[1 0 0]';
>> C=A'*A;
>> b1=A'*b;
>> format long
>> x=C\b1
x =
   0.33333333333333
  -0.00000000000000
>> % METODO QR
>> [Q,R]=qr(A);
>> size(Q)
ans =
     3     3
>> size(R)
ans =
     3     2
>> % COMMENTO: NON E'LA FATTORIZZAZIONE CERCATA
>> % "R" DEVE ESSERE QUADRATA.
>> [Q,R]=qr(A,0);
>> size(Q)
ans =
     3     2
>> size(R)
ans =
     2     2
>> b1=Q'*b;
>> size(b1)
ans =
     2     1
>> x=R\b1
x =
   0.33333333333333
   0.00000000000000
>> % SOLUZIONE VIA SVD.
>> [U,S,V]=svd(A);
>> size(U)
ans =

3     3
>> size(V)
ans =
2     2
>> size(b)
ans =
3     1
>> c1=U'*b;
>> size(c1)
ans =
3     1
>> c2=c1(1:size(V,1));
>> s=c2./diag(S);
>> x=V*s
x =
0.33333333333333
0.00000000000000
>> gamma_value=norm(A*x-b,2)
gamma_value =
0.81649658092773
>> abs((U(:,3))'*b)
ans =
0.81649658092773
>>

Commentiamo i risultati a partire da quanto visto per il metodo QR, visto che la risoluzione via equazioni normali è del tutto ovvia.

1. L’implementazione in Matlab di QR per matrici rettangolari richiede una precisazione. Sia ^{}la soluzione del problema sovradeterminato. Il comando [Q,R]=qr(A) comporta che sia quadrata, mentre [Q,R]=qr(A,0) implica che è quadrata, come da noi richiesto.

2. Per quanto riguarda la soluzione via SVD, l’unico punto da osservare è che in un precedente teorema si dice che

_{1}^{}\frac{_{}^{}}{_{}}_{}

dove _{}, _{}sono le righe di e . La sommatoria è limitata alle prime righe e per questo abbiamo effettuato

>> c1=U'*b;
>> c2=c1(1:size(V,1));
>> s=c2./diag(S);

Si osservi che il rango per la matrice in questione è uguale a cioè 2, che altri non è che ( 1). Se il rango non fosse stato massimo, cioè , prima sarebbe stato opportuno calcolare il rango e poi sostituire c1(1:size(V,1),:);

con c1(1:r,:);

Ricordiamo che 1:r genera il vettore riga [12 ] e che c2=c1(1:r); assegna a c2 le righe dalla prima all’ -sima del vettore c1

>> r=5;
>> 1:r
ans =
     1     2     3     4     5
>> c1=[1 4 2 6 4 8 4]
c1 =
     1     4     2     6     4     8     4
>> c1(1:r)
ans =
     1     4     2     6     4

Essendo _{}le righe di , posto _{}_{}^{}/ _{}risulta che

_{1}^{}\frac{_{}^{}}{_{}}_{}

Nell’ultima parte si mostra che effettivamente la norma 2 del residuo cioè coincide con

2 _{1}^{}_{}^{}^{2}

In merito osserviamo che il rango della matrice è 2 ed 3 e quindi

2 _{1}^{}_{}^{}^{2} _{}^{}^{2}

da cui

2  _{3}^{}^{2}  _{3}^{}

3. Può restare il dubbio di come si faccia a sapere per via numerica che il rango di è 2. Questo corrisponde al numero di coefficienti diagonali di non nulli e quindi da

>> A=[1 1; 1 -1; 1 3];
>> [U,S,V]=svd(A);
>> diag(S)
ans =
    3.4641
    1.4142
>> t=find(abs(diag(S))>0);
>> rango=length(t)
rango =
     2
>>

deduciamo che il rango è 2.

6 Un esempio: compressione immagini.

Consideriamo i siti web [6], [2], [4]. L’argomento esposto è la compressione di immagini via SVD. In particolare si accenna ad una interessante implementazione in Matlab (non funziona in Octave!) relativa alla compressione di immagini. Entriamo nei dettagli. Sia ^{}la fattorizzazione SVD della matrice . Se _{}sono gli elementi diagonali della matrice diagonale e rettangolare , _{}la -sima colonna di , _{}la -sima colonna di , allora

_{1}^{}_{}_{}_{}^{}

Infatti, se _{}�^{}è il vettore colonna avente tutte componenti nulle eccetto la sima che è uguale a 1, indicato con _{}^{}�^{}la -sima riga di ^{}e con _{}la -sima colonna di abbiamo dopo qualche conto (provare a farli, anche a modo proprio!) che

_{1}^{}_{}_{}_{}^{}

da cui evidenziando il prodotto tra matrici , ricordando la sua proprietà associativa e ricordando la linearità delle operazioni

( ^{})
_{1}^{}_{}_{}_{}^{}
_{1}^{}_{}_{}_{}^{}
_{1}^{}_{}( _{}) _{}^{}
_{1}^{}_{}_{}_{}^{}

(6)

per cui essendo _{}la -sima colonna di , _{}la -sima colonna di abbiamo _{} e _{}^{}_{}^{}, nonchè

_{1}^{}_{}_{}_{}^{}

Osserviamo che le matrici _{}: _{}_{}^{}hanno rango 1 in quanto la loro immagine è uno spazio di dimensione 1, visto che

_{}^{}( ) _{}

dove al solito ( ) è il prodotto scalare in �^{}.

Se l’indice viene ordinato cosicchè i valori singolari siano in ordine decrescente, cioè

1 _{2} _{}

e tronchiamo la somma dopo termini otteniamo una approssimazione di rango della matrice . Se in particolare la matrice è una matrice ottenuta codificando i bit di un’immagine, si capisce che si può comprimere la foto utilizzando per invece di la matrice

_{1}^{}_{}_{}_{}^{}

Vediamone un implementazione in Matlab. Scarichiamo dal sito [4] il file compression.zip e una volta decompresso, utilizziamo la routine imcompr.m per comprimere un’immagine di tipo bitmap, come ad esempio il gatto salvato in gatto.bmp. Osserviamo che il gatto è codificato tramite una matrice di interi di dimensione 416 325.

Lanciando dalla shell di Matlab (non funziona in Octave!) il comando

>> [im] = imcompr ('gatto.bmp',20,'gatto20.bmp');
>> [im] = imcompr ('gatto.bmp',50,'gatto50.bmp');
>> [im] = imcompr ('gatto.bmp',100,'gatto100.bmp');

Il risultato sono 3 foto dello stesso gatto, con diversi fattori di compressione, in quanto abbiamo usato valori di diversi (r=20, 50, 100).

Per ulteriori dettagli si confronti [3], p.180, esempio 6.9. Nota. Siò notare che dopo la compressione, la dimensione del file compresso è maggiore di quello originale. La ragione è che il secondo è relativo ad una matrice avente componenti intere mentre quelle del primo sono reali. Nonostante questo problema, l’idea mostra una possibile strada per comprimere immagini. Inoltre non è chiaro perchè la compressione funzioni in pratica e quindi .

7 Un esempio: approssimazione polinomiale.

Sia una funzione reale e continua e siano _{}punti a due a due distinti appartenenti al dominio di . Si cerca il polinomio

1 ( ) _{0} _{1} _{1} ^{1} _{}

per cui sia minimo lo scarto quadratico

_{1}^{}[ _{1} ( _{}) ( _{})]^{2}

Tale polinomio _{1} esiste unico ed è noto come polinomio di miglior approssimazione . Si prova che ció corrisponde a risolvere il problema sovradeterminato ¯ dove

; ; ; \begin{matrix}1 & _{1} & & _{1}^{1} \\ 1 & _{2} & & _{2}^{1} \\ & & & \\ 1 & _{} & & _{}^{1} \end{matrix} ; ; ;

Sistemi lineari sovradeterminati e SVD

Figura 1: Compressione di una foto via SVD. e ¯ ( ( _{1}) ( _{})), ^{}, , ( _{}) ^{}. Questa idea è alla base dell’approssimazione fornita dalla toolbox polyfit. Vediamo i dettagli. Notiamo subito che

1 ( _{}) _{0} _{1} _{}_{1} _{}^{1} _{}

e quindi per definizione di ¯

_{1}^{}[ _{1} ( _{}) ( _{})]^{2} _{1}^{}[ _{}_{}¯ ]^{2}

Ricordiamo ora che la soluzione di un sistema sovradeterminato (con ) è quella che minimizza _{2} e quindi pure

_{2}^{2} _{}[( )_{}_{}]^{2}

Di conseguenza il vettore avente quali coefficienti i termini _{}del polinomio di miglior approssimazione del problema di approssimazione è il vettore che risolve il

sistema sovradeterminato  ¯ .

References

[1] D. Bini, M. Capovani, O. Menchi, Metodi numerici per l’algebra lineare, Zanichelli, 1988.

[2] Cleve’s corner, Professor SVD http://www.mathworks.com/company/newsletters/news_notes/oct06/clevescorner.html .

[3] A. Quarteroni ed F. Saleri, Introduzione al Calcolo Scientifico, Springer, 2006.

[4] L. Rosa, Image compression http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4772&objectType=File .

[5] K. Sandberg, The Singular Value Decomposition of an Image http://amath.colorado.edu/courses/4720/2000Spr/Labs/SVD/svd.html

[6] G.W. Steward, On the Early History of the Singular Value Decomposition http://www.ima.umn.edu/preprints/April92/952.pdf .

[7] Wikipedia (Singular value decomposition)

http://en.wikipedia.org/wiki/Singular_value_decomposition .