Эффективная обработка списка малых векторов с помощью компиляции

Ответ 1

Ваше приложение фактически достаточно, чтобы понять, в чем проблема. Для первой версии вы вызываете CopyTensor во внутренний цикл, и это основная причина неэффективности, так как в куче должно быть выделено множество небольших буферов, а затем выпущено. Чтобы проиллюстрировать, вот версия, которая не копирует:

fun3 =
 Compile[{{vec, _Real, 2}},
   Module[{sum = vec[[1]], len = Length[vec[[1]]]},   
     Do[sum[[j]] += vec[[i, j]], {j, 1, len}, {i, 2, Length[vec]}];
     sum], CompilationTarget -> "C"]

(кстати, я считаю, что сравнение скорости более справедливо при компиляции на C, так как виртуальная машина Mathematica, например, гораздо более сильно препятствует вложенным циклам). Эта функция все еще медленнее, чем ваша, но примерно в 3 раза быстрее, чем fun1, для таких маленьких векторов.

Остальная часть неэффективности, я считаю, присуща этому подходу. Тот факт, что вы можете разложить проблему на решение для сумм отдельных компонентов, - это то, что делает вашу вторую функцию эффективной, потому что вы используете структурные операции, такие как Transpose, и, что самое главное, это позволяет выжать больше инструкций из внутреннего цикла, Поскольку это самое главное - вы должны иметь как можно меньше инструкций во внутреннем цикле. Вы можете видеть из CompilePrint, что это действительно так для fun1 vs fun3. В некотором роде вы нашли (для этой проблемы) эффективный высокоуровневый способ вручную развернуть внешний цикл (один по индексу координат). Альтернатива, которую вы предлагаете, попросит компилятор автоматически развернуть внешний цикл, основываясь на дополнительной информации о векторной размерности. Это похоже на правдоподобную оптимизацию, но еще не реализовано для виртуальной машины Mathematica.

Отметим также, что для больших длин векторов (скажем, 20) разница между fun1 и fun2 исчезает, поскольку стоимость выделения/освобождения памяти при тензорном копировании становится несущественной по сравнению со стоимостью массивного присвоения (которая все еще выполняется более эффективно, когда вы назначаете вектор в вектор - возможно, потому, что в этом случае вы можете использовать такие вещи, как memcpy).

В заключение я считаю, что, хотя было бы неплохо, если бы эта автоматизация оптимизации была, по крайней мере, в этом конкретном случае, это своего рода низкоуровневая оптимизация, которую трудно ожидать полностью автоматическим - даже оптимизация компиляторов C не всегда выполняют его. Одна вещь, которую вы могли бы попробовать, - жестко закодировать длину вектора в скомпилированную функцию, затем использовать SymbolicCGenerate (из пакета CCodeGenerator`) для создания символа C, затем использовать ToCCodeString для генерации кода C (или, каким бы то ни было образом вы используете для получения кода C для скомпилированной функции), а затем попытайтесь создать и загрузить библиотеку вручную, включив все оптимизации для компилятора C с помощью параметров CreateLibrary. Будет ли это работать, я не знаю. EDIT Я действительно сомневаюсь, что это вообще поможет, поскольку петли уже реализованы с помощью goto -s для скорости при генерации кода C, и это, вероятно, не позволит компилятору попытаться развернуть цикл.

Ответ 2

Это всегда хороший вариант для поиска функции, которая делает именно то, что вы хотите сделать.

In[50]:= fun3=Compile[{{vec,_Real,2}},Total[vec]]

Out[50]= CompiledFunction[{vec},Total[vec],-CompiledCode-]

In[51]:= Do[fun3[data],{10}]//Timing

Out[51]= {0.121982,Null}

In[52]:= fun3[data]===fun1[data]

Out[52]= True

Другой вариант, менее эффективный (* из-за транспонирования *), заключается в использовании Listable

fun4 = Compile[{{vec, _Real, 1}}, Total[vec], 
  RuntimeAttributes -> {Listable}]

In[63]:= Do[fun4[Transpose[data]],{10}]//Timing

Out[63]= {0.235964,Null}

In[64]:= Do[Transpose[data],{10}]//Timing

Out[64]= {0.133979,Null}

In[65]:= fun4[Transpose[data]]===fun1[data]

Out[65]= True

Ответ 3

Как я могу сказать Compile, что один из его аргументов не является произвольной матрицей n * m, а специальным n * 2 или n * 3, поэтому он может делать эту оптимизацию автоматически, а не использовать сложение общего вектора функция для добавления крошечных векторов длины-2 или длины-3?

Вы пытаетесь выручить "Титаник" с помощью ложки!

На моей машине с Mathematica 7.0.1 ваш первый пример занимает 4 с, вторая занимает 2 секунды, а Total - 0,1 с. Таким образом, у вас есть объективно крайне неэффективное решение (Total на 40 × быстрее!) И правильно определили и рассмотрели один из наименее важных вкладов в плохую производительность (что делает его на 2 раза быстрее).

Низкая производительность Mathematica связана с тем, как оценивается код Mathematica. С точки зрения языка программирования термин переписывающий термин Mathematica может быть мощным решением для компьютерной алгебры, но он безумно неэффективен для оценки программ общего назначения. Следовательно, оптимизации, такие как специализирование для 2D-векторов, не будут рассчитываться в Mathematica, как они могут на скомпилированных языках.

Я быстро перешел к оптимизации в Mathematica, написав код для добавления двух векторных элементов за раз:

fun3 = Compile[{{vec, _Real, 2}}, 
  Module[{x = vec[[1]][[1]], y = vec[[1]][[2]]}, 
   Do[x += vec[[i]][[1]]; y += vec[[i]][[2]], {i, 2, Length[vec]}];
   {x, y}]]

Это на самом деле еще медленнее, принимая 5.4s.

Но это математика в худшем случае. Mathematica полезен только тогда, когда:

• Вы действительно не заботитесь о производительности, или

• В стандартной библиотеке Mathematica уже есть функции (например, Total в этом конкретном случае), которые позволяют эффективно решать вашу проблему либо с помощью одного вызова, либо путем написания script для совершения нескольких вызовов.

Как сказал Руэбенко, Mathematica предоставляет встроенную функцию для решения этой проблемы для вас с помощью одного вызова (Total), но это бесполезно в общем случае, о котором вы спрашиваете. Объективно лучшим решением является предотвращение неэффективности ядра Mathematica путем переноса вашей программы на язык, который оценивается более эффективно.

Например, самое наивное возможное решение в F # (скомпилированный язык, который я должен передать, но почти любой другой будет делать) - использовать 2D-массив:

let xs =
  let rand = System.Random()
  Array2D.init 1000000 2 (fun _ _ -> rand.NextDouble())

let sum (xs: float [,]) =
  let total = Array.zeroCreate (xs.GetLength 1)
  for i=0 to xs.GetLength 0 - 1 do
    for j=0 to xs.GetLength 1 - 1 do
      total.[j] <- total.[j] + xs.[i,j]
  total

for i=1 to 10 do
  sum xs |> ignore

Это сразу на 8 раз быстрее, чем ваше быстрое решение! Но подождите, вы можете сделать все возможное, используя систему статического типа через свой собственный 2D-векторный тип:

[<Struct>]
type Vec =
  val x : float
  val y : float
  new(x, y) = {x=x; y=y}

  static member (+) (u: Vec, v: Vec) =
    Vec(u.x+v.x, u.y+v.y)

let xs =
  let rand = System.Random()
  Array.init 1000000 (fun _ -> Vec(rand.NextDouble(), rand.NextDouble()))

let sum (xs: Vec []) =
  let mutable u = Vec(0.0, 0.0)
  for i=1 to xs.Length-1 do
    u <- u + xs.[i]
  u

Это решение занимает всего 0,057, что на 70% выше, чем у вашего оригинала и значительно быстрее, чем любое решение на базе Mathematica, размещенное здесь до сих пор! Язык, скомпилированный для эффективного кода SSE, скорее всего, будет намного лучше.

Вы можете подумать, что 70 × - особый случай, но я видел много примеров, когда перенос кода Mathematica на другие языки давал огромные ускорения:

• Саль Мангано "производительность критическая" pricer для американских опций получила ускорение 960 × от портирования на F #.

• Sal Mangano красно-черные деревья в Mathematica получил 100-кратное ускорение от портирования на F #.

• Мой трассировщик лучей в Mathematica получил 100 000 × ускорение от портирования на OCaml. Даниэль Лихтблау из Wolfram Research позже оптимизировал Mathematica, но даже его версия была на 1000 × медленнее, чем мой OCaml.