Альтернатива использованию квадратной формы (Matlab)

Question

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

% calculate distances between all points
distances = pdist(m);

Но поскольку pdist возвращает одномерный массив расстояний, нет простого способа выяснить, какая точка имеет наименьшее среднее расстояние (напрямую).Вот почему я использую squareform, а затем вычисляю наименьшее среднее расстояние, например так:

% convert found distances to matrix of distances
distanceMatrix = squareform(distances);

% find index of point with smallest average distance
[~,j] = min(mean(distanceMatrix,2));

Расстояния усредняются для каждого столбца, а переменная j является индексом для столбца (и точка) с наименьшим средним расстоянием.

Это работает, но квадратная форма занимает много времени (этот фрагмент кода повторяется тысячи раз), поэтому я ищу способ оптимизировать его. Кто-нибудь знает более быстрый способ вывести точку с наименьшим средним расстоянием от результатов pdist?

Answer 1

Я думаю, что для вашей задачи использование SQUAREFORM - лучший способ с точки зрения векторизации. Если вы посмотрите на содержание этой функции по

edit squareform

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

[r, c] = size(m);
distanceMatrix = zeros(r);
distanceMatrix(tril(true(r),-1)) = distances;
distanceMatrix = distanceMatrix + distanceMatrix';

Затем запустите тот же код, что и вы, чтобы найти medioid.

Answer 2

Когда pdist вычисляет расстояния между парами наблюдений (1,2, ..., n), расстояния располагаются в следующем порядке:

(2,1), (3,1), ..., (м, 1), (3,2), ..., (м, 2), ..., (м, м – 1))

Чтобы продемонстрировать это, попробуйте следующее:

> X = [.2 .1 .7 .5]';
> D = pdist(X)
.1  .5  .3   .6  .4  .2

В этом примере X хранит n = 4 наблюдения.Результат, D, представляет собой вектор расстояний между наблюдениями (2,1), (3,1), (4,1), (3,2), (4,2), (5,4).Это расположение соответствует элементам нижней треугольной части следующей n-n-матрицы:

M =
0 0 0 0
.1 0 0 0
.5.6 0 0
.3 .4 .2 0

Обратите внимание, что D ( 1 ) = M ( 2,1 ), D (2 ) = ( 3,1 ) и так далее.Таким образом, один из способов получить пару индексов в M, которые соответствуют D (k), состоит в том, чтобы вычислить линейный индекс для D (k) в M. Это можно сделать следующим образом:

% matrix size
n = 4;
% r(j) is the no. of elements in cols 1..j, belonging to the upper triangular part 
r = cumsum(1:n-1);       
% p(j) is the no. elements in cols 1..j, belonging to the lower triangular part 
p = cumsum(n-1:-1:1);
% The linear index of value D(k)
q = find(p >= k, 1);
% The subscript indices of value D(k)
[i j] = ind2sub([n n], k + r(q));

Обратите внимание, что n, r и p нужно устанавливать только один раз.С этого момента вы можете найти индекс для любого заданного k, используя последние две строки.Давайте проверим это:

for k = 1:6
   q = find(p >= k, 1);
   [i, j] = ind2sub([n n], k + r(q));
   fprintf('D(%d) is the distance between observations (%d %d)\n', k, i, j);
end

Вот вывод:
D (1) - расстояние между наблюдениями (2 1)
D (2) - расстояние между наблюдениями (3 1)
D (3) - расстояние между наблюдениями (4 1)
D (4) - расстояние между наблюдениями (3 2)
D (5) - расстояние между наблюдениями (4 2)
D (6) - расстояние между наблюдениями (4 3)

Answer 3

Вот реализация, которая не требует вызова метода квадратной формы:

N1 = 10;
dim = 5;

% generate points
X = randn(N1, dim);

% find mean distance
for iter=N1:-1:1
    d_mean(iter) = mean(pdist2(X(iter,:),X([1:(iter-1) (iter+1):end],:),'euclidean'));
    % D(iter,:) = pdist2(X(iter,:),X([1:(iter-1) (iter+1):end],:),'euclidean');
end

[val ind] = min(d_mean);

Но, не зная больше о вашей проблеме, я понятия не имею, будет ли она быстрее.

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

Удачи.

Answer 4

Нет необходимости использовать squareform:

distances = pdist(m);
l=length(distances);
n=(1+sqrt(1+4*l))/2;
m=[];
for i=1:n
  idx=[1+i:n:length(distances)];
  m(i)=mean(distances(idx));
end

j=min(m);

Я не уверен, но, возможно, это можно также векторизовать, но сейчас уже поздно.