text:
\section{Robust node generation on irregular domains and surfaces}
\label{sec:node_gen}

Having tested our geometric modeling technique, we present our algorithm (Algorithm \ref{alg:node_gen_1}) for generating scattered nodes on irregular domains $\Omega \subset \mathbb{R}^d$, $d=2,3$ using a combination of Poisson disk sampling and the new geometric modeling technique. Algorithm \ref{alg:node_gen_1} has the important feature that it ensures that interior nodes maintain a user-specified separation distance of $h$ from the boundary and between each other. In the remainder of this section, we will discuss the salient features of Algorithm \ref{alg:node_gen_1}, present some modifications to allow for time-varying embedded boundaries, explore the quality of the node sets it produces (using histograms), and test its suitability for RBF-FD discretizations.
\begin{algorithm}
\caption{Node generation for domains with smooth boundaries}
\label{alg:node_gen_1}
\begin{algorithmic}[1]
  \Statex{\bf Given}: $\chid = \{(\vX_d)_{\ell}\}_{\ell=1}^{N_d}$, a set of ``seed'' nodes  on domain boundary.
	\Statex{\bf Given}: $h$, the average separation distance between nodes.
	\Statex{\bf Generate}: $\chib = \{(\vX_b)_j\}_{j=1}^{N_b}$, a set of boundary nodes with spacing $h$.
	\Statex{\bf Generate}: $\eta_b$, the set of outward unit normals on the boundary.
	\Statex{\bf Generate}: $\chii = \{\vx_k\}_{k=1}^{N_i}$, a set of interior nodes with spacing $h$.
	\Statex{\bf Generate}: $\rchi = \chii \cup \chib$.
	\State Obtain $\chib$ using Algorithm \ref{alg:ss_decim}.
	\State Evaluate derivatives of \eqref{eq:sbf_poly_model} at the parametric evaluation points obtained from Algorithm \ref{alg:ss_decim} to obtain $\eta_b$.
	\State Use the normals $\eta_b$ to project $\chib$ inwards a distance $h$, giving a set $\hat{\chib}$; this set defines an inner boundary and an inner domain.
	\State Build a kd-tree on the set $\hat{\chib}$.
	\State Generate the oriented bounding volume (OBV) corresponding to $\hat{\chib}$ using Algorithm \ref{alg:pca_obb}.
	\State Fill the OBV with points of (approximate) spacing $h$ using Poisson disk sampling in Algorithm \ref{alg:pds}~\cite{Bridson07}.
	\State Test each sample against the outward normal at its closest inner boundary point (found using the kd-tree). If it is outside the inner domain, discard.
	\State All remaining Poisson disk samples (not including $\hat{\chib}$) form the set $\chii$.
	\State Set $\rchi = \chii \cup \chib$.
\end{algorithmic}
\end{algorithm}

\subsection{Sampling on surfaces}
\label{sec:node_gen_surf}

A node generation algorithm for RBF-FD methods must generate boundary nodes that are approximately a user-specified distance $h$ apart. Since all domain boundaries are submanifolds of $\mathbb{R}^d$, this problem involves generating a set of parametric evaluation nodes that result in approximately quasi-uniform Cartesian node sets on surfaces. We present an algorithm based on two simple ideas: supersampling of our geometric model, and decimation. As implemented in Algorithm \ref{alg:ss_decim}, we first supersample the parametric interpolant \eqref{eq:sbf_poly_model} to the surface at a large set of parametric evaluation points in $[-\pi,\pi)$ (curves) and $[-\pi,\pi)\times [-\pi/2, \pi/2)$ (surfaces). The resulting Cartesian nodes are not quasi-uniformly spaced, but tend to ``bunch'' according to the parametric map. Thus, the Cartesian node set is now thinned (or \emph{decimated}) to approximately enforce that all nodes be no closer than some separation distance $h$. This requires only the implementation of a ball query (range search) algorithm-- in our case provided by our kd-tree implementation-- in conjunction with a simple depth-first traversal of the node set (in any order). This approach is often referred to as \emph{sample elimination}. A strength of our approach is that the resulting set of samples comes from the geometric model, and retains the high-order accuracy conferred by the model.
\begin{algorithm}
\caption{Sampling on surfaces}
\label{alg:ss_decim}
\begin{algorithmic}[1]
  \Statex{\bf Given}: $\chid = \{(\vX_d)_{\ell}\}_{\ell=1}^{N_d}$, a set of ``seed'' nodes  on domain boundary.
	\Statex{\bf Given}: $h$, the average separation distance between nodes.
	\Statex{\bf Given}: $\tau$, the supersampling parameter.
	\Statex{\bf Generate}: $\chib = \{(\vX_b)_j\}_{j=1}^{N_b}$, a set of boundary nodes  with spacing $h$.
	\Statex{\bf Generate}: $\Lambda^e$, a set of parametric evaluation nodes for \eqref{eq:sbf_poly_model}
	\Statex{\bf Generate}: the SBF geometric model in \eqref{eq:sbf_poly_model}.
	\State Using \eqref{eq:sbf_poly_model}, fit a geometric model to the nodes in $\chid$.
	\State To find $N_b$ corresponding to a value of $h$, estimate the surface area (or perimeter) $a_d$ of the OBV corresponding to $\chid$ (generated by Algorithm \ref{alg:pca_obb}). Then, $N_b = a_d h^{-(d-1)}$.
	\State Evaluate \eqref{eq:sbf_poly_model} at $\hat{N}_b = \tau N_b$ parametric evaluation points either in $[-\pi,\pi)$ or $[-\pi,\pi)\times [-\pi/2, \pi/2)$. Let $\hat{\chib}$ be the resulting candidate set of Cartesian evaluation points.
	\State Build a kd-tree on $\hat{\chib}$.
	\State Initialize $g$, an $\hat{N}_b$ array of flags, to 1.
	\State Store $\hat{\chib}$ in $\hat{N}_b \times d$ matrix $\hat{X}_b$.
	\For {$k=1,\hat{N}_b$}
		\If {g(k)$ \neq 0$}
				\State Set idxs = indices of points within distance $h$ of $\hat{X}_b$ (using, for instance, a kd-tree ball query).
				\State Set g(idxs $\neq k$) = 0.
		\EndIf
	\EndFor
	\State Set $\hat{g} =$ find(g=1), the vector of indices corresponding to flags of 1.
	\State Collect all parametric evaluation nodes corresponding to indices in $\hat{g}$ into array $\Lambda^e$.
	\State Evaluate t\eqref{eq:sbf_poly_model} at $\Lambda^e$ to obtain $\chib$, the $N_b$ Cartesian points on the boundary.
\end{algorithmic}
\end{algorithm}
Algorithm \ref{alg:ss_decim} clearly only approximately enforces that points be a distance $h$ apart in the \emph{Euclidean} norm, much like Bridson's algorithm for Poisson disk sampling~\cite{Bridson07}. We will explore the spatial distributions of the node sets produced by Algorithm \ref{alg:ss_decim} in Section \ref{sec:res_hist1}.

\subsection{Oriented bounding boxes (OBBs) from PCA}
\begin{algorithm}
\caption{OBB generation using PCA}
\label{alg:pca_obb}
\begin{algorithmic}[1]
  \Statex{\bf Given}: $X_b$, the $N_b \times d$ matrix containing a set of points for which we want an OBB.
	\Statex{\bf Generate}: $B$, the matrix of OBB vertices.
	\State Compute $\vx_c = \frac{1}{N_b} \begin{bmatrix} \sum\limits_{i=1}^{N_b} X_b(i,1) & \sum\limits_{i=1}^{N_b} X_b(i,2) & \sum\limits_{i=1}^{N_b} X_b(i,3) \end{bmatrix}$, the centroid of $X_b$.
	\State Compute $M = \mathbf{1} \otimes \vx_c$, the $N_b \times d$ matrix with $\vx_c$ as each of its rows.
	\State Compute $C = \frac{1}{N_b}(X_b-M)^T (X_b-M)$, the $d \times d$ normalized covariance matrix corresponding to $X_b$.
	\State Decompose $C$ as $C =  VDV^T$, where $V$ contains eigenvectors (rotations about origin), $D$ contains eigenvalues (scaling).
	\State Compute unrotated boundary points $\hat{X}_b = X_b V$.
	\State Find column minimum and maximum values for $\hat{X_b}$, \emph{i.e.,} two vertices of the unrotated bounding box.
	\State Using these two vertices, find the side lengths of the box along each coordinate direction.
	\State Using two vertices and side lengths, find the other $2^d - 2$ vertices of the bounding box.
	\State Store all unrotated bounding box vertices in $\hat{B}$, a $2^d \times d$ matrix.
	\State Compute the OBB vertices $B = \hat{B} V^T$.
\end{algorithmic}
\end{algorithm}
Algorithm \ref{alg:node_gen_1} requires an oriented-bounding volume (OBV) for node generation. Almost any simple shape will serve as a bounding volume, but we restrict ourselves to oriented bounding boxes (OBBs), since intersection tests of vectors against boxes simply involve testing a set of linear inequalities. Further, the Poisson disk sampling algorithm is trivial to implement for a rectangular/cuboidal domain. Our goal is to use a technique that has low computational complexity, is easy to implement, and is meshfree. Further, a technique that produces approximately minimal volume OBBs is desirable, since this will result in fewer inside/outside tests. Principal Component Analysis (PCA) fits all our requirements~\cite{Dimitrov2006}, and can easily be used to generate an OBB for our domain in $O(N_b)$ operations. This procedure is summarized in Algorithm \ref{alg:pca_obb}.

\subsection{Poisson disk sampling}

\comment{Algorithm \ref{alg:node_gen_1} is primarily structured around Poisson disk sampling inside the OBB generated by Algorithm \ref{alg:pca_obb}. Our approach for Poisson disk sampling in the OBB is a straightforward implementation of the algorithm described in~\cite{Bridson07}. For completeness, we present this procedure in Algorithm \ref{alg:pds}, adapted to our notation.
\begin{algorithm}
\caption{Fast Poisson disk sampling}
\label{alg:pds}
\begin{algorithmic}[1]
  \Statex{\bf Given}: $h$, the minimum distance between nodes.
	\Statex{\bf Given}: $d$, the spatial dimension.
	\Statex{\bf Given}: $\hat{k}$, the number of samples to choose before rejection (the Poisson neighborhood size).
	\Statex{\bf Given}: $B$, the matrix of vertices of the domain OBB.
	\Statex{\bf Generate}: The set of $N_{obb}$ Poisson disk samples in the domain OBB.
	\State Initialize a $d$-dimensional Cartesian background grid $\calG$ with cell size $\frac{h}{\sqrt{d}}$. A value of $-1$ in $\calG$ indicates no sample, any non-negative integer is the \emph{linear index} of the sample in a cell.
	\State Generate the first uniform random node $\vx_0$ within the domain OBB. Insert this node into $\calG$.
	\State Initialize the array of node indices (active list) $\calI$ with the index of $0$ (corresponding to $\vx_0$).
	\While{$\calI$ is not empty}
	\State Choose a random index $i$ from $\calI$.
	\State Generate $\hat{k}$ uniform random nodes in the spherical annulus centered at $\vx_i$ with inner radius $h$ and outer radius $2h$. Place these nodes in set $\mathcal{P}$.
	\State Compare all of $\mathcal{P}$ against existing nodes within distance $h$ (using $\calG$ to facilitate comparisons).
	\State If none of the nodes in $\mathcal{P}$ are sufficiently far from existing samples, remove $i$ from the active list.
	\State If any of the nodes in $\mathcal{P}$ are sufficiently far from existing samples, add their indices to the active list, store these \emph{valid samples} in the set $\mathcal{V}$.
	\State If any of the nodes in $\mathcal{P}$ are outside the domain OBB, reject them.
	\EndWhile
	\State The set of $N_{obb}$ Poisson disk samples is $\mathcal{V}$.
\end{algorithmic}
\end{algorithm}
The complexity of this algorithm is $O(\hat{k} N_{obb})$. Since $\hat{k}$ is a constant, the cost scales as $O(N_{obb})$. In Section \ref{sec:results}, we will explore the impact of $\hat{k}$ on the uniformity of node distributions and on the accuracy for RBF-FD discretizations of PDEs. We explore the node distributions obtained from the use of Algorithm \ref{alg:pds} within Algorithm \ref{alg:node_gen_1} in Section \ref{sec:res_hist2}.}

\subsection{Domains with time-varying embedded boundaries}

Biological problems often involve domains with not only irregular outer boundaries, but irregular embedded inner boundaries, \emph{e.g.,} platelets and red blood cells in a blood vessel~\cite{SWFKIJNMF2015}. While the authors have begun developing numerical methods for simulations in such complex domains~\cite{SWFKIJNMF2014,ShankarJCP2017,SFJCP2018}, our methods currently lack a robust node generation algorithm that can modify its node sets \emph{locally} when the embedded inner boundaries deform, and are added or removed.

\begin{figure}[h!]
\centering
\subfloat[Star domain]
{
	\includegraphics[scale=0.3]{"Figures/unemb_bdry2D".pdf}
	\label{fig:ea}
}
\subfloat[Star with embedded ellipse]
{
	\includegraphics[scale=0.3]{"Figures/emb_bdry2D".pdf}

	\label{fig:eb}
}
\caption{Local node modification for embedded boundaries using Algorithm \ref{alg:node_mod}. The embedded ellipse is shown using circles.}
\label{fig:node_mod}
\end{figure}

Fortunately, Algorithm \ref{alg:node_gen_1} is easily modified to tackle this problem. Consider the scenario where Algorithm \ref{alg:node_gen_1} has generated node sets $\chib$ and $\chii$ on the boundary and interior of an irregular domain $\Omega_0$ without any embedded boundaries. Now consider the inclusion of irregular boundaries $\Gamma = \{\Gamma_j\}_{j=1}^{N_p}$ enclosing domains $\{\Omega_j\}_{j=1}^{N_p}$ into $\Omega_0$ in such a way that these boundaries define a new domain $\Omega = \Omega_0 \backslash \bigcup\limits_{j=1}^{N_p}\Omega_j$. Our node generator must automatically and efficiently generate a node set in $\Omega$ that \emph{deletes} the nodes contained in $\tilde{\Omega} = \bigcup\limits_{j=1}^{N_p}\Omega_j$. In fact, our algorithm also deletes any nodes that are within an $h$ distance of the embedded boundaries, again ensuring that the global separation distance is approximately $h$. This is detailed in Algorithm \ref{alg:node_mod}.
\begin{algorithm}
\caption{Node set modification for (closed) embedded boundaries}
\label{alg:node_mod}
\begin{algorithmic}[1]
	\Statex{\bf Given}: $h$, the average separation distance between nodes in domain $\Omega_0$.
	\Statex{\bf Given}: $\alpha$, the ratio of the numbers of inner boundary nodes to outer boundary nodes (see Section \ref{sec:nodemod_complex}).
	\Statex{\bf Given}: $\chib = \{(\vX_b)_j\}_{j=1}^{N_b}$, a set of boundary nodes on $\partial \Omega_0$.
	\Statex{\bf Given}: $\eta_b$, the set of outward unit normals on the boundary $\partial \Omega_0$.
	\Statex{\bf Given}: $\rchi = \{\vx_k\}_{k=1}^{N}$, the set of nodes on $\Omega_0$ (including the boundary).
	\Statex{\bf Given}: $N_d^{\Gamma}$ seed nodes on each of the embedded boundaries $\Gamma_j, j=1,\hdots,N_p$.
	\Statex{\bf Generate}: $Z$, the $N \times 2$ array whose first column indicates whether a point $\vx_k \in \Omega$, and second column is $j$ if $\vx_k \in \Omega_j$.
	\Statex{\bf Generate}: $\tilde{\rchi}$, the (modified) set of $\tilde{N}$ nodes on the irregular domain $\Omega = \Omega_0 \backslash \tilde{\Omega}$.
	\State Using \eqref{eq:sbf_poly_model}, fit a geometric model to the seed nodes on each embedded boundary $\Gamma_j,j=1,\hdots,p$.
	\State Set $N_b^{\Gamma} = \alpha N_b$.
	\State For $j=1,\hdots,N_p$, obtain $\chibgj$, the $N_b^{\Gamma}$ boundary points on the $j$th embedded boundary using Algorithm \ref{alg:ss_decim}.
	\State Evaluate derivatives of the $N_p$ boundary interpolants at $N_b^{\Gamma}$ parametric points to obtain $\eta_b^{\Gamma_j}$, the set of unit normals on each boundary (pointing into $\Omega_0$).
	\State Extend each embedded boundary $\Gamma_j$ by distance $h$ in their normal directions to obtain $\tilde{\Gamma}_j$ and the corresponding points $\tilde{\chibgj}$.
	\State Let $\tilde{\chibg} = \bigcup\limits_{j=1}^{N_p} \tilde{\chibgj}$.
	\State Add $\tilde{\chibg}$ to the domain kd-tree currently containing $\chib$.
	\State Generate the OBB for each $\tilde{\Gamma}_j$ using Algorithm \ref{alg:pca_obb} on the $N_b^{\Gamma}$ boundary points.
	\State Initialize the 2D array (map) Z with dimensions $N \times 2$ to zeros.
	\For {$k=1,N$}
		\For{$j=1,N_p$}
				\If {$\vx_k$ is within the $j$\textsuperscript{th} OBB}
					\State Find the closest point $\zeta$ on $\tilde{\Gamma}_j$ using the domain kd-tree.
					\State Use $\zeta$ to test $\vx_k$ against $\tilde{\Gamma}_j$.
					\If {$\vx_k$ is inside $\tilde{\Gamma}_j$}
						\State Set $Z(k,1)$ = 1, $Z(k,2) = j$.
						\State Break and go to next $k$ value.
					\EndIf
				\EndIf
		\EndFor
	\EndFor
	\State Let g = find(Z(:,1)=0 be the list of indices of points with $\Omega$.
	\State Set $\tilde{\rchi} = \rchi(g=0)$, $\tilde{N}$ = length($\tilde{\rchi}$).
\end{algorithmic}
\end{algorithm}
This procedure involves only locally changing the original node set as new objects are introduced, and also accounts for cases where an embedded boundary may touch the domain outer boundary. It is important to note that this algorithm features potential early termination: a point is only tested against bounding boxes and embedded boundaries until it is found to be contained within an embedded boundary. While other efficiency improvements are possible, we forgo them for simplicity. Figure \ref{fig:node_mod} shows an example of the node sets obtained using Algorithm \ref{alg:node_mod}. In Figure \ref{fig:ea}, we show the interior nodes on the star domain. Figure \ref{fig:eb} shows the node set obtained when an ellipse is embedded into the star domain. It is easy to see that Algorithm \ref{alg:node_mod} only locally modifies node sets, as desired.

Algorithm \ref{alg:node_mod} does not take into account the case where an embedded boundary is removed. There are two options within the above framework for modifying the node sets in such a situation. The first option is to fully generate a new node set on $\Omega$ after an embedded boundary is removed; this would involve using Algorithm \ref{alg:node_mod} to modify the original node set on $\Omega_0$ to account for the embedded boundaries that were \emph{not removed}. This option is rather wasteful. Instead, Algorithm \ref{alg:node_mod} associates each grid point $\vx_k$ with an embedded boundary $\Omega_j$. Thus, to handle a \emph{removed} embedded boundary $\Omega_j$, we need only find all nodes in the original node set that are contained within $\Omega_j$ by searching the second dimension of the array $Z$ for the value $j$. Of course, this second option relies on knowing which of the embedded boundaries were removed. If this is not known, the first option discussed above is preferable.

\subsection{Boundary refinement and ghost nodes}

The numerical solution of PDEs with RBF-FD sometimes requires denser node sets near domain boundaries~\cite{FlyerNS,FlyerPHS,ShankarJCP2017}, and/or ghost nodes outside the domain boundary to enforce boundary conditions~\cite{FlyerNS,BarnettPHS,SFJCP2018}. Our node generator must therefore possess these capabilities. Generating ghost nodes is straightforward: given a set of boundary nodes and unit outward normals, simply copy the nodes a distance $h$ outside the domain along the outward normals. Further, Algorithm \ref{alg:node_gen_1} can be easily modified to provide a node set refined near a boundary: simply copy the boundary nodes at some user-specified distance inward from the boundary. This places a layer of nodes between $\chib$ and the inner boundary $\hat{\chib}$. This approach generalizes straightforwardly to multiple layers of boundary refinement. However, a drawback of this approach is that it does not generate \emph{graded} boundary-refined node sets,~\emph{i.e.}, node sets that smoothly vary in space. We leave this generalization for future work.

\input{Results}

\section{Complexity Analysis}
\label{sec:complex}

We now derive the computational complexity of the Algorithms \ref{alg:node_gen_1} and \ref{alg:node_mod}.

\subsection{Complexity of one-time node generation}
\label{sec:nodegen_complex}

The total computational cost of our node generator is a sum of the individual costs of the following steps:
\begin{itemize}
\item Forming a kd-tree on the domain boundary: The kd-tree is formed on $N_b$ points for a cost of $C_1 = O(d N_b \log N_b)$. Each subsequent nearest neighbor search costs $O(\log N_b)$ operations.

\item SBF boundary representation: Since this involves $N_d$ seed nodes, it requires the solution of the $N_d \times N_d$ dense linear system \eqref{eq:sbf_linsys}, which has an asymptotic cost of $O(N_d^3)$. This interpolant is evaluated to obtain $\tau N_b$ boundary points for a cost of $O(N_b N_d)$, then thinned for a worst-case cost of $O(N_b \log N_b)$ using the kd-tree. If we assume $N_d = \gamma N_b$, where $\gamma \in (0,1]$ is some small number, the evaluation cost can be written as $O(\gamma N_b^2)$. Cost: $C_2 = O(\gamma^3 N_b^3) + O(\gamma N_b^2) + O(N_b \log N_b)$.

\item OBB: The primary cost in the PCA computation of the OBB is the $O(d^2 N_b)$ cost of forming the $d \times d$ covariance matrix. The subsequent eigendecomposition of the $d \times d$ matrix is computed in $d^3$ flops. Then, the rotations and max/min computations to find the unrotated bounding box cost $O(N_b)$. The costs of rotations of the bounding boxes can be neglected when $d=2,3$. Cost: $C_3 = O(N_b)$.

\item Poisson disk sampling in the OBB: Using Algorithm \ref{alg:pds}, Poisson disk samples for the OBB can be computed for a cost of $C_4 = O(N_{obb})$.

\item Eliminate OBB nodes outside the domain boundary: This requires finding the closest boundary point to each of the $N_{obb}$ Poisson disk samples in the OBB using the kd-tree. The total cost is therefore given by $C_5 = O(N_{obb} \log N_b)$.
\end{itemize}

The total cost is given by $C_{tot} = \sum\limits_{i=1}^5 C_i$. Dropping the big-$O$ notation for convenience, we have
\begin{align}
C_{tot} = d N_b \log N_b + \gamma^3 N_b^3 + \gamma N_b^2 + N_b \log N_b + N_b + N_{obb} + N_{obb} \log N_b.
\label{eq:complexity}
\end{align}
All our costs are currently in terms of $N_b$ and $N_{obb}$. However, it is more appropriate to express $C_{tot}$ in terms of the total number of points in the domain $N = N_i + N_b$. This requires dimension-specific simplifications.

\subsubsection{2D complexity estimates}

We now present the complete 2D complexity estimates in terms of $N$. Assuming the nodes are approximately uniformly spaced on both the boundary and the interior, we have
\begin{align}
h_i = \lf(\frac{a}{N_i}\rt)^{\frac{1}{2}}, h_b = \frac{q}{N_b}.
\end{align}
If we wish the nodes to be spaced comparably on the boundary and interior, it is reasonable to assume that $h_i = h_b = h$. Solving for $N_b$ then gives
\begin{align}
N_b &= q a^{-\frac{1}{2}} N_i^{\frac{1}{2}}.\label{eq:Nb2d}
\end{align}
Substituting \eqref{eq:Nb2d} in \eqref{eq:complexity} and neglecting sublinear terms in $N_i$ gives us
\begin{align}
C_{tot} = \gamma^3 q^3 a^{-\frac{3}{2}} N_i^{\frac{3}{2}} + \gamma \frac{q^2}{a} N_i + 2q a^{-\frac{1}{2}} N_i^{\frac{1}{2}} + 2N_{obb}.
\label{eq:ctot1}
\end{align}
We must now relate $N_{obb}$ to $N_i$ to get the total cost in terms of $N_i$. We know that
\begin{align}
h_i &= \lf(\frac{a}{N_i}\rt)^{\frac{1}{2}} = \lf(\frac{a_{obb}}{N_{obb}}\rt)^{\frac{1}{2}}, \\
\implies N_{obb} &= \frac{a_{obb}}{a} N_i.
\end{align}
Substituting this new expression into \eqref{eq:ctot1} and retaining only linear and higher-order terms in $N_i$, we obtain
\begin{align}
C_{tot} = \gamma^3 q^3 a^{-1.5} N_i^{1.5} + \gamma \frac{q^2}{a} N_i + 2\frac{a_{obb}}{a} N_i.
\end{align}
Since $N = N_i + N_b$ and $N_b << N_i$, it is reasonable to replace $N_i$ by $N$, giving us
\begin{align}
C_{tot} = \gamma^3 q^3 a^{-1.5} N^{1.5} + \gamma \frac{q^2}{a} N + 2\frac{a_{obb}}{a} N.
\end{align}
The constant in front of the leading-order $N^{1.5}$ term is typically small in practice for large $N$ and smooth domain boundaries. In such a case, the $O(N)$ terms dominate. Thus, for small $\gamma$ and reasonably compact domains, we have
\begin{align}
C_{tot} \approx \lf(\gamma \frac{q^2}{a} + 2\frac{a_{obb}}{a}\rt) N.
\end{align}
Our node generator thus has a theoretical computational complexity of $O(N)$ in 2D, despite the use of a global RBF interpolant on the boundary. As we show in Section \ref{sec:scale_res}, we recover the $O(N)$ complexity in practice.

\subsubsection{3D complexity estimates}

We proceed much in the same way as in the 2D case, but obtain a different estimate. First, assuming the nodes are approximately uniformly spaced on both the boundary and the interior, we now have
\begin{align}
h_i &= \lf(\frac{v}{N_i} \rt)^{\frac{1}{3}}, h_b = \lf(\frac{s}{N_b} \rt)^{\frac{1}{2}}.
\end{align}
Now, letting $h_i = h_b = h$, solving for $N_b$ gives
\begin{align}
N_b & = s v^{-\frac{2}{3}} N_i^{\frac{2}{3}}. \label{eq:Nb3d}
\end{align}
Using \eqref{eq:Nb3d} in \eqref{eq:complexity} and neglecting sublinear terms in $N_i$, we have
\begin{align}
C_{tot} = \gamma^3 s^3 v^{-2} N_i^2 + \gamma s^2 v^{-\frac{4}{3}} N_i^{\frac{4}{3}} + 2 N_{obb}.
\end{align}
Once again relating $N_{obb}$ to $N_i$ using the volume of the OBB $v_{obb}$, we have
\begin{align}
C_{tot} = \gamma^3 s^3 v^{-2} N_i^2 + \gamma s^2 v^{-\frac{4}{3}} N_i^{\frac{4}{3}} + 2 \frac{v_{obb}}{v} N_i.
\end{align}
For a sufficiently large number of $N_b$ (small $h$) and fixed $N_d$, $\gamma$ can be as low as $O(10^{-3})$. Thus, both the quadratic and superlinear terms in the above expansion may be small for a compact domain, leaving us with
\begin{align}
C_{tot} \approx 2 \frac{v_{obb}}{v} N_i \approx 2 \frac{v_{obb}}{v} N.
\end{align}
Even if the superlinear term is retained, this represents a complexity of $O(N^{1.\bar{33}})$. However, in practice, we observe close to $O(N)$ complexity, as we show in Section \ref{sec:scale_res}.

\subsection{Complexity of node set modification}
\label{sec:nodemod_complex}

We now derive the complexity associated with Algorithm \ref{alg:node_mod}. The total \emph{worst-case} computational cost for removing nodes from the interior of $N_p$ embedded boundaries can be written as
\begin{align}
C_{mod} = O( N_p(N_d^{\Gamma})^3) + O(N_p N_b^{\Gamma} (N_d^{\Gamma})^2) + O(N_p N_b^{\Gamma}) + O(N N_p),
\end{align}
where the first term corresponds to fitting the SBF geometric model to each embedded boundary, the second term to evaluations of the models to obtain normals, the third term to the computation of the OBBs of the embedded boundaries, and the last term to testing each of the $N$ domain nodes against the $N_p$ OBBs. Now, let
\begin{align}
N_b^{\Gamma} = \alpha N_b, N_d^{\Gamma} = \alpha N_d.
\end{align}
If approximately the same spacing between nodes is maintained on all boundaries, we have $h_{\Gamma} = h_b$, which in 2D gives us $\alpha = \frac{q^{\Gamma}}{q}$ (rounded to the nearest integer). Similarly, in 3D, we have $\alpha = \frac{s^{\Gamma}}{s}$. Dropping the big-$O$ notation for convenience, $C_{mod}$ can be written as
\begin{align}
C_{mod} = \alpha\gamma^2(\alpha^2\gamma + 1) N_p N_b^3 + \alpha N_p N_b + N N_p.
\label{eq:cmod}
\end{align}
Using the fact that $N = N_i + N_b$, \eqref{eq:cmod} can be written as
\begin{align}
C_{mod}^{2\rm D} = \alpha\gamma^2(\alpha^2\gamma + 1)  N_p q^3 a^{-1.5} N_i^{1.5} + \alpha N_p q a^{-\frac{1}{2}} N^{\frac{1}{2}} + N N_p.
\end{align}
By a similar argument as in Section \ref{sec:nodegen_complex}, the $N_i^{1.5}$ term is small, leaving an estimate of
\begin{align}
C_{mod}^{2\rm D} \approx N N_p,
\end{align}
where sublinear terms in $N$ and $N_i$ have been neglected. Since the number of embedded boundaries is always much smaller than the number of nodes, \emph{i.e.,} $N_p << N$, we have
$C_{mod}^{2\rm D} = O(N)$. The 3D derivation is similar. Using the definition of $N_b$ in terms of $N_i$ in 3D, \eqref{eq:cmod} can be re-written as
\begin{align}
C_{mod}^{3\rm D} =  \alpha\gamma^2(\alpha^2\gamma + 1) N_p s^3 v^{-2} N_i^{2} + \alpha N_p s v^{-\frac{2}{3}} N^{\frac{2}{3}} + N N_p.
\end{align}
Again, in practice, the linear term in $N$ dominates.  Thus, we have
\begin{align}
C_{mod}^{3\rm D} \approx N N_p,
\end{align}
which is in practice $O(N)$ since $N_p << N$. If a bounding box hierarchy is used to keep track of the embedded boundaries, the $O(N N_p)$ term can be further shrunk to $O(N \log N_p)$, but we leave this approach for future work. It is important to note that our estimates for node generation only constitute the worst-case. In general, it is unlikely that a Poisson disk sample needs to be tested against \emph{all} OBBs of the embedded inner boundaries. Indeed, the disk sample is deleted as soon as it is found to lie within any embedded boundary, requiring no additional tests. The true complexity estimate in this scenario is non-deterministic. We leave this analysis for future work also.

\subsection{Scaling of the node generator}
\label{sec:scale_res}

\begin{figure}[h!]
\centering
\subfloat[Star, $N_d = 128$]
{
	\includegraphics[scale=0.5]{"Figures/2DCostsNd128".pdf}

	\label{fig:s1a}
}
\subfloat[Star, $N_d = 256$]
{
	\includegraphics[scale=0.5]{"Figures/2DCostsNd256".pdf}
	\label{fig:s1b}
}

\subfloat[Bumpy sphere, $N_d = 400$]
{
	\includegraphics[scale=0.5]{"Figures/3DCostsNd400".pdf}

	\label{fig:s1c}
}
\subfloat[Bumpy sphere, $N_d = 800$]
{
	\includegraphics[scale=0.5]{"Figures/3DCostsNd800".pdf}
	\label{fig:s1d}
}
\caption{Cost of node generation using Algorithm \ref{alg:node_gen_1} on 2D (top row) and 3D (bottom row) irregular domains (on a loglog plot). The timings are shown as a function of the Poisson neighborhood size $\hat{k}$, the supersampling parameter $\tau$, the number of data sites $N_d$, and the total number of nodes $N$. The dashed lines are lines of best fit used to measure slopes.}
\label{fig:scaling1}
\end{figure}

We now present scaling results (in the form of timings) for our node generator. All timings were done using a C/C++ code on the six-core Intel Coffee Lake 8700-K (12 logical cores with hyperthreading) clocked at 4.56 GHz, with 16 GB of 2600 MHz DDR4 RAM. We focus on the star domain from Section 3 and the bumpy sphere domain from~\cite{FuselierWright2013,SWFKJSC2014}. We perform two experiments: the first measures the cost of node generation in these two irregular domains, while the second measures the cost of modifying these node sets with embedded inner boundaries. The results of the first experiment are shown in Figure \ref{fig:scaling1}, and the results from the second experiment are shown in Figure \ref{fig:scaling2}.

Figure \ref{fig:scaling1} verifies that the cost of node generation is indeed approximately $O(N)$ in both 2D (Figure \ref{fig:s1a}) and 3D (Figure \ref{fig:s1b}. The cost in each case goes up slightly as $N_d$, $\hat{k}$, and $\tau$ are increased (though this is hard to see from the graphs), but the slopes remain approximately linear. Code profiling shows that the bulk of the time is taken up by kd-tree operations for large $N$. We also verified that the geometric modeling costs scaled as $O(N_b N_d)$ (not shown).  We return to this in the discussion section.
\begin{figure}[h!]
\centering
\subfloat[Star with embedded ellipse]
{
	\includegraphics[scale=0.5]{"Figures/2DModCostsNd256".pdf}

	\label{fig:s2a}
}
\subfloat[Bumpy sphere with embedded red blood cell]
{
	\includegraphics[scale=0.5]{"Figures/3DModCostsNd800".pdf}
	\label{fig:s2b}
}
\caption{Cost of node set \emph{modification} using Algorithm \ref{alg:node_mod} on 2D (left) and 3D (right) irregular domains (on a loglog plot). The number of inner boundary data nodes is fixed to $N_d^{\Gamma} = 24$ in 2D and $N_d^{\Gamma} = 200$ in 3D. The timings are shown as a function of the Poisson neighborhood size $\hat{k}$, the supersampling parameter $\tau$, and the total number of nodes $N$ (measured prior to modification). The dashed lines are lines of best fit used to measure slopes. }
\label{fig:scaling2}
\end{figure}
For the next experiment, we place a single ellipse tilted at an angle of $\pi/4$ with respect to the $x$-axis at the center of the star domain; for the 3D analogue of this experiment, we place a single red blood cell (RBC) in the interior of the bumpy sphere domain. We then measure the costs of modifying the node sets as the node spacing $h$ is decreased and $N$ is increased. It is important to note that this also has the effect of increasing the number of nodes on the boundaries of these inner embedded objects. The results are shown in Figure \ref{fig:scaling2} for different values of $\hat{k}$ and $\tau$. Both the 2D and 3D experiments (Figures \ref{fig:s2a} and \ref{fig:s2b}, respectively) show once again that our asymptotic estimates hold true: we see that the cost of node set modification is approximately $O(N)$ (see figure captions). Obviously, the precise cost depends on the shape, volume, and number of the inner embedded boundaries. Nevertheless, the experiment verifies that one can indeed modify these node sets much more quickly than one can generate them, reflecting the local nature of our node set modification algorithm. From these results, the results from Sections \ref{sec:res_hist1}--\ref{sec:res_err}, we conclude that it is reasonable to set $\tau = 2$ and $\hat{k} = 15$ for our applications. It may be reasonable to adapt $\tau$ according to the local curvature of the domain boundary. We leave this for future work.


text:
\chapter[Schoenberg-von Neumann construction of \HE]{Schoenberg-von Neumann construction of the energy space \HE}
\label{sec:Construction-of-HE}

\headerquote{If people do not believe that mathematics is simple, it is only because they do not realize how complicated life is.}{---~John~von~Neumann}

Studying the geometry of state space $X$ through vector spaces of functions on $X$ is a fundamental idea and variations of it can be traced back in several areas of mathematics. In the setting of Hilbert space, it originates with a suggestion of B. O. Koopman \cite{Koo27,Koo36,Koo57} in the early days of ``modern'' dynamical systems, ergodic theory, and the systematic study of representations of groups. A separate impetus in 1932 were the two ergodic theorems, the $L^2$ variant due to von Neumann \cite{vN32c} and the pointwise variant due to G. D. Birkhoff. While Birkhoff's version is deeper, von Neumann's version really started a whole trend: mathematical physics, quantization \cite{vN32c}, and operator theory; especially the use of the adjoint operator and the deficiency indices which we find useful in \S\ref{sec:the-boundary}--\S\ref{sec:Lap-on-HE}; cf.~\cite{vN32a,vN32b}. Further, there is an interplay between Hilbert space on the one side, and pointwise results in function theory on the other: In fact, the $L^2$-mean ergodic theorem of von Neumann is really is a corollary to the spectral theorem in its deeper version (spectral resolution via projection-valued measures) as developed in by M. H. Stone and J. von Neumann in the period 1928-1932; cf.~\cite{vN32b} and \cite[Ch. 2]{Arveson:invitation-to-Cstar}. This legacy motivates the material in this section, as well as our overall approach.

\section{Schoenberg and von Neumann's embedding theorem}
\label{sec:vonNeumann's-embedding-thm}

In Theorem~\ref{thm:R^F-embed-ERN-in-Hilbert} we show that an resistance network equipped with resistance metric may be embedded in a Hilbert space in such a way that $R$ is induced from the inner product of the Hilbert space. As a consequence, we obtain an alternative and independent construction of the Hilbert space \HE of finite-energy functions. This provides further justification for \HE as the \emph{natural} Hilbert space for studying the metric space $(\Graph,R^F) = ((\Graph,\cond),R^F)$ and \Fin as the natural Hilbert space for studying the metric space $(\Graph,R^W)$. Although we will be interested in both $(\Graph, R^F)$ and $(\Graph, R^W)$, for brevity, we sometimes refer to both as $(\Graph, R)$ when the distinction is not important.
  \glossary{name={$(\Graph,R)$},description={the network $(\Graph,\cond)$ as a metric space under $R$},sort=G}

It is a natural question to ask whether or not a metric space may be naturally represented as a Hilbert space, and von Neumann proved a general result which provides an answer. The reader may wish to consult \S\ref{sec:von-Neumann's-embedding-theorem} for the statement of this result (Theorem~\ref{thm:vonNeumann's-embedding-thm}) in the form it is applied below, as well as the relevant definitions. We apply Theorem~\ref{thm:vonNeumann's-embedding-thm} to the metric space $(\Graph,R)$ and to obtain a Hilbert space \sH and a natural embedding $(\Graph,R) \to \sH$. It turns out that the Hilbert space is \HE when the embedding is applied with $R=R^F$ and \Fin when applied with $R=R^W$; see Remark~\ref{rem:uniqueness-of-HE}! Therefore, the Hilbert space \HE is the natural place to study $(\Graph,R)$. The reader may find the references \cite{vN32a,Ber84,Ber96,Schoe38a} to be helpful; see also Theorem~\ref{thm:Schoenberg's-Thm}.

The following theorem is inspired by the work of von Neumann and Schoenberg \cite{Schoe38b,Ber84}, but is a completely new result. One aspect of this result that contrasts sharply with the classical theory is that the embedding is applied to the metric $R^{1/2}$ instead of $R$, for each of $R=R^F$ and $R=R^W$.

\begin{theorem}\label{thm:R^F-embed-ERN-in-Hilbert}
  $(\Graph,R^F)$ may be isometrically embedded in a Hilbert space.
  \begin{proof}
    According to Theorem~\ref{thm:vonNeumann's-embedding-thm}, we need only to check that $R^F$ is negative semidefinite (see Definition~\ref{def:negative-semidefinite}). Let $f:\verts \to \bR$ \saty $\sum_{x \in \verts} f(x) = 0$. We must show
that $\sum_{x,y \in F} \cj{f(x)} R^F(x,y) f(y) \leq 0$, for any finite subset $F \ci \verts$. From \eqref{eqn:def:R^F(x,y)-energy}, we have
    \linenopax
    \begin{align*}
      \sum_{x,y \in F} \cj{f(x)} R^F(x,y) f(y)
      &= \sum_{x,y \in F} \cj{f(x)} \energy(v_x - v_y) f(y) \\
      &= \sum_{x,y \in F} \cj{f(x)} \energy(v_x) f(y) - 2 \sum_{x,y \in F} \cj{f(x)} \la v_x, v_y \ra_\energy f(y) \\
          &\hstr[10]+ \sum_{x,y \in F} \cj{f(x)} \energy(v_x) f(y) \\
      &= -2\left\la \sum_{x \in F} f(x) v_x, \sum_{y \in F} f(y) v_y \right\ra_\energy \\
      &= -2\left\|\sum_{x \in F} f(x) v_x \right\|_\energy^2
      \leq 0.
    \end{align*}
    For the second equality, note that the first two sums vanish by the assumption on $f$.
  \end{proof}
\end{theorem}

\begin{cor}\label{thm:R^W-embed-ERN-in-Hilbert}
  $(\Graph,R^W)$ may be isometrically embedded in a Hilbert space.
  \begin{proof}
    Because the energy-minimizer in \eqref{eqn:def:R^W(x,y)-energy} is $f_x = \Pfin v_x$, we can repeat the proof of Theorem~\ref{thm:R^F-embed-ERN-in-Hilbert} with $f_x$ in place of $v_x$ to obtain the result.
  \end{proof}
\end{cor}

\begin{remark}\label{rem:uniqueness-of-HE}
  Since $R^F(x,y) =\|v_x-v_y\|_\energy^2$ by \eqref{eqn:def:R^F(x,y)-energy}, Theorem~\ref{thm:uniqueness-of-vNeu-embedding} shows that the embedded image of $(\Graph,R^F)$ is unitarily equivalent to the \energy-closure of $\spn\{v_x\}$, which is \HE. Similarly, $R^W(x,y) =\|f_x-f_y\|_\energy^2$, where $f_x := \Pfin v_x$, by \eqref{eqn:def:R^W(x,y)-energy}, whence the embedded image of $(\Graph,R^W)$ is unitarily equivalent to the \energy-closure of $\spn\{f_x\}$, which is \Fin.
\end{remark}

\begin{remark}\label{rem:norm-vs-quasinorm}
  One can choose any vertex $o \in \verts$ to act as the ``origin'' and it becomes the origin of the new Hilbert space during the construction outlined in \S\ref{sec:von-Neumann's-embedding-theorem}. As a quadratic form defined on the space of all functions $v:\verts \to \bC$, the energy is indefinite and hence allows one to define only a quasinorm. There are ways to deal with the fact that \energy does not ``see constant functions''. One possibility is to adjust the energy so as to obtain a true norm, as follows:
  \begin{equation}\label{eqn:energy-o}
    \energy^o(u,v) := \energy(u,v) + u(o)v(o).
  \end{equation}
  The corresponding quadratic form is immediately seen to be a norm; this approach is carried out in \cite{FOT94}, for example, and also occasionally in the work of Kigami. This is discussed in \S\ref{sec:grounded-energy-space} under the name ``grounded energy form'.

  We have instead elected to work ``modulo constants''; the kernel of \energy is the set of constant functions, and inspection of von Neumann's embedding theorem (cf.~\eqref{eqn:vNeu's-quotient-completion}) shows that it is precisely these functions which are ``modded out'' in von Neumann's construction. However, the constant functions resurface as multiples of the vacuum vector in the Fock space representation of \S\ref{sec:bdG-as-equivalence-classes-of-paths}.
\end{remark}

\begin{comment}\label{thm:unitary-vN-to-H_E}
  For $\HE = \dom \energy/\{constants\}$, we have a unitary equivalence $U:H_{vN} \to \HE$ where
  \begin{equation}\label{eqn:unitary-vN-to-H_E}
    U(w) = v \q\text{by}\q
    U:\sum_x \gh_x w_x \mapsto \sum_x \gh_x v_x,
  \end{equation}
  where for each $x \in \verts$, $\Lap v_x = \gd_x - \gd_o$ and $w_x := \tfrac1{\sqrt2}(\gd_x - \gd_o)$.
  \begin{proof}
    First, note that for each $x \in X$, the existence of $\{v_x\}$ is ensured by Theorem~\ref{thm:Pot-is-nonempty-by-current-flows}. Since $\spn\{w_x\}$ is dense in $H_{vN}$, let $
      w = \sum_{x \neq o} \gh_x w_x$
    be any finite linear combination from $\{w_x\}$ so that $\gh_o := -\sum_{x \neq o} \gh_x$ is finite. Since $w_o = 0$, we have excluded this term from the sum.

    Now we show that \eqref{eqn:unitary-vN-to-H_E} is isometry and that it is well-defined:\footnote{In general, a formula like \eqref{eqn:unitary-vN-to-H_E} may not give a well-defined function due to linear dependency.}
    \linenopax
    \begin{align*}
      \|w\|_{vN}^2
      = \left\|\sum_{x \neq o} \gh_x w_x \right\|_{vN}^2
      &= \frac12 \left\|\sum_{x \neq o} \gh_x \gd_x
        - \sum_{x \neq o} \gh_x \gd_o \right\|_{vN}^2
      = \frac12 \left\|\gh_o \gd_o + \sum_{x \neq o} \gh_x \gd_x\right\|_{vN}^2.
    \end{align*}
    A tedious application of \eqref{eqn:vNeu-inner-prod} in the form $\|f\|_{vN}^2 = -\sum_{x,y} f(x) R(x,y) f(y)$ shows that $\|w\|_{vN}^2$ is equal to
    \linenopax
    \begin{align*}
      \frac12 \sum_{x \neq o} \gh_x R(o,x) \gh_o
        + \frac12 \sum_{x,y} \gh_x R(x,y) \gh_y
      &= \frac{\gh_o}2 \sum_{x \neq o} \gh_x \energy(v_x)
        + \frac12 \sum_{x,y} \gh_x \energy(v_x-v_y) \gh_y,
    \end{align*}
    since $R(x,o) = \energy(v_x)$, and $\energy(v_x-v_y) = R(x,y)$ by \eqref{eqn:def:R(x,y)-energy}; note that $v_x \in \Pot(x,o)$ and $v_x-v_y \in \Pot(x,y)$. By expanding $\energy(v_x-v_y) = \la v_x-v_y,v_x-v_y \ra_\energy$ we finally see that $\|w\|_{vN}^2$ is equal to
    \linenopax
    \begin{align}\label{eqn:energy-as-quadratic-form-condensed}
      \sum_{x,y \neq o} \gh_x \energy(v_x,v_y)\gh_y
      = \energy\left(\sum_{x \neq o} \gh_x v_x, \sum_{y \neq o} \gh_y v_y \right)
      = \left\| \sum_{x \neq o} \gh_x v_x\right\|_\energy^2
      = \left\| v \right\|_\energy^2.
    \end{align}
    This shows we have an isometry from a dense subset of $H_{vN}$ to the span of $\{v_x\}$, which is dense in \HE by Theorem~\ref{thm:vx-dense-in-HE}.
  \end{proof}
\end{comment}

\begin{comment}
\begin{remark}\label{rem:appearance-of-Greens-oper-in-vNeu-constr}
  Note that the unitary $U$ in Theorem~\ref{thm:unitary-vN-to-H_E} is essentially $\sG = \Lap^{-1}$, the Green's operator; see \S\ref{sec:Green-function-and-operator}.
\end{remark}

\begin{cor}\label{thm:resistance-distance-as-norm-difference}
  For any $x,y \in \verts$, one has $R(o,x) = \|v_x\|_\energy^2$ and more generally
  \begin{equation}\label{eqn:resistance-distance-as-norm-difference}
    R(x,y) = \|v_x-v_y\|_\energy^2.
  \end{equation}
  \begin{proof}
    From the von Neumann construction we already had $R(x,y) = \|w_x - w_y\|_{vN}^2$, and the unitary $U$ of Theorem~\ref{thm:unitary-vN-to-H_E} is an isometry.
  \end{proof}
\end{cor}

\begin{cor}\label{thm:R(x,y)=sup<vx-vy,u>}
  The resistance distance may also be computed as
  \linenopax
  \begin{equation}\label{eqn:R(x,y)=sup<vx-vy,u>}
    R(x,y)^{1/2} = \sup_{\|u\|_\energy \leq 1} \la v_x-v_y,u\ra_\energy.
  \end{equation}
  \begin{proof}
    This follows immediately from Corollary~\ref{thm:resistance-distance-as-norm-difference}, since the inner product between two vectors always attains its maximum when the vectors are parallel:
    \linenopax
    \begin{align*}
      \sup_{\|u\|_\energy \leq 1} \la v_x-v_y,u\ra_\energy
      &= \left\la v_x-v_y,\frac{v_x-v_y}{\|v_x-v_y\|_\energy}\right\ra_\energy
      = \la v_x-v_y,u\ra_\energy = R(x,y)^{1/2}.
      &&\qedhere
    \end{align*}
  \end{proof}
\end{cor}

\begin{remark}\label{rem:sup-versions-of-R}
  It is interesting to compare \eqref{eqn:R(x,y)=sup<vx-vy,u>} to \eqref{eqn:def:R(x,y)-sup}.
\end{remark}

\begin{cor}\label{thm:vNeu-metric}
  The von Neumann metric (that is, the metric of $\sH_{vN}$) is given by $d_{vN}(x,y) = \sqrt{R(x,y)} = \|v_x-v_y\|_{vN}$.
\end{cor}

\begin{cor}\label{thm:H_energy-as-functions-of-finite-distance}
  $v \in \HE$ if and only if $\sum_{x,y \in \verts} v(x) R(x,y) v(y) < \iy$.
  \begin{proof}
    This is immediate from the preceding construction and
    \linenopax
    \begin{align*}
      \sum_{x,y \in \verts} v(x) R(x,y) v(y)
      &= \left\|\sum_{x \in \verts} v(x)(\gd_x-\gd_o)\right\|_\energy^2.
      \qedhere
    \end{align*}
  \end{proof}
\end{cor}

\begin{remark}\label{rem:Fin(verts)-not-dense-in-HE}
  Note that $\Fin(\verts)$, the finitely supported functions on \verts, may not be dense in \HE; the binary tree in Example~\ref{exm:binary-tree:nonapproximable} illustrates the reason: if $v$ is the function defined to be 1 on one branch of the tree and $-1$ on the other (and 0 at 0), then $\energy(v)=2$ but one can show $\energy(f-v) \geq \frac12$ for any finitely supported function on \verts.
  \glossary{name={$\Fin(X)$},description={the functions with finite support in $X$},sort=f,format=textbf}
\end{remark}

\begin{remark}\label{rem:define-HE-as-domE}
  One could also have simply defined $\HE := \dom \energy/\{constants\}$ from the beginning. However, the easiest way to show that \HE is complete requires the operators developed in the next section; see Lemma~\ref{thm:HE=domE-is-complete}.
\end{remark}

\begin{remark}\label{rem:Hvn-equiv-to-HE}
  Note that the uniqueness of von Neumann's construction (Theorem~\ref{thm:uniqueness-of-vNeu-embedding}) shows there is nothing special about the embedding $x \mapsto w_x$ described just above; up to unitary equivalence, \HE is the only Hilbert space into which $(\Graph,R)$ may be metrically embedded. This motivates the next section.
\end{remark}

\version{}{\marginpar{Answer this}
{\Large \textbf If you use $R^W$ instead of $R^F$, does the embedding give you $\Fin$???}
}
\end{comment}

\section{\HE as an invariant of \Graph}
\label{sec:HE-as-an-invariant}
In this section, we show that \HE may be considered as an invariant of the underlying graph.

\begin{defn}\label{def:ERN-morphism}
  Let \Graph and $H$ be resistance networks with respective conductances $\cond^{\Graph}$ and $\cond^H $. A \emph{morphism of resistance networks} is a function $\gf:(\Graph, \cond^{\Graph}) \to (H, \cond^{H})$ between the vertices of the two underlying graphs for which

  \begin{equation}\label{eqn:ERN-morphism}
    \cond^H_{\gf(x) \gf(y)} = r \cond^{\Graph}_{xy}, \qq 0 < r < \iy,
  \end{equation}
  for some fixed $r$ and all $x,y \in \verts$.

  Two resistance networks are \emph{isomorphic} if there is a bijective morphism between them.
  Note that this implies
  \begin{equation}\label{eqn:edges-in-H}
    \edges[H] = \{(\gf(x),\gf(y)) \suth (x,y) \in \edges\}.
  \end{equation}
\end{defn}

\begin{defn}\label{def:metric-space-morphism}
  A \emph{morphism of metric spaces} is a homothetic map, that is, an isometry composed with a dilation:
  \begin{equation}\label{eqn:metric-space-morphism}
    \gf:(X,d_X) \to (Y,d_Y),
    \qq d_Y(\gf(a),\gf(b)) = r d_X(a,b), \qq 0 < r < \iy,
  \end{equation}
  for some fixed $r$ and all $a,b \in X$. An \emph{isomorphism} is, of course, an invertible morphism.
\end{defn}

We allow for a scaling factor $r$ in each of the previous definitions because an isomorphism amounts to a relabeling, and rescaling is just a relabeling of lengths. More formally, an isomorphism in any category is an invertible mapping, and dilations are certainly invertible for $0 < r < \iy$.

\begin{theorem}\label{thm:R-is-a-functor}
  For each of $R=R^F,R^W$, there is a functor $\sR:(\Graph, \cond) \to ((\Graph, \cond),R)$ from the category of resistance networks to the category of metric spaces.
  \begin{proof}
    One must check that an isomorphism $\gf:(\Graph,\cond_{\Graph}) \to (H,\cond_H)$ of resistance networks induces an isomorphism of the corresponding metric spaces, so check that \gf preserves \energy. We use $x,y$ to denote vertices in \Graph and $s,t$ to denote vertices in $H$.
    \linenopax
    \begin{align}
      \la u \comp \gf, v \comp \gf\ra_\energy
      &= \sum_{x,y} \cond_{xy} (\cj{u \comp \gf}(x) - \cj{u \comp \gf}(y)) (v \comp \gf(x) - v \comp \gf(y)) \notag \\
      &= r^{-1} \sum_{x,y} \cond_{\gf(x)\gf(y)} (\cj{u(\gf(x))} - \cj{u(\gf(y))}) (v(\gf(x)) - v(\gf(y))) \notag \\
      &= r^{-1} \sum_{s,t} \cond_{st} (\cj{u}(s) - \cj{u}(t)) (v(s) - v(t)) \notag \\
      &= r^{-1} \la u \comp \gf, v \comp \gf\ra_\energy,
      \label{eqn:scaling-of-energy-under-isom}
    \end{align}
    where we can change to summing over $s,t$ because \gf is a bijection. Therefore, the reproducing kernel $\{v_x\}$ of $(\Graph,R^F)$ (or $\{\Pfin v_x\}$ of $(\Graph,R^W)$) is preserved, and hence so is the metric.
  \end{proof}
\end{theorem}

\begin{cor}\label{thm:scaling-of-Lap-under-isomorphism}
  If \gp is an isomorphism of resistance networks with scaling ratio $r$, then
  \linenopax
  \begin{align}
    \Lap(v \comp \gf) = r^{-1} \Lap(v) \comp \gf.
  \end{align}
  \begin{proof}
    Compute $\Lap(v \comp \gf)(x)$ exactly as in \eqref{eqn:scaling-of-energy-under-isom}.
  \end{proof}
\end{cor}

\begin{cor}\label{thm:isomorphic-ERNs-induce-isomorphic-metric-spaces}
  An isomorphism $\gf:(\Graph, \cond^{\Graph}) \to (H, \cond^{H})$ of resistance networks induces an isomorphism of metric spaces (where the resistance networks are equipped with their respective effective resistance metrics).
\end{cor}

We use the notation $\clspn{S}$ to denote the closure of the span of a set of vectors $S$ in a Hilbert space, where the closure is taken with respect to the norm of the Hilbert space. The following theorem is just an application of Theorem~\ref{thm:uniqueness-of-vNeu-embedding} with the quadratic form $\tilde Q = \la \cdot\,,\,\cdot\ra_\sK$.

\begin{theorem}\label{thm:unitary-invariance-of-HE}
  If there is a Hilbert space $\sK = \clspn{k_x}$ for some $k:X \to \sK$ with $d(x,y)=\|k_x-k_y\|_\sK^2$, then there is a unique unitary isomorphism $U:\sH \to \sK$ and it is given by $U:\sum_x \gx_x w_x \mapsto \sum_x \gx_x k_x$.
  \begin{comment}
    Let $\sum \gx_x = 0$. It is conceivable that $U$ fails to be well-defined because of linear dependency; we show this is not the case:
    \linenopax
    \begin{align*}
      \left\| \sum \gx_x w_x\right\|^2
      &= \sum_{x,y} \gx_x \la w_x, w_y\ra_{vN} \gx_y \\
      &= \sum_{x,y} \gx_x \gx_y (d^2(x,y) - d^2(x,o) - d^2(y,o))
        \qq \text{by \eqref{eqn:crossterm-sums-of-madness}}\\
      &= \sum_{x,y} \gx_x \gx_y d^2(x,y) - \sum_{x} \gx_x d^2(x,o) \cancel{\sum_y \gx_y} - \sum_{y} \gx_y d^2(y,o) \cancel{\sum_x \gx_x}.
    \end{align*}
    But $\left\|\sum \gx_x k_x\right\|^2 = \sum_{x,y} \gx_x \gx_y d^2(x,y)$, too, by the same argument. This shows the existence of $U$, as well as the fact that it is well-defined.
  \end{comment}
\end{theorem}

\begin{remark}\label{rem:HE-is-an-invariant}
  Theorem~\ref{thm:unitary-invariance-of-HE} may be interpreted as the statement that there is a functor from the category of metric spaces (with negative semidefinite metrics) into the category of Hilbert spaces. However, we have avoided this formulation because the functor is not defined for the entire category of metric spaces.
  For us, it suffices to note that the composition is a functor from resistance networks to Hilbert spaces, so that $\HE = \HE(\Graph)$ is an invariant of \Graph.
\end{remark}

\begin{remark}\label{rem:first-quantization}
  To obtain a first quantization, one would need to prove that a contractive morphism between resistance networks induces a contraction between the corresponding Hilbert spaces. In other words,
  \linenopax
  \begin{align*}
    f:\Graph_1 \to \Graph_2
    \q\implies\q
    T_f:\HE(\Graph_1) \to \HE(\Graph_2)
  \end{align*}
  with $\|T_f v\|_\energy \leq \|v\|_\energy$ whenever $f$ is contractive. The authors are currently working on this endeavour in \cite{JoPea10c}. The second quantization is discussed in Remark~\ref{rem:Wiener-improves-Minlos}.
\end{remark}

\version{}{\marginpar{Add stuff about rough isometry \& rough equivalence, someday.}}

\section{Remarks and references}
\label{sec:Remarks-and-References-construction-of-HE}

Of the results in the literature of relevant to the present chapter, the references \cite{Ba04, Ale75, PaSc72, Schoe38a, Schoe38b} are especially relevant. Some of the cited references for this chapter are more specialized, but for prerequisite material (if needed), the reader may find the books \cite{Gui72} by Guichardet, \cite{Hida80} by Hida, and \cite{PaSc72} by Parthasarathy \& Schmidt especially relevant. Standard applications of a negative definite function include either the construction of an abstract Hilbert space \cite{vN55,vN32c} or else the construction of measures on a path space \cite{PaSc72,Minlos63}.

We use the terminology ``Schoenberg-von Neumann embedding'' to denote a set of general principles, both classical and modern:
\linenopax
\begin{align*}
  \xymatrix{
  \text{geometry and physics} \ar[d]^{\text{\cite{Koo57}}} \\
  \text{Hilbert space} \ar[d]^{\text{\cite{Nel73a}}} \\
  \text{$L^2$ space of random variables}
  }
\end{align*}

Some examples of the Schoenberg-von Neumann embedding include: Brownian motion \cite{Nelson64, Nel69}, second quantization and quantum fields \cite{Minlos63, Gro70, Hida80}, stochastic integrals \cite{Malliavin}, spin models \cite{Lig93}, quantum spin models \cite{Pow67, Pow76a, Pow76b}.

See Chapter~\ref{sec:Magnetism-and-long-range-order} below for further details, especially Theorem~\ref{thm:Powers'-correlation-estimate}.

config: arxiv


TITLE: Show that $\int_a^b xdx = \frac{b^2-a^2}{2}$ using upper and lower Riemann sums
QUESTION [1 upvotes]: Here's my work so far:
Let $f(x) = x$ and suppose that $P = \{x_0, x_1, x_2,...,x_n\}$ is the partition for interval $[a,b]$, creating $n$ number of subintervals with equal lengths.
$L_f(P)$ denotes the lower Riemann sum associated with partition P and $U_f(P)$ denotes the upper Riemann sum associated with partition P. So from what I know, $U_f(P)$ will always be greater than or equal to the true definite integral. Likewise, $L_f(P)$ will always be less than or equal to the true definite integral. Like so...
$$L_f(P) \leq \int_a^b xdx \leq U_f(P)$$
To calculate $L_f(P)$, for each integer $i$ in $[1,n]$, let $m_i = \textrm{min}\{f(x) \,|\, x \  \textrm{in} \ [x_{i-1}, x_i]\}$. Then...
\begin{align}
L_f(P) &= \sum_{i=1}^n m_i(x_i - x_{i-1}) = (\frac{b-a}{n})(a+ ...+m_{n})
\end{align}
how can I continue? My goal is to show the upper and lower sums are equal.

REPLY [2 votes]: You are close to being done. Note that $x_i = a + \frac{i}{n}(b-a)$, so
\begin{align*}
\sum_{i=1}^{n}m_i &= \sum_{i=0}^{n-1}(a+ \frac{i}{n}(b-a)) \\
&= na + \frac{b-a}{n}\sum_{i=0}^{n-1} i \\
&= na + \frac{b-a}{n}\frac{n(n-1)}{2}
\end{align*}
The lower Riemann sum then becomes
\begin{align*}
L_f(P)=\frac{b-a}{n}\cdot\sum_{i=1}^nm_i &= a(b-a) + \frac{n-1}{n} \frac{(b-a)^2}{2}
\end{align*}
A very similar computation gives us that
$$U_f(P) = a(b-a)+\frac{n+1}{n} \frac{(b-a)^2}{2}.$$
Letting $n$ approach $\infty$ we can see that we get the same limit namely
\begin{align*}
a(b-a) + \frac{(b-a)^2}{2} &= \frac{2ab - 2a^2 + b^2 + a^2 - 2ab}{2} = \frac{b^2 - a^2}{2}
\end{align*}


TITLE: Identifying closed form of two-parameter sequence
QUESTION [2 upvotes]: I've come across a two-parameter sequence $a_{nk}$ with $n=1,2,\ldots$ and $k=1,2,\ldots,n$, and I would like to identify a closed expression for it. So far I have the elements
$a_{1k}=\{1\}\\ a_{2k}=\{1,3\}\\ a_{3k}=\{1,8,12\}\\ a_{4k}=\{1,15,60,60\}\\ a_{5k}=\{1,24,180,480,360\}\\ a_{6k}=\{1,35,420,2100,4200,2520\}\\ a_{7k}=\{1,48,840,6720,25200,40320,20160\}\ ,$
with $k$ increasing to the right.
There seems to be a lot of order in the elements, and I think there's probably a not too complicated expression for it - but I haven't been able to find one, despite trying the integer sequence encyclopedia on various subsets and other ways to look for a system.
So far I've identified the following:
Looking at the second $k$ for each $n>1$ they can be described with a simple formula giving $\{3,8,15,24,35,48\}=\{1\times3,2\times4,3\times5,4\times6,5\times7,6\times8\}$.
Also each row $n$ is divisible by $n-1$ except for $n=1$ and $k=1$. Performing this division gives (neglecting all $a_{nk}$ for $k=1$)
$a_{2k}=\{3\}\\ a_{3k}=\{4,6\}\\ a_{4k}=\{5,20,20\}\\ a_{5k}=\{6,45,120,90\}\\ a_{6k}=\{7,84,420,840,504\}\\ a_{7k}=\{8,140,1120,4200,6720,3360\}\ .$
Now the previously mentioned order in $k=2$ is obvious.
But there I'm stuck. Any ideas?
UPDATE: I've realised the elements have a factor ${n+1 \choose k-1}$. Dividing out this factor leaves the elements
1
1 1
1 2 2
1 3 6 6
1 4 12 24 24
1 5 20 60 120 120
1 6 30 120 360 720 720 ,
which in itself looks pretty binomial-ish..

REPLY [1 votes]: Hint: You can find the numbers in your updated section as A008279 in OEIS.


TITLE: Tensor with two indices
QUESTION [0 upvotes]: A tensor with two indices can be represented by a $3\times3$ matrix.
\begin{equation}
A= \left( \begin{array}{ccc}
A_{11} & A_{12} & A_{13} \\
A_{21} & A_{22} & A_{23} \\
A_{31} & A_{32} & A_{33} \end{array} \right)
\end{equation}
It make sense why is a $2$D matrix but why it's a $3\times3$ and not a $2\times2$, or $4\times4$?

REPLY [1 votes]: A rank-$n$ tensor is a linear map from a sequence of vector spaces to the reals(/complexes, but let's keep things simple). So it eats vectors, which are normally all in $d$ dimensions. Therefore a rank-$2$ tensor can be written as a two-index object $a_{ij}$, which acts linearly on $d$-dimensional vectors $v_i,w_j$ as
$$ \sum_{i=1}^d \sum_{j=1}^d a_{ij} v_i w_j. $$
But this is the same as the summation from the matrix product $v^T A w$, so $A$ is a $d \times d$ matrix.