\subsection{Proof: Linear Functions}

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda \textrm{geo}(m, m')] \; \; \textrm{d} m' \textrm{d} m $$

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma ||x - m||^2 - \gamma ||y - m'||^2 - \lambda ||m - m'||^2] \; \; \textrm{d} m' \textrm{d} m $$

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma (x - m)^\intercal (x - m) - \gamma (y - m')^\intercal (y - m') - \lambda (m - m')^\intercal (m - m')] \; \; \textrm{d} m' \textrm{d} m $$

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma x^\intercal x + 2 \gamma x^\intercal m - \gamma m^\intercal m     -\gamma y^\intercal y + 2 \gamma y^\intercal m' - \gamma m'^\intercal m'     -\lambda m^\intercal m + 2 \lambda m^\intercal m' - \lambda m'^\intercal m'  ] \; \; \textrm{d} m' \textrm{d} m $$

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma x^\intercal x -\gamma y^\intercal y  + 2 \gamma x^\intercal m + 2 \gamma y^\intercal m' + 2 \lambda m^\intercal m' - (\gamma + \lambda) m^\intercal m  - (\gamma+\lambda) m'^\intercal m'  ] \; \; \textrm{d} m' \textrm{d} m $$

$$k(x,y) =  \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\gamma x^\intercal x -\gamma y^\intercal y  + 2 \gamma x^\intercal m + 2 \gamma y^\intercal m'] \exp[ 2 \lambda m^\intercal m' - (\gamma + \lambda) m^\intercal m  - (\gamma+\lambda) m'^\intercal m'  ] \; \; \textrm{d} m' \textrm{d} m $$

~

~

------------------------------
Working Backwards..

$$k(x,y) =  \exp[-\gamma ||x-y||^2] \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}}   \exp[- 2 \gamma x^\intercal y + 2 \gamma x^\intercal m + 2 \gamma y^\intercal m'] \exp[ 2 \lambda m^\intercal m' - (\gamma + \lambda) m^\intercal m  - (\gamma+\lambda) m'^\intercal m'  ] \; \; \textrm{d} m' \textrm{d} m $$

\subsection{Proof: Triangle Proof, simple case}

Consider a straight line $\mathcal{M}$, which we take to be a 1-dimensional manifold embedded in $\mathbb{R}^2$, and two points on the line $A, B$. We take a second straight line $\mathcal{G}$ which intersects $\mathcal{M}$ at $A$, and a third line perpendicular to $\mathcal{M}$ that intersects $\mathcal{M}$ at $B$ and intersects $\mathcal{G}$ at some point $C$. We define $\epsilon = ||C - B||$ and $\delta = ||B - A||$.

Assume we that we can take a finite number of samples from each manifold, i.e. $M = \{ x:x \in \mathbb{R}^2, x \in \mathcal{M} \}$ and $G = \{ x:x \in \mathbb{R}^2, x \in \mathcal{G} \}$. We define the two kernels $k_{M}(x,y)$ and $k_{G}(x,y)$:

$$k_{M}(x,y) =  \frac{1}{|M|^2}\sum_{m \in M} \sum_{m' \in M} \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda \textrm{geo}(m, m')]$$

We want to show that as $\epsilon$ increases, $\frac{k_M(x,y)}{k_G(x,y)}$ decreases for any $x,y \in \mathbb{R}^2$ s.t. $x,y$ are on the same side of the line.

$$\frac{k_M(x,y)}{k_G(x,y)} = \frac{\sum_{m \in M} \sum_{m' \in M} \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda \textrm{geo}(m, m')]}
{\sum_{g \in G} \sum_{g' \in G} \exp[-\gamma (||x - g||^2 + ||y - g'||^2) - \lambda \textrm{geo}(g, g')]}$$

Since $\mathcal{M}$ and $\mathcal{G}$ are straight lines, $\textrm{geo}(x,y) = ||x - y||$:

$$= \frac{\sum_{m \in M} \sum_{m' \in M} \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda ||m-m'||]}
{\sum_{g \in G} \sum_{g' \in G} \exp[-\gamma (||x - g||^2 + ||y - g'||^2) - \lambda ||g-g'||]}$$

Define $a(m) = \frac{||m-A||}{||B-A||}(C-B)$ and $\alpha(m) = ||a(m)||$. We can then reparametrize $g = m + a(m)$:

$$ = \sum_{m \in M} \sum_{m' \in M} \frac{ \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda ||m - m'||]}
{\exp[-\gamma (||x - m - a(m)||^2 + ||y - m' - a(m')||^2) - \lambda ||m + a(m)-m' - a(m')||]}$$

\begin{multline*}
    = \sum_{m \in M} \sum_{m' \in M} \exp[-\gamma (||x - m||^2 - ||x - m - a(m)||^2] \\
        \times \exp[-\gamma(||y - m'||^2- ||y - m' - a(m')||^2)]\\
        \times \exp[- \lambda(||m - m'|| - ||m + a(m)-m' - a(m')||) ]
\end{multline*}

Third line:

$$\exp[- \lambda(||m - m'|| - ||m + a(m)-m' - a(m')||)]$$
$$=\exp[- \lambda(||m - m'|| - \sqrt{||m-m'||^2 + (\frac{\epsilon}{\delta}||m-m'||)^2})]$$
$$=\exp[- \lambda(||m - m'|| - \sqrt{||m-m'||^2 (1+ \frac{\epsilon^2}{\delta^2})})]$$
$$=\exp[- \lambda||m - m'|| (1- \sqrt{1+ \frac{\epsilon^2}{\delta^2}})]$$

First line:
$$\exp[-\gamma (||x - m||^2 - ||x - m - a(m)||^2)]$$
$$=\exp[-\gamma (||x - m||^2 - ||x-m||^2 + a(m)^\intercal a(m) - 2a(m)^\intercal x  - 2 a(m)^\intercal m)]$$
$$=\exp[-\gamma (a(m)^\intercal a(m) - 2a(m)^\intercal x  - 2 a(m)^\intercal m)]$$
$$=\exp[-\gamma a(m)^\intercal a(m) -  2 \gamma a(m)^\intercal (x  - m)]$$
$$\leq \exp[-\gamma ||a(m) || -  2 \gamma ||a(m)|| \, || x  - m ||]$$
$$\leq \exp[-2 \gamma ||a(m) || (\frac12 + || x  - m ||)]$$

\iffalse
We can approximate the manifold using finite samples. Here is the discrete form of our kernel function. We also normalize the weight vectors with their L2 norm.

$$k(x,y) := \sum_{m \in \mathcal{M}} \sum_{m' \in \mathcal{M}} \frac{w(x, m)}{\sqrt{\sum_{m \in \mathcal{M}} w(x,m)^2}} \, \frac{w(x, m')}{\sqrt{\sum_{m' \in \mathcal{M}} w(x,m')^2}}\, g(m, m')$$

\fi
\iffalse
\subsection{Method: Weighted Geodesic between the closest projection}

Given a point $x$, we define it's closest projection on the manifold as $z_x$.

We define a kernel function $K(x,y)$ as follows:
$$K(x,y) := \int_{z \in S^2} w(x, z) g(z, z_y) + \int_{z \in S^2} w(y, z) g(z, z_x)$$
\fi

\subsection{Proof: Euler-Lagrange}

$$k(x,y) =  \int^{A}_{B}  \int^{A}_{B}   \exp[-\gamma (||x - m||^2 + ||y - m'||^2) - \lambda \textrm{geo}(m, m')] \; \textrm{d} m \textrm{d} m' $$

$$ =  \int^{A}_{B}  \exp[-\gamma ||y - m'||^2] \int^{A}_{B}   \exp[-\gamma ||x - m||^2  - \lambda \textrm{geo}(m, m')] \; \textrm{d} m \textrm{d} m' $$

Consider the inner integral over $m \in \mathcal{M}$ which is taken w.r.t. some fixed $m'$. Using the implicit function theorem, define a function $f$ s.t. $F(t,f(t)) = 0$.

$$J[f] = \int^{A}_{B} L(t, f(t), f'(t)) $$
$$s.t. \quad L(t, f(t), f'(t)) = \exp[-\gamma ||x - f(t)||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s] \sqrt{1+f'(t)^2 }$$

We consider $J$ a functional operating on $f$. Using Euler-Lagrange: $\frac{\partial L}{\partial f} = \frac{\textrm{d}}{\textrm{d}t} \frac{\partial L}{\partial f'}$

$$\frac{\partial L}{\partial f} =  \exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s] \sqrt{1+f'(t)^2 } [-\lambda \sqrt{1+f'(f(t))^2} -2\delta (x_2 - f(t))]$$

$$=  -[\lambda (1+f'(f(t))^2)^{\frac12} + 2\delta(x_2-f(t))](1+f'(t)^2)^{\frac12}\exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s]$$

$$\frac{\partial L}{\partial f'} = \exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s]\frac{f'(t)}{\sqrt{1+f'(t)^2}}$$

\begin{multline*}
    \frac{\textrm{d}}{\textrm{d}t} \frac{\partial L}{\partial f'} = \exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s]\frac{f'(t)}{\sqrt{1+f'(t)^2}}(-\lambda f'(t) \sqrt{1+f'(t)^2}) \\
    + \exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s] \frac{f''(t) f'(t)^2}{(1+f'(t)^2)^{\frac{3}{2}}}
\end{multline*}

$$=\exp[-\gamma ||x - \{t,f(t)\}||^2  - \lambda \int^{f(t)}_{m'} \sqrt{1+f'(s)^2} \; \textrm{d}s][-\lambda f'(t)^2 + \frac{f''(t) f'(t)^2}{(1+f'(t)^2)^{\frac{3}{2}}}]$$

Setting $\frac{\partial L}{\partial f} = \frac{\textrm{d}}{\textrm{d}t} \frac{\partial L}{\partial f'}$:

$$ -[\lambda (1+f'(f(t))^2)^{\frac12} + 2\delta(x_2-f(t))](1+f'(t)^2)^{\frac12} = -\lambda f'(t)^2 + \frac{f''(t) f'(t)^2}{(1+f'(t)^2)^{\frac{3}{2}}}$$

$$ -[\lambda (1+f'(f(t))^2)^{\frac12} + 2\delta(x_2-f(t))](1+f'(t)^2)^{2} = -\lambda f'(t)^2 (1+f'(t)^2)^{\frac{3}{2}}  + f''(t) f'(t)^2$$

\clearpage
\subsection{Proof: Euler-Lagrange v2}

$$k(x,y) =  \int_{m \in \mathcal{M}}  \int_{m' \in \mathcal{M}}  \exp[-\rho (||x - m||^2 + ||y - m'||^2) - \lambda d_{geo}(m, m')] \; \textrm{d} m \textrm{d} m' $$

Assume that $\mathcal{M}$ is a 1-dimensional manifold embedded in $\mathbb{R}^2$. We define an interval $(a,b) \subset \mathbb{R}$ and a smooth function $\gamma:(a,b) \rightarrow \mathcal{M}$. We reparametrize the kernel function:

$$\int^{a}_{b}  \int^{a}_{b}   \exp[-\rho (||x - \gamma(t)||^2 + ||y - \gamma(t')||^2) - \lambda \int_{t'}^{t} \sqrt{g(\dot \gamma(s), \dot \gamma(s))} \; \textrm{d}s ]\; ||\dot \gamma(t)|| \; ||\dot \gamma(t')|| \; \textrm{d} t \textrm{d} t' $$

$$ =  \int^{a}_{b} \exp[-\rho ||y - \gamma(t')||^2 ] ||\dot \gamma(t')|| \int^{a}_{b}   \exp[-\rho ||x - \gamma(t)||^2 - \lambda \int_{t'}^{t} \sqrt{\sum_{i,j} g_{ij}\dot \gamma(s)_i \dot \gamma(s)_j} \; \textrm{d}s ] ||\dot \gamma(t)||\; \textrm{d} t \textrm{d} t' $$

$\gamma^i$ indicates the $i^{th}$ component of vector $\gamma$. Since we assume that the decision boundary is an embedding in $\mathbb{R}^2$, the Riemannian Metric is Euclidean; i.e.

\[ g_{ij} =
\begin{cases}
    1 & i = j \\
    0 & i \neq j
 \end{cases}
\]

We rewrite the geodesic distance:

$$ =  \int^{a}_{b} \exp[-\rho ||y - \gamma(t')||^2 ] \int^{a}_{b}   \exp[-\rho ||x - \gamma(t)||^2 - \lambda \int_{t'}^{t} ||\dot \gamma (s)|| \; \textrm{d}s ] \; ||\dot \gamma(t)|| \; \textrm{d} t \textrm{d} t' $$

Consider the inner integral with respect to $t$, with some fixed $t'$:

$$J[\gamma] = \int^{a}_{b} L(t, \gamma(t), \dot \gamma(t)) $$

$$s.t. \quad L(t, \gamma(t), \dot \gamma(t)) = \exp[-\rho ||x - \gamma(t)||^2 - \lambda \int_{t'}^{t} ||\dot \gamma(s)|| \; \textrm{d}s ] \; || \dot \gamma(t) ||$$

For simplicity, define:

\begin{equation}
    L_1 =  \exp[-\rho ||x - \gamma(t)||^2] \quad \quad  L_2 =  \exp[- \lambda \int_{t'}^{t} ||\dot \gamma(s)|| \; \textrm{d}s ]  \quad \quad L_3 =  ||\dot \gamma(t) ||
\end{equation}

This results in the Lagrangian:

\begin{equation}
    L(t, \gamma(t), \dot \gamma(t)) = L_1 L_2 L_3
\end{equation}

Taking derivative $\frac{\partial L}{\partial \gamma_i}$:
$$\frac{\partial L}{\partial \gamma} = L_2 L_3 \frac{\partial L_1}{\partial \gamma} $$
\begin{equation} \label{eqn:EulerLagrangeLHS}
    =2 \rho L_1 L_2 L_3  (x_i-\gamma_i(t))
\end{equation}

\begin{comment}
Taking derivative $\frac{\partial L_2}{\partial \dot \gamma_i}$:
$$\frac{\partial L_2}{\partial \dot \gamma_i}$$

Taking derivative $\frac{\partial L_3}{\partial \dot \gamma_i}$:
\begin{equation}
    \frac{\partial L_3}{\partial \dot \gamma_i} = \frac{\sum_{i=1}^{d-1} \dot \gamma_i (t) \ddot \gamma_i  (t)}{||\dot \gamma(t)||}
\end{equation}
\end{comment}

Taking derivative $\frac{\partial L}{\partial \dot \gamma_i}$:

\begin{equation}
    \frac{\partial L}{\partial \dot \gamma_i} = L_1 L_3 \frac{\partial L_2}{\partial \dot \gamma_i} + L_1 L_2 \frac{\partial L_3}{\partial \dot \gamma_i}
\end{equation}

Taking derivative $\frac{\textrm{d}}{\textrm{d} t} (\frac{\partial L}{\partial \dot \gamma_i})$:

$$= \frac{\textrm{d} L_1}{\textrm{d} t} L_3 \frac{\partial L_2}{\partial \dot \gamma_i}
+ L_1 \frac{\textrm{d} L_3}{\textrm{d} t} \frac{\partial L_2}{\partial \dot \gamma_i}
+ L_1 L_3 \frac{\textrm{d}}{\textrm{d} t}(\frac{\partial L_2}{\partial \dot \gamma_i})
+ \frac{\textrm{d} L_1}{\textrm{d} t} L_2 \frac{\partial L_3}{\partial \dot \gamma_i}
+ L_1 \frac{\textrm{d} L_2}{\textrm{d} t}  \frac{\partial L_3}{\partial \dot \gamma_i}
+ L_1 L_2 \frac{\textrm{d}}{\textrm{d} t}(\frac{\partial L_3}{\partial \dot \gamma_i})
$$

Using Euler-Lagrange for $L_2$, we get  $\frac{\textrm{d}}{\textrm{d} t} (\frac{\partial L_2}{\partial \dot \gamma_i}) = \frac{\partial L_2}{\partial \gamma} \; \Rightarrow \; \frac{\textrm{d}}{\textrm{d} t} (\frac{\partial L_2}{\partial \dot \gamma_i}) = 0$. Similarly,  $\frac{\textrm{d}}{\textrm{d} t} (\frac{\partial L_3}{\partial \dot \gamma_i}) = 0$.

\begin{equation} \label{eqn:EulerLagrangeRHS}
    = \frac{\textrm{d} L_1}{\textrm{d} t} L_3 \frac{\partial L_2}{\partial \dot \gamma_i}
    + L_1 \frac{\textrm{d} L_3}{\textrm{d} t} \frac{\partial L_2}{\partial \dot \gamma_i}
    + \frac{\textrm{d} L_1}{\textrm{d} t} L_2 \frac{\partial L_3}{\partial \dot \gamma_i}
    + L_1 \frac{\textrm{d} L_2}{\textrm{d} t}  \frac{\partial L_3}{\partial \dot \gamma_i}
\end{equation}

We use Euler-Lagrange and combine Eq. \ref{eqn:EulerLagrangeRHS} with Eq. \ref{eqn:EulerLagrangeLHS}:
\begin{equation} \label{eqn:EulerLagrangeCombined}
    2 \rho L_1 L_2 L_3  (x-\gamma(t)) = \frac{\textrm{d} L_1}{\textrm{d} t} L_3 \frac{\partial L_2}{\partial \dot \gamma_i}
    + L_1 \frac{\textrm{d} L_3}{\textrm{d} t} \frac{\partial L_2}{\partial \dot \gamma_i}
    + \frac{\textrm{d} L_1}{\textrm{d} t} L_2 \frac{\partial L_3}{\partial \dot \gamma_i}
    + L_1 \frac{\textrm{d} L_2}{\textrm{d} t}  \frac{\partial L_3}{\partial \dot \gamma_i}
\end{equation}
---------------------------------------------

Taking derivative $\frac{\partial L_2}{\partial \dot \gamma_i}:$

\begin{equation}
    \frac{\partial L_2}{\partial \dot \gamma_i} = - \lambda L_2 \int_{t'}^{t} \frac{\partial}{\partial \dot \gamma_i} (||\dot \gamma(s)||)  \; \textrm{d}s
\end{equation}

\begin{equation}
    = - \lambda L_2 \int_{t'}^{t} \frac{ \dot \gamma_i (s)}{||\dot \gamma (s)||}  \; \textrm{d}s
\end{equation}
\begin{equation}
    := -\lambda L_2 A
\end{equation}

Taking derivative $\frac{\partial L_3}{\partial \dot \gamma_i}:$
\begin{equation}
    \frac{\partial L_3}{\partial \dot \gamma_i} = \frac{ \dot \gamma_i (t)}{||\dot \gamma (t)||}
\end{equation}

Taking derivative $\frac{\textrm{d} L_1}{\textrm{d} t}:$
\begin{equation}
    \frac{\textrm{d} L_1}{\textrm{d} t} = 2 \rho L_1  \sum_{i=1}^{2} \dot \gamma_i(t) (x_i-\gamma_i(t))
\end{equation}

Taking derivative $\frac{\textrm{d} L_2}{\textrm{d} t}:$
\begin{equation}
    \frac{\textrm{d} L_2}{\textrm{d} t} = -\lambda L_2 ||\dot \gamma(t) ||
\end{equation}

Taking derivative $\frac{\textrm{d} L_3}{\textrm{d} t}:$
\begin{equation} \label{eqn:EulerLagrangeDerivativeLast}
    \frac{\textrm{d} L_3}{\textrm{d} t} = \frac{ \sum_{i=1}^{2} \dot \gamma_i (t) \ddot \gamma_i (t) }{||\dot \gamma (t)||}
\end{equation}
\begin{comment}
\begin{equation}
    = \frac{ \dot \gamma(t)^{\intercal} \ddot \gamma(t)  }{||\dot \gamma (t)||}
\end{equation}
\end{comment}
---------------------------------------------

Combining Eq. \ref{eqn:EulerLagrangeCombined} - Eq. \ref{eqn:EulerLagrangeDerivativeLast}:

\begin{align*}
    2 \rho L_1 L_2 L_3  (x_i-\gamma_i(t)) = [2 \rho L_1  \sum_{i=1}^{2} \dot \gamma_i(t) (x_i-\gamma_i(t))] L_3 [- \lambda L_2 A] \\
    + L_1 \frac{ \sum_{i=1}^{2} \dot \gamma_i (t) \ddot \gamma_i (t) }{||\dot \gamma (t)||} [- \lambda L_2 A] \\
    + [2 \rho L_1  \sum_{i=1}^{2} \dot \gamma_i(t) (x_i-\gamma_i(t))] L_2 [\frac{ \dot \gamma_i (t)}{||\dot \gamma (t)||}] \\
    + L_1 [-\lambda L_2 ||\dot \gamma(t) ||] [\frac{ \dot \gamma_i (t)}{||\dot \gamma (t)||}]
\end{align*}

\begin{align*}
    2 \rho L_3  (x_i-\gamma_i(t)) = - 2 \lambda \rho L_3 A  \sum_{i=1}^{2}[ \dot \gamma_i(t) (x_i-\gamma_i(t))] \\
    - \lambda \frac{ \sum_{i=1}^{2} \dot \gamma_i (t) \ddot \gamma_i (t)  }{||\dot \gamma (t)||} A \\
    + 2 \rho \frac{ \dot \gamma_i (t)}{||\dot \gamma (t)||}  \sum_{i=1}^{2} [\dot \gamma_i(t) (x_i-\gamma_i(t))] \\
    -\lambda \dot \gamma_i (t)
\end{align*}

\begin{align*}
    \sum_{i=1}^{2} \dot \gamma_i (t) \ddot \gamma_i (t) = -2 \rho  (x_i-\gamma_i(t)) ||\dot \gamma(t)||^2 \frac{1}{\lambda A}\\
    - 2 \rho ||\dot \gamma(t)||^2 \sum_{i=1}^{2}[ \dot \gamma_i(t) (x_i-\gamma_i(t))] \\
    + 2 \rho \dot \gamma_i (t) \frac{1}{\lambda A} \sum_{i=1}^{2}[ \dot \gamma_i(t) (x_i-\gamma_i(t))] \\
    - \dot \gamma_i (t) ||\dot \gamma(t)|| \frac{1}{ A}
\end{align*}

\begin{align*}
    \sum_{i=1}^{2} \dot \gamma_i (t) \ddot \gamma_i (t) = \frac{\rho}{\lambda}[ -2  (x_i-\gamma_i(t)) ||\dot \gamma(t)||^2 \frac{1}{ A}
    +2  \dot \gamma_i (t) \frac{1}{ A} \sum_{i=1}^{2}[ \dot \gamma_i(t) (x_i-\gamma_i(t))] ]\\
    - \rho[ 2 ||\dot \gamma(t)||^2 \sum_{i=1}^{2}[ \dot \gamma_i(t) (x_i-\gamma_i(t))]] \\
    - \dot \gamma_i (t) ||\dot \gamma(t)|| \frac{1}{ A}
\end{align*}

\emph{Notes.} A small $\rho$ (weighting) or $\lambda$ (geodesic) indicates that points far away from each other can still have a large impact. I.E. the gaussian distribution has high variance. As these parameters increase, the "Neighborhood" of relevat points shrink, and at $\infty$, the only neighborhood is effectively zero. When these parameters are either zero or infinity, the actual distance between the two points don't really matter since it's the same for any two points.

\begin{comment}
Taking derivatives:

$$\frac{\partial L}{\partial \dot \gamma^i} =  \exp[-\rho ||x - \gamma(t)||^2 - \lambda \int_{t'}^{t} ||\dot \gamma(s)||  \; \textrm{d}s] $$

$$\frac{\partial L}{\partial \gamma} =  2 \rho ( x-\gamma(t)) \exp[-\rho ||x - \gamma(t)||^2 - \lambda \int_{t'}^{t} ||\dot \gamma(s)||  \; \textrm{d}s] \dot \gamma$$
\end{comment}

\clearpage
\subsection{Proof: Euler-Lagrange v3}
We define the manifold $\mathcal{M}$ as the decision boundary for our black-box function. Assume that we can draw uniform samples from the decision boundary $m \sim p(m)$. We weight $p(m)$ according to the $L_2$ distance between each drawn sample $m$ and the fixed points $x, y$. Define the weighted distribution

\begin{comment}
\begin{equation}
    q(m | x, \rho) = \frac{1}{Z_m} \exp[-\rho || x-m||_2^2] p(m) \quad s.t. \; \;  Z_m = \int \exp[-\rho || x-m||_2^2] p(m) \; \textrm{d}m
\end{equation}
\end{comment}
\begin{equation}
    q(m | x, \rho) \propto  \exp[-\rho || x-m||_2^2] p(m)
\end{equation}

We evaluate the kernel function $k(x,y)$ by taking the expected value over the weighted distribution of samples from the manifold.

\begin{equation}
k(x,y) = \mathbb{E}_{m \sim q(m|x,\rho), m' \sim q(m'|y,\rho)}[\exp[-\lambda \textrm{geo}(m, m')]]
\end{equation}

\begin{equation}
= \int \int \exp[-\lambda \textrm{geo}(m, m')] \; q(m|x, \rho) \; q(m'|y, \rho) \; \textrm{d} m \textrm{d} m'
\end{equation}

In practice we can evaluate the integrals using Monte Carlo sampling, drawing $K$ samples from $p(m)$:
\begin{equation}
    k(x,y) \approx \frac{1}{Z_m Z_{m'} K^2}\sum_{i=1}^K \sum_{j=1}^K \exp[-\lambda \textrm{geo}(m_i, m_j)] \exp[-\rho (|| x-m_i||_2^2 + || y-m_j||_2^2)]
\end{equation}

$Z_m$, $Z_{m'}$ are the normalization constants for $q(m|x,\rho)$ and $q(m'|y,\rho)$, respectively. We can similarly estimate these values using $K$ Monte Carlo samples:

\begin{equation}
    Z_m = \int \exp[-\rho || x-m||_2^2] p(m) \; \textrm{d}m  \approx \frac{1}{K}\sum_{i=1}^K \exp[-\rho || x-m_i||_2^2]
\end{equation}

Putting it all together:
\begin{equation}
    \tilde k(x,y) = \frac{\sum_{i=1}^K \sum_{j=1}^K \exp[-\lambda \textrm{geo}(m_i, m_j)] \exp[-\rho (|| x-m_i||_2^2 + || y-m_j||_2^2)]}{(\sum_{i=1}^K \exp[-\rho || x-m_i||_2^2])(\sum_{j=1}^K \exp[-\rho || y-m_j||_2^2])}
\end{equation}
\begin{theorem}
Given a 1-dimensional manifold $\mathcal{M}$ and two points $x,y \in \mathbb{R}^2$. For any two points $a,b \in \mathcal{M}$, if we replace the manifold between $a$ and $b$ with some function $g$, the function $g$ that minimizes $k_{\mathcal{M}}(x,y)$ is a straight line between $a$ and $b$.
\end{theorem}
\begin{remark}
    This represents model complexity.
\end{remark}
\begin{corollary}
Any function other than the straight line will increase $k(x,y)$
\end{corollary}

\begin{proof}

Assume that $\mathcal{M}$ is a 1-dimensional manifold embedded in $\mathbb{R}^2$. Without loss of generality, assume that point $a = (0,0)$ and point $b = (0,c)$ for some constant $c$. Define a smooth, bijective function $\gamma:[0,1] \rightarrow \mathcal{M}$ such that $\gamma(0) = a$ and $\gamma(1) = b$. \textbf{We also define $\pi(t)$ (and similarly $\pi(t')$) as the truncacted distribution of $T = \gamma^{-1}(M)$, where $M \sim q(M|x,\rho)$.}\footnote{this needs to be fixed -- not sure how to define this} We reparametrize the kernel function:

\begin{equation}
k(x,y) = \int \int \exp[-\lambda \textrm{geo}(m, m')] \; q(m|x, \rho) \; q(m'|y, \rho) \; \textrm{d} m \textrm{d} m'
\end{equation}

\begin{equation}
=\int^{0}_{1}  \int^{0}_{1}   \exp[- \lambda \int_{t'}^{t} \sqrt{g(\dot \gamma(s), \dot \gamma(s))} \; \textrm{d}s ]\;  \pi(t) \; \pi(t')  \; \textrm{d} t \textrm{d} t'
\end{equation}

\begin{equation}
=\int^{0}_{1}  \pi(t') \int^{0}_{1}   \pi(t)  \; \exp[- \lambda \int_{t'}^{t} \sqrt{\sum_{i,j} g_{ij}\dot \gamma(s)_i \dot \gamma(s)_j} \; \textrm{d}s ]\;    \; \textrm{d} t \textrm{d} t'
\end{equation}

$\gamma^i$ indicates the $i^{th}$ component of vector $\gamma$. Since we assume that the decision boundary is an embedding in $\mathbb{R}^2$, the Riemannian Metric is Euclidean; i.e.

\[ g_{ij} =
\begin{cases}
    1 & i = j \\
    0 & i \neq j
 \end{cases}
\]

We rewrite the geodesic distance:

\begin{equation}
=\int^{0}_{1}  \pi(t') \int^{0}_{1}   \pi(t)  \; \exp[- \lambda \int_{t'}^{t} ||\dot \gamma (s)|| \; \textrm{d}s ]\;    \; \textrm{d} t \textrm{d} t'
\end{equation}

First, we prove that for any given $t, t'$, the inner integral $\int_{t'}^{t} ||\dot \gamma (s)|| \textrm{d}s$ is minimized by a straight line. This is trivial (using Euler-Lagrange).

Next, consider the middle integral $\int^{0}_{1}   \pi(t)  \exp[- \lambda \int_{t'_0}^{t} ||\dot \gamma (s)|| \; \textrm{d}s ] \textrm{d} t$ where $t'$ is fixed to be some value, denoted $t'_0$.

\end{proof}

\scratchline

We next take the expected value of the LHS and RHS of Eq. \ref{eq:inequality} with respect to the distributions $q(g_\epsilon(t))|x, \rho) q(g_\epsilon(t'))|y, \rho)$ and $q((t, g_0(t))|x, \rho) q((t', g_0(t'))|y, \rho)$, respectively. Note that these two distributions are different weightings of the same underlying distribution $p(m)$ (Eq. \ref{eq:q_weighting}).

\begin{equation}
    \lim_{\rho \rightarrow 0} q(g_\epsilon(t))|x, \rho) = p(g_\epsilon(t))) = p((t, g_0(t))) = \frac{1}{c} \mathbbm{1}_{[0,c]}(t)
\end{equation}

\begin{comment}
** But the problem is that $p((t, g_\epsilon(t))) \neq p((t, g_0(t)))$ so we have the same problem as before.. we need a way to go from $g_\epsilon(t)$ to $g_0(t)$.

Maybe the key is that we put a uniform distribution over $t$, then compare $p((t, g_\epsilon(t)))$ and $p((t, g_0(t)))$. There should be some mapping between $g_\epsilon(t)$ and $g_0(t)$. Because $\ddot g_\epsilon(t) \geq \ddot g_0(t)$, the scaling factor must be $\in [0,1]$? Therefore the density at any given interval will be lower? Since you're essentially stretching a uniform distribution over a wider support.
\end{comment}

\begin{comment}

\begin{align*}\numberthis
    \int_0^c \int_0^c \exp[-\lambda \textrm{geo}((t, g_\epsilon(t)), g_\epsilon(t')))] p((t, g_\epsilon(t))) p(g_\epsilon(t'))) \; \textrm{d} t \textrm{d} t' \\
    \leq \int_0^c \int_0^c \exp[-\lambda \textrm{geo}(g_0(t), g_0(t'))] p((t, g_0(t))) p((t', g_0(t'))) \; \textrm{d} t \textrm{d} t'
\end{align*}
\end{comment}

\begin{align*}\numberthis
    \frac{1}{c^2} \int_0^c \int_0^c \exp[-\lambda \textrm{geo}((t, g_\epsilon(t)), g_\epsilon(t')))] \; \textrm{d} t \textrm{d} t' \\
    \leq \frac{1}{c^2} \int_0^c \int_0^c \exp[-\lambda \textrm{geo}(g_0(t), g_0(t'))]  \; \textrm{d} t \textrm{d} t'
\end{align*}

\begin{align*}\numberthis
    k_\epsilon(x,y) \leq k_0(x,y)
\end{align*}

We therefore establish that the inequality is true in the limit when $\rho$ approaches zero. Since $k(x,y)$ is continuous with respect to $\rho$, it follows that there exists a $\rho>0$ such that the inequality still holds.

\section{Proof of Theorems and Multiclass Extension}

\subsection{Theorem \ref{thm:kernel_expgeodesic}: Relation to Exponential Geodesic Kernel} \label{proof:points_on_manifold}

\begin{align*}
    k(x,y) =  \int \int \exp[-\lambda d_{\textrm{geo}}(m, m')] q(m | x, \rho) q(m' | y, \rho)  \; \; \textrm{d} m' \textrm{d} m \\
    s.t. \quad q(m | x, \rho) \propto  \exp[-\rho || x-m||_2^2] p(m) \\
\end{align*}

Note that $\rho$ controls how to weight manifold samples close to $x, y$.
We take $\lim_{\rho \rightarrow \infty}$:

\[ \lim_{\rho \rightarrow \infty}  q(m | x, \rho)q(m' | y, \rho)=
\begin{cases}
    1 & x = m \; \textrm{and}\;  y = m' \\
    0 & \textrm{Otherwise}
 \end{cases}
\]

\iffalse
= \[ \begin{cases}
    \int_{m \in \mathcal{M}} \int_{m' \in \mathcal{M}} \exp[-\lambda d_{\textrm{geo}}(m, m')] \; \; \textrm{d} m' \textrm{d}m & x = m \; \textrm{and}\;  y = m' \\
    0 & \textrm{Otherwise}
 \end{cases}
\]
\fi

Therefore the function within the integral of $k(x,y)$ evaluates to zero at all points except $x = m$ and $y=m'$. Since $x,y \in \mathcal{M}$ we can evaluate the integral:
$$k(x,y) = \exp[-\lambda d_{\textrm{geo}}(x, y)] $$

\subsection{Theorem \ref{thm:kernel_complexity}: Kernel Similarity and Decision Boundary Complexity} \label{proof:kernel_complexity}

\begin{comment}
Consider $h(\rho) - h(c)$:

\begin{equation}
= \int \int \mathcal{A}[\frac{\exp[-\rho (|| x-m||_2^2 + || y-m'||_2^2)]}{Z_m(\rho) Z_{m'}(\rho)} - \frac{\exp[-c (|| x-m||_2^2 + || y-m'||_2^2)]}{Z_m(c) Z_{m'}(c)}]  \textrm{d} m \textrm{d} m'
\end{equation}

where $\mathcal{A} = \exp[-\lambda d_{\textrm{geo}}(m, m')] p(m, m')$.
\end{comment}

\begin{comment}
We are given a binary classifier $F:\mathbb{R}^d \rightarrow [0,1]$ and a classification threshold $c$, where $F(x)\geq c$ is class 1, and $F(x)<c$ is class 0. We define the manifold $\mathcal{M} = \{m:F(m) = c\}$.
\end{comment}

\begin{comment}
Based on the definition of perturbation, r is identity for all the charts except for chart i. Based on the definition of the WEG kernel, we only investigate the effect of perturbation
on a compact subset K_i of U_i.

Define a piecewise linear manifold $\mathcal{M}$ of dimension $(d-1)$ embedded in $\mathbb{R}^d$; i.e. $\mathcal{M}$ is a manifold such that for every point $x \in \mathcal{M}$ there exists an open set $U$ where $x \in U$ as well as a piecewise linear homeomorphism $f:U \rightarrow \mathbb{R}^{d-1}$.

f is piecewise linear. each one of them can be covered by a ball. you cover the whole space with the cover. it has a finite subcover. Each one corresponds to a linear region.

We fix a point $x_0 \in \mathcal{K_i}$ and take an open neighborhood $U_0$ of $x_0$ s.t. there exists $f:U_0 \rightarrow \mathbb{R}^{d-1}$. We construct a bounded set $U_1 \subset U_0$, and define the boundary of $U_1$ as $\partial U_1 = \bar U_1 \setminus \textrm{int}_{\mathbb{R}^{d-1}} (U_1)$, where int indicates interior and $\bar U_1$ indicates the closure of $U_1$. We also define $\mathcal{T} = [0,1]^{d-1}$.
\end{comment}

From definition \ref{def:1}, given any perturbation $\widetilde{\mathcal{P}}$ on $\mathcal{P}$, there must exist a compact subset $K_i \subset U_i$ s.t.  $\left.R\right|_{\mathcal{P} \setminus \textrm{int}(K_i)} = \left.\textrm{id}\right|_{\mathcal{P} \setminus \textrm{int} (K_i)}$ and $\left.R\right|_{\textrm{int}(K_i)} \neq \left.\textrm{id}\right|_{\textrm{int}(K_i)}$. Furthermore there exists a linear homeomorphism between an open subset $\widetilde{U_i} \subseteq U_i$ with $\mathbb{R}^{d-1}$ which contains $K_i$.

We parametrize $K_i$ using a smooth function $g: \mathcal{T} \rightarrow  K_i \:$ s.t. $g(t) \in \partial K_i \; \forall t \in \partial \mathcal{T}$.

We further define $g_\epsilon(t) = g(t) + \epsilon \eta (t)$, for some perturbation $\epsilon \in \mathbb{R}$ and a smooth function $\eta: \mathcal{T} \rightarrow  \mathbb{R}^{d-1}$
We also restrict $\eta$ such that $\eta(t) = \mathbf{0} \; \forall t \in \partial \mathcal{T}$ and $\exists \; t_0 \in \mathcal{T}$ s.t. $\eta(t_0) \neq g(t_0)$.
In other words, $\eta$ is a smooth function where $g_\epsilon(t) = g(t) \; \forall \epsilon>0, \forall t \in \partial \mathcal{T}$, but is not identical to $g$ for all $t \in \mathcal{T}$. Using $g_\epsilon(t)$, we define the manifold $\mathcal{P}_\epsilon = \{g_\epsilon(t):t \in \mathcal{T}\}$.

To complete the proof, we want to show that the kernel similarity between any two given points $x,y \in \mathbb{R}^d$ is lower when using the manifold $\mathcal{P}_\epsilon$ for $\epsilon>0$ as opposed to the manifold $\mathcal{P}_0$.
We therefore want to compare the two respective kernels $k_{\epsilon}(x,y)$ and $k_0(x,y)$.
Note that in this proof we consider the local effects of $\mathcal{P}$ on the kernel similarity through $\mathcal{P}_0$ and $\mathcal{P}_\epsilon$ exclusively, ignoring the manifold $\mathcal{P} \setminus U_0$
Using Euler-Lagrange, we can calculate a lower bound for $d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))$. In particular, for any $t,t' \in \mathcal{T}$, $d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t'))) \geq d_{\textrm{geo}}(g_0(t), g_0(t'))$.

\begin{equation}
    d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t')) \geq d_{\textrm{geo}}(g_0(t), g_0(t'))
\end{equation}

\begin{equation}  \label{eq:inequality}
    \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \leq  \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))]
\end{equation}

\begin{equation} \label{eq:inequality2}
    \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \;  \textrm{d}t \textrm{d}t' \leq \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] \;  \textrm{d}t \textrm{d}t'
\end{equation}

Note that in Eq. \ref{eq:inequality2} we are integrating over all possible values of $t,t'$, therefore the inequality is tight iff $g_\epsilon(t) = g_0(t)$ $\forall t \in \mathcal{T}$; i.e. $\epsilon=0$ (see proof in \ref{proof:integral}). The case of $\epsilon = 0$ is trivial; we instead assume $\epsilon > 0$, in which case we can establish the following strict inequality:

\begin{equation} \label{eq:inequality3}
    \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \;  \textrm{d}t \textrm{d}t' < \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] \;  \textrm{d}t \textrm{d}t'
\end{equation}

Define uniform random variables $T$, $T'$ over the domain of $g$, i.e. $T, T' \sim \mathcal{U}_{\mathcal{T}}$. Then we have:

\begin{equation} \label{eq:inequality_int_over_t}
    \mathbb{E}_{T,T' \sim \mathcal{U}_{[0,1]}}[ \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(T), g_\epsilon(T'))]] < \mathbb{E}_{T,T' \sim \mathcal{U}_{[0,1]}}[ \exp[-\lambda d_{\textrm{geo}}(g_0(T), g_0(T'))]]
\end{equation}

\begin{comment}
Under the transformations $g_\epsilon$ and $g_0$, we obtain distributions $p_\epsilon$ and $p_0$ over the modified manifolds $\mathcal{P}_\epsilon$ and $\mathcal{P}_0$, respectively.
\textbf{this needs to be explained better}
\end{comment}

\begin{equation} \label{eq:inequalityfinal}
    \mathbb{E}_{M, M' \sim p_\epsilon(M)}[ \exp[-\lambda d_{\textrm{geo}}(M, M')]] < \mathbb{E}_{M, M' \sim p_0(M)}[ \exp[-\lambda d_{\textrm{geo}}(M, M')]]
\end{equation}

\begin{comment}
\begin{align*} \numberthis \label{eq:inequalityfinal}
    \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))] \frac{||\dot g_\epsilon(t)|| \; ||\dot g_\epsilon(t)||}{\int \int ||\dot g_\epsilon(t)|| \;  ||\dot g_\epsilon(t')|| \textrm{d}t \textrm{d}t'} \; \textrm{d} t \textrm{d} t' \\
    \leq  \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))]  \frac{||\dot g_0(t)|| \; ||\dot g_0(t)||}{\int \int ||\dot g_0(t)|| \; || \dot g_0(t')||\textrm{d} t \textrm{d} t'} \; \textrm{d} t \textrm{d} t'
\end{align*}
\end{comment}

\begin{comment}
\begin{align*} \numberthis \label{eq:inequality2}
    \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))] p(g_\epsilon(t))) p(g_\epsilon(t'))) \; \textrm{d} t \textrm{d} t' \\
    <  \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] p((t, g_0(t))) p((t', g_0(t'))) \; \textrm{d} t \textrm{d} t'
\end{align*}
\end{comment}
We define the random variable $M = g_\epsilon(T)$ with distribution $p_\epsilon(M)$. The distribution $p_\epsilon(M)$ represents the uniform distribution $\mathcal{U}_{\mathcal{T}}$ mapped to the manifold $\mathcal{P}_\epsilon$ using $g_\epsilon (T)$. The step from Eq. \ref{eq:inequality_int_over_t} to Eq. \ref{eq:inequalityfinal} uses a property of distribution transformations (Eq. 2.2.5 in \cite{casella_berger}).

Next, compare either side of Eq. \ref{eq:inequalityfinal} to our kernel formulation shown below in Eq. \ref{eq:inequalitykernel}. The kernel $k_\epsilon(x,y | \rho, \lambda)$ takes an expected value over $q_\epsilon(M|x,\rho)$ and $q_\epsilon(M'|y,\rho)$, which are equivalent to $p_\epsilon(M)$ and $p_\epsilon(M')$ weighted with respect to $x$, $y$, and a hyperparameter $\rho \geq 0$.

\begin{align*} \numberthis \label{eq:inequalitykernel}
    k_\epsilon(x,y | \rho, \lambda) = \mathbb{E}_{M \sim q_\epsilon(M|x,\rho), M' \sim q_\epsilon(M'|y,\rho)}[\exp[-\lambda d_{\textrm{geo}}(M, M')]] \\
    s.t. \quad q_\epsilon(M | x, \rho) \propto  \exp[-\rho || x-M||_2^2] p_\epsilon(M)  \; \; \\
    s.t. \quad q_\epsilon(M' | y, \rho) \propto  \exp[-\rho || y-M'||_2^2] p_\epsilon(M') \\
\end{align*}

Note that when $\rho$ is set to zero, $q(M|x,0) = p(M)$ and $q(M'|y,0) = p(M')$. Therefore Eq. \ref{eq:inequalityfinal} is equivalent to the inequality $k_\epsilon(x,y|0, \lambda) < k_0(x,y | 0, \lambda)$.

We next want to prove that the inequality $k_\epsilon(x,y | \rho, \lambda) < k_0(x,y|\rho, \lambda)$ also holds for non-zero values of $\rho$. For convenience, define

\begin{equation}
    f(\rho) = k_0(x,y|\rho, \lambda) - k_\epsilon(x,y | \rho, \lambda)
\end{equation}

Under this definition, we want to prove there exists $\rho_0>0$ such that $f(\rho) > 0$ $\forall \rho < \rho_0$. From Eq. \ref{eq:inequalityfinal}, we established that $f(0)>0$. Assume that

\begin{equation}
\lim_{\rho \rightarrow 0} f(\rho) = c
\end{equation}

It therefore follows that $c>0$.
In addition, note that $f(\rho)$ is continuous with respect to $\rho$ (see proof in section \ref{proof:kernel_continuity}). Therefore for any $\epsilon > 0$ there exists $\delta > 0$ s.t. $\rho < \delta$ implies $|f(\rho) - c| < \epsilon$.

We choose $\epsilon = c$ and the define the corresponding $\delta$ to be $\rho_0$. Therefore:

\begin{equation}
    \rho < \rho_0 \Rightarrow |f(\rho) - c| < c
\end{equation}

\begin{equation}
    \rho < \rho_0 \Rightarrow  0 < f(\rho) < 2c
\end{equation}

Since this result holds for any $i$, it follows that the piecewise linear manifold $\mathcal{P}$ is a local minimum under any perturbation along a specific chart or combination of charts with respect to the kernel similarity $k(x,y) \; \forall x,y \in \mathbb{R}^d$.

\subsubsection{Proof: Eq. \ref{eq:inequality3}} \label{proof:integral}
We want to prove:

\begin{align*} \numberthis \label{eq:integral_proof1}
    \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \;  \textrm{d}t \textrm{d}t' = \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] \;  \textrm{d}t \textrm{d}t' \\
    \Rightarrow g_\epsilon(t) = g_0(t) \quad \forall t \in \mathcal{T}
\end{align*}

Consider the LHS of Eq. \ref{eq:integral_proof1}:

\begin{equation}
    \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \;  \textrm{d}t \textrm{d}t' = \int_\mathcal{T} \int_\mathcal{T} \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] \;  \textrm{d}t \textrm{d}t'
\end{equation}

\begin{equation} \label{eq:integral_proof2}
    \int_\mathcal{T} \int_\mathcal{T} \underbrace{\exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] -\exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))]}_{h(t,t')} \;  \textrm{d}t \textrm{d}t' = 0
\end{equation}

Define $h(t,t')$ as the function inside the integrals in Eq. \ref{eq:integral_proof2}. From Eq. \ref{eq:inequality}, $h(t,t') \geq 0$ $\forall t,t' \in \mathcal{T}$. Since $h$ is continuous (see proof in \ref{proof:h_continuity}) and $\int_\mathcal{T} \int_\mathcal{T} h(t,t') \textrm{d}t \textrm{d}t' = 0$, it follows that $h(t,t') = 0$ $\forall t,t' \in \mathcal{T}$ (Ch.6 \cite{baby_rudin}).

It therefore follows that:

\begin{equation}\label{eq:integral_proof3}
    \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] = \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] \quad \forall t,t' \in \mathcal{T}
\end{equation}

From the definition of $\eta(t)$ in $g_\epsilon(t) = g(t) + \epsilon \eta(t)$, there must exist $t \in \mathcal{T}$ s.t. $\eta(t) \neq 0$. Therefore $\epsilon$ must be zero for Eq. \ref{eq:integral_proof3} to hold. It follows that $g_\epsilon(t) = g_0(t)$ $\forall t \in \mathcal{T}$.

\subsubsection{Proof: Continuity of \texorpdfstring{$h(t,t')$}{h(t,t')}} \label{proof:h_continuity}

We prove that $h(t,t')$ is continuous with respect to $t,t'$. First note that by definition, $g_\epsilon(t)$ is a continuous parametrization of the manifold $\mathcal{P}_\epsilon$. From \cite{burago2022course}, it follows that for any two points $g_\epsilon(t), g_\epsilon(t') \in \mathcal{P}_\epsilon$,  $d_{\textrm{geo}}(g_\epsilon(t),g_\epsilon(t'))$ is continuous. Since the exponential functional preserves continuity and the sum of continuous functions are also continuous, it follows that $h(t,t')$ is continuous.

\subsubsection{Proof: Continuity of \texorpdfstring{$k(x,y)$}{k(x,y)} With Respect To \texorpdfstring{$\rho$}{Rho}} \label{proof:kernel_continuity}

We prove that $k(x,y)$ is continuous with respect to $\rho$.
\begin{equation}
k(x,y) = \int \int \exp[-\lambda d_{\textrm{geo}}(m, m')] \; q(m|x, \rho) \; q(m'|y, \rho) \; \textrm{d} m \textrm{d} m'
\end{equation}

\begin{align*} \numberthis \label{eq:continuousproof}
    =\frac{1}{Z_m(\rho) Z_{m'}(\rho)}  \int \int \mathcal{A} \underbrace{\exp[-\rho (|| x-m||_2^2 + || y-m'||_2^2)]}_{Z(\rho)} \; \textrm{d} m \textrm{d} m' \\
    s.t. \quad Z_m(\rho) = \int \exp[-\rho || x-m||_2^2] p(m) \; \textrm{d}m \; \; \; \;  \\
    \quad Z_{m'}(\rho) = \int \exp[-\rho || y-m'||_2^2] p(m') \; \textrm{d}m' \\
    \mathcal{A} = \exp[-\lambda d_{\textrm{geo}}(m, m')] p(m)p(m') \; \; \; \; \; \; \; \quad
\end{align*}

Define $h(\rho) = \rho \mathcal{B}$, where $\mathcal{B}$ is a constant. Consider $h(\rho) - h(\rho_0)$, where $\rho_0$ is a fixed positive constant:

\begin{equation}
    |h(\rho) - h(\rho_0)| = |\rho \mathcal{B} - \rho_0 \mathcal{B}|
\end{equation}

\begin{equation}
    = |(\rho - \rho_0)\mathcal{B}| < \delta |\mathcal{B}|
\end{equation}

It follows that $\forall \; \epsilon>0$, $\exists \; \delta = \frac{\epsilon}{|\mathcal{B}|} > 0$ such that $|\rho - \rho_0| < \delta \; \Rightarrow \; |h(\rho) - h(\rho_0)| < \epsilon$. Therefore $h$ is continuous for all $\rho \in \mathbb{R}^+$.

We set $\mathcal{B}$ to be $|| x-m||_2^2$, $|| y-m'||_2^2$, and $|| x-m||_2^2 + || y-m'||_2^2$, which shows that $Z_m (\rho)$, $Z_{m'}(\rho)$, and $Z(\rho)$ are also continuous, respectively. It then follows that the entirety of Eq. \ref{eq:continuousproof} is continuous.

\begin{comment}
\scratchline
Assume that $\mathcal{P}$ is a 1-dimensional manifold embedded in $\mathbb{R}^2$. We choose two points $a, b \in \mathcal{P}$. We parametrize the manifold using a continuous function $g: [0,1] \rightarrow  \mathbb{R}^2$, where $g(0) = a$, $g(1) = b$, and $\ddot g(t) = 0 \; \forall t \in [0,1]$.
We further define $g_\epsilon(t) = g(t) + \epsilon \eta (t)$, for some $\epsilon \in \mathbb{R}$ and where $\eta: [0,1] \rightarrow  \mathbb{R}^2$ is some function such that $\eta(0) = 0$ and $\eta(1) = 0$. We also restrict $\eta$ to be a continuous function of $t$ such that $\exists \; t_0 \in [0,1]$ s.t. $\eta(t_0) \neq g(t_0)$, i.e. $\eta$ is not identical to $g$.
Using this parametrization, define the modified manifold $\mathcal{M_\epsilon}$ as the set $\{g_\epsilon(t): t \in [0,1]\}$
\hl{Becreful, for instance you are already using it here with a lot of underneath assumption why this is a manifold?}
.
\end{comment}

\begin{comment}
We assume there exists some distribution $p(t)$ over the support $t \in [0,c]$ from which samples can be drawn from; without loss of generality assume $p(t)$ is uniform. We can then define $p(m)$ as follows:

\begin{equation} \label{eq:pdf_m}
    p(m) = p((t, g_\epsilon(t))) = p(t) = \frac{1}{c} \mathbbm{1}_{[0,c]}(t)
\end{equation}
\end{comment}

\begin{comment}
We rewrite the kernel function using the defined parametrization. For convenience, denote $k_\epsilon(x,y)$ as the kernel $k(x,y)$ defined under the manifold $\mathcal{M_\epsilon}$.
\begin{equation}
k_\mathcal{\epsilon}(x,y) = \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t), g_\epsilon(t'))] q(g_\epsilon(t)|x, \rho) q(g_\epsilon(t')|y, \rho) \; \textrm{d} t \textrm{d} t'
\end{equation}
\end{comment}

\begin{comment}
\begin{align*}
    k_{\epsilon}(x,y) - k_{0}(x,y) = \\
    \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))] q(g_\epsilon(t))|x, \rho) q(g_\epsilon(t'))|y, \rho) \; \textrm{d} t \textrm{d} t' \\
    - \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] q((t, g_0(t))|x, \rho) q((t', g_0(t'))|y, \rho) \; \textrm{d} t \textrm{d} t'
\end{align*}

Combining the two integrals:
\begin{align*}
    = \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))] q(g_\epsilon(t))|x, \rho) q(g_\epsilon(t'))|y, \rho)  \\
    -  \exp[-\lambda d_{\textrm{geo}}(g_0(t), g_0(t'))] q((t, g_0(t))|x, \rho) q((t', g_0(t'))|y, \rho) \; \textrm{d} t \textrm{d} t'
\end{align*}
\end{comment}

\begin{comment}
\textbf{NOTE: THIS NEEDS TO BE GENERALIZED TO $\mathbb{R}^d$}. Specifically, we need the $t_0$ to be within the endpoints.

To create a contradiction, assume that $g_\epsilon(t) \neq g_0(t)$ $\forall t \in [0,1]$. Therefore there must exist $t_0 \in (0,1)$ such that $g_\epsilon(t_0) \neq g_0(t_0)$. Note that from our definition  of $g_\epsilon$ we have $g_\epsilon(t) = g_0(t)$ $\forall t \in \{0,1\}$. Consider the geodesic distance between two boundary points $g_\epsilon(0) = g_0(0)$ and $g_\epsilon(1)=g_0(1)$. Per Euler-Lagrange, the minimum geodesic distance between the two points is

\begin{equation}
    d_{\textrm{geo}}(g_0(0), g_0(1))
\end{equation}

Since $g_\epsilon(t_0) \neq g_0(t_0)$ and $g_\epsilon$ is continuous on (0,1), it follows that

\begin{equation}
    d_{\textrm{geo}}(g_\epsilon(0), g_\epsilon(1)) > d_{\textrm{geo}}(g_0(0), g_0(1))
\end{equation}

\begin{equation}
    \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(0), g_\epsilon(1))] <  \exp[-\lambda d_{\textrm{geo}}(g_0(0), g_0(1))]
\end{equation}

\begin{equation}
    h(0,1) > 0
\end{equation}

This causes a contradiction and concludes the proof.
\end{comment}

\begin{comment}
\subsubsection{Proof: \texorpdfstring{$h(t,t') = 0 \; \forall t,t' \in [0,1]$}{h(t,t')=0}} \label{proof:h_equal_zero}
We prove that if $h(t,t') \geq 0$ and $h$ is continuous then this implies $h(t,t') = 0$ $\forall t,t' \in [0,1]$.

Proof by contradiction. Assume that $h(t,t') \neq 0$ $\forall t,t' \in [0,1]$, therefore there must exist $t_0,t_0'$ such that $h(t_0,t_0') > 0$.
From continuity, there must exist $\delta>0$ such that $\forall |t - t_0| < \delta$ and $\forall |t' - t_0'| < \delta$
\end{comment}

\begin{comment}
Note that $f(\rho)$ is continuous with respect to $\rho$ (see proof in Sec. \ref{proof:kernel_continuity}). We choose a $\rho_1$ such that $f(\rho_1) \leq 0$ (otherwise the inequality holds for all $\rho$, which trivially satisfies the proof). Since $f(0) > 0$ and $f(\rho_1) \leq 0$, there must exist $\rho_0 \in (0, \rho_1]$ s.t. $f(\rho_0) = 0$ (Intermediate Value Theorem, 4.23 in \cite{baby_rudin}). Therefore for any $\rho < \rho_0$ and any $x,y \in \mathbb{R}^2$, it follows that $f(\rho) > 0$.

Consider the set $\Omega = \{\rho: \rho \in \mathbb{R}^+, f(\rho) \leq 0 \}$. We take $\rho_0 = \inf \Omega$. Since $f(0)>0$, it follows that $0 \notin \Omega$ and $\rho_0 >0$. Therefore $f(\rho) > 0$ $ \forall \rho \in (0,\rho_0)$

 Therefore there exists $\rho_0$ such that $\rho < \rho_0$ $\Rightarrow$ $|f(\rho_0) - f(0)| < \epsilon\;$ $\forall \epsilon>0$. We choose $\rho_0$ and $\rho_\epsilon$ such that:

\begin{equation}
    |f_\epsilon(\rho_\epsilon) - f_\epsilon(0)| \leq \frac{\delta}{2} \quad |f_0(\rho_0) - f_0(0)| \leq \frac{\delta}{2}
\end{equation}
\end{comment}

\begin{comment}
More formally, $\forall \; \epsilon>0$, $\exists \; \delta > 0$ such that $|\rho - a| < \delta \; \Rightarrow \; |h(\rho) - h(a)| < \epsilon$. We take $\epsilon = \frac{1}{n}$ $n \in \mathbb{N}$. Therefore there must exist a $\delta_n$ s.t. $q_n < \delta_n$ $\Rightarrow$ $|f(q_n) - f(0)| < \frac{1}{n}\end{comment}

\begin{comment}
Note that the inequality in Eq. \ref{eq:inequality} holds for any $t, t' \in [0,c]$. Therefore it follows:

\begin{equation}
    \sup_{t, t' \in [0,c]} -\lambda \exp[d_{\textrm{geo}}((t, g_\epsilon(t)), g_\epsilon(t')))] \leq \inf_{t, t' \in [0,c]} -\lambda \exp[d_{\textrm{geo}}(g_0(t), g_0(t'))]
\end{equation}

\begin{align*}\numberthis
    \mathbb{E}[-\lambda \exp[d_{\textrm{geo}}((t, g_\epsilon(t)), g_\epsilon(t')))]]\leq \sup_{t, t' \in [0,c]} -\lambda \exp[d_{\textrm{geo}}((t, g_\epsilon(t)), g_\epsilon(t')))] \\
    \leq \inf_{t, t' \in [0,c]} -\lambda \exp[d_{\textrm{geo}}(g_0(t), g_0(t'))] \leq \mathbb{E}[-\lambda \exp[d_{\textrm{geo}}(g_0(t), g_0(t'))]]
\end{align*}

\begin{align*}\numberthis
    k_\epsilon(x,y) \leq k_0(x,y)
\end{align*}
\end{comment}

\begin{comment}
 Consider $\lim_{\rho \rightarrow \infty} f (\rho)$:

\begin{align*} \numberthis
    \lim_{\rho \rightarrow \infty} f (\rho) = \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_0(t)), g_0(t')))] \delta(x - g_0(t)))  \delta(y - g_0(t'))) \; \textrm{d} t \textrm{d} t' \\
    - \int_0^c \int_0^c \exp[-\lambda d_{\textrm{geo}}(g_\epsilon(t)), g_\epsilon(t')))] \delta(x - g_\epsilon(t)))  \delta(y - g_\epsilon(t'))) \; \textrm{d} t \textrm{d} t'
\end{align*}

where $\delta$ is the dirac delta function. It is clear that the inequality $\lim_{\rho \rightarrow \infty} f(\rho) >0 $ does not necessarily hold for certain values of $x$ and $y$, specifically when $x,y \in \mathcal{P}_\epsilon$ and $x,y \notin \mathcal{P}_0$. However, we will next show that that there exists a $\rho_0>0$ such that $f(\rho) >0$ $\forall \rho < \rho_0$ $\forall x,y \in \mathbb{R}^2$.
\end{comment}

\subsection{Extending to Multiclass Classifiers} \label{sec:multiclass}
In the multiclass case we define a black-box prediction model $F:\mathcal{X} \to \mathbb{R}^c$. We consider the one-vs-all DB for every class $y \in \mathcal{Y} = \{1, \ldots, c\}$, defined as $\mathcal{M}_y = \{x \in \mathbb{R}^{d} : F_y(x) = \max_{i \in \mathcal{Y}} \: F_i(x) = \max_{j \neq y \in \mathcal{Y}} F_j(x) \}$, where $F_k$ indicates the model output for class $k$. We then apply the GPEC framework separately to each class using the respective DB. The uncertainty estimate of the GP model would be of dimension $d \times s$.

\section{Implementation Details} \label{app:implementation_details}

\subsection{Adversarial Sample Filtering Multi-class Models}
\label{app:adv}
We elect to sample from multi-class neural network decision boundaries by running binary search on train-point adversarial example pairs. Specifically, given a test-point $x_0 \in \mathbb{R}^d$ and model prediction $y = \text{argmax}_{k \in \mathcal{Y}} F(x_0)$, decision boundary points may be generated by the following procedure:

First, for each class $v \in \mathcal{Y}$ a set of $M_v$ points is randomly sampled from the set of train points on which the model predicts class $v$:
\begin{equation}
 \mathcal{X}_v \subseteq
\{ x : \text{argmax}_{k \in \mathcal{Y}} \ F(x) = v, \ x \in \mathcal{X}_{tr} \}, \ |\mathcal{X}_v| = M_v
\end{equation}
$\forall v \in \mathcal{Y}$. An untargeted adversarial attack using some $l_p$ norm and radius $\epsilon$ is run on each point in $\mathcal{X}_y$, the set of points with the same class prediction as $x_0$. Each attack output $Attack_{un}(x, \epsilon) \in \mathbb{R}^d$ is paired with its corresponding input, resulting in the set
\begin{equation}
\mathcal{X}_{y'} = \{(x, Attack_{Un}(x, \epsilon)): x \in \mathcal{X}_y \},
\end{equation}
where for an element $(a,b) \in \mathcal{X}_{y'}$ we have $ \text{argmax}_{k \in \mathcal{Y}} \ F(a) = y, \text{argmax}_{k \in \mathcal{Y}} \ F(b) = v \neq y$, where $v$ is an unspecified class.

Likewise a targeted adversarial attack, with target class $y$, is run on each point in each of the sets of points that are not predicted as class $y$. Each  attack output $Attack_y(x, \epsilon) \in \mathbb{R}^d$ may be paired with its input $x$ resulting in sets
\begin{equation}
\mathcal{X}_{v'} = \{ (x, Attack_y(x, \epsilon)) : x \in \mathcal{X}_v \}
\end{equation}
$\forall v \neq y \in \mathcal{Y}$. Here, for an element $(a,b) \in \mathcal{X}_{v'}$ we have $\text{argmax}_{k \in \mathcal{Y}} \ F(a) = v, \text{argmax}_{k \in \mathcal{Y}} \ F(b) = y$.

Thus, we have generated a diverse set of $\sum_{v \in \mathcal{Y}}M_v$ pairs of points that lie on opposite sides of the decision boundary for class $y$. The segment between any pair from a given set $\mathcal{X}_{v'} \ v \neq y$ will necessarily contain a point on the class $v$ v.s. class $y$ decision boundary. Likewise, in the interest of further diversity, segments between any pair from the set $\mathcal{X}_{y'}$ will contain a point on the class $v$ v.s. class $y$ decision boundary, where $v \neq y \in \mathcal{Y}$ is unspecified. A binary search may be run on each pair to find the boundary point in the middle.

In practice the entire procedure may be amortized for each class, and ran for all classes as a single post-processing step immediately after training. This results in a dictionary of boundary points which may be efficiently queried on demand via the model predicted class of any given test point.

Each adversarial attack is attempted multiple times, once using each radius value $\epsilon$ in the list: $[0.0, 2e^{-4}, 5e^{-4}, 8e^{-4}, 1e^{-3}, 1e^{-3},1.5e^{-3},2e^{-3},3e^{-3},1e^{-2},1e^{-1},3e^{-1},5e^{-1},1.0]$. For a given input, the output of the successful attack with smallest $\epsilon$ is used. If no attack is successful at any radius, the input is discarded from further consideration.

In this work the Foolbox \cite{rauber2017foolbox, rauber2017foolboxnative} implementation  of the Projected Gradient Descent (PGD) \cite{madry2018towards} attack with the $l_\infty$ norm was used for both targeted and untargeted attacks.  The  $M_c$ values used for the relevant datasets are indicated below in  Appendix \ref{app:dataset_details}.

\section{Societal Impacts}
As machine learning models are increasingly relied upon in a diverse set of high-impact domains ranging from health-care to financial lending \citep{esteva2019guide, kose2021segmentation, doshi2017towards, 9155614, ashutosh}, it is crucial that users of these models can accurately interpret why predictions are made.
An understanding of why a model is making a certain prediction is important for users to trust it -- for instance a doctor may wish to know if a skin-cancer classifier's high test-set accuracy comes from the leveraging of truly diagnostic features, or a specific imaging device artifact.
However, further spurred by the advent of deep learning's increasing popularity \citep{krizhevsky2017imagenet}, many of the models deployed in these high-stakes fields are complex black box's; producing predictions which are non-trivial to explain the reasoning behind.

The development of many methods for explaining black-box predictions has arisen from this situation \citep{ribeiroLIME, shap_lundberg, covertFeatureRemovalUnifying2020, masoomi2022explanations}, but explanations may have varying quality and consistency. Before utilizing explanations in practice, it is essential that users know when, and when not, to trust them. Explanation uncertainty is one proxy for this notion of trust, in which more uncertain explanations may be deemed less trustworthy. In this work, we explore a new way to model explanation uncertainty, in terms of local decision-boundary complexity.
In tandem with the careful consideration of domain experts, our methodology may be used to assist in determining when explanations are reliable. Our theoretical results provide new insights towards what explanation uncertainty entails, and open the door for future methods expounding upon our formulation.

\section{Gaussian Process Background}
\label{app:GPReview}
A single-output Gaussian Process represents a distribution over \textit{functions} $f : \mathcal{X} \rightarrow \mathbb{R}$
\begin{equation}
f(x) \sim \mathcal{GP}(m(x), k(x,x')).
\end{equation}
Here $m : \mathcal{X} \rightarrow \mathbb{R} $ and $k : (\mathcal{X}, \mathcal{X}) \rightarrow \mathbb{R}$ are the mean and kernel (or covariance) functions respectively, which are chosen \textit{a priori} to encode the users assumptions about the data. The kernel function $k(x,x')$ reflects a notion of similarity between data points for which predictive distributions over $f(x),f(x')$ respect.
The prior $m(x)$ -- frequently considered to be less important -- is commonly chosen to be the constant $m(x) = 0$.

Specifically, a GP is an infinite collection of R.V's $f(x)$, each indexed by an element $x \in \mathcal{X}$. Importantly, any finite sub-collection of these R.V's
\begin{equation}
f(X_{tr}) = (f(x_1) \ldots, f(x_n)) \in \mathbb{R}^d,
\end{equation}
corresponding to some index set $X_{tr} = \{x_i\}_{i=1}^n \subset \mathcal{X}$,  follows the multivariate normal (MVN) distribution, i.e.
\begin{equation}
    f(X_{tr}) \sim \mathcal{N}(m(X_{tr}),K(X_{tr}, X_{tr})).
\end{equation}
Here the mean vector $m(X_{tr}) = (m(x_1), \ldots, m(x_n)) \in \mathbb{R}^n$ represents the mean function applied on each $x \in X_{tr}$ and the covariance matrix $K \in \mathbb{R}^{n\times n}$, also known as the gram matrix, contains each pairwise kernel-based similarity value $K_{ij} = k(x_i, x_j)$. Kernel function outputs correspond to dot products in potentially infinite dimensional expanded feature space, which allows for the encoding of nuanced notions of similarity; e.g. the exponential geodesic kernel referenced in this work \citep{feragen_geodesic_exponential_kernel}.

Making predictions with a GP is analogous to simply conditioning this normal distribution on our data. Considering a set of input,noise-free label pairs
\begin{equation}
\mathcal{D} = \{ (x_i, f(x_i)\}_{i=1}^n
\end{equation}
we may update our posterior over \textit{any subset} of the R.V's $f(x)$ by considering the joint normal over the subset and $\mathcal{D}$ and conditioning on $\mathcal{D}$. For instance, when choosing a singleton index set $\{x_0\}$, the posterior over $f(x_0) | D$ is another normal distribution which may be written
as\footnote{assuming prior $m(x) =0$}
\begin{equation}
    f(x_0) \sim \mathcal{N}(\bar{f}(x_0), \mathbb{V}[f(x_0)])
\end{equation}
where
\begin{align}
    \bar{f}(x_0) &= K(x_0,X_{tr})^T K(X_{tr},X_{tr})^{-1} f(X_{tr})\\
    \mathbb{V}[f(x_0)] &= k(x_0,x_0) - K(x_0,X_{tr})^TK(X_{tr},X_{tr})^{-1}K(x_0,X_{tr}) \label{noise-free-var}
\end{align}
and $K(x_0,X_{tr}) \in \mathbb{R}^{d}$ is defined element-wise by $K(x_0,X_{tr})_i = k(x_0, x_i)$.

Now we may consider the situation where our labels are noisy:
\begin{equation}
\mathcal{D} = \{(x_i,y_i)\}_{i=1}^n , \  y_i = f(x_i) + \epsilon, \ \epsilon \sim \mathcal{N}(0,\sigma^2), \ \sigma^2 \in \mathbb{R}_+.
\end{equation}
Here $y_i$ is equal to the quantity we wish to model, $f(x_i)$, with the addition of noise variable $\epsilon$. The conditional is still a MVN, but the mean and variance equations are slightly modified
\begin{align}
    \bar{f}(x_0) &= K(x_0,X_{tr})^T (K(X_{tr},X_{tr}) + \sigma^2I)^{-1} Y\\
    \mathbb{V}[f(x_0)] &= k(x_0,x_0) - K(x_0,X_{tr})^T(K(X_{tr},X_{tr}) + \sigma^2I)^{-1}K(x_0,X_{tr}) \label{noise-var},
\end{align}
where $Y \in \mathbb{R}^{n}$ has elements $Y_i = y_i$.

Notice how the variance $\sigma^2 I$ is added to $K(X_{tr},X_{tr})$ in the quadratic form in Eq. \ref{noise-var}, resulting in smaller eigenvalues after matrix inversion. Since this quadratic form is subtracted, the decision to model labels as noisy increases the uncertainty (variance) of estimates the GP posterior provides. This agrees with the intuition that noisy labels should result in more uncertain predictions.

While GPs may also be defined over vector valued functions, in this work the independence of each output component is assumed, allowing for modeling with $c\geq 1$ independent GPs. For more details see Ch.2 of \cite{GPML}, from which the notation and content of this section were inspired.

\section{Experiment Setup}

\subsection{Datasets and Models} \label{app:dataset_details}

\textbf{Census.} The UCI Census dataset consists of 32,561 samples from the 1994 census dataset. Each sample is a single person's response to the census questionaire. An XGBoost model is trained using the 12 features to predict whether the individual has income $\geq$ \$50k.

 \textbf{Online Shopper.} The UCI Online Shoppers dataset consists of clickstream data from 12,330 web sessions. Each session is generated from a different individual and specifies whether a revenue-generating transaction takes place. There are 17 other features including device information, types of pages accessed during the session, and date information. An XGBoost model is trained to predict whether a purchase occurs.

\textbf{German Credit.} The German Credit dataset consists of 1,000 samples; each sample represents an individual who takes credit from a bank. The classification task is to predict whether an individual is considered a good or bad risk. Features include demographic information, credit history, and information about existing loans. Categorical features are converted using a one-hot encoding, resulting in 24 total features.

\textbf{MNIST.} The MNIST dataset \citep{lecun-mnisthandwrittendigit-2010} consists of 70k grayscale images of dimension 28x28. Each image has a single handwritten numeral, from 0-9. A fully connected network with layer sizes 784-700-400-200-100-10 and ReLU activation functions was trained and validated on on $50,000$ and $10,000$ image label pairs, respectively. Training lasted for $30$ epochs with initial learning rate of $2$ and a learning rate decay of $\gamma = 0.5$ when training loss is plateaued. During adversarial example generation we used $M_y = 500$ and $M_c = 50 \ \forall c \neq y$.

\begin{comment}
\textbf{CIFAR10} An X-Y-Z layer fully connected network, with ReLU activation functions after all but the last layer was trained and validated on on $40,000$ and $10,000$ image label pairs respectively. \hl{Training lasted ... }

During adversarial example generation we used $M_y = ...$ and $M_c = ... \forall c \neq y$.

\end{comment}

\textbf{Fashion MNIST} The Fashion MNIST dataset \citep{fmnist} contains 70k grayscale images of dimension 28x28. There are 10 classes, each indicating a different article of clothing. We train a MLP model with the same architecture used for the MNIST dataset, however we increase training to $100$ epochs and increase the initial learning rate to $3$. During adversarial example generation we used $M_y = 500$ and $M_c = 50 \ \forall c \neq y$.

\subsection{Competitor Implementation Details} \label{app:competitor_details}

\textbf{BayesLIME and BayesSHAP.} \cite{reliable:neurips21} extend the methods LIME and KernelSHAP to use a Bayesian Framework. BayesLIME and BayesSHAP are fit using Bayesian linear regression models on perturbed outputs of the black-box model. The posterior distribution of the model weights are taken as the feature attributions instead of the frequentist estimate that characterizes LIME and KernelSHAP. We take the expected value of the posterior distribution as the point estimate for feature attributions, and the 95\% credible interval as the estimate of uncertainty. To implement BayesLIME and BayesSHAP we use the public implementation\footnote{https://github.com/dylan-slack/Modeling-Uncertainty-Local-Explainability}. We set the number of samples to 200, disable discretization for continuous variables, and calculate the explanations over all features. Otherwise, we use the default parameters for the implementation.

\textbf{CXPlain.} \cite{cxplain} introduces the explanation method CXPlain, which trains a surrogate explanation model based on a causal loss function. After training the surrogate model, the authors propose using a bootstrap resampling technique to estimate the variance of the predictions. In our experiments we implement the publicly available code\footnote{https://github.com/d909b/cxplain}. We use the default parameters, which include using a 2-layer UNet model \cite{unet} for the image datasets and a 2-layer MLP model for the tabular datasets. We take a 95\% confidence interval from the bootstrapped results as the estimate of uncertainty.
\section{Additional Results} \label{app:additional_results}

\subsection{Execution Time Comparison} \label{app:execution_time}
\input{tables/table_time.tex}

In Table \ref{table:time} we include an execution time comparison between the methods implemented in this paper. Results are averaged over 100 test samples. For MNIST and f-MNIST datasets, results evaluate the time to calculate uncertainty estimates with respect to all classes. All experiments were run on an internal cluster using AMD EPYC 7302 16-Core processors.
We observe from the results that the methods that amortization methods (GPEC-WEG, GPEC-RBF, CXPlain) are significantly faster than perturbation methods BayesLIME and BayesSHAP.

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1\linewidth]{figures/figure_fullfat_appendix.pdf}
   \end{center}

   \caption{Comparison of the change in quantified uncertainty of explanations as we change the number of samples for BayesSHAP and SSV. Row (A) visualizes the combined uncertainty estimate using GPEC and either BayesSHAP or SSV, using 200 samples for approximating the BayesSHAP / SSV explanation. In Row (B) we decrease the number of samples to 5 and recalculate the estimated uncertainty. Row ($\Delta$) represents the change in uncertainty estimate between (A) and (B). We see that the average uncertainty changes as we decrease the number of samples, which indicates that GPEC is able to capture the uncertainty arising from BayesSHAP / SSV approximation.}
   \label{fig:fullfat_appendix}
\end{figure}
\subsection{Evaluating Explainer Uncertainty through GPEC} \label{app:addnl_results_fullfat}

In section \ref{sec:experiment_fullfat} we evaluate GPEC's ability to combine uncertainty from the black-box decision boundary and the uncertainty estimate from BayesSHAP and SSV explainers. In Figure \ref{fig:fullfat_appendix} we extend this experiment to evaluate how well GPEC can capture the explainer uncertainty. We calculate the combined GPEC+explainer estimate using different numbers of approximation samples.

Both BayesSHAP and SSV depend on sampling to generate their explanations; having fewer samples increases the variance of their estimates. As we decrease the number of samples from 200 (Row A) to 5 (Row B) we would expect that the explainer uncertainty, and consequently the combined GPEC uncertainty, would increase. We see in Row $\Delta$ that the results follow our intuition; uncertainty increases for most of the plotted test points and uncertainty does not decrease for any points.

\subsection{Additional Results for Uncertainty Visualization Experiment} \label{app:addnl_results_uncertainty_visualization}

\textbf{Results for Y-Axis Feature} In Figure \ref{fig:uncertaintyfigure_feat1} we visualize the estimated explanation uncertainty as a heatmap for a grid of explanations. The generated plots only visualize the uncertainty for the feature on the x-axis. Due to space constraints, we list the results for the y-axis feature in the appendix, in Figure \ref{fig:uncertaintyfigure_feat1}. We can see that the results are in line with those from the x-axis figure.

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1.0\linewidth]{figures/UncertaintyFigure_appendix.pdf}
   \end{center}

   \caption{Complement to Figure \ref{fig:uncertaintyfigure}. Visualization of estimated explanation uncertainty where the heatmap represents level of uncertainty for the feature on the \textbf{y-axis}.}
   \label{fig:uncertaintyfigure_feat1}
\end{figure}

\textbf{Black-box model output} For reference, in Figure \ref{fig:uncertaintyfigure_bb_output} we plot the probability output of the XGBoost models used in the visualization experiment (Figure \ref{fig:uncertaintyfigure}).

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1.0\linewidth]{figures/UncertaintyFigure_bb_output.pdf}
   \end{center}

   \caption{Output of binary classifiers used in the experiments. The color indicates the predicted probability of each point. Points with probability $\geq$ 0.5 are classified as class 1, and points with probability $<0.5$ are classified as class 0. The decision boundary, which is the set of points $\{m: m \in \mathbb{R}^2, f(m) = 0.5 \}$ where $f$ is the classifier, is represented by the black line.}
   \label{fig:uncertaintyfigure_bb_output}
\end{figure}

\subsection{Sensitivity Analysis of WEG Kernel Parameters} \label{app:sensitivity_analysis}

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1.0\linewidth]{figures/sensitivity_1.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-10pt}
   \caption{Hyperparameter sensitivity analysis for the German Credit Dataset. Heatmap of estimated uncertainty for the x-axis variable under different $\rho$ and $\lambda$ parameter choices.}
   \label{fig:sensitivity_germancredit}
\end{figure}

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1.0\linewidth]{figures/sensitivity_2.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-10pt}
   \caption{Hyperparameter sensitivity analysis for the Census Dataset. Heatmap of estimated uncertainty for the x-axis variable under different $\rho$ and $\lambda$ parameter choices.}
   \label{fig:sensitivity_census1}
\end{figure}

\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1.0\linewidth]{figures/sensitivity_3.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-10pt}
   \caption{Hyperparameter sensitivity analysis for the Census Dataset. Heatmap of estimated uncertainty for the x-axis variable under different $\rho$ and $\lambda$ parameter choices.}
   \label{fig:sensitivity_census2}
\end{figure}

The WEG kernel formulation uses two parameters, $\rho$ and $\lambda$. The parameter $\rho$ controls the weighting between each datapoint and the manifold samples. As $\rho$ increases, the WEG kernel places more weight on manifold samples close in $\ell_2$ distance to the given datapoint. The parameter $\lambda$ acts as a bandwidth parameter for the exponential geodesic kernel. Increasing $\lambda$ increases the effect of the geodesic distance along the manifold. Therefore decision boundaries with higher complexity will have an increased effect on the WEG kernel similarity. In Figures \ref{fig:sensitivity_germancredit}, \ref{fig:sensitivity_census1}, and \ref{fig:sensitivity_census2} we plot heatmaps for various combinations of $\rho$ and $\lambda$ parameters to evaluate the change in the uncertainty estimate. The black line is the decision boundary and the red points are the samples used for training GPEC. Please note that the heatmap scales are not necessarily the same for each plot.

\section{Limitations and Conclusion}
\vspace{-2mm}
Generating uncertainty estimates for feature attribution explanations are essential for building reliable explanations.
We introduce a novel GP-based approach that can be used with any black-box classifier and feature attribution method to generate uncertainty estimates that capture the complexity of the black-box model.
Experiments show that capturing this type of uncertainty improves understanding of the explanations and the black-box model itself.
Regarding limitations, GPEC relies on the ability to approximate the black-box DB, which increases the training time. The effect of this is minimized during inference, as the DB only needs to be sampled during training. Additionally, in its current implementation GPEC is limited to black-box classifiers; we leave the extension to regression models as future work.
\section{Experiments}
\vspace{-2mm}
\vspace{-1mm}
\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1\linewidth]{figures/UncertaintyFigure.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-2pt}
   \caption{Visualization of estimated uncertainty of explanations for different models and competing methods. The heatmap represents level of uncertainty for a grid of explanations for the feature on the x-axis; dark areas represent high uncertainty whereas lighter areas represent low uncertainty. Red points represent training samples for GPEC and CXPlain, and the reference samples for BayesSHAP and BayesLIME. The black line represents the decision boundary for the XGBoost binary classifier. }
   \label{fig:uncertaintyfigure}
\end{figure}
We evaluate the uncertainty estimates from GPEC on a variety of datasets and classifiers. In section \ref{sec:experiment_uncertainty_visualization} we compare the DB-aware uncertainty and approximation uncertainty components of GPEC.
Section \ref{sec:experiment_regularizationtest} evaluates how GPEC captures DB complexity.
Section \ref{sec:experiment_fullfat} evaluates how well GPEC combines sources of uncertainty.
Due to space constraints we have additional results in the appendix, including sensitivity analysis for GPEC parameters (\ref{app:sensitivity_analysis}), execution time comparison \ref{app:execution_time}, and additional experiments on combining uncertainty (\ref{app:addnl_results_fullfat}).
All experiments were run on an internal cluster using AMD EPYC 7302 16-Core processors, and all source code will be made public.
\vspace{-2mm}
\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1\linewidth]{figures/figure_fullfat.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-5pt}
   \caption{Comparison of the change in uncertainty estimate of explanations when combining sources of uncertainty with GPEC. Row (A) visualizes the estimate using only the WEG Kernel, which captures the DB-based uncertainty assuming noiseless explanations. Row (B) combines the estimate in (A) with an estimate of approximation uncertainty from the given explainer, BayesSHAP and SSV. The third row ($\Delta$) represents the change in uncertainty estimate between (A) and (B).}
   \label{fig:fullfat}
\end{figure}
\vspace{-1mm}
\subsection{Experimental Setup}
\vspace{-1mm}
\textbf{Datasets and Models.}
Experiments are performed on three tabular datasets (Census, Online Shoppers \citep{sakarRealtimePredictionOnline2019}, German Credit) from the UCI data repository \citep{uci_datasets} and two image datasets (MNIST \citep{lecun-mnisthandwrittendigit-2010} and fashion-MNIST (f-MNIST) \citep{fmnist}).
GPEC can be used with any black-box model; for our experiments we use XGBoost \citep{xgboost} with log-loss for the tabular datasets and a 4-layer Multi-Layered Perception (MLP) model for the image datasets. Additional dataset details are outlined in App. \ref{app:dataset_details}.

\begin{comment}
For the tabular datasets we convert categorical variables to numeric using one-hot-encoding and train XGBoost \citep{xgboost} models with log-loss. For MNIST we train a Multi-Layered Perceptron (MLP) model with the ReLU activation function, 3 hidden layers, and 700 / 400 / 200 nodes in each hidden layer, respectively.
\end{comment}

\textbf{Implementation Details.}
We consider two versions of GPEC: GPEC-WEG which is fit with the WEG kernel, and GPEC-RBF which is fit with the RBF Kernel. Unless otherwise stated, we choose $\lambda = 1.0$ and $\rho = 0.1$ (see App. \ref{app:sensitivity_analysis} for experiments on parameter sensitivity). For the Uncertainty Visualization (Sec. \ref{sec:experiment_uncertainty_visualization}) and Regularization Test (Sec. \ref{sec:experiment_regularizationtest}) we use GPEC with the KernelSHAP explainer.
The DBPS binary search algorithm \citep{DBPS_decision_boundary} is used to sample from the decision boundaries of the black-box models. For tabular datasets we apply DBPS to a grid of samples. For image datasets we first apply adversarial attacks to reduce search space size \citep{karimi2019characterizing}. We then apply ISOMAP \cite{tenenbaumGlobalGeometricFramework2000}, which uses a graph-based approach to estimating geodesic distances. Additional implementation details are available in App. \ref{app:implementation_details}.

\textbf{Competing Methods.}
We compare GPEC to three competing explanation uncertainty estimation methods. BayesSHAP and BayesLIME \citep{reliable:neurips21} are extensions of KernelSHAP and LIME, respectively, that fit Bayesian linear regression models to perturbed data samples. After fitting, $95\%$ credible intervals can be estimated by sampling the posterior distribution of the feature attributions. CXPlain \citep{cxplain} trains a separate explanation model using a causal-based objective, and applies a bootstrap resampling approach to estimate explanation uncertainty. We report the $95\%$ confidence interval from the set of bootstrapped explanations.
Unless otherwise stated, we use the default settings in the provided implementation for all three of these methods.

\begin{comment}
\textbf{Algorithm - GPEC (global)}
1. Calculate Exponential Geodesic Kernel (This step is the global step).
a. Get samples from the decision boundary
b. Get Geodesic Distance between Samples to get EG Kernel

sample training distribution to get {x'}
Using {x'}:
1. Get WEG Kernel:
    a. Calculate EG kernel (M x M matrix)
        aa. Get samples from decision boundary
        ab. Get geodesic distances between samples
    b. Calculate weighting
    c. Calculate WEG kernel (n x n matrix)

2. Get explanations E(x')
3. Train a GP using WEG kernel and (Potentially Noisy) Explanations
4. Approximate predictive distribution for x*. Get the associated variance.

The main difference is in where you get your samples from. For local, you can save computation by sampling locally around x*. For global, you want to sample from the training distribution.

Doing this locally is pretty pointless. We should instead focus on the amortization aspect.

\end{comment}

\vspace{-1mm}
\subsection{Uncertainty Visualization} \label{sec:experiment_uncertainty_visualization}
\vspace{-2mm}
In order to visualize explanation uncertainty, we train XGBoost binary classifiers using two selected features, which we take to be the black-box model.
Both GPEC and competing methods are used to quantify the uncertainty for the feature attributions of the test samples, which we choose to be a grid of 10,000 samples over the data domain. Uncertainty estimates for the x-axis feature are plotted in Figure \ref{fig:uncertaintyfigure} (results for the y-axis variable are included in App. \ref{app:addnl_results_uncertainty_visualization}). Test points are plotted using a heatmap, where darker values indicate high uncertainty.
GPEC-WEG, GPEC-RBF, and CXPlain require training samples for their amortized uncertainty estimates; these samples, also used as the reference distribution for BayesSHAP and BayesLIME, are plotted in red. The DB is plotted in black for reference.
We observe that the uncertainty heatmap for GPEC-WEG attributes higher uncertainty to samples near nonlinearities in the DB.
In contrast, GPEC-RBF provides uncertainty estimates that relate only to the training samples; test sample uncertainty is proportional to distance from the training samples.
The methods BayesSHAP, BayesLIME, and CXPlain results in relatively uniform uncertainty estimates over the test samples. CXPlain shows areas of higher uncertainty for Census, however the magnitude of these estimates are small. The uncertainty estimates produced by these competing methods are unable to capture the properties of the black-box model.

\vspace{-1mm}
\subsection{Regularization Test} \label{sec:experiment_regularizationtest}
\vspace{-2mm}
\input{tables/table_regularization.tex}

In this experiment we compare the average uncertainty of explanations as the model is increasingly regularized in order to assess the impact of restricting model complexity. For XGBoost models, we vary the parameter $\gamma$, which penalizes the number of leaves in the regression tree functions (Eq. 2 in \cite{xgboost}).
For neural networks, we use two types of regularization. First, we add an $\ell_2$ penalty to the weights; this penalty increases as the parameter $\lambda$ increases. Second, we change the ReLU activation functions to Softplus; a smooth approximation of ReLU with smoothness inversely proportional to a parameter $\beta$, which we vary \citep{dombrowskiExplanationsCanBe2019}.
We observe in Table \ref{table:regularize} that the average uncertainty estimate generated by GPEC-WEG decreases as regularization increases, suggesting that the uncertainty estimates reflect the overall complexity of the underlying black-box model. For the tabular datasets, the estimates for BayesSHAP, BayesLIME, and CXPlain stay relatively flat. Interestingly, the estimates from these methods decrease for the image datasets. We hypothesize that the regularization for the neural network model also increases overall stability of the explanations. GPEC can capture both the uncertainty from WEG kernel and also the estimated uncertainty from the function approximation for the explainer, which we demonstrate in Section \ref{sec:experiment_fullfat}.

\vspace{-2mm}
\subsection{Combining Approximation Uncertainty} \label{sec:experiment_fullfat}
\vspace{-2mm}

Using the WEG kernel with noisy explanation labels enables GPEC to incorporate uncertainty estimates from the explainer. Explainers such as BayesSHAP provide an estimate of uncertainty for their explanations, which can be used as fixed noise in the Gaussian likelihood for GPEC. For explainers with no such estimate, we can empirically resample explanations for the same test sample and take the explanation variance as an estimate of uncertainty. Here, we compare GPEC uncertainty results before and after applying explanation noise estimates. We combine GPEC with two different explainers: BayesSHAP and Shapley Sampling Values (SSV) \citep{shapley_sampling}. The two selected explainers are different methods of approximating SHAP values, however the former has an inherent method for quantifying its uncertainty while the latter does not. In Figure \ref{fig:fullfat} row (A) we plot the heatmap for uncertainty estimates using GPEC using only the WEG Kernel and no added noise from the explainers. In row (B) we add the additional uncertainty estimates from BayesSHAP and SSV. The difference in quantified uncertainty is visualized in the third row ($\Delta$). We observe that combining the uncertainty estimates from BayesSHAP and SSV in the GPEC Gaussian Process formulation increases the overall uncertainty estimate.

\section{Uncertainty Quantification for Black-Box Explainers} \label{sec:GP_basic}
\vspace{-2mm}
\begin{comment}
I need to better describe the different steps in the process -- overexplain everything.
How does this capture different sources of uncertainty? Why are we using a GP?
Why do we use our kernel? Geodesic? Projections?
How do we determine what the parameters are?
How do we get uncertainty estimates?
Why don't we just use astuteness or local lipschitzness?
\end{comment}

\begin{comment}
In order to estimate uncertainty of the explanation, we want to understand the behavior of the nearby explanations. Given access to the the explainer $E$, a natural option is to randomly sample the data space around $x$ and observe the behavior of the corresponding explanations. \hl{However, this in itself introduces a source of uncertainty given that only a finite number of explanations can be sampled.} Since we are approximating a potentially high-dimensional space with a finite number of explanations, this task lends itself to the GP framework which allows us to interpolate between sampled explanations while capturing uncertainty as a function of similarity between samples.
\end{comment}

\begin{figure}[t]
    \begin{center}
    \includegraphics[width=0.9\linewidth]{figures/figure1.pdf}
    \end{center}
    \setlength{\belowcaptionskip}{-15pt}
    \caption{Schematic for the GPEC uncertainty estimation process. GPEC can be used in conjunction with a black-box classifier and explainer to derive an estimate of explanation uncertainty. GPEC takes samples from classifier's decision boundary plus (possibly noisy) explanations from the explainer and fits a Gaussian Progress Regression model with Weighted Exponential Geodesic Kernel. The variance of the predictive distribution combines uncertainty from the black-box classifier complexity and the approximation uncertainty from explainer.}
    \label{fig:fig1}
\end{figure}

We now outline the GPEC framework, which we parameterize with a Gaussian Process (GP) regression model.
We define a vector-valued GP which is trained on data samples as input and explanations as labels. More concretely, consider a data sample $x \in \mathcal{X} \subset \mathbb{R}$ that we want to explain in the context of a black-box prediction model $F:\mathcal{X} \to [0,1]$, in which the output corresponds to the probability for the positive class.
For simplicity, we consider the binary classifier case; we extend to the multiclass case in Section \ref{sec:multiclass}. We apply any given feature attribution method $E:\mathcal{X} \to \mathbb{R}^s$, where $s$ is the dimension of the explanations $e = E(x)$.

Let us draw samples $x_1, \dots, x_M$ from the data distribution and generate their respective explanations $e_1, \dots, e_M$.
We assume that each explanation $e_m$ is generated from an unobserved latent function $\mathcal{E}$ plus some independent Gaussian noise $\eta_m$.
\begin{equation*}
\setlength{\abovedisplayskip}{10pt}
   e_m = \mathcal{E}(x_m) + \eta_m \quad s.t. \quad  \quad \underbrace{\mathcal{E}(x_m) \sim \mathcal{GP}(0,k(x, x'))}_{\textrm{Decision Boundary-Aware Uncertainty}} \quad \underbrace{\eta_m \sim \mathcal{N}(0,\tau_m^{-1}  )}_{\textrm{Function Approximation Uncertainty}}\numberthis \label{eqn:GP}
\end{equation*}
where $k(\cdot, \cdot)$ is the specified kernel function and $\tau_m >0 $.
We decompose each explanation into two components, $\mathcal{E}(x_m)$ and $\eta_m$, which represent two separate sources of uncertainty: 1) a \emph{decision boundary-aware }uncertainty which we capture using the kernel similarity, and 2) a \emph{function approximation }uncertainty from the explainer. After specifying $\mathcal{E}(x_m)$ and $\eta_m$, we can combine the two sources by calculating the posterior predictive distribution for the GP model; we take the 95\% confidence interval of this distribution to be the GPEC uncertainty estimate. Due to space constraints, we provide an overview of GP regression in App. \ref{app:GPReview}.

\begin{comment}
We see that the two components  of Eq. \ref{eqn:GP} capture two different sources of uncertainty. $\eta_m$ represents the uncertainty associated with generating the specific explanation for $x_m$. Many explanation methods rely on stochastic methods for approximation which results in noisy explanations; this noise can be captured by the parameter $\tau$ and incorporated into the GP uncertainty estimate.
\end{comment}

\textbf{Function Approximation Uncertainty.} The $\eta_m$ component of Eq. \ref{eqn:GP} represents the uncertainty stemming from explainer specification.
For example, this uncertainty can be captured from  increased variance from undersampling (perturbation-based explainers) or from the increased bootstrap resampling variance from model misspecification (surrogate-based explainers).

We can therefore define $\tau$ based on the chosen explainer. Explainers that include some estimate of uncertainty (e.g. BayesLIME, BayesSHAP, CXPlain) can be directly used to estimate $\tau$. For example, CXPlain returns a distribution of feature importance values for a sample $x_m$; the variance of this distribution can be directly used as the variance of $\eta_m$. For other stochastic explanation methods that do not explicitly estimate uncertainty, we can estimate $\tau$ empirically by resampling explanations for the same data sample:
\vspace{-2mm}
\begin{equation}
   \hat \tau_m = [\frac{1}{|K|}\sum_{i=1}^{K} (E_i(x_{m}) - \bar E(x_m))^2]^{-1} \quad s.t. \quad \bar E(x_m) = \frac{1}{|K|} \sum_{i=1}^K E_i(x_{m})
\end{equation}
where each $E_i(x_m)$ is a sampled explanation for $x_m$. Alternatively, for deterministic explanation methods we can omit the noise and assume that the generated explanations are noiseless.

\textbf{Decision Boundary-Aware Uncertainty.} In contrast, the $\mathcal{E}(x_m)$ component of Eq. \ref{eqn:GP} draws possible functions from the GP prior that could have generated the observed explanations, with the function space defined by the choice of kernel.
Intuitively, the kernel encodes our \emph{a priori} definition of how similar two explanations should be given the similarity of the inputs.
In other words, how much information do we expect a given point $x$ to provide for a nearby point $x'$ with respect to their explanations?
We use this kernel specification to define a source of uncertainty dependent on the behavior of nearby explanations.
In particular, we consider a novel kernel formulation that reflects the complexity of the DB in a local neighborhood of the samples; this is detailed in Section \ref{sec:weighted_geodesic_kernel}.
For any given kernel, we can interpret the distribution of possible functions $\mathcal{E}(x_m)$ as an estimate of uncertainty based on the set of previously observed or sampled explanations.

\begin{comment}
\textbf{Combined GPEC Uncertainty Estimate.}
Once the kernel and noise parameter $\tau_m$ are specified, the two sources of uncertainty can be combined during GP inference. Since the explanations are $s-$dimensional, the specified GP regression model is trained to have vector-valued output; for simplicity we assume that each output dimension is independent. We can then calculate the posterior predictive distribution for any given samples $x^*$, which are conditioned on the previously sampled explanations and their respective uncertainty. We take the 95\% confidence interval as the uncertainty estimate for GPEC.
\end{comment}

\vspace{-2mm}
\section{Weighted Exponential Geodesic Kernel} \label{sec:weighted_geodesic_kernel}
\vspace{-2mm}

Intuitively, the GP kernel encodes the assumption that each individual explanation gives some information about the explanations around it; the choice of kernel defines the neighborhood and magnitude of this shared information.
\begin{comment}
Given a sample $x$, kernel $k(x,x')$ and $\alpha\in (0,1)$, we can define this \emph{explanation neighborhood} as:
\begin{equation}
\mathcal{N}_{\alpha}(x) =\{x' : k(x,x') \geq \alpha \}
\end{equation}
\end{comment}
In the GP framework introduced in Section \ref{sec:GP_basic}, the kernels define the relationship or similarity between two explanations $\mathcal{E}(x)$ and $\mathcal{E}(x')$ based soley on some function of their inputs $x$ and $x'$.
For example, stationary kernels assign the same similarity for any two inputs regardless of where in the data manifold they are located.
Instead, we want to encode the assumption that when the DB is complex, knowing an explanation $\mathcal{E}(x)$ gives limited information about other nearby explanations (see Fig. \ref{fig:kernel_comparison}).
\vspace{-1mm}
\subsection{Geometry and Geodesics of the Decision Boundary}
\vspace{-2mm}
We want to relate kernel similarity to the behavior of the black-box model; specifically, the complexity or smoothness of the DB. Given any two points on the DB, the relative complexity of the boundary segment between them can be approximated by the segment length. The simplest form that the DB can take is a linear boundary connecting the two points; this is exactly the minimum distance between the points. As the complexity of the DB grows, there is a general corresponding increase in segment length.

Concretely, we define the DB as a $(d-1)$-dimensional Riemannian manifold embedded in $\mathbb{R}^d$. Given the binary classifier $F$, we define $\mathcal{M} = \{x \in \mathbb{R}^{d} : F(x) = \frac{1}{2} \}$ representing the DB\footnote{Without loss of generality, we assume that the decision rule for the classifier is set to be $\frac{1}{2}$} of $F$. Given this interpretation, we can define distances along the DB as geodesic distances in $\mathcal{M}$.
\begin{equation}
\setlength{\abovedisplayskip}{3pt}
\setlength{\belowdisplayskip}{3pt}
   d_{geo}(m,m') = \min_{\gamma} \int_{0}^{1} ||\gamma(t)|| \textrm{d}t \quad \;  \forall m,m' \in \mathcal{M}
\end{equation}
The mapping $\gamma:[0,1] \rightarrow \mathcal{M}$, defined such that $\gamma(0) = m$ and $\gamma(1) = m'$, is a parametric representation of a 1-dimensional curve on $\mathcal{M}$.
We can adapt geodesic distance in our kernel selection through the exponential geodesic (EG) kernel \citep{feragen_geodesic_exponential_kernel}, which is a generalization of the Radial Basis Function (RBF) kernel substituting $\ell_2$ distance with geodesic distance:
\begin{equation}
\setlength{\abovedisplayskip}{8pt}
\setlength{\belowdisplayskip}{8pt}
   k_{EG}(x, x') = \exp(-\lambda  d_{geo}(x,x')) \quad \quad k_{RBF}(x, x') = \exp(-\lambda  || x - x'||_2^2)
\end{equation}
The EG kernel has been previously investigated in the context of Riemannian manifolds \citep{feragen_geodesic_exponential_kernel, feragen_open_problem_geodesic_kernel}. In particular, while prior work shows that the EG kernel fails to be positive definite for all values of $\lambda$ in non-Euclidean space, there exists large intervals of $\lambda>0$ for the EG kernel to be positive definite. Appropriate values can be selected through grid search and cross validation; we assume that a valid value of $\lambda$ has been selected.

Note that the manifold we are interested in is the DB; applying the EG kernel on the data manifold (i.e. directly using $k_{EG}$ in the GPEC formulation) would not capture model complexity. In addition, na{\"\i}vely using the EG kernel on the DB manifold would not capture the local complexity with respect to a given explanation; the similarity would be invariant to observed explanations. We therefore need to relate the geodesic distances on $\mathcal{M}$ to samples in the data space $\mathcal{X}$.
\subsection{Weighting the Decision Boundary}
\vspace{-2mm}
\begin{figure}[t]
   \begin{center}
   \includegraphics[width=1\linewidth]{figures/kernel_comparison.pdf}
   \end{center}
   \setlength{\belowcaptionskip}{-5pt}
   \caption{Comparison of RBF and WEG kernel similarity neighborhoods. The gray highlighted region $\mathcal{N}_{0.8}(x) = \{x' : k(x,x') \geq 0.8 \}$ indicate points where the kernel similarity is greater than $0.8$ with respect to the red point $x$. The black line is the decision boundary for the classifier $f(x_1, x_2) = 2 \cos(\frac{10}{x_1}) - x_2$, where $f(x_1, x_2) \geq 0$ indicates class 1, and $f(x_1, x_2) < 0$ indicates class 0. We see that for GPEC, $\mathcal{N}_{0.8}(x)$ shrinks as $x$ moves closer to segments of the decision boundary which are more complex.}
   \label{fig:kernel_comparison}
\end{figure}

Consider a probability distribution $p(M)$ with its support defined over $\mathcal{M}$.
We weight $p(M)$ according to the $\ell_2$ distance between $M$ and a fixed data sample $x$:
\begin{comment}
\begin{equation}
    q(m | x, \rho) = \frac{1}{Z_m} \exp[-\rho || x-m||_2^2] p(m) \quad s.t. \; \;  Z_m = \int \exp[-\rho || x-m||_2^2] p(m) \; \textrm{d}m
\end{equation}
\end{comment}
\begin{equation} \label{eq:q_weighting}
    q(M | x, \rho) \propto  \exp[-\rho || x-M||_2^2] p(M)
\end{equation}
We evaluate the kernel function $k(x,x')$ by taking the expected value over the weighted distribution.
\begin{equation} \label{eq:kernel_weight}
k_{WEG}(x,x') = \int \int \exp[-\lambda d_{geo}(m, m')] \; q(m|x, \rho) \; q(m'|x', \rho) \; \textrm{d} m \textrm{d} m'
\end{equation}
\begin{comment}
\begin{equation}
= \mathbb{E}_{M \sim q(M|x,\rho), M' \sim q(M'|y,\rho)}[\exp[-\lambda d_{geo}(M, M')]]
\end{equation}
\end{comment}
Our formulation is an example of a marginalized kernel \citep{tsuda_marginalizedkernel}: a kernel defined on two observed samples $x,x'$ and taking the expected value over some hidden variable. Given that the underlying EG kernel is positive definite, it follows that the WEG kernel forms a valid kernel.

To evaluate the WEG kernel, we theoretically investigate two properties of the kernel. Theorem \ref{thm:kernel_expgeodesic} shows that the WEG kernel is an extension of the EG kernel for data samples not directly on the DB.
\begin{theorem} \label{thm:kernel_expgeodesic}
   Given two points $x,x' \in \mathcal{M}$, then $\lim_{\lambda \rightarrow \infty} k_{WEG}(x,x') = k_{EG}(x, x')$
\end{theorem}
Proof details are available in App. \ref{proof:points_on_manifold}. Intuitively, as we increase $\lambda$, the manifold distribution closest to the points $x,x'$ becomes weighted increasingly heavily. At the limit, the weighting becomes concentrated entirely on the points $x,x'$ themselves, which recovers the EG kernel. Therefore we see that the WEG kernel is simply a weighting of the EG kernel, which is controlled by $\lambda$.

Theorem \ref{thm:kernel_complexity} establishes the negative relationship between WEG kernel similarity and DB complexity. Note that the notation $\textrm{int}(S)$ indicates the interior of a set $S$, $\left.\textrm{f}\right|_{\mathcal{S}}$ represents the function $f$ restricted to $S$, and $\textrm{id}$ indicates the identity mapping.
\begin{comment}
\begin{definition}
    Let $\cup_{\alpha \in A}\{U_\alpha\}$ be an open cover for a manifold $\mathcal{P} \subset \mathbb{R}^d$. Let $R:\mathcal{P} \subset \mathbb{R}^d \rightarrow \widetilde{\mathcal{P}}$ be a diffeomorphism  s.t. $\forall i \neq \alpha \;$ $\left.R\right|_{\mathcal{P} \setminus K_i} = \left.\textrm{id}\right|_{\mathcal{P} \setminus K_i}$ and $\left.R\right|_{K_i} \neq \left.\textrm{id}\right|_{K_i}$ where $K_i
    \subseteq U_i$ and is compact. We define the manifold $\widetilde{\mathcal{P}} = R(\mathcal{P}) \subset \mathbb{R}^d$ as a perturbation of $\mathcal{P}$.
\end{definition}
\end{comment}
\begin{definition} \label{def:1}
Let $\{U_\alpha\}_{\alpha\in I}$ be charts of an atlas for a manifold $\mathcal{P} \subset \mathbb{R}^d$, where $I$ is a set of indices. Let $\mathcal{P}$ and $\widetilde{\mathcal{P}}$ be differentiable manifolds embedded in $\mathbb{R}^d$, where $\mathcal{P}$ is a Piecewise Linear (PL) manifold. Let  $R: \mathcal{P} \rightarrow \widetilde{\mathcal{P}}$ be a diffeomorphism. We say $\widetilde{\mathcal{P}}$ is a perturbation of $\mathcal{P}$ on the $i^{\textrm{th}}$ chart if $R$ satisfies the following conditions. There exists a compact subset $K_i \subset U_i$ s.t.  $\left.R\right|_{\mathcal{P} \setminus \textrm{int}(K_i)} = \left.\textrm{id}\right|_{\mathcal{P} \setminus \textrm{int} (K_i)}$ and $\left.R\right|_{\textrm{int}(K_i)} \neq \left.\textrm{id}\right|_{\textrm{int}(K_i)}$. Furthermore there exists a linear homeomorphism between an open subset $\widetilde{U_i} \subseteq U_i$ with $\mathbb{R}^{d-1}$ which contains $K_i$.
\end{definition}

\begin{theorem} \label{thm:kernel_complexity}
Let $\mathcal{P}$ be a $(d$--$1)$-dimension PL manifold embedded in $\mathbb{R}^d$. Let $\widetilde{\mathcal{P}}$ be a perturbation of $\mathcal{P}$ and define $\tilde k(x,x')$ and $k(x,x')$ as the WEG kernel defined on $\widetilde{\mathcal{P}}$ and $\mathcal{P}$ respectively. Then $\tilde k(x,x') < k(x,x') \; \forall x,x' \in \mathbb{R}^d$.

\end{theorem}
\begin{comment}
\begin{corollary}
    . Under the assumption that the local perturbations do not intersect, then $\mathcal{P}$
\end{corollary}

\begin{theorem} \label{thm:kernel_complexity}
Given a 1-dimensional manifold $\mathcal{M}$ and two points $x,x' \in \mathbb{R}^2$. For any two points $a,b \in \mathcal{M}$, if we replace the manifold between $a$ and $b$ with some function $g$, the function $g$ that maximizes $k_{\mathcal{M}}(x,x')$ is a straight line between $a$ and $b$.
\end{theorem}
\end{comment}

Proof details are
in App. \ref{proof:kernel_complexity}. Theorem \ref{thm:kernel_complexity} implies that the linearization of the DB maximizes the kernel similarity; any deviation from the linearization will correspondingly decrease kernel similarity.

\begin{comment}

\begin{remark}
    Replacing a given segment of the decision boundary with a straight line results in a linearization. Other papers have shown that linearization of the decision boundary results in lower model complexity.
\end{remark}

\begin{corollary}
Any function other than the straight line will decrease $k(x,y)$
\end{corollary}

\begin{theorem}
   Given a manifold $\mathcal{M}$, two points $a,b \in \mathcal{M}$, and two points $x,y \in \mathbb{R}^d$. We define $\mathcal{M}$ as the manifold $\mathcal{M}$ with the interval between $[a,b]$ replaced with a straight line. Then $k_\mathcal{M'}(x,y) \geq k_{W}(x,y)$ with equality when $\mathcal{M} = \mathcal{M}'$.
\end{theorem}
\end{comment}

To improve the interpretation of the WEG kernel, we apply a normalization to scale the similarity values to be between $[0,1]$. We construct the normalized kernel $k_{WEG}^*$ as follows:
\vspace{-2mm}

\begin{equation}
    k_{WEG}^*(x,x') = \frac{k_{WEG}(x,x')}{\sqrt{k_{WEG}(x,x) k_{WEG}(x',x')}}
\end{equation}

\vspace{-2mm}
\subsection{WEG Kernel Approximation}
\vspace{-1mm}
In practice, the integral in Eq. \ref{eq:kernel_weight} is intractable; we approximate the expected value using Monte Carlo (MC) sampling with $K$ samples:
\begin{equation}
    k_{WEG}(x,x') \approx \frac{1}{Z_m Z_{m'} K^2}\sum_{i=1}^K \sum_{j=1}^K \exp[-\lambda d_{geo}(m_i, m_j)] \exp[-\rho (|| x-m_i||_2^2 + || x'-m_j||_2^2)]
\end{equation}
$Z_m$, $Z_{m'}$ are the normalization constants for $q(M|x,\rho)$ and $q(M'|y,\rho)$, respectively. We can similarly estimate these values using MC sampling:
\vspace{-3mm}
\begin{equation}
    Z_m = \int \exp[-\rho || x-m||_2^2] p(m) \; \textrm{d}m  \approx \frac{1}{K}\sum_{i=1}^K \exp[-\rho || x-m_i||_2^2]
\end{equation}
\begin{comment}
Putting it all together:
\begin{equation}
    \tilde k(x,y) = \frac{\sum_{i=1}^K \sum_{j=1}^K \exp[-\lambda d_{geo}(m_i, m_j)] \exp[-\rho (|| x-m_i||_2^2 + || y-m_j||_2^2)]}{(\sum_{i=1}^K \exp[-\rho || x-m_i||_2^2])(\sum_{j=1}^K \exp[-\rho || y-m_j||_2^2])}
\end{equation}
\end{comment}
\textbf{GPEC Implementation Details.}
To use the GPEC framework, we first train the GP regression model with a set of training samples. Once the GP is trained, we can derive uncertainty estimates for test samples from the 95\% confidence interval of the posterior predictive distribution.
The WEG kernel requires additional setup time to use in the GP. We first calculate the EG kernel matrix by sampling the black-box DB, then weigh each element of the matrix per eq. \ref{eq:q_weighting}. The estimation of a classifier's DB is an ongoing area of research; a number of methods have been proposed in the literature (e.g. \cite{DBPS_decision_boundary, karimi2019characterizing}).
Once GPEC is trained, training cost is amortized during inference; estimating uncertainty for test samples using a GP generally has time complexity of $\mathcal{O}(n^3)$, which can be reduced to $\mathcal{O}(n^2)$ using BBMM \citep{gpytorch}, and further with variational methods (e.g., \cite{gal_variational_GP}).

\begin{comment}
We can simplify the decision boundary estimation in $k_{\mathcal{M}}$ in two ways.
First, by assuming a two-class prediction from the black-box model $F$. This assumption aligns with many feature attribution methods, which simplify $F$ to a one-vs-all classifier in the multiclass case.
Second, because the weight components of $k_{\mathcal{M}}$ (Eq. \ref{eq:q_weighting}) decreases proportionally with distance to the manifold, the decision boundary only needs to be estimated locally. Given a threshold $\beta \approx 0$, we can limit the decision boundary search to a set $\mathcal{B}_{\beta}(x, x') = \{y: \exp[-\gamma ( ||x-y||^2 + ||x'-y||^2)] \geq \beta \}$.
\end{comment}

\section{Related Works}
\vspace{-2mm}
\textbf{Explanation Methods.}
A wide variety of methods have been proposed for the purpose of improving transparency for pre-trained black-box prediction models \citep{guidotti_survey, arrietta_survey}.
Within this category of \emph{post-hoc} methods, many methods focus on \emph{local} explanations, that is, explaining individual predictions rather than the entire model.
Some of these methods generate explanations through local \emph{feature selection} \citep{l2x_chen, IFG}. In this work, we focus primarily on \emph{feature attribution} methods, which return a real-valued score for each feature in the sample. For example, LIME \citep{ribeiroLIME} trains a local linear regression model to approximate the black-box model.
\cite{shap_lundberg} generalizes LIME and five other feature attribution methods using the SHAP framework, which fulfill a number of desirable axioms.
While some methods such as LIME and its SHAP variant, KernelSHAP, are model-agnostic, others are designed for specific model architectures, such as neural networks \citep{LRP_saliency, pmlr-v70-shrikumar17a, sundararajanAxiomaticAttributionDeep2017, erionImprovingPerformanceDeep2021}, tree ensembles \citep{lundbergLocalExplanationsGlobal2020}, or Bayesian neural networks \citep{bykov2020much}, taking advantage of those specific architectures. Another class of method involves the training of  separate surrogate models to  explain the black-box model \citep{dabkowski_gal, l2x_chen, cxplain, guo_explaining_nonparametric_approach, jethani2022fastshap}. Once trained, surrogate-based methods are typically very fast; explanation generation generally only requires a single inference step from the surrogate model.

\textbf{Reliability of Explanations.}
While feature attribution methods have gained wide popularity, a number of issues relating to the reliability of such methods have been uncovered. \cite{alvarez-melis2018} investigate the notion of robustness and show that many feature attribution methods are sensitive to small changes in input. This has been further investigated in the adversarial setting for perturbation-based methods \citep{slack_foolingLIMESHAP} and neural network-based methods. \citep{ghorbaniInterpretationNeuralNetworks2019}. \cite{Kindermans2019} show that many feature attribution methods are affected by distribution transformations such as those common in preprocessing. The generated explanations can also be very sensitive to hyperparameter choice \cite{bansalSAMSensitivityAttribution2020}.
A number of metrics have been proposed for evaluating explainer reliability, such as with respect to adversarial attack \citep{hsiehEvaluationsMethodsExplanation2021}, local perturbations \citep{alvarez-melis2018, VSI_visani2022}, black-box smoothness \citep{zulqarnain_lipschitzness}, fidelity to the black-box model \citep{yehFidelitySensitivityExplanations2019}, or combinations of these metrics \citep{bhattEvaluatingAggregatingFeaturebased2020}.

\textbf{Explanation Uncertainty.}
One option for improving the trustworthiness of explainers is to quantify the associated explanation uncertainty. Bootstrap resampling techniques have been proposed as a way to get uncertainty estimates from surrogate-based explainers \citep{cxplain, schulz2022uncertainty}. \cite{guo_explaining_nonparametric_approach} also proposes a surrogate explainer parameterized with a Bayesian mixture model.
Alternatively, \cite{bykov2020much} introduces a method for explaining Bayesian neural networks, which can be transferred to non-Bayesian neural networks.
\cite{improving_kernelshap} derives an unbiased version of KernelSHAP and investigates an efficient way of estimating its uncertainty.
Several methods investigate LIME and KernelSHAP in a Bayesian context; for example calculating a posterior over attributions \citep{reliable:neurips21}, investigating the use of priors for explanations \citep{bayelime_zhao21}, or using active learning to improve sampling \citep{select_wisely_and_explain_Saina2022}.

\begin{comment}
Our proposed method complements these other approaches by taking a chosen explainer, possibly with an associated uncertainty estimate, and combining them with a novel decision boundary-based notion of uncertainty. Using GPEC with another explainer therefore improves understanding of the black-box model.
\end{comment}
However, all existing methods for quantifying explanation uncertainty only consider the uncertainty of the explainer as a function approximator.
This work introduces an additional notion of uncertainty for explainers that takes into account the uncertainty of the explainer due to the DB.

\section{Introduction}
\vspace{-2mm}

Machine learning models are becoming increasingly prevalent in a wide variety of industries and applications.
In many such applications, the best performing model is opaque; post-hoc explainability methods are one of the crucial tools by which we understand and diagnose the model's predictions.
Recently, many explainability methods, termed \emph{explainers}, have been introduced in the category of local feature attribution methods.
That is, methods that return a real-valued score for each feature of a given data sample, representing the feature's relative importance with respect to the sample prediction.
These explanations are local in that each data sample may have a different explanation.
Using local feature attribution methods therefore helps users better understand nonlinear and complex black-box models, since these models are not limited to using the same decision rules throughout the data distribution.

Recent works have shown that existing explainers can be inconsistent or unstable.  For example, given similar samples, explainers might provide different explanations \citep{alvarez-melis2018,slack_foolingLIMESHAP}.  When working in high-stakes applications, it is imperative to provide the user with an understanding of whether an explanation is reliable, potentially problematic, or even misleading.  A way to guide users regarding an explainer's reliability is to provide corresponding \emph{uncertainty quantification} estimates.

One can consider explainers as function approximators; as such, standard techniques for quantifying the uncertainty of estimators can be utilized to quantify the uncertainty of explainers.  This is the strategy utilized by existing methods for producing uncertainty estimates of explainers \citep{reliable:neurips21,cxplain}.
However, we observe that for explainers, this is not sufficient; because in addition to uncertainty due to the function approximation of explainers, explainers also have to deal with the uncertainty due to the complexity of the decision boundary (DB) of the black-box model in the local region being explained.

\begin{comment}
A number of explainability methods have been proposed, which capture different aspects of an opaque model's inner workings.
In particular, several methods \cite{ribeiroLIME} use a simplified linear model to locally approximate the complex decision boundaries of the opaque model for a given data sample.
These methods generate perturbed test samples and solve a weighted least-squares approximation on the corresponding model outputs. Followup works have evaluated the sensitivity of this linear approximation to sample size, kernel weighting, and reference baseline parameters \cite{}. \cite{reliable:neurips21} has extended LIME to a Bayesian setting, with the approximation of a Bayesian posterior for the parameters of the simplified linear model. This Bayesian posterior is interpreted as an approximation uncertainty given the specified number of perturbed samples. which they use to establish guidelines on sample size parameter.

However, existing linear methods fail to capture the complexity of the opaque model in the locality of the test sample, as defined by the neighborhood of sampled perturbations around the test sample. While LIME allows the user to specify neighborhood size through the sampling distribution and kernel weighting function, the linearization effectively averages the opaque model's nonlinearities. In addition, without \emph{a priori} knowledge of the model's topology, it is unclear how the neighborhood size should be chosen.

\end{comment}

Consider the following example: we are using a prediction model for a medical diagnosis using two features, level of physical activity and body mass index (BMI) (Fig. \ref{fig:example1}).
In order to understand the prediction and give actionable recommendations to the patient, we use a feature attribution method to evaluate the relative importance of each feature.
Because of the nonlinearity of the prediction model, patients A and B show very similar symptoms, but are given very different recommendations.
Note that while this issue is related to the notion of explainer uncertainty, measures of uncertainty that only consider the explainer would not capture this phenomenon.
This suggests that any notion of uncertainty is incomplete without capturing information related to the local behavior of the model.
Therefore, the ability to quantify uncertainty for DB-related explanation instability is desirable.

\begin{figure}[t]
    \begin{center}
    \includegraphics[width=0.8\linewidth]{figures/example1.pdf}
    \end{center}
    \setlength{\belowcaptionskip}{-15pt}
    \caption{Illustrative example of how similar data samples can result in very different feature importance scores given a black-box model with nonlinear decision boundary. Here, two similar patients with similar predictions are given opposing feature importance scores which could result in misguided recommendations. We define a similarity based on the geometry of the decision boundary between any two given samples (red line). While the two patients are close together in the Euclidean sense, they are dissimilar under the proposed WEG kernel similarity. Using GPEC would return a high uncertainty measure for the explanations, which would flag the results for further investigation.}
    \label{fig:example1}

\end{figure}

We approach this problem from the perspective of similarity: given two samples and their respective explanations, how closely related should the explanations be? From the previous intuition, we define this similarity based on a geometric perspective of the DB complexity between these two points.
Specifically, we propose a novel geodesic-based kernel similarity metric, which we call the Weighted Exponential Geodesic (WEG) kernel.
The WEG kernel encodes our expectation that two samples close in Euclidean space may not actually be similar if the DB within a local neighborhood of the samples is highly complex.

Using our similarity formulation, we propose the \textbf{G}aussian \textbf{P}rocess \textbf{E}xplanation Un\textbf{C}ertainty (GPEC) framework, which is an instance-wise, model-agnostic, and explainer-agnostic method to quantify the uncertainty of explanations.
The proposed notion of uncertainty is complementary to existing quantification methods. Existing methods primarily estimate the uncertainty related to the choice in model parameters and fitting the explainer, which we call \emph{function approximation uncertainty}, and does not capture uncertainty related to the DB.
GPEC can combine the DB-based uncertainty with function approximation uncertainty derived from any local feature attribution method.

In summary, we make the following contributions:
\vspace{-0.5mm}
\setlist{nolistsep}
\begin{itemize} [leftmargin=0.3cm]
    \itemsep0.4em
    \item We introduce a geometric perspective on capturing explanation uncertainty and define a novel geodesic-based similarity between explanations. We prove theoretically that the proposed similarity captures the complexity of the decision boundary from a given black-box classifier.
    \item We propose a novel Gaussian Process-based framework that combines A) uncertainty from decision boundary complexity and B) explainer-specific uncertainty to generate uncertainty estimates for any given feature attribution method and black box model.
    \item Empirical results show GPEC uncertainty improves understanding of feature attribution methods.
\end{itemize}