统计学习 | 矩阵正态分布 (matrix normal distribution)

前言:在机器学习和统计学习中,正态分布的身影无处不在,最为常见的是标准正态分布和多元正态分布 (multivariate normal distribution),两者分别作用于标量 (scalar) 和向量 (vector)。实际上,也存在一种正态分布的形式,它作用于矩阵,并广泛地应用于贝叶斯向量自回归模型 (Bayesian vector autoregression) 中。本文接下来将从大家所熟知的正态分布出发,先介绍矩阵正态分布,然后讨论矩阵正态分布在贝叶斯方法中的应用。

1 从标准正态分布到矩阵正态分布

1.1 标准正态分布

正态分布在机器学习和统计学习中随处可见,对于单一的随机变量 x\in\mathbb{R} ,其正态分布的形式为

\mathcal{N}\left(x\mid \mu,\sigma\right)=\left(2\pi\right)^{-1/2}\sigma^{-1}\exp\left(-\frac{1}{2}\sigma^{-2}\left(x-\mu\right)^2\right)

其中, \mu 表示均值, \sigma 表示标准差 ( \sigma^2 表示方差)。当然,标准正态分布在我们中学时代就已经接触了。

1.2 多元正态分布

随着我们学习线性代数、概率论等相关课程,我们又认识到:若向量 \boldsymbol{x}\in\mathbb{R}^{m} 服从正态分布,则存在一个多元正态分布 (multivariate normal distribution),形式为

\mathcal{N}\left(\boldsymbol{x}\mid \boldsymbol{\mu},\Sigma\right)=\left(2\pi\right)^{-m/2}\left|\Sigma\right|^{-1/2}\exp\left(-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^\top\Sigma^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right)

=\left(2\pi\right)^{-m/2}\left|\Sigma\right|^{-1/2}\exp\left(-\frac{1}{2}\text{tr}\left[\left(\boldsymbol{x}-\boldsymbol{\mu}\right)\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^\top\Sigma^{-1}\right]\right)

其中, \boldsymbol{\mu}\in\mathbb{R}^{m} 对应着正态分布的均值, \Sigma 则表示协方差矩阵。

需要说明的是,这里将多元正态分布的指数项写成矩阵迹 (trace) 的形式是为了方面后续认识矩阵正态分布,其中,在多元正态分布的写法中, \left(\boldsymbol{x}-\boldsymbol{\mu}\right)^\top\Sigma^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)=\text{tr}\left[\left(\boldsymbol{x}-\boldsymbol{\mu}\right)\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^\top\Sigma^{-1}\right] 是恒成立的。

1.3 矩阵正态分布

除了作用于向量的正态分布,实际上还存在一种正态分布,它作用于矩阵,常见于贝叶斯向量自回归模型 (Bayesian vector autoregression model)。一般而言,对于随机矩阵 X\in\mathbb{R}^{m\times n} ,若其服从矩阵正态分布,则形式为

\mathcal{MN}_{m\times n}\left(X\mid M,U,V\right)

=(2 \pi)^{-mn / 2}|{V}|^{-m / 2}|{U}|^{-n / 2}\exp \left(-\frac{1}{2} \operatorname{tr}\left[{V}^{-1}({X}-{M})^{\top} {U}^{-1}({X}-{M})\right]\right)

其中,符号 \mathcal{MN}_{m\times n}\left(\cdot\right) 来自于矩阵正态分布 (matrix normal distribution) 英文首字母的简写,下标指代随机矩阵的大小;矩阵 M\in\mathbb{R}^{m\times n} ,与随机矩阵 X 大小相同,对应于均值项;矩阵 U\in\mathbb{R}^{m\times m}V\in\mathbb{R}^{n\times n} 对应于协方差矩阵。

到这里,大家可能会发现:矩阵正态分布与多元正态分布在概率密度函数形式上存在较大差异。但实际上,我们依然可以从多元正态分布的角度理解矩阵正态分布,只要将矩阵正态分布进行向量化处理便可以得到多元正态分布形式,即

X\sim \mathcal{MN}_{m\times n}\left( M,U,V\right)\Leftrightarrow \text{vec}\left(X\right)\sim\mathcal{N}\left(\text{vec}\left(M\right),V\otimes U\right)

这两者完全等价 (证明过程可以参考wiki),公式中的符号 \otimes 表示Kronecker积 (这个运算规则可参考外积、Kronecker积和张量积一文,本文不做详述), V\otimes U\in\mathbb{R}^{(mn)\times (mn)} 表示多元正态分布的协方差矩阵;符号 \text{vec}\left(\cdot\right) 表示将给定矩阵按列组织成一个向量。

不过,这里有一个值得思考的问题:既然矩阵正态分布和多元正态分布可以等价,那为何还要构造出如此复杂的矩阵正态分布呢?

可能的原因要归结到协方差矩阵 V\otimes U\in\mathbb{R}^{(mn)\times (mn)} 上。实际上,假设我们想生成一个大小为 100\times 200 的随机矩阵 X ,并要求矩阵 X 在概率上服从矩阵正态分布。此时,若利用多元正态分布进行生成,则需要协方差矩阵 V\otimes U 的大小为 (100\times 200)\times (100\times 200)=20000\times 20000 ,元素数量为 4\times 10^8 ,显然,这个数字很惊人,毕竟存储这么大的矩阵就需要消耗计算机比较多的内存了。

因此,在很多时候,将矩阵正态分布用等价的多元正态分布替换有点不切实际。

2 贝叶斯向量自回归模型

2.1 向量自回归模型的数学表达式

在多元时间序列 (multivariate time series) 预测问题中,如果我们想借助回归模型对未来的序列进行合理估计,可采用简单的向量回归模型:

\boldsymbol{y}_{t}=A_{1}\boldsymbol{y}_{t-1}+\cdots+A_{d}\boldsymbol{y}_{t-d}+\boldsymbol{\epsilon}_{t},~t=d+1,...,T

其中, \boldsymbol{y}_{t}\in\mathbb{R}^{m} ​表示​时刻的时间序列变量,是一个大小为​ m 的向量;​ A_{1},...,A_{d} 表示大小 m\times m 为​的自回归系数矩阵;向量 \boldsymbol{\epsilon}_{t}\in\mathbb{R}^{m} 表示误差,常采用高斯噪声进行刻画。

一般而言,这条表达式可以被称为 \text{VAR}(d) ​,其中,​ d 表示向量自回归模型的阶数。现在,假设​ m=5,~d=4 ,则我们需要估计的参数为 \left\{A_1,A_2,A_3,A_4\right\} ​,这些参数所包含的参数数量为​ 5^2\times 4=100 .

为方便书写,向量自回归模型也可以写成如下形式:

\boldsymbol{y}_{t}=A^\top{\boldsymbol{z}}_{t}+\boldsymbol{\epsilon}_{t},~t=d+1,...,T

其中, A=[A_1,...,A_d]^\top\in\mathbb{R}^{(md)\times m} 是由系数矩阵构成;向量 {\boldsymbol{z}}_{t}=\left[\begin{array}{c} \boldsymbol{y}_{t-1} \\ \vdots\\ \boldsymbol{y}_{t-d} \\ \end{array}\right]\in\mathbb{R}^{(md)\times 1}

这条公式也可以进一步写成如下形式:

Y\approx {Z}A

其中,矩阵 Y=\left[\begin{array}{c} \boldsymbol{y}_{d+1}^\top \\ \vdots \\ \boldsymbol{y}_{T}^\top \\ \end{array}\right]\in\mathbb{R}^{(T-d)\times m} ;矩阵 {Z}=\left[\begin{array}{c} {\boldsymbol{z}}_{d+1}^\top \\ \vdots \\ {\boldsymbol{z}}_{T}^\top \\ \end{array}\right]\in\mathbb{R}^{(T-d)\times (md)} .

2.2 构建贝叶斯模型

假设向量自回归模型中的误差项 \boldsymbol{\epsilon}_{t} 为高斯噪声,即时间序列的观测值服从正态分布,则可以构建如下的贝叶斯向量自回归模型:

\begin{aligned} \boldsymbol{y}_{t}&\sim\mathcal{N}\left(A^\top{\boldsymbol{z}}_{t},\Sigma\right),~t=d+1,...,T \\ A&\sim\mathcal{MN}_{(md)\times m}\left(M_0,\Psi_0,\Sigma\right) \\ \Sigma &\sim\mathcal{IW}(S_0,\nu_0) \end{aligned}

其中,符号 \mathcal{IW}\left(\cdot\right) 表示inverse Wishart分布, \mathcal{IW}(\Sigma\mid S_0,\nu_0) 表达式为

\mathcal{IW}(\Sigma\mid S_0,\nu_0)\propto\left|\Sigma\right|^{-(\nu_0+m+1)/2}\exp\left(-\frac{1}{2}\text{tr}\left(S_0\Sigma^{-1}\right)\right)

2.3 推导后验分布

构建贝叶斯模型之后,我们需要对贝叶斯模型进行后验估计,简而言之,依据贝叶斯准则(后验分布正比于似然函数和先验分布的乘积,即 \text{posterior}\propto \text{likelihood}\times \text{prior} ),我们需要写出后验分布的具体表达式。

  • (1) 似然函数

在上述贝叶斯向量自回归模型中,似然来自于时间序列的观测值 \boldsymbol{y}_{t}\in\mathbb{R}^{m},~t=d+1,...,T ,即

\mathcal{L}\left(Y\mid A,\Sigma\right)=\left(2\pi\right)^{-m(T-d)/2}\left|\Sigma\right|^{-(T-d)/2}

\times \exp\left(-\frac{1}{2}\sum_{t=d+1}^{T}\left(\boldsymbol{y}_{t}-A^\top {\boldsymbol{z}}_{t}\right)^\top\Sigma^{-1}\left(\boldsymbol{y}_{t}-A^\top {\boldsymbol{z}}_{t}\right)\right)

其中,指数项中的求和部分可以被改写成如下形式:

\sum_{t=d+1}^{T}\left(\boldsymbol{y}_{t}-A^\top {\boldsymbol{z}}_{t}\right)^\top\Sigma^{-1}\left(\boldsymbol{y}_{t}-A^\top {\boldsymbol{z}}_{t}\right)

=\text{tr}\left[\left(Y-{Z}A\right)\Sigma^{-1}\left(Y-{Z}A\right)^\top\right]=\text{tr}\left[\Sigma^{-1}\left(Y-{Z}A\right)^\top\left(Y-{Z}A\right)\right]

综上所述,似然函数为

\mathcal{L}\left(Y\mid A,\Sigma\right)=\left(2\pi\right)^{-m(T-d)/2}\left|\Sigma\right|^{-(T-d)/2}

\times \exp\left(-\frac{1}{2}\text{tr}\left[\Sigma^{-1}\left(Y-{Z}A\right)^\top\left(Y-{Z}A\right)\right]\right) .

  • (2) 先验分布

\Sigma \sim\mathcal{IW}(S_0,\nu_0) 得,

\mathcal{IW}(\Sigma\mid S_0,\nu_0)\propto\left|\Sigma\right|^{-(\nu_0+m+1)/2}\exp\left(-\frac{1}{2}\text{tr}\left(S_0\Sigma^{-1}\right)\right)

A\sim\mathcal{MN}_{(md)\times m}\left(M_0,\Psi_0,\Sigma\right) 得,

\mathcal{MN}_{(md)\times m}\left(A\mid M_0,\Psi_0,\Sigma\right)

\propto\left|\Sigma\right|^{-md/2}\exp\left(-\frac{1}{2}\text{tr}\left[\Sigma^{-1}\left(A-M_0\right)^\top\Psi_0^{-1}\left(A-M_0\right)\right]\right)

(3) 后验分布

依据贝叶斯准则,从

\begin{aligned} p\left(A,\Sigma\mid Y\right)\propto &\mathcal{L}\left(Y\mid A,\Sigma\right)\times \mathcal{IW}\left(\Sigma\mid S_0,\nu_0\right) \\ &\times \mathcal{MN}_{(md)\times m}\left(A\mid M_0,\Psi_0,\Sigma\right) \end{aligned}

中可以推导出关于 A\in\mathbb{R}^{(md)\times m},\Sigma\in\mathbb{R}^{m\times m} 的后验分布为

\begin{aligned} A \sim\mathcal{MN}_{(md)\times m}\left(M^*,\Psi^*,\Sigma\right), ~~~~\Sigma \sim\mathcal{IW}\left(S^*,\nu^*\right) \\ \end{aligned}

后验分布中的参数为

\begin{cases} \Psi^*=\left(\Psi_0^{-1}+{Z}^\top{Z}\right)^{-1}, \\ M^*= \Psi^*\left(\Psi_0^{-1}M_0+{Z}^\top Y\right), \\ S^*=S_0+YY^\top+M_0^\top\Psi_0^{-1}M_0-\left(M^{*}\right)^\top\left(\Psi^{*}\right)^{-1}M^{*}, \\ \nu^*=\nu_0+T-d. \end{cases}

3 如何生成服从矩阵正态分布的随机矩阵?

现假设矩阵 X\in\mathbb{R}^{m\times n} 服从矩阵正态分布 \mathcal{MN}_{m\times n}\left(M,U,V\right) ,若已知 M\in\mathbb{R}^{m\times n},~U\in\mathbb{R}^{m\times m},~V\in\mathbb{R}^{n\times n} ,则生成随机矩阵 X 的原理如下:

第一步:对矩阵 U,V 作Cholesky分解,即 U=PP^\top,~V=QQ^\top ,分别得到矩阵 PQ
第二步:用标准正态分布生成一个大小为 m\times n 的矩阵 X_0
第三步:计算 X=M+PX_0Q^\top ,其中, X\sim \mathcal{MN}_{m\times n}\left(M,U,V\right) .

在Python中,可以利用numpy写出随机矩阵 X 的生成函数:

def mnrnd(M, U, V):
    """
    Generate matrix normal distributed random matrix.
    M is a m-by-n matrix, U is a m-by-m matrix, and V is a n-by-n matrix.
    """
    import numpy as np
    
    X0 = np.random.rand(M.shape)
    P = np.linalg.cholesky(U)
    Q = np.linalg.cholesky(V)
    return M + np.matmul(np.matmul(P, X0), Q.T)

4 参考

[1] matrix normal distribution - Wiki

[2] slide | Bayesian Vector Autoregressions

[3] Forecasting with Bayesian Vector Autoregressions

[4] Bayesian Autoregressive Time Series Models

[5] slide | Bayesian VARs

编辑于 2019-07-18

文章被以下专栏收录