線形代数の練習問題をPythonでやってみた

「データサイエンスのための数学」の練習問題を解いていて、手計算がとにかく苦手でミスしまくりなのでPythonで確認しながらやってみました。

この本には練習問題の解説がついてないので、独自の解説っぽいものも書いていきます。

※この記事は随時追記していきます。

練習問題 1.2 基本変形で逆行列を求める

問: 行列の基本変形によって以下の行列の逆行列をもとめよ

$$ \begin {pmatrix}
1 & 2 \\
2 & -1 \\
\end {pmatrix} $$

解法: $ \textbf{A}\textbf{A}^{-1} = \textbf{E} $より、$\textbf{A}$を$\textbf{E}$に変形する基本変形と同じ基本変形を$\textbf{E}$にすれば、それは$\textbf{A}^{-1}$になる。($\textbf{E}$は単位行列)

pythonで掃き出し法をやっていきます。一番上の行から、順に対角成分を1にして同じ列を0にする、という処理を行うのが効率が良いとよびのり先生が言っておりました。なのでその方法でやっていきます。

a = np.array([[1,2],[2,-1]], dtype=float)
e = np.identity(2)
combine = [a,e]

for mat in combine:
  print(mat)
[[ 1.  2.]
 [ 2. -1.]]
[[1. 0.]
 [0. 1.]]

とりあえずnp.linalg.invでaの逆行列を見てみます

np.linalg.inv(a)
array([[ 0.2,  0.4],
       [ 0.4, -0.2]])

これが正解です。では掃き出し法をやっていきます。

# (2,1)成分を0にする
for mat in combine:
  mat[1] = mat[1] + mat[0]*(-2)
  print(mat)
[[ 1.  2.]
 [ 0. -5.]]
[[ 1.  0.]
 [-2.  1.]]
# (2,2)成分を1にする
for mat in combine:
  mat[1] = mat[1] * (-1/5)
  print(mat)
[[ 1.  2.]
 [-0.  1.]]
[[ 1.   0. ]
 [ 0.4 -0.2]]
# (1,2)成分を0にする
for mat in combine:
  mat[0] = mat[0] - mat[1]*2
  print(mat)
[[ 1.  0.]
 [-0.  1.]]
[[ 0.2  0.4]
 [ 0.4 -0.2]]

ということで逆行列を掃き出し法で求めることができました。

ちなみに、

5*e
array([[ 1.,  2.],
       [ 2., -1.]])

となるので、

$ \frac{1}{5}
\begin {pmatrix}
1 & 2 \\
2 & -1 \\
\end {pmatrix} $

とスカラーとの積で綺麗にできます。

練習問題1.3 逆行列で連立方程式を解く

問: 逆行列を求めることで以下の連立方程式を解け

$$ (1) \begin{equation}
\left\{ \,
\begin{aligned}
& 2x + 3y = 350 \\
& 5x + 6y = 800
\end{aligned}
\right.
\end{equation}$$

解法: この連立方程式を$ \textbf{A}x = b $の形にすると、

$$ \begin {pmatrix}
2 & 3 \\
5 & 6 \\
\end {pmatrix}
\begin {pmatrix}
x \\
y \\
\end {pmatrix}
=
\begin {pmatrix}
350 \\
800 \\
\end {pmatrix} $$

また、

$ \textbf{A}^{-1} \textbf{A} \textbf{x} = \textbf{A}^{-1} \textbf{b} $
$ \textbf{E} \textbf{x} = \textbf{A}^{-1} \textbf{b} $
$ \textbf{x} = \textbf{A}^{-1} \textbf{b} $

なので、$ \textbf{A} $の逆行列と$ \textbf{b} $の積が$ \textbf{x} $になることがわかる。

ではPythonで逆行列を求めていく。

a = np.array([[2,3],[5,6]], dtype=float)
e = np.identity(2)
combine = [a,e]

print(np.linalg.inv(a))
# out:
# array([[-2.        ,  1.        ],
#        [ 1.66666667, -0.66666667]])

for mat in combine:
  print(mat)
# out:
# [[2. 3.]
#  [5. 6.]]
# [[1. 0.]
#  [0. 1.]]

for mat in combine:
  mat[0] = mat[0] * (1/2)  # (1,1)成分を1にする
  mat[1] = mat[1] - mat[0]*5  # (2,1)成分を0にする
  mat[1] = mat[1]*(-2/3)  # (2,2)成分を1にする
  mat[0] = mat[0] + mat[1]*(-1.5)  # (1,2)成分を0にする
  print(mat)
# out:
# [[ 1.  0.]
#  [-0.  1.]]
# [[-2.          1.        ]
#  [ 1.66666667 -0.66666667]]

これで逆行列を求めることができた。

$ \textbf{x} = \textbf{A}^{-1} \textbf{b} $より、

x = e.dot(np.array([350, 800]))
print(x)
# out:
# array([100.,  50.])

ということで、$x = 100, y = 50$が解である。

では次の問題

$$ (2) \begin{equation}
\left\{ \,
\begin{aligned}
& x + 2y – z = 20 \\
& 2x – y + 2z = 110 \\
& -x – y + z = -10
\end{aligned}
\right.
\end{equation}$$

a = np.array([[1,2,-1],[2,-1,2],[-1,-1,1]], dtype=float)
e = np.identity(3)
combine = [a,e]

for mat in combine:
  mat[1] = mat[1] + mat[2]*2
  mat[2] = mat[2] + mat[0]

  mat[1] = mat[1] * (-1/3)
  mat[0] = mat[0] + mat[1]*(-2)
  mat[2] = mat[2] + mat[1]*(-1)

  mat[2] = mat[2]*(3/4)
  mat[0] = mat[0] + mat[2]*(-5/3)
  mat[1] = mat[1] + mat[2]*(4/3)

  print(mat)
# out: 
# [[ 1.00000000e+00  0.00000000e+00 -2.22044605e-16]
#  [ 0.00000000e+00  1.00000000e+00  0.00000000e+00]
#  [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
# [[-0.25  0.25 -0.75]
#  [ 1.    0.    1.  ]
#  [ 0.75  0.25  1.25]]

print(4 * e)
# out: 
# array([[-1.,  1., -3.],
#        [ 4.,  0.,  4.],
#        [ 3.,  1.,  5.]])

print(e.dot(np.array([20,110,-10])))
# out:
# array([30., 10., 30.])

よって逆行列は

$ \frac{1}{4}
\begin {pmatrix}
-1 & 1 & -3 \\
4 & 0 & 4 \\
3 & 1 & 5
\end {pmatrix} $

であり、この連立方程式の解は

$ x=30, y=10,z=30 $となる

練習問題1.11 左右から行列を掛けて標準化

問: $p$個の観測項目についての$n$人のデータである$n\times p$行列$\textbf{X}$に対し、$\textbf{A}\textbf{X}\textbf{B}$のように左右から掛けることで行列$\textbf{X}$を標準化する行列$\textbf{A}, \textbf{B}$を求めよ。ただし、各観測項目の平均値と分散はそれぞれ$ (\bar{x_1}, \sigma_1^2), … ,(\bar{x_n}, \sigma_n^2) $とする。

解法: 標準化は$ \frac{x_{ij} – \mu_j}{\sigma_j} $で表されるように、それぞれのデータから各項目ごとの平均値を引いて、それを各項目ごとの標準偏差で割れば良い。

よって、$\textbf{A},\textbf{B}$を以下のような行列にする。

解:$\textbf{A} = \textbf{E}_n – \frac{1}{n} \textbf{1}_n $, $\textbf{B} = diag(\frac{1}{\sigma_i})$

この$\textbf{A}$のような行列を、中心化作用素行列といいます。$\frac{1}{n} \textbf{1}_n$は全ての要素が$\frac{1}{n}$である$n$次元ベクトルで、もしこれだけを行列に左から掛けると以下のように、それぞれの列の平均値が要素であるm次元ベクトルになることがわかります。

$$ \begin{eqnarray}
\frac{1}{n} \textbf{1}_n^T \textbf{X} &=& \frac{1}{n} (1, … 1)
\begin{pmatrix}
x_{11} & \cdots & x_{1m} \\
\vdots & & \vdots \\
x_{n1} & \cdots & x_{nm} \\
\end{pmatrix} \\
&=& ( \sum_{i=1}^n \frac{1}{n} x_{i1}, … ,\sum_{i=1}^n \frac{1}{n} x_{im} ) \\
&=& (\frac{1}{n} \sum_{i=1}^n x_{i1}, … ,\frac{1}{n} \sum_{i=1}^n x_{im} )
\end{eqnarray} $$

$\textbf{A} \textbf{X}$を考えると、これは

$ \begin{eqnarray}
\textbf{A} \textbf{X} &=& \textbf{E}_n \textbf{X} – \frac{1}{n} \textbf{1}_n \textbf{X} \\
&=& \textbf{X} – \frac{1}{n} \textbf{1}_n \textbf{X}
\end{eqnarray}
$

となり、これは行列$\textbf{X}$の全ての要素から、その要素の列全体の平均値を引いた行列になります。

また$\textbf{B} = diag(\frac{1}{\sigma_i})$について考えると、これは、

$ \begin{pmatrix}
\frac{1}{\sigma_{1}} & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0 & \cdots & \frac{1}{\sigma_{p}} \\
\end{pmatrix} $

ですが、これを右から掛けることにより、行列$\textbf{X}$の全ての要素を、要素の列全体の標準偏差で割るということになります。

Pythonで実際にやってみます。データはおなじみirisを使っていきます。

import seaborn as sns
import pandas as pd
iris = sns.load_dataset('iris')

iris.describe()
sepal_lengthsepal_widthpetal_lengthpetal_width
count150.000000150.000000150.000000150.000000
mean5.8433333.0573333.7580001.199333
std0.8280660.4358661.7652980.762238
min4.3000002.0000001.0000000.100000
25%5.1000002.8000001.6000000.300000
50%5.8000003.0000004.3500001.300000
75%6.4000003.3000005.1000001.800000
max7.9000004.4000006.9000002.500000

これを標準化するのは以下のようにscikit-learnを使えば簡単にできます。

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
data_std = sc.fit_transform(iris.iloc[:,:-1])

pd_data_std = pd.DataFrame(data_std)
pd_data_std.describe()
0123
count1.500000e+021.500000e+021.500000e+021.500000e+02
mean-4.736952e-16-7.815970e-16-4.263256e-16-4.736952e-16
std1.003350e+001.003350e+001.003350e+001.003350e+00
min-1.870024e+00-2.433947e+00-1.567576e+00-1.447076e+00
25%-9.006812e-01-5.923730e-01-1.226552e+00-1.183812e+00
50%-5.250608e-02-1.319795e-013.364776e-011.325097e-01
75%6.745011e-015.586108e-017.627583e-017.906707e-01
max2.492019e+003.090775e+001.785832e+001.712096e+00


ではこれを、左右から行列を掛けることで再現していきます。

X = iris.drop(columns="species")

n = 150
p = 4
A = np.identity(n) - (1/n)*np.ones(n)

pd.DataFrame(A.dot(X.values)).describe()
0123
count1.500000e+021.500000e+02150.0000001.500000e+02
mean4.263256e-16-1.136868e-150.000000-4.263256e-16
std8.280661e-014.358663e-011.7652987.622377e-01
min-1.543333e+00-1.057333e+00-2.758000-1.099333e+00
25%-7.433333e-01-2.573333e-01-2.158000-8.993333e-01
50%-4.333333e-02-5.733333e-020.5920001.006667e-01
75%5.566667e-012.426667e-011.3420006.006667e-01
max2.056667e+001.342667e+003.1420001.300667e+00

$\textbf{A}\textbf{X}$で、平均値が0になりました。標準偏差は変化していません。

sigmas = [1/(X[c].std()) for c in X.columns]
B = np.diag(sigmas)

df_axb = pd.DataFrame(A.dot(X.values).dot(B))
df_axb
0123
0-0.8976741.015602-1.335752-1.311052
1-1.139200-0.131539-1.335752-1.311052
2-1.3807270.327318-1.392399-1.311052
3-1.5014900.097889-1.279104-1.311052
4-1.0184371.245030-1.335752-1.311052
1451.034539-0.1315390.8168591.443994
1460.551486-1.2786800.7035640.919223
1470.793012-0.1315390.8168591.050416
1480.4307220.7861740.9301541.443994
1490.068433-0.1315390.7602110.788031

統計量を調べてみます。

df_axb.describe()
0123
count1.500000e+021.500000e+021.500000e+021.500000e+02
mean5.684342e-16-2.676378e-15-9.473903e-17-5.684342e-16
std1.000000e+001.000000e+001.000000e+001.000000e+00
min-1.863780e+00-2.425820e+00-1.562342e+00-1.442245e+00
25%-8.976739e-01-5.903951e-01-1.222456e+00-1.179859e+00
50%-5.233076e-02-1.315388e-013.353541e-011.320673e-01
75%6.722490e-015.567457e-017.602115e-017.880307e-01
max2.483699e+003.080455e+001.779869e+001.706379e+00

平均値が0、標準偏差が1になっているので、標準化できています。