線形代数の練習問題を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*earray([[ 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} $になることがわかる。
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_length sepal_width petal_length petal_width count 150.000000 150.000000 150.000000 150.000000 mean 5.843333 3.057333 3.758000 1.199333 std 0.828066 0.435866 1.765298 0.762238 min 4.300000 2.000000 1.000000 0.100000 25% 5.100000 2.800000 1.600000 0.300000 50% 5.800000 3.000000 4.350000 1.300000 75% 6.400000 3.300000 5.100000 1.800000 max 7.900000 4.400000 6.900000 2.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()0 1 2 3 count 1.500000e+02 1.500000e+02 1.500000e+02 1.500000e+02 mean -4.736952e-16 -7.815970e-16 -4.263256e-16 -4.736952e-16 std 1.003350e+00 1.003350e+00 1.003350e+00 1.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-01 3.364776e-01 1.325097e-01 75% 6.745011e-01 5.586108e-01 7.627583e-01 7.906707e-01 max 2.492019e+00 3.090775e+00 1.785832e+00 1.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()0 1 2 3 count 1.500000e+02 1.500000e+02 150.000000 1.500000e+02 mean 4.263256e-16 -1.136868e-15 0.000000 -4.263256e-16 std 8.280661e-01 4.358663e-01 1.765298 7.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-02 0.592000 1.006667e-01 75% 5.566667e-01 2.426667e-01 1.342000 6.006667e-01 max 2.056667e+00 1.342667e+00 3.142000 1.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_axb0 1 2 3 0 -0.897674 1.015602 -1.335752 -1.311052 1 -1.139200 -0.131539 -1.335752 -1.311052 2 -1.380727 0.327318 -1.392399 -1.311052 3 -1.501490 0.097889 -1.279104 -1.311052 4 -1.018437 1.245030 -1.335752 -1.311052 ... ... ... ... ... 145 1.034539 -0.131539 0.816859 1.443994 146 0.551486 -1.278680 0.703564 0.919223 147 0.793012 -0.131539 0.816859 1.050416 148 0.430722 0.786174 0.930154 1.443994 149 0.068433 -0.131539 0.760211 0.788031
統計量を調べてみます。
df_axb.describe()0 1 2 3 count 1.500000e+02 1.500000e+02 1.500000e+02 1.500000e+02 mean 5.684342e-16 -2.676378e-15 -9.473903e-17 -5.684342e-16 std 1.000000e+00 1.000000e+00 1.000000e+00 1.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-01 3.353541e-01 1.320673e-01 75% 6.722490e-01 5.567457e-01 7.602115e-01 7.880307e-01 max 2.483699e+00 3.080455e+00 1.779869e+00 1.706379e+00
平均値が0、標準偏差が1になっているので、標準化できています。
