行列計算 - Intro to Python

目次

LU分解を行う

scipyパッケージのlinalgサブモジュールのlu関数を使う。

n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = LUと分解することを、行列AのLU分解という。lu関数は置換行列Pを用いてA = P L Uと分解する。戻り値はLが下三角行列、Uが上三角行列、Pが置換行列。以下は4次の正方行列のパスカル行列をLU分解した例。最後に、元の行列に戻るかどうか計算している。

>>> import numpy as np
>>> from scipy.linalg import pascal
>>> from scipy.linalg import lu
>>> aa = pascal(4)
>>> print(aa)
[[ 1  1  1  1]
 [ 1  2  3  4]
 [ 1  3  6 10]
 [ 1  4 10 20]]
>>> pp, ll, uu = lu(aa)
>>> print(ll)
[[1.         0.         0.         0.        ]
 [1.         1.         0.         0.        ]
 [1.         0.33333333 1.         0.        ]
 [1.         0.66666667 1.         1.        ]]
>>> print(pp)
[[1. 0. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]]
>>> print(uu)
[[ 1.          1.          1.          1.        ]
 [ 0.          3.          9.         19.        ]
 [ 0.          0.         -1.         -3.33333333]
 [ 0.          0.          0.         -0.33333333]]
>>> print(pp @ ll @ uu)
[[ 1.  1.  1.  1.]
 [ 1.  2.  3.  4.]
 [ 1.  3.  6. 10.]
 [ 1.  4. 10. 20.]]

p_indicesオプションにTrueを指定すると、置換行列ではなく下三角行列Lの行を並び替えるためのベクトルが得られる。

>>> p, ll, uu = lu(aa, p_indices = True)
>>> print(p)
[0 2 3 1]
>>> print(ll[p, ] @ uu)
[[ 1.  1.  1.  1.]
 [ 1.  2.  3.  4.]
 [ 1.  3.  6. 10.]
 [ 1.  4. 10. 20.]]

コレスキー分解を行う

numpyパッケージのlinalg.cholesky関数を使う。コレスキー分解とは、既知の正定値対称行列Aに対してA = U^T Uを満たす上三角行列U（上三角コレスキー因子）、またはA = L L^Tを満たす下三角行列L（下三角コレスキー因子）を求めること。デフォルトでは下三角コレスキー因子が求まる。upperオプションにTrueを指定すると、上三角コレスキー因子が求まる。

>>> import numpy as np
>>> aa = np.array([1, 1, 1, 1, 2, 3, 1, 3, 6]).reshape(3, 3)
>>> print(aa)
[[1 1 1]
 [1 2 3]
 [1 3 6]]
>>> np.linalg.det(aa)
np.float64(1.0)
>>> uu = np.linalg.cholesky(aa, upper = True)
>>> print(uu)
[[1. 1. 1.]
 [0. 1. 2.]
 [0. 0. 1.]]
>>> print(uu.T.dot(uu))
[[1. 1. 1.]
 [1. 2. 3.]
 [1. 3. 6.]]
>>> ll = np.linalg.cholesky(aa)
>>> print(ll)
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 2. 1.]]
>>> print(ll.dot(ll.T))
[[1. 1. 1.]
 [1. 2. 3.]
 [1. 3. 6.]]

コレスキー分解を行う

scipyパッケージのlinalgサブパッケージのcholesky関数を使う。コレスキー分解とは、既知の正定値対称行列Aに対してA = U^T Uを満たす上三角行列U（上三角コレスキー因子）、またはA = L L^Tを満たす下三角行列L（下三角コレスキー因子）を求めること。下三角行列Lを求める場合は、lowerオプションをTrueにすること。

>>> import numpy as np
>>> from scipy.linalg import cholesky
>>> aa = np.array([1, 1, 1, 1, 2, 3, 1, 3, 6]).reshape(3, 3)
>>> print(aa)
[[1 1 1]
 [1 2 3]
 [1 3 6]]
>>> np.linalg.det(aa)
1.0
>>> uu = cholesky(aa)
>>> print(uu)
[[1. 1. 1.]
 [0. 1. 2.]
 [0. 0. 1.]]
>>> print(uu.T.dot(uu))
[[1. 1. 1.]
 [1. 2. 3.]
 [1. 3. 6.]]
>>> ll = cholesky(aa, lower = True)
>>> print(ll)
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 2. 1.]]
>>> print(ll.dot(ll.T))
[[1. 1. 1.]
 [1. 2. 3.]
 [1. 3. 6.]]

QR分解を行う

numpyパッケージのlinalgサブモジュールのqr関数を使う。

m行n列の行列A(m >= n)が与えられたとき、Q^T Q = Iを満たす行列Qと上三角行列（対角成分より下の成分はすべて0の行列）Rを用いてA = Q Rと分解することを、行列のQR分解という。QR分解はcomplete型とreduced型の2種類がある。それぞれ分解した行列の行数と列数が異なる

complete型 A(m×n) = Q(m×m) R(m×n)
reduced型 A(m×n) = Q(m×n) R(n×n)

以下は適当な行列を作成してQR分解した例。Q^T Q = Iを満たしていることと、Q R = Aと元に戻ることも計算している。modeオプションに'complete'を指定すると、complete型の計算を行う。'reduced'を指定するか何も指定しなければreduced型の計算結果が返される。

>>> import numpy as np
>>> import math
>>> m = 4
>>> n = 3
>>> aa = np.zeros((m, n))
>>> for i in range(m):
...     for j in range(n):
...         aa[i, j] = (-1) ** (i + 1) * math.comb(2 * (i + 1) + j + 1, i + 1)
...
>>> print(aa)  # m×n
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]
>>>
>>> # complete型 A(m×n) = Q(m×m) R(m×n)
>>> qq, rr = np.linalg.qr(aa, mode = 'complete')
>>> print(qq)
[[-0.02286814  0.33447694  0.78430749  0.52198083]
 [ 0.07622713 -0.54743741 -0.28540799  0.78297125]
 [-0.26679495  0.72431094 -0.53999463  0.33555911]
 [ 0.96046183  0.25260863 -0.10867309  0.0434984 ]]
>>> print(rr)
[[131.18688959 217.87238107 341.07829022]
 [  0.           2.93693148   9.35015978]
 [  0.           0.          -0.41767692]
 [  0.           0.           0.        ]]
>>> print(qq.T @ qq)
[[ 1.00000000e+00 -1.52611076e-16 -1.23191320e-16 -1.16499927e-16]
 [-1.52611076e-16  1.00000000e+00  1.17431915e-17  1.39258140e-17]
 [-1.23191320e-16  1.17431915e-17  1.00000000e+00 -4.75351482e-17]
 [-1.16499927e-16  1.39258140e-17 -4.75351482e-17  1.00000000e+00]]
>>> print(qq @ rr)
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]
>>>
>>> # reduced型 A(m×n) = Q(m×n) R(n×n)
>>> qq, rr = np.linalg.qr(aa, mode = 'reduced')
>>> print(qq)
[[-0.02286814  0.33447694  0.78430749]
 [ 0.07622713 -0.54743741 -0.28540799]
 [-0.26679495  0.72431094 -0.53999463]
 [ 0.96046183  0.25260863 -0.10867309]]
>>> print(rr)
[[131.18688959 217.87238107 341.07829022]
 [  0.           2.93693148   9.35015978]
 [  0.           0.          -0.41767692]]
>>> print(qq.T @ qq)
[[ 1.00000000e+00 -1.52611076e-16 -1.23191320e-16]
 [-1.52611076e-16  1.00000000e+00  1.17431915e-17]
 [-1.23191320e-16  1.17431915e-17  1.00000000e+00]]
>>> print(qq @ rr)
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]

QR分解を行う

scipyパッケージのlinalgサブモジュールのqr関数を使う。

m行n列の行列A(m >= n)が与えられたとき、Q^T Q = Iを満たす行列Qと上三角行列（対角成分より下の成分はすべて0の行列）Rを用いてA = Q Rと分解することを、行列のQR分解という。QR分解はcomplete型とreduced型の2種類がある。それぞれ分解した行列の行数と列数が異なる

complete型 A(m×n) = Q(m×m) R(m×n)
reduced型 A(m×n) = Q(m×n) R(n×n)

以下は適当な行列を作成してQR分解した例。Q^T Q = Iを満たしていることと、Q R = Aと元に戻ることも計算している。modeオプションに'full'を指定すると、complete型の計算を行う。'economic'を指定するか何も指定しなければreduced型の計算結果が返される。

>>> import numpy as np
>>> import math
>>> from scipy import linalg
>>> m = 4
>>> n = 3
>>> aa = np.zeros((m, n))
>>> for i in range(m):
...     for j in range(n):
...         aa[i, j] = (-1) ** (i + 1) * math.comb(2 * (i + 1) + j + 1, i + 1)
...
>>> print(aa)  # m×n
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]
>>>
>>> # complete型 A(m×n) = Q(m×m) R(m×n)
>>> qq, rr = linalg.qr(aa, mode = 'full')
>>> print(qq)
[[-0.02286814  0.33447694  0.78430749  0.52198083]
 [ 0.07622713 -0.54743741 -0.28540799  0.78297125]
 [-0.26679495  0.72431094 -0.53999463  0.33555911]
 [ 0.96046183  0.25260863 -0.10867309  0.0434984 ]]
>>> print(rr)
[[131.18688959 217.87238107 341.07829022]
 [  0.           2.93693148   9.35015978]
 [  0.           0.          -0.41767692]
 [  0.           0.           0.        ]]
>>> print(qq.T @ qq)
[[ 1.00000000e+00 -1.52611076e-16 -1.23191320e-16 -1.16499927e-16]
 [-1.52611076e-16  1.00000000e+00  1.17431915e-17  1.39258140e-17]
 [-1.23191320e-16  1.17431915e-17  1.00000000e+00 -4.75351482e-17]
 [-1.16499927e-16  1.39258140e-17 -4.75351482e-17  1.00000000e+00]]
>>> print(qq @ rr)
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]
>>>
>>> # reduced型 A(m×n) = Q(m×n) R(n×n)
>>> qq, rr = linalg.qr(aa, mode = 'economic')
>>> print(qq)
[[-0.02286814  0.33447694  0.78430749]
 [ 0.07622713 -0.54743741 -0.28540799]
 [-0.26679495  0.72431094 -0.53999463]
 [ 0.96046183  0.25260863 -0.10867309]]
>>> print(rr)
[[131.18688959 217.87238107 341.07829022]
 [  0.           2.93693148   9.35015978]
 [  0.           0.          -0.41767692]]
>>> print(qq.T @ qq)
[[ 1.00000000e+00 -1.52611076e-16 -1.23191320e-16]
 [-1.52611076e-16  1.00000000e+00  1.17431915e-17]
 [-1.23191320e-16  1.17431915e-17  1.00000000e+00]]
>>> print(qq @ rr)
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]

pivotingオプションにTrueを指定すると、ピボット選択が行われる。このピボットの情報は戻り値の3番目に含まれている。

>>> print(aa)
[[ -3.  -4.  -5.]
 [ 10.  15.  21.]
 [-35. -56. -84.]
 [126. 210. 330.]]
>>> qq, rr, p = linalg.qr(aa, mode = 'economic', pivoting = True)
>>> print(p)
[2 0 1]
>>> print(qq)
[[-0.01465387  0.29965789  0.79845252]
 [ 0.06154627 -0.53604535 -0.30955363]
 [-0.2461851   0.75472417 -0.50713352]
 [ 0.96715574  0.23076386 -0.0972919 ]]
>>> print(rr)
[[ 3.41206682e+02  1.31137526e+02  2.17870880e+02]
 [ 0.00000000e+00 -3.59852742e+00 -3.04345583e+00]
 [ 0.00000000e+00  0.00000000e+00  1.31063688e-01]]
>>> print(qq.T @ qq)
[[ 1.00000000e+00  1.81298658e-17 -3.43596805e-18]
 [ 1.81298658e-17  1.00000000e+00 -3.38842160e-17]
 [-3.43596805e-18 -3.38842160e-17  1.00000000e+00]]
>>> print(aa[:, p])
[[ -5.  -3.  -4.]
 [ 21.  10.  15.]
 [-84. -35. -56.]
 [330. 126. 210.]]
>>> print(qq @ rr)
[[ -5.  -3.  -4.]
 [ 21.  10.  15.]
 [-84. -35. -56.]
 [330. 126. 210.]]