#  _ADS_ : Approximate Densest Subgraph for Novel Image Discovery ††thanks:
Citation: Authors. Title. Pages…. DOI:000000/11111.

Shanfeng Hu
Northumbria University
Newcastle upon Tyne, UK
<EMAIL_ADDRESS>

###### Abstract

The volume of image repositories continues to grow. Despite the availability
of content-based addressing, we still lack a lightweight tool that allows us
to discover images of distinct characteristics from a large collection. In
this paper, we propose a fast and training-free algorithm for novel image
discovery. The key of our algorithm is formulating a collection of images as a
perceptual distance-weighted graph, within which our task is to locate the
$K$-densest subgraph that corresponds to a subset of the most unique images.
While solving this problem is not just NP-hard but also requires a full
computation of the potentially huge distance matrix, we propose to relax it
into a $K$-sparse eigenvector problem that we can efficiently solve using
stochastic gradient descent (SGD) without explicitly computing the distance
matrix. We compare our algorithm against state-of-the-arts on both synthetic
and real datasets, showing that it is considerably faster to run with a
smaller memory footprint while able to mine novel images more accurately.

_K_ eywords novel image, K-densest subgraph, sparse eigenvector, stochastic
gradient descent

##  1 Introduction

We are naturally attracted by visually novel scenes in our daily lives [1, 2].
While various photo management softwares can provide the functionalities of
chronologically and semantically indexing these visual assets for easy access,
they still lack the basic ability to let a user navigate some of the most
memorable images she has taken before. This feature may not just satisfy our
own curiosity of recalling "old surprises" but can also prompt sufficiently
novel images to be shared on social media for views attraction. Posting less
interesting images, in contrast, can lead to lower user engagement on social
media [3]. Despite the availability of various content-based retrieval of
images using state-of-the-art deep learning technologies [4, 5], we still lack
a lightweight and efficient tool to automatically mine images of distinct
nature.

Related Work. While assessing whether an image is of great interest or not is
subject to personal perception to some extend, there are work on image
memorability showing that the judgements of human participants are actually
quite consistent among each other [6, 7]. This has since then sparked a range
of efforts that attempt to predict the memorability of images using machine
learning and deep learning techniques [6, 8, 9, 10, 11, 12, 13, 14]. Despite
the success, these learning-based methods can suffer a few limitations. First,
they require a ground-truth set of human-annotated memorability scores for
training, which can be costly to collect. Second, the predictive power of the
trained models depends on whether test images are sufficiently similar to
those in the training set. As new images arrive, distributional shifts can
happen and as a result incur prediction errors. And third, the existing
methods make predictions on an image-by-image case, without comparing the
relative novelty of images as a whole. However, this problem would be more
natural to solve as a ranking rather than a regression task.

Visual saliency detection is another area of research that could help identify
images that significantly deviate from others within a collection. However,
most saliency detection methods only look for unique regions or objects within
individual images [15, 16, 17, 18, 19], which lack the objective of assessing
the uniqueness of whole images. While a variety of group-level saliency
detection approaches have been proposed [20, 21, 22], they are mainly designed
to search for co-occurring prominent objects of an image collection. Salient
features have also been used for image memorability prediction [23, 9], but it
remains unclear how to formalize the relationship between local saliency and
global ranking of image novelties.

Our Contributions. In view of the aforementioned challenges, we propose a
lightweight algorithm for novel image discovery that is mathematically clear
to understand, practically simple to implement, fast to run, and free of
training (hence requiring no human annotated labels). Different from the
existing work on image memorability prediction, we formulate the problem as
searching for a global ranking of images in a collection that maximizes a
total novelty score. The key idea of our algorithm is representing a
collection of images as an undirected complete graph, with the nodes being the
images and each edge weighted using the perceptual distance of the two
corresponding images. This allows us to identify a subset of the $K$ most
distinct images by finding the $K$-densest subgraph within the complete graph
[24]. The total novelty score is the sum of the mutual distances among
selected images.

Due to the NP-hardness of the $K$-densest subgraph problem [24], a wide range
of approximate algorithms have been proposed to solve it efficiently [25, 26,
27]. However, a main limitation of these methods is that they require a full
graph representation of data, which in our case would mean that we need to
pre-compute the pairwise distance matrix of all images in a collection for
novelty finding. If we have 10,000 images in the collection, then we will have
a $10000\times 10000$ dense matrix that is trivially intractable to handle on
a resources-limited device such as our phone. To address this challenge, we
propose to relax the original problem into a half-discrete, half-continuous
version, which looks for a $K$-sparse eigenvector that maximizes a weighted
sum of the mutual distances among the selected images. The higher a weight (an
element in the eigenvector), the more novel the corresponding image. A zero
weight indicates that the image is not selected.

The proposed relaxation allows us to solve the problem much more efficiently
using SGD, without explicitly computing or storing the full distance matrix of
images. The key step is that we can approximate the full gradient evaluation
at each step with a Monte-Carlo estimate, which only requires sampling a mini-
batch of images for comparison with each image to update its weight. While SGD
and its variants are de facto optimizers for deep learning [28], we further
equip it with a clipping operator that prunes unpromising images early on
while ensuring the $K$-sparsity constraint of the solved eigenvector. Due to
the scalability of SGD, we are able to efficiently discover a required amount
of novel images from a potentially large image repository.

##  2 Novel Image Discovery

###  2.1 Problem Formulation

Suppose we are given a set of $N$ images,
$\mathcal{C}=\\{\bm{x}_{i}\\}_{i=1}^{N}$. For each image $\bm{x}_{i}$ we are
able to extract a $M$-dimensional feature vector,
$\phi(\bm{x}_{i})\in\mathbb{R}^{M}$, which represents the visual and semantic
contents of the image. The problem of novel image discovery can then be
formulated as finding a subset of at most $K$ images
$\mathcal{S}(\mathcal{C})=\\{\bm{x}_{(k)}\\}_{k=1}^{K}$ within the collection
$\mathcal{C}$ that are the most distinct from each other as well as from the
remaining ones in the feature space. As explained in Section 1, the existing
work on image memorability prediction treats this problem as an individual-
based regression task [6, 10, 11, 12, 13], without optimizing the global
ranking of images within a collection. Instead, we formalize it as a
mathematical optimization problem:

$\displaystyle\mathcal{S}(\mathcal{C})=\underbrace{\underset{\mathcal{S}}{\text{arg
max\ }}J(\mathcal{S})\text{\ s.t.
$\mathcal{S}\subset\mathcal{C}\land|\mathcal{S}|\leq K$}}_{\text{novel image
discovery problem}}$ (1)

with $|\mathcal{S}|$ being the cardinality of the solution set and the
objective function $J(\mathcal{S})$ is defined as

$J(\mathcal{S})=\underbrace{\sum_{\bm{x}_{i}\in\mathcal{S}}{\sum_{\bm{x}_{j}\in\mathcal{S}}{d(\bm{x}_{i},\bm{x}_{j})}}}_{\text{total
novelty score}}$ (2)

in which we use the Euclidean distance in the feature space to measure the
closeness for every pair of images

$d(\bm{x}_{i},\bm{x}_{j})=\|\phi(\bm{x}_{i})-\phi(\bm{x}_{j})\|$ (3)

We refer to the objective (2) as the _total novelty score_ as it measures the
sum of the mutual feature distances among all selected images. The higher the
score, the more separated the images in $\mathcal{S}$ from each other. We are
searching for the $K$ most unique images within $\mathcal{C}$ to maximize the
score.

If we view the given collection of images $\mathcal{C}$ as an undirected
complete graph in which each image is a node and the edge between every pair
of images is weighted using their feature distance (3), then the problem (1)
amounts to solving for the $K$-densest subgraph of the complete graph [24].
This problem is well known for its NP-hardness, for which various approximate
solutions have been proposed in the past [25, 29, 26, 30]. However, these
methods generally require that every pair of data points need to be compared
in the solution process, which may not be feasible when the underlying set
$\mathcal{C}$ is very large (e.g., containing over 10,000 images). This
presents a particular computational hurdle on resources-limited devices such
as our phone. We address this challenge in Section 2.2 using our proposed
algorithm.

###  2.2 The Proposed Algorithm

Here, we present a simple and scalable algorithm for approximately solving
problem (1) to find the subset of most unique images within a collection. We
start with relaxing the problem into a half-discrete, half-continuous version
in Section 2.2.1, which we then efficiently solve using SGD and sparsity
clipping in Section 2.2.2. We analyze the computational complexity of our
algorithm in Section 2.3. We give the pseudocode in Algorithm (1).

Algorithm 1 Approximate Densest Subgraph (ADS)

Input: $\mathcal{C}=\\{\bm{x}_{i}\\}_{i=1}^{N}$, a collection of $N$ images.

Parameters: $K$, number of novel images to be found from $\mathcal{C}$;
$\phi$, feature extractor; $E$, epoch size (default: 20); $J$, batch size
(default: 16); $\alpha$, learning rate (default: 0.001); $\beta$, momentum
(default: 0.9).

Output: $\bm{s}(\mathcal{C})$, a $K$-sparse vector of novelty weights.

1: $\\{\phi(\bm{x}_{i})\\}_{i=1}^{N}\xleftarrow{\phi}\mathcal{C}$

2: Let $t=1$, $\bm{s}^{0}=\frac{1}{\sqrt{N}}[1,1,\cdots,1]^{T}$,
$\Delta^{0}_{i}=0$ for $1\leq i\leq N$

3: while $t\leq E$ do

4: Let $i=1$

5: while $i\leq N$ do

6:
$\mathcal{J}^{t}_{i}\leftarrow\text{random\\_choice}(\\{i|\bm{s}_{i}^{t-1}>0\\},J)$

7: Let $\Delta_{i}^{t}=0$

8: for $j\in\mathcal{J}^{t}_{i}$ do

9:
$\Delta_{i}^{t}\leftarrow\Delta_{i}^{t}+\|\phi(\bm{x}_{i})-\phi(\bm{x}_{j})\|\bm{s}^{t-1}_{j}$

10: end for

11: if $t\geq 2$ then

12: $\Delta_{i}^{t}\leftarrow(1-\beta)\Delta_{i}^{t}+\beta\Delta_{i}^{t-1}$

13: end if

14:
$\bm{s}^{t}_{i}\leftarrow\bm{s}^{t-1}_{i}+2\alpha\frac{N}{J}\Delta_{i}^{t}$

15: $\Delta_{i}^{t-1}\leftarrow\Delta_{i}^{t}$

16: $i\leftarrow i+1$

17: end while

18: Let $K^{t}=\lfloor N-\frac{N-K}{E}t\rfloor$

19:
$\bm{s}^{t}\leftarrow\bm{s}^{t}\circ[I(1\in\text{top-}K^{t}(\bm{s}^{t})),\cdots,I(N\in\text{top-}K^{t}(\bm{s}^{t}))]^{T}$

20: $\bm{s}^{t}\leftarrow\frac{\bm{s}^{t}}{\|\bm{s}^{t}\|}$

21: $t\leftarrow t+1$

22: end while

23: $\bm{s}(\mathcal{C})\leftarrow\bm{s}^{T}$

24: Return $\bm{s}(\mathcal{C})$

####  2.2.1 Sparse Continuous Relaxation

Let us assume for the moment that we can afford to compute and store the full
distance matrix of images in the collection $\mathcal{C}$,
$D\in\mathbb{R}_{\geq 0}^{N\times N}$, with each element
$D_{ij}=d(\phi(\bm{x}_{i}),\phi(\bm{x}_{j}))$ as computed in equation (3). It
is easy to see that problem (1) can be equivalently expressed as the following
one:

$\bm{s}(\mathcal{C})=\underbrace{\underset{\bm{s}}{\text{arg max\
}}\bm{s}^{T}D\bm{s}\text{\ s.t. $\bm{s}\in\\{0,1\\}^{N}\land|\bm{s}|_{0}\leq
K$}}_{\text{reformulation of (\ref{eqn:discrete-problem}) using matrix form}}$
(4)

in which the optimization variable $\bm{s}$ is restricted to be a vector of 0
and 1 of length $N$, with $\bm{s}_{i}=1$ indicating
$\bm{x}_{i}\in\mathcal{S}(\mathcal{C})$ and $\bm{s}_{i}=0$ otherwise. The
objective
$J(\bm{s})=\bm{s}^{T}D\bm{s}=\sum_{i=1}^{N}{\sum_{j=1}^{N}{D_{ij}\bm{s}_{i}\bm{s}_{j}}}$
hence equals to the original one in (2). The L0 norm of the vector,
$|\bm{s}|_{0}$, namely the number of non-zero elements in the vector, is
required to be no greater than the maximal number of selected images $K$.

While the new formulation (4) is still NP-hard to solve, we now have the
opportunity to relax it into a half-discrete, half-continuous version by
allowing each element of the solution vector to take on non-negative real
values:

$\bm{s}(\mathcal{C})=\underbrace{\underset{\bm{s}}{\text{arg max\
}}\bm{s}^{T}D\bm{s}\text{\ s.t. $\bm{s}\in\mathbb{R}_{\geq
0}^{N}\land\|\bm{s}\|=1\land|\bm{s}|_{0}\leq K$}}_{\text{sparse continuous
relaxation of (\ref{eqn:discrete-problem-matrix})}}$ (5)

In this case, $\bm{s}_{i}>0$ would mean that the $i$-th image is selected into
the subset $\mathcal{S}(\mathcal{C})$ with a certain weight while
$\bm{s}_{i}=0$ indicates that the image is not selected. Since $K$ is normally
chosen to be much smaller than $N$, the solution is therefore enforced to be
sparse while maximizing the weighted sum of the mutual feature distances among
the selected images. Intuitively speaking, the more novel an image
$\bm{x}_{i}$ given that it is selected, the higher its score
$\bm{s}_{i}(\mathcal{C})$, providing a fine-grained ranking of images within
$\mathcal{S}(\mathcal{C})$ according to their relative uniqueness. Solving the
original problem (1) or (4), however, does not enjoy this feature as the
solution would be binary with no ranking information generated.

It is interesting to note that when $K=N$ (i.e., the sparsity constraint is
dropped), the relaxed problem (5) is further reduced to solving for the unit-
length principal eigenvector of the distance matrix $D$ that is associated
with the largest eigenvalue. The Perron–Frobenius theorem [31] states that due
to the symmetry and non-negativity of $D$, the solution is guaranteed to be
non-negative hence automatically satisfying the constraint
$\bm{s}\in\mathbb{R}_{\geq 0}^{N}$. The principal eigenvector can be found
through the classical power iteration:
$\bm{s}^{t}=\frac{D\bm{s}^{t-1}}{\|D\bm{s}^{t-1}\|}$, for which $\bm{s}^{0}$
can be initialized to be the unit-length uniform vector and the iteration can
be stopped at step $t$ when $\|\bm{s}^{t}-\bm{s}^{t-1}\|<\epsilon$, with
$\epsilon$ being a small constant for convergence checking. However, this
requires the $N\times N$ dense distance matrix to be fully instantiated for
iteration, which can be costly when $N$ is large. Also, the method cannot
guarantee the $K$-sparsity of the sought solution.

####  2.2.2 Stochastic Approximation

Here, we address the aforementioned challenges of power iteration (Section
2.2.1) with the following steps. First, we notice that by temporarily dropping
the constraints in problem (5), we can locally maximize the objective by
following the gradient direction:

$\displaystyle\bm{s}^{t}$
$\displaystyle=\bm{s}^{t-1}+\alpha\frac{\partial(\bm{s}^{T}D\bm{s})}{\partial\bm{s}}|_{\bm{s}^{t-1}}$
(6) $\displaystyle=\bm{s}^{t-1}+2\alpha D\bm{s}^{t-1}$ (7)

where $\alpha>0$ is the learning rate (or step size). When $\alpha$ is
sufficiently large, the update rule (7) has a roughly similar effect of
running one step of power iteration. While the gradient evaluation at each
step is exact, it involves the full distance matrix and can be computationally
heavy. We choose to approximate it using an Monte-Carlo estimate of the
gradient for each component of the solution vector:

$\displaystyle\bm{s}^{t}_{i}$
$\displaystyle=\bm{s}^{t-1}_{i}+2\alpha\sum_{j=1}^{N}{D_{ij}\bm{s}^{t-1}_{j}}$
(8)
$\displaystyle\approx\underbrace{\bm{s}^{t-1}_{i}+2\alpha\left[(1-\beta)\frac{N}{|\mathcal{J}_{i}^{t}|}\sum_{j\in\mathcal{J}_{i}^{t}}{D_{ij}\bm{s}^{t-1}_{j}}+\beta\Delta^{t-1}_{i}\right]}_{\text{stochastic
gradient descent with momentum}}$ (9)

where the step from (8) to (9) is done by replacing the costly exact sum with
the rescaled sum over a small mini-batch of randomly picked images
$\mathcal{J}_{i}^{t}\subset\\{1,2,\cdots,N\\}$ for image $\bm{x}_{i}$ at step
$t$. We also maintain an exponential moving average of the historical sums
$\Delta^{t-1}_{i}$ for each image $\bm{x}_{i}$ to stabilize the SGD learning
process. $\beta$ controls the strength of the momentum.

The intuition of (9) is that at each step $t$ when we update the $i$-th
image’s novelty weight $\bm{s}^{t-1}_{i}$, we randomly draw a small sample of
images from the collection $\mathcal{C}$ and compare each of them with
$\bm{x}_{i}$. If an image $\bm{x}_{j}$ for $j\in\mathcal{J}_{i}^{t}$ has a
large distance from $\bm{x}_{i}$ and its weight $\bm{s}^{t-1}_{j}$ is also
high, then we increase $\bm{s}^{t-1}_{i}$ by a large amount in the update
because the image $\bm{x}_{i}$ is found to be distinct from $\bm{x}_{j}$ that
is also likely to be novel. Otherwise, when either the two images are very
close in the feature space or if $\bm{x}_{j}$ is not sufficiently distinct
from the previous step, we only increase the weight by a comparatively small
size to indicate that $\bm{x}_{i}$ may not be a good candidate of novelty
given the evidence of the current step.

To impose the $K$-sparsity constraint, we apply the following clipping
operation after executing each step of (9):

$\bm{s}^{t}_{i}=\begin{cases}\bm{s}^{t}_{i}\text{\ if
$i\in\text{top-}K^{t}(\bm{s}^{t})$}\\\ 0\text{\ otherwise}\end{cases}$ (10)

where $\text{top-}k^{t}(\cdot)$ returns the indices of the $K^{t}$ largest
elements of a given vector. After truncating the remaining least weights of
$\bm{s}^{t}$ to zeros, we normalize it to become a unit-length vector for the
next iteration: $\bm{s}^{t}\leftarrow\frac{\bm{s}^{t}}{\|\bm{s}^{t}\|}$. It
can be seen that by gradually decreasing $K^{t}$ from $N$ down to the required
number of selected images $K$ during iteration, we can effectively filter out
images that are not sufficiently novel in the SGD learning process. This also
suggests that when we sample $\mathcal{J}_{i}^{t}$ for image $\bm{x}_{i}$ at
step $t$, we only need to draw those images with a non-zero weight (i.e.,
$\bm{s}_{j}^{t-1}>0$) as those with a zero weight do not contribute to the
update (9).

Our algorithm as described in (1) has a few meta parameters to specify for
each problem instance. The foremost one is $K$ \- the number of novel images
to be found and returned. A feature extractor $\phi(\bm{x})$ is also needed to
produce dense embeddings of images for perceptual distance calculation. The
SGD learning rate $\alpha$ controls the step size of each update, the momentum
parameter $\beta$ stabilizes the learning process, and
$J=|\mathcal{J}_{i}^{t}|$ is the batch size that determines the number of
random samples for gradient estimation. Following the practice of deep
learning, we also introduce the epoch size $E$ as the maximal number of
learning steps to be executed. Given a fixed $E$, we can schedule the sparsity
$K^{t}$ to gradually decrease from $N$ to $K$ through the learning process.
One simple heuristic is $K^{t}=\lfloor N-\frac{N-K}{E}t\rfloor$ for $1\leq
t\leq E$.

###  2.3 Complexity Analysis

It is clear from the pseudocode (1) that our algorithm does not need to store
the $N\times N$ distance matrix for novel image discovery. When $N$ is large
(e.g., over 10,000), the run-time space saving is huge compared with any
algorithms that require a full realization of the matrix (or the underlying
graph structure) to work. The time complexity of our algorithm is
$\Theta(E(NJM+N\log N))$, where $E$ is the epoch size, $J$ is the batch size,
and $M$ is the feature dimension. To sense check this, the time complexity of
the classical power iteration is $\Theta((\frac{M}{2}+L)N^{2})$, where $L$ is
the required number of steps towards convergence. It can be seen that our
algorithm scales well with the growth of $N$.

##  3 Results and Discussion

We put our algorithm to the test of synthetic datasets in Section 3.1 and
demonstrate its novelty detection accuracy by identifying images of novel
classes in Section 3.2, while analysing its gradient approximation behaviour
in these cases. We also compare it against the state-of-the-art methods on
increasingly larger datasets in Section LABEL:subsec:comparison. Without
mentioned otherwise, we run all experiments using the default parameters as
stated in the pseudocode (1).

(a) a

(b) b

(c) c

(d) d

Figure 1: A 2-dimensional synthetic dataset. The dataset contains 100
uniformly generated 2D points with a few of them perturbed to become
anomalies. The colour of each data point is from cold to warm indicating the
corresponding value of the sparse eigenvector computed by our algorithm (1).
The sparsity parameter $K$ is set to 9 (the exact number of anomaly points).
The epoch size $E$ is displayed under each sub figure from (a) to (h). The
remaining parameters are set to their default values.

(a) e

(b) f

(c) g

(d) h

Figure 2: A 50-dimensional synthetic dataset. This dataset has 1,000 points
generated from a mixture of two 50-dimensional Gaussian distributions (we only
plot the first two features for visualization). The colour of each data point
is from cold to warm indicating the corresponding value of the sparse
eigenvector computed by our algorithm (1). The sparsity parameter $K$ is set
to 120 (the exact number of anomaly points). The epoch size $E$ is displayed
under each sub figure from (a) to (h). The remaining parameters are set to
their default values.

###  3.1 Detecting Novelties in Synthetic Data

We first test the reaction of our algorithm to novel data points in both a low
and a high dimensional setting in Figure 1 and 2 respectively. For the low-
dimensional case we generate 100 uniformly distributed 2D points and make some
of them anomalies by adding random noises to their coordinates. The high-
dimensional dataset contains 1,000 points drawn from a mixture of two 50D
Gaussian distributions, and the anomalies are synthesized in the same way as
for the 2D case. It is interesting to see that our algorithm can locate a
majority of deviating data points in both datasets after just one gradient
step (9). For the 2D dataset, the third step manages to find all ground-truth
anomalies while leaving one false positive, which is then corrected by the
following step. Similar iterative corrections can be observed for the 50D
dataset, with a default number of 20 steps successfully identifying over 95%
of true anomalies. The default parameters of our algorithm work well in both
settings, leaving just the sparsity value $K$ to choose by a user.

(a) a

(b) b

Figure 3: Tiny-ImageNet train set (500 images from the class Egyptian
cat-n02124075 and 1 randomly chosen image from each of the remaining 119
classes) results. (a) shows the accuracy of our computed sparse eigenvector
used for differentiating normal and novel images with respect to each value of
$K$ \- an image is classified to be novel if its computed weight is non-zero.
(b) shows the weight of each image along with its t-SNE embeddings derived
from the MobileNet-V2 extracted 1,000-dimensional features when $K$ is set to
119 as the ground-truth number of anomalies.

(a) a

(b) b

(c) c

(d) d

(e) e

(f) f

(g) g

(h) h

Figure 4: Tiny-ImageNet test set (10,000 images) results. (a) and (b)
visualize our computed novelty weight for each image along with its PCA and
t-SNE embeddings derived from the MobileNet-V2 extracted 1,000-dimensional
features. (c)-(h) show the 6 images with the highest weight (displayed under
each corresponding sub figure) found by our algorithm with the sparsity
$K=100$.

###  3.2 Identifying Novel Images from Real Datasets

We now apply our algorithm to the highly non-trivial Tiny-ImageNet dataset
[32] and show the results in Figure 3 and 4. We use the popular lightweight
MobileNet-V2 [33] pre-trained neural network to extract a 1,000-dimensional
feature vector for each image. The first experiment we do is mixing the 500
training images of one particular category with 1 randomly drawn image from
each of the remaining 119 classes that serve as true novelties. By varying the
sparsity $K$ from 1 to 250 while keeping the other parameters as default, we
can run our algorithm to produce a series of anomaly predictions and measure
their accuracy using the class-balanced precision, recall, and f-measure
metrics. Figure 3a shows that the precision of our algorithm for novel class
identification remains high under a varying number of images found. The recall
keeps increasing when $K$ is increased from 1 to around 150 and then drops
slightly after that, which shows that our algorithm works well when the
sparsity is set just a bit higher than the actual number of novel images in
the dataset. The f-measure curve also reflects the overall monotonic
performance increase of our algorithm when $K$ is not set too high. The local
fluctuations of these metrics are due to the non-convexity of the problem (5)
and the inherent randomness of SGD 9, which can be simply smoothed out by
averaging multiple runs with different seeds. Figure 3b visualizes our
computed sparse eigenvector along with a t-SNE embedding [34] of the dataset.
It is interesting to note that the identified image anomalies concentrate
around the periphery as well as low-density regions of the empirical data
distribution in the feature space. To reproduce this on a held-out test set,
we also run our algorithm on the 10,000 Tiny-ImageNet test images (with no
labels available) with $K=100$ and show the results in Figure 4. While
principal component analysis (PCA) does not seem able to delineate the dataset
structures as t-SNE do, both visualizations show that our algorithm can
accurately identify images of low probability from the given collection, with
some examples displayed in Figure 4c to 4h.

##  4 Conclusion

We have proposed a simple yet effective algorithm for finding images of
distinct characteristics within a large collection. The key of our algorithm
is formulating the problem as solving for the $K$-densest subgraph within a
perceptual distance-weighted complete graph of images. While the original
problem is NP-hard, we have successfully relaxed it to be a sparse eigenvector
problem, which we have effectively solved using stochastic gradient descent
with sparsity clipping. Unlike any supervised learning methods, our algorithm
is training-free and therefore does not require human-annotated labels to
work. It also scales well with the number of the images in a collection in
terms of both space and time complexity. The experimental results show the
accuracy of our algorithm for picking up novel images.

An important future work would be establishing a theoretical analysis of the
convergence properties of our algorithm (1). The SGD method has previously
been applied to PCA problems [35] with certain convergence guarantee of eigen-
decomposition [36]. However, the analysis cannot extend to our case as we work
with distance rather than covariance matrices.

## References

  * [1] X. Yang, X. Qian, T. Mei, Learning salient visual word for scalable mobile image retrieval, Pattern Recognition 48 (10) (2015) 3093–3101.
  * [2] A. Clark, A nice surprise? predictive processing and the active pursuit of novelty, Phenomenology and the Cognitive Sciences 17 (3) (2018) 521–534.
  * [3] Y. Li, Y. Xie, Is a picture worth a thousand words? an empirical study of image content and social media engagement, Journal of Marketing Research 57 (1) (2020) 1–19.
  * [4] M. Zhang, X. Zhe, S. Chen, H. Yan, Deep center-based dual-constrained hashing for discriminative face image retrieval, Pattern Recognition 117 (2021) 107976\.
  * [5] X. Zhang, M. Shen, X. Li, F. Feng, A deformable cnn-based triplet model for fine-grained sketch-based image retrieval, Pattern Recognition (2021) 108508.
  * [6] P. Isola, J. Xiao, A. Torralba, A. Oliva, What makes an image memorable?, in: CVPR 2011, IEEE, 2011, pp. 145–152.
  * [7] A. Jaegle, V. Mehrpour, Y. Mohsenzadeh, T. Meyer, A. Oliva, N. Rust, Population response magnitude variation in inferotemporal cortex predicts image memorability, Elife 8 (2019) e47596.
  * [8] J. Kim, S. Yoon, V. Pavlovic, Relative spatial features for image memorability, in: Proceedings of the 21st ACM international conference on Multimedia, 2013, pp. 761–764.
  * [9] A. Khosla, A. S. Raju, A. Torralba, A. Oliva, Understanding and predicting image memorability at a large scale, in: Proceedings of the IEEE international conference on computer vision, 2015, pp. 2390–2398.
  * [10] P. Jing, Y. Su, L. Nie, H. Gu, J. Liu, M. Wang, A framework of joint low-rank and sparse regression for image memorability prediction, IEEE Transactions on Circuits and Systems for Video Technology 29 (5) (2018) 1296–1309.
  * [11] H. Squalli-Houssaini, N. Q. Duong, M. Gwenaëlle, C.-H. Demarty, Deep learning for predicting image memorability, in: 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2018, pp. 2371–2375.
  * [12] E. Akagunduz, A. G. Bors, K. K. Evans, Defining image memorability using the visual memory schema, IEEE transactions on pattern analysis and machine intelligence 42 (9) (2019) 2165–2178.
  * [13] S. Perera, A. Tal, L. Zelnik-Manor, Is image memorability prediction solved?, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2019, pp. 0–0.
  * [14] S. Basavaraju, A. Sur, Image memorability prediction using depth and motion cues, IEEE Transactions on Computational Social Systems 7 (3) (2020) 600–609.
  * [15] L. Itti, C. Koch, E. Niebur, A model of saliency-based visual attention for rapid scene analysis, IEEE Transactions on pattern analysis and machine intelligence 20 (11) (1998) 1254–1259.
  * [16] X. Hou, L. Zhang, Saliency detection: A spectral residual approach, in: 2007 IEEE Conference on computer vision and pattern recognition, Ieee, 2007, pp. 1–8.
  * [17] Y. Yan, J. Ren, G. Sun, H. Zhao, J. Han, X. Li, S. Marshall, J. Zhan, Unsupervised image saliency detection with gestalt-laws guided optimization and visual attention based refinement, Pattern Recognition 79 (2018) 65–78.
  * [18] W. Wang, Q. Lai, H. Fu, J. Shen, H. Ling, R. Yang, Salient object detection in the deep learning era: An in-depth survey, IEEE Transactions on Pattern Analysis and Machine Intelligence (2021).
  * [19] Y. Lv, B. Liu, J. Zhang, Y. Dai, A. Li, T. Zhang, Semi-supervised active salient object detection, Pattern Recognition 123 (2022) 108364.
  * [20] D. Zhang, J. Han, C. Li, J. Wang, Co-saliency detection via looking deep and wide, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 2994–3002.
  * [21] D. Zhang, H. Fu, J. Han, A. Borji, X. Li, A review of co-saliency detection algorithms: Fundamentals, applications, and challenges, ACM Transactions on Intelligent Systems and Technology (TIST) 9 (4) (2018) 1–31.
  * [22] C. Wang, Z.-J. Zha, D. Liu, H. Xie, Robust deep co-saliency detection with group semantic, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 2019, pp. 8917–8924.
  * [23] M. Mancas, O. Le Meur, Memorability of natural scenes: The role of attention, in: 2013 IEEE International Conference on Image Processing, IEEE, 2013, pp. 196–200.
  * [24] U. Feige, D. Peleg, G. Kortsarz, The dense k-subgraph problem, Algorithmica 29 (3) (2001) 410–421.
  * [25] C. Tsourakakis, F. Bonchi, A. Gionis, F. Gullo, M. Tsiarli, Denser than the densest subgraph: extracting optimal quasi-cliques with quality guarantees, in: Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, 2013, pp. 104–112.
  * [26] C. Tsourakakis, The k-clique densest subgraph problem, in: Proceedings of the 24th international conference on world wide web, 2015, pp. 1122–1132.
  * [27] A. Jazayeri, C. Yang, Frequent subgraph mining algorithms in static and temporal graph-transaction settings: A survey, IEEE Transactions on Big Data (2021).
  * [28] L. Bottou, Large-scale machine learning with stochastic gradient descent, in: Proceedings of COMPSTAT’2010, Springer, 2010, pp. 177–186.
  * [29] L. Qin, R.-H. Li, L. Chang, C. Zhang, Locally densest subgraph discovery, in: Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2015, pp. 965–974.
  * [30] C. Chekuri, K. Quanrud, M. R. Torres, Densest subgraph: Supermodularity, iterative peeling, and flow, in: Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2022, pp. 1531–1555.
  * [31] S. U. Pillai, T. Suel, S. Cha, The perron-frobenius theorem: some of its applications, IEEE Signal Processing Magazine 22 (2) (2005) 62–75.
  * [32] Y. Le, X. Yang, Tiny imagenet visual recognition challenge, CS 231N 7 (7) (2015) 3.
  * [33] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, L.-C. Chen, Mobilenetv2: Inverted residuals and linear bottlenecks, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 4510–4520.
  * [34] L. Van der Maaten, G. Hinton, Visualizing data using t-sne., Journal of machine learning research 9 (11) (2008).
  * [35] E. Oja, J. Karhunen, On stochastic approximation of the eigenvectors and eigenvalues of the expectation of a random matrix, Journal of mathematical analysis and applications 106 (1) (1985) 69–84.
  * [36] O. Shamir, Convergence of stochastic gradient descent for pca, in: International Conference on Machine Learning, PMLR, 2016, pp. 257–265.


# Spectral invariants and playing hide-and-seek on surfaces

Peter G. Doyle  Jean Steiner

(Original dated 27 September 2005
minor revisions dated 2009
No Copyright††thanks:  The authors hereby waive all copyright and related or
neighboring rights to this work, and dedicate it to the public domain. This
applies worldwide. )

###### Abstract

We prove the expected duration of a game of hide-and-seek played on a
Riemannian manifold under the laws of Brownian Motion is a spectral invariant:
it is a zeta-regularized version of the ‘trace’ of the Laplacian. An analogous
hide-and-seek game may be played on Markov chains, where the spectral
invariant that emerges is a classical quantity known as Kemeny’s constant. We
develop the analogies between the two settings in order to highlight the
connections between the regularized trace and Kemeny’s constant. Our proof
relies on the connections between Green’s functions and expected hitting
times, and the fact that the regularized trace may be approached via the
Green’s function.

AMS 2000 Mathematics Subject Classification  58J65, 60J05, 58J50

_Why bother about journals? The arXiv is where math research lives._ _— Oded
Schramm, 12 April 2008_

##  1 Introduction

The (regularized) trace of the Laplacian is a geometric spectral invariant
which has been shown to have interesting extremal properties. In [14],
Morpurgo showed that among all metrics conformally equivalent to the standard
round metric on the sphere, the regularized trace is minimized at the standard
round metric. More recently, in [17], Okikiolu has shown that the infimum of
the trace within the conformal class of any metric with a fixed volume is
bounded above by the trace of the standard round sphere with the same fixed
volume. The trace is closely related to a sharp Sobolev inequality, as is laid
out in both [14] and [17], and in some sense these results suggest that the
trace ‘hears’ the sphericalness of the surface. Analogous statements hold on
higher dimensional manifolds for the trace of the Paneitz operator, a higher
order elliptic operator, and the connection to a Sobolev inequality persists
in the higher dimensional setting.

In this paper we will give a probabilistic interpretation of the trace. We
will describe two hide-and-seek games in two different settings. The first
setting is that of a finite state Markov chain, the operator of giving rise to
our spectral invariant will be $I-P$, where $P$ is the transition matrix for
the Markov chain, and the second setting is that of a surface, in which case
the analogous operator is the Laplacian, $\Delta$. In both settings, the
operators $(I-P)$ and $\Delta$ have discrete eigenvalues, and the operators
capture the difference between the value of the function at a point and a
weighted average at the neighbors, and in both cases we will need to consider
the appropriate Green’s functions and approximate inverses to $(I-P)$ and
$\Delta$, in order to analyze the hide-and-seek games. In both settings, we
will see that the expected duration of the first of our hide-and-seek games is
a spectral invariant, and the expected duration of the second game is the
density for the spectral invariant.

Now, for the remainder of the introduction, we will describe the two hide-and-
seek games, and then interpret a spectral invariant on surfaces and Markov
chains in terms of the expected durations of the games . In each of the hide-
and-seek games we consider, one player’s location is chosen deterministically
and the other is chosen at random, according to the equilibrium measure. In
order to play this game on a finite state, discrete time Markov chain with $n$
states, given by $\mathcal{S}=\\{1,2,\ldots,n\\}$, the seeker will begin at a
state $i\in\mathcal{S}$, and the hider will be at a state $j\in\mathcal{S}$.
On a surface, $M$, the seeker will begin at a point $x\in M$, and the hider
will be at a point $y\in M$ On a surface, the probability of hitting a
particular point is $0$, so the seeker can’t find $y$. However, the seeker
can, with probability one, hit a ball of radius $\epsilon$ around a point, so
the hide-and-seek games will end when the seeker, having started out at the
point $x$, comes within a distance of $\epsilon$ from the hider at $y$ (i.e.
the seeker must reach the epsilon ball, $B_{\epsilon}(y)$). The two games that
will be considered are as follows:

Game I

The seeker is deterministically located at point $x\in M$, or at state
$i\in\mathcal{S}$. The hider is at a random location, $y\in M$ or
$j\in\mathcal{S}$.

Game II

The seeker is randomly located at $x\in M$, or at state $i\in\mathcal{S}$ .
The hider is at deterministic location, $y\in M$ or at state
$j\in\mathcal{S}$.

A priori, the expected duration of both games should depend on the
deterministically chosen data (i.e. the starting point for Game I, and the
ending point for Game II), but we will see that Game I actually turns out to
be independent of the deterministic data.

In order to give a more precise statement of our main results (and to give a
bit of context for the interest in a spectral invariant such as the
regularized trace), let’s think spectrally about each of the two settings. In
the case of surfaces, the Laplacian, $\Delta$ has a discrete spectral
resolution with eigenvalues $\\{\lambda_{j}\\}_{j=0}^{\infty}$, and the
eigenvalues are unbounded, with the lowest eigenvalue of $0$ corresponding to
the constant functions, and asymptotic growth $\lambda_{j}\sim j$ for large
$j$. Through a regularization process, one can make sense of both the
determinant and the trace of $\Delta$. The determinant is a well-known
spectral invariant that picks out geometric features such as constant
curvature, since it is extremized at constant curvature as stated by Osgood,
Phillips and Sarnak, in [18]. While the story of the log-determinant on
surfaces, is reasonably well understood, investigations of the regularized
trace have been undertaken more recently. The behavior of the regularized
trace on the sphere is well-understood, but there is much to be understood
about its behavior on other surfaces. The regularized trace is defined below
in (11), but heuristically one may think of it as
$Tr^{\prime}(\Delta^{-1})=\sum_{j\geq 1}^{\prime}\frac{1}{\lambda_{j}}$, where
we point out that $\Delta^{-1}$ only makes sense once one specifies its
behavior on the null-space, and the prime reminds us that the sum does not
actually converge. In [14], Morpurgo considers the regularized trace and its
behavior under conformal change, and he finds that it is extremized on the
round sphere. In [20], Steiner gives a characterization of the density for the
regularized trace. In work of Chiu, in [3], the regularized trace is shown to
be linearly related to the log-determinant for constant curvature tori, though
the ‘trace’ is not identified as such. Additional work towards understanding
the behavior of the regularized trace appears in work of Steiner and Doyle in
[6], and also in forthcoming work by Okikiolu in [17]. Additionally the
density for the regularized trace in connection with a condition for
criticality of the log-determinant in work of Okikiolu [16] and Richardson
[19].

In this paper, our discussion of spectral invariants on surfaces will be
complemented by a parallel discussion which takes place on a finite state,
discrete time Markov chain. If the Markov chain is irreducible and ergodic,
the transition matrix $P=(p_{ij})$ will have eigenvalues
$1>\nu_{1}\geq\nu_{2}\geq\cdots\nu_{n-1}$. Spectral invariants have also been
considered on Markov chains, see for example the survey of Aldous and Fill in
[1], and also work of Diaconis and Stroock in [4], Levene and Loizou in [11],
111There are many important references absent from this list.. The spectral
invariant in the present discussion is the ‘trace’ for the operator
$(I-P)^{-1}$, which will need to be suitably regularized, since the operator
$(I-P)$ has a one-dimensional null-space. It turns out that a quantity called
Kemeny’s constant, which is denoted by $K^{i}=K$ and defined below in (21) is
exactly this regularized trace. The quantity $K^{i}$, given below in (21) is
an expression given in terms of expected first hitting times, and, at first
glance, it appears that $K^{i}$ depends on $i$. However, $K^{i}$ was observed
by Kemeny to be constant in [10]. It turns out that Kemeny’s constant is
essentially the regularized trace of the inverse, i.e.
$Tr(I-P)^{-1}=\sum_{j=1}^{n-1}\frac{1}{1-\nu_{j}}$ and this perspective is put
forth by Meyer in [13]. Since this quantity appears in a number of places, we
point out that in his notes, [1] Aldous discusses the constancy of $K^{i}$,
which is referred to in the ‘Random target identity’, a phrase coined by
Lovász and Winkler in [12]. Doyle won a $ 50 prize for providing an
‘intuitively plausible’ explanation for why $K^{i}$ is constant in [5], and
Levene and Loizou [11] give a physical interpretation of Kemeny’s constant as
the expected length of time it will take a random (web) surfer to reach the
final destination web-page.

Given the suggestive notation we’re using here to describe these two
regularized traces on Markov chains and on surfaces, one might guess that
there should be a connection between the regularized trace on surfaces,
$Tr^{\prime}(\Delta^{-1})=\sum_{j\geq 1}^{\prime}\frac{1}{\lambda_{j}}$, and
Kemeny’s constant on Markov chains ,
$Tr^{\prime}(I-P^{-1})=\sum_{j=1}^{n-1}\frac{1}{1-\nu_{j}}$. Indeed, the hide-
and-seek games give us a way to related these two quantities, by considering
the the expected durations of the two games.

Our main result is the following theorem, which is a paraphrase of Theorem
4.2, and Corollaries 3.5 and 3.7 .

######  Theorem 1.1.

Consider the two hide-and-seek games described above on either a surface
$(M,g)$ with volume $1$, or on a Markov Chain with states $\mathcal{S}$ and
transition matrix $P$.

  1. 1.

In both settings, the expected duration for Game I, which starts
deterministically at point $x\in M$ or state $i\in\mathcal{S}$ , is
independent of the point $x\in M$ and $i\in\mathcal{S}$. Moreover, on a Markov
chain, the expected duration of the game is the spectral invariant Kemeny’s
constant, where $K=`Tr((I-P)^{-1})=\sum_{j=1}^{n-1}\frac{1}{1-\nu_{j}}$. On
surfaces, the expected duration of the game is essentially the regularized
trace, up to a universal correction depending on $\epsilon$.

  2. 2.

The expected duration for Game II, which ends deterministically at a point $y$
or a state $j$ is a density for the corresponding quantities in Game I, and in
general, the expected duration of Game II does depend on the deterministically
chosen data.

The density for the two regularized traces on surfaces and Markov chains may
be given as a mass-like term in the Green’s function, so Theorem 1.1 gives a
probabilistic interpretation of the regularized trace on surfaces, as well as
a probabilistic interpretation of the mass-like quantity in the Green’s
function. In order to prove Theorem 1.1, we will relate the expected hitting
times to the Green’s function for $\Delta$ on surfaces, and we will consider
the Green’s function for the operator $(I-P)$ on Markov chains.

#### Organization

In Section 2, the relevant background for Markov chains and surfaces is
detailed. Throughout our exposition, we emphasize the analogous statements on
Markov chains and on surfaces. In each case, we consider an operator that
gives the difference between the value of a function at a point, and a
weighted average value of the function on it’s neighbors, which is $(I-P)$ and
$\Delta$. We discuss the spectral resolutions and Green’s functions for the
relevant operators in each of the two settings. While we are largely recording
well-known background in Section 2, we believe that the carefully laid out
analogies do not appear explicitly in the literature, it seems worthwhile to
carefully develop the analogies between the two settings. In Section 3, we
analyze the expected duration of the hide-and-seek games on Markov chains, and
in Section 4 we consider the expected duration for surfaces. In each of these
two sections, we give a formulation of the expected first hitting times in
terms of Green’s functions, which we apply in computing the expected duration
of each of the hide-and-seek games. The main results are in Theorem 4.2, and
Corollaries 3.5 and 3.7.

##  2 Background

In this section we give the necessary background in order to set-up our hide-
and-seek games and their probabilistic interpretations on Markov chains and on
surfaces. We draw parallels between the analogous statements for the two
settings throughout this discussion.

###  2.1 Markov Chain background

A discrete, time-homogeneous Markov chain is a sequence of random variables,
$X_{0},X_{1},X_{2},\ldots$, valued in $S=\\{1,2,3,\ldots n\\}$. From the time-
homogeneity and Markov properties, we know that the Markov chain is
characterized by initial distribution, $p_{0}$, and transition matrix $P$ with
$(P)_{ij}=P(i,j)\geq 0$ where:

$(P)_{ij}=P(i,j)=\mathbf{P}(X_{k}=j|X_{0}=i_{0},X_{1}=i_{1},\ldots
X_{k-1}=i)=\mathbf{P}(X_{k}=j|X_{k-1}=i)$

In words, the $ij$ entry of the transition matrix, $P_{ij}$ gives the
probability that the Markov chain will go from state $i$ to state $j$ at any
given time step.

#### Difference and averaging operator, $(I-P)$

We will consider the operator $I-P$, which gives the difference between the
value of the function at a point, and the weighted average of the values at
the neighboring states (here, a state $j$ is a neighbor if $P_{ij}\neq 0$).
Consider the vector $\mathbf{f}=(f_{j})_{j=1}^{n}$ to be a function, and
compare the following equation to the analogous statement for the Laplacian on
surfaces, which appears in (9) in the next section.

$[(I-P)f]_{i}=f_{i}-\sum_{j}P_{ij}f_{j}$ (1)

In terms of electrical sources and potentials, the operator $(I-P)$ maps
potentials to source densities.222A careful treatment of the connections
between Markov Chains and electrical networks and potential theory is given by
Doyle and Snell in [7]

#### The spectrum of $(I-P)$

The matrix $P$ is row stochastic, so $\sum_{j=1}^{n}P_{ij}=1$ for all $1\leq
i\leq n$, and the constant vector $\mathbf{1}$ is a right eigenvector for $P$,
with eigenvalue $1$:

$P\mathbf{1}=\mathbf{1}$

From linear algebra, we know that there will also be a right eigenvector for
$P^{T}$ corresponding to the eigenvalue $1$, where $P^{T}w=w$, and taking the
transpose will give rise to $w^{T}P=w^{T}$.

######  Definition 2.1.

If the Markov chain is irreducible and ergodic, then the eigenvalue $1$ has
multiplicity $1$, and there is a unique positive left eigen-measure, $w$,
which is called the equilibrium measure and satisfies

$w^{T}P=w^{T},\hskip 10.0ptw_{i}>0,\hskip 10.0pt\sum_{i}w_{i}=1$

We recall two interpretations of the equilibrium measure that we will exploit
in the sequel:

######  Proposition 2.2 (See Grinstead and Snell, Ch. 11 [8] ).

For an irreducible, ergodic Markov chain $P$ with equilibrium measure $w^{T}$,
we have

  1. 1.

For large $n$, $w_{j}$ is the probability that $X_{n}=j$ (by the law of large
numbers), and.

$P^{n}\rightarrow W\text{ as }n\rightarrow\infty\text{ where }W=(W_{i}),\text{
with }W_{i}=(w^{T})$ (2)

  2. 2.

The expected first return time for $j$ is $\frac{1}{w_{j}}$.

Also, if $P$ is irreducible and ergodic, then $P$ has eigenvalues
$1>\nu_{1}\geq\nu_{2}\cdots\nu_{n-1}>-1$.

#### Spectrum of $(I-P)$.

Now, we will be interested in the operator $(I-P)$, which has eigenvalues
$0,1-\nu_{1},1-\nu_{2}\ldots 1-\nu_{n-1}$. For an irreducible ergodic Markov
chain, we know that the operator $(I-P)$ has a $1$-dimensional null-space, and
its range is $(n-1)$ dimensional:

$(I-P)\mathbf{1}=0;\hskip 30.0pt\text{ and }w^{t}(I-P)=0,\hskip 5.0pt\text{ so
}\forall v\in\mathbb{R}^{n},w^{t}(I-P)v=0\hskip 10.0pt\forall
v\in\mathbb{R}^{n}$ (3)

Thinking electrically, the fact that $\mathbf{1}$ is in the null-space says
that if we start with a potential that is identical at each state, there will
be no net current flowing anywhere at all, and so the source strength at each
state will be $0$. The second statement says that the source density
corresponding to any potential $v$ must have total source strength $0$, where
the source strength is measured with respect to the equilibrium measure $w$.

#### The Green’s function and inverting $(I-P)$

We are interested in inverting the operator $(I-P)$, but now we can see that
we can only hope to invert the operator on its range, namely the $(n-1)$
dimensional subspace, $\\{v\in\mathbb{R}^{n}:v\perp w\\}$. Since $\mathbf{1}$
is in the null-space for $(I-P)$, the inverse will only be unique up to a
choice of an additive constant. In order to specify the Green’s function,
which is given by the matrix $G=(G_{ij})$, we take a frame of $n$ source
densities, $S_{j}$ with $S_{j}\in\mathbb{R}^{n}\perp w$ corresponding to
having a source of strength $1$ at $j$, and a distributed sink of total
strength $1$. The source density, $S_{j}$, is the analog to the source density
on surfaces which is $\delta_{x}(y)-\frac{1}{V}$. On Markov chains, $S_{j}$
may be given in terms of the standard Euclidean basis vector $\mathbf{e}_{j}$
by:

$\mathbf{S}_{j}=\frac{1}{w_{j}}\mathbf{e_{j}}-\mathbf{1},\text{ which
satsifies }\sum w_{i}S_{j}(i)=w_{j}\frac{1}{w_{j}}-\sum w_{i}=0.$ (4)

If we can find the potential corresponding to each of the source densities,
$S_{j}$, then we can solve the equation $(I-P)u=f$ as long as $f\perp w$. For
each $j$, we would like to solve the following equation, which should be
compared with the distributional equation for the Green’s function on
surfaces, given below in (13):

$[(I-P)G]_{j}=\mathbf{S}_{j}=\frac{1}{w_{j}}\mathbf{e_{j}}-\mathbf{1}\text{ or
}(I-P)G=(\mathbf{S}_{1}\mathbf{S}_{2}\cdots\mathbf{S}_{n})$ (5)

Since the constants are in the null-space for $(I-P)$, the Green’s function
will only be unique up to a choice of an additive constant. We will demand
that $WG=0$, and Proposition 2.5 will demonstrate that our choice is
probabilistically natural. Moreover, this choice of normalization is the
analogous choice to taking $\int G(x,y)\mathrm{d}V(x)=0$ on surfaces, which
ensures that $\Delta^{-1}:1\rightarrow 0$. In the Markov chain setting, $G$,
the Green’s function exists, and the following Proposition gives a useful
concrete representation for the Green’s function.

######  Proposition 2.3.

[Meyer, [13]] Consider an irreducible and ergodic Markov chain with transition
matrix $P$, and (unique positive) equilibrium measure $w^{T}$, and let $W$
denote the rank one matrix with rows $W_{j}=w^{T}$, and consider the diagonal
matrix $D=Diag(\frac{1}{w_{i}})$ with entries $D_{ii}=\frac{1}{w_{i}}$, and
$\mathbf{1}$ denote the $n\times n$ matrix where each entry is a $1$. There is
a unique Green’s function, denoted by $G=(G_{ij})$ which satisfies the
following two requirements:

$\left\\{\begin{array}[]{ll}(I-P)G=D-\mathbf{1}&\text{ or }(I-P)G=(I-W)D\\\
WG=0&\end{array}\right.$ (6)

The Green’s function is given by:

$G=\left(\left(I-(P-W)\right)^{-1}-W\right)D$ (7)

######  Remark 2.4.

The matrix $D$ may be thought of as playing the same role as the inverse
metric plays in Riemannian geometry. The inverse metric, raises an index by
turning a covector into a vector. In the Markov chain setting, a co-vector is
a (possibly signed) measure, and by acting on the right, the matrix $D$ turns
a measure into a potential, which is a vector. The following interpretation of
the measure was also given by Meyer in [13].

######  Proposition 2.5.

[Compare Theorem 3.2 of Meyer, [13]] The normalization that we have chosen for
the Green’s function, given by $WG=0$, gives rise to a signed measure obtained
from the Green’s function:

$Z=GD^{-1}=\left(\left(I-(P-W)\right)^{-1}-W\right)$ (8)

The expected number of excess visits to a state $j$, for a Markov chain with
an initial distribution of $p_{0}=\mathbf{e}_{i}$, versus a Markov chain with
an initial probability distribution given by $p_{o}=w$ is given by $Z_{ij}$.

######  Remark 2.6.

Our notation for $Z=\left(I-(P-W)\right)^{-1}-W$ as the fundamental matrix is
consistent with Aldous’ notation in [1], however, Grinstead and Snell in [8]
use a different normalization and their fundamental matrix is just
$\left(I-(P-W)\right)^{-1}$.

######  Remark 2.7.

Speaking informally, if you spend a little extra time in a state $j_{1}$, you
won’t be able to spend quite as much time in the other states, so the total
number of excess visits to all of the states must be $0$. Therefore, the
probabilistic interpretation of Proposition 2.5 makes it clear that
$\sum_{j}Z_{ij}=0$.

###### Proof.

In order to verify the probabilistic interpretation of $Z_{ij}$, we need to
consider two Markov processes, both governed by the original transition matrix
$P$, but having two different initial distributions: $\mathcal{P}_{i\bullet}$
corresponds to having the initial distribution $p_{o}=\mathbf{e}_{i}$, which
means the chain is started at state $i$ with probability $1$, and
$\mathcal{P}_{w\bullet}$ corresponds to having the initial distribution
$p_{o}=w$, which means the chain is started at random with respect to the
equilibrium measure $w$. So, $\mathcal{P}_{ij}=P_{ij}$, which is the
probability that the chain goes from state $i$ to $j$ in one step, and since
$(w^{T}P)_{j}=w_{j}$, we know that $\mathcal{P}_{wj}=w_{j}$, which is the
probabillity that the state will be in state $j$ after one step, given that it
started at a random state chosen with respect to the equilibrium measure.

Now, we introduce a $0-1$ random variable $\chi^{k}_{\alpha j}$, where
$\alpha$ is the initial distribution of the Markov chain, and

$\chi^{k}_{\alpha j}=\left\\{\begin{array}[]{ll}1&X_{k}=j\text{ given
}X_{0}\text{ determined by }p_{0}=\alpha\\\ 0&\text{ else }\end{array}\right.$

Now, since $P_{ij}^{k}$ is the probability that the Markov chain is in state
$j$ at step $k$, given that it started in state $i$ (i.e. given that
$\alpha=\mathbf{e}_{i}$), we know that $\chi^{k}_{\mathbf{e}_{i}j}=1$ with
probability $(P^{k})_{ij}$. Since $w^{T}P^{k}=w^{T}$, we know that
$\chi^{k}_{wj}=1$ with probability $w_{j}$. The total number of visits within
the first $N$ steps for a process with initial distribution $\alpha$ would be
$\sum_{k=1}^{N}\chi^{k}_{\alpha j}$, so the expected excess number of visits
to a site $j$ within the first $N$ steps for the two processes is given by:

$\mathbb{E}[\sum_{k=0}^{N}\chi^{k}_{\mathbf{e}_{i}j}-\sum_{k=0}^{N}\chi^{k}_{wj}]=1\cdot\delta_{ij}+1\cdot
P_{ij}+1\cdot P^{2}_{ij}+\cdots 1\cdot P^{N}_{ij}-\stackrel{{\scriptstyle
n}}{{\overbrace{\left(1\cdot w_{j}+1\cdot w_{j}+\ldots+1\cdot
w_{j}\right)}}}\\\ =\delta_{ij}+P_{ij}+P^{2}_{ij}+\cdots
P^{N}_{ij}-(N+1)w_{j}$

Now, since this expression holds for any $j$, we may write it in terms of
matrices, and since the sum is finite, we may rearrange the sum as follows:

$\mathbb{E}[\sum_{k=0}^{N}\chi^{k}_{\mathbf{e}_{i}}-\sum_{k=0}^{N}\chi^{k}_{w}]=[I-W+(P-W)+(P^{2}-W)+\cdots(P^{N}-W)]_{ij}$

Now, by using the fact that $w$ is the stationary measure for $P$ and the
related fact that $PW=WP=W$, if coefficients are collected, one can show that
$(P^{k}-W)=(P-W)^{k}$. Now, we may re-write our expression for the expected
excess visits within the first $N$ steps, and we can recognize a truncation of
the Neumann series for the fundamental matrix, $Z$, defined in (8):

$\displaystyle\mathbb{E}[\sum_{k=0}^{N}\chi^{k}_{\mathbf{e}_{i}}-\sum_{k=0}^{N}\chi^{k}_{w}]$
$\displaystyle=$ $\displaystyle[I-W+(P-W)+(P^{2}-W)+\cdots(P^{N}-W)]_{ij}$
$\displaystyle=$ $\displaystyle[I-W+(P-W)+(P-W)^{2}+\cdots(P-W)^{N}]_{ij}$
$\displaystyle=$
$\displaystyle\left[\left(\sum_{l=0}^{N}(P-W)^{k}\right)-W\right]_{ij}$

Of course, in order to compute the expected number of excess visits between
the two processes, we must sum over all $l$, and so we have

$\mathbb{E}[\text{ Excess visits to }j\text{ for }\mathcal{P}_{ij}\text{
versus
}\mathcal{P}_{wj}]=\left[\left(\sum_{l=0}^{\infty}(P-W)^{k}\right)-W\right]_{ij}=Z_{ij}$

In case the reader is worrying about issues of convergence for this infinite
series, we should observe here that the Law of Large Numbers for irreducible
and ergodic Markov chains implies that $P^{n}\rightarrow W$, and since the
other eigenvalues for $P$ will be less than $1$, we can estimate
$|P^{n}-W|\leq|\lambda_{2}|^{n}$, where $|\lambda_{2}|<1$, and so the Neumann
series for $Z$ converges. ∎

###  2.2 Surface background

The second story presented in this paper takes place on a compact surface
without boundary $(M,g)$, so we now introduce the same cast of characters that
was just introduced in the context of Markov chains. In order to set-up the
notation, let $\mathrm{d}V$ denote the volume element, $\nabla$ denote the
gradient, and $\Delta=-div\cdot\nabla$ denotes the positive Laplacian operator
corresponding to $g$. If $(M,g)$ has volume $V$, then the corresponding
probability measure will be $\frac{\mathrm{d}V}{V}$.

#### Difference and averaging operator, $\Delta$

The operator $\Delta$ is analogous to the operator $(I-P)$ on Markov chains,
since both operators give the difference between the value of the function at
a point and the weighted average over neighboring points. In the surface
setting, we recall that the Laplacian operator gives the $\epsilon^{2}$
difference between a function $f$ at a point $x_{0}$, and the average of the
neighbors on a sphere of radius $\epsilon$ by considering Taylor’s theorem and
writing the metric in normal coordinates at the point $x_{0}$. This statement
should be familiar from the behavior in $\mathbb{R}^{n}$, where the most
familiar evidence of this averaging behavior is the well-known ‘Mean Value
Property’ for a harmonic function: if $\Delta f=0$ then the value of the
function at a point is equal to the average value on a sphere of radius
$\epsilon$. Note that we can compare the surface formulation of this averaging
property to that of $(I-P)$, which was stated in (1) :

$\Delta f(x_{0})\epsilon^{2}\approx
f(x_{0})-\frac{1}{Vol(S_{\epsilon}^{1})}\int_{S_{\epsilon}^{1}(x_{0})}f(x)\mathrm{d}S\text{
for }f\in C^{\infty}(M)$ (9)

#### Spectrum of $\Delta$ and the regularized trace

Since the Laplacian is a self-adjoint elliptic operator of order $2$, it has
discrete eigenvalues $\\{\lambda_{j}\\}_{j=0}^{\infty}$ that are bounded
below, with $0<\lambda_{1}\leq\lambda_{2}\cdots$, and Weyl’s law describes the
growth of the eigenvalues:

$\text{ Weyl's Law: }\hskip
5.0pt\lim_{j\rightarrow\infty}\frac{\lambda_{j}}{j}=\frac{4\pi}{Vol(M)}$ (10)

Now, in order to make connections with a spectral invariant precise, we
introduce the operator $\Delta^{-s}$, which is a linear operator on
$L^{2}(M,g)$. A branch of the logarithm is fixed so that
$\lambda^{s}=e^{s\log\lambda}$, with $1^{s}=1$. In order to define the
regularized trace, we introduce the operator $\Delta^{-s}$. From Weyl’s law
and the functional calculus, we see that $\Delta^{-s}$ is a Hilbert-Schmidt
operator for $\Re s\geq 1$. If $\Re(s)>1$, then $\Delta^{-s}$ is trace-class,
and the spectral Zeta function is defined to be the trace of $\Delta^{-s}$:

$Z(s)\stackrel{{\scriptstyle
def}}{{=}}Trace(\Delta^{-s})=\sum_{j=1}^{\infty}\lambda_{j}^{-s}\hskip
30.0pt\Re s>1$

The zeta function has a simple pole at $s=1$, and it is well known that the
zeta function may be continued to a meromorphic function defined on the
complex plane. In [14], Morpurgo considers the following spectral invariant,
which is the regularized trace of $\Delta^{-1}$:

######  Definition 2.8.

Assume $(M,g)$ is a closed compact surface with volume normalized to be
$4\pi$, then the regularized trace of $\Delta$, or the regularized Zeta
function at 1 is defined as follows:

$\tilde{Z}_{g}(1)=\lim_{s\rightarrow 1}\left(Z(s)-\frac{1}{s-1}\right)$ (11)

The regularized trace may also be approached via the Green’s function for
$\Delta$, as will be stated in Proposition 2.10.

#### The Green’s function and $\Delta^{-1}$

As is the case for Markov chains, the Green’s function allows us to invert the
Laplacian, $\Delta$, a linear operator on the Sobolev space $H_{2}(M)\subset
L^{2}(M,g)$. The null-space and range are $\Delta$ are characterized by the
following analog to the Markov chain situation described in (3).

$\Delta 1=0;\hskip 30.0pt\int_{M}\Delta f\mathrm{d}V=0\forall f\in H_{2}(M)$
(12)

Again, the inverse can not be defined on all of $L^{2}(M,g)$, since the image
of $\Delta$ consists of the functions that are orthogonal to the constants.
Additionally, the inverse for $\Delta$ will only be unique up to the choice of
a constant. As was the case on Markov chains, the Green’s function, $G(x,y)$
is the potential corresponding to a source density that has a source strength
of $1$ at a point $x$ and a distributed sink that has a total strength of $1$,
which corresponds to solving the distributional equation (13). We will choose
the constant so that the $\int G(x,y)\mathrm{d}V(x)=0$, or in otherwords
$\Delta^{-1}$ takes constant functions to $0$. The following proposition
describes the Green’s function on a surface. While we cannnot write down the
Green’s function explicitly on a surface in the way that we wrote down the
Green’s function on a Markov chain in Proposition 2.3, the parametrics in (15)
is the next best thing. Since the surface behaves locally like a piece of the
plane, we expect to see the same logarithmic singularity that one sees for the
Green’s function on Euclidean space.

######  Proposition 2.9.

[See, for example, Aubin [2]] On a surface $(M,g)$ with volume $V$, there
exists a unique Green’s function, $G(x,y)$ which is the integral kernel for
the operator $\Delta^{-1}$, it is symmetric in $x$ and $y$, and satisfies the
following two equations:

$\displaystyle\Delta_{y}^{distr}G(x,y)\stackrel{{\scriptstyle\mathrm{d}V}}{{=}}\delta_{x}(y)-\frac{1}{Vol(M)}$
(13)
$\displaystyle\int_{M}G(x,y)\mathrm{d}V(y)=\int_{M}G(x,y)\mathrm{d}V(x)=0$
(14)

The Green’s function is smooth away from the diagonal, and for small geodesic
distance, $d(x,y)$, there is a logarithmic singularity, with:

$G(x,y)=\frac{-1}{2\pi}\log(d(x,y))+m(x)+O(d(x,y))$ (15)

The term $m(x)=\lim_{y\rightarrow x}[G(x,y)+\frac{1}{2\pi}\log d(x,y)]$ is
called the Robin’s mass.

#### Comparing the two methods of regularizing the trace for $\Delta$

In the context of this paper, we will be considering expected hitting times,
which means we will be working with the Green’s function, so it is worth
recording the fact that the regularized trace may either be defined through
the zeta function as in (11), or by integrating the Robin’s mass, which
appears in (15). The following proposition gives the precise relation between
the two quantities:

######  Proposition 2.10.

[Morpurgo [15], S. [20]] Assume $(M,g)$ is a closed compact surface with
volume $V$. Then, the Robin’s mass $m(x)$ appearing in (15) is a density for
the spectral invariant $\tilde{Z}(1)$ defined in (11), and the regularized
trace may be written in terms of the integral of the Robin’s mass as follows:

$\tilde{Z}(1)=\int_{M}m(x)\mathrm{d}V(x)+\frac{V}{4\pi}(-2\log
2+2\gamma)\text{ with Euler gamma: }\gamma\approx.5772$ (16)

##  3 Hide-and-seek and Kemeny’s Constant on Markov Chains

#### Expected first hitting times via the Green’s function

In order to analyze the expected duration of the hide-and-seek games on a
Markov chain that were explained in Items Game I and Game II of the
Introduction, we need to know what the expected first hitting times are to get
from a state $i$ to a state $j$. After defining the hitting times, for
expository purposes, we will formulate the system of equations that the
hitting times must satisfy, we will then compute these hitting times in terms
of the Green’s function. A good reference for the relation between the
expected first hitting times and the Green’s function on Markov chains is the
Grinstead and Snell’s text on Probability, [8].

######  Definition 3.1.

Let $m_{ij}$ denote the expected first hitting time for a Markov process to
reach state $j$, given that it started in state $i$. In other words, if $T$ is
the random variable satisfying $T=\inf\\{n\geq 1:X_{n}=j\\}$, then

$m_{ij}\stackrel{{\scriptstyle def}}{{=}}\mathbb{E}[T|X_{0}=i],\text{ the
expected first hitting time. \hskip 10.0pt Take }m_{ii}=0$ (17)

Now, in order to see the relation between the Green’s function and the first
hitting times, we derive (18), which is satisfied by the expected hitting time
$m_{ij}$. We use a ‘first-step’ analysis, i.e. we will exploit the adage that
states, ‘A journey of a thousand miles begins with a single step’: if you are
wandering on the Markov chain at random and you find yourself at a state $i$,
then, you will step to any state $k$ with probability $p_{ik}$, and the
expected length of the remaining journey to state $j$ will be the sum of that
first step and the expected time it will take you to then get from $k$ to $j$,
i.e.

$m_{ij}=1+\sum_{k}p_{ik}m_{kj}$

Consequently, the matrix of expected hitting times must satisfy:

$\displaystyle[(I-P)M]_{ij}$ $\displaystyle=$ $\displaystyle 1\text{ for
}j\neq i$ (18) $\displaystyle m_{ii}$ $\displaystyle=$ $\displaystyle 0$

By a similar first-step analysis for a return trip to $i$, one sees that we
must have $1+\sum_{k=1}^{n}p_{ik}m_{ki}=\frac{1}{w_{i}}$, since if we take one
step from state $i$, we will go to any other state $k$ with probability $k$,
and then we expect the remainder of the journey to take $m_{ki}$. Of course,
from Proposition 2.2, we know that the expected return time to $i$ is exactly
$\frac{1}{w_{i}}$. Consolidating all of the requirements on $m_{ij}$, we find
that the expected hitting times must solve the following equations:

$\displaystyle(I-P)M$ $\displaystyle=$
$\displaystyle\mathbf{\hat{1}}-(\frac{1}{w_{i}}e_{i})=\mathbf{\hat{1}}-D$ (19)
$\displaystyle m_{ii}$ $\displaystyle=$ $\displaystyle 0$

Now, these equations should look familiar to us from the discussion of the
Green’s function in Section 2.

######  Lemma 3.2.

[Grinstead and Snell [8]] For an irreducible ergodic Markov chain where $G$
denotes the Green’s function defined in Proposition 2.3, and $M$ denotes the
matrix of expected hitting times defined in 3.1, then $M$ may be represented
in terms of the Green’s function:

$m_{ij}=G_{jj}-G_{ij}$ (20)

#### Kemeny’s constant and density

Now, we are finally in a position to analyze the expected duration of the
hide-and-seek games on a Markov chain.

#### Game I

Observe that for the first hide-and-seek game, in which a seeker begins at a
deterministically chosen state $i$ and must find a hider located at a randomly
chosen state $j$ chosen with probability $w_{j}$, the expected duration of the
game is given by the following quantity:

$K^{i}\stackrel{{\scriptstyle def}}{{=}}\sum_{j}m_{ij}w_{j}$ (21)

A priori, the expected duration of the hide-and-seek game should depend on the
starting point $i$, as indicated by the $i$ dependence in (21).

######  Theorem 3.3 (Kemeny, [10]).

The quantity $K_{i}$ is a constant $K$, independent of $i$.

We will refer to the constant quantity, $K$ defined in Equation 21 as Kemeny’s
constant.

One can prove this theorem by showing that $K_{i}$ is harmonic, and therefore
constant, as is done by Doyle in [5]. Or, one can prove the theorem explicitly
computing $K_{i}$ in terms of the Green’s function. By substituting the
expression for $m_{ij}$ that was given in terms of the Green’s function in
(20), and observing that one is left with a quantity that is independent of
$i$, namely, for the matrix $Z$ addressed in Proposition 2.5, which is the
signed measure corresponding to $G$, one finds:

$K=K^{i}=\sum_{j}Z_{jj}=Tr(Z)$

Of course, the trace of a matrix is a spectral invariant, which is the content
of the next theorem.

######  Theorem 3.4 (Levene-Loizou, [11], Aldous [1], Meyer [13]).

Kemeny’s constant is a spectral invariant. If the eigenvalues for $P$ are
given by $1,\nu_{1},\ldots\nu_{n-1}$, then,

$K=K^{i}=Tr(Z)=\sum_{j=1}^{n-1}\frac{1}{1-\nu_{j}}.$

Heuristically, Kemeny’s constant can be thought of as $Tr((I-P)^{-1})$, but of
course the honest approach is to use the expression for the Green’s function
given in (7) and observe that $K=Tr(Z)=Tr((I-(P-W))^{-1}-W)$, and $Z$ has
eigenvalues $0,1-\nu_{1},1-\nu_{2},\ldots 1-\nu_{n-1}$.

Now, we may give the following interpretation of the Game I, in which the
seeker’s state is chosen deterministically at a state $i$, and the hider’s
state is chosen at random.

######  Corollary 3.5.

On a Markov Chain, the expected duration of Game I, detailed in Item Game I of
the Introduction, is a spectral invariant that is independent of the starting
point $i$.

#### Game II

Now, in Game II, the seeker is located at a randomly chosen state, $i$ with
probability $w_{i}$, and a determinstically located hider must be found at
$j$. The expected duration of this game is given by:

$K_{j}=\sum_{i}w_{i}m_{ij}$ (22)

######  Lemma 3.6.

The quantity $K_{j}$ is a density for Kemeny’s constant, and it is the
diagonal term of $G$, the Green’s function defined in 2.3:

$K_{j}=G_{jj}$

Now, we may connect this density to our second hide-and-seek Game:

######  Corollary 3.7.

The expected duration of hide-and-seek Game II, detailed in Item Game II, is
given by the quantity $K_{j}$ of (22), and in general it will depend on the
point $j$. The expected duration of Game II is a density for Kemeny’s constant
(here the measure is the stationary measure, $w$), and thus $K_{j}$ is a
density for a spectral invariant.

##  4 Hide-and-seek and the regularized trace on Surfaces

Now, we would like to consider the same two hide-and-seek Games on a surface,
but we must begin by discussing a stochastic process. The seeker’s path will
follow the laws for Brownian motion.

#### Brownian Motion on manifolds and expected first hitting times via Green’s
functions

(See [9] for details concerning the construction of Brownian Motion on a
manifold). In order to compute the expected duration of the hide-and-seek
games, we must consider the behavior of a Brownian ‘seeker’. On a surface,
$(M,g)$, one can construct a Brownian motion $\\{\Sigma_{t}^{x}\\}_{t\geq 0}$,
which begins at a point $x\in M$, i.e. $\Sigma_{0}^{x}=x$. Since a closed
compact surface is stochastically complete, the Brownian motion will be
recurrent. While the probability of hitting any particular point $y$ is $0$
almost everywhere, the probability of hitting the boundary of an
$\epsilon$-ball around $y$ will be $1$. We let the random variable
$\tau^{x}(B_{\epsilon}(y))$ denote the first hitting time for a Brownian
motion beginning at $x$ to reach an $\epsilon$ ball around $y$, i.e.,

$\tau^{x}(B_{\epsilon}(y))=\inf\\{t>0:\Sigma_{t}^{x}\in\partial
B_{\epsilon}(y)\\},$

and we write the expected value of $\tau^{x}$ for a given epsilon ball as
$H_{\epsilon}(x,y)=\mathbb{E}[\tau^{x}(B_{\epsilon}(y)]$. Recall that $\Delta$
is the analg to the operator $(I-P)$ in the Markov chain setting, so the
quantities $H(x,y)$ and $m_{ij}$ are the first hitting times for surfaces and
Markov chains respectively, and it should come as no surprise that $H(x,y)$
solves a partial differential equation which is the continuous space analog of
(18). In particular, it is well-known that $H_{\epsilon}(x,y)$ solves the
following boundary value problem333Here we consider a re-scaled Brownian
motion to follow the geometer’s convention for the Laplacian:

$\left\\{\begin{array}[]{r@{\quad\quad\quad}l}\Delta H_{\epsilon}(x,y)=1&x\in
M\setminus B_{\epsilon}(y)\\\ H_{\epsilon}(x,y)=0&x\in\partial[M\setminus
B_{\epsilon}(y)]\end{array}\right.$ (23)

Now, the expected hitting time, $H_{\epsilon}(x,y)$ can be written in terms of
the Green’s function, as stated in the following lemma, which is the
continuous analog to Lemma 3.2 on Markov chains. Notice that on surfaces we
need to introduce a $\log\epsilon$ regularization, which reflects the fact
that it would take infinitely long to actually hit the point $y$.

######  Lemma 4.1.

On a surface $(M,g)$, if $H_{\epsilon}(x,y)$ is the expected value of the
first hitting time for a Brownian motion beginning at $x$ to reach
$B_{\epsilon}(y)$ satisfying (23), and $G(x,y)$ is the Green’s function for
the Laplacian characterized in Proposition 2.9, then $H_{\epsilon}(x,y)$ may
be written uniquely in terms of $G(x,y)$ and the robin’s mass $m(y)$ appearing
in (15) as follows:

$\left\\{\begin{array}[]{r@{\quad\quad\quad}l}H_{\epsilon}(x,y)=-VG(x,y)-\frac{V}{2\pi}\log\epsilon+Vm(y)+u_{\epsilon}(x)&x\in
M\setminus B_{\epsilon}(y)\\\
H_{\epsilon}(x,y)=0&x\in\bar{B_{\epsilon}(y)}\end{array}\right.$ (24)

Here, $u^{\epsilon}(x)$ is harmonic on $M\setminus{B_{\epsilon}(y)}$ and
$\|u^{\epsilon}\|_{C^{0}}=O(\epsilon)$.

###### Proof.

While one only needs to verify that $H_{\epsilon}(x,y)$ defined in Equation 24
satisfies (23), we believe it is more illuminating to prove the lemma by
demonstrating how to construct the correct function, $H_{\epsilon}(x,y)$. We
begin by observing that a good candidate for solving the hitting time problem
in (23) would be $-G(x,y)$, since we know from (13) that if $x\neq y$, we have
(for $V=Vol(M,g)$)

$\Delta_{x}^{distr}G(x,y)\stackrel{{\scriptstyle\mathrm{d}V}}{{=}}-\frac{1}{V}\text{
for }x\neq y$

Of course, away from $y$, the right-hand side of the equation is smooth, so it
holds point-wise for $x\neq y$, and we must re-scale the equation so that we
are working with the probability measure, which is given by
$\frac{\mathrm{d}V}{Vol(M)}$:

$\Delta_{y}\left(-Vol(M)G(x,y)\right)\stackrel{{\scriptstyle\frac{\mathrm{d}V}{Vol(M)}}}{{=}}1\text{
on }M\setminus\overline{B_{\epsilon}(y)}$

Now, while we have arranged for the Laplacian to behave properly on
$M\setminus\overline{B_{\epsilon}(y)}$, but of course, we need to normalize
our candidate it so that the hitting time is $0$ at the boundary, so we will
have $H^{\epsilon}_{y}(x)=-Vol(M)G(x,y)+(correction)$. From the parametrics of
the Green’s function given in (15), we know that for $x\in\partial
B_{\epsilon}(y)$, we have:

$-Vol(M)G(x,y)=\frac{V}{2\pi}\log
d(x,y)-Vm(y)+O(d(x,y))=\frac{V}{2\pi}\log\epsilon-Vm(y)+O(\epsilon)$

Now, by correcting $-Vol(M)G(x,y)$ with the logarithmic term and the Robin’s
mass term, we have the hitting time that solves the correct differential
equation on the interior of our region, and it will be correct up to a
function that is $O(\epsilon)$ on the boundary of the ball. In order to obtain
an exact solution we need to deal with the $O(\epsilon)$ correction term,
without affecting the differential equation on the interior. In particular, we
must find a function $u_{\epsilon}(x)$ satisfying:

$\displaystyle u_{\epsilon}(x)=-Vol(M)G(x,y)-(\frac{V}{2\pi}\log\epsilon-
Vm(y))$ $\displaystyle x\in\partial B_{\epsilon}(y)$ $\displaystyle\Delta
u_{\epsilon}(x)=0$ $\displaystyle x\in M\setminus\overline{B_{\epsilon}(y)}$

From standard PDE theory on a surface (see, for example, [21]), it is known
that we can solve this equation for the harmonic extension $u_{\epsilon}(x)$,
and $u_{\epsilon}(x)$ will be unique. From the maximum principle, we know that
$u_{\epsilon}$ will attain its maximum on the boundary, $\partial
B_{\epsilon}(y)$, so we know that

$\|u_{\epsilon}\|_{C^{0}(M\setminus\overline{B_{\epsilon}(y)}}=O(\epsilon)$

Now, since $\Delta u_{\epsilon}(x)=0$, we may add it to the modified Green’s
function without disrupting the differential equation on the interior, so we
now know that the function $H_{\epsilon}(x,y)$ defined in Equation 24 does
satisfy the Dirichlet problem for the expected hitting times. The function
$H_{\epsilon}(x,y)$ is unique and positive because $H_{\epsilon}(x,y)$ is sub-
harmonic and $H_{\epsilon}(x,y)=0$ on the boundary, so by the maximum
principle, since it is not constant, it must be positive on
$M\setminus\overline{B_{\epsilon}(y)}$. ∎

Now that we know the expected hitting time for an individual point, we may
consider the two hide-and-seek games.

#### Hide-and-seek on surfaces

We summarize the behavior of the two hide-and-seek games in the following
theorem.

######  Theorem 4.2.

On a surface $(M,g)$ with volume $V$ and Robin’s mass $m(x)$ defined in (15),
the two hide-and-seek games described in the Introduction exhibit the
following behavior on $M$.

  1. 1.

The expected duration for Game I, which starts deterministically at $x$, is

$\int_{M\setminus
B_{\epsilon}(x)}H_{\epsilon}(x,y)\mathrm{d}V(y)=-\frac{V}{2\pi}\log\epsilon+\int_{M}m(y)\mathrm{d}V(y)+O(\epsilon),$

and the expected duration for Game I is independent of the starting point,
$x$.

  2. 2.

The expected duration for Game II, which ends deterministically at $y$ is
$\int_{M\setminus
B_{\epsilon}(y)}H_{\epsilon}(x,y)\mathrm{d}V(x)=-\frac{V}{2\pi}\log\epsilon+Vm(y)+O(\epsilon)$.

###### Proof.

In order to verify the first statement, we must allow the point $y$ to be
randomly chosen with respect to the probability measure
$\frac{\mathrm{d}V(y)}{V}$, so we will integrate the $\epsilon$-regularized
hitting time, $H_{\epsilon}(x,y)$ given in Lemma 4.1 with respect to $y$. Note
that if according to the rules of our game, we have already found any point
$y\in B_{\epsilon}(x)$, so we need only integrate over $M\setminus
B_{\epsilon}(x)$:

$\mathbb{E}[\text{duration Game
I}]=\int\limits_{M}H^{\epsilon}(x,y)\frac{\mathrm{d}V(y)}{V}\\\
=\int\limits_{M\setminus
B(x,\epsilon)}H^{\epsilon}(x,y)\frac{\mathrm{d}V(y)}{V}=\int\limits_{M\setminus
B(x,\epsilon)}\left(-G(x,y)-\frac{V}{2\pi}\log\epsilon+Vm(y)+u_{\epsilon}(y)\right)\frac{\mathrm{d}V(y)}{V}$
(25)

Now, since we have normalized the Green’s function to integrate to $0$, the
contribution of the integral of $G(x,y)$ over $M\setminus B(x,\epsilon)$ will
be $O(\epsilon)$, and we can lump this contribution together with that of the
harmonic function $u_{\epsilon}(y)$, for which the $C^{0}$ norm is also
$O(\epsilon)$:

$\mathbb{E}[\text{ duration Game
I}]=O(\epsilon)-\frac{V}{2\pi}\log\epsilon+\int\limits_{M\setminus
B(x,\epsilon)}m(y)\mathrm{d}V(y)=-\frac{V}{2\pi}\log\epsilon+\int\limits_{M}m(y)\mathrm{d}V(y)+O(\epsilon)$

Note that this term is clearly constant with respect to $x$, since anything
involving $x$ has been integrated away (and we can estimate the $O(\epsilon)$
term uniformly).

Now, in order to obtain the expected duration of Game II, we do almost the
same calculation, however since the seeker is randomly chosen, we must
integrate with respect to $\frac{\mathrm{d}V(x)}{V}$.

$\displaystyle\int\limits_{M}H^{\epsilon}(x,y)\frac{\mathrm{d}V(x)}{V}$
$\displaystyle=$ $\displaystyle\int\limits_{M\setminus
B(y,\epsilon)}H^{\epsilon}(x,y)\frac{\mathrm{d}V(x)}{V}$ $\displaystyle=$
$\displaystyle\int\limits_{M\setminus
B(y,\epsilon)}\left(-G(x,y)-\frac{V}{2\pi}\log\epsilon+Vm(y)+u_{\epsilon}(y)\right)\frac{\mathrm{d}V(x)}{V}$
$\displaystyle=$ $\displaystyle-\frac{V}{2\pi}\log\epsilon+Vm(y)+O(\epsilon)$

∎

Now, recalling that the integral of the Robin’s mass is related to the
regularized trace on a surface, as stated in Proposition 2.10, so we see that
we have developed a probabilistic interpretation for the regularized trace and
its density by way of the two hide-and-seek games.

######  Corollary 4.3.

On surfaces, the expected duration of Game I is the regularized trace for
$\Delta^{-1}$, up to an $\epsilon$-regularization. The expected duration of
Game II on surfaces is the density for the $\epsilon$-regularized trace.

#### Application

Now, we can use the probabilistic intuition to consider how the regularized
trace for different surfaces might compare. At a strictly intuitive level, we
indicate how one might apply this to a family of rectangular tori with a flat
metric, where area is fixed, and the lengths are getting longer and longer. If
one thinks about the regularized trace in terms of the hide-and-seek game, it
seems quite plausible that as the rectangle is becoming longer (and therefore
skinnier since the area is fixed), there will be more and more area that is
farther away from a given point. Therefore, it will take longer to reach the
regions that are increasingly far away, so as the fundamental rectangles
become longer and skinnier, the regularized trace becomes larger. As
demonstrated in [6], by connecting the regularized trace with the Dedekind eta
function, one can verify that this behavior holds. Moreover, the probabilistic
intuition suggests that if one ‘blows a bubble’ into a long-skinny torus, than
most of the area will be closer to a given point, so the regularized trace
should be smaller, and in particular it should be close to that of the sphere.
This phenomenon does indeed occur, and it is described in [6].

## References

  * [1] David Aldous and Jim Fill.  Reversible Markov Chains and Random Walks on Graphs.  2002\.  http://www.stat.berkeley.edu/users/aldous/RWG/book.html.
  * [2] Thierry Aubin.  Some nonlinear problems in Riemannian geometry.  Springer-Verlag, Berlin, 1998.
  * [3] Patrick Chiu.  Height of flat tori.  Proc. Amer. Math. Soc., 125(3):723–730, 1997.
  * [4] Persi Diaconis and Daniel Stroock.  Geometric bounds for eigenvalues of Markov chains.  Ann. Appl. Probab., 1(1):36–61, 1991.
  * [5] Peter Doyle.  The Kemeny constant of a Markov chain.  2009, arXiv:0909.2636v1 [math.PR].
  * [6] Peter Doyle and Jean Steiner.  Blowing bubbles on the torus.  2005\.
  * [7] Peter G. Doyle and J. Laurie Snell.  Random Walks and Electric Networks.  The Mathematical Association of America, 1984, arXiv:math/0001057v1 [math.PR].
  * [8] Charles M. Grinstead and J. Laurie Snell.  Introduction to Probability: Second Revised Edition.  American Mathematical Society, Providence, Rhode Island, United States, 1997.
  * [9] Elton P. Hsu.  Stochastic analysis on manifolds, volume 38 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2002.
  * [10] John G. Kemeny and J. Laurie Snell.  Finite Markov chains.  The University Series in Undergraduate Mathematics. D. Van Nostrand Co., Inc., Princeton, N.J.-Toronto-London-New York, 1960.
  * [11] Mark Levene and George Loizou.  Kemeny’s constant and the random surfer.  Amer. Math. Monthly, 109(8):741–745, 2002.
  * [12] Laslo Lovasz and Peter Winkler.  Efficient stopping rules for Markov chains.  In STOC ’95: Proceedings of the twenty-seventh annual ACM symposium on Theory of computing, pages 76–82, New York, NY, USA, 1995. ACM Press.
  * [13] Carl D. Meyer, Jr.  The role of the group generalized inverse in the theory of finite Markov chains.  SIAM Rev., 17:443–464, 1975.
  * [14] C. Morpurgo.  The logarithmic Hardy-Littlewood-Sobolev inequality and extremals of zeta functions on ${S}^{n}$.  Geom. Funct. Anal., 6(1):146–171, 1996.
  * [15] Carlo Morpurgo.  Sharp inequalities for functional integrals and traces of conformally invariant operators.  Duke Math. J., 114(3):477–553, 2002.
  * [16] K. Okikiolu.  Critical metrics for the determinant of the Laplacian in odd dimensions.  Ann. of Math. (2), 153(2):471–531, 2001.
  * [17] K. Okikiolu.  Extremals for logarithmic Hardy-Littlewood-Sobolev inequalities on compact manifolds.  2006, arXiv:math/0603717v1 [math.SP].
  * [18] B. Osgood, R. Phillips, and P. Sarnak.  Extremals of determinants of Laplacians.  J. Funct. Anal., 80(1):148–211, 1988.
  * [19] Ken Richardson.  Critical points of the determinant of the Laplace operator.  J. Funct. Anal., 122(1):52–83, 1994.
  * [20] Jean Steiner.  A geometrical mass and its extremal properties for metrics on $S^{2}$.  Duke Math. J., 129(1):63–86, 2005.
  * [21] Michael E. Taylor.  Partial differential equations. I, volume 115 of Applied Mathematical Sciences.  Springer-Verlag, New York, 1996.  Basic theory.