Достаточно эффективный чистофункциональный матричный продукт в Haskell?

Я знаю Haskell немного, и мне интересно, можно ли написать что-то вроде матрично-матричного продукта в Haskell, что является следующим:

Чистофункциональные: нет IO или State монады в своей сигнатуре типа (мне все равно, что происходит в теле функции. То есть меня не волнует, функция body использует монады, пока вся функция чиста). Я могу использовать этот матрично-матричный продукт в чистой функции.
Память: нет malloc или указателей. Я знаю, что можно "написать C" в Haskell, но вы теряете безопасность памяти. Фактически запись этого кода на C и его взаимодействие с Haskell также лишена безопасности памяти.
Насколько эффективен, например, Java.. Для конкретности позвольте предположить, что я говорю о простой тройной петле, одинарной точности, непрерывной компоновке столбцов (float[], not float[][]) и матрицы размером 1000x1000 и одноядерный процессор. (Если вы получаете 0,5-2 операции с плавающей запятой за цикл, вы, вероятно, находитесь в шаге.)

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

Я уже знаю, что

Реализация тройного цикла не самая эффективная. Это довольно кэш-невообразимо. Лучше использовать хорошо написанную BLAS в этом конкретном случае. Однако нельзя всегда рассчитывать на доступность библиотеки C для того, что вы пытаетесь сделать. Интересно, может ли быть достаточно эффективный код в обычном Haskell.
Некоторые люди написали целые исследовательские работы, демонстрирующие № 3. Однако я не научный сотрудник. Интересно, возможно ли сохранить простые вещи в Haskell.
Нежное введение в Haskell имеет реализацию матричного продукта. Однако он не удовлетворяет вышеуказанным требованиям.

Адресация комментариев:

У меня есть три причины: во-первых, требование "no malloc or pointers" пока еще не определен (я призываю вас написать любую часть Haskell код, который не использует указатели);

Я видел много программ Haskell, не использующих Ptr. Возможно, это относится к тому, что на уровне машинной инструкции будут использоваться указатели? Это не то, что я имел в виду. Я имел в виду уровень абстракции исходного кода Haskell.

во-вторых, атака на исследования CS неуместна (и, кроме того, я не может представить ничего проще, чем использовать код, который кто-то еще уже написано для вас); в-третьих, существует множество матричных пакетов на Hackage (и подготовка к заданию этого вопроса должна включать рассмотрение и отклонение каждого).

Кажется, что ваши # 2 и # 3 одинаковы ( "использовать существующие библиотеки" ). Меня интересует матричный продукт как простой тест того, что Haskell может сделать сам по себе, и позволяет ли он "упростить простые вещи". Я мог бы легко найти числовую проблему, у которой нет готовых библиотек, но тогда я должен был бы объяснить проблему, тогда как все уже знают, что такое матричный продукт.

Как Java может удовлетворять 1.? Любой Java-метод по существу :: IORef Arg -> ... -> IORef This -> IO Ret

Это относится к корню моего вопроса, фактически (+1). Хотя Java не претендует на отслеживание чистоты, Haskell делает. В Java, является ли функция чистой или нет, указывается в комментариях. Я могу утверждать, что матричный продукт чист, хотя я и мутация в теле функции. Вопрос в том, совместим ли подход Haskell (чистота, закодированный в системе типов) с эффективностью, безопасностью и простотой.

Ответ 1

Насколько эффективен, например, Java. Для конкретности позвольте предположить, что я говорю о простой тройной петле, одинарной точности, непрерывной компоновке столбцов (float [], not float [] []) и матрицах размером 1000x1000 и одноядерном процессоре. (Если вы получаете 0,5-2 операции с плавающей запятой за цикл, вы, вероятно, находитесь в шаге)

Так что-то вроде

public class MatrixProd {
    static float[] matProd(float[] a, int ra, int ca, float[] b, int rb, int cb) {
        if (ca != rb) {
            throw new IllegalArgumentException("Matrices not fitting");
        }
        float[] c = new float[ra*cb];
        for(int i = 0; i < ra; ++i) {
            for(int j = 0; j < cb; ++j) {
                float sum = 0;
                for(int k = 0; k < ca; ++k) {
                    sum += a[i*ca+k]*b[k*cb+j];
                }
                c[i*cb+j] = sum;
            }
        }
        return c;
    }

    static float[] mkMat(int rs, int cs, float x, float d) {
        float[] arr = new float[rs*cs];
        for(int i = 0; i < rs; ++i) {
            for(int j = 0; j < cs; ++j) {
                arr[i*cs+j] = x;
                x += d;
            }
        }
        return arr;
    }

    public static void main(String[] args) {
        int sz = 100;
        float strt = -32, del = 0.0625f;
        if (args.length > 0) {
            sz = Integer.parseInt(args[0]);
        }
        if (args.length > 1) {
            strt = Float.parseFloat(args[1]);
        }
        if (args.length > 2) {
            del = Float.parseFloat(args[2]);
        }
        float[] a = mkMat(sz,sz,strt,del);
        float[] b = mkMat(sz,sz,strt-16,del);
        System.out.println(a[sz*sz-1]);
        System.out.println(b[sz*sz-1]);
        long t0 = System.currentTimeMillis();
        float[] c = matProd(a,sz,sz,b,sz,sz);
        System.out.println(c[sz*sz-1]);
        long t1 = System.currentTimeMillis();
        double dur = (t1-t0)*1e-3;
        System.out.println(dur);
    }
}

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

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

Я не хочу, чтобы это звучало как вызов, , но обратите внимание, что Java может полностью удовлетворить все перечисленные выше

И то, что Java дает легко, поэтому я возьму это.

(Если вы получаете 0,5-2 операции с плавающей запятой за цикл, вы, вероятно, находитесь в шаге)

Нигде рядом, я боюсь, ни на Java, ни в Haskell. Слишком много недостатков кэша, чтобы достичь этой пропускной способности с помощью простого тройного цикла.

Выполняя то же самое в Haskell, снова не думая о том, чтобы быть умным, простой простой тройной цикл:

{-# LANGUAGE BangPatterns #-}
module MatProd where

import Data.Array.ST
import Data.Array.Unboxed

matProd :: UArray Int Float -> Int -> Int -> UArray Int Float -> Int -> Int -> UArray Int Float
matProd a ra ca b rb cb =
    let (al,ah)     = bounds a
        (bl,bh)     = bounds b
        {-# INLINE getA #-}
        getA i j    = a!(i*ca + j)
        {-# INLINE getB #-}
        getB i j    = b!(i*cb + j)
        {-# INLINE idx #-}
        idx i j     = i*cb + j
    in if al /= 0 || ah+1 /= ra*ca || bl /= 0 || bh+1 /= rb*cb || ca /= rb
         then error $ "Matrices not fitting: " ++ show (ra,ca,al,ah,rb,cb,bl,bh)
         else runSTUArray $ do
            arr <- newArray (0,ra*cb-1) 0
            let outer i j
                    | ra <= i   = return arr
                    | cb <= j   = outer (i+1) 0
                    | otherwise = do
                        !x <- inner i j 0 0
                        writeArray arr (idx i j) x
                        outer i (j+1)
                inner i j k !y
                    | ca <= k   = return y
                    | otherwise = inner i j (k+1) (y + getA i k * getB k j)
            outer 0 0

mkMat :: Int -> Int -> Float -> Float -> UArray Int Float
mkMat rs cs x d = runSTUArray $ do
    let !r = rs - 1
        !c = cs - 1
        {-# INLINE idx #-}
        idx i j = cs*i + j
    arr <- newArray (0,rs*cs-1) 0
    let outer i j y
            | r < i     = return arr
            | c < j     = outer (i+1) 0 y
            | otherwise = do
                writeArray arr (idx i j) y
                outer i (j+1) (y + d)
    outer 0 0 x

и вызывающий модуль

module Main (main) where

import System.Environment (getArgs)
import Data.Array.Unboxed

import System.CPUTime
import Text.Printf

import MatProd

main :: IO ()
main = do
    args <- getArgs
    let (sz, strt, del) = case args of
                            (a:b:c:_) -> (read a, read b, read c)
                            (a:b:_)   -> (read a, read b, 0.0625)
                            (a:_)     -> (read a, -32, 0.0625)
                            _         -> (100, -32, 0.0625)
        a = mkMat sz sz strt del
        b = mkMat sz sz (strt - 16) del
    print (a!(sz*sz-1))
    print (b!(sz*sz-1))
    t0 <- getCPUTime
    let c = matProd a sz sz b sz sz
    print $ c!(sz*sz-1)
    t1 <- getCPUTime
    printf "%.6f\n" (fromInteger (t1-t0)*1e-12 :: Double)

Итак, мы делаем почти то же самое на обоих языках. Скомпилируйте Haskell с помощью -O2, Java с javac

$ java MatrixProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592897E10
8.193
$ ./vmmult 1000 "-13.7" 0.013
12915.623
12899.999
8.35929e10
8.558699

И получившиеся времена довольно близки.

И если мы скомпилируем Java-код для native, с gcj -O3 -Wall -Wextra --main=MatrixProd -fno-bounds-check -fno-store-check -o jmatProd MatrixProd.java,

$ ./jmatProd 1000 "-13.7" 0.013
12915.623
12899.999
8.3592896512E10
8.215

до сих пор нет большой разницы.

В качестве специального бонуса тот же алгоритм в C (gcc-O3):

$ ./cmatProd 1000 "-13.7" 0.013
12915.623047
12899.999023
8.35929e+10
8.079759

Таким образом, это не показывает принципиальной разницы между простой Java и простым Haskell, когда дело доходит до вычислительно-интенсивных задач с использованием чисел с плавающей запятой (при работе с целочисленной арифметикой на средних и больших числах использование GMP GHC делает Haskell опережающим Java BigInteger огромный запас для многих задач, но это, конечно, проблема библиотеки, а не язык), и оба они близки к C с помощью этого алгоритма.

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

Если мы улучшаем шаблон доступа, умножив матрицу основных строк на матрицу со столбцом, все станет быстрее, gcc-скомпилированный C заканчивает его 1.18s, java занимает 1.23s, а собранный ghc Haskell занимает около 5.8 с, что может быть сокращено до 3 секунд, используя бэкэнд llvm.

Здесь проверка диапазона библиотекой массивов действительно болит. Используя доступ к неконтролируемому массиву (как следует, после проверки ошибок, поскольку проверки уже выполняются в коде, управляющем циклами), родной бэкенд GHC заканчивается в 2.4s, идя через бэкенд llvm, позволяет вычислять завершение в 1.55s, который является достойным, хотя и значительно медленнее, чем C и Java. Используя примитивы из GHC.Prim вместо библиотеки array, бэкэнд llvm создает код, который работает в 1.16 (опять же без проверки границ на каждом доступе, но при этом в процессе вычисления могут быть получены только действительные индексы, которые в этом случае могут быть легко доказаны, поэтому здесь не жертвует безопасность памяти¹; проверка каждого доступа приводит к времени до 1,96 с, что значительно лучше, чем проверка границ array).

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

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

Ответ 2

Есть два угла, чтобы атаковать эту проблему.

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

Однако я не научный сотрудник. Интересно, возможно ли сохранить простые вещи в Haskell.

Если вы укажете эти документы, возможно, мы сможем их расшифровать.
Разработка программного обеспечения, в этих строках, хорошо понятна, понятна и даже проста. Разумный кодер Haskell будет использовать тонкую оболочку вокруг BLAS или искать такую оболочку в Hackage.

Расшифровка новейших исследований - это непрерывный процесс, который переводит знания от исследователей к инженерам. Это был научный сотрудник по информатике C.A.R. Хоар, который сначала обнаружил быструю сортировку и опубликовал статью об этом. Сегодня это редкий выпускник компьютерных наук, который не может лично реализовать quicksort из памяти (по крайней мере, тех, которые недавно закончили).

Бит истории

Почти этот точный вопрос задавался в истории несколько раз раньше.

Можно ли написать матричную арифметику в Fortran, которая выполняется так же быстро, как сборка?
Можно ли написать матричную арифметику в C так же быстро, как Fortran?
Можно ли написать матричную арифметику в Java, которая выполняется так же быстро, как C?
Можно ли написать матричную арифметику в Haskell, которая так же быстро, как Java?

До сих пор ответ всегда был "еще не", а затем "достаточно близко". Достижения, которые делают это возможным, связаны с улучшением написания кода, улучшением компиляторов и улучшением самого языка программирования.

В качестве конкретного примера C не смог превзойти Fortran во многих реальных приложениях, пока компиляторы C99 не получили широкого распространения в последнее десятилетие. Предполагается, что в Fortran разные массивы имеют различное хранилище друг от друга, тогда как в C это обычно не так. Поэтому компиляторам Fortran было разрешено делать оптимизации, которые компиляторы C не могли. Ну, пока C99 не выйдет, и вы можете добавить квалификатор restrict к вашему коду.

Компиляторы Fortran ждали. В конце концов, процессоры стали достаточно сложными, что хорошая сборка была сложнее, и компиляторы стали достаточно сложными, чтобы Fortran был быстрым.

Затем программисты C ждали до 2000-х годов для возможности писать код, сопоставляемый с Fortran. До этого момента они использовали библиотеки, написанные на Fortran или ассемблере (или и то, и другое), или с уменьшенной скоростью.

Java-программистам также приходилось ждать компиляторов JIT и приходилось ждать появления определенных оптимизаций. Компиляторы JIT были первоначально эзотерической концепцией исследования, пока они не стали частью повседневной жизни. Оптимизация проверки границ также необходима для того, чтобы избежать проверки и ветвления для каждого доступа к массиву.

Вернуться к Haskell

Итак, понятно, что программисты Haskell "ждут", как и программисты Java, C и Fortran перед ними. Чего мы ждем?

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

И, возможно, мы ждем некоторой комбинации выше.

О чистоте

Чистота и монады много объединяются в Haskell. Причина этого в том, что в Haskell нечистые функции всегда используют монаду IO. Например, монада State на 100% чиста. Поэтому, когда вы говорите: "чистая" и "сигнатура типа не использует монаду State", это фактически полностью независимые и отдельные требования.

Однако вы также можете использовать монаду IO в реализации чистых функций, и на самом деле это довольно просто:

addSix :: Int -> Int
addSix n = unsafePerformIO $ return (n + 6)

Хорошо, да, эта глупая функция, но она чиста. Это даже очевидно чисто. Тест на чистоту двоякий:

Оказывает ли тот же результат одинаковые данные? Да.
Получают ли какие-либо семантически значимые побочные эффекты? Нет.

Причина, по которой нам нравится чистота, состоит в том, что чистые функции легче создавать и манипулировать, чем нечистые функции. То, как они реализованы, не имеет значения. Я не знаю, знаете ли вы об этом, но Integer и ByteString являются в основном оболочками вокруг нечистых функций C, хотя интерфейс чист. (Там работает над новой реализацией Integer, я не знаю, насколько это далеко.)

Окончательный ответ

Вопрос в том, совместим ли подход Haskell (чистота, закодированный в системе типов) с эффективностью, безопасностью и простотой.

Ответ на эту часть - "да", так как мы можем использовать простые функции из BLAS и помещать их в чистую, безопасную для типа оболочку. Тип оболочки кодирует безопасность функции, хотя компилятор Haskell не может доказать, что реализация функции чиста. Наше использование unsafePerformIO в его реализации является признанием того, что мы доказали чистоту функции, и это также уступка, что мы не могли найти способ выразить это доказательство в системе типа Haskell.

Но ответ еще не "пока", так как я не знаю, как полностью реализовать эту функцию в Haskell как таковой.

Исследования в этой области продолжаются. Люди смотрят на такие системы, как Coq и новые языки, такие как Agda, а также события в самой GHC. Чтобы увидеть, какой тип системы нам нужно будет доказать, что высокопроизводительные процедуры BLAS можно безопасно использовать. Эти инструменты также могут использоваться с другими языками, такими как Java. Например, вы можете написать доказательство в Coq, что ваша реализация Java чиста.

Я приношу свои извинения за ответ "да и нет", но ни один другой ответ не будет признавать как вклад инженеров (которые заботятся о "да" ), так и исследователи (которые еще не заботятся ").

P.S. Просьба привести статьи.

Ответ 3

Как и Java, Haskell не лучший язык для написания числового кода.

Haskell численно-тяжелая когенерация - это... средняя. У него не было много исследований за этим, как у Intel и GCC.

Что Haskell дает вам вместо этого, это способ четко связать ваш "быстрый" код с остальной частью вашего приложения. Помните, что 3% кода отвечает за 97% времени выполнения вашего приложения. 1

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

Но я отвлекаюсь.

Благодаря тому, что эти высоко оптимизированные функции изолированы и с аналогичным интерфейсом с остальной частью вашего кода Haskell, вы можете выполнять оптимизацию на высоком уровне с помощью очень мощных правил перезаписи Haskell, которые позволяют вам писать такие правила, как reverse . reverse == id, который автоматически уменьшить сложные выражения во время компиляции 2. Это приводит к чрезвычайно быстрым, чисто функциональным и простым в использовании библиотекам, таким как Data.Text 3 и Data.Vector [4].

Объединяя высокие и низкие уровни оптимизации, мы получаем гораздо более оптимизированную реализацию, причем каждая половина ( "C/asm" и "Haskell" ) относительно легко читается. Оптимизация низкого уровня выполняется на его родном языке (C или сборке), оптимизация высокого уровня получает специальный DSL (правила перезаписи Haskell), и остальная часть кода не обращает на это внимания.

В заключение, да, Haskell может быть быстрее, чем Java. Но он обманывает, перейдя через C для сырого FLOPS. Это намного сложнее сделать на Java (а также иметь гораздо более высокие накладные расходы для Java FFI), поэтому этого избежать. В Хаскелле это естественно. Если ваше приложение тратит чрезмерное количество времени на числовые вычисления, то, возможно, вместо того, чтобы смотреть на Haskell или Java, вы смотрите на Fortran для своих нужд. Если ваше приложение тратит большую часть своего времени на крошечную часть кода, чувствительного к производительности, тогда ваш FFI - ваш лучший выбор. Если ваше приложение не проводит какое-либо время в числовом коде... тогда используйте все, что вам нравится. =)

Haskell (и Java, если на то пошло) не является Fortran.

1 Эти цифры были составлены, но вы получите мою точку зрения.

2 http://www.cse.unsw.edu.au/~dons/papers/CLS07.html

3 http://hackage.haskell.org/package/text

[4] http://hackage.haskell.org/package/vector

Теперь, чтобы это было в стороне, чтобы ответить на ваш реальный вопрос:

Нет, в настоящее время нет умных писать умножения на матрицу в Haskell. В настоящий момент REPA - это канонический способ сделать это [5]. Реализация частично нарушает безопасность памяти (они используют unsafeSlice), но "поврежденная безопасность памяти" изолирована от этой функции, на самом деле очень безопасна (но не легко проверена компилятором), и ее легко удалить, если что-то пойдет не так (замените "unsafeSlice" с "slice" ).

Но это Хаскелл! Очень редко характеристики производительности функции выполняются изолированно. Это может быть плохо (в случае утечки пространства) или очень-очень хорошая вещь.

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

Что, если бы вы были ученым с миллионами точек данных и хотели бы размножать огромные матрицы? [7]

Для этих людей у нас есть mmultP [6]. Это выполняет матричное умножение, но является параллельным данными и подвержено вложенным данным REPA parallelism. Также обратите внимание, что код практически не изменяется от последовательной версии.

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

Подобно ByteString, Text и Vector, массивы REPA подлежат слиянию. 2 Кстати, вы должны читать 2 - это очень хорошо написанная статья. Это, в сочетании с агрессивной вставкой соответствующего кода и высокопараллельной природы REPA, позволяет нам выразить эти высокоуровневые математические концепции с очень продвинутыми оптимизациями высокого уровня за кулисами.

Хотя метод написания эффективного матричного умножения в чистых функциональных языках в настоящее время не известен, мы можем приблизиться (нет автоматической векторизации, нескольких чрезмерных разломов, чтобы добраться до фактических данных и т.д.), но ничего близкого что IFORT или GCC. Но программа не существует на острове, и создание острова в целом хорошо работает, намного проще в Haskell, чем в Java.

[5] http://hackage.haskell.org/packages/archive/repa-algorithms/3.2.1.1/doc/html/src/Data-Array-Repa-Algorithms-Matrix.html#mmultS

[6] http://hackage.haskell.org/packages/archive/repa-algorithms/3.2.1.1/doc/html/src/Data-Array-Repa-Algorithms-Matrix.html#mmultP

[7] По сути, лучший способ сделать это - использовать графический процессор. Для Haskell доступно несколько GPU DSL, которые делают это возможным изначально. Они очень аккуратные!