Aplicatii ale descompunerii SVD la compresia imaginii

Factorizarea SVD a unei matrici

Pentru orice matrice $A\in\mathbb{R}^{m\times n}$ de rang $r\leq \min(m,n)$ exista doua matrici ortogonale $U\in\mathbb{R}^{m\times m}$, $V\in\mathbb{R}^{n\times n}$si numerele reale pozitive $\sigma_1\geq \sigma_2\geq\cdots\geq \sigma_r$ astfel incat $A$ se descompune in produsul $A=U\Sigma V^T$, adica:

$$\underbrace{A}_{m\times n}=\underbrace{U}_{m\times m}\underbrace{\left[\begin{array}{ccccccc} \sigma_1&0&\ldots& 0&0&\ldots&0\\ 0&\sigma_2&\ldots&0&0&\ldots&0\\ \vdots&\vdots&\ldots&\vdots&\vdots&\ldots&\vdots\\ 0&0&\ldots&\sigma_r&0&\ldots&0\\ 0&0&\ldots&0&0&\ldots& 0\\ \vdots&\vdots&\ldots&\vdots&\vdots&\ldots&\vdots\\ 0&0&\ldots&0&0&\ldots&0 \end{array}\right]}_{m\times n}\underbrace{V^T}_{n\times n}$$

Aceasta descompunere se numeste descompunerea SVD a matricii $A$, iar $\sigma_1, \ldots, \sigma_r$, valorile singulare ale matricii $A$.

Notand cu $u_1, u_2, \ldots, u_m$ coloanele matricii $U$ si respectiv cu $v_1, v_2, \ldots, v_n$, coloanele matricii $V$, descompunerea SVD se exprima ca o combinatie liniara cu coeficientii $\sigma_i$ a matricilor de rang $1$, $u_i v_i^T$, $i=\overline{1,r}$:

$$A=U\Sigma V^T=\sigma_1 u_1v_1^T+\sigma_2 u_2 v_2^T+\cdots+\sigma_r u_r v_r^T$$

Daca sirul valorilor singulare:

$\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_r>0$.

este rapid descrescator, adica incepand de la un anumit rang $k$, valorile $\sigma_{k+1}, \ldots, \sigma_r$ sunt foarte mici comparativ cu primele valori, atunci primii $k$ termeni din descompunerea SVD permit reconstructia matricii $A$. Mai precis, din primii $k$ termeni construim o aproximatie a matricii $A$: $$A_k=\sigma_1(u_1v_1^T)+\cdots+\sigma_k(u_kv_k^T)$$

Matricea aproximare, $A_k$, factorizeaza astfel:

$$A_k =\underbrace{[u_1|u_2|\ldots|u_k]}_{m\times k} \underbrace{\left[\begin{array}{cccc} \sigma_1&0&\ldots&0\\ 0&\sigma_2&\ldots&0\\ \vdots&\vdots&\ldots&\vdots\\ 0&0&\ldots&\sigma_k\end{array}\right]}_{k\times k}\underbrace{\left[\begin{array}{c} v_1^T\\ v_2^T\\\vdots\\v_k^T\end{array}\right]}_{k\times n}$$

si deci are rangul $k$.

Conform Teoremei Eckart-Young matricea $A_k$ este cea mai buna aproximatie de rang $k$ a matricii $A$, adica dintre toate matricile $B\in\mathbb{R}^{m\times n}$ de rang $k$, distanta de la $A$ la $B$ este minimizata de $A_k$:

$$||A-A_k||_F\leq ||A-B||_F, \quad \forall\:\: B$$

Precizare: Daca $M=(m_{ij})\in\mathbb{R}^{m\times n}$, atunci $||M||_F=\sqrt{\sum_{i,j}m_{ij}^2}$ noteaza norma Frobenius a matricii $M$.

Eroarea aproximarii lui $A$ prin $A_k$ este $Er=||A-A_k||_F=\sqrt{\sigma_{k+1}^2+\cdots+\sigma_r^2}$.

Calculul factorilor din descompunerea SVD in Python

Declarand matricea $A\in\mathbb{R}^{m\times n}$ ca un numpy.array, functia numpy.linalg.svd apelata astfel:

[U,s, Vt]=numpy.linalg.svd(A) sau U, s, Vt=numpy.linalg.svd(A)

calculeaza elementele descompunerii SVD a matricii $A$. U si Vt sunt array-uri 2D, reprezentand matricile $U$ si $V^T$ din formula teoretica $A=U\Sigma V^T$, iar s este un array 1D (un vector) ce contine $\min{(m,n)}$ elemente s=$[\sigma_0,\sigma_1, \ldots, \sigma_{r-1}, 0, \ldots, 0]$, adica valorile singulare in ordine descrescatoare completate cu zerouri, daca este cazul, pana la $\min{(m,n)}$ elemente.

Numarul valorilor singulare este egal cu rangul matricii $A$.

Reconstructia matricii $A$ din factorii SVD, $U$, $\Sigma$, $V^T$

Pentru a genera din vectorul valorilor singulare $s$, matricea $\Sigma$, care are acelasi shape ca matricea $A$, se identifica numarul valorilor singulare (rangul matricii A):

Testul s>1e-7 returneaza un vector de booleeni, avand acelasi nr de elemente ca s:

Cu alte cuvinte, b[i]=True, daca s[i]>1e-7 si b[i]=False, in caz contrar, $i=\overline{0,len(s)-1}$

Matricea $\Sigma$ se constituie din vectorul s astfel:

Linia Sigma[:r,:r]= np.diag(s[:r]) atribuie sublocului matricii Sigma format din primele $r$ linii si $r$ coloane, matricea diagonala generata din vectorul s[:r], adica din vectorul ce contine valorile singulare.

Produsul celor trei matrici din descompunerea SVD, comparat cu matricea $A$:

Compresia imaginilor

O imagine in tonuri de gri este reprezentata de o matrice $A=(a_{ij})$. Elementul $a_{ij}$ al matricii este codul $c\in\{0,1,2,\ldots,255\}$, al nivelului de gri al pixelului din linia $i$, coloana $j$.

O imagine color poate fi reprezentata de 3 matrici $R$, $G$, $B$, de aceeasi rezolutie $m\times n$ si in pozitia $(i,j)$ a fiecarei matrici este stocata respectiv "cantitatea" de red, green, blue, (r, g, b), din care este sintetizata culoarea pixelului din pozitia $(i,j)$ a imaginii. $r, g, b\in\{0,1, \ldots, 255\}$ (pentru a gasi tripletul (r,g,b) ce genereaza o culoare sau pentru a afisa culoarea dintr-un triplet dat se pot folosi instrumentele de aici).

Imaginile in tonuri de gri si cele color, in format png, sunt citite de functia plt.imread in matrici ale caror elemente sunt numere de tip float din intervalul $[0,1]$ (codurile de culoare, $c$, sunt scalate la $c/255.0\in[0,1]$).

Compresia unei imagini digitale consta in reducerea volumului de date necesare pentru a reprezenta imaginea.

Descompunerea SVD a matricii reprezentative a unei imagini in tonuri de gri ofera un instrument de comprimare a imaginii.

Si anume, determinand descompunerea SVD a matricii, $A=\sum_{i=1}^r\sigma_i u_i v_i^T$, si renuntand la termenii ce au coeficientii $\sigma_{k+1}, \ldots, \sigma_r$ suficient de mici in comparatie cu $\sigma_1, \sigma_2,\ldots, \sigma_k$, obtinem o aproximare $A_k$ a imaginii $A$. Cu alte cuvinte imaginea $A_k$ se obtine prin compresia imaginii $A$. Imaginea initiala este constituita din $r$ vectori $u_i$ din $\mathbb{R}^m$, $r$ vectori $v_i$ din $\mathbb{R}^n$ si $r$ valori singulare, in timp ce pentru compresie sunt necesari $k<<r$ vectori $u_i$, respectiv $v_i$, si $k$ valori singulare.

Sa aplicam acum descompunerea SVD a unei imagini in tonuri de gri. Imaginea o selectam din cele trei imagini reprezentative ale unei imagini color, in format png:

• A se observa sistemul stramb de coordonate relativ la care este raportata imaginea. Punctul de coordonate $(0,0)$ este plasat in coltul stanga-sus.

• Imaginea color este reprezentata de un aray 3D, adica 3 array-uri 2D. Interogand imag.shape deducem ca fiecare array 2D are 295 linii si 284 coloane. Cele trei array-uri se acceseaza respectiv, ca imag[:,:,0], imag[:,:,1], imag[:,:,2]. In primul array un element este setat pe cantitatea de red continuta de pixelul respectiv din imaginea color, al doilea array pe cantitatea de green si al treielea pe cantitatea de blue.

Afisam cele trei array-uri ca imagini in tonuri de gri:

Copiem array-ul imag[:,:,0] ce contine cantitatile de red din imaginea initiala, in array-ul A caruia ii aplicam descompunerea SVD:

Sa ilustram nivelul de marime al valorilor singulare consecutive, vizualizand punctele de coordonate (i, sig[i]), $i=\overline{0,r-1}$:

Se observa din grafic o descrestere rapida a valorilor singulare.

Afisand primele 30 valori singulare si respectiv un set de valori consecutive $\sigma_{80}, \ldots, \sigma_{105}$, se observa descresterea brusca de la $\sigma_0=175.36$ la $\sigma_1=35.5$. Incepand cu $\sigma_{105}$ valorile singulare devin subunitare, deci neglijabile.

Rezulta atunci ca orice aproximatie $A_k$ a imaginii $A$, cu $k$ putin mai mare decat $105$ conduce la o imagine ce aproximeaza imaginea initiala foarte bine, din punct de vedere vizual.

Functia urmatoare realizeaza compresia imaginii de rang r printr-una de rang k<<r si o afiseaza:

Comprimam imaginea A, aproximand-o prin imagini de rang $k<<r$, cu $k$ crescand progresiv.

Observam ca pe masura ce se adauga mai multi termeni in descompunerea aproximativa, se amplifica informatia vizuala si imaginea $A_{105}=\sum_{k=0}^{104}\sigma_k u_k v_k^T$ nu se distinge de imaginea originala.

Teoretic, eroarea de aproximare a matricii $A$ prin $A_{105}$ este:

Imaginea $A_{30}$, care este neclara, aproximeaza imaginea $A$ cu o eroare mai mare:

Acest experiment ilustreaza ca daca o matrice imagine are doar cateva valori singulare mari si restul au un nivel de marime cu mult mai mic, atunci valorile singulare neglijabile in raport cu primele sunt coeficienti ai unor termeni din descompunerea SVD, care contin zgomot sau informatie redundanta. Acesti termeni pot fi neglijati fara pierdere de informatie.

Daca imaginea este insa in format jpg atunci ea este citita de plt.imread intr-un array cu elemente de tip uint8, adica din multimea $\{0,1,2,3,\ldots, 255\}$. In acest caz, valorile singulare sunt cu mult mai mari si e dificil sa decidem aproximatia care conduce la eroare redusa.

Observam fie din colorbar, fie afisand tipul de date din array-ul $B$, ca $B$ are elementele in multimea $\{0,1, \ldots, 255\}$.

Din graficul din stanga, care vizualizeaza punctele $(i,\sigma_i)$ se vede ca valorile singulare sunt descrescatoare, dar nu si cat de mult descresc sub 10000.

In dreapta sunt reprezentate punctele $(i, \log_{10}(\sigma_i))$, cu scara logaritmica pe axa y. Din aceasta se observa ca cea mai mica valoare singulara este apropiata de $10^2$.

Afisand subsiruri de valori singulare:

remarcam ca spre deosebire de exemplul precedent, in care de la un anumit rang valorile singulare erau subunitare, in acest caz al unei matrici cu elemente de 255 ori mai mari decat in cazul precedent, este dificil sa decidem care matrice $B_k$ ar fi o aproximatie buna pentru imaginea $B$. Sa afisam de exemplu aproximatia $B_{100}$ si sa calculam eroarea de aproximare:

Imaginea aproximatie este vizual aproape identica cu imaginea $B$. Eroarea de aproximare este:

Aceasta situatie bizara - aproximatie vizual buna, dar eroarea teoretica mare, se datoreaza faptului ca valorile matricii $B$ nu sunt normalizate la $[0,1]$.

De aceea pentru a lua o decizie potrivita privind alegerea ordinului $k$ al aproximarii, care nu conduce la pierdere de informatie relevanta din imagine si nici nu pastreaza informatie redundanta, se poate proceda in doua moduri care intr-un fel sunt echivalente:

1. Matricea $B$ sa fie convertita la dtype=float si elementele ei rescalate la elemente din $[0,1]$, adica $B\mapsto B/255.0$ si apoi sa fie afisata imaginea $C=B/255.0$ si calculata descompunerea ei SVD si o aproximatie $C_k$ adecvata (practic asa procedeaza functiile ce opereaza asupra array-urilor reprezentative ale imaginilor png).

2. Pastram matricea $B$, dar pentru a decide care ordin $k$ conduce la o aproximatia buna, interpretam matricea $B$ ca fiind $B=255.0 C$ si deducem relatia dintre valorile singulare ale unei matrici $C$ si ale scalatei sale $B=aC$, precum si relatia dintre eroarea $||B-B_k||_F$ si $||C-C_k||$.

Din relatia dedusa pentru erori vom observa ca putem decide mai usor care aproximatie $B_k$ este suficient de buna.

Tema: Deduceti relatiile mentionate la 2 si experimentati intr-un Ipython Notebook ambele modalitati de compresie si analiza a erorii pentru o imagine jpg, ce contine poza proprie.