SVD, pseudoinverse en PCA: SVD, pseudoinverse en PCA in R
SVD en beeldcompressie
R heeft de functie svd aan boord om numeriek de singulierewaardenontbinding te bepalen. Eerst rekenen we het voorbeeld uit de lestekst na:
> A <- matrix(c(2, 0, 0, -3), nrow=2, byrow=TRUE); A
[,1] [,2]
[1,] 2 0
[2,] 0 -3
> svdA <- svd(A) # singulierewaardendeontbinding (SVD) van A
> U <- svdA$u; U # matrix U uit singulierewaardendeontbinding
[,1] [,2]
[1,] 0 1
[2,] 1 0
> D <- svdA$d; D
[1] 3 2
> S <- diag(D) # matrix S uit singulierewaardenontbinding
> V <- svdA$v; V # matrix V uit singulierewaardenontbinding
[,1] [,2]
[1,] 0 1
[2,] -1 0
> U %*% S %*% t(V) - A # controle van svd(A)
[,1] [,2]
[1,] 0 0
[2,] 0 0
Het volgende voorbeeld van een singuliere matrix illustreert het verschil tussen exact en numeriek rekenen. \[\begin{aligned}\matrix{1 & -1\\ -1 & 1} &= \matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ \frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2}}\matrix{2 & 0\\ 0 & 0}\matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ -\frac{1}{2}\sqrt{2} & -\frac{1}{2}\sqrt{2}}\\ \\ &=\frac{1}{2}\sqrt{2} \matrix{-1 & 1\\ 1 & 1}\matrix{2 & 0\\ 0 & 0}\matrix{-1 & 1\\ -1 & -1}\frac{1}{2}\sqrt{2} \end{aligned} \]
> B <- matrix(c(1, -1, -1, 1), nrow=2, byrow=TRUE); B
[,1] [,2]
[1,] 1 -1
[2,] -1 1
> svdB <- svd(B) # singulierewaardenontbinding (SVD)
> U <- svdB$u; U # matrix U uit SVD
[,1] [,2]
[1,] -0.7071068 0.7071068
[2,] 0.7071068 0.7071068
> S <- diag(svdB$d); S # matrix S uit SVD
[,1] [,2]
[1,] 2 0
[2,] 0 0
> V <- svdB$v; V # matrix V uit SVD
[,1] [,2]
[1,] -0.7071068 0.7071068
[2,] 0.7071068 0.7071068
> U %*% S %*% t(V) - B
[,1] [,2]
[1,] -3.330669e-16 4.440892e-16
[2,] 0.000000e+00 -2.220446e-16
> round(.Last.value, 3) # afronden
[,1] [,2]
[1,] 0 0
[2,] 0 0
Nu een voorbeeld van beeldcompressie aan de hand van een patroon bestaande uit witte en zwarte hokjes op een \(9\times 9\) rooster. We specificeren de figuur met een matrix bestaande uit nullen (voor witte hokjes) en enen (voor zwarte hokjes): \[M = \matrix{ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0}\] In R kunnen we dit patroon visualiseren.
> M <- matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 1, 0, 1, 1, 1, 0, 1, 0,
1, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 0, 1, 1, 1, 0, 1, 0,
0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0), nrow=9, byrow=TRUE) # patroon
> dev.new(width=240, height=240) % vierkante grafische output
> # patroon als zwart-wit figuur (donker=255, licht=0)
> image(255*M, col=gray.colors(2, rev=TRUE), axes=FALSE)
Gebruik col=gray.colors(2, rev=TRUE, start=0, end=1) wanneer je een zwart-wit plaatje wenst.
Stel nu dat we iemand dit patroon willen sturen. Dan kunnen we de matrix \(M\) sturen, oftewel 81 getallen nul en één. Maar met de singulierewaardenontbinding kun je minder getallen sturen, terwijl het plaatje toch gereconstrueerd kan worden.
De eigenwaarden van de matrix \(M\) zijn in stijgende volgorde als volgt:
> ev <- eigen(M)
> round(sort(ev$values), 3)
[1] -2.000 -0.903 0.000 0.000 0.000 0.000 1.194
[8] 2.000 3.709
Dus zijn de singuliere waarden van \(M\) in dalende volgorde, afgerond op 1 decimaal, gelijk aan \(3.7\), \(2.0\), \(1.2\), \(0.9\) en \(0\).
> svdM <- svd(M)
> svdM$d
> round(svdM$d, 4)
[1] 3.7093 2.0000 2.0000 1.1939 0.9032 0.0000 0.0000
[8] 0.0000 0.0000
We benaderen \(M\) nu met de eerste termen van de singulierewaardenontbinding. Maar in plaats van de numerieke waarden van de ontbinding te gebruiken voor beeldweergave ronden we matrixelementen af op gehele getallen \(0\) en \(1\) om het patroon te reconstrueren.
Stel dat we alleen de hoogste singuliere waarde \(3.7\) gebruiken. Dan lijkt het plaatje al wat op de oorspronkelijke figuur:
> S1 <- diag(c(S[1,1],0,0,0,0,0,0,0,0))
> M1 <- U %*% S1 %*% t(V)
> dev.new(width=240, height=240)
NULL
> image(255*M1, col=gray.colors(2, rev=TRUE), axes=FALSE)Hiervoor hoeven we, in de verkorte notatie, maar \(1+2\times 9 = 19\) numerieke waarden op te sturen, namelijk de singuliere waarde en de twee vectoren \(\vec{u}_1\) en \(\vec{v}_1\) .
Als we ook de singuliere waarde \(2\) gebruiken, dan krijgen we de oorspronkelijke figuur terug.
> M3 <- U[,1:3] %*% S[1:3,1:3] %*% t(V[,1:3])
> dev.new(width=240, height=240)
> image(255*M3, col=gray.colors(2, rev=TRUE), axes=FALSE)
Nu volstaat het opsturen van \(3+6\times 9 = 57\) numerieke waarden, namelijk \(3\) eigenwaarden en \(2\times 3= 6\) vectoren. Dit lijkt niet zo veel minder dan \(81\), het gaat toch al om een reductie van ongeveer \(30\%\)
\(\phantom{xx}\)
De winst is spectaculairder bij het sturen van grotere digitale beelden, bijvoorbeeld op de in computer vision veel gebruikte foto van \(512\times 512\) pixels.
We lezen de foto in, berekenen de singulierewaardenontbinding van de bijpassende (\(512\times 512\))-matrix met grijswaarden, en benaderen de originele foto met 8, 16, 32, 64, en 128 termen uit de ontbinding.
> library(imager) # laad R pakket voor image processing
> Lena<- load.image("G:/Lena512.bmp") # lees bitmap foto in
> dev.new(height=512, width=512)
NULL
> plot(Lena, axes=FALSE) # toon de foto
> svdLena <- svd(scale(Lena)) # singulierewaardenontbinding van geschaalde matrix
> # definitie van functie om hoogste singuliere waarden alleen te gebruiken
> compressLena <- function(dim) {
U <- svdLena$u[,1:dim];
S <-diag(svdLena$d)[1:dim,1:dim];
V <- svdLena$v[,1:dim];
plot(as.cimg(U %*% S %*% t(V)), axes=FALSE)
}
> # toon verschillende gecomprimeerde foto's
> compressLena(8)
> compressLena(16)
> compressLena(32)
> compressLena(64)
> compressLena(128)
Hieronder zie je de beelden voor verschillende waarden van het aantal termen \(k\) die uit de singulierewaardenontbinding gehaald worden. Duidelijk is dat je minder informatie nodig hebt om toch een foto te hebben die op het zicht weinig afwijkt van het origineel (bij \(k=64\) is de compressiefactor ongeveer \(\frac{1}{4}\)). Bij kleurenfoto's is de winst zelfs groter.
| \(k=8\) | \(k=16\) |
 |
 |
| \(k=32\) | \(k=64\) |
 |
 |
| \(k=128\) | origineel: \(k=512\) |
 |
 |
