atmosphere_gm. Российский анонс библиотеки быстрых сеточных методов разделения смесей вероятностных распределений.

Всем добрый день и с Наступившим!

Библиотека

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

Библиотека разделена на пакеты:

1. kernel - весь BLAS и подобный API, который можно оптимизировать. Также утилиты работы с многомерной кластеризацией на гиперторах и многомерными индексами, поиском корней функции одной переменной.
2. fractal - описание фракталов.
3. cyclone - алгоритмы независимые от платформы
1. phiallocator - максимально оптимизированный и гибкий модуль управления памятью при поточном получении данных. Оптимально работает в том числе с данными, которые поступают неравномерными порциями. Данный пункт очень важен: количество перемещений не должно превышать количество новых данных, и при этом нужно иметь в памяти плоскую матрицу для оптимальной работы в том числе с ее транспонированной версией. Управление детерминировано по времени и может использоваться в программах реального времени.
2. generator - итератор по фракталам
3. eye - модуль структуры одной итерации. Работает с информацией с предыдущих шагов при принятии решений (то, что нам известно с предыдущего шага мы не пересчитываем, что во многом зависит от обстоятельств работы и аргументов итерации). В нем сосредоточена вся логика алгоритма. В зависимости от того, как им пользоваться можно получить совершенно разные алгоритмы, в том числе родительский алгоритм Алексея Назарова.
4. cyclone - примеры управления структурой одной итерации. В настоящий момент сделан двухпараметрический пример для смесей нормальных законов. Всего 150 строчек задают все параметры, весь алгоритм от чтения файла до полного управления структурой одной итерации и параметров сетки, и описания отображения собственно нормального распределения на параметрическую сеть.

BLAS API стандартный. Возможна работы без BLAS вовсе.

Технические возможности оптимизации и использования API

1. Библиотека написана на языке программирования D (компиляторы dmd (dmc-32 или win64), gdc (с версии 4.9 часть gcc), ldc (llvm компилятор, включен в анонс llvm 3.4)). Язык имеет С прямой API - можно интегрировать с openBLAS, cuBLAS, Intel MKL, MATLAB, Julia, R и т.д.
2. Возможно использовать библиотеку в системах реального времени. Память управляется самим алгоритмом и не участвует в сборке мусора, поэтому сборщик можно вообще отключить.
3. В плане оптимизации кода компиляторами библиотека полностью эквивалента аналогичной программе написанной на чистом С. Возможно использовать SIMD или встроеный ассемблер. Весь код, который может быть подвержен оптимизации на уровне железа вынесен в пакет kernel.
4. Библиотека платформнонезависима. Сейчас работает на всех основных ОС для ПК или серверов. Потом добавлю поддержку смартфонов :-)
5. Вы можете задать свою функцию правдоподобия в модуле kernel, для этого также нужно сделать согласованную версию кластеризации. Изменение алгоритма при этом не потребуется. По умолчанию алгоритм можно запускать с двумя версиями функции правдоподобия: логарифм произведения сумм (быстрая), или сумма логарифмов сумм (всегда точная).

Стабильность

Программа стабильна и протестирована до уровня RC (release candidate). Сама структура библиотеки располагает к быстрому обнаружению ошибок. Для тестирования использовалась специально созданная система визуализации. В прямом смысле алгоритм тестировался и визуализировался на каждом шаге внутри итерации. При этом API вероятно будет немного меняться согласно пожеланиям пользователей.

Сроки

1. Препринт по алгоритму будет в середине февраля. С этого момента с удовольствием начну совместную работу (если нужно раньше, прошу связаться со мной). В публикации будут использоваться как можно более заезженные данные, ее основная цель - описание алгоритма, а не выводов. Был бы рад иметь после ряд совместных статей с актуальными данными и уже детальным анализом.
2. Исходный код будет опубликован под свободной лицензией на github (рассчитываю на совместную поддержку) как только Виктор Юрьевич примет препринт для публикации.
3. В течении двух недель с момента публикации кода, он будет детально прокомментирован на английском.

Дополнительно

• Разработана библиотека на основе openGL для визуализации данных в реальном времени. Сейчас 2D, будет и 3D к концу учебного года (к примеру для обобщенных гамма, картинки обещают быть красивыми).
• Биндинг blas к D.

TODO

1. Поддержка для kernel различных драйверов:
1. Intel MKL (в том числе SIMD реализация функции логарифма для ФП и других)
2. cuBLAS
3. обобщенная реализация SIMD, в том числе AVX2 и AVX3.
2. 3D визуализация

Оптимизация EM-алгоритмов для вероятностного тематического моделирования на примере PLSA-EM-алгоритма

Как правило тестирование EM алгоритмов затратно по времени. Данная статья призвана обозначить основные способы оптимизации и готовые инструменты, применимые к широкому спектру алгоритмов, работающих с матрицами на примере PLSA-EM алгоритма. 

Матричная форма

Последующее изложение происходит в терминологии "Модификации EM-алгоритма для вероятностного тематического моделирования" К. В. Воронцов, А. А. Потапенко. В ней вы можете ознакомиться с постановкой задачи и с исходным алгоритмом, который далее будет представлен в матричной форме.

Матричное представление удобно по следующим причинам:

1. Простота формулировки, записи и программной реализации. 
2. Алгоритмы работы с матрицами хорошо изучены и оптимизированы.

Основные оптимизации

SIMD

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

Multithreading

Использование преимущества многоядерных процессоров само по себе довольно просто. В тестах оптимизированного тривиального алгоритма перемножения матриц, приведенных в данной статье, используется группировка перемножения матриц по внешнему циклу. Это позволяет создавать в 2 раза меньше новых потоков.

BLAS 

BLAS (Basic Linear Algebra Subprograms) является наиболее распространенным API для работы с матрицами, его реализации доступны как под свободными, так и под коммерческими лицензиями. BLAS разработано и оптимизировано для широкого спектра устройств, начиная со смартфонов и заканчивая суперкомпьютерами. Существуют реализации для графических ускорителей, их тесты можно найти в работе "BLAS Comparison on FPGA, CPU and GPU" Srinidhi Kestury, John D. Davis, Oliver Williams.

OpenBLAS

Наиболее удачным и производительным вариантом для настольных компьютеров (и ноутбуков) является реализация OpenBLAS, являющаяся потомком gotoBLAS2, разработанной Kazushige Goto . В основе реализации, кроме алгоритмов быстрого умножения матриц, использования SIMD инструкций и многопоточности, лежит оптимизация алгоритмов с учетом особенностей кеширования памяти для конкретной процесcорной архитектуры "Anatomy of high-performance matrix multiplication", "High-performance implementation of the level-3 BLAS" Kazushige Goto. Кроме того, ядро OpenBLAS имеет различные драйвера в зависимости от версии процессора. OpenBLAS опубликовано под лицензией BSD.

Оболочки для BLAS

MATLAB

MATLAB по умолчанию использует BLAS, представленный в Intel MKL. Существует реализация BLAS от AMD - ACML, ее также возможно подключить к MATLAB. При этом нужно обратить особое внимание на механизмы реализации многопоточности в конкретной версии MATLAB и ACML (зависит от компилятора).

R

Вероятно самой гибкой оболочкой для BLAS является язык R. В работе "Benchmarking Single- and Multi-Core BLAS Implementations and GPUs for use with R" Dirk Eddelbuettel можно найти сравнение реализаций GotoBLAS, Intel MKL, ATLAS и других. Тут стоит отметить, что, чем производительнее BLAS, тем больше оптимизации в ней сделано под конкретные процессоры.

Julia

Julia - новый язык с JIT компилятором на основе LLVM. В качестве ядра использует OpenBLAS. Пройдя по ссылке можно увидеть таблицу (приведена ниже) сравнения производительноcти, в том числе с MATLAB.

Тестирование

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

1. Тривиальная реализация на основе матричной формулировки.
2. Многопоточная тривиальная реализация на основе матричной формулировки, с оптимизированным скалярным произведением, суммой (используется при нормировании).
3. BLAS
4. BLAS c использованием оптимизированной суммы.

Платформа

Processor	Intel i5-4570, Haswell
Instruction set	AVX2
System	Ubuntu 13.10
Compiler	GDC-4.8.1
Compiler flags	-march=native -fno-bounds-check -frename-registers -frelease -O3

Реализация

Структура матрицы (язык D)

Для чистоты эксперимента здесь представлена простая плоская реализация двумерной матрицы. Производительность данной реализации равна аналогичной на C++.

Тривиальная реализация

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

Оптимизированная тривиальная реализация

Оптимизированная тривиальная реализация абсолютна аналогична тривиальной, за исключением:

1. Используется оптимизированное скалярное произведение
2. Используется распараллеливание по каждому внешнему циклу.

BLAS реализация

Используется многопоточный OpenBLAS 2.8, оптимизированный под архитектуру Intel Haswell.

BLAS + fast sum реализация

Аналогична BLAS, но используется оптимизированное суммирование.

Результаты

Было проведено 64 теста с различным количеством тем, документов и слов для каждой реализации. Результаты приведены в следующей таблице:
Отдельно вынесены 2 графика

Заключение

BLAS является достаточно мощным иструментом, с которым можно работать как нативно через компилируемые языки, так и через оболочки, такие как MATLAB, R и Julia. Однако, при оптимизации тривиального алгоритма можно достичь сопостовимой производительности, что позволяет оптимизировать смежные версии алгоритмов.

SIMD implementation of dot-product.Benchmarks

Dot product can be auto-vectorized only with fast-math or similar compiler option. fast-math allows compiler to use associative floating point transformations. Other math functions like exponent can be damaged consequently.
Possible solutions:

1. Compile fast_math code from other program separately and then link it. This is easy solution. However this is a step back to C.
2. To introduce a @fast_math attribute. This is hard to realize. But I hope this will be done for future compilers.
3. Do vectorization yourself. In that case you need to realize SIMD accessory functions like unaligned load.

So I wrote implementations of dot product for real and complex numbers.

Code

Dot product for real numbers:

Complex version is similar. The main loop:

There are gdc and ldc implementations, but ldc implementation was not tested. Source code is available at GitHub.

Benchmarks

Processor	Intel i5-4570, Haswell
Instruction set	AVX2
System	Ubuntu 13.10
Compiler	GDC-4.8.1
Compiler flags	-march=native -fno-bounds-check -frename-registers -frelease -O3