5.2. 数値計算#
NumPy には、基本的な算術演算に加えて、統計量の計算や固有値分解などの行列演算といった、数値解析に必要なさまざまな機能が用意されています。本節では、NumPy を使った代表的な数値計算の具体例を紹介していきます。
5.2.1. 算術計算#
配列同士の計算を行うときは、計算対象となる配列の形が同じである必要があります。形が異なる配列同士で計算をすると、ValueError が発生します。
x = np.array([1, 2, 3, 4])
y = np.array([1, 0, 1, 0])
x + y
array([2, 2, 4, 4])
x = np.array([1, 2, 3, 4])
y = np.array([1, 0])
x + y
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[3], line 4
1 x = np.array([1, 2, 3, 4])
2 y = np.array([1, 0])
----> 4 x + y
ValueError: operands could not be broadcast together with shapes (4,) (2,)
配列の計算対象が 1 つの値の場合、ブロードキャスト（broadcasting）の枠組みにより計算ができるようになります。
x = np.array([1, 2, 3, 4])
y = 1
x + y
array([2, 3, 4, 5])
二次元配列の場合も、基本的には 計算対象となる配列の形が同じ でなければなりません。形が一致していれば、各要素ごとの計算が行われます。
x = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
y = np.array([[1, 1, 1],
[0, 1, 1],
[0, 0, 1]])
x + y
array([[ 2, 3, 4],
[ 4, 6, 7],
[ 7, 8, 10]])
一方で、形が異なる配列同士では、ブロードキャストにより、自動的に形をそろえて計算ができる場合があります。例えば、次のように 行数が一致している場合には、列方向に値が繰り返されて計算されます。
x = np.array([[2, 4, 6, 8],
[1, 3, 4, 5]])
y = np.array([[1],
[0]])
x + y
array([[3, 5, 7, 9],
[1, 3, 4, 5]])
同様に、列数が一致している場合は、行方向に値が繰り返されて計算できます。
x = np.array([[2, 4, 6, 8],
[1, 3, 4, 5]])
y = np.array([1, 0, 1, 0])
x - y
array([[1, 4, 5, 8],
[0, 3, 3, 5]])
ブロードキャストは一定のルールに従って動作しますが、その仕組みはやや複雑で、最初は理解しづらいかもしれません。「最初は」というか、今の筆者は、理解しようとすることさえすでに諦めています。なので、実務では簡単な配列で試し計算しながら進めるのが確実です。
5.2.2. 要約統計量#
NumPy には、データの特徴を把握するためによく使われる要約統計量を計算する関数が用意されています。以下にいくつかの例を示します。なお、これらの機能は NumPy だけでなく、Python の標準ライブラリに含まれる statistics モジュールや、データ分析に特化した Pandas でも利用できます。同じ処理を複数の方法で書くとコードが読みづらくなるため、プログラム全体ではどれか一つのパッケージに統一して使うのが望ましいです。
x = np.array([1.49, 1.54, 1.53, 1.29, 1.40])
np.mean(x)
np.float64(1.45)
np.var((x))
np.float64(0.00884)
np.std(x)
np.float64(0.09402127418834527)
np.median(x)
np.float64(1.49)
配列内の最大値や最小値を調べることもでき、さらに、それらの値が配列のどの位置（インデックス）にあるかを知ることも可能です。
np.min(x)
np.float64(1.29)
np.max(x)
np.float64(1.54)
np.argmin(x)
np.int64(3)
np.argmax(x)
np.int64(1)
これらの関数は、デフォルトでは配列全体を対象にして集計を行います。そのため、多次元配列に対して最大値や最小値などを計算すると、出力は 1 つの値のみになります。
x = np.array([[4.83, 2.80, 1.56, 0.83],
[5.44, 2.54, 1.88, 0.97],
[5.42, 2.33, 1.76, 0.74]])
x.shape
(3, 4)
np.max(x)
np.float64(5.44)
配列全体ではなく、特定の次元（軸）ごとに集計したい場合は、axis オプションを指定します。たとえば、二次元配列に対して axis=0 を指定すると、1 次元方向（行をまたぐ方向）に向かって集計が行われ、列ごとの集計結果が得られます。
np.mean(x, axis=0)
array([5.23 , 2.55666667, 1.73333333, 0.84666667])
axis=1 を指定すると、2 次元方向（列をまたぐ方向）に向かって集計が行われ、行ごとの集計結果が得られます。
np.mean(x, axis=1)
array([2.505 , 2.7075, 2.5625])
5.2.3. 行列計算#
5.2.3.1. 基本計算#
NumPy では、二次元配列を数学の行列として扱うことができ、行列計算のためのさまざまな関数が用意されています。行列の形が同じであれば、加減乗除などの四則演算は、同じ位置にある要素同士で自動的に行われます。
x = np.array([[0, 1],
[1, 1]])
y = np.array([[1, 2],
[3, 4]])
x + y
array([[1, 3],
[4, 5]])
x - y
array([[-1, -1],
[-2, -3]])
二つの行列を * 演算子で掛け算すると、対応する位置の要素同士が掛け合わされます。この積のことをアダマール積（Hadamard product）と呼びます。計算結果は行列であり、入力行列と同じ形になります。
x * y
array([[0, 2],
[3, 4]])
x / y
array([[0. , 0.5 ],
[0.33333333, 0.25 ]])
5.2.3.2. 行列の積#
行列同士の積には複数の種類があり、それぞれ意味と計算方法が異なります。* 演算子で計算するアダマール積のほかに、内積（inner product）やクロス積、テンソル積などがあります。内積は、np.matmul で計算できます[1]。
x = np.array([[ 2, 3],
[ 1, 2]])
y = np.array([[ 2, -3],
[-1, 2]])
np.matmul(x, y)
array([[1, 0],
[0, 1]])
一方、クロス積（cross product）は、物理学などで回転や力の方向を求める際に使われる演算です。クロス積はベクトル積（vector product）とも呼ばれます。また、テンソル積（tensor product）は線形代数で扱われる数学的な概念です。日本語ではこれらの演算をまとめて外積と呼ぶこともありますが、意味や用途は異なります。Python では、クロス積は np.cross、テンソル積は np.outer を使って計算できます。
np.cross(x, y)
array([-12, 4])
np.outer(x, y)
array([[ 4, -6, -2, 4],
[ 6, -9, -3, 6],
[ 2, -3, -1, 2],
[ 4, -6, -2, 4]])
問題 NC-1
平面直行座標系 XY 上の点 \(P(3, 4)\) を、原点を中心に反時計回りに \(\frac{\pi}{2}\) だけ回転させたとき得られる点を \(P'\) とします。点 \(P'\) の座標を計算してください。なお、点 \((x, y)\) を原点を中心に反時計回りに \(\theta\) だけ回転させた時に得られる新しい座標は次のように計算できます。
5.2.3.3. 逆行列#
ある行列 \(\mathbf{A}\) に対して、\(\mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I}\) となるような行列 \(\mathbf{A}^{-1}\) が存在する時、\(\mathbf{A}^{-1}\) を \(\mathbf{A}\) の逆行列といいます。NumPy では np.linalg.inv 関数を利用して逆行列を計算できます。
x = np.array([[2, 3],
[1, 2]])
np.linalg.inv(x)
array([[ 2., -3.],
[-1., 2.]])
なお、逆行列はすべての行列に存在するわけではなく、たとえば次のような正則でない行列では逆行列が存在しません。正則でない行列に対して np.linalg.inv を使うと、LinAlgError が発生します。
x = np.array([[1, 2],
[2, 4]])
np.linalg.inv(x)
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
Cell In[28], line 3
1 x = np.array([[1, 2],
2 [2, 4]])
----> 3 np.linalg.inv(x)
File ~/.pyenv/versions/3.11.4/envs/ebook/lib/python3.11/site-packages/numpy/linalg/_linalg.py:669, in inv(a)
666 signature = 'D->D' if isComplexType(t) else 'd->d'
667 with errstate(call=_raise_linalgerror_singular, invalid='call',
668 over='ignore', divide='ignore', under='ignore'):
--> 669 ainv = _umath_linalg.inv(a, signature=signature)
670 return wrap(ainv.astype(result_t, copy=False))
File ~/.pyenv/versions/3.11.4/envs/ebook/lib/python3.11/site-packages/numpy/linalg/_linalg.py:163, in _raise_linalgerror_singular(err, flag)
162 def _raise_linalgerror_singular(err, flag):
--> 163 raise LinAlgError("Singular matrix")
LinAlgError: Singular matrix
逆行列が存在しない場合でも、その代わりに用いられるのがムーア・ペンローズ逆行列（Moore–Penrose inverse）です。NumPy では、np.linalg.pinv を使ってムーア・ペンローズ逆行列を計算することができます。
np.linalg.pinv(x)
array([[0.04, 0.08],
[0.08, 0.16]])
逆行列は、行列式（determinant）が 0 でないときに存在します。逆に、行列式が 0 のときには逆行列は存在しません。この行列式は np.linalg.det 関数で計算できます。
np.linalg.det(x)
np.float64(0.0)
問題 NC-2
入力された 2 行 2 列の行列（二次元配列）に対して、行列式が 0 のときにムーア・ペンローズ逆行列を計算し、そうでなければ逆行列を計算するようなプログラムを作成してください。
x = np.array([[ 2, 3],
[ 1, 2]])
行列を扱う上で、固有値と固有ベクトルの存在も重要です。固有値と固有ベクトルを計算するには np.linalg.eig 関数を利用します。
x = np.array([[5, 4],
[2, 1]])
eigen_values, eigen_vectors = np.linalg.eig(x)
eigen_values
array([ 6.46410162, -0.46410162])
eigen_vectors
array([[ 0.9390708 , -0.59069049],
[ 0.34372377, 0.80689822]])
試しに、ここで計算した固有値および固有ベクトルが、\(A \boldsymbol{x} = \lambda \boldsymbol{x}\) という関係（固有方程式）を満たしているかどうか、実際に計算して確認してみましょう。
eigen_values * eigen_vectors
array([[ 6.07024909, 0.27414041],
[ 2.22186537, -0.37448277]])
np.matmul(x, eigen_vectors)
array([[ 6.07024909, 0.27414041],
[ 2.22186537, -0.37448277]])