Можно ли оптимизировать этот код Matlab для выполнения векторного квантования с помощью центроидов из k-средних?

Я создал кодовую книгу с использованием k-средних размеров 4000x300 (4000 центроидов, каждый с 300 функциями). Затем, используя кодовую книгу, я хочу пометить входной вектор (для целей сортировки позже). Входной вектор имеет размер Nx300, где N - общее количество входных экземпляров, которые я получаю.

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

Мой текущий код Matlab выглядит так:

function labels = assign_labels(centroids, X)
labels = zeros(size(X, 1), 1);

% for each X, calculate the distance from each centroid
for i = 1:size(X, 1)
    % distance of X_i from all j centroids is: sum((X_i - centroid_j)^2)
    % note: we leave off the sqrt as an optimization
    distances = sum(bsxfun(@minus, centroids, X(i, :)) .^ 2, 2);
    [value, label] = min(distances);
    labels(i) = label;
end     

Однако этот код все еще довольно медленный (для моих целей), и я надеялся, что есть способ его дальнейшей оптимизации.

Одна очевидная проблема заключается в том, что существует цикл for, который является проклятием для хорошей производительности в Matlab. Я пытался придумать способ избавиться от этого, но безуспешно (я изучал использование arrayfun в сочетании с bsxfun, но не получил, чтобы это работало). В качестве альтернативы, если кто-то знает какой-либо другой способ ускорить это, я был бы очень признателен.

Обновить

После некоторого поиска я не смог найти отличного решения с использованием Matlab, поэтому решил посмотреть, что используется в пакете Python scikits.learn для 'euclidean_distance' (сокращенно):

 XX = sum(X * X, axis=1)[:, newaxis]
 YY = Y.copy()
 YY **= 2
 YY = sum(YY, axis=1)[newaxis, :]
 distances = XX + YY
 distances -= 2 * dot(X, Y.T)
 distances = maximum(distances, 0)

который использует биномиальную форму евклидова расстояния ((x-y) ^ 2 -> x ^ 2 + y ^ 2 - 2xy), которое из того, что я читал, обычно выполняется быстрее. Мой полностью непроверенный перевод Matlab:

 XX = sum(data .* data, 2);
 YY = sum(center .^ 2, 2);
 [val, ~] = max(XX + YY - 2*data*center');

person Abe Schneider    schedule 20.04.2011    source источник
comment
связанные: эквивалент pdist2 в MATLAB версии 7   -  person Amro    schedule 11.07.2012


Ответы (4)


Вы можете векторизовать его, преобразовав в ячейки и используя cellfun:

[nRows,nCols]=size(X);
XCell=num2cell(X,2);
dist=reshape(cell2mat(cellfun(@(x)(sum(bsxfun(@minus,centroids,x).^2,2)),XCell,'UniformOutput',false)),nRows,nRows);
[~,labels]=min(dist);

Объяснение:

  • Мы назначаем каждую строку X отдельной ячейке во второй строке
  • Эта часть @(x)(sum(bsxfun(@minus,centroids,x).^2,2)) является анонимной функцией, которая аналогична вашей строке distances=..., и, используя cell2mat, мы применяем ее к каждой строке X.
  • В этом случае метки являются индексами минимальной строки каждого столбца.
person abcd    schedule 20.04.2011
comment
Похоже, что num2cell (X) превращает каждый элемент X в свой собственный вектор. Однако для (x-центроидов) мы хотим вычесть каждую строку X против каждого центроида. Так следует ли вместо этого читать: XCell = num2cell (X, 2)? - person Abe Schneider; 21.04.2011
comment
@ Абэ, ты права. Хотя он, безусловно, возвращает тот же ответ, при его применении к каждому элементу возникают ненужные накладные расходы на вызов функции. Я это исправил. Но помните, что bsxfun и cellfun обычно являются просто краткими способами написания цикла, и не обязательно должны быть быстрее (иногда это так, но не всегда). Если рассчитать время вашего цикла и кода cellfun для матрицы того же размера, что и у вас, они были довольно даже на отметке 93 секунды, отличаясь только десятым местом. - person abcd; 21.04.2011
comment
Большое спасибо! Да, я надеялся, что внутренний цикл Matlab может быть лучше, чем их общие циклы for, но я также не вижу больших улучшений. Также цикл for позволяет использовать parfor. Первоначально хотел использовать графический процессор, но bsxfun не позволяет использовать анонимные функции и, следовательно, не может передавать «центроиды». Одно улучшение, которое мне удалось найти, - это чья-то публикация на SO, которая зацикливалась на центроидах вместо данных. У меня 4000 центроидов и от 1e5 до 1e6 векторов признаков в моих данных. Так что, перебирая центроиды, я могу ускорить процесс с помощью матричной математики. - person Abe Schneider; 21.04.2011
comment
В моем предыдущем комментарии не было места, но я использовал код назначения метки отсюда: stackoverflow.com/questions/1373516/matlabk-means-clustering/ - person Abe Schneider; 21.04.2011

Используйте следующую функцию для расчета расстояния. Вы должны увидеть ускорение на порядок

Две матрицы A и B имеют столбцы как размеры и строки как каждую точку. A - это ваша матрица центроидов. B - ваша матрица точек данных.

function D=getSim(A,B)
    Qa=repmat(dot(A,A,2),1,size(B,1));
    Qb=repmat(dot(B,B,2),1,size(A,1));
    D=Qa+Qb'-2*A*B';
person twerdster    schedule 24.04.2011
comment
+1 вам. Это действительно быстрее. Если подумать, я не знаю, почему я этого не сделал, вместо того, чтобы использовать cellfun. В последнее время я слишком много использую cellfun, и мне следует убрать его из своего набора инструментов, по крайней мере, временно :) - person abcd; 24.04.2011

Для истинной матричной реализации вы можете попробовать что-то вроде:

  P2 = kron(centroids, ones(size(X,1),1));
  Q2 = kron(ones(size(centroids,1),1), X);

  distances = reshape(sum((Q2-P2).^2,2), size(X,1), size(centroids,1));

Примечание. Предполагается, что данные организованы как [x1 y1 ...; x2 y2 ...; ...]

person jmetz    schedule 22.04.2011
comment
Это использует много памяти - как и любая действительно векторизованная версия. (Более медленная) альтернатива уже предоставлена ​​@yoda. Также вы можете заменить kron на repmat - это была просто версия, которую я уже использовал в одном из моих собственных проектов. - person jmetz; 22.04.2011
comment
Хм. Я собираюсь попробовать это, так как мой текущий подход займет ›4 дня. Тем не менее, я не знаком с kron, поэтому я попытаюсь переписать его, используя описанные выше повторные маты. Если у вас есть шанс, не могли бы вы проверить? - person Abe Schneider; 22.04.2011
comment
Ничего. На бумаге это выглядит так, как будто я делаю обе матрицы 3-х мерными (повторяю матрицу данных для размера (центроиды, 1) раз и повторяю каждый размер центроида (данные, 1) раз для каждого центроида, я могу просто вычесть, квадрат, sum, а затем минимизировать их.Однако я не могу придумать хороший способ убедиться, что вторая матрица. - person Abe Schneider; 22.04.2011
comment
Конечно - не могли бы вы дать мне представление о том, насколько большим будет N? Как я уже упоминал, использование памяти может стать большим довольно быстро - я не думаю, что он будет обрабатывать N намного больше, чем несколько сотен или около того максимума, учитывая 300-дневное пространство измерения. - person jmetz; 22.04.2011
comment
Хорошо, полезно знать. Я сейчас читаю о кронах, но раньше не видел, чтобы они использовались таким образом. Что касается N, он довольно большой, так что это может быть проблемой. Мои данные обычно составляют 100000, а мои центроиды около 4000. Я могу изменить размер данных, так как я их загружаю, но если он слишком мал, я боюсь, что это повлечет за собой большие затраты. - person Abe Schneider; 23.04.2011
comment
Это излишне сохраняет данные, что становится слишком дорогостоящим для больших массивов. Каждая точка в матрице расстояний - это просто скалярное произведение, которое можно записать как | x | ^ 2 -2 * x '* y + | y ​​| ^ 2, что, в свою очередь, довольно легко превратить в матричное вычисление. - person twerdster; 25.04.2011
comment
@mutzmatron - Ему вообще не нужно перебирать наборы данных. Проверьте мой ответ на код. - person twerdster; 26.04.2011
comment
@mutzmatron Ваш метод хранит данные без надобности. Для центроидов размером 4000x300 и X размером 10000x300 мой метод использует около 650 МБ и несколько секунд. Вашему методу немедленно не хватает памяти. - person twerdster; 27.04.2011
comment
@twerdster - искренние извинения - вы правы, ваш код не хранит промежуточные матрицы перед суммами, что позволяет избежать накладных расходов, о которых я упоминал. +1! - person jmetz; 28.04.2011

Вы можете использовать более эффективный алгоритм поиска ближайшего соседа, чем грубая сила. Самый популярный подход - Kd-Tree. Среднее время запроса O (log (n)) вместо сложности грубой силы O (n). Что касается реализации Kd-Trees на Maltab, вы можете посмотреть здесь

person Quant    schedule 07.11.2013