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

# はじめに

テンソルを扱っていると、「あれ？いまどの足がどこにあるんだっけ？」とごちゃごちゃになることが多い。というわけで、[HOSVDによる画像の低ランク近似](http://qiita.com/kaityo256/items/0e8438c6c05897dadc2c)を題材に、一つ一つ確認しながらステップを進めてみる。

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

# ステップ・バイ・ステップ
## ステップ1　画像の読み込みまで

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

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

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

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

![uma.jpg](https://qiita-image-store.s3.amazonaws.com/0/79744/1fef2e44-69b6-92a3-0541-bbc63fd044cc.jpeg)

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

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

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

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

```pycon
>>> w, h = img.size
```

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## ステップ4 情報の復元

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

まずは次元の確認。

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

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

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

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

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

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

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

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

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

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

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

こうしてできた復元画像がこんな感じ。
![approximated.jpg](https://qiita-image-store.s3.amazonaws.com/0/79744/3c4b763d-5389-f323-77de-a60353c2d65d.jpeg)

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

# まとめ

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

しっかしPythonは便利だね・・・[^ruby]。


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

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