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

Question

我的任务是编写计算矩阵 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)（这意味着程序计算的数字不应大于或小于计算的数字 "perfectly" 比 4*epsilon).

100*epsilon表示准确率为86%。

你对如何让它更快更准确有什么想法吗？我正在寻找数学技巧进一步转换给定的公式，以及可以加速程序的一般 MATLAB 技巧？

编辑：使用 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));

您对如何改进它有任何进一步的想法吗？

Answer 1

程序似乎适用于大范围的输入：

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

截止日期 round-off errors，对于非常小的 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')

Answer 2

您可以用矢量化代码替换循环。这通常比循环更有效，因为循环中有一个条件，这对 branch prediction:

不利

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 s.

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

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

Optimizing algorithm calculating (sin(x)-x)*x^{-3} (in matlab)

algorithm

optimization

matlab

numerical-methods