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

Posted 2021-03-19

我的任务是编写计算矩阵的最佳程序。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));

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

答案

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

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

另一答案

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

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

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')