SciPyによる線形代数(scipy.linalg)①｜SciPy入門 #1

f:id:lib-arts:20190225182115g:plain

SciPyについて色々と話題になったのでまとめていければと思います。

SciPy — SciPy v1.2.1 Reference Guide
まとめるにあたっては上記の公式チュートリアルが良さそうだったのでこちらをベースにまとめていきます。内容に関してはまずは線形代数(Linear Algebra)について取り扱えればということで、scipy.linalgについてまとめられればと思います。

Linear Algebra (scipy.linalg) — SciPy v1.2.1 Reference Guide
#1では上記の内容をまとめますが、1回で取り扱うには若干分量が多かったのでBasic routinesまでの内容としたいと思います。
以下目次になります。

1. numpy.matrix vs 2D numpy.ndarray
2. Basic routines
3. まとめ

1. numpy.matrix vs 2D numpy.ndarray
行列を扱う上で基本的に用いられるNumPyのnumpy.matrixやnumpy.ndarrayなどについて言及されています。numpy.matrixは便利である反面、下記のようにあまりよくない指摘もあがっています。

Despite its convenience, the use of the numpy.matrix class is discouraged, since it adds nothing that cannot be accomplished with 2D numpy.ndarray objects, and may lead to a confusion of which class is being used.

上記を受けて、ScyPyでの実装についてまとめられています。(numpy.matrixは省略します)

f:id:lib-arts:20190225183917p:plain
チュートリアルでは上図のような動作例についてまとまっています。>>>のついている行を貼り付ければ実際に実行することができます。計算の中身としては、単に行列を生成し、積を計算したり逆行列を計算したりしているだけです。

動作環境は今回の目的が動作確認のため、Jupyterを推奨します。

Jupyterの導入に関しては上記記事にまとめましたのでわからない方はこちらを確認いただけたらと思います。（環境構築に関してはいろいろなやり方がまとめられていますが、ライブラリの動作確認や理解が目的の際は基本的にコードが動けば良いです。そのため、別のやり方で構築されている場合は上記は気にしないでいただければと思います。）

2. Basic routines
・逆行列の計算(Finding Inverse)

f:id:lib-arts:20190225185219p:plain
上図のような計算をSciPyを用いて行います。コードは下記を貼り付ければ実行できます。

import numpy as np
from scipy import linalg
A = np.array([[1,3,5],[2,5,1],[2,3,8]])

print(A)
print(linalg.inv(A))
print(A.dot(linalg.inv(A))) #double check

上記のコードはチュートリアルのコードを動作確認用に編集したものです。ライブラリのインポートはすでに行っていれば不要です。機能の説明としては、linalg.invで逆行列を計算することができます。

・連立一次方程式を解く(Solving linear system)

f:id:lib-arts:20190225190516p:plain
linalg.invの出力結果を新しいベクトルにかけることで、上記のように連立一次方程式の解を求めることができます。
が、linalg.solveを利用する方が速度面や解の精度面で良いということでlinalg.solveが連立一次方程式を解くにあたっては推奨されています。下記のコードで実際に検証がされています。

import numpy as np #すでに行っていれば不要
from scipy import linalg #すでに行っていれば不要
A = np.array([[1, 2], [3, 4]])
print(A)
b = np.array([[5], [6]])
print(b)

print(linalg.inv(A).dot(b)) # slow
print(A.dot(linalg.inv(A).dot(b)) - b) # check1
print(np.linalg.solve(A, b)) # fast
print(A.dot(np.linalg.solve(A, b)) - b) # check2

上記の結果を見ると、check1はごく小さな小数の値を持つのに対し、check2では0になっています。このことが示すように、連立一次方程式を解く際はlinalg.solveを用いる方が良さそうです。ちなみにここでnp.linalg.solveが使われている理由はざっと見た感じだと見つからなかったですが、scipy.linalg.solveで実行しても同じ結果が出たので一旦はどちらでも良いという理解をしました。（後日必要があれば詳しく調べてみたいと思います。）

・行列式の計算(Finding Determinant)
行列式の計算にあたってはlinalg.detを用います。

import numpy as np #すでに行っていれば不要
from scipy import linalg #すでに行っていれば不要
A = np.array([[1,2],[3,4]])
print(A)
print(linalg.det(A))

上記のように引数として行列を与えるだけで行列式の計算が可能です。

・最小二乗法を解く(Solving linear least-squares problems)

f:id:lib-arts:20190225193208p:plain
最小二乗法においては上記の関数の最小化を行います。

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
c1, c2 = 5.0, 2.0
i = np.r_[1:11]
xi = 0.1*i
yi = c1*np.exp(-xi) + c2*xi
zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi))

A = np.c_[np.exp(-xi)[:, np.newaxis], xi[:, np.newaxis]]
c, resid, rank, sigma = linalg.lstsq(A, zi)
xi2 = np.r_[0.1:1.0:100j]
yi2 = c[0]*np.exp(-xi2) + c[1]*xi2

plt.plot(xi,zi,'x',xi2,yi2)
plt.axis([0,1.1,3.0,5.5])
plt.xlabel('$x_i$')
plt.title('Data fitting with linalg.lstsq')
plt.show()