行列の低ランク近似と近似度

特異値分解を使った「行列の低ランク近似」に関するメモ。

まずは、視覚的に確認しやすい画像データで試してみる。

高さ200px、幅320px、インデックスカラーのRAWデータを用意する。
f:id:ichou1:20181010153239p:plain

1pxを1バイトで表現、0〜255の範囲の値をとる200行 x 320列の行列ができる。

以下は左上を拡大した画像。
先頭6バイト分のみ値を添えてみた。値が小さいほど黒に近く、値が大きいほど白に近くなる。
f:id:ichou1:20181010160014p:plain

CSVに落としたものをRに読み込み、

A <- data.matrix(read.table("200_320_gray_scale.csv", header=F, sep=","))

ランク数を「r」とおいて、特異値分解する。

A.svd <- svd(A)

u <- A.svd$u[,1:r, drop=F]
d <- diag(A.svdA$d)[1:r, 1:r, drop=F]
v <- A.svd$v[,1:r, drop=F]

そして、再度、csvに書き戻す。

write.table(round(u %*% d %*% t(v)),
                 str_c("200_320_gray_scale_low_rank_approximate_", r), ".csv"),
                 sep=",", row.names=FALSE, col.names=FALSE)
r = 1

f:id:ichou1:20181010161502p:plain

r = 1 (拡大)

f:id:ichou1:20181010162151p:plain

r = 10

f:id:ichou1:20181010162446p:plain

r = 10 (拡大)

f:id:ichou1:20181010163421p:plain

r = 20

f:id:ichou1:20181010162541p:plain

r = 20 (拡大)

f:id:ichou1:20181010163239p:plain


続けて、フロベニウスノルム(全成分の二乗和の平方根)を確認してみる。

フロベニウスノルム(元々の行列)
norm(A, "F")  # 26607.34 

続けて、近似した行列のフロベニウスノルム

フロベニウスノルム(ランク1、5、10、15、20、25、50、75、100で近似した行列)
norm(A_rank1,   "F")  # 25696.45
norm(A_rank5,   "F")  # 26186.34
norm(A_rank10,  "F")  # 26333.71
norm(A_rank15,  "F")  # 26412.39
norm(A_rank20,  "F")  # 26460.39
norm(A_rank25,  "F")  # 26493.76
norm(A_rank50,  "F")  # 26570.16
norm(A_rank75,  "F")  # 26594.95
norm(A_rank100, "F")  # 26603.22

近似に使うランク数を大きくするほど、フロベニウスノルムも大きくなっていることが確認できる。
尚、特異値の上位20は以下のとおりであった。

> A.svd$d[1:20]
 [1] 25695.8853  3754.6842  2374.1272  1782.2953  1588.0870  1466.8724
 [7]  1308.7928  1248.0294  1112.3381  1042.8967   983.9771   963.9730
[13]   916.1656   841.3105   837.3214   794.9816   735.2430   722.9311
[19]   685.1815   634.1618

下位20は以下のとおり。

> A.svd$d[180:200]
 [1] 15.898742 15.417982 14.987853 14.435300 14.174012 13.721447 13.289164
 [8] 12.747505 12.436280 11.677728 11.473563 11.297093 10.792612 10.623457
[15]  9.714841  9.523776  9.005887  8.843195  8.200154  8.091183  7.484766
各ランクの特異値(ランク1〜200)

f:id:ichou1:20181012210530p:plain

各ランクの特異値(ランク11〜200)

f:id:ichou1:20181012213726p:plain


(参考)特異値を全部足すと以下のとおり。

> sum(A.svd$d)
[1] 75484.84
近似に使うランク数とフロベニウスノルム

f:id:ichou1:20181012210506p:plain


フロベニウスノルムは以下のやり方でも求まる。

列方向または行方向の内積を求めて(分散共分散行列になる)、

AT_A <- t(A) %*% A # column-wise 320row x 320col
A_AT <- A %*% t(A) # row-wise    200row x 200col    

トレース(対角成分の和)の平方根を求める。

> sqrt(sum(diag(AT_A)))
[1] 26607.34
> sqrt(sum(diag(A_AT)))
[1] 26607.34