いぬおさんのおもしろ数学実験室

おいしい紅茶でも飲みながら数学、物理、工学、プログラミング、そして読書を楽しみましょう

写真から立体を再現(11) 一般逆行列を使って求めた解はどのくらい正確か

 すみません、数学的に誤差を評価……ということではありません。一般逆行列を使って求めた「解」はある意味誤差が最小になりますが、実際に「解」を求めてそれがどんな感じなのか見てみよう、ということです。
 「ある意味」というのはこういうことです。x,yに何か値を入れると連立方程式の左辺の値が決まります。連立方程式なので3本の等式が成立しなければなりませんが、無理なこともあります。実際、下の例では3つの式を満たすx,yは存在しません。そこで次善の策です。各式でx,yを代入し、左辺との差をとって2乗します。それらを加えたものを残差(残差平方和)と呼びます。残差が最小になるようなx,yが見つかれば、それは本来の意味の解ではないですが「まあまあの解」と言ってもよいでしょう。
 次のような連立方程式を考えます。
f:id:Inuosann:20200418205453p:plain:w250
①、②の解はx=2、y=1で、これは③を満たしません。だからこの連立方程式の解はありません。でも、上で説明した意味の「解」は存在します。行列で表現すると次のようになります。
f:id:Inuosann:20200418205624p:plain:w250
これをAと表したとき、の左からAの一般逆行列というものをかけるとその「まあまあの解」が求まります。具体的にAの一般逆行列がどういう形なのか、それは前にブログで紹介しました。
www.omoshiro-suugaku.com
www.omoshiro-suugaku.com
ここではまたPythonのライブラリを使いましょう。

#=========================
#一般逆行列を試す
#
import numpy as np
import math
import sys
#行列をセットする
A = np.array([[1, 1], [1, -1], [0.9, -1.2]])
ginv = np.linalg.pinv(A) #Aの一般逆行列
print(ginv)
arr = np.dot(ginv, [3, 1, 0.8]) #一般逆行列を左からかける
print('\n', arr) 

print("\n計算された「解」をAにかけてみると……")
arr1 = np.dot(A, arr)
print(arr1)
sys.exit()

実行の結果は次の通り。
-----------------
[[ 0.53176471 0.27764706 0.21176471]
[ 0.45764706 -0.20352941 -0.28235294]]

[2.04235294 0.94352941]

計算された「解」をAにかけてみると……
[2.98588235 1.09882353 0.70588235]
-----------------
なかなかいい感じです! 見ると分かりますが、③は②を少しだけ変えた式です。だからなるべくよい解は本来の①、②の解にかなり近いのですね。同じようなことをすでにブログでやりましたが、今回は「おお、まあまあよい解だ!」と納得できそうな例を選びました。
 少し前からPythonでときどきプログラムを書いたりしています。行列の一部にベクトルを代入したり、転置行列の扱いなど、分かっている人には分かるんでしょうがぼくは使いづらい気もしています。それでも簡単にいろいろできるし、面白いと思います。