Оптимизация алгоритма расчета (sin (x) -x) * x ^ {- 3} (в matlab) - PullRequest
Новая
       28

Оптимизация алгоритма расчета (sin (x) -x) * x ^ {- 3} (в matlab)

2 голосов
/ 22 апреля 2020

Моя задача - написать оптимальную программу, которая вычисляет матрицу 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 раза медленнее, чем оптимальное решение).

В тестах участвуют около 13 миллионов образцов, из которых около 6 миллионов имеют значение меньше 0,1. Диапазон выборок (-8π, 8π).

Целевая точность (100%) равна 4*epsilon, где epsilon равно 2^(-52) (это означает, что числа, рассчитанные программой, не должны быть больше или меньше, чем числа, рассчитанные «отлично», чем 4*epsilon).

100*epsilon означает точность 86%.

У вас есть какие-либо идеи о том, как сделать это быстрее и точнее? Я ищу как математические приемы о том, как далее преобразовать данную формулу, так и общие советы MATLAB, которые могут ускорить программы?

РЕДАКТИРОВАТЬ: Используя метод Хорнера, я удалось повысить эффективность до 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));

Есть ли у вас какие-либо дальнейшие идеи о том, как улучшить ее?

Ответы [ 2 ]

3 голосов
/ 22 апреля 2020

Программа работает нормально для большого диапазона ввода: 

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

Из-за ошибок округления , у вас могут возникнуть проблемы с вычислением его для очень малых значений 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

1 голос
/ 22 апреля 2020

Вы можете заменить свой l oop векторизованным кодом. Это обычно более эффективно, чем l oop, потому что в l oop есть условное выражение, что плохо для предсказания ветвления : 

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;

Переписать основное уравнение, чтобы избежать cubi c root дает 3-кратное ускорение для этого вычисления: 

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

Сравнение скорости:

Следующий скрипт сравнивает время для этого метод по сравнению с кодом l oop ОП. Я использую данные, которые имеют 7 миллионов значений, равномерно распределенных в (-8π, 8π), и еще 6 миллионов значений, равномерно распределенных в (-0,1,0.1).

Код OP * l oop занимает 2,4412 с, и векторизованное решение занимает 0,7224 с. Используя метод Хорнера в OP и переписанное выражение 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
...