HOSVDによる画像の低ランク近似 ステップ・バイ・ステップ

はじめに

テンソルを扱っていると、「あれ?いまどの足がどこにあるんだっけ?」とごちゃごちゃになることが多い。というわけで、HOSVDによる画像の低ランク近似を題材に、一つ一つ確認しながらステップを進めてみる。

以下、実習を想定し、過程をしつこいくらい詳細に書く。

ステップ・バイ・ステップ

ステップ1 画像の読み込みまで

以下、表示の都合上Pythonの対話モードで書いているが、実際にはipythonでやったほうがいろいろ便利。

まずは必要となるライブラリをインポートしよう。

>>> import numpy as np
>>> from scipy import linalg
>>> from PIL import Image

次に適当な画像を読み込む。ちょうど手元に手頃なサイズの画像があったので、それを使うことにしよう。

image0.jpeg

これはItanum2の精霊馬である。

これを読み込んで、サイズを確認してみよう。

>>> img = Image.open("uma.jpg")
>>> img.size
(480, 340)

480X340ピクセルの画像データだ。後のため、高さと幅を変数に落としておく。

>>> w, h = img.size

さて、イメージデータをNumpyの配列にしてみる。

>>> X = np.asarray(img)
>>> X.shape
(340, 480, 3)

順番が「高さ、幅、3(RGB)」となっていることに注意。

ステップ2 足を入れ替えてSVDの準備

これからSVDをするのだが、SVDをする足を一番右に持ってきて、さらに3階のテンソルを行列の形にする必要がある。足の入れ替えにはtransposeを使う。

再度Xの次元を確認する。

>>> X.shape
(340, 480, 3)

これはいま「高さ、幅、色」となっているが、まず幅について足を潰したいので、「高さ、色、幅」の順番にしなければならない。「高さ、幅、色」の順番が(0,1,2)なので、目標とする順番は(0,2,1)である。

>>> X.transpose(0,2,1).shape
(340, 3, 480)

順番が入れ替わった。これをさらに、左2つの足をまとめて行列にする。

>>> X.transpose(0,2,1).reshape(h*3,w).shape
(1020, 480)

無事に、(高さと色, 幅)の次元を持つ行列になった。これをX1として保存しておこう。

>>> X1 = X.transpose(0,2,1).reshape(h*3,w)

同様に、一番右に高さの次元を持つ行列を作りたい。もう一度Xの次元を見る。

>>> X.shape
(340, 480, 3)

これを(480, 3, 340)にしたい1。そのためには(1,2,0)の順番にすれば良い。

>>> X.transpose(1,2,0).shape
(480, 3, 340)

同様に、左の足を2本まとめて行列にして、これをX2に保存しておこう。

>>> X2 = X.transpose(1,2,0).reshape(w*3,h)
>>> X2.shape
(1440, 340)

右の足の次元が高さになっていますね。ここまででSVDの準備完了。

ステップ3 SVDしてコアテンソルを作る

まずX1をSVDする。さらにそれぞれのサイズを確認してみる。

>>> U, s, A1 = linalg.svd(X1)
>>> U.shape
(1020, 1020)
>>> s.shape
(480,)
>>> A1.shape
(480, 480)

さて、A1の情報を5%に圧縮する。もともと480行だったので、上位24行だけ取ってくる。これをa1としよう。

>>> a1 = A1[:24, :]
>>> a1.shape
(24, 480)

同様にX2をSVDして、A2を上位17行(340の5%)だけとってくる。

>>> U, s, A2 = linalg.svd(X2)
>>> A2.shape
(340, 340)
>>> a2 = A2[:17,:]
>>> a2.shape
(17, 340)

ここでえられたa1, a2の転置を取ると、幅、高さをつぶす圧縮行列になる。

>>> a1.T.shape
(480, 24)
>>> a2.T.shape
(340, 17)

それぞれ480→24、340→17に潰す行列となっていることがわかる。これをXの適当な足と内積を取ればコアテンソルを得る。

まず、Xにa1.Tをかける。どの足とどの足をつぶせば良いのかは、shapeを見ればわかる。

>>> X.shape
(340, 480, 3)
>>> a1.T.shape
(480, 24)

一番左を0番目として、Xの1番目とa1.Tの0番目の足を潰せば良い。これをX3に代入する。

>>> X3 = np.tensordot(X, a1.T, (1,0))
>>> X3.shape
(340, 3, 24)

無事に480→24に圧縮された。ここで注意すべきは、tensordot後の足の順番。慣れるまでは新しいテンソルのどの足が何になっているかを毎回shapeで確認したほうが良いと思う。

次に、X3の高さを表す足にa2.Tをかけて圧縮する。あらためてどれとどれをつぶすべきか調べる。

>>> X3.shape
(340, 3, 24)
>>> a2.T.shape
(340, 17)

足の大きさが全部違う場合は、同じ次元のところをつぶせばよいとすぐわかって便利。この場合は0番同士をつぶせば良いので、

>>> x = np.tensordot(X3,a2.T,(0,0))
>>> x.shape
(3, 24, 17)

無事に、Xの次元を削減したコアテンソルxが得られた。

最終的に、Xから、幅の圧縮/復元行列a1, 高さの圧縮/復元行列a2, コアテンソルxに圧縮された。この圧縮率を調べてみよう。

>>> float(x.size + a1.size + a2.size)/X.size
0.03783496732026144

ということで、もともとのデータの3.8%くらいまで圧縮できたことになる。

ステップ4 情報の復元

いま、手元にはコアテンソルxと、幅、高さの復元行列a1,a2がある。コアテンソルにa1,a2をかければ、元のテンソルの近似テンソルYが得られるので、それを作って見よう。

まずは次元の確認。

>>> x.shape
(3, 24, 17)
>>> a1.shape
(24, 480)

xの1番とa1の0番をつぶせばよいから、

>>> x2 = np.tensordot(x,a1,(1,0))
>>> x2.shape
(3, 17, 480)

幅が24から480に復元された。同様に高さを復元したものをx3とする。

>>> x3 = np.tensordot(x2,a2,(1,0))
>>> x3.shape
(3, 480, 340)

ただし、このままでは足の順番が違う。もう一度Xを見てみる。

>>> X.shape
(340, 480, 3)

Xと見比べれば、足を(2,1,0)の順番にすればよいことがわかるので、

>>> Y = x3.transpose(2,1,0)
>>> Y.shape
(340, 480, 3)

これが、コアテンソルxから、復元行列a1, a2を用いて復元されたテンソルである。 これで元のデータと同じ次元になったので、Imageオブジェクトを作ることができる。ただし、Image.fromarrayにつっこむ前にuint8型にしておく必要がある。

>>> img2 = Image.fromarray(np.uint8(Y))
>>> img2.save("approximated.jpg")

こうしてできた復元画像がこんな感じ。 image1.jpeg

まぁ、「元が画像である」という情報をまったく使わずに3.8%まで圧縮したデータから復元したわりにはそこそこ復元できたかな、という感じ。

まとめ

HOSVD(Higher Order Singular Value Decomposition)を使って画像を圧縮、復元してみた。IPython使うと、タブ補完とか効いてすごく便利。テンソルの足を潰したあと、何番目が何を表す足だったか混乱しがちだが、いちいちshapeで次元を調べながら作業すると間違い/混乱も減る気がする。

しっかしPythonは便利だね・・・2

  1. ここで、高さと色の足を入れ替えて(3,480,340)の順番にしないのは、まとめた足のデータが(R,G,B,R,G,B…)と並んでほしいから。(3,480,340)の順番だと、データが(RRR…,GGG….,BBB…)と、X1の時と異なる並び方になってしまう。 

  2. 一応筆者の常用言語はRubyなのだが、ライブラリの充実度はやっぱりPythonに一日の長がある感じ・・・。