Оптимизация 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++.

	struct Matrix(T)
	{
	this(size_t height, size_t width)
	{
	ptr = (new T[height * width]).ptr;
	this.height = height;
	this.width = width;
	}
	T[] opIndex(size_t heightI) @safe pure
	{
	assert(heightI < height);
	auto l = heightI * width;
	return ptr[l .. l + width];
	}
	auto ref opIndex(size_t heightI, size_t widthI)
	{
	assert(heightI < height);
	assert(widthI < width);
	return ptr[heightI * width + widthI];
	}
	T[] array()
	{
	return ptr[0..width * height];
	}
	T[] front()
	{
	assert(height);
	return ptr[0..width];
	}
	void popFront()
	{
	assert(height);
	ptr+= width;
	height--;
	}
	bool empty()
	{
	return length == 0;
	}
	size_t length() const
	{
	return height;
	}
	T[][] arrays()
	{
	auto r = new T[][height];
	foreach(i, ref e; r)
	{
	const l = i * width;
	e = ptr[l..l+width];
	}
	return r;
	}
	T* ptr;
	size_t height;
	size_t width;
	}
	unittest
	{
	auto m = Matrix!float(3, 4);
	foreach(i; 0..3)
	{
	auto r = m[i];
	foreach(j; 0..4)
	{
	const q = (i+3) * (j-5);
	m[i, j] = q;
	m[i, j]++;
	assert(m[i, j] == q+1);
	assert(r[j] == q+1);
	}
	}
	}

view raw matrix.d hosted with ❤ by GitHub

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

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

	foreach(d; 0..dl)
	{
	auto pt = pdt[d];
	auto nw = ndw[d];
	auto qw = qdw[d];
	foreach(w; 0..wl)
	qwd[w, d] = qw[w] = nw[w] / dotProduct(pt, pwt[w]);
	}
	foreach(t; 0..tl)
	{
	auto pd = ptd[t];
	auto pw = ptw[t];
	auto rw = rtw[t];
	foreach(w; 0..wl)
	rw[w] = pw[w] * dotProduct(pd, qwd[w]);
	foreach(d; 0..dl)
	rdt[d, t] = pdt[d, t] * dotProduct(pw, qdw[d]);
	}
	foreach(t; 0..tl)
	{
	auto rw = rtw[t];
	auto pw = ptw[t];
	const s = rw.sum();
	foreach(w; 0..wl)
	pwt[w, t] = pw[w] = rw[w] / s;
	}
	foreach(d; 0..dl)
	{
	auto rt = rdt[d];
	auto pt = pdt[d];
	const s = rt.sum();
	foreach(t; 0..tl)
	ptd[t, d] = pt[t] = rt[t] / s;
	}

view raw commonVersion.d hosted with ❤ by GitHub

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

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

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

BLAS реализация

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

	sgemm(Order.RowMajor, Transpose.NoTrans, Transpose.NoTrans,
	dl, wl, tl, 1.0, pdt.ptr, tl, ptw.ptr, wl, 0.0, qdw.ptr, wl);
	qdw.array[] /= ndw.array[];
	sgemm(Order.RowMajor, Transpose. Trans, Transpose.NoTrans,
	tl, wl, dl, 1.0, pdt.ptr, tl, qdw.ptr, wl, 0.0, rtw.ptr, wl);
	rtw.array[] *= ptw.array[];
	sgemm(Order.RowMajor, Transpose.NoTrans, Transpose. Trans,
	dl, tl, wl, 1.0, qdw.ptr, wl, ptw.ptr, wl, 0.0, rdt.ptr, tl);
	rdt.array[] *= pdt.array[];
	for(auto a = ptw, b = rtw; !a.empty; a.popFront, b.popFront)
	a.front[] = b.front[] / b.front.sum;
	for(auto a = pdt, b = rdt; !a.empty; a.popFront, b.popFront)
	a.front[] = b.front[] / b.front.sum;

view raw blasVersion.d hosted with ❤ by GitHub

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

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

Результаты

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

Заключение

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