计算(sin(x)-x)*x^{-3}的优化算法(matlab中)

问题描述 投票:0回答:1

我的任务是编写计算矩阵的最佳程序。Y,给定矩阵 X,其中

y = (sin(x)-x) x-3

这是我目前写的代码。

n = size(X, 1);
m = size(X, 2);
Y = zeros(n, m);
d = n*m; 

for i = 1:d
    x = X(i);
    if abs(x)<0.1
        Y(i) = -1/6+x.^2/120-x.^4/5040+x.^6/362880;
    else
        Y(i) = (sin(x)-x).*(x.^(-3));
    end
end

所以,一般来说,公式在0附近是不准确的, 所以我用泰勒定理来逼近它。

遗憾的是这个程序的准确率是91%,效率只有24%(所以比最优解慢了4倍)。

测试的样本约为1300万个,其中约600万个样本的值小于0.1。样本的范围是(-8π , 8π)。

目标精度(100%)为 4*epsilon 哪儿 epsilon 等于 2^(-52) (也就是说,程序计算的数字不应该比 "完美 "计算的数字大或小,而非 4*epsilon).

100*epsilon 意味着准确率达到86%。

你有什么办法可以让它更快更准确吗? 我既要寻找如何进一步变换给定公式的数学技巧,又要寻找可以加速程序的一般MATLAB技巧?

EDIT:使用Horner方法,我已经成功地将这个程序的效率提高到81%(准确率仍然是91%)。

function Y = main(X)

Y = (sin(X)-X).*(X.^(-3));
i = abs(X) < 0.1;
Y(i) = horner(X(i));

function y = horner (x)

pow = x.*x;
y = -1/6+pow.*(1/120+pow.*(-1/5040+pow./362880));

你有什么进一步的想法来改进它吗?

algorithm matlab optimization numerical-methods
1个回答
1
投票

你可以用向量化的代码代替你的循环。这通常比循环更有效率,因为循环里有一个条件,这对于 分支预测:

Y = (sin(X)-X).*(X.^(-3));
i = abs(X) < 0.1;
Y(i) = -1/6+X(i).^2/120-X(i).^4/5040+X(i).^6/362880;

重写一次方程以避免立方根,使该计算速度提高了3倍。

Y = (sin(X)./X - 1) ./ (X.*X);

速度比较。

下面的脚本比较了这种方法和OP的循环代码的时间。我使用的数据有700万个值均匀分布在(-8π,8π),另外600万个值均匀分布在(-0.1,0.1)。

OP的循环代码需要2.4412秒,而矢量化解法需要0.7224秒。使用OP的Horner方法和重写的 sin 表达式需要0.1437秒。

X = [linspace(-8*pi,8*pi,7e6), linspace(-0.1,0.1,6e6)];
timeit(@()method1(X))
timeit(@()method2(X))

function Y = method1(X)
n = size(X, 1);
m = size(X, 2);
Y = zeros(n, m);
d = n*m; 

for i = 1:d
    x = X(i);
    if abs(x)<0.1
        Y(i) = -1/6+x.^2/120-x.^4/5040+x.^6/362880;
    else
        Y(i) = (sin(x)-x).*(x.^(-3));
    end
end
end

function Y = method2(X)
Y = (sin(X)-X).*(X.^(-3));
i = abs(X) < 0.1;
Y(i) = -1/6+X(i).^2/120-X(i).^4/5040+X(i).^6/362880;
end

function Y = method3(X)
Y = (sin(X)./X - 1) ./ (X.*X);
i = abs(X) < 0.1;
Y(i) = horner(X(i));
end

function y = horner (x)
pow = x.*x;
y = -1/6+pow.*(1/120+pow.*(-1/5040+pow./362880));
end

3
投票

程序似乎在很大的输入范围内都能正常工作。

x = linspace(-8*pi,8*pi,13e6); % 13 million samples in the desired range
y = (sin(x)-x)./x.^3;
plot(x,y)

enter image description here

Due due due 四舍五入对于很小的x值,你可能会有计算上的问题。

x = 0
y = (sin(x)-x)./x.^3
y =

   NaN

你已经有了0附近的函数的泰勒数列展开图 由于泰勒展开图不包括除法 x在这个区域附近,你可以期待泰勒函数有更好的表现。

x = -1e-6:1e-9:1e-6;
y = (sin(x)-x)./x.^3;
y_taylor = -1/6 + x.^2/120 - x.^4/5040 + x.^6/362880;
plot(x,y,x,y_taylor); legend('y','taylor expansion','location','best')

enter image description here

© www.soinside.com 2019 - 2024. All rights reserved.