图信号处理学习笔记（3）：基于GMRF的图估计

上一篇文章中提及了图信号插值的一种方法，这是基于图已给定的情况。在实际的推荐系统等应用中，不仅要针对未知节点进行预测，还要获得一个良好的，切合训练数据的graph。

Reference：
[1] E. Pavez and A. Ortega, “Generalized Laplacian precision matrix estimation for graph signal processing,” 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Shanghai, 2016, pp. 6350-6354.
[2] H. E. Egilmez, E. Pavez and A. Ortega, “Graph Learning From Data Under Laplacian and Structural Constraints,” in IEEE Journal of Selected Topics in Signal Processing, vol. 11, no. 6, pp. 825-841, Sept. 2017.
[3] C. Zhang, D. Florencio, and P. A. Chou, “Graph signal processing—A probabilistic framework,” Microsoft Research, WA, USA, Tech. Rep. MSR-TR-2015-31, 2015.

一、优化问题概述

本质上，图学习可以归结为一个矩阵优化求解问题：给定节点集合$V$和对应的信号集合（训练数据集），估计存在的节点连接$E$及其权重$w$，或者换句话说，估计其图描述矩阵$Q$。在文[1]中，以任意的主对角线为正的对称矩阵（例如协方差矩阵$K$），来输出一个泛化拉普拉斯矩阵（Generalized Laplacian，简称为GL）$Q$。GL矩阵本身只需要满足对称且非对角线的元素为负数或0，对主对角线的元素没有要求。因此，在笔记（1）中提到的拉普拉斯矩阵，标准化拉普拉斯矩阵等，都是GL矩阵的特例。

所以，事实上GL矩阵可以由一个满足拉普拉斯矩阵特性的矩阵和一个对角矩阵相加。前者表达了节点与节点之间的关系，后者表达了单个节点自身的隐藏信息。

GL矩阵在这个问题中相比Laplacian和Normalized Laplacian有如下优势：1. 减少约束数量，便于优化。2. 使输出的矩阵更加泛化，能够包含更多节点自身的信息。3. 相比于普通矩阵，GL矩阵依然能满足节点域定理（Nodal Domain Theorem），使得该矩阵能够有实际意义上的频谱分解。节点域定理简单地说，就是能够保证在特征值变大的情况下，特征向量的过零点变多。（也就是实际意义上的高低频的区分）

先来讨论怎样的GL矩阵能够满足要求。以推荐系统为例，我们需要通过图的连接来表述各个用户间的相似程度。两个用户间的连接权值越大，他们的偏好越相似，则他们对特定事物的爱好程度也越相似。由此可见，最后生成的图，要做到在连接权值大的两个节点之间，对数据集达到平滑（Smoothness） 的效果。举例：用户A对物品1的偏好为5，用户B对物品1的偏好为4，用户C的偏好为3，显然A与B之间的相似度要更大，即连接权值要更大，此时可以发现，从偏好5到偏好4没有急剧的变化，这就是达到了平滑的效果。

因此，优化函数的建立，就是要以找到一个使数据集在上面的表现尽可能平滑的图为目标。在笔记（1）中提到，对于单个图信号，要度量其在图上的平滑程度，可以通过瑞利熵（Rayleigh Quotient） 来分析。假定图的Laplacian矩阵为$\mathbf{L}$，则图信号${\mathbf{x}}_i$在$\mathbf{L}$中的平滑程度为：
$${\sigma }_i={\mathbf{x}}_i^T\mathbf{L}{\mathbf{x}}_i$$

该式的值的范围在$[{\lambda }_{min},{\lambda }_{max}]$之间，即最小特征值与最大特征值之间。因此，当${\sigma }_i$越小，则图信号${\mathbf{x}}_i$在$\mathbf{L}$上越平滑。

那么很明显地，要寻找最优的$\mathbf{L}$在满足GL矩阵的约束的条件下，使得目标函数$f$最小化。$f$可定义如下：
$$f=\sum _{i=1}^N{\sigma }_i=\sum _{i=1}^N{\mathbf{x}}_i^T\mathbf{L}{\mathbf{x}}_i=tr({\mathbf{X}}^T\mathbf{L}\mathbf{X})=tr(\mathbf{L}\mathbf{X}{\mathbf{X}}^T)=Ntr(\mathbf{L}\mathbf{K})$$

由上式，$\mathbf{K}=\mathbf{X}{\mathbf{X}}^T/N$即为数据集$\mathbf{X}$在零均值条件下的协方差矩阵。

本文中，整个算法可以归结为一个基于高斯马尔科夫随机场（Gaussian Markov Random Field，GMRF） 的优化问题：一个随机向量$\mathbf{x}=(x_1,...,x_n{)}^T$被称为关于图 G = ( V = { 1 , . . . , n } , E ) \mathcal G=(\mathcal V=\{1,...,n\},\mathcal E) G=(V={ 1,...,n},E)的GMRF，假定均值向量为$\mu$，精度矩阵（协方差矩阵$\mathbf{K}$的逆矩阵，满足半正定特性）为$\mathbf{Q}$，则随机向量$\mathbf{x}$的概率密度函数为：
$$\mathbf{p}(\mathbf{x}|\mathbf{Q})=\frac{|\det (\mathbf{Q}){|}^{1/2}}{(2\pi {)}^{N/2}}\exp (-\frac{1}{2}({\mathbf{x}}^T-\mu )\mathbf{Q}(\mathbf{x}-\mu ))$$
并且${\mathbf{Q}}_{ij}$仅在节点$i$和节点$j$的边缘存在时（即$E_{ij}\in \mathcal{E}$）才不为零。

事实上，对于无向图$\mathcal{G}$，令其邻接矩阵为$\mathcal{A}\ge 0$，在高斯随机场的前提条件下，我们可以用泛化的精度矩阵$\mathbf{Q}$来表示图的GL矩阵$\mathcal{L}$。这是因为$\mathcal{Q}$可以与图的邻接矩阵$\mathcal{A}$形成一一映射的关系。或者说，一个半正定矩阵可以和一个非负对称矩阵形成一一映射的关系。 其映射方式如下：

定义一个基于$n\ast n$的矩阵$\mathbf{W}$和$n\ast n$的矩阵$\mathbf{Q}$的映射，使得：
$$\begin{gathered} {\mathbf{W}}_{ij}=-{\mathbf{Q}}_{ij}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}i̸=j \\ {\mathbf{W}}_{ii}= \end{gathered}$$