社内向けのに用意した機械学習入門の一部抜粋です。
（単純な連立方程式の解法から、最新の機械学習で使われる手法を紹介しようと試みたものの、その趣旨があまり伝わらず...不評でしたが。）
$pinv$ 関数の答えを計算してみようというもの。

例題

次の連立方程式を一般化逆行列を用いて解け。（回答は下にあります。）

(1)  \ \ \ \left\{
\begin{array}{ll}
x=1   \\
2x=0 
\end{array}
\right.\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
(2)  \ \ \ \left\{
\begin{array}{ll}
x=1   \\
3x=0 
\end{array}
\right. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(3) \ \ \ \
\begin{array}{ll}
x+y=1   
\end{array}

(1),(2)は解の不能な連立方程式,(3)は不定な方程式になります。(1)と(2)は係数が違うだけですが係数の違いは何か意味があるでしょうか？

連立方程式の解き方

例えば、連立方程式が

 A \dot\  x=b    \ \ \ \ \ \ \ \ \ \ \ \ \ \             (1)

であれば、その解は

A^t \dot\ A \dot\ x= A^t \dot\ b

として、正方行列$(A^t\dot\ A)$を作る。もしこれに逆行列が存在する場合は、下記のように

x=(A^t \dot\ A)^{-1}\dot\ A^t\dot\ b \ \ \ \ \ \ \ \ \ \ \ \  (2)

として、$x$を求めることが可能になるが、もし$A^t\dot\ A$に逆行列が無い場合は、どうすればよいだろう？

この場合は考え方を改める必要がある。
もとの方程式を(1)、次のように書いてみよう。

 A \dot\ x-b=0

そして、この方程式を解く代わりに、左辺の2乗を取って、

  ||Ax-b||^2 =\epsilon^2 \ \ \ \ (4)

とおいて、右辺の$\epsilon^2$が0になるか、または最小になる$x$を方程式の解に採用するように解$x$の条件を緩める。
もし、$x$が(1)の厳密解の場合は$\epsilon^2=0$となり最小になる条件を完全に満たすので、この考え方を最小しても何ら問題はない。
（１）の厳密解が分からない場合は、試した範囲で最小の$\epsilon^2$になる$x$をとりあえずの解として妥協するわけである。

式(1)の厳密解が複数ある場合でも、(4)の$\epsilon$を最小にする解は唯一であることに注意しよう。このような近似解$x$を最小二乗解と呼ぶ。

これは一般化逆行列$A^{+}$を用いて

x=A^+\dot\ b

として求めることができる。
$A^{+}$の求め方は複数あるが、説明しやすいものとしては行列$A$の特異値分解を使う方法と、正則化を用いる方法である。

特異値分解

特異値分解とは、任意の行列Aについて直交行列$S,V$と,$i$番目の対角成分のみに$\lambda_i$を持つ特異値行列$D$を用いて

A=S\dot\ D(\lambda_i)\dot\ V^T  \\
st. \ \ \ S\cdot\ S^T=I \\ 
\ \ \ \ \ \ \ V\cdot\ V^T=I

のように表すことができることである。
$D(\lambda_i)$は正方行列ではないが、対角行列になるので、その逆行列は対角成分を逆数にして、転置したものを新たに定義する。

D^{-1}(\lambda_i)\equiv D^T(\lambda_i^{-1})

例えば、$A$が$5\times 3$でランク３をもつ行列ならば、特異値行列$D$は

\left(\begin{matrix}
\lambda_{1} & 0 & 0\\
0 & \lambda_{2} & 0\\
0 & 0 & \lambda_{3} \\
0 & 0 &  0 \\
0 & 0 & 0
\end{matrix}
\right)

となるとすると、$D^{-1}$は

\left(\begin{matrix}
\lambda_{1}^{-1} & 0 & 0　& 0 & 0\\
0 & \lambda_{2}^{-1} & 0 & 0 & 0\\
0 & 0 & \lambda_{3}^{-1} & 0 & 0\\
\end{matrix}
\right)

A^{+}=(S\dot\ D(\lambda_i) \dot\ V^T)^{-1} =V\dot\ D^{-1}(\lambda_i)\dot\ S^{T}

正則化による方法

式(2)で$(A^t \dot\ A)$の逆行列が無い場合の対処方法として、正則化を用いることができる。
$I$を単位行列、$\lambda >0$を正則化パラメータとして、行列の方程式に$(A^t \dot\ A　\lambda I)$
のように加える。するとこの行列は対角成分の$\lambda$のおかげで必ず逆行列を持つ正則行列になる。
すると、

x=(A^t \dot\ A+\lambda\dot I)^{-1}\dot\ A^t\dot\ b  \ \ \ \ \ \ \ \ \ \ \ \  (3)

となるが、ここで$\lim \lambda \rightarrow 0$とすれば、$x$の厳密解に収束する。
つまり、正則化は逆行列を持たない行列に対角成分に値$\lambda$を加え、正則行列して、逆行列を求めたあとで、
$\lambda$を0の極限に持ち込むことで$x$の解を求める方法である。

これは有力な方法で、$\lambda$を完全にゼロにはせずに、ゼロからわずかに近い値に留めることで、
厳密解からやや離れた解を求めることができる。この解は厳密解ではないが、局所厳密解
から離れたロバストな解を求めることができる。このため正則化は解の安定性（ロバスト性）を上げることができる。

では問題を解いてみよう。

例題

次の連立方程式を一般化逆行列を用いて解け。

(1)  \ \ \ \left\{
\begin{array}{ll}
x=1   \\
2x=0 
\end{array}
\right.\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
(2)  \ \ \ \left\{
\begin{array}{ll}
x=1   \\
3x=0 
\end{array}
\right. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(3) \ \ \ \
\begin{array}{ll}
x+y=1   
\end{array}

(1)の解き方

行列で書くと次のようになります。

\left(
\begin{matrix}
1\\
2 
\end{matrix}
\right)\cdot x=
\left(
\begin{matrix}
1\\
0 
\end{matrix}
\right)\ \ \ \ \ \ \ \ \ \ \ (1)

左辺について特異値分解を$A=S\cdot D \cdot V^T$に対応させて書くと次のようになる。

\left(
\begin{matrix}
1\\
2 
\end{matrix}
\right)=
\left(
\begin{matrix}
\frac{1}{\sqrt{5}} & -\frac{2}{\sqrt{5}} \\
\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}}
\end{matrix}
\right)\cdot
\left(
\begin{matrix}
\sqrt{5} \\ 0
\end{matrix}
\right)\cdot 1

となります。
左辺のランクは１でフロベニウスノルムは$\sqrt{5}$ですから、特異値行列$D$の成分が分かります。行列$S$は$2\times 2$の直交行列なので容易に見つけられるはずです。
したがって一般逆行列は$A^{+}=V\cdot D^{-1}\cdot S^T$で計算できるから、

\left(
\begin{matrix}
1\\
2 
\end{matrix}
\right)^{+}=
1 \cdot
\left(
\begin{matrix}
\frac{1}{\sqrt{5}} & 0
\end{matrix}
\right)\cdot \left(
\begin{matrix}
\frac{1}{\sqrt{5}} & \frac{2}{\sqrt{5}} \\
-\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}}
\end{matrix}
\right)

よって、

\left(
\begin{matrix}
1\\
2 
\end{matrix}
\right)^{+}=\left(\frac{1}{5},\frac{2}{5}\right)

となる。pinv実際疑似逆行列は以下になる。

import numpy as np

np.linalg.pinv([[1],[2]])
#return: array([[0.2, 0.4]])

(1)式より、$x$は、

x=A^{+}\cdot \left(\begin{matrix}1 \\ 0\end{matrix}\right)

よって、

x=\frac{1}{5}

が求まります。

問題(2)の解き方

問題(1)と係数が違うだけで、問題(2)と同じですね。
にもかからわず、求まる解が異なることに注意してください。

前問と同様に行列で書くと

\left(
\begin{matrix}
1\\
3 
\end{matrix}
\right)\cdot x=
\left(
\begin{matrix}
1\\
0 
\end{matrix}
\right)\ \ \ \ \ \ \ \ \ \ \ (2)

左辺について特異値分解を$A=S\cdot D \cdot V^T$に対応させて書くと次のようになる。

\left(
\begin{matrix}
1\\
3 
\end{matrix}
\right)=
\left(
\begin{matrix}
\frac{1}{\sqrt{10}} & -\frac{3}{\sqrt{10}} \\
\frac{3}{\sqrt{10}} & \frac{1}{\sqrt{10}}
\end{matrix}
\right)\cdot
\left(
\begin{matrix}
\sqrt{10} \\ 0
\end{matrix}
\right)\cdot 1

となるので、

\left(
\begin{matrix}
1\\
3 
\end{matrix}
\right)^{+}=
1 \cdot
\left(
\begin{matrix}
\frac{1}{\sqrt{10}} & 0
\end{matrix}
\right)\cdot \left(
\begin{matrix}
\frac{1}{\sqrt{10}} & \frac{3}{\sqrt{10}} \\
-\frac{3}{\sqrt{10}} & \frac{1}{\sqrt{10}}
\end{matrix}
\right)=\left(\frac{1}{10},\frac{3}{10}\right)

となる。よって、

x=\left(\frac{1}{10},\frac{3}{10}\right)\cdot\ \left(\begin{matrix}1\\ 0\end{matrix}\right)=\frac{1}{10}

import numpy as np

np.linalg.pinv([[1],[3]])
# array([[0.1, 0.3]])

となります。
問題(2)は問題(1)の第二の両辺を3/2しただけだが異なる解が得られることに注意しよう。このような現象が起こることは記憶にとどめておいて損はないはず。
$x=0$という情報の重みが（1）と(2)で3/2倍異なって解釈しているために、求まる解が異なるのである。

問題(3)の解き方

問題(1)と同様に行列で書くと次のようになります。

\left( \begin{matrix}1 &1\end{matrix}\right)
\cdot \left(
\begin{matrix}
x\\
y 
\end{matrix}
\right)=1 \ \ \ \ \ \ \ \ \ \ \ (3)

左辺について特異値分解を$A=S\cdot D \cdot V^T$に対応させて書くと次のようになる。

\left(
\begin{matrix}
1&
1 
\end{matrix}
\right)=
1\cdot
\left(
\begin{matrix}
\sqrt{2} & 0
\end{matrix}
\right)\cdot 
\left(
\begin{matrix}
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\
\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}
\end{matrix}
\right)

となるので、

\left(
\begin{matrix}
1 &
1 
\end{matrix}
\right)^{+}=
 \left(
\begin{matrix}
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\
\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}
\end{matrix}
\right)
\cdot
\left(
\begin{matrix}
\frac{1}{\sqrt{2}} \\
 0
\end{matrix}
\right)\cdot 1=\left(\begin{matrix}\frac{1}{2}\\ \frac{1}{2}\end{matrix}\right)

となる。よって、

\left(\begin{matrix}x \\y\end{matrix}\right)=\left(\begin{matrix}\frac{1}{2}\\ \frac{1}{2}\end{matrix}\right)

となる。xy平面内の直線$x+y=1$内で、$x,y$のL2ノルムつまり原点からの２乗距離が最小なものが解になるのが確認できる。

import numpy as np

np.linalg.pinv([[1,1],[0 ,0]])  # 方程式が足りないので、第二式として0x+0y=0 を用いている。。
# return: array([[0.5, 0. ],
#        [0.5, 0. ]])

正則化による方法$\lambda \rightarrow 0$でも同様の解が得られる。試してください。

まとめ

上記例題から明らかになるが、手計算では、(A)と(B)はともに(C)に等しい。
しかし、実は実際に $pinv$で計算すると明らかなのだが、(A)の最小二乗解と(B)の最小二乗解と(C)の最小二乗解は異なることに注意しよう。

更に、(A)の第二式の係数２と(B)の係数３００は何を意味するだろうか？これはその式の重みを表している。例えば(B)の場合、$$x+y+z=1$$ よりも、第二式の$$x+y+z=0$$ の方に３００倍の重みがつく。（式(4)において、第二式の解の誤差が３００倍になる為）。

(A)  \ \ \ \left\{
\begin{array}{ll}
x+y+z=1   \\
2(x+y+z)=0 
\end{array}
\right.\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
(B)  \ \ \ \left\{
\begin{array}{ll}
x+y+z=1   \\
300(x+y+z)=0 
\end{array}
\right. \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(C) \ \ \ \
\begin{array}{ll}
x+y+z=1   \\
x+y+z=0 
\end{array}