Эффективный способ получить полномочия вектора

Question

Я написал код, который численно использует полиномы Лежандра до некоторого высокого n-го порядка. Например:

....
case 8 
p = (6435*x.^8-12012*x.^6+6930*x.^4-1260*x.^2+35)/128; return
case 9 
...

Если векторx долго это может стать медленным. Я видел, что есть разница в производительности между скажем x.^4 а также x.*x.*x.*x и думал, что смогу использовать это для улучшения своего кода. Я использовал timeit и обнаружил, что для:

x=linspace(0,10,1e6);
f1= @() power(x,4)
f2= @() x.4;
f3= @() x.^2.^2
f4= @() x.*x.*x.*x

f4 быстрее в 2 раза, чем остальные. Однако, когда я иду в x.^6 разница очень мала (x.*x.*x).^2 а также x.*x.*x.*x.*x.*x (в то время как все другие варианты медленнее).

Можно ли сказать, что будет наиболее эффективным способом получить степень вектора? Можете ли вы объяснить, почему такая большая разница в производительности?

Answer 1

Это не совсем ответ на ваш вопрос, но он может решить вашу проблему:

x2 = x.*x; % or x.^2 or power(x,2), whichever is most efficient
p = ((((6435*x2-12012)*x2+6930)*x2-1260)*x2+35)/128

Таким образом, вы делаете мощность только один раз, и только с показателем 2. Этот прием может быть применен ко всем многочленам Лежандра (в многочленах нечетной степени один x2 заменяется x).

Answer 2

Кажется, что Mathworks имеет специальные квадраты в своей степенной функции (к сожалению, это все встроенный закрытый источник, который мы не можем видеть). В моем тестировании на R2013b кажется, что .^, power, а также realpow использовать тот же алгоритм. Для квадратов, я полагаю, они имеют специальный случай, чтобы это было x.*x,

1.0x (4.4ms):   @()x.^2
1.0x (4.4ms):   @()power(x,2)
1.0x (4.5ms):   @()x.*x
1.0x (4.5ms):   @()realpow(x,2)
6.1x (27.1ms):  @()exp(2*log(x))

Для кубов история другая. Они больше не в специальном корпусе. Снова, .^, power, а также realpow все похожи, но гораздо медленнее на этот раз:

1.0x (4.5ms):   @()x.*x.*x
1.0x (4.6ms):   @()x.*x.^2
5.9x (26.9ms):  @()exp(3*log(x))
13.8x (62.3ms): @()power(x,3)
14.0x (63.2ms): @()x.^3
14.1x (63.7ms): @()realpow(x,3)

Давайте перейдем к 16-й степени, чтобы увидеть, как масштабируются эти алгоритмы:

1.0x (8.1ms):   @()x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x
2.2x (17.4ms):  @()x.^2.^2.^2.^2
3.5x (27.9ms):  @()exp(16*log(x))
7.9x (63.8ms):  @()power(x,16)
7.9x (63.9ms):  @()realpow(x,16)
8.3x (66.9ms):  @()x.^16

Так: .^, power а также realpow все выполняются в постоянном времени относительно показателя степени, если только он не был в специальном регистре (-1 также, кажется, был в специальном регистре). С использованием exp(n*log(x)) трюк также постоянное время по отношению к показателю степени и быстрее. Единственный результат, который я не совсем понимаю, почему повторное возведение в квадрат медленнее, чем умножение.

Как и следовало ожидать, увеличение размера x в 100 раз увеличивает время одинаково для всех алгоритмов.

Итак, мораль этой истории? При использовании скалярных целочисленных показателей всегда делайте умножение самостоятельно. Там много умов в power и друзья (экспонента может быть плавающей точкой, вектором и т. д.). Единственным исключением является ситуация, когда Mathworks выполнил оптимизацию для вас. В 2013b, похоже, x^2 а также x^(-1), Надеюсь, они будут добавлять больше со временем. Но, как правило, возведение в степень сложно, а умножение легко. В коде, чувствительном к производительности, я не думаю, что вы можете ошибиться, всегда набирая x.*x.*x.*x, (Конечно, в вашем случае, следуйте совету Луиса и используйте промежуточные результаты в течение каждого семестра!)

function powerTest(x)

f{1} = @() x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x.*x;
f{2} = @() x.^2.^2.^2.^2;
f{3} = @() exp(16.*log(x));
f{4} = @() x.^16;
f{5} = @() power(x,16);
f{6} = @() realpow(x,16);

for i = 1:length(f)
    t(i) = timeit(f{i});
end

[t,idxs] = sort(t);
fcns = f(idxs);

for i = 1:length(fcns)
    fprintf('%.1fx (%.1fms):\t%s\n',t(i)/t(1),t(i)*1e3,func2str(fcns{i}));
end

Answer 3

Вот несколько мыслей:

power(x,4) а также x.^4 эквивалентны (просто прочитайте документацию).

x.*x.*x.*x вероятно, оптимизирован для чего-то вроде x.^2.^2

---

x.^2.^2 вероятно, оценивается как: Возьмите квадрат каждого элемента (быстро) и снова возьмите квадрат этого (снова быстро).

x.^4 вероятно, непосредственно оценивается как: Возьмите четвертую степень каждого элемента (медленно).

Не так странно видеть, что 2 быстрые операции занимают меньше времени, чем 1 медленная. Очень жаль, что оптимизация не выполняется в случае с Power 4, но, возможно, она не всегда будет работать или стоить дорого (проверка ввода, память?).

---

О времени: На самом деле разница намного больше, чем фактор 2!

Поскольку вы вызываете их в функции сейчас, накладные расходы на функцию добавляются в каждом случае, уменьшая относительные различия:

y=x;tic,power(x,4);toc
y=x;tic,x.^4;toc
y=x;tic,x.^2.^2;toc
y=x;tic,x.*x.*x.*x;toc

дам:

Elapsed time is 0.034826 seconds.
Elapsed time is 0.029186 seconds.
Elapsed time is 0.003891 seconds.
Elapsed time is 0.003840 seconds.

Таким образом, это разница почти в 10 раз. Однако обратите внимание, что разница во времени в секундах все еще незначительна, поэтому для большинства практических приложений я бы просто использовал простой синтаксис.