Singulärvärdesuppdelning (SVD)

Varje matris — oavsett form — kan skrivas som tre enklare operationer efter varandra:

$A = U\Sigma V^\top$

En matrismultiplikation är alltså alltid: vrid, sträck, vrid. Inget annat.

Singulärvärdena i $\Sigma$ är sorterade i storleksordning och säger hur mycket matrisen sträcker i varje riktning. Är det femte värdet litet betyder det att matrisen knappt gör något i den femte riktningen — och den kan då kastas nästan gratis.

Det är hela idén med lågrangsapproximation: behåll de $k$ största singulärvärdena, släng resten.

För $A \in \mathbb{R}^{m\times n}$ med rang $r$:

$A = \sum_{i=1}^{r} \sigma_i\, u_i v_i^\top, \qquad \sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_r > 0$

Eckart–Young-satsen: den bästa rang-$k$-approximationen av $A$ (i Frobenius- och spektralnorm) fås genom att helt enkelt trunkera summan:

$A_k = \sum_{i=1}^{k} \sigma_i u_i v_i^\top, \qquad \|A - A_k\|_F^2 = \sum_{i>k}\sigma_i^2$

Det är ett anmärkningsvärt resultat: den optimala approximationen kräver ingen sökning, den läses av direkt.

Kompressionen: $A$ har $mn$ tal, $A_k$ har $k(m + n + 1)$. För en 1000×1000-matris med $k = 50$ är det 100 050 mot 1 000 000 — en tiondel.

Fyra användningar:

Kopplingen till egenvärden: $\sigma_i^2$ är egenvärdena till $A^\top A$, och $v_i$ dess egenvektorer. Men SVD finns för alla matriser, även icke-kvadratiska och singulära — till skillnad från egenvärdesuppdelning. Det är därför den är arbetshästen.

Beräkningen kostar $O(mn\min(m,n))$ för full SVD. Behöver du bara de $k$ största finns randomiserade metoder som är dramatiskt snabbare — sklearn.utils.extmath.randomized_svd eller scipy.sparse.linalg.svds.

import numpy as np

rng = np.random.default_rng(0)

# En matris som EGENTLIGEN har rang 3, plus lite brus
A = rng.normal(size=(200, 5)) @ rng.normal(size=(5, 150))
A = A[:, :3] @ rng.normal(size=(3, 150)) + 0.1 * rng.normal(size=(200, 150))

U, s, Vt = np.linalg.svd(A, full_matrices=False)
print(np.round(s[:8], 2))
# [419.4  406.63 242.5    2.61   2.5    2.47   2.44   2.42]
#   ↑ tre stora, sedan ett hopp ner till brusnivån

# Förklarad varians per komponent
andel = s**2 / (s**2).sum()
print(np.round(np.cumsum(andel)[:5], 4))     # [0.4394 0.8524 0.9993 0.9993 0.9993]

def trunkera(U, s, Vt, k):
    return U[:, :k] * s[:k] @ Vt[:k]

for k in (1, 3, 10, 50):
    Ak = trunkera(U, s, Vt, k)
    fel = np.linalg.norm(A - Ak) / np.linalg.norm(A)
    lagrat = k * (A.shape[0] + A.shape[1] + 1)
    print(f"k={k:>3}  relativt fel {fel:.4f}  lagring {lagrat / A.size:.1%} av originalet")
# k=  1  relativt fel 0.7487  lagring 1.2% av originalet
# k=  3  relativt fel 0.0268  lagring 3.5% av originalet
# k= 10  relativt fel 0.0248  lagring 11.7% av originalet

# Eckart–Young: felet är exakt summan av de bortkastade kvadrerade singulärvärdena
k = 3
print(round(float(np.linalg.norm(A - trunkera(U, s, Vt, k))**2), 4),
      round(float((s[k:]**2).sum()), 4))            # samma tal

# LoRA-idén: en rang-r uppdatering av ett stort lager
d, r = 4096, 8
print(f"fullt lager {d*d:,} parametrar, LoRA-rang {r}: {2*d*r:,} "
      f"({2*d*r/(d*d):.2%})")
# fullt lager 16,777,216 parametrar, LoRA-rang 8: 65,536 (0.39%)