ゲルシュゴリンの定理というものがある。名前がなんというか、印象的な定理で、内容も分かりやすくて気に入っている。 いくつかの行列の例に対してPythonで検証してみる。
まずテスト用の行列を簡単に準備できるように、ランダム行列を生成する関数を用意しておく。
import numpy as np def generate_random_square_matrix_of_size(n): mat_real = np.matrix(np.random.normal(0, 1, (n, n))) mat_imag = np.matrix(np.random.normal(0, 1, (n, n))) mat = mat_real + 1j * mat_imag return mat
関数名からわかるだろうが、n x nのランダム行列を作っている。各乱数は中心0、標準偏差1の正規分布とした。
Pythonでは虚数をjであらわすが、ここで、単にjと書くだけではだめで、1jという風に、実数の表現の直後にjを後置することに注意する。j単体では変数名と解釈されてしまう。
上の関数が正規分布でデータを作っていることを確認しておく。
import numpy as np import tqdm def check_distribution(pngfile): values = [np.random.choice(np.ravel( generate_random_square_matrix_of_size(1))).real for _ in tqdm.tqdm(range(50_000), desc='Generating')] plt.rcParams['font.size'] = 20 fig = plt.figure() ax = fig.add_subplot(111) ax.hist(values, bins=200, density=True) ax.set_xlabel('value') ax.set_ylabel('density') ax.grid(True) fig.tight_layout() plt.savefig(str(pngfile)) plt.cla() plt.clf() plt.close()
確かに正規分布らしい分布になっている。 念のため、複素平面上の分布も確認しておく。
import numpy as np import tqdm def check_with_scatter(pngfile): values = [np.random.choice(np.ravel( generate_random_square_matrix_of_size(1))) for _ in tqdm.tqdm(range(10_000), desc='Generating')] plt.rcParams['font.size'] = 18 fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(np.real(values), np.imag(values), s=4, c='tab:orange') ax.set_aspect('equal') ax.grid(True) ax.spines['left'].set_position('zero') ax.spines['right'].set_color('none') ax.spines['bottom'].set_position('zero') ax.spines['top'].set_color('none') ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') ax.plot(1, 0, '>k', transform=ax.get_yaxis_transform(), clip_on=False) ax.plot(0, 1, '^k', transform=ax.get_xaxis_transform(), clip_on=False) ax.set_xlabel('Re', position=(1, 0), ha='right') ax.set_ylabel('Im', position=(0, 1), ha='right', rotation=0) plt.savefig(str(pngfile)) plt.cla() plt.clf() plt.close()
大丈夫そうだ。
それでは早速、行列を作って検証してみる。 行列の固有値はnp.linalg.eigvals()関数で求まる。これを、先ほどの例と同じように、実部と虚部に分けてプロットする。
eigvals = np.linalg.eigvals(mat)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(np.real(eigvals), np.imag(eigvals))
続いてゲルシュゴリンの定理に出てくる円を描画する。ゲルシュゴリンの定理の円は、対角成分を中心として、半径を、その行の非対角成分の絶対値の和とする円である。定義に従って計算して、先ほどのaxに描画してみる。
center_list = []
radius_list = []
for i in range(len(mat)):
center_list.append(mat[i, i])
# i-th radius is the sum of absolute values of
# i-th row values except diagonal element
radius = np.sum(np.abs(
np.where(np.arange(len(mat)) == i, 0, mat[i, :])))
radius_list.append(radius)
theta_data = np.linspace(0, 2 * np.pi, 50)
dr_data = np.exp(1j * theta_data)
for c, r in zip(center_list, radius_list):
z_data = c + r * dr_data
ax.plot(np.real(z_data), np.imag(z_data))
あとは軸のデザインなどを適宜調整して描画する。
無事描画できた。上図の場合、円が全て連結してしまっているので、連結成分に関する主張までは確認できないが、３つの円の中に合計で３つの固有値が入っている、ということは確認できた。
乱数のseedを変えていくつか試していて気付いたが、正規分布だとどうしても円が重なりやすい。 対角成分のスケールを思い切って大きくすればばらけるだろう。あるいは非対角成分を小さくすることで個別の半径を小さくするという考え方もできる（その極限が対角行列で、対角成分が固有値になる）。
三重対角行列のように少しだけ非対角要素に染み出していたり、行方向/列方向に規則が見られる状況だと面白くなるように思う。
連結性の自動判定や、対角成分の分布を変える検証は、また気が向いた際に行うことにする。
参考：
