有没有办法对这个matlab代码进行矢量化？

【问题标题】：Is there any way to vectorize this matlab code?有没有办法对这个matlab代码进行矢量化？
【发布时间】：2016-03-04 19:12:05
【问题描述】：

我有以下代码：

function [Ps,Pd,Pv] = Arii2010_Modified_1Pixel(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
MeanOrientationAngleStep = 1/100;
StandardDeviationStep = 1/100;
NumberOfAngles = floor(2*pi/MeanOrientationAngleStep);
NumberOfStandandardDeviations = 1/StandardDeviationStep;
for i = 0:NumberOfAngles
    for j = 0:NumberOfStandandardDeviations
        theta = -pi+i*MeanOrientationAngleStep;
        sigma = j*StandardDeviationStep;
        Fv = FindFv(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33);
    end
end
end

function Fv = FindFv(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
coef1 = 1/2-(((1-2*sigma^2)*(1-sigma^2)*cos(2*theta))/(2*(1+sigma^2)));
coef2 = ((1-4*sigma^2)*(1-3*sigma^2)*(1-2*sigma^2)*(1-sigma^2)*cos(4*theta))/(8*(1+sigma^2)*(1+2*sigma^2)*(1+3*sigma^2));
coef3 = ((1-2*sigma^2)*(1-sigma^2)*sin(2*theta))/(2*sqrt(2)*(1+sigma^2));
coef4 = ((1-4*sigma^2)*(1-3*sigma^2)*(1-2*sigma^2)*(1-sigma^2)*sin(4*theta))/(4*sqrt(2)*(1+sigma^2)*(1+2*sigma^2)*(1+3*sigma^2));
Fv1 = C11/(coef1+coef2);
aQuadratic1 = -2*coef2*(coef1+coef2);
aQuadratic2 = (coef3-coef4)^2;
aQuadratic = aQuadratic1 - aQuadratic2;
bQuadratic1 = 2*C11*coef2-C22*(coef1+coef2);
bQuadratic2 = -2*C12_real*(coef3-coef4);
bQuadratic = bQuadratic1 - bQuadratic2;
cQuadratic1 = C11*C22;
cQuadratic2 = C12_real^2+C12_imag^2;
cQuadratic = cQuadratic1 - cQuadratic2;
rQuadratic = roots([aQuadratic bQuadratic cQuadratic]);
sel1 = rQuadratic == real(rQuadratic);
rQuadratic = rQuadratic(sel1);
sel2 = rQuadratic == abs(rQuadratic);
rQuadratic = rQuadratic(sel2);
Fv2 = max(rQuadratic);
aCubic1 = 2*coef2*(1-coef1+coef2)*(coef1+coef2);
aCubic2 = 2*coef2*(coef3-coef4)*(coef3+coef4);
aCubic3 = -(coef1+coef2)*(coef3+coef4)^2;
aCubic4 = 2*coef2^3;
aCubic5 = -(1-coef1+coef2)*(coef3-coef4)^2;
aCubic = aCubic1+aCubic2-aCubic3-aCubic4-aCubic5;
bCubic1 = -2*coef2*(C11*(1-coef1+coef2)+C33*(coef1+coef2))+C22*(1-coef1+coef2)*(coef1+coef2);
bCubic2 = -2*coef2*(C23_real*(coef3-coef4)+C12_real*(coef3+coef4))+2*C13_real*(coef3-coef4)*(coef3+coef4);
bCubic3 = 2*C23_real*(coef1+coef2)*(coef3+coef4)+C11*(coef3+coef4)^2;
bCubic4 = (4*C13_real+C22)*coef2^2;
bCubic5 = 2*C12_real*(1-coef1+coef2)*(coef3-coef4)+C33*(coef3-coef4)^2;
bCubic = bCubic1+bCubic2-bCubic3-bCubic4-bCubic5;
cCubic1 = 2*C11*C33*coef2-C11*C22*(1-coef1+coef2)-C22*C33*(coef1+coef2);
cCubic2 = 2*(C12_real*C23_real-C12_imag*C23_imag)*coef2-2*(C13_real*C23_real+C13_imag*C23_imag)*...
    (coef3-coef4)-2*(C12_real*C13_real+C12_imag*C13_imag)*(coef3+coef4);
cCubic3 = -(C23_real^2+C23_imag^2)*(coef1+coef2)-2*C11*C23_real*(coef3+coef4);
cCubic4 = 2*(C13_real^2+C13_imag^2+C13_real*C22)*coef2;
cCubic5 = -(C12_real^2+C12_imag^2)*(1-coef1+coef2)-2*C33*C12_real*(coef3-coef4);
cCubic = cCubic1+cCubic2-cCubic3-cCubic4-cCubic5;
dCubic1 = C11*C22*C33;
dCubic2 = 2*(C12_real*(C13_real*C23_real+C13_imag*C23_imag)+C12_imag*(C13_imag*C23_real-C13_real*C23_imag));
dCubic3 = C11*(C23_real^2+C23_imag^2);
dCubic4 = C22*(C13_real^2+C13_imag^2);
dCubic5 = C33*(C12_real^2+C12_imag^2);
dCubic = dCubic1+dCubic2-dCubic3-dCubic4-dCubic5;
rCubic = roots([aCubic bCubic cCubic dCubic]);
sel3 = rCubic == real(rCubic);
rCubic = rCubic(sel3);
sel4 = rCubic == abs(rCubic);
rCubic = rCubic(sel4);
Fv3 = max(rCubic);
Fv = min([Fv1 Fv2 Fv3]);
end

我想知道是否有像bsxfun、blockproc 或arrayfun 这样的函数来矢量化这部分代码？ 

for i = 0:NumberOfAngles
    for j = 0:NumberOfStandandardDeviations
        theta = -pi+i*MeanOrientationAngleStep;
        sigma = j*StandardDeviationStep;
        Fv = FindFv(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33);
    end
end

【问题讨论】：

  • 由于FindFv 似乎包含roots，它一次接受一个多项式，我认为您不能进行太多矢量化：您可能需要调用FindFvNumberOfAngles*NumberOfStandandardDeviations 次。如果不是为了这个，您可以对FindFv 中的每个操作进行矢量化处理，并使用max/min 和感兴趣的维度上的可选输入参数。它可以帮助向量化roots 之前的所有内容，并使用双循环来覆盖每个多项式。它会更快，但我不知道多少。
  • @AndrasDeak 事实上我不想向量化 FindFv 函数中的代码。我想对调用 FindFv 函数的部分 for i = 0:NumberOfAngles for j = 0:NumberOfStandandardDeviations theta = -pi+i*MeanOrientationAngleStep; sigma = j*StandardDeviationStep; Fv = FindFv(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33); end end 进行矢量化
  • 不幸的是，矢量化不是魔术棒和精灵尘的秘诀：) 如果您调用一个函数 1000 次，而该函数一次只能工作一个，那么您无能为力.有一件事：在 MEX 中重写函数并预编译它。
  • @AndrasDeak 能不能多解释一下最后一句话，介绍好的资源入手？
  • MEX 是一种编写 C 代码的方法，它不太可能与 matlab 接口。 Some docs and examples are here。由于 C 和很容易遇到的微妙问题，使用 MEX 是一件很痛苦的事情。但它使用编译后的代码，因此比高级 MATLAB 更快。无论如何，我对 MEX 和正确的矢量化都不太熟悉。我会尝试首先对您的函数进行矢量化并将循环放在roots 位周围，这可能更简单/更快。我真的不知道你的胜算。

标签： performance matlab vectorization nested-loops mex


【解决方案1】：

我已经评论说，您可能应该将循环放在矢量化工作函数中，或/并使用已编译的 MEX 代码。这是我对第一种方法的看法：

function [Ps,Pd,Pv] = Arii2010_Modified_1Pixel(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
MeanOrientationAngleStep = 1/100;
StandardDeviationStep = 1/100;
NumberOfAngles = floor(2*pi/MeanOrientationAngleStep);
NumberOfStandandardDeviations = 1/StandardDeviationStep;

%// set up vectorized input
thetav = -pi + (0:NumberOfAngles)*MeanOrientationAngleStep;
sigmav = (0:NumberOfStandandardDeviations)*StandardDeviationStep;

%// grid for parameter pairs
[theta, sigma] = meshgrid(thetav,sigmav);

%// run vectorized function
Fv_vec = FindFv_vec(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33);
end


function Fv_vec = FindFv_vec(sigma,theta,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
coef1 = 1./2-(((1-2.*sigma.^2).*(1-sigma.^2).*cos(2.*theta))./(2.*(1+sigma.^2)));
coef2 = ((1-4.*sigma.^2).*(1-3.*sigma.^2).*(1-2.*sigma.^2).*(1-sigma.^2).*cos(4.*theta))./(8.*(1+sigma.^2).*(1+2.*sigma.^2).*(1+3.*sigma.^2));
coef3 = ((1-2.*sigma.^2).*(1-sigma.^2).*sin(2.*theta))./(2.*sqrt(2).*(1+sigma.^2));
coef4 = ((1-4.*sigma.^2).*(1-3.*sigma.^2).*(1-2.*sigma.^2).*(1-sigma.^2).*sin(4.*theta))./(4.*sqrt(2).*(1+sigma.^2).*(1+2.*sigma.^2).*(1+3.*sigma.^2));
%// this is OK:
Fv1 = C11./(coef1+coef2);

aQuadratic1 = -2.*coef2.*(coef1+coef2);
aQuadratic2 = (coef3-coef4).^2;
aQuadratic = aQuadratic1 - aQuadratic2;
bQuadratic1 = 2.*C11.*coef2-C22.*(coef1+coef2);
bQuadratic2 = -2.*C12_real.*(coef3-coef4);
bQuadratic = bQuadratic1 - bQuadratic2;
cQuadratic1 = C11.*C22;
cQuadratic2 = C12_real.^2+C12_imag.^2;
cQuadratic = cQuadratic1 - cQuadratic2;
%// move this down to one block for looping:
%//rQuadratic = roots([aQuadratic bQuadratic cQuadratic]);
%//sel1 = rQuadratic == real(rQuadratic);
%//rQuadratic = rQuadratic(sel1);
%//sel2 = rQuadratic == abs(rQuadratic);
%//rQuadratic = rQuadratic(sel2);
%//Fv2 = max(rQuadratic);

aCubic1 = 2.*coef2.*(1-coef1+coef2).*(coef1+coef2);
aCubic2 = 2.*coef2.*(coef3-coef4).*(coef3+coef4);
aCubic3 = -(coef1+coef2).*(coef3+coef4).^2;
aCubic4 = 2.*coef2.^3;
aCubic5 = -(1-coef1+coef2).*(coef3-coef4).^2;
aCubic = aCubic1+aCubic2-aCubic3-aCubic4-aCubic5;
bCubic1 = -2.*coef2.*(C11.*(1-coef1+coef2)+C33.*(coef1+coef2))+C22.*(1-coef1+coef2).*(coef1+coef2);
bCubic2 = -2.*coef2.*(C23_real.*(coef3-coef4)+C12_real.*(coef3+coef4))+2.*C13_real.*(coef3-coef4).*(coef3+coef4);
bCubic3 = 2.*C23_real.*(coef1+coef2).*(coef3+coef4)+C11.*(coef3+coef4).^2;
bCubic4 = (4.*C13_real+C22).*coef2.^2;
bCubic5 = 2.*C12_real.*(1-coef1+coef2).*(coef3-coef4)+C33.*(coef3-coef4).^2;
bCubic = bCubic1+bCubic2-bCubic3-bCubic4-bCubic5;
cCubic1 = 2.*C11.*C33.*coef2-C11.*C22.*(1-coef1+coef2)-C22.*C33.*(coef1+coef2);
cCubic2 = 2.*(C12_real.*C23_real-C12_imag.*C23_imag).*coef2-2.*(C13_real.*C23_real+C13_imag.*C23_imag).*...
    (coef3-coef4)-2.*(C12_real.*C13_real+C12_imag.*C13_imag).*(coef3+coef4);
cCubic3 = -(C23_real.^2+C23_imag.^2).*(coef1+coef2)-2.*C11.*C23_real.*(coef3+coef4);
cCubic4 = 2.*(C13_real.^2+C13_imag.^2+C13_real.*C22).*coef2;
cCubic5 = -(C12_real.^2+C12_imag.^2).*(1-coef1+coef2)-2.*C33.*C12_real.*(coef3-coef4);
cCubic = cCubic1+cCubic2-cCubic3-cCubic4-cCubic5;
dCubic1 = C11.*C22.*C33;
dCubic2 = 2.*(C12_real.*(C13_real.*C23_real+C13_imag.*C23_imag)+C12_imag.*(C13_imag.*C23_real-C13_real.*C23_imag));
dCubic3 = C11.*(C23_real.^2+C23_imag.^2);
dCubic4 = C22.*(C13_real.^2+C13_imag.^2);
dCubic5 = C33.*(C12_real.^2+C12_imag.^2);
dCubic = dCubic1+dCubic2-dCubic3-dCubic4-dCubic5;
%// move this down to one block for looping:
%//rCubic = roots([aCubic bCubic cCubic dCubic]);
%//sel3 = rCubic == real(rCubic);
%//rCubic = rCubic(sel3);
%//sel4 = rCubic == abs(rCubic);
%//rCubic = rCubic(sel4);
%//Fv3 = max(rCubic);

Fv_vec = zeros(size(theta));
for k=1:numel(theta)
   rQuadratic = roots([aQuadratic(k) bQuadratic(k) cQuadratic(k)]);
   %//sel1 = rQuadratic == real(rQuadratic);
   %//rQuadratic = rQuadratic(sel1);
   %//sel2 = rQuadratic == abs(rQuadratic);
   %//rQuadratic = rQuadratic(sel2);

   %// suggestion instead: (with optional tol=1e-10 or something else earlier)
   rQuadratic = rQuadratic(abs(rQuadratic-abs(rQuadratic))<1e-10);  % allow some tolerance due to floating-point
   Fv2 = max(rQuadratic);

   rCubic = roots([aCubic(k) bCubic(k) cCubic(k) dCubic(k)]);
   %//sel3 = rCubic == real(rCubic);
   %//rCubic = rCubic(sel3);
   %//sel4 = rCubic == abs(rCubic);
   %//rCubic = rCubic(sel4);
   rCubic = rCubic(abs(rCubic-abs(rCubic))<1e-10);
   Fv3 = max(rCubic);
   Fv_vec(k) = min([Fv1(k); Fv2; Fv3]);  %// Fv1 is a matrix!
end 

end

由于您没有提供任何示例输入，因此我没有尝试提供一个，所以我没有测试它。

我所做的是将每个算术运算更改为它们的矢量化对 (.*,./,.^)，并重组代码以使用单个循环来计算 roots。输入theta 和sigma 现在可以是矩阵，对应于输入向量对（用meshgrid 构造）。我的解决方案假定所有输入参数都是标量。我还更改了查找最大绝对值的部分，允许由于浮点运算导致的一些小错误。

【讨论】：

  • 我正在做这样的事情。尝试使用 for 循环来获取根，然后尝试使用 mex 文件更改那一小部分代码。这适用吗？
  • @sepideh 根查找可能会占用运行时的一小部分。您必须运行分析器才能查看占用最多运行时间的内容，但我怀疑所有这些算术运算。所以你应该像上面那样做，如果太慢（你找不到瓶颈），然后在 MEX 中重写整个函数（完成算术）。
  • 我已经运行了分析器并解释了这个瓶颈。请看Is there anyway to accelerate the process through MEX files? Will this help significantly?
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2019-03-10
  • 1970-01-01
  • 1970-01-01
  • 2016-03-26
  • 1970-01-01
  • 1970-01-01
  • 2021-02-11
相关资源
最近更新 更多
热门标签
Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode