<p>I don't really have an answer for you, but I think these pointers should help you find an answer:</p>

<ol>
<li><p>You claim to have memory problems. Are you sure your affinity matrix is sparse? It seems like only the diagonal degree matrix is sparse in your code. Usually when running spectral clustering on pixels/voxels one designs the affinity matrix to be very sparse (8 connect or 26 connect).</p></li>
<li><p>You describe your clusters as "they are small". Spectral clustering has <a href="http://wpre.weizmann.ac.il/math/Nadler/sites/math.Nadler/files/publications/fundamental_limitations_of_spectral_clustering.pdf" rel="nofollow noreferrer">known issues</a> with clusters in very different scales. Are you sure you are getting satisfactory results?</p></li>
<li><p>How do you compute the affinity (similarity) between neighboring voxels? Can you measure <em>dissimilarity</em> as well? That is, can you say for some voxels that they should <em>not</em> belong to the same cluster? If so, have you considered using <a href="https://stackoverflow.com/a/19510366/1714410">correlation clustering</a>? This method is more robust to different cluster scales and can <em>automatically</em> detect the number of clusters.</p></li>
<li><p>Have you considered using multiscale/<a href="http://arxiv.org/abs/1204.4867" rel="nofollow noreferrer">multigrid</a> methods to coarsen your data instead of brutally slicing it into "slabs"? </p></li>
<li><p>Have you looked at <a href="https://arxiv.org/abs/1801.01587" rel="nofollow noreferrer">spectralNet</a>? If I am not mistaken, this method should enable you to "learn" the spectral clustering on part of the points and then use the net to "extrapolate" the clustering to new points.</p></li>
</ol>

<hr>

<p><strong>Upadate:</strong><br>
In light of <a href="https://stackoverflow.com/questions/50446375/how-to-combine-split-runs-of-spectral-clustering-for-a-huge-affinity-matrix/50460927?noredirect=1#comment87999026_50460927">Leo's comment</a>, I would say that when it comes to spectral clustering of very large data, brutally slicing the data into "slabs" and then trying to "stitch" the results together might not be the best coarse of action (not that I think it is not possible). A better way to approach the problem is by significantly sparsifying the affinity matrix: compute pair-wise affinities for each point only to its neighbors, resulting with affinity matrix that is mostly sparse. This way one can process all the points at once without the need to "slice" and "stitch".  </p>

<p>As for the difference between spectral clustering and correlation clustering:<br>
<em>Why spectral clustering is able to cluster <strong>all</strong> points even when the input affinity matrix is so sparse?</em> how can it tell that point <code>a</code> and a far away point <code>c</code> should belong in the same cluster even when no affinity was computed between them?<br>
The simple answer is <a href="https://en.wikipedia.org/wiki/Transitive_relation" rel="nofollow noreferrer">transitivity</a> of affinities: if <code>a</code> is similar to <code>b</code> and <code>b</code> is similar to <code>c</code> then <code>a</code> and <code>c</code> should be clustered together.<br>
<em>Where's the catch?</em> In spectral clustering all entries in the affinity matrix are non-negative, which means that unless there is absolutely no path connecting <code>a</code> and <code>c</code> (slim chance) there is some "transitive affinity" suggesting <code>a</code> and <code>c</code> should belong to the same cluster. Therefore, if look at the math of spectral clustering you'll notice that the "trivial solution", i.e., placing all points in the same cluster, provides a global optimum to the problem. One must artificially force the solution to have <code>k</code> clusters to avoid the trivial solution.<br>
<em>What can be done?</em> If you only consider positive affinities the value 0 is ambiguous: it means either "I didn't bother to compute the affinities between these points", but it can also mean "I think these two points should <em>not</em> be in the same cluster". To overcome this ambiguity we can introduce <strong>negative affinities</strong> this way if <code>A(i, j) &gt; 0</code> means point <code>i</code> and point <code>j</code> should be in the same cluster with certainty <code>A(i, j)</code>, while if <code>A(i, j) &lt; 0</code> means <code>i</code> and <code>j</code> should <strong>not</strong> be in the same cluster (with certainty <code>|A(i, j)|</code>). Introducing negative affinities breaks the "transitivity chains" that may link far away points, no it is no longer trivial to place all points in the same cluster.<br>
<em>How to take advantage of negative affinities?</em> When your affinity matrix has both positive (attraction) and negative (repulsion) values, you can cluster the points using correlation clustering which basically tries to maximize the affinities/attraction between points within each cluster and simultaneously maximize the repulsion between points in different clusters. A nice property of correlation clustering is that it "automatically" discover the underlying number of clusters, see sec. 2 of <a href="http://www.wisdom.weizmann.ac.il/~bagon/pub/LargeScaleCorrClust_2011.pdf" rel="nofollow noreferrer">this paper</a>.  </p>
How to combine split runs of spectral clustering for a huge affinity matrix
<p><strong><em>Leading up to the question</em></strong></p>

<p>I have a 2D complex valued image with a short series of values. I want to cluster similar pixels / segment the image. There is one more or less static image with a superimposed image that has some blobs in it that have a changing value (mostly the angle of the complex number) over the short series. They are also slightly discernable in the norm of the image.</p>

<p>My first attempt was k-means, but that really clustered according to the means (there is a distinction in mean values, especially compared to surrounding pixels, but the temporal and angular information is greater). My second attempt was ICA and then look at the k components with the largest magnitude, and that did successfully identify certain regions in my image as being different, but did not identify the group of pixels I was interested in (visually it is not hard to recognize them, but they are small).</p>

<p><strong><em>Current situation</em></strong></p>

<p>So because my first two tries did not work out, I looked around with google and it seemed <a href="http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf" rel="nofollow noreferrer">spectral clustering</a> might be appropriate. But I have some serious issues when using the method, mostly to do with limited available memory. I then thought, since I have so many pixels, I can just apply spectral clustering to seperate slabs of the data.</p>

<p>Someone <a href="https://www.researchgate.net/post/How_to_perform_spectral_clustering_for_large_size_imagesDoes_parallel_computing_can_solve_the_issue_of_scalability" rel="nofollow noreferrer">here</a> 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 <a href="http://file:///C:/Users/Lennart/Downloads/SpectralClusteringforaLargeDataSetbyReducingtheSimilarityMatrixSize%20(1).pdf" rel="nofollow noreferrer">this</a> paper, but that method does not process all the data in slabs. It rather excludes vectors that are not close to a principal component.</p>

<p><strong><em>Question</em></strong></p>

<p>My question has 2 parts:</p>

<p><strong>1</strong>. How do I combine the results for the seperate segments? The eigenvectors are different and the cluster numbers are different. The result looks like it worked in the seperate slabs.</p>

<p><strong>2</strong>. No distance / affinity between pixels in seperate slabs is taken into account. Can I make 'slabs between slabs'? For those slabs L and A are not symmetric, no clue how to perform the method then. Perhaps I can somehow compare / merge all eigenvectors at the end?</p>

<p>(3. Is there a similar or better method that does not need so much memory. Computation time is also borderline acceptable, easily exploding)</p>

<p><strong><em>Matlab code example</em></strong></p>

<pre><code>%% 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')
</code></pre>

<p>This is what the generated data looks like:</p>

<p><a href="https://i.stack.imgur.com/CFFXg.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/CFFXg.jpg" alt="mean data"></a></p>

<p>The bright blobs in the right image are what Im after. They have the strongest variation over time as well (and are correlated in time).</p>

<p>Then to set up the clustering:</p>

<pre><code>%% 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);
</code></pre>

<p>Perform spectral clustering in slabs that fit into memory:</p>

<pre><code>% 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&gt;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
</code></pre>

<p>The resulting clustering:</p>

<p><a href="https://i.stack.imgur.com/hLJ4m.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/hLJ4m.png" alt="k means clusters on spectral slabs"></a></p>

<p>The result looks like the method is working to some degree within each slab; the goal was to cluster all blobs with slightly higher norm and stronger angle variation (high saturation blobs from <code>tempim2</code>), which seem recognizable in the result. Now its mostly the seperate slabs that are the issue and that there are no cross-slab clusters. This took my computer about 15 minutes. I reduced the number of eigenvalues and the image size for this example so it runs in an acceptable amount of time. I think that illustrates part of my problem.</p>