第6章 行列分解

6.1 特異値分解

例53

Z <- matrix(c(0, 5, -1, -2, -4, 1), nrow = 3)
Z

##      [,1] [,2]
## [1,]    0   -2
## [2,]    5   -4
## [3,]   -1    1

svd(Z)

## $d
## [1] 6.681937 1.533530
## 
## $u
##            [,1]        [,2]
## [1,] -0.1987442  0.97518046
## [2,] -0.9570074 -0.21457290
## [3,]  0.2112759 -0.05460347
## 
## $v
##            [,1]       [,2]
## [1,] -0.7477342 -0.6639982
## [2,]  0.6639982 -0.7477342

svd(t(Z))

## $d
## [1] 6.681937 1.533530
## 
## $u
##            [,1]      [,2]
## [1,] -0.7477342 0.6639982
## [2,]  0.6639982 0.7477342
## 
## $v
##            [,1]        [,2]
## [1,] -0.1987442 -0.97518046
## [2,] -0.9570074  0.21457290
## [3,]  0.2112759  0.05460347

例54

Z <- matrix(c(0, 5, 5, -1), nrow = 2)
Z

##      [,1] [,2]
## [1,]    0    5
## [2,]    5   -1

svd(Z)

## $d
## [1] 5.524938 4.524938
## 
## $u
##            [,1]       [,2]
## [1,]  0.6710053 -0.7414525
## [2,] -0.7414525 -0.6710053
## 
## $v
##            [,1]       [,2]
## [1,] -0.6710053 -0.7414525
## [2,]  0.7414525 -0.6710053

eigen(Z)

## eigen() decomposition
## $values
## [1]  4.524938 -5.524938
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7414525 -0.6710053
## [2,] -0.6710053  0.7414525

6.2 Eckart-Youngの定理

svd.r <- function(z, r) {
  n <- min(nrow(z), ncol(z))
  ss <- svd(z)
  tt <- ss$u %*% diag(c(ss$d[1:r], rep(0, n - r))) %*% t(ss$v)
  return(tt)
}

例55

m <- 200
n <- 150
z <- matrix(rnorm(m * n), nrow = m)
F.norm <- NULL
for (r in 1:n) {
  m <- svd.r(z, r)
  F.norm <- c(F.norm, norm(z - m, "F") ^ 2)
}
plot(1:n, F.norm, type = "l", xlab = "階数", ylab = "Frobeniusノルムの二乗")

例56

library(jpeg)

## Warning: package 'jpeg' was built under R version 4.0.3

image <- readJPEG('lion.jpg')
rank.seq <- c(2, 5, 10, 20, 50, 100)
mat <- array(0, dim = c(nrow(image), ncol(image), 3))
for (j in rank.seq) {
  for (i in 1:3)
    mat[, , i] <- svd.r(image[, , i], j)
  writeJPEG(mat, paste("compressed/lion_compressed", "_mat_rank_", j, ".jpg", sep = ""))
}

mat.r <- function(z, mask, r) {
  z <- as.matrix(z)
  min <- Inf
  m <- nrow(z)
  n <- ncol(z)
  for (j in 1:5) {
    guess <- matrix(rnorm(m * n), nrow = m)
    for (i in 1:10)
      guess <- svd.r(mask * z + (1 - mask) * guess, r)
    value <- norm(mask * (z - guess), "F")
    if (value < min) {
      min.mat <- guess
      min <- value
    }
  }
  return(min.mat)
}

例57

library(jpeg)
image <- readJPEG('lion.jpg')
m <- nrow(image)
n <- ncol(image)
mask <- matrix(rbinom(m * n, 1, 0.5), nrow = m)
rank.seq <- c(2, 5, 10, 20, 50, 100)
mat <- array(0, dim = c(nrow(image), ncol(image), 3))
for (j in rank.seq) {
  for (i in 1:3)
    mat[, , i] <- mat.r(image[, , i], mask, j)
  writeJPEG(mat, paste("compressed/lion_compressed", "_mat_rank_", j, ".jpg", sep = ""))
}

6.3 ノルム

6.4 低階数近似のスパースの適用

soft.svd <- function(lambda, z) {
  n <- ncol(z)
  ss <- svd(z)
  dd <- pmax(ss$d - lambda, 0)
  return(ss$u %*% diag(dd) %*% t(ss$v))
}

mat.lasso <- function(lambda, z, mask) {
  z <- as.matrix(z)
  m <- nrow(z)
  n <- ncol(z)
  guess <- matrix(rnorm(m * n), nrow = m)
  for (i in 1:20)
    guess <- soft.svd(lambda, mask * z + (1 - mask) * guess)
  return(guess)
}

例59

library(jpeg)
image <- readJPEG('lion.jpg')
m <- nrow(image[, , 1])
n <- ncol(image[, , 1])
p <- 0.5
lambda <- 0.5
mat <- array(0, dim = c(m, n, 3))
mask <- matrix(rbinom(m * n, 1, p), ncol = n)
for (i in 1:3)
  mat[, , i] <- mat.lasso(lambda, image[, , i], mask)
writeJPEG(mat, paste("compressed/lion_compressed", "_mat_soft.jpg", sep = ""))