【转载请注明出处】http://www.cnblogs.com/mashiqi
通过杨老师的细心点拨，我获得了以前没有的看问题的视角，重新审视了以前看的文章。通过这些新的视角，我从这些内容中看到了以前没有注意到的东西，看到了更本质的东西。希望多多交流！
$ extrm{JSE in my view:}$
杨老师已经把JSE和EM讲的很明白了，现在把它们都写下来，免得以后又忘了。在JSE中，$mu$从fixed number变成random variable，从而引入了本来没有的uncertainty。经过积分$int p (z|mu) cdot p(mu; heta) dmu$，使得这种多出来的uncertainty被消除掉，同时引入一个新的麻烦：参数$ heta$。但是我们要注意到的是，新麻烦“$ heta$”可比就麻烦“$mu$的uncertainty”要好处理得多，这或许就是statistics的作用吧。现在$mu$暂时消失了，我们要做的事就是根据observed data去估计这个新麻烦$ heta$。在JSE中是用了一个统计量$Q = (1+sigma^2) / sum_i z_i^2$来进行估计的，这个统计量$Q$符合inverse-Chi distribution。但我觉得也可以使用其他的方法进行估计。
因为$z_i sim N(0,1+sigma^2)$，所以也可以用如下式子来估计$sigma^2$：$$1+sigma^2 = frac{1}{N-1} sum_{i=1}^{N} (z_i - ar{z})^2,~ar{z} = frac{1}{N} sum_{i=1}^{N} z_i$$但是，经过Matlab的比较后发现，JSE的performance要优于采用样本方差来估计的方法。下面是Matlab代码：
function JSE_sigma_compare(N) if nargin < 1 N = 20; end rand('state',1000); % rng(1000,'v5uniform'); randn('state',1000); % rng(1000,'v5normal'); Times = 1000; mu_MLE_err = zeros(Times,1); mu_JSE_err_type1 = zeros(Times,1); mu_JSE_err_type2 = zeros(Times,1); for j = 1:Times % In order to fully check the effectiveness, mu should be sampled from % various type of distributions, such as uniform or Gaussian % distribution. mu = 3 + randn(N,1); % mu = randn(N,2); z = mu + randn(N,1); % MLE mu_MLE = z; mu_MLE_err(j) = norm(mu_MLE - mu)^2/N; % JSE type1, estimate sigma^2 using inverse-Chi expectation temp = mean(z,2); mu_JSE_type1 = ( 1 - (N-2)/norm(z)^2 ) * z; mu_JSE_err_type1(j) = norm(mu_JSE_type1 - mu)^2/N; % JSE type2, estimate sigma^2 using sampling variance temp = norm(z-mean(z))^2; mu_JSE_type2 = ( 1 - (N-1)/temp ) * z; mu_JSE_err_type2(j) = norm(mu_JSE_type2 - mu)^2/N; end mu_MLE_err_mean = mean(mu_MLE_err); mu_JSE_err_type1_mean = mean(mu_JSE_err_type1); mu_JSE_err_type2_mean = mean(mu_JSE_err_type2); figure(1),boxplot([mu_JSE_err_type1 mu_JSE_err_type2]); str = strcat('N =',32,num2str(N),'.',' ');fprintf(str); str = strcat('The mean error of MLE, with N=',num2str(N),', is',32,num2str(mu_MLE_err_mean),' '); fprintf(str); str = strcat('The mean error of JSE_type1, with N=',num2str(N),', is',32,num2str(mu_JSE_err_type1_mean),' '); fprintf(str); str = strcat('The mean error of JSE_type2, with N=',num2str(N),', is',32,num2str(mu_JSE_err_type2_mean),' '); fprintf(str);
另外，可以比较一下JSE和加正则化化因子的方法：$$min_w |t-Phi w|_2^2 + lambda|w|_2^2$$这种加上因子$lambda|w|_2^2$的方法，究其本质，是在一个确定的高斯分布上选择合适的$w$的方法。因为系数$lambda$是固定的，或者是通过人为调整的，总之不是从数据中优化得来的，那么$w$的高斯分布就一定是确定的。这个方法可以理解为：从一个确定的random variable的分布里面，选择一个fixed number出来，使得这个fixed number能使目标函数达到最值(或极值)。这个方法和JSE看起来很像，然而，它们真的真的是有本质的区别的！！！JSE算法，是从一堆random variable($left{ X|X sim(0,sigma^2),sigma geq 0
ight}$)中，选出一个最适合的random variable，而不是最适合的fixed number，来使得目标函数达到最值，这就是区别。这个对random variable的选择是通过积分算子这个桥梁才得以实现的，也就是通过marginal实现。当然，在实际应用中，我们通常不是把这个最合适的random variable给出来最优解，因为它不是一个数嘛，而是给出它的期望值，作为最优解。
如果真正看明白了上面的例子，就能对random variable的理解加深一些。不光是JSE，其实很多其他的方法也是这样来处理随机变量的。
$ extrm{EM in my view:}$
EM算法通过增加隐变量，引入一些uncertainty，同时使问题得到简化。然后通过对$ln p(oldsymbol{x},oldsymbol{z}; heta)$求期望的方法将这个uncertainty转化为比较好处理的certainty，也就是确定的参数$ heta$。然后求$ heta$的值，使得$oldsymbol{E}_{oldsymbol{z}|oldsymbol{x}; heta_{old}} left[ ln p(oldsymbol{x},oldsymbol{z}; heta) ight]$达到极大值。不停的交替求期望和极大值。就是这样！
总结一下，JSE和EM大致都可以分为三个步骤：第一，引入一个新的uncertainty，从而使得问题变得更加make sense，更加直观，反正加入uncertainty能带来一些好处；第二，要通过一些方法消除这个uncertainty带来的byproduct，把这种byproduct再转化为一些certain的东西，JSE通过积分，EM通过求期望都达到了这样的目的；第三，去估计这个这个certain的东西，JSE使用统计量，EM使用Evidence Maximization达到这样的目的。在每一个步骤中使用不同的方法，就可以解决不同类型的问题！
$ extrm{Relevance Vector Machine:}$
RVM是从线性回归模型派生出来的，所以首次从线性回归模型的角度来看RVM。RVM的线性回归模型是：$$t_n = sum_{k} phi_k(x_n)w_n + epsilon_n$$，其中${x_n,t_n}$是input-target pairs，它们是已知的; $phi_k(cdot)$是basis function; $epsilon_n sim N(0,sigma^2)$; 而$oldsymbol{w}$和$sigma^2$就是需要被求出的参数。现在可以将$oldsymbol{w}$和$sigma^2$看作fixed number，然后通过MLE直接求出；也可以将他们看作random variable通过概率的方法求出它们的最优估计值。RVM是选择的后一种方式。
通过给每个$w_k$加上高斯分布：$w_k sim N(0,alpha_i^{-1})$，使得$oldsymbol{w}$成为random ariable。但是成为random variable之后，自然引入了uncertainty了，这时候不能像处理fixed number一样处理$oldsymbol{w}$了。这时候为了消除这个uncertainty，通常有几种手段，RVM的方法是对$oldsymbol{w}$做积分将其消除掉：$$p(oldsymbol{t};oldsymbol{alpha},sigma^2) = int_{oldsymbol{w}} p(oldsymbol{t}|oldsymbol{w};sigma^2) p(oldsymbol{w};oldsymbol{alpha}) doldsymbol{w}.$$到了这里，$oldsymbol{w}$暂时消失了，只剩下参数$oldsymbol{alpha}$和$sigma^2$有待求解了。求解这些参数的时候，可以像JSE那样用估计量的方法，也可以直接像MLE那样求最大似然值，方法多种多样了。
我们在JSE里面看到了对take home message的使用，其实在RVM中，也可以把其他样本的信息拿来大家分享，使得得到的结果更加精确一些。在RVM模型中，其实$alpha_k$也是有一个先验的，并且$alpha_k,~(k=1,cdots,K)$的先验拥有共同的参数$a,b$：$$alpha_k sim Gamma(alpha_k;a,b)$$，其中$Gamma(alpha;a,b) = Gamma(a)^{-1} b^a alpha^{a-1} e^{-balpha}$。可惜的是，RVM中将$a,b$的值设定成了固定值，如果$a,b$的值也能像JSE里面的$sigma^2$一样可以从数据中学习的话，那么这样是不是$alpha_k$的值会更好一点呢？进而使得算法在性能上提高一些？
$ extrm{Bayesian Compressive Sensing:}$
下面再回味“Bayesian Compressive Sensing”一下这篇文章，发现文章所描绘的模型一下子在我眼前清晰了许多。这篇文章是将RVM运用到CS上面，求解CS的问题。不是一般性，CS的模型可以写成$$t = Aw + e$$其中$x$是待恢复的数据，$A$包括观测矩阵和稀疏矩阵，$t$是观测到的值，$e$是误差或者理解成噪声。Compressive Sensing说：$$argmin_w left{ |t-Aw|_2^2 + ho|w|_0 ight} approx argmin_w left{ |t-Aw|_2^2 + ho|w|_1 ight}$$而$$argmin_w left{ |t-Aw|_2^2 + ho|w|_1 ight} Leftrightarrow argmax_w N(t|Aw,sigma^2) cdot e^{- ho|w|_1}$$于是这篇文章后面就开始了使用RVM模型来求解问题。但是Laplace先验$e^{- ho|w|_1}$又不是很好处理，于是这个Laplace先验就被作者换成了高斯先验+Gamma先验，这就和RVM完全一样的。到这里，这篇文章基本上就清楚了。剩下的RVM的内容就不重复赘述了。
$ extrm{Bayesian Compressive Sensing using Laplace Priors:}$
在这篇文章中，作者使用了多重先验，最开始我完全搞不懂这样做的意义，不过现在明白一些了，把它们写下来。第一层先验和第二层先验要结合起来一起看，因为它们结合起来构成了Laplace先验，这正是“Bayesian Compressive Sensing”这篇文章所没有完成的事，而在这篇文章中得到了解决。第一层是系数先验：$$p(w_i|gamma_i) = N(w_i|0,gamma_i)$$其中$w_i$是系数，具体可以参见上一个section。第二层是超参数$gamma_i$的先验：$$p(gamma_i|lambda) = frac
{lambda}{2} exp left( -frac{lambda gamma_i}{2}
ight),~gamma_i geq 0,lambda geq 0$$这两个先验合起来对$gamma_i$求一个积分，就可以变成Laplace先验：$$p(w_i|lambda) = int_{gamma_i} p(w_i|gamma_i) p(gamma_i|lambda) dgamma_i = frac{lambda}{2} exp left( -sqrt{lambda} |w_i|
ight)$$这就完成了对$| cdot |_1$的实现。
第三层先验是对超参数$lambda$的先验：$$p(lambda;v) = Gamma(lambda;v/2,v/2)$$根据文章的原话，这最后一层先验的由来是因为：“the last stage of $(18)$ is embedded to calculate $lambda$.”，其中这个公式$(18)$就是：$p(lambda;v) = Gamma(lambda;v/2,v/2)$。不知道我解读是否正确，加上这个先验的目的就是为了计算$lambda$，但其实不加这先验也可以计算$lambda$啊。这就像是JSE里面那样，$sigma^2$不需要先验分布也可以确定啊。如果加这个先验的原因是因为take home message的话，在$p(gamma_i|lambda),~(i=1,cdots,N)$的时候已经通过$lambda$达到了borrow from others的目的，那干嘛要再通过一个$v$来borrow from others？发现这些问题，或许是因为我还没看懂，或许是因为我已经看懂一些了，从而能从别人的paper中发现一些问题了。我希望是后者，不过如果是我没看懂，也没关系，在多看几次呗。
模型分析了之后，再来分析一下它的求解过程。文章中所优化的目标函数是：$$mathcal{L} = ln p(y,gamma,eta,lambda) = ln int_{w}
p(y|w,eta)p(w|gamma)p(gamma|lambda)p(lambda)p(eta) dw$$但是根据我前面的分析，我觉得是不是应该对如下这个函数进行优化呢？：$$mathcal{L} = ln p(y,eta,lambda) = ln int_{w} int_{gamma} p(y|w,eta)p(w|gamma)p(gamma|lambda)p(lambda) dgamma p(eta) dw$$这是因为$$int_{gamma} p(y|w,eta)p(w|gamma)p(gamma|lambda)p(lambda) dgamma Leftrightarrow N(y|w,eta) cdot p_{ extrm{Laplace}}(w)$$要对$gamma$进行积分才是真正实现了Laplace先验嘛！作者不这样做的原因，我想有可能是因为两个积分符号处理起来太复杂了吧。把对$gamma$的积分去掉这一点，让我联想到了EM算法中是如何处理uncertainty的，我来试着分析一下。由于对$gamma$的积分不好处理，但是又要cancel掉这个超参数，所以就取了一个能使$mathcal{L}$达到极值的$hat{gamma}$，一旦random variable的概率分布坍缩到某一个确定的数上面时，uncertainty就消失了，问题得到了简便。当然这样处理是一种近似处理。
文章后面还有不小的一部分是讲如何加速计算过程的，做的也很好，但是由于和本次总结的主要目标关系不大，所以略去不分析了。