Pythonで線形代数！〜行列・応用編（行列式・固有値）：数学×Pythonプログラミング入門（5/5 ページ）

[羽山博，著]

練習問題

　では、練習問題に取り組みましょう。問題の考え方について解説した動画も用意してあります。ぜひご視聴ください。

動画4　行列応用の練習問題

（1）行列式を使ってハート形の面積を求める

　この練習問題は目標2の延長線上にあるものです。目標2では単純な図形を使いましたが、ここでは複雑な図形でやってみましょう。図6のようなハート形の内部の面積を求めてください。

図6　ハート形の面積を求める
この図形の頂点は565個ある。データは、リスト12のコードに記載されているように、CSVファイルとして提供しているので、それを利用するとよい。

　このハート形の頂点データはリスト12のコードで取得できます。

import numpy as np

data = np.loadtxt('https://raw.githubusercontent.com/Gessys/math/main/data/heart.csv', skiprows=1,delimiter=',')

print(data)
# 出力例：
# [[ 0.         -1.        ]
#  [ 0.01       -0.97073482]
#  [ 0.02       -0.95348411]
#  ...
#  [-0.02       -0.95348411]
#  [-0.01       -0.97073482]
#  [ 0.         -1.        ]]

リスト12　ハート形の頂点データを取得するためのコード
0列目がx座標、1列目がy座標となっている。最初の三角形は原点と(0,−1)、(0.01,−0.97073482)で作られ、次の三角形は原点と(0.01,−0.97073482)、(0.02,−0.95348411)で作られる。

　ちなみに、このデータは以下の式に基づいて作成したものです*2。macOSに付属のGrapherで積分の近似値を求めると4.443ぐらいになるので、それに近い値が得られれば正解です。

*2　このデータは、コードを書きやすくするため、座標の絶対位置を並べたものにしてありますが、先頭だけを絶対的な位置にし、それ以降は移動量を記録しておいた方が取り扱いが便利です（先頭の位置を変えるだけで図形が移動できます）。また、このCSVファイルの内容をSVGファイルに取り込めば、ベクトルグラフィックスを取り扱うためのソフトウェアやブラウザーで表示できるようになります。ただし、このデータはy軸の下の方が小さな値になっているので、グラフィックを取り扱う環境でy軸の上の方が小さな値になっている場合には、上下が逆さまになります。……というのは全くの余談ですが。

（ヒント）行列を転置しても行列式の値は変わりません。これまでのお話では下の式の左辺の形で行列式を求めていましたが、右辺の形でも結果は同じになります。

（2）分散・共分散行列の固有値と固有ベクトルを求めて主成分分析を行う

　固有値を利用すると、主成分分析と呼ばれる分析を行うことができます。そこで、以下のようなデータを使って主成分分析を行ってみましょう。これは、連載第5回で散布図を作成するために使った成績データです。散布図も併せて描いておきます（リスト13、図7）。

import numpy as np
import matplotlib.pyplot as plt

english = np.array([98, 77, 64, 45, 58, 62, 89, 30])
mathmatics = np.array([88, 54, 72, 85, 70, 81, 63, 27])

plt.scatter(english, mathmatics, marker='*')
plt.xlabel('English')
plt.ylabel('Math')
plt.show()

リスト13　主成分分析に使うためのデータ
8人の英語と数学の成績。これらのデータの特徴をうまく表すような直線を求めようというのが主成分分析。

図7　英語と数学の成績をプロットした散布図
散布図を見ると、なんとなく「勉強が得意か不得意か」といった特徴と、「文科系・理数系」といった特徴が見て取れそうではある。

　またまたいきなりですね。主成分分析とはいったい何か、というところからスタートしないと手の付けようがないので、簡単に説明します。主成分分析とは、データが持つ情報をできるだけ失うことなく、

という式を作り、a₁, a₂, ... a_nを求めることです。

　主成分分析は、多次元のデータの次元を削減するために使われます（図7の例は2次元なので、次元を削減してもあまり意味がなさそうですが、例えば、このデータの特徴を「勉強が得意か不得意か」だけで表そうということです）。

　ここでは、変数が英語の成績と数学の成績だけなので、英語の成績をx、数学の成績をyとして、上の式を以下のように表すことにします（添え字が多いと煩雑になるのでx₁を単にxとし、x₂をyとしました）。

　主成分の求め方と、主成分を解釈するために必要な主成分得点を求める手順を示します。

• ステップ1： numpyモジュールのcov関数を使って分散・共分散行列を求める
• ステップ2： 分散・共分散行列の固有値を求める
• ステップ3： 固有値の大きいものに対する固有ベクトルを求める（これが第1主成分となる）
• ステップ4： 次に大きい固有値に対する固有ベクトルを求める（これが第2主成分となる）
• ステップ5： 主成分と各データを使って主成分得点を求める

　n次元のデータであれば最大n個の主成分が求められますが、重要なのは大きな固有値に対する主成分です。ステップ1〜4まではこれまでの知識でできると思います。ステップ5については、例えば第1主成分を表す固有ベクトルが(a₁,a₂)であるとき、各データを（4）式に代入し、

を求め、z₁〜z_nの平均z̄を計算して、

とします。主成分得点は、固有ベクトルを座標軸としたときの各データの位置を表すものと考えられます。

　ここでは、主成分分析の仕組みと計算方法について詳しく説明する余裕がないので、理屈抜きで、情報損失量を最小にするような直線（主成分）を求めると、その係数a₁,a₂,..., a_nが、分散・共分散行列の固有ベクトルになっていると考えて下さい*3。

*3　主成分分析の方法は、重回帰分析と似ています（が、少し違います）。重回帰分析では実測値と回帰式の「高さ」の差を最小にしますが、主成分分析では、直線との距離、つまり各データから直線に下ろした垂線の長さを最小にするような計算になります。そうすれば、各データの持つ情報が最大限、その直線に反映されるというわけです。斜めから日が差し込むより、真上から（垂直に）差し込む方が日光のエネルギーを最大限受けられます。あくまで例えですが、それと同じようなことです。

　取りあえず、これでコードが書けると思います。プログラムの作成に取り組んでみましょう。

練習問題の解答例

　以下、解答とプログラムの作成例です。もちろん、異なるやり方もあるので、これらが唯一の答えというわけではありません。

（1）行列式を使ってハート形の面積を求める

　これは簡単です。読み込んだデータの0行目と1行目の行列式を求め、1行目と2行目の行列式を求め、2行目と3行目の行列式を求め……という具合に進めて、それらの値を合計すればいいですね。合計した値は平行四辺形の面積なので、最後に2で割ることを忘れずに。

import numpy as np

data = np.loadtxt('https://raw.githubusercontent.com/Gessys/math/main/data/heart.csv', skiprows=1,delimiter=',')

total = 0
rows = data.shape[0] # 行数を得る
for i in range(rows-1):
  total += np.linalg.det(data[i:i+2])
print(total / 2)
# 出力例：
# 4.4401816396000084

リスト14　ハート形の頂点データから内部の面積を求めるコード
data[i:i+2]はi行目からi+2行目の手前までという意味。行列式の値を求めてtotalに加えていく。出力例を見ると、4.443に近い値が得られたことが分かる。

　実は、図6を見ても分かるのですが、この例では行列式が負になる部分が1つもありません。が、それは問題ではありません。負になる部分があっても、目標2で見たようにちゃんと面積が求められます。

（2）分散・共分散行列の固有値と固有ベクトルを求めて主成分分析を行う

　主成分分析という聞き慣れない手法が登場したので、身構えてしまった（固まってしまった）方もおられるかもしれませんが、手順通りに進めればできます。

import numpy as np

english = np.array([98, 77, 64, 45, 58, 62, 89, 30])
mathmatics = np.array([88, 54, 72, 85, 70, 81, 63, 27])

V = np.cov(english, mathmatics)  # 分散・共分散行列を求める
print("分散・共分散行列：")
print(V)

DP = np.linalg.eig(V)  # 固有値／固有ベクトルを求める
print("固有値／固有ベクトル：")
print(DP)
# 出力例：
# 分散・共分散行列：
# [[498.83928571 201.64285714]
#  [201.64285714 396.85714286]]
# 固有値／固有ベクトル：
# (array([655.83843329, 239.85799528]), array([[ 0.78903768, -0.61434481],
#        [ 0.61434481,  0.78903768]]))

リスト15　主成分分析を行う
分散・共分散行列はnumpyモジュールのcov関数で求める。出力例の対角要素が分散、左下と右上が共分散となる。固有値／固有ベクトルの求め方はこれまでと同じ。

　固有値の大きな方は655.838...です。それに対する固有ベクトルは(0.789...,0.614...)となっています。これが第1主成分のa₁とa₂に当たります。第2主成分は、固有値239.857...に対する固有ベクトル(−0.614...,0.789...)で求められます。これらの固有ベクトルを見ると、第1主成分と第2主成分は直交することが分かります。

　次に、主成分得点です。最大の固有値を第1主成分、最小の固有値を第2主成分とするのは固有値が2個しかない場合にしか使えないのであまりスマートなやり方ではありませんが、この例ではうまくいきます。リスト15に続けてリスト16のコードを入力してください。

# 第1主成分
idx = DP[0].argmax()  # 固有値が2つなので最大の方を第1主成分とする
a1= DP[1][0, idx]
a2= DP[1][1, idx]
# 全てのzを求める（暫定）
z = a1 * english + a2 * mathmatics
# zの平均を求める
zmean = z.mean()
# zを求める
print("第1主成分")
print(a1, a2)
print("第1主成分の主成分得点")
print(z - zmean)

# 第2主成分
idx = DP[0].argmin()  # 固有値が2つなので最小の方を第2主成分とする
a1= DP[1][0, idx]
a2= DP[1][1, idx]
# 全てのzを求める（暫定）
z = a1 * english + a2 * mathmatics
# zの平均を求める
zmean = z.mean()
# zを求める
print("第2主成分")
print(a1, a2)
print("第2主成分の主成分得点")
print(z - zmean)
# 出力例：
# 第1主成分
# 0.7890376766407048 0.6143448094038383
# 第1主成分の主成分得点
# [ 38.33642279   0.87890806   1.67962484  -5.3256085   -4.28329084
#    5.63065277  15.87646347 -52.79317259]
# 第2主成分
# -0.6143448094038383 0.7890376766407048
# 第2主成分の主成分得点
# [ -3.86772704 -17.79376704   4.39539366  26.32543483   6.50338716
#   12.72542237 -18.06456567 -10.22357827]

リスト16 主成分得点を求める
argmaxメソッドで固有値の最大値のインデックスを求め、それに対応する固有ベクトルの要素を第1主成分としてa1、a2に代入する。あとは手順通りにzの値を求める。第2主成分については、argminメソッドで固有値の最小値のインデックスを求め、それに対応する固有ベクトルの要素をa1、a2に代入する。こちらも手順通りにzの値を求める。

　第1主成分の主成分得点を見ると、点数の高い人は主成分得点が大きく、点数の低い人は主成分得点が小さいことが分かります。例えば、最初の英語98点、数学88点の人の主成分得点は38.336...です。一方、最後の英語30点、数学27点の人の主成分得点は-52.973...です。従って、第1主成分は、やはり「勉強が得意か不得意か」という特徴を表しているようです。

　一方、第2主成分の主成分得点は、英語よりも数学の成績がいい人の主成分得点が大きく、数学よりも英語の成績がいい人の主成分得点が小さくなっています。例えば、英語が45点、数学が85点の人の主成分得点は26.325...で、逆に、英語が89点、数学が63点の人の主成分得点は-18.064...という値でした。こちらも当初の予想通り「文科系・理数系」といった特徴を表しているようです。

---

　……というわけで、今回は、AI／機械学習に必要とされる線形代数の知識のうち、行列式や固有値／固有ベクトルなどの取り扱いについて、基本的な考え方と関数を使って簡単に計算する方法、幾つかの応用例を見てきました。手作業での計算が大変なので難しいとか面倒だという印象はあるかもしれませんが、関数を使ってとにかく計算してみると、意外に簡単だということが分かるのではないでしょうか。

　さて、次回ですが、しばらく線形代数の話題が続いたので、気分を変えて三角関数の利用について見ていくことにします。

「数学×Pythonプログラミング入門」