如何为巨大的亲和力矩阵组合谱聚类的拆分运行
How to combine split runs of spectral clustering for a huge affinity matrix
引出问题
我有一个二维复值图像，其中包含一系列短值。我想聚类相似的像素/分割图像。有一个或多或少的静态图像和一个叠加图像，其中有一些斑点，这些斑点在短系列上具有变化的值（主要是复数的角度）。它们在图像的规范中也略显可辨。
我的第一次尝试是 k-means，但它确实根据均值进行了聚类（平均值存在差异，尤其是与周围像素相比，但时间和 angular 信息更大）。我的第二次尝试是 ICA，然后查看幅度最大的 k 个分量，它确实成功地将图像中的某些区域识别为不同，但没有识别出我感兴趣的像素组（视觉上不难识别）认识他们，但他们很小）。
现状
所以因为我的前两次尝试没有成功，所以我用 google 环顾四周，似乎 spectral clustering 可能是合适的。但是我在使用该方法时遇到了一些严重的问题，主要与可用内存有限有关。然后我想，因为我有这么多像素，我可以将光谱聚类应用于单独的数据块。
有人 here suggests clustering slabs first and then combine them, he then says 'at the end you will have the problem of recombining them and this problem can be solved easily'. The bits designated as 'easy' in explanations are hardly ever easy of course. He links to this 论文，但是那个方法并没有处理 slabs 中的所有数据。它排除了不接近主成分的向量。
问题
我的问题分为两部分：
1。如何合并单独细分的结果？特征向量不同，簇数不同。结果看起来它在单独的板中工作。
2。不考虑单独平板中像素之间的距离/亲和力。我可以制作'slabs between slabs'吗？对于那些 L 和 A 不对称的板，不知道如何执行该方法。也许我可以以某种方式比较/合并最后的所有特征向量？
(3.有没有类似或者更好的方法，不需要那么多内存，计算时间也能接受，容易爆)
Matlab代码示例
%% generate data
% get some outer region without data
tempdisk = strel('disk',922/2); tempdisk = double(repmat((1+sqrt(-1)).*tempdisk.Neighborhood,[1 1 15]));
% make some noise
tempnoise = (rand(921,921,15)+sqrt(-1).*rand(921,921,15))./10;
% 'background signal'
tempim1 = double(imresize(mean(imread('cameraman.tif'),3),[921,921])); tempim1 = repmat(tempim1./max(tempim1(:)),[1 1 15]);
% 'target signal'
tempim2 = double(rgb2hsv(imread('fabric.png'))); tempim2 = imresize(tempim2(:,:,2),[921,921]); tempim2 = repmat(tempim2./max(tempim2(:)),[1 1 15]);
sin1 = repmat(permute(sin(2.*pi.*(0:14)./15),[1 3 2]),[921 921 1]);
% combine into data
complexdata = (sin1.*(1.0.*tempim1+0.5.*tempim2.^2).*exp(-sqrt(-1).*2.*pi.*sin1.*(tempim2.^2)).*tempdisk+tempnoise)./1.5;
%this is what the mean data looks like
meannorm = mean(abs(complexdata),3);
meanangle = mean(angle(complexdata),3);
figure; subplot(1,2,1); imshow(meannorm,[]); title('mean norm'); subplot(1,2,2); imshow(meanangle,[]); title('mean angle')
这是生成的数据的样子：
右图中明亮的斑点就是我想要的。它们随时间的变化也最强（并且在时间上相关）。
然后设置集群：
%% perform spectral clustering in seperate slabs
% method from http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf
%get all pixel vectors in a single matrix
complexrows = reshape(permute(complexdata, [3,1,2]), [15, 921*921]);
%k means and eigs dont accept complex, so convert to real here?
complexrowsTranspose = [real(complexrows);imag(complexrows)]';
%lets say 10000 by 10000 matrices are still ok
npix = 10000;
nslabpix = floor(length(complexrowsTranspose)/npix);
nrestpix = rem(length(complexrowsTranspose), npix);
在适合内存的平板中执行光谱聚类：
% spectral clustering
keig = 50;%how many eigenvectors needed? more is better
affinity_sigma = 1;% i dont understand how to calculate this from the paper
tic
% start with last slab (dynamically preallocate)
for islabpix = (nslabpix+1):-1:1;
%print progress
islabpix/nslabpix
toc
if islabpix>nslabpix
pixrange = (1:nrestpix) + ((islabpix-1)*npix);;
else
pixrange = (1:npix) + ((islabpix-1)*npix);
end
%calculate affinity between all voxels in slab
Aff = exp( -squareform(pdist(complexrowsTranspose(pixrange,:))).^2/(2*affinity_sigma^2) ); % affinity matrix
%calculate degree matrix for normalization
Dsq = sparse(size(Aff,1),size(Aff,2)); %degree matrix
for idiag=1:size(Aff,1)
Dsq(idiag,idiag) = sum(Aff(idiag,:))^(1/2);
end
%normalize affinity matrix
Lap = Dsq * Aff * Dsq; %normalize affinity matrix
%calculate eigenvectors of affinity matrix
[eigVectors(pixrange,1:keig), eigValues] = eigs(Lap, keig); %eigenvectors of normalized aff mat
normEigVectors(pixrange,1:keig) = eigVectors(pixrange,1:keig)./repmat(sqrt(sum(abs(eigVectors(pixrange,1:keig)).^2,2)), [1 keig]); %normalize rows of eigen vectors, normr only works on real numbers
% perform k means clustering on weights for eigenvectors
[idx,C,sumd,D] = kmeans([real(normEigVectors(pixrange,1:keig)),imag(normEigVectors(pixrange,1:keig))], 5); %k means on normalized eigenvecotrs
idxval(pixrange) = idx;
end
%reshape into image
idxim = reshape(idxval, [921, 921]);
figure; imshow(idxim,[])
toc
生成的聚类：
结果看起来该方法在每个平板中都在某种程度上起作用；目标是将所有具有稍高范数和更强角度变化的斑点（来自 tempim2 的高饱和度斑点）聚集在一起，这在结果中似乎是可识别的。现在问题主要是单独的平板，并且没有交叉平板集群。这花了我的电脑大约 15 分钟。我减少了本示例的特征值数量和图像大小，因此它可以在可接受的时间内运行。我认为这说明了我的部分问题。
我真的没有答案，但我认为这些指示应该可以帮助您找到答案：
您声称有记忆问题。你确定你的亲和力矩阵是稀疏的吗？您的代码中似乎只有对角线度矩阵是稀疏的。通常当 运行 在 pixels/voxels 上进行谱聚类时，将亲和矩阵设计得非常稀疏（8 连接或 26 连接）。
您将集群描述为 "they are small"。光谱聚类有 known issues 个不同尺度的聚类。你确定你得到了满意的结果吗？
如何计算相邻体素之间的亲和力（相似度）？你也可以测量 差异性 吗？也就是说，你能说一些体素应该不属于同一个簇吗？如果是这样，您是否考虑过使用 correlation clustering？该方法对不同的簇尺度鲁棒性更强，可以自动检测簇数
您是否考虑过使用多尺度/multigrid 方法来粗化您的数据，而不是将其粗暴地分割成 "slabs"？
你看过spectralNet了吗？如果我没记错的话，这个方法应该可以让你"learn"对部分点进行谱聚类，然后使用网络"extrapolate"对新点进行聚类
更新：
鉴于 ，我想说的是，当涉及到非常大的数据的光谱聚类时，将数据粗暴地切成 "slabs" 然后尝试 "stitch" 结果可能不是最好的粗略操作（不是我认为这是不可能的）。解决该问题的更好方法是显着稀疏亲和矩阵：仅计算每个点与其相邻点的成对亲和度，从而得到大部分稀疏的亲和矩阵。这样就可以一次处理所有的点，而不需要 "slice" 和 "stitch"。
关于谱聚类和相关聚类的区别：
为什么谱聚类能够聚类所有点，即使输入亲和力矩阵如此稀疏？它如何分辨出那个点a和一个很远的点 c 应该属于同一个集群，即使它们之间没有计算亲和力？
简单的答案是亲和力的 transitivity：如果 a 类似于 b 并且 b 类似于 c 那么 a 和 c 应该聚集在一起。
问题在哪里？ 在谱聚类中，亲和矩阵中的所有条目都是非负的，这意味着除非绝对没有连接 a 和 c 的路径（机会很小）有一些 "transitive affinity" 建议 a 和 c 应该属于同一个集群。因此，如果查看谱聚类的数学，您会注意到 "trivial solution"，即将所有点放在同一个聚类中，为问题提供了全局最优解。必须人为地强制解决方案具有 k 个集群以避免琐碎的解决方案。
可以做什么？ 如果只考虑正相关性，值 0 是不明确的：它表示 "I didn't bother to compute the affinities between these points"，但也可以表示 "I think these two points should not be in the same cluster"。为了克服这种歧义，我们可以这样引入 negative affinities 如果 A(i, j) > 0 意味着点 i 和点 j 应该在同一个集群中并且确定 A(i, j)，而如果 A(i, j) < 0 意味着 i 和 j 应该 而不是 在同一个集群中（肯定 |A(i, j)|） .引入负亲和力打破了可能 link 远点的 "transitivity chains"，不，将所有点放在同一个集群中不再是微不足道的。
如何利用负亲和力？ 当你的亲和力矩阵同时具有正（吸引力）和负（排斥）值时，你可以使用相关聚类来聚类这些点，它基本上试图最大化affinities/attraction 在每个簇内的点之间，同时最大化不同簇中的点之间的排斥力。一个很好的 属性 相关聚类是它 "automatically" 发现聚类的基础数量，见 sec。 this paper 中的 2 个。
在 I switched to only calculating the affinity between pixels that are adjacent (within radius of 4 pixels) and using sparse matrices. I then get equal clustering for the entire image. To make a sparse adjacency matrix I use Shai's function sparse_adj_matrix 之后，否则内存仅由邻接矩阵填充。
tic
complexrowsTranspose = [real(complexrows);imag(complexrows)]';
indexnonzero = find(mean(tempdisk,3) > 0);
idxval = zeros(size(complexrowsTranspose, 1),1);
[irow jcol] = sparse_adj_matrix([size(meannorm,1), size(meannorm,2)], 4, 2);
keep = ismember(irow, indexnonzero);
irow(~keep) = [];
jcol(~keep) = [];
clear keep
sigma = 1;
Avect = exp(-sum((complexrowsTranspose(irow,:)-complexrowsTranspose(jcol,:)).^2,2)/(2*sigma^2));
iAval = find([Avect>0].*[Avect<Inf]);
Aff = sparse(irow(iAval),jcol(iAval),Avect(iAval),length(complexrowsTranspose),length(complexrowsTranspose));
Dvector = sum(Aff,2).^(1/2);
Dindex = find(Dvector);
D = sparse(Dindex, Dindex, Dvector(Dindex),size(Aff,1),size(Aff,2));
L = D * Aff * D;
keig = 25;
[Vect, Val] = eigs(L, keig);
normVect = Vect./repmat(sqrt(sum(abs(Vect).^2,2)), [1 keig]);
[kidx,kC,ksumd,kD] = kmeans(normVect, 5);
kmeansim = reshape(kidx, [921, 921]);
figure; imshow(kmeansim,[])
toc
这是生成的集群的样子：
看起来好多了。然而，我感兴趣的集群没有出现（'angle' 图像中具有高值的斑点，在摄影师的深色外套内）。特别是具有相似平均范数值的像素被聚类，在短系列上没有相似的变化，也没有相似的角度（复值）。
我会尝试调整输入值和邻接半径来计算相关性。
更新
我尝试只输入角度、整个复数值（并调整代码以适用于复数值）、更改计算亲和力的半径、输入 1 相关性而不是距离，但我没有将摄影师外套中的小亮角斑点聚集在一起。
然后我下载了 this 代码并尝试如下：
complexrowsTranspose = complexrows';
[icol irow] = sparse_adj_matrix([921, 921], 1, Inf);
complexrowsTminmean = complexrowsTranspose -repmat(mean(complexrowsTranspose , 2), [1, 15]);
complexrowsTstd = sqrt( std(real(complexrowsTranspose), [], 2).^2+std(imag(complexrowsT), [], 2).^2 );
complexrowsTcorr = sum(real(complexrowsTminmean(icol).*complexrowsTminmean(irow)), 2)./complexrowsTstd(irow)./complexrowsTstd(icol)./(15-1);
Asparse = sparse(irow, icol, 1-complexrowsTcorr, 921*921, 921*921);
Asparse(isnan(Asparse)) = 0;
K = AL_ICM(Asparse);
但我没有让它超越第一次迭代。我计算这些复杂向量的相关性的方式可能无法满足函数的要求。
引出问题
我有一个二维复值图像，其中包含一系列短值。我想聚类相似的像素/分割图像。有一个或多或少的静态图像和一个叠加图像，其中有一些斑点，这些斑点在短系列上具有变化的值（主要是复数的角度）。它们在图像的规范中也略显可辨。
我的第一次尝试是 k-means，但它确实根据均值进行了聚类（平均值存在差异，尤其是与周围像素相比，但时间和 angular 信息更大）。我的第二次尝试是 ICA，然后查看幅度最大的 k 个分量，它确实成功地将图像中的某些区域识别为不同，但没有识别出我感兴趣的像素组（视觉上不难识别）认识他们，但他们很小）。
现状
所以因为我的前两次尝试没有成功，所以我用 google 环顾四周，似乎 spectral clustering 可能是合适的。但是我在使用该方法时遇到了一些严重的问题，主要与可用内存有限有关。然后我想，因为我有这么多像素，我可以将光谱聚类应用于单独的数据块。
有人 here suggests clustering slabs first and then combine them, he then says 'at the end you will have the problem of recombining them and this problem can be solved easily'. The bits designated as 'easy' in explanations are hardly ever easy of course. He links to this 论文，但是那个方法并没有处理 slabs 中的所有数据。它排除了不接近主成分的向量。
问题
我的问题分为两部分：
1。如何合并单独细分的结果？特征向量不同，簇数不同。结果看起来它在单独的板中工作。
2。不考虑单独平板中像素之间的距离/亲和力。我可以制作'slabs between slabs'吗？对于那些 L 和 A 不对称的板，不知道如何执行该方法。也许我可以以某种方式比较/合并最后的所有特征向量？
(3.有没有类似或者更好的方法，不需要那么多内存，计算时间也能接受，容易爆)
Matlab代码示例
%% generate data
% get some outer region without data
tempdisk = strel('disk',922/2); tempdisk = double(repmat((1+sqrt(-1)).*tempdisk.Neighborhood,[1 1 15]));
% make some noise
tempnoise = (rand(921,921,15)+sqrt(-1).*rand(921,921,15))./10;
% 'background signal'
tempim1 = double(imresize(mean(imread('cameraman.tif'),3),[921,921])); tempim1 = repmat(tempim1./max(tempim1(:)),[1 1 15]);
% 'target signal'
tempim2 = double(rgb2hsv(imread('fabric.png'))); tempim2 = imresize(tempim2(:,:,2),[921,921]); tempim2 = repmat(tempim2./max(tempim2(:)),[1 1 15]);
sin1 = repmat(permute(sin(2.*pi.*(0:14)./15),[1 3 2]),[921 921 1]);
% combine into data
complexdata = (sin1.*(1.0.*tempim1+0.5.*tempim2.^2).*exp(-sqrt(-1).*2.*pi.*sin1.*(tempim2.^2)).*tempdisk+tempnoise)./1.5;
%this is what the mean data looks like
meannorm = mean(abs(complexdata),3);
meanangle = mean(angle(complexdata),3);
figure; subplot(1,2,1); imshow(meannorm,[]); title('mean norm'); subplot(1,2,2); imshow(meanangle,[]); title('mean angle')
这是生成的数据的样子：
右图中明亮的斑点就是我想要的。它们随时间的变化也最强（并且在时间上相关）。
然后设置集群：
%% perform spectral clustering in seperate slabs
% method from http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf
%get all pixel vectors in a single matrix
complexrows = reshape(permute(complexdata, [3,1,2]), [15, 921*921]);
%k means and eigs dont accept complex, so convert to real here?
complexrowsTranspose = [real(complexrows);imag(complexrows)]';
%lets say 10000 by 10000 matrices are still ok
npix = 10000;
nslabpix = floor(length(complexrowsTranspose)/npix);
nrestpix = rem(length(complexrowsTranspose), npix);
在适合内存的平板中执行光谱聚类：
% spectral clustering
keig = 50;%how many eigenvectors needed? more is better
affinity_sigma = 1;% i dont understand how to calculate this from the paper
tic
% start with last slab (dynamically preallocate)
for islabpix = (nslabpix+1):-1:1;
%print progress
islabpix/nslabpix
toc
if islabpix>nslabpix
pixrange = (1:nrestpix) + ((islabpix-1)*npix);;
else
pixrange = (1:npix) + ((islabpix-1)*npix);
end
%calculate affinity between all voxels in slab
Aff = exp( -squareform(pdist(complexrowsTranspose(pixrange,:))).^2/(2*affinity_sigma^2) ); % affinity matrix
%calculate degree matrix for normalization
Dsq = sparse(size(Aff,1),size(Aff,2)); %degree matrix
for idiag=1:size(Aff,1)
Dsq(idiag,idiag) = sum(Aff(idiag,:))^(1/2);
end
%normalize affinity matrix
Lap = Dsq * Aff * Dsq; %normalize affinity matrix
%calculate eigenvectors of affinity matrix
[eigVectors(pixrange,1:keig), eigValues] = eigs(Lap, keig); %eigenvectors of normalized aff mat
normEigVectors(pixrange,1:keig) = eigVectors(pixrange,1:keig)./repmat(sqrt(sum(abs(eigVectors(pixrange,1:keig)).^2,2)), [1 keig]); %normalize rows of eigen vectors, normr only works on real numbers
% perform k means clustering on weights for eigenvectors
[idx,C,sumd,D] = kmeans([real(normEigVectors(pixrange,1:keig)),imag(normEigVectors(pixrange,1:keig))], 5); %k means on normalized eigenvecotrs
idxval(pixrange) = idx;
end
%reshape into image
idxim = reshape(idxval, [921, 921]);
figure; imshow(idxim,[])
toc
生成的聚类：
结果看起来该方法在每个平板中都在某种程度上起作用；目标是将所有具有稍高范数和更强角度变化的斑点（来自 tempim2 的高饱和度斑点）聚集在一起，这在结果中似乎是可识别的。现在问题主要是单独的平板，并且没有交叉平板集群。这花了我的电脑大约 15 分钟。我减少了本示例的特征值数量和图像大小，因此它可以在可接受的时间内运行。我认为这说明了我的部分问题。
我真的没有答案，但我认为这些指示应该可以帮助您找到答案：
您声称有记忆问题。你确定你的亲和力矩阵是稀疏的吗？您的代码中似乎只有对角线度矩阵是稀疏的。通常当 运行 在 pixels/voxels 上进行谱聚类时，将亲和矩阵设计得非常稀疏（8 连接或 26 连接）。
您将集群描述为 "they are small"。光谱聚类有 known issues 个不同尺度的聚类。你确定你得到了满意的结果吗？
如何计算相邻体素之间的亲和力（相似度）？你也可以测量 差异性 吗？也就是说，你能说一些体素应该不属于同一个簇吗？如果是这样，您是否考虑过使用 correlation clustering？该方法对不同的簇尺度鲁棒性更强，可以自动检测簇数
您是否考虑过使用多尺度/multigrid 方法来粗化您的数据，而不是将其粗暴地分割成 "slabs"？
你看过spectralNet了吗？如果我没记错的话，这个方法应该可以让你"learn"对部分点进行谱聚类，然后使用网络"extrapolate"对新点进行聚类
更新：
鉴于
关于谱聚类和相关聚类的区别：
为什么谱聚类能够聚类所有点，即使输入亲和力矩阵如此稀疏？它如何分辨出那个点a和一个很远的点 c 应该属于同一个集群，即使它们之间没有计算亲和力？
简单的答案是亲和力的 transitivity：如果 a 类似于 b 并且 b 类似于 c 那么 a 和 c 应该聚集在一起。
问题在哪里？ 在谱聚类中，亲和矩阵中的所有条目都是非负的，这意味着除非绝对没有连接 a 和 c 的路径（机会很小）有一些 "transitive affinity" 建议 a 和 c 应该属于同一个集群。因此，如果查看谱聚类的数学，您会注意到 "trivial solution"，即将所有点放在同一个聚类中，为问题提供了全局最优解。必须人为地强制解决方案具有 k 个集群以避免琐碎的解决方案。
可以做什么？ 如果只考虑正相关性，值 0 是不明确的：它表示 "I didn't bother to compute the affinities between these points"，但也可以表示 "I think these two points should not be in the same cluster"。为了克服这种歧义，我们可以这样引入 negative affinities 如果 A(i, j) > 0 意味着点 i 和点 j 应该在同一个集群中并且确定 A(i, j)，而如果 A(i, j) < 0 意味着 i 和 j 应该 而不是 在同一个集群中（肯定 |A(i, j)|） .引入负亲和力打破了可能 link 远点的 "transitivity chains"，不，将所有点放在同一个集群中不再是微不足道的。
如何利用负亲和力？ 当你的亲和力矩阵同时具有正（吸引力）和负（排斥）值时，你可以使用相关聚类来聚类这些点，它基本上试图最大化affinities/attraction 在每个簇内的点之间，同时最大化不同簇中的点之间的排斥力。一个很好的 属性 相关聚类是它 "automatically" 发现聚类的基础数量，见 sec。 this paper 中的 2 个。
在 sparse_adj_matrix 之后，否则内存仅由邻接矩阵填充。
tic
complexrowsTranspose = [real(complexrows);imag(complexrows)]';
indexnonzero = find(mean(tempdisk,3) > 0);
idxval = zeros(size(complexrowsTranspose, 1),1);
[irow jcol] = sparse_adj_matrix([size(meannorm,1), size(meannorm,2)], 4, 2);
keep = ismember(irow, indexnonzero);
irow(~keep) = [];
jcol(~keep) = [];
clear keep
sigma = 1;
Avect = exp(-sum((complexrowsTranspose(irow,:)-complexrowsTranspose(jcol,:)).^2,2)/(2*sigma^2));
iAval = find([Avect>0].*[Avect<Inf]);
Aff = sparse(irow(iAval),jcol(iAval),Avect(iAval),length(complexrowsTranspose),length(complexrowsTranspose));
Dvector = sum(Aff,2).^(1/2);
Dindex = find(Dvector);
D = sparse(Dindex, Dindex, Dvector(Dindex),size(Aff,1),size(Aff,2));
L = D * Aff * D;
keig = 25;
[Vect, Val] = eigs(L, keig);
normVect = Vect./repmat(sqrt(sum(abs(Vect).^2,2)), [1 keig]);
[kidx,kC,ksumd,kD] = kmeans(normVect, 5);
kmeansim = reshape(kidx, [921, 921]);
figure; imshow(kmeansim,[])
toc
这是生成的集群的样子：
看起来好多了。然而，我感兴趣的集群没有出现（'angle' 图像中具有高值的斑点，在摄影师的深色外套内）。特别是具有相似平均范数值的像素被聚类，在短系列上没有相似的变化，也没有相似的角度（复值）。
我会尝试调整输入值和邻接半径来计算相关性。
更新
我尝试只输入角度、整个复数值（并调整代码以适用于复数值）、更改计算亲和力的半径、输入 1 相关性而不是距离，但我没有将摄影师外套中的小亮角斑点聚集在一起。
然后我下载了 this 代码并尝试如下：
complexrowsTranspose = complexrows';
[icol irow] = sparse_adj_matrix([921, 921], 1, Inf);
complexrowsTminmean = complexrowsTranspose -repmat(mean(complexrowsTranspose , 2), [1, 15]);
complexrowsTstd = sqrt( std(real(complexrowsTranspose), [], 2).^2+std(imag(complexrowsT), [], 2).^2 );
complexrowsTcorr = sum(real(complexrowsTminmean(icol).*complexrowsTminmean(irow)), 2)./complexrowsTstd(irow)./complexrowsTstd(icol)./(15-1);
Asparse = sparse(irow, icol, 1-complexrowsTcorr, 921*921, 921*921);
Asparse(isnan(Asparse)) = 0;
K = AL_ICM(Asparse);
但我没有让它超越第一次迭代。我计算这些复杂向量的相关性的方式可能无法满足函数的要求。