なんとなく需要がありそうなので書くことにしました。一般化逆行列とか疑似逆行列とか聞くけど具体的にどんなものなのか、どうやって計算するのか分からない、雰囲気を知りたい、という方向けの記事です。とても基本的なことしか書いていません。

方程式の数と変数の数が一致しない線形連立方程式

次の線形連立方程式を考えます．

\boldsymbol{Ax}=\boldsymbol{b}

ただし，$\boldsymbol{A}$は$m\times n$既知定数行列，$\boldsymbol{x}$は$n$次の未知ベクトル，$\boldsymbol{b}$は$m$次の既知定数ベクトルです．もし$\boldsymbol{A}$が正方行列（$n=m$）かつ正則ならば，解は

\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{b}

と一意に求まります．

工学上の諸問題において$n$と$m$は常に等しいとは限らず，$\boldsymbol{A}$も正則とは限りません．それでもなんらかの$\boldsymbol{x}$を求めたいという時がしばしばあります．一般化逆行列はそのような際に用いられるものです．

例えば$m<n$，つまり$\boldsymbol{A}$が横長である状況を考えましょう．ただし簡単のため$\boldsymbol{A}$は行フルランク（つまり$\boldsymbol{A}$の全ての行が線形独立）であるとします．このとき，元の線形連立方程式の解$\boldsymbol{x}$は無数に存在しますが，そのうちノルム（2乗ノルム）が最小となるものを求めてみましょう．これは次の最小化問題を解くことで得られます．

\frac{1}{2}\boldsymbol{x}^{\mathrm{T}}\boldsymbol{x} \rightarrow \mathrm{min.}
\qquad
\mathrm{subject~to}
\quad
\boldsymbol{Ax}-\boldsymbol{b}=\boldsymbol{0}

等式制約条件付きの凸計画問題なので，Lagrange乗数$\boldsymbol{\lambda}$を用いて

L\overset{\mathrm{def}}{=}
\frac{1}{2}\boldsymbol{x}^{\mathrm{T}}\boldsymbol{x}
-\boldsymbol{\lambda}^{\mathrm{T}}\left(\boldsymbol{Ax}-\boldsymbol{b}\right)

を定義すれば，$\boldsymbol{x}$は

\displaylines{
\left(\frac{\partial L}{\partial\boldsymbol{x}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{x}=\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\lambda}
\\
\left(\frac{\partial L}{\partial\boldsymbol{\lambda}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{Ax}=\boldsymbol{b}
}

の解です．上式を下式に代入すれば，

\boldsymbol{AA}^{\mathrm{T}}\boldsymbol{\lambda}=\boldsymbol{b}

となります．$\boldsymbol{A}$は行フルランクなので$\boldsymbol{AA}^{\mathrm{T}}$は正則な正方行列です．よって

\boldsymbol{\lambda}=\left(\boldsymbol{AA}^{\mathrm{T}}\right)^{-1}\boldsymbol{b}

より

\boldsymbol{x}=\boldsymbol{A}^{\mathrm{T}}\left(\boldsymbol{AA}^{\mathrm{T}}\right)^{-1}\boldsymbol{b}

を得ます．これをノルム最小型一般化逆行列解と呼びます．また，行列$\boldsymbol{A}^{+}=\boldsymbol{A}^{\mathrm{T}}\left(\boldsymbol{AA}^{\mathrm{T}}\right)^{-1}$をノルム最小型一般化逆行列と呼びます．

$\boldsymbol{A}^{+}$は次の性質を満たします．

\displaylines{
\boldsymbol{A}\boldsymbol{A}^{+}\boldsymbol{A}=\boldsymbol{A}
\\
\boldsymbol{A}^{+}\boldsymbol{A}\boldsymbol{A}^{+}=\boldsymbol{A}^{+}
}

こんどは$m>n$，つまり$\boldsymbol{A}$が縦長である状況を考えましょう．簡単のため$\boldsymbol{A}$は列フルランク（つまり$\boldsymbol{A}$の全ての列が線形独立）であるとします．このとき，元の線形連立方程式の解$\boldsymbol{x}$は一般的には存在しないので，次の最小二乗問題の解をもって代替することにしましょう．

E=\frac{1}{2}\left(\boldsymbol{Ax}-\boldsymbol{b}\right)^{\mathrm{T}}\left(\boldsymbol{Ax}-\boldsymbol{b}\right) \rightarrow \mathrm{min.}

これは無制約凸計画問題なので，$\boldsymbol{x}$は

\left(\frac{\partial E}{\partial\boldsymbol{x}}\right)^{\mathrm{T}}=\boldsymbol{0}
\quad\Leftrightarrow\quad
\boldsymbol{x}=\left(\boldsymbol{A}^{\mathrm{T}}\boldsymbol{A}\right)^{-1}\boldsymbol{A}^{\mathrm{T}}\boldsymbol{b}

と直ちに得られます．$\boldsymbol{A}$が列フルランクなので，$\boldsymbol{A}^{\mathrm{T}}\boldsymbol{A}$が正則であることにご注意下さい．これを最小二乗型一般化逆行列解と呼びます．また，行列$\boldsymbol{A}^{-}=\left(\boldsymbol{A}^{\mathrm{T}}\boldsymbol{A}\right)^{-1}\boldsymbol{A}^{\mathrm{T}}$を最小二乗型一般化逆行列と呼びます．

$\boldsymbol{A}^{-}$は次の性質を満たします．

\displaylines{
\boldsymbol{A}\boldsymbol{A}^{-}\boldsymbol{A}=\boldsymbol{A}
\\
\boldsymbol{A}^{-}\boldsymbol{A}\boldsymbol{A}^{-}=\boldsymbol{A}^{-}
}

以上のように，行列$\boldsymbol{A}$に対して

\displaylines{
\boldsymbol{A}\boldsymbol{A}^{\#}\boldsymbol{A}=\boldsymbol{A}
\\
\boldsymbol{A}^{\#}\boldsymbol{A}\boldsymbol{A}^{\#}=\boldsymbol{A}^{\#}
}

を満たす$\boldsymbol{A}^{\#}$を$\boldsymbol{A}$の一般化逆行列，またこれを用いて得た元の方程式の解$\boldsymbol{x}=\boldsymbol{A}^{\#}\boldsymbol{b}$を一般化逆行列解とそれぞれ呼びます．

Moore-Penroseの一般化逆行列

重要な事実として，任意の非零行列$\boldsymbol{A}$は必ずある列フルランク行列$\boldsymbol{B}$（$m\times r$行列）およびある行フルランク行列$\boldsymbol{C}^{\mathrm{T}}$（$r\times n$行列）の積として表せます．

\boldsymbol{A}=\boldsymbol{B}\boldsymbol{C}^{\mathrm{T}}

ただし，$r$は$\boldsymbol{A}$のランクであり，$r\leq m$かつ$r\leq n$を満たす自然数です．このとき元の方程式は，

\displaylines{
\boldsymbol{B}\boldsymbol{y}=\boldsymbol{b}
\\
\boldsymbol{C}^{\mathrm{T}}\boldsymbol{x}=\boldsymbol{y}
}

と2段階に分解できるので，まず一つめの方程式の最小二乗型一般化逆行列解が

\boldsymbol{y}=\left(\boldsymbol{B}^{\mathrm{T}}\boldsymbol{B}\right)^{-1}
\boldsymbol{B}^{\mathrm{T}}\boldsymbol{b}

と求まり，次いでこの$\boldsymbol{y}$を用いて二つめの方程式のノルム最小型一般化逆行列解が

\boldsymbol{x}=\boldsymbol{C}\left(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{C}\right)^{-1}\boldsymbol{y}

と求まります．これらをまとめた

\boldsymbol{x}=\boldsymbol{C}\left(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{C}\right)^{-1}
\left(\boldsymbol{B}^{\mathrm{T}}\boldsymbol{B}\right)^{-1}\boldsymbol{B}^{\mathrm{T}}\boldsymbol{b}

は，「元の方程式の両辺の誤差二乗を最小にするもののうちノルムが最小な$\boldsymbol{x}$」となっています．これをMoore-Penroseの一般化逆行列解と呼びます．また，行列$\boldsymbol{A}^{\dagger}=\boldsymbol{C}\left(\boldsymbol{C}^{\mathrm{T}}\boldsymbol{C}\right)^{-1}\left(\boldsymbol{B}^{\mathrm{T}}\boldsymbol{B}\right)^{-1}\boldsymbol{B}^{\mathrm{T}}$をMoore-Penroseの一般化逆行列もしくは疑似逆行列と呼びます．

念のため，$\boldsymbol{A}^{\dagger}$も次の性質を満たします．

\displaylines{
\boldsymbol{A}\boldsymbol{A}^{\dagger}\boldsymbol{A}=\boldsymbol{A}
\\
\boldsymbol{A}^{\dagger}\boldsymbol{A}\boldsymbol{A}^{\dagger}=\boldsymbol{A}^{\dagger}
}

Moore-Penroseの一般化逆行列の求め方

$\boldsymbol{A}$を$\boldsymbol{B}$と$\boldsymbol{C}^{\mathrm{T}}$に分解する方法は無数にあります．例えば特異値分解を用いれば，ある二つのユニタリ行列$\boldsymbol{U}$（$m\times r$行列）と$\boldsymbol{V}^{\mathrm{T}}$（$r\times n$行列），それにある対角行列$\boldsymbol{\Sigma}=\mathrm{diag}\{\sigma_{k}\}$（$r\times r$行列，全ての$k=1,\cdots,r$について$\sigma_{k}>0$）を用いて

\boldsymbol{A}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\mathrm{T}}

と表わすことができます．ユニタリ行列なので$\boldsymbol{U}^{\mathrm{T}}\boldsymbol{U}=\boldsymbol{1}$かつ$\boldsymbol{V}^{\mathrm{T}}\boldsymbol{V}=\boldsymbol{1}$（ただし$\boldsymbol{1}$は任意のサイズの単位行列）となることに注意すれば，Moore-Penroseの一般化逆行列は

\boldsymbol{A}^{\dagger}=
\boldsymbol{V}\left(\boldsymbol{V}^{\mathrm{T}}\boldsymbol{V}\right)^{-1}
\left(\boldsymbol{\Sigma}\boldsymbol{U}^{\mathrm{T}}\boldsymbol{U}\boldsymbol{\Sigma}\right)^{-1}\boldsymbol{\Sigma}\boldsymbol{U}^{\mathrm{T}}
=\boldsymbol{V}\boldsymbol{\Sigma}^{-1}\boldsymbol{U}^{\mathrm{T}}

となります（$\boldsymbol{B}=\boldsymbol{U}\boldsymbol{\Sigma}$、$\boldsymbol{C}^{\mathrm{T}}=\boldsymbol{V}^{\mathrm{T}}$とおいたということです）．

特異値分解はコストの大きい計算ですので，実用には適しません．よく使われるのはLQ分解（またはQR分解）です．これは$\boldsymbol{A}$を$n\times r$下三角行列$\boldsymbol{L}$および$r\times m$直交行列$\boldsymbol{Q}^{\mathrm{T}}$の積$\boldsymbol{A}=\boldsymbol{LQ}^{\mathrm{T}}$に分解するもので，Gram=Schmidtの直交化によって比較的高速に行えます（割愛します）．このときMoore-Penroseの一般化逆行列は

\boldsymbol{A}^{\dagger}=\boldsymbol{Q}(\boldsymbol{Q}^{\mathrm{T}}\boldsymbol{Q})^{-1}(\boldsymbol{L}^{\mathrm{T}}\boldsymbol{L})^{-1}\boldsymbol{L}^{\mathrm{T}}

ですが，$\boldsymbol{Q}^{\mathrm{T}}$は直交行列なので$\boldsymbol{Q}^{\mathrm{T}}\boldsymbol{Q}=\boldsymbol{1}$が成り立ちます．したがって，

\boldsymbol{A}^{\dagger}=\boldsymbol{Q}(\boldsymbol{L}^{\mathrm{T}}\boldsymbol{L})^{-1}\boldsymbol{L}^{\mathrm{T}}

と得られます．なお、LQ分解でなくQR分解$\boldsymbol{A}=\boldsymbol{QR}^{\mathrm{T}}$を用いた場合には

\boldsymbol{A}^{\dagger}=\boldsymbol{R}(\boldsymbol{R}^{\mathrm{T}}\boldsymbol{R})^{-1}\boldsymbol{Q}^{\mathrm{T}}

となります（導出は省きます）．

像空間と核空間

Moore-Penroseの一般化逆行列解は一意に存在しますが，$r<n$ならば元の方程式の解は無数に存在します．あらためて，元の方程式を2段階に分けた

\displaylines{
\boldsymbol{B}\boldsymbol{y}=\boldsymbol{b}
\\
\boldsymbol{C}^{\mathrm{T}}\boldsymbol{x}=\boldsymbol{y}
}

を考えましょう．これらのうち，無数の解を作り出しているのは後者です．より具体的には，$r<n$ならば

\boldsymbol{C}^{\mathrm{T}}\boldsymbol{x}=\boldsymbol{0}

が無数の非自明解$\bar{\boldsymbol{x}}$を持ちます．これらの$\boldsymbol{A}$による写像は（常に$\boldsymbol{0}$になるので）表に現れません．このような$\bar{\boldsymbol{x}}$の集合は$\boldsymbol{A}$の核空間（零空間）と呼ばれます．また，$\boldsymbol{A}$による写像は$\boldsymbol{B}$を基底とする空間の元となります．$\boldsymbol{B}$が張る空間は$\boldsymbol{A}$の像空間と呼ばれます．

もし$\boldsymbol{C}^{\mathrm{T}}$がユニタリ行列$\boldsymbol{V}^{\mathrm{T}}$ならば，核空間はその直交補空間$\boldsymbol{1}-\boldsymbol{V}\boldsymbol{V}^{\mathrm{T}}$と一致します．$\boldsymbol{A}$の特異値分解を利用すると、上述の通り

\boldsymbol{A}^{\dagger}=\boldsymbol{V}\boldsymbol{\Sigma}^{-1}\boldsymbol{U}^{\mathrm{T}}

ですから

\boldsymbol{A}^{\dagger}\boldsymbol{A}
=\boldsymbol{V}\boldsymbol{\Sigma}^{-1}\boldsymbol{U}^{\mathrm{T}}
\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\mathrm{T}}
=\boldsymbol{V}\boldsymbol{V}^{\mathrm{T}}

となります．したがって核空間の任意の元$\bar{\boldsymbol{x}}$は，任意の$n$次ベクトル$\boldsymbol{z}$を用いて

\bar{\boldsymbol{x}}=\left(\boldsymbol{1}-\boldsymbol{A}^{\dagger}\boldsymbol{A}\right)\boldsymbol{z}

と表すことができます．すなわち元の方程式の一般解は

\boldsymbol{x}=\boldsymbol{A}^{\dagger}\boldsymbol{b}+\left(\boldsymbol{1}-\boldsymbol{A}^{\dagger}\boldsymbol{A}\right)\boldsymbol{z}

となります．

なお，行列$\boldsymbol{1}-\boldsymbol{A}^{\dagger}\boldsymbol{A}$は，ランクが$n-r$の$n\times n$対称行列となりますので，

\boldsymbol{1}-\boldsymbol{A}^{\dagger}\boldsymbol{A}=\boldsymbol{D}\boldsymbol{D}^{\mathrm{T}}

を満たす$n\times(n-r)$行列$\boldsymbol{D}$が存在します．このような$\boldsymbol{D}$は，たとえばCholesky分解により求めることができます．これを用いれば一般解は

\boldsymbol{x}=\boldsymbol{A}^{\dagger}\boldsymbol{b}+\boldsymbol{D}\boldsymbol{v}

と書くこともできます．ただし$\boldsymbol{v}$は任意の$r$次ベクトルです．
$\boldsymbol{B}$の像空間と$\boldsymbol{D}$の像空間は、当然ながら互いに直交関係にあり、これらを合わせることで$n$次元の空間全体が張られます。

もし$\boldsymbol{A}$が列フルランクならば，$\boldsymbol{V}^{\mathrm{T}}$は正方行列となるので

\boldsymbol{1}-\boldsymbol{A}^{\dagger}\boldsymbol{A}\equiv\boldsymbol{O}

となります．つまりこのとき$\boldsymbol{A}$に核空間は存在しません．

ZMでの実装

ZMは、与えられた線形方程式のMoore-Penroseの一般化逆行列解を求める関数zLESolveMP()と、与えられた行列のMoore-Penroseの一般化逆行列を求める関数zMPInv()、与えられた行列のMoore-Penroseの一般化逆行列をもう一つ与えられた行列に左から掛ける関数zMulMPInvMatMat()を用意しています。方程式を解く計算は逆行列を求める計算よりも軽量にでき、かつ頻用されるので別関数にしています。コードは次のようになっています。

zLESolveMP()（関連する内部関数も含む）

/* weighted-norm-minimizing redundant linear equation solver
 * without checking size consistency. */
zVec zLESolveNormMinDST(zMat a, zVec b, zVec w, zVec ans, zMat m, zVec v, zIndex idx, zVec s)
{
  w ? zMatQuadNC( a, w, m ) : zMulMatMatTNC( a, a, m );
  if( !zLESolveGaussDST( m, b, v, idx, s ) ) return NULL;
  zMulMatTVecNC( a, v, ans );
  return w ? zVecAmpNCDRC( ans, w ) : ans;
}

/* error-minimizing inferior linear equation solver
 * without checking size consistency. */
zVec zLESolveErrorMinDST(zMat a, zVec b, zVec w, zVec ans, zMat m, zVec v, zIndex idx, zVec s)
{
  if( w ) zVecAmpNCDRC( b, w );
  zMulMatTVecNC( a, b, v );
  zMatTQuadNC( a, w, m );
  return zLESolveGaussDST( m, v, ans, idx, s );
}

/* allocate working memory for a lienar equation solver with Moore=Penrose's inverse matrix. */
static bool _zLEAllocWorkMP(zMat a, zVec b, zMat *l, zMat *u, zVec *bcp, zIndex *idx)
{
  *bcp = zVecClone( b );
  *l = zMatAllocSqr( zMatRowSizeNC(a) );
  *u = zMatAlloc( zMatRowSizeNC(a), zMatColSizeNC(a) );
  *idx = zIndexCreate( zMatRowSizeNC(a) );
  return *bcp && *l && *u && *idx ? true : false;
}

/* free working memory for a lienar equation solver with Moore=Penrose's inverse matrix. */
static void _zLEFreeWorkMP(zMat *l, zMat *u, zVec *bcp, zIndex *idx)
{
  zVecFree( *bcp );
  zMatFreeAO( 2, *l, *u );
  zIndexFree( *idx );
}

/* compute left-lower part of the linear equation. */
static void _zLESolveMP1(zMat l, zMat u, zVec b, zVec c, zVec we, zMat m, zVec v, zVec s, zIndex idx1, zIndex idx2, int rank)
{
  if( rank < zMatColSizeNC(l) ){
    zMatColReg( l, rank );
    zMatRowReg( u, rank );
    zLESolveErrorMinDST( l, b, we, c, m, v, idx2, s );
  } else
    zLESolveL( l, b, c, idx1 );
}

/* generalized linear equation solver using Moore-Penrose's
 * inverse (MP-inverse, pseudoinverse) based on LQ decomposition. */
zVec zLESolveMP(zMat a, zVec b, zVec wn, zVec we, zVec ans)
{
  int rank;
  zMat l, q, m;
  zVec bcp, c, v, s;
  zIndex idx1, idx2;

  if( !_zLEAllocWorkMP( a, b, &l, &q, &bcp, &idx1 ) ) goto TERMINATE2;
  if( ( rank = zLQDecomp( a, l, q, idx1 ) ) == 0 )
    goto TERMINATE2; /* extremely irregular case */
  c = zVecAlloc( rank );
  zLEAllocWork( &m, &v, &s, &idx2, rank );
  if( !c || !m || !v || !s || !idx2 ) goto TERMINATE1;

  _zLESolveMP1( l, q, bcp, c, we, m, v, s, idx1, idx2, rank );
  zMatIsSqr(q) ?
    zMulMatTVec( q, c, ans ) :
    zLESolveNormMinDST( q, c, wn, ans, m, v, idx2, s );

 TERMINATE1:
  zVecFree( c );
  zLEFreeWork( m, v, s, idx2 );
 TERMINATE2:
  _zLEFreeWorkMP( &l, &q, &bcp, &idx1 );
  return ans;
}

zMPInv()とzMulMPInvMatMat()（関連する内部関数も含む）

/* initialization: LQ decomposition */
static int _zMPInvAllocWork1(zMat m, zMat *l, zMat *q, zIndex *idx)
{
  int rank;

  if( ( rank = zLQDecompAlloc( m, l, q, idx ) ) == -1 ){
    zMatFreeAO( 2, *l, *q );
    zIndexFree( *idx );
    return -1;
  }
  return rank;
}

/* initialization: workspace for matrix computation */
static int _zMPInvAllocWork2(int rank, int row, zMat *tmp1, zMat *tmp2, zMat *tmp3)
{
  *tmp1 = zMatAlloc( rank, row );
  *tmp2 = zMatAllocSqr( rank );
  *tmp3 = zMatAlloc( rank, row );
  if( !*tmp1 || !*tmp2 || !*tmp3 ){
    zMatFreeAO( 3, *tmp1, *tmp2, *tmp3 );
    return -1;
  }
  return rank;
}

/* Moore-Penrose's inverse matrix. */
int zMPInv(zMat m, zMat mp)
{
  int rank;
  zMat l, q, tmp1, tmp2, tmp3;
  zIndex idx;

  if( ( rank = _zMPInvAllocWork1( m, &l, &q, &idx ) ) == -1 )
    return -1;
  if( _zMPInvAllocWork2( rank, zMatRowSizeNC(l), &tmp1, &tmp2, &tmp3 ) == -1 )
    return -1;

  if( zMatIsSqr(l) )
    zMatInv( l, tmp3 );
  else{
    zMatT( l, tmp1 );
    zMulMatTMat( l, l, tmp2 );
    zMulInvMatMat( tmp2, tmp1, tmp3 );
  }
  zMulMatTMat( q, tmp3, mp );

  zMatFreeAO( 5, l, q, tmp1, tmp2, tmp3 );
  zIndexFree( idx );
  return rank;
}

/* multiply Moore-Penrose's inverse matrix from the left side to another matrix. */
zMat zMulMPInvMatMat(zMat m1, zMat m2, zMat m)
{
  int rank;
  zMat l, q, tmp1, tmp2, tmp3;
  zIndex idx;

  if( ( rank = _zMPInvAllocWork1( m1, &l, &q, &idx ) ) == -1 )
    return NULL;
  if( _zMPInvAllocWork2( rank, zMatColSizeNC(m2), &tmp1, &tmp2, &tmp3 ) == -1 )
    return NULL;

  if( zMatIsSqr(l) )
    zMulInvMatMat( l, m2, tmp3 );
  else{
    zMulMatTMat( l, m2, tmp1 );
    zMulMatTMat( l, l, tmp2 );
    zMulInvMatMat( tmp2, tmp1, tmp3 );
  }
  zMulMatTMat( q, tmp3, m );

  zMatFreeAO( 5, l, q, tmp1, tmp2, tmp3 );
  zIndexFree( idx );
  return m;
}

ほか、核空間を用いた計算を行う関数も幾つか用意していますが、それらは割愛します。