通过向量化符号优化程序

Question

大家好，我正在从事图像处理工作，并在 MATLAB 中编写了一小段代码。代码很慢。

我在这里给出我的代码 sn-p

for i=1:10
    //find c1,c2,c3
    //c1 c2 and c3 change at each iteration

    u = (1./((abs(P-c1))^m) + 1./((abs(P-c2))^m) + 1./((abs(P-c3))^m));
    u1 = 1./((abs(P-c1))^m)./u;
    u2 = 1./((abs(P-c2))^m)./u;    
    u3 = 1./((abs(P-c3))^m)./u;
end

让我在这里解释一下变量：

P,u,u1,u2 and u3 are all matrices of size 512x512
c1,c2 and c3 are constants of dimension 1x1
m is a constant with value = 2

我想在循环中重复此操作（比如 10 次）。但是我的代码很慢。

分析器的结果如下：

程序的总运行时间为 4.6 秒。然而，上面列出的四个步骤本身需要大约 80% 的时间。

所以我想让我的代码运行更快。 我的第一次编辑

我修改后的代码sn-p

for i=1:10

    //find c1 and c2
    //c1 and c2 changes at each iteration

    a=((abs(P-c1))^m); 
    b=((abs(P-c2))^m); 
    c=((abs(P-c3))^m);

    x=1./a; y=1./b; z=1./c;
    u = (x + y + z);
    u1 = x./u;
    u2 = y./u;    
    u3 = z./u;
end

现在程序计算上述步骤的计算时间为 2.47 秒：

所以这比我的第一种方法快得多。

第二次修改 对于 i=1:10 //找到c1,c2,c3 //c1 c2 和 c3 在每次迭代时改变

    a=(P-c1).*(P-c1); 
    b=(P-c2).*(P-c2); 
    c=(P-c3).*(P-c3);

    x=1./a; y=1./b; z=1./c;
    u = (x + y + z);
    u1 = x./u;
    u2 = y./u;    
    u3 = z./u;
end

现在程序计算时间为 0.808 秒。

上面描述的四个步骤计算上面的速度非常快。

我相信它可以做得更快。你们能帮我进一步优化 我的代码。 这对于大于 512 的矩阵（例如 1024 、 2048 等）非常有用。

提前致谢。

Answer 1

您当前的代码是：

a=((abs(P-c1))^m); 
b=((abs(P-c2))^m); 
c=((abs(P-c3))^m);

x=1./a; y=1./b; z=1./c;
u = (x + y + z);
u1 = x./u;
u2 = y./u;    
u3 = z./u;

首先，认识到绝对值函数是乘法的。所以|AB| = |A|x|B|。现在，abs(P-C1)^m 等价于 abs( (P-C1)^m )。

初步看一眼就表明瓶颈中的一些计算可以重用。具体来说，由于c1、c2 和c3 是常量，因此如果您尝试重用它们（以增加内存为代价），计算速度会稍微加快。

temp_P2 = P*P;
temp_PCA = P*ones(size(P));
temp_PCB = ones(size(P))*P;

a = abs(temp_P2 - c1*temp_PCA - c1*temp_PCB + c1^2 * length(P))

temp_PCA 和 temp_PCB 的计算也可以避免，因为乘以常数矩阵总是等于构造具有常数行或列的秩 1 矩阵。

我并不声称这些修改中的任何一个都会加快您的代码速度，但它们绝对值得一试。

Answer 2

第一个建议是：

如果 m = 2 并且它没有改变，为什么不尝试这些替代方案：

A*A

如果 m = 2 那么你真的需要 abs 吗？

你正在做的这部分

1./a

比

快

a.^(-1)

所以我在这部分没有看到更好的选择。

您可以尝试的另一件事是这个。而不是：

x=1./a; y=1./b; z=1./c;
    u = (x + y + z);
    u1 = x./u;
    u2 = y./u;    
    u3 = z./u;

你可以有这个：

    u = (x + y + z);
    u1 = 1./(a.*u);
    u2 = 1./(b.*u);    
    u3 = 1./(c.*u);

这样我想通过删除 3 个变量会快一点。但是代码变得不那么可读了。