On Nonlinear Dimensionality Reduction, Linear Smoothing and Autoencoding

TL;DR

This paper introduces a new theoretical framework linking nonlinear dimensionality reduction (NLDR) to linear smoothing, enabling better understanding, comparison, and development of NLDR algorithms through differential equations.

Contribution

It establishes a novel connection between NLDR methods and linear smoothers, providing a unified analysis and deriving new algorithms based on differential operators.

Findings

Existing NLDR methods asymptotically approximate solutions to differential equations.

Many NLDR techniques can be characterized by three limiting differential operators.

Local Tangent Space Alignment is shown to be superior among certain methods.

Abstract

We develop theory for nonlinear dimensionality reduction (NLDR). A number of NLDR methods have been developed, but there is limited understanding of how these methods work and the relationships between them. There is limited basis for using existing NLDR theory for deriving new algorithms. We provide a novel framework for analysis of NLDR via a connection to the statistical theory of linear smoothers. This allows us to both understand existing methods and derive new ones. We use this connection to smoothing to show that asymptotically, existing NLDR methods correspond to discrete approximations of the solutions of sets of differential equations given a boundary condition. In particular, we can characterize many existing methods in terms of just three limiting differential operators and boundary conditions. Our theory also provides a way to assert that one method is preferable to another; indeed, we show Local Tangent Space Alignment is superior within a class of methods that assume a global coordinate chart defines an isometric embedding of the manifold.

Tables

Table 1. Table 1: Categorization of NLDR methods by their limit operator and boundary conditions. Within each category, methods with better properties appear closer to the bottom. Newly identified boundary conditions are highlighted in grey. LLE and existing variants have a dependence on the manifold curvature which makes the induced second-order penalty different from 𝔼 ​ T ​ r ​ ( ℋ ​ f ) 2 𝔼 𝑇 𝑟 superscript ℋ 𝑓 2 {\mathbb{E}}\;Tr(\mathcal{H}f)^{2} .

Boundary Condition	2nd order penalty
	$𝔼Tr{\left(\mathscr{H}f\right)}^2$	$𝔼{\Vert \mathscr{H}f\Vert }^2$
None	Coefficient Laplacian
$f=0$	LLR-Laplacian
$\frac{\partial f}{\partial \eta }=0$	Diffusion maps
	Laplacian Eigenmaps
	(LDR-, m-)LLE∗	HLLE
	LDR-LLE+	LTSA
$\frac{{\partial }^{2}f}{\partial {\eta }^2}=\beta \Delta f$	LLR-Laplacian

Table 2. Table 2: Properties of NLDR methods. Non-convergent or unstable methods may converge if regularization is added. Lower order smoothers generally have less variance. New methods and convergence results are highlighted in grey.

Method	$L⪰0$	Order of smoother	Converges & stable
Diffusion Maps	✗	$0^{th}$	✓
Laplacian Eigenmaps	✓	$0^{th}$	✓
LLE	✗	$2^{nd}$	✗
LDR-LLE	✗	$2^{nd}$	✗
LDR-LLE+	✗	$1^{st}$	✓
LLR-Laplacian	✗	$1^{st}$	✓
HLLE	✓	$2^{nd}$	✓
LTSA	✓	$1^{st}$	✓
Coefficient Laplacian	✓	$1^{st}$	✓

Table 3. Table 3: Table of mean Hessian values and the prediction of the 2 n ​ d superscript 2 𝑛 𝑑 2^{nd} derivative given the Laplacian showing fair correspondence of the simulated values to the predicted ones on a modestly fine grid.

$2^{nd}$ Derivative	Estimated value
$\frac{{\partial }^{2}f}{\partial u_1^2}$	0.935
$\frac{{\partial }^{2}f}{\partial u_2^2}$	-0.346
$\frac{{\partial }^{2}f}{\partial u_1\partial u_2}$	0.00
$\Delta f$	0.588
Predicted $\frac{{\partial }^{2}f}{\partial u_1^2}$	0.882

Equations

$$\hat{\mathbf{y}}$$

$${\tilde{X}}_{\mathcal{N}}$$

$${\int }_{{\mathcal{N}}_x}\acute{q}(u)d{V}_g(u)$$

$$(S_{h}f)(x)$$

$$+Tr\left(\mathcal{H}f(x)S_h(y-x)(y-x{)}^T\right)(x)+\cdots$$

$$+O(h^r),$$

$$\Delta u$$

$$X_i$$

$$w_i$$

$$\propto 1_{{\mathcal{N}}_i}-U(I-\lambda (D^2+\lambda I{)}^{-1})U^T{1}_{{\mathcal{N}}_i}.$$

$$\sum _i\Vert \mathcal{H}f(x_i){\Vert }^2$$

$$\underset{\mathbf{A}}{\min }\sum _i{J}_i(\mathbf{Y},\mathbf{A})$$

$$=Y^T\left(\sum _i{Q}_i^{LTSA}\right)Y.$$

$$L_h^{LTSA}f$$

$$\frac{{\partial }^{2}f}{\partial {\eta }_x^2}(x)$$

$$J_{{\Delta }^2}(f)=⟨f,{\Delta }^{2}f⟩=⟨\Delta f,\Delta f⟩=\int Tr(\mathcal{H}f(z){)}^{2}d{V}_g(z)$$

$$J_{HLLE}(f)=\int Tr\left((\mathcal{H}f(z){)}^2\right)d{V}_g(z),$$

$$w_i$$

$$\beta \frac{{\partial }^{2}f}{\partial {\eta }_x^2}(x)$$

$$S_{{\mathcal{N}}_i}$$

$$L$$

$$Q_x^{HLLE}$$

$$={\alpha }_h(Q_x^{LTSA}+(H_x^{(2)}-I_{{\mathcal{N}}_x})).$$

$${\alpha }_h{g}^T{R}_{x}f$$

$$\frac{{\partial }^{2}f}{\partial {\eta }_x^2}(x)$$

$$E(u_i^2|u_1)=(m+1{)}^{-1}(1-u_1^2).$$

$$h^2(L_h{\delta }_x)(y)$$

$$-\sum _{i=2}^m\Vert {\tilde{u}}_i(\cdot ,y){\Vert }^{-2}{\tilde{u}}_i(y,y)⟨{\tilde{u}}_i(\cdot ,y),{\delta }_x⟩$$

$$-\Vert v_1(\cdot ,y){\Vert }^{-2}{v}_1(y,y)⟨v_1(\cdot ,y),{\delta }_x⟩$$

$$={\delta }_x(y)-Vol({\mathcal{N}}_y{)}^{-1}I(y\in {\mathcal{N}}_x)$$

$$-Vol({\mathcal{N}}_y{)}^{-1}\Vert v_1(\cdot ,y){\Vert }^{-2}⟨{\tilde{u}}_1(\cdot ,y),I(\cdot \in {\mathcal{N}}_y)⟩v_1(x,y)$$

$$=-\Vert v_1(\cdot ,y){\Vert }^{-2}({\tilde{u}}_1(x,y){\mu }_y-{\mu }_y^2)$$

$$=\Vert v_1(\cdot ,y){\Vert }^{-2}(Vol({\mathcal{N}}_y{)}^{-1}{h}^{-1}{u}_1(y){\mu }_y+{\mu }_y^2)$$

$$⟨L_h{\delta }_x,I(\cdot \in {\mathcal{N}}_x)⟩$$

$$=\Omega (h^{-2})$$

$$\displaystyle\to\infty\quad\mbox{as $h\to 0$}$$

$$(L_{h}f)(x)$$

Full text

1 Introduction

One of the major open problems in machine learning involves the development of theoretically-sound methods for identifying and exploiting the low-dimensional structures that are often present in high-dimensional data. Such methodology would have applications not only in supervised learning, but also in visualization and nonlinear dimensionality reduction, semi-supervised learning, and manifold regularization.

Some initial steps have been made in this direction over the years under the rubric of “manifold learning methods.” These nonlinear dimension reduction (NLDR) methods have permitted interesting theoretical analysis and allowed the field to move beyond linear dimension reduction. But the theoretical results have fallen short of providing characterizations of the overall scope of the problem, including the similarities and differences of existing methods and their respective advantages. Consider, for example, the classical problem of finding nonlinear embeddings of a Swiss roll with hole, shown in Figure 1. Of the methods shown, Local Tangent Space Alignment is clearly the best at recovering the underlying planar structure of the manifold, but there is no existing theoretical explanation of this fact. Nor are there answers to the natural question of whether there are scenarios in which other methods would be better, why other methods perform worse, or whether the deficiencies can be corrected. The current paper aims to tackle some of these problems. Not only do we provide new characterizations of NLDR methods, but by correcting deficiencies of existing methods we are able to propose new methods in Laplace-Beltrami approximation theory.

We analyze the general class of manifold learning methods that we refer to as local, spectral methods. These methods construct a matrix using only information in local neighborhoods and take a spectral decomposition to find a nonlinear embedding. This framework includes the commonly used methods: Laplacian Eigenmaps (LE) Belkin & Niyogi (2003), Local Linear Embedding (LLE) Roweis & Saul (2000), Hessian LLE (HLLE) Donoho & Grimes (2003), Local Tangent Space Alignment (LTSA) Zhang & Zha (2004), and Diffusion Maps (DM) Coifman & Lafon (2006). We also consider several recent improvements to these classical methods, including Low-Dimensional Representation-LLE Goldberg & Ritov (2008), Modified-LLE Zhang & Wang (2007), and MALLER Cheng & Wu (2013). Outside of our scope are global methods that construct dense matrices that encode global pairwise information; these include multidimensional scaling, principal components analysis, isomap Tenenbaum et al. (2000), and maximum variance unfolding Weinberger & Saul (2006). We note in passing that these global methods have serious practical limitations, either in terms of a strong linearity assumption or computational intractability.

Our general approach proceeds by showing that the embeddings for convergent local methods are solutions of differential equations, with a set of boundary conditions induced by the method. The treatment of the boundary conditions is the critical novelty of our approach. This approach allows us to categorize the methods considered into three classes of differential operators and accompanying boundary conditions, as shown in Table 1. We are also able to delineate properties that allow comparisons of methods within a class; see Table 2. These theoretical results allow us to, for example, conclude that when the goal is to find an isometry-preserving, global coordinate chart, LTSA is best among the classical methods considered. We can also show that HLLE belongs to the same class and converges to the same limit. Hence, they can be used exchangeably as smoothness penalties. We also improve existing Laplace-Beltrami approximations. In particular, we give a Laplace-Beltrami approximation that is consistent as a smoothing penalty even when boundary conditions are not satisfied.

Our analysis is based on the following two insights. First, matrices obtained from local, spectral methods can be seen as an operator that returns the bias from smoothing a function. This allows us to make a connection to the theory of linear smoothing in statistics. Second, as the neighborhood sizes for each local method decrease, we obtain a linear operator that evaluates infinitesimal differences in that neighborhood. In other words, we obtain convergence to a differential operator. Thus, the asymptotic behavior of a procedure can be characterized by the corresponding differential operator in the interior of the manifold and (crucially) by the boundary bias in the linear smoother.

2 Preliminaries

We begin by introducing some basic mathematical concepts and providing a high-level overview of the relationship between NLDR and linear smoothing.

Local, spectral methods share a construction paradigm:

Choose a neighborhood for each point $x_{i}$; 2. 2.

Construct a matrix $L$ with $L_{ij}\neq 0$ only if $x_{j}$ is a neighbor of $x_{i}$; 3. 3.

Obtain a nonlinear embedding from an eigen- or singular value decomposition of $L$.

Step 2 for constructing the matrix $L$ can be seen to be equivalent to constructing a linear smoother. We use this equivalence to compare existing NLDR methods and generate new ones by examining the asymptotic bias of the smoother on the interior and boundary of the manifold. In particular, we show that whenever the neighborhoods shrink in Step 1 as the number of points $n\to\infty$, the matrix $L$ converges to a differential operator on the interior of the manifold. Nonlinear dimensionality reduction methods thus solve a eigenproblem involving this differential operator. The solutions of the eigenproblem on a compact manifold are typically only well-separated and identifiable in the presence of boundary conditions. If $L$ is symmetric and the smoother’s bias has a different asymptotic rate on the boundary than the interior, then it imposes a corresponding boundary condition for the eigenproblem. If $L$ is not symmetric, then the boundary conditions can be determined from the boundary bias for $L$ and its adjoint $L^{T}$.

2.1 Linear smoothers

Consider the standard regression problem for data generated by the model $Y_{i}=f(X_{i})+\epsilon_{i}$ where $\epsilon_{i}$ is a zero-mean error term. A regression model is a linear smoother if

[TABLE]

and the smoothing matrix $S(\mathbf{x})$ does not depend on the outcomes $\mathbf{y}$. Examples of linear smoothers include linear least squares regression on a matrix $X$, kernel and local polynomial regression, and ridge regression.

In the regression problem, one wishes to infer the conditional mean function $f$ given noisy observations $Y$. For manifold learning problems, the goal is to learn noiseless global coordinates $\Phi_{i}=f(X_{i})$ but the outcomes $\Phi$ are never observed. In both linear smoothing and manifold learning, the smoothing matrix makes no use of the response and simply maps rough functions into a space of smoother functions. The embedding that is learned by NLDR methods is a space of functions that are best reconstructed by the smoother. In other words, the smoother defines a form of autoencoder and the learned space is a space of likely outputs from the autoencoder.

2.2 NLDR constructions of linear smoothers

It is a simple but crucial fact that each of the NLDR methods that we consider construct or approximate a matrix $L=G(I-S)$ where $G$ is a diagonal matrix and $S$ is a linear smoother. Thus, $L$ measures the bias of the prediction $Sf$ weighted by $G$.

For example, Diffusion Maps, with a Gaussian kernel constructs the Nadaraya-Watson smoother $S=D^{-1}K$, where $D$ is the diagonal degree matrix and $K$ is the kernel matrix. The constructed embedding consists of the right singular vectors corresponding to the smallest singular values of $L^{DM}=I-S$. Laplacian Eigenmaps, using an unnormalized graph Laplacian, differ only by the reweighting $L^{LE}=DL^{DM}=D(I-S)$, which ensures that the matrix $L^{LE}$ is positive semidefinite.

2.3 Assumptions

The usual setting for linear smoothing is Euclidean space. To account for the manifold setting, we must make some regularity assumptions and demonstrate that calculations in Euclidean space approximate calculations on the manifold sufficiently closely.

We consider a smooth, compact $m$-dimensional Riemannian manifold $\mathcal{M}$ with smooth boundary embedded in $\mathbb{R}^{d}$. From this manifold, points $\{X_{1},\ldots,X_{n}\}\in\mathbb{R}^{d}$ are drawn from a uniform distribution on the manifold. This uniformity assumption can be weakened to admit a smooth density and obtain density weighted asymptotic results.

In the continuous case, we consider neighborhoods $\mathcal{N}_{x}=B(x,h)$ for a given bandwidth $h$ where distance is Euclidean distance in the ambient space. When appropriate, we will take $h\to 0$ as $n\to\infty$. In the discrete case, denote by $\mathcal{N}_{i}$ the indices for neighbors of $X_{i}$ in the point cloud. This neighborhood construction is only mildly restrictive since kNN neighborhoods are asymptotically equivalent to $\epsilon$ neighborhoods in the interior of the manifold when points are sampled uniformly and the radius of the neighborhoods is $\Theta(h)$ on the boundary.

2.4 Local coordinates

Most existing methods and the ones we construct require access to a local coordinate system with manifold dimension $m$ for each point. This coordinate system is estimated using a local PCA or SVD at each point. Let $X_{\mathcal{N}}$ be the $|\mathcal{N}|\times d$ matrix of ambient space coordinates for points in neighborhood $\mathcal{N}$ centered on $x$ and $\tilde{X}_{\mathcal{N}}$ be the corresponding matrix centered on $x$. The $m$ top right singular vectors

[TABLE]

give estimated tangent vectors at $x$. The top rescaled left singular vectors, $\tau_{i}=\lambda_{i}U_{\mathcal{N},i}$, project points in $\mathcal{N}_{i}$ to the tangent space. The normal coordinates $u_{i}$, corresponding to the geodesics traced out by the tangent vectors $V_{\mathcal{N},i}$, agree closely with these tangent space coordinates. Specifically, by Lemma 6 in Coifman & Lafon (2006), $u_{i}=\tau_{i}+q_{3}(\mathbf{t})+O(h^{4})$, where $q_{3}$ is a homogeneous polynomial of degree three whenever the coordinates are in a unit ball of radius $h$. We adopt the convention that $u_{i}(y)$ refers to the normal coordinates at $x\in\mathcal{M}$ for $y\in\mathcal{N}_{x}$.

Our results rely on integration in normal coordinates for a ball of radius $h$ in the distance metric of the ambient space. To account for the manifold curvature, the volume form and neighborhood sizes must be accounted for in the integral. Lemma 7 in Coifman & Lafon (2006) further provides a Taylor expansion for the volume form for the Riemannian metric $dV_{g}(u)=1+q_{2}(u)+O(h^{3})$, where $q_{2}$ is a homogeneous polynomial of degree two. Likewise, distances in the ambient space and normal coordinates differ by $\|y-x\|^{2}=\|u(y)\|^{2}+\tilde{q}_{2}(u(y))+O(h^{3})$ where $u(y)$ denotes normal coordinates for $y$ about $x$, $\tilde{q}_{2}$ is homogenous degree 2, and the distance on the left-hand side is with respect to the ambient space. Consequently, integrals for any homogeneous polynomial $\acute{q}$ of degree one or two satisfy

[TABLE]

where the integral on the right represents integration in $\mathbb{R}^{m}$ and the $O(h^{4})$ term hides a smooth function that depends on the curvature of the manifold at $x$.

For the purposes of this paper, we do not account for error from estimation of the tangent space. Thus we assume that we have a sufficiently accurate estimate of the normal coordinates. Accounting for this estimation error is a natural direction for further work.

3 Analysis of existing methods

Each of the existing NLDR methods construct a matrix $L_{n}$ from a point cloud of $n$ points with bandwidth $h$. The resulting nonlinear embeddings are obtained from the bottom eigenvectors of $L_{n}$ or $L_{n}^{T}L_{n}$. We consider $L_{n}$ as a discrete approximation to an operator $\mathbf{L}_{h}:\mathcal{F}\to\mathcal{F}$ on functions $\mathcal{F}\subset C^{\infty}(\mathcal{M})$. We examine the limit operator $\mathbf{L}_{h}$ constructed by each of the NLDR methods for a fixed bandwidth $h$. We show that as the bandwidth $h\to 0$, $\mathbf{L}_{h}\to\mathbf{L}_{0}$ where $\mathbf{L}_{0}$ is a differential operator. The exception is LLE where there is no well defined limit. The stochastic convergence of the empirical constructions $L_{n}\stackrel{{\scriptstyle n\to\infty}}{{\to}}L_{h}$ is not considered in this paper.

3.1 Taylor expansions and local polynomial bases

As described in Section 2.2, most existing NLDR methods can be expressed as $L_{h}=D_{h}(I-S_{h})$, where $D_{h}$ is a multiplication operator corresponding to a diagonal matrix and $S_{h}$ is a linear smoother. Denote by $S_{h}(x,\cdot)$ the linear function such that $\langle S_{h}(x,\cdot),f\rangle=(S_{h}f)(x)$. This function exists by the Riesz representation theorem. We say $S_{h}$ is local when for all $x\in\mathcal{M}$, the support of $S_{h}(x,\cdot)$ is contained in a ball of radius $h$ centered on $x$. We assume that $S_{h}$ is a bounded operator (as unbounded operators are necessarily poor smoothers).

Consider $S_{h}$ applied to a function $f\in C^{\infty}(\mathcal{M})$. A Taylor expansion of $f$ in a normal coordinates at $x$ gives

[TABLE]

where the error term holds since $f$ has bounded $r^{th}$ derivatives due to the compactness of $\mathcal{M}$. Here $y-x=u$ is a function denoting the normal coordinate map at $x$. From this it is clear that the asymptotic behavior of $S_{h}$ in a neighborhood of $x$ is determined by its behavior on a basis of local polynomials. Furthermore, it can well approximated by examining the behavior on only low-degree polynomials.

3.2 Convergence to a differential operator

Of particular interest is the case where $S_{n}$ operates locally on shrinking neighborhoods. In this case, the following theorem (proven in the supplementary material) is an immediate consequence of applying the sequence of smoothers to the Taylor expansion.

Theorem 1**.**

Let $S_{n}$ be a sequence of linear smoothers where the support of $S_{n}(x,\cdot)$ is contained in a ball of radius $h_{n}$. Further assume that the bias in the residuals $h_{n}^{-k}(I-S_{n})f=o(1)$ for some integer $k$ and all $f\in\mathcal{C}\subset C^{\infty}(\mathcal{M})$. Then if $h_{n}\to 0$ and $h_{n}^{-k}(I-S_{n})$ converges to a bounded linear operator as $n\to\infty$, then it must converge to a differential operator of order at most $k$ on the domain $\mathcal{C}$.

As a result, each of the existing NLDR methods considered can be seen as discrete approximations for the solutions of differential equations $Df=\lambda f$ given some boundary conditions, where $D$ is some differential operator and $\lambda\in\mathbb{R}$.

3.3 Laplacian and boundary behavior

The convergence of Laplacians constructed from point clouds to a weighted Laplace-Beltrami operator has been well-studied by Hein et al. (2007) and Ting et al. (2010). In particular, the spectral convergence of Laplacians constructed using kernels has been studied by Belkin & Niyogi (2007), von Luxburg et al. (2008), and Berry & Sauer (2016). In the presence of a smooth boundary, the eigenproblem for the Diffusion Maps Laplacian that converges to the unweighted Laplace-Beltrami operator has been shown by Coifman & Lafon (2006) and Singer & Wu (2016) to impose Neumann boundary conditions. Specifically, it solves the differential equation

[TABLE]

where $\eta_{x}$ is the vector normal to the boundary at $x$. This result is easily extended to other Laplacian constructions.

As the Diffusion Maps operator is the bias operator for a Nadaraya-Watson kernel smoother, and the unnormalized Laplacian is a rescaling of this bias, the Neumann boundary conditions can be derived from existing analyses for kernel smoothers. The boundary bias of the kernel smoother is $c\cdot h\cdot\partial f/\partial\eta_{x}(x)+o(h)$ for some constant $c$ when $x\in\partial\mathcal{M}$ while in the interior it is of order $O(h^{2})$ when points are sampled uniformly. This trivially follows from noting that the first moment of a symmetric kernel is zero when the underlying distribution is uniform and that the boundary introduces asymmetry in the direction $\eta_{x}$. Hence, the operator must be scaled by $\Theta(h^{-2})$ to obtain a non-trivial limit in the interior, but this scaling yields divergent behavior at the boundary when the Neumann boundary condition $\partial f/\partial\eta_{x}(x)=0$ is not met. Since eigenfunctions cannot exhibit this divergent behavior, they must satisfy the boundary conditions. The significance of this result is given by the following simple theorem.

Theorem 2**.**

Let $L$ be an operator imposing the Neumann boundary condition $\partial f/\partial\eta_{x}(x)=0$. Then, if there exists an isometry-preserving global coordinate chart, a spectral decomposition of $L$ cannot recover the global coordinates, as they are not in the span of the eigenfunctions of $L$.

3.4 Local Linear Embedding

Local Linear Embedding (LLE) solves the boundary bias problem by explicitly canceling out the first two terms in the Taylor expansion. For a point cloud $\mathbf{X}$, LLE constructs a weight matrix $W$ satisfying $W1=1$ and $W\mathbf{X}=\mathbf{X}$. Goldberg & Ritov (2008) and Ting et al. (2010) showed that this condition is not sufficient to identify a limit operator. The behavior of LLE depends on a regularization term. We give a more refined analysis of the constrained ridge regression procedure used to generate the weights.

Let $\tilde{\mathcal{N}}_{i}=\mathcal{N}_{i}\backslash\{i\}$ be the neighborhood that excludes $i$ and $\tilde{X}_{\mathcal{N}_{i}}=X_{\mathcal{N}_{i}}-1_{|\mathcal{N}_{i}|}X_{i}$ be the points in $\tilde{\mathcal{N}}_{i}$ centered on $X_{i}$. To reconstruct $X_{i}$ from its neighbors, one solves the regression problem

[TABLE]

under the constraint that the weights sum to one, $\|w_{i}\|_{1}=1$. Adding a ridge regression penalty $\lambda$ and applying Lagrange multipliers gives $w_{i}\propto(\tilde{X}_{\tilde{\mathcal{N}}_{i}}\tilde{X}_{\tilde{\mathcal{N}}_{i}}^{T}+\lambda I)^{-1}1_{\mathcal{N}_{i}}$.

Consider the singular value decomposition $\tilde{X}_{\tilde{\mathcal{N}}_{i}}=UDV^{T}$, and let $U_{c}$ be the orthogonal complement of $U$. This allows us to rewrite the weights

[TABLE]

In other words, the weights are the same as the constant weights used by a $h$-ball Laplacian but with a correction term. This correction term in Eq. 8 depends only on the curvature of the manifold. One can see this by noting that the top $m$ left singular vectors correspond to the points expressed in normal coordinates around $X_{i}$. As the neighborhood is symmetric around $X_{i}$, it follows that $U_{1:m}^{T}1\approx 0$. The remaining left singular vectors correspond to directions of curvature of the manifold. From this, is it easy to see that LLE with $h$-ball neighborhoods is equivalent to $\Delta_{h}-C$ in the interior of the manifold where $\Delta_{h}$ is the $h$-ball Laplacian and $C$ is some correction term that reduces the penalty on the curvature of a function when it matches the curvature of the manifold.

3.5 Hessian LLE

Hessian LLE performs a local quadratic regression in each neighborhood. It then discards the locally linear portion in order to construct a functional that penalizes second derivatives. The Taylor expansion in Eq. 3.1 gives that the entries of the Hessian are coefficients of a local quadratic polynomial regression.

Given a basis $Z_{\mathcal{N}_{i}}$ of monomials of degree $\leq 2$ in local coordinates for each point in neighborhood $\mathcal{N}_{i}$ where the first $m+1$ basis elements span linear and constant functions, the entries in the estimated Hessian $\mathcal{H}(f)_{X_{i}}$ are $\Pi(Z_{\mathcal{N}_{i}}^{T}Z_{\mathcal{N}_{i}})^{-1}Z_{\mathcal{N}_{i}}f$ where $\Pi$ is the projection onto the $m+2$ through $1+m+m(m-1)/2$ coordinates.

Rather than directly estimating the Hessian, HLLE performs the following approximation. Orthonormalize the basis of monomials using Gram-Schmidt and discard the first $m+1$ vectors to obtain $\tilde{Z}_{\mathcal{N}_{i}}$. By orthonormality, the regression coefficients in this basis are simply the inner products $\tilde{Z}_{\mathcal{N}_{i}}^{T}f$. We note that the approximation only recovers the Hessian in the interior of the manifold and when the underlying sampling distribution is uniform.

Given Hessian estimates $\hat{\beta}_{i}=\tilde{Z}_{\mathcal{N}_{i}}f$ for each point, the sum of their Frobenius norms is easily expressed as

[TABLE]

3.6 Local Tangent Space Alignment

Local tangent space alignment (LTSA) computes a set of global coordinates that can be aligned to local normal coordinates via local affine transformations.

Given neighbors of $X_{i}$ expressed in local coordinates $U_{\mathcal{N}_{i}}\in\mathbb{R}^{|\mathcal{N}_{i}|\times m}$, LTSA finds global coordinates $\mathbf{Y}$ and affine alignment transformations $A_{i}$ that minimize the difference between the aligned global and local coordinates, $J_{i}(\mathbf{Y},\mathbf{A})=\|CY_{\mathcal{N}_{i}}-U_{\mathcal{N}_{i}}A_{i}\|_{F}^{2}$, where $C=1-11^{T}/|\mathcal{N}_{i}|$ is a centering matrix.

Given the global coordinates $Y$, this is a least-squares regression problem for a fixed set of covariates $U_{\mathcal{N}_{i}}$. Thus, the best linear predictors for $CY_{\mathcal{N}_{i}}$ are given by $CY_{\mathcal{N}_{i}}H^{(1)}_{\mathcal{N}_{i}}$ where $H^{(1)}_{\mathcal{N}_{i}}$ is the hat matrix for the covariates $U_{\mathcal{N}_{i}}$. The objective can be expressed using local operators $Q^{LTSA}_{i}$ as

[TABLE]

4 Equivalence of HLLE and LTSA

Although LTSA is derived from a very different objective than HLLE, we will show that they are asymptotically equivalent. This greatly strengthens a result in Zhang et al. (2018) which showed the equivalence of a modified version of HLLE to LTSA under the restrictive condition that there are exactly $m(m+1)/2+m-1$ neighbors.

We show the asymptotic equivalence in two steps. First, we show $L_{h}^{HLLE}$ and $L_{h}^{LTSA}$ converge to the same differential operator by showing convergence in the weak operator topology. Then, we give an argument that derives the boundary condition for both methods.

4.1 Continuous extension

Both HLLE and LTSA are composed of the sum of local projection operators $Q_{i}$ on the point cloud. As we wish to study their limit behavior on smooth functions, we consider the continuous extension to the manifold. To form the continuous extension, the sum is replaced by an integral

[TABLE]

where $V_{g}$ is the natural volume measure on $\mathcal{M}$, $\alpha_{h}$ is an appropriate normalizing constant, $Vol_{h}$ is the volume of a ball of radius $h$ in $\mathbb{R}^{m}$, and $Q_{y}$ are local operators. In the case of LTSA, $Q^{LTSA}_{y}=I_{\mathcal{N}_{y}}-H_{y}$ where $H_{y}$ is the projection onto linear functions in the neighborhood $\mathcal{N}_{y}$ and on a local normal coordinate system rather than linear functions on a discrete set of points.

To identify the order of the differential operator that LTSA converges to, we simply need to find the scaling that yields a non-trivial limit. For the functional $f^{T}L_{h}^{LTSA}f$ to non-trivially converge as $h\to 0$, the component terms $f^{T}Q_{y}f$ should also non-trivially converge. Since $Q_{y}$ is a projection onto the space orthogonal to linear functions, this is equivalent to $f^{T}Q_{y}^{T}Q_{y}f=\|Q_{y}f\|^{2}=O(\alpha_{h}h^{4}(\|\mathcal{H}f_{y}\|^{2}+h))$. For this to converge non-trivially, one must have $\alpha_{h}=O(h^{-4})$. The same argument holds for HLLE. Thus, both LTSA and HLLE yield fourth-order differential operators in the interior of $\mathcal{M}$.

4.2 HLLE and LTSA difference is neglible

To show the equivalence of HLLE and LTSA, we show that the remainder term resulting from their difference is the bias term of an even higher order smoother. As such, this bias term goes to zero faster than the HLLE and LTSA bias terms, which yields convergence in the weak operator topology. This gives following theorem. Proof details are given in the supplementary material.

Theorem 3**.**

$L_{h}^{LTSA}-L_{h}^{HLLE}\to 0$* as $h\to 0$ in the weak operator topology of $C^{\infty}(\mathcal{M})$ equipped with the $L_{2}$ norm.*

4.3 HLLE and LTSA boundary behavior

To establish that the limit operator has the same eigenfunctions, we must show that the boundary conditions that those eigenfunctions satisfy also match. We provide a theorem where a single proof identifies the boundary condition for both methods. This boundary condition applies to the second derivatives of a function and admits linear functions, unlike the Neumann boundary condition imposed by graph-Laplacian-based methods.

Theorem 4**.**

$L^{LTSA}_{h}f(x)\to\infty$* and $L^{HLLE}_{h}f(x)\to\infty$ as $h\to 0$ for any $x\in\partial\mathcal{M}$ and $f\in C^{\infty}(\mathcal{M})$ unless*

[TABLE]

where $\eta_{x}$ is the tangent vector orthogonal to the boundary.

We outline the proof, providing the details in the supplementary material. The proof shows that boundary bias of the LTSA smoother is of order $\Omega(h^{2}Tr(\mathcal{H}f(x)\Sigma))$ for some matrix $\Sigma$. This is of lower order than the required scaling $h^{-4}$ for the functional to converge.

Since locally linear and constant functions are in the null space by construction, only the behavior of the smoother on quadratic polynomials needs to be considered. The difficulty in analyzing the behavior arises from needing to average over multiple local regressions on different neighborhoods and dealing with a non-product measure. By exploiting symmetry, the problem of a non-product measure is removed, and the multivariate regression problem can be reduced to a univariate linear regression with respect to $u_{1}$, where $u_{1}$ is the coordinate function in the direction $\eta_{x}$ orthogonal to the boundary. Due to the shape of the asymmetric neighborhood, the quadratic functions $u_{i}^{2}$, $i>1$, induce a nonzero conditional mean $E(u_{i}^{2}|u_{1})=(1-u_{1}^{2})/(m+1)$. Since $u_{1}^{2}$ has a negative sign in this expression, the coefficients for $u_{1}^{2}$ and each $u_{i}^{2}$ can be set to cancel out the boundary bias to obtain the result.

4.4 Partial equivalence of HLLE and Laplacian

Although the Laplace-Beltrami operator induces a penalty on the gradient and appears fundamentally different from the Hessian penalty induced by HLLE and LTSA, they can be compared by squaring the Laplace-Beltrami operator to create a fourth-order differential operator that induces a slightly different penalty on the Hessian.

Consider the smoothness penalty obtained by squaring the Laplace-Beltrami operator compared to the HLLE penalty.

[TABLE]

where the second equality follows from the self-adjointness of the Laplace-Beltrami operator. These penalties are identical for one-dimensional manifolds. However, a slightly different penalty is induced on multidimensional manifolds. Given global coordinates $g_{1},\ldots,g_{m}$, the quadratic polynomial $g_{1}g_{2}$ is in the null space of the twice-iterated Laplacian but not in the null space of the HLLE penalty.

5 Alternate constructions

We provide a few examples to illustrate how alternative constructions can address undesirable properties of existing methods and the implications of the changes. In particular, we propose a convergent and more stable variation of LLE. We find LLE and its variants behave similarly to a local linear regression smoother, replacing the Neumann boundary condition for Laplacians with a second-order boundary condition that admits linear functions. Furthermore, we generate a new Laplace-Beltrami approximation that generates a smoothing penalty that does not impose boundary conditions. The resulting operator is shown to not have well-separated eigenfunctions.

5.1 Low-dimensional representation LLE+

LLE has the deficiency that the curvature of the manifold affects the resulting operator and smoothness penalty. In the worst case when the regularization parameter is very small, the embedding $\Phi:\mathcal{M}\to\mathbb{R}^{d}$ of the manifold in the ambient space lies close to the null space of the resulting LLE operator. In other words, a spectral decomposition of the LLE operator may simply recover a linear transformation of the ambient space.

Low dimensional representation Goldberg & Ritov (2008) LLE (LDR-LLE) is a modified version of LLE that removes some of the effect of curvature in the manifold by reconstructing each point based on its tangent space coordinates only. As this still cancels out the first two terms of a Taylor expansion, it is an approximation to the Laplace-Beltrami operator. The weights at $X_{i}$ are chosen to be in the subspace spanned by points in $\mathcal{N}_{i}$ in the ambient space and orthogonal to the tangent space, and thus, explicitly penalizes functions with curvature that matches the manifold. However, as the singular values for directions orthogonal to the tangent space decrease as $O(h^{2})$ compared to $O(h)$ for the tangent space, the method is subject to numerical instability. As such we classify LDR-LLE as well as LLE as being based on a second-order smoother. To compensate for this instability, it adds a regularization term.

We propose a further modification, LDR-LLE+, which removes the artificial restriction on the weights to the span of the points in the neighborhood and generates a first-order smoother. The weights are given by

[TABLE]

where $U_{1:m}$ are the top left singular vectors of $X_{\mathcal{N}_{i}}$ after centering on $X_{i}$.

Another variant of LLE is modified LLE (mLLE) Zhang & Wang (2007). While LDR-LLE explicity constructs a local smoother that can be expressed by linear functions in the tangent coordinates and a bias operator in which linear functions are in its null space, mLLE achieves a similar effect by adding vectors orthogonal to the tangent space to the bias operator in order to remove non-tangent directions from the null space.

5.2 Local Linear Regression

A simple alternative construction is to directly use a local linear regression (LLR) smoother to approximate the Laplace-Beltrami operator. Since only functions linear in the normal coordinates should be included in the null space of the bias operator, each local regression is performed on the projection of the neighborhood into its tangent space. Cheng & Wu (2013) propose this linear smoother in the context of local linear regression on a manifold and note its connection to the Laplace-Beltrami operator without analyzing the convergence of the associated operator. Our contribution is determining the boundary conditions.

In particular, we study the boundary behavior for the operator $L^{LLR*}L^{LLR}$ where $L^{LLR*}$ is the adjoint of $L^{LLR}$. The resulting boundary conditions have the same form as those for LTSA, albeit with different constants, and admit linear functions as eigenfunctions.

Theorem 5**.**

Let $S_{h}$ be the expected continuous extension for a local linear regression smoother with bandwidth $h$ and $L_{h}=h^{-2}Vol_{h}^{-1}(I-S_{h})$. Then there is some $\beta>0$ such that for any $f\in C^{\infty}(\mathcal{M})$ and $x\in\partial\mathcal{M}$, $L_{h}^{*}L_{h}f(x)$ converges as $h\to 0$ only if

[TABLE]

where $\eta_{x}$ is the vector normal to the boundary at $x$.

This result is obtained by examining the behavior of the adjoint $L_{h}^{*}$ by reducing it to computing the average residual over a set of univariate linear regressions. The result is that the adjoint has boundary bias $\Omega(g(x))$ when applied to a function $g$. Combined with existing analyses for the boundary bias of $L^{(LLR)}$ Ruppert & Wand (1994) gives $(L_{h}^{*}L_{h}f)(x)=\Omega(h^{-2})\Omega(1)\to\infty$ if the boundary condition is not met. Interestingly, this yields a Dirichlet boundary condition when swapping the order the operators are applied $L_{h}L_{h}^{*}$ as $h\to 0$.

We note that this procedure is similar to LDR-LLE. The row weights must sum to one to preserve constant functions, $S_{n}1=1$, and each point’s tangent coordinates are perfectly reconstructed since the coordinate functions are linear functions of themselves. Empirically we found these perform very similarly as shown in Figure 1, but Local Linear Regression was more likely to exclude the quadratic interaction polynomial from its bottommost eigenfunctions.

5.3 Laplacian without boundary conditions

For existing Laplacian approximations, the smoothness penalty $f^{T}L_{n}f\to\|\nabla f\|^{2}$ is only guaranteed to converge when $f$ satisfies the boundary conditions imposed by $L$. Furthermore, the only existing Laplace-Beltrami approximation that yields a positive semidefinite matrix and non-negative smoothness penalty is the graph Laplacian.

We derive a construction, the Coefficient Laplacian, that is both self-adjoint and guarantees convergence of the smoothness penalty for all $C^{\infty}(\mathcal{M})$ functions. This is achieved by explicitly estimating the gradient with a local linear regression. For each neighborhood $\mathcal{N}_{i}$, estimate normal coordinates $U_{\mathcal{N}_{i}}$ with respect to $X_{i}$. Define

[TABLE]

Thus, $S_{\mathcal{N}_{i}}Y_{\mathcal{N}_{i}}$ yields the coefficients for a linear regression of the centered $Y_{\mathcal{N}_{i}}-1Y_{i}$ on $U_{\mathcal{N}_{i}}$ without the intercept term. Thus, $f^{T}Lf=\sum_{i}f^{T}Q_{\mathcal{N}_{i}}f\to\sum_{i}\|\nabla f\|^{2}$ as $n\to\infty$ and $h\to 0$ sufficiently slowly.

The resulting operator has desirable behavior when used as a smoothing penalty. However, its eigenfunctions are useless for nonlinear dimensionality reduction. This occurs for two reasons. First, removal of boundary conditions yields an operator without a pure point spectrum. For example, in one dimension, the resulting operator solves the differential equation $y^{\prime\prime}=-\lambda y$ which has solutions $y=cos(\sqrt{\lambda}x)$ for all $\lambda>0$. Thus, it is unclear what solutions are picked out by an eigendecomposition of the discrete approximations. Second, one can show that $L_{h}f$ does not converge uniformly at the boundary even when the functional $f^{T}L_{h}f$ converges. As a result, the eigenvectors of the discrete approximations are uninformative. The resulting deficiencies are illustrated in Figure 2.

6 Discussion

This paper explores only one aspect of the connection between manifold learning and linear smoothing, namely the analysis of manifold learning methods and ways to use that analysis to improve existing NLDR methods.

However, this connection can be exploited in many additional ways. For instance, the methods developed for manifold learning can be used for linear smoothing. In particular, the linear smoother associated with LTSA can be used as a smoother for non-parametric regression problems. A simple consequence of our analysis is that the smoother can achieve $O(h^{4})$ bias in the interior even though it only uses a first-order local linear regression. By contrast, local linear regression yields a bias of $O(h^{2})$ in the interior. The connection can also be exploited for multi-scale analyses. In the case of the Laplacian, the weighted bias operator $L$ is the infinitesimal generator of a diffusion with $P_{t}=exp(-Lt)$ defining its transition kernels and can be used to generate Diffusion wavelets Coifman & Maggioni (2006). Other bias operators can be substituted to generate orthogonal bases with local, compact support. Another possibility is to combine existing linear smoothers with manifold learning penalties to generate new smoothers and NLDR techniques by using the theory of generalized additive models Buja et al. (1989).

Appendix A Convergence to differential operator

Theorem 6**.**

Let $S_{n}$ be a sequence of linear smoothers where the support of $S_{n}(x,\cdot)$ is contained in a ball of radius $h_{n}$. Further assume that the bias in the residuals $h_{n}^{-k}(I-S_{n})f=o(1)$ for some integer $k$ and all $f\in\mathcal{C}\subset C^{\infty}(\mathcal{M})$. Then if $h_{n}\to 0$ and $h_{n}^{-k}(I-S_{n})$ converges to a bounded linear operator as $n\to\infty$, then it must converge to a differential operator of order at most $k$ on the domain $\mathcal{C}$.

Proof.

Let $S_{\infty}$ is be the limit of $h^{-k}(I-S_{n})$. Let $q_{\ell}$ be any monomial with degree $\ell>k$ and $I_{x,n}$ be some smooth function that is $1$ on a ball of radius $h_{n}$ around $x$. Then the pointwise product $\|q_{\ell}\cdot I_{x,n}\|=O(h_{n}^{\ell})$. Since $S_{\infty}$ is bounded, $(S_{n}q_{\ell})(x)=(S_{n}(I_{x,n}\cdot q_{\ell}))(x)=O(h_{n}^{\ell-k})\to 0$ . Furthermore, the convergence is uniform over all $x$. Thus, the behavior is determined on a basis of polynomials of degree at most $k$. Applying $S_{\infty}$ to a Taylor expansion with remainder gives the desired result. ∎

Appendix B Equivalence of HLLE and LTSA

Theorem 7**.**

$L_{h}^{LTSA}-L_{h}^{HLLE}\to 0$* as $h\to 0$ in the weak operator topology of $C^{\infty}(\mathcal{M})$ equipped with the $L_{2}$ norm.*

Proof.

The local operators $Q^{HLLE}_{i}$ for HLLE can be described more succinctly in terms of the difference of two linear smoothers. Let $H_{x}^{(2)}$ be hat matrix for a local quadratic regression in the neighborhood $\mathcal{N}_{x}$.

[TABLE]

Let $R_{x}=H^{(2)}_{x}-I_{\mathcal{N}_{x}}$. For any $f,g\in C^{\infty}(\mathcal{M})$,

[TABLE]

since quadratic terms and lower can be removed from the Taylor expansion. ∎

Appendix C Boundary behavior of HLLE and LTSA

Theorem 8**.**

For any $x\in\partial\mathcal{M}$ and function $f\in C^{\infty}(\mathcal{M})$, $L^{LTSA}_{h}f(x)\to\infty$ as $h\to 0$ unless

[TABLE]

where $\eta_{x}$ is the tangent vector orthogonal to the boundary.

To do this, we show that the required scaling by $h^{-4}$ for the functional to converge causes the value at the boundary to go to $\infty$ when this condition is not met.

Let $x\in\mathcal{M}$. Assume $h$ is sufficiently small so that there exists a normal coordinate chart at $x$ such that $B(x,2h)$ is contained in the neighborhood for the chart. Furthermore, choose the normal coordinates such that the first coordinate $u_{1}$ corresponds to the direction $\eta_{x}$ orthogonal to the boundary and pointing inwards.

Consider the behavior of $L_{h}$ on the basis of polynomials expressed in this coordinate chart. Denote by $Q_{y}$ the hat operator for a linear regression on linear and constant functions restricted to a ball of radius $h$ centered on $y$. Note that any polynomial of degree at most 1 is in the null space of any $Q_{y}$ by construction. Likewise, by symmetry, any polynomial $u_{i}u_{j}$ with $j<i$ or of odd degree is in the null space of $Q_{y}$. This leaves only the polynomials $u_{i}^{2}$.

Each $(Q_{y}f)(z)$ computes the residual at $z$ from regressing $f$ against linear functions in the neighborhood at $y$. The residual $r_{i}=(Q_{y}u_{i}^{2})(0)=(Q_{y}(u_{i}-y_{i})^{2})(0)$ for all $y\in B(0,h)$ since linear functions are in the null space. By symmetry of the neighborhoods in the interior, $r_{y}=(Q_{0}u_{i}^{2})(-y)$ whenever $i>1$. In other words, in the interior of the manifold, averaging over the multiple regressions is equivalent to averaging over the residuals of a single regression. Since the expected residual is always [math] for least squares regression, $(L_{h}f)(x)=O(h^{4})$ for any $x\in Int(\mathcal{M})$.

More intuitively, one wishes to evaluate the residual at $z$ over all neighborhoods $\mathcal{N}_{y}$ that include $z$. Because linear functions are in the null space, this depends only on the curvature of the response $f$ at $x$. Furthermore, since all neighborhoods have essentially identical shape, the hat matrices for every neighborhood are also identical, so the residual only depends on the offset $z-y$. Averaging over the different $y$ is equivalent to fixing a single regression at $B(z,h)$ and averaging over the residuals at all offsets $y\in B(z,h)$.

However, the boundary behavior is very different since the shape of the neighborhoods can change. First, consider the quadratic polynomial $u_{1}^{2}$ in the direction orthogonal to the tangent space on the boundary. Since the residual is always evaluated at the boundary for a linear regression and $u_{1}^{2}$ is convex, the residual is necessarily positive. Figure 4 provides an intuitive illustration of this. Furthermore, since evaluating $(Q^{(\alpha h)}_{y}(\alpha h)^{-2}u_{1}^{2})(x)=O((\alpha h)^{-4})$ It follows that $(L_{h}f)(x)=\alpha h^{2}\partial^{2}f/\partial u_{1}^{2}(x)+O(h^{3})$ for some non-zero constant $\alpha$.

For polynomials $u_{i}^{2}$ where $i>1$, we note that any linear regression on linear functions $u_{j}$ has corresponding coefficients $\beta_{j}=0$ for all $j>1$ by symmetry of the neighborhood. Thus, the prediction can be reduced to a linear regression where the response is the conditional mean $E(u_{i}^{2}|u_{1})$ given the neighborhood. Integrating over slices of an $m$-sphere gives

[TABLE]

Since constant terms are removed by the regression and since summing over the residual function for a given value of $u_{1}$ gives ${\mathbb{E}}(u_{i}^{2}-\hat{u}_{i}^{2}|u_{1})={\mathbb{E}}(u_{i}^{2}-{\mathbb{E}}(u_{i}^{2}|u_{1})|u_{1})+{\mathbb{E}}({\mathbb{E}}(u_{i}^{2}|u_{1})-\hat{u}_{i}^{2}|u_{1})={\mathbb{E}}({\mathbb{E}}(u_{i}^{2}|u_{1})-\hat{u}_{i}^{2}|u_{1})$, the summed residuals are the same as those obtained for the function $-u_{1}^{2}/(m+1)$.

Thus, for any $x\in\partial\mathcal{M}$, $(L_{h}f)(x)\to\infty$ unless $0=\partial^{2}f/\partial u_{1}^{2}+\frac{m+1}{m-1}\sum_{i=2}^{m}\partial^{2}f/\partial u_{i}^{2}$ which is equivalent to equation 19.

We checked this boundary condition using simulation on a manifold isomorphic to a rectangle. We took 10 estimated eigenfunctions and computed their Hessian at a point on the boundary. This generates a $10\times 3$ matrix consisting of the estimates $\frac{\partial^{2}f}{\partial u_{1}^{2}},\frac{\partial^{2}f}{\partial u_{2}^{2}},\frac{\partial^{2}f}{\partial u_{1}\partial u_{2}}$. We take the svd of this matrix. The distribution of the top, middle, and bottom singular value is shown in figure 5. There is one eigenvalue clearly close to 0 that represents the boundary condition. The average bottom right singular vector is given in table 3 and show fairly good correspondence to the theoretical calculations on a modestly fine grid.

Appendix D Boundary Behavior of Local linear regression

As the boundary bias of local linear regression is already well studied, existing results can help determine the boundary conditions for the resulting operator. However, as the local linear regression smoother is not self-adjoint, the behavior of the adjoint must also be determined. Let $x\in\partial\mathcal{M}$ and $\mathbf{u}$ be normal coordinates at $x$ in a neighborhood $\mathcal{B}$ containing $B(x,2h)\cap\mathcal{M}$. We again choose $u_{1}$ to correspond to the tangent vector that is normal to the boundary and pointing inwards so that $u_{1}\geq 0$. We wish to evaluate the boundary behavior of $L^{*}_{h}L_{h}f(x)=\langle L_{h}\delta_{x},L_{h}f\rangle$ where $L^{*}_{h}$ is the adjoint of $L_{h}$ and $\delta_{x}$ is the Dirac delta. The main part of our proof is to evaluate $L_{h}\delta_{x}$.

Similar to the proof for theorem 8, we first show we can reduce the problem to a univariate linear regression. We do this through an orthogonal basis. Then we use the usual regression equations to actually compute the value. For each $y\in B(x,h)\cap\mathcal{M}$, denote $\mathcal{N}_{y}=B(y,h)\cap\mathcal{M}$ and $\tilde{u}^{(0)}_{i}(\cdot,y)=h^{-1}Vol(\mathcal{N}_{y})^{-1}(u_{i}(\cdot)1(\cdot\in\mathcal{N}_{y})-u_{i}(y))$. It is easy to see that by symmetry that there are functions $\tilde{u}_{i}(\cdot,y)=\tilde{u}^{(0)}_{i}(\cdot,y)+\epsilon_{i}$ which are orthogonal to each other and to the constant function in $L_{2}(\mathcal{N}_{y})$ for $i>1$ and where $\epsilon_{i}=O(h^{2})$. Likewise, let $\tilde{u}_{1}(\cdot,y)$ be similarly defined to be orthogonal to $\tilde{u}_{i}(\cdot,y)$ for $i>1$. To generate an orthonormal basis, we must orthogonalize it with respect to the constant function as well. Gram-Schmidt gives $v_{1}(\cdot,y)=\tilde{u}_{1}(\cdot,y)-Vol(\mathcal{N}_{y})^{-1}\langle\tilde{u}_{1}(\cdot,y),I(\cdot\in\mathcal{N}_{y})\rangle I(\cdot\in\mathcal{N}_{y})$ is orthogonal to the constant function. Thus, by evaluating

[TABLE]

where the inner product is taken with respect to $L_{2}(\mathcal{B})$. Integrating over $\mathcal{B}$, the first two terms each have magnitude 1 and cancel each other out. The third term may be rewritten

[TABLE]

where $\mu_{y}=Vol(\mathcal{N}_{y})^{-1}\langle\tilde{u}_{1}(\cdot,y),I(\cdot\in\mathcal{N}_{y})\rangle>0$ unless $u_{1}(y)(1+O(h^{2}))\geq h$.

It follows that

[TABLE]

By applying this to a Taylor expansion, one concludes that for any continuously differentiable function which is non-zero at $x\in\partial\mathcal{M}$, $|(L_{h}^{*}f)(x)|\to\infty$ as $h\to 0$.

By Theorem 2.2 in Ruppert & Wand (1994), a point $x$ with distance $dist(x,\partial\mathcal{M})<h$ has boundary bias

[TABLE]

where $\alpha>0$ is some constant and $D_{x,h}$ is the unit $m$-sphere cut along the plane orthogonal to the first coordinate at $u_{1}=dist(x,\partial\mathcal{M})/h$. This is a linear function in the Hessian. Furthermore, all odd moments of $u_{i}$ for $i>1$ are 0, and their second moments are all equal. It is a linear function of the diagonal of the Hessian and of the form $(L_{h}f)(x)=-\beta\frac{\partial^{2}f}{\partial u_{1}^{2}}(x)+(\Delta f)(x)+o(1)$ for some $\beta\neq 0$.

Thus, the boundary condition that functions must satisfy is

[TABLE]

Although it takes the same form as the HLLE boundary condition, the constants are different. However, we do not know of a reason to prefer one boundary condition over the other.