Python 绘制 Wigner Matrices

Python 绘制 Wigner Matrices

上随机矩阵暑期课程的时候，老师讲到Wigner Matrices的内容，刚好自己编程水平有限，就试着编程实现了一下。

Wigner Matrices定义：

令$X_N=(X_{ij})$是$N\times N$的对称方阵，对角线及其上元素均为i.i.d. 满足
$\mathbb{E}{X}_{ij}=0$ 且 $\mathbb{E}|X_{ij}{|}^2=1$
$X_{ij}=X_{ji}$ (对称性).

• Wigner Matrices的谱为：
  $Y_N=\frac{X_N}{\sqrt{N}}$

• semi-circular 分布为：
  $x\to \frac{\sqrt{4-x^2}}{2\pi }$

接下来我们绘图来验证半圆率，通俗的理解就是当$N$取足够大时，绘制上述方阵的特征值会在$[-2,2]$区间形成半圆，具体理论证明待补.

Wigner’s theorem (1948):
The histogram of a Wigner matrix converges to the semi-circular distribution！

问题:生成期望为0，方差为1的对称矩阵

网上很多关于对称矩阵的程序是简单用

n = 100 #任意整数
A = np.random.randn(n)
B = (A + A^T) / 2

完成的，但实际上，这样得到的样本方差并不是1，样本量足够大的时候是0.5.这样会导致直方图不是在$[-2,2]$中.

解决方案1：

所以我先生成与方阵上三角元素数量相等的服从N(0,1)的随机数，再通过for循环将元素写进预先生成的零矩阵中，具体程序如下

import numpy as np
import matplotlib.pyplot as plt
n = 1500  # 矩阵的维度
m = (n*n-n)/2+n
X = np.random.randn(int(m))
# print(X)
print('原始生成X的方差是：',np.var(X))

matrix = np.zeros((n, n))  # 创建一个n x n的全零矩阵
k = 0
for i in range(n):
    for j in range(i, n):    #从i开始到n结束
        matrix[i][j] = X[k] # 上三角部分赋值
        k = k+1
matrix += matrix.T - np.diag(matrix.diagonal())
print('生成对称矩阵的方差是：',np.var(matrix))

解决方案2：

一个更为简便的写法，直接进行矩阵分解，待补

绘图部分

鉴于自己总是忘记plt的用法，今天记录详细一些：
首先用

plt.figure(figuresize=(6,6))

确定画布尺寸，这句命令需要在绘图的最开始，否则产生空白画布.
其次，使用

plt.hist(eigenvalue, bins=50, density=True, label='spectrum', edgecolor = 'black')

来绘制直方图，hist是histogram的缩写，第一项参数是需要统计频率的值，bins代表需要绘制的条数，normed = true被移除了，改成density = True，使用edgecolor来指定条与条之间的分割线.
最后，使用

plt.savefig()

来保存图片，该命令需要在plt.show()前，否则会得到一张空白的画布，这是因为show本身产生一张空白画布.

完整程序

完整程序如下，在电脑上命名为semi-circle_law.py

import numpy as np
import matplotlib.pyplot as plt
n = 1500  # 矩阵的维度
m = (n*n-n)/2+n
X = np.random.randn(int(m))
# print(X)
print('原始生成X的方差是：',np.var(X))

matrix = np.zeros((n, n))  # 创建一个n x n的全零矩阵
k = 0
for i in range(n):
    for j in range(i, n):    #从i开始到n结束
        matrix[i][j] = X[k] # 上三角部分赋值
        k = k+1
matrix += matrix.T - np.diag(matrix.diagonal())
print('生成对称矩阵的方差是：',np.var(matrix))
#绘图
plt.figure(figsize=(6, 6))
Z = np.random.uniform(low=-2,high=2,size=n)
M = np.sqrt(4 - np.power(Z,2)) / (2*np.pi)
plt.scatter(Z,M,c='r')

Y = matrix / np.sqrt(n)
eigenvalue, eigenvector = np.linalg.eig(Y)
plt.hist(eigenvalue, bins=50, density=True, label='spectrum', edgecolor = 'black')
plt.legend()
plt.savefig('/Users/XXX/Desktop/matlab_test/Python/002.png') #需要在plt.show前
plt.show()