とあるランダム行列の固有値は円板状に分布するらしい。正しくはこう
円則
n×n 実正方行列(または複素正方行列)において各行列要素を独立同一分布で平均ゼロ $E(X_j,k)=0$、分散 $E(|X_{j,k}|^{2})=1/n$ のように規格化すると、行列のサイズを非常に大きくしていく(n → ∞)に従い固有値は複素平面上の単位円盤上で一様に分布するようになるというもの。
引用元: ランダム行列#円則
早速実装して実験みた。
使用言語はHaskell、線形代数ライブラリはhmatrix、可視化ライブラリにはChartを使った。
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Diagrams
import Numeric.LinearAlgebra
-- | 複素数x+iyを(x, y)型のタプルに変換
toPoint :: Complex a -> (a, a)
toPoint c = (realPart c, imagPart c)
main = do
let n = 1000 -- 正方行列のサイズ
m <- (sqrt (1 / fromIntegral n) *) <$> randn n n -- 行列要素が平均値0分散1/nの正規分布に従う行列をサンプリング
let es = map toPoint . toList . fst $ eig m -- 行列の固有値を求めて[(Double, Double)]型に変換
toFile (def {_fo_size = (600, 600)}) "output.svg" $ -- プロットしてファイルに書き出す
plot $ points "Eigenvalues" es
実験に使う関数は以上。
hmatrixでランダム行列を計算する方法について解説しておくと、まずは
randomVector :: Seed -- シード値
-> RandDist -- 値が従う分布
-> Int -- 生成するベクトルの長さ
-> Vector Double
というランダムなベクトルを計算する関数が最も基本的な関数になる。
type Seed = Int
data RandDist = Uniform
| Gaussian
となっており randomVector は値が一様分布もしくは正規分布に従うランダムなベクトルを生成できることが分かる。
ランダム行列を生成する場合は randomVector を利用した簡便な関数が用意されている。
-- | 行列要素が[0,1]区間の一様分布に従うようなランダム行列を生成する
rand :: Int -> Int -> IO (Matrix Double)
-- | 行列要素が標準正規分布に従うようなランダム行列を生成する
randn :: Int -> Int -> IO (Matrix Double)
さらにこれらの行列を変形してより一般的な行列を得られる関数も存在する。
uniformSample :: Seed -- シード値
-> Int -- 行数
-> [(Double, Double)] -- 各列が従う一様分布の区間
-> Matrix Double
gaussianSample :: Seed -- シード値
-> Int -- 行数
-> Vector Double -- 平均値ベクトル
-> Herm Double -- 分散共分散行列
-> Matrix Double
uniformSampleは[0,1]区間に限らない一様分布に従う行列要素をもつランダム行列、 gaussianSample は各列が自明でない平均値と分散共分散構造をもつ多変量正規分布に従うようなランダム行列を計算する。
実験
行列のサイズを大きくしていきながら固有値の振る舞いを観察する。
n = 10
n = 100
n = 500
n = 1000
固有値が綺麗に円板状に分布していく様子が見て取れる。すごい。今回は実数を行列要素に持つ正方行列で実験したので、よく見ると固有値の分布がx軸に対して線対称になっていることが分かる。円則の証明は Terence Tao, Van Vu "Random Matrices: The circular Law" に載っているらしい。