0811 Affinity Propagation

message passing 방식

“Affinity Propagation attempts to find the exemplars that maximize the net similarity, i.e. the overall sum of similarities between all exemplars and their member data points.”^[1]

한문장으로 AP를 잘 설명한것 같아서 퍼왔다. Sparse matrix를 다루는 페이지에도 AP의 간략한 설명이 있으니 참고.

해당 논문^[2]을 보면, 내용 자체는 어렵지 않고, 실제로 사용할 수 있는 MATLAB코드도 Frey Lab에 공개되어 있다. 아래 코드와 수식은 대부분 해당 논문과 Frey Lab이 출처이고 이해를 돕기 위해 약간 수정하거나 간략한 주석을 추가했다. (링크한 것은 이해하기 쉬운 코드이고, Frey Lab에 가보면 더 깔끔하게 동작하는 코드를 얻을 수 있다.)

AP는 다음 세 식에 따라 동작한다.

$$\begin{equation} r(i,k) \leftarrow s(i,k) - \max_{k'\neq k} \{ a(i, k') + s(i, k') \} \end{equation}$$ $$\begin{equation} a(i,k) \leftarrow \min \Big\{0, r(k,k) + \sum_{i' \not\in \{i, k\}} \max\{0, r(i', k)\} \Big\} \end{equation}$$ $$\begin{equation} a(k,k) \leftarrow \sum_{i' \neq k} \max \{0, r(i', k) \} \end{equation}$$

	% test data generation
	snr = 7;
	xy = [awgn(repmat(1, 100, 1), snr) awgn(repmat(2, 100, 1), snr);
	awgn(repmat(3, 100, 1), snr) awgn(repmat(4, 100, 1), snr)];
	S = squareform(pdist(xy)); % distance matrix
	S = -S + max(S(:))+9;
	S(logical(eye(200))) = 0;

	% affinity propagation
	% http://www.psi.toronto.edu/affinitypropagation/FreyDueckScience07.pdf

	N=size(S,1); A=zeros(N,N); R=zeros(N,N); % Initialize messages
	S=S+1e-12randn(N,N)(max(S(:))-min(S(:))); % Remove degeneracies
	lam=0.5; % Set damping factor
	for iter=1:1000000
	%
	% Compute responsibilities (formula (1))
	%
	Rold=R;
	AS=A+S;
	[Y,I]=max(AS,[],2); %row max. I is index.
	for i=1:N
	AS(i,I(i))=-realmax; % remove the largest. ref.realmax=-1.7977e+308
	end;
	[Y2,I2]=max(AS,[],2); %second max
	R=S-repmat(Y,[1,N]); % S - {first max}
	for i=1:N
	R(i,I(i))=S(i,I(i))-Y2(i); % {first max} - {second max}
	end;
	R=(1-lam)R+lamRold; % Dampen responsibilities

	%
	% Compute availabilities (formula (2), (3))
	%
	Aold=A;
	Rp=max(R,0);
	for k=1:N
	Rp(k,k)=R(k,k); % set diagnal to 0. and give R(k,k)-term of formula (2)
	end;
	A=repmat(sum(Rp,1),[N,1])-Rp; %'sigma column except me'
	dA=diag(A);
	A=min(A,0);
	for k=1:N
	A(k,k)=dA(k);
	end;
	A=(1-lam)A+lamAold; % Dampen availabilities
	end;
	fprintf('end. iter = %d, K = %d \n', iter, K);
	E=R+A; % Pseudomarginals
	I=find(diag(E)>0);
	K=length(I); % Indices of exemplars

	fprintf('K = %d\n', K);
	[tmp, c]=max(S(:,I),[],2);
	c(I)=1:K;
	idx=I(c); % Assignments


	scatter(xy(:,1), xy(:,2), [], c);
	title(K);

view raw ap.mat hosted with ❤ by GitHub

	% [idx,netsim,dpsim,expref]=apclusterLeveraged(s,p,frac,sweeps)
	%
	% frac is the fraction of data points to be considered.
	% sweeps is the number of times to run affinity propagation.
	%

	function [idx,netsim,dpsim,expref]=apclusterLeveraged(S,p,frac,sweeps);
	% p(i) indicates the preference that data point i be chosen as an exemplar.
	% Often a good choice is to set all preferences to median(s);

	%
	% dpsim = the sum of the similarities of the data points to their exemplars
	% expref = the sum of the preferences of the identified exemplars
	% netsim = dpsim + expref
	%

	% C = zeros(size(S),'uint8'); % for counting S-computations
	if isa(S,'function_handle'), N=S(0); else N=size(S,1); end; M=ceil(frac*N);
	s=zeros(M*N+N,3);
	s(1:M*N,1)=repmat([1:N]',[M,1]);

	% diagonal part
	s(M*N+(1:N),1)=(1:N)';
	s(M*N+(1:N),2)=(1:N)';
	s(M*N+(1:N),3)=realmin('single'); % 1.1755e-38 : noise


	ii=randperm(N); ii=ii(1:M);
	netsimbest=-1e300; netsimbest=-Inf;
	for sw=1:sweeps
	% make sparse distance matrix.
	% rows = use all(=N) points.
	% cols = only picks M points out of N points
	for m=1:M
	s((m-1)N+1:mN,2)=ii(m);
	s((m-1)N+1:mN,3)=S(1:N,ii(m))+realmin('single');
	% C(1:N,ii(m)) = C(1:N,ii(m))+1; % counting S-computations
	end;

	% run sparse.
	[idx,netsim,dpsim,expref]=apclustermex(sparse(s(:,1),s(:,2),s(:,3),N,N),p);
	fprintf('.');
	if netsim > netsimbest \| isinf(netsimbest), % should allow updates if netsim=-Inf;
	netsimbest=netsim;
	idxbest=idx; dpsimbest=dpsim; exprefbest=expref;
	tmp1=unique(idx)';
	tmp2=setdiff([1:N],tmp1);
	ii=[tmp1,tmp2(randperm(length(tmp2)))]; % sampling and shuffle
	ii=ii(1:M);
	else
	tmp1=unique(idxbest)';
	tmp2=setdiff([1:N],tmp1);
	ii=[tmp1,tmp2(randperm(length(tmp2)))];
	ii=ii(1:M);
	end;
	end;
	idx=idxbest; netsim=netsimbest; dpsim=dpsimbest; expref=exprefbest;

view raw leveraged_ap.mat hosted with ❤ by GitHub

	N=100; x=rand(N,2); % Create N, 2-D data points
	M=N*N-N; s=zeros(M,3); % Make ALL N^2-N similarities
	j=1;
	for i=1:N
	for k=[1:i-1,i+1:N]
	s(j,1)=i; s(j,2)=k; s(j,3)=-sum((x(i,:)-x(k,:)).^2);
	j=j+1;
	end;
	end;
	p=median(s(:,3)); % Set preference to median similarity
	[idx,netsim,dpsim,expref]=apclusterSparse(s,p,'plot');
	fprintf('Number of clusters: %d\n',length(unique(idx)));
	fprintf('Fitness (net similarity): %f\n',netsim);
	figure; % Make a figures showing the data and the clusters
	for i=unique(idx)'
	ii=find(idx==i); h=plot(x(ii,1),x(ii,2),'o'); hold on;
	col=rand(1,3); set(h,'Color',col,'MarkerFaceColor',col);
	xi1=x(i,1)ones(size(ii)); xi2=x(i,2)ones(size(ii));
	line([x(ii,1),xi1]',[x(ii,2),xi2]','Color',col);
	end;
	axis equal tight;

view raw sparse_eg1.mat hosted with ❤ by GitHub

	function [idx,netsim,dpsim,expref]=apclusterSparse(s,p,varargin);

	% Handle arguments to function
	if nargin<2 error('Too few input arguments');
	else
	maxits=1000; convits=100; lam=0.9; plt=0; details=0; nonoise=0;
	i=1;
	while i<=length(varargin)
	if strcmp(varargin{i},'plot')
	plt=1; i=i+1;
	elseif strcmp(varargin{i},'details')
	details=1; i=i+1;
	elseif strcmp(varargin{i},'nonoise')
	nonoise=1; i=i+1;
	elseif strcmp(varargin{i},'maxits')
	maxits=varargin{i+1};
	i=i+2;
	if maxits<=0 error('maxits must be a positive integer'); end;
	elseif strcmp(varargin{i},'convits')
	convits=varargin{i+1};
	i=i+2;
	if convits<=0 error('convits must be a positive integer'); end;
	elseif strcmp(varargin{i},'dampfact')
	lam=varargin{i+1};
	i=i+2;
	if (lam<0.5)\|\|(lam>=1)
	error('dampfact must be >= 0.5 and < 1');
	end;
	else i=i+1;
	end;
	end;
	end;
	if lam>0.9
	fprintf('\n*** Warning: Large damping factor in use. Turn on plotting\n');
	fprintf(' to monitor the net similarity. The algorithm will\n');
	fprintf(' change decisions slowly, so consider using a larger value\n');
	fprintf(' of convits.\n\n');
	end;

	% Check that standard arguments are consistent in size
	if length(size(s))~=2 error('s should be a 2D matrix');
	elseif length(size(p))>2 error('p should be a vector or a scalar');
	elseif size(s,2)==3
	tmp=max(max(s(:,1)),max(s(:,2)));
	if length(p)==1 N=tmp; else N=length(p); end;
	if tmp>N
	error('data point index exceeds number of data points');
	elseif min(min(s(:,1)),min(s(:,2)))<=0
	error('data point indices must be >= 1');
	end;
	elseif size(s,1)==size(s,2)
	N=size(s,1);
	if (length(p)~=N)&&(length(p)~=1)
	error('p should be scalar or a vector of size N');
	end;
	else error('s must have 3 columns or be square'); end;

	% Make vector of preferences
	if length(p)==1
	p=p*ones(N,1);
	end;

	% Append self-similarities (preferences) to s-matrix
	tmps=[repmat([1:N]',[1,2]),p];
	s=[s;tmps];
	M=size(s,1);

	% In case user did not remove degeneracies from the input similarities,
	% avoid degenerate solutions by adding a small amount of noise to the
	% input similarities
	if ~nonoise
	rns=randn('state');
	randn('state',0);
	s(:,3)=s(:,3)+(epss(:,3)+realmin100).*rand(M,1);
	randn('state',rns);
	end;

	% Construct indices of neighbors % refer my comment below.
	ind1e=zeros(N,1);
	for j=1:M
	k=s(j,1);
	ind1e(k)=ind1e(k)+1;
	end;
	ind1e=cumsum(ind1e);
	ind1s=[1;ind1e(1:end-1)+1];
	ind1=zeros(M,1);

	for j=1:M
	k=s(j,1);
	ind1(ind1s(k))=j;
	ind1s(k)=ind1s(k)+1;
	end;
	ind1s=[1;ind1e(1:end-1)+1];
	ind2e=zeros(N,1);

	for j=1:M
	k=s(j,2);
	ind2e(k)=ind2e(k)+1;
	end;
	ind2e=cumsum(ind2e);
	ind2s=[1;ind2e(1:end-1)+1];
	ind2=zeros(M,1);

	for j=1:M
	k=s(j,2);
	ind2(ind2s(k))=j;
	ind2s(k)=ind2s(k)+1;
	end;

	ind2s=[1;ind2e(1:end-1)+1];

	% Allocate space for messages, etc
	A=zeros(M,1);
	R=zeros(M,1);
	t=1;
	if plt
	netsim=zeros(1,maxits+1);
	end;
	if details
	idx=zeros(N,maxits+1);
	netsim=zeros(1,maxits+1);
	dpsim=zeros(1,maxits+1);
	expref=zeros(1,maxits+1);
	end;

	% Execute parallel affinity propagation updates
	e=zeros(N,convits);
	dn=0;
	i=0;
	while ~dn
	i=i+1;

	% Compute responsibilities
	for j=1:N
	ss=s(ind1(ind1s(j):ind1e(j)),3);
	as=A(ind1(ind1s(j):ind1e(j)))+ss;

	[Y,I]=max(as);
	as(I)=-realmax;
	[Y2,I2]=max(as);

	r=ss-Y;
	r(I)=ss(I)-Y2;

	R(ind1(ind1s(j):ind1e(j)))=(1-lam)*r+ ...
	lam*R(ind1(ind1s(j):ind1e(j)));
	end;

	% Compute availabilities
	for j=1:N
	rp=R(ind2(ind2s(j):ind2e(j)));
	rp(1:end-1)=max(rp(1:end-1),0);
	a=sum(rp)-rp;
	a(1:end-1)=min(a(1:end-1),0);
	A(ind2(ind2s(j):ind2e(j)))=(1-lam)*a+ ...
	lam*A(ind2(ind2s(j):ind2e(j)));
	end;

	% Check for convergence
	E=((A(M-N+1:M)+R(M-N+1:M))>0);
	e(:,mod(i-1,convits)+1)=E;
	K=sum(E);
	if i>=convits \|\| i>=maxits
	se=sum(e,2);
	unconverged=(sum((se==convits)+(se==0))~=N);
	if (~unconverged&&(K>0))\|\|(i==maxits)
	dn=1;
	end;
	end;

	% Handle plotting and storage of details, if requested
	if plt\|\|details
	if K==0
	tmpnetsim=nan; tmpdpsim=nan; tmpexpref=nan; tmpidx=nan;
	else
	tmpidx=zeros(N,1); tmpdpsim=0;
	tmpidx(find(E))=find(E); tmpexpref=sum(p(find(E)));
	discon=0;
	for j=find(E==0)'
	ss=s(ind1(ind1s(j):ind1e(j)),3);
	ii=s(ind1(ind1s(j):ind1e(j)),2);
	ee=find(E(ii));
	if length(ee)==0 discon=1;
	else
	[smx imx]=max(ss(ee));
	tmpidx(j)=ii(ee(imx));
	tmpdpsim=tmpdpsim+smx;
	end;
	end;
	if discon
	tmpnetsim=nan; tmpdpsim=nan; tmpexpref=nan; tmpidx=nan;
	else tmpnetsim=tmpdpsim+tmpexpref;
	end;
	end;
	end;
	if details
	netsim(i)=tmpnetsim; dpsim(i)=tmpdpsim; expref(i)=tmpexpref;
	idx(:,i)=tmpidx;
	end;
	if plt
	netsim(i)=tmpnetsim;
	figure(234);
	tmp=1:i; tmpi=find(~isnan(netsim(1:i)));
	plot(tmp(tmpi),netsim(tmpi),'r-');
	xlabel('# Iterations');
	ylabel('Net similarity of quantized intermediate solution');
	drawnow;
	end;
	end;
	% Identify exemplars
	E=((A(M-N+1:M)+R(M-N+1:M))>0); K=sum(E);
	if K>0
	tmpidx=zeros(N,1); tmpidx(find(E))=find(E); % Identify clusters
	for j=find(E==0)'
	ss=s(ind1(ind1s(j):ind1e(j)),3);
	ii=s(ind1(ind1s(j):ind1e(j)),2);
	ee=find(E(ii));
	[smx imx]=max(ss(ee));
	tmpidx(j)=ii(ee(imx));
	end;
	EE=zeros(N,1);
	for j=find(E)'
	jj=find(tmpidx==j); mx=-Inf;
	ns=zeros(N,1); msk=zeros(N,1);
	for m=jj'
	mm=s(ind1(ind1s(m):ind1e(m)),2);
	msk(mm)=msk(mm)+1;
	ns(mm)=ns(mm)+s(ind1(ind1s(m):ind1e(m)),3);
	end;
	ii=jj(find(msk(jj)==length(jj)));
	[smx imx]=max(ns(ii));
	EE(ii(imx))=1;
	end;
	E=EE;
	tmpidx=zeros(N,1); tmpdpsim=0;
	tmpidx(find(E))=find(E); tmpexpref=sum(p(find(E)));
	for j=find(E==0)'
	ss=s(ind1(ind1s(j):ind1e(j)),3);
	ii=s(ind1(ind1s(j):ind1e(j)),2);
	ee=find(E(ii));
	[smx imx]=max(ss(ee));
	tmpidx(j)=ii(ee(imx));
	tmpdpsim=tmpdpsim+smx;
	end;
	tmpnetsim=tmpdpsim+tmpexpref;
	else
	tmpidx=nan*ones(N,1); tmpnetsim=nan; tmpexpref=nan;
	end;
	if details
	netsim(i+1)=tmpnetsim; netsim=netsim(1:i+1);
	dpsim(i+1)=tmpnetsim-tmpexpref; dpsim=dpsim(1:i+1);
	expref(i+1)=tmpexpref; expref=expref(1:i+1);
	idx(:,i+1)=tmpidx; idx=idx(:,1:i+1);
	else
	netsim=tmpnetsim; dpsim=tmpnetsim-tmpexpref;
	expref=tmpexpref; idx=tmpidx;
	end;
	if plt\|\|details
	fprintf('\nNumber of identified clusters: %d\n',K);
	fprintf('Fitness (net similarity): %f\n',tmpnetsim);
	fprintf(' Similarities of data points to exemplars: %f\n',dpsim(end));
	fprintf(' Preferences of selected exemplars: %f\n',tmpexpref);
	fprintf('Number of iterations: %d\n\n',i);
	end;
	if unconverged
	fprintf('\n*** Warning: Algorithm did not converge. The similarities\n');
	fprintf(' may contain degeneracies - add noise to the similarities\n');
	fprintf(' to remove degeneracies. To monitor the net similarity,\n');
	fprintf(' activate plotting. Also, consider increasing maxits and\n');
	fprintf(' if necessary dampfact.\n\n');
	end;

view raw sparse_eg2.mat hosted with ❤ by GitHub

	snr = 7;
	xy = [awgn(repmat(1, 5, 1), snr) awgn(repmat(2, 5, 1), snr);
	awgn(repmat(3, 5, 1), snr) awgn(repmat(4, 5, 1), snr)];
	S = zeros(10*9/2, 3);
	C = combnk(1:10, 2);

	%>> xy
	%
	%xy =
	%
	% 0.4872 2.0837
	% 1.0468 1.9632
	% 1.3226 1.1366
	% 2.1549 1.8039
	% 0.7021 1.1983
	% 3.3754 3.7318
	% 2.6033 4.2189
	% 3.0447 4.3303
	% 2.7568 4.7647
	% 3.1356 3.9133

	for i=1:length(C)
	xi = C(i,1);
	yi = C(i,2);
	d = norm(xy(xi, :)-xy(yi, :)); %this is falsy. d should be -norm(...)
	S(i, :) = [C(i,1), C(i,2), d];
	end

	%>> S
	%
	%S =
	%
	% 1 2 0.57245
	% 1 3 1.2629
	% 1 4 1.691
	% 1 5 0.91104
	% 1 6 3.3253
	% 1 7 3.0061
	% 1 8 3.4041
	% 1 9 3.5126
	% 1 10 3.2189
	% 2 3 0.87139
	% 2 4 1.1194
	% 2 5 0.83891
	% 2 6 2.9241
	% 2 7 2.7406
	% 2 8 3.0975
	% 2 9 3.2821
	% 2 10 2.8576
	% 3 4 1.0668
	% 3 5 0.62358
	% 3 6 3.309
	% 3 7 3.3378
	% 3 8 3.6284
	% 3 9 3.9013
	% 3 10 3.3162
	% 4 5 1.5739
	% 4 6 2.2818
	% 4 7 2.4562
	% 4 8 2.6785
	% 4 9 3.0213
	% 4 10 2.3262
	% 5 6 3.6831
	% 5 7 3.5691
	% 5 8 3.9111
	% 5 9 4.1159
	% 5 10 3.6459
	% 6 7 0.91282
	% 6 8 0.6837
	% 6 9 1.2039
	% 6 10 0.30071
	% 7 8 0.45522
	% 7 9 0.56697
	% 7 10 0.61373
	% 8 9 0.52117
	% 8 10 0.42676
	% 9 10 0.93185

	apclusterSparse(S, 0);

view raw sparse_eg3 hosted with ❤ by GitHub

	[(1:55)' s ind1 ind2]
	ans =
	1 1 2 1.2908 1 46
	2 1 3 0.98971 2 1
	3 1 4 1.3646 3 47
	4 1 5 0.53551 4 2
	5 1 6 3.0272 5 10
	6 1 7 3.1811 6 48
	7 1 8 2.6518 7 3
	8 1 9 2.6847 8 11
	9 1 10 2.9787 9 18
	10 2 3 1.4746 46 49
	11 2 4 0.19874 10 4
	12 2 5 1.2703 11 12
	13 2 6 3.1465 12 19
	14 2 7 3.4175 13 25
	15 2 8 2.9135 14 50
	16 2 9 2.7438 15 5
	17 2 10 3.1162 16 13
	18 3 4 1.4072 17 20
	19 3 5 0.45646 47 26
	20 3 6 2.0407 18 31
	21 3 7 2.2126 19 51
	22 3 8 1.6809 20 6
	23 3 9 1.6951 21 14
	24 3 10 1.9934 22 21
	25 4 5 1.2638 23 27
	26 4 6 2.9858 24 32
	27 4 7 3.2659 48 36
	28 4 8 2.7697 25 52
	29 4 9 2.5816 26 7
	30 4 10 2.9576 27 15
	31 5 6 2.4971 28 22
	32 5 7 2.6632 29 28
	33 5 8 2.132 30 33
	34 5 9 2.1501 49 37
	35 5 10 2.4495 31 40
	36 6 7 0.36396 32 53
	37 6 8 0.45161 33 8
	38 6 9 0.40532 34 16
	39 6 10 0.059503 35 23
	40 7 8 0.53177 50 29
	41 7 9 0.72777 36 34
	42 7 10 0.35529 37 38
	43 8 9 0.44005 38 41
	44 8 10 0.39225 39 43
	45 9 10 0.38415 51 54
	46 1 1 3.4908e-307 40 9
	47 2 2 7.0924e-307 41 17
	48 3 3 1.0571e-307 42 24
	49 4 4 5.8453e-307 52 30
	50 5 5 1.4614e-306 43 35
	51 6 6 9.2654e-307 44 39
	52 7 7 5.9746e-307 53 42
	53 8 8 1.1337e-306 45 44
	54 9 9 6.0503e-307 54 45
	55 10 10 1.7434e-306 55 55

view raw sparse_eg4 hosted with ❤ by GitHub

(similiarity를 실수로 distance를 넣었을 때 결과가 상당히 특이했다. 좀 더 실험해봐야 할것 같다.)

AP는 $N^2$ size의 availability, responsibility matrix를 모두 메모리에 들고 있어야 해서, 아주 큰 데이터에 대해서는 돌려볼 수가 없다($N=$1M이고 similiarity하나당 4byte라고 했을 때, 1M×1M×4byte$\approx$232TB가 필요하다). 메모리문제를 약간 완화시키는 방법으로 <a href='leveraged_ap.html'>leveraged ap</a>가 있는데, 완화효과가 그리 크지 않고, cluster center의 개수에 영향을 줄 수 있어서 별로인것 같다. 예를 들어, frac을 작게 잡아서 column을 10개밖에 쓰지 못하게 되면 center의 개수도 10개 이하가 된다. row가 수백만인 경우도 흔하기 때문에 결국 좋은 방법이라고 하기는 힘들다. 이렇게 했을 때 결국 AP에 수렴하는지 보인 문서도 없는 듯 보인다.(대강 봐서는 수렴할것 같기는 하다)

Similiarity matirx가 sparse할 때도 사용할 수 있다. <a href='apcluster_sparse.html'>별도로 정리한 것</a>이 있다.

Fujiwara¹⁾에 위 식들의 간략한 버전이 나온다. 이해하기 더 쉬운것 같다. 표기만 약간 바꾸어 다시 적으면 다음과 같다.

$$ \begin{align} &r(i,j) = s(i,j) - \max_{k\neq j}\{a(i,k) + s(i,k)\} & (i \neq j) \\ &r(i,i) = s(i,i) - \max_{k\neq i}\{s(i,k)\} & \\ &a(i,j) = \min\left\{ 0, r(j,j) + \sum_{k\neq i,j} \max \{ 0, r(k,j) \} \right\} & (i \neq j ) \\ &a(i,i) = \sum_{k\neq i} \max\{ 0, r(k,i) \} \end{align} $$

좌변의 $a$, $r$과 우변의 $a$, $r$은 의미가 달라서(damping을 고려하지 않는다면 그냥 iteration전/후차이) Fujiwara에는 좌변이 $\rho$, $\alpha$로 표기되어 있다.

‘responsibility’와 ‘availability’에 관한 상세한 설명이 Frey에 나오지만, Fujiwara에 나온 간략한 설명도 좋은것 같다. $r(i,j)$는 $i$가 $j$에게 보내는 메시지이고, $a(i,j)$는 $j$가 $i$에게 보내는 메시지인데 의미는 모두 ‘$j$가 exemplar로 얼마나 적당하냐’ 하는 것이다($ \definecolor{blue}{RGB}{0,0,200} \color{blue} r(i\rightarrow j), a(i\leftarrow j) $로 적어도 될듯. s, r로 해서 send, receive가 더 나으려나 생각도 해봤는데 send, receive는 누구 기준해서 주고받는거냐가 헷갈려서 별로 도움 안될것 같다. $\overset{i\rightarrow j}M, \overset{i\leftarrow j}M$도 괜찮을것 같다). 그래서 나중에 exemplar를 고를 때도 둘의 합이 최대가 되는 곳($\arg\max \{ r(i,j)+a(i,j):j=1,2,\dots,N\} $)이 선택된다.

식(5)는 Frey 읽을때도 조금 이해가 안갔었다. Frey에는 식(1)아래에 다음과 같이 실려있다.

For $k = i$, the responsibility $r(k,k)$ is set to the input preference that point $k$ be chosen as an exemplar, $s(k,k)$, minus the largest of the similarities between point $i$ and all other candidate exemplars. This “self-responsibility” reflects accumulated evidence that point $k$ is an exemplar, based on its input preference tempered by how ill-suited it is to be assigned to another exemplar.

식(5)는 이 말을 식으로 옮겨둔 것으로 보이는데, 정작 Frey Lab의 코드(<a href='ap_code.html'>위의 링크</a>에 있는 코드)를 보면 이 부분이 없다. 즉 $r(i,k)$는 $i$와 $k$가 같든 다르든 식(1)에 따라 갱신된다. 만약 $r$이 매번 저렇게 갱신되면 $\lambda$가 0일 때 $r(k,k)$가 변화가 없어야 하지만 실제로는 그렇지 않다. 따라서 (5)는 $a(i,j)=0$인 초기조건에서 식(4)의 $k=i$일 때인것 같다. Fray의 논문에 나온 부분도 초기조건 설명이 아닌가 한다. (억지로 이렇게 끼워맞춰 생각해보지만 속시원한 느낌이 없기는 매한가지)

(작성중)

spiral에 대한 성능은?

↑ Fujiwara, Yasuhiro, Go Irie, and Tomoe Kitahara. “Fast algorithm for affinity propagation.” IJCAI Proceedings-International Joint Conference on Artificial Intelligence. Vol. 22. No. 3. 2011.
↑ Frey, Brendan J., and Delbert Dueck. “Clustering by passing messages between data points.” science 315.5814 (2007): 972-976.