性能 – matlab代码的矢量化

我对矢量化有点新意.已经尝试过但不能.有人可以帮助我对这段代码进行矢量化,并简要说明你是如何做到这一点的,这样我就可以适应思维过程.谢谢.

function [result] = newHitTest (point,Polygon,r,tol,stepSize)
%This function calculates whether a point is allowed.

%First is a quick test is done by calculating the distance from point to 
%each point of the polygon. If that distance is smaller than range "r", 
%the point is not allowed. This will slow down the algorithm at some 
%points, but will greatly speed it up in others because less calls to the 
%circleTest routine are needed.
polySize=size(Polygon,1);
testCounter=0;

for i=1:polySize
d = sqrt(sum((Polygon(i,:)-point).^2));

if d < tol*r
    testCounter=1;
    break
end
end

if testCounter == 0
circleTestResult = circleTest (point,Polygon,r,tol,stepSize);
testCounter = circleTestResult;
end

result = testCounter;

鉴于Polygon是2维的信息,point是行向量而其他变量是标量,这是新函数的第一个版本(向下滚动以查看有很多方法可以为这只猫设置皮肤)：

function [result] = newHitTest (point,Polygon,r,tol,stepSize)
result = 0;
linDiff = Polygon-repmat(point,size(Polygon,1),1);
testLogicals = sqrt( sum( ( linDiff ).^2 ,2 )) < tol*r;    
if any(testLogicals); result = circleTest (point,Polygon,r,tol,stepSize); end

Matlab中矢量化的思考过程涉及尝试使用单个命令尽可能多地操作数据.大多数基本的内置Matlab函数在多维数据上运行非常有效.使用for循环与此相反,因为您将数据分解为较小的段以进行处理,每个段都必须单独解释.通过使用for循环进行数据分解,您可能会失去与Matlab内置函数背后的高度优化代码相关的一些巨大性能优势.

在您的示例中要考虑的第一件事是主循环中的条件中断.你无法摆脱矢量化过程.相反,计算所有可能性,为数据的每一行创建结果数组,然后使用any关键字查看是否有任何行已发出信号表示应调用circleTest函数.

注意：在Matlab中有效地有条件地打破计算是不容易的.但是,由于您只是在循环中计算欧几里德距离的形式,您可能会通过使用矢量化版本并计算所有可能性来看到性能提升.如果你的循环中的计算更昂贵,输入数据很大,并且你想在遇到某个条件时立即爆发,那么用编译语言编写的matlab扩展可能比矢量化版本快得多.你可能正在进行不必要的计算.但是,假设您知道如何使用编译为本机代码的语言编写与Matlab内置函数的性能相匹配的代码.

回到主题……

首先要做的是在Polygon和行向量点之间取线性差异(代码示例中为linDiff).要以矢量化方式执行此操作,2个变量的维度必须相同.实现此目的的一种方法是使用repmat复制每一行的点,使其与Polygon相同.但是,bsxfun通常是repmat(as described in this recent SO question)的优秀替代品,使代码成为……

function [result] = newHitTest (point,Polygon,r,tol,stepSize)
result = 0;
linDiff = bsxfun(@minus, Polygon, point);
testLogicals = sqrt( sum( ( linDiff ).^2 ,2 )) < tol*r;    
if any(testLogicals); result = circleTest (point,Polygon,r,tol,stepSize); end

我通过在第二轴上求和将您的d值转换为d列(注意从Polygon中删除数组索引并在sum命令中添加2).然后我进一步评估逻辑数组testLogicals内联计算距离度量.您很快就会发现重载矢量化的缺点是它可以使代码对于不熟悉Matlab的人来说更不易读,但性能提升是值得的.评论是非常必要的.

现在,如果你想完全疯狂,你可以说测试函数现在非常简单,它保证使用’匿名函数’或’lambda’而不是完整的函数定义.是否值得执行circleTest的测试也不需要stepSize参数,这也许是使用匿名函数的另一个原因.您可以将测试转换为匿名函数,然后在调用脚本中使用circleTest,使代码在某种程度上自我记录. . .

doCircleTest = @(point,Polygon,r,tol) any(sqrt( sum( bsxfun(@minus, Polygon, point).^2, 2 )) < tol*r);

if doCircleTest(point,Polygon,r,tol)
    result = circleTest (point,Polygon,r,tol,stepSize); 
else
    result = 0;
end

现在一切都是矢量化的,函数句柄的使用给了我另一个想法. . .

如果您计划在代码中的多个点执行此操作,那么重复if语句会有点难看.为了保持dry,将条件函数的测试放入单个函数似乎是明智的,就像在原始帖子中所做的那样.但是,该函数的效用将非常狭窄 – 它只测试是否应该执行circleTest函数,然后在需要时执行它.

现在想象一下,经过一段时间,你有一些其他的条件函数,就像circleTest一样,有自己的doCircleTest.很可能重用条件切换代码.为此,创建一个类似于原始函数的函数,它采用默认值,计算上廉价的测试函数的布尔结果,以及昂贵的条件函数及其相关参数的函数句柄……

function result = conditionalFun( default, cheapFunResult, expensiveFun, varargin )
if cheapFunResult
    result = expensiveFun(varargin{:});
else
    result = default;
end
end %//of function

您可以使用以下命令从主脚本中调用此函数. . .

result = conditionalFun(0, doCircleTest(point,Polygon,r,tol), @circleTest, point,Polygon,r,tol,stepSize);

…它的美妙之处在于你可以使用任何测试,默认值和昂贵的功能.对于这个简单的例子来说,或许有点矫枉过正,但是当我提出使用函数句柄的想法时,我的思绪就会徘徊.