Este post es la PARTE-1 de dos posts. La PARTE-1 repasará la descomposición de matrices utilizando las técnicas de Singular Vector Decomposition (SVD), randomized Singular Vector Decomposition (rSVD) y Nonnegative Matrix Factorization (NMF) principalmente para la compresión de imágenes.
La PARTE-2 estudiará la eliminación del fondo de los vídeos utilizando randomized SVD y robust PCA.
Singular Vector Decomposition (ejemplo simple)
Singular Vector Decomposition (SVD) es una técnica matemática para descomponer una matriz en tres matrices separadas, cada una con propiedades especiales. En esencia, SVD transforma las variables correlacionadas en un conjunto de variables no correlacionadas.
SVD se puede utilizar para comprender mejor un conjunto de datos porque te indica el número de dimensiones importantes, así como también te puede simplificar el conjunto de datos reduciendo el número de columnas eliminando -desde el punto de vista del álgebra lineal- las que tienen dependencia lineal. Por lo tanto, SVD puede utilizarse para reducir un conjunto de datos de muchas dimensiones a un número menor de dimensiones, conservando al mismo tiempo la información importante.
Supongamos que A es una matriz con cada variable en una columna, y las observaciones en filas, entonces SVD está compuesto por:
A = U · d · VT
Donde:
- A: matriz original (m x n)
- U: matriz ortogonal (m × n)
- d: matriz diagonal (n × n)
- V: matriz ortogonal (n × n)
Empecemos con un pequeño ejemplo utilizando una matriz A: 10 x 3, con 3 variables (o columnas) y 10 registros (o filas).
|
1 2 3 4 5 6 7 |
set.seed(0) A <- matrix(runif(30, min=0, max=100), nrow=10) A <- data.frame(A) colnames(A) <- paste0("feature_",1:NCOL(A)) rownames(A) <- paste0("register",1:NROW(A)) A |
|
1 2 3 4 5 6 7 8 9 10 11 |
## feature_1 feature_2 feature_3 ## register1 89.66972 6.178627 77.744522 ## register2 26.55087 20.597457 93.470523 ## register3 37.21239 17.655675 21.214252 ## register4 57.28534 68.702285 65.167377 ## register5 90.82078 38.410372 12.555510 ## register6 20.16819 76.984142 26.722067 ## register7 89.83897 49.769924 38.611409 ## register8 94.46753 71.761851 1.339033 ## register9 66.07978 99.190609 38.238796 ## register10 62.91140 38.003518 86.969085 |
El siguiente código transformará la matriz A en 3 matrices (matriz diagonal: d, y las matrices u y v) de acuerdo con la descomposición del SVD y usando la función svd(). A continuación se pueden ver estas matrices:
|
1 2 3 |
A_svd <- svd(A) A_svd |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
## $d ## [1] 300.06826 113.25773 84.80736 ## ## $u ## [,1] [,2] [,3] ## [1,] -0.3450932 -0.39375868 -0.44297281 ## [2,] -0.2499063 -0.54713537 0.23660896 ## [3,] -0.1517606 -0.01260518 -0.09916967 ## [4,] -0.3595852 -0.06634029 0.27641129 ## [5,] -0.2982925 0.23555537 -0.40124974 ## [6,] -0.2256024 0.19268910 0.52579073 ## [7,] -0.3584251 0.09509590 -0.22232164 ## [8,] -0.3471264 0.47793821 -0.20284986 ## [9,] -0.3897184 0.28669365 0.36108378 ## [10,] -0.3541244 -0.35883174 0.05603342 ## ## $v ## [,1] [,2] [,3] ## [1,] -0.6948113 0.1875190 -0.6943154 ## [2,] -0.5277211 0.5229674 0.6693397 ## [3,] -0.4886182 -0.8314696 0.2644060 |
Hay que tener en cuenta algunos aspectos de los resultados de las matrices SVD:
- Cada fila de la matriz u se corresponde con cada fila de la matriz original.
- Las columnas de la matriz v coinciden con las columnas de la matriz original.
- Las columnas de la matriz u deben ser consideradas como topic_1, topic_2 y topic_3 y no coinciden directamente con las columnas originales de la matriz (feature_1, feature_2, feature_3).
- La matriz diagonal (d) se corresponde con cada tema, y está ordenada por la importancia del tema (el primer elemento es más importante y así sucesivamente).
De acuerdo con eso, podemos poner nombres a cada columna de modo aclaratorio:
|
1 2 3 4 5 |
U <- data.frame(A_svd$u) colnames(U) <- paste0("topic",1:NCOL(U)) rownames(U) <- paste0("register",1:NROW(A)) U |
|
1 2 3 4 5 6 7 8 9 10 11 |
## topic1 topic2 topic3 ## register1 -0.3450932 -0.39375868 -0.44297281 ## register2 -0.2499063 -0.54713537 0.23660896 ## register3 -0.1517606 -0.01260518 -0.09916967 ## register4 -0.3595852 -0.06634029 0.27641129 ## register5 -0.2982925 0.23555537 -0.40124974 ## register6 -0.2256024 0.19268910 0.52579073 ## register7 -0.3584251 0.09509590 -0.22232164 ## register8 -0.3471264 0.47793821 -0.20284986 ## register9 -0.3897184 0.28669365 0.36108378 ## register10 -0.3541244 -0.35883174 0.05603342 |
|
1 2 3 4 5 |
D <- data.frame(diag(A_svd$d)) colnames(D) <- paste0("topic",1:NCOL(U)) rownames(D) <- paste0("topic",1:NCOL(U)) D |
|
1 2 3 4 |
## topic1 topic2 topic3 ## topic1 300.0683 0.0000 0.00000 ## topic2 0.0000 113.2577 0.00000 ## topic3 0.0000 0.0000 84.80736 |
|
1 2 3 4 5 |
V <- data.frame(A_svd$v) colnames(V) <- colnames(A) rownames(V) <- paste0("topic",1:NCOL(U)) V |
|
1 2 3 4 |
## feature_1 feature_2 feature_3 ## topic1 -0.6948113 0.1875190 -0.6943154 ## topic2 -0.5277211 0.5229674 0.6693397 ## topic3 -0.4886182 -0.8314696 0.2644060 |
Observa la correlación entre las columnas de la matriz u:
|
1 2 |
cor(A_svd$u) |
|
1 2 3 4 |
## [,1] [,2] [,3] ## [1,] 1.0000000 -0.1230321288 0.1185052849 ## [2,] -0.1230321 1.0000000000 0.0007930161 ## [3,] 0.1185053 0.0007930161 1.0000000000 |
Como ejemplo, y sólo para ver las diferencias entre las características o features (de la matriz original) y los topics, podemos replicar la matriz A en el objeto B, y luego poner ceros en la feature_3:
|
1 2 3 |
B <- A B[,3] <- 0 |
Ahora podemos calcular la SVD de la matriz B. Como se puede ver abajo, los valores cero son visibles sólo en la matriz diagonal, y en la matriz v, pero no en la matriz u.
|
1 2 3 |
SVD_B <- svd(B) SVD_B |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
## $d ## [1] 267.38125 87.74893 0.00000 ## ## $u ## [,1] [,2] [,3] ## [1,] -0.2787209 0.5726378541 -0.172499723 ## [2,] -0.1256765 0.0008567322 -0.012160150 ## [3,] -0.1503636 0.1019880483 0.979107903 ## [4,] -0.3269347 -0.2162331156 -0.031176728 ## [5,] -0.3562184 0.2909811492 -0.051100385 ## [6,] -0.2364834 -0.5506977253 -0.009723958 ## [7,] -0.3794382 0.1819959212 -0.050268679 ## [8,] -0.4436522 0.0167429159 -0.052397739 ## [9,] -0.4229711 -0.4286735662 -0.035489988 ## [10,] -0.2729534 0.0991208907 -0.035126831 ## ## $v ## [,1] [,2] [,3] ## [1,] -0.7887458 0.6147195 0 ## [2,] -0.6147195 -0.7887458 0 ## [3,] 0.0000000 0.0000000 1 |
Ahora podemos reconstruir la matriz B original usando las matrices u,v y d y calculando B=u·d·vT.
|
1 2 3 |
B_ <- SVD_B$u %*% diag(SVD_B$d) %*% t(SVD_B$v) B_ |
|
1 2 3 4 5 6 7 8 9 10 11 |
## [,1] [,2] [,3] ## [1,] 89.66972 6.178627 0 ## [2,] 26.55087 20.597457 0 ## [3,] 37.21239 17.655675 0 ## [4,] 57.28534 68.702285 0 ## [5,] 90.82078 38.410372 0 ## [6,] 20.16819 76.984142 0 ## [7,] 89.83897 49.769924 0 ## [8,] 94.46753 71.761851 0 ## [9,] 66.07978 99.190609 0 ## [10,] 62.91140 38.003518 0 |
Ahora podemos comprobar si la matriz B original y la matriz B_ reconstruida son iguales (o casi iguales) simplemente restando las matrices. Se puede ver que las dos matrices coinciden muy bien.
|
1 2 3 |
options(scipen=99) B - B_ |
|
1 2 3 4 5 6 7 8 9 10 11 |
## feature_1 feature_2 feature_3 ## register1 0.000000000000000000000 -0.000000000000056843419 0 ## register2 -0.000000000000003552714 -0.000000000000007105427 0 ## register3 -0.000000000000007105427 -0.000000000000007105427 0 ## register4 -0.000000000000021316282 -0.000000000000028421709 0 ## register5 0.000000000000000000000 0.000000000000000000000 0 ## register6 0.000000000000014210855 -0.000000000000028421709 0 ## register7 -0.000000000000014210855 -0.000000000000007105427 0 ## register8 -0.000000000000014210855 -0.000000000000014210855 0 ## register9 -0.000000000000014210855 0.000000000000000000000 0 ## register10 -0.000000000000014210855 -0.000000000000014210855 0 |
También podemos calcular el RMSE de todos los valores:
|
1 2 |
sqrt(sum((B - B_)^2)/NROW(B)) |
|
1 |
## [1] 0.00000000000002632367 |
Singular Vector Decomposition (ejemplo con imagen)
Ahora supongamos que queremos reducir el tamaño de una imagen comprimiéndola. ¡Podemos usar SVD para este trabajo!
En primer lugar, necesitamos instalar/cargar el paquete rsvd, que además contiene una imagen de ejemplo de un tigre que utilizaremos para nuestro propósito, las dimensiones de la imagen son en escala de grises de 1600 x 1200 (8 bits[0-255]/255).
|
1 2 3 4 5 |
# loading example image library(rsvd) data(tiger) dim(tiger) |
|
1 |
## [1] 1600 1200 |
El siguiente paso es utilizar svd() a la imagen tiger, por lo que obtenemos nuestras tres matrices: u, d, y v.
|
1 2 3 4 |
# Image compression using SVD SVD <- svd(tiger) str(SVD) |
|
1 2 3 4 |
## List of 3 ## $ d: num [1:1200] 528 119.8 86.1 77 66.1 ... ## $ u: num [1:1600, 1:1200] -0.0134 -0.0133 -0.0133 -0.0134 -0.0133 ... ## $ v: num [1:1200, 1:1200] -0.0357 -0.0357 -0.0355 -0.0356 -0.0357 ... |
A continuación, podemos representar en una gráfica la importancia de los “topics”, recordemos que están ordenados por orden de importancia, por lo que en el lado izquierdo de la gráfica se sitúan los más importantes y el lado derecho los menos importantes:
|
1 2 3 |
# topics importance plot(SVD$d^2/sum(SVD$d^2), type="b", pch=16, xlab="topics", ylab="variance explained") |
Para este ejemplo, seleccionaremos los 100 “topics” más importantes. Recuerde que los topics están ordenados, así que es tan sencillo como seleccionar las 100 primeras columnas de las matrices.
Así que de 1200 topics, sólo escogeremos los primeros 100 y luego reconstruiremos la imagen usando las matrices reducidas:
|
1 2 3 |
levels <- 100 # the number of topics we will use tiger.SVD = SVD$u[,1:levels] %*% diag(SVD$d[1:levels]) %*% t(SVD$v[,1:levels]) # reconstructing the image |
En esencia, estamos usando sólo la parte naranja de las matrices para reconstruir la matriz A. Ver el diagrama de abajo.
Finalmente, podemos mostrar la imagen original y la imagen reducida:
|
1 2 3 4 5 6 7 8 9 10 |
# Display orginal and reconstrucuted image # imager library plots faster than image() funtion # imager vignette: https://dahtah.github.io/imager/imager.html library(imager) IM_1 <- as.cimg(tiger[,1200:1]) IM_2 <- as.cimg(tiger.SVD[,1200:1]) plot(imappend(list(IM_1,IM_2),axis="x")) text(10,50,"original", pos=4, col="white") text(1610,50,"SVD reduced [100]", pos=4, col="white") |
También podemos estimar el tamaño de las matriz reducida, ésta se ha reducido hasta ocupar un 9,52% de la original.
|
1 2 3 |
#Size of original matrices: print(object.size(SVD$u)+object.size(SVD$d)+object.size(SVD$v), units="Mb") |
|
1 |
## 25.6 Mb |
|
1 2 3 |
#Size of reduced matrices (SVD [100]): print(object.size(SVD$u[,1:levels])+object.size(SVD$d[1:levels])+object.size(SVD$u[,1:levels]), units="Mb") |
|
1 |
## 2.4 Mb |
Randomized SVD (ejemplo con imagen)
En el apartado anterior, utilizamos el método SVD para comprimir las imágenes. También podemos usar el método randomized SVD (usando el paquete rsvd) que es equivalente al SVD clásico pero más rápido. Puedes echarle un vistazo en PDF para detalles del benchmark tiempos de proceso.
Una simple comparación de velocidad entre los dos métodos se muestra a continuación.
|
1 2 3 4 5 |
SVD_time <- system.time({ X <- svd(tiger) }) rSVD_time <- system.time({ X <- rsvd(tiger, k=100) }) # classical SVD time calculation: SVD_time |
|
1 2 |
## user system elapsed ## 8.984 0.050 9.041 |
|
1 2 3 |
# randomized SVD [100] time calculation: rSVD_time |
|
1 2 |
## user system elapsed ## 1.594 0.030 1.624 |
Como puedes ver, en este caso, el randomize SVD es mucho más rápida que el SVD clásico. A continuación, veremos cómo funciona rSVD, ten en cuenta que tienes que indicarle a la función el número de topics o el low-rank decomposition (el argumento k de la función) de antemano.
|
1 2 3 4 |
# Image compression using randomized SVD rSVD <- rsvd(tiger, k=100) tiger.rSVD = rSVD$u %*% diag(rSVD$d) %*% t(rSVD$v) # reconstruct image |
Comparando los resultados gráficos del SVD y rSVD, se ven similares.
|
1 2 3 4 5 6 7 |
# Display orginal, and reconstrucuted images using SVD and rSVD IM_3 <- as.cimg(tiger.rSVD[,1200:1]) plot(imappend(list(IM_1,IM_2,IM_3),axis="x"), axes=FALSE) text(10,60,"original", pos=4, col="white", cex=0.8) text(1610,60,"SVD [100]", pos=4, col="white", cex=0.8) text(3210,60,"rSVD [100]", pos=4, col="white", cex=0.8) |
Calcularemos ahora el error (RMSE) de la descomposición en función del valor de k:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
rSVD2 <- list() tiger.rSVD_2 <- list() result <- list() for(k in c(2,5,10,20,50,100,200,500,1200)){ cat("\nk:",k,"/",1200) rSVD2[[k]] <- rsvd(tiger, k=k) size <- object.size(rSVD2[[k]]$u[,1:k])+object.size(rSVD2[[k]]$d[1:k])+object.size(rSVD2[[k]]$u[,1:k]) # object size ## tiger.rSVD_2[[k]] <- rSVD2[[k]]$u %*% diag(rSVD2[[k]]$d) %*% t(rSVD2[[k]]$v) # reconstructing image error <- sqrt(sum((tiger - tiger.rSVD_2[[k]])^2)/NROW(as.vector(tiger))) # calculating error result[[k]] <- data.frame(k=k, error=round(error,2), size_Mb=as.numeric(as.character(round(size/1024/1024,2)))) } |
|
1 2 3 4 5 6 7 8 9 10 |
## ## k: 2 / 1200 ## k: 5 / 1200 ## k: 10 / 1200 ## k: 20 / 1200 ## k: 50 / 1200 ## k: 100 / 1200 ## k: 200 / 1200 ## k: 500 / 1200 ## k: 1200 / 1200 |
|
1 2 |
result <- do.call("rbind",result) |
Como se puede observar a continuación, incrementando el valor de k (low-rank decomposition) disminuye el RMSE y aumenta el tamaño que ocupan en memoria las matrices:
|
1 2 3 4 |
par(mfrow=c(1,2)) plot(result$k, result$error, type="b", pch=16, ylab="error (RMSE)", xlab="k") plot(result$k, result$size_Mb, type="b", pch=16, ylab="matrices mem size (Mb)", xlab="k", col="blue") |
A continuación observaremos las diferentes imágenes resultantes de diferentes valores de k:
|
1 2 3 4 5 6 7 8 9 10 |
# Display orginal and reconstrucuted image IM <- list() for(k in result$k) IM[[k]] <- as.cimg(tiger.rSVD_2[[k]][,1200:1]) ## par(mfrow=c(3,3), mar=c(0,0,0,0)) for(i in 1:NROW(result)){ plot(IM[[result[i,"k"]]], axes=FALSE) text(10,60,paste0("k = ",result[i,"k"],". mem size = ",result[i,"size_Mb"], " Mb"), pos=4, col="white") } |
Non-negative Matrix Factorization (ejemplo con imagen)
Non-negative matrix factorization (NMF) es otro método de descomposición matricial. Descompone una matriz V en dos matrices W y H:
V = W · H
Las columnas de W se llaman vectores de base y las columnas de H representan las codificaciones, que están en correspondencia uno a uno con los vectores de base de W.
Una gran diferencia con respecto al SVD es que las tres matrices no tienen elementos negativos. SDV puede tener valores negativos, en NMF no es posible tener valores negativos, los valores no negativos son más fáciles de interpretar que los valores negativos, por lo que ésta es la razón principal de la popularidad de NMF sobre SVD. En contrapartida, usando NMF no se podrá recuperar la matriz original si ésta originalmente tiene números negativos.
Puede obtener más información sobre el paquete NMF en PDF.
|
1 2 3 4 5 6 |
# Image compression using NMF r <- 100 library(NMF) NMFd_time <- system.time({ NMFd <- nmf(tiger, rank=r) }) tiger.NMF <- NMFd@fit@W %*% NMFd@fit@H # reconstructing image |
Hay que tener en cuenta que el NMF es computacionalmente intensivo en el cálculo y se vuelve inviable para datos masivos. Podemos ver abajo que el tiempo consumido para procesar nuestra imagen tiger es enorme, unas 900 veces más lento que rSVD.
|
1 2 3 |
# NMF (k=100) time calculation: NMFd_time |
|
1 2 |
## user system elapsed ## 2157.356 19.122 2178.343 |
A continuación podemos encontrar la imagen original, y las tres transformaciones utilizadas en este post con una “low-rank decomposition” de 100.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
# Display orginal, and reconstrucuted images using SVD and rSVD IM_4 <- as.cimg(tiger.NMF[,1200:1]) par(mfrow=c(1,1), mar=c(0,0,0,0)) plot(imappend(list(IM_1,IM_2,IM_3,IM_4),axis="x"), axes=FALSE) text(10,60,"original", pos=4, col="white", cex=0.8) text(1610,60,"SVD [100]", pos=4, col="white", cex=0.8) text(3210,60,"rSVD [100]", pos=4, col="white", cex=0.8) text(4810,60,"NMF [100]", pos=4, col="white", cex=0.8) # calculating error text(1610,1000,paste0("SVD RMSE error = ",round(sqrt(sum((tiger - tiger.SVD)^2)/NROW(as.vector(tiger))),6)), pos=4, col="white", cex=0.8) text(3210,1000,paste0("rSVD RMSE error = ",round(sqrt(sum((tiger - tiger.rSVD)^2)/NROW(as.vector(tiger))),6)), pos=4, col="white", cex=0.8) text(4810,1000,paste0("NMF RMSE error = ",round(sqrt(sum((tiger - tiger.NMF)^2)/NROW(as.vector(tiger))),6)), pos=4, col="white", cex=0.8) # calculating size SVD_size <- object.size(SVD$u[,1:levels])+object.size(SVD$d[1:levels])+object.size(SVD$u[,1:levels]) rSVD_size <- object.size(rSVD$u[,1:levels])+object.size(rSVD$d[1:levels])+object.size(rSVD$u[,1:levels]) NMF_size <- object.size(NMFd@fit@W)+object.size(NMFd@fit@H) text(1610,1100,paste0("size = ",round(as.numeric(SVD_size)/1024/1024,2)," Mb"), pos=4, col="white", cex=0.8) text(3210,1100,paste0("size = ",round(as.numeric(rSVD_size)/1024/1024,2)," Mb"), pos=4, col="white", cex=0.8) text(4810,1100,paste0("size = ",round(as.numeric(NMF_size)/1024/1024,2)," Mb"), pos=4, col="white", cex=0.8) |
En este post hemos visto los métodos de descomposición de matrices de SVD, rSVD y NMF aplicados a imágenes. En la PARTE-2 utilizaremos estas técnicas a videos de videovigilancia para la eliminación del fondo de los vídeos y así separar el fondo de los objetos en movimiento.
Nos vemos en el próximo post!
Session Info:
|
1 2 3 |
------------------------------------ Total R execution time: 39 mins ------------------------------------ |
|
1 2 3 4 5 6 7 8 9 |
setting value version R version 3.4.3 (2017-11-30) os macOS High Sierra 10.13.5 system x86_64, darwin15.6.0 ui RStudio language (EN) collate es_ES.UTF-8 tz Europe/Madrid date 2018-08-18 |
|
1 2 |
------------------------------------ Packages: |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
[1] "bitops - 1.0-6 - 2013-08-17 - CRAN (R 3.4.0)" [2] "cluster - 2.0.6 - 2017-03-10 - CRAN (R 3.4.3)" [3] "imager - 0.40.2 - 2017-04-24 - CRAN (R 3.4.0)" [4] "magrittr - 1.5 - 2014-11-22 - CRAN (R 3.4.0)" [5] "NMF - 0.21.0 - 2018-03-06 - CRAN (R 3.4.4)" [6] "pkgmaker - 0.22 - 2014-05-14 - CRAN (R 3.4.0)" [7] "plyr - 1.8.4 - 2016-06-08 - CRAN (R 3.4.0)" [8] "RCurl - 1.95-4.11 - 2018-07-15 - cran (@1.95-4.)" [9] "registry - 0.3 - 2015-07-08 - CRAN (R 3.4.0)" [10] "reshape2 - 1.4.2 - 2016-10-22 - CRAN (R 3.4.0)" [11] "rngtools - 1.2.4 - 2014-03-06 - CRAN (R 3.4.0)" [12] "rsvd - 0.9 - 2017-12-08 - CRAN (R 3.4.3)" [13] "RWordPress - 0.2-3 - 2018-03-04 - Github (duncantl/RWordPress@ce6d2d6)" [14] "stringr - 1.2.0 - 2017-02-18 - CRAN (R 3.4.0)" [15] "XMLRPC - 0.3-1 - 2018-08-17 - Github (duncantl/XMLRPC@add9496)" |
Appendix, all the code:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 |
set.seed(0) A <- matrix(runif(30, min = 0, max = 100), nrow = 10) A <- data.frame(A) colnames(A) <- paste0("feature_", 1:NCOL(A)) rownames(A) <- paste0("register", 1:NROW(A)) A A_svd <- svd(A) A_svd U <- data.frame(A_svd$u) colnames(U) <- paste0("topic", 1:NCOL(U)) rownames(U) <- paste0("register", 1:NROW(A)) U D <- data.frame(diag(A_svd$d)) colnames(D) <- paste0("topic", 1:NCOL(U)) rownames(D) <- paste0("topic", 1:NCOL(U)) D V <- data.frame(A_svd$v) colnames(V) <- colnames(A) rownames(V) <- paste0("topic", 1:NCOL(U)) V cor(A_svd$u) B <- A B[, 3] <- 0 SVD_B <- svd(B) SVD_B B_ <- SVD_B$u %*% diag(SVD_B$d) %*% t(SVD_B$v) B_ options(scipen = 99) B - B_ sqrt(sum((B - B_)^2)/NROW(B)) # loading example image library(rsvd) data(tiger) dim(tiger) # Image compression using SVD SVD <- svd(tiger) str(SVD) # topics importance plot(SVD$d^2/sum(SVD$d^2), type = "b", pch = 16, xlab = "topics", ylab = "variance explained") levels <- 100 # the number of topics we will use tiger.SVD = SVD$u[, 1:levels] %*% diag(SVD$d[1:levels]) %*% t(SVD$v[, 1:levels]) # reconstructing the image # Display orginal and reconstrucuted image imager library plots faster than # image() funtion imager vignette: # https://dahtah.github.io/imager/imager.html library(imager) IM_1 <- as.cimg(tiger[, 1200:1]) IM_2 <- as.cimg(tiger.SVD[, 1200:1]) plot(imappend(list(IM_1, IM_2), axis = "x")) text(10, 50, "original", pos = 4, col = "white") text(1610, 50, "SVD reduced [100]", pos = 4, col = "white") # Size of original matrices: print(object.size(SVD$u) + object.size(SVD$d) + object.size(SVD$v), units = "Mb") # Size of reduced matrices (SVD [100]): print(object.size(SVD$u[, 1:levels]) + object.size(SVD$d[1:levels]) + object.size(SVD$u[, 1:levels]), units = "Mb") SVD_time <- system.time({ X <- svd(tiger) }) rSVD_time <- system.time({ X <- rsvd(tiger, k = 100) }) # classical SVD time calculation: SVD_time # randomized SVD [100] time calculation: rSVD_time # Image compression using randomized SVD rSVD <- rsvd(tiger, k = 100) tiger.rSVD = rSVD$u %*% diag(rSVD$d) %*% t(rSVD$v) # reconstruct image # Display orginal, and reconstrucuted images using SVD and rSVD IM_3 <- as.cimg(tiger.rSVD[, 1200:1]) plot(imappend(list(IM_1, IM_2, IM_3), axis = "x"), axes = FALSE) text(10, 60, "original", pos = 4, col = "white", cex = 0.8) text(1610, 60, "SVD [100]", pos = 4, col = "white", cex = 0.8) text(3210, 60, "rSVD [100]", pos = 4, col = "white", cex = 0.8) rSVD2 <- list() tiger.rSVD_2 <- list() result <- list() for (k in c(2, 5, 10, 20, 50, 100, 200, 500, 1200)) { cat("\nk:", k, "/", 1200) rSVD2[[k]] <- rsvd(tiger, k = k) size <- object.size(rSVD2[[k]]$u[, 1:k]) + object.size(rSVD2[[k]]$d[1:k]) + object.size(rSVD2[[k]]$u[, 1:k]) # object size ## tiger.rSVD_2[[k]] <- rSVD2[[k]]$u %*% diag(rSVD2[[k]]$d) %*% t(rSVD2[[k]]$v) # reconstructing image error <- sqrt(sum((tiger - tiger.rSVD_2[[k]])^2)/NROW(as.vector(tiger))) # calculating error result[[k]] <- data.frame(k = k, error = round(error, 2), size_Mb = as.numeric(as.character(round(size/1024/1024, 2)))) } result <- do.call("rbind", result) par(mfrow = c(1, 2)) plot(result$k, result$error, type = "b", pch = 16, ylab = "error (RMSE)", xlab = "k") plot(result$k, result$size_Mb, type = "b", pch = 16, ylab = "matrices mem size (Mb)", xlab = "k", col = "blue") # Display orginal and reconstrucuted image IM <- list() for (k in result$k) IM[[k]] <- as.cimg(tiger.rSVD_2[[k]][, 1200:1]) ## par(mfrow = c(3, 3), mar = c(0, 0, 0, 0)) for (i in 1:NROW(result)) { plot(IM[[result[i, "k"]]], axes = FALSE) text(10, 60, paste0("k = ", result[i, "k"], ". mem size = ", result[i, "size_Mb"], " Mb"), pos = 4, col = "white") } # Image compression using NMF r <- 100 library(NMF) NMFd_time <- system.time({ NMFd <- nmf(tiger, rank = r) }) tiger.NMF <- NMFd@fit@W %*% NMFd@fit@H # reconstructing image # NMF (k=100) time calculation: NMFd_time # Display orginal, and reconstrucuted images using SVD and rSVD IM_4 <- as.cimg(tiger.NMF[, 1200:1]) par(mfrow = c(1, 1), mar = c(0, 0, 0, 0)) plot(imappend(list(IM_1, IM_2, IM_3, IM_4), axis = "x"), axes = FALSE) text(10, 60, "original", pos = 4, col = "white", cex = 0.8) text(1610, 60, "SVD [100]", pos = 4, col = "white", cex = 0.8) text(3210, 60, "rSVD [100]", pos = 4, col = "white", cex = 0.8) text(4810, 60, "NMF [100]", pos = 4, col = "white", cex = 0.8) # calculating error text(1610, 1000, paste0("SVD RMSE error = ", round(sqrt(sum((tiger - tiger.SVD)^2)/NROW(as.vector(tiger))), 6)), pos = 4, col = "white", cex = 0.8) text(3210, 1000, paste0("rSVD RMSE error = ", round(sqrt(sum((tiger - tiger.rSVD)^2)/NROW(as.vector(tiger))), 6)), pos = 4, col = "white", cex = 0.8) text(4810, 1000, paste0("NMF RMSE error = ", round(sqrt(sum((tiger - tiger.NMF)^2)/NROW(as.vector(tiger))), 6)), pos = 4, col = "white", cex = 0.8) # calculating size SVD_size <- object.size(SVD$u[, 1:levels]) + object.size(SVD$d[1:levels]) + object.size(SVD$u[, 1:levels]) rSVD_size <- object.size(rSVD$u[, 1:levels]) + object.size(rSVD$d[1:levels]) + object.size(rSVD$u[, 1:levels]) NMF_size <- object.size(NMFd@fit@W) + object.size(NMFd@fit@H) text(1610, 1100, paste0("size = ", round(as.numeric(SVD_size)/1024/1024, 2), " Mb"), pos = 4, col = "white", cex = 0.8) text(3210, 1100, paste0("size = ", round(as.numeric(rSVD_size)/1024/1024, 2), " Mb"), pos = 4, col = "white", cex = 0.8) text(4810, 1100, paste0("size = ", round(as.numeric(NMF_size)/1024/1024, 2), " Mb"), pos = 4, col = "white", cex = 0.8) |