高效的Matlab：mex & profile

刚一条句子让wordpress bug了，拆成两篇算了

for的终极取代：多线程混合编程

我们也常会遇到并不能简答并行化的情形，例如针对数组元素的复杂的操作。这时的一个选择就是混合编程: mex, 它可以实现matlab与C/C++/Fortan之间的相互调用。matlab自带的user guide里有详细介绍。

在这里我们讨论两个小trick. 首先注意到向量化的高斯核代码中我们用到了exp(A)，A可能是一个非常大的矩阵。而matlab却会不谙风情的先将这个大矩阵在内存里复制一次，然后再逐一算exp. 有时我们明明算好内存能存下，但matlab却报out of memory了。

我们可以用C来对内存进行直接操作。例如for(i...){A(i)=exp(A(i));}. 不过还有另外一个问题：现在我们计算机大多是多核多线程，matlab也很早就支持了多线程加速，而简单的for只会单线程，所以mex后的代码可能会慢很多。这时的我们的一个选择就是openmp，它为共享内存的多线程提供了非常方便的实现。

下面附上一计算exp(A)的多线程代码，它无需额外内存:

#include "mex.h"
#include
#include
#include
using namespace std;

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
  double *A, *pp;
  mwSize n, m, p, i, j;

  A = mxGetPr(prhs[0]);
  n = mxGetN(prhs[0]); // # of rows in C
  m = mxGetM(prhs[0]); // # of columns in C
  p = omp_get_num_procs();

#pragma omp parallel for private (j)    \
    default(shared) num_threads(p)
  for (i = 0; i < n; i ++) {
    for (j = 0; j < m; j ++)
      A[i*m+j] = exp(A[i*m+j]);
  }
  return;
}

编译时需要稍修改下编译参数来开启openmp的支持。例如linux下gcc需要加上flag-fopenmp，详见此说明，Windows下的visual studio的相应的flag是/openmp.

但用这段代码时需要很小心。matlab使用了lazy copy的技术，既当我们运行命令A=B时，matlab不会立即将B复制一份给A, 只有当A或者B被修改的时候才执行复制。如果我们用上面代码对B做指数运算，例如mex_exp(B), 这时matlab不知道B变了，所以复制还是不会执行，那么A就会和B一起被修改了，这也许会给后面的代码造成麻烦。

代码究竟慢在哪里

我们可以方便使用tic,toc来知道每段代码的执行时间。但要知道每条语句的详细运行情况的话，非profile莫属。

先运行profile on，再运行需要测试的代码，然后使用profile viewer来查看报告。例如，上面两段不同的计算高斯核的实现的profile报告如下（为了报告更有信息，将其中的一条长语句拆开成了三句）。

根据报告我们就知道每条语句的执行情况，从而可以对其中时间耗费最多的部分进行优化。

无止尽的优化

大多人类活动就是一优化过程。无论是金钱，利益，地位，还是算法，代码，编译，硬件。但除了最优化外，我们还追求一些其他。例如我们用matlab是为了实现和调试上的方便，虽然优化能让计算机做得更快，但是却花去了我们的工作时间，而且使得代码难读懂。糟糕的实现是大忌，过度的优化也不可取。所以，我们需要在计算时间，开发时间，代码的可读性之间做一个权衡。

高效的Matlab：向量化

前面一直在谈理论，为了避免成为纯理论blog, 这次聊聊最实验的部分：coding.

Matlab是数值计算和仿真中使用最为广泛的软件。它便捷，高效，通常可以直接将伪代码算法转换成matlab代码。但在方便的同时，matlab里面也有很多性能的陷阱。例如令人诟病的for语句。

在matlab里，至少有两种方法将一个数组每个元素加1. 一种来自一般的程序语言for i=1:length(a), a(i)=a(i)+1; end；另一种则借鉴数学记法a=a+1。

但在matlab中，这两实现的性能却有天壤之别。虽然matlab在不断优化for的执行，不过即使是matlab 2010b, 在我的笔记本上（Core 2 Duo 2.4 GHz, 4GB, Win7），对于一个10万元素的一维数组，第一种实现需5微秒，而第二种只用0.4微妙。后者快上10倍。

事实上，用过matlab就会发现，matlab自带的函数通常都有很高的性能。例如刚提到的a=a+1，或是排序 sort(a)，以及令人称道的矩阵乘法，SVD分解等等。事实上，对于两矩阵a和b相乘，如果简单用三层for循环来实现（用C），性能肯定会比matlab差很多。

虽然matlab的界面是用java实现以方便跨平台，在性能上不尽人意，不过内核是用C实现的。可以试试命令matlab -nojvm。对于矩阵运算，matlab是调用高度优化的BLAS/LAPACK库。例如BLAS level3提供的一般矩阵乘法，它将矩阵拆成很多小块，充分利用内存的连续性，cpu的cache架构，多线程等来优化执行。matlab做的只是提供一个方便的人机接口。当然，其他函数例如sort matlab则提供自己优化过的实现。

但为什么用for时就会慢呢。主要原因是matlab是一种脚本语言，通过解释执行。每次我们运行一段代码时，matlab先预解析一遍，生成p-code, 然后一句一句执行。预解析当然是会有代价。在调用for时，matlab会每循环一次就会重新解析下for循环内部的代码，所以自然就慢了。

提高性能的办法是尽量使用matlab提供的函数，而非自己写。例如all（是否元素全非0）, cumsum（累加）, diff（相邻元素的差）,prod（累乘）。另外一些小技巧包括预先通过zeros等操作分配内存。在这里我们着重讨论for的优化，通常称其为向量化（vectorizing）。

用bsxfun来取代for

bsxfun是一个matlab自版本R2007a来就提供的一个函数，作用是”applies an element-by-element binary operation to arrays a and b, with singleton expansion enabled.”

举个例子。假设我们有一列向量和一行向量。

a = randn(3,1), b = randn(1,3)

a =
   -0.2453
   -0.2766
   -0.1913

b =
    0.6062    0.5655    0.9057

我们可以很简单的使用matlab的外乘c=a*b来得到. 但如果我们想用”外加”呢？既. 这时我们可以用c=bsxfun(@plus,a,b)来实现。

bsxfun的执行是这样的，如果a和b的大小相同，那么c=a+b. 但如果有某维不同，且a或b必须有一个在这一维的维数为1, 那么bsxfun就将少的这个虚拟的复制一些来使与多的维数一样。在我们这里，b的第一维只有1（只一行），所以bsxfun将b复制3次形成一个3×3的矩阵，同样也将a复制成3×3的矩阵。这个等价于c=repmat(a,1,3)+repmat(b,3,1)。这里

repmat(a,1,3)

ans =
   -0.2453   -0.2453   -0.2453
   -0.2766   -0.2766   -0.2766
   -0.1913   -0.1913   -0.1913

repmat是显式的复制，当然带来内存的消耗。而bsxfun是虚拟的复制，实际上通过for来实现，等效于for(i=1:3),for(j=1:3),c(i,j)=a(i)+b(j);end,end。但bsxfun不会有使用matlab的for所带来额外时间。实际验证下这三种方式

>> c = bsxfun(@plus,a,b)
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

>> c = repmat(a,1,3)+repmat(b,3,1)
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

>> for(i=1:3),for(j=1:3),c(i,j)=a(i)+b(j);end,end,c
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

从计算时间上来说前两种实现差不多，远高于for的实现。但如果数据很大，第二种实现可能会有内存上的问题。所以bsxfun最好。

这里@plus是加法的函数数柄，相应的有减法@minus, 乘法@times, 左右除等，具体可见 doc bsxfun.

下面看一个更为实际的情况。假设我们有数据A和B, 每行是一个样本，每列是一个特征。我们要计算高斯核，既. 当然可以用双重for实现（如果第一直觉是用三重for的话…）。

K1 = zeros(size(A,1),size(B,1));
for i = 1 : size(A,1)
    for j = 1 : size(B,1)
        K1(i,j) = exp(-sum((A(i,:)-B(j,:)).^2)/beta);
    end
end

使用2,000×1,000大小的A和B, 运行时间为88秒。

考虑下面向量化后的版本：

sA = (sum(A.^2, 2));
sB = (sum(B.^2, 2));
K2 = exp(bsxfun(@minus,bsxfun(@minus,2*A*B', sA), sB')/beta);

使用同样数据，运行时间仅0.85秒，加速超过100倍。

如要判断两者结果是不是一样，可以如下

assert(all(all(abs(K1-K2)<1e-12)))

在高斯核中，我们通常可以取是样本之间的平均距离，既, 此时我们可以使用下面的代码：

n = size(X,1);
S = sum(X);
beta = 2*norm(X,'fro')^2/(n-1) - (S*S')*2/n/(n-1);

各位有兴趣可以比比与for实现的性能差异，呵呵。