一些关于随机矩阵的算法

©作者 | 刘敏

学校 | Universität Bonn

研究方向 | Probability

本文介绍一下我硕士论文中用到的关于随机矩阵 GUE 的算法，真的超级好使，谁用谁知道！关于 GUE 的简单介绍，可以看下：

https://zhuanlan.zhihu.com/p/161375201

这篇文章的主要参考文献是  [1][2][3]  。所有代码都是使用 Matlab 编写。

那我们首先来回顾一下，GUE 的定义：

DEFINITION 1.1（Gaussian unitary ensemble）假设 是独立同分布的标准高斯随机变量（期望为 0，方差为 1），那么 的 GUE 就被定义为

或者展开写的话就是

然后呢，我们比较关心的东西是他的最大的那个特征值，我们表示为 。那最简单的方法来实现这个矩阵，就是按照定义，也就是

function GUE = GUE_matrix_MC_create_GUE(size,seed)
    %set random seed
    rng(seed);
    tempMat=randn(size)+1i*randn(size);
    GUE=(tempMat+tempMat')/2;
end

但这个方法其实很不好用，主要有下面这两个理由：

1. 对存储的要求非常大，也就是 。比如说我们需要大概 80G 去存储一个 1w 乘 1w 的矩阵。
2. 构造出来的是一个 dense 的矩阵，也就是大多数分量都不是零！那当我们要去算 的时候，我们基本只能使用最基本的算特征值的方法，复杂度就是 !
   

那我们需要找点其他方法搞搞，也就是：是否可以找到一个矩阵

1. 他对存储的要求比较低。

   
   

2. 他有点特殊，可以用一些算法复杂度比较低的方法来算他最大的特征值。 

   
   

3. 他最大的特征值的分布是等于 的分布的。

那在 [1] 里面， 两位作者证明了下面这个矩阵是满足这三个要求的：

这里面 是一个高斯分布，他的期望是 0 方差是 2， 是 chi square distribution with freedom 的平凡根， 。这里要注意的是：

1. 和 随机变量都是两两互为独立的。

2. sub-digonal 和 super-digonal 上是相等的！

那我们就可以通过下面这段代码来实现他的构造：

function triMat = GUE_matrix_MC_create_TriMat(size,seed)
    %set random seed
    rng(seed);
    %set subdiagonal/superdigonal as chi-distributed
    d=sqrt(1/2)*sqrt(chi2rnd(beta*[size:-1:1]))';
    %set up digonal
    d1=(randn(size,1));
    triMat=spdiags(d,1,size,size)+spdiags(d1,0,size,size)+spdiags(d,1,size,size)';
end

这个方法确实好，通过观察（2.1）我们可以发现：

1. 我们只需要 的存储空间。
   
   
2. 他具有 tridigonal 和 irreducible 的结构（因为他的 sub-digonal 上的元素 a.s. 不等于 0），那我们就可以用一些比较厉害的算法来计算他最大的特征值了！比如说 bisection method（这个方法真的不错，感兴趣的可以看看这本书的 [4]  lecture 30），他的算法复杂度只有 。

当然也可以用 Matlab 自带的算最大特征值的函数 ，但这个的复杂度是 。在下面这个 图里面，我比较了一下他们三者的算法复杂度，也就是最原始的 GUE + ，（2.1）+ 以及（2.1）+ bisection method，然后矩阵的大小 ，测时的方法就是 Matlab 自带的 。

▲ 从上到下依次为GUE+eigs, (2.1)+eigs以及(2.1)+bisection，我们可以看到他们的算法复杂度分别为n^3, n^2以及n.

关于 bisection method 的代码我就不贴了吧，毕竟我也是从别人那里下载的，如果大家想下载的话，可以去 [2] 的作者主页下载（ http://www.mit.edu/ ）。

但是上面三个方法本质上都是对 Monte Carlo 方法的修修补补，并不能克服 Monte Carlo 方法自身的 的趋近速度，我当时想要得到一个 2000 乘 2000 矩阵的靠谱数值，花了大概七八个小时。所以我们需要点新东西，那接下来要介绍的方法就有点厉害了，完全换了一个思路！

首先，我们其实已经知道 的分布函数，我们只是想研究他的一些其他特质，那我们为什么不能直接从他的分布函数入手呢？如果这个可行的话，那我们完全可以不再用 Monte Carlo 方法啊，那我们回顾一下， 的分布函数可以写成 Fredholm determinant：

这里

其中

是第 个 oscillator wave-function， 是第 个 Hermite polynomial。那我们观察 发现，右边的那个积分是一个 维的积分，那对于这种积分，我们可是有很有效的方法来模拟的！比如说 Gauss-Legendre 或者 r Curtis-Clenshaw，也就是说，我们可以把式子 右边近似为

那现在的问题就是，这个误差有多少，趋近的有多快啊？那在 [3] 中，Bornemann 证明了 （其实他证明了一个更一般的情形，这里为了表述方便我就取一个特殊形式了）：

THEOREM 2.1 假设 是 analytic 且 exponential decay，也就是存在 使得

那么 趋近于 exponentially fast。

那对于定义在 中的 ，他是满足这个方法的，所以我们可以用这种方法来算他的分布！进而可以算他的期望或者其他的一些性质！这个方法真的超级快，算一个 2000*2000 矩阵的最大特征值的期望可能不需要两秒吧！并且这个方法不仅仅适用于 random matrix，在 KPZ-model 里面，大部分的 kernel 都满足这个性质，比如说对于 TASEP 的分布，我们可以通过下面这个代码来实现：

function [result] = step_TASEP_cdf(sigma,t,s)

s=step_TASEP_proper_interval(t,sigma,s);
c2=sigma^(-1/6)*(1-sigma^(1/2))^(2/3);
delta_t=c2^(-1)*t^(-1/3);
n=sigma*t;
MAX=(t+n-2*(sigma)^(1/2)*t-1/2)/(c2*t^(1/3));

for k=1:length(s)   
    if s(k)> MAX
        result(k)=1;
    else
        s_resc=s(k)+delta_t;
        x=s_resc:delta_t:MAX;

        x=x';
        result(k)=det(eye(length(x))-step_TASEP_kernel(t,sigma,x,x)*delta_t);%Bornemann Method
     end
end

end

代码中标注为 Bornemann method 的地方就是用的上面说的方法，这里面我们不需要选取 是因为这个分布函数本身就是个离散的。强烈建议大家读一下 [3] ，一定会有很大的收获！并且一定会有意外收获！可以点击这里 阅读：

https://arxiv.org/pdf/0804.2543.pdf

这篇文章就是简单的介绍了一下有关于 random matrix 的算法，之后可能会陆续介绍一下 KPZ-universality 相关的东西，也就是我自己的方向，真的超级有趣！

参考文献

[1] Dumitriu I, Edelman A. Matrix models for beta ensembles[J]. Journal of Mathematical Physics, 2002, 43(11): 5830-5847.

[2] ersson P O. Random matrices. Numerical methods for random matrices[J]. 2002.

[3] Bornemann F. On the numerical evaluation of Fredholm determinants[J]. Mathematics of Computation, 2010, 79(270): 871-915.

[4] Trefethen L N, Bau III D. Numerical linear algebra[M]. Siam, 1997.

更多阅读

#投 稿 通 道#

 让你的文字被更多人看到 

如何才能让更多的优质内容以更短路径到达读者群体，缩短读者寻找优质内容的成本呢？答案就是：你不认识的人。

总有一些你不认识的人，知道你想知道的东西。PaperWeekly 或许可以成为一座桥梁，促使不同背景、不同方向的学者和学术灵感相互碰撞，迸发出更多的可能性。 

PaperWeekly 鼓励高校实验室或个人，在我们的平台上分享各类优质内容，可以是最新论文解读，也可以是学术热点剖析、科研心得或竞赛经验讲解等。我们的目的只有一个，让知识真正流动起来。

📝 稿件基本要求：

• 文章确系个人原创作品，未曾在公开渠道发表，如为其他平台已发表或待发表的文章，请明确标注 

• 稿件建议以 markdown 格式撰写，文中配图以附件形式发送，要求图片清晰，无版权问题

• PaperWeekly 尊重原作者署名权，并将为每篇被采纳的原创首发稿件，提供业内具有竞争力稿酬，具体依据文章阅读量和文章质量阶梯制结算

📬 投稿通道：

• 投稿邮箱：hr@paperweekly.site 

• 来稿请备注即时联系方式（微信），以便我们在稿件选用的第一时间联系作者

• 您也可以直接添加小编微信（pwbot02）快速投稿，备注：姓名-投稿

△长按添加PaperWeekly小编

🔍

现在，在「知乎」也能找到我们了

进入知乎首页搜索「PaperWeekly」

点击「关注」订阅我们的专栏吧

·