MATLAB Implementation: light-weight vs. heavy-weight

转载自https://dahuasky.wordpress.com/2008/11/07/matlab-implementation-light-weight-vs-heavy-weight/

这次，回到一个实际一点的问题，关于matlab的实现。

同一个数学问题，在实际计算中，往往是可能有多种途径的。虽然殊途同归，但是效率很可能大相径庭（即使这些不同途径在理论上有相同的复杂度）。对于小规模计算，这种差别也许无关重要，但是在大规模的simulation或者experiment中，注意不同方式的效率差异可能会让你在跑程序上花费的时间从一天变成一个小时。而且由于计算机字长和精度的限制，选择一个不恰当（但是理论仍旧是正确的）的途径，很可能导致数值溢出。

在这篇文章里，只是以高斯分布的概率密度计算为例子，说明不同计算途径可能导致的不同后果。

在实际问题中，有两种常见的使用高斯模型的情况：

1. 在一个高维空间，使用少数（几个或者几十个）高斯模型；
2. 在一个大型网络中（比如一个图像的每个像素是一个node），在每个节点使用一个（或几个）低维高斯模型。

这两种问题在matlab中的实现需要循不同的路径。

高维高斯模型的pdf计算

在log-scale计算pdf

我们接触许多重要分布（包括Gaussian），在数学上都属于一个类别：Exponential family。当我们需要计算这类分布的概率密度(pdf)时，直接计算概率密度是不明智的，尤其在高维空间。因为 log p(x) 往往会偏离零点甚远的距离，直接计算很可能会导致overflow或者underflow。对于，这类分布，计算log p(x)而不是p(x)本身，应该成为一种惯例。

很多时候，我们的最终目标不是p(x)或者log p(x)，而是由此得到后验概率q(x)。那么给出log likelihood之后，如何计算posteriori呢？最直接的方法是，用exp函数转为p(x)，除以这些p(x)的总和。对于单个样本，MATLAB的实现如下：

% Input loglik:  m x 1 vector with loglik(k) giving log_k p(x)% Input logpri:  m x 1 vector with logpri(k) giving log-prior of k-th model% Ouput q:       m x 1 vector with q(k) being the posteriori that the sample is from the k-th modelL = loglik + logpri;p = exp(L);q = p / sum(p);

通常在matlab实际应用中，我们不会只计算单个样本，而是大群样本一起计算。那么，就可以用下面的实现(注意bsxfun的运用，它出现在2007a，从此以后被大量运用于矩阵和向量之间的运算，比如把一个向量加到或者乘到矩阵中的每一行或者列，这是新版matlab中最标准和高效的写法)：

% Input loglik: m x n matrix, where loglik(k, i) gives log p_k(x_i),  % Input logpri: m x 1 column vector, where logpri(k) gives the log-prior of the k-th model% Output q:     m x n matrix, where q(k, i) is the posteriori that the i-th sample is from the k-th modelL = bsxfun(@plus, loglik, logpri);      % L <- log(prior) + log(likelihood)p = exp(L);                        % p <- prior * likelihoodq = bsxfun(@times, p, 1 ./ sum(p, 1));  % posterior: q(k,i) = p(k,i) / sum_{l=1}^m p(l,i) 

这在理论上是正确的，但是exp(L)这一步很可能会导致溢出，这也是我们用log-scale的原因。要解决这个问题，关键在于posteriori只取决于p的数值之间的相对比例，而无关于它们的绝对数值。因此，我们可以给各个p同时放大或者缩小到一定的水平，然后再进行计算。对p同时乘以一个因子，相当于给它们的对数同时加上一个常数。这里，我的策略是，给L的每列同时加上一个数，使得它的最大值恰好是零，然后在求posteriori。

L = bsxfun(@plus, loglik, logpri);L = bsxfun(@minus, L, max(L, [], 1));   % shift each column of L to a safe levelp = exp(L);q = bsxfun(@times, p, 1 ./ sum(p, 1));

为什么是要控制最大值，而不是最小值呢。道理是这样的：如果最大值是1（i.e. 它的对数是0），那么最小值可能下溢出变成0。但这在posteriori是没关系的，一个模型它的posteriori是0，还是1e-500，都是一回事——这个样本不是从这个模型出来的，在后验概率中，我们只关心对样本有意义的模型。而如果，我们反过来，把最小值调到1，那么最大值可能出现overflow，在matlab里面就是表达成inf，那么后验概率就根本求不出来了。

计算Mahalanobis distance

高斯分布的log-pdf可以写成下面的形式

log p(x) = (1/2) * (d * log(2 * pi) + log det( sigma ) +  (x – mu)^T * inv(sigma) * (x – mu))

这里，d 是维数，mu是均值向量，sigma是协方差矩阵。这个式子有三项，d * log(2 * pi)是一个常数项，在matlab里面计算很直接，这里不讨论了。而后面两项在计算上都有值得注意的问题。

我们首先来看第三项（这一项在数学上叫Mahalanobis distance）：

(x – mu)^T * inv(sigma) * (x – mu)

对于单个样本，它的直接实现是：

(x - mu)' * inv(sigma) * (x - mu)

对于多个样本，用for-loop显然不是高效的方法。我们可以通过向量化写成：

% X: d x n matrix, with each column being a sample% D: the difference between the input samples and the mean% md: 1 x n vector of Mahalanobis distancesD = bsxfun(@minus, X, mu);md = sum(D .* (inv(sigma) * D), 1);

这里面最重要的一步是inv(sigma) * D的计算，这也是最花时间的地方。对于高维空间，inv(sigma)是很耗时的。这里面有多种策略可以运用于提高运算效率。

1. 如果这种log-pdf的计算要被调用很多次，那么在第一次计算前预先计算好inv(sigma)存放起来，可以有效减少后续计算所需要的时间；
2. 如果只是对整批样本计算一次，那么在MATLAB里面更高效简洁的写法是sigma D。一些刚接触matlab的朋友可能没注意到，matlab里面对于inv(A) * B有着一种功能上等效，但是性能上更好的写法：A B。
3. 如果sigma具有某种特殊结构，那么计算时可以充分利用来提高效率。比如sigma是对角阵的话，那么inv(sigma) * D可以这样计算：

   bsxfun(@times, D, 1 ./ diag(sigma));

   如果一开始就把sigma存成d x 1的向量，只存放对角线上的值，那么diag(sigma)也省了，而且存储的空间复杂度也从O(d^2)降到O(d)。

对于高维空间，还需要注意的问题是奇异问题。在实际问题中，sigma可能是从样本中估计出来的，在样本量有限的情况下，很可能出现sigma是singular的情况，或者poor condition。这时就需要考虑regularize或者dimension reduction，先把sigma变成non-singular的，然后在用上述实现去求log-pdf。

计算log (det (sigma))

在log-pdf中还有一项是log(det(sigma))。这句话可以直接写在matlab里面，理论上是正确的，但是你很可能得到的结果是inf或者-inf，尤其是维数上千的情况。这里面的问题是，一个大矩阵的determinant可能是很大或者很小，以至overflow或者underflow，所以我们只能计算它的对数，而计算路径又不能经过它本身。方法有很多，基于线性代数的原理，我们可以有以下几种策略：

1. 对于一个矩阵，它的determinant就是它的所有eigenvalue的乘积，那么log-determinant就是所有eigenvalue的对数和。给予这个原理，我们可以给予特征值分解实现：

   eigvals = eig(sigma);v = sum(log(eigvals));

   类似的，也可以使用奇异值分解SVD。

2. 不过特征值分解或者奇异值分解的复杂度相当高。一种更高效的思路是使用下面的定理: 三角阵的determinant等于对角线上元素的乘积；两个矩阵乘积的determinant等于它们determinant的积。另外，任何一个正定阵，都可以分解为两个三角阵的乘积，这就是矩阵理论中的Cholesky分解。在MATLAB提供了Cholesky分解的非常高效的实现。基于Cholesky factorization，求log-determinant的过程可以表示为

   v = 2 * sum(log(diag(chol(sigma))));

   这是一种安全高效的实现。即使在det(sigma)保证不溢出的情款下，基于chol分解的实现也不比log(det(sigma))来得慢，比起基于eignevalue的实现则快了很多。对于一般的矩阵（可能非对称），类似的，可以使用LU分解或者QR分解。

大量低维高斯模型的pdf计算

在很多应用里面（比如image modeling)，一种通常的做法是，有大量的低维模型。这里，我们讨论二维的情况。这时我们会有很多不同的covariance matrix，它们都很小。而matlab并没有提供内建的函数，同时对大批矩阵求inverse，求cholesky factorization，或者求determinant。

但是，低维的好处是，它的矩阵的inverse, determinant都有很简单的解析表达式。可以利用它们进行计算：

求一批2 x 2矩阵的determinant

我们知道，一个2 x 2矩阵C的determinant的表达式是：C(1, 1) * C(2, 2) – C(1, 2) * C(2, 1)。 基于这点，我们可以对大批矩阵的求determinant，进行向量化：

% Input S: 2 x 2 x m array, with each page giving a covariance% Output detv: 1 x m array, with v(i) = det(S(:,:,i))mS = reshape(S, 4, m);detv = mS(1, :) .* mS(4, :) - mS(2, :) .* mS(3, :);

这种做法，比起用for-loop对每个矩阵逐一计算快很多。

对一批矩阵求inverse

同样的，我们可以给予解析表达式进行计算

invS = [mS(4,:); -mS(2,:); -mS(3,:); mS(1,:)] ./ detv;invS = reshape(invS, 2, 2, m);

我们看到inv要用到detv的值，所以在实际计算中，应该安排先计算determinant，然后再利用detv去计算inverse。举一反三，对于整个pdf的其它步骤，也可以类似的充分利用解析表达式直接向量化。

对于二维和三维，同过解析表达式向量化是可能的，但是四维以上，虽然仍旧是低维空间，但是解析表达就变得愈发困难。但是，仍旧可以通过各种方式来提高整体的计算效率。而在这里面，采用Expoential family representation比起采用传统的基于mu, sigma的表达更加高效，可以有效实现从estimation到evaluation的全面向量化，并且少数高维模型和大量低维模型为了适应效率要求而产生的在implementation上的差别，可以被重新弥合。关于这个内容，暂时不探讨了。

总的来说，和C/C++不同，MATLAB的实现和数学从来都是密不可分的，高效的实现离不开对数学和数值计算过程的理解。把数学思路充分体现于实现过程中，这是我一直以来特别喜欢matlab的一个重要原因。而反过来，解决计算过程中的效率苦难，也是众多计算科学中的活跃课题。无论是Machine learning，还是optimization，大量的努力其实也就是为了降低一些问题的计算复杂度和提高逼近精度，或者把intractable的问题变得tractable（至少在approximate的意义下）。