©作者 | 刘敏
学校 | 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
但这个方法其实很不好用,主要有下面这两个理由:
对存储的要求非常大,也就是 。比如说我们需要大概 80G 去存储一个 1w 乘 1w 的矩阵。
构造出来的是一个 dense 的矩阵,也就是大多数分量都不是零!那当我们要去算 的时候,我们基本只能使用最基本的算特征值的方法,复杂度就是 !
那我们需要找点其他方法搞搞,也就是:是否可以找到一个矩阵
他对存储的要求比较低。
他有点特殊,可以用一些算法复杂度比较低的方法来算他最大的特征值。
他最大的特征值的分布是等于 的分布的。
那在 [1] 里面, 两位作者证明了下面这个矩阵是满足这三个要求的:
这里面 是一个高斯分布,他的期望是 0 方差是 2, 是 chi square distribution with freedom 的平凡根,。这里要注意的是:
和 随机变量都是两两互为独立的。
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)我们可以发现:
我们只需要 的存储空间。
他具有 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」
点击「关注」订阅我们的专栏吧
·