# Decision Making for Human-in-the-loop Robotic Agents via
Uncertainty-Aware Reinforcement Learning

Siddharth Singi12, Zhanpeng He13, Alvin Pan3, Sandip Patel2, Gunnar A.
Sigurdsson4,
Robinson Piramuthu4, Shuran Song3 and Matei Ciocarlie2
https://roamlab.github.io/hula/ 1 denotes joint first authorship.2Dept. of
Mechanical Engineering, 3Dept. of Computer Science, Columbia University, New
York, USA4Amazon Alexa AIThis work was supported in part by an award from the
Columbia Center of Artificial Intelligence Technology.

###### Abstract

In a Human-in-the-Loop paradigm, a robotic agent is able to act mostly
autonomously in solving a task, but can request help from an external expert
when needed. However, knowing when to request such assistance is critical: too
few requests can lead to the robot making mistakes, but too many requests can
overload the expert. In this paper, we present a Reinforcement Learning based
approach to this problem, where a semi-autonomous agent asks for external
assistance when it has low confidence in the eventual success of the task. The
confidence level is computed by estimating the variance of the return from the
current state. We show that this estimate can be iteratively improved during
training using a Bellman-like recursion. On discrete navigation problems with
both fully- and partially-observable state information, we show that our
method makes effective use of a limited budget of expert calls at run-time,
despite having no access to the expert at training time.

##  I Introduction

Deep Reinforcement Learning (DRL) has shown great progress in learning
decision-making for complex robotic skills [1, 2, 3, 4] using experiences
collected by a robotic agent exploring an environment and receiving reward
signals. Traditional RL agents use a policy learned during training in order
to act autonomously at deployment. However, even a well-trained agent can
encounter situations when deployed that are hard to make decisions for, for
reasons such as partial state observability, uncertain dynamics, changes in
state distributions between training and testing, etc.

The Human-in-the-Loop (HitL) paradigm has been developed in robotics precisely
for situations where an agent can act autonomously most of the time, but would
still benefit from receiving assistance from an available (tele-)operator,
usually assumed to be a human expert. This paradigm is particularly powerful
if the agent itself makes the decision of when to request assistance, thus
freeing the operator from having to monitor task progress. However, this
approach gives rise to a critical decision-making problem on the agent’s part:
when to request help? Too few such requests can lead to the robot making
mistakes, but too many requests will overload the expert and lose the benefit
of semi-autonomous operation.

In this paper, we propose a DRL-based method for a HitL agent to make the
critical determination of when to request expert assistance. We posit that the
best moment to request such assistance is when the agent is highly uncertain
in the successful outcome of the task. From an RL perspective, we relate this
uncertainty to the variance of the return from the current state, as perceived
by the agent. Numerous RL algorithms provide methods for the agent to estimate
the expected return from a given state. We show that similar methods can be
used to estimate the variance of return as well during the training process.
At deployment time, the agent can then request expert assistance when its
estimate of the return variance from the current state falls below a given
threshold. We dub our method HULA, for Human-in-the-loop Uncertainty-aware
Learning Agent, and illustrate its operation in Fig. 1.

Figure 1: An illustration of HULA, the method we propose in this paper. In a
partially-observable environment, an agent without the help of an expert (A)
cannot localize itself accurately due to partial observability, goes down the
wrong passage and fails to reach the target. A HULA agent (B) decides to
request assistance from an available external expert in the states marked with
a red E and achieves the goal (denoted by a green star $\bigstar$). In both
cases, the agent can only observe a 5x5 grid around its current location;
shadow areas represent unobserved regions throughout the navigation.

Critically, HULA does not need to make any calls to the expert during
training. This stands in contrast to a standard RL approach, where an agent
could learn how to use an expert simply by making numerous such calls at train
time. Nevertheless, our method is able to make effective use of a limited
budget of expert calls in order to improve task performance. We summarize our
key contributions as follows:

  * •

To the best of our knowledge, we are the first to propose a method for an HitL
RL agent to learn how to effectively budget its interactions with a human
expert at deployment time, without needing any expert calls during training.

  * •

We show that the variance of return, which can be estimated during training
using Bellman-like equations, is an effective measure for agent uncertainty
and can be used at deployment to determine when to request assistance.

  * •

Our experiments on discrete navigation problems with both fully- and
partially-observable state information show that our method is as effective or
better in managing its budget of expert calls compared to a standard learning
approach that also makes thousands of expert calls during training.

##  II Related Work

Our work extends standard reinforcement learning algorithms to human-in-the-
loop policy that solves robotics tasks with the help of humans. Researchers
have investigated how to learn such a policy. For example, Arakawa et al.
propose DQN-TAMER [5], which is a method that incorporates a human observer
model by real-time human feedback during training. Expected Local Improvement
(ELI) [6] trains a state selector that suggests states requiring to query new
actions from human experts. PAINT [7] learns a classifier to identify
irreversible states and query the expert in case entering such states cannot
be avoided. Hug-DRL [8] leverages a control transfer mechanism between humans
and automation that corrects the agent’s unreasonable actions during training
by human intervention, where it has demonstrated significant potential in
autonomous driving applications. While these methods exploit human
intervention during training, our work relies on the internal uncertainty of
an agent and only uses an expert in test time.

HULA explicitly estimates the uncertainty of an RL agent and use it to make
decisions about requesting an expert’s assistance. Prior work has explored
representing the uncertainty of outcomes during RL training [9, 10, 11]. For
instance, Kahn et al. use the variance of the dynamics model to represent the
uncertainty of collision to avoid damage from high-speed collisions [12].
Ensemble Quantile Networks (EQN) [13] estimates both aleatoric and epistemic
uncertainties via the combination of quantile networks(IQN) [14] and
randomized prior function(RPF) [15]. These works use the uncertainty of an
agent to perform task-specific decision-making, while our work differs in
using it to collaborate with a human. Moreover, our work represents an agent’s
uncertainty using the variance of returns, which is a direct indicator of task
performance.

The closest work to ours is RCMP, which learns a HitL policy that queries the
expert when the epistemic uncertainty estimated by the variance of multiple
value functions is high [16]. However, RCMP requires expert queries during
training, whereas our work learns an uncertainty model based on the variance
of return and does not query an expert in train time.

##  III Method

Consider a problem for an autonomous agent formulated as a standard Markov
decision process (MDP). A general MDP is defined as a tuple of $(S,A,r,p)$,
where $S$ represents the state space, $A$ represents the action space,
$r(s,a)$ is a reward function that evaluates immediate reward of action $a\in
A$ in a state $s\in S$, and $p$ represent a transition distribution
$p(s_{t+1}|s_{t},a_{t})$. The goal of solving an MDP is to learn a policy
$\pi(a|s)$ that maximizes its expected return $\mathbb{E}_{\pi}[R]$, where
$R=\sum\gamma^{t}r(s_{t},a_{t})$ and $\gamma$ is a discount factor.

Going beyond the standard MDP formulation, we now assume the availability of
an expert that can give instructions to the robotic agent. When queried from a
given state $s_{t}$, this expert can directly provide an action
$a_{exp}(s_{t})$ for the robot to execute. We assume that the expert always
provides high-quality advice, i.e. a policy of always following expert-
provided actions from all states will achieve a satisfactorily high return on
the MDP problem above. However, such a policy would be impractical: we assume
that the expert has limited bandwidth or availability, thus it is desirable
for the agent to balance the goal of achieving high returns with the goal of
not overloading the expert with requests. Such a scenario is typical in HitL
robotics applications.

Our method aims to determine when the agent should request expert assistance
while making the most of a limited number of such calls available during
deployment. Furthermore, we will like to limit (or eliminate) the number of
calls made to the expert during its training phase.

To achieve this goal, our method leverages the uncertainty of the outcomes of
the agent during exploration. Specifically, we use the variance of the return
from a given state as a measure of uncertainty. Intuitively, a state with high
variance of the return is one from which the agent, acting alone, could
achieve a range of outcomes, ranging from very poor to very effective. We
posit that these are the best states in which to request assistance. In
contrast, states where the variance is low are those in which the agent is
certain of the outcome of the task, and there are fewer ways in which external
assistance can be of help.

Most RL algorithms work by estimating the expected return from a given state;
this estimate is continuously refined during training time. We modify this
training procedure also to maintain and refine an estimate of
$\texttt{Var}_{\pi}(R)$, the variance of the returns under the learned policy,
as an indicator of the agent’s uncertainty. At test time, our agent can use
this estimate of variance to decide when to request help. We detail these
processes next.

###  III-A Estimation for the variance of return

The expected return from a given state is encapsulated by the state value
functions, which are explicitly relied on by most RL algorithms. However,
while most RL algorithms are not concerned with the variance of the return, we
would like to compute and maintain similar estimates for it during training.
By definition, the variance of the return is defined by:

$\texttt{Var}(R)=\mathbb{E}[R^{2}]-\mathbb{E}[R]^{2}$

The expected return of taking action $a_{t}$ in state $s_{t}$ is given by the
Q-function under policy $\pi$:

$\displaystyle\mathbb{E}[R]=Q^{\pi}(s_{t},a_{t})$ (1)

Since the transition function is not available, it is usually approximated via
Bellman update:

$\displaystyle Q^{\pi}(s_{t},a_{t})\leftarrow(1-\alpha)$ $\displaystyle
Q^{\pi}(s_{t},a_{t})+\alpha(r(s_{t},a_{t})$
$\displaystyle+\gamma\max_{a_{t+1}\in A}Q(a_{t+1},s_{t+1}))$ (2)

The second moment of return $\mathbb{E}[R^{2}]$, which we call
$M^{\pi}(s_{t},a_{t})$, can be approximated using a sampling-based method by:

$\displaystyle M^{\pi}$ $\displaystyle(s_{t},a_{t})=\mathbb{E}[R^{2}]$
$\displaystyle=\mathbb{E}[[r(s_{t},a_{t})+\gamma\sum_{k=t+1}^{k=N}r(s_{k},a_{k})]^{2}]$
$\displaystyle=\mathbb{E}[r^{2}(s_{t},a_{t})+2\gamma\sum_{k=t}^{k=N}r(s_{k},a_{k})+(\sum_{k=t+1}^{k=N}r(s_{k},a_{k}))^{2}]$
$\displaystyle=r^{2}(s_{t},a_{t})+2\gamma\sum_{s_{t+1}\in
S}p(s_{t+1}|s_{t},a_{t})Q^{\pi}(s_{t+1},a_{t+1})$ $\displaystyle\hskip
10.00002pt+\gamma^{2}\sum_{s_{t+1}\in
S}p(s_{t+1}|s_{t},a_{t})M^{\pi}(s_{t+1},a_{t+1})$ (3)

Since the transition function is not available for a model-free agent, we can
estimate $M^{\pi}$ using a Bellman-like formula with sampled data:

$\displaystyle\vspace{2em}M^{\pi}(s_{t},a_{t})$
$\displaystyle\leftarrow(1-\alpha)M^{\pi}(s_{t},a_{t})+\alpha
M^{\prime}(s_{t},a_{t})$ (4) where $\displaystyle M^{\prime}(s_{t},a_{t})$
$\displaystyle=r^{2}(s_{t},a_{t})+2\gamma Q^{\pi}(s_{t+1},a_{t+1})$
$\displaystyle+\gamma^{2}M(s_{t+1},a_{t+1})$ (5)

Here, $a_{t+1}=\operatorname*{argmax}_{a_{t+1}\in A}Q(s_{t+1},a_{t+1})$, since
we assume a greedy agent that will take the actions that lead to max expected
return during deployment.

Finally, combining $M$ and $Q$, we can approximate the variance of returns in
a state by:

$\displaystyle\texttt{Var}(R)=M(s_{t},a_{t})-Q^{2}(s_{t},a_{t})$ (6)

Hence, the variance of the return given a state and action pair only depends
on the $Q$ function as well as the $M$ (second function), which in turn
depends on $Q$ and can be updated via a Bellman-like recursion.

###  III-B HULA: Complete method

Algorithm 1 HULA Training

1: while returns have not converged do

2: Collect $H$ trajectories using the current control policy $\pi$ with an MDP

3: Update estimator of the Q-function $Q_{\theta_{1}}$ using Eq. (2)

4: Update estimator of the M-function $M_{\theta_{2}}$ using Eq. (5)

5: With updated $M$ and $Q$, compute state-dependent variance using Eq. (6)

6: end while

We are now ready to integrate the variance estimation into a complete learning
algorithm. We use algorithms from the Q-learning family, as these already
build an explicit estimate of the $Q$ function at train time, and also use a
greedy policy w.r.t. $Q$ at run-time.

During training, in addition to the function estimator for the $Q$ function
(referred to as $Q_{\theta_{1}}$), we build and update an estimator for $M$
(referred to as $Q_{\theta_{2}}$). This estimator can be of the same type as
used for $Q$: for tabular Q-learning, we can use a table, while for Deep
Q-Networks (DQN) [17] we use a deep neural network. At every iteration, we use
Eqs. (4-5) to update the estimator for $M$. In the case of tabular Q-learning,
Eq. (4) can be used directly, while in the case of DQN it can be transformed
in a loss function analogous to the one used for $Q$. The complete procedure
is shown in Alg. 1.

Finally, during deployment, the agent can use the trained estimators for the
$Q$ and $M$ functions to compute the variance of the return for the current
state. If this variance exceeds a threshold, the agent will request assistance
and execute the action provided by the expert. Otherwise, the agent will
execute the action prescribed by its own policy (in this case, a greedy
selection based on the trained $Q$-function).

We note a key feature of our method: it does not require the presence of an
expert during train time. The only change at train time compared to
traditional, expert-less Q-learning is the addition of the function estimator
for $M$, which will be used at run-time to help provide a variance estimator.
In turn, the variance estimator is used to decide when to request assistance
from the expert, which is only needed at deployment.

Algorithm 2 HULA Deployment

1: while not done do

2: Compute variance $v_{t}$ in current state using Eq. (6)

3: if $v_{t}\geq\epsilon$ then

4: Execute action $a_{exp}(s_{t})$ provided by expert

5: else

6: Execute action from own policy $\pi(a_{t}|s_{t})$

7: end if

8: end while

##  IV Evaluation

Figure 2: Fully-observable grid worlds. The trap world (A) requires the agent
to navigate around traps. In the shortcut world (B), the agent can leverage
the expert to use the shortcut alley to achieve higher returns. Figure 3:
Partially-observable grid world. In this environment, an agent may not
localize itself by observing a surrounding 5x5 region. For example, the blue
state and the red state result in identical observations.

###  IV-A Environments

We evaluate HULA in a discrete navigation scenario, where the agent must reach
a goal while avoiding obstacles and/or traps. Furthermore, we test our
approach both on problems where the agent has exact knowledge of its current
state (which is tractable via tabular Q-learning) and problems where the agent
only has access to sensor data of limited range, creating ambiguity (which
requires more powerful function estimators, and the use of DQN).

####  IV-A1 Fully-observable MDPs

These are discrete navigation environments where the agent is provided its
exact location as observation. In addition to the goal cell, the environment
also contains obstacles and traps: colliding with a trap would terminate the
episode, while colliding with an obstacle results in a failure of action, and
the agent would stay in its original cell. The agent also receives a small
penalty for each step it takes. Available actions for an agent are moving up,
down, left, and right, and the observations are the coordinates $(x,y)$ of the
agent’s location. We use the maps shown in Fig. 2, and train and test
individually on each map.

For this simple class of problems, an unsophisticated agent can easily learn
to act optimally on its own. However, we introduce uncertainty in the form of
a stochastic transition function: at every step, after selecting an action,
the agent moves in the desired direction with probability $\psi$, and moves in
a random direction with probability $1-\psi$ (in practice, we use
$\psi=0.45$). We think of this setting as a ”slippery world”. Thus, even with
an optimal policy, the robot may not reach its goal.

To help, we introduce an expert: when called, the expert provides the optimal
action in the direction of the goal, and the action provided by the expert is
not subject to transition function stochasticity. In this setting, the expert
can thus be thought of as an “action corrector” with a better understanding of
the environment dynamics, and whose actions always produce the expected
results. Intuitively, we would expect an agent to call the expert in high-risk
situations where it needs to avoid a wrong move, such as very close to traps.
In this experiment, we test HULA on two maps: 1. trap world that requires an
agent to navigate around traps with an expert; 2. shortcut grid where the
agent can leverage the expert to take the shortcut for higher returns.

####  IV-A2 Partially-observable MDP’s

To test how our algorithm performs in more complex environments, we extend to
a case of partially observable MDP’s where the agent does not have complete
information about the system. Instead of observing state information as the
$(x,y)$ position in the grid, the observation now only includes a finite patch
of grid cells around the agent’s current location (we use 5x5 grids in our
implementation).

As shown in Figure 3, due to partial observability, it is not always possible
for the agent to uniquely identify its own location in the map from
observations. In such ambiguous regions, the same action taken based on
identical observations can lead to utterly different results. However, other
parts of the map are uniquely identifiable based on sensor observations, and
thus an autonomous agent can always select the optimal action.

In this experiment, the expert observes the full state and always provides
optimal actions to the agent. Intuitively, we would expect the agent to make
use of the expert when traversing ambiguous areas of the map.

Figure 4: A learned variance map for the trap world environment. (A) shows the
return variance computed by Monte-Carlo Sampling with the same policy. (B)
shows the estimated variance by our method (darker colors represent higher
variances).

###  IV-B Evaluation Approach

The main evaluation metric for our approach focuses on its ability to make
effective use of the expert: we would like to see how performance (measured as
average episodic return) changes as a function of the number of expert calls
made during deployment. We can measure this performance by varying the value
of the variance threshold $\epsilon$ used in Alg. 2: for large values of
$\epsilon$, the agent will never make use of the expert; conversely, if
$\epsilon$ is very small, the agent will call the expert in every state. By
sweeping the value of $\epsilon$ between these extremes, we obtain more or
less autonomous agents, and can plot performance as a function of the number
of expert calls that result in each case.

As a baseline, we use a standard RL approach that simply integrates the expert
into the training procedure. We refer to this baseline as the penalty-based
agent. Specifically, the penalty-based agent treats calling expert as an extra
action $a_{call}$ that can be called alongside the other actions at both
training and deployment time.

However, since the expert is optimal, the penalty-based agent will learn to
always call for help. To learn a nontrivial policy that does not call the
expert excessively, we employ a penalty of calling an expert in the reward
function at train time: $r^{\prime}(s,a)=r(s,a)+c$, where $c<0$ is a penalty
assessed only when calling an expert. By training this method with varying
values of $c$, we again obtain agents that are more or less autonomous: the
penalty-based agents trained with high $c$ will call the expert less often,
while those trained will low $c$ will call more often. Again, we sweep the
value of $c$ between these extremes, and plot performance as a function of the
number of expert calls that result in each case. However, for a fair
comparison against HULA, we do not assess the expert penalty at deployment
time.

##  V Results

Figure 5: Comparisons between a variance map learned by HULA and the ground
truth for the shortcut world. Figure 6: Comparisons between a variance map
learned by HULA and the ground truth for the partially-observable environment.

###  V-A Uncertainty Estimation

TABLE I: Quantitative results of the expert call efficiency. Accuracy of top-$N$ variance states estimation. PO represents the partially-observable experiments. Experiment | Trap world | Shortcut world | PO
---|---|---|---
Top-5 | 0.8 | 0.2 | 0.8
Top-10 | 0.8 | 0.6 | 0.7

The first question we would like to answer is whether HULA successfully
captures the uncertainty of the agent and allocate the expert calls
efficiently. We compare the learned variances by our method with the ground
truth variance of returns by Monte Carlo sampling the environment with the
same policy without the assistance of an expert.

As shown in Figure 4, in the fully-observable experiment, HULA successfully
estimates a variance map that is similar to the ground truth. Specifically, it
is uncertain about its outcome in states close to traps and confident when it
is far from traps.

In the partially observable case, HULA identifies ambiguous states when the
agent cannot localize itself in the map. In these states, the agent fails to
produce an action faithfully in that the same action taken can lead to
different results (e.g., hitting the boundary or reaching the goal). Our
results align with the ground-truth variance map, which also has high
uncertainty in a similar region.

To further investigate HULA’s efficiency in expert allocation, we compare the
states with the highest $N$ variances with those in the ground truth variance
map. These states are usually where the agent calls the expert, especially
when the agent has a budget of $N$ expert calls.

Our results show that HULA can estimate the uncertainty of an RL agent
accurately and allocate its expert call based on the return variance
efficiently. As shown in Table I, in trap world, $80\%$ of the top-$5$ and
top-$10$ variance states estimated by HULA are also the top-$5$ and $10$
accordingly in the ground-truth variance. In the shortcut world, although our
method only correctly estimates $1$ of $5$ highest variance states, its
performance is higher and reaches $60\%$ in the top-$10$ case. This indicates
that it allocates expert calls efficiently if given more budget. In the
partially-observable experiment, our agent also recovers most of the high
variance states in both the top-$5$ and top-$10$ evaluation.

Figure 7: # of expert calls vs Average returns for the fully-observable
experiments. Curves are generated with the rolling mean method (window size of
4) for visualization purposes.

###  V-B Task Performance

In all experiments, our results show that HULA agents achieve higher
performance when requesting help from experts compared to the ones that does
not use an expert. In the trap world, as shown in Fig. 7, HULA achieves
similar performance as the penalty-based agent while using the same number of
expert calls. Interestingly, in the shortcut world, we find that the penalty-
based method has a higher average return than HULA. This is because our method
does not evaluate the effect of calling an expert during training. As shown in
Fig. 8, the penalty-based agent learns to move towards states that ask for
help from an expert since it knows it will call an expert in a future state,
whereas HULA navigates around the uncertain states to avoid them. However, in
the highest-variance states (e.g., states between traps), both agents are able
to call the expert and complete the task.

In this experiment, we find that HULA is robust to different expert call
budgets compared to the penalty-based agent. For example, in the shortcut
world, the penalty-based agent cannot learn a policy that uses 2 expert calls
to complete the task even if we perform hyper-parameter sweeping with the
expert penalty $c$. This can be caused by the high expert penalty that
modifies the original reward structure and hinders the learning of the task.
In contrast, our method can flexibly incorporate different expert call budgets
and maintain reasonable performance.

In the partially-observed grid world, the HULA agent outperforms the penalty-
based agent when the allowed expert call is limited in the range of $(2,8]$.
When both agents are given enough expert assistance, they achieve similar
performance. This indicates that HULA efficiently utilizes the expert to
localize itself and improve the task performance.

An important feature of our method is that it does not introduce extra
complexity in training an RL agent. Practically, we stop training when the $Q$
function converges. Therefore, the training efficiency is similar to its
underlying RL algorithm. In contrast, the penalty-based agent requires access
to an expert during training time. Our experiments show that training each
penalty-based agent for the fully-observable environment results in about
$70000$ expert calls. The requirement of an expert during train time limits
its application on learning policies that interact with humans.

Figure 8: Example rollouts of the HULA agent and the penalty-based agent in
the shortcut world.  Figure 9: # of expert calls vs Average returns for the
partially-observable experiments. Curves are generated with the rolling mean
method (window size of 4) for visualization purposes.

##  VI Conclusion

We introduced HULA, a method to learn human-in-the-loop policies by estimating
an RL agent’s uncertainty. We proposed a return variance estimation method
that captures an RL agent’s uncertainty. Our experimental results demonstrate
that HULA can capture an RL agent’s uncertainty and use it to request
assistance from experts to achieve high task performance. An important feature
of our method is that it does not require the presence of an expert during
training. We envision that our approach can be applied to more complex
problems. For example, we plan to extend HULA to continuous RL algorithms
(e.g. DDPG [18]) and learns HitL policies to solve continuous control
problems. We also would like to explore the direction of learning uncertain
agents in other domains, such as language-guided navigation.

## References

  * [1] Oliver Kroemer, Scott Niekum and George Dimitri Konidaris  “A Review of Robot Learning for Manipulation: Challenges, Representations, and Algorithms”  In _CoRR_ abs/1907.03146, 2019  arXiv: http://arxiv.org/abs/1907.03146
  * [2] Christian Daniel, Gerhard Neumann, Oliver Kroemer and Jan Peters  “Hierarchical Relative Entropy Policy Search”  In _Journal of Machine Learning Research_ 17.93, 2016, pp. 1–50  URL: http://jmlr.org/papers/v17/15-188.html
  * [3] Jan Peters and Stefan Schaal  “Reinforcement learning of motor skills with policy gradients”  In _Neural Networks_ 21, 2008, pp. 682–  DOI: 10.1016/j.neunet.2008.02.003
  * [4] E. Theodorou, J. Buchli and S. Schaal  “Reinforcement learning of motor skills in high dimensions: A path integral approach”  In _2010 IEEE International Conference on Robotics and Automation_ , 2010, pp. 2397–2403
  * [5] Riku Arakawa, Sosuke Kobayashi, Yuya Unno, Yuta Tsuboi and Shin-ichi Maeda  “Dqn-tamer: Human-in-the-loop reinforcement learning with intractable feedback”  In _arXiv preprint arXiv:1810.11748_ , 2018
  * [6] Travis Mandel, Yun-En Liu, Emma Brunskill and Zoran Popović  “Where to add actions in human-in-the-loop reinforcement learning”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 31.1, 2017
  * [7] Annie Xie, Fahim Tajwar, Archit Sharma and Chelsea Finn  “When to Ask for Help: Proactive Interventions in Autonomous Reinforcement Learning”  In _arXiv preprint arXiv:2210.10765_ , 2022
  * [8] Jingda Wu, Zhiyu Huang, Chao Huang, Zhongxu Hu, Peng Hang, Yang Xing and Chen Lv  “Human-in-the-loop deep reinforcement learning with application to autonomous driving”  In _arXiv preprint arXiv:2104.07246_ , 2021
  * [9] Max Jaderberg, Volodymyr Mnih, Wojciech M. Czarnecki, Tom Schaul, Joel Z. Leibo, David Silver and Koray Kavukcuoglu  “Reinforcement Learning with Unsupervised Auxiliary Tasks”  In _ArXiv_ abs/1611.05397, 2016
  * [10] Shakir Mohamed and Danilo Jimenez Rezende  “Variational Information Maximisation for Intrinsically Motivated Reinforcement Learning”  In _NIPS_ , 2015
  * [11] Deepak Pathak, Pulkit Agrawal, Alexei A. Efros and Trevor Darrell  “Curiosity-Driven Exploration by Self-Supervised Prediction”  In _2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , 2017, pp. 488–489
  * [12] Gregory Kahn, Adam Villaflor, Vitchyr Pong, Pieter Abbeel and Sergey Levine  “Uncertainty-aware reinforcement learning for collision avoidance”  In _arXiv preprint arXiv:1702.01182_ , 2017
  * [13] Carl-Johan Hoel, Krister Wolff and Leo Laine  “Ensemble quantile networks: Uncertainty-aware reinforcement learning with applications in autonomous driving”  In _arXiv preprint arXiv:2105.10266_ , 2021
  * [14] Will Dabney, Georg Ostrovski, David Silver and Rémi Munos  “Implicit quantile networks for distributional reinforcement learning”  In _International conference on machine learning_ , 2018, pp. 1096–1105  PMLR
  * [15] Ian Osband, John Aslanides and Albin Cassirer  “Randomized prior functions for deep reinforcement learning”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [16] Felipe Leno Da Silva, Pablo Hernandez-Leal, Bilal Kartal and Matthew E Taylor  “Uncertainty-aware action advising for deep reinforcement learning agents”  In _Proceedings of the AAAI conference on artificial intelligence_ 34.04, 2020, pp. 5792–5799
  * [17] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra and Martin A. Riedmiller  “Playing Atari with Deep Reinforcement Learning”  In _ArXiv_ abs/1312.5602, 2013
  * [18] Timothy P. Lillicrap, Jonathan J. Hunt, Alexander Pritzel, Nicolas Manfred Otto Heess, Tom Erez, Yuval Tassa, David Silver and Daan Wierstra  “Continuous control with deep reinforcement learning”  In _CoRR_ abs/1509.02971, 2015


# Estimating expectation values using approximate quantum states

Marco Paini<EMAIL_ADDRESS>Rigetti Computing, 138 Holborn, London, EC1N
2SW, UK.  Amir Kalev<EMAIL_ADDRESS>Information Sciences Institute, University
of Southern California, Arlington, VA 22203, USA.  Dan Padilha  Rigetti
Computing, 138 Holborn, London, EC1N 2SW, UK.  Brendan Ruck  Rigetti
Computing, 2919 Seventh St, Berkeley, CA 94710, USA.

###### Abstract

We introduce an approximate description of an $N$–qubit state, which contains
sufficient information to estimate the expectation value of any observable to
a precision that is upper bounded by the ratio of a suitably–defined seminorm
of the observable to the square root of the number of the system’s identical
preparations $M$, with no explicit dependence on $N$. We describe an
operational procedure for constructing the approximate description of the
state that requires, besides the quantum state preparation, only single–qubit
rotations followed by single–qubit measurements. We show that following this
procedure, the cardinality of the resulting description of the state grows as
$3MN$. We test the proposed method on Rigetti’s quantum processor unit with
12, 16 and 25 qubits for random states and random observables, and find an
excellent agreement with the theory, despite experimental errors.

##  1 Introduction

Quantum density operators represent our knowledge of the state of quantum
systems and give us a way to calculate expectation values and predict
experimental results via the Born rule. Our knowledge of a system’s density
operator is equivalent to our ability to calculate the expectation value of
any observable of the system. However, for an $N$–body system, reconstructing
the system’s density operator requires the knowledge of exponentially many
numbers in $N$. This in turn translates into practical difficulties in
estimating and storing density operators, in the form of density matrices,
even for systems composed of just a few tens of qubits.

Over the last few years, there has been a surge of interest in devising
protocols to reduce the resources required to estimate and classically store
quantum states [1, 2, 3, 4, 5, 6, 7, 8, 9]. There have been numerous
approaches to reduce the cost of performing quantum state tomography by
assuming that the state of the system has a specific structure. These
approaches include, for example: matrix product states tomography [3], where
the state is assumed to have an efficient matrix product state representation
as a matrix; neural network tomography [6, 7], where the state is assumed to
have an efficient deep neural network representation; and compressed sensing
tomography [4, 5], where the state is assumed to have a low-rank density
matrix representation. The theory of Probably Approximately Correct (PAC)
learning was utilised in [2] to show that if, instead of full tomography, we
are interested in approximating with high probability the expectation values
of observables drawn from a distribution, then only polynomially–many samples
in $N$ of the observables’ expectation values are required. It was later also
shown [8] that estimating with high probability the expectation values of a
set of observables can be achieved using only polynomially–many copies in $N$
of the system. As for storage, it is well established [1] that quantum states
that can be represented as low bond–dimension matrix product states can be
efficiently stored on classical computers. Recently it was shown [9] that a
general quantum state of $N$ qubits can be approximately described using its
inner product with an order of $\sqrt{2^{N}}$ stabiliser states. This result
thus provides a polynomial improvement for classically storing quantum states
over a brute–force approach of storing density matrices. In Sec. 4 we discuss
these and other works in connection to the results presented in this paper.

The growth of the number of copies of the system in $N$ for estimation and the
considerable requirements in $N$ for classically storing the information
associated with a generic quantum state lead to the question of whether a
protocol that significantly reduces or eliminates the dependence on $N$ is
admissible. In this work, we give an affirmative answer to the question.
Instead of focusing on the estimation and storage of the density matrix, we
start by considering the equivalent problem of calculating the expectation
value of any observable of the system up to a certain precision. In this
article we show that, for a system of $N$ qubits, the expectation value of any
observable can be calculated with arbitrary precision with a number of copies
of the system independent of $N$, but dependent on a properly–defined seminorm
of the observable. The information obtained through the measurements on the
copies of the system will be used to form an approximate description of the
quantum state – a description with a cardinality that grows linearly in $N$
and thus can be constructed and stored for large systems. The method for
estimating expectation values, and the exact definition of an observable
seminorm and of the storage approximate quantum state will be given in the
next section.

##  2 Expectation value estimation and the approximate quantum state

We will consider specifically a system of $N$ qubits. However, the
considerations that will be made are applicable in general to other systems,
with finite and infinite–dimensional Hilbert spaces, since they are based on
the general theory of irreducible representations of groups which provide
explicit forms for a basis in the space of linear bounded operators defined on
a Hilbert space. We refer the reader to Refs. [10, 11] where these
considerations have been developed.

We start by presenting the method to estimate the expectation value of any
observable and notion of the approximate state for a single–qubit system. We
then generalise these ideas for the more general $N$–qubit case.

###  2.1 One qubit

A single–qubit density operator can be written as [10, 11]

$\displaystyle\rho=\int_{\Sigma}\frac{d\vec{n}}{4\pi}\sum_{m=\pm
1}p(m,\vec{n})K_{1}(m,\vec{n}),$ (1)

with the unitary spherical surface ${\Sigma}$ as integration domain,
$\vec{n}=(\cos\varphi\sin\vartheta,\sin\varphi\sin\vartheta,\cos\vartheta)$,
with $\vartheta\in[0,\pi],\varphi\in[0,2\pi)$, is a unit vector in
$\mathds{R}^{3}$, $p(m,\vec{n})$ is the probability of obtaining the result
$m\in\\{-1,+1\\}$ when measuring the qubit in the eigenbasis of
$\vec{\sigma}\cdot\vec{n}$, where
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of the Pauli
operators, and $K_{1}(m,\vec{n})$ is the 1–qubit kernel operator given by

$\displaystyle
K_{1}(m,\vec{n})=\frac{1}{2}(\mathbb{1}+3\,m\,\vec{\sigma}\cdot\vec{n}),$ (2)

see App. A for more details.

Let us define a single–qubit (unbiased) estimator function

$\displaystyle{R}_{1}[\;\cdot\;](m,\vec{n})=\text{Tr}\Bigl{[}\;\cdot\;K_{1}(m,\vec{n})\Bigr{]},$
(3)

where the $\cdot$ represents the argument of the function ${R}_{1}$, i.e. a
single–qubit observable (for a lighter notation we sometime omit the explicit
dependence of the estimator on $(m,\vec{n})$). In particular, the expectation
value of the Pauli operator $\sigma_{i}$, with $i=0,1,2,3$ and
$\sigma_{0}:=\mathbb{1},$ $\sigma_{1}:=\sigma_{x},\,\sigma_{2}:=\sigma_{y},$
and $\sigma_{3}:=\sigma_{z}$, can be written in terms of ${R}_{1}[\sigma_{i}]$
as

$\displaystyle\langle\sigma_{i}\rangle=\int_{\Sigma}\frac{d\vec{n}}{4\pi}\sum_{m=\pm
1}p(m,\vec{n}){R}_{1}[\sigma_{i}](m,\vec{n}).$ (4)

Calculating the estimator of the Pauli operators we arrive at

$\displaystyle{R}_{1}[\mathbb{1}](m,\vec{n})$ $\displaystyle=1,$ (5)
$\displaystyle{R}_{1}[\sigma_{\alpha}](m,\vec{n})$
$\displaystyle=3\,m\,n_{\alpha},\,\,\text{for}\,\,\alpha=x,y,z.$ (6)

The estimator ${R}_{1}[\sigma_{\alpha}](m,\vec{n})$ is a random variable due
to its dependence on the stochastic measurement outcome $m$ and the random
unit vector $\vec{n}$. Hereafter, unless noted otherwise, we use the subscript
$\alpha$ to index $x,y,z$ and subscript $i$ to index $0,1,2,3$. In App. A we
show that the variance of the estimator ${R}_{1}[\sigma_{\alpha}]$ is upper
bounded by $3$,

$\displaystyle\text{Var}({R}_{1}[\sigma_{\alpha}](m,\vec{n}))$
$\displaystyle\leqslant 3.$ (7)

Since ${R}_{1}[\mathbb{1}](m,\vec{n})$ is independent of $m$ and $\vec{n}$, it
is deterministic and thus has a null variance.

Equation (4) implies that we can construct an estimate of
$\langle\sigma_{i}\rangle$ using the following Monte Carlo sampling procedure:

  1. 1.

Sample a point $\vec{n}$ uniformly from the surface of a unit sphere,

  2. 2.

prepare the qubit in the state $\rho$ and measure it in the eigenbasis of
$\vec{\sigma}\cdot\vec{n}$,

  3. 3.

record the outcome $m\in\\{-1,1\\}$,

  4. 4.

repeat steps 1–3 $M$ times.

The estimate of $\langle\sigma_{i}\rangle$, which we denote by
$\langle\sigma_{i}\rangle_{M}$, is simply given by the average of
corresponding estimator ${R}_{1}[\sigma_{i}](m,\vec{n})$ over the data set:

$\displaystyle\langle\sigma_{i}\rangle_{M}=\frac{1}{M}\sum_{j=1}^{M}{R}_{1}[\sigma_{i}](m_{j},\vec{n}_{j}),$
(8)

where
$\vec{n}_{j}=(\cos\varphi_{j}\sin\vartheta_{j},\sin\varphi_{j}\sin\vartheta_{j},\cos\vartheta_{j})$
is the unit vector sampled at the $j$th iteration and $m_{j}$ is the
corresponding measurement outcome. Note that the estimate
$\langle{\mathbb{1}}\rangle_{M}=1$ for any value of $M$ and it has zero
variance. Moreover, by construction, in the limit $M{\to}\infty$,
$\langle\sigma_{\alpha}\rangle_{M}{\to}\langle\sigma_{\alpha}\rangle$, for
$\alpha=x,y,z$. The value of $M$ determines the precision of the estimation.
For $M$ sufficiently large, $\langle\sigma_{\alpha}\rangle_{M}$ has a normal
statistical error, which due to Eq. (7) is smaller than or equal to
$\sqrt{3/M}$.

The measurement in the eigenbasis of $\vec{\sigma}\cdot\vec{n}$, at step 2
above, can be implemented by first applying a unitary operator $U$, that
rotates this basis to the eigenbasis of $\sigma_{z}$, and then measuring in
the eigenbasis $\sigma_{z}$ (the computational basis). This unitary operator
is given by $U=e^{i\frac{\vartheta}{2}\vec{\sigma}\cdot\vec{n}_{\perp}}$,
where $\vec{n}_{\perp}=(-\sin\varphi,\cos\varphi,0)$. We call the vector
$s_{j}=(m_{j},\vartheta_{j},\varphi_{j})$ the $j$-th snapshot of the qubit
system. We refer to the collection of the $M$ snapshots, ${\cal
S}=\\{s_{j}\\}_{j=1}^{M}$, that we obtain following the sampling procedure
above, as the approximate quantum state, as it enables us to estimate the
expectation value of any single–qubit observable to a certain precision as we
now describe.

Let $O$ be a single–qubit observable, expanded in the Pauli basis
$\\{\sigma_{i}\\}_{i=0}^{3}$,

$\displaystyle O=\sum_{i=0}^{3}a_{i}\,\sigma_{i}.$ (9)

Given the estimated expectation values of $\\{\sigma_{i}\\}_{i=0}^{3}$, by
linearity of the estimator function (3), we can construct an estimate of
$\langle O\rangle$ as

$\displaystyle\langle{O}\rangle_{M}=\sum_{i=0}^{3}a_{i}\,\langle\sigma_{i}\rangle_{M}.$
(10)

We prove in App. A that for $M$ sufficiently large, $\langle{O}\rangle_{M}$
has a normal statistical error smaller than or equal to $\|O\|/{\sqrt{M}}$,
that is

$\displaystyle\text{Std}(\langle{O}\rangle_{M})\leqslant\frac{\|O\|}{\sqrt{M}}\,,$
(11)

where Std stands for ‘standard deviation’ and

$\displaystyle\|O\|=\sqrt{3\sum_{i=1}^{3}a_{i}^{2}}$ (12)

is a seminorm of $O$ that depends on its expansion in the Pauli operator
basis.

###  2.2 $N$ qubits

The method developed for one qubit in the previous subsection is easily
generalised to $N$ qubits, in the assumptions that the qubits are
distinguishable and that a qubit–specific measurement capability is available.
The space of the $N$–qubit system is the tensor product of the single–qubit
spaces and thus we can write a general $N$–qubit density operator as [10, 11]

$\displaystyle\rho=$
$\displaystyle\bigg{(}{\frac{1}{(4\pi)^{N}}\prod_{k=1}^{N}\int_{{\Sigma}_{k}}{d{\vec{n}}_{k}}\sum_{m_{k}}}\bigg{)}p(\\{m,\vec{n}\\})K(\\{m,\vec{n}\\}),$
(13)

where $\\{m,\vec{n}\\}$ is a shorthand for
$(m_{1},{\vec{n}}_{1},\ldots,m_{N},{\vec{n}}_{N})$, the $N$–qubit kernel
operator $K(\\{m,\vec{n}\\})$ equals to the tensor product of single–qubit
kernels

$\displaystyle
K(\\{m,\vec{n}\\})=\bigotimes_{k=1}^{N}K_{1}(m_{k},{\vec{n}}_{k}),$ (14)

and $p(\\{m,\vec{n}\\})$ represents the joint probability of finding the
outcome $(m_{1},\ldots,m_{N})\in\\{-1,1\\}^{N}$ when measuring the qubits in
the product eigenbasis of
${\vec{\sigma}}_{1}\cdot{\vec{n}}_{1}\otimes\cdots\otimes{\vec{\sigma}}_{N}\cdot{\vec{n}}_{N}$.
Similarly to the single qubit case, this measurement can be implemented with
single–qubit rotations followed by a measurement in the computational basis on
each qubit.

We define the $N$–qubit (unbiased) estimator function

$\displaystyle{R}[\;\cdot\;](\\{m,\vec{n}\\})$
$\displaystyle=\text{Tr}\Bigl{[}\;\cdot\;K(\\{m,\vec{n}\\})\Bigr{]}.$ (15)

As a result of the tensor product structure of the kernel operator,
${R}[P_{i}]$, where $P_{i}=\sigma_{i_{1}}\otimes\cdots\otimes\sigma_{i_{N}}$,
$i=(i_{1},\ldots\,i_{N})$, is a monomial of Pauli operators (a ‘Pauli
monomial’ or a ‘Pauli string’ for short), and hereafter $i_{k}=0,\ldots,3$, is
simply the product of the corresponding single–qubit random variables,

$\displaystyle{R}[P_{i}](\\{m,\vec{n}\\})$
$\displaystyle=\prod_{k=1}^{N}{R}_{1}[\sigma_{i_{k}}](m_{k},\vec{n}_{k}),$
(16)

where $R_{1}[\sigma_{i_{k}}](m_{k},\vec{n}_{k})$ is given in Eqs. (5) and (6).

Similarly to the single qubit case of the previous section, we can write the
expectation value of $P_{i}$ with respect to the state of the system in terms
of ${R}[P_{i}]$ as

$\displaystyle\langle P_{i}\rangle=$
$\displaystyle\bigg{(}{\frac{1}{(4\pi)^{N}}\prod_{k=1}^{N}\int_{{\Sigma}_{k}}{d{\vec{n}}_{k}}\sum_{m_{k}}}\bigg{)}$
(17) $\displaystyle\times p(\\{m,\vec{n}\\}){R}[P_{i}](\\{m,\vec{n}\\}).$

Eq. (17) implies that an estimate $\langle P_{i}\rangle_{M}$ of $\langle
P_{i}\rangle$ can be calculated by performing the single–qubit tomographic
procedure introduced in the above section, in parallel on all of the $N$
qubits, as in Fig. 1. We call the vector
$s_{j}=({m_{1}}_{j},{{}{\vec{n}}_{1}}_{j},\ldots,{m_{N}}_{j},{{}{\vec{n}}_{N}}_{j})$
the $j$th snapshot of a system of $N$ qubits, and we refer to a set $\cal S$
of such $M$ independent snapshots as the approximate state of the $N$–qubit
system.

The use of single–qubit Haar–random basis measurements was recently proposed,
and experimentally implemented, to probe certain properties of an $N$–qubit
system, such as purity [12, 13] and the Frobenius inner product of two states
[14]. In light of these results, as we now show, the measurement protocol we
introduced above is a universal protocol for estimating the expectation value
of any $N$–qubit observable. The protocol is efficient in terms of time and
space resources, as long as the observable admits an efficient representation
in the Pauli basis and has a small seminorm (e.g. not growing exponentially
with $N$).

To start with, we can use the approximate state $\cal S$ to calculate an
estimate $\langle P_{i}\rangle_{M}$ of $\langle{P_{i}}\rangle$ as

$\displaystyle\langle P_{i}\rangle_{M}$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}{R}[P_{i}](\\{m,\vec{n}\\})$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}\prod_{k=1}^{N}R_{1}[\sigma_{i_{k}}]({m_{k}}_{j},{{}{\vec{n}}_{k}}_{j}),$
(18)

where the first equality follows from Eq. (16). We show in App. B that for $M$
sufficiently large, $\langle P_{i}\rangle_{M}$ ($i\neq
0_{v}\doteq(0,\ldots,0)$) has a normal statistical error that is smaller than
or equal to $\sqrt{3^{r_{i}}/M}$, where $r_{i}$ represents the weight of
$P_{i}$, i.e., the number of qubits on which it acts non–trivially. Note that
the estimated expectation value for the $N$–qubit identity operator
$\langle\mathbb{1}^{\otimes N}\rangle_{M}=1$ and it has zero variance.

With the ability to efficiently calculate the estimated expectation value of
any given Pauli monomial operator, we can estimate the expectation value of
$N$–qubit observables. Let $O$ be an $N$–qubit observable such that

$\displaystyle O=\sum_{i}\,a_{i}\,P_{i}\,,$ (19)

From linearity of the estimator function (15) and from Eq. (18), the estimated
expectation value of $O$ is given by

$\displaystyle\langle{O}\rangle_{M}$
$\displaystyle=\sum_{i}a_{i}\langle{P_{i}}\rangle_{M}$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}\sum_{i}a_{i}\prod_{k=1}^{N}R_{1}[\sigma_{i_{k}}]({m_{k}}_{j},{{}{\vec{n}}_{k}}_{j}).$
(20)

Similar to $\langle{P_{i}}\rangle_{M}$, $\langle{O}\rangle_{M}$ is an unbiased
estimate, i.e., $\langle{O}\rangle_{M}{\to}\langle{O}\rangle$ in the limit
$M{\to}\infty$. Practically, for the efficient calculation of (20), the sum on
$i$ should have at most a polynomial number of terms in $N$.

While we assumed in Eq. (19) that the observable admits an efficient
representation in terms of Pauli monomials, as a result of the tensor product
structure of the kernel operator, Eq. (14), our procedure is efficient to
calculate the expectation value of any operator of the form

$\displaystyle O=\sum_{i}\,a_{i}\,O_{i_{1}}\otimes O_{i_{2}}\cdots\otimes
O_{i_{N}}\,,$ (21)

where $O_{i_{k}}$ is a single-qubit operator (on the $k$th qubit space). If
this is the case then,

$\displaystyle\langle
O\rangle_{M}=\frac{1}{M}\sum_{i}a_{i}\sum_{j=1}^{M}\prod_{k=1}^{N}R_{1}[O_{i_{k}}].$
(22)

Note that any Pauli monomial $P_{i}$ is a special case of $O_{i_{1}}\otimes
O_{i_{2}}\cdots\otimes O_{i_{k}}$.

We show in App. B that the standard deviation of $\langle{O}\rangle_{M}$ has
the upper bound

$\displaystyle{\rm
Std}(\langle{O}\rangle_{M})\leqslant\frac{{\|O\|}}{\sqrt{M}},$ (23)

with the seminorm of $O$ defined as

$\displaystyle\|O\|=\sqrt{\sum_{i,j(\neq
0_{v})}3^{r_{ij}}\Delta_{ij}|a_{i}||a_{j}|},$ (24)

where the sum is extended to all values of $i$ and $j$ except
$0_{v}\doteq(0,0,\ldots,0)$, $r_{ij}$ is the number of indices $k$ for which
$i_{k}\neq 0\wedge j_{k}\neq 0$, $\Delta_{ij}=0$ if there exists a $k$ such
that $i_{k}\neq 0\wedge j_{k}\neq 0\wedge i_{k}\neq j_{k}$ and $\Delta_{ij}=1$
otherwise. Eq. (24) reduces to Eq. (12) when $O$ is an observable on a
single–qubit space.

Importantly, in this procedure, the expectation values estimation error has no
global dependence on the system size $N$. The error only depends on the
observable seminorm and the number of snapshots. One could argue that the
presence in the seminorm of $3^{r_{ij}}$ reflects an exponential growth of the
estimates’ variance in $N$, as, for example, for an observable given by the
tensor product of $N$ Pauli operators, which has seminorm ${3}^{N/2}$.
However, the seminorm is a property of the observable and not of the dimension
of the space: an $N$–qubit observable $O$ and its extension to a larger
$(N{+}L)$–qubit space, $O\otimes\mathbb{1}^{\otimes L}$, have the same
seminorm (24), independent of $L$. Irrespective of $N$, because of Eq. (23),
the expectation values of all observables with unit seminorm can be obtained
with statistical error bounded by $1/\sqrt{M}$. Any observable with seminorm
different from zero can be obtained from an observable with unit seminorm
through multiplication by a real number, representing in fact the observable’s
seminorm. For a given number of identical preparations and measurements, the
statistical error in the expectation value of the generic observable has a
bound only dependent on how “large” the observable is. In App. B, we also
explain how the seminorm ${\|O\|}_{2}$, defined by

$\displaystyle{\|O\|}_{2}=\sqrt{\sum_{i\neq 0_{v}}3^{r_{i}}\,a_{i}^{2}},$ (25)

divided by $\sqrt{M}$, with $r_{i}$ representing the weight of $P_{i}$, is
generally a good approximation for the standard deviation of
$\langle{O}\rangle_{M}$. This result may be important for applications of the
method since ${\|O\|}_{2}$ can be significantly smaller than $\|O\|$ and
since, for a chosen statistical error, the number of identical preparations
required is inversely proportional to the variance of $\langle{O}\rangle_{M}$.
See App. C for further details.

Figure 1: Obtaining the approximate quantum state. (a) To construct the
approximate description of an $N$–qubit state we repeat the following steps
$M$ times: (1) prepare a copy of the $N$–qubit state, (2) measure each qubit
in a random direction
${\vec{n}}_{j}=(\cos\varphi_{j}\sin\vartheta_{j},\sin\varphi_{j}\sin\vartheta_{j},\cos\vartheta_{j})$,
(3) record the measurement outcome of each qubit. (b) A circuit realization of
the protocol. The single–qubit rotations are defined by $\vartheta_{j}$ and
$\vec{n}_{j_{\perp}}=(-\sin\varphi_{j},\cos\varphi_{j},0)$, $j=1,\ldots,N$.

We close this part with a remark on the notion of the approximate quantum
state, $\cal S$. We have shown that the approximate description of the state,
as we have defined it, comes with an explicit operational procedure to
estimate the expectation value of any generic observable with a rigorous bound
on the estimation statistical error. Therefore, the approximate state together
with the rule (20) provides a description of a physical system in terms of
statistical inferences based on observed experimental data and, in principle,
could be assumed as primitive concepts instead, or, at least, as not requiring
the existence of density operators and the Born rule.

##  3 Experimental results

As a proof–of–concept, we have tested the theory on Rigetti Computing quantum
processor unit (QPU), Aspen7-28Q-A, with three systems, involving $N=12,16$
and $25$ qubits, respectively. In each experiment the system is prepared in a
random state. The state preparation circuit is composed of two layers of gates
applied to the computational state $|0\rangle^{\otimes N}$. In the first
layer, the gates are chosen at random from a flat distribution over the
single–qubit gate set $\\{X,Y,Z,H,S,T\\}$ and applied to $N/2$ qubits, chosen
at random. While in the second layer the two–qubit gate
$XY(\alpha)=e^{-i\alpha(X\otimes X+Y\otimes Y)}$ is applied to $50\%$ of pairs
of qubits, chosen at random, with different values of $\alpha$ sampled
uniformly from $[0,2\pi)$. This gate composition results in shallow quantum
circuits with $N/2$ random single–qubit gates and $N/4$ random two–qubit
gates, to ensure that the experimentally prepared state has a high fidelity
with the target state.

For each system, an approximate description of the prepared state was obtained
using $M=10^{4}$ snapshots (to check the results to $0.01$ accuracy). The
approximate state was then used to estimate the expectation value of two sets
of observables: (1) 20 random observables, (2) 20 random projectors onto
computational basis states (we have selected 20 for the number of states so
that it would be visually ‘easy’ to check the results in Fig. 2 against the
second standard deviation given it is at $\sim$95%). Each random observable in
the first set is a linear combination of 20 Pauli monomials. Each Pauli
monomial is constructed by randomly sampling from the set
$\\{\sigma_{i}\\}_{i=0}^{3}$ for each qubit, whereas the linear combination
coefficients (the $a_{i}$’s) are sampled uniformly from the $(0,1]$. The
projectors onto computational basis states were sampled uniformly from the set
$\\{0,1\\}^{N}$. For simplicity in the analysis, the 20 random observables are
normalised to have a unit seminorm (the projectors already have unit seminorm
with probability exponentially close to 1 as a function of $N$, see App. C).

Fig. 2 shows the results for the 25–qubit experiment. In this figure, we
compare the estimated expectation values of random observables and of random
computational basis projectors to their expectation values with respect to the
25–qubit target state. The results are plotted as a function of the number of
snapshots. Table 1 summarises the experimental results of the $12$–, $16$– and
$25$–qubit systems for $10^{4}$ snapshots. The figure and the table show that
the experimental results generally mirror the theory, with the fraction of
estimates within one and two standard deviations from their true values not
far from $68\%$ and $95\%$ respectively. Other than statistical
considerations, the cases with a fraction of the estimates greater than the
expected percentages from a normal distribution can be explained in terms of
how tightly the observables’ seminorms bound the standard deviation of the
estimators, in particular for the case of the projectors.

Figure 2: The error in the estimated expectation value of (a) 20 random
observable and (b) 20 random projectors onto computational basis states in a
25–qubit system as a function of the number of snapshots. Table 1:  The
fraction of the estimated expectation values, of (a) random observables and
(b) random projectors onto computational basis states, that are within one and
two standard deviations from their true expectation values for $10^{4}$
snapshots.

To complete the picture, including an explanation for the cases of a lower
fraction of the estimates than predicted, we need to start examining the
effects of hardware errors. It is worth noting that we obtain agreement with
the theory despite an order of $p_{\rm err}=5\%$ (averaged) readout error per
qubit, which amounts to flipping the measurement outcome from $+1$ to $-1$ and
vice versa. The resilience of our protocol to single–qubit readout errors can
be attributed to the randomness of the measurement bases, as was discussed in
[15]. Moreover, the following considerations provide additional insights in
the robustness of the estimation protocol to readout errors.

Firstly, let us consider the estimation of the expectation value of
single–qubit observable. Let $\cal S$ be the approximate state description
composed of $M$ snapshots. Suppose that in $M_{\rm err}$ of the snapshots we
obtained the erroneous measurement outcome, $M_{\rm err}/M\approx p_{\rm
err}$. The estimation of the expectation value of any single–qubit Pauli
operator $\sigma_{\alpha}$ ($\alpha=x,y,z$) is

$\displaystyle\langle{\sigma_{\alpha}}\rangle_{M}$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}3m_{j}{n_{\alpha}}_{j}$
$\displaystyle=\frac{M-M_{\rm err}}{M}\Big{(}\frac{1}{M-M_{\rm
err}}\sum_{j=1}^{M-M_{\rm err}}3m_{j}{n_{\alpha}}_{j}\Big{)}$
$\displaystyle+\frac{M_{\rm err}}{M}\Big{(}\frac{1}{M_{\rm
err}}\sum_{j=1}^{M_{\rm err}}3m_{j}^{\rm err}{n_{\alpha}}_{j}\Big{)},$ (26)

where $m_{j}^{\rm err}$ is the erroneous outcome. Since $m_{j}^{\rm
err}=-m_{j}$ we arrive at

$\displaystyle\langle{\sigma_{\alpha}}\rangle_{M}$
$\displaystyle=\frac{M-M_{\rm err}}{M}\Big{(}\frac{1}{M-M_{\rm
err}}\sum_{j=1}^{M-M_{\rm err}}3m_{j}{n_{\alpha}}_{j}\Big{)}$
$\displaystyle-\frac{M_{\rm err}}{M}\Big{(}\frac{1}{M_{\rm
err}}\sum_{j=1}^{M_{\rm err}}3m_{j}{n_{\alpha}}_{j}\Big{)},$ (27)

which, for sufficiently large $M$ and $M_{\rm err}$, can be approximated by

$\displaystyle\langle{\sigma_{\alpha}}\rangle_{M}$
$\displaystyle\approx(1-p_{\rm err})\langle{\sigma_{\alpha}}\rangle-p_{\rm
err}\langle{\sigma_{\alpha}}\rangle$ $\displaystyle=(1-2p_{\rm
err})\langle{\sigma_{\alpha}}\rangle.$ (28)

Therefore, for a single–qubit observable $O=\sum_{i}a_{i}\sigma_{i}$ we arrive
at

$\displaystyle\langle{O}\rangle_{M}$ $\displaystyle\approx
a_{0}\mathbb{1}+\sum_{i=1}^{3}a_{i}(1-2p_{\rm
err})\langle{\sigma_{i}}\rangle.$ (29)

Thus, the readout error introduces an attenuator $(1-2p_{\rm err})$ to each
coefficient different than $a_{0}$ in the expectation value.

Moving to the $N$–qubit case, we first note that only odd numbers of readout
errors per snapshot contribute to the estimation error of an $N$–qubit
expectation value. This is because the estimated expectation value is
calculated via a product of the measurement outcomes which take values $\pm
1$, c.f. Eq. (18). Therefore, if a given snapshot includes readout errors on
an even number of qubits, there are even number of erroneous outcomes
$m_{j}^{\rm err}$ such that $m_{j}^{\rm err}=-m_{j}$, and the estimated
expectation value remains unaffected by these errors. Let $p_{\rm odd}(r)$ be
the probability that in a given snapshot there is an odd number of readout
errors on $r$ qubits. Given a readout error $p_{\rm err}$ per qubit,

$\displaystyle p_{\rm odd}(r)$ $\displaystyle=\sum_{k\leavevmode\nobreak\
\rm{odd}}{r\choose{k}}p_{\rm err}^{k}(1-p_{\rm err})^{r-k}$
$\displaystyle=\frac{1}{2}\Big{(}1-(1-2p_{\rm err})^{r}\Big{)}.$ (30)

Similarly to the single qubit analysis, one can then show that the expectation
value of a Pauli monomial $P_{i}$ with weight $r_{i}$ is

$\displaystyle\langle{P_{i}}\rangle_{M}$ $\displaystyle\approx(1-2p_{\rm
odd}(r_{i}))\langle{P_{i}}\rangle=(1-2p_{\rm
err})^{r_{i}}\langle{P_{i}}\rangle.$ (31)

We therefore obtain that the estimated expectation value of an $N$–qubit
observable $O=\sum_{i}\,a_{i}\,P_{i}$ is given by

$\displaystyle\langle{O}\rangle_{M}$
$\displaystyle\approx\sum_{i}a_{i}(1-2p_{\rm
err})^{r_{i}}\langle{P_{i}}\rangle.$ (32)

This equation indicates that the readout error does not depend on $N$. It
introduces an attenuator $(1-2p_{\rm err})^{r_{i}}$ to each term in the
expectation value. Since the attenuator is an exponentially decreasing
function in the weight $r_{i}$, only low–weight Pauli monomials contribute to
the expectation value $\langle{O}\rangle_{M}$. Thus, the proposed method is
robust against readout error when estimating expectation values that have
contributions dominantly from low–weight Pauli monomials (this does not
necessarily exclude observables, such as projectors, that have contributions
from all Pauli monomials: the contribution of higher weight Pauli monomials
could be small on a class of states of interest, as in the case of the low
entanglement states of this section’s random projectors experiment). Finally,
we remark that other errors, e.g. gate operation errors, could be included in
this framework in a similar fashion and that, clearly, (32) is easily
generalisable to the case of different readout errors per qubit.

##  4 Related work

Quantum state tomography. In general, estimating an $N$–qubit state to error
$\epsilon$ in trace distance requires ${\cal O}(2^{N}/\epsilon^{2})$ copies of
the state [16, 17]. The estimation of the density matrix with our tomographic
procedure does not break this lower bound. If we choose the basis of $N$–qubit
Pauli operators, Eq. (23) tells us that a few matrix elements will require $M$
to be exponentially large in $N$. A different basis, such as the standard
basis $|i\rangle\langle j|$, has exponentially many matrix elements bounded by
1. However, calculating the trace distance requires adding exponentially many
of them, creating again the need for an exponentially large $M$. The key point
of this paper, however, is to be able to estimate approximately and
efficiently expectation values without resorting to the notion of density
matrix.

Shadow tomography. In shadow tomography the task is to estimate with high
probability the expectation values of a set of $K$ observables, to a certain
precision $\epsilon$. It was shown that [8, 18] this can be done using
$M={\cal O}(\log(K)^{2}N^{2}\epsilon^{-8})$ copies of the state, by
implementing a non–local measurement on all copies of the state. Recently,
these results were improved [19]. The authors of [19] proposed an algorithm
where $M={\cal O}(B\log(K)/\epsilon^{2})$ copies of the state, with
$B\geqslant\max_{i}\|{O_{i}}\|^{2}_{\rm shadow}$, allow one to predict the
expectation values of $K$ observables, $\\{O_{i}\\}_{i=1}^{K}$, to a precision
$\epsilon$. The norm $\|{\cdot}\|_{\rm shadow}$ depends on the particular
measurement procedure used in the algorithm. For example, for random Pauli
measurements the authors show (Proposition 3 there) that for a $k$-local
observable $\|{O}\|_{\rm shadow}\leqslant 2^{k}\|{O}\|_{\infty}$, where
$\|\cdot\|_{\infty}$ denotes the spectral norm.

The algorithm of [19] has a similar structure to the algorithm proposed in
this work. Similarly to [19], the choice of randomised single–qubit
measurement directions of Eq. (13) is not the only possible one. The framework
of [10, 11] generates valid formulas for any irreducible projective
representation in the space of the system, including for example the Clifford
group and the representation of $SU(2)$ in the space of a spin
$s=(2^{N}-1)/2$. Furthermore, one could also easily derive a general
expression for the norm, bound of the standard deviation of the estimator, as
$\sqrt{\text{Tr}[O^{2}]}$ multiplied by a factor only dependent on the group
representation chosen. It is important to note, however, that the use of the
Cauchy-Schwarz inequality in the derivation would generate suboptimal norms,
often exponentially larger than the true standard deviation of the estimator.
A tight norm (or seminorm) is key, as it gives a guarantee of performance with
a lower number of preparations, especially if the separation from a suboptimal
one is exponential. The closeness to the standard deviation of the estimator
is what justifies the importance of $\|{O}\|_{2}$ of (25): $\|{O}\|_{2}$ can
be exponentially smaller than norms depending exponentially on the locality or
on the dimension of the system (App. C for the case of projectors onto the
computational basis).

Finally, while the task in shadow tomography is to estimate with high
probability the expectation values of a set of $K$ observables, to a certain
precision $\epsilon$, we focus here on the estimation with high probability
the expectation values of a given, yet arbitrary, observable to a certain
precision $\epsilon$. We show that accomplishing this task requires $M={\cal
O}(\|{O}\|^{2}/\epsilon^{2})$ copies of the state. If a certain precision is
required for the combined evaluation of expectation values of several
observables, one would simply calculate the joint probability of several
normal variables of being smaller than a certain value which would lead to a
log contribution to the error in the number of observables.

##  5 Conclusions and outlook

We have introduced a concept of approximate description of a quantum state for
a system of $N$ qubits. We have shown that the approximate state comes with an
explicit operational procedure to estimate the expectation value of any
observable that admits an efficient representation in the Pauli basis with
statistical error only dependent on the cardinality of the approximate quantum
state and on the observable’s seminorm, independent of $N$ and growing at most
exponentially with the locality of the observable. The operational procedure
to obtain the approximate state, other than the initial state preparations,
only requires single–qubit operations and measurement. The natural robustness
of random single–qubit operations to experimental noise, as we saw in the
experimental results and also discussed in [15], suggests that the proposed
protocol could be appropriate for NISQ computers.

We expect the notion of the approximate quantum state to enable improvement
over existing variational optimisation methods, as it would not require
separate iterations for each set of values of the parameters. Machine learning
optimisations, in particular for generative models known as Born Machines
[20], may be an interesting case, since the explicit dependence of the
expectation values of the projectors onto computational basis states on the
variational parameters [21] may be obtained even for a large number of qubits,
since, as shown in App. B, the estimators for the projectors onto
computational basis state have, with high probability, variance bounded by 1
for all values of $N$.

The definition of the snapshots, constituting an approximate quantum state,
can be generalised, too. The snapshots can in fact be derived from a group
different from $SU(2)$. This generalisation has interest that goes beyond the
purely mathematical exercise. As we have seen, $SU(2)$ tomography works well
for observables expressed as linear combinations of products of Pauli
operators, especially if the weights of the Pauli monomials are low.
Electronic Hamiltonians can be represented as linear combinations of products
of Pauli operators, but at the cost, in general, of introducing Pauli
monomials with weights increasing with $N$. In second quantisation, electronic
Hamiltonians contain linear combinations of products of only two or four
ladder operators for any $N$. The bosonic tomographic formulas of [10] can be
extended to the fermionic case, choosing the additive group of Grassmann
numbers [22], instead of the additive group of complex numbers, as tomographic
group and expressing the snapshots as vectors of sampled real numbers
descending from Grassmann numbers and of measured presence or absence of
electrons in orbitals. The estimators for ladder operators and for their
linear combinations can then be calculated directly without the need for the
mapping to Pauli operators and without the consequent appearance of weights of
Pauli monomials growing with $N$.

##  6 Acknowledgments

MP would like to thank Chad Rigetti and Michael Brett for the support, Amy
Brown, Colm Ryan, David Garvin, Duncan Fletcher, Eric Peterson, Genya
Crossman, John Lapeyre, John Macaulay, Juan Bello-Rivas, Mark Hodson, Max
Henderson, Riccardo Manenti, Robert Smith and Sohaib Alam for the reviews and
the feedback. AK acknowledges the support from the National Science Foundation
Grant No. 2037301.

## References

  * Vidal [2003] G. Vidal.  Efficient classical simulation of slightly entangled quantum computations.  _Phys. Rev. Lett._ , 91:147902, Oct 2003.  URL https://doi.org/10.1103/PhysRevLett.91.147902.
  * Aaronson [2007] S. Aaronson.  The learnability of quantum states.  _Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences_ , 463(2088):3089–3114, 2007\.  URL https://doi.org/10.1098/rspa.2007.0113.
  * Cramer et al. [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu.  Efficient quantum state tomography.  _Nat. Comm._ , 1:149, 2010.  URL https://doi.org/10.1038/ncomms1147.
  * Gross et al. [2010] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert.  Quantum state tomography via compressed sensing.  _Phys. Rev. Lett._ , 105(15):150401, 2010.  URL https://doi.org/10.1103/PhysRevLett.105.150401.
  * Kalev et al. [2015] A. Kalev, R. Kosut, and I. Deutsch.  Quantum tomography protocols with positivity are compressed sensing protocols.  _NPJ Quantum Information_ , 1:15018, 2015.  URL https://doi.org/10.1038/npjqi.2015.18.
  * Gao and Duan [2017] X. Gao and L.-M. Duan.  Efficient representation of quantum many-body states with deep neural networks.  _Nat. Comm._ , 8:662, 2017.  URL https://doi.org/10.1038/s41467-017-00705-2.
  * Torlai et al. [2017] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo.  Neural-network quantum state tomography.  _Nat. Phys._ , 18:447, 2017.  URL https://doi.org/10.1038/s41567-018-0048-5.
  * Aaronson [2018] S. Aaronson.  Shadow tomography of quantum states.  In _Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing_ , STOC 2018, pages 325–338, New York, NY, USA, 2018. ACM.  ISBN 978-1-4503-5559-9.  URL https://doi.org/10.1145/3188745.3188802.
  * Gosset and Smolin [2019] D. Gosset and J. Smolin.  A compressed classical description of quantum states.  In Wim van Dam and Laura Mancinska, editors, _14th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2019)_ , volume 135 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 8:1–8:9, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-112-2.  URL https://doi.org/10.4230/LIPIcs.TQC.2019.8.
  * Paini [2000] M. Paini.  Quantum tomography via group theory.  arXiv:quant-ph/0002078, 2000.  URL https://arxiv.org/pdf/quant-ph/0002078.pdf.
  * D’Ariano et al. [2003] G. M. D’Ariano, L. Maccone, and M. Paini.  Spin tomography.  _J. Opt. B: Quantum and Semiclassical Optics_ , 5(1):77–84, jan 2003.  URL https://doi.org/10.1088/1464-4266/5/1/311.
  * Brydges et al. [2019] T. Brydges, A. Elben, P. Jurcevic, C. Vermersch, B.and Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos.  Probing Rényi entanglement entropy via randomized measurements.  _Science_ , 364(6437):260–263, 2019.  URL https://doi.org/10.1126/science.aau4963.
  * Elben et al. [2020a] A. Elben, J. Yu, G. Zhu, M. Hafezi, F. Pollmann, P. Zoller, and B. Vermersch.  Many-body topological invariants from randomized measurements in synthetic quantum matter.  _Science Advances_ , 6(15), 2020a.  URL https://doi.org/10.1126/sciadv.aaz3666.
  * Elben et al. [2020b] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller.  Cross-platform verification of intermediate scale quantum devices.  _Phys. Rev. Lett._ , 124:010504, Jan 2020b.  URL https://doi.org/10.1103/PhysRevLett.124.010504.
  * Vermersch et al. [2019] B. Vermersch, A. Elben, L. M. Sieberer, N. Y. Yao, and P. Zoller.  Probing scrambling using statistical correlations between randomized measurements.  _Phys. Rev. X_ , 9:021061, Jun 2019.  URL https://doi.org/10.1103/PhysRevX.9.021061.
  * O’Donnell and Wright [2016] R. O’Donnell and J. Wright.  Efficient quantum tomography.  In _Proceedings of the Forty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’16, page 899–912, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450341325.  URL https://doi.org/10.1145/2897518.2897544.
  * Haah et al. [2016] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu.  Sample-optimal tomography of quantum states.  In _Proceedings of the Forty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’16, page 913–925, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450341325.  URL https://doi.org/10.1145/2897518.2897585.
  * Aaronson and Rothblum [2019] S. Aaronson and G. N. Rothblum.  Gentle measurement of quantum states and differential privacy.  In _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , STOC 2019, page 322–333, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450367059.  URL https://doi.org/10.1145/3313276.3316378.
  * Huang et al. [2020] H.-Y. Huang, R. Kueng, and J. Preskill.  Predicting many properties of a quantum system from very few measurements.  _Nature Physics_ , Jun 2020.  ISSN 1745-2481.  URL https://doi.org/10.1038/s41567-020-0932-7.
  * Cheng et al. [2018] S. Cheng, J. Chen, and L. Wang.  Information perspective to probabilistic modeling: Boltzmann machines versus born machines.  _Entropy_ , 20(8), 2018.  URL https://doi.org/10.3390/e20080583.
  * Benedetti et al. [2019] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini.  Parameterized quantum circuits as machine learning models.  _Quantum Science and Technology_ , 4(4):043001, nov 2019.  URL https://doi.org/10.1088/2058-9565/ab4eb5.
  * Cahill and Glauber [1999] K. E. Cahill and R. J. Glauber.  Density operators for fermions.  _Phys. Rev. A_ , 59:1538–1555, Feb 1999.  URL https://doi.org/10.1103/PhysRevA.59.1538.
  * Li [2013] Y. Li.  Talk 5 : Averages over the unitary group.  In _Selected topics in Mathematical Physics : Quantum Information Theory_ , 2013.
  * Ambainis and Emerson [2007] A. Ambainis and J. Emerson.  Quantum t-designs: t-wise independence in the quantum world.  In _Proceedings of the Twenty-Second Annual IEEE Conference on Computational Complexity_ , CCC ’07, pages 129–140, Washington, DC, USA, 2007\. IEEE Computer Society.  ISBN 0-7695-2780-9.  URL https://doi.org/10.1109/CCC.2007.26.

##  Appendix A One qubit

Based on group-theoretical considerations, it was shown in Refs. [10, 11] that
the density operator on a Hilbert space of spin-$s$ particle (a $2s+1$
dimensional space) can be written as

$\displaystyle\rho=\sum_{m_{s}=-s}^{s}\int_{\Sigma}\frac{d\vec{n}}{4\pi}p(m_{s},\vec{n})K(m_{s},\vec{n}),$
(33)

with the unitary spherical surface ${\Sigma}$ as integration domain,
$\vec{n}=(\cos\varphi\sin\vartheta,\sin\varphi\sin\vartheta,\cos\vartheta)$,
with $\vartheta\in[0,\pi],\varphi\in[0,2\pi)$, $p(m,\vec{n})$ probability of
obtaining the result $m_{s}\in[-s,s]$ when measuring the spin in the
eigenbasis of $\vec{s}\cdot\vec{n}$, where $\vec{s}=(s_{x},s_{y},s_{z})$ is
spin of the system, and $K(m_{s},\vec{n})$ is a kernel function given by

$\displaystyle
K(m_{s},\vec{n})=\frac{2s+1}{\pi}\int_{0}^{2\pi}d\psi\sin^{2}\frac{\psi}{2}e^{i\psi(m_{s}-\vec{s}\cdot\vec{n})}.$
(34)

Specialising Eq. (34) for a spin-$1/2$ particle where $m_{s}=\pm 1/2\doteq
m/2$ and $\vec{s}=(\sigma_{x},\sigma_{y},\sigma_{z})/2\doteq\vec{\sigma}/2$ we
can write

$K_{1}(m,\vec{n})=\frac{2}{\pi}\int_{0}^{2\pi}d\psi\sin^{2}\frac{\psi}{2}e^{i\frac{\psi}{2}(m-\vec{\sigma}\cdot\vec{n})},$
(35)

where the subscript of $K_{1}$ is used to emphasise that is the kernel
operator on a single qubit. We calculate (35) utilising the identity
$e^{-i\frac{\psi}{2}\vec{\sigma}\cdot\vec{n}}=\mathbb{1}\cos{\frac{\psi}{2}}-i\,\vec{\sigma}\cdot\vec{n}\,\sin{\frac{\psi}{2}}$
and directly obtain

$K_{1}(m,\vec{n})=\frac{1}{2}(\mathbb{1}+3\,m\,\vec{\sigma}\cdot\vec{n}).$
(36)

In the main text we have defined

$R_{1}[\sigma_{\alpha}](m,\vec{n})\doteq\text{Tr}[\sigma_{\alpha}K_{1}(m,\vec{n})]=3\,m\,n_{\alpha},$
(37)

an unbiased estimator of $\sigma_{\alpha}$, with $\alpha=x,y,z$ and
$n_{\alpha}$ defined by $\vec{n}=(n_{x},n_{y},n_{z})$. The variance of
$R_{1}[\sigma_{\alpha}]$ is

$\displaystyle\text{Var}(R_{1}[\sigma_{\alpha}](m,\vec{n}))=\langle{R}_{1}^{2}[\sigma_{\alpha}](m,\vec{n})\rangle-\langle{R}_{1}[\sigma_{\alpha}](m,\vec{n})\rangle^{2}\leqslant\langle{R}_{1}^{2}[\sigma_{\alpha}](m,\vec{n})\rangle.$
(38)

We calculate the last term of (38) utilising (37):

$\displaystyle\langle{R}_{1}[\sigma_{\alpha}](m,\vec{n}){R}_{1}[\sigma_{\beta}](m,\vec{n})\rangle$
$\displaystyle=\int_{\Sigma}\frac{d\vec{n}}{4\pi}\sum_{m=\pm
1}p(m,\vec{n}){R}_{1}[\sigma_{\alpha}](m,\vec{n}){R}_{1}[\sigma_{\beta}](m,\vec{n})$
(39)
$\displaystyle=9\int_{\Sigma}\frac{d\vec{n}}{4\pi}\,n_{\alpha}\,n_{\beta}\sum_{m=\pm
1}p(m,\vec{n})=9\int_{\Sigma}\frac{d\vec{n}}{4\pi}\,n_{\alpha}\,n_{\beta}=3\,\delta_{\alpha\beta},$

with $\delta_{\alpha\beta}$ is Kronecker delta. To calculate the variance of
$R_{1}[O](m,\vec{n})$, we expand $O$ in the Pauli operator basis
$O=\sum_{i=0}^{3}a_{i}\,\sigma_{i}$, so that

$\displaystyle\text{Var}(R_{1}[O](m,\vec{n}))$
$\displaystyle=\text{Var}\bigg{(}a_{o}\,R_{1}[\mathbb{1}](m,\vec{n})+\sum_{i=1}^{3}a_{i}\,R_{1}[\sigma_{i}](m,\vec{n})\bigg{)}$
(40)
$\displaystyle=a_{0}^{2}\,\underbrace{\text{Var}(R_{1}[\mathbb{1}](m,\vec{n}))}_{=0}+\text{Var}\bigg{(}\sum_{i=1}^{3}a_{i}\,R_{1}[\sigma_{i}](m,\vec{n})\bigg{)}$
$\displaystyle+2\,\underbrace{\text{Cov}\bigg{(}R_{1}[\mathbb{1}](m,\vec{n}),\sum_{i=1}^{3}a_{i}\,R_{1}[\sigma_{i}](m,\vec{n})\bigg{)}}_{=0}$
$\displaystyle=\text{Var}\bigg{(}\sum_{i=1}^{3}a_{i}\,R_{1}[\sigma_{i}](m,\vec{n})\bigg{)}\leqslant\bigg{\langle}\bigg{(}\sum_{i=1}^{3}a_{i}\,R_{1}[\sigma_{i}](m,\vec{n})\bigg{)}^{\\!2}\,\bigg{\rangle}$
$\displaystyle=\sum_{i=1}^{3}a_{i}^{2}\,\langle
R_{1}^{2}[\sigma_{i}](m,\vec{n})\rangle+\sum_{i\neq j(\neq
0)}\,a_{i}\,a_{j}\langle
R_{1}[\sigma_{i}](m,\vec{n})R_{1}[\sigma_{j}](m,\vec{n})\rangle$
$\displaystyle=3\sum_{i=1}^{3}a_{i}^{2}+3\sum_{i\neq j(\neq
0)}\,a_{i}\,a_{j}\,\delta_{ij}=3\sum_{i=1}^{3}a_{i}^{2}.$

We define a seminorm of $O$ as

$\displaystyle\|O\|=\sqrt{3\sum_{i=1}^{3}a_{i}^{2}}$ (41)

and obtain

$\displaystyle\text{Var}({R}_{1}[O](m,\vec{n}))\leqslant{\|O\|}^{2}.$ (42)

Therefore, for $M$ sufficiently large, $\langle{O}\rangle_{M}$ has a normal
statistical error smaller than or equal to $\|O\|/{\sqrt{M}}$.

##  Appendix B $N$ qubits

For a system of $N$ qubits, we expand the generic observable $O$ as in (19)

$\displaystyle O=\sum_{i}\,a_{i}\,P_{i}$ (43)

where $P_{i}=\sigma_{i_{1}}\otimes\cdots\otimes\sigma_{i_{N}}$,
$i=(i_{1},\ldots\,i_{N})$, is a Pauli monomial operator and
$i_{k}=0,\ldots,3$, and proceed in the calculation of the variance of $R[O]$
in a similar way to the one qubit case (the dependence of $R[\cdot]$ on
$(m_{1},\vec{n}_{1},\ldots,m_{N},\vec{n}_{N})$ is not explicitly indicated for
a lighter notation):

$\displaystyle\text{Var}(R[O])$
$\displaystyle=a_{0}^{2}\,\text{Var}(R[\mathbb{1}])+\text{Var}\bigg{(}\sum_{i\neq
0_{v}}a_{i}\,R[P_{i}]\bigg{)}+2\,\text{Cov}\bigg{(}R[\mathbb{1}],\sum_{i\neq
0_{v}}a_{i}\,R[P_{i}]\bigg{)}$ (44)
$\displaystyle=\text{Var}\bigg{(}\sum_{i\neq
0_{v}}a_{i}\,R[P_{i}]\bigg{)}\leqslant\bigg{\langle}\bigg{(}\sum_{i\neq
0_{v}}a_{i}\,R[P_{i}]\bigg{)}^{\\!2}\,\bigg{\rangle}=\sum_{i,j(\neq
0_{v})}\,a_{i}\,a_{j}\langle R[P_{i}]R[P_{j}]\rangle,$

where the last sum on $i$ and $j$ is extended to all values of $i$ and $j$
except $0_{v}\doteq(0,0,\ldots,0)$. In general, for every $k$, the indices
$i_{k}$ and $j_{k}$ appearing in each $\langle R[P_{i}]R[P_{j}]\rangle$ of
(44) can assume any integer value between 0 and 3. For given $i$ and $j$, let
$S$ be the set of indices $k$ for which $i_{k}\neq 0\wedge j_{k}\neq 0$, $s$
the cardinality of $S$ and $\overline{S}$ the set containing the remaining
indices $k$ of ${1,\ldots,N}$ with cardinality ${\bar{s}}=N-s$. We use (16)
for a Pauli monomial to arrive at

$\displaystyle\langle R[P_{i}]R[P_{j}]\rangle$
$\displaystyle=\bigg{(}\prod_{{h}\in{\overline{S}}}\int_{{\Sigma}_{h}}\frac{d{\vec{n}}_{h}}{4\pi}\sum_{m_{h}}{R}_{1}[\sigma_{i_{h}}](m_{i_{h}},\vec{n}_{i_{h}}){R}_{1}[\sigma_{j_{h}}](m_{j_{h}},\vec{n}_{j_{h}})\bigg{)}$
(45)
$\displaystyle\times\bigg{(}\prod_{k\in{S}}\int_{{\Sigma}_{k}}\frac{d{\vec{n}}_{k}}{4\pi}\sum_{m_{k}}{R}_{1}[\sigma_{i_{k}}](m_{i_{k}},\vec{n}_{i_{k}}){R}_{1}[\sigma_{j_{k}}](m_{j_{k}},\vec{n}_{j_{k}})\bigg{)}p(\\{m,\vec{n}\\})$
$\displaystyle=\bigg{(}\prod_{{h}\in{\overline{S}}}\int_{{\Sigma}_{h}}\frac{d{\vec{n}}_{h}}{4\pi}\sum_{m_{h}}{R}_{1}[\sigma_{i_{h}}](m_{i_{h}},\vec{n}_{i_{h}}){R}_{1}[\sigma_{j_{h}}](m_{j_{h}},\vec{n}_{j_{h}})\bigg{)}$
$\displaystyle\times\bigg{(}9^{s}\prod_{k\in{S}}\int_{{\Sigma}_{k}}\frac{d{\vec{n}}_{k}}{4\pi}\,{n}_{i_{k}}{n}_{j_{k}}\bigg{)}\bigg{(}\prod_{k\in{S}}\sum_{m_{k}}p(\\{m,\vec{n}\\})\bigg{)}$
$\displaystyle=\bigg{(}\prod_{{h}\in{\overline{S}}}\int_{{\Sigma}_{h}}\frac{d{\vec{n}}_{h}}{4\pi}\sum_{m_{h}}{R}_{1}[\sigma_{i_{h}}](m_{i_{h}},\vec{n}_{i_{h}}){R}_{1}[\sigma_{j_{h}}](m_{j_{h}},\vec{n}_{j_{h}})\bigg{)}$
$\displaystyle\times\bigg{(}3^{s}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\bigg{)}p(m_{h_{1}},\vec{n}_{h_{1}},\ldots,m_{h_{\bar{s}}},\vec{n}_{h_{\bar{s}}})$
$\displaystyle=\bigg{(}\prod_{{h}\in{\overline{S}}}\int_{{\Sigma}_{h}}\frac{d{\vec{n}}_{h}}{4\pi}\sum_{m_{h}}{R}_{1}[\sigma_{i_{h}}](m_{i_{h}},\vec{n}_{i_{h}}){R}_{1}[\sigma_{j_{h}}](m_{j_{h}},\vec{n}_{j_{h}})\bigg{)}$
$\displaystyle\times\bigg{(}3^{s}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\int_{{\Sigma}_{k}}\frac{d{\vec{n}}_{k}}{4\pi}\sum_{m_{k}}{R}_{1}[\sigma_{i_{k}}^{2}]\bigg{)}p(\\{m,\vec{n}\\})$
$\displaystyle=3^{s}\bigg{(}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\bigg{)}\bigg{(}\prod_{{h}\in{\overline{S}}}\int_{{\Sigma}_{h}}\frac{d{\vec{n}}_{h}}{4\pi}\sum_{m_{h}}{R}_{1}[\sigma_{i_{h}}\sigma_{j_{h}}](m_{i_{h}},\vec{n}_{i_{h}},m_{j_{h}},\vec{n}_{j_{h}})\bigg{)}$
$\displaystyle\times\bigg{(}\prod_{k\in{S}}\int_{{\Sigma}_{k}}\frac{d{\vec{n}}_{k}}{4\pi}\sum_{m_{k}}{R}_{1}[\sigma_{i_{k}}\sigma_{j_{k}}](m_{i_{k}},\vec{n}_{i_{k}},m_{j_{k}},\vec{n}_{j_{k}})\bigg{)}p(\\{m,\vec{n}\\})$
$\displaystyle=3^{s}\bigg{(}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\bigg{)}\langle
R[P_{i}P_{j}]\rangle=3^{s}\bigg{(}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\bigg{)}\text{Tr}[P_{i}P_{j}\rho]\leqslant
3^{s}\bigg{(}\prod_{k\in{S}}\delta_{{i_{k}}{j_{k}}}\bigg{)}\doteq
3^{r_{ij}}\Delta_{ij}.$

If we introduce a seminorm of $O$ as

$\displaystyle\|O\|=\sqrt{\sum_{i,j(\neq
0_{v})}3^{r_{ij}}\Delta_{ij}|a_{i}||a_{j}|}\,,$ (46)

Eq. (44) becomes

$\displaystyle\text{Var}({R}[O](\\{m,\vec{n}\\}))\leqslant{\|O\|}^{2}.$ (47)

The seminorm (46) is a tighter upper bound for the variance of $R[O]$ than the
formally simpler seminorm

$\displaystyle{\|O\|}_{1}=\sum_{i\neq 0_{v}}\sqrt{3^{r_{i}}}|a_{i}|,$ (48)

with $r_{i}$ representing the weight of $P_{i}$, since

$\displaystyle\|O\|^{2}=\sum_{i,j(\neq
0_{v})}3^{r_{ij}}\Delta_{ij}|a_{i}||a_{j}|\leqslant\sum_{i,j(\neq
0_{v})}\sqrt{3^{r_{i}}3^{r_{j}}}|a_{i}||a_{j}|=\Bigl{(}\sum_{i\neq
0_{v}}\sqrt{3^{r_{i}}}|a_{i}|\Bigr{)}^{2}={\|O\|}_{1}^{2}$ (49)

This is not surprising, since the seminorm (48) can be shown to be a bound for
the variance of $R[O]$ by expanding $\text{Var}\bigl{(}\sum_{i\neq
0_{v}}a_{i}\,R[P_{i}]\bigr{)}$ of (44) and using Schwarz inequality without
resorting to explicit properties of $R$ as in (45). The seminorm (46) can also
be expressed as

$\displaystyle\|O\|^{2}={\|O\|}_{2}^{2}+\sum_{i\neq j(\neq
0_{v})}3^{r_{ij}}\Delta_{ij}|a_{i}||a_{j}|,$ (50)

with the seminorm ${\|O\|}_{2}$ defined by

$\displaystyle{\|O\|}_{2}=\sqrt{\sum_{i\neq 0_{v}}3^{r_{i}}\,a_{i}^{2}}.$ (51)

We use (45) and rewrite (44) as

$\displaystyle\text{Var}(R[O])\leqslant\sum_{i,j(\neq
0_{v})}\,3^{r_{ij}}\,\Delta_{ij}\,\text{Tr}[P_{i}P_{j}\rho]\,a_{i}\,a_{j}={\|O\|}_{2}^{2}+\sum_{i\neq
j(\neq
0_{v})}\,3^{r_{ij}}\,\Delta_{ij}\,\text{Tr}[P_{i}P_{j}\rho]\,a_{i}\,a_{j}.$
(52)

Eq. (52) shows how the mixed terms of (50) are a weak upper bound for the last
term of (52). For a randomly chosen $O$ or, in the same way, for a random
$\rho$, the addends in the last term of (52) have mixed positive and negative
signs. Furthermore, for most states $\rho$, the expectation values of the
expressions $P_{i}P_{j}$ tend to become smaller for higher values of $r_{ij}$.
As a result, ${\|O\|}_{2}$ is generally a good approximation of the standard
deviation of $R[O]$.

To formalise this statement, we start by showing that the expectation value of
the mixed terms of (52) over the space of the density operators of the system
is equal to zero. The probability of choosing a density operator from the
space of density operators is the probability of choosing a set of eigenvalues
multiplied by the probability of choosing a density operator with specific
eigenvalues conditional to the choice of the eigenvalues. If for any choice of
the eigenvalues, the expectation value over the density operators with the
chosen eigenvalues is zero, then the expectation value over all density
operators is zero. All density operators with the same eigenvalues can be
obtained from one with a unitary transformation and, vice versa, all unitary
transformations applied to a density operator preserve the eigenvalues.
Therefore, we can calculate the expectation value of the mixed terms of (52)
over the density operators with given eigenvalues by averaging over the
unitary group [23] with the unitary transformations applied to an arbitrary
density operator $\rho_{0}$

$\displaystyle\int_{U}dU\sum_{i\neq j(\neq
0_{v})}3^{r_{ij}}\,\Delta_{ij}\,a_{i}\,a_{j}\,\text{Tr}[P_{i}P_{j}\,U\rho_{0}\,U^{\dagger}]=\sum_{i\neq
j(\neq
0_{v})}3^{r_{ij}}\,\Delta_{ij}\,a_{i}\,a_{j}\,\text{Tr}\bigg{[}P_{i}P_{j}\int_{U}dU\,U\rho_{0}\,U^{\dagger}\bigg{]}$
(53) $\displaystyle=\sum_{i\neq j(\neq
0_{v})}3^{r_{ij}}\,\Delta_{ij}\,a_{i}\,a_{j}\,\text{Tr}\Bigl{[}\frac{P_{i}P_{j}}{2^{N}}\Bigr{]}=0,$

which implies

$\displaystyle\bigg{\langle}\sum_{i\neq j(\neq
0_{v})}\,3^{r_{ij}}\,\Delta_{ij}\,\text{Tr}[P_{i}P_{j}\rho]\,a_{i}\,a_{j}\bigg{\rangle}_{\\!\\!\rho}=0.$
(54)

To gain further information on how close ${\|O\|}_{2}$ is likely to be to
${\|O\|}$ for a random state, we now consider the variance of the mixed terms
of (50) over the space of the density operators of the system. For simplicity,
we specialise the analysis to pure states. From the unitary group notation of
(53), we switch to the average over the states expressed as in [24] and recall
the following identities:

$\displaystyle\int_{\psi}{d\psi}|\psi\rangle\\!\langle\psi|^{\otimes
2}=\frac{2}{2^{N}(2^{N}+1)}\,\mathbb{1}_{\rm sym},$ (55)

with $\mathbb{1}_{\rm sym}$ projector onto the symmetric subspace of ${\cal
H}^{\otimes 2}$, and

$\displaystyle\text{Tr}[A\otimes B\,\mathbb{1}_{\rm
sym}]=\frac{1}{2}\big{(}\text{Tr}[A]\text{Tr}[B]+\text{Tr}[AB]\big{)},$ (56)

for any operator $A$ and $B$. The variance of the mixed terms of (52) over the
states of the system can be expressed as

$\displaystyle\text{Var}\bigg{(}\sum_{i\neq j(\neq
0_{v})}\,3^{r_{ij}}\,\Delta_{ij}\,\text{Tr}\big{[}P_{i}P_{j}|\psi\rangle\\!\langle\psi|\big{]}\,a_{i}\,a_{j}\bigg{)}_{\\!\psi}=\int_{\psi}{d\psi}\bigg{(}\sum_{i\neq
j(\neq
0_{v})}a_{i}\,a_{j}\,3^{r_{ij}}\,\Delta_{ij}\,\text{Tr}\big{[}P_{i}P_{j}|\psi\rangle\\!\langle\psi|\big{]}\bigg{)}^{\\!2}$
(57) $\displaystyle=\sum_{i\neq j(\neq 0_{v})}\sum_{i^{\prime}\neq
j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\int_{\psi}{d\psi}\,\big{(}\text{Tr}\big{[}P_{i}P_{j}|\psi\rangle\\!\langle\psi|\big{]}\,\text{Tr}\big{[}P_{i^{\prime}}P_{j^{\prime}}|\psi\rangle\\!\langle\psi|\big{]}\big{)}$
$\displaystyle=\sum_{i\neq j(\neq 0_{v})}\sum_{i^{\prime}\neq j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\,\text{Tr}\bigg{[}(P_{i}P_{j})\otimes(P_{i^{\prime}}P_{j^{\prime}})\int_{\psi}{d\psi}|\psi\rangle\\!\langle\psi|^{\otimes
2}\bigg{]}$ $\displaystyle=\frac{2}{2^{N}(2^{N}+1)}\sum_{i\neq j(\neq
0_{v})}\sum_{i^{\prime}\neq j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\,\text{Tr}\big{[}(P_{i}P_{j})\otimes(P_{i^{\prime}}P_{j^{\prime}})\,\mathbb{1}_{\rm
sym}\big{]}$ $\displaystyle=\frac{1}{2^{N}(2^{N}+1)}\sum_{i\neq j(\neq
0_{v})}\sum_{i^{\prime}\neq j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\,\text{Tr}[P_{i}P_{j}P_{i^{\prime}}P_{j^{\prime}}]$
$\displaystyle=\frac{1}{(2^{N}+1)}\sum_{i\neq j(\neq
0_{v})}\sum_{i^{\prime}\neq j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\,\overline{\Delta}_{iji^{\prime}j^{\prime}},$

with $\overline{\Delta}_{iji^{\prime}j^{\prime}}=1$ when
$i,j,i^{\prime},j^{\prime}$ are such that
$P_{i}P_{j}P_{i^{\prime}}P_{j^{\prime}}=\mathbb{1}$, otherwise
$\overline{\Delta}_{iji^{\prime}j^{\prime}}=0$. A general upper bound of the
last expression is likely to be weak, given its strong dependence on $O$. The
specific case of $O$ given by a projector onto a computational basis state is
examined in App. C.

##  Appendix C Projectors onto the computational basis states

In this section we examine properties of the seminorm for an important special
case, namely, the case of projectors onto a generic element of the
computational basis state,

$\displaystyle|x\rangle\langle
x|=\bigotimes_{k=1}^{N}\,\frac{1}{2}(\mathbb{1}\pm\sigma_{z_{k}}),$ (58)

where the $\pm$ of each factor is defined by the choice of the element of the
computational basis $|x\rangle$. This case is important not only for its
practical significance in a machine learning context, as briefly indicated in
Sec. 5, but also because it is representative of the cases with bounded
${\|O\|}_{2}$ and diverging ${\|O\|}$ for increasing $N$, which makes the
determination of whether the statistical error of the estimator can be
approximated by ${\|O\|}_{2}$ of particular interest. To calculate the
seminorms of any projector onto a computational basis state, we expand the
tensor product of (58) and use the seminorm definitions. For
${\||x\rangle\langle x|\|}_{2}$ this simply means

$\displaystyle{\||x\rangle\langle
x|\|}_{2}^{2}=\frac{1}{4^{N}}\bigg{(}{{N}\choose{N-1}}\,3^{1}+{{N}\choose{N-2}}\,3^{2}+\ldots+{{N}\choose{0}}\,3^{N}\bigg{)}=\frac{1}{4^{N}}\big{(}(1+3)^{N}\big{)}-\frac{1}{4^{N}}=1-\frac{1}{4^{N}}.$
(59)

A more involved combinatorial calculation would show that ${\||x\rangle\langle
x|\|}^{2}\leqslant(3/2)^{N}$ and ${\||x\rangle\langle x|\|}^{2}=O((3/2)^{N})$.
We come back to (57) for any projector onto computational basis state and
observe that, for given $i$ and $j$, there are $2^{N}$ ordered pairs
$i^{\prime},j^{\prime}$ such that
$\overline{\Delta}_{iji^{\prime}j^{\prime}}=1$. The operator $P_{i}P_{j}$ is
in fact the tensor product of single–qubit identities and operators
$\sigma_{z}$. Each of them becomes the identity when multiplied by itself and
this can be obtained with two choices for every qubit of $P_{i^{\prime}}$ and
$P_{j^{\prime}}$. Unfortunately, to complicate things,
${r_{i^{\prime}j^{\prime}}}$ is not independent of the choice of $i^{\prime}$
and $j^{\prime}$ for given $i$ and $j$; however, we can simplify the
evaluation of (57) by considering a limit case and determining an upper bound.
In particular, for every $i$ and $j$, with $i\neq j$ and $i$ and $j$ different
from $0_{v}$,

$\displaystyle\sum_{i^{\prime},j^{\prime}}\,3^{r_{i^{\prime}j^{\prime}}}\overline{\Delta}_{iji^{\prime}j^{\prime}}<\sum_{i^{\prime},j^{\prime}}\,3^{r_{i^{\prime}j^{\prime}}}\overline{\Delta}_{iii^{\prime}j^{\prime}}=\sum_{l=0}^{N}{{N}\choose{l}}\,3^{l}=4^{N}.$
(60)

With the substitutions $a_{i}=(1/2)^{N}$ and $\Delta_{ij}=1$ valid for any
projector $|x\rangle\langle x|$ and any value of the indexes, the last term of
Eq. (57) becomes

$\displaystyle\frac{1}{(2^{N}+1)}\sum_{i\neq j(\neq
0_{v})}\sum_{i^{\prime}\neq j^{\prime}(\neq
0_{v})}a_{i}\,a_{j}\,a_{i^{\prime}}a_{j^{\prime}}\,3^{r_{ij}}\,3^{r_{i^{\prime}j^{\prime}}}\,\Delta_{ij}\,\Delta_{i^{\prime}j^{\prime}}\,\overline{\Delta}_{iji^{\prime}j^{\prime}}$
(61) $\displaystyle<\frac{1}{(2^{N}+1)\,4^{N}}\sum_{i\neq j(\neq
0_{v})}3^{r_{ij}}=\frac{{\||x\rangle\langle
x|\|}^{2}}{2^{N}+1}<\bigg{(}\frac{3}{4}\bigg{)}^{\\!N}.$

Eq. (61) shows that, for any projector $|x\rangle\langle x|$,
$\text{Var}(R[|x\rangle\langle x|])$ for a random state is bounded by
${\||x\rangle\langle x|\|}_{2}^{2}<1$ with probability exponentially close to
1 as a function of $N$.