Вычисление нормы векторного выражения || aW + bX + cY ||

Question

Я аспирант. Во введении моей диссертации меня заинтересовал компромисс между выразительностью и характеристиками инструментов линейной алгебры.

В качестве простого примера я использую вычисление нормы векторного выражения. Код C для моего примера:

float normExpression3(float a, float *W, float b, float *X, float c, float*Y){
double norm = 0;
for (int i=0; i<n; ++i) // n in [3e6; 2e8]
{
    float tmp = a*W[i]+b*X[i]+c*Y[i];
    norm+=tmp*tmp;
}
return sqrtf(norm);

}

Я сравниваю достигнутые характеристики с разными техниками. Поскольку векторы велики (несколько миллионов элементов), производительность ограничена пропускной способностью памяти. Однако между этими подходами существуют огромные различия.

Оптимизированная версия C, которую я написал, не выразительна (новая функция должна быть записана как 4-й вектор) и очень уродлива (с резьбой и векторизацией), но достигла 6,4 GFlops. С другой стороны, код MATLAB очень хорош:

result = norm(a*W+b*X+c*Y)

, но достигает только 0,28 GFlops.

C ++ шаблоны выражений à la Blitz ++ обеспечивает пользователю выразительность и производительность (6,5 GFlops).

В рамках моего анализа я хотел бы знать, как функциональные языки могут сравниваться с этими подходами. Я думал о том, чтобы показать пример в Haskell или в OCaml (AFAIK, оба, как считается, хорошо подходят для такого рода операций).

Я не знаю ни одного из этих языков. Я мог бы научиться одному из них, чтобы привести мой пример, но это не будет честным сравнением: я не уверен, что смогу обеспечить реализацию, позволяющую как выразительность, так и производительность.

Итак, два моих вопроса: 1) какой язык лучше всего подходит? 2) как можно эффективно вычислить норму векторных выражений без ущерба для общности реализации?

Заранее спасибо!

Уилфрид К.

---

Редактировать: исправлен тип аккумулятора norm для float до double

Answer 1

Для чего стоит, следующая версия вашей функции в OCaml:

let normExpression3 a w b x c y =
    let n = Array.length w in
    if not (n = Array.length x && n = Array.length y)
        then invalid_arg "normExpression3";
    let (@) = (Array.unsafe_get : float array -> int -> float) in
    let rec accum a w b x c y n i norm =
        if i = n then sqrt norm else
        let t = a *. (w @ i) +. b *. (x @ i) +. c *. (y @ i) in
        accum a w b x c y n (i + 1) (norm +. t)
    in accum a w b x c y n 0 0.

Это делает некоторые допуски для производительности, а именно:

• Не проверенный доступ к массиву (илискорее, проверки границ массива вручную выводятся из цикла)
• Доступ к мономорфному массиву
• Рекурсивный внутренний цикл, чтобы избежать упаковки и распаковки поплавкового аккумулятора
• Лямбда-лифтинг внутреннегоцикл, чтобы избежать обращения к закрытым значениям

Последняя оптимизация должна быть проверена по закрытому внутреннему циклу, так как при большом количестве параметров разлив регистров может доминировать над стоимостью ссылки на закрытые параметры.

Обратите внимание, что обычно такого рода оптимизации не будут беспокоить, если только вы не будете соревноваться в бенчмарке ;-) Обратите внимание, что вам понадобится , чтобы проверить это с 64-бит OCaml, так как массивы иначе ограничены 4 мегаэлементами.

Answer 2

1) какой язык лучше всего подходит?

Любой из них используется для такого рода задач.Основной проблемой будет наличие необходимых библиотек (например, для векторов или матриц) и необходимость параллелизма.

Библиотеки, такие как vector и repa в Haskell хорошо подходят.И в случае с repa вы также получаете параллелизм.

2) как можно эффективно вычислить норму векторных выражений без ущерба для общности реализации?

Один из подходов заключается в использовании методов метапрограммирования для генерации специализированных реализаций вычислительных ядер из высокоуровневых описаний.В функциональных языках это довольно распространенный метод, основанный на небольших доменных языках, с пользовательскими генераторами кода.

См., Например, Специализированные генераторы симуляторов для высокопроизводительных методов Монте-Карло илиработа в OCaml на FFTW.

Answer 3

Не так много, как ответ функционального языка как таковой, но, пожалуйста, обратите внимание, что ваша реализация для вычисления norm (в C) сильно отличается от того, как matlab фактически вычисляет его.

И да, для этого действительно есть очень веские причины.Скорее всего, ваше приближение norm в значительной степени бесполезно (как то, как оно реализовано в настоящее время) для любого реального случая использования.Пожалуйста, не стоит недооценивать «трудности», связанные с вычислением численно превосходных приближений norm с.

Answer 4

Как сказал дон, рассмотрите Repa.Вот простой код, с которого можно начать.

import Data.Array.Repa as R

len :: Int
len = 50000

main = do
    let ws = R.fromList (Z :. len) [0..len-1]
        xs = R.fromList (Z :. len) [10498..10498 + len - 1]
        ys = R.fromList (Z :. len) [8422..8422 + len - 1]
    print (multSum 52 73 81 ws xs ys)

multSum a b c ws xs ys = R.map (a*) ws +^ R.map (b*) xs +^ R.map (c*) ys

Это все еще оставляет вам возможность найти хороший способ получить данные с диска и в массив Repa.Я думаю, что все это следует читать как Llazy ByteString и использовать Repa.fromFunction, возможно, кто-то вмешается более умным способом.