OpenMV MicroPython
6.13. 線形代数¶
カメラ上での線形代数は小さな行列の作業です。IMU サンプルをワールドフレームに融合する 3x3 回転、レンズを補正するキャリブレーション行列、カルマンフィルタの状態共分散更新、正規方程式が小さな線形求解として現れる多項式フィッティングなどです。numpy.linalg と scipy.linalg のサブモジュールは、まさにこのスケールをカバーします。
これらの関数は2つのモジュールに存在します。行列積の演算は numpy のトップレベルにあり、分解と行列の逆行列は numpy.linalg の下に、専用の線形システムソルバーは scipy.linalg の下にあります。例えば 2x2 の行列乗算は次のようになります:
from ulab import numpy as np
A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
# [43.0, 50.0]])
6.13.1. 利用できるもの¶
dot()-- 行列積またはベクトル積。cross()-- 3次元ベクトルの外積。trace()-- 対角成分の和。inv()-- 行列の逆行列。det()-- 行列式。cholesky()-- コレスキー分解（対称正定値入力）。eig()-- 実対称行列の固有値と固有ベクトル。norm()-- ベクトルまたは行列の2ノルム。qr()--mode='reduced'（デフォルト）またはmode='complete'による QR 分解。solve_triangular()--Aが三角行列の場合にA @ x = bを解く。cho_solve()--Aのコレスキー因子が与えられた場合にA @ x = bを解く。
6.13.2. dot、cross、trace¶
dot() は、カメラ上で行列乗算を記述する方法です。デスクトップ版 numpy が同じ目的で提供する @ 演算子は ndarray には実装されていないため、行列・ベクトル積、行列・行列積、ベクトルのドット積はすべて dot() 呼び出しを通します:
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.dot(a, b) # 32.0 (scalar product)
np.cross(a, b) # array([-3.0, 6.0, -3.0])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a) # matrix-vector product
np.dot(m, m) # matrix-matrix product
np.trace(m) # 1 + 5 + 9 = 15.0
dot() の結果は常に dtype float です。cross() は2つの3要素ベクトルの外積、trace() は正方行列の主対角成分の和です。
6.13.3. inv と det¶
m = np.array([[1, 2, 3, 4],
[4, 5, 6, 4],
[7, 9, 9, 4],
[3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))
逆行列はガウス・ジョルダン消去法で計算されるため、行列が特異な場合（消去中に対角成分がゼロになる場合）には inv() は ValueError を送出します。RAM コストはおよそ入力サイズの2倍です。
行列式は同じ消去法を再利用します。実行時間は本質的に逆行列と同じです。
アプリケーションの目的が線形システムを解くことである場合は、逆行列を計算して乗算するのではなく、以下の専用ソルバーを使ってください。どちらも高速で、数値的にも振る舞いが良好です。
6.13.4. cholesky¶
対称正定値行列 A に対して、cholesky() は A = L @ L.T となる下三角行列 L を返します:
a = np.array([[25, 15, -5],
[15, 18, 0],
[-5, 0, 11]])
L = np.linalg.cholesky(a)
入力が正定値でないか、または対称でない場合は ValueError が送出されます。
コレスキー因子は LU 分解の半分の作業量であり、行列が対称正定値であることが分かっているあらゆる問題（共分散更新、最小二乗フィッティングからの正規方程式）に適した出発点です。
6.13.5. eig¶
eig() は実対称行列に対してのみ動作します。非対称行列は ValueError を送出します。(eigenvalues, eigenvectors) の2要素タプルを返します:
a = np.array([[1, 2, 1, 4],
[2, 5, 3, 5],
[1, 3, 6, 1],
[4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)
注意:
6.13.6. norm¶
ベクトルまたは行列のユークリッド（フロベニウス）ノルム:
v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.linalg.norm(v) # 7.416...
np.linalg.norm(m) # 16.881...
オプションの axis= キーワードは、配列全体ではなく単一の軸に沿ってノルムを取ります。
6.13.7. qr¶
qr() は、矩形行列 A（形状 (M, N)）を A == Q @ R となる正規直交行列 Q と上三角行列 R に分解します:
A = np.arange(6).reshape((3, 2))
q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)
q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)
この分解は連続するギブンス回転によって実装されています。行列が対称でない最小二乗問題に適した選択です。
6.13.8. システムを解く¶
ulab.scipy.linalg の下にある2つの専用ソルバーは、どちらも np.dot(np.linalg.inv(A), b) よりも高速で正確です。
solve_triangular(a, b, lower=False)()--aが三角行列であると仮定してa @ x = bを解く:A = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 2, 1, 8]]) b = np.array([4, 2, 4, 2]) x = sp.linalg.solve_triangular(A, b, lower=True)
cho_solve(L, b)()-- コレスキー因子Lが与えられた場合に、A = L @ L.TとなるA @ x = bを解く:L = np.linalg.cholesky(A) x = sp.linalg.cho_solve(L, b)
A の構造が許す限り、逆行列を計算するのではなくこれらを活用してください。これらは消去作業と明示的な逆行列計算の両方を節約します。
6.13.9. 小さな線形システムを解く¶
A = np.array([[3, 0, 1, 1],
[0, 1, 0, 2],
[1, 0, 1, 1],
[1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x)) # should equal b
同じ問題は、A を分解して適切なソルバーを呼び出すことでも表現できます。A が適切な構造を持つ場合、こちらの方が高速で正確です。
完全な引数レベルのリファレンスについては、numpy.linalg --- 線形代数ルーチン と scipy.linalg --- 線形代数ルーチン を参照してください。