\documentclass[11pt,twocolumn]{article}
\usepackage{times,amsmath,amsfonts}
\usepackage[dvips]{graphicx,graphics}



\setlength{\textheight}{10in} \setlength{\textwidth}{7.1in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.3in} \setlength{\headsep}{0in}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.


% Add the period after section numbers.  Adjust spacing.
%\newcommand{\Section}[1]{\vspace{-5pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
%\newcommand{\SubSection}[1]{\vspace{-2pt}\subsection{\hskip -1em.~~#1}
%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray*}}
\newcommand{\eqn}{\end{eqnarray*}}
\newcommand{\bq}{\begin{eqnarray}}
\newcommand{\eq}{\end{eqnarray}}

\begin{document}

\title{Stat 992: Lecture 17\\
Laplace Operator.}
\author{Moo K. Chung \tt{mchung@stat.wisc.edu}}
\maketitle \thispagestyle{empty}

%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.8]{homework05-1.eps}
%\end{figure}

\begin{enumerate}
\item {\em Problem statement.} In the previous lecture we showed
how to estimate Laplacian in regular grid $\mathbb{Z}^n$. Now we
show how to estimate Laplacian in irregular grid such as polygonal
surfaces in $\mathbb{R}^2$ for the formulation of diffusion
smoothing (Figure 1). The question is how one estimate Laplacian
or any other differential operators in the polygonal surface.
Assume we have observations $Y_i$ at each point $p_i$, which is
assumed to follow additive model $$Y_i = \mu(p_i) + \epsilon(p_i),
p_i \in \mathbb{R}^2$$ where $\mu$ is a smooth continuous function
and $\epsilon$ zero mean Gaussian random fields. We want to
estimate
$$\Delta \mu (p_0)= \frac{\partial^2 \mu}{\partial x^2}\Big|_{p_0} +
\frac{\partial^2 \mu}{\partial y^2}\Big|_{p_0}.$$ Unfortunately,
the geometry of polygonal surfaces forbid direct application of
finite difference scheme. There are a couple of radically
different approaches answering this problem (See one of my
publication that uses finite element method).

\item {\em Local polynomial regression.} Let $p_i = (x_i,y_i)$ be
the coordinate of vertices of polygons. We estimate the Laplacian
at $p_0=(0,0)$ by fitting a quadratic polynomial of the form
 \bq \mu(u,v) =\beta_0 + \beta_1u + \beta_2v + \beta_3u^2 +
\beta_4uv + \beta_5v^2.\label{eq:quadratic}\eq Let
$(u,v)=(x_i,y_i)$ and $\mu(x_i,y_i)=Y_i$ and set up a normal
equation, i.e. \bqn Y_i =\beta_0 + \beta_1x_i + \beta_2y_i +
\beta_3x_i^2 + \beta_4x_iy_i + \beta_5y_i^2.\eqn

Let $Y=(Y_1,\cdots,Y_m)'$, $\beta = (\beta_0,\cdots,\beta_5)'$ and
design matrix
$$\mathbb{X}= \begin{pmatrix}
1& x_1 & y_1 & x_1^2 & x_1y_1 & y_1^2\\
1& x_2 & y_2 & x_2^2 & x_2y_2 & y_2^2\\
  \hdotsfor[2]{6} \\
1&   x_m & y_m & x_m^2 & x_my_m & y_m^2
\end{pmatrix}.$$
Then we have the following matrix equation
$$Y =\mathbb{X}\beta.$$ The unknown coefficients $\beta_i$ are
estimated by the least-squares approximation such that
$$\hat \beta = (\hat\beta_0,\dots,\hat \beta_5)'
=(\mathbb{X}'\mathbb{X})^{-}\mathbb{X}'Y,$$ where $^{-}$ denotes
generalized inverse, which can be obtained through the singular
value decomposition (SVD). Note that $\mathbb{X}'\mathbb{X}$ is
nonsingular if $m > 6$.


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.3]{Figure1.eps}
\caption{Polygonal mesh in $\mathbb{R}^2$.}
\end{figure}


\item {\em Moore-Penrose Generalized inverse.} Generalized inverse
usually used is that of Moore-Penrose. See C.R. Rao and S.K.
Mitra, Generalized inverse of matrices and its applications,
Wiley, New Work, 1971; Penrose, R., A Generalized Inverse for
Matrices. Proc. Cambridge Phil. Soc. 51, 406-413, 1955. It is
usually defined as matrix $\mathbb{X}$ satisfying four conditions
$$\mathbb{X}\mathbb{X}^{-}\mathbb{X}=\mathbb{X}, \; \mathbb{X}^{-}\mathbb{X}\mathbb{X}^{-}
 = \mathbb{X}^{-},$$
$$(\mathbb{X}\mathbb{X}^{-})'=\mathbb{X}\mathbb{X}^{-}, \; (\mathbb{X}^{-}\mathbb{X})'=
\mathbb{X}^{-}\mathbb{X}.$$

Let $\mathbb{X}$ be $m \times p$ matrix with $m \geq p$. Then SVD
of $\mathbb{X}$ is
$$\mathbb{X}=UDV',$$
where $U_{m \times p}$ has orthonormal columns, $D_{p \times
p}=Diag(d_1,\cdots,d_p)$ is diagonal with non-negative elements
and $V_{p \times p}$ is orthogonal. Let
$D^{-}=Diag(d_1^{-},\cdots,d_p^{-}),$ $d_i^{-} = 1/d_i$ if $d_i
\neq 0$ and $d_i^{-} = 0$ if $d_i=0$. Then it can be shown that
the Moore-Penrose generalized inverse is given by
$$\mathbb{X}^{-}=VD^{-}U'.$$

Then the Laplacian of is
$$\widehat {\Delta \mu}(p_0)= 2\hat \beta_3 + 2\hat
\beta_5.$$ Note that $\hat \beta_3$ and $\hat \beta_5$ are
expressed as the weighted averaging of $Y_i$'s. Hence
$$\widehat {\Delta \mu}(p_0) = \sum_{i=1}^m w_i Y_i$$
 where weight $w_i$ are equivalent to the the sum of
the $i$-th component of the 3rd and 5th rows of $\mathbb{X}^{-}$.
Note that the weights $w_i$ are functions of the coordinates of
vertices $p_1,\cdots,p_m$.

{\bf Problem 28.} Express weights $w_i$ in terms of $x_i$'s and
$y_i$'s. {\tt MAPLE} or {\tt MATHEMATICA} may help you. :)




\end{enumerate}

\end{document}