\documentclass[12pt]{article}
\usepackage{amsmath,amsthm,amsfonts,amssymb}
\usepackage{euscript,graphicx,subfigure,times,version}
\include{amsmath_def}
\include{symbols}
\excludeversion{hide}
\parindent=0in

\title{MAT 280: Laplacian Eigenfunctions: Theory, Applications, and Computations\\
Lecture 15: Use of Laplacian Eigenfunctions and Eigenvalues for
Analyzing Data on a Domain of Complicated Shape}
\author{Lecturer: Naoki Saito\\
        Scribe: Lawrence Austria/Allen Xue}
\date{May 23, 2007}

% SECTION 1 -------------------------------------------------------

\begin{document}
\maketitle

This lecture is based on \cite{NS1} and \cite{NS2}.

\section{Motivation}
Consider a bounded domain of general shape $\varOmega \subset\Rf^d$.
The region $\varOmega$ maybe as simple as a rectangular domain in
$\Rf^2$, or as complicated as a map of an island or a retinal
ganglion cell of a mouse; see Figure~\ref{fig:complicated-domain}.

%%% Insert 1-1.eps
%%% Insert 1-2.eps

\begin{figure}[!h]
\begin{center}
\subfigure[]{\includegraphics[width=2in]{jpn5.eps}}
\subfigure[]{\includegraphics[width=2in]{chalupa_page8.eps}}
\caption{Complicated Domains.}\label{fig:complicated-domain}
\end{center}
\end{figure}

\ \\
We wish to analyze the spatial frequency inside of the object
defined in $\varOmega$, at the same time, we need to avoid the Gibbs
phenomenon due to the boundary $\varGamma=\partial\varOmega$.  We
also want to represent the object information efficiently for
analysis, interpretation, and discrimination; this includes
examining fast decaying expansion coefficients relative to a
meaningful basis. Further, we want to extract geometric information
about the domain $\varOmega$. For more details, see \cite{NS1},
\cite{NS2}.

% SECTION 2 -------------------------------------------------------

\section{Laplacian Eigenfunctions}
Our previous attempt was to extend the object to the outside
smoothly and then bound it nicely with a rectangular box followed by
the ordinary Fourier analysis.  Because various domains give rise to
different sets of Laplacian eigenfunctions, it makes sense to
analyze (and synthesize) the object using genuine basis functions
tailored to the domain.  For instance, {\em sines} (and {\em
cosines}) are eigenfunctions of the Laplacian on the rectangular
domain with Dirichlet (and Neumann) boundary condition.  Also {\em
Spherical harmonics, Bessel functions}, and {\em Prolate Spheroidal
wave functions}, are part of the eigenfunctions of the Laplacian
(via separation of variables) for the {\em spherical, cylindrical},
and {\em spheroidal} domains, respectively.

% Section 2.1 ---

\subsection{Difficulties}
Consider  the operator $\mathcal{L} =-\varDelta$ in $L^2(\varOmega)$
with an appropriate boundary condition.  In the case when $\Omega$
is simple, e.g., a bounded interval $(a,b) \subset \Rf$, or a
rectangular region in $\Rf^d$, analysis of $\mathcal{L}$ is
relatively straightforward.  In general, however, analysis of
$\mathcal{L}$ is difficult mainly due to its unboundedness.  On the
other hand, it is much better to analyze its inverse, i.e., the
Green's operator because it is {\bf compact} and {\bf self-adjoint}.
Thus $\mathcal{L}^{-1}$ has discrete spectra (i.e., countable number
of eigenvalues with finite multiplicity) except for the 0 spectrum.
We use the fact that $\mathcal{L}$ induces a complete orthonormal
basis for $L^2(\varOmega)$ to allow us to perform {\bf eigenfunction
expansion} in $L^2(\varOmega)$.  Another difficulty is computing
such eigenfunctions; directly solving the Helmholtz equation (or the
Laplacian eigenvalue problem) on a general domain $\varOmega$ is
tough. In addition, computing the Green's function for a general
$\varOmega$ satisfying the usual boundary conditions (e.g.,
Dirichlet or Neumann) is very difficult too.

% SECTION 3 ---------------------------------------------------------

\section{Integral Operators Commuting with the Laplacian}
Instead of computing the eigenfunctions of $\mathcal{L}$ on a general
domain, we look at certain integral operators {\bf commuting} with
$\mathcal{L}$.  The key idea is to find an integral operator commuting
with the Laplacian without imposing the strict boundary condition a
priori.  Then we know that the eigenfunctions of the Laplacian is the
same as those of the integral operator, which is much easier to deal
with--thanks to the following fact:

\begin{thm}[G. Frobenius 1878?; B. Friedman 1956]
Suppose $\mathcal{K}$ and $\mathcal{L}$ commute and one of them has an
eigenvalue with finite multiplicity.  Then $\mathcal{K}$ and
$\mathcal{L}$ share the same eigenfunction corresponding to that
eigenvalue.  That is, $\mathcal{L}\varphi = \lambda\varphi$ and
$\mathcal{K}\varphi = \mu \varphi$.
\end{thm}

\ \\
Now, let's replace the Green's function $G(\bx,\by)$ by the {\bf
fundamental solution of the Laplacian}:

\begin{equation}\nonumber
K(\bx,\by)= \left\{ \begin{array}{lcl} -\frac{1}{2}|x-y|
  & \mbox{for} & d=1 \\[8pt] -\frac{1}{2\pi}\log|\bx-\by| &
  \mbox{for} & d=2 \\[8pt]
  \frac{1}{(d-2)\omega_d}|\bx-\by|^{2-d} & \mbox{for} &
  d>2 \end{array}\right.
\end{equation}

where $\omega_d$ is the surface area of the $d$-dimensional unit
ball. The price we pay is to have a rather implicit, {\bf non-local}
boundary condition (although we do not need to deal with this
condition directly).  Let $\mathcal{K}$ be the integral operator
whose kernel is $K$, i.e., for all $f \in L^2(\Omega)$, define:

\begin{equation*}
\mathcal{K}(f(\bx)) \triangleq \int_{\Omega}
K(\bx,\by)f(\by)\dd{\by}.
\end{equation*}

On the other hand, what we gain are the following:

\begin{thm}
The integral operator $\mathcal{K}$ commutes with the Laplacian
$\mathcal{L} = -\varDelta$ with the following non-local boundary
condition:

\begin{equation*}
\Int{\varGamma}{}K(\bx,\by)\normalderivs{\varphi}{\by}(\by)
\dd{s}(\by)=-\frac{1}{2}\varphi(\bx)+ \mathrm{pv}\Int{\varGamma}{}
\normalderivs{K(\bx,\by)}{\by}\varphi(\by)\dd{s}(\by)
\end{equation*}

for all $\bx \in \varGamma = \partial\varOmega$, where $\varphi$ is
an eigenfunction common for both operators.
\end{thm}

\begin{cor}
The integral operator $\mathcal{K}$ is compact and self-adjoint in
$L^2(\varOmega)$.  Thus the kernel $K(\bx,\by)$ had the following
eigenfunction expansion (in the sense of mean convergence):

\begin{equation*}
K(\bx,\by)\sim \sum_{j=1}^{\infty}
\mu_j\varphi_j(\bx)\overline{\varphi_j(\by)}
\end{equation*}

and $\{\varphi_j\}_{j=1}^{\infty}$ forms an orthonormal basis of
$L^2(\varOmega)$
\end{cor}

% Section 3.1 ---

\subsection{Examples}

\begin{ex}
Consider the simplest example, when $\varOmega = (0,1) \subset
\mathbb{R}$.  Then our integral operator $\mathcal{K}$ with the
kernel $K(x,y) = -\frac{1}{2}|x-y|$ gives rise to the eigenvalue
problem:

\begin{displaymath}
  \left\{ \begin{array}{ll}
    -\varphi'' = \lambda \varphi \text{ ,  for } x \in (0,1) &
    \text{(DE)} \\
    \varphi(0)+\varphi(1)=-\varphi'(0)=\varphi'(1) & \text{(BC).}
\end{array} \right.
\end{displaymath}

%%% Insert Figure 3-1.eps
\begin{figure}[!h]
\begin{center}
\includegraphics[width=4in]{line1dtop5basis.eps}
\caption{First 5 eigenfunctions.}\label{fig:topeig}
\end{center}
\end{figure}

Note that the kernel $K(x,y)$ is of {\bf Toeplitz} form; thus, the
eigenvectors must have even and odd symmetry ~\cite{CANTONI-BUTLER}.
In this case, we have the following explicit solution ( the first 5
eigenfunctions are shown in Figure~\ref{fig:topeig}):

\begin{itemize}
\item $\lambda_0 \approx -5.756915$, which is a solution of the equation:
$\tanh(\frac{\sqrt{-\lambda_0}}{2})=\frac{2}{\sqrt{-\lambda_0}}$.
And the corresponding eigenfunction is
$$
\varphi_0(x)=A_0\cosh(\sqrt{-\lambda_0})\(x-\frac{1}{2}\);
$$
\item $\lambda_{2m-1}=(2m-1)^{2}\pi^{2}$ for all $m \in\N$, and the corresponding eigenfunction is
$$
\varphi_{2m-1}(x)= \sqrt{2}\cos(2m-1)\pi x;
$$
\item $\lambda_{2m}$ for all $m \in\N$ is a solution of the
  equation: $\tan
  \frac{\sqrt{\lambda_{2m}}}{2}=-\frac{2}{\sqrt{\lambda_{2m}}}$, and the corresponding eigenfunction is
$$
\varphi_{2m}(x)=A_{2m} \cos\sqrt{\lambda_{2m}}\(x-\frac{1}{2}\).
$$
The constants  $A_{2m}$, $m=0,1,2,\ldots$ are normalization
constants to have $\Vert\varphi_{2m}\Vert=1$.
\end{itemize}
\end{ex}

\ \\
\begin{ex}
Now let $\varOmega = B(0,1) \subset \Rf^2$ be the unit disc. Our
operator $\mathcal{K}$ with the kernel
$K(\bx,\by)=-\frac{1}{2\pi}\log|\bx-\by|$, gives rise to the
eigenvalue problem:

\begin{displaymath}
  \left\{ \begin{array}{ll}
    -\varDelta \varphi = \lambda \varphi \text{ ,  in $\varOmega$ } &
    \text{(DE)}  \\
    \frac{\partial\varphi}{\partial\nu}\Big\vert_{\Gamma}=\frac{\partial \varphi}{\partial r} \Big\vert_{\Gamma} = -\frac{\partial\mathcal{H}\varphi}{\partial\theta}\Big\vert_{\Gamma}, & \text{(BC)} \end{array} \right.
\end{displaymath}
where $\mathcal{H}$ is the Hilbert transform for the circle:

$$
\mathcal{H}f(\theta) \triangleq
\frac{1}{2\pi}\text{pv}\int_{-\pi}^{\pi}
f(\eta)\cot\frac{\theta-\eta}{2}\dd{\eta}
$$
for all $\theta \in [-\pi,\pi]$. Now, let $\beta_{k,\ell}$ be the
$\ell$-th zero of the Bessel function of the first kind of order
$k$, then


\begin{displaymath}
\varphi_{m,n}(r,\theta) = \left\{ \begin{array}{ll}
  J_{m}(\beta_{m-1,n}r)\cos(m\theta) & \mbox{ for all } m,n \in\N   \\
  J_{m}(\beta_{m-1,n}r)\sin(m\theta) & \mbox{ for all } m,n \in\N   \\
  J_0(\beta_{0,n}r) & \mbox { if } m=0, n \in\N
  \end{array}\right.
\end{displaymath}

\begin{displaymath}
\lambda_{m,n} = \left\{ \begin{array}{ll}
  \beta_{m-1,n}^{2} & \mbox{ if  } m,n \in\N  \\
  \beta_{0,n}^{2} & \mbox { if } m=0, n \in\N
\end{array}\right.
\end{displaymath}

\end{ex}

%%% Insert Figure 3-4.eps
%%% Insert Figure 3-5.eps
\begin{figure}[!h]
\begin{center}
\subfigure[]{\includegraphics[width=2in]{circbasis_top25.eps}}
\subfigure[]{\includegraphics[width=2in]{circtrefethen.eps}}
\caption{Left: First 25 eigenfunctions of Our integral operator
$\mathcal{K}$; Right: First 25 eigenfunctions of the Dichlet-Laplace
via separation of variables.} \label{fig:compare-eign}
\end{center}
\end{figure}

\begin{ex}
Consider the unit ball $\varOmega = B(0,1) \subset \Rf^3$. Then our
integral operator uses the kernel:
\begin{equation*}
K(\bx,\by)= \frac{1}{4\pi |\bx-\by|}.
\end{equation*}

The first 9 eigenfunctions cut at the equator viewed from the south
is depicted in Figure~\ref{fig:hemispehere}.

%%% Insert Figure 3-6.eps
\begin{figure}[h]\centering
  \includegraphics[width=4in]{northernhemispeheretop9_noaxes.eps}
  \caption{First 9 eigenfunctions cut at the equator viewed from the south.}
  \label{fig:hemispehere}
\end{figure}
\end{ex}

% SECTION 4 ----------------------------------------------------------

\section{Discretization of the Problem}

Assume that the whole data set consists of a collection of data
sampled on a rectangular grid, and that each sampling cell is a box
of size $\prod_{i=1}^{d} \Delta x_i$. Further, assume that an object
of our interest $\varOmega$ consists of a subset of these boxes
whose centers are $\{\bx_i\}_{i=1}^{N}$.  Under these assumptions,
we can approximate the integral eigenvalue problem
$\mathcal{K}\varphi=\mu\varphi$ with a simple quadrature rule with
node-weight pairs $(\bx_j,w_j)$ as follows:

\begin{equation*}
\sum_{j=1}^{N} w_j K(\bx_i,\bx_j) \varphi(\bx_j)=\mu\varphi(\bx_i)
\end{equation*}

\begin{equation*}
w_j=\prod_{k=1}^{d} \Delta x_k
\end{equation*}

for all $1 \leq j \leq N$.  Let $K_{i,j} \triangleq w_j
K(\bx_i,\bx_j)$, $\varphi_i \triangleq \varphi(\bx_i)$, and $\bvphi
\triangleq (\varphi_1, \cdots, \varphi_N)^T \in \mathbb{R}^N$.  The
above equation can be written in a matrix-vector form
$K\bvphi=\mu\bvphi$, where $K=(K_{i,j}) \in \Rf^{N \times N}$. Under
our assumptions, the weight $w_j$ does not depend on $j$; therefore,
$K$ is {\em symmetric}.

% SECTION 5 ----------------------------------------------------------

\section{Applications}

% Section 5.1 ---

\subsection{Image Approximation}

%%% Insert Figure 5-1-1.eps
%%% Insert Figure 5-1-2.eps

\begin{figure}[!h]
\begin{center}
\subfigure[$\chi_{\varOmega}$]{\includegraphics[width=2.25in]{jpn5.eps}}
\subfigure[$\chi_{\varOmega}\cdot\text{Barbara}$]{\includegraphics[width=2.25in]{barbxjpn5.eps}}
\caption{The characteristic function of the Japanese Islands (a) and
the Barbara image overlaid over the islands.} \label{fig:barb}
\end{center}
\end{figure}

In this section, we provide a brief summary of the approximation
capability of the Laplacian eigenfunction for given image data on an
irregularly-shaped domain, and compare the performance with that of
the standard wavelet-based methods.

%%% Insert Figure 5-1-3.eps
%%% Insert Figure 5-1-4.eps
\begin{figure}[!h]
\begin{center}
\subfigure[$100$-term
Approx.]{\includegraphics[width=2.25in]{barbxjpn5_nonlin100.eps}}
\subfigure[Error]{\includegraphics[width=2.25in]{barbxjpn5_nonlin100_err.eps}}
\caption{The $100$-term approximation and the residual error using
the Laplacian eigenfunctions.} \label{fig:barbjpn-error}
\end{center}
\end{figure}

\begin{figure}[!h]
\begin{center}
\subfigure[$100$-term
Approx.]{\includegraphics[width=2.25in]{barbxjpn5_wsym4_recon100.eps}}
\subfigure[Error]{\includegraphics[width=2.25in]{barbxjpn5_wsym4_nonlin100_err.eps}}
\caption{The $100$-term approximation and the residual error using
the 2DWavelets (Symmlet 8).} \label{fig:bj-w2-err}
\end{center}
\end{figure}

%%% Insert Figure 5-1-5.eps
\begin{figure}[!h]
\begin{center}
\subfigure[$100$-term
Approx.]{\includegraphics[width=2.25in]{barbxjpn5_1dwsym4_recon100.eps}}
\subfigure[Error]{\includegraphics[width=2.25in]{barbxjpn5_1dwsym4_nonlin100_err.eps}}
\caption{The $100$-term approximation and the residual error using
the 1DWavelets (Symmlet 8).} \label{fig:bj-w1-err}
\end{center}
\end{figure}

\ \\
Let $\varOmega$ be the digitized image of the islands of Japan. We
define the characteristic function $\chi_\varOmega(\bx)$ to indicate
the space of the island. As an example, we form a function (or data)
on $\Omega$ by multiplying the standard Barbara image with
$\chi_\varOmega$; see Figure~\ref{fig:barb}.

%%% Insert Figure 5-1-6.eps
\begin{figure}[!h]
\begin{center}
\includegraphics[width=4in]{barbxjpn5_lap_vs_1dwsym4_wsym4_wsym4tensor.eps}
\caption{Comparison of the approximation error in $\ell^2$ using
these transforms.} \label{fig:compare-decay}
\end{center}
\end{figure}

\ \\
For this particular example, the number of samples forming the
data on the island is 1625; therefore, computing the Laplacian
eigenfunctions defined on $\varOmega$ involves a kernel matrix of
size 1625 $\times$ 1625.

\ \\
Figure~\ref{fig:barbjpn-error} shows a reconstruction of the data
set using the largest 100 coefficients. Figure~\ref{fig:bj-w2-err}
shows another reconstruction using the largest 100 coefficients
computed via the standard 2D wavelet basis called {\em Symmlet 8}.
Notice that the 100 wavelet coefficients cannot even capture the
boundary of the domain, not to mention the data on the domain.

\ \\
To be fairer, though, we organize the 1625 data points in a
one-dimensional array in a column scanning order.  After applying
the 1D wavelet transform to this 1D array using the Symmlet 8
filter, we recompute the 100-term approximation and its residual
errors, as shown in Figure~\ref{fig:bj-w1-err}. Observe the
stripe-shaped artifacts in the approximation due to the
deconstruction of the 2D spatial coherency of the original data by
arranging them into a 1D array.

\ \\
We conclude this section with the following conjecture, which
essentially says that the approximation performance of our
eigenfunctions for a general domain is the same order as that of the
Fourier cosine series (or DCT) for a rectangular domain:

\begin{conj}[NS 2005]
For $f \in C(\overline{\varOmega})$ with $\nabla f \in
BV(\overline{\varOmega})$ defined on a $C^2$-domain $\varOmega$, the
expansion coefficients $\langle f,\varphi_k \rangle$ with respect to
the Laplacian eigenbasis decay $O(k^{-1})$.  Thus, the $N$-th
approximation error measured in $L^2$-norm should have a decay rate of $O(N^{-1.5})$.
\end{conj}

% Section 5.2 ---

\subsection{Statistical Image Analysis; Comparison with PCA}
See \cite[Sec. 5]{NS2} for the details.

%%% Insert Figure 5-1-6.eps
%\begin{figure}[!h]\centering
%\includegraphics[width=5in]{5-2-1.eps}
%\caption{Comparison with PCA: Energy Distribution over Coordinates}
%\end{figure}
%%% Insert Figure 5-1-6.eps
%\begin{figure}[!h]\centering
%  \includegraphics[width=3in]{5-2-2.eps}
%  \caption{Comparison with PCA: 13th Basis Vector}
%\end{figure}

% Section 5.3 ---

\subsection{Solving the Heat Equation on a Complicated Domain}
It is well known that the semigroup $\e^{t\varDelta}$ can be
diagonalized using the Laplacian eigenbasis.  More precisely, for
any initial heat distribution $u_0(\bx)\in L^2(\varOmega)$, we have
the heat distribution at time $t$ formally as:

\begin{equation*}
u(\bx, t)=\e^{t\varDelta}u_0=\sum_{j=1}^{\infty}
\e^{-t\lambda_j}\langle u_0,\varphi_j \rangle \varphi_j(\bx)
\end{equation*}
which is based on the expansion of the Green's function for the heat
equation $\rho_t(\bx,\by)$ via the Laplacian eigenfunctions as
follows:

\begin{equation*}
\rho_t(\bx,\by)=\sum_{j=1}^{\infty} \e^{-\lambda_j t} \varphi_j(\bx)
\overline{\varphi_j(\by)},\quad (t,\bx,\by) \in (0,\infty) \times
\overline{\varOmega} \times \overline{\varOmega}.
\end{equation*}

In practice, the domain $\Omega$ is discretized by a finite number
(i.e., $N\in\N$) of sample points (or pixels), thus the Laplacian
eigenfunctions become the Laplacian eigenvectors of length $N$.
Therefore, we can write $\e^{t\varDelta}$ in the matrix-vector
notation as

\begin{equation*}
\varPhi \e^{-t\varLambda} \varPhi^T = \varPhi \text{ diag}(\e^{-\lambda_1
  t},\cdots, \e^{-\lambda_N t})\varPhi^T = \sum_{j=1}^{N}
  \e^{-\lambda_j t} \bvphi_j \bvphi_j^T
\end{equation*}

where $\varPhi = (\bvphi_1, \cdots, \bvphi_N)$ is the Laplacian
eigenbasis matrix of size $N \times N$, and $\varLambda$ is the
diagonal matrix consisting of eigenvalues of the Laplacian, which
are the inverse of the eigenvalues of the kernel matrix
$K=(K_{i,j})$, i.e., $\varLambda_{k,k}=\lambda_k=\frac{1}{\mu_k}$.

\ \\
Given an initial heat distribution over the domain, $\bu_0 \in
\mathbb{R}^N$, we can compute the heat distribution at time $t$ as:
$$
\bu(t)=\varPhi \e^{-t \varLambda}\varPhi^T \bu_0.
$$

Figure~\ref{fig:jpn-heat} shows the result of a simple numerical
experiment on the Japanese islands data sets.

%%% Insert Figure 5-3-1.eps
\begin{figure}[!h]
\begin{center}
\includegraphics[width=5in]{heatjpn5t0_t500.eps}
\caption{The propagation of the heat distribution over the Japanese
Islands. The initial point heat source was put at the location of
Mt. Fuji.} \label{fig:jpn-heat}
\end{center}
\end{figure}


To learn more relationships between heat equation and spectral
geometry, see \cite{DODZIUK}.
% Section 5.4 ---

\subsection{Clustering Mouse Retinal Ganglion Cells}

In this section, we examine how Laplacian eigenfunctions help us
understand how the structural/geometric properties of mouse {\bf
retinal ganglion cells} (RGCs) relate to the cell types and their
functionality. For the background of mouse retinal ganglion cell
datasets, see \cite{COOMBS} and \cite{KONG}. Our data consists of 3D
images of dendrites/axons of RGCs. The process (often arduous)
involves examining each images via specialized software to extract
geometric/morphological parameters (totally 14 parameters) followed
by a conventional (often arduous) clustering algorithm.

\ \\
The parameters include: somal size, dendric field size, total
dendrite length, branch order, and mean internal branch length.
Typically, it takes half a day to process each cell with a lot of
human interaction.

%%% Insert Figure 5-4-1.eps
\begin{figure}[!h]
\begin{center}
\includegraphics[width=3in]{chalupa_page8.eps}
\caption{Mouse Retinal Ganglion Cells (plan and side views).}
\label{fig:RGC}
\end{center}
\end{figure}

% Section 5.4.1 -

\subsubsection{Very Preliminary Studies}
Our preliminary study on mouse RGCs consists of 1) using 2D plane
projection of the data instead of the full 3D, computing the
smallest $k$ Laplacian eigenvalues using our method (i.e., the
largest $k$ eigenvalues of $K$) for each image; 2) constructing a
{\bf feature vector} per image; and 3) performing visualization and
clustering (see Figure~\ref{fig:RGC-features}).  The possible
feature vectors reflecting geometric information are:

$$
\mathbf{F}_1 = (\lambda_1, \cdots, \lambda_k)^T
$$
$$
\mathbf{F}_2 = (\mu_1, \cdots, \mu_k)^T
$$
$$
\mathbf{F}_3 = \(\frac{\lambda_1}{\lambda_2}, \cdots ,
\frac{\lambda_1}{\lambda_k}\)^T
$$
$$
\mathbf{F}_4 = \(\frac{\mu_1}{\mu_2}, \cdots,\frac{\mu_1}{\mu_k}\)^T
$$

\begin{figure}[!h]
\begin{center}
\subfigure[$\bF_1$]{\includegraphics[width=2in]{lapeigvaltop3.eps}}
\subfigure[$\bF_3$]{\includegraphics[width=2in]{lapeiginvfea1.eps}}
\caption{Visualization anc clustering using $\bF_1$ (a) and $\bF_3$
(b).} \label{fig:RGC-features}
\end{center}
\end{figure}

% Section 5.4.2 -

\subsubsection{Challenges of the Mouse Retinal Ganglion Cells Problem}
The shapes of the RGCs are very complicated.  In addition,
interpretation of our eigenvalues are not yet fully understood
compared to the usual Dirichlet-Laplacian case that have been well
studied. Perhaps using the original 3D data may yield a more
meaningful result instead of the projected 2D data.  Another way to
approach this problem is to construct actual graphs based on the
connectivity and analyze them directly via spectral graph theory and
diffusion maps (see Lectures 18-20).  Lastly, because we are dealing
with very complicated domains with many sample points (see
Figure~\ref{fig:challenge}), we need to develop a faster algorithm
to reduce computational burden.

% SECTION 6 ------------------------------------------------------------

\section{A Possible Fast Algorithm For Computing
  $\mathbf{\varphi}_j$'s}
We take advantage of the fact that our kernel function $K(\bx,\by)$
is of special form.  More precisely, it is the fundamental solution
of the Laplacian used in Potential theory. One way to construct a
faster algorithm is to accelerate the matrix-vector product
$K\bvphi$ using the {\bf Fast Multipole Method
(FMM)}~\cite{ROKHLIN-GREENGARD}. We convert the kernel matrix to the
tree-structured matrix via FMM whose submatrices are nicely
organized in terms of their ranks.


%%% Insert Figure 6-1.eps
%%% Insert Figure 6-2.eps

\begin{figure}[h]
\begin{center}
\subfigure[Whole
islands]{\includegraphics[width=2.5in]{Kislands_interact.eps}}
\subfigure[Separated
islands]{\includegraphics[width=2.5in]{Kislands.eps}}
\caption{The
kernel matrix where the communications between and within the
islands are kept (a) and the one without in-between communication
(b).} \label{fig:whole-split}
\end{center}
\end{figure}

Our current implementation cost $O(N\log N)$; Another approach is to
construct an $O(N)$ matrix-vector product module fully utilizing
{\bf rank} information.  For more detailed treatment, see the {\em
HSS} algorithm of Chandrasekaran et al. (2006)~\cite{GUMING-HSS-1}.
We then embed the matrix-vector product module in the Krylov
subspace method, e.g., Lanczos iteration.  As a result, the
computational cost for each eigenvalue and eigenvector is $O(N)$.

\ \\
Lastly, if a domain of data support consists of disconnected
components, then we can reduce the original problem into a set of
smaller problem.  For instance, we can separate the Japanese island
into four major islands, as shown in Figure~\ref{fig:whole-split}.


% SECTION 7 -----------------------------------------------------------

\section{Conclusions}

%%% Insert Figure 7-1.eps
\begin{figure}[!h]
\begin{center}
  \includegraphics[width=3.5in]{geom_map.eps}
  \caption{Real Challenge: The number of sampling points in this map
  is 387,394 over the Japanese Islands.}
  \label{fig:challenge}
  \end{center}
\end{figure}

We demonstrated that our Laplacian eigenfunctions may be useful for
object oriented image analysis and synthesis in which the user can
define the image domain freely and explicitly with the help of
interactive device (e.g., pointer/mouse) or some automatic
segmentation algorithm. We also demonstrated that our method based
on the eigenanalysis of the commuting integral operator leads to
unconventional non-local boundary condition for the Laplacian
eigenvalue problem, but that the numerical implementation is
straightforward and is amenable to modern fast algorithms. Our
experiments and analogy with the analytic examples suggest that we
should be able to get fast-decaying expansion coefficients if the
images are in $C(\overline{\Omega})$ and $\nabla f\in
\text{BV}(\overline{\Omega})$, where the boundary of $\Omega$ is
smooth. In essence, our method can be viewed as a replacement of DCT
for the general shape domain. This means that our eigenfunctions
have a variety of potential applications e.g., interpolation,
extrapolation, local feature computation, and perhaps compression.
In addition, it connects several interesting mathematics, including
spectral geometry, spectral graph theory, isoperimetric
inequalities, Toeplitz operators, PDEs, potential theory, and
almost-periodic functions.

\begin{thebibliography}{10}

\bibitem{NS1}
{\sc N. Saito}, ``Geometric Harmonics as a Statistical Image Processing
Tool for Images Defined on Irregularly Shaped Domains,'' in {\em
  Proc. IEEE Workshop on Statistical Image Processing}, Bordeaux,
France, Jul. 2005

\bibitem{NS2}
{\sc N. Saito}, ``Data Analysis and Representation Using
Eigenfunctions of Laplacians on a General Domain,'' Submitted to
{\em Applied and
  Computational Harmonic Analysis}, Mar. 2007

\bibitem{DODZIUK}
{\sc J. Dodziuk}, "Eigenvalues of the Laplacian and the Heat
Equation," {\em Amer. Math. Monthly}, vol.~88, no.~9, pp.~686--695,
1981.

 \bibitem{COOMBS}
 {\sc J. Coombs, D. van der List, G.-Y. Wang, and L. M. Chalupa}, ``Morphological Properties of Mouse Retinal Ganglion Cells'',  {\em Neuroscience}, vol.~140, no.~1, pp.~123--136, 2006.

\bibitem{KONG}
{\sc J.-H. Kong, D. R. Fish, R. L. Rockhill, and R. H. Masland},
``Diversity of Ganglion Cells in the Mouse Retina: Unsupervised
Morphological Classification and its Limits'', {\em The Journal of
Comparative Neurology}, vol.~489, no.~3, pp.~293--310, 2005.

\bibitem{CANTONI-BUTLER}
{\sc A. Cantoni, P. Butler}, ``Eigenvalues and eigenvectors of
symmetric centrosymmetric matrices'', {\em Linear Algebra Appl.},
vol.~13, pp.~275--288, 1976.

\bibitem{ROKHLIN-GREENGARD}
{\sc L. Greengard, V. Rokhlin}, ``A fast algorithm for particle
simulations'', {\em J. Comput. Phys.}, vol.~73, pp.~325--348, 1987.

\bibitem{GUMING-HSS-1}
{\sc S. Chandrasekaran, M. Gu, and T. Pals}, ``A fast ULV
decomposition solver for hierarchically semi-separable
representations'', {\em SIAM J. Matrix Anal. Appl.}, vol.~28, no.~3, pp.~603--622, 2006.

\end{thebibliography}


\end{document}
