锐单电子商城 ， 一站式电子元器件采购平台！

  • 电话：400-990-0325

购物车

一些关于随机矩阵的算法

时间：2022-08-03 17:19:00 二极管kpz

cbbbc2a09be3c0de8b9cb332f53d0114.png

来源：PaperWeekly 本文约1500字，建议阅读5分钟 random matrix 的算法。

本文介绍了我在硕士论文中使用的随机矩阵 GUE 算法，真的超级好用，谁用谁知道！ GUE 简单介绍，可见：

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

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

让我们先回顾一下，GUE 的定义：

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

本文介绍了我在硕士论文中使用的随机矩阵 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.

编辑：于腾凯

相关文章