让MATLAB更快

让MATLAB更快

MATLAB是一门在数值运算方面非常高效的语言，对于提高MATLAB的效率，有一些大家熟知的方法，比如避免使用for循环，向量化，使用C-mex写核心，等等。这些原则大体是没错的，但是教条地运用有时反而会降低效率。

根据问题的规模，适当选择向量化的策略

向量化的基本思想是以空间换时间，通过把要处理的对象转化成一个矩阵，集中处理。这比起非向量化的处理方式往往需要消耗更多的内存。在内存紧张的时候，内存分配的开销可能是很大的。

比如，给定一个列向量v1，和一个行向量v2, 计算矩阵M，使得M(i, j) = v1(i) + v2(j)。一种常用的向量化实现：

m = length(v1); 

n = length(v2);

M = repmat(v1, [1, n]) + repmat(v2, [m, 1]);

这个过程中，除了M以外，还要创建两个大小为m x n的临时矩阵。repmat创建临时矩阵有一定的overhead，主要花费在元素索引上。由于元素加法的向量化收益本身并不显著，扣除repmat的额 外开销后，得益就更小了。另外，当这些矩阵很大，接近物理内存容纳能力的时候，必然导致和虚拟内存进行频繁的交换，耗费大量的IO时间，向量化可能得不偿 失。这种情况下，下面这种向量化程度略低的方式可能是更好的选择：

M = zeros(m, n);

if m < n

    for i = 1 : m

        M(i, :) = v1(i) + v2;

    end

else

    for i = 1 : n

        M(:, i) = v1 + v2(i);

    end

end

这个例子是要说明在设计向量化的时候，要对其它潜在的影响因素也进行考虑，权衡得失。对于这个例子本身，在2007版以后的matlab，有一个通用的简洁的写法。

M = bsxfun(@plus, v1, v2);

MATLAB内部对于bsxfun有很专业的优化，这个写法在时间和空间上都比较节省。

对一组d x d矩阵求逆，假设这些矩阵放在大小为d x d x n的数组里。在一般条件下，多个矩阵求逆没有特别的向量化方法，通常就是逐个求

n = size(A, 3);

B = zeros(size(A));

for i = 1 : n

    B(:, :, i) = inv(A(:, :, i));

end

不过，当d很小而n很大的时候，比如对大量2x2矩阵求逆（这在设计几何的问题中很常见），那么直接利用2x2矩阵的求逆共识，进行向量化了：

a = reshape(A(1, 1, :), [1, n]);

b = reshape(A(1, 2, :), [1, n]);

c = reshape(A(2, 1, :), [1, n]);

d = reshape(A(2, 2, :), [1, n]);



g = 1 ./ (a .* d - b .* c);

r11 = d .* g;

r12 = -b .* g;

r21 = -c .* g;

r22 = a .* g;



B = reshape([r11; r21; r12; r22], [2, 2, n]); 

在正式的代码中，还需要处理可能为奇异(a*d-b*c=0)的情况。

某些非数值计算问题的向量化

有些问题，表面上不像是数值计算，但是通过适当的转化，也可以通过向量化解决。

非等量复制: 比如有一组数[10, 20, 30], 各自要复制成[3, 2, 5]份，最终产生出 [10 10 10 20 20 30 30 30 30 30]。这个问题，虽然简单，但是因为不是那么“整齐”，因此向量化不是特别直接。下面给出一种思路：

首先产生出起始位置标志向量[1 0 0 1 0 1 0 0 0 0]，然后通过cumsum，产生出[1 1 1 2 2 3 3 3 3 3]，最后以此为索引，可以得到结果。

% vs: the values to be copied

% ns: the numbers of copies 



vs = vs(ns > 0);

ns = ns(ns > 0);



sp = [1, cumsum(ns(1:end-1)) + 1];

result = vs(cumsum(sp));

对于图像处理或者二维信号问题，善用filter。

比如，要对每个像素，基于其局部邻域进行一些计算。即使整个计算本身不是线性的，但是只要能分解成线性局部运算的组合，就可以利用filter。

比如，对于图像I, 产生V，使得V(i, j)是像素I(i, j)周围w x w邻域的像素的方差。

I = im2double(I);

h = ones(w, w) / (w * w);

M = imfilter(I, h, 'symmetric');

M2 = imfilter(I.^2, h, 'symmetric');

V = M2 - M.^2;