Эффективное вычисление массива в Haskell

Потому что всем нам нужен высокопроизводительный и красивый код.

Цель этого сообщения в блоге - объяснить мои результаты в исследовании интерфейса линейной алгебры в Haskell. Это исследование было проведено с целью проинструктировать нас в theam, где искать высокопроизводительные операции линейной алгебры, и дать нам набросок пространства решений в этой области. Это ни в коем случае не полное расследование, и я полностью открыт для мысли, что, возможно, неправильно оценил некоторые из выявленных мной проблем, на самом деле это то, что я надеюсь спровоцировать, разместив это там, поскольку это просто означало бы, что эта проблема решить проще, чем я предполагал.

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

Некоторые прошлые исследования

Когда говорят об эффективных алгоритмах линейной алгебры, часто говорят о Базовых программах линейной алгебры (BLAS), и в некотором смысле этот набор дополнительных интерфейсов стал синонимом из-за количества часов, которые были потрачены вложить в оптимизацию различных процедур и алгоритмов, составляющих технический интерфейс BLAS. Первоначальная идея BLAS состоит в том, что существует небольшой набор операций, которые можно в значительной степени оптимизировать, используя знания о современных компьютерах, и большинство других полезных программ могут быть написаны на этих операциях. Эти исходные операции называются уровнем 1 [2] и состоят из векторных операций, таких как суммирование, норма, максимум, скалярное произведение и т. Д. Эти операции оптимизируются в первую очередь с помощью векторных процессоров. Со временем, когда многопроцессорность стала очень важной, интерфейс уровня 2 [5], который семантически мог быть реализован в терминах уровня 1, был создан для матрично-векторной такие операции, как матрично-векторные произведения, обновления, решение линейных систем. Наконец, был создан уровень 3 [6], использующий преимущества современных архитектур памяти для реализации матрично-матричных операций, таких как матричное умножение и обновления. Все эти оптимизации должны использоваться любым серьезным конкурентом для систем, основанных на BLAS или реализации BLAS, используемой под капотом. Один из принципов дизайна, который меня особенно вдохновил, - это идея денотационной экономии, поскольку, как только вы поймете n-й уровень операций, вы сможете понять (n + 1) -й, реализовав (n + 1) ) th в примитивах n-го уровня, даже несмотря на то, что (n + 1) th реализованы эффективными способами, которые вам не нужно понимать как пользователю.

Одна из точек разработки в области реализации BLAS называется Автоматически настраиваемое программное обеспечение линейной алгебры (ATLAS), в котором используется метод под названием Автоматическая эмпирическая оптимизация программного обеспечения. По сути, поиск выполняется после установки по ряду параметров, связанных с производительностью, обоснованием этого является то, что ряд применяемых оптимизаций по существу будет работать только для некоторых параметров, которые зависят от компьютера пользователя. Если какое-либо программное обеспечение может найти эти параметры для себя, вместо того, чтобы полагаться на пользователя, чтобы указать, какие из них являются лучшими, программист может написать программное обеспечение, которое может работать быстро даже на будущем оборудовании, поскольку константы будут меняться в зависимости от машины. система скомпилирована для. Требования для поддержки этого программного метода изложены в [7] и заключаются в следующем:

  1. Изоляция подпрограмм, критичных к производительности
  2. Метод адаптации программного обеспечения. Текущие подходы к разным условиям
  3. Надежные контекстно-зависимые таймеры
  4. Соответствующая эвристика поиска

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

Eigen - это библиотека C ++, которая обеспечивает гораздо лучшую нотацию, чем BLAS, поскольку ее интерфейс абстрагирован от модели памяти, и можно просто написать такой код, как :

Matrix2d x = a * (b + c) - d * a;

Прелесть Eigen заключается в том, что приведенный выше код будет переведен в форму, в которой каждая операция в правой части присваивания не будет выделять новую матрицу, даже если это их семантика. Это возможно в C ++ относительно простым способом, но, как мы увидим в Haskell, мы должны приложить немного больше усилий для слияния, поскольку у нас нет такого прославленного понятия присваивания.

Некоторые современные подходы

Чтобы говорить о текущих подходах к этой проблеме, мы должны понимать, что мы хотим делать. Цель этого документа - изучить пространство решений для эффективных численных вычислений в Haskell, и поэтому нам необходимо изучить три области:

  1. Центральные процессоры (ЦП)
  2. Графические процессоры (ГП)
  3. Вычислительные кластеры

У каждой из перечисленных выше парадигм есть свои применения, и они могут использоваться в тандеме для большего эффекта, но тот, который предлагает самую высокую пропускную способность для таких операций, как те, которые нам нужны, - это GPU, в виде сотен или Тысячи ядер трудно превзойти с восемью или шестнадцатью. При этом я бы сказал, что идеальная ситуация - это та ситуация, в которой мы оставались безразличными к этому выбору, поэтому я обратился к предметно-ориентированным языкам.

Первое встречное решение для численных вычислений в Haskell - это библиотека vector. Vector - это, по его собственным словам, эффективная реализация массивов с индексом Int (как изменяемых, так и неизменяемых) с мощной структурой оптимизации цикла. Пакет структурирован в терминах класса Vector и экземпляров класса в распакованных, сохраненных извне, упакованных в коробку, а также в примитивных векторах. Эта библиотека может использоваться для численных вычислений и является подходящим решением для многих более простых задач, особенно с учетом того, что библиотека статистика предлагает множество дополнительных функций, которые можно использовать для вычисления различных фактов. о ваших данных. Здесь я нахожу коэффициент корреляции Пирсона между двумя популяциями образцов:

import Data.Vector.Generic
import Statistics.Correlation 
import DataFetchingLib (getData)
main = do
  xs <- getData "Population1"
  ys <- getData "Population2"
  print (pearson (zip xs ys))

Вот еще один пример использования пакета math-functions, от которого зависит статистика.

import GHC.Exts (fromList)
import Data.Vector.Generic (Vector(..))
import qualified Data.Vector.Unboxed as Unboxed
import Numeric.RootFinding
import Numeric.Polynomial
import Numeric.MathFunctions.Constants (m_epsilon)
wasFound :: Root a -> Bool
wasFound (Root _) = True
wasFound _ = False
main :: IO ()
main = do
  -- f(x) = 1 + x^1 + x^2 + ... + x^100
  let coeffs :: Unboxed.Vector Double = fromList [x | x <- [1..100]]
  let poly x = evaluatePolynomial x coeffs
  let roots = [ridders (m_epsilon * 4) (x, x + 0.2) poly
              | x <- [-1000.0, -999.8 .. 1000.0]]
  print (wasFound `filter` roots)

Другое решение с некоторыми очевидными преимуществами - hmatrix, которое просто предоставляет набор операций с поддержкой BLAS / LAPACK на плотных матрицах. Одним из больших преимуществ этой библиотеки является пакет Numeric.LinearAlgebra.Static, который предлагает векторы и матрицы статического размера, так что вы можете исключить большой класс ошибок при работе с такими программами. Эта библиотека - явный выигрыш, если у вас есть множество сложных задач линейной алгебры в вашем приложении, и вы можете себе это позволить на процессоре.

Еще одно рассмотренное мной решение, гораздо более подходящее для произвольного анализа данных, - repa, оно было написано как результат [9]. Библиотека, опять же, своими словами, обеспечивает высокопроизводительные, регулярные, многомерные, полиморфные параллельные массивы формы, гарантируя, что числовые данные хранятся в распакованном виде. Он имеет структуру распараллеливания, которая по сути является автоматической и по этой причине может создавать невероятно эффективный код. Существует плагин GHC, repa-plugin, который обеспечивает слияние потоков данных, что, по сути, нам понравилось в Eigen. библиотека, возможность не выделять столько памяти, но при этом использовать удобную нотацию. Это достигается только с помощью repa , а также путем использования отложенного представления наряду с обычным представлением массива, поскольку это позволяет вам отложить создание массивов для потребителя любого красивого массива, который вы Оформил в отложенном конструкторе. Эта библиотека, похоже, специализируется на обработке изображений и имеет функции, которые будут поддерживать это очень хорошо.

Последнее решение, которое я рассмотрел, было ускорение, предметно-ориентированный язык (DSL) для высокопроизводительных вычислений в Haskell, который стремится стать хорошо подходящим для GPU исполнение [10]. Он использует многие из тех же методов, что и repa, но он подходит для выполнения в очень многих контекстах, а accelerate-llvm предлагает многоядерные, а также серверные части на GPU. Явным преимуществом ускорения является отсутствие контекста: это глубокое внедрение, а repa - неглубокое. Это означает, что, в то время как массивы repa представлены непосредственно своей семантикой, вычисления ускорения хранятся в форме абстрактного синтаксического дерева (AST), что означает их можно запускать в разных контекстах. Для меня это кажется большим преимуществом, так как можно написать код анализа данных и выполнить статический анализ, а также запустить его очень удобным способом.

Обсуждение

Сначала я посмотрел на старые решения на C, Fortran и C ++ и на то, что люди веками использовали для решения этих проблем, затем я посмотрел на современные решения Haskell и попытался понять, что изменилось, чтобы сделать эти новые решения более желательны сейчас, и как мы можем их улучшить или использовать как есть. Первое, что я заметил, это то, что в BLAS, который на самом деле все еще лежит в основе многих функций линейной алгебры, которые вы собираетесь выполнять на компьютере, был написан в то время, когда ограничения были очень сильными. другой. Первоначальная статья о BLAS была
опубликована в 70-х годах, когда память можно было измерять в килобайтах или мегабайтах, тогда как сейчас ее измеряют в гигабайтах и ​​терабайтах. Понятно, что модель изменилась и что есть больше возможностей для кода на более высоком уровне, но сколько там места?

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

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main() 
{
  Matrix2d a, b, c;
a << 1, 2,
     3, 4;
b << 2, 3,
     1, 4;
c = a + b - a - b + a;
}

Этот код очень удобен для чтения и, согласно документации Eigen, должен компилироваться до кода, который выделяет только три матрицы. Это то, что нам нужно имитировать в любом решении, которое мы выберем, поскольку мы должны до некоторой степени забыть о памяти, если мы хотим работать на высоком уровне, но мы, безусловно, хотим, чтобы за кулисами эти вещи достигли хоть немного оптимизирован. Проблема, конечно, в Haskell заключается в том, что здесь нет прославленного понятия присваивания, и вместо этого нам приходится иметь дело с чистыми значениями и функциями. Решение, как заявлено в repa и ускорении, состоит в том, чтобы адаптировать своего рода DSL, который скрывает сложность этого позади вас. В repa это делается с помощью «отложенного» матричного представления, представленного функцией от индекса к значению, так что вы можете отложить построение сложных матриц до потребителя матрицы. На мой взгляд, путь, который выбирает ускорение, намного проще и расширяемо, по сути, просто используя AST за кулисами для разумной оптимизации времени и пространства.

Когда мы говорим о необходимых параметрах памяти и времени, легко думать, что эти ограничения можно решить с помощью экспериментов на одной машине, но мне ясно, что это уже не так. Методы разработки алгоритмов часто проектируют один алгоритм для одной машины, но ATLAS немного отличался, у него было пространство реализаций, которое, как было известно, содержало достойные реализации на некоторых машинах, и ему было разрешено покачиваться в этом пространстве и рисовать. автоматически выводит то, что ему лучше всего делать [7]. Существуют интересные оптимизации, которые могут быть выполнены на основе прогнозирования производительности кеша [8], а также константы, доступные во время выполнения [9], которые, возможно, мы сможем изучить с такой платформы, как ускорение или другой глубоко встроенному DSL нравится, так как в этом параметре у нас есть доступ к синтаксическому представлению во время выполнения. Это может быть интересным начальным приложением машинного обучения в Haskell и мотивирующим случаем для предстоящих исследований интерфейсов машинного обучения.

Цели

Я думаю, что в краткосрочной перспективе будет выгодно использовать ускорение и набор бэкэндов, представленных в акселерат-llvm. Одно очевидное преимущество заключается в том, что она разрабатывается членом сообщества DataHaskell, а это означает, что автор заинтересован в расширении этой библиотеки для случаев использования науки о данных. Эта стратегия позволит нам развернуть высокопроизводительный анализ данных с ускорением на GPU на облачных платформах, таких как AWS и (довольно скоро) Google Cloud Platform.

Ясно, что в долгосрочной перспективе нам понадобится создать карту этой территории, и я не претендую на понимание всего пространства дизайна каким-либо образом. Я думаю, что было бы здорово иметь исчерпывающие и действенные знания о производительности vector, hmatrix, repa, ускорения, а также других привязок, таких как Eigen, выпущенные для Haskell. При этом ускорение и repa - единственные платформы, которые поддерживают произвольные тензоры. Мне ясно, что на данный момент ускорение превосходит repa, если выполняется на графическом процессоре, хотя repa имеет преимущество. что он полагается на неглубокое встраивание и, таким образом, позволяет оптимизатору GHC Haskell творить с ним чудеса и, таким образом, может быть смешан с произвольным кодом с хорошими результатами, тогда как ускорение полагается на более глубокое встраивание и, таким образом, сохраняет больше информации при переводе и компиляции исходного кода в различные серверные части.

Я также хочу упомянуть, что введение линейных типов в язык программирования Haskell, реализованное в GHC, дало бы нам отличные способы устранения некоторых основных проблем, возникающих из-за сложности памяти. Похоже, это может произойти в будущем [3]. Если кто-то тем временем хочет выяснить, как это может происходить, я рекомендую изучить ржавчину.

Источники

0. Веб-сайт репы: http://repa.ouroborus.net/
1. Репозиторий ускоренного github: https://github.com/AccelerateHS/accelerate
2. Базовый линейный Подпрограммы алгебры для использования Fortran. C.L. Лоусон, Р.Дж. Хэнсон, Д. Кинкейд, Ф. Крог: http://www.cs.utexas.edu/users/kincaid/blas.pdf
3. Описание линейных типов GHC: https://ghc.haskell.org/trac/ghc/wiki / LinearTypes
4. repa-plugin: https://hackage.haskell.org/package/repa-plugin
5. Расширенный набор подпрограмм базовой линейной алгебры Fortran. Джек Дж. Донгарра, Джерми Дю Кроз, Свен Хаммарлинг, Ричард Дж. Хэнсон.
6. Набор подпрограмм базовой линейной алгебры уровня 3. Джек Дж. Донгарра, Джереми Дю Кроз, Свен Хаммарлинг, Иэн Дафф.
7. Автоматически настраиваемое программное обеспечение линейной алгебры. Р. Клинт Уэйли, Джек Дж. Донгарра.
8. Оптимизация программы на основе прогнозирования производительности кэша во время компиляции. Уэсли К. Каплов, Болеслав К. Шимански.
9. Регулярные, шейп-полиморфные, параллельные массивы в Haskell. Габриэле Келлер, Мануэль М. Т. Чакраварти, Роман Лещинский, Саймон Пейтон Джонс, Бен Липпмейр.
10. Оптимизация чисто функциональных программ на GPU. Тревор Л. Макдоннелл.