# Online Detection of Golden Circuit Cutting Points
††thanks: This research was supported in part by NSF Awards 2216923, 2117439,
and 2230111.

Daniel T. Chen  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Ethan H. Hansen  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Xinpeng Li  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Aaron Orenstein  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Vinooth Kulkarni  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Vipin Chaudhary  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>Qiang Guan  Kent State University
Kent, OH, USA
<EMAIL_ADDRESS>Ji Liu  Argonne National Laboratory
Lemont, IL, USA
<EMAIL_ADDRESS>Yang Zhang  University of Illinois Urbana-Champaign
Champaign, IL, USA
<EMAIL_ADDRESS>Shuai Xu  Case Western Reserve University
Cleveland, OH, USA
<EMAIL_ADDRESS>
###### Abstract

Quantum circuit cutting has emerged as a promising method for simulating large
quantum circuits using a collection of small quantum machines. Running low-
qubit circuit “fragments” not only overcomes the size limitation of near-term
hardware, but it also increases the fidelity of the simulation. However,
reconstructing measurement statistics requires computational resources—both
classical and quantum—that grow exponentially with the number of cuts. In this
manuscript, we introduce the concept of a golden cutting point, which
identifies unnecessary basis components during reconstruction and avoids
related downstream computation. We propose a hypothesis-testing scheme for
identifying golden cutting points, and provide robustness results in the case
of the test failing with low probability. Lastly, we demonstrate the
applicability of our method on Qiskit’s Aer simulator and observe a reduced
wall time from identifying and avoiding obsolete measurements.

###### Index Terms:

quantum circuit cutting, circuit cutting, circuit knitting, circuit
reconstruction, hypothesis-testing, golden cutting point

##  I Introduction

Quantum circuit cutting refers to the method of splitting quantum circuits
into a set of small independent circuit fragments [1]. Using circuit cutting
methods, large quantum circuits can be simulated by a collection of smaller
machines, barring some addition classical computing resources. Moreover, it
was also empirically shown that cutting the circuit reduces the affect of
noise [2, 3] and can be used for error mitigation [4]. There has also been
work on properly accounting for statistical shot noise [5] and adaptation to
specific problems such as combinatorial optimization [6]. Thus, this technique
holds great promise for resolving many practical issues with utilizing quantum
hardware, particularly in the NISQ era [7].

However, circuit cutting suffers greatly in practice as the runtime grows
exponentially with the number of cuts. Akin to quantum state tomography,
circuit cutting works by classically tracking all quantum degrees of freedom
at the cut locations. Thus, exponential scaling comes naturally as the quantum
state of interest grows. There have been many efforts to reduce this cost
through randomized measurements [8, 9], classical sampling [10], and
variational optimization [11]. Nonetheless, the exponential growth in runtime
is unlikely to vanish without imposing structural assumptions on the circuit.
Alternatively, finding applications of circuit cutting that avoid the scaling
issue, as demonstrated in [4], also remains a problem of interest.

In this manuscript, we build upon our previous work [12] and propose an
algorithm for online detection of neglectable basis elements during circuit
cutting. In [12], we showed that some reconstruction procedures can be sped up
if we impose extra assumptions on the circuit—namely, whether a basis element
can be neglected. However, such assumptions cannot be easily detected a
priori. Thus, we propose a hypothesis testing scheme at each cut location
that, at no additional cost in run time, identifies whether there is
statistically significant evidence against the assumption being true. We
empirically demonstrate the viability of our method on the Qiskit simulator
[13], and examine scaling effects with respect to important algorithm
parameters.

The outline of the paper is as follows. Section II re-derives the circuit
cutting in the general bipartition case and introduces the concept of a
“golden cut”—a cut location that has basis elements which can be neglected. We
also show the algorithm for detecting golden cuts as well as their statistical
properties. In Section III, we demonstrate the applicability of our method in
a simple, one-cut case. Meanwhile, we explore additional properties of the
proposed algorithm under varying parameters. Lastly, we discuss some future
directions in Section IV.

##  II Theory

Figure 1: An example of circuit bipartition using $K$ cuts. The circuit is
comprised of arbitrary unitary gates $U_{f_{1}}$ and $U_{f_{2}}$ such that the
middle $K$ qubits are possibly entangled. The circuit can be split into two
fragments by performing $K$ cuts between the two gates, and classically
recombining measurement outcomes in basis $\bm{M}$ and preparing eigenstates
of $\bm{M}$ for a collection of operators that forms a basis over density
matrices.

We start the section by briefly introducing circuit cutting only for the
bipartite case. Readers are encouraged to consult [1, 5, 9] for more detailed
derivations and/or general treatment of circuit cutting. Section II-A ends
with the definition of golden cuts—circuit structure that induces neglectable
basis elements. Then, Section II-B derives an online detection algorithm for
these structures.

###  II-A Circuit Bipartition

Suppose an $N$-qubit quantum circuit induces a state $\rho$. Further suppose
we want to perform $K$ cuts on a quantum circuit to divide it into two
fragments, $f_{1}$ and $f_{2}$ (see Figure 1). Indexing each cut by an integer
in $[K]$, then we can write the cutting scheme as an injective function
$c:[K]\to[N]$ that maps a cut to the respective qubit being cut. The premise
of circuit cutting is to rewrite $\rho$ in terms of the states induced by the
circuit fragments, $\rho_{f_{1}}$ and $\rho_{f_{2}}$, albeit with some
parameterization.

As first shown in Peng et al. [1], this can be done via expanding $\rho$ at
the location of the cut. That is, for a basis set over $2\times 2$ Hermitian
matrices, $\mathcal{B}=\\{I,X,Y,Z\\}$, the following decomposition holds:

$\displaystyle\rho=\frac{1}{2^{K}}\sum_{\bm{M}\in\mathcal{B}^{K}}\rho_{f_{1}}(\bm{M})\otimes\rho_{f_{2}}(\bm{M})$
(1)

where

$\displaystyle\bm{M}=\begin{pmatrix}M_{c(1)}&M_{c(2)}&\dots&M_{c(K)}\end{pmatrix},~{}~{}M_{c(i)}\in\mathcal{B}.$
(2)

The state of each fragment $\rho_{f_{i}}$, $i=1,2,\dots$ is parameterized by
an operator $\bm{M}$ and depends on the particular gates contained in the
circuit. Let $U_{f_{i}}$ denote the unitary operation induced by the quantum
gates on each fragment, and let $|\bm{0}\rangle$ be an $N_{f_{1}}$-qubit
“zero” state. Then, we can write $\rho_{f_{i}}$ as the following (up to
appropriate qubit permutation):

$\displaystyle\rho_{f_{1}}(\bm{M})$
$\displaystyle=\textnormal{tr}_{c(1),\dots,c(K)}\left(\bigotimes_{i\in[K]}M_{c(i)}U_{f_{1}}|\bm{0}\rangle\langle\bm{0}|U_{f_{1}}^{\dagger}\right),$
(3) $\displaystyle\rho_{f_{2}}(\bm{M})$
$\displaystyle=U_{f_{2}}\left(\bigotimes_{i\in[K]}M_{c(i)}\otimes|\bm{0}\rangle\langle\bm{0}|\right)U_{f_{2}}^{\dagger}.$
(4)

The choice of basis is arbitrary, and we chose the normalized Pauli basis for
simplicity. Note that the above equation lacks a physical interpretation as
elements in $\mathcal{B}$ are traceless (except for $I$) and hence, are not
quantum states.

To resolve this issue, we note that each operator $M$ admits spectral
decomposition. Letting

$\displaystyle\bm{r}=\begin{pmatrix}r_{c(1)}&r_{c(2)}&\dots&r_{c(K)}\end{pmatrix}\in\\{\pm
1\\}^{K}$ (5)

be a tuple of eigenvalues, we define

$\displaystyle\bm{M}^{\bm{r}}=\begin{pmatrix}M_{c(1)}^{r(1)}&M_{c(2)}^{r(2)}&\dots&M_{c(K)}^{r(K)}\end{pmatrix}$
(6)

to be the $\bm{r}$-th eigenstate of operator $\bm{M}$. Let $\bm{s}\in\\{\pm
1\\}^{K}$ and $\bm{M}^{\bm{s}}$ be similarly defined. Applying this
decomposition gives the reconstruction formula in the bipartition case:

$\displaystyle\rho=\frac{1}{2^{K}}\sum_{\begin{subarray}{c}\bm{M}\in\mathcal{B}^{K},\\\
\bm{r},\bm{s}\in\\{\pm
1\\}^{K}\end{subarray}}\textnormal{Par}(\bm{r})~{}\textnormal{Par}(\bm{s})~{}\rho_{f_{1}}(\bm{M}^{\bm{r}})\otimes\rho_{f_{2}}(\bm{M}^{\bm{s}})$
(7)

where $\textnormal{Par}(\bm{r})$ denotes the parity of a string of
eigenvalues, i.e., $\textnormal{Par}(\bm{r})=\prod_{i}r_{i}$. The formula
above lends itself to a measure-and-prepare scheme for realizing quantum
circuit cutting: for each basis element $\bm{M}$, we measure the upstream
circuit in the basis, prepare the downstream circuit into the eigenstates of
the same basis, and reweight the outcome of the downstream circuit by the
probability of observing the respective eigenstate upon measuring the upstream
fragment.

Alternatively, for any desired quantum observable $O$, suppose the operator
can be decomposed to accommodate the two fragments, i.e., $O=O_{f_{1}}\otimes
O_{f_{2}}$ up to appropriate permutation of qubit indices. Then, we can arrive
at an expression for the expectation of the uncut circuit in terms of the
fragments $\rho_{f_{i}}$ and their respective observables $O_{f_{i}}$:

$\displaystyle\textnormal{tr}(O\rho)$
$\displaystyle=\frac{1}{2^{K}}\sum_{\bm{M},\bm{r},\bm{s}}\textnormal{Par}(\bm{r})\textnormal{Par}(\bm{s})\textnormal{tr}\left(\left(O_{f_{1}}\otimes
O_{f_{2}}\right)\left(\rho_{f_{1}}\otimes\rho_{f_{2}}\right)\right)$ (8)
$\displaystyle=\frac{1}{2^{K}}\sum_{\bm{M},\bm{r},\bm{s}}\textnormal{Par}(\bm{r})\textnormal{Par}(\bm{s})~{}\textnormal{tr}\left(O_{f_{1}}\rho_{f_{1}}\right)\textnormal{tr}\left(O_{f_{2}}\rho_{f_{2}}\right)$
(9)

where we implicitly apply $\bm{M}$ to $\rho_{f_{1}}$ and $\rho_{f_{2}}$.

Note that the decomposition assumption $O=O_{f_{1}}\otimes O_{f_{2}}$ is
without loss of generality. For any choice of Hermitian operator $O$, one can
expand it with respect to Pauli strings, i.e.,

$\displaystyle O=\sum_{\bm{S}\in\mathcal{B}^{N}}a_{S}~{}\bm{S}$ (10)

for some set of real coefficients $\\{a_{S}\\}$. So, by linearity of the trace
operator, we obtain a generalized expression for the expectation:

$\displaystyle\textnormal{tr}(O\rho)$
$\displaystyle=\frac{1}{2^{K}}\sum_{\bm{S},\bm{M},\bm{r},\bm{s}}a_{S}\textnormal{Par}(\bm{r})\textnormal{Par}(\bm{s})\textnormal{tr}\left(\bm{S}_{f_{1}}\rho_{f_{1}}\right)\textnormal{tr}\left(\bm{S}_{f_{2}}\rho_{f_{2}}\right)$
(11)

where $\bm{S}_{f_{1}}$ and $\bm{S}_{f_{2}}$ are the Pauli strings separated
according to the circuit cutting scheme, i.e.,
$\bm{S}=\bm{S}_{f_{1}}\otimes\bm{S}_{f_{2}}$ under appropriate qubit
permutations.

We now formally define the golden circuit cutting point.

######  Definition 1.

Consider an $N$-qubit circuit amenable to bipartition with $K$ cuts. We’re
interested in the expectation of the circuit-induced state with respect to
some quantum observable $O=O_{f_{1}}\otimes O_{f_{2}}$. The cutting scheme
admits a _golden cutting point_ if there exists $\bm{M}_{*}\in\mathcal{B}^{K}$
such that

$\displaystyle\sum_{\bm{r}\in\\{\pm
1\\}^{K}}\textnormal{Par}(\bm{r})~{}\textnormal{tr}\left(O_{f_{1}}\rho_{f_{1}}(\bm{M}_{*}^{\bm{r}})\right)=0$
(12)

More simply put, a golden cutting point refers to the existence of a basis
element that leads to systematic cancellations. Golden cutting points neither
necessarily exists nor are unique. For each such basis element, one does not
need to execute the circuit downstream of the cut with initialization
corresponding to the neglected basis.

Golden cutting points can be constructed via circuit design by restricting the
set of rotations allowed prior to cutting so long as the structure of the
quantum circuit permits. However, one should not expect such a property to
hold for an arbitrary algorithm. Thus, we propose an “online” scheme for
detecting the existence of golden cutting points and establish robustness of
misidentifying golden cuts.

###  II-B Identifying Golden Cutting Points

With no knowledge of the existence of golden cutting points, one must execute
each of the $4^{K}$ upstream circuits and another $4^{K}$ downstream circuits
($f_{1}$ and $f_{2}$ respectively in the bipartition case). To detect golden
cutting points in the absence of a priori knowledge, we propose to conduct a
hypothesis test for each of the $4^{K}$ upstream circuits, determine whether
there is statistically significant evidence for the existence of a golden
cutting point, then run the corresponding downstream circuit.

Denote

$\displaystyle\tau=\sum_{\bm{r}\in\\{\pm
1\\}^{K}}\textnormal{Par}(\bm{r})~{}\textnormal{tr}\left(O_{f_{1}}\rho_{f_{1}}(\bm{M}_{*}^{\bm{r}})\right)$
(13)

as the quantity we want to verify magnitude of. Inheriting the bipartition
assumption from the previous section, we can rewrite $\tau$ as estimating the
expectation of a larger observable

$\displaystyle\tau=\textnormal{tr}\left(\left(O_{f_{1}}\otimes\bm{M}_{*}\right)U_{f_{1}}|\bm{0}\rangle\langle\bm{0}|U_{f_{1}}^{\dagger}\right).$
(14)

where, again, $U_{f_{1}}|\bm{0}\rangle$ is the state induced by the upstream
fragment. Writing $\tau$ in this form allows us to employ standard techniques
for estimating quantum observables.

Assume for convenience that $O_{f_{1}}$ is a Pauli-string. To estimate the
expectation, we measure each qubit in the respective Pauli basis (by
performing a rotation $V$) $m$ times and obtain an ensemble of bitstring
samples $\\{\hat{b}_{i}\\}_{i=1}^{m}$. Therefore, we can estimate $\tau$ be
constructing

$\displaystyle\hat{\tau}=\frac{1}{m}\sum_{i=1}^{m}\langle\hat{b}_{i}|V^{\dagger}(O_{f_{1}}\otimes\bm{M})V|\hat{b}_{i}\rangle.$
(15)

Alternatively, one can think of estimating the distribution of strings of
eigenvalues (which we’ll call eigenstrings for short) induced by the
measurements. Write $p_{b}$ for the probability of obtaining eigenstring
$\bm{b}$, and $\hat{p}_{b}$ for the empirical probability. Moreover, let
$\bm{p}$ and $\hat{\bm{p}}$ denote the vector of probabilities. Then, we can
sum the parity of each eigenstring weighted by the (empirical) probability to
arrive at an alternative expression for the estimator:

$\displaystyle\hat{\tau}=\sum_{\bm{b}\in\\{\pm
1\\}^{N_{f_{1}}}}\textnormal{Par}(\bm{b})~{}\hat{p}_{b}.$ (16)

The proposition below establishes the standard error and asymptotic normality,
which are convenient for hypothesis testing.

######  Proposition 1 (Asymptotic Normality).

Given a circuit amenable to the bipartite circuit cutting scheme (cf. II-A),
let $\hat{\tau}$ be the estimator of $\tau$ expressed in Equation 15 and let
$O_{f_{1}}$ admit decomposition as in Equation 10. Then, $\hat{\tau}$ is
asymptotically normal, i.e.,

$\displaystyle\frac{\hat{\tau}-\tau}{\textnormal{std}(\hat{\tau})}\to\mathcal{N}(0,1)$
(17)

where the standard deviation of the estimator is expressed as

$\displaystyle\textnormal{std}(\hat{\tau})=\left(\sum_{S\in\mathcal{B}^{N_{f_{1}}}}\frac{a_{S}^{2}}{N}~{}\bm{\chi}^{\intercal}(\textnormal{diag}(\hat{\bm{p}_{S}})-\hat{\bm{p}_{S}}\hat{\bm{p}_{S}}^{\intercal})\bm{\chi}\right)^{1/2}$
(18)

and $\bm{\chi}$ is the vector of parities, i.e.,
$\bm{\chi}_{b}=\textnormal{Par}(b).$

###### Proof.

We first consider the case were $O_{f_{1}}$ is a Pauli string $\bm{S}$, then
proceed to generalize to arbitrary quantum observables.

Estimating the expectation of the Pauli string $\bm{S}$ with $m$ shots, using
the formalism presented in Equation 16, gives the standard error

$\displaystyle\textnormal{Var}(\hat{\tau})$
$\displaystyle=\textnormal{Cov}\left(\sum_{b}\hat{p}_{b}~{}\textnormal{Par}(b),\sum_{b^{\prime}}\hat{p}_{b^{\prime}}~{}\textnormal{Par}(b^{\prime})\right)$
(19)
$\displaystyle=\sum_{b,b^{\prime}}\textnormal{Cov}(\hat{p}_{b},\hat{p}_{b^{\prime}})\textnormal{Par}(b)\textnormal{Par}(b^{\prime})$
(20)
$\displaystyle=\frac{1}{m}\left(\sum_{b=b^{\prime}}p_{b}(1-p_{b})-\sum_{b\neq
b^{\prime}}p_{b}p_{b^{\prime}}\textnormal{Par}(b)\textnormal{Par}(b^{\prime})\right)$
(21)

where the third equality follows from the covariance of multinomial
distributions. Using the empirical quantities for each $p_{b}$ and writing in
matrix form gives the estimated standard error

$\displaystyle\textnormal{std}(\hat{\tau})=\sqrt{\frac{1}{m}\bm{\chi}^{\intercal}(\textnormal{diag}(\hat{\bm{p}_{S}})-\hat{\bm{p}_{S}}\hat{\bm{p}_{S}}^{\intercal})\bm{\chi}}.$
(22)

Consider the decomposition in Equation 10. As the estimation procedure runs
independently, variances add. Hence, we arrive at the form in Equation 18.

Lastly, asymptotic normality is established by the equivalence formulation
presented in Equation 15 and Equation 16. In the form of Equation 15, we can
express $\hat{\tau}$ as the average over independent samples. In combination
with finiteness of $O_{f_{1}}$, the Central Limit Theorem holds, implying
asymptotic normality of $\hat{\tau}$. ∎

Using the above proposition, we can deduce an algorithm for detecting golden
cutting points. For each basis element in $\mathcal{B}^{K}$, we will compute
$\hat{\tau}$ and perform a statistical test for whether $\tau\neq 0$. If we’ve
gathered statistically significant evidence for $\tau$ being non-zero, then we
would run the downstream fragment parameterized by the respective basis
element. On the other hand, if $\hat{\tau}$ is sufficiently close to zero,
then we classify the cut as a golden cutting point and proceed without running
the corresponding downstream fragment. Using $\Phi$ to denote the CDF of a
standard Gaussian, we summarize the above procedure in Algorithm 1.

Input: fragments $f_{1}$ and $f_{2}$, observable $O_{f_{1}}$, $O_{f_{2}}$,
significance level $\alpha\in(0,1)$

Output: Expectation $\textnormal{tr}(O\rho)$

1 for  _$\bm{M}\in\mathcal{B}^{K}$_ do

2 Compute $\hat{\tau}$ using Eqn. 16

3 if  _$|\hat{\tau}| >\Phi^{-1}(1-\alpha)\cdot\textnormal{std}(\hat{\tau})$_
then

4 Reject the hypothesis and compute
$\textnormal{Par}(\bm{s})~{}\textnormal{tr}(O_{f_{2}}\rho_{f_{2}}(\bm{M}^{\bm{s}})$
for all $s\in\\{-1,+1\\}^{K}$

5 else

6 Fail to reject the hypothesis and set quantities related to $\bm{M}$ to zero

7

Reconstruct the full expectation from fragment data using Eqn. 9

Algorithm 1 Online detection of golden cutting points Figure 2: Behavior of
Algorithm 1 for varied shots and $\alpha$ levels averaged over 1000
independent trials. The rate of true positives (A) is consistent with the
specified $\alpha$, and the true negative rate (B) converges to one as the
number of shots increases. The reconstruction error monotonically vanishes
both in the presence (C) and the absence (D) of golden cuts.

While we can control the rate of correctly identifying golden cutting points
by the significance level $\alpha$, we would also like to derive a way of
controlling the rate of false negatives. In fact, falsely identifying a non-
golden cutting point as golden is more problematic than falsely identifying
golden as non-golden. This is because in the latter case the reconstructed
bitstring distribution will likely still be within acceptance error ranges,
but that’s not the case for the former. Thus, we hope to lower the probability
of false negatives by taking a sufficient number of shots during the
estimation procedure.

######  Remark 1.

As we have showed that the estimator converges weakly to a normal
distribution, we will assume Gaussianity to facilitate analysis. Suppose
$\hat{\tau}\sim\mathcal{N}(\tau,b^{2}/m)$ where
$b=\sqrt{m}\cdot\textnormal{std}(\hat{\tau})$. We hope that for a basis
element where $\tau>\epsilon>0$, the hypothesis testing scheme would reject
it, perhaps with a small probability of error $\delta$. By the Chernoff bound,
we know that

$\displaystyle\Pr(|\hat{\tau}-\tau|>\epsilon)\leq
2e^{-m\epsilon^{2}/2b^{2}}=\delta.$ (23)

Thus, to estimate $\tau$ to any desired $\epsilon$ accuracy with probability
$1-\delta$, we need

$\displaystyle m\geq\frac{2b^{2}}{\epsilon^{2}}\log\frac{2}{\delta}$ (24)

measurements. Since $b$ is not known a priori, we can upper bound it by

$\displaystyle\bm{1}^{\intercal}(\textnormal{diag}(\bm{q})+\bm{q}\bm{q}^{\intercal})\bm{1}.$
(25)

The distribution $\bm{q}$ that maximizes the above quantity is the uniform
distribution. Thus, we can arrive at a definitive upper bound
$b\leq\frac{3}{2}(1-2^{-N_{f_{1}}})$.

Alternatively, one can interpret the proposed sample complexity as accepting
an $\epsilon$ margin for Equation 12 in the sense that we declare a cut is
golden if $|\hat{\tau}|<\epsilon$, thereby accepting an additional additive
error of magnitude $\mathcal{O}(\epsilon)$ to the result of the circuit
reconstruction. On the other hand, if $\tau>\epsilon$, we wish to identify it
with probability $1-\delta$. Note that in the limit of $m\to\infty$, variance
vanishes and the true positive and negative rates approach one.

##  III Experiments

In this section, we numerically demonstrate the applicability of Algorithm 1.
The algorithm was implemented in Qiskit and executed on the Aer simulator. We
examine its statistical (Section III-A) and runtime (Section III-B) properties
through studying its dependency on the number of shots and the significance
level $\alpha$.

###  III-A Statistical Analysis

As is standard in analyzing binary decisions, we analyzed our statistical test
by providing instances when the null hypothesis is true and when it’s false.
Specifically, we provided circuits either with or without golden cuts, and
observed the probability of correctly identifying the existence (or lack) of
golden cuts.

We consider a simple three-qubit circuit amenable to cutting on the second
qubit. First, we generated a circuit containing a golden cut by appending two
$R_{X}$ gates on the first and second qubits, then an $R_{Y}$ gate only on the
first qubit—this is $U_{f_{1}}$ from Figure 1, and we let $K=1$. Rotation
angle $\theta$ was set to a value far from zero ($\theta=0.5$) to ensure $X$
would be the only golden cutting axis and to focus only on statistical shot
noise. To generate a circuit known to not contain a golden cutting point, we
applied the same procedure then appended an additional $R_{Y}$ gate on the
second qubit (the qubit being cut) as well. Finally, we generated $U_{f_{2}}$
randomly across qubits 2 and 3 using Qiskit’s `random_circuit` function. Once
the circuit was constructed, for each shot count-$\alpha$ pair, we repeated
1000 independent executions of Algorithm 1 and collected the frequency at
which the algorithm correctly identified the circuit structure. Results are
displayed in subplots A and B of Figure 2.

Subplot A shows the probability of failing to reject the null hypothesis given
the null hypothesis is true, which should be exactly $\alpha$. The numerics
aligned with the theoretical value with the exception of cases with low shot
counts. This can be understood as our estimator is built upon asymptotic
statements on the sampling distribution. Subplot B demonstrates the rate of
true negatives. We can see that, given sufficient samples, we always
identified non-golden cuts correctly. For lower significance values, we are
more prone to rejecting the null hypothesis, explaining the faster rate of
convergence towards 1 for smaller $\alpha$-values.

We also examined the quality of the reconstruction by calculating the distance
between the empirical and theoretical bitstring distributions. The theoretical
distribution is obtained by taking large number of shots without circuit
cutting. We employed the $\ell^{2}$-distance to quantify how far apart two
distributions are, i.e., for discrete distributions $\bm{p}$ and $\bm{q}$,

$\displaystyle d(\bm{p},\bm{q})=\sqrt{\sum_{i}(p_{i}-q_{i})^{2}}.$ (26)

Again, we executed Algorithm 1 independently for 1000 trials and collected the
$\ell^{2}$ distance between the empirical, reconstructed bitstring
distribution and the respective theoretical distribution at varying numbers of
circuit execution shots and alpha levels. Results are found in subplots C and
D of Figure 2.

In general, the reconstruction error decreases monotonically with the number
of shots, and there was not a significant difference among choices of
$\alpha$. In the case of low shot count and no existing golden cut, there
seems to be more statistical fluctuation when reconstructing. Considering
subplot B above, we know that this region is prone to false negatives, and
thus neglecting bases that should not be neglected.

###  III-B Runtime Analysis

| run time (sec)
---|---
| w/ optimization | w/o optimization
$\alpha$ | $10^{-1}$ | 0.0771$\pm$0.0006 | 0.0959$\pm$0.0004
$10^{-2}$ | 0.0749$\pm$0.0004 | 0.0961$\pm$0.0004
$10^{-3}$ | 0.0747$\pm$0.0003 | 0.0962$\pm$0.0004
TABLE I: Runtime comparison between circuit cutting procedures with and
without optimization from Algorithm 1. Independent trials were repeated 1000
times. We can see that neglecting basis elements consistently run faster
despite spending computing overhead on hypothesis testing.

To obtain timed runtime values, we generated a circuit with a golden cut and
executed Algorithm 1. Recall that, depending on the results of the hypothesis
test on the upstream circuit, the downstream circuit might not be executed for
certain bases. Then, we ran the same cut circuit and performed the usual
reconstruction routine without hypothesis testing or golden cutting
optimization. Both of these processes were timed individually over 1000 trials
and at varying alpha levels. Results for this can be found in Table I. In
general, we see roughly a $20\%$ decrease in runtime upon performing the
optimization. As $\alpha$ decreases, we tend to reject the null hypothesis
more often, thereby executing the downstream circuit more often.

##  IV Conclusion

In this manuscript, we proposed an online detection algorithm for finding
golden cutting points—circuit structures that induce neglectable basis
elements during circuit reconstruction. The detection algorithm was built on
performing a hypothesis test for each basis element, and executing the
respective downstream circuit only if the null hypothesis is rejected. The
detection does not require additional circuit executions. We showed
numerically that, under sufficient number of shots, golden cuts will always be
detected and there is no drastic difference among choices of significance
levels.

Empirically testing the algorithm on quantum hardware and at large scale—both
of which contribute additional noise that can affect the quality of
estimators—are left to future work. Another immediate open question is the
prevalence of circuit structure amenable to golden cuts in applicable
circuits. For instance, we found that the SupermarQ [14] and QASMBench [15]
benchmark suites both feature a handful of benchmarks that exhibit this
circuit structure. Variational circuits whose ansatz can be flexible might
also be a candidate to apply golden cutting point restrictions for better
scalability.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers. The
colors used in figures were based on Nord Theme [16].

## References

  * [1] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, “Simulating large quantum circuits on a small quantum computer,” _Physical Review Letters_ , vol. 125, no. 15, oct 2020. [Online]. Available: https://doi.org/10.1103/PhysRevLett.125.150504
  * [2] T. Ayral, F.-M. Le Régent, Z. Saleem, Y. Alexeev, and M. Suchara, “Quantum divide and compute: Hardware demonstrations and noisy simulations,” in _2020 IEEE Computer Society Annual Symposium on VLSI (ISVLSI)_. IEEE, 2020, pp. 138–140. [Online]. Available: https://ieeexplore.ieee.org/document/9155024
  * [3] T. Ayral, F.-M. L. Régent, Z. Saleem, Y. Alexeev, and M. Suchara, “Quantum divide and compute: exploring the effect of different noise sources,” _SN Computer Science_ , vol. 2, no. 3, pp. 1–14, 2021. [Online]. Available: https://doi.org/10.1007/s42979-021-00508-9
  * [4] J. Liu, A. Gonzales, and Z. H. Saleem, “Classical simulators as quantum error mitigators via circuit cutting,” _arXiv preprint arXiv:2212.07335_ , 2022\. [Online]. Available: https://arxiv.org/abs/2212.07335
  * [5] M. A. Perlin, Z. H. Saleem, M. Suchara, and J. C. Osborn, “Quantum circuit cutting with maximum-likelihood tomography,” _npj Quantum Information_ , vol. 7, no. 1, p. 64, 2021. [Online]. Available: https://doi.org/10.1038/s41534-021-00390-6
  * [6] Z. H. Saleem, T. Tomesh, M. A. Perlin, P. Gokhale, and M. Suchara, “Quantum Divide and Conquer for Combinatorial Optimization and Distributed Computing,” _arXiv preprint_ , Jul. 2021. [Online]. Available: https://arxiv.org/abs/2107.07532
  * [7] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, 2018.
  * [8] A. Lowe, M. Medvidović, A. Hayes, L. J. O’Riordan, T. R. Bromley, J. M. Arrazola, and N. Killoran, “Fast quantum circuit cutting with randomized measurements,” 2022. [Online]. Available: https://arxiv.org/abs/2207.14734
  * [9] D. T. Chen, Z. H. Saleem, and M. A. Perlin, “Quantum divide and conquer for classical shadows,” _arXiv preprint arXiv:2212.00761_ , 2022. [Online]. Available: https://arxiv.org/abs/2212.00761
  * [10] D. Chen, B. Baheri, V. Chaudhary, Q. Guan, N. Xie, and S. Xu, “Approximate quantum circuit reconstruction,” in _2022 IEEE International Conference on Quantum Computing and Engineering (QCE)_. IEEE, 2022, pp. 509–515.
  * [11] G. Uchehara, T. M. Aamodt, and O. Di Matteo, “Rotation-inspired circuit cut optimization,” _arXiv preprint arXiv:2211.07358_ , 2022. [Online]. Available: https://arxiv.org/abs/2211.07358
  * [12] D. T. Chen, E. H. Hansen, X. Li, V. Kulkarni, V. Chaudhary, B. Ren, Q. Guan, S. Kuppannagari, J. Liu, and S. Xu, “Efficient quantum circuit cutting by neglecting basis elements,” _arXiv preprint arXiv:2304.04093_ , 2023.
  * [13] A. tA-v et al., “Qiskit: An open-source framework for quantum computing,” 2021.
  * [14] T. Tomesh, P. Gokhale, V. Omole, G. S. Ravi, K. N. Smith, J. Viszlai, X.-C. Wu, N. Hardavellas, M. R. Martonosi, and F. T. Chong, “Supermarq: A scalable quantum benchmark suite,” in _2022 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_ , 2022, pp. 587–603.
  * [15] A. Li, S. Stein, S. Krishnamoorthy, and J. Ang, “Qasmbench: A low-level quantum benchmark suite for nisq evaluation and simulation,” _ACM Transactions on Quantum Computing_ , vol. 4, no. 2, feb 2023. [Online]. Available: https://doi.org/10.1145/3550488
  * [16] S. Greb, “Nord theme.” [Online]. Available: https://www.nordtheme.com/


# Using cGANs for Anomaly Detection: Identifying Astronomical Anomalies in
JWST Imaging

Ruby Pearce-Casey School of Physical Sciences
The Open University, Milton Keynes, MK7 6AA UK  Hugh Dickinson  School of
Physical Sciences
The Open University, Milton Keynes, MK7 6AA UK  Stephen Serjeant  School of
Physical Sciences
The Open University, Milton Keynes, MK7 6AA UK  Jane Bromley  School of
Computing
The Open University, Milton Keynes MK7 6AA UK

††journal: AAS Journal††software: Astropy (Astropy Collaboration et al., 2013,
2018, 2022)††software: Sep (Barbary, 2016; Bertin & Arnouts, 1996)

[latex]

##  1 Abstract

We present a proof of concept for mining JWST imaging data for anomalous
galaxy populations using a _conditional Generative Adversarial Network_
(cGAN). We train our model to predict long wavelength NIRcam fluxes
($LW:F277W,F356W,F444W$ between $2.4\mu m-5.0\mu m$) from short wavelength
fluxes ($SW:F115W,F150W,F200W$ between $0.6\mu m-2.3\mu m$) in $\sim 2000$
galaxies. We test the cGAN on a population of 37 Extremely Red Objects (EROs)
discovered by the CEERS JWST Team (Barro et al., 2023). Despite their red long
wavelength colours, the EROs have blue short wavelength colours
($F150W-F200W\sim 0\,\text{mag}$) indicative of bimodal SEDs. Surprisingly,
given their unusual SEDs, we find that the cGAN accurately predicts the _LW_
NIRcam fluxes of the EROs. However, it fails to predict _LW_ fluxes for other
rare astronomical objects, such as a merger between two galaxies, suggesting
that the cGAN can be used to detect some anomalies.111The source code and
trained model are available at
https://github.com/RubyPC/Anomaly_Detection_with_cGANs.

##  2 Method

We use public domain JWST NIRcam data in the CEERS field, which covers
$100sq.arcmin$ of the EGS field, and our network is based on the _pix2pix_
architecture described in Isola et al. (2017). The data for the model were
obtained from the CEERS data release $0.5$. Due to the different sensitivities
of Module A and B of the JWST imaging, this study focussed on the subset of
galaxies that were imaged using Module A. The images for each NIRcam filter
were cropped into two regions of the NIRcam imaging field and a PSF-matching
was done to match the PSF of each filter to that of F444W. Sources were
detected using Python’s _sep_
package222https://sep.readthedocs.io/en/v1.1.x/(Barbary, 2016; Bertin &
Arnouts, 1996) using $F277W$ as the detection filter. A pixel threshold value
$\geq 2.0\sigma$ was used with $\sigma$, the global background RMS. This
detection threshold is the minimum signal-to-noise ($S/N$) ratio of a pixel to
be regarded as a detection. All detected sources were extracted as $90\times
90$ pixel cutout images using Python’s _Astropy_
package333https://readthedocs.org/projects/astropy/(Astropy Collaboration et
al., 2013, 2018, 2022).

The cGAN was pre-trained on _SW_ NIRcam filters for a subset of sources
detected to predict _LW_ filters using $90\%$ of the $\sim 2000$ galaxies for
training and $10\%$ for validation. After training the cGAN for $100$ epochs,
the validation dataset was used to assess the generalisation capability of the
model. To assess the performance of the cGAN, photometry was measured on the
predicted _LW_ fluxes to compare the SED of the objects with those predicted
by the model. Figure 1a shows the distribution of the absolute difference
between the ground truth SED, $y$, and the SED predicted by the cGAN,
$\hat{y}$, against a subset of the validation data set across the three _LW_
NIRcam filters. The absolute difference is predominantly situated at $\sim 0$
for the majority of sources for all of the _LW_ filters. The distribution is
roughly uniform for the $F356W$ filter where there is a small difference
between the ground truth SED and the predicted SED, respectively.
Nevertheless, this shows that the model can accurately predict NIRcam _LW_
filters such that we can reliably construct an SED.

##  3 Results

We tested the cGAN on the $37$ EROs identified by Barro et al. (2023) with
$F444W<28\,\text{mag}$ and photometric redshifts $5<z<7$. The $37$ EROs were
deliberately excluded from the model's training dataset. All of the EROs are
unresolved, point-like sources in all observed filters and are typically faint
in all _SW_ NIRcam filters. Barro et al. (2023) found that the EROs exhibit a
distinctive, bimodal SED with blue colours for $\lambda>2\mu m$ and they
suggest that the flux in some of the _LW_ filters is partially boosted by
strong emission lines, making the colours redder than the stellar continuum. A
defining feature of these EROs is that _all_ of them have blue colours in the
_SW_ NIRcam filters ($F115W-F200W\sim 0$). Figure 1b shows the cGAN's
prediction on $3$ of the EROs with the top row being the _SW_ input to the
model and the middle and bottom rows being the cGAN's prediction of the _LW_
and the ground truth _LW_ , respectively. Despite the peculiar, bimodal SEDs
of the EROs, the cGAN performs just as well in predicting the _LW_ filters for
a sample of EROs, predicting red _LW_ NIRcam fluxes. Additionally, Figure 1c
shows the comparison between the ground truth SED of an identified ERO and the
corresponding prediction from the cGAN. Although EROs are rare, and with a
peculiar SED, the cGAN accurately predicts the steep rise in the SED at longer
NIR wavelengths.

Figure 1: _The results of the $cGAN$. 1a shows the distribution of the
absolute difference in SED between the ground truth and the
$cGAN\textquotesingle s$ prediction, respectively, after testing on a subset
of the validation data set of $\sim 2000$ objects. 1b shows the
$cGAN\textquotesingle s$ prediction of the _LW_ channel for a set of $3$ EROs
in the middle row comparing to the ground truth in the bottom row with the top
row being the _SW_ input to the model, and 1c illustrates the
$cGAN\textquotesingle s$ prediction for one ERO, showing its generalisation
power of predicting the _LW_ fluxes of objects with peculiar SEDs by
excellently constructing the characteristic steep rise in the $SED$. Finally,
1d shows the $cGAN\textquotesingle s$ poor prediction on a pair of merging
galaxies with 1e showing the poorly constructed SED, where the cGAN can be
used as an anomaly detector._

However, we find that the cGAN cannot accurately predict _LW_ filters for
_all_ galaxy types. In particular, we find its predictions for galaxies with
unusual or disturbed morphologies are less accurate, such as can be seen in
Figure 1d where the prediction of the _LW_ in the middle row is not accurate
in comparison to the ground truth in the bottom row for a pair of merging
galaxies. Figure 1e shows the construction of the SED for the galaxy merger
and it is clear that the model fails to predict the _LW_ NIRcam fluxes. Thus,
we suggest that the cGAN can be used as an anomaly detector; it can accurately
predict _LW_ filters for the majority of galaxies in the CEERS field,
including EROs, but its performance suffers when predicting on rare objects
such as a merger between galaxies.

##  4 Conclusion

In conclusion, our results demonstrate our cGAN's capability of detecting some
rare astronomical objects. Therefore, for this specific case using JWST NIRcam
filters, not only can cGANs be used for image generation and image
segmentation tasks but also for anomaly detection. This technology can be
deployed on the rapidly increasing archive of JWST imaging survey observations
as well as forthcoming giant survey data sets from e.g. Roman and Euclid;
anomalous systems such as strong gravitational lenses have already been seen
in JWST surveys (e.g. Pearson et al., 2023; van Dokkum et al., 2023).

##  5 Acknowledgements

Ruby Pearce-Casey thanks the Science and Technology Facilities Council for
support under grant ST/W006839/1. This work made use of Astropy: a community-
developed core Python package and an ecosystem of tools and resources for
astronomy, and Sep. .

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barbary (2016) Barbary, K. 2016, Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058
  * Barro et al. (2023) Barro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, Extremely red galaxies at $z=5-9$ with MIRI and NIRSpec: dusty galaxies or obscured AGNs?  https://arxiv.org/abs/2305.14418
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Isola et al. (2017) Isola, P., Zhu, J.-Y., Zhou, T., & Efros, A. A. 2017, in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 5967–5976, doi: 10.1109/CVPR.2017.632
  * Pearson et al. (2023) Pearson, J., Serjeant, S., Wang, W.-H., et al. 2023, arXiv e-prints, arXiv:2309.00888, doi: 10.48550/arXiv.2309.00888
  * van Dokkum et al. (2023) van Dokkum, P., Brammer, G., Wang, B., Leja, J., & Conroy, C. 2023, arXiv e-prints, arXiv:2309.07969, doi: 10.48550/arXiv.2309.07969


# Exponential Utility Maximization in a Discrete Time Gaussian Framework

Yan Dolinsky  Department of Statistics, Hebrew University of Jerusalem, Israel
<EMAIL_ADDRESS>and  Or Zuk  Department of Statistics, Hebrew
University of Jerusalem, Israel<EMAIL_ADDRESS>
###### Abstract.

The aim of this short note is to present a solution to the discrete time
exponential utility maximization problem in a case where the underlying asset
has a multivariate normal distribution. In addition to the usual setting
considered in Mathematical Finance, we also consider an investor who is
informed about the risky asset’s price changes with a delay ${D}$. Our method
of solution is based on the theory developed in [4] and guessing the optimal
portfolio.

Y. Dolinsky is supported in part by the GIF Grant 1489-304.6/2019 and the ISF
grant 230/21.
O. Zuk is supported in part by the ISF grant 2392/22.

Mathematical Subject Classification (2010):

91B16, 91G10

Keywords:

Utility Maximization, Hedging with Delay, Banded Matrices

##  1\. Introduction and the Main Result

Taking into account frictions is an important challenge in financial
modelling. In this note, we focus on the friction arising from the fact that
investment decisions may be based only on delayed information, and the actual
present market price is unknown at the time of decision making. This
corresponds to the case where there is a time delay in receiving market
information (or in applying it), which causes the trader’s filtration to be
delayed with respect to the price filtration.

We start by briefly reviewing some of the relevant literature. In [18] the
author solved the mean-variance hedging problem with partial information for
European contingent claims in a setup where the asset price process is a
martingale. Subsequently, the study of the mean-variance hedging problem was
extended beyond the martingale setup (see for instance [7, 13]). The work [15]
applied Malliavin calculus to the utility maximization problem, and in [17]
the authors introduce and study a general framework of optimal stochastic
control with delayed information. Another related work is [2] which deals with
scaling limits of utility indifference prices with vanishing delay and risk
aversion inversely proportional to the delay.

To formulate our main result, let $(\Omega,\mathcal{F},\mathbb{P})$ be a
complete probability space carrying one risky asset which we denote by
$S=(S_{k})_{0\leq k\leq n}$ where $n\in\mathbb{N}$ is a fixed finite time
horizon. We assume that the investor has a bank account that, for simplicity,
bears no interest. We assume that the initial stock price $S_{0}$ is a known
constant and the random vector $(S_{1},...,S_{n})$ has a non-degenerate
multivariate normal distribution. Hence,

$(S_{1}-S_{0},S_{2}-S_{1},...,S_{n}-S_{n-1})\sim\mathcal{N}(\mu,{\Sigma})$

where $\mu=(\mu_{1},...,\mu_{n})\in\mathbb{R}^{n}$ is a constant (row) vector
and ${\Sigma}\in M_{n}(\mathbb{R})$ (as usual, $M_{n}(\mathbb{R})$ denotes the
set of all $n\times n$ real matrices) is a constant positive definite matrix.
Let ${\Lambda}:={\Sigma}^{-1}$ be the precision matrix.

We fix a nonnegative integer number ${D}\in\mathbb{Z}_{+}$ and consider a
situation where there is a delay of ${D}$ trading times. Hence, the investor’s
flow of information is given by the filtration
$\mathcal{G}^{{D}}_{k}:=\sigma\left\\{S_{0},...,S_{\left(k-{D}\right)^{+}}\right\\}$,
$k=0,1,...,n$. The case ${D}=0$ means that there is no delay and corresponds
to the usual setting.

A trading strategy is a random vector $\gamma=(\gamma_{1},...,\gamma_{n})$
which is predictable with respect to the above filtration. Namely, for any
$k$, $\gamma_{k}$ is a $\mathcal{G}^{{D}}_{k-1}$ measurable random variable
($\gamma_{k}$ is the number of shares at time $k-1$). Denote by
$\mathcal{A}_{{D}}$ the set of all trading strategies. For
$\gamma\in\mathcal{A}_{{D}}$ the corresponding portfolio value at the maturity
date is given by

$V^{\gamma}_{n}=\sum_{i=1}^{n}\gamma_{i}(S_{i}-S_{i-1}).$

The investor’s preferences are described by an exponential utility function

$u(x)=-\exp(-\alpha x),\ \ x\in\mathbb{R},$

with absolute risk aversion parameter $\alpha>0$, and her goal is to

(1.1) $\mbox{Maximize}\ \ \mathbb{E}_{\mathbb{P}}\left[-\exp\left(-\alpha
V^{\gamma}_{n}\right)\right]\ \ \mbox{over}\ \ \gamma\in\mathcal{A}_{{D}}$

where $\mathbb{E}_{\mathbb{P}}$ denotes the expectation with respect to the
market probability measure $\mathbb{P}$. Clearly, for any portfolio strategy
$\gamma$ and a constant $\lambda\in\mathbb{R}$, $V^{\lambda\gamma}_{n}=\lambda
V^{\gamma}_{n}$. Thus, without loss of generality, we take the risk aversion
$\alpha=1$.

Next, we introduce some notations. Let ${\mathcal{S}}_{{D}}\subset
M_{n}(\mathbb{R})$ be the set of all positive definite matrices $Q\in
M_{n}(\mathbb{R})$ that satisfy $Q_{ij}=0$ for $|i-j|>{D}$. Namely,
${\mathcal{S}}_{{D}}$ is the set of all banded positive definite matrices with
lower bandwidth and upper bandwidth equal to ${D}$. Let
${\mathcal{T}}_{{D}}\subset M_{n}(\mathbb{R})$ be the set of all symmetric
matrices $\Gamma\in M_{n}(\mathbb{R})$ that satisfy $\Gamma_{ij}=0$ for
$|i-j|\leq{D}.$

We arrive at the main result of the paper.

######  Theorem 1.1.

There exists a unique decomposition

(1.2) ${\Lambda}=\hat{Q}^{-1}+\hat{\Gamma}$

where $\hat{Q}\in{\mathcal{S}}_{{D}}$ and
$\hat{\Gamma}\in{\mathcal{T}}_{{D}}$. The maximizer
$\hat{\gamma}=(\hat{\gamma}_{1},...,\hat{\gamma}_{n})$ for the optimization
problem (1.1) is unique and is given by the linear form

(1.3)
$\hat{\gamma}_{i}=\sum_{j=1}^{n}{\Lambda}_{ij}\mu_{j}-\sum_{j=1}^{i-1}\hat{\Gamma}_{ij}(S_{j}-S_{j-1}),\
\ i=1,...,n.$

The corresponding value is given by

(1.4)
$\mathbb{E}_{\mathbb{P}}\left[-\exp\left(-V^{\hat{\gamma}}_{n}\right)\right]=-\sqrt{\frac{|\hat{Q}|}{|\Sigma|}}\exp\left(-\frac{1}{2}\mu{\Lambda}\mu^{\prime}\right)$

where ′ denotes transposition and $|\cdot|$ is the determinant of the matrix
$\cdot$.

The proof of Theorem 1.1 is presented in the next section.

######  Remark 1.2.

There are two trivial cases where the optimal trading strategy is
deterministic. The first one is the case of no information ${D}=n-1$. The
second case is where ${\Sigma}$ is a diagonal matrix, i.e. the stock
increments $S_{1}-S_{0},...,S_{n}-S_{n-1}$ are independent. In both of these
cases we have $\hat{Q}=\Sigma$, $\hat{\Gamma}=0$ and the value is
$\exp\left(-\frac{1}{2}\mu{\Lambda}\mu^{\prime}\right)$. In particular, if the
stock increments $S_{1}-S_{0},...,S_{n}-S_{n-1}$ are independent the optimal
(deterministic) strategy does not depend on ${D}$ and given by
$\hat{\gamma}_{i}=\frac{\mu_{i}}{\Sigma_{ii}}$, $i=1,...,n$.

For the case ${D}\in\\{0,1\\}$ the decomposition (1.2) is explicit provided
that the inverse matrix ${\Lambda}$ is known.

######  Corollary 1.3.

(I) For ${D}=0$ (no delay) the optimal strategy is given by

$\hat{\gamma}_{i}=\sum_{j=1}^{n}{\Lambda}_{ij}\mu_{j}-\sum_{j=1}^{i-1}{\Lambda}_{ij}(S_{j}-S_{j-1}),\
\ i=1,...,n.$

The value is

$\mathbb{E}_{\mathbb{P}}\left[-\exp\left(-V^{\hat{\gamma}}_{n}\right)\right]=-\frac{1}{\left(|\Sigma|\prod\limits_{i=1}^{n}{\Lambda}_{ii}\right)^{1/2}}\exp\left(-\frac{1}{2}\mu{\Lambda}\mu^{\prime}\right)$

(II) For ${D}=1$ the optimal strategy is given by

$\hat{\gamma}_{i}=\sum_{j=1}^{n}{\Lambda}_{ij}\mu_{j}-\sum_{j=1}^{i-2}\left({\Lambda}_{ij}-\frac{\prod\limits_{k=j}^{i-1}{\Lambda}_{kk+1}}{\prod\limits_{k=j+1}^{i-1}{\Lambda}_{kk}}\right)(S_{j}-S_{j-1}),\
\ i=1,...,n.$

The value is

$\displaystyle\mathbb{E}_{\mathbb{P}}\left[-\exp\left(-V^{\hat{\gamma}}_{n}\right)\right]$
$\displaystyle=-\left(\frac{\prod\limits_{i=2}^{n-1}{\Lambda}_{ii}}{|\Sigma|\prod\limits_{i=1}^{n-1}\left({\Lambda}_{ii}{\Lambda}_{i+1i+1}-{\Lambda}^{2}_{ii+1}\right)}\right)^{1/2}\exp\left(-\frac{1}{2}\mu{\Lambda}\mu^{\prime}\right).$

###### Proof.

(I) Observe that for ${D}=0$ we have the following properties: The matrix
$\hat{Q}$ is diagonal with $\hat{Q}^{-1}_{ii}={\Lambda}_{ii}$ for all $i$; The
matrix $\hat{\Gamma}$ is given by $\hat{\Gamma}_{ij}={\Lambda}_{ij}$ for
$i\neq j$ and $\hat{\Gamma}_{ij}=0$ otherwise. This together with (1.3)–(1.4)
completes the proof.
(II) The matrix $\hat{Q}$ satisfies $\hat{Q}_{ij}=0$ for $|i-j|>1$ (i.e.
$\hat{Q}$ is tridiagonal). Hence, from (2.2) in [3]

$[\hat{Q}^{-1}]_{ij}=[\hat{Q}^{-1}]_{ji}=\frac{\prod\limits_{k=j}^{i-1}[\hat{Q}^{-1}]_{kk+1}}{\prod\limits_{k=j+1}^{i-1}[\hat{Q}^{-1}]_{kk}}\
\ \mbox{if}\ \ i>j+1$

and from (3.1) in [3]

$|\hat{Q}^{-1}|=\frac{\prod\limits_{i=1}^{n-1}\left([\hat{Q}^{-1}]_{ii}[\hat{Q}^{-1}]_{i+1i+1}-[\hat{Q}^{-1}]^{2}_{ii+1}\right)}{\prod\limits_{i=2}^{n-1}[\hat{Q}^{-1}]_{ii}}.$

Since $\hat{\Gamma}\in{\mathcal{T}}_{1}$ then

$\displaystyle[\hat{Q}^{-1}]_{ij}={\Lambda}_{ij}\ \ \mbox{if}\ \ |i-j|\leq 1.$

By combining the above equalities with (1.3)–(1.4) we complete the proof. ∎

######  Remark 1.4.

Let us emphasize that we do not assume any Markov structure for the underlying
asset, hence even for the “usual” case ${D}=0$ Theorem 1.1 gives an efficient
way to solve a non-Markovian optimization problem.

We end this section with the following example.

######  Example 1.5.

Consider a simple model where the increments have zero mean (i.e. $\mu=0$) and
the covariance matrix is the Kac-Murdock-Szegö matrix

$\Sigma_{ij}:=\rho^{|i-j|}\ \ \forall i,j$

for some constant $\rho\in(0,1)$. Assume that $n\geq 3$.

The matrix $\Sigma$ has a determinant equal to
$|\Sigma|=\left(1-\rho^{2}\right)^{n-1}$, and a simple tridiagonal inverse
that is given by (see [6])

${\Lambda}_{ij}=\begin{cases}0,&|i-j|>1\\\
-\frac{\rho}{1-\rho^{2}},&|i-j|=1\\\ \frac{1+\rho^{2}}{1-\rho^{2}},&1<i=j<n\\\
\frac{1}{1-\rho^{2}},&\mbox{otherwise}.\end{cases}$

Case I: ${D}=0$. From Corollary 1.3 we obtain that the optimal strategy is
$\hat{\gamma}_{1}=0$, and
$\hat{\gamma}_{i}=\frac{\rho}{1-\rho^{2}}(S_{i-1}-S_{i-2})$ for $i>1$. The
value is $-\sqrt{\frac{1-\rho^{2}}{(1+\rho^{2})^{n-2}}}$.
Case II: ${D}=1$. Corollary 1.3 gives the optimal strategy

$\hat{\gamma}_{1}=\hat{\gamma}_{2}=0\ \ \mbox{and}\ \ \
\hat{\gamma}_{i}=\frac{1+\rho^{2}}{1-\rho^{2}}\sum_{j=1}^{i-2}\left(-\frac{\rho}{1+\rho^{2}}\right)^{i-j}(S_{j}-S_{j-1})\
\ \mbox{for}\ \ i>2.$

The value is
$-\sqrt{\frac{\left(1-\rho^{2}\right)\left(1+\rho^{2}\right)^{n-2}}{\left(1+\rho^{2}+\rho^{4}\right)^{n-3}}}.$

##  2\. Proof of Theorem 1.1

In this section we prove Theorem 1.1.

###### Proof.

The proof will be done in three steps.

Step I: In this step we provide an existence and uniqueness for the
decomposition (1.2). We apply the results from [4] which perfectly fit our
purposes. First we introduce some notations. Let $C\in M_{n}(\mathbb{R})$ and
$I,J\in\\{1,..,n\\}^{k}$ be two integer vectors of length $k\leq n$ denoting
row and column indices, satisfying $I=(i_{1},..,i_{k}),J=(j_{1},..,j_{k})$
with $1\leq i_{1}<i_{2}<..<i_{k}\leq n$, $1\leq j_{1}<j_{2}<..<j_{k}\leq n$.
We denote by $C_{J}^{I}$ the minor specified by these indices, i.e. the
determinant of the $k\times k$ sub-matrix determined by taking the elements
$C_{ij}$ for $i=i_{1},..i_{k};j=j_{1},..j_{k}$. For any two natural numbers
$b\geq a$ we denote the vector $(a,a+1,...,b)$ by $[a:b]$.

Since ${\Lambda}$ is positive definite it follows by Sylvester’s criterion
(see e.g. [8]) that its principal minors
${\Lambda}^{[k+1:k+{D}]}_{[k+1:k+{D}]}>0$ for all $k\leq n-{D}-1$ and
${\Lambda}^{[k:k+{D}]}_{[k:k+{D}]}>0$ for all $k\leq n-{D}$. Hence, by
applying Theorem 5.5 in [4] for $r=s={D}+1$ we obtain that there is a unique
invertible matrix $R$ such that $R^{-1}\in{\mathcal{S}}_{{D}}$ and
${\Lambda}-R\in{\mathcal{T}}_{{D}}$. Set $\hat{Q}:=R^{-1}$,
$\hat{\Gamma}={\Lambda}-R$. Since ${\Sigma}$ is symmetric then from uniqueness
it follows that $R$ is also symmetric and the first step is completed.

Step II:  Denote by $\mathcal{Q}_{{D}}$ the set of all equivalent probability
measures $\mathbb{Q}\sim\mathbb{P}$ with finite entropy
$\mathbb{E}_{\mathbb{Q}}\left[\log\left(\frac{d\mathbb{Q}}{d\mathbb{P}}\right)\right]<\infty$
relative to $\mathbb{P}$ that satisfy

(2.1) $\mathbb{E}_{\mathbb{Q}}[S_{t}-S_{s}|\mathcal{G}^{{D}}_{s}]=0\ \ \forall
t\geq s.$

In this step we prove the following verification result: If a triplet
$(\tilde{\gamma},\tilde{\mathbb{Q}},C)\in\mathcal{A}_{{D}}\times\mathcal{Q}_{{D}}\times\mathbb{R}$
satisfies

(2.2)
$V^{\tilde{\gamma}}_{n}+\log\left(\frac{d{\tilde{\mathbb{Q}}}}{d\mathbb{P}}\right)=C$

then $\tilde{\gamma}\in\mathcal{A}_{{D}}$ is the unique optimal portfolio for
the optimization problem (1.1) and the corresponding value is
$\mathbb{E}_{\mathbb{P}}\left[-\exp\left(-V^{\tilde{\gamma}}_{n}\right)\right]=-e^{-C}.$

Indeed, using similar arguments as in the proof of Lemma 2.1 in [2], it
follows that for any $\gamma\in\mathcal{A}_{{D}}$ and
$\mathbb{Q}\in\mathcal{Q}_{{D}}$ we have
$\mathbb{E}_{\mathbb{Q}}[V^{\gamma}_{n}]=0$ provided that
$\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\gamma}_{n}\right)\right]<\infty$.
From (2.2)
$\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\tilde{\gamma}}_{n}\right)\right]=e^{-C}.$
Hence, $\mathbb{E}_{\tilde{\mathbb{Q}}}[V^{\tilde{\gamma}}_{n}]=0$ and by
applying (2.2) again

(2.3)
$\log\left(\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\tilde{\gamma}}_{n}\right)\right]\right)=-C=-\mathbb{E}_{\tilde{\mathbb{Q}}}\left[\log\left(\frac{d\tilde{\mathbb{Q}}}{d\mathbb{P}}\right)\right].$

Next, it is well known that for a strictly concave utility maximization
problem, the optimizer is unique (if exists). Thus, in view of (2.3), in order
to complete the proof of the second step it remains to show that in general we
have the inequality

(2.4)
$\log\left(\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\gamma}_{n}\right)\right]\right)\geq-\mathbb{E}_{\mathbb{Q}}\left[\log\left(\frac{d\mathbb{Q}}{d\mathbb{P}}\right)\right]\
\ \forall(\gamma,\mathbb{Q})\in\mathcal{A}_{{D}}\times\mathcal{Q}_{{D}}.$

Let us establish (2.4). Without loss of generality assume that
$\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\gamma}_{n}\right)\right]<\infty$
(otherwise (2.4) is trivial). Then for any $z\in\mathbb{R}$

$\displaystyle\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\gamma}_{n}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{P}}\left[\exp\left(-V^{\gamma}_{n}\right)+z\frac{d{\mathbb{Q}}}{d\mathbb{P}}V^{\gamma}_{n}\right]$
$\displaystyle\geq\mathbb{E}_{\mathbb{P}}\left[z\frac{d{\mathbb{Q}}}{d\mathbb{P}}\left(1-\log\left(z\frac{d{\mathbb{Q}}}{d\mathbb{P}}\right)\right)\right]$
$\displaystyle=z-z\log
z-z\mathbb{E}_{\mathbb{Q}}\left[\log\left(\frac{d{\mathbb{Q}}}{d\mathbb{P}}\right)\right].$

The first equality is due to $\mathbb{E}_{{\mathbb{Q}}}[V^{\gamma}_{n}]=0$.
The inequality follows from the Legendre-Fenchel duality inequality $xy\leq
e^{x}+y(\log y-1)$ for all $x,y\in\mathbb{R}$ by setting $x=-V^{\gamma}_{n}$
and $y=z\frac{d\mathbb{Q}}{d\mathbb{P}}$. The last equality is
straightforward. From simple calculus it follows that the concave function
$z\rightarrow z-z\log
z-z\mathbb{E}_{\mathbb{Q}}\left[\log\left(\frac{d{\mathbb{Q}}}{d\mathbb{P}}\right)\right]$,
$z>0$ attains its maximum at
$z^{*}:=\exp\left(-\mathbb{E}_{\mathbb{Q}}\left[\log\left(\frac{d{\mathbb{Q}}}{d\mathbb{P}}\right)\right]\right)$
and the corresponding maximal value is also $z^{*}$. This completes the proof
of (2.4).

Step III: In view of Step I the portfolio $\hat{\gamma}$ from (1.3) is well
defined. From the fact that $\hat{\Gamma}_{ij}=0$ for $|i-j|\leq{D}$ we obtain
that $\hat{\gamma}\in\mathcal{A}_{{D}}$. Set

$C:=\frac{1}{2}\left(\log|\Sigma|-\log|\hat{Q}|+\mu{\Lambda}\mu^{\prime}\right)$

and define the measure $\hat{\mathbb{Q}}$ by the Radon-Nikodym derivative

(2.5)
$\frac{d\hat{\mathbb{Q}}}{d\mathbb{P}}:=\exp\left(C-V^{\hat{\gamma}}_{n}\right).$

From Step II it follows that in order to complete the proof of Theorem 1.1 it
remains to establish that the measure $\hat{\mathbb{Q}}$ is a probability
measure which satisfies $\hat{\mathbb{Q}}\in\mathcal{Q}_{{D}}$.

To that end define $X=(X_{1},...,X_{n})$ by $X_{i}=S_{i}-S_{i-1}$,
$i=1,...,n$. From (1.3) we have (recall that $\hat{\Gamma}$ is symmetric)
$V^{\hat{\gamma}}_{n}=\hat{\gamma}X^{\prime}=\mu{\Lambda}X^{\prime}-\frac{1}{2}X\hat{\Gamma}X^{\prime}$.
This together with (1.2) and (2.5) gives

$\frac{d\hat{\mathbb{Q}}}{d\mathbb{P}}=\frac{\frac{\exp\left(-\frac{1}{2}X\hat{Q}^{-1}X^{\prime}\right)}{\sqrt{(2\pi)^{n}|\hat{Q}|}}}{\frac{\exp\left(-\frac{1}{2}(X-\mu){\Lambda}(X-\mu)^{\prime}\right)}{\sqrt{(2\pi)^{n}|\Sigma|}}}.$

Hence, the relation $\left(X;{\mathbb{P}}\right)\sim\mathcal{N}(\mu,\Sigma)$
implies that $\hat{\mathbb{Q}}$ is a probability measure and
$\left(X;\hat{\mathbb{Q}}\right)\sim\mathcal{N}(0,\hat{Q})$. In particular
$\mathbb{E}_{\hat{\mathbb{Q}}}\left[\log\left(\frac{d\hat{\mathbb{Q}}}{d\mathbb{P}}\right)\right]<\infty$.
Finally, from the property $\hat{Q}_{ij}=0$ for $|i-j|>{D}$ we obtain that
(under the probability measure $\hat{\mathbb{Q}}$) $X_{k}$ is independent of
$(X_{1},...,X_{(k-1-{D})^{+}})$ for all $k$. This yields (2.1) and completes
the proof. ∎

We end this section with the following remark.

######  Remark 2.1.

Observe that for any positive definite matrix $Q\in M_{n}(\mathbb{R})$, the
probability measure $\mathbb{Q}\sim\mathbb{P}$ which is given by

$\frac{d{\mathbb{Q}}}{d\mathbb{P}}=\frac{\frac{\exp\left(-\frac{1}{2}XQ^{-1}X^{\prime}\right)}{\sqrt{(2\pi)^{n}|Q|}}}{\frac{\exp\left(-\frac{1}{2}(X-\mu){\Lambda}(X-\mu)^{\prime}\right)}{\sqrt{(2\pi)^{n}|\Sigma|}}}$

satisfies $\left(X;{\mathbb{Q}}\right)\sim\mathcal{N}(0,Q)$. In particular for
a diagonal $Q$ the measure $\mathbb{Q}$ is an equivalent martingale measure,
and so the market is arbitrage-free in the classical sense (no delay).
Clearly, for any ${D}$ the set $\mathcal{Q}_{{D}}$ contains the set of all
equivalent martingale measures.

##  3\. Computational Results

We start this section by providing an explicit algorithm for the computation
of $\hat{Q}$ and $\hat{\Gamma}$ from (1.2). Recall the notations from Sections
1,2. In view of Corollary 3.2 in [4] and the fact that
$\hat{Q}\in{\mathcal{S}}_{{D}}$, we look for a matrix $R:=\hat{Q}^{-1}$ that
satisfies $R_{ij}={\Lambda}_{ij}$ for $i,j<{D}+1$ and has vanishing super and
sub ${D}+1$ minors, i.e. minors of the form $R^{I}_{J}$ where
$I=(i_{1},...,i_{{D}+1})$, $J=(j_{1},...,j_{{D}+1})$ that satisfy one of the
inequalities $j_{1}>i_{{D}+1}-{D}$ (super minors) or $i_{1}>j_{{D}+1}-{D}$
(sub minors). The matrix $R$ is symmetric hence we can restrict ourselves to
looking only at super minors. Since ${\Lambda}$ is positive definite and
$R_{ij}={\Lambda}_{ij}$ for $i,j<{D}+1$, it follows that the principal minors
$R^{[i+1:i+{D}]}_{[i+1:i+{D}]}>0$ for all $i\leq n-{D}-1$. Thus, from Theorem
3.4 in [4] we conclude that the matrix $R:=\hat{Q}^{-1}$ is the unique
symmetric matrix that satisfies $R_{ij}={\Lambda}_{ij}$ for $|i-j|\leq{D}$ and

(3.1) $R^{[i:i+{D}]}_{(i+1,i+2,...,i+{D},i+m)}=0,\ \ i\in[1:n-{D}-1],\ \
m\in[1+{D}:n-i].$

Observe that it is possible to order the $({D}+1)$-minors in (3.1) such that
each minor is a function of $({D}+1)^{2}-1$ known elements of $R$, and one
unknown element, hence providing a linear equation for this unknown element.
For example, the first minor $R^{[1:1+{D}]}_{[2:2+{D}]}$ is determined by
$(R_{ij}),i=1,..,{D}+1;j=2,..,{D}+2$ which are all set to the corresponding
elements of ${\Lambda}$ because $|i-j|\leq{D}$ except for $R_{1{D}+2}$. For
four integers $a\leq i\leq b\leq j$, denote by $[a:b]^{i\to j}$ the vector
obtained by omitting the index $i$ from $[a:b]$, and adding in the last
coordinate the index $j$, i.e. $[a:b]^{i\to j}=(a,a+1,..,i-1,i+1,..,b,j)$.
Hence, in the Laplace expansion for the $({D}+1)$-minor

$0=R^{[1:1+{D}]}_{[2:2+{D}]}=(-1)^{{D}}R_{12+{D}}R^{[2:1+{D}]}_{[2:1+{D}]}+\sum_{j=1}^{{D}}(-1)^{1+j}R_{11+j}R^{[2:1+{D}]}_{[2:2+{D}]^{1+j\to
2+{D}}}$

all the ${D}$-minors are known, and the unknown term $R_{12+{D}}$ appears only
in the first term. It is therefore possible to express $R_{12+{D}}$ uniquely
as a linear function of the known $({D}+1)^{2}-1$ $R_{ij}$ values.

Once $R_{1{D}+2}$ is determined, we can express the next $({D}+1)$-minor
$R^{[2:2+{D}]}_{[3:3+{D}]}$ in terms of $({D}+1)^{2}-1$ known $R_{ij}$
elements and the unknown element $R_{2{D}+3}$, and determine it in a similar
manner. Algorithm 1 describes how to proceed and determine all the unknown
elements $R_{ij}$ one by one, hence it provides the entire decomposition of
${\Lambda}$ into $R+\hat{\Gamma}$ (recall that $R=\hat{Q}^{-1}$).

Algorithm 1 Banded Matrix Decomposition

1:Input: $\Sigma\in M_{n}(\mathbb{R})$ \- a real Positive Definite matrix,
${D}\leq n-1$ \- band width.

2:Output: Two matrices $R+\hat{\Gamma}={\Sigma}^{-1}$ such that
$R^{-1}\in{\mathcal{S}}_{{D}},\hat{\Gamma}\in{\mathcal{T}}_{{D}}$.

3:Compute the inverse ${\Lambda}=\Sigma^{-1}$ and set
$R_{ij}={\Lambda}_{ij},\forall|i-j|\leq{D}$.

4:for $m={D}+1$ to $n-1$ do

5: for $i=1$ to $n-m$ do

6: Set
$R_{i+mi},R_{ii+m}=\frac{\sum\limits_{j=1}^{{D}}(-1)^{j+{D}}R_{ii+j}R_{[i+1:i+{D}]^{i+j\to
i+m}}^{[i+1:i+{D}]}}{R_{[i+1:i+{D}]}^{[i+1:i+{D}]}}.$

7:Set $\hat{\Gamma}={\Lambda}-R$.

8:Output the matrices $R,\hat{\Gamma}$.

######  Remark 3.1.

As matrix inversion can be performed in $O\left(n^{2.373}\right)$ floating
points operations [20], the total number of floating points operations for
Algorithm 1 is $O\left(n^{2.373}\right)$ (for matrix inversion), plus
$O\left(n^{2}\times{D}^{2.373}\right)$ for computing the minors (see [1]),
hence yielding an overall computational complexity of
$O\left(n^{2}(n^{0.373}+{D}^{2.373})\right)$ floating point operations. When
$\Sigma$ has a special structure, specialized faster algorithms can be used.
For example, if $\Sigma$ is a Toeplitz matrix, inversion can be performed in
$O\left(n^{2}\right)$ floating point operations using Trench’s algorithm [19],
hence the overall complexity is reduced to $O\left(n^{2}{D}^{2.373}\right)$.

Next, we implemented Algorithm 1 in $R$ code, freely available at github.
Using this implementation, we present some numerical results for the simple
model in Example 1.5, and for a more realistic model of a discretized
fractional Brownian motion. A fractional Brownian motion (fBm) [12] with Hurst
parameter ${\mathcal{H}}\in(0,1)$ is a centered Gaussian process
$B^{{\mathcal{H}}}_{t}$, $t\geq 0$ with the covariance function

(3.2)
$Cov\left(B^{\mathcal{H}}_{s},B^{\mathcal{H}}_{t}\right)=\frac{1}{2}\left(t^{2{\mathcal{H}}}+s^{2{\mathcal{H}}}-|t-s|^{2{\mathcal{H}}}\right).$

If ${\mathcal{H}}=\frac{1}{2}$ we recover the standard Brownian motion. Models
of asset prices based on fBm have long attracted the interest of researchers
for their properties of long-range dependence, see for instance, [11, 21]. For
some explicit computations for portfolio optimization in the fBm framework see
[9, 10].

We consider a discretized version of the fBm framework. Namely, we fix
${\mathcal{H}}\in(0,1)$, a time-resolution $\Delta t>0$ and a final time
$T=n\Delta t$ for some $n\in\mathbb{N}$, and consider the fBm market which is
active at times $0,\Delta t,2\Delta t,\dots,n\Delta t=T$. The stock price at
time $k\Delta t$ is given by $S_{k}=S_{0}+B^{{\mathcal{H}}}_{k\Delta t}$,
$k=0,1,\dots,n$. In this case $\mu=0$ and by (3.2) the covariance matrix is
the Toeplitz matrix given by

$\displaystyle\Sigma_{ij}$
$\displaystyle=Cov\left(B^{{\mathcal{H}}}_{(i+1)\Delta
t}-B^{{\mathcal{H}}}_{i\Delta t},B^{{\mathcal{H}}}_{(j+1)\Delta
t}-B^{{\mathcal{H}}}_{j\Delta t}\right)$ (3.3) $\displaystyle=\frac{\Delta
t^{2{\mathcal{H}}}}{2}\left(|i-j-1|^{2{\mathcal{H}}}+|i-j+1|^{2{\mathcal{H}}}-2|i-j|^{2{\mathcal{H}}}\right).$

By following carefully the proof of Proposition 1.6 in [14] we obtain that
$\Sigma$ is positive definite. Moreover, when $\Delta t$ is changed $\Sigma$
is multiplied by a constant scalar, hence for fixed values of $n$, the delay
${D}$ and the ${\mathcal{H}}$ parameter, the value is invariant to $\Delta t$.

We computed numerically the value as a function of the delay ${D}$ for both
models, for $n=64$ and $\mu=0$. The results, displayed in Fig. 1, show the
decrease in the value as ${D}$ in increasing for both models. We first
examined the Kac-Murdock-Szegö covariance matrix $\Sigma_{ij}:=\rho^{|i-j|}$
from Example 1.5. As expected, for any delay ${D}$ the value is monotonically
increasing with the correlations parameter $\rho$, hence longer delays can be
tolerated with only a negligible drop in the value as $\rho$ increases. For
example, when there is no delay the value is significantly reduced compared to
the maximum $0$ only for $\rho<0.3$, while for a delay of ${D}=4$ a
significant drop in the value is observed as soon as $\rho<0.9$.

Next, we examined the covariance matrix representing the discretized
fractional Brownian motion from (3.3). For high ${\mathcal{H}}$ values
(${\mathcal{H}}>0.5$), the increments are positively correlated hence the
value increases with ${D}$ towards the maximum $0$. For ${\mathcal{H}}=0.5$ we
reduce to the independent Brownian motion case, hence the value drops to $-1$
regardless of the delay. For low ${\mathcal{H}}$ values (${\mathcal{H}}<0.5$),
the increments are negatively correlated hence the value increases, but the
increase is non-monotonic in ${\mathcal{H}}$. The reason is that for very low
${\mathcal{H}}$ values the correlation between the increments decays faster to
$0$ with their time distance, hence a delay results in almost complete loss of
information regarding the current price.

The entire computations for each model at a resolution of $0.01$ of $\rho$ or
${\mathcal{H}}$ (in total computing the value for $100$ parameter values
$\times$ $8$ values of ${D}$, i.e. $800$ different $\Sigma$ matrices) took
$\sim\\!200$ seconds on a standard PC laptop in our implementation. While the
value is invariant to the total time $T$, changing the discretization
parameter $n$ may alter the value. Increasing $n$ can be interpreted as either
keeping a fixed time horizon $T$ and decreasing $\Delta t$, or as fixing
$\Delta t$ and increasing the horizon $T$. For a fixed delay ${D}$, the value
increases towards zero with $n$, this make sense since the continuous time fBm
model leads to arbitrage (see for instance [16, 5]). For example, as shown in
Figure 1, for ${\mathcal{H}}=0.2$ and $D=1$ the value is roughly $-0.426$. The
value increases to $-0.123$ and $-0.007$ for $n=128$ and $256$, respectively.
Increasing the delay ${D}_{n}$ with $n$ while keeping the ratio
$\frac{{D}_{n}}{n}$ fixed amounts to refining the discretization towards the
continuous limit with a fixed continuous delay $0<{D}<T$. The value appears to
converge to a negative number for this case (i.e. for
$\frac{{D}_{n}}{n}=\frac{1}{64}$), with $-0.283$ for $n=256$ and $-0.228$ for
$n=1024$. A more thorough theoretical and numeric investigation of the scaling
of the value as a function of the problem’s parameters for this and other
models remains for future work.

Figure 1. The value for different delays ${D}$ (shown in different colors) for
$n=64$, $\mu=0$. Top: The value is shown as a function of the $\rho$ parameter
(x-axis) for the Kac-Murdock-Szegö covariance matrix
$\Sigma_{ij}:=\rho^{|i-j|}$ from Example 1.5. The results for ${D}=0,1$ can be
deduced from the closed-form derivations shown in the example. Bottom: The
value is shown as a function of the Hurst index ${\mathcal{H}}$ (x-axis) for a
covariance matrix representing the discretized fractional Brownian motion from
(3.3).

## Acknowledgments

We thank the anonymous referees whose comments improved the quality of the
paper.

## References

  * Aho and Hopcroft [1974] A.V. Aho and J.E. Hopcroft.  _The design and analysis of computer algorithms_.  Pearson Education India, 1974.
  * Bank and Dolinsky [2021] P. Bank and Y. Dolinsky.  A note on utility indifference pricing with delayed information.  _SIAM Journal on Financial Mathematics_ , 12:SC–31–SC–43, 2021.
  * Barrett [1979] W.W. Barrett.  A theorem on inverse of tridiagonal matrices.  _Linear Algebra and its Applications_ , 26:211–217, 1979\.
  * Barrett and Feinsilver [1981] W.W. Barrett and P.J. Feinsilver.  Inverses of banded matrices.  _Linear Algebra and its Applications_ , 41:111–130, 1981\.
  * Cheredito [2003] P. Cheredito.  Arbitrage in fractional brownian motion models.  _Finance and Stochastics_ , 7(4):533–553, 2003\.
  * Dow [2003] M. Dow.  Explicit inverses of Toeplitz and associated matrices.  _ANZIAM_ , 44:185–215, 2003.
  * Frey [2000] R. Frey.  Risk minimization with incomplete information in a model for high-frequency data.  _Mathematical Finance_ , 10:215–225, 2000.
  * Gilbert [1991] G.T. Gilbert.  Positive definite matrices and Sylvester’s criterion.  _The American Mathematical Monthly_ , 98(1):44–46, 1991.
  * Guasoni et al. [2019] P. Guasoni, Z.Nika, and M.Rasonyi.  Trading fractional brownian motion.  _SIAM Journal on Financial Mathematics_ , 10:769–789, 2019\.
  * Guasoni et al. [2021] P. Guasoni, Y.Mishura, and M.Rasonyi.  High-frequency trading with fractional brownian motion.  _Finance and Stochastics_ , 25:277–310, 2021.
  * Lo [1991] A.W. Lo.  Long-term memory in stock market prices.  _Econometrica_ , 59:1279–1313, 1991.
  * Mandelbrot and Van Ness [1968] B.B. Mandelbrot and J.W. Van Ness.  Fractional brownian motions, fractional noises and applications.  _SIAM review_ , 10(4):422–437, 1968.
  * Mania et al. [2008] M. Mania, R. Tevzadze, and T. Toronjadze.  Mean–variance hedging under partial information.  _SIAM Journal on Control and Optimization_ , 47:2381–2409, 2008.
  * Nourdin [2012] I. Nourdin.  _Selected Aspects of Fractional Brownian Motion_.  Springer-Verlag: Italia, 2012.
  * Nunno and Øksendal [2009] G.D. Nunno and B. Øksendal.  Optimal portfolio, partial information and Malliavin calculus.  _Stochastic Processes and their Applications_ , 81:303–322, 2009.
  * Rogers [1997] L.C.G. Rogers.  Arbitrage with fractional brownian motion.  _Mathematical Finance_ , 7(1):95–105, 1997.
  * Saporito and Zhang [2019] Y.F. Saporito and J. Zhang.  Stochastic control with delayed information and related nonlinear master equation.  _SIAM Journal on Control and Optimization_ , 57:693–717, 2019.
  * Schweizer [1994] M. Schweizer.  Risk–minimizing hedging strategies under restricted information.  _Mathematical Finance_ , 4:327–342, 1994.
  * Trench [1964] W.F. Trench.  An algorithm for the inversion of finite Toeplitz matrices.  _Journal of the Society for Industrial and Applied Mathematics_ , 12(3):515–522, 1964.
  * Williams [2012] V.V. Williams.  Multiplying matrices faster than Coppersmith-Winograd.  In _Proceedings of the forty-fourth annual ACM symposium on Theory of computing_ , pages 887–898, 2012.
  * Willinger et al. [1999] W. Willinger, M.S. Taqqu, and V. Teverovsky.  Stock market prices and long–range dependence.  _Finance and Stochastics_ , 3:1–13, 1999.


# Studying Confirmation Bias in Hashtag Usage on Twitter

Dominik Kowald  Know-Center GmbHGraz, Austria<EMAIL_ADDRESS>and
Elisabeth Lex  Graz University of Technology & Know-Center GmbHGraz, Austria
<EMAIL_ADDRESS>
###### Abstract.

The micro-blogging platform Twitter allows its nearly 320 million monthly
active users to build a network of follower connections to other Twitter users
(i.e., followees) in order to subscribe to content posted by these users
(Myers and Leskovec, 2014; Kwak et al., 2010). With this feature, Twitter has
become one of the most popular social networks on the Web and was also the
first platform that offered the concept of hashtags (see post by Chris Messina
in 2007111https://twitter.com/chrismessina/status/223115412). Hashtags are
freely-chosen keywords, which start with the hash character “#”, to annotate,
categorize and contextualize Twitter posts (i.e., tweets) (Romero et al.,
2011; Huang et al., 2010).

Problem. Although hashtags are widely accepted and used by the Twitter
community, the heavy reuse of hashtags that are popular in the personal
Twitter networks (i.e., own hashtags and hashtags used by followees) can lead
to filter bubble effects (Pariser, 2011) and thus, to situations, in which
only content associated with these hashtags are presented to the user. These
filter bubble effects are also highly associated with the concept of
confirmation bias, which is the tendency to favor and reuse information that
confirms personal preferences (Plous, 1993). One example would be a Twitter
user who is interested in political tweets of US president Donald Trump.
Depending on the hashtags used, the user could either be stuck in a pro-Trump
(e.g., “#MAGA”) or contra-Trump (e.g., “#fakepresident”) filter bubble.

Therefore, the goal of this paper is to study confirmation bias and filter
bubble effects in hashtag usage on Twitter by treating the reuse of hashtags
as a phenomenon that fosters confirmation bias.

Method. To reach this goal, we crawl two datasets from Twitter. The first one
(i.e., CompSci dataset) consists of researchers from the field of computer
science (= seed users) and their followees, while the second one (i.e., Random
dataset) consists of random people (= seed users) and their followees. The
statistics of these datasets are depicted in Table 1 and details of the
crawling strategy are described in (Kowald et al., 2017). For all the seed
users (i.e., $|U_{S}|$) in these datasets, we analyze (i) individual hashtag
reuse (i.e., reusing own hashtags), and (ii) social hashtag reuse (i.e.,
reusing hashtags, which has been previously used by a followee) with respect
to hashtag usage types and temporal effects.

Finally, we discuss ways to overcome confirmation bias and filter bubble
effects in Twitter by means of recommendation mechanisms that focus on beyond-
accuracy metrics.

Code. For reasons of reproducibility, we conduct this study using our open-
source tag recommendation and benchmarking framework TagRec (Kowald et al.,
2017; Trattner et al., 2015), which is available via our GitHub
repository222https://github.com/learning-layers/TagRec.

Dataset | $|U_{S}|$ | $|U|$ | $|T|$ | $|HT|$ | $|HTAS|$
---|---|---|---|---|---
CompSci | 2,551 | 91,776 | 5,649,359 | 1,081,403 | 9,161,842
Random | 3,466 | 127,112 | 8,157,702 | 1,507,773 | 13,628,750
Table 1. Statistics of our CompSci and Random Twitter datasets. Here,
$|U_{S}|$ is the number of seed users, $|U|$ is the number of distinct users,
$|T|$ is the number of tweets, $|HT|$ is the number of distinct hashtags and
$|HTAS|$ is the number of hashtag assignments. We aim to study confirmation
bias in hashtag usage on Twitter by analyzing (i) the reuse rate of these
hashtag assignments, and (ii) temporal effects of hashtag usage.  Figure 1.
Analysis of hashtag reuse types as a proxy for studying confirmation bias. For
each hashtag assignment, we study whether the corresponding hashtag has been
used by the same user before in time (“individual”), by some of the followed
users (“social”), by both (“individual/social”), by anyone else (“network”) or
neither of them (“external”). We find that between 66% and 81% of hashtag
assignments in both datasets can be explained by individual or social hashtag
usage (i.e., the sum of “individual”, “social” and “individual/social”). This
is a clear indication of confirmation bias in hashtag usage on Twitter.

Results. In our datasets, we first analyze hashtag assignments as well as
hashtag reuse practices with the aim of studying confirmation bias in hashtag
usage on Twitter. Specifically, for each hashtag assignment, we study whether
the corresponding hashtag has either been used by the same user before
(“individual”), by some of her followees (“social”), by both
(“individual/social”), by anyone else in the dataset (“network”) or by neither
of them (“external”). The results of this study are shown in Figure 1. We find
that 66% of hashtag assignments in the CompSci dataset and 81% in the Random
dataset can be explained by individual or social hashtag reuse, which is a
clear indication of confirmation bias and filter bubble effects.

(a)  (b)  (c)  (d)

Figure 2. Temporal effects on individual and social hashtag reuse for the
CompSci and Random datasets (plots are in log-log scale). Plots (a) and (b)
show that the more recently a hashtag $ht$ was used by a user $u$, the higher
its individual reuse count (i.e., people tend to reuse hashtags that have been
used very recently by their own). Plots (c) and (d) show that the more
recently a user $u$ was exposed to a hashtag $ht$, which was used by the
followees $F_{u}$, the higher its social reuse count (i.e., people tend to
reuse hashtags that have been used recently in the social network). We find
that temporal effects play an important role in individual and social hashtag
reuse, which is our measure for studying the confirmation bias in hashtag
usage on Twitter.

Next, we study to what extent temporal effects play a role in the reuse of
individual and social hashtags on Twitter. Specifically, we analyze the
recency of hashtags assignments (i.e., the time since the last hashtag usage).
The results for this study are shown in the plots of Figure 2 and we find that
temporal effects have an important influence on both individual as well as
social hashtag reuse. People tend to reuse hashtags that were used very
recently by their own and/or by their Twitter followees. Interestingly, there
is a clear peak after 24 hours in all four plots, which further indicates that
users typically (re)use the same set of hashtags in this time span and thus,
tend to tweet about similar topics on a daily basis. This shows that also
temporal effects are an important influence for confirmation bias and
filtering bubble effects in hashtag usage on Twitter.

Conclusion & Discussion. Our findings of this paper are two-fold:

  1. (1)

We find that between 66% and 81% of hashtag assignments in our datasets can be
explained by individual or social hashtag reuse, which is a clear indication
of confirmation bias in hashtag usage on Twitter.

  2. (2)

We find that people tend to reuse hashtags that were used very recently by
their own and/or by their Twitter followees, which means that confirmation
bias in hashtag usage is also strongly influenced by temporal effects.

Based on this, we now discuss ways to overcome confirmation bias and filter
bubble effects by using recommendation mechanisms. In our previous work (i.e.,
(Kowald et al., 2017)), we presented a novel hashtag recommendation algorithm,
which utilizes the Base-Level Learning (BLL) equation of the cognitive
architecture ACT-R (Anderson et al., 2004) for suggesting hashtags to users.
Specifically, this algorithm ranks the hashtags of a user and her followees
based on the factors of frequency and recency (see also (Kowald and Lex, 2016;
Kowald et al., 2015)). Although this approach achieves good results with
respect to recommendation accuracy (see BLLI,S in Figure 3), it fosters
hashtag reuse and thus, confirmation bias and filter bubbles effects in our
two Twitter datasets.

(a) CompSci dataset (b) Random dataset

Figure 3. Precision / Recall plots showing the accuracy of BLLI (i.e., only
individual hashtags), BLLS (i.e., only social hashtags), CF (i.e.,
Collaborative Filtering), BLLI,S (i.e., individual and social hashtags) for
$k$ = 1 - 10 recommended hashtags. BLLI,S provides the best results with
respect to predication accuracy but also fosters the reuse of hashtags and
thus, confirmation bias and filter bubbles effects in our two Twitter
datasets.

In order to overcome these effects, we propose to focus on beyond-accuracy
metrics of recommender systems such as diversity and serendipity (Kaminskas
and Bridge, 2017; Kowald and Lex, 2015). If we again take the example of a
Twitter user who is interested in political tweets of US president Donald
Trump, we could think of a hybrid recommendation approach with the goal to
achieve a trade-off between accuracy, serendipity and diversity. For example,
if we think of a pro-Trump Twitter user, a typical accuracy-focused
recommender system would only recommend pro-Trump tweets and hashtags (e.g.,
“#MAGA”). Thus, our idea is to mix these pro-Trump recommendations also with
contra-Trump recommendations (e.g., “#fakepresident”) to boost serendipity and
diversity, and to give the user the chance to break the filter bubble.

Future Work. Finally, for future work, we also plan to verify our findings in
larger Twitter data samples than the ones used in this paper as well as in
other online social networks that feature the concept of hashtags, such as
Instagram and Facebook.

Additionally, it would be interesting to evaluate the influence of social
hashtags exposure, for example, by investigating reweet or mention networks in
Twitter.

Acknowledgments. This work is funded by the Know-Center GmbH, which is part of
the Austrian COMET program, and the H2020 project AFEL under grant agreement
687916.

Keywords. Twitter; Hashtag Usage; Confirmation Bias; Filter Bubble; Tag
Recommendations; ACT-R; TagRec

## References

  * (1)
  * Anderson et al. (2004) John R Anderson, Daniel Bothell, Michael D Byrne, Scott Douglass, Christian Lebiere, and Yulin Qin. 2004.  An integrated theory of the mind.  Psychological review 111, 4 (2004), 1036\.
  * Huang et al. (2010) Jeff Huang, Katherine M. Thornton, and Efthimis N. Efthimiadis. 2010\.  Conversational Tagging in Twitter. In Proc. of HT ’10. ACM.
  * Kaminskas and Bridge (2017) Marius Kaminskas and Derek Bridge. 2017.  Diversity, serendipity, novelty, and coverage: a survey and empirical analysis of beyond-accuracy objectives in recommender systems.  ACM Transactions on Interactive Intelligent Systems (TiiS) 7, 1 (2017), 2.
  * Kowald et al. (2017) Dominik Kowald, Simone Kopeinik, and Elisabeth Lex. 2017\.  The tagrec framework as a toolkit for the development of tag-based recommender systems. In Proc. of UMAP’17. ACM, 23–28.
  * Kowald et al. (2015) Dominik Kowald, Simone Kopeinik, Paul Seitlinger, Tobias Ley, Dietrich Albert, and Christoph Trattner. 2015.  Refining Frequency-Based Tag Reuse Predictions by Means of Time and Semantic Context.  In Mining, Modeling, and Recommending’Things’ in Social Media. Springer, 55–74.
  * Kowald and Lex (2015) Dominik Kowald and Elisabeth Lex. 2015.  Evaluating Tag Recommender Algorithms in Real-World Folksonomies: A Comparative Study. In Proc. of RecSys’15. ACM, New York, NY, USA, 265–268.
  * Kowald and Lex (2016) Dominik Kowald and Elisabeth Lex. 2016.  The Influence of Frequency, Recency and Semantic Context on the Reuse of Tags in Social Tagging Systems. In Proc. of HT’16. ACM, 237–242.  DOI:http://dx.doi.org/10.1145/2914586.2914617
  * Kowald et al. (2017) Dominik Kowald, Subhash Chandra Pujari, and Elisabeth Lex. 2017.  Temporal effects on hashtag reuse in twitter: A cognitive-inspired hashtag recommendation approach. In Proc. of WWW’17. International World Wide Web Conferences Steering Committee, 1401–1410.
  * Kwak et al. (2010) Haewoon Kwak, Changhyun Lee, Hosung Park, and Sue Moon. 2010\.  What is Twitter, a Social Network or a News Media?. In Proc. of WWW ’10. ACM, New York, NY, USA, 591–600.  DOI:http://dx.doi.org/10.1145/1772690.1772751
  * Myers and Leskovec (2014) Seth A. Myers and Jure Leskovec. 2014.  The Bursty Dynamics of the Twitter Information Network. In Proc. of WWW’14. ACM, New York, NY, USA, 913–924.  DOI:http://dx.doi.org/10.1145/2566486.2568043
  * Pariser (2011) Eli Pariser. 2011\.  The Filter Bubble: What the Internet Is Hiding from You.  (2011).
  * Plous (1993) Scott Plous. 1993\.  The psychology of judgment and decision making. Mcgraw-Hill Book Company.
  * Romero et al. (2011) Daniel M. Romero, Brendan Meeder, and Jon Kleinberg. 2011\.  Differences in the Mechanics of Information Diffusion Across Topics: Idioms, Political Hashtags, and Complex Contagion on Twitter. In Proc. of WWW ’11. ACM, New York, NY, USA, 695–704.  DOI:http://dx.doi.org/10.1145/1963405.1963503
  * Trattner et al. (2015) Christoph Trattner, Dominik Kowald, and Emanuel Lacic. 2015\.  TagRec: towards a toolkit for reproducible evaluation and development of tag-based recommender algorithms.  ACM SIGWEB Newsletter Winter (2015), 3.