画像圧縮〜少ないデータ量で近似する
０．問題
＜演習１＞[br]画像の濃淡を数値であらわすとき、識別力を下げずに次元を下げたい。[br]8☓8で「０」を表す画像データ行列を加工して、行列の階数を下げて圧縮の程度を確認しよう。[br]
１．演習１
「0」画像の濃淡を８行８列の行列に表したものをXとする。[br]この行列を３行列U、S、Vの積に分解するSVD分解をしよう。[br]すると、Sは対角成分に特異値s1,s2,s3,.....に、他は０になる。UとVが直交行列。[br]XはΣs[sub]i[/sub]U[sub]i[/sub]V[sub]i[/sub][sup]t[/sup]で近似できる。iはその行列のi列目を切り出した１列の列ベクトルの取り出しです。[br]X1＝s[sub]1[/sub]U[sub]1[/sub]V[sub]1[/sub][sup]t[br][/sup]X2＝s[sub]1[/sub]U[sub]1[/sub]V[sub]1[/sub][sup]t[/sup]+s[sub]2[/sub]U[sub]2[/sub]V[sub]2[/sub][sup]t[br][/sup]X3＝s[sub]1[/sub]U[sub]1[/sub]V[sub]1[/sub][sup]t[/sup]+s[sub]2[/sub]U[sub]2[/sub]V[sub]2[/sub][sup]t[/sup]+s[sub]3[/sub]U[sub]3[/sub]V[sub]3[/sub][sup]t[br]。。。。。。。。。[br][/sup]X６ ≒ X[sup][br][/sup]行列のサイズが８であっても、特異値s6まであり、S7=s8=0だとしたら、[br]X6を求めれば、もとのXとほぼ同じ画像に復元できる。[br][br]これを低ランク近似というデータ量がどのくらいへるかを計算してみる。[br]X1はs1で１個、U1列ベクトルに8個、V1列ベクトルに8個、合計1+8+8=17個[br]X2は２列目の特異値に対応するデータが追加されて、17＊2=34個[br]X3は３列目の特異値に対応するデータが追加されて、17*3=51個[br]X4は4列目の特異値に対応するデータが追加されて、17*4=68個[br].......[br]X6は17*6=102個[br]もとのデータ量は8*8=64個だから、３列目まではデータ節約になる。[br]もし800*600=48万個のデータでできた画像が、rank100に圧縮できれば、[br]100(1+800+600)≒14万個のデータで済むので、保持するデータ量の節約になる。[br]つまり、S,U,Vの行列データを端からすべて保持せず、端から100番目まであれば済むことになる。
[b][size=150]＜julia＞[/size][/b][color=#38761d][br][/color][color=#38761d]using [/color]LinearAlgebra[br]x0=[ 0. 0. 5. 13. 9. 1. 0. 0.;[br] 0. 0. 13. 15. 10. 15. 5. 0.;[br] 0. 3. 15. 2. 0. 11. 8. 0.;[br] 0. 4. 12. 0. 0. 8. 8. 0.;[br] 0. 5. 8. 0. 0. 9. 8. 0.;[br] 0. 4. 11. 0. 1. 12. 7. 0.;[br] 0. 2. 14. 5. 10. 12. 0. 0.;[br] 0. 0. 6. 13. 10. 0. 0. 0.][br]using PyPlot[br][color=#38761d]#もとの画像の表示[br][/color]plt.imshow(x0)[br]#画像の行列をU,S,Vに分解する。[br]U,S,V=svd(x0)[br][color=#1e84cc]#==================================================================[br]SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}[br]U factor:[br]8×8 Matrix{Float64}:[br] -0.227147 0.489927 0.275667 … -0.556251 0.0 0.531409[br] -0.53908 0.273281 0.0188344 0.251859 -0.12505 0.0302385[br] -0.387732 -0.318589 0.193466 -0.408655 -0.0777116 -0.31044[br] -0.302179 -0.323856 0.267671 0.557167 -0.134868 0.517763[br] -0.263389 -0.314291 0.26639 -0.177026 -0.43388 -0.277758[br] -0.335106 -0.329088 -0.0536914 … -0.100526 0.828133 0.058899[br] -0.423318 0.0821909 -0.815776 -0.0375244 -0.185262 0.00673247[br] -0.235119 0.514853 0.274298 0.326904 0.227499 -0.521129[br]singular values:[br]8-element Vector{Float64}:[br] 48.30784500260761[br] 24.95585263978704[br] 8.020753260580447[br] 6.029373007622524[br] 3.534479321858529[br] 0.6157632795726049[br] 0.0[br] 0.0[br]Vt factor:[br]8×8 Matrix{Float64}:[br] -0.0 -0.121635 -0.635846 -0.351655 … -0.547891 -0.262225 -0.0[br] -0.0 -0.199337 -0.182616 0.678604 -0.292419 -0.344252 -0.0[br] 0.0 0.141722 -0.0620057 0.466305 -0.400015 0.690527 0.0[br] -0.0 -0.289877 -0.598803 0.25419 0.659918 0.108522 -0.0[br] 0.0 -0.707353 0.443133 0.1381 -0.0853768 -0.203097 0.0[br] 0.0 -0.583958 -0.0597936 -0.33869 … -0.107305 0.53186 0.0[br] 0.0 0.0 0.0 0.0 0.0 0.0 1.0[br] 1.0 0.0 0.0 0.0 0.0 0.0 0.0[/color]
[color=#38761d]using[/color] LinearAlgebra[br][color=#38761d]using[/color] PyPlot[br][color=#38761d]#表示用の関数を作る。kを指定して合計するΣS[sub]k[/sub]U[sub]k[/sub]V[sub]k[/sub][sup]t[/sup][br]function [/color][b]disp[/b](k)[br] x=sum([s[k]*U[:,k]*transpose(V[:,k]) for k in 1:k])[br] print("rank",k)[br] plt.imshow(x)[br][color=#38761d]end[br][br]#１つ実行するたびに、１つの画像が表示されます。[br][/color][b]disp(1)[br]disp(2)[br]......[br]disp(6)[/b]
[b]＜Python＞[br][/b][color=#38761d]#juliaと同様にできる。[br][/color][color=#38761d]import[/color][b] numpy as np[br]x0=np.array([[ 0., 0., 5., 13., 9., 1., 0., 0.],[br] [ 0., 0., 13., 15., 10., 15., 5., 0.],[br] [ 0., 3., 15., 2., 0., 11., 8., 0.],[br] [ 0., 4., 12., 0., 0., 8., 8., 0.],[br] [ 0., 5., 8., 0., 0., 9., 8., 0.],[br] [ 0., 4., 11., 0., 1., 12., 7., 0.],[br] [ 0., 2., 14., 5., 10., 12., 0., 0.],[br] [ 0., 0., 6., 13., 10., 0., 0., 0.]])[br][/b][color=#38761d]import[/color][b] numpy as np[br][/b][color=#38761d]from[/color][b] numpy.linalg [/b][color=#38761d]import[/color][b] svd[br]U, S, Vt = svd(x0)[br]print(U)[br]print(S)[br]print(Vt)[br][/b][size=85][color=#38761d]#===============================================================[br][[-0.22714678 0.48992739 0.27566678 -0.15755125 -0.12550139 -0.55625083 0. 0.53140882][br] [-0.53908005 0.27328068 0.01883441 0.71916474 0.19280258 0.25185888 -0.12505044 0.03023848][br] [-0.38773218 -0.31858921 0.19346641 -0.20168304 0.63296791 -0.40865501 -0.07771161 -0.31043999][br] [-0.30217909 -0.32385608 0.26767054 -0.36448633 0.05103709 0.5571669 -0.1348682 0.51776255][br] [-0.26338915 -0.31429094 0.26639012 0.09414371 -0.67474504 -0.1770263 -0.43387994 -0.27775823][br] [-0.33510624 -0.32908792 -0.05369139 0.11924292 -0.250848 -0.10052563 0.8281329 0.05889901][br] [-0.42331758 0.08219089 -0.81577588 -0.31623618 -0.11303018 -0.03752439 -0.18526211 0.00673247][br] [-0.23511864 0.51485336 0.27429832 -0.40170401 -0.11332252 0.32690435 0.22749926 -0.52112878]][br][48.307845 24.95585264 8.02075326 6.02937301 3.53447932 0.61576328 0. 0. ][br][[-0. -0.12163488 -0.63584592 -0.35165517 -0.29714822 -0.54789077 -0.26222547 -0. ][br] [-0. -0.19933667 -0.18261557 0.67860378 0.51224488 -0.29241939 -0.34425199 -0. ][br] [ 0. 0.14172169 -0.06200574 0.46630482 -0.34898491 -0.40001459 0.69052728 0. ][br] [-0. -0.28987701 -0.59880272 0.25418992 -0.21336759 0.65991771 0.10852187 -0. ][br] [ 0. -0.70735326 0.44313298 0.13810032 -0.48546381 -0.08537679 -0.20309735 0. ][br] [ 0. -0.58395848 -0.05979355 -0.33869034 0.49630315 -0.10730494 0.53185986 0. ][br] [ 0. 0. 0. 0. 0. 0. 0. 1. ][br] [ 1. 0. 0. 0. 0. 0. 0. 0. ]][br][/color][color=#38761d][br]# 取り出したベクトルを１列ベクトルとして整形してから、転置して行列の積を計算する。[br][/color][size=150][color=#38761d]from[/color] matplotlib [color=#38761d]import[/color] pyplot[br][color=#38761d]def[/color] [b]disp[/b](k):[br] x=sum([S[k]*U[:,k].reshape(8,1)@((Vt[k].reshape(8,1)).transpose()) for k in range(k)])[br] print("rank",k)[br] fig, ax = pyplot.subplots()[br] ax.imshow(x, cmap = 'binary')[br][br][/size][/size][color=#38761d]#１つ実行するたびに、１つの画像が表示されます。[/color][br]disp(1)[br]disp(2)[br]....[br]disp(6)[br][br]
[b]＜geogebra＞[br][/b][b]x0={{ 0., 0., 5., 13., 9., 1., 0., 0.},[br] {0., 0., 13., 15., 10., 15., 5., 0.},[br] { 0., 3., 15., 2., 0., 11., 8., 0.},[br] { 0., 4., 12., 0., 0., 8., 8., 0.},[br] { 0., 5., 8., 0., 0., 9., 8., 0.},[br] { 0., 4., 11., 0., 1., 12., 7., 0.},[br] { 0., 2., 14., 5., 10., 12., 0., 0.},[br] { 0., 0., 6., 13., 10., 0., 0., 0.}}[br]Ms=SVD(A) [br][/b][color=#0000ff]#３つの行列まとめて返すから、１，２，３番と番号をつけて表示してみる。[/color][b][br]U=Ms(1)[br]S=Ms(2)[br]V=Ms(3)[br][/b][color=#38761d]# geogebraはベクトルと行列を区別している。行ベクトルは{}を２重にしないと行列と見ない。[br][/color][color=#38761d]# Sequence(Element(U(1),k),k,1,8) でUの１行ベクトルを８行１列行列に変換する。[br]# {V(1)}でVの１行ベクトルを１行８列行列に変換する。[br]# s1U1V1tだけでも以下のように長くなる。[br][/color]x1=S(1,1) Sequence(Element(U(1),k),k,1,8) * {V(1)}[br]x2=Sum(Sequence(S(j,j) Sequence({Element(U(j),k)},k,1,8) {V(j)},j,1,2))[br]x3=Sum(Sequence(S(j,j) Sequence({Element(U(j),k)},k,1,8) {V(j)},j,1,3))[br]x4=Sum(Sequence(S(j,j) Sequence({Element(U(j),k)},k,1,8) {V(j)},j,1,4))[br]x5=Sum(Sequence(S(j,j) Sequence({Element(U(j),k)},k,1,8) {V(j)},j,1,5))[br][b]x6=Sum(Sequence(S(j,j) Sequence({Element(U(j),k)},k,1,8) {V(j)},j,1,6))[br][color=#38761d]#このように近似を続けると、[br][/color][/b][size=200][b][br]x6 ≒ x0[/b][/size][size=200]となる。[/size]