固有値問題

(3) 画像の固有空間

前回、カルーネン・レーベ展開によるデータの統計解析手法である「主成分分析」を紹介しました。直交基底による線形変換を利用した手法は画像の解析や変換などの処理にも応用されています。この章では、画像の基底とその線形変換がどのように利用されているかを紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 基底画像

画像は平面上に格子状に配置された「画素」と呼ばれる点の集まりによって表現されています。一つの画素が 1 か 0 のみで表現できるとき、つまり 1 ビットしかない場合は「二値画像」といい、その画像は二種類の色 ( 例えば白と黒 ) しか表現できません。現在、一般的に使われる画素は、赤・青・緑 ( これを光の三原色といいます ) の三種類の値を持っています。これらがそれぞれ 1 ビットの値を持つ場合、全部で 8 種類の色 ( 赤・青・緑・黄・シアン・マゼンタ・白・黒 ) を表現できます。さらに、三原色が 8 ビット ( 1 バイト ) の値を持てば、約 1677 種類の色を表現できるようになります ( これを「トゥルーカラー ( True Color ) といいます」)。明るさだけを値として持ち、0 を黒、最大値を白とするような表現も可能で、この場合は白黒画像になります。

画像は、ある値を持った点の集まりなので、これをベクトルとして表現することができます。画像が w x h の大きさの画素から構成されている場合、それは w x h 次元空間の一点とみなすことができるので、画像に対して内積とノルムを次のように定義することができます。

( f, g ) = Σy{1→h}( Σx{1→w}( fxygxy ) )

|| f ||2 = ( f, f )1/2 = Σy{1→h}( Σx{1→w}( fxy2 ) )

内積が 0 になる画像は互いに直交すると考えることができます。

m 個の画素を持つ任意の画像が、ある m 個の画像 { ek } ( k = 1, 2, ... m ) の線形結合で一通りに表せるとき、{ ek } は基底になり、各 ek を「基底画像」といいます。各基底画像が直交すれば、その基底は直交基底であり、さらに全基底画像のノルムが 1 であれば正規直交基底になります。

ここで、m 個の基底画像を使うのではなく、その一部である N 個の基底画像 ( N < m ) を使って画像を近似することを考えます。

f ≅ Σk{1→N}( ckek )

これを、画像 f の { ek } による「展開」といいます。「最小二乗法」に従えば、各係数 ck を次のようにすると、最もよく近似することができます。

ck = ( f, ek )

2) 固有空間法 ( Eigenspace Method )

画像に対して「カルーネン・レーベ展開 ( KLT )」を利用することを検討します。N 枚の画像 fi ( i = 1, 2, ... N ) があり、各画像の大きさは m であったとします。{ fi } の平均をとり、それを平均画像 g0 としたとき、

g0 = (1/N)Σi{1→N}( fi )

と表すことができます。また、各画像の平均画像からの差を { xi } とすると、

xi = fi - g0

{ xi } を使い、行列 M を以下のように計算します。

M = Σi{1→N}( xixiT )

行列 M は m x m の対称行列になります。この行列の固有値を λj、対応する単位固有ベクトルを gj ( 但し j = 1, 2, ... m ) とすると、xigj の線形結合として、次のように表すことができます。

xi = Σj{1→m}( cijgj )

但し、cij = ( xi, gj )

このように、各ベクトルの平均との差 { xi } ( i = 1, 2, ... N ) に対する和を使って対称行列 M を求め、その単位固有ベクトル { gj } ( j = 1, 2, ... m ) を正規直交基底として線形結合の形に表す展開方法を KLT というのでした。M の固有値が、大きい順に λ1 > λ2 > ... > λm と並べ替えられていたとすると、固有ベクトル g1 は { xi } を最もよく近似していることになり、二番めは g2、三番めは g3 という具合に近似性はだんだんと小さくなっていきます。M の各成分を N で割った行列は共分散行列と呼ばれ、その固有値はデータの主成分に対する分散を表しているので、これを利用してデータを分析することができ、それを「主成分分析」というのでした。{ gj } は言わば、画像の主成分であるとも言えます。

KLT によって展開された画像は、次のように表すことができます。

fi=g0 + xi
=g0 + Σj{1→m}( ( xi, gj )gj )
=g0 + ( xi, g1 )g1 + ( xi, g2 )g2 + ... + ( xi, gm )gm

固有値が小さい ( つまり分散の小さい ) 基底画像は、他の基底画像に比べて寄与率が小さいので無視することができます。よって、それを利用してデータを圧縮することができます。また、各係数はそれぞれの画像を表すパラメータであるため、画像認識などに応用することも可能です。このような利用法を総称して「固有空間法 ( Eigenspace Method )」といいます。

結局は KLT を行なうことと同じなので、主成分分析の場合と同様に処理を行なえばいいように見えますが、実際にはそのままで処理することはできません。画像を扱う場合は、さらに別の処理を検討する必要があります。


3) 特異値分解 ( Singular Value Decomposition )

問題になるのは、行列 M の大きさです。M の行列数 m は画像の大きさを表すので、例えば 1024 x 768 の大きさの画像を扱う場合、m = 786432 になります。これほどのサイズの行列に対して固有値と固有ベクトルを求めることは処理に時間がかかりすぎて不可能です。そこで、

M = Σi{1→n}( xixiT )

の固有値と固有ベクトルを別の形で計算することを検討します。

行列 M の r 行 c 列の成分 mrc は、xi = ( xi1, xi2, ... xim )T としたとき

mrc = Σi{1→n}( xirxic )

の形に表すことができます。ここで xi を列成分とする行列 P = ( x1, x2, ... xn ) を定義すると、これは m x n の大きさになります。このとき、PPT の r 行 c 列の成分 prc は、各 xi の第 r 番目の成分と、同じく各 xi の第 c 番目の成分の積について和を取ったものになるので、

prc = Σi{1→n}( xirxic )

になります。よって、mrc = prc であり、

M = PPT

が成り立ちます。ここで、N = PTP を定義すると、N の r 行 c 列の成分 nrc は、xrxc の内積そのものであり、

nrc = ( xr, xc ) = Σi{1→n}( xrixci )

になります。M は半正値対称行列でしたが、N も同様に半正値対称行列であることが証明できます。まず、

nrc = ( xr, xc ) = ( xc, xr ) = ncr

より、対称行列であることが示され、

( x, Nx ) = ( x, PTPx ) = ( Px, Px ) = ||Px||2 ≥ 0

ここで、( x, ATAx ) = ( Ax, Ax ) であることを利用しています ( 補足 1 )。従って、N も半正値対称行列になります。

M のゼロでない任意の固有値を λ ( > 0 ) とします。λ に対する固有ベクトルを u ( ≠ 0 ) としたとき、

Mu = PPTu = λu

が成り立ちます。左側からPTを掛けると、

PTPPTu = λPTu

PTP = Nより

N(PTu) = λ(PTu)

よって、λ は N の固有値でもあり、その固有ベクトルは PTu になります。これを逆に、N のゼロでない任意の固有値を λ ( > 0 )、λ に対する固有ベクトルを v ( ≠ 0 ) としたとき、

M(Pv) = λ(Pv)

となって、λ は M の固有値でもあり、その固有ベクトルは Pv になります。M の固有ベクトル u が単位ベクトルならば、

||PTu||2 = ( PTu, PTu ) = ( u, PPTu ) = ( u, Mu ) = ( u, λu ) = λ||u||2 = λ

なので、N の固有ベクトル v を単位ベクトルにするためには

v = ± PTu / √λ

とすれば求められます。逆に、M の固有ベクトル u を、v を使って単位ベクトルの形に表すためには

u = ± Pv / √λ

とすればよいことになります。固有ベクトルは符号を変えても成り立つので、符号は不定になります。

固有値がゼロになる場合、Mu = 0, Nv = 0 です。この時、

||PTu||2 = ( PTu, PTu ) = ( u, PPTu ) = ( u, Mu ) = ( u, 0 ) = 0

||Pv||2 = ( Pv, Pv ) = ( v, PTPv ) = ( v, Nv ) = ( v, 0 ) = 0

となり、PTu = Pv = 0 が成り立ちます。


以上の結果から、M のゼロでない固有値とその固有ベクトルを求めるには、M そのものを使う代わりに N を使ってもよいことになります。M のサイズは画像の大きさに依存しますが、N は扱う画像の枚数に依存します。主成分分析で扱ったデータは、データの種類に比べてデータ数の方が多かったのですが、今回はその逆になるわけです。そのようなときは、N を利用した方が有利です。扱う画像の枚数が 100 枚程度であったとすれば、1024 x 768 の大きさの画像を扱う場合、M のサイズが 786432 なのに対し、Nのサイズは 100 であり、固有値を求めることが十分可能になります。

行列 M = PPT の固有値を大きい順に並べ、それを λ1, λ2, ... λm とします。各固有値に対する単位固有ベクトルを u1, u2, ... um として、これらを列ベクトルとする m 行 m 列の行列 U を定義します。すなわち、

U = ( u1 u2 ... um )

になります。N = PTP についても同様に、固有値を大きい順に並べ、その単位固有ベクトルを v1, v2, ... vn として、これらを列ベクトルとする n 行 n 列の行列 V を定義します。すなわち、

V = ( v1 v2 ... vn )

です。それぞれの行列は、ゼロでない固有値の数が等しいので、その数を r とすると、i > r で λi = 0 が成り立ちます。

P は m x n の大きさの行列であり、n 行 n 列の行列 V との積を計算することができます。PV を計算すると、次のようになります。

PV=P( v1 v2 ... vn )
=( Pv1 Pv2 ... Pvn )

ゼロでない固有値に対する固有ベクトルについては、± Pv / √λ = u の関係が成り立つので、i ≤ r では Pvi = ± √λiui になります。固有値がゼロになる固有ベクトルに対しては、Pv = 0 になります。従って、ui の係数が全て正数となるように ui の符号を選べば、

PV=( √λ1u1 √λ2u2 ... √λrur 0 ... 0 )
=
( u1 u2 ... um )|√λ1,...|
|√λ2,...|
|:|
|...√λr, |
|:|
|...|

上の書き方ではわかりにくいですが、√λi を成分とする行列は m x n 行列であり、左上側には、r x r の大きさで対角成分上に √λ1から √λr まで並んだ対角行列があって、その他の成分は全てゼロになります。この行列を D とすれば、( u1 u2 ... um ) = Uなので

PV = UD

従って、

P = UDVT

が成り立つことになり、P を固有値分解と似たような形に分解することができます。D は、行列ともに r より後の成分が全てゼロなので、U の r 番目より右側の列ベクトルと、VT の r 番目より下側の行ベクトル ( つまり Vの r 番目より右側の列ベクトル ) がどんな値を持っていても結果は変化しません。よって、

Ur = ( u1 u2 ... ur )

Vr = ( v1 v2 ... vr )

とすれば、Ur は m x r、Vr は n x r の大きさの行列になり、これを使って

P=Ur|√λ1,|VrT = UrDrVrT
|√λ2,|
|:|
|√λr|

と変形することができます。但し、Dr は √λi を対角成分に持つ r x r の対角行列を表しています。このような分解手法を、m x n 行列 Pの「特異値分解(Singular Value Decomposition)」といいます。


KLT では、ベクトルの平均との差 { xi } を求めた上で、Σ( xixiT ) を計算して対称行列を求めていました。また、これは、xi を列ベクトルとする行列を P としたとき、PPT を計算しても同じことであることを先ほど示しました。さらに、PTP を計算して行列を求め、そこから固有値と固有ベクトルを求めても ( 固有ベクトルは変換をしなければなりませんが )、同じ結果が得られることが上記の内容から分かります。

データの種類が m 個あるのに対して、データは一つしかなかったと仮定します。この時、データは m 次元空間上の一点で表されます。よって、分布は、原点を含む一次元空間上 ( つまり直線上 ) のみに広がることになり、ゼロではない固有値はただ一つのみです。データ数が二つになると、分布は原点を含む二次元空間上 ( 平面上 ) に広がることになります。この時、ゼロでない固有値は二つになり、残りはゼロです。データ数が m 個になれば、分布は m 次元上に広がるので、全ての固有値がゼロではなくなります。もちろん、データは互いに線形独立であることが条件になります ( 例えば、データが全て直線上に並んでしまえば、一次元で表すことができてしまいます )。
KLT の場合、全データの平均が原点になります。データ数が一つなら、それは原点と一致するので、固有値は全てゼロです ( 共分散行列は零行列になります )。データ数が二つなら、原点を通る直線上にデータが並び、一次元上のみに分散が生じるので、ゼロでない固有値は一つのみです。データ数が m 個になっても、分布は、原点を含む m - 1 次元の部分空間上に存在することになって、ゼロの固有値が少なくとも一つ存在することになります。データ数が m 個を超えて初めて、全ての固有値がゼロでない可能性が生じるわけです。

画素数 m の画像をベクトルとして、N 枚の画像に KLT を行なった場合、空間の次元が m なのに対して点の数は N です。通常は m の方が圧倒的に大きいので、画像ベクトルからなる N 個の点は高々 N 次元の部分空間上に存在することになり、それ以上の次元に対しては分散がゼロ、すなわち固有値がゼロになるわけです。共分散行列として今まで使っていた行列 M のサイズが点の数より大きい場合、その差の分だけゼロの固有値があることが分かっています。よって、行列 N を使って固有値と固有ベクトルを求めた方が無駄なく処理できるというわけです。

実際に、次の 2 x 3 行列 PMN を求めてみます。

P =|1,2,3|
|3,2,1|

まずは MN を計算します。

M = PPT =|1,2,3||1,3|=|14,10|
|3,2,1||2,2||10,14|
|3,1|
N = PTP =|1,3||1,2,3|=|10,8,6|
|2,2||3,2,1||8,88|
|3,1||6,810|

それぞれの固有値を求めます。

M の固有方程式は

( λ - 14 )2 - (-10)2 = λ2 - 28λ + 96
= (λ - 4)(λ - 24) = 0

で、λ = 4, 24

N の固有方程式は

( λ - 8 )( λ - 10 )2 + (-8)・(-8)・(-6) + (-6)・(-8)・(-8) - (-6)・(-6)・( λ - 8 ) - (-8)・(-8)・( λ - 10 ) - (-8)・(-8)・( λ - 10 )
= λ3 - 28λ2 + 96λ
= λ(λ - 4)(λ - 24) = 0

で、λ = 0, 4, 24

固有ベクトルは、M の場合

λ = 4 のとき

|-10,-10||u11|=|0|
|-10,-10||u12||0| より u1 = ( -1/√2, 1/√2 )T

λ = 24 のとき

|10,-10||u21|=|0|
|-10,10||u22||0| より u2 = ( 1/√2, 1/√2 )T

N の場合

λ = 4 のとき

|-6,-8-6||v11|=|0|
|-8,-4-8||v12||0|
|-6,-8-6||v13||0| より v1 = ( 1/√2, 0, -1/√2 )T

λ = 24 のとき

|14,-8-6||v21|=|0|
|-8,16-8||v22||0|
|-6,-814||v23||0| より v2 = ( 1/√3, 1/√3, 1/√3 )T

λ = 0 のとき

|-10,-8-6||v31|=|0|
|-8,-8-8||v32||0|
|-6,-8-10||v33||0| より v3 = ( 1/√6, -2/√6, 1/√6 )T

このとき、MN の固有ベクトルの関係は次のようになります。

Pv1 / √4 =|1,2,3||1 / 2√2|=|-1/√2| = u1
|3,2,1||0||1/√2|
|-1 / 2√2|
Pv2 / √24 =|1,2,3||1 / 6√2|=|1/√2| = u2
|3,2,1||1 / 6√2||1/√2|
|1 / 6√2|
Pv3 =|1,2,3||1/√6|=|0| = 0
|3,2,1||-2/√6||0|
|1/√6|
PTu1 / √4 =|1,3||-1 / 2√2|=|1/√2| = v1
|2,2||1 / 2√2||0|
|3,1||-1/√2|
PTu2 / √24 =|1,3||1 / 4√3|=|1/√3| = v2
|2,2||1 / 4√3||1/√3|
|3,1||1/√3|

U = ( u1 u2 ), V = ( v1 v2 v3 ) とすれば、

PV=( Pv1, Pv2, Pv3 ) = ( √4u1, √24u2, 0 )
=
|u1,u2||√4,0,0|
|0,√24,0|
=UD

但し、

D =|√4,0,0|
|0,√24,0|

よって、行列 P は次のように分解することができます。

P=UDVT
=
|-1/√2,1/√2||√4,0,0||-1/√2,0,1/√2|
|1/√2,1/√2||0,√24,0||1/√3,1/√3,1/√3|
|1/√6,-2/√6,1/√6|

また、Ur = ( u1, u2 ), Vr = ( v1, v2 ) とし、D の第三列めの列ベクトルを除いた 2 x 2 の対角行列を Dr とすれば、

P=UrDrVrT
=
|-1/√2,1/√2||√4,||-1/√2,0,1/√2|
|1/√2,1/√2||√24||1/√3,1/√3,1/√3|

が成り立ちます。


特異値分解を行なうためのサンプル・プログラムを以下に示します。

/*
  CalcPPt : M = P(Pt) を求める
*/
void CalcPPt( const Matrix< stat_type >& p, SymmetricMatrix< stat_type >* m )
{
  typedef Matrix< ProbDist::stat_type >::size_type size_type;
  typedef Matrix< ProbDist::stat_type >::const_iterator const_iterator;

  m->assign( p.rows() );
  for ( size_type r = 0 ; r < p.rows() ; ++r ) {
    for ( size_type c = r ; c < p.rows() ; ++c ) {
      const_iterator row1 = p[r];
      const_iterator row2 = p[c];
      (*m)[r][c] = std::inner_product( row1, row1.end(), row2, ProbDist::stat_type() );
    }
  }
}

/*
  CalcPtP : N = (Pt)P を求める
*/
void CalcPtP( const Matrix< stat_type >& p, SymmetricMatrix< stat_type >* n )
{
  typedef Matrix< ProbDist::stat_type >::size_type size_type;
  typedef Matrix< ProbDist::stat_type >::const_iterator const_iterator;

  n->assign( p.cols() );
  for ( size_type r = 0 ; r < p.cols() ; ++r ) {
    for ( size_type c = r ; c < p.cols() ; ++c ) {
      const_iterator col1 = p.cbegin( r );
      const_iterator col2 = p.cbegin( c );
      (*n)[r][c] = std::inner_product( col1, col1.end(), col2, ProbDist::stat_type() );
    }
  }
}

/*
  SingularValueDecomposition : 特異値分解 ( P = UD(Vt) ) を行う

  p : 特異値分解を行う対象の行列 P
  u, d, v : 行列 U, D, V ( U, V は直交行列、D は対角行列 )
  e : 対角化を行うときに、非対角成分をゼロと判断するためのしきい値
  maxCnt : 対角化を行うときの反復処理最大回数
*/
void SingularValueDecomposition( const Matrix< stat_type >& p,
                                 Matrix< stat_type >* u, vector< stat_type >* d, Matrix< stat_type >* v,
                                 stat_type e, unsigned int maxCnt )
{
  SquareMatrix< stat_type > r;    // 上三角行列 R (対角成分が固有値になる)

  // 横長の行列
  if ( p.rows() <= p.cols() ) {
    SymmetricMatrix< stat_type > m; // M = P(Pt)
    SquareMatrix< stat_type > ev;    // M の固有ベクトル

    CalcPPt( p, &m ); // M = P(Pt) を求める

    // 対象行列の対角化 ( 対角行列の対角成分と直交行列 U を求める )
    Householder_DoubleShiftQR< stat_type >( m, &r, &ev, e, maxCnt );
    ev.clone( u );

    // szN : P の列数 : N = (Pt)P のサイズ
    Matrix< stat_type >::size_type szN = p.cols();
    // szM : P の行数 : M = P(Pt), R, U のサイズ
    Matrix< stat_type >::size_type szM = p.rows();

    // V = (Pt)U / √λ
    v->assign( szN, szM );
    for ( Matrix< stat_type >::size_type k = 0 ; k < szN ; ++k ) {
      for ( Matrix< stat_type >::size_type j = 0 ; j < szM ; ++j ) {
        Matrix< stat_type >::const_iterator ptRow = p.cbegin( k );
        Matrix< stat_type >::const_iterator uCol = u->cbegin( j );
        (*v)[k][j] = std::inner_product( ptRow, ptRow.end(), uCol, stat_type() );
        if ( r[j][j] != 0 ) (*v)[k][j] /= std::sqrt( std::abs( r[j][j] ) );
        if ( r[j][j] < 0 ) (*v)[k][j] *= -1;
      }
    }
  // 縦長の行列
  } else {
    SymmetricMatrix< stat_type > n; // N = (Pt)P
    SquareMatrix< stat_type > ev;    // N の固有ベクトル

    CalcPtP( p, &n ); // N = (Pt)P を求める

    // 対象行列の対角化 ( 対角行列の対角成分と直交行列 V を求める )
    Householder_DoubleShiftQR< stat_type >( n, &r, &ev, e, maxCnt );
    ev.clone( v );

    // szN : P の列数 : N = (Pt)P, R, V のサイズ
    Matrix< stat_type >::size_type szN = p.cols();
    // szM : P の行数 : M = P(Pt) のサイズ
    Matrix< stat_type >::size_type szM = p.rows();

    // U = PV / √λ
    u->assign( szM, szN );
    for ( Matrix< stat_type >::size_type k = 0 ; k < szM ; ++k ) {
      for ( Matrix< stat_type >::size_type j = 0 ; j < szN ; ++j ) {
        Matrix< stat_type >::const_iterator pRow = p[k];
        Matrix< stat_type >::const_iterator vCol = v->cbegin( j );
        (*u)[k][j] = std::inner_product( pRow, pRow.end(), vCol, stat_type() );
        if ( r[j][j] != 0 ) (*u)[k][j] /= std::sqrt( std::abs( r[j][j] ) );
        if ( r[j][j] < 0 ) (*u)[k][j] *= -1;
      }
    }
  }

  // 固有値の抽出
  d->assign( r.size(), stat_type() );
  for ( unsigned int i = 0 ; i < r.size() ; ++i )
    (*d)[i] = std::sqrt( std::abs( r[i][i] ) );
}

CalcPPt は行列 M = PPT を、CalcPtP は行列 N = PTP を求めるためのルーチンです。それぞれ、P の各行・列の組み合わせで内積を計算するだけの単純な処理ですが、行列クラス Matrix のメンバ関数である operator[]cbegin を利用して各行・列の先頭を指す反復子を得ています。operator[] が指定した行番号の行の先頭 ( 左端の要素 ) 位置を返すようになっているのは、反復子も operator[] を使って任意の位置の要素へアクセスできることを利用して、p[r][c] と表すことで r 行 c 列の要素へアクセスできるようにするためです。一般に、添字の一番目は行番号が用いられるので、p[c][r] ではなく p[r][c] と表せるようにしてあります。なお、画像などの場合、x, y 座標で表すときは x を第一引数とすることが多いので、例えば operator() では列番号を第一引数とする表し方 ( p( c, r ) の形式 ) を別に用意するという手もあります。

SingularValueDecomposition 関数が特異値分解のメイン・ルーチンとなります。まず行列 P の列と行のサイズを比較して場合分けを行います。もし列数の方が大きい ( 横長の行列の ) ときは、M = PPT を求め、その固有値と固有ベクトルを前章で紹介した Householder_DoubleShiftQR を使って得ます。

P = UDVT

より

D-1UTP = VT

が成り立つので

V = PTUD-1

となります。したがって、V の r 行 c 列目の要素は、PT の r 行目 ( すなわち P の r 列目 ) と U の c 列目の内積を D の c 行 c 列目の要素 ( 固有値の平方根 ) で割れば求められます。また、横長の行列のときは、このような処理をすればゼロの固有値は発生しなくなります。

縦長の行列の場合は、N = PTP を求めてその固有値と固有ベクトルを計算し、

U = PVD-1

を利用して U の要素を求めています。


4) 固有空間法の応用例

固有空間法を使った最も有名な応用例の一つに、顔認識技術があります。顔の画像をベクトル fi として、複数の顔の画像を元に平均顔画像 g0 と各平均差ベクトル xi = fi - g0 を求め、カルーネン・レーベ変換によって固有値と固有ベクトルを計算します。固有値の大きいものから順番に、その固有値に対する固有ベクトルと平均差ベクトルの内積 ( gj, xi ) を計算して係数 cij を求めると、各画像は平均画像と固有ベクトルによる線形結合で表すことができます。

fi = g0 + ci1g1 + ci2g2 + ...

各画像に対する係数は画像によって異なるので、それぞれの画像を認識するためのパラメータとして係数を利用することができます。こういった技術は「パターン認識」の一分野であり、音声の認識や文字認識などにも応用されています。パターン認識は、対象がどのパターン ( これをクラスといいます ) に含まれるかを識別する手法であり、クラスを分類するために利用されるパラメータを求めるために固有空間法が利用できるわけです。このパラメータは「特徴ベクトル ( Feature Vector )」、特徴ベクトルの張る空間は「特徴空間 ( Feature Space )」と呼ばれます。パターン認識とは、簡単に言うと、特徴空間を複数に分割してそれぞれをクラスとし、各特徴ベクトルがどのクラスに属するかを判定する手法であり、どのような特徴ベクトルを利用するかは非常に重要な要素になります。顔認識についても、画像全体をベクトルと見なすのではなく、例えば目や鼻・口の位置や、顎の輪郭など、個人差の大きなパラメータを選んだ方が認識率が向上するでしょう。現在の顔認識システムでは、誤認識が低く、また髪型などの変化に対してもきちんと認識できるような工夫がされているものと思います。

他にも、データ圧縮などに利用することができます。例えば、Aさんの顔は g0 + 25g1 + 30g2 - 5g3 までを計算すればかなり近似できるとすれば、g0 と、g1 から g3 までの画像と、それぞれの係数だけあれば A さんの顔を表現することができます。複数の顔画像を保存しなければならないとき、半分程度の基底画像で顔を認識できるだけの近似画像が得られるのであれば、基底画像の半数と係数だけを保存しておけば充分であり、保存量を小さくすることができます。

下の画像は、10 枚の顔画像を固有空間法によって基底画像に分解し、その線形結合によって近似画像を作成した結果を並べたものです。顔画像は、The Japanese Female Facial Expression (JAFFE) Database から拝借しました。10 枚の顔画像を使ったので固有値は 10 個になり、それに対する基底画像も 10 枚になります。その中の一名の顔写真に対して、全体の平均画像を一番左側、その右側に、近似画像の中から四枚分を並べたのが下の画像になります。

画像展開

下の画像は、元の顔画像です。上の画像と見比べると、第三近似画像あたりからほとんど変化がないことが分かると思います。

元の画像

下の表は、10 枚の画像を展開したときの係数の一覧です。No.7 と No.9 以外は最初の係数だけで大きな差があるので、この係数だけでほとんどの画像が判別可能です。二つ目の係数も含めれば、全ての画像が判別できることになります。

表 4-1. 基底画像に対する係数の算出結果
No.係数
1-463.2332535.54943.675-1936.901317.37318.3821037.16-91.5742668.533-136.080
2215.016-1229.081946.03970.026-205.25774.7955-850.624-927.7861495.93139.811
3-724.5305.6329150.5458-1067.77643.221316.307-1686.52-975.023-1259.98-186.390
4679.7322310.86-957.674683.606-1278.22-1301.22703.959-1228.99-326.301346.544
5306.088-1255.37-815.593918.9821783.64-1787.6523.0515576.16532.2521-538.359
6-413.485-1145.45-890.708-1678.39-890.357-643.506-574.2681081.00371.9311040.26
7-2286.92-1644.44-1229.95996.332970.6381474.92945.989-316.201-329.017900.482
81641.983147.87-3.463711545.48-214.505830.437-829.8471332.34-50.9558-65.1921
9-2222.54-1242.05-1923.4-488.942-1228.59823.590300.379-42.2748434.008-1273.45
103267.89-1483.532880.5357.5750-897.948-106.059930.712592.342-1036.40-227.629

今回のテストは画像の枚数も少なく、単純に画像全体を使って基底を求めているだけのものですが、それでも面白い結果が得られました。人の顔は通常左右対称とはならないそうですが、平均顔を見ると左右対称の画像になっています。非常に整った顔立ちに見えますが、個性がないような印象を受けるのは全体の平均を取っているからでしょうか。少ない枚数でも平均顔になるそうですが、もう少し枚数を増やすとこれが変化するかどうかは試してみる価値がありそうです。また、今回はサンプルが日本人女性の顔のみでしたが、これを多国籍にするとか、男性も加えてみると、また違った結果が得られるかもしれません。


補足1) ( x, ATAy ) = ( Ax, Ay ) の証明

任意の正方行列に対して ( x, ATy ) = ( Ax, y ) が成り立つことは「固有値問題 (1) 対称行列の固有値」の「3) 直交行列 ( Orthogonal Matrix )」で証明していましたが、行列 A が正方行列ではない場合、転置行列の行列数も反転するので、この式は成り立たなくなります ( その前に、内積自体が成り立ちません )。しかし、( x, ATAy ) と ( Ax, Ay ) は、ベクトルのサイズと行列の列数が等しければ内積の計算が可能であり、両式が等しくなることも証明できます。

行列 Aの c 列目の列ベクトルを ac とし、ac の r 番目の要素を arc とします。
ATA の r 行 c 列目の要素は ( ar, ac ) となるので、

( x, ATAy )=( x, ( Σi( ( a1, ai )yi ), Σi( ( a2, ai )yi ), ... ) )
=x1( Σi( ( a1, ai )yi ) ) + x2( Σi( ( a2, ai )yi ) ) + ...
=x1( ( a1, a1 )y1 + ( a1, a2 )y2 + ... ) + x2( ( a2, a1 )y1 + ( a2, a2 )y2 + ... ) + ...

また、

( Ax, Ay )=( ( Σi( a1ixi ), Σi( a2ixi ), ... ), ( Σi( a1iyi ), Σi( a2iyi ), ... ) )
=Σk( ( Σi( akixi ) )( Σj( akjyj ) ) )
=Σk( ( ak1x1 + ak2x2 + ... )( ak1y1 + ak2y2 + ... ) )
=Σk( x1( ak1ak1y1 + ak1ak2y2 + ... ) + x2( ak2ak1y1 + ak2ak2y2 + ... ) + ... )
=x1( Σk( ak1ak1 )y1 + Σk( ak1ak2 )y2 + ... ) + x2( Σk( ak2ak1 )y1 + Σk( ak2ak2 )y2 + ... ) + ...
=x1( ( a1, a1 )y1 + ( a1, a2 )y2 + ... ) + x2( ( a2, a1 )y1 + ( a2, a2 )y2 + ... ) + ...

よって、( x, ATAy ) = ( Ax, Ay ) になります。


<参考文献>
  1. 「これなら分かる応用数学教室」 金谷健一著 (共立出版)
  2. 「わかりやすい パターン認識」 石井健一郎他 共著 (オーム社)
  3. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system