\section{Introduction and Background}

Let $H^2$ denote the classical Hardy--Hilbert space on the unit disk $\mathbb{D}$, the space of analytic functions on $\mathbb{D}$ that have square-summable Maclaurin coefficients. In 1949, Arne Beurling \cite{beurling1949two} proved several seminal results regarding function and operator theory in $H^2$. One result showed that every function in $H^2$ may be factored as $\theta U$, where $|\theta| =1$ almost everywhere on the unit circle $\mathbb{T}$ (coined as \textit{inner}), and $U$ is such that $\log|U(0)| = \int_0^{2\pi} \log|U(e^{i\theta})|$ (coined \textit{outer}). This result allowed him to characterize shift-invariant subspaces of $H^2$; those subspaces $M \subseteq H^2$ with the property that $SM \subseteq M$, where $S: H^2 \to H^2$ is given by $f(z) \mapsto zf(z)$. Namely, these subspaces must equal $\theta H^2$, for some inner function $\theta$.

The simplest $H^2$-inner functions are called \textit{finite Blaschke products}, given for $\beta_1, \ldots, \beta_n \in \mathbb{D}$ as
\[
B(z) = \lambda \prod_{j = 1}^n \frac{z - \beta_j}{1 - \overline{\beta_j}z},
\]
where $\lambda \in \mathbb{T}$.
One may check that $|B| = 1$ on $\mathbb{T}$, and so $B$ is in fact inner. We recall that \textit{Blaschke factors}, given when $n=1$ in $B$ above, define the automorphisms of $\mathbb{D}$, up to a unimodular constant. See the recent book \cite{GarciaSC2018} for a nice treatment of finite Blaschke products and their applications.

Let us start with an observation. Applying a partial fractions decomposition, and assuming simple zeros different from the origin, one can check that any finite Blaschke product can be expressed as
\[
B(z) = c_0 - \sum_{j = 1}^n \frac{c_j}{1 - \overline{\beta_j}z}
\]
for some constants $c_0, c_1, \ldots, c_n \in \mathbb{C}$.
Further, each term in the sum can be seen as a scalar multiple of the Szeg\H{o} kernel,  $s_{\lambda}(z) = 1/(1-\overline{\lambda}z)$, $\lambda \in \mathbb{D}$.
Noting that $s_0(z) = 1$, we have
\[
B(z) = c_0s_0(z) - \sum_{j = 1}^n c_j s_{\beta_j}(z).
\]
Further, the Szeg\H{o} kernel is the \textit{reproducing kernel} for $H^2$. Namely, for every $\lambda \in \mathbb{D}$ and every $f \in H^2$, $\langle f, s_\lambda \rangle_{H^2} = f(\lambda)$. Consequently, every Blaschke product with simple zeros can be seen as a linear combination of reproducing kernels. If $B$ also has repeated zeros, certain derivatives of kernel functions will also be needed in the linear combination. Nonetheless, it turns out this observation actually characterizes finite Blaschke products among inner functions: an $H^2$-inner function $f$ is a finite Blaschke product if and only if $f$ is a linear combination of reproducing kernels and their derivatives.

We will prove this result in a more general setting in Theorem \ref{modified_le}. Toward this generalization, we note that the Hardy space $H^2$ is just one example of a function space where point evaluation can be recovered via its inner product. Spaces of this type are known as \textit{reproducing kernel Hilbert spaces} (RKHS).

The aim of this paper is to show, working in a general class of RKHSs, and using a generalized definition of inner, that there are analogous inner functions characterized as certain linear combinations of reproducing kernels. We call functions characterized in this way \textit{analogues of finite Blaschke products}. Further, we show precisely how these functions arise as certain Gram determinants, or as certain projections onto shift-invariant subspaces generated by polynomials.

Throughout this paper, $\mathcal{H}$ will denote a reproducing kernel Hilbert space of analytic functions on some planar domain $\Omega \subset \mathbb{C}$ with $0 \in \Omega$. Formally, $\mathcal{H}$ is a complete inner product space comprised of functions mapping $\Omega \to \mathbb{C}$ that has the reproducing kernel property, i.e. for each $w \in \Omega$, the linear functional given by $f \mapsto f(w)$ is bounded. Equivalently, by the Riesz representation theorem, for each $w \in \Omega$, there exists a unique function $k_w \in \mathcal{H}$ such that for all $f \in \mathcal{H}$, $\langle f, k_w \rangle = f(w)$. We also ask that $\mathcal{H}$ has the following properties:
\begin{enumerate}
\item The forward shift operator $S$, given by $f(z) \mapsto zf(z)$, is bounded on $\mathcal{H}$.
\item The analytic polynomials $\mathcal{P}$ form a dense subset of $\mathcal{H}$.
\end{enumerate}

The following notation will be used throughout the paper:
\begin{itemize}
    \item When $V \subseteq \mathcal{H}$ is a closed subspace, we will use $\Pi_V: \mathcal{H} \to V$ to denote the orthogonal projection from $\mathcal{H}$ onto $V$.
    \item When $X \subseteq \mathcal{H}$ is a subset, we will use $\operatorname{cl}_{\mathcal{H}}(X)$ to denote the closure of $X$ in $\mathcal{H}$.
    \item For $f \in \mathcal{H}$, we will use the standard notation $\operatorname{ord}_0(f)$ to denote the order of the zero of $f$ at the origin.
    \item For $f \in \mathcal{H}$, we will let $Z(f)$ be the multiset containing the zeros $f$, i.e. $Z(f)$ is the zero set of $f$, each zero listed with its multiplicity.
\end{itemize}

The primary examples of such spaces we keep in mind are those where $\Omega = \mathbb{D}$ and $\mathcal{H}$ is a subset of the analytic functions on $\mathbb{D}$. A well-known family of such spaces are the so-called weighted Hardy spaces. Given a sequence of positive numbers $w = \{w_k\}_{k \ge 0}$, with $\lim_{k \to \infty} w_k/w_{k+1} = 1$, define the weighted Hardy space as
\[
H^2_w : = \left\{ f(z) = \sum_{k \ge 0} a_k z^k \in \textrm{Hol}(\mathbb{D}) : \sum_{k \ge 0} w_k |a_k|^2 < \infty \right\}.
\]
For $f(z) = \sum_{k \ge 0} a_k z^k$ and $g(z) = \sum_{k \ge 0} b_k z^k$ in $H^2_w$, their inner product in $H^2_w$ is given by
\[
\langle f, g \rangle_w = \sum_{k\ge 0}w_k a_k \overline{b_k}.
\]
One may also verify (e.g. see \cite[Section~2.1]{Cowen1995}) that these spaces are reproducing kernel Hilbert spaces with reproducing kernel
\[
k^w_\beta(z) = \sum_{n\ge0} \frac{1}{w_k}\left(\bar{\beta} z\right)^n
\]
for, a priori, $\beta, z \in \mathbb{D}$. It turns out that under appropriate conditions on the weight sequence, the above kernel function extends to points on the unit circle as well.

If we let $\alpha \in \mathbb{R}$ and $w = \{(k+1)^\alpha\}_{k\ge 0}$, we recover the so-called Dirichlet-type space $\mathcal{D}_\alpha$. Of these spaces, $\alpha = -1$ coincides with the Bergman space $A^2$, $\alpha = 0$ corresponds to the Hardy space $H^2$, and $\alpha = 1$ gives the Dirichlet space $\mathcal{D}$.

\section{Inner functions}
Although analogous definitions of Beurling's inner and outer functions in spaces other than $H^2$ have been made, the objects they describe are not as well understood.
We give these definitions here.

\begin{definition}[Cyclic function]
Say that $f \in \mathcal{H}$ is cyclic (for $S$ in $\mathcal{H}$) if
\[
[f] : = \operatorname{cl}_\mathcal{H}\left(\textrm{span}\{ z^kf : k \ge 0 \}\right)
\]
is equal to all of $\mathcal{H}$.
\end{definition}
Note that $[f]$ is always a shift-invariant subspace, i.e. $S[f] \subseteq [f]$, and is the smallest such containing $f$. The space $[f]$ is read ``bracket $f$" or ``the (shift-invariant) subspace generated by $f$." In $H^2$, it is well known that cyclic functions are precisely the outer functions.

Recall that a bounded analytic function $f$ on $\mathbb{D}$ is called an inner function if $|f(\zeta)|=1$ for a.e. $|\zeta|=1$. Inner functions play an important role in function theory and operator theory on Hardy spaces. See \cite{ChalendarNWEJM2015} for a recent survey of classical and new results linking inner functions and operator theory. Besides the above definition, inner functions can also be characterized via the inner product in $H^2$. Indeed, it can be checked that a function $f\in H^2$ is inner if and only if $\|f\|_{H^2}=1$ and $\inner{z^mf}{f}=0$ for all integers $m\geq 1$.

In the case of the Dirichlet space $\mathcal{D}$, Richter \cite{RichterJRAM1988} showed that any shift-invariant subspace is also generated by a single function that satisfies the same orthogonality properties above.
 Aleman, Richter and Sundberg \cite{aleman1996beurling} proved an analogue of Beurling's Theorem for the Bergman space $A^2$: any invariant subspace $\mathcal{M}$ of $A^2$ is generated by the so-called wandering subspace $\mathcal{M}\ominus z\mathcal{M}$. Any unit norm function in this subspace satisfies $\|f\|_{A^2}=1$ and $z^mf\perp f$ for all $m\geq 1$ and is called an $A^2$-inner function. Prior to this work, Hedenmalm \cite{HedenmalmJRAM1991} showed the existence of so-called contractive zero-divisors, which play the role of Blaschke products in the Bergman space. In certain cases, explicit formulas for these functions have been given, e.g. see MacGregor and Stessin \cite{macgregor1994weighted} and Hansbo \cite{hansbo1998reproducing}. These results are phrased in the language of \textit{extremal functions}. Although our work here will not explicitly cover this aspect, it is well known that (normalized) inner functions are solutions to the extremal problem
\[
\sup\left\{ \operatorname{Re}(g^{(d)}(0)) : g \in \mathcal{M}, \|g\| \le 1 \right\},
\]
where $\mathcal{M}$ is a shift-invariant subspace and $d$ is the smallest integer so that $z^d \notin \mathcal{M}^\perp$.
 See \cite[Chapters~5 and 9]{DurenAMS2004} and \cite[Chapter~3]{HedenmalmSpringer2000} for a detailed discussion of inner functions on Bergman spaces $A^p$.
Researchers have thus defined the notion of inner functions in more general reproducing kernel Hilbert spaces.

\begin{definition}[Inner function]
Say that $f \in \mathcal{H} \setminus \{0\}$ is $\mathcal{H}$-inner if, for all $k \ge 1$,
\[
\langle f, z^k f \rangle = 0.
\]
\end{definition}

This definition was originally considered in \cite{aleman1996beurling} for the Bergman space. The authors there also require an inner function to be of unit norm, as well as other authors. In a recent paper, Cheng, Mashreghi and Ross \cite{cheng2019inner} introduced and studied the notion of inner functions with respect to a bounded linear operator. However, they do not require unit norm, which turns out more convenient in several situations. We will follow their approach in this work.
Although no function-theoretic description of inner functions is known in general reproducing Hilbert spaces, there are known constructions of certain types of inner functions. We will introduce one of these constructions at the end of Section \ref{RepPoints}.

B\'en\'eteau et al. \cite{BeneteauJLMS2016, BeneteauCMB2018} studied inner functions and examined the connections between them and optimal polynomial approximants on weighted Hardy spaces. They also described a method to construct inner functions that are analogues of finite Blaschke products with simple zeroes. In \cite{SecoCAOT2019}, Seco discussed inner functions on Dirichlet-type spaces and characterized such functions as those whose norm and multiplier norm are equal.
In \cite{LeAMP2020}, the second author studied inner functions on weighted Hardy spaces and obtained generalizations of several results from \cite{BeneteauCMB2018, SecoCAOT2019}.
In a recent paper \cite{BeneteauAMP2019}, B\'en\'eteau et al. investigated inner functions on general simply connected domains in the complex plane. We would like to mention that operator-valued inner functions on vector-valued weighted Hardy spaces have also been defined and studied \cite{BallIEOT2013, BallASM2013, OlofssonAiA2007}. In particular, Ball and Bolotnikov \cite{BallASM2013} obtained a realization of inner functions on vector-valued weighted Hardy spaces. In \cite{BallCM2016}, they investigated the expansive multiplier property of inner functions. They obtained a sufficient condition on the weight sequences for which any inner function has the expansive multiplier property. Recently, Cheng, Mashreghi, and Ross considered inner vectors for Toeplitz operators \cite{MR4061942} and for the shift operator in the Banach space setting of $\ell^p_A$ \cite{MR3976584}.

We point out here that if $f$ is inner, then $\langle pf, f \rangle = p(0)\|f\|^2$ for all $p \in \mathcal{P}$, and equivalently, up to a constant, $f$ is the reproducing kernel at the origin for $[f]$ (we ask for $0 \in \Omega$ for this kernel to make sense). The converse is also true, and we record this result below.

\begin{prop}
Let $f \in \mathcal{H}$ and suppose $d = \operatorname{ord}_0(f)$. Then $f$ is $\mathcal{H}$-inner if and only if $f$ is a non-zero constant multiple of $\Pi_{[f]}(k_0^{(d)})$.
Moreover, $\Pi_{[f]}(k_0^{(d)})$ is always an $\mathcal{H}$-inner function.
\end{prop}
\begin{proof}
Consider the forward implication. For all analytic polynomials $p$, we have
\begin{align*}
    \left\langle pf, \overline{f^{(d)}(0)}f -  \Pi_{[f]}(k_0^{(d)}) \right\rangle
    &= p(0)f^{(d)}(0) - \sum_{k = 0}^d {d \choose k}p^{(k)}(0) f^{(d - k)}(0) \\
    & = p(0)f^{(d)}(0) - p(0)f^{(d)}(0) = 0.
\end{align*}
Thus, $\overline{f^{(d)}(0)}f =  \Pi_{[f]}(k_0^{(d_f)})$.

Conversely, suppose $f = c \Pi_{[f]}(k_0^{(d)})$ for some non-zero constant $c$. Then, for $k \ge 1$, we have
\begin{align*}
    \langle z^kf, f \rangle &= \langle z^kf,  c \Pi_{[f]}(k_0^{(d)})\rangle = \overline{c} \langle z^kf, k_0^{(d)}\rangle = 0
\end{align*}
since $z^kf(z)$ vanishes at the origin with order at least $d+1$.
Thus, $f$ is $\mathcal{H}$-inner. Further, this shows that $\Pi_{[f]}(k_0^{(d)})$ is alway $\mathcal{H}$-inner.
\end{proof}

In \cite{cheng2019inner}, the authors conducted a robust exploration of inner functions. It was shown there \cite[Proposition~3.1]{cheng2019inner} that every inner function is given by
$\Pi_{[Sf]^{\perp}}(f)$. We show here that, up to a constant, this function is the same as a projection of a kernel onto a shift invariant subspace.
\begin{prop}
Let $f \in \mathcal{H}$ and let $d = \operatorname{ord}_0(f)$. Put \[J = f - \Pi_{[Sf]}(f) = \Pi_{[Sf]^{\perp}}(f)\] and $v=\Pi_{[f]}(k^{(d)}_0)$.
Then we have
\[v = \frac{v^{(d)}(0)}{f^{(d)}(0)}J.\]
\end{prop}

\begin{proof}
Note that any element of $[Sf]$ vanishes at the origin with multiplicities at least $d+1$. It follows that $k^{(d)}_0\perp [Sf]$ and hence $v\perp [Sf]$, which implies that $v\in [f]\ominus [Sf]$. On the other hand, we also have $J\in [f]\ominus[Sf]$. Since $[f]\ominus[Sf]$ is a one dimensional space, we conclude that $v = \lambda J$ for some constant $\lambda$. To find the constant $\lambda$, let us compute the inner product
\[J^{(d)}(0) = \inner{J}{k^{(d)}_0} = \inner{f}{k^{(d)}_0} - \inner{\Pi_{[Sf]}(f)}{k^{(d)}_0} = f^{(d)}(0)\]
because $k_0^{(d)}\perp [Sf]$. It then follows that
\[\lambda = \frac{\inner{v}{k^{(d)}_0}}{\inner{J}{k^{(d)}_0}} = \frac{v^{(d)}(0)}{f^{(d)}(0)}\]
and the conclusion follows.
\end{proof}

In \cite{shapiro1962zeros}, Shapiro and Shields used Gram determinants to produce linear combinations of reproducing kernels that are Dirichlet-inner functions. This was then generalized by B\'en\'eteau et al. in \cite{BeneteauCMB2018}, and further by the second author \cite{LeAMP2020} to weighted Hardy spaces over the unit disk. Surprisingly, as we will see later, even in general RKHSs in which monomials are not necessarily orthogonal, such a construction (Definition \ref{SS-function}) is the only way to produce inner functions that are linear combinations of kernels (see Theorem \ref{modified_le}).

\subsection{Gram Determinants}\label{Gram}
Let $v_1, \ldots, v_n$ be vectors in an inner product space. We will denote the associated \textit{Gram matrix} by
\[
G(v_1, \ldots, v_n) = (\langle v_i, v_j\rangle)_{1 \le i, j \le n} =
\begin{bmatrix}
\langle v_1, v_1\rangle & \dots & \langle v_1, v_n\rangle\\
\vdots & \ddots & \vdots \\
\langle v_n, v_1\rangle & \dots & \langle v_n, v_n\rangle
\end{bmatrix}.
\]
The \textit{Gram determinant} is then $\det(G(v_1, \ldots, v_n))$.
Note that for any $x = (x_1, \ldots, x_n) \in \mathbb{C}^n$,
\[
\langle G(v_1, \ldots, v_n)x, x \rangle = \Big\langle \sum_{i = 1}^n x_i v_i, \sum_{j = 1}^n x_j v_j \Big\rangle \ge 0.
\]
Hence, every Gram matrix is positive semidefinite.
Moreover, the vectors $\{ v_1, \ldots v_n\}$ are linearly independent if and only if $G(v_1, \ldots, v_n)$ has full rank, and equivalently, if and only if $\det(G(v_1, \ldots, v_n)) > 0$.

Similarly, for any vector $u$, we define
\[
D(u; v_1, \ldots, v_n) := \det
\begin{bmatrix}
u & \langle u, v_1\rangle & \dots & \langle u, v_n\rangle\\
v_1 & \langle v_1, v_1\rangle & \dots & \langle v_1, v_n\rangle\\
\vdots & \vdots & \ddots & \vdots \\
v_n &\langle v_n, v_1\rangle & \dots & \langle v_n, v_n\rangle
\end{bmatrix}.
\]
Note that $D(u; v_1, \ldots, v_n)$ is a linear combination of $u, v_1, \ldots, v_n$.

For our purposes, it is critical to note that $D(u; v_1, \ldots, v_n)$ is orthogonal to each $v_j, 1 \le j \le n$. When $\delta := \det(G(v_1, \ldots, v_n)) > 0$, we have that
\[
u - \delta^{-1} D(u; v_1, \ldots, v_n)
\]
is a linear combination of the vectors $\{v_1, \ldots, v_n \}$, since the coefficient of $u$ in $D(u; v_1, \ldots, v_n)$ is $\delta$.
Consequently, we have the following well-known lemma.
\begin{lemma}\label{SS-proj}
Let $v_1, \ldots, v_n$ be linearly independent vectors in an inner product space $\mathcal{V}$. Then for any $u\in\mathcal{V}$,
\[
\frac{D(u; v_1, \ldots, v_n)}{\det(G(v_1, \ldots, v_n))}
\]
is the orthogonal projection of $u$ onto $(\operatorname{span}\{v_1, \ldots, v_n \})^\perp$.
\end{lemma}

We note here that for any vectors $u, v_1,\ldots, v_n$ with $n\geq 2$, the coefficient of $v_{n}$ in $D(u; v_1,\ldots, v_n)$ is exactly
$-\big\langle D(u; v_1,\ldots, v_{n-1}), v_n\big\rangle$.

For a set of distinct points $\beta_1, \ldots, \beta_n \subset \mathbb{D}$, the authors in \cite{BeneteauCMB2018} coined the \textit{Shapiro--Shields function} as the normalization of $D(1; k_{\beta_1}, \ldots, k_{\beta_n})$, where the inner product is taken in some $H^2_w$. This follows the construction of Shapiro and Shields in \cite{shapiro1962zeros}. They also showed that this function is always inner.

\section{Reproducible Points}\label{RepPoints}
It is important to note that the spaces in which we are working are RKHSs on $\Omega$ but may also have kernels that extend to points outside of $\Omega$. For example, the linear functional of point evaluation extends boundedly to the unit circle in the Dirichlet-type spaces when $\alpha >1$ (more on this below in Example \ref{D-kernels}).

It is our aim to uncover the relationship between $\mathcal{H}$-inner functions, linear combinations of reproducing kernels, generalized Shapiro--Shields functions, and projections of kernels onto shift-invariant subspaces generated by polynomials. If $\beta \notin \Omega$, but point evaluation is bounded at $\beta$, we can still make sense of, for example, the function $D(u; k_\beta)$ when $u \in \mathcal{H}$.
We introduce here a framework to handle such points.

\begin{definition}[Reproducible point/order]
Say $\beta \in \mathbb{C}$ is a reproducible point of order $m$ in $\mathcal{H}$ if the linear functional
\[
p \mapsto p^{(m)}(\beta)
\]
extends from $\mathcal{P}$ to a bounded linear functional on $\mathcal{H}$.
\end{definition}

It is evident that any $\beta\in\Omega$ is reproducible of order $m$ for all $m\geq 1$. On the other hand, points outside of $\Omega$ may only be reproducible up to a certain finite order. First we establish a simple fact about the orders of a reproducible point.

\begin{lemma}\label{reprod}
If $\beta$ is reproducible of order $m$, then it is also reproducible of all orders $0\leq j\leq m$.
\end{lemma}

\begin{proof}
By assumption, there exists a constant $C_m>0$ such that for all polynomials $q$, one has
\[
|q^{(m)}(\beta)| \leq C_m\|q\|.
\]
For any $0\leq j\leq m$, since multiplication by $z$ is bounded, there exists a constant $B_{j}>0$ such that $\|(z-\alpha)^{m-j}p\|\leq B_{j}\|p\|$ for all polynomials $p$. On the other hand,
\[
\binom{m}{j}\, p^{(j)}(\beta) = \frac{d^{m}}{dz^{m}}\left((z-\beta)^{m-j}p(z)\right)\Big\rvert_{z=\beta}.
\]
It then follows that
\[
|p^{(j)}(\beta)| \leq \frac{C_{m}B_{j}}{\binom{m}{j}}\|p\|
\]
Therefore the map $p\mapsto p^{(j)}(\beta)$ extends to a bounded linear functional on $\mathcal{H}$.
\end{proof}

\begin{definition}
Let $\beta \in \mathbb{C}$ be reproducible of some order. Define the reproducible order of $\beta$ (in $\mathcal{H}$) as
\[
\operatorname{ro}(\beta) := \sup_m\{ p \mapsto p^{(m)}(\beta) \textrm{ extends as a bounded functional from $\mathcal{P}$ to $\mathcal{H}$}\}.
\]
\end{definition}

\begin{example}[$D_\alpha$, $\alpha >1$]\label{D-kernels}
Consider the Dirichlet-type spaces $\mathcal{D}_\alpha$. Recall that monomials have norm $\|z^n\|_\alpha = (n+1)^{\alpha/2}$ and that for $|\beta|<1$,
\[
k_\beta(z) = \sum_{n\ge 0} \frac{\bar{\beta}^n}{(n+1)^{\alpha/2}}\frac{z^n}{(n+1)^{\alpha/2}} = \sum_{n\ge 0}\frac{(\bar{\beta}z)^n}{(n+1)^{\alpha}}.
\]
For $\alpha>1$, it is evident that $k_{\beta}$ is a function in $D_{\alpha}$ even for $|\beta|=1$, which implies that all points on the unit circle are reproducible points.

In addition, it is well known that for $|\beta|<1$, the linear functional given by $f \mapsto f^{(m)}(\beta)$ has a reproducing kernel given by
\[
k_\beta^{(m)}(z) = \frac{\partial^m}{\partial \bar{\beta}^m} k_\beta(z) = \sum_{n\geq m}\frac{n(n-1)\cdots (n-m+1)}{(n+1)^{\alpha}}\bar{\beta}^{n-m}\,z^{n}.
\]
The above series expansion shows that $k_{\beta}^{(m)}$ belongs to $D_{\alpha}$ for some (and hence all) $|\beta|=1$ if and only if $\alpha > 2m+1$. As a consequence, all points on the unit circle are reproducible of order $r$ in $D_{\alpha}$, where $r$ is the largest natural number strictly less than $\frac{\alpha-1}{2}$.
\end{example}

\begin{example}[Local Dirichlet spaces]
\label{Ex:localD}
Let $\zeta \in \mathbb{T}$ and let $\delta_\zeta$ denote the Dirac measure at $\zeta$. One may form the \textit{local Dirichlet space at $\zeta$}
\[
\mathcal{D}_{\delta_\zeta}:= \left\{ f \in \textrm{Hol}(\mathbb{D}) : \int_{\mathbb{D}}|f'(z)|^2 \frac{1 - |z|^2}{|z- \zeta|^2}dA(z) < \infty \right\}
\]
where $dA$ is normalized area measure on $\mathbb{D}$. It is well known that the polynomials are dense in $\mathcal{D}_{\delta_\zeta}$ and that it is a reproducing kernel Hilbert space on $\mathbb{D}$. Additionally, $\mathcal{D}_{\delta_\zeta}$ has the property that point evaluation is bounded at $\zeta \in \mathbb{T}$  but not at any other point on $\mathbb{T}$ (see \cite[Theorems 7.2.1 and 8.1.2 (ii)]{el2014primer}).

\end{example}

\begin{example}[$L$ Regions]
Let $\Delta$ be an infinite sequence of disjoint closed discs whose centers lie on the positive real axis and decrease monotonically to zero. By deleting $\Delta$ from $\mathbb{D}$, one obtains an infinitely connected region, known as an $L$ region (see \cite{zalcman1969bounded, mccarthy1994bounded}).
The origin is a boundary point of the region, and in \cite{mccarthy1994bounded} it was shown that for certain reproducing kernel Hilbert spaces of analytic functions on the region, where the polynomials are dense, point evaluation is bounded at the origin, dependent on the rate of decay of the radii of the disks in $\Delta$. Uniformly perturbing the disks of $\Delta$ to the right by a fixed positive amount $\epsilon > 0$, we obtain an infinitely connected $L$ like region containing zero, where the results describing bounded point evaluation at the origin hold then for the boundary point $\epsilon$.
We communicate this example to highlight that we need no hypotheses on the connectedness of $\Omega$ and that there is interesting reproducible behavior in this case.
\end{example}

The previous examples show that spaces with bounded point evaluation can behave very differently, and point evaluation can extend outside of the domain $\Omega$ in various ways.

As we will see, Shapiro--Shields functions can be viewed as projections of a kernel at zero onto certain shift invariant subspaces generated by polynomials. Consequently, we would like to connect reproducibility with the zeros of a polynomial.
\begin{definition}[Reproducible zeros]
Let $p \in \mathcal{P}$. The multiset of reproducible zeros of $p$ (in $\mathcal{H}$) is
\[
R(p) : = \{ \beta \in \mathbb{C} : p^{(m)}(\beta) = 0 \textrm{ and $\beta$ is reproducible of order $m$ in $\mathcal{H}$} \},
\]
listed with multiplicity.
\end{definition}
Namely, by multiset and ``listed with multiplicity,'' we require that if $p$ has a zero of order $m$ at $\beta$, then $\beta$ appears in $R(p)$ with multiplicity $\min\{\operatorname{ord}(\beta), m\}$. For example, if the point 1 is reproducible of order 2, but not of order 3 (i.e. $\operatorname{ro}(1) = 2$), and the point $\pi i$ is not reproducible in $\mathcal{H}$, then for $p(z) = z(z - 1)^3(z - \pi i)$, we have $R(p) = \{ 0, 1, 1\}$.
Although the multiset of reproducible zeros of a polynomial depends on $\mathcal{H}$, for convenience, we have not included $\mathcal{H}$ in our notation. However, all use of this notation will be clear.

We would also like to study the multisets of reproducible points that coincide with reproducible zero multisets of polynomials. The following definition allows us to do that.

\begin{definition}[Reproducible multiset]
A finite multiset $Z$ is a reproducible multiset if it can be written as
\[
Z = \Big\{ \underbrace{0 = \beta_0, \ldots, \beta_0}_{m_0 \textrm{ times}}, \underbrace{\beta_1, \ldots, \beta_1}_{m_1 \textrm{ times}}, \ldots, \underbrace{\beta_s, \ldots, \beta_s}_{m_s \textrm{ times}} \Big\}
\]
where each $\beta_j$ is a distinct reproducible point, $\beta_0 = 0$ appears with multiplicity $m_0$ (possibly zero), and for each $1 \le j \le s$, $\beta_j$ appears with at least multiplicity 1, but no multiplicity higher than its reproducible order, i.e. $1 \le m_j \le \operatorname{ro}(\beta_j)$.
\end{definition}

Note that for any polynomial $p$, $R(p)$ is a reproducible multiset. We can now make a full generalization of the Shapiro--Shields function.

\begin{definition}[Shapiro--Shields function]\label{SS-function}
Let $Z$ be a reproducible multiset and put
\[
Z = \Big\{ \underbrace{0 = \beta_0, \ldots, \beta_0}_{m_0 \textrm{ times}}, \underbrace{\beta_1, \ldots, \beta_1}_{m_1 \textrm{ times}}, \ldots, \underbrace{\beta_s, \ldots, \beta_s}_{m_s \textrm{ times}} \Big\}.
\]
The \textit{Shapiro--Shields function} associated to $Z$ is then defined as
\[
\S_Z = D(k_0^{(m_0)}; k_0^{(m_0 - 1)},\ldots, k_0, k_{\beta_1}^{(m_1 - 1)},\ldots, k_{\beta_1},k_{\beta_s}^{(m_s - 1)},\ldots, k_{\beta_s}).
\]
\end{definition}

It is imperative to note that $\S_Z$ vanishes at each $\beta_j$ with multiplicity $m_j$ when $0 \le j \le s$.
Also, since $\mathcal{P}$ is contained in $\mathcal{H}$, it can be shown that reproducing kernels are linearly independent. Indeed, for any linear combination of kernel functions $f=\sum_{m,n}c_{mn}k_{\beta_n}^{(m)}$ in which not all coefficients $c_{mn}$ are zero, one can always find a polynomial $p$ such that
\[\inner{p}{f} = \langle p, \sum_{m,n}c_{mn}k_{\beta_n}^{(m)} \rangle = \sum_{m,n}\overline{c}_{mn}p^{(m)}(\beta_n) \neq 0.\]
This implies that $f$ is not zero.
In light of Lemma \ref{SS-proj}, $\S_{Z}$ is a nonzero constant multiple of
the projection of $k_0^{(m_0)}$ onto the orthogonal complement of
\[
\operatorname{span}\{k_0^{(m_0 - 1)},\ldots, k_0, k_{\beta_1}^{(m_1 - 1)},\ldots, k_{\beta_1},k_{\beta_s}^{(m_s - 1)},\ldots, k_{\beta_s}\}.
\]

\section{Analogues of Finite Blaschke products}

In this section, we will show that for $p,q \in \mathcal{P}$, $[p] = [q]$ if and only if $R(p) = R(q)$. We will then unify the perspective of Shapiro--Shields functions and  projections of kernels at the origin onto shift-invariant subspaces generated by polynomials, giving \textit{analogues of finite Blaschke products}.

\subsection{Shift Invariant Subspaces}
Note that $k_0^{(d)} - \Pi_{[p]}(k_0^{(d)}) \perp [p]$, so in order to understand $\Pi_{[p]}(k_0^{(d)})$, it is useful to have a characterization of $[p]^\perp$. We do this first when all the zeros of $p$ are reproducible.

\begin{proposition}
\label{P:invariant_f}
Let $f \in \mathcal{P}$ and suppose that
\[
R(f) = Z(f) = \Big\{ \underbrace{\beta_1, \ldots, \beta_1}_{r_1 \textrm{ times}}, \underbrace{\beta_2, \ldots, \beta_2}_{r_2 \textrm{ times}}, \ldots, \underbrace{\beta_n, \ldots, \beta_n}_{r_n \textrm{ times}} \Big\}.
\]
Then
\[
[f] = \Big(\operatorname{span}\big\{k_{\beta_j}^{(\ell)}: 0\leq\ell\leq r_j-1,\ 1\leq j\leq n\big\}\Big)^{\perp}.
\]
\end{proposition}

\begin{proof}
Let $\mathcal{M}$ denote the right hand-side. Then $\mathcal{M}$ consists of all functions $h$ in $\mathcal{H}$ for which $h^{(\ell)}(\beta_j)=0$ for all $0\leq \ell \leq r_j-1$ and $1\leq j\leq n$. Note that for each $p\in\mathcal{P}$, the polynomial $fp$ belongs to $\mathcal{M}$. It follows that $f\mathcal{P}\subseteq\mathcal{M}$ and hence $[f]\subseteq\mathcal{M}$, which implies $\mathcal{M}^{\perp}\subseteq [f]^{\perp}$. On the other hand, since kernel functions are linearly independent (due to the fact that polynomials belong to $\mathcal{H}$), the space $\mathcal{M}^{\perp}$ is of dimension $d=r_1+\cdots+r_n$. To prove the equality, we only need to show that the dimension of $[f]^{\perp}$ is at most $d$.

We have $f\mathcal{P} + \mathcal{P}_{d-1} = \mathcal{P}$, where $\mathcal{P}_{d-1}$ denotes the space of all polynomials of degree at most $d-1$. Taking closure and using the fact that the sum of a closed subspace with a finite dimensional subspace is closed, we have
\[\mathcal{H} = \operatorname{cl}_{\mathcal{H}}(\mathcal{P}) = \operatorname{cl}_{\mathcal{H}}(f\mathcal{P}+\mathcal{P}_{d-1}) = [f] + \mathcal{P}_{d-1}.\]
As a result, the dimension of $[f]^{\perp}$ is at most that of $\mathcal{P}_{d-1}$, which is $d$. Therefore, we have $\mathcal{M}^{\perp}=[f]^{\perp}$, which implies $[f]=\mathcal{M}$ as required.
\end{proof}

This proposition generalizes \cite[Lemma~4.7]{cheng2019inner} where the authors require the zeros of the polynomial to be contained in $\Omega$ and additional properties are imposed on the space.

We will also show that if $f$ has zeros that are not reproducible, this does not change the structure of $[f]$. First though, we need a proposition.

\begin{proposition}
\label{P:unbd_eval_der}
 Let $\beta$ be a complex number and $m$ be a non-negative integer. Then the following statements hold.
 \begin{enumerate}[(a)]
    \item $\beta$ is not a reproducible point if and only if $(z-\beta)$ is cyclic.
   \item $\beta$ is a reproducible point with $\operatorname{ro}(\beta)\leq m$ if and only if \[\operatorname{cl}_{\mathcal{H}}((z-\beta)^{m+2}\mathcal{P})=\operatorname{cl}_{\mathcal{H}}((z-\beta)^{m+1}\mathcal{P}).\]
 \end{enumerate}
\end{proposition}

\begin{proof}
For any integer $k\geq 0$, define $\mathcal{X}_k = (z-\beta)^k\mathcal{P}$. It is clear that $\mathcal{X}_{k+1}\subset\mathcal{X}_{k}$, which shows that the identity $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{k+1})=\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{k})$ holds if and only if $\mathcal{X}_{k+1}$ is dense in $\mathcal{X}_{k}$ with respect to the norm induced from $\mathcal{H}$.

On the other hand, define $\Lambda_{k}:\mathcal{P}\rightarrow\mathbb{C}$ by $\Lambda_{k}(p) = p^{(k)}(\beta)$ for $p\in\mathcal{P}$.
Observe that
\begin{align*}
    \ker(\Lambda_{k}|_{\mathcal{X}_k})
    & = \Big\{(z-\beta)^{k}q(z):\ q\in\mathcal{P}\text{ such that } \Lambda_k\big((z-\beta)^kq(z)\big) = 0\Big\}\\
    & = \Big\{(z-\beta)^{k}q(z):\ q\in\mathcal{P}\text{ such that } q(\beta)=0\Big\}\\
    & = \mathcal{X}_{k+1}.
\end{align*}
It follows from a well-known result in functional analysis (e.g., see Proposition 5.2 and Theorem 5.3 in \cite[Chapter~III]{conway1990course}) that $\Lambda_k|_{\mathcal{X}_k}$ (being a nonzero functional) is unbounded if and only if $\mathcal{X}_{k+1}=\ker(\Lambda|_{\mathcal{X}_k})$ is dense in $\mathcal{X}_k$.

Therefore, we have just showed that for any $k\geq 0$, $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{k+1}) = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{k})$ if and only if the linear function $\Lambda_k|_{\mathcal{X}_k}$ is unbounded.

(a) Since $\mathcal{X}_0=\mathcal{P}$ and $\mathcal{X}_1=(z-\beta)\mathcal{P}$, the function $(z-\beta)$ is cyclic if and only if $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_1) = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_0)$, which, from the argument above, is equivalent to the fact that $\Lambda_0|_{\mathcal{X}_0}$ is unbounded. Since $\Lambda_0(h)=h(\beta)$ for all $h\in\mathcal{P}$, the unboundedness of $\Lambda_0$ means exactly that $\beta$ is not a reproducible point.

(b) Suppose first $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{m+2})=\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{m+1})$. Then the linear functional $\Lambda_{m+1}$ is unbounded on $\mathcal{X}_{m+1}$, hence, unbounded on $\mathcal{P}$. This implies that $\beta$ is not reproducible of order $m+1$. That is, $\operatorname{ro}(\beta)\leq m$.

Let us now prove the converse. Suppose that $\operatorname{ro}(\beta)\leq m$. To simplify the notation, define $n=\operatorname{ro}(\beta)$. Note that for each $0\leq k\leq n$, the linear functional $\Lambda_{k}$ is bounded but $\Lambda_{n+1}$ is unbounded on $\mathcal{P}$. We show that actually $\Lambda_{n+1}|_{\mathcal{X}_{n+1}}$ is unbounded, which then implies that $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+2}) = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+1})$.

Suppose, for the purpose of obtaining a contradiction, that $\Lambda_{n+1}|_{\mathcal{X}_{n+1}}$ is bounded. For any $h\in\mathcal{P}$, we write
\[h = \sum_{0\leq j\leq n}\frac{h^{(j)}(\beta)}{j!}(z-\beta)^{j} + \tilde{h},\] where $\tilde{h}\in\mathcal{X}_{n+1}$. Then
\[\Lambda_{n+1}(h) = \Lambda_{n+1}|_{\mathcal{X}_{n+1}}(\tilde{h})\]
and by triangle inequality,
\begin{align*}
    \|\tilde{h}\|
    & = \Big\|h -\sum_{0\leq j \leq n}\frac{h^{(j)}(\beta)}{j!}(z-\beta)^{j}\Big\| \leq \|h\| + \sum_{0\leq j\leq n}\frac{|h^{(j)}(\beta)|}{j!}\|(z-\beta)^j\|\\
    & \leq \|h\| + \sum_{0\leq j\leq n}\frac{\|\Lambda_j\|}{j!}\|(z-\beta)^j\|\cdot \|h\|.
\end{align*}
Therefore,
\[|\Lambda_{n+1}(h)| = \|\Lambda_{n+1}|_{\mathcal{X}_{n+1}}(\tilde{h})\| \leq \|\Lambda_{n+1}|_{\mathcal{X}_{n+1}}\|\cdot\|\tilde{h}\| \leq C\|h\|,\]
which implies that $\Lambda_{n+1}$ is bounded on $\mathcal{P}$, a contradiction.

We have thus showed that $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+2}) = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+1})$. Now,
\begin{align*}
    \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+3})  & = \operatorname{cl}_{\mathcal{H}}\Big((z-\beta)\cdot\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+2})\Big)\\
    & = \operatorname{cl}_{\mathcal{H}}\big((z-\beta)\cdot\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+1})\big)\\
    & = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{n+2}).
\end{align*}
It then follows inductively that  $\operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{m+2}) = \operatorname{cl}_{\mathcal{H}}(\mathcal{X}_{m+1})$.
\end{proof}

The propositions above allows us to provide a complete description of $[f]$ whenever $f$ is a polynomial.

\begin{theorem}\label{feqfr}
Let $f \in \mathcal{P}$. For each distinct $\beta_j \in R(f)$, let $r_j$ be the multiplicity of $\beta_j$, i.e.
\[
R(f) = \big\{ \underbrace{\beta_1, \ldots, \beta_1}_{r_1 \textrm{ times}}, \underbrace{\beta_2, \ldots, \beta_2}_{r_2 \textrm{ times}}, \ldots, \underbrace{\beta_n, \ldots, \beta_n}_{r_n \textrm{ times}}\big\}.
\]
Then
\[
[f] = \left[ \prod_{\beta \in R(f)} (z - \beta) \right] =  \Big(\operatorname{span}\big\{k_{\beta_j}^{\ell}: 0 \leq \ell \le r_j - 1,\  1 \leq j \leq n \big\}\Big)^{\perp}.
\]
\end{theorem}
\begin{proof}

We first recall the fact that for any multipliers $g,h$ of $\mathcal{H}$, we have
\[
[g\cdot h] = \operatorname{cl}_{\mathcal{H}}(g\cdot[h]).
\]
If $\beta$ is a non-reproducible zero of $f$, then by Lemma \ref{P:unbd_eval_der}, we have $[f]=[f/(z-\beta)]$. So it suffices to consider only the zeros of $f$ with some reproducible order. Put $f(z) = p(z) \prod_{j = 1}^n (z - \beta_j)^{d_j}$, each $\beta_j$ distinct with $d_j \ge r_j$, $1 \le j \le n$, and with the zeros of $p \in \mathcal{P}$ being all of the non-reproducible zeros of $f$ (i.e. $p$ is cyclic). Then $[f] = [f/p]$. So without loss of generality, we may assume that $p(z)$ is identically one. Let $h(z) = (z - \beta_1)^{d_1}$. Then by Lemma \ref{P:unbd_eval_der}, we have $[h] = [(z - \beta_1)^{d_1}] = [(z - \beta_1)^{r_1}]$. Letting $g = f/h$, we have
\begin{align*}
    [f] & = [g\cdot h] = \operatorname{cl}_{\mathcal{H}}(g\cdot[h])\\
& = \operatorname{cl}_{\mathcal{H}}\big(g\cdot[(z-\beta_1)^{d_1}]\big) = [g\cdot (z-\beta_1)^{r_1}]\\
& = \left[\frac{f}{(z-\beta_1)^{d_1 - r_1}}\right].
\end{align*}
Repeating this argument for each $\beta_j$, $2 \le j \le n$, we have
\begin{align*}
    [f] & = \left[\frac{f}{\prod_{j=1}^{n}(z - \beta_j)^{d_j - r_j}}\right] = \Big[\prod_{j=1}^{n}(z-\beta_j)^{r_j}\Big]  = \left[ \prod_{\beta \in R(f)} (z - \beta) \right].
\end{align*}
Applying Proposition \ref{P:invariant_f} then gives the result.
\end{proof}

This theorem, to which was previously alluded, is a generalization proved for the $H^2_w$ spaces by the first author in \cite[Theorem~5.4]{felder2020general}. Let us illustrate this theorem by applying it to an example in various spaces.

\begin{example}
Let $f(z)=z^2(z-\frac{i}{2})(z^2-1)^2$. Then the zero multiset of $f$ is
\[Z(f) = \Big\{0,0,\frac{i}{2}, -1, -1, 1, 1\Big\}.\]
Interestingly, by Theorem \ref{feqfr}, the shift-invariant subspace $[f]$ depends greatly on the underlying Hilbert space.

If $\mathcal{H}=D_{\alpha}$ for $\alpha\leq 1$ (which includes the Hardy, Bergman and Dirichlet spaces), then $R(f)=\{0,0,i/2\}$ and
\[[f] = \Big(\operatorname{span}\big\{k_{0}, k_{0}^{(1)}, k_{i/2}\big\}\Big)^{\perp}.\]
If $\mathcal{H}=D_{\alpha}$ for $1<\alpha\leq 3$, then by Example \ref{D-kernels}, we have $R(f)=\{0,0,i/2,-1,1\}$, which then implies that
\[[f] =\Big(\operatorname{span}\big\{k_{0}, k_{0}^{(1)}, k_{i/2}, k_{-1}, k_{1}\big\}\Big)^{\perp}.\]
If $\mathcal{H}=D_{\alpha}$ for $\alpha>3$, then by Example \ref{D-kernels} again, $R(f)=Z(f)$ and so
\[[f] = \Big(\operatorname{span}\big\{k_{0}, k_{0}^{(1)}, k_{i/2}, k_{-1}, k_{-1}^{(1)}, k_{1}, k_{1}^{(1)}\big\}\Big)^{\perp}.\]
On the other hand, if $\mathcal{H} = \mathcal{D}_{\delta_{1}}$, the local Dirichlet space at $1$, then by Example \ref{Ex:localD}, we have $R(f)=\{0,0,i/2,1\}$ and hence,
\[[f] = \Big(\operatorname{span}\big\{k_{0}, k_{0}^{(1)}, k_{i/2}, k_{1}\big\}\Big)^{\perp}.\]
\end{example}

Theorem \ref{feqfr} also has an immediate and useful corollary.
\begin{corollary}\label{equal-perp}
Let $p, q \in \mathcal{P}$. Then $[p] = [q]$ if and only if $R(p) = R(q)$.
\end{corollary}
\begin{proof}
The backward implication is given directly by Theorem \ref{feqfr}.
So let $[p] = [q]$ and suppose for contradiction that $R(p) \neq R(q)$. WLOG, there exists $\beta \in R(p)$ with $\beta \notin R(q)$ or $\beta$ having greater multiplicity in $R(p)$ than in $R(q)$. In either case, Theorem \ref{feqfr} implies there is some $n \ge 0$ so that $k_\beta^{(n)} \perp [p] = [q]$. But then $\langle p, k_\beta^{(n)} \rangle = \langle q, k_\beta^{(n)} \rangle = 0$, which is a contradiction, since $\beta \notin R(q)$ or $\beta$ has multiplicity strictly less than $n+1$ in $R(q)$.
\end{proof}

\subsection{Inner Functions and Linear Combinations of Kernels}
We now show that each inner function that arises as a certain linear combination of reproducing kernels can be identified with a shift invariant subspace and a Shapiro--Shields function. The following theorem generalizes a result of the second author in \cite[Theorem~3.7]{LeAMP2020}, proved initially in the $H^2_w$ spaces. The significance of our contribution is that we do not require monomials be orthogonal and require very little geometric assumption on the underlying set for which $\mathcal{H}$ is an RKHS, providing a very general setting for which these results hold. When $\mathcal{H} = H^2$, this result describes classical Blaschke products and in general gives our \textit{analogues of finite Blaschke products}.

\begin{theorem}\label{modified_le}
Suppose that
\[
B = \sum_{j=0}^{s}\sum_{\ell=0}^{m_{j}}c_{j,\ell}\,k_{\lambda_j}^{(\ell)}
\]
is an $\mathcal{H}$-inner function, with $c_{j,m_{j}}\neq 0$.
Then $B$ is a constant multiple of $\Pi_{[f]}(k_0^{(d)})$ for some $d$ and some polynomial $f$. Further, $B$ must also be a Shapiro--Shields function.
\end{theorem}

The function $B$ here is what we call an \textit{analogue of a finite Blaschke product}.

\begin{proof}
Take $B$ as above and without loss of generality, assume that $\|B\|=1$. Let $g \in \mathcal{P}$. Then
\begin{align*}
    g(0) & = \inner{gB}{B} = \sum_{j=0}^{s}\sum_{\ell=0}^{m_{j}}\bar{c}_{j,\ell}\inner{gB}{k_{\lambda_j}^{(\ell)}}\\
    & = \sum_{j=0}^{s}\sum_{\ell=0}^{m_{j}}\bar{c}_{j,\ell}\sum_{m=0}^{\ell}\binom{\ell}{m}g^{(m)}(\lambda_j)\cdot B^{(\ell-m)}(\lambda_j)\\
    & = \sum_{j=0}^{s}\sum_{m=0}^{m_j}\Big\{\sum_{\ell=m}^{m_j}\bar{c}_{j,\ell}\binom{\ell}{m}B^{\ell-m}(\lambda_j)\Big\}g^{(m)}(\lambda_j).
\end{align*}
It follows that $0\in\{\lambda_0,\ldots,\lambda_s\}$ and for all $j$, we have
\begin{align}
\label{Eqn:inner_condition_B}
    \sum_{m=0}^{m_j}\Big\{\sum_{\ell=m}^{m_j}\bar{c}_{j,\ell}\binom{\ell}{m}B^{(\ell-m)}(\lambda_j)\Big\}g^{(m)}(\lambda_j) =
    \begin{cases}
    0, & \text{ if } \lambda_j \neq 0\\
    g(0), & \text{ if } \lambda_j = 0
    \end{cases}
\end{align}
for all polynomials $g$. Without loss of generality, we shall always assume that $\lambda_0=0$. Then for $j=0$ and $0\leq m\leq m_0$, identity \eqref{Eqn:inner_condition_B} gives
\[\sum_{\ell=m}^{m_0}\bar{c}_{0,\ell}\binom{\ell}{m}B^{(\ell-m)}(0) =
\begin{cases}
1, & \text{ if } m=0,\\
0, & \text{ if } m\geq 1.
\end{cases}
\]
Since $c_{0,m_0}\neq 0$, we conclude that either $m_0=0$, or $m_0\geq 1$ and $B^{(\ell)}(0)=0$ for all $0\leq\ell\leq m_0-1$. That is, $B\perp\{k_{0}^{(\ell)}: 0\leq\ell\leq m_0-1\}$.

On the other hand, for $j\neq 0$ and $0\leq m\leq m_j$, by identity \eqref{Eqn:inner_condition_B},
\[\sum_{\ell=m}^{m_j}\bar{c}_{j,\ell}\binom{\ell}{m}B^{(\ell-m)}(\lambda_j) = 0.\]
Since $c_{j,m_j}\neq 0$, it follows that $B^{(\ell)}(\lambda_j)=0$ for all $0\leq\ell\leq m_j$. That is, $B\perp\{k_{\lambda_j}^{(\ell)}: 0\leq\ell\leq m_j\}$. Let $\mathcal{M}$ be the subspace spanned by the functions
\[\big\{k_{0}^{(\ell)}: 0\leq\ell\leq m_0-1\big\}\cup\big\{k_{\lambda_j}^{(\ell)}: 1\leq j\leq s, 0\leq\ell\leq m_j\big\},\] where the first set is considered to be empty if $m_0=0$. Then we have $B = c_{0,m_0}\Pi_{\mathcal{M}^{\perp}}(k_{0}^{(m_0)})$. By Proposition \ref{P:invariant_f}, $\mathcal{M}^{\perp}$ can be recognized as $[f]$, where \[f(z)=z^{m_0}\prod_{j=1}^{s}(z-\lambda_j)^{m_j+1}.\]
As a result, we have shown that $B$ is a projection of $k_0^{(d)}$ onto $[f]$ for some $d$ and some polynomial $f$. Further, by Lemma \ref{SS-proj}, we know that $B$ must also equal the Shapiro-Shield function $\S_{Z(f)}$.
\end{proof}

\begin{remark}
Recall that finite Blaschke products are the only rational inner functions on the Hardy space (e.g. see \cite[Section~3.5]{GarciaSC2018}). In fact, since any rational function with poles outside of the closed unit disk can be written as a linear combination of $H^2$-kernel functions, this also follows from our Theorem \ref{feqfr}.

We can also apply our result to investigate rational inner functions on the Bergman space $A^2$. Suppose $\lambda_1,\ldots,\lambda_s$ are distinct nonzero points on the open unit disk and a function of the form
\[B(z) = c_0 + \sum_{j=1}^{s}c_{j}k_{\lambda_j} = c_0 + \sum_{j=1}^{s}\frac{c_j}{(1-\overline{\lambda}_jz)^2}\]
is $A^2$-inner. Theorem \ref{feqfr} shows that $B$ must vanish at all these points. Therefore, there exists a polynomial $q$ of degree at most $s$ such that
\[B(z) = \frac{q(z)\prod_{j=1}^{s}(z-\lambda_j)}{\prod_{j=1}^{s}(1-\overline{\lambda}_jz)^2}.\qquad \]
On the other hand, it follows from the formula for $B$ as a linear combination of Bergman kernels, all residues of $B$ at $1/\overline{\lambda}_1,\ldots, 1/\overline{\lambda}_s$ are zero. As a consequence, $q$ is uniquely determined (up to a constant) from this condition. This approach provides a different way to construct analogues of finite Blaschke products in $A^2$, which does not involve determinants. The above argument can be adapted for the case of repeated points.
\end{remark}

\begin{corollary}
\label{C:rational_inner_functions}
Assume that any function that is analytic and non-vanishing on a neighborhood of $\overline{\Omega}$ is cyclic in $\mathcal{H}$ (this is true for all weighted Hardy spaces $\mathcal{H}_{w}$, including the Hardy, Bergman and Dirichlet-type spaces). Let $f$ be an $\mathcal{H}$-inner function that is analytic on a neighborhood of $\overline{\Omega}$.  Then $f$ must be a linear combination of kernel functions and hence, $f$ is determined by Theorem \ref{modified_le}. In particular, all rational $A^2$-inner functions are Shaprio--Shields functions.
\end{corollary}

\begin{proof}
Write $f(z) = p(z)g(z)$, where $p$ is a polynomial and $g$ is analytic and non-vanishing on a neighborhood of $\overline{\Omega}$. By the definition of inner functions, we see that $f$ belongs to the orthogonal complement of $[Sf] = [qg]$, where $q(z)=zp(z)$. Since $g$ is cyclic, it follows that $[Sf]=[q]$. On the other hand, Theorem \ref{feqfr} shows that $[q]^{\perp}$ is spanned by a finite collection of kernel functions. As a result, $f$ is of the form described in Theorem \ref{modified_le}.
\end{proof}

We can also provide a description of the various ways in which we view analogues of finite Blaschke products.
\begin{prop}\label{projections}
Let $f \in \mathcal{P}$ and write
\[
R(f) = \Big\{\underbrace{0 = \beta_0, \ldots, \beta_0}_{r_0 \textrm{ times}}, \underbrace{\beta_1, \ldots, \beta_1}_{r_1 \textrm{ times}}, \ldots, \underbrace{\beta_n, \ldots, \beta_n}_{r_n \textrm{ times}} \Big\}
\]
and let $\tilde{f}(z) = \prod_{\beta \in R(f)} (z - \beta)$.
Then, up to a constant, the following are equal:
\begin{enumerate}
\item $\Pi_{[f]}(k_0^{(r_0)})$
\item $\Pi_{[\tilde{f}]}(k_0^{(r_0)})$
\item The Shapiro--Shields function $\S_{R(f)}$, associated to the reproducible zeros of $f$.
\item $\varphi := k_0^{(r_0)} - \sum_{j = 0}^{n} \sum_{\ell=0}^{r_j-1} c_{\ell,j}k^{(\ell)}_{\beta_j}$, where the constants $c_{\ell,j}$ are given by $\langle \varphi, k_{\beta_j}^{(\ell)} \rangle = 0$ for $0 \le \ell \le r_j-1$ and $0 \le j \le n$.
\end{enumerate}

\end{prop}
\begin{proof}
The fact that $\Pi_{[f]}(k_0^{(r_0)}) = \Pi_{[\tilde{f}]}(k_0^{(r_0)})$ is given by Theorem \ref{feqfr}.
Now, let $\mathcal{M} = \textrm{span}\left\{k_{\beta_j}^{(\ell)}: 0\leq\ell\leq r_j-1,\ 0\leq j\leq n \right\}^{\perp}$. Again, by Theorem \ref{feqfr}, $\mathcal{M} = [f] = [\tilde{f}]$. So by Lemma \ref{SS-proj}, we have $\S_{R(f)}$ and $\Pi_{[\tilde{f}]}(k_0^{(r_0)})$ are equal up to a constant multiple.
Following the discussion after Definition \ref{SS-function}, we know that $\langle \Pi_{[\tilde{f}]}(k_0^{(r_0)}), k_{\beta_j}^{(\ell)} \rangle=0$ for $0 \le \ell \le r_j-1$ and $0 \le j \le n$, which is also required of $\varphi$. Since $\dim(\mathcal{M}^\perp) = r_0 + \cdots + r_n$, this uniquely determines both $\Pi_{[\tilde{f}]}(k_0^{(r_0)})$ and $\varphi$, which then must be equal up to a constant multiple.
\end{proof}

We have included the equivalence $(iv)$ above in order to provide a more computationally explicit description. Conditions $(i), (ii)$, and $(iv)$ generalize results of the first author in \cite{felder2020general}, initially proven for $H^2_w$ spaces.

\section{Extraneous Zeros}

So far, we have seen various descriptions of analogues of finite Blaschke products, namely those in Proposition \ref{projections}. We also saw that for $p, q \in \mathcal{P}$, $R(p) = R(q)$ if and only if $[p] = [q]$.
However, in certain settings, we will see that we have a surprising strengthening of this fact. Namely, Corollary \ref{single_vec_proj} will show, under certain assumptions, that $[p] = [q]$ if and only if $\Pi_{[p]}(k_0^{(d)}) = \Pi_{[q]}(k_0^{(d)})$, for appropriate $d$.

In a similar fashion, one may ask, for different reproducible multisets $A$ and $B$, is it possible that $\S_A$ is a constant multiple of $\S_B$?
We start by pointing out that for a reproducible multiset $Z$, the function $\S_Z$ can vanish off of $Z$ (e.g. see discussion in  \cite[Section~5]{cheng2019inner}). In this case, we naturally say that $\S_Z$ has an extraneous zero. On the other hand, the following lemma asserts that any such extraneous zero must be different from the origin.

\begin{lemma}
\label{L:SS_extraneous_zero}
Let $Z$ be a reproducible multiset with $0$ appearing $m_0$ times. Then $\S_Z$ vanishes at the origin with order precisely $m_0$. Namely, any extraneous zero of a Shapiro--Shields function cannot be at the origin.
\end{lemma}

\begin{proof}
Let $p(z) = \prod_{\beta \in Z}(z - \beta)$. Then the inner product $\inner{\S_{Z}}{k_0^{(m_0)}}$ is a nonzero constant multiple of
\[
    \binner{\Pi_{[p]}(k_0^{(m_0)})}{k_0^{(m_0)}} = \big\|\Pi_{[p]}(k_0^{(m_0)})\big\|^2,
\]
which is not zero since, by hypothesis, $k_0^{(m_0)}$ is not orthogonal to $p$. As a result, any extraneous zero of $\S_Z$ must be different from the origin.
\end{proof}

We will also need the following simple lemma about projection on orthogonal complements of sets of vectors.
\begin{lemma}
\label{L:orthogonal_proj}
Let $M$ and $N$ be two sets of vectors in a Hilbert space. Then for any $v$, we have $\Pi_{(M\cup N)^{\perp}}(v)=\Pi_{M^{\perp}}(v)$ if and only if $\Pi_{M^{\perp}}(v)\perp N$.
\end{lemma}

\iffalse
\begin{proof}
The necessary direction is clear. For the sufficient direction, assume that $\Pi_{M^{\perp}}(v)\perp N$. Then we have $\Pi_{M^{\perp}}(v)$ belongs to $(M\cup N)^{\perp}$. In addition, note that
\[v-\Pi_{M^{\perp}}(v)\in (M^{\perp})^{\perp} \subseteq ((M\cup N)^{\perp})^{\perp}.\]
It then follows that $\Pi_{M^{\perp}}(v)$ is the orthogonal projection of $v$ onto $(M\cup N)^{\perp}$, as required.
\end{proof}
\fi

In the following proposition, we show that shift-invariant subspaces are characterized uniquely by the projection of a kernel function at the origin, as long as Shapiro--Shields functions do not possess extraneous reproducible zeroes.

\begin{proposition}\label{extraneous}
The following two statements are equivalent.
\begin{itemize}
    \item[(a)] All Shapiro--Shields functions do not possess extraneous zeros that are reproducible points.
    \item[(b)] For any two polynomials $p$ and $q$, if $\Pi_{[p]}(k_0^{(d_p)})=\Pi_{[q]}(k_0^{(d_q)})$, then $[p]=[q]$. Here $d_p$ (respectively, $d_q$) denotes the multiplicity of the zero of $p$ (respectively, $q$) at the origin.
\end{itemize}
\end{proposition}

\begin{proof}
Assume that for some multiset $Z$ of reproducible points with $\beta_0 = 0$ appearing $m_0 \ge 0$ times and non-zero values $\beta_1, \ldots, \beta_s$ each appearing with multiplicity $m_j$, $1 \le m_j \le \textrm{ord}(\beta_j)$, the corresponding Shapiro--Shields function
\[\S_{Z} = D\big(k_0^{(m_0)}; k_0^{(m_0 - 1)},\ldots, k_0, k_{\beta_1}^{(m_1-1)},\ldots, k_{\beta_1},k_{\beta_s}^{(m_s-1)},\ldots, k_{\beta_s}\big)\]
has an extraneous zero. This means that either $\S_{Z}$ has an additional zero at a reproducible point $\beta\notin Z$ or it has a zero of multiplicity at least $m_j+1$ at some $\beta=\beta_j$. By Lemma \ref{L:SS_extraneous_zero}, $\beta\neq 0$. Note that $\S_{Z}$ is a non-zero constant multiple of $f=\Pi_{[p]}(k_0^{(m_0)})$, where
\[p(z) = \prod_{j=0}^{s}(z-\beta_j)^{m_j}.\]
It follows that $f(\beta)=0$ (so $f\perp k_{\beta}$) in the case $\beta\notin Z$, or $f^{(m_j)}(\beta_j)=0$ (so $f\perp k_{\beta_j}^{(m_j)}$) in the case $\beta=\beta_j$ for some $j$. Put $q(z) = (z-\beta)p(z)$. Then $f$ coincides with $\Pi_{[q]}(k_0^{(m_0)})$, by Proposition \ref{P:invariant_f} and Lemma \ref{L:orthogonal_proj}. Since $\beta\neq 0$, we see that $d_q = d_p = m_0$ and hence $\Pi_{[p]}(k_0^{d_p}) = \Pi_{[q]}(k_0^{d_q})$ even though $[p]\neq [q]$.

Conversely, assume that (a) holds. Let $p, q$ be two polynomials. Using Theorem \ref{feqfr}, we write
\[[p] = \Big\{k_{\beta_j}^{(\ell)}:  0\leq\ell < m_j,\ 0\leq j\leq s\Big\}^{\perp},\]
and
\[[q] = \Big\{k_{\alpha_j}^{(\ell)}:  0\leq\ell < n_j,\ 0\leq j\leq t\Big\}^{\perp}.\]
We know that $\Pi_{[p]}(k_0^{(d_p)})$ is a non-zero constant multiple of the Shapiro-Shield function corresponding to a multiset $Z_1$ consisting of $\beta_j$ with multiplicity $m_j$ for $0\leq j\leq s$. In the expansion of $\Pi_{[p]}(k_0^{(d_q)})$, for each $1\leq j\leq s$, the coefficient of $k_{\beta_j}^{(m_j-1)}$ is plus or minus of $\inner{\S_{\tilde{Z}_1}}{k_{\beta_j}^{(m_j-1)}}$, where $\tilde{Z}_1$ is obtained from $Z_1$ by removing one appearance of $\beta_j$. By the hypothesis, this coefficient is nonzero. By the same argument, the coefficient of any $k_{\alpha_j}^{(n_j-1)}$ in the expansion of $\Pi_{[q]}(k_0^{(d_q)})$ is nonzero. It follows that if $\Pi_{[p]}(k_0^{(d_p)})=\Pi_{[q]}(k_0^{(d_q)})$, then due to the linear independence of the kernel functions, each $\beta_j$ equals some $\alpha_{\ell}$ and $m_j = n_{\ell}$, and vice versa. Therefore, $[p]=[q]$.
\end{proof}

In general, the projection of a single vector onto a subspace does not characterize the subspace itself. However, in the setting of Theorem \ref{extraneous}, we have a surprising immediately corollary.

\begin{corollary}\label{single_vec_proj}
Let $p,q \in \mathcal{P}$ with $\operatorname{ord}_0(p) = d_p$ and $\operatorname{ord}_0(q) = d_q$.  Suppose that $\Pi_{[p]}(k_0^{(d_p)})$ and $\Pi_{[q]}(k_0^{(d_q)})$ vanish only on $R(p)$ and $R(q)$, respectively. Then, $[p] = [q]$ if and only if $\Pi_{[p]}(k_0^{(d_p)}) = \Pi_{[q]}(k_0^{(d_q)})$.
\end{corollary}

When Shapiro--Shields functions do possess extraneous zeros, we show one way in which two different reproducible multisets can generate the same Shapiro--Shields function.

\begin{theorem}
Let $A$ be a reproducible multiset. Then for any reproducible multiset $B$ such that $A\subset B\subset Z(\S_A)$, the function $\S_B$ is a constant multiple of $\S_{A}$.
\end{theorem}

\begin{proof}
It suffices to prove the case $B$ having exactly one element more than $A$. Write
\[A = \Big\{\underbrace{0 = \beta_0, \ldots, \beta_0}_{r_0 \textrm{ times}}, \underbrace{\beta_1, \ldots, \beta_1}_{r_1 \textrm{ times}}, \ldots, \underbrace{\beta_n, \ldots, \beta_n}_{r_n \textrm{ times}} \Big\},
\]
where $r_0\geq 0$ and $r_1,\ldots, r_n\geq 1$. Then $\S_A$ is a nonzero constant multiple of the projection $\Pi_{M^{\perp}}(k_0^{(r_0)})$, where
\[M = \{k_{\beta_j}^{(\ell)}: 0\leq\ell< r_j,\ 0\leq j\leq n\}.\]
Let $\beta$ be the extra element that $B$ has more than $A$. Note that $\beta\neq 0$ by Lemma \ref{L:SS_extraneous_zero}. If $\beta=\beta_{s}$ for some $1\leq s\leq n$, set $u = k_{\beta_s}^{(r_s)}$. If $\beta\neq\beta_j$ for all $j$, let $u=k_{\beta}$. Since $B\subset Z(\S_A)$, we have that $\Pi_{M^{\perp}}(k_0^{(r_0)})\perp u$. As a consequence,
\[\Pi_{M^{\perp}}(k_0^{(r_0)}) = \Pi_{(M\cup\{u\})^{\perp}}(k_0^{(r_0)}).\]
On the other hand, $\S_B$ is a nonzero constant multiple of  $\Pi_{(M\cup\{u\})^{\perp}}(k_0^{(r_0)})$. The conclusion of the theorem then follows.
\end{proof}

We conclude with a few remarks. For $f$ a polynomial, we have described $[f]$ in general reproducing kernel Hilbert space of analytic functions. As we have seen, $[f]$ depends closely on reproducible points of the underlying space. It would be interesting to go beyond the polynomial case but it seems to be a more difficult problem.

Blaschke products (and more generally, inner functions) on the Hardy space possess a remarkable multiplicative property: if $B_1$ and $B_2$ are finite Blaschke products, then $B_1B_2$ is also a finite Blaschke product. On the other hand, this property does not hold on other spaces, such as the Bergman and Dirichlet spaces, which can be seen from  direct calculation. While the multiplicative property of classical Blaschke products trivially follows from the definition, it is quite curious why such a property even holds true, in light of Theorem  \ref{modified_le}.
In terms of extremal problems, it is somehow remarkable that the solution to the linear $n$-point extremal problem is given by the \textit{product} of the $n$ individual 1-point problems.

Lastly, we would like to mention that analyticity is not actually needed to establish the results throughout this paper. Instead, one could make the weaker assumption that functions in $\mathcal{H}$ are simply smooth at the origin. We have chosen to require analyticity because it is unclear if there are any (interesting) RKHSs for which the polynomials are dense and the shift is bounded, yet not comprised of analytic functions.

\subsection*{Acknowledgements} The first author would like to thank John McCarthy and Dima Khavinson for helpful discussion.


\section{Introduction}\label{section: introduction}
Estimating the covariance function has an essential role in various fields involving spatial data analysis, such as in climategoraphy, ecology and neurosceince \citep{hohn1998geostatistics,zhong2008classifying,mihoub2016modeling}. A common example for estimating the covariance function is in fitting Gaussian process regression (GPR) for predicting spatial data \citep{rasmussen2003gaussian}. However, estimating the covariance function is also required in many other settings, including unsupervised learning tools, such as in dimensional reduction using Gaussian process latent variable models \citep{li2016review}.

Commonly, the covariance function of a sample $\boldsymbol{z}(S)=[z(\boldsymbol{s}_1),...,z(\boldsymbol{s}_n)]$ --- where $z(\boldsymbol{s}_i)\coloneqq z_i\in\mathbb{R}$ is sampled from the random field $\mathcal{Z}(\mathcal{S})$ at the coordinate values $\boldsymbol{s}_i\in\mathcal{S}\subset\mathbb{R}^{q}$ ---
is estimated using a \textit{kernel function},  $\mathcal{K}_{\boldsymbol{\theta}}(\boldsymbol{s}_i-\boldsymbol{s}_j):\; \mathbb{R}^q\to \mathbb{R},$ where $\boldsymbol{\theta}\in\Theta$ is the kernel's parameters. A popular kernel function is the exponential kernel:
\begin{align}\label{Kernel: exponential, isotropy}
\sigma_{s}^2\times\exp{\big(-\frac{\|\boldsymbol{h}_{i,j}\|}{\lambda}\big)+\sigma^2_{\epsilon}},
\end{align}
where $\boldsymbol{h}_{i,j}=\boldsymbol{s}_i-\boldsymbol{s}_j,$
$\sigma_{s}^2\in \mathbb{R}^{+}$ is the signal parameter, $\lambda\in \mathbb{R}^{+}$ is the length-scale parameter and $\sigma^2_{\epsilon}\in \mathbb{R}^{+}$ is the variance of the independent error term. In this case, $\boldsymbol{\theta}=[\text{exponential kernel},\sigma_{s}^2,\lambda,\sigma^2_{\epsilon}].$\footnote{Commonly, the variance of the error term, $\sigma^2_{\epsilon},$ which does not depend on the distance is not included in the kernel and therefore the covariance function is the kernel plus the error term, however in order to simplify writing the error term is included in the kernel function.}

The exponential kernel, as well as many other kernels, assumes weak stationarity, i.e.,
\begin{itemize}
    \item $\mathbb{E} z(\boldsymbol{s})=\mathbb{E} z(\boldsymbol{s}+\boldsymbol{s}_{\tau})=\mu,$ where $\mu\in \mathbb{R}$ and $\boldsymbol{s}_{\tau}\in \mathbb{R}^{q},$ such that $\boldsymbol{s}+\boldsymbol{s}_{\tau}\in\mathcal{S}.$
    \item $\mathrm{Cov}\big(z(\boldsymbol{s}_i),z(\boldsymbol{s}_j)\big)=\mathrm{Cov}\big(z(\boldsymbol{s}_i+\boldsymbol{s}_{\tau}),z(\boldsymbol{s}_j+\boldsymbol{s}_{\tau})\big).$
\end{itemize}
Under the stationary assumption the covariance function can be denoted by $C(\boldsymbol{h}),$ which emphasizes that the directed distance between the locations is the sufficient argument for the covariance function (rather than the coordinate values themselves). In many use cases, the raw data should be preprocessed in order to be stationary.

A stronger assumption than stationarity is rotational invariance, $C(\boldsymbol{h})=C(\|\boldsymbol{h}\|),$ which is called \textit{isotropy}. The isotropy assumption is taken frequently, also in cases it does not hold \citep{rajala2018review}, i.e., in \textit{anisotropic} settings. There are several anisotropy types, the most common is range anisotropy (also called geometric anisotropy), where $C(\boldsymbol{h})$ decreases differently in different directions, but the sill, $\text{lim}_{\|\boldsymbol{h}\|\to\infty}C(\boldsymbol{h}),$ and the nugget effect, $\text{lim}_{\|\boldsymbol{h}\|\to 0}C(\boldsymbol{h}),$ do not depend on the direction. Figure \ref{Figure: Demonstrating Settings} demonstrates anisotropy and isotropy settings.

\begin{figure}[h!]
    \centering
        \includegraphics[width=1\linewidth]{Figures/scatterplot2021-06-05n5000.pdf}
         \caption{Simulated data of isotropy and anisotropy settings. Anisotropy setting: The left figure presents covariance functions in two directional axes that decay differently as a function of the distance (with $log_{10}$ scale). The green line is the covariance function for the longitudinal axis, and the blue line is for the latitudinal axis. The middle figure presents a simulated sample of an anisotropy setting. As once can see the variance in the latitudinal axis is larger than in the longitudinal axis. This is due to the higher correlation in the longitudinal axis than in the latitudinal axis. Isotropy setting: The right figure presents a simulated sample of an anisotropy setting.}\label{Figure: Demonstrating Settings}
    \end{figure}

A specific type of range anisotropy is elliptic anisotropy, when a linear transformation of $\boldsymbol{h}$ induces isotropy. For example, the standard exponential kernel function, which assumes isotropy (see expression (\ref{Kernel: exponential, isotropy})), can be generalized using rotation and scaling matrices in order to capture elliptic anisotropy:
\begin{align}\label{Kernel: exponential, anisotropy}
   \sigma_{s}^2\times\exp{\big(-\|A\boldsymbol{h}\|}\big)+\sigma^2_{\epsilon}=\sigma_{s}^2\times\exp{\big(-\sqrt{\boldsymbol{h}^{t}A^{t}A\boldsymbol{h}}\big)}+\sigma^2_{\epsilon},
\end{align}
where for $q=2:$
\[
A\coloneqq A(\lambda_1,\lambda_2,\eta)=
\begin{bmatrix}
\frac{1}{\lambda_1} & 0 \\
0 & \frac{1}{\lambda_2}
\end{bmatrix}
\begin{bmatrix}
\cos(\eta) & -\sin(\eta) \\
\sin(\eta) & \cos(\eta)
\end{bmatrix}
,
\]
$\eta\in[0,\pi]$  and $\eta+\pi/2$ are the anisotropy direction axes and $\lambda_i\in\mathbb{R}^{+}$ are the anisotropic scales. When $\lambda_1=\lambda_2$ expression (\ref{Kernel: exponential, anisotropy}) is reduced to expression (\ref{Kernel: exponential, isotropy}).

Several hypothesis testing approaches for detecting anisotropy were proposed in recent years (for a literature review see Section \ref{section: Literature Review}), most of them focus on testing elliptic anisotropy. In many cases the tests are implemented on the variogram function which is a transformation of the covariance function:
\[
\gamma(\boldsymbol{h}_{i,j})=\frac{1}{2}\mathrm{Var}\big(z(\boldsymbol{s}_i)-z(\boldsymbol{s}_j)\big).
\]
It is easy to see that $\gamma(\boldsymbol{h}_{i,j})=C(0)-C(\boldsymbol{h}_{i,j}).$

Frequently, the variogram is estimated by the empirical variogram. The empirical variogram in the range $\boldsymbol{h}^*\pm\delta,$ where $\delta\in\mathbb{R}^d$ is:
\[
\widehat{\gamma}(\boldsymbol{h})=\frac{1}{2\times|\mathcal{H}^*|}\sum_{(i,j)\in \mathcal{H}^*}\big(z(\boldsymbol{s}_i)-z(\boldsymbol{s}_j)\big)^2,\;\boldsymbol{h}\in\boldsymbol{h}^*\pm\delta\]
where $\mathcal{H}^*=\{(i,j)|\boldsymbol{s}_i-\boldsymbol{s}_j=\boldsymbol{h}^*\pm\delta\}$ and $|\mathcal{H}^*|$ is the size of the set $\mathcal{H}^*.$
In this paper we mainly refer to the covariance function.

\section{Our Approach}
 In this section we present our algorithms for testing anisotropy using re-sampling mechanisms. The first approach is based on parametric bootstrap hypothesis testing \citep{mackinnon2009bootstrap} and therefore requires parametric assumptions, the second approach is a non-parametric approach and therefore requires fewer assumptions.
\subsection{Parametric Bootstrap Based Test}\label{section: Parametric Bootstrap}
As in a standard hypothesis testing framework, two hypotheses are compared --- the null hypothesis ($H_0$) assumes that $\boldsymbol{z}(S)$ was sampled from a distribution with an isotropic covariance function, while the alternative hypothesis ($H_1$) assumes that $\boldsymbol{z}(S)$ was sampled from a distribution with an anisotropic covariance function. $H_1$ can specify the suspected anisotropic directional  axes, $\{\eta_i\}_{i\in[1,...,R]},$ or avoid specifying the anisotropic directional axes, but only posit that $R$ anisotropic directional axes exist. Another setting is when $H_1$ specifies ranges of $\{\eta_i\}_{i\in[1,...,R]}.$ For example, using the exponential kernel function family, elliptical anisotropy can be tested when $H_0$ assumes that the covariance function is the standard exponential kernel (expression (\ref{Kernel: exponential, isotropy})), while $H_1$ assumes that the covariance function is the elliptical exponential kernel function, presented in expression (\ref{Kernel: exponential, anisotropy}).
Other properties, such as stationarity and normality of $\boldsymbol{z}(S),$ are assumed equally by both hypotheses. The normality assumption will be relaxed in Section \ref{section: Non-Parametric Rotational Sampling Test}.

Algorithm \ref{Algorithm: parametric} presents our parametric bootstrap hypothesis testing for anisotropy detection.
\begin{algorithm}[ht!]
\caption{Parametric Bootstrap Hypothesis Testing for Anisotropy Detection}
\label{Algorithm: parametric}
\begin{algorithmic}[1]
\STATE {\bfseries Input:} $\{S,\boldsymbol{z}(S)\}.$
\STATE {\bfseries Output:} P-value.
\STATE  Estimate $\mathcal{K}_{\boldsymbol{\theta}|H_0}(\cdot),\;\mathcal{K}_{\boldsymbol{\theta}|H_1}(\cdot)$ and denote their covariance matrices scores at $S$ as $\widehat{\Sigma}_{\boldsymbol{\theta}|H_0},\;\widehat{\Sigma}_{\boldsymbol{\theta}|H_1},$ respectively.
\STATE Calculate the following anisotropic discrepancy measure:
    \[\phi=\ell\big(\boldsymbol{z};\widehat{\Sigma}_{\boldsymbol{\theta}|H_1}\big)-\ell\big(\boldsymbol{z};\widehat{\Sigma}_{\boldsymbol{\theta}|H_0}\big),\]
    where $\ell\big(\boldsymbol{z};\widehat{\Sigma}_{\boldsymbol{\theta}|H_0}\big)$ and $\ell\big(\boldsymbol{z};\widehat{\Sigma}_{\boldsymbol{\theta}|H_1}\big)$ are the log-likelihood of $\boldsymbol{z}$ under the two hypotheses.
\FOR{$b=1$ {\bfseries to} $B\in \mathbb{N}$}
 \STATE Sample one set of observations from $N_n(\mu\mathbbm{1},\widehat{\Sigma}_{\boldsymbol{\theta}|H_0}),$ and denote the sample as $\boldsymbol{z}^{(b)}.$
\STATE Estimate $\mathcal{K}_{\boldsymbol{\theta}|H_0}(\cdot),\;\mathcal{K}_{\boldsymbol{\theta}|H_1}(\cdot)$ using $\boldsymbol{z}^{(b)}$ in the same way as in line 3 and denote their covariance matrices scores at $S$ as $\widehat{\Sigma}^{(b)}_{\boldsymbol{\theta}|H_0},\;\widehat{\Sigma}^{(b)}_{\boldsymbol{\theta}|H_1},$ respectively.
\STATE Calculate
    \[\phi^{(b)}=\ell\big(\boldsymbol{z}^{(b)};\widehat{\Sigma}^{(b)}_{\boldsymbol{\theta}|H_1}\big)-\ell\big(\boldsymbol{z}^{(b)};\widehat{\Sigma}^{(b)}_{\boldsymbol{\theta}|H_0}\big).\]
  \ENDFOR
\STATE  \[\text{P-value}=|\{\phi\leq\phi^{(b)}|b\in[1,...,B]\}|/B,\] where $|\cdot|$ is the set size.
\end{algorithmic}
\end{algorithm}

Technical details for Algorithm \ref{Algorithm: parametric}:
\begin{itemize}
    \item Line 3: The kernel parameters can be estimated using various approaches, such as maximum likelihood and restricted maximum likelihood (REML) of  $\boldsymbol{z}.$ More details about kernel parameters estimation can be found in \citet{rasmussen2003gaussian}.
    \item Line 4: $\phi,$ the anisotropic discrepancy measure, can also be written as follows:
  \[\phi=\Big(-\ell\big(\boldsymbol{z};\widehat{\Sigma}_{\boldsymbol{\theta}|H_0}\big)\Big)-\Big(-\ell\big(\widehat{\Sigma}_{\boldsymbol{\theta}|H_1}\big)\Big).\]
  Therefore $\phi$ is the loss function of the null hypothesis minus the loss function of the alternative hypothesis (where the loss function in this case is minus log-likelihood).
  \item Line 5: $B$ controls the P-value resolution. For example, when $B=200$ the P-value resolution is $0.005.$ Also, $\mathrm{Var}(\text{P-value})$ decreases with $B$ (still, commonly the main factor affecting on $\mathrm{Var}(\text{P-value})$ is $\mathrm{Var}(\boldsymbol{z}(S)),$ rather than $B$).
  \item Line 6.: $\mu$ can be estimated by the mean. Alternatively, the data can be normalized, such that $\mu=0.$
  \item Line 10: By definition, the P-value is the probability of rejecting $H_0$ when $H_0$ is true, i.e., $P(\text{rejecting } H_0|H_0).$  Here, the distribution of $\phi$ under the null, $P(\phi|H_0),$ is numerically estimated (rather than derived analytically) and therefore the algorithm's output is $P(\text{rejecting } H_0|\widehat{H_0}).$ This is an inherent property in sampling-based hypothesis testing approaches.
  \end{itemize}

The logic behind Algorithm \ref{Algorithm: parametric} is that $P(\phi|H_0)$ is estimated using parametric bootstrap approach, i.e., by simulating $\{\phi^{(b)}\}_{b\in[1,...,B]}$ using the best empirical parametric representation of the data under the null hypothesis, $N_n(\mu\mathbbm{1},\widehat{\Sigma}_{\boldsymbol{\theta}|H_0}).$ Then, $\phi,$ which is calculated using the data itself is compared to $\{\phi^{(b)}\}_{b\in[1,...,B]}$ in order to derive the P-value. Obviously, this framework involves significance parametric assumptions (normality and a specific kernel function family), however as will be presented in Section \ref{section: Literature Review}, these assumptions are relatively mild, and do not involve arbitrary hyperparameters and asymptotic approximations. Also, the parametric bootstrap hypothesis testing approach allows flexibility in different aspects. First, the statistic $\phi$ can be modified to other loss functions measuring the anisotropic discrepancy magnitude, such as test set error or even prediction errors, e.g., AIC \citep{akaike1974new}, Cp \citep{mallows1973some} and cross-validation \citep{stone1974cross} error types.
Secondly, Algorithm \ref{Algorithm: parametric} allows controlling $B,$ which tradeoffs between the resultant P-value resolution and $\mathrm{Var}(\text{P-value})$ on one hand and the computational cost on the other hand.

We note that our proposed approach is a straight forward application of parametric bootstrap to this problem, however to our knowledge this has not been previously proposed in the context of testing anisotropy. As we show below, it yields strong results even when the assumptions do not strictly apply.

\subsubsection{Parametric Bootstrap with a Non-Parametric Covariance Function}\label{section: Parametric bootstrap with non-parametric covariance}
Besides normality, the main parametric assumption in Algorithm \ref{Algorithm: parametric} is the assumption of kernel function family. This section presents how this assumption can be potentially avoided by estimating the covariance function using a non-parametric regression with monotonicity constraints instead of a kernel function.

Since $y_{i,j}\coloneqq(z_i-\mu)\times (z_j-\mu)$ is an estimator of $C(\boldsymbol{h}_{i,j}),$ then $C(\boldsymbol{h})$ can be estimated by modeling $\boldsymbol{y}=\{y_{i,j}\}_{(i,j)\in([1,...,n],[i,...,n])}$ (which contains $(n+1)\times n/2$ points) as a function of its corresponding location differences, $H=\{\boldsymbol{h}_{i,j}\}_{(i,j)\in([1,...,n],[i,...,n])}.$ The model should be constrained such that it yields non-negative values and the function is monotonically decreasing with respect to the distance (Tobler's law, \cite{tobler1970computer}). An example of a framework that facilitates enforcing monotonic constraints is isotonic regression \citep{stylianou2002dose,luss2014generalized}. Other popular tools that facilitate enforcing monotonic constraints are \href{https://stat.ethz.ch/R-manual/R-patched/library/mgcv/html/mono.con.html}{cubic regression splines model in mgcv package in R software} and \href{https://xgboost.readthedocs.io/en/latest/}{xgboost package in Python software}. Moreover, the idea of replacing the kernel function framework by monotonic regression in other applications is not new, for example see \cite{bjornstad2001nonparametric,wu2003nonparametric,choi2014modeling}.\footnote{Unfortunately, these studies do not supply code for implementation of their algorithm.}

Replacing the kernel function in Algorithm \ref{Algorithm: parametric} by a monotonic regression model requires specifying the suspected anisotropic directional axes $\{\eta_i\}_{i=1}^{R}.$ For example, when implementing isotonic regression in Algorithm \ref{Algorithm: parametric}, the monotonicity constraint in the null hypothesis model is forced on $\{\|\boldsymbol{h}_{i,j}\|\}_{(i,j)\in([1,...,n],[i,...,n])},$ while in the alternative hypothesis is forced on the projections of $H$ on  $\{\eta_i\}_{i=1}^{R}.$
Besides the covariance function estimation, all other parts in Algorithm \ref{Algorithm: parametric} remain the same.

We have invested substantial effort in implementing this approach with non-parametric covariance estimation, trying out the different tools mentioned above with no success. The main reason for this failure seems to be the difficulty of estimating the covariance matrices well. They are consistently poorly estimated and the resulting tests have very poor performance.
Other issues which come up are:
\begin{itemize}
\item The fitted monotonic regression model might be a non strictly positive definite function, therefore the covariance matrix that is plugged-in the likelihood function can be non-positive definite (which make it as a illegitimate covariance matrix). A simple workaround for solving this issue is removing the  eigenvectors with the non-positive eigenvalues from the estimated covariance matrices.
\item As opposed to Section \ref{section: Parametric Bootstrap}, here the anisotropic directional axes must be prespecified.
\item Models enforcing monotonicity are computationally expensive, especially models with multiple covariates, as in the anisotropic model. Taking into account the large sample size in our application, $(n+1)\times n/2,$ then when $n$ is not very small, the running time is long. Therefore, fitting the model using a sub-sample of $\{\boldsymbol{y},H\}$ can be necessary sometimes. Note that since the elements in $\{\boldsymbol{y},H\}$ are correlated, due to the overlap between the $z(\boldsymbol{s}_i)$'s elements  composing $\boldsymbol{y}$ and due to the spatial correlation between $z(\boldsymbol{s}_i)$'s themselves, then the effective sample size is smaller than the nominal one.
\end{itemize}

\subsection{Non-Parametric Rotational Sampling Test}\label{section: Non-Parametric Rotational Sampling Test}
Algorithm \ref{Algorithm: parametric} is based on the normality assumption, which is very common in spatial regression of a continuous dependent variable, however in many cases the deviation from normality is prominent and it is not possible to assume normality \citep{horvath2020testing}.
Here we present a different approach which relies on the rotational symmetry that exists under the null hypothesis, to test it with no distributional assumptions. However, the suspected anisotropic directions have to be specified in advance (exactly or at least approximately). This approach is presented in Algorithm \ref{Algorithm: non-parametric}. For improved readability, we present Algorithm \ref{Algorithm: non-parametric}
where two perpendicular anisotropic directional axes are specified, however as will be demonstrated in Section \ref{section: Real Data}, it can be easily generalized for more than two and non-perpendicular anisotropic directional axes.
\begin{algorithm}[ht!]
\caption{Non-Parametric Rotational Sampling Hypothesis Testing for Anisotropy Detection}
\label{Algorithm: non-parametric}
\begin{algorithmic}[1]
\STATE {\bfseries Input:} $\{S,\boldsymbol{z}(S)\},\;\eta,\;\alpha\in[0,\pi/4).$
\STATE {\bfseries Output:} P-value.
\STATE Derive $\{\boldsymbol{y},H\}$ using $\{S,\boldsymbol{z}(S)\}.$
\STATE  Calculate
\begin{align*}
\phi_{isotropy}&=\underset{\boldsymbol{\theta}|H_0}{\min} \sum_{i=1}^{n}\sum_{j=i}^{n}\big(y_{i,j}-\mathcal{K}_{\boldsymbol{\theta}|H_0}(\|\boldsymbol{h}_{i,j}\|)\big)^2\\
\phi_{anisotropy}&=\underset{\boldsymbol{\theta}|H_1}{\min} \sum_{i=1}^{n}\sum_{j=i}^{n}\big(y_{i,j}-\mathcal{K}_{\boldsymbol{\theta}|H_1}(\boldsymbol{h}_{i,j})\big)^2
\end{align*}
\[\phi=\phi_{isotropy}-\phi_{anisotropy}.\]
\FOR{$b=1$ {\bfseries to} $B\in \mathbb{N}$}
 \STATE Sample random directional axis $\eta^{(b)}$ from  $[\eta+\alpha,\eta+\pi/2-\alpha].$
 \STATE Calculate
\[\phi^{(b)}_{anisotropy}=\underset{\boldsymbol{\theta}|H_1,\eta^{(b)}}{\min} \sum_{i=1}^{n}\sum_{j=i}^{n}\big(y_{i,j}-\mathcal{K}_{\boldsymbol{\theta}|H_1,\eta^{(b)}}(\boldsymbol{h}_{i,j})\big)^2,\]
where $\mathcal{K}_{\boldsymbol{\theta}|H_1,\eta^{(b)}}$ is the anisotropic kernel with the anisotropic directional axes $\{\eta^{(b)},\eta^{(b)}+\pi/2\}$ instead of $\{\eta,\eta+\pi/2\}.$
\[\phi^{(b)}=\phi_{isotropy}-\phi^{(b)}_{anisotropy}.\]
  \ENDFOR
\STATE  \[\text{P-value}=|\{\phi\leq\phi^{(b)}|b\in[1,...,B]\}|/B.\]
\end{algorithmic}
\end{algorithm}

Technical details for Algorithm \ref{Algorithm: non-parametric}:
\begin{itemize}

    \item Line 4.: Since Algorithm \ref{Algorithm: non-parametric} enables avoiding assuming normality, it is more reasonable to use squared errors loss function than a likelihood-based loss function.
    \item Line 6.: In order to increase power, $\alpha$ prevents sampling axes that are close to the specified anisotropic directional axes.

    When the specified anisotropic directional axes are not perpendicular or the number of axes is larger than two, then the space of $\alpha$ is changed accordingly.
    \item Line 7.: Unlike in Algorithm \ref{Algorithm: parametric}, since a new data is not simulated, $\phi_{isotropy}$ can be reused in $\phi^{(b)}.$
\end{itemize}

As in Algorithm \ref{Algorithm: parametric}, Algorithm \ref{Algorithm: non-parametric}
also estimates $P(\phi|H_0)$ using simulation, however here the \textit{directional axes are simulated}, rather than parametric bootstrap samples. In that way the original data can also be used for calculating $\{\phi^{(b)}\}_{b\in[1,...,B]}$ and
distribution assumptions can be avoided. On the other hand, sampling the directional axes forces specifying the suspected anisotropic directional axes (unlike in Algorithm \ref{Algorithm: parametric}).

The requirement for specifying suspected anisotropic directional axes can be partially relaxed, such that ranges that include the suspected anisotropic directional axes are specified instead of the exact values. It can be done by utilizing $\alpha$ for constructing non-overlapping domains of $\{\eta_i\}_{i\in[1,...,R]}$ and $\{\eta^{(b)}\}_{b\in[1,...,B]}.$ For example, in case the suspected anisotropic directional ranges are $\{\eta_1\in[-\alpha,\alpha],\eta_2=\eta_1+\pi/2\},$ then the estimation of $\phi_{anisotropy}$ should also include optimization of $\eta_1.$ Correspondingly, the sampling space in line 6 is $[2\times\alpha,\pi/2-2\times\alpha],$ and the optimization in line 7 is also over $\eta^{(b)*}\in[\eta^{(b)}-\alpha,\eta^{(b)}+\alpha].$

Similarly to Algorithm \ref{Algorithm: parametric}, the kernel can be potentially replaced by non-parametric monotonic regression. In that way, both main parametric assumptions --- normality and kernel structure --- are avoided.

\section{Literature Review}\label{section: Literature Review}
The most common approach for detecting anisotropy is by analyzing popular graphs, such as directional variograms plot and rose diagram \citep{weller2016review}, which compare empirical variograms with respect to different directional exes. Although graphs can be very informative, analyzing them is open to subjective interpretations and therefore cannot be used for deriving objective scientific conclusions. Another drawback in using these graphs, which also appears in many hypothesis testings methods for detecting anisotropy, is their reliance on hyperparameters (e.g., $\delta,$ see Section \ref{section: introduction}) which may affect the conclusion.

A popular hypothesis testing approach for detecting anisotropy is utilizing the asymptotic distribution of the empirical variogram (which was derived by \citet{baczkowski1987approximate}) in order to calculate the P-value of tests whose statistic is based on the empirical variograms. For example, in \cite{guan2004nonparametric}, lags, $[\boldsymbol{h}_1,....,\boldsymbol{h}_K],$ that relate to the suspected anisotropic directional axes are selected, such that under $H_0$ (which assumes isotropy) $AG=0,$ where $G=[\gamma(\boldsymbol{h}_1),....,\gamma(\boldsymbol{h}_K)]$ and A is the contrast matrix, while under $H_1,$ $AG\neq0.$ The statistic of the test is $c(A\widehat{G})^{t}(A\Sigma_GA^{t})^{-1}A\widehat{G},$ where $\widehat{G}=[\widehat{\gamma}(\boldsymbol{h}_1),....,\widehat{\gamma}(\boldsymbol{h}_K)],$ $\Sigma_G$ is the variance of the asymptotic distribution of $A\widehat{G}$ under $H_0$ and $c\in \mathbb{R}.$ The P-value can be derived by the asymptotic distribution of the statistic (which is  $\chi^2$).
The main drawbacks in this method is deriving $\Sigma_G$ and the selection of the smoothing parameter $\delta$ when calculating $\widehat{G}.$ Also, the asymptotic derivations in \cite{guan2004nonparametric} assume specific sampling mechanisms of $S,$ and therefore the test is not valid in many applications which do not follow their sampling mechanisms. \citet{maity2012testing} modify \citeauthor{guan2004nonparametric}'s method for testing isotropy in the covariance function framework by estimating the covariance function at $[\boldsymbol{h}_1,...,\boldsymbol{h}_K]$ using kernel regression. Also, they suggest to estimate the variance matrix of the asymptotic distribution using block bootstrap. Both methods, \cite{guan2004nonparametric} and \cite{maity2012testing}, do not rely on specification of a kernel covariance/variogram function, however their tests rely on selection of hyperparameters --- $\delta$ in \cite{guan2004nonparametric} and the kernel regression function and its bandwidth in \cite{maity2012testing}. \cite{weller2017nonparametric} claims that the resulting tests are very sensitive to the hyperparameters selection.

Detecting anisotorpy correlation structure in the spectral domain, where the  asymptotic derivations are simpler than in the spatial domain, is also common \citep{weller2016review}.
\cite{van2020general} propose a unified framework for detecting different properties of the covariance function, including anisotropy. Similarly to \cite{guan2004nonparametric} and \cite{maity2012testing}, \citeauthor{van2020general}'s test requires pre-specification of the directional axes, however the statistic assessing the deviation from isotropy is based on the periodogram, which is the corresponding object to variogram/covariance function in the spectral domain. Still, some challenges are common with \citeauthor{maity2012testing}'s approach. First, \citeauthor{van2020general} provide an asymptotic result, and it is unclear how their method performs for small sample sizes (the smallest sample size that is used in their simulations is $n=1200$). Also, this method requires selecting a set of lags and number and spacing of frequencies for calculating the statistic.
Besides \cite{van2020general}, there are other papers that analyze covariance function properties, however many of them focus on different types of symmetry \citep{weller2017nonparametric}, which is closely related to isotropy (isotropy does imply symmetry, but symmetry does not imply isotropy).

\section{Simulation Experiments}\label{section: simulation}
This section presents simulation experiments which compare between Algorithm \ref{Algorithm: parametric}, Algorithm \ref{Algorithm: non-parametric} and the algorithm presented by \citet{maity2012testing}, which we refer to as MS.

\paragraph{Distributional setting}
The sample, $\{z(\boldsymbol{s}_i)\}_{i=1}^{n},$ was drawn from a normal distribution with zero mean and the elliptic  exponential kernel covariance function (see expression \ref{Kernel: exponential, anisotropy}) with the parameters $\sigma_{\epsilon}^2=1,\sigma_s^2=1,\eta=0,\lambda_1=1$ and different $\lambda_2$ values. The coordinate values, $\{\boldsymbol{s}_i\}_{i=1}^{n},$ were drawn from the uniform distribution $U(0,1).$ The simulation experiment was executed for different sample sizes, $n=200/500/1000,$ and different anisotropic intensity, $\lambda_2=1/2/5/10$ ($\lambda_2=1$ means an isotropic kernel). Also, the number of resamples $B$ was set to $200.$ For the simulation code see at \textbf{(anonymized for review)}.

\paragraph{Hypotheses definition} In all the three algorithms the null hypothesis assumes isotropy (see expression \ref{Kernel: exponential, isotropy}) and the alternative hypothesis specifies the anisotropic directional axes $\{\eta_1=0,\eta_2=\pi/2\}.$
As presented in Section \ref{section: Parametric Bootstrap}, Algorithm \ref{Algorithm: parametric} enables avoiding specifying anisotropic directional axes, however in order to present a fair comparison between all the algorithms, the hypotheses are defined here in the same way for all the algorithms. For other settings, where the anisotropic directional axes are not specified for Algorithm \ref{Algorithm: parametric}, see Section \ref{section: Real Data}.

\paragraph{Technical details}
For Algorithm \ref{Algorithm: parametric}, the \href{https://scikit-learn.org/stable/modules/gaussian_process.html}{Gaussian Process module} from the sklearn package in python was used to estimate the kernel parameters. For Algorithm \ref{Algorithm: non-parametric}, in order to reduce running time, the test was implemented only on $10,000$ observations which were drawn without replacement from $\{\boldsymbol{y},H\}.$ Also, $\alpha$ was set to $\pi/36$ (which is $5^o$). For MS, SpTest package in R software  \citep{weller2015sptest} was used with the parameters:
$\text{nBoot}=75,\;\text{blockdims}=\{1,1\},\;\text{grid}=\{0.1,0.1\},$ defaulted lags array of $\Lambda=\{\boldsymbol{h}_1=(0.1,0),\;\boldsymbol{h}_2=(0,0.1),\;\boldsymbol{h}_1=(0.1,0.1),\;\boldsymbol{h}_1=(-0.1,0.1)\}$
and the contrast matrix $A=\begin{bmatrix}
1 & -1 & 0&0\\
0 & 0 & 1&-1
\end{bmatrix}.$ Each simulation setting was repeated about $200$ times in order to estimate the P-value distribution.

The following table presents the estimated Type I error and power under significance level of $0.05,$ (i.e., $\widehat{P}(\text{P-value}<0.05)$) for the different settings.

\begin{table}\footnotesize\label{Table: power}
\caption{The table presents empirical power for different settings for significance level of $0.05.$ The $\lambda_2=1$ column is the Type I error estimates.}
\fbox{%
\begin{tabular}{|l|ccc|ccc|ccc|ccc|}
                           & \multicolumn{3}{c|}{$\lambda_2=1$}                              & \multicolumn{3}{c|}{$\lambda_2=2$}                           & \multicolumn{3}{c|}{$\lambda_2=5$}                            & \multicolumn{3}{c|}{$\lambda_2=10$}       \\
Algorithm\textbackslash{}n & 200                        & 500                        & 1000  & 200                       & 500                       & 1000 & 200                       & 500                       & 1000  & 200                       & 500                       & 1000  \\
Algorithm \ref{Algorithm: parametric}                         & 0.05 & 0.04 & 0.05 & 0.11	&0.31	&0.43
& 0.47	&0.70	&0.94&0.65&0.89&	0.99 \\
Algorithm \ref{Algorithm: non-parametric}                         &0.14&	0.07&	0.07 & 0.05&0.15	&0.19& 0.16&	0.29&	0.24
& 0.27&	0.30&	0.33
  \\
MS  &0.02&	0.04&	0 & 0.06	&0.04 &0.05 & 0.02	&0.10 &0.10&  0.07&	0.10&	0.17  \\
\end{tabular}}
\end{table}

As we can see, Algorithm \ref{Algorithm: parametric} has the highest power. It is not surprising, since Algorithm \ref{Algorithm: parametric} makes a stronger distributional assumption --- normality of $\boldsymbol{z}$ --- which holds in our setting. Also, as expected, the power of all the three algorithms increases with the sample size and $\lambda_2,$ however the power of MS is relatively low in our settings. Also, as we can see in the $\lambda_2=1$ column, Algorithm \ref{Algorithm: parametric} controls the Type I error properly for significance level of $0.05,$ Algorithm \ref{Algorithm: non-parametric} also controls the Type I error properly for moderate and large sample (but not for $n=200$) and MS does not control Type I error properly.

Figure \ref{Figure: power for n=500} presents the empirical cumulative distribution of the P-value of the three algorithms for $n=500$ and for $\lambda_2=1/5/10.$
\begin{figure}[h]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/allCurvesdate2021-04-28.pdf}
         \caption{P-value empirical cumulative distribution for Algorithms \ref{Algorithm: parametric}, \ref{Algorithm: non-parametric} and MS for $n=500$ and different  $\lambda_2$ values.}\label{Figure: power for n=500}
    \end{figure}

As we can see, Algorithms \ref{Algorithm: parametric} has a higher power than the other algorithms for any significance level. Still, Algorithms \ref{Algorithm: non-parametric} also has a much higher power than MS.

\section{Real Data Analysis}\label{section: Real Data}
This section analyzes and demonstrates implementation of algorithms \ref{Algorithm: parametric}, \ref{Algorithm: non-parametric} and MS in various use cases using four datasets that were introduced by \cite{hohn1998geostatistics}. The datasets can be found in the \href{http://www.wvgs.wvnet.edu/www/geostat/GeostatPetGeol.html}{book's website}. The code can be found at \textbf{(anonymized for review)}.

\paragraph{Granny Creek Field Data}
The Granny Creek Field dataset contains $181$ measurements of a sandstone base elevation in Granny Creek Field, central West Virginia.
Figure \ref{Figure: Granny Creek Field} presents the analyzed standard score elevation measurements, $\boldsymbol{z}(S),$ and their standardized coordinates, $S.$ Based on prior geographical knowledge, \citeauthor{hohn1998geostatistics} claims that the expected anisotropic directional axes are $\{0,\pi/2\}.$ After analyzing the directional variograms graph, \citeauthor{hohn1998geostatistics} suggests that the anisotropic directional axes are $\{\eta_1=10\times\pi/18,\eta_2=\pi/9\}$ (10 and 100 degrees). Here we use Algorithms \ref{Algorithm: parametric}, \ref{Algorithm: non-parametric} and MS for testing whether $\{\eta_1=0,\eta_2=\pi/2\}$ are anisotropic directional axes. Using exponential kernel and setting $B=200,$ the P-value of Algorithm \ref{Algorithm: parametric} is smaller than $0.005$ and the P-value of Algorithm \ref{Algorithm: non-parametric} with $\alpha=\pi/36$ is $0.025.$ Figure \ref{Figure: Granny Creek Field} compares the $\phi$ value with $\{\phi^{(b)}\}_{b=1}^{200}$ for both algorithms. Implementing MS algorithm using the same setting as in Section \ref{section: simulation}, yields $\text{P-value}=0.011.$ Thus, all three testing approaches support anisotropy.

\begin{figure}[h]
    \centering
        \includegraphics[width=1\linewidth]{Figures/GrannyResultsAll2021-06-06.pdf}
         \caption{Granny Creek Field. The left figure is the scatter plot of the Granny Creek Field dataset after scaling, where the color intensity indicates the elevation value. The middle and the right figures compare the $\phi$ value (in black) with $\{\phi^{(b)}\}_{b=1}^{200}$ for Algorithm \ref{Algorithm: parametric} and Algorithm \ref{Algorithm: non-parametric}, respectively.}\label{Figure: Granny Creek Field}
    \end{figure}

\paragraph{Mississippian Sandstone Data}
The Mississippian Sandstone dataset contains
$348$ measurements of subsea depth of a mississippian-age reservoir sandstone base in Ritchie County, West Virginia. Besides standard scaling, an outlier observation was removed from the analysis. The data is presented in Figure \ref{Figure: Mississippian Sandstone}).

Based on prior geographical knowledge, \citeauthor{hohn1998geostatistics} claims that the correlation in the $\{\pi/4,\;3\times\pi/4\}$ directional axes is suspected to be different than the correlation in the $\{\pi,\;\pi/2\}$ directional axes. Therefore, the elliptical transformation of the anisotropic kernel (see equation \ref{Kernel: exponential, anisotropy}) is:
\[
A=
\begin{bmatrix}
\frac{1}{\lambda_1} & 0& 0& 0 \\
0 & \frac{1}{\lambda_1} & 0& 0\\
0 & 0 &\frac{1}{\lambda_2}& 0\\
0 & 0 &0&\frac{1}{\lambda_2}\\
\end{bmatrix}
\begin{bmatrix}
\cos(0) & \sin(0) \\
\cos(\frac{\pi}{2}) & \cos(\frac{\pi}{2}) \\
\cos(\frac{\pi}{4}) & \sin(\frac{\pi}{4}) \\
\cos(\frac{3\pi}{4}) & \cos(\frac{3\pi}{4}) \\
\end{bmatrix},
\]
i.e., $\lambda_1$ is the length-scale parameter of $\{\eta_3=0,\eta_4=\pi/2\}$ and $\lambda_2$ is the length-scale parameter of $\{\eta_1=\pi/4,\eta_2=3\times\pi/4\}.$

The P-value for Algorithm \ref{Algorithm: parametric} is $0.93$ and for Algorithm \ref{Algorithm: non-parametric}  with $\alpha=\pi/72$ is $0.7.$
Figure \ref{Figure: Mississippian Sandstone} compares the $\phi$ value with $\{\phi^{(b)}\}_{b=1}^{200}.$ Implementing MS algorithm with
contrast matrix, $A=\begin{bmatrix}
1 & 0 &-1&0\\
0 & 1 & 0&-1\\
1 & 1 & -1&-1\\
\end{bmatrix},$ (other parameters are the same as in Section \ref{section: simulation}), gives $\text{P-value}=0.19.$

\begin{figure}[h]
    \centering
        \includegraphics[width=1\linewidth]{Figures/MississippianResultsAll2021-06-06.pdf}
         \caption{Mississippian Sandstone. The left figure is the scatter plot of the Mississippian Sandstone dataset after scaling, where the color intensity indicates the depth value. The middle and the right figures compare the $\phi$ value (in black) with $\{\phi^{(b)}\}_{b=1}^{200}$ for Algorithm \ref{Algorithm: parametric} and Algorithm \ref{Algorithm: non-parametric}, respectively.}\label{Figure: Mississippian Sandstone}
    \end{figure}

\paragraph{Gas Potential in Barbour County Data} This data contains  $674$ measurements of gas initial potential, producing from upper Devonian sandstones and siltstones in Barbour County, West Virginia. Based on \citeauthor{hohn1998geostatistics}'s suggestion, in order to satisfy the normality assumption (in Algorithm \ref{Algorithm: parametric}), the standard score of the log of the gas initial potential is analyzed. Figure \ref{Figure: Barbour County Data} presents the data before and after pre-processing. In this example there is no prior information about suspected anisotropic directional axes. Implementing Algorithm \ref{Algorithm: parametric} without specifying suspected anisotropic directional axes (however assuming that there are two anisotropic directional axes) gives $\text{P-value}=0.015.$ Implementing Algorithms \ref{Algorithm: non-parametric} and specifying defaulted anisotropic directional axes of $\{\eta_1=0,\eta_2=\pi/2\}$ with $\alpha=\pi/36$ gives $\text{P-value}<0.005.$ A comparison of the $\phi$ value with $\{\phi^{(b)}\}_{b\in[1,...,200]}$ is presented in Figure \ref{Figure: Barbour County Data}. Implementing MS algorithm using the same setting given in Section \ref{section: simulation}, gives $\text{P-value}=0.01.$

\begin{figure}[h]
    \centering
        \includegraphics[width=1\linewidth]{Figures/BarbourResultsAll2021-06-06.pdf}
         \caption{Gas potential in Barbour County. The left figure is the scatter plot of the gas potential in Barbour County dataset after scaling, where the color intensity indicates the gas potential. The middle and the right figures compare the $\phi$ value (in black) with $\{\phi^{(b)}\}_{b=1}^{200}$ for Algorithm \ref{Algorithm: parametric} and Algorithm \ref{Algorithm: non-parametric}, respectively.}\label{Figure: Barbour County Data}
    \end{figure}

\paragraph{Greenbrier Limestone}
The Greenbrier Limestone dataset contains measurements at $2,335$ tops of the Greenbrier
Limestone. Based on graphical analysis,
\citeauthor{hohn1998geostatistics} claims that the anisotropic directional axes are \{$\eta_1=\pi/3,\;\eta_2=2\times\pi/3$\}. Since using graphical analyses for selecting the suspected anisotorpic directional axes of the alternative hypothesis violates our hypothesis testing scheme (as well as many other standard hypothesis testing schemes), here we use the default anisotropic directional axes, $\{\eta_1=0,\eta_2=\pi/2\}.$ Also, unlike in the analysis of the Gas Potential in Barbour Country dataset, these axes will also be specified for Algorithms \ref{Algorithm: parametric}. The P-value of Algorithm \ref{Algorithm: parametric} is $0.005$ and the P-value of Algorithm \ref{Algorithm: non-parametric} with $\alpha=\pi/36$ is $0.3.$ Figure \ref{Figure: Greenbrier Limestone} compares the $\phi$ value with $\{\phi^{(b)}\}_{b\in[1,...,B]}.$ Implementing MS algorithm using the same setting given in Section \ref{section: simulation} but with $grid=[0.5,0.5]$ (due to convergence issue), gives $\text{P-value}=0.12.$ Therefore, as we can see, while Algorithm \ref{Algorithm: parametric} rejects the null hypothesis with significance level of $0.05,$ Algorithm \ref{Algorithm: parametric} and MS do not reject.

\begin{figure}[h]
    \centering
        \includegraphics[width=1\linewidth]{Figures/GreenbrierResultsAll2021-06-06.pdf}
         \caption{Greenbrier limestone. The left figure is the scatter plot of the greenbrier limestone dataset after scaling, where the color intensity indicates the elevation values. The middle and the right figures compare the $\phi$ value (in black) with $\{\phi^{(b)}\}_{b=1}^{200}$ for Algorithm \ref{Algorithm: parametric} and Algorithm \ref{Algorithm: non-parametric}, respectively.}\label{Figure: Greenbrier Limestone}
    \end{figure}

\section{Summary and Conclusions}
In this paper we introduce and analyze parametric and non-parametric hypothesis testing algorithms for detecting anisotropic covariance structure. In both algorithms, the statistic distribution under the null hypothesis is obtained using resampling mechanisms, while asymptotic parametric approximations and hyperparameters selection, which are the Achilles' heel of other methods, are avoided.

Both algorithms propose general frameworks that can also be applied when the covariance function is estimated using a non-parametric monotonic regression model, which enables avoiding assuming a parametric kernel function. However, as explained in Section \ref{section: Parametric bootstrap with non-parametric covariance}, modeling the covariance function using a non-parametric monotonic regression model is a difficult task and requires further research.

The numerical results demonstrate the high power that is obtained by the proposed algorithms, and also demonstrates how the algorithms can easily be adjusted in various settings. In our view, the simple resampling-based methods we propose here should be considered as leading practical approaches for testing anisotropy.