白黒画像からそれらしいカラー画像を生成する技術Colorizationに関する論文を理解しようとすると、習ったはずの線形代数を復習するエッセンスがたっぷり含まれていることに気づいた。習ったはずだからと教科書に書かれている公式を駆使するのではなくて、なるべく原点（高校卒くらいの知識＋α）にもどって考えてみる。
Colorizationとは
Colorizationの問題は，グレースケール画像$Y(r) (r \in \mathbb{R}^2)$とユーザにより複数$(i=1 \dots s)$の画素$r_i$の色が$U(r_i)=U_i$, $V(r_i)=V_i$で与えられたとき、以下に示すエネルギー関数が最小となる$U(r),V(r)$の全画素値を見つけること。
\[
J(U) = \sum_r \left|U(r) - \sum_{s \in N(r)}w_{sr}U(s)\right|^2
\]
\[
w_{sr} = e^{-(Y(r) - Y(s))^2/2\sigma_r^2}
\]
ただし、$U(r)$と$V(r)$は全く対称なので、$U(r)$のみ代表して記載している。
また、$N(r)$は画素$r$の近傍画素の集合、$\sigma_r$は画素$r$の近傍画素値$\{Y(r)|r \in N(r)\}$の分散である。
$w_{sr}$は、色画像の近傍画素どうしの性質を拘束する重みで、輝度画像の$Y(r)$と$Y(s)$が近い値、つまりのっぺりとしたら$U(r)$と$U(s)$の値ものっぺりでなければならなく、$Y(r)$と$Y(s)$に大きな差があれば、つまりエッジの両側であれば、$U(r)$と$U(s)$もエッジで良い。
このエネルギー$J$が$\sum(測定値-予測値)^2$の形をしていたらLevenberg-Marquardt法が使えるので簡単だけど、この場合は$\sum(予測値-予測値 \times 測定値)^2$のような形をしている。どうやって最適化するのが正攻法なのか。
エネルギー関数をベクトルと行列（二次形式）で表す
論文中にも書かれているが$J=x^TAx$の形にして考えてみる。まず、全画素数を$n$とし、色画像$U(r)$を$n$次元ベクトル$u=[U(r_1),U(r_2),\dots,U(r_n)]^T$で表現するとエネルギー関数は次のようにベクトルの内積で表し、求めるパラメータである色画像の画素値$U$の要素だけを分離する。
\[
J(U) =
\left|
\left[
\begin{array}{c}
U(r_1) - (w_{11}U(r_1) + \cdots + w_{1?}U(r_?) + \cdots)\\
U(r_2) - (\cdots + w_{22}U(r_2) + w_{2?}U(r_?) + \cdots)\\
\vdots\\
U(r_n) - (\cdots + w_{n?}U(r_?) + \cdots + w_{nn}U(r_n))\\
\end{array}
\right]
\right|^2
\]
\[
=
{
\left|
\left[
\begin{array}{cccc}
1 - w_{11} & \cdots & - w_{1?} & \cdots\\
\cdots & 1 - w_{22} & - w_{2?} & \cdots\\
\vdots&& \ddots\\
\cdots & -w_{n?} & \cdots & 1 - w_{nn}\\
\end{array}
\right]
\left[
\begin{array}{c}
U(r_1)\\
U(r_2)\\
\vdots\\
U(r_n)
\end{array}
\right]
\right|^2
}
\]
\[
=
\left[
\begin{array}{cccc}
U(r_1) & U(r_2) & \cdots & U(r_n)
\end{array}
\right]
(I-W)^T (I-W)
\left[
\begin{array}{c}
U(r_1)\\
U(r_2)\\
\vdots\\
U(r_n)
\end{array}
\right]
= u^T A u
\]
ただし、記号?は近傍画素に対応するインデクス。
ここで$I$は単位行列、$W$は以下に示す任意の2組の画素の間の重みを網羅的に表す行列である。
\[
W =
\left[
\begin{array}{cccc}
w_{11} & w_{12} & \cdots & w_{1n}\\
w_{21} & w_{22} & & w_{2n}\\
\vdots & & \ddots & \vdots\\
w_{n1} & w_{n2} & \cdots & w_{nn}\\
\end{array}
\right]
\]
しかも、画素$r_i$が$r_j$の近傍でない場合、すなわち$r_i \notin N(r_j)$のとき、$w_{ij} = 0$となる。
要するに行列$W$は至る所の要素が0の行列である。また、$(I-W)^T (I-W)$の形から行列$A$は実対称行列である。
いずれにせよ、これでエネルギー関数$J(U)$が二次形式の形を取ることが示された。
\[
J(U)=u^TAu
\]
エネルギー関数の性質を考える
エネルギー関数が二次形式になっていることは分かった。二次形式は端的に言えば二次関数の多変数バージョン。ただ方向によって上向きに凸、下向きに凸、平坦な場合の組み合わせが存在するということ。行列$A$が正定値実対称行列であれば、この関数は下記のような形になって、図のように0を頂点にどの方向に向かっても上向きに反り上がる関数となる（図は$x^2+y^2+xy$）この関数は、全ての変数$U(r_i) = 0\ (i=1,\cdots,n)$のときのみ、最小値0をとる。
\[
a_{1}x_{1}^2 + a_{2}x_{1}^2 + \cdots + a_{n}x_{n}^2
\]
[Gnuplotコード表示／非表示]
reset
set xrange[-5:5]
set yrange[-5:5]
set zrange[0:40]
set pm3d
set ticslevel 0
set isosample 40
set hidden3d
set terminal png
set output 'gnuplot.png'
splot x**2+y**2+x*y
set terminal x11
splot x**2+y**2+x*y
エネルギー関数は確かに$x_i=0(i=1,\cdots,n)$で、もともと二乗和なので0以上である。ところがこの関数は、上のようなおわん型をしていない。おわん型で常に0以上なら$x_i=0(i=1,\cdots,n)$以外では0にならないはず。ところがこの関数に輝度一定($U(r_1) = U(r_2) = \cdots = U(r_n)$)の画像を入れると関数値は0になる。画素値$U(r)$と近傍画素値の重み付き平均$\sum_{s \in N(r)}w_{sr}U(s)$が等しくなってしまうから。
このことから少なくともエネルギー関数は、直線$U(r_1) = U(r_2) = \cdots = U(r_n)$の方向は一定値0をとる関数であることが分かる。これは一部の固有値が0をとることに等しい。この場合、下図のようにある方向には二次関数のように上向きに反り返り、ある方向には谷のようになだらかな形の関数になる。この図は、関数$x^2+y^2-2xy$のグラフで、$x$方向および$y$方向は上向きに反り上がっているが、$x=y$方向は0になっている。
[Gnuplotコード表示／非表示]
reset
set xrange[-5:5]
set yrange[-5:5]
set zrange[-20:40]
set pm3d
set ticslevel 0
set isosample 80
set hidden3d
set view 66,338
set terminal png
set output 'gnuplot.png'
splot (x-y)**2
set terminal x11
splot (x-y)**2
全ての画素の輝度値が0というのでも困るし、おわんの底が無いので最適な（エネルギー関数が最小となる）画素値が求まらないじゃないか。そこで次に、拘束条件としてユーザにより与えられた色$U(r_i)=U_i$をこの関数に導入することを考える。
拘束条件を加えたエネルギー関数
拘束条件は、ユーザにより与えられる複数$(i=1 \dots s)$の画素$r_i$の色$U(r_i)=U_i$, $V(r_i)=V_i$である。拘束条件が与えられれば、すぐに思いつくのはラグランジュの未定定数法だけど、この場合は非常に簡単な拘束条件なので、もとのエネルギー関数に代入する方法が良さそうだ。
\[
J(U) = u^T A u
=
\left[
\begin{array}{c}
U(r_{s+1})\\
U(r_{s+2})\\
\vdots\\
U(r_n)\\
U_1\\
\vdots\\
U_s
\end{array}
\right]^T
A
\left[
\begin{array}{c}
U(r_{s+1})\\
U(r_{s+2})\\
\vdots\\
U(r_n)\\
U_1\\
\vdots\\
U_s
\end{array}
\right]
\]
このような形になると変数$u$の場所に定数が含まれてしまって、この式の最小値を求めるなど一見難しそうに見える。しかし、全部掛け合わせてパラしてしまうと、以下のような形になるはず。
\[
J = [U(r_i)^2の項] + [U(r_i)U(r_j)の項] + [U(r_i)の項] + [定数の項]
\]
さらにこの式は、$U'(r_i) = f(U(r_1), \cdots, U(r_n))$などと変数変換すれば、以下のように二乗の項のみになる。
\[
J = [U'(r_i)^2の項] + [定数の項]
\]
これは$y=x^2 + 2x + 2$という二次式が与えられても、$x' = x + 1$と変数変換（一次変換）することで、$y=(x+1)^2+1=x'^2+1$と書きなおすことができて、$x = -1$のときに最小値1をとることが分かることと同じ。これの多変数バージョンだと考えればいい。多変数の場合も同じように、全ての変数の二乗の項＋定数となるように、適切に一次変換してくれる行列$L$を見つける問題となる。
\[
\left[
\begin{array}{c}
x'_1\\
\vdots\\
x'_n
\end{array}
\right] = L
\left[
\begin{array}{c}
x_1\\\
\vdots\\
x_n
\end{array}
\right]
\]
上で述べた、エネルギー関数を変数の二乗の項＋定数に変換できる、という根拠は行列の対角化可能性である。実対称行列は常に対角化可能だから。それがどう関係あるかというと、
もし、エネルギー関数を（変数変換により）変形し、下記の式のような形に出来るなら、二乗の項＋定数の形にできる。つまり行列Aを下記のように対角化できるなら変数変換だけで最適解$x_i = 0(i=1,\cdots,n-1)$と最小値$a_n c^2$が計算できることになる。
\[
\left[
\begin{array}{c}
x_1\\
x_2\\
\vdots\\
x_{n-1}\\
c
\end{array}
\right]^T
\left[
\begin{array}{cccc}
a_1 & 0 & & 0 & 0\\
0 & a_2 & & 0 & 0\\
& & \ddots& \vdots&\vdots\\
0 & 0 & \cdots& a_{n-1} & 0\\
0 & 0 & \cdots& 0 & a_n\\
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2\\
\vdots\\
x_{n-1}\\
c
\end{array}
\right]
=
\sum_{i=1}^{n-1} a_i x_i^2 + a_n c^2\ \ \ (a_i \geqq 0)
\]
実際のエネルギー関数を見ると制約条件が複数($i=1,\cdots,s$)あるので定数部分が複数ある。
複数あったとしても、下記のように行列のこの定数部分を１つにまとめることができる。
\[
J(U) = u^T A u =
\left[
\begin{array}{c|c}
変数 & 定数\\
\end{array}
\right]
\left[
\begin{array}{c|c}
2次式の係数 & 1次式の係数\\
\hline
1次式の係数 & 定数項\\
\end{array}
\right]
\left[
\begin{array}{c}
変数\\
\hline
定数\\
\end{array}
\right]
\]
\[
=
\left[
\begin{array}{cc}
{\bf x}^T & {\bf c}^T
\end{array}
\right]
\left[
\begin{array}{cc}
{\bf a^{(2)}} & {\bf a^{(1)}}\\
{\bf a^{(1)}}^T & {\bf a^{(0)}}
\end{array}
\right]
\left[
\begin{array}{c}
{\bf x}\\
{\bf c}
\end{array}
\right]
=
\left[
\begin{array}{cc}
{\bf x}^T & 1
\end{array}
\right]
\left[
\begin{array}{cc}
{\bf a^{(2)}} & {\bf a}^{(1)}{\bf c}\\
{\bf c}^T{\bf a}^{(1)T} & {\bf c}^T {\bf a}^{(2)} {\bf c}
\end{array}
\right]
\left[
\begin{array}{c}
{\bf x}\\
1
\end{array}
\right]
=
u'^TA'u'
\]
ただし、${\bf a^{(2)}}$は二次の係数、${\bf a}^{(1)}$は一次の係数、${\bf c}^T {\bf a}^{(2)} {\bf c}$は定数項である。これでようやく定数の要素が一箇所になった。
拘束条件付きの最適化
この関数を最小化する輝度$U(r_i)$を求めるには、既に述べたように$U(r_i)$を変数変換して、行列$A'$が対角化されるようにしなければならない。変数変換は一次変換でいいのでこの変換を$Lu'$とする。行列$A'$を対角化する行列$L$は、ひと通りではないけれど少なくとも、$L$の最終行、つまり$n-s+1$行目が対角成分以外全部0でないといけない。なぜなら$u'$の最後の要素は定数項に対応する1であり、$L$の最終行の要素に余計な成分が入っていると、線形変換により変数が混じってしまうことになる。これでは折角定数だけ下にもってきて1要素にまとめた意味がなくなってしまう。つまり、変換後も定数部分は定数のままでいるように保証されるような変換がよく、以下のような行列$L$で$u'$を$u''$に変換して対角化しなければならない。
\[
u'' = Lu' =
\left[
\begin{array}{cccc}
l_{1,1} & \cdots & l_{1,(n-s)} & l_{1,(n-s+1)}\\
\vdots & \ddots & \vdots & \vdots\\
l_{n,1} & \cdots & l_{(n-s),(n-s)} & l_{(n-s),(n-s+1)}\\
0 & \cdots & 0 & l_{(n-s+1),(n-s+1)}\\
\end{array}
\right]
u'
\]
この条件を満たす対角化方法が少なくとも1つあり、コレスキー分解である。ただし、行列$A$が正定値でないので行列$A'$も正定値でない可能性が高い。コレスキー分解は正定値限定なので、正定値の制約のない修正コレスキー分解を用いる。この分解では、$L$が上三角行列になるが、上の条件を満たし、以下のように元の変数を一次変換してくれる。
\[
J = u'^TA'u' = u'^T (L^T D L) u' = (Lu')^T D (Lu')
\]
\[
D=
\left[
\begin{array}{ccc}
d_1 & & \\
& \ddots & \\
& & d_{n-s+1}
\end{array}
\right]
\]
$D$は対角行列で行列$A'$が対角化されたということである。
一旦、コレスキー分解ができ$u''=[u''_1, \cdots, u''_{(n-s)}, u''_{(n-s+1)}]^T$と求まれば、最適解は$u''_i = 0 (i=1, \cdots, n-s)$でその時のエネルギー関数の値は$d_{(n-s+1),(n-s+1)}u''^2_{(n-s+1)}$である。このままでは一次変換したままの変数なので以下の式により元の画素値に逆変換する。
\[
\tilde{u} = L^{-1}u''
\]
参考
