Compressed sensing: CS for networked data, Diffusion wavelets and Analog CS measurements.

The Special Issue on Compressive Sampling of the IEEE Signal Processing Magazine will also feature Compressed sensing for networked data by Jarvis Haupt, Waheed Bajwa, Mike Rabbat, and Robert Nowak. In the article, one can learn how one can exploit the connectivity in the network "as the key to obtaining effective sparsifying transformations for networked data". Hence because Network traffic patterns are sparse in the diffusion wavelet framework, the idea is to use this framework to enable Compressed Sensing measurements from networks.


They also discuss the specific schemes one could implement with a particular attention to data from wireless networks and the analog production of the compressed sensing data through the use of the network. Very nice.

Compressed Sensing: How could Hardware implementations become disruptive technologies ?

I am still wondering how technology developed with Compressed Sensing could be construed at Disruptive Technology, here is reminder of what disruptive technology is:

"...Disruptive technologies do not destroy existing market leaders overnight. They do not get adopted by the entire market at the same time. They do not initially seem to be "better" products (in fact, in the early going, they are often distinctly "worse.") They are not initially a viable option for mainstream users. They do not win head-to-head feature tests. Initially, they do not even seem to be a threat.

Disruptive technologies take advantage of a new manufacturing/business process or technology to provide a cheaper, more convenient, simpler solution that meets the needs of the low end of the market. Low-end users don't need all the features in the Incumbent's product, so they rapidly adopt the simpler solution. Meanwhile, the Incumbent canvasses its mainstream customers, reassures itself that they want the feature-rich products, and dismisses the Disruptor as a niche player in an undesirable market segment. The Incumbent may dabble with the new technology or process, but only in service of its existing business model.

Then the Disruptor improves its products, adding more features while keeping the convenience and low cost. Now the product appeals to more mainstream users, who adopt it not because it's "better" but because it's simpler and cheaper. Seeing this, the Incumbent continues adding ever more features and functionality to its core product to try to maintain its value proposition for higher end customers. And so on. Eventually, the Incumbent's product overshoots the needs of the mass market, the Disruptor grabs the mainstream customers, and, lo and behold, the technology has been "disrupted."

In a previous entry, I was wondering How could Compressed Sensing be seen as a disruptive technology ?
and used some guidance given by Todd Proebsting when he was giving a talk on innovation in programming languages at LL2.

First of all, Compressed Sensing is really a technique not a technology, so I will try to address the technologies that are beginning to use CS at and see if those can be considered disruptive technologies.

1. In the case of MRI, Compressed Sensing should clearly not be labeled as disruptive as it has clearly been adopted by the large MRI companies to increase the throughput of their machines.
2. In the Rice single pixel camera case, there is currently a clear advantage in the realm of imaging parts of the light spectrum where the development many sensitive pixels is clearly too expensive. The other interesting aspect of the single pixel imager is the ability to remove the compression scheme which takes most of the battery. But in an age where batteries are seen as the technology to improve, should an imaging hardware be developed just for the purpose of being less power hungry ?
   
3. In the case of hyperspectral imagers, any current hyperspectral imager cost between 70K$ to 100K$. A low cost product while not being very efficient (accurate in certain bands) would clearly make the CS hyperspectral imager as a disruptive technology. One of the major underlying reason as to why CS hyperspectral imager have a reasonable possibility of creating a niche market is the ability to provide a large field of view and dynamic range compared to current technology.
   

For each of these technologies I will try to address the following in the future:

• Opportunity: what’s the issue?
  
• Current solutions: what’s done now
• Proposal: New disruptive technology
  
• Disadvantages why some (many?) will scoff
• Unmet needs benefits to adopters:

I am sure most of you can make out what the "many will scoff" will be made off ?

Monday Morning Algorithm 8: An Example of Compressed Measurements Speeding Up Diffusion Maps

In the previous installment of the Monday Morning Algorithm, we featured the diffusion method. I also mentioned that one could use compressed sensing and obtain the same result: i.e the ability to sort several high dimensional images because only one parameter is being varied. When going through this exercise, one can note that of the most expensive operation is the inter-distance measured between all the elements of the set. If every picture is 10 MP, one has to evaluate the difference between elements made of ten million numbers. In Compressed Sensing, one learns that projecting a large image on the right kind of random projections, one can keep the amount of information in the original image. We are going to reuse our last algorithm but we will be projecting the original images of Buzz Lightyear on a smaller space. We will be using a Gaussian ensemble and will be projecting from an image of dimension 57600 down to 100 Compressed Measurements. Let us note that this not an efficient operation unless we have directly access to a compressed measurements or if one does extra operations with the compressed measurements.

clear
% Stephane Lafon, “Diffusion maps and geometric harmonics”,
% Ph.D. dissertation, Yale University (May
% 2004)
% written by Cable Kurwitz and Igor Carron
% This algorithm shows how random projections
% used in compressed sensing provides the same
% results as the original images.
%
% Number of Random Projections
number_RP =100;
% get all the files in jpg
d = dir('buzz*.jpg')
ntotal = size(d,1)
for i=1:ntotal
d(i).name
end
Img = imread(d(1).name);
s=size(Img,1);
w=size(Img,2);
hd=size(Img,3);
swhd=s*w*hd;
%reshape these images into one vector each
for ii=1:ntotal
Img = imread(d(ii).name);
b1(:,ii) = double(reshape(Img, swhd,1));
end
%
% Compressed sensing part
sa=randn(number_RP,swhd);
for i=1:swhd
sa(:,i)=sa(:,i)./norm(sa(:,i));
end
b=sa*b1;
%b=b1;
%
%
for i=1:ntotal
for j=1:ntotal
D1(i,j)=norm(b(:,i)-b(:,j),2);
end
end
%
n1=ntotal;
D_sum=0;
for l = 1:n1;
D_sum = D_sum + min(nonzeros(D1(l,:)));
end;
epsilon = D_sum/n1/10; %
%
% Step 1
%
kernel_1 = exp(-(D1./epsilon).^1);
%
% Step 2
%
one_V = ones(1,n1);
%
% Step 3
%
p = kernel_1*one_V';
kernel_2 = kernel_1./((p*p').^1);
%
% Step 4
%
v = sqrt(kernel_2*one_V');
%
% Step 5
%
K = kernel_2./(v*v');
%
% Step 6
%
%% Singular Value Decomposition
[U S V] = svd(K);
% Normalization
for i=1:ntotal-1
phi(:,i) = U(:,i)./U(:,1); %
end
eig_funct = phi(:,1:6);
%
% The eigenvalues of \delta are approximated by those of K, and its
% eigenfunctions Ái are approximated by phi(:,i) = U(:; i):=U(:; 1)
%
vv=eig(K);

% Initial order of Buzz images
figure(1)
for i=1:ntotal
subplot(3,4,i)
imshow(d(i).name)
end

X = vv(ntotal-1).*phi(:,2);
% sorting the Buzz images using the first
% eigenfunctions the Laplace Beltrami operator
[Y,I] = sort(X);
% d(I).name
%
figure(2)
for i=1:ntotal
subplot(3,4,i)
imshow(d(I(i)).name)
end
%
figure(3)
plot(vv(ntotal-1).*phi(:,2),'.')
title(' Buzz Lightyear')
xlabel('Index')
ylabel('Lambda_1 phi_1')
%
figure(4)
plot(vv(ntotal-1).*phi(:,2),vv(ntotal-2).*phi(:,3),'.')
title(' Buzz Lightyear')
xlabel('Lambda_1 phi_1')
ylabel('Lambda_2 phi_2')
%
figure(5)
plot(eig(K),'*')
xlabel('Index of eigenvalue')
ylabel('Value of eigenvalue')
title(' Eingenvalues of K')
vv=eig(K);
%

It looks like 100 projections are enough to recover the same results as the original analysis. One can note that this process of figuring out what is the number of random projections needed would require the type of analysis evaluating the intrinsic dimensionality of the set as mentioned here by the folks at Rice.

All the Monday Morning Algorithms (and attendant codes/files) are listed here.

Compressed Sensing: Fast Compressive Sampling with Structurally Random Matrices

Thong Do mentions that the paper presented here is already on the Rice website. Thong Do, Trac Tran and Lu Gan are co-authors of the fascinating Fast Compressive Sampling with Structurally Random Matrices which seems to lead to highly sparse measurement matrices. The abstract reads:

This paper presents a novel framework of fast and efficient compressive sampling based on the new concept of structurally random matrices. The proposed framework provides four important features. (i) It is universal with a variety of sparse signals. (ii) The number of measurements required for exact reconstruction is nearly optimal. (iii) It has very low complexity and fast computation based on block processing and linear filtering. (iv) It is developed on the provable mathematical model from which we are able to quantify trade-offs among streaming capability, computation/memory requirement and quality of reconstruction. All currently existing methods only have at most three out of these four highly desired features. Simulation results with several interesting structurally random matrices under various practical settings are also presented to verify the validity of the theory as well as to illustrate the promising potential of the proposed framework.

I cannot wait to see how it will fare with regards to reconstruction capabilities, i.e the ability of having very sparse measurement matrices and still allow for few measurements and good reconstruction. I need to look into this deeper.

Compressed Sensing: MCA news and Sublinear Recovery of Sparse Wavelet Signals

Coming back from NIPS 07, Jort Gemmeke mentioned the following two items:

Morphological Component Analysis : an Adaptive Thresholding Strategy by Jerome Bobin, Jean-Luc Starck, Jalal Fadili, Yassir Moudden, and David Donoho, the abstract reads:

In a recent paper, a method called MCA (Morphological Component Analysis) has been proposed to separate the texture from the natural part in images. MCA relies on an iterative thresholding algorithm, using a threshold which decreases linearly towards zero along the iterations. This paper shows how the MCA convergence can be drastically improved using the mutual incoherence of the dictionaries associated to the different components. This modified MCA algorithm is then compared to Basis Pursuit, and experiments show that MCA and BP solutions are similar in terms of sparsity, as measured by l1 norm, but MCA is much faster and gives us the possibility to handle large scale data sets.

The MOM/matched filtering and thresholding bears some similarity with another paper mentioned in this entry from some of the same co-authors. I haven't looked at it deeper yet. Always on the same topic of Morphological Component Analysis, Jort also mentions the release of the MCA 802 C++ toolbox from Jalal Fadili's webpage.

Toward a dedicated Compressed Sensing based on the signal signature? It looks like Ray Maleh and Anna Gilbert are beginning to do that in Sublinear recovery of sparse wavelet signals. The abstract reads:

There are two main classes of decoding algorithms for “compressed sensing,” those which run time time polynomial in the signal length and those which use sublinear resources. Most of the sublinear algorithms focus on signals which are compressible in either the Euclidean domain or the Fourier domain. Unfortunately, most practical signals are not sparse in either one of these domains. Instead, they are sparse (or nearly so) in the Haar wavelet system. We present a modified sublinear recovery algorithm which utilizes the recursive structure of Reed-Muller codes to recover a wavelet-sparse signal from a small set of pseudo-random measurements. We also discuss an implementation of the algorithm to illustrate proof-of-concept and empirical analysis.

Compressed Sensing, Compressive Sampling: A Customized Search Engine

One of the reason I link to the web pages of people involved in Compressed Sensing is twofold. It is designed to:

• see what their other interests are and,
• help Google customized search engine focus on the subject of Compressed Sensing.

If you use the search box below (also on the right side) you will specifically search through all the pages of these researchers and the webpages they link to. To the students writing papers on the subject, please make a web page (even a tiny placeholder), it will go a long toward establishing an hopefully long lasting presence on the web in the future.

Compressed Sensing: Intrinsic Dimensionality of a Dataset, StOMP for large scale problems and 2007 WD5 is aiming at Mars

In Random Projections for Manifold Learning, Chinmay Hegde, Mike Wakin, and Richard Baraniuk try to devise a means by which one can determined the intrinsic dimensionality of a dataset. This is important as it has a direct impact on the size of the measurement matrix that one wants to be as small as possible. It uses the Grassberger-Procaccia algorithm to estimate the intrinsic dimension of the manifold using random projections of the original dataset. The abstract reads:

We derive theoretical bounds on the performance of manifold learning algorithms, given access to a small number of random projections of the input dataset. We prove that with the number of projections only logarithmic in the size of the original space, we may reliably learn the structure of the nonlinear manifold, as compared to performing conventional manifold learning on the full dataset.

In the NIPS version of that technical report, the abstract reads:

We propose a novel method for linear dimensionality reduction of manifold modeled
data. First, we show that with a small number M of random projections of sample points in RN belonging to an unknown K-dimensional Euclidean manifold, the intrinsic dimension (ID) of the sample set can be estimated to high accuracy. Second, we rigorously prove that using only this set of random projections, we can estimate the structure of the underlying manifold. In both cases, the number of random projections required is linear in K and logarithmic in N, meaning that K less than M much less than N. To handle practical situations, we develop a greedy algorithm to estimate the smallest size of the projection space required to perform manifold learning. Our method is particularly relevant in distributed sensing systems and leads to significant potential savings in data acquisition, storage and transmission costs.

in the paper the algorithm is given in an explicit way and requires ISOMAP. These two papers are important for Machine Learning as they now give a faster mean of estimating the underlying manifold using random projections.


David Donoho, Yaakov Tsaig, Iddo Drori, and Jean-Luc Starck, Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit, the abstract reads:

Finding the sparsest solution to underdetermined systems of linear equations y = x is NP-hard in general. We show here that for systems with ‘typical’/‘random’ , a good approximation to the sparsest solution is obtained by applying a fixed number of standard operations from linear algebra. Our proposal, Stagewise Orthogonal Matching Pursuit (StOMP), successively transforms the signal into a negligible residual. Starting with initial residual r0 = y, at the s-th stage it forms the ‘matched filter’ T rs−1, identifies all coordinates with amplitudes exceeding a specially-chosen threshold, solves a least-squares problem using the selected coordinates, and subtracts the least squares fit, producing a new residual. After a fixed number of stages (e.g. 10), it stops. In contrast to Orthogonal Matching Pursuit (OMP), many coefficients can enter the model at each stage in StOMP while only one enters per stage in OMP; and StOMPtakes a fixed number of stages (e.g. 10), while OMP can take many (e.g. n). StOMP runs much faster than competing proposals for sparse solutions, such as l1 minimization and OMP, and so is attractive for solving large-scale problems.
We use phase diagrams to compare algorithm performance. The problem of recovering a k-sparse vector x0 from (y, ) where  is random n × N and y = x0 is represented by a point (n/N, k/n) in this diagram; here the interesting range is k less than n less than N. For n large, StOMP correctly recovers (an approximation to) the sparsest solution of y = x over a region of the sparsity/indeterminacy plane comparable to the region where l1 minimization is successful. In fact, StOMPoutperforms both l1 minimization and OMP for extremely underdetermined problems. We rigorously derive a conditioned Gaussian distribution for the matched filtering coefficients at each stage of the procedure and rigorously establish a large-system limit for the performance variables of StOMP. We precisely calculate large-sample phase transitions; these provide asymptotically precise limits on the number of samples needed for approximate recovery of a sparse vector by StOMP. We give numerical examples showing that StOMPrapidly and reliably finds sparse solutions in compressed sensing, decoding of error-correcting codes, and overcomplete representation.

Always pay attention when Dave Donoho says something. I also found the qualifying exams slide of Deanna Needell interesting, especially slide 116/145 an information I did not know. I am also impatient to see what Thong Do a student of Trac Tran is talking about in this seminar:

In this work, we propose a novel fast compressive sampling algorithm that is based on a new concept of structurally random ensembles. The new proposed algorithm not only enables the system to compressively sample signals very fast but also requires very low complexity with other highly desired features for practical applications such as: integer friendliness, streaming capability. Surprisingly, this novel algorithm is also proposed to solve a closely related problem of high dimensional reduction in theoretical computer sciences.

I mentioned earlier the Titan flyby of today, but guess what kids, this one trumps it: 2007 WD5, an asteroid is going to hit Mars where we currently have five robots ( Spirit and Opportunity, Odyssey, Mars Express and the Mars Reconnaissance Orbiter). Woohoo, where is my popcorn ?.

Compressed Sensing: Compressive sampling vs conventional imaging, the Birth of the Everlasting Dataset

I just found this presentation entitled Compressive Sensing: A New Framework for Imaging by Richard Baraniuk, Justin Romberg and Robert Nowak which seems to be a tutorial presentation for ICIP 2006. Much of this tutorial appear in previous presentations except for the presentation part starting at slide 93 to slide 120, I have found this part of the presentation very refreshing as I am beginning to see the connection between results of spectral theory of Random Matrices and Compressed Sensing. At some point the subject will merge with this one but I am not sure how for the time being.

Another aspect of the image reconstruction, which I had not seen before, is also explicitly presented. Once you have retrieved random measurements/projections, Compressed Sensing techniques enable you to reconstruct an image with different bases. If, in 20 years, we have better bases, the same measurements can be used and will reproduce a better image than the one we can currently reconstruct. Case in point, if we had done compressed sensing in the mid-1990's, we would have reconstructed images with edges using wavelets bases. Then curvelets came up , using the same measurements we would have been able to get a clearer imageright when Emmanuel Candes discovered the basis!

The underlying reason being that curvelets provide better approximation properties to edges, hence providing a sparse representation of images (most natural images are made of edges).

The presentation (excerpted from reference [1]) then goes on to show a comparison between recovery using Compressed Sensing with different bases and a normal pixel sampling.


The wedgelets are coming from this implementation. The conclusions are nice summary:

"Compressive sampling techniques can yield accurate reconstructions even when the signal dimension greatly exceeds the number of samples, and even when the samples themselves are contaminated with significant levels of noise. CS can be advantageous in noisy imaging problems if the underlying image is highly compressible or if the SNR is sufficiently large"

I can see how there is going to be similar activity in the field of data mining where the underlying data dimension is larger than 2 and for which there is not currently a reasonable basis to decompose datasets on. Images are 10 MB large at most and get reduced to about 1 MB jpeg (this gives an idea of the compression using the DCT). Compressed Sensing with current bases yield about 4 MB of random measurements a reduction of about 60 %. With this in mind, one wonders how many random measurements would be needed to eventually provide a good reconstruction on different types of datasets. Does taking 47 random measurements of Shakespeare's Hamlet provide enough information to reconstruct it ? Maybe not today, if not, when ?

Reference: [1] Jarvis Haupt and Robert Nowak, Compressive sampling vs conventional imaging

Compressed Sensing Thoughts, Diffusion maps methods for Brain Imaging and for Computer Generated Graphics

Muthu Muthukrishnan is wondering aloud about the substring matching problem and how it can be used to find the reconstruction of a signal made of subsignals that are sparse using the compressed sensing framework.

I also found items of interest and relevant to dimensionality reduction applied to different domains. The first one is an article on Recovering cerebral white matter structures with Spectral Clustering of Diffusion MRI Data by Demian Wassermann, Maxime Descoteaux and Rachid Deriche.
The abstract reads:

White matter fiber clustering allows to get insight about anatomical structures in order to generate atlases, perform clear visualizations and compute statistics across subjects, all important and current neuroimaging problems. In this work, we present a Diffusion Maps clustering method applied to diffusion MRI in order to cluster and segment complex white matter fiber bundles. It is well-known that Diffusion Tensor Imaging (DTI) is restricted in complex fiber regions with crossings and this is why recent High Angular Resolution Diffusion Imaging (HARDI) such has Q-Ball Imaging (QBI) have been introduced to overcome these limitations. QBI reconstructs the diffusion orientation distribution function (ODF), a spherical function that has its maxima agreeing with the underlying fiber populations. In this paper, we introduce the usage of the Diffusion Maps technique and show how it can be used to directly cluster set of fiber tracts, that could be obtained through a streamline tractography for instance, and how it can also help in segmenting fields of ODF images, obtained through a linear and regularized ODF estimation algorithm based on a spherical harmonics representation of the Q-Ball data. We first show the advantage of using Diffusion Maps clustering over classical methods such as N-Cuts and Laplacian Eigenmaps in both cases. In particular, our Diffusion Maps requires a smaller number of hypothesis from the input data, reduces the number of artifacts in fiber tract clustering and ODF image segmentation and automatically exhibits the number of clusters in both cases by using an adaptive scale-space parameter. We also show that our ODF Diffusion Maps clustering can reproduce published results using the diffusion tensor (DT) clustering with N-Cuts on simple synthetic images without crossings. On more complex data with crossings, we show that our ODF-based method succeeds to separate fiber bundles and crossing regions whereas the DT-based methods generate artifacts and exhibit wrong number of clusters. Finally, we illustrate the potential of our approach on a real brain dataset where we successfully segment well-known fiber bundles.

It is pretty long but worth reading.

Also, Gabriel Peyre makes available his course note and some matlab programs to help out on the topic of Data Processing Using Manifold Methods. This course seems specific to the use of dimensionality reduction technique like diffusion maps to do some manipulation of 3-D graphics. The presentation of his course reads:

This course is an introduction to the computational theory of manifolds. Manifold models arise in various area of mathematics, image processing, data mining or computer science. Surfaces of arbitrary dimension can be used to model non-linear datasets that one encounters in modern data processing. Numerical methods allow to exploit this geometric non-linear prior in order to extract relevant information from the data. These methods include in particular local differential computations (related to the Laplacian operator and its variants) and global distance methods (related to geodesic computations). In this course, you will learn how to perform differential and geodesic computations on images, volumes, surfaces and high dimensional graphs.

The course includes a set of Matlab experiments. These experiments give an overview of various tasks in computer vision, image processing, learning theory and mesh processing. This includes computation of shortest paths, Voronoi segmentations, geodesic Delaunay triangulations, surface flattening, dimensionality reduction and mesh processing.

One can notice that the flattening of meshes has a parallel with methods like the MVU, except that in the MVU case, additional interdistance constraints have to be met.

Then again in CG, people are ok with "it's good enough" :-)

Monday Morning Algorithm 7: Diffusion Maps for Dimensionality Reduction and Manifold Parametrization.

There aren't that many dimensionality reduction approaches that also try to provide some parametrization of the manifold. The diffusion map approach is one of them. It tries to approximate a Laplace-Beltrami operator on the manifold. This is computed through the algorithm of page 33 of Stephane Lafon's dissertation [1]. Extensions can be found here and interactive animated examples can be found here. For today's algorithm, we will use images of Buzz Lightyear rotating around a single axis. The file is here. Obviously this is an unoptimized code. For better implementation of diffusion maps and other dimensionality reduction technique, this site is a good choice .

What do we do here ? we first compute the inter distance between each frame, then plug this into a kernel which eventually provides an idea of how connected each frame is to each other. Using a thermal heat transfer analogy, a frame is a node and the new distance between any one node to another one is inversely proportional the number of nodes between the two nodes. The new metric induced by this new distance enables a new parametrization of the manifold. In our case, the manifold is made of one object rotating around one axis. Parametrization of the underlying manifold should yield a one parameter family since only one dimension is being varied between each frame (one rotation). In order to show this, we use the first eigenfunction of the Laplace Beltrami and sort each frame according to its projection on that eigenfunction. The sorting along that dimension permits one to sort Buzz Lightyear rotating continuously from one end to the other.

clear
% Stephane Lafon, “Diffusion maps and geometric harmonics”,
% Ph.D. dissertation, Yale University (May
% 2004)
% written by Cable Kurwitz and Igor Carron
%
% get all the files in jpg
d = dir('buzz*.jpg')
ntotal = size(d,1)
for i=1:ntotal
d(i).name
end
Img = imread(d(1).name);
s=size(Img,1);
w=size(Img,2);
hd=size(Img,3);
swhd=s*w*hd;
%reshape these images into one vector each
for ii=1:ntotal
Img = imread(d(ii).name);
b1(:,ii) = double(reshape(Img, swhd,1));
end
%
% Compressed sensing part
%sa=randn(20,swhd);
%b=sa*b1;
b=b1;
%
%
for i=1:ntotal
for j=1:ntotal
D1(i,j)=norm(b(:,i)-b(:,j),2);
end
end
%
n1=ntotal;
D_sum=0;
for l = 1:n1;
D_sum = D_sum + min(nonzeros(D1(l,:)));
end;
epsilon = D_sum/n1/10; %
%
% Step 1
%
kernel_1 = exp(-(D1./epsilon).^1);
%
% Step 2
%
one_V = ones(1,n1);
%
% Step 3
%
p = kernel_1*one_V';
kernel_2 = kernel_1./((p*p').^1);
%
% Step 4
%
v = sqrt(kernel_2*one_V');
%
% Step 5
%
K = kernel_2./(v*v');
%
% Step 6
%
%% Singular Value Decomposition
[U S V] = svd(K);
% Normalization
for i=1:ntotal-1
phi(:,i) = U(:,i)./U(:,1); %
end
eig_funct = phi(:,1:6);
%
% The eigenvalues of \delta are approximated by those of K, and its
% eigenfunctions Ái are approximated by phi(:,i) = U(:; i):=U(:; 1)
%
vv=eig(K);

% Initial order of Buzz images
figure(1)
for i=1:ntotal
subplot(3,4,i)
imshow(d(i).name)
end

X = vv(ntotal-1).*phi(:,2);
% sorting the Buzz images using the first
% eigenfunctions the Laplace Beltrami operator
[Y,I] = sort(X);
% d(I).name
%
figure(2)
for i=1:ntotal
subplot(3,4,i)
imshow(d(I(i)).name)
end
%
figure(3)
plot(vv(ntotal-1).*phi(:,2),'.')
title(' Buzz Lightyear')
xlabel('Index')
ylabel('Lambda_1 phi_1')
%
figure(4)
plot(vv(ntotal-1).*phi(:,2),vv(ntotal-2).*phi(:,3),'.')
title(' Buzz Lightyear')
xlabel('Lambda_1 phi_1')
ylabel('Lambda_2 phi_2')
%
figure(5)
plot(eig(K),'*')
xlabel('Index of eigenvalue')
ylabel('Value of eigenvalue')
title(' Eingenvalues of K')
vv=eig(K);
%

Please note the ability to test the same algorithm by reducing first the frames through random projections (one needs to delete some % ). It is pretty obvious that one would need a robust metric in the first place instead of L2, one could foresee the use of a specifically designed metric (Improving Embeddings by Flexible Exploitation of Side Information) on top of the reduction operated with the random projections. The specifically designed metric would alleviate the need to come up with the ad-hoc epsilon of this method.

The diffusion mapping has been taken to another level by Mauro Maggioni with diffusion wavelets and Sridhar Mahadevan with Value Function Approximation with Diffusion Wavelets and Laplacian Eigenfunctions.


Reference: [1] S. Lafon, “Diffusion maps and geometric harmonics”, Ph.D. dissertation, Yale University (May 2004)