Friday, September 13, 2013

Оптимизация 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 использованием оптимизированной суммы.

Платформа

ProcessorIntel i5-4570, Haswell
Instruction set   AVX2
SystemUbuntu 13.10
CompilerGDC-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. Однако, при оптимизации тривиального алгоритма можно достичь сопостовимой производительности, что позволяет оптимизировать смежные версии алгоритмов.

Saturday, August 17, 2013

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