Ctrl+K

SVD

Hide code cell source 
import numpy as np;
import scipy as sc;
import matplotlib.pyplot as plt;
from scipy.optimize import curve_fit;

%matplotlib inline

4.3. SVD#

Pokažimo još jedan način kako riješiti diskretan problem najmanjih kvadrata. Alat za to je SVD (singularna dekompozicija ili dekompozicija singularnih vrijednosti, eng. singular value decomposition). Na kraju poglavlja pokazat ćemo da problem najmanjih kvadrata nije jedino mjesto gdje se može primijeniti SVD.

Definicija 4.7

Neka je \(A \in \R^{n\times m}\), uz \(n\geq m\), matrica ranga \(r\leq m\). Njezina singularna dekompozicija je rastav oblika

\[A = U \Sigma V^T,\]

pri čemu su \(U \in \R^{n \times n}\) i \(V \in \R^{m \times m}\) ortogonalne, a \(\Sigma \in \R^{n\times m}\) matrica koja zadovoljava \(\Sigma_{i,j} = 0\) za sve \(i\not= j\), te vrijednosti \(\sigma_i := \Sigma_{i,i}\) zadovoljavaju

\[\sigma_1 \geq \dots \geq \sigma_r > 0 = \sigma_{r+1} =\dots = \sigma_{m}.\]

Skraćeni zapis SVD je

\[A = U_r \Sigma_r V_r^T,\]

pri čemu su \(U_r\) i \(V_r\) matrice dobivene od prvih \(r\) stupaca matrica \(U\) i \(V\) redom, a \(\Sigma_r = \Sigma (1:r,1:r)\).

SVD se analogno definira za pravokutne matrice kod kojih je \(n \leq m\).

Vrijednosti \(\sigma_i\) nazivaju se singularnim vrijednostima (eng. singular values). Pokazuje se da su singularne vrijednosti za matricu jedinstvene, iako same matrice \(U\) i \(V\) ne moraju biti (slično kao kod svojstvenih vrijednosti za hermitske matrice). Stupci matrica \(U\) i \(V\) redom se nazivaju lijevi i desni singularni vektori.

Označimo s \(U_r^\perp\) i \(V_r^\perp\) matrice dobivene od preostalih stupaca matrica \(U\) i \(V\). Ako imamo potpuni SVD matrice \(A\), onda vrijedi

\[\begin{split}A=\begin{bmatrix} U_{r} & U_{r}^{\perp} \end{bmatrix} \begin{bmatrix} \Sigma_{r} & 0\\ 0 & 0 \end{bmatrix} \begin{bmatrix} V_{r}^{T}\\ (V_{r}^{\perp})^{T} \end{bmatrix} = U_{r}\Sigma_{r}V_{r}^{T},\end{split}\]

što opravdava skraćeni zapis.

Teorem 4.8

Svaka matrica \(A \in \R^{n\times m}\) dopušta SVD.

Dokaz.

Bez smanjenja općenitosti pretpostavimo da je \(n \geq m\) (u protivnom gledamo \(A^T\)). Tvrdnju dokazujemo indukcijom po broju stupaca \(m\).

Za \(m=1\), neka je \(a\) jedini stupac matrice \(A \in \R^{n \times 1}\). Označimo \(\sigma_1 = \|a\|_2\), \(u_1 = a/ \|a\|_2 \in \R^n\). Tada je

\[\begin{split} A = \underbrace{\mb{cccc} u_1 & u_2 & \ldots & u_n \me}_{U} \cdot \underbrace{\mb{c} \sigma_1 \\ 0 \\ \vdots \\ 0 \me}_{\Sigma} \cdot \underbrace{[1]}_{V^T} \end{split}\]

SVD matrice \(A\), pri čemu su \(u_2, \ldots, u_n\) proizvoljni vektori koji nadopunjavaju \(\{u_1\}\) do ortonormirane baze od \(\R^n\). U matrici \(\Sigma\) ima \(n-1\) nula.

Pretpostavimo sada da za svaku matricu s \(m-1\) stupaca postoji SVD. Neka je \(A \in \R^{n \times m}\) matrica s \(m\) stupaca. Označimo \(\sigma_1 = \|A\|_2\) te neka je \(v_1 \in \R^m\) vektor takav da vrijedi \(\|Av_1\|_2 = \sigma_1\). Stavimo i \(u_1 = Av_1 / \|Av_1\|_2 \in \R^n\), stoga vrijedi \(Av_1 = \sigma_1 u_1\).

Nadopunimo matrice \([u_1]\) i \([v_1]\) na proizvoljan način do ortogonalnih: \(\tilde{U} = [u_1, X] \in \R^{n \times n}\), \(\tilde{V} = [v_1, Y] \in \R^{m \times m}\). Promotrimo matricu

\[\begin{split}\begin{align*} \tilde{A} := \tilde{U}^T \cdot A \cdot \tilde{V} &= \mb{c} u_1^T \\ X^T \me \cdot A \cdot \mb{cc} v_1 & Y \me \\ &= \mb{cc} u_1^T A v_1 & u_1^T A Y \\ X^T A v_1 & X^T A Y \me \\ &= \mb{cc} \sigma_1 & u_1^T A Y \\ \sigma_1 X^T u_1 & X^T A Y \me \\ &= \mb{cc} \sigma_1 & w^T \\ 0 & A_1 \me . \end{align*}\end{split}\]

Ovdje smo u zadnjem koraku iskoristili ortogonalnost matrice \(\tilde{U}\), te označili \(w^T = u_1^T A Y \in \R^{1 \times (m-1)}\) i \(A_1 = X^T A Y \in \R^{(n-1) \times (m-1)}\). Kako su matrice \(\tilde{U}\) i \(\tilde{V}\) ortogonalne, vrijedi \(\|\tilde{A}\|_2 = \|A\|_2\). Međutim,

\[\begin{split} \left\| \tilde{A} \mb{c} \sigma_1 \\ w \me \right\|_2 = \left\| \mb{cc} \sigma_1 & w^T \\ 0 & A_1 \me \mb{c} \sigma_1 \\ w \me \right\|_2 = \left\| \mb{c} \sigma_1^2 + \|w\|^2 \\ A_1 w \me \right\|_2 \geq \sigma_1^2 + \|w\|_2^2 = \sqrt{\sigma_1^2 + \|w\|_2^2} \cdot \left\| \mb{c} \sigma_1 \\ w \me \right\|_2\end{split}\]

pa je \(\|\tilde{A}\|_2 \geq \sqrt{\sigma_1^2 + \|w\|^2}\). Zbog \(\|\tilde{A}\|_2 = \sigma_1\) to je moguće samo u slučaju \(\|w\|_2 = 0\), tj. \(w = 0\). Stoga imamo

\[\begin{split} \tilde{A} = \mb{cc} \sigma_1 & 0 \\ 0 & A_1 \me .\end{split}\]

Kako je \(A_1 \in \R^{(n-1) \times (m-1)}\) matrica s \(m-1\) stupaca, po pretpostavci indukcije za nju postoji SVD: \(A_1 = U_1 \Sigma_1 V_1^T\), pa slijedi

\[\begin{split} \tilde{A} = \tilde{U}^T A \tilde{V} = \mb{cc} \sigma_1 & 0 \\ 0 & U_1 \Sigma_1 V_1^T \me = \mb{cc} 1 & 0 \\ 0 & U_1 \me \cdot \mb{cc} \sigma_1 & 0 \\ 0 & \Sigma_1 \me \cdot \mb{cc} 1 & 0 \\ 0 & V_1 \me^T .\end{split}\]

Konačno, množenjem s \(\tilde{U}\) slijeva i s \(\tilde{V}^T\) zdesna imamo:

\[\begin{split} A = \underbrace{\left(\tilde{U} \mb{cc} 1 & 0 \\ 0 & U_1 \me \right)}_{U} \cdot \underbrace{ \mb{cc} \sigma_1 & 0 \\ 0 & \Sigma_1 \me }_{\Sigma} \cdot \underbrace{ \left(\tilde{V} \mb{cc} 1 & 0 \\ 0 & V_1 \me \right)^T}_{V^T},\end{split}\]

što je SVD matrice \(A\) jer je umnožak ortogonalnih matrica ponovno ortogonalna matrica, a matrice \(\mb{cc} 1 & 0 \\ 0 & U_1 \me\) i \(\mb{cc} 1 & 0 \\ 0 & V_1 \me\) su ortogonalne. Uočimo još da je \(\|A_1\|_2 \leq \|A\|_2 = \sigma_1\) (dokažite to!), pa će, „odmotamo” li induktivni dokaz, na dijagonali od \(\Sigma\) brojevi biti monotono opadajući.

Rješavanje diskretnog problema najmanjih kvadrata SVD-om#

Riješimo sada problem najmanjih kvadrata koristeći SVD. Neka su \(A\) i \(b\) matrica i vektor iz problema najmanjih kvadrata i neka je

\[A = U \Sigma V^T = U_r \Sigma_r V_r^T,\]

SVD matrice \(A\). Ponovno koristeći unitarnu invarijantnost, imamo

\[\begin{split}\| Ax-b\|^2_2 &= \left\| U^T (Ax-b) \right\|^2_2 \\ &= \left\| U^T U \Sigma V^T x- U^T b \right\|^2_2 \\ &= \left\| \begin{bmatrix} \Sigma_r V_r^T \\ 0\end{bmatrix} x- \begin{bmatrix} U_r^T \\ (U_r^\perp)^T \end{bmatrix} b \right\|^2_2 \\ &= \left\| \begin{bmatrix} \Sigma_r V_r^Tx - U_r^Tb \\ 0 - (U_r^\perp)^T b\end{bmatrix} \right\|^2_2 \\ &= \left\| \Sigma_r V_r^Tx - U_r^Tb \right\|_2^2 + \left\| (U_r^\perp)^T b \right\|^2_2.\end{split}\]

Slično kao za u slučaju rješavanja QR faktorizacijom, zaključujemo da je

\[\min_{x \in \R^m} \| Ax-b\|_2 = \left\| (U_r^\perp)^T b \right\|_2,\]

te se postiže za \(x\) koji je rješenje sustava \(\Sigma_r V_r^Tx = U_r^Tb\). Taj je sustav ponovno jednostavno za riješiti, budući da se regularna dijagonalna matrica \(\Sigma_r\) i ortogonalna matrica \(V_r^T\) lagano invertiraju.

Pokažimo i kod u Pythonu za rješenje Primjera 4.1 korištenjem singularne dekompozicije.

Hide code cell source 
def SVDLinearniNajmanjiKvadrati (x,y):
    # Izrada sustava jednadžbi
    n = len(x)
    A = np.zeros((n, 2))
    A[:, 0] =x[:]
    A[:, 1] = np.ones(n)
    B = y

    # Rješavanje sustava
    U,S,Vt = np.linalg.svd(A, full_matrices=False)
    # koristimo skraćeni zapis, U i Vt imaju nisu kvadratne nego imaju 2 stupca/retka
    # Vt je transponirana matrica V iz definicije SVD-a
    # S je dijagonala dijagonalne matrice 2x2
    rhs = U.T @ B
    rhs[:] = rhs[:] / S[:]
    # djelovanje s S^(-1)
    rj = Vt.T @ rhs
    a = rj[0]
    b = rj[1]
    return a, b

# Zadani podaci
x = np.array([1, 2, 3, 4, 5])
y = np.array([2.5, 3.7, 3.5, 4.5, 4.9])

a, b = SVDLinearniNajmanjiKvadrati(x, y)

# Ispis rezultata
print("Vrijednost a:", a)
print("Vrijednost b:", b)

Vrijednost a: 0.5600000000000002
Vrijednost b: 2.139999999999999
Hide code cell source 
# Dio koda za vizualizaciju
def model(x, a, b):
  return a * x + b
x_line = np.linspace(min(x), max(x), 100)
y_line = model(x_line, a, b)
plt.plot(x, y, 'o', label='Podaci')
plt.plot(x_line, y_line, label='Pravac')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Skup podataka i optimalni pravac')
plt.legend()
plt.grid(True)
plt.show()
_images/0f34dd6ac72fe4d57fd3e9be6f711f73452693d1567b4628aa1c0e6f833f3738.png

Dodatna svojstva SVD-a#

Singularna dekompozicija ima puno veću primjenu nego u rješavanju problema najmanjih kvadrata. Pokažimo prvo teorijska svojstva SVD-a, a kasnije ćemo pokazati i neke primjere.

Sljedeća propozicija slijedi direktno iz definicije SVD.

Propozicija 4.9

Dana je matrica \(A \in \R^{n\times m}\). Tada su netrivijalne singularne vrijednosti matrice \(A\) upravo drugi korijeni netrivijalnih svojstvenih vrijednosti matrica \(AA^T\) i \(A^TA\).

Iz singularne dekompozicije matrice jednostavno čitamo jezgru i sliku matrice, te unitarno invarijantne norme matrice \(A\).

Propozicija 4.10

Neka je \(A = U \Sigma V^T\) singularna dekompozicija matrice \(A\in \R^{n\times m}\). Neka je točno \(r\) njenih singularnih vrijednosti različito od nula. Neka su \(v_i\) i \(u_i\) oznake za \(i\)-te stupce matrica \(V\) i \(U\). Tada vrijedi:

  • \(r(A) = r\);

  • \(\mathop{\mathrm{Ker}}A = \text{span}\{ v_{r+1}, \ldots, v_{m} \}\);

  • \(\mathop{\mathrm{Im}} A = \text{span}\{ u_{1}, \ldots, u_{r} \}\);

  • \(\|A\|_2 = \sigma_1\);

  • \(\| A \|_F^2 = \sigma_1^2 + \ldots + \sigma_r^2\).

Dokaz.

Matrica \(\Sigma_r\) očito ima rang \(r\), dok množenjem te matrice regularnim matricama ne mijenja njezin rang. Zato je \(r(A) = r\).

Pomnožimo jednakost \(A = U \Sigma V^T\) zdesna matricom \(V\). Dobivamo

\[A V = U \Sigma.\]

Odavde čitamo da za svaki indeks \(1 \leq i \leq r\) vrijedi \(A v_i = \sigma_i u_i \not= 0\), dok za sve indekse \(i > r\) vrijedi \(A v_i = 0\). Sada tvrdnje druga i treća tvrdnja slijede iz činjenice da je \(\{ v_1, v_2, \ldots, v_n \}\) baza za \(\mathbb R^n\), a \(\{ u_1, u_2, \ldots, u_m \}\) baza za \(\mathbb R^m\).

2-norma matrice je unitarno invarijantna: množenjem matrice slijeva ili zdesna unitarnom matricom ne mijenja normu te matrice. Zato vrijedi

\[\| A \|_2 = \| U \Sigma V^T \|_2 = \| \Sigma \|_2 = \sigma_1;\]

uvjerite se iz definicije 2-norme da vrijedi zadnja jednakost. Uočimo i da ova tvrdnja zapravo direktno slijedi i iz našeg dokaza egzistencije SVD-a.

Frobeniusova norma matrice je također unitarno invarijantna pa imamo

\[\| A \|_F^2 = \| U \Sigma V^T \|_F^2 = \| \Sigma \|_F^2 = \sigma_1^2 + \ldots + \sigma_r^2,\]

čime je dokazana i zadnja tvrdnja.

Na kraju dajemo još rezultat koji govori o aproksimacijama matrice nižeg ranga. Kako ćemo vidjeti kasnije, ovo svojstvo SVD-a je jedan od glavnih razloga za njezinu široku upotrebu u raznim primjenama.

Teorem 4.11 (Eckart, Young, Mirsky)

Neka je \(A = U \Sigma V^T\) singularna dekompozicija matrice \(A \in \R^{n\times m}\) te neka je \(1 \leq k \leq \min\{m, n\}\). Tada je

\[\min_{r(X) = k} \| A - X \|_2 = \sigma_{k+1}.\]

Minimum se uzima po svim matricama \(X \in \R^{n \times m}\) ranga \(k\), a ostvaruje za matricu \(X = U_k \Sigma_k V_k^T\), gdje je \(\Sigma_k\) gornja lijeva podmatrica matrice \(\Sigma\) reda \(k\), a \(U_k\) i \(V_k\) matrice sastavljene od prvih \(k\) stupaca matrica \(U\) i \(V\).

Dokaz.

Neka je \(X\) proizvoljna matrica ranga \(k\). Zbog teorema o rangu i defektu je \(\dim \text{Ker}(X) = m-k\). Promotrimo potprostor \(\mathcal{V} = \text{span}\{v_{1}, \ldots, v_{k+1}\}\) razapet odgovarajućim desnim singularnim vektorima. Kako je \(\dim \text{Ker}(X) + \dim \mathcal{V} > m\), postoji jedinični vektor \(x \in \text{Ker}(X) \cap \mathcal{V}\). Stoga \(x = \sum_{j=1}^{k+1} \xi_j v_j\) za neke \(\xi_1, \ldots, \xi_{k+1} \in \R\) takve da je \(\sum_{j=1}^{k+1} \xi_j^2 = 1\). Slijedi:

\[\begin{split}\begin{align*} \|A-X\|_2^2 &\geq \|(A-X)x\|_2^2 = \|Ax\|_2^2 = \|U\Sigma V^T x\|_2^2 = \|\Sigma V^T x\|_2^2 \\ &= \left\| \mb{c} \sigma_1 v_1^T x \\ \sigma_2 v_2^T x \\ \vdots \\ \sigma_m v_m^T x \me \right\|_2^2 = \sum_{j=1}^m (\sigma_j v_j^T x)^2 = \sum_{j=1}^{k+1} \sigma_j^2 \xi_j^2 \\ &\geq \sum_{j=1}^{k+1} \sigma_{k+1}^2 \xi_j^2 \\ &= \sigma_{k+1}^2. \end{align*}\end{split}\]

Lako se vidi da se jednakost postiže za \(X\) iz iskaza teorema – rastavimo SVD od \(A\) ovako:

\[\begin{split} A = \mb{cc} U_k & U_{\perp} \me \cdot \mb{cc} \Sigma_k & 0 \\ 0 & \Sigma_{\perp} \me \cdot \mb{cc} V_k & V_{\perp} \me^T = U_k \Sigma_k V_k^T + U_{\perp} \Sigma_{\perp} V_{\perp}^T. \end{split}\]

Stoga je

\[ \|A - U_k \Sigma_k V_k^T\|_2 = \|U_{\perp} \Sigma_{\perp} V_{\perp}^T\|_2 = \|\Sigma_{\perp}\|_2 = \sigma_{k+1}. \]

Ponovno smo iskoristili da se 2-norma ne mijenja množimo li matricu slijeva ortonormalnom (ili zdesna transponiranom ortonormalnoj), te da je 2-norma dijagonalne matrice jednaka po modulu maksimalnom elementu dijagonale.

Sadržaj