【L】on dynamical isometry
本篇博文作为a note on the learning of deep architectures的续篇研究了一般网络的dynamical isometry理论。
Jefferey Penington在2017年发表的论文[1]中进一步对非线性网络的dynamical isometry进行了研究。论文开篇提出了三个疑问：1.深度网络的Input-Output 雅可比矩阵的奇异值的分布与网络深度，初始化权值的统计结构以及非线性特性有什么关系? 2.什么样的网络配置可以达到dynamical isometry?3.那些不受梯度弥散或梯度爆炸影响的非线性网络中，达成dynamical isometry是否能使网络学习更迅速?
要给出这些问题的答案，我们先来研究非线性神经网络中，各层preactivation的分布以及不同example输入下preactivation的协方差的分布（实际上研究activation也能给出一样的结果）。这里我们先约定神经网络里的一层由上一层的postactivation $a^{(l-1)}$,权值矩阵$W^{(l)}$,偏置向量$b^{(l)}$组成，其中preactivation $h^{(l)}=W^{(l)}a^{(l-1)}+b^{(l)}$。激活函数为$\phi$。
现在我们定义preactivation 的方差的定义，第l层的preactivation为$q^{(l)}$,第l层关于example中a,b的协方差为$q^{(l)}_{ab}$。并且权值服从均值为0,方差为$\frac{{\sigma_{w}}^2}{N_{l-1}}$的分布,偏置服从均值为0,方差为${\sigma_{b}}^2$的分布。
\begin{align*}
q^{(l)}=\frac{1}{N_l}\sum_{i=0}^{N_l}{h^{(l)}_i}^2 \tag{1}\\
q^{(l)}_{ab}=\frac{1}{N_l}\sum_{i_0}^{N_l}h^{(l)}_a h^{(l)}_b \tag{2}
\end{align*}
另一方面，有
\[\begin{align*}
VAR(h^{(l)}_i)&=E({h^{(l)}_i}^2)-E({h^{(l)}}_i)^2\\&=E({h^{(l)}_i}^2)\\&=E(({W_{i,}^{(l)}}^Ta^{(l-1)})^T{W_{i,}^{(l)}}^Ta^{(l-1)})+E(b_i^2)\\&=N_{l-1}VAR(w)VAR(a^{(l-1)})+VAR(b)\\&=\frac{\sigma_{w}^2}{N_{l-1}}N_{l-1}VAR(a^{(l-1)})+\sigma_{b}^2\\&=\sigma_{w}^2VAR(a^{(l-1)})+\sigma_{b}^2
\end{align*}\]
又$a^{(l-1)}=\phi(h^{(l-1)})$,由于$h^{(l-1)}$等价于由同样分布独立随机变量相加，所以由中心极限定理，应该是服从均值为0,方差为$q^{(l-1)}$的高斯分布。于是有
\[
\begin{align*}
q^{(l)}&=\sigma_{w}^2\int d(z){\phi(\sqrt{q^{(l-1)}}z)}^2dz+\sigma_{b}^2 \tag{3}
\end{align*}
\]
其中，z为服从标准正太分布的随机变量，$d(z)=\frac{1}{\sqrt{2\pi}}e^{-z^2}$为标准正态分布的密度函数。
类似地，我们有
\[
\begin{align*}
COV(h^{(l)}_a,h^{(l)}_b)&=E(h^{(l)}_a h^{(l)}_b)-E(h^{(l)}_a)E(h^{(l)}_b)\\&=E(({W_{i,}^{(l)}}^Ta^{(l-1)}_a)^T{W_{i,}^{(l)}}^Ta^{(l-1)}_b)+E(b_i^2)\\&=N_{l-1}VAR(w)COV(a^{(l-1)}_a,a^{(l-1)}_b)+VAR(b)\\&=\sigma_{w}^2COV(a^{(l-1)}_a,a^{(l-1)}_b)+\sigma_{b}^2
\end{align*}
\]
其中由于$a^{(l-1)}=\phi(h^{(l-1)})$,而$h^{(l-1)}_a,h^{(l-1)}_b$各自是均值为0，方差为$q^{(l-1)}_a,q^{(l-1)}_b$的正态分布，协方差为$q^{(l-1)_{ab}}$。
不妨令
\[
\begin{align*}
h^{(l-1)}_a=mz_1+nz_2\\
h^{(l-1)}_b=gz_1+pz_2
\end{align*}
\]
其中,$z_1,z_2$均为相互独立的服从标准正态分布的随机变量（独立正态随机变量的和依然服从高斯分布），(m,n,g,p)为使得下边方程成立的元组:
\[
\begin{align*}
VAR(h^{(l-1)}_a)=q^{(l-1)}_a\\
VAR(h^{(l-1)}_b)=q^{(l-1)}_b\\
COV(h^{(l-1)}_a,h^{(l-1)}_b)=q^{(l-1)}_{ab}
\end{align*}
\]
得到其中一个解,$m=\sqrt{q^{(l-1)}_a},n=0,g=\sqrt{q^{(l-1)}_b}C^{(l-1)}_{ab},p=\sqrt{q^{(l-1)}_b}\sqrt{1-{C^{(l-1)}_{ab}}^2}$,其中$C^{(l-1)}_{ab}=\frac{q^{(l-1)}_{ab}}{\sqrt{q^{(l-1)}_a q^{(l-1)}_b}}$.而且由于对称性，$q^{(l-1)}_a=q^{(l-1)}_b$,因此定义$q^{*}=q^{(l-1)}_a=q^{(l-1)}_b$
于是
\[
\begin{align*}
q^{(l)}_{ab}&=\sigma_{w}^2COV(a^{(l-1)}_a,a^{(l-1)}_b)+\sigma_{b}^2\\&=\sigma_{w}^2\int\int d(z_1)d(z_2)\phi(u_1) \phi(u_2)dz_2dz_1+\sigma_{b}^2 \tag{4}
\end{align*}
\]
其中$d(z_1)=\frac{1}{\sqrt{2\pi}}e^{-{z_1}^2},d(z_2)=\frac{1}{\sqrt{2\pi}}e^{-{z_2}^2},u_1=\sqrt{q^{(l-1)}_a}z_1,u_2=\sqrt{q^{(l-1)}_b}C^{(l-1)}_{ab}z_1+\sqrt{q^{(l-1)}_b}\sqrt{1-{C^{(l-1)}_{ab}}^2}z_2$。特别地,有
\[
C^{(l)}_{ab}=\frac{1}{q^{*}}(\sigma_{w}^2\int\int d(z_1)d(z_2)\phi(u_1) \phi(u_2)dz_2dz_1+\sigma_{b}^2) \tag{5}
\]
并记$\chi=\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}}$
可见无论是方差还是协方差都可以写成按层迭代的形式。那么什么时候会形成一个不动点呢？对于固定的$\sigma_{w},\sigma_{b}$来说，当$\frac{\partial{q^{(l)}}}{\partial{q^{(l-1)}}}<1$时,$q_{(l)}$稳定收敛;类似地,当$\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}}<1$时,$C^{(l)}_{ab}$收敛。那么收敛到什么值呢？实验表明，通过调节$\sigma_{w},\sigma_{b}$的大小可以影响收敛的情况。根据实验,当$\sigma_{b}=0.3$时,$\sigma_{w}$越大，$q^{(l)}$收敛的值越小。而且收敛速度很快，一般到达3层到4层就可以达到收敛，所以研究协方差的迭代过程中直接就把方差看成是收敛的。而$C^{(l)}_{ab}$的收敛速度慢一些,一般到达20层左右收敛，而且总有不动点$c^{*}=1$.固定$\sigma_{b}$并不断地增加$\sigma_{w}$，当$\sigma_{w}$较小，$c^{*}=1$是唯一的不动点，且$\chi_1<1$，不动点处于稳定状态；当$\sigma_{w}$继续增大，$\chi_1$也增大而且逐渐超过1,$c^{*}=1$的不动点不再稳定，开始产生新的不动点，随着$\sigma_{w}$增大,新的不动点就越小，最极端的时候达到0。当$\chi_1>1$时，网络倾向于对输入空间的流形进行折叠，呈现出混沌的状态；当$\chi_1<1$时，网络对输入空间的流形进行捋顺，呈现出平整有序的状态。$\chi_1(\sigma_{w},\sigma_{b})=1$构成了在第一个不动点时的一个分割混沌和有序的边界。
上边的分析都是针对固定的$\sigma_{w}$和$\sigma_{b}$，这只能说明不同的权值配置下网络的效果的差异，和dynamical isometry有什么关系呢。dynamical isometry主要影响的是网络训练时的速度。答案就与雅可比矩阵有关了。
雅可比矩阵为什么会影响网络训练的速度？先看神经网络任意两层间的雅可比矩阵的定义:
\[J=\frac{\partial{h^{(l)}}}{\partial{h^{(l-1)}}}=DW\]
具体意义就是，l-1层preactivation对l层的preactivation的贡献矩阵，其中D为对角矩阵，对角元素为激活函数在对应输入下的梯度，W为l层权重矩阵。
考虑误差传递,$\delta_l=\frac{\partial{L}}{\partial{h^{l}}}$，假设不同维度对误差的贡献互不相关且均值为0方差为1,有$COV(\delta_l,\delta_l)=I$,要使误差传递到第l-1层时对l-1层的分量依然保持类似的性质，必然有$\delta_{l-1}=J\delta_l,COV(\delta_{l-1},\delta_{l-1})=(J\delta_l)((J\delta_l)^T=JJ^T\approx I$,从而要求J的所有奇异值都接近O(1)。满足这个条件的J使得深层误差的corvariance较好地传递到浅层，保证了误差贡献的均衡性，从而网络学习得较快。而从正向传播来看，给第$l-1$层preactivation添加扰动u,也就是$h^{(l-1)}+u$传递到第l层就成了$h^{(l)}+Ju$(由泰勒一阶展开可得),$\frac{\|Ju\|^2}{\|u\|^2}$表征了扰动经过网络的缩张效应，从数学上来说这表征了J的平方奇异值谱的分布，分布越大越不接近O(1)缩张程度越大网络越不稳定,对于误差传递的情况也缓解了梯度的弥散或者爆炸。这就是dynamical isometry的本质。
进一步地，将这个任意相邻层的雅可比矩阵推广到最后一层与输入之间的雅可比矩阵：
\[J=\prod_{l=0}^{D-1}D^{(l)}W^{(l)}\]
同样要实现dynamical isometry也要求有J的奇异值接近O(1)。
然后当我们需要分析一个网络的dynamic isometry，就要分析对应的雅可比矩阵的奇异值的谱。这就和上边的迭代有一定的关系了，因为
\[
\begin{align*}
\chi_1&=\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}|_{1}}\\&=\frac{1}{N_l}tr(JJ^T) \tag{6}
\end{align*}
\]
上式是怎么得来的呢？联系（1）式和（3）式，（2）和（4）式之间的联系，有
\[
\begin{align*}
\frac{1}{N_l}tr(J^{(l)}{J^{(l)}}^T)&=VAR(W^{(l)}\phi^{'}(h^{(l-1)}))\\&=\sigma_{w}^2\int d(z){\phi^{'}(\sqrt{q^{(l-1)}}z)}^2dz\\&=\sigma_{w}^2\int d(z){\phi^{'}(\sqrt{q^{*}}z)}^2dz \tag{7}
\end{align*}
\]
而由（5）式,
\[
\begin{align*}
\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}}&=\frac{1}{q^*}\sigma_{w}^2\int\int d(z_1)d(z_2)\phi(u_1)\frac{\partial{\phi(u_2)}}{\partial{C^{(l-1)}_{ab}}}dz_1dz_2\\
&=\frac{1}{q^*}\sigma_{w}^2\int\int d(z_1)d(z_2)\phi(u_1)\phi^{'}(u_2)(\sqrt{q^{*}}z_1+\sqrt{q^{*}}\frac{-C^{(l-1)}_{ab}}{\sqrt{1-{C^{(l-1)}_{ab}}^2}}z_2)dz_1dz_2 \tag{8}\\
\end{align*}
\]
又容易有$\int_{-\infty}^{+\infty} d(z)F^{'}(z)dz=\int_{-\infty}^{+\infty} d(z)dF(z)=d(z)F(z)|_{-\infty}^{+\infty}-\int F(z)dd(z)=\int d(z) F(z)zdz$，所以式（8）变成
\[
\begin{align*}
\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}}&=\frac{1}{q^*}\sigma_{w}^2(\int d(z_2)(\int d(z_1)\phi(u_1)\phi^{'}(u_2)\sqrt{q^{*}}z_1dz_1)dz_2\\&+\int d(z_1)\phi(u_1)(\int d(z_2)\phi^{'}(u_2)\sqrt{q^{*}}\frac{-C^{(l-1)}_{ab}}{\sqrt{1-{C^{(l-1)}_{ab}}^2}}z_2dz_2)dz_1)\\
&=\frac{1}{q^*}\sigma_{w}^2(\int\int d(z_2)d(z_1)q^{*}\phi^{'}(u_1)\phi^{'}(u_2)dz_1dz_2+\int\int d(z_1)d(z_2)\phi(u_1)\phi^{''}(u_2)q^{*}C^{(l-1)}_{ab}dz_2dz_1\\
&-\int\int d(z_1)d(z_2)\phi(u_1)\phi^{''}(u_2)q^{*}C^{(l-1)}_{ab}dz_2dz_1)\\
&=\sigma_{w}^2\int\int d(z_2)d(z_1)\phi^{'}(u_1)\phi^{'}(u_2)dz_1dz_2
\end{align*}
\]
当$C^{(l-1)}_{ab}=1$,便有
\[\begin{align*}\frac{\partial{C^{(l)}_{ab}}}{\partial{C^{(l-1)}_{ab}}}&=\sigma_{w}^2\int d(z_1){\phi^{'}(\sqrt{q^{(l-1)}}z_1)}^2dz_1\\&=\sigma_{w}^2\int d(z_1){\phi^{'}(\sqrt{q^{*}}z_1)}^2dz_1 \tag{9}\end{align*}\]
（7）式右边等于（9）式右边,从而当$C^{(l-1)}_{ab}=c^{*}=1$，（6）式成立。
也就是说我们通过合理设置$(\sigma_{w}^2,\sigma_{b})$的配置（在混沌与有序的边界$\chi_1(\sigma_{w},\sigma_{b})=1$上寻找），使得$\chi_1=1$,便能使得雅可比矩阵的奇异值的平方均值为1，实现一定程度上的dynamical isometry，防止梯度弥散或者爆炸。
P.S：鉴于一开始网络至少在二十层以内$C^{(l)}_{ab}$都不会收敛，我们可以一开始使$q^{(0)}=q^{*},q^{(0)}_{ab}=q^{*}$，从而使得它一开始就在不动点。另外对于更深层网络，比如上百层的网络，这二十多层的收敛延迟可以说是小意思了，不影响整体的dynamical isometry。
本质上来说使得$\frac{1}{N_l}tr(JJ^T)=\chi_1=1$描述的是雅可比矩阵奇异值的方差(统计学第二动量)同时也是$JJ^T$的特征值的平均值(统计学第一动量)。这其实还不能很充分地证明雅可比矩阵的奇异值均匀地接近O(1)，因为存在$\lambda_{max}和\lambda_{min}$极大和极小的情况，大部分奇异值也许是被平均的。所以$\chi_1=1$只能是一个实现dynamical isometry的必要但不充分条件。
更进一步地说，要提高dynamical isometry的置信度，就要使用更高阶的统计学动量，比如$JJ^T$的特征值方差，以及最大特征值的范围。
[1]中引入了free probablity的相关理论，实现了对$JJ^T$的特征谱估计。下边介绍一下这种分析工具的数学基础：
定义矩阵特征值的谱密度:
\[\rho_{X}(\lambda)=\left \langle\frac{1}{N}\sum_{i=1}^{N}\delta(\lambda-\lambda_i) \right \rangle \tag{10}\]
表示对任意随机矩阵X的平均特征频谱密度。
定义$\rho_X$的Stieltjes变换：
\[G_X(z)=\int_{R}\frac{\rho_X(t)}{z-t}dt \tag{11}\]
以及逆变换:
\[\rho_X(\lambda)=-\frac{1}{\pi}lim_{\varepsilon\to 0^{+}} ImG_X(\lambda+i\varepsilon) \tag{12}\]
以及与$G_X(z)$相关的动量生成函数$M_X(z)$:
\[M_X(z)=zG_X(z)-1=\sum_{k=1}^{\infty}\frac{m_k}{z^k} \tag{13}\]
而且$m_k=\int d\lambda \rho_X(\lambda)\lambda^{k}=\frac{1}{N}<trX^{k}>_X$
以及S变换:
\[S_X(z)=\frac{1+z}{z{M^{-1}}_{X}(z)} \tag{14}\]
S变换满足一个特性:
\[S_{AB}(z)=S_{A}(z)S_{B}(z)\]
所以$S_{Jj^T}=\prod_{l=1}^{L}S_{W_{l}{W_{l}}^T}S_{{D_l}^2}=S_{WW^T}^LS_{D^2}^L$,这说明我们通过研究分别研究权重分布以及不同的激活函数，使用（11）,（13）,（14）式可以得到$S_{WW^T}$和$S_{D^2}$,再反过来一步步使用逆变换从$S_{JJ^T}$求得$M_{JJ^T}$,然后根据$M_{JJ^T}$得到高阶动量的估计（用到了Lagrange Inversion Theorem）。
论文中通过理论分析(推导略复杂，推导不来了，读者自己看论文吧)以及实验验证，得到了下边的结论：
1.线性神经网络中，使用正交权重初始化直接就能获得dynamical isometry,能够获得深度无关的学习时间;高斯初始化，使得$JJ^T$的$\lambda_{max}\sim eL$,第二动量$\sigma_{JJ^T}=L$,特征谱随网络深度而越来越扩散。
2.非线性网络中,Relu网络无论如何都不能获得dynamical isometry的特性，$\lambda_{max}$与$\sigma_{JJ^T}$与深度L呈正相关；tanh激活的网络，高斯初始化$\lambda_{max}$与$\sigma_{JJ^T}$与深度L呈正相关，正交权重初始化当合适选择$(\sigma_{w},b)$使得对应的$q^{*}$较小（隐含了$\chi_1=1$）的情况下(实验中的值是$(1.0,10^{-5})$)，$\lambda_{max}$与$\sigma_{JJ^T}$与深度L也成弱正相关，但是系数相当小，达到dynamical isometry。
3.非线性网络中因为$\lambda_{max}$在达到dynamical isometry情况下也与L成正相关,所以网络最佳学习率总是随着网络的深度增加而减少,$\eta_{opt} \sim O(\frac{1}{L})$,这和线性网络的情况一样。另外在最佳学习率的情况下的学习时间$t\sim O(\sqrt{L})$，受深度影响较小；而线性网络中，最优学习率训练时间是与深度无关的。dynamical isometry也使得网络可以容忍更大的学习率.
4.达到dynamical isometry后深度网络的训练算法RMSPROP，ADAM，SGD,SGD+MOTUMN等都获得了显著更好的结果。
5.合理选择$(\sigma_{w},\sigma_{b})$（使$q^{*}=0$的不一定最好）可以使得训练过程中dynamical isometry也能一定程度保持，而且达到更好的泛化效果。
[3]S.S. Schoenholz, J.Gilmer, S.Ganguli, J.S.Dickstein.Deep Information Propagation.arXiv:1611.01232