到目前为止,我已搜查,我知道有几种方法(1,2,3,4)到目前为止,我用下面的代码:
Fv_calc(:,2) = arrayfun(@(n) MaxPositiveRoot2DegreePolynomial(QuadraticCoefficients(n,:)), 1:size(QuadraticCoefficients,1));
其中MaxPositiveRoot2DegreePolynomial是以下函数:
function root = MaxPositiveRoot2DegreePolynomial(c)
d = c./c(1);
root = eig([0 -d(3);1 -d(2)]);
condition = root == abs(root);
root = root(condition);
if isempty(root)
root = 0;
end
root = max(root);
end
其中QuadraticCoefficients是62271x3矩阵包含的系数的每一行a
,b
,c
一般二次方程的。ax^2+bx+c
关于我要解决的问题,所有根都是真实的,因此我使用了一种修改过的函数来查找根,以免浪费时间检查根是否为真。
通过对代码进行性能分析,我发现该函数的每次运行MaxPositiveRoot2DegreePolynomial
大约需要0.000047秒,对于62271次运行来说大约是2.914925371秒。但是代码
Fv_calc(:,2) = arrayfun(@(n) MaxPositiveRoot2DegreePolynomial(QuadraticCoefficients(n,:)), 1:size(QuadraticCoefficients,1));
需要3.198597803才能运行。所以大约0.283672431只是取arrayfun
。我想看看这次是否有减少的方法?
这是一个完美的示例,它显示了您应该继续使用循环的地方。我打算使用的向量化解决方案最终比循环方法要慢。
但这可能取决于您的Matlab版本,因为我已经将Matlab R2015b与优化的JIT-Compiler一起使用,在较旧的版本中,矢量化方法可能会更快。
[n,~] = size(C);
%// division
d = bsxfun(@rdivide,C,C(:,1));
%// create blocks for eigen matrix block procession
X = ( [zeros(n,1) ones(n,1) -d(:,3) -d(:,2)] ).'; %'
Y = reshape(X,2,[]);
%// use block processing to get eigenvalues of blocks
rt = blockproc(Y,[2 2],@(x) eig(x.data) ).'; %'
%// check if eigen values are real positives
condition = rt == abs(rt);
%// output filter
out = rt(condition);
function [t] = bench()
%// load your data
DATA = load('matlab');
C = DATA.QuadraticCoefficients;
% functions to compare
fcns = {
@() thewaywewalk(C);
@() sepideh(C);
};
% timeit
t = zeros(2,1);
for ii = 1:10;
t = t + cellfun(@timeit, fcns);
end
end
function out = thewaywewalk(C) %thewaywewalk
[n,~] = size(C);
d = bsxfun(@rdivide,C,C(:,1));
X = ( [zeros(n,1) ones(n,1) -d(:,3) -d(:,2)] ).'; %'
Y = reshape(X,2,[]);
rt = blockproc(Y,[2 2],@(x) eig(x.data) ).'; %'
condition = rt == abs(rt);
out = rt(condition);
end
function out = sepideh(C) %sepideh
[n,~] = size(C);
out = zeros(n,1);
for ii = 1:n
c = C(ii,:);
d = c./c(1);
rt = eig([0 -d(3);1 -d(2)]);
condition = rt == abs(rt);
rt = rt(condition);
if isempty(rt)
rt = 0;
end
rt = max(rt);
out(ii) = rt;
end
end
ans =
12.2086 %// thewaywewalk
9.2176 %// sepideh
我仍然看不到if条件的意义。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句