remarkRemark hypothesisHypothesis claimClaim Operator Augmentation for General
Noisy Matrix SystemsP.A. Etter, L. Ying

# Operator Augmentation for General Noisy Matrix Systems

Philip A. Etter () Institute for Computational and Mathematical Engineering,
Stanford University<EMAIL_ADDRESS>Lexing Ying () Department of
Mathematics, Stanford University<EMAIL_ADDRESS>
###### Abstract

In the computational sciences, one must often estimate model parameters from
data subject to noise and uncertainty, leading to inaccurate results. In order
to improve the accuracy of models with noisy parameters, we consider the
problem of reducing error in a linear system with the operator corrupted by
noise. To address this problem, we extend the elliptic _operator augmentation_
framework [4] to the general nonsymmetric matrix case. We show that under the
conditions of right-hand-side isotropy and noise symmetry that the optimal
operator augmentation factor for the residual error is always positive,
thereby making the framework amenable to a necessary bootstrapping step. While
the above conditions are unnecessary for positive optimal augmentation factor
in the elliptic case, we provide counter-examples that illustrate their
necessity when applied to general matrices. When the noise in the operator is
small, however, we show that the condition of noise symmetry is unnecessary.
Finally, we demonstrate through numerical experiments on Markov chain problems
that operator augmentation can significantly reduce error in noisy matrix
systems — even when the aforementioned conditions are not met.

###### keywords:

Operator Augmentation, Random Matrices, Monte Carlo, Polynomial Expansion,
Nonsymmetric Matrices.

Linear and multilinear algebra; matrix theory. Statistics. Computer Science.

##  1 Introduction

There are a plethora of different situations in the natural, mathematical, and
computer sciences that necessitate computing the solution to a linear system
of equations given by

(1) $\mathbf{A}\mathbf{x}=\mathbf{b}\,,$

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ and
$\mathbf{x},\mathbf{b}\in\mathbb{R}^{n}$ for $n\in\mathbb{N}$. When both the
matrix $\mathbf{A}$ and $\mathbf{b}$ are known, there are many decades of
research on how to solve the system Eq. 1 efficiently. Unfortunately, for a
variety of reasons, it is often the case that the true matrix $\mathbf{A}$ is
not known exactly, and must be estimated from data (see [9, 8]). In this
situation, there is a sampling error between the unobserved true matrix
$\mathbf{A}$ and the matrix $\hat{\mathbf{A}}$ one constructs from data. The
discrepancy between $\mathbf{A}$ and $\hat{\mathbf{A}}$ is often referred to
as _model uncertainty_ , as it stems from incomplete or inaccurate information
about the underlying system. This model uncertainty means that with naive
application of the inverse of the observed matrix $\hat{\mathbf{A}}$, one is
not solving the desired system Eq. 1, but rather, the system

(2) $\hat{\mathbf{A}}\hat{\mathbf{x}}=\mathbf{b}\,,$

where $\hat{\mathbf{x}}=\hat{\mathbf{A}}^{-1}\mathbf{b}\in\mathbb{R}^{n}$ is
the solution we observe when we solving the observed system naively. Often, we
will write

(3) $\hat{\mathbf{A}}=\mathbf{A}+\hat{\mathbf{Z}}\,,$

where one can think of the matrix $\hat{\mathbf{Z}}$ as constituting the noise
or sampling error in our measurements of the system Eq. 1. Hence, the sampling
error $\hat{\mathbf{Z}}$ between $\mathbf{A}$ and $\hat{\mathbf{A}}$
translates into sampling error between the true solution $\mathbf{x}$ and the
naively estimated solution $\hat{\mathbf{x}}$.

The question of interest in this paper is whether one can increase the
accuracy of the solution $\hat{\mathbf{x}}$ vis-á-vis the true solution
$\mathbf{x}$ by modifying the sampled operator $\hat{\mathbf{A}}$. The paper
[4] proposes the operator augmentation framework; whereby one _augments_ the
operator $\hat{\mathbf{A}}^{-1}$ by an appropriately chosen augmentation
matrix $\hat{\mathbf{K}}$,

(4) $\hat{\mathbf{A}}^{-1}-\beta\hat{\mathbf{K}}\,.$

Solving the system of equations with this augmented inverse operator instead
of the naive operator $\hat{\mathbf{A}}^{-1}$, yields an augmented estimate
$\tilde{\mathbf{x}}_{\beta}$ of $\mathbf{x}$,

(5)
$\tilde{\mathbf{x}}_{\beta}=(\hat{\mathbf{A}}^{-1}-\beta\hat{\mathbf{K}})\mathbf{b}\,.$

Once $\hat{\mathbf{K}}$ has been chosen, obtaining an optimal value
$\beta^{*}$ for the augmentation factor $\beta$ involves solving the quadratic
optimization problem,

(6)
$\beta^{*}=\operatorname*{arg\,min}_{\beta}\mathbb{E}\|\tilde{\mathbf{x}}_{\beta}-\mathbf{x}\|^{2}_{\mathbf{B}}=\operatorname*{arg\,min}_{\beta}\mathbb{E}\|(\hat{\mathbf{A}}^{-1}-\beta\hat{\mathbf{K}})\mathbf{b}-\mathbf{A}^{-1}\mathbf{b}\|^{2}_{\mathbf{B}}\,,$

for some positive definite norm matrix $\mathbf{B}$. Of course, the true value
of $\beta^{*}$ is not possible to obtain, by virtue of the fact that we do not
have access to the matrix $\mathbf{A}$. Therefore, operator augmentation
attempts to bootstrap $\beta^{*}$ by substituting $\hat{\mathbf{A}}$ for
$\mathbf{A}$ and instead solving

(7)
$\tilde{\beta}^{*}=\operatorname*{arg\,min}_{\beta}\mathbb{E}\|(\hat{\mathbf{A}}_{b}^{-1}-\beta\hat{\mathbf{K}})\mathbf{b}-\hat{\mathbf{A}}^{-1}\mathbf{b}\|^{2}_{\mathbf{B}}\,,$

where now the expectation $\mathbb{E}$ is over a bootstrapped noisy matrix
$\hat{\mathbf{A}}_{b}$, sampled from a distribution we estimate from
$\hat{\mathbf{A}}$. This minimization problem can now be solved with Monte
Carlo, although one might opt to use specially truncated Taylor expansions of
$\hat{\mathbf{A}}_{b}^{-1}$ about $\hat{\mathbf{A}}^{-1}$ to accelerate
computation [4].

When $\mathbf{A}$ is a positive definite symmetric operator, one can prove for
several choices of augmentation matrix $\hat{\mathbf{K}}$, that, under a wide
range of norms and noise distributions satisfying
$\mathbb{E}[\hat{\mathbf{A}}]\preceq\mathbf{A}$, the optimal augmentation
factor $\beta^{*}$ will always be positive [4]. In effect, optimal operator
augmentation always _shrinks_ the noisy matrix $\hat{\mathbf{A}}^{-1}$ towards
some appropriate base point. The intuition behind this comes from Jensen’s
inequality (see Fig. 1). Matrix inversion operation is a convex function (wrt.
the Löwner Order) on the set of positive definite matrices
$S(\mathbb{R}^{n})_{+}$. This means that naively inverting an unbiased
estimator $\hat{\mathbf{A}}$ for $\mathbf{A}$ will always overshoot the true
goal $\mathbf{A}^{-1}$. One can think of operator augmentation as correcting
for this overshooting effect. This is similar to the intuition behind Stein
estimation [5], which makes the observation that in high dimensions, the
standard Pitman estimator is always inadmissible because appropriately
shrinking the Pitman estimate towards the origin produces an estimator with
strictly better risk. The guarantee that $\beta^{*}>0$ is desirable because
amplifying the noise in the matrix $\hat{\mathbf{A}}$ could potentially
produce unstable results, particularly because $\tilde{\beta}^{*}$ is
estimated from a bootstrap.

Figure 1: An example of the overshooting effect. If we take a single sample of
the of the scalar random variable $X\sim\Gamma(2,1/2)$, and invert it, the pdf
of the inverted $1/X$ has an expectation that is significantly larger (x2)
than the inversion $1/\mathbb{E}[X]$. This means that naively trying to
estimate $1/\mathbb{E}[X]$ with only a single sample will likely give a
significant overestimate. The same principle can apply when $X$ is a random
matrix.

However, while [4] provides theory and numerical experiments for the space of
elliptic operators $S(\mathbb{R}^{n})_{+}$, it does not provide any theory or
experiments for general matrix systems in $GL(\mathbb{R}^{n})$111The general
linear group over $\mathbb{R}^{n}$. In this paper, we attempt to fill this gap
by studying the operator augmentation algorithm on the full generality of
nonsingular matrices in $GL(\mathbb{R}^{n})$. Immediately, the theoretical
task of analyzing the algorithm is made more difficult by the fact that
$GL(\mathbb{R}^{n})$ doesn’t have the same convexity structure relative to the
matrix inversion operation $\mathbf{A}^{-1}$ that the SPD matrices
$S(\mathbb{R}^{n})_{+}$ do. Despite this obstacle, we are able to provide
error reduction guarantees in the residual norm under the conditions of noise
symmetry and right-hand-side isotropy — and we provide proofs of these
guarantees in Section 4. Furthermore, we show the existence of adversarial
examples that demonstrate the necessity of the aforementioned conditions and
thereby lead us to a better understanding of the phenomena that may cause
operator augmentation to fail in the nonsymmetric case. We provide these case
studies in Section 5. Fortunately, the noise symmetry assumption can be
discarded if one is willing to assume a small noise regime, as we demonstrate
in Section 6. Finally, we examine more general forms of operator augmentation
and anisotropic right-hand-side in Section 7.

To complete our discussion about nonsymmetric operator augmentation, we
benchmark the technique on value function estimation problems for Markov
Decision Processes (MDPs), and show that despite the limitations of the
nonsymmetric setting, operator augmentation still achieves a substantial
reduction in error. We give an overview of our benchmark problems and
numerical results in Section 9.

##  2 Related Work

As we discussed earlier, the _operator augmentation_ framework for elliptic
problems was originally proposed and analyzed in [4]. However, the paper [4]
only considers the elliptic case and does not study the general matrix case.

The spirit of operator augmentation derives from Stein Estimation [5]. In
their seminal paper, James and Stein consider the high dimensional problem of
estimating a parameter corrupted by noise and show that the standard Pitman
estimator is always inadmissible because appropriately shrinking the Pitman
estimate towards the origin produces an estimator with strictly better risk.
At a fundamental level, one can frame our work as taking this idea and
applying it to the novel setting of matrices corrupted by noise.

A pertinent area of the literature related to our work is _model uncertainty_.
Quantifying and representing model uncertainty is important in many different
fields of computational science, ranging from structural dynamics [9] to
weather and climate prediction [8]. However, work relating to model or
parameter uncertainty is usually domain specific and focuses more on
establishing a model for uncertainty than it does on trying to reduce error in
the resulting predictions. In contrast, our work focuses entirely on reducing
error, rather than quantifying it. Our work is also not restricted to a
particular domain, class of problems, or randomness structure, as long as
those problems are linear.

We note that our setting shares some similarities with the problem of
uncertainty quantification (UQ). However, the problem we face here is
different from the standard uncertainty quantification setting in a subtle but
very important way. In UQ, one is usually given a distribution $\mathcal{P}$
and a map $T$ and asked to estimate statistics about the pushforward
distribution $T_{*}\mathcal{P}$ (i.e., expectation, standard error, etc.).
Practitioners typically accomplish this task via Monte Carlo techniques [7] or
some form of stochastic Galerkin projection [12] or collocation method [11].
In contrast, in our setting we have a distribution
$\mathcal{P}=D_{\omega^{*}}$ and our goal is to accurately compute the image
under $T$ of a statistic of $\mathcal{P}$, rather than a statistic of the
push-forward of $\mathcal{P}$ under $T$. Here $T(\cdot)$ is matrix inversion
and the statistic we are interested in is the true matrix
$s(D_{\omega^{*}})=\mathbf{A}$. Indeed, the issue that
$\mathbb{E}[\hat{\mathbf{A}}^{-1}]$ and $\mathbb{E}[\hat{\mathbf{A}}]^{-1}$
are significantly different is at the core of our paper, as we are (in the
unbiased case) interested in $\mathbb{E}[\hat{\mathbf{A}}]^{-1}$, and not
$\mathbb{E}[\hat{\mathbf{A}}^{-1}]$. What does, however, fall into the purview
of uncertainty quantification is the question of evaluating how well one does
of estimating $\mathbb{E}[\hat{\mathbf{A}}^{-1}]$; but this is ultimately an
orthogonal discussion.

The central problem in this paper is also not dissimilar to the setting of
matrix completion seen in [2, 6]. In matrix completion, one usually seeks to
recover a low-rank ground truth matrix $\mathbf{M}_{ij}$ from observations
that have been corrupted by additive noise, e.g.,
$\mathbf{N}=\mathbf{M}+\mathbf{Z}$. If $\mathcal{P}_{\Omega}$ denotes the
subset sampling operator on matrix space, then one is trying to recover
$\mathbf{M}$ from

(8)
$\mathcal{P}_{\Omega}(\mathbf{N})=\mathcal{P}_{\Omega}(\mathbf{M})+\mathcal{P}_{\Omega}(\mathbf{Z})\,.$

However, the operator augmentation and matrix completion settings are subtly
different. The matrix completion analogue of $\mathbf{A}$ is the actual linear
operator $\mathcal{P}_{\Omega}$, and not the matrix $\mathbf{M}$. Morally, one
may think of the matrix competition problem as solving the under-determined
linear system

(9) $\mathcal{P}_{\Omega}(\mathbf{M})=\mathcal{P}_{\Omega}(\mathbf{N})\,$

by assuming a low-rank regularity on $\mathbf{M}$. The randomness in this
problem lies completely in the right hand side $\mathbf{N}$, and not in the
actual linear operator $\mathcal{P}_{\Omega}$.

We also draw attention to the related field of perturbation matrix analysis.
In this setting, one is usually interested in proving results about how
various properties of matrices change under a perturbation to the elements of
the matrix. A seminal example of work in this field is the Davis-Kahan Theorem
[3], which quantifies the extent to which the invariant sub-spaces of a matrix
change under perturbations. In a similar vein, work in backward stability
analysis revolves around understanding the behavior of the solution of a
linear system under perturbations to the matrix. However, backward stability
analysis typically adopts a worst-case mentality in analysis. In contrast, we
care about average case error — and more importantly, how one can reduce it.

We should briefly mention that the mathematical branch of random matrix theory
(RMT) studies the spectral properties of random matrix ensembles [1, 10].
However, RMT results usually apply only when the entries of the random
matrices are independent and in the large matrix limit. We find these
assumptions to be too stringent for the problem at hand.

In conclusion, we do not believe that the setting we introduce in this paper,
where the operator is noisy, has been studied in the proposed fashion before.
There is little precedent in the literature for the operator augmentation
method beyond the original paper [4].

##  3 Probabilistic Model and Problem Formalism

To provide mathematical clarity to underlying setting, we briefly devote time
to explaining the formalism of the operator augmentation setting. Formally, we
model this process by assuming that there exists an underlying parameter space
$\Omega$ (with sigma algebra $\Sigma$), where the parameters $\omega\in\Omega$
contain a description of the system that produces the matrices above (e.x.,
$\omega$ may be measurements of a scattering background, edge weights, vertex
positions, etc.). The relationship between parameters and matrices is given by
a Borel-measurable mapping $\mathcal{M}:\Omega\longrightarrow
GL(\mathbb{R}^{n})$, where $GL(\mathbb{R}^{n})$ denotes the group of non-
singular matrices in $\mathbb{R}^{n\times n}$. For example, $\omega\in\Omega$
may be a weighted graph, and $\mathcal{M}(\omega)$ may denote a minor of its
Laplacian. We suppose that there exist some _unobserved_ true system
parameters $\omega^{*}\in\Omega$ that produce the true matrix
$\mathbf{A}=\mathcal{M}(\omega^{*})$. We also suppose that there exists a
known family of distributions $\mathbb{P}_{\omega}$ over $\Omega$ indexed by
$\omega\in\Omega$ that describes the observed randomness in the system if
$\omega$ were to be the true system parameters. It is this the relationship
between $\omega^{*}$ and the distribution $\mathbb{P}_{\omega^{*}}$ that we
assume is known as part of the model (but not the true system parameters
$\omega^{*}$ themselves). Once this family has been specified, the
distribution of $\hat{\mathbf{A}}$ is given by
$\mathcal{M}_{*}\mathbb{P}_{\omega^{*}}$, where $\mathcal{M}_{*}$ denotes the
pushforward. For shorthand, we will denote
$D_{\omega^{*}}\equiv\mathcal{M}_{*}\mathbb{P}_{\omega^{*}}$.

Finally, we formally state our problem as follows: suppose that we observe
exactly _one_ sample $\hat{\omega}\sim\mathbb{P}_{\omega^{*}}$. Given
$\mathbf{b}$ or a distribution from which $\mathbf{b}$ is sampled, we wish to
build an inverse operator $\tilde{\mathbf{A}}_{\beta}^{-1}$ such that
$\tilde{\mathbf{x}}=\tilde{\mathbf{A}}_{\beta}^{-1}\mathbf{b}$ produces a
better estimate of $\mathbf{x}=\mathbf{A}^{-1}\mathbf{b}$ than
$\hat{\mathbf{x}}=\hat{\mathbf{A}}^{-1}\mathbf{b}$. In our case, we take

(10)
$\tilde{\mathbf{A}}_{\beta}^{-1}\equiv\hat{\mathbf{A}}^{-1}-\beta\hat{\mathbf{K}}\,,$

for some $\beta\in\mathbb{R}$, $\hat{\mathbf{K}}\in\mathbb{R}^{n\times n}$ to
be specified later. For ease of understanding, we present flow-chart of the
above probabilistic setting together with the operator augmentation framework
in Fig. 2.

Figure 2: A flow-chart of the probabilistic setting of operator augmentation
as well as the algorithm itself. Operator augmentation aims to find a $\beta$
that gives an optimal reduction in error given an augmentation matrix
$\hat{\mathbf{K}}$.

##  4 Provable Error Reduction Guarantees

In order to provide a theory for nonsymmetric operator augmentation that
mirrors the theory for symmetric operator augmentation, we focus on the semi-
Bayesian error as defined for symmetric operator augmentation. For semi-
Bayesian error, we assume that the right-hand-side $\mathbf{b}$ in the
equation

(11) $\mathbf{A}\mathbf{x}=\mathbf{b}\,,$

is drawn at random from a prior distribution $P$. For a given value of
augmentation factor $\beta$, we are interested in the error between the
operator augmentation $\tilde{\mathbf{x}}_{\beta}$ and the true solution
$\mathbf{x}$, averaged over both the matrix noise distribution
$\hat{\mathbf{A}}\sim D_{\omega^{*}}$ as well as the prior $P$ for
$\mathbf{b}$.

(12)
$\mathcal{E}_{\mathbf{B}}^{\text{Bayes}}(\tilde{\mathbf{x}}_{\beta})\equiv\mathbb{E}_{P}\mathbb{E}_{D_{\omega^{*}}}\left[\|\tilde{\mathbf{x}}_{\beta}-\mathbf{x}\|_{\mathbf{B}}^{2}\right]\,.$

Note that the formulation of this error is meant to include cases where $P$ is
supported at a single point (i.e., $\mathbf{b}$ is deterministic). Viewing the
problem in terms of operators, the error can equivalently be written in terms
of the operators $\hat{\mathbf{A}}$, $\mathbf{A}$ and right-hand-side
$\mathbf{b}$,

(13)
$\mathcal{E}_{\mathbf{B}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathbb{E}_{P}\mathbb{E}_{D_{\omega^{*}}}\left[\|\tilde{\mathbf{A}}_{\beta}^{-1}\mathbf{b}-\mathbf{A}^{-1}\mathbf{b}\|_{\mathbf{B}}^{2}\right]\,.$

Letting $\mathbf{R}$ be the auto-correlation matrix of $\mathbf{b}$, i.e.,
$\mathbf{R}=\mathbb{E}_{P}[\mathbf{b}\mathbf{b}^{T}]$, one can transform this
into a definition purely in terms of operators. Define the semi-norm

(14)
$\|\mathbf{A}\|_{\mathbf{B},\mathbf{R}}^{2}\equiv\langle\mathbf{A},\mathbf{A}\rangle_{\mathbf{B},\mathbf{R}}\equiv\operatorname*{tr}(\mathbf{R}^{1/2}\mathbf{A}^{T}\mathbf{B}\mathbf{A}\mathbf{R}^{1/2})\geq
0\,.$

Note that, in the isotropic-$L_{2}$ case where
$\mathbf{B}=\mathbf{R}=\mathbf{I}$, $\|\cdot\|_{\mathbf{B},\mathbf{R}}$ is the
well-known Frobenius norm. When the prior $P$ for $\mathbf{b}$ has non-
degenerate auto-correlation $\mathbf{R}$, the semi-norm
$\|\cdot\|_{\mathbf{B},\mathbf{R}}$ is a proper norm, as
$\mathbf{R}^{1/2}\mathbf{A}^{T}\mathbf{B}\mathbf{A}\mathbf{R}^{1/2}$ is a
symmetric positive definite matrix. Nonetheless, using the linearity and
cyclic properties of the trace on Eq. 13, one obtains that

(15)
$\mathcal{E}_{\mathbf{B}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathbb{E}_{D_{\omega^{*}}}\left[\|\tilde{\mathbf{A}}_{\beta}^{-1}-\mathbf{A}^{-1}\|_{\mathbf{B},\mathbf{R}}^{2}\right]\,.$

Therefore, the problem of finding a more accurate estimate
$\tilde{\mathbf{x}}_{\beta}$ of the solution $\mathbf{x}$ in the
$\|\cdot\|_{\mathbf{B}}$ norm can equivalently be seen as finding an estimator
$\tilde{\mathbf{A}}_{\beta}^{-1}$ so that the inverse
$\tilde{\mathbf{A}}_{\beta}^{-1}$ is closer to $\mathbf{A}^{-1}$ in the
$\|\cdot\|_{\mathbf{B},\mathbf{R}}$ semi-norm.

In the elliptic case, the choice of augmentation matrix $\hat{\mathbf{K}}$ is
usually a function of $\hat{\mathbf{A}}$ – as the intent is to move the
sampled matrix $\hat{\mathbf{A}}$ closer to some base-point matrix. For
example, the choice $\hat{\mathbf{K}}=\hat{\mathbf{A}}$ moves our estimated
operator $\tilde{\mathbf{A}}_{\beta}$ closer to the zero matrix $\mathbf{0}$.
Because of this, it is important to establish that the optimal $\beta^{*}$ is
greater than zero – if this is not the case, then our method is essentially
amplifying the noise in $\hat{\mathbf{A}}$ rather than diminishing it. This is
potentially dangerous when paired with the need to bootstrap the quantity
$\beta^{*}$.

A simple calculation gives that

(16)
$\beta^{*}=\frac{\mathbb{E}_{D_{\omega^{*}}}\langle\hat{\mathbf{K}},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{B},\mathbf{R}}}{\mathbb{E}_{D_{\omega^{*}}}\|\hat{\mathbf{K}}\|^{2}_{\mathbf{B},\mathbf{R}}}\,,$

hence the sign of $\beta^{*}$’s sign is equivalent to the sign of
$\mathbb{E}_{D_{\omega^{*}}}\langle\hat{\mathbf{K}},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{B},\mathbf{R}}$,
and we would therefore like an augmentation method to exhibit
$\mathbb{E}_{D_{\omega^{*}}}\langle\hat{\mathbf{K}},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{B},\mathbf{R}}\geq
0$. We will begin by studying the simplest choice of operator augmentation,

(17) $\hat{\mathbf{K}}=\hat{\mathbf{A}}^{-1}\,.$

The overshooting effect demonstrated in Fig. 1 gives one reason to believe
that this choice of augmentation is a reasonable one, as it involves shrinks
the inverse operator towards zero.

###  4.1 Conditions

Unlike in the symmetric case, to prove a rigorous statement about the sign of
$\beta^{*}$ for Eq. 17 in the nonsymmetric case, one must place a number of
additional conditions on the constituent components of the model. We will
discuss the necessity of these conditions in more detail in Section 5 – in
short, each of these conditions has an adversarial example associated with it
that causes the theory to fail when the conditions are not assumed.

In order to prove positivity of $\beta^{*}$, we introduce the following
constraints:

  1. 1.

Unbiased Noise: We assume that the noise matrix $\hat{\mathbf{Z}}$ is
unbiased, i.e., $\mathbb{E}[\hat{\mathbf{A}}]=\mathbf{A}$.

  2. 2.

Isotropy: We assume that
$\mathbf{R}=\mathbb{E}_{P}[\mathbf{b}\mathbf{b}^{T}]=\mathbf{I}$. This means
that there is no preferred direction in which we care about the accuracy of
the estimator $\tilde{\mathbf{A}}_{\beta}^{-1}$.

  3. 3.

Noise Symmetry: We assume that the distribution of the matrix
$\hat{\mathbf{A}}$ is symmetric about its mean, namely that $\hat{\mathbf{Z}}$
has the same distribution as $-\hat{\mathbf{Z}}$.

  4. 4.

Residual Norm: We specifically choose our norm of interest $\mathbf{B}$ to be
the residual norm $\mathbf{B}=\mathbf{A}^{T}\mathbf{A}$. The residual norm is
often used as an objective in nonsymmetric iterative methods.

We note that for the theory of elliptic operator augmentation, items (2)
through (4) are not necessary – and so their apparent necessity in the
nonsymmetric case is a curious mathematical phenomenon of the group
$GL(\mathbb{R}^{n})$.

###  4.2 Proof

With the conditions outlined above, we proceed to prove the positivity of the
augmentation factor $\beta^{*}$,

######  Theorem 4.1.

Let $\hat{\mathbf{A}}$ be a random matrix, invertible almost everywhere, such
that $\mathbb{E}[\hat{\mathbf{A}}^{-2}]$ exists. Under the conditions outlined
in Section 4.1, the optimal augmentation factor is always nonnegative.

######  Proof 4.2.

We must verify:

(18)
$\mathbb{E}\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}\geq
0\,.$

Expanding,

(19)
$\mathbb{E}\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}=\operatorname*{tr}\mathbb{E}\left[\hat{\mathbf{A}}^{-T}\mathbf{A}^{T}\mathbf{A}(\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1})\right]\,,$

Since $\mathbf{A}$ is invertible and $\hat{\mathbf{A}}$ is invertible almost
everywhere, let us define the matrix $\hat{\mathbf{Y}}$,

(20)
$\hat{\mathbf{Y}}=\hat{\mathbf{A}}\mathbf{A}^{-1}-\mathbf{I}=\hat{\mathbf{Z}}\mathbf{A}^{-1}\,.$

Note that this definition implies that
$\mathbb{E}[\hat{\mathbf{Y}}]=\mathbf{0}$ as well as that the distribution of
$\hat{\mathbf{Y}}$ is symmetric, since the distribution of $\hat{\mathbf{Z}}$
is symmetric by condition (3).

One can rearrange to obtain an expression for $\hat{\mathbf{A}}^{-1}$,

(21)
$\hat{\mathbf{A}}^{-1}=\mathbf{A}^{-1}(\mathbf{I}+\hat{\mathbf{Y}})^{-1}\,.$

And we substitute the above expression into Eq. 19,

(22)
$\begin{split}\mathbb{E}\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}&=\operatorname*{tr}\mathbb{E}\left[(\mathbf{I}+\hat{\mathbf{Y}})^{-T}(\mathbf{I}+\hat{\mathbf{Y}})^{-1}-(\mathbf{I}+\hat{\mathbf{Y}})^{-T}\right]\,.\\\
\end{split}$

Since the distribution of $\hat{\mathbf{Y}}$ is symmetric, it suffices to
verify that

(23)
$\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}+\mathbf{Y})^{-1}-(\mathbf{I}+\mathbf{Y})^{-T}\right]+\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y})^{-1}-(\mathbf{I}-\mathbf{Y})^{-T}\right]\geq
0\,,$

or alternatively,

(24)
$\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}+\mathbf{Y})^{-1}+(\mathbf{I}-\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y})^{-1}\right]\geq\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}+(\mathbf{I}-\mathbf{Y})^{-T}\right]\,,$

for all matrices $\mathbf{Y}$ for which $\mathbf{I}+\mathbf{Y}$ and
$\mathbf{I}-\mathbf{Y}$ are nonsingular.

In order to verify Eq. 24, we begin by considering the matrix,

(25)
$((\mathbf{I}-\mathbf{Y})^{-T}-(\mathbf{I}+\mathbf{Y})^{-1})((\mathbf{I}-\mathbf{Y})^{-1}-(\mathbf{I}+\mathbf{Y})^{-T})\succeq\mathbf{0}\,.$

Since this matrix is positive definite (by virtue of the fact that it has the
form $\mathbf{M}^{T}\mathbf{M}$), it follows that the trace of the above
matrix is positive,

(26)
$\operatorname*{tr}\left[((\mathbf{I}+\mathbf{Y})^{-T}-(\mathbf{I}-\mathbf{Y})^{-1})((\mathbf{I}+\mathbf{Y})^{-1}-(\mathbf{I}-\mathbf{Y})^{-T})\right]\geq
0\,.$

The above inequality can be rearranged,

(27)
$\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}+\mathbf{Y})^{-1}+(\mathbf{I}-\mathbf{Y})^{-1}(\mathbf{I}-\mathbf{Y})^{-T}\right]\geq\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y}^{2})^{-T}+(\mathbf{I}-\mathbf{Y}^{2})^{-1}\right]\,.$

Note that, by the cyclic property of the trace, that the left hand sides of
both Eq. 24 and Eq. 27 are identical. Therefore, it suffices to prove that

(28)
$\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y}^{2})^{-T}+(\mathbf{I}-\mathbf{Y}^{2})^{-1}\right]=\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}+(\mathbf{I}-\mathbf{Y})^{-T}\right]\,,$

from which Eq. 24 will follow.

To prove Eq. 28, we note that

(29)
$\begin{split}&\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}+(\mathbf{I}-\mathbf{Y})^{-T}\right]\\\
&\quad=\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y})^{T}(\mathbf{I}-\mathbf{Y})^{-T}+(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}+\mathbf{Y})^{T}(\mathbf{I}-\mathbf{Y})^{-T}\right]\\\
&\quad=\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y}+\mathbf{I}+\mathbf{Y})^{T}(\mathbf{I}-\mathbf{Y})^{-T}\right]\\\
&\quad=2\operatorname*{tr}\left[(\mathbf{I}+\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y})^{-T}\right]\\\
&\quad=2\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y}^{2})^{-T}\right]\\\
&\quad=\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y}^{2})^{-T}+(\mathbf{I}-\mathbf{Y}^{2})^{-1}\right]\,.\end{split}$

We remark that the critical step that requires isotropy is the assertion that
the left hand sides of both Eq. 24 and Eq. 27 are equal, namely that

(30)
$\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y})^{-1}(\mathbf{I}-\mathbf{Y})^{-T}\right]=\operatorname*{tr}\left[(\mathbf{I}-\mathbf{Y})^{-T}(\mathbf{I}-\mathbf{Y})^{-1}\right]\,.$

This is no longer necessarily true if we replace the isotropic trace operator
$\operatorname*{tr}(\cdot)$ with an anisotropic operator
$\operatorname*{tr}(\mathbf{R}^{1/2}(\cdot)\mathbf{R}^{1/2})$ for general
$\mathbf{R}$.

##  5 The Necessity of Conditions

In this section, we provide adversarial examples of how violating the
conditions in Section 4.1 can lead to situations where the optimal
augmentation factor $\beta^{*}$ for $\hat{\mathbf{K}}=\hat{\mathbf{A}}^{-1}$
is negative. This provides a concrete lens of how the group
$GL(\mathbb{R}^{n})$ differs from the symmetric positive definite cone
$S(\mathbb{R}^{n})_{+}$, as such situations are not possible in the SPD case.
These examples demonstrate how one might use the structure of
$GL(\mathbb{R}^{n})$ to construct situations where increasing the “noise” in
our estimator can actually lead to more accurate results.

###  5.1 The Necessity of Isotropy

First, we investigate the necessity of the isotropy condition. Suppose that
$\mathbf{A}=\mathbf{I}$ and $\hat{\mathbf{Z}}$ has two-atom distribution with
atoms:

(31) $\mathbf{Z}_{1}=\left[\begin{array}[]{cc}0&k\\\
0&-k\end{array}\right],\qquad\mathbf{Z}_{2}=\left[\begin{array}[]{cc}0&-k\\\
0&k\end{array}\right]\,,$

where each atom occurs with equal probability, and $k\gg 1$. Clearly in this
situation, the error distribution of $\hat{\mathbf{Z}}$ is symmetric and has
mean zero. Therefore, all of the other conditions in Section 4.1 are met. The
augmented operator $\tilde{\mathbf{A}}_{\beta}^{-1}$ is given by
$\tilde{\mathbf{A}}_{\beta}^{-1}=(1-\beta)\hat{\mathbf{A}}^{-1}$ and the
random matrix $\hat{\mathbf{A}}$ has two outcomes with equal probability,

(32) $\mathbf{A}_{1}=\left[\begin{array}[]{cc}1&k\\\
0&-k+1\end{array}\right],\qquad\mathbf{A}_{2}=\left[\begin{array}[]{cc}1&-k\\\
0&k+1\end{array}\right]\,.$

These outcomes have inverses

(33) $\mathbf{A}_{1}^{-1}=\left[\begin{array}[]{cc}1&k/(k+1)\\\
0&-1/(k-1)\end{array}\right],\qquad\mathbf{A}_{2}^{-1}=\left[\begin{array}[]{cc}1&k/(k+1)\\\
0&1/(k-1)\end{array}\right],.$

This means that in either case, we have

(34) $\hat{\mathbf{A}}^{-1}=\left[\begin{array}[]{cc}1&1\\\
0&0\end{array}\right]+O(1/k)\,.$

With this matrix ensemble, we can take the prior $P$ to be deterministic, such
that

(35) $\mathbf{b}=\left[\begin{array}[]{c}2\\\ -1\end{array}\right]\,.$

This immediately makes the problem with this setup evident, as

(36) $\hat{\mathbf{A}}^{-1}\mathbf{b}=\left[\begin{array}[]{c}1\\\
0\end{array}\right]+O(1/k),\qquad\mathbf{A}^{-1}\mathbf{b}=\left[\begin{array}[]{c}2\\\
-1\end{array}\right]\,.$

It is therefore clear that the minimum of the objective

(37)
$\mathcal{E}_{\mathbf{A}^{T}\mathbf{A}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathcal{E}_{\mathbf{I}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathbb{E}\|(1-\beta)\hat{\mathbf{A}}^{-1}\mathbf{b}-\mathbf{A}^{-1}\mathbf{b}\|_{2}^{2}$

must achieve its minimum at

(38) $\beta^{*}=-1+O(1/k)\,.$

One might initially conclude that the ability of this example to undermine
Theorem 4.1 has something to do with the assumption that $k$ is very large.
This is not the case. The largeness of $k$ is assumed only for illustrative
purposes. One can verify numerically (see Fig. 3) for the above choice of
$\hat{\mathbf{A}}$, that
$\mathbb{E}[\hat{\mathbf{A}}^{-T}\hat{\mathbf{A}}^{-1}-\hat{\mathbf{A}}^{-T}]$
has a negative eigenvalue for all $k\neq\pm 1,0$. Therefore, for any $k\neq\pm
1,0$, one can find a $\mathbf{b}$ such that $\beta^{*}$ is negative.

Figure 3: The lowest eigenvalue of
$\mathbb{E}[\hat{\mathbf{A}}^{-T}\hat{\mathbf{A}}^{-1}-\hat{\mathbf{A}}^{-T}]$
as defined in Section 5.1.

In particular, one should note that without a requirement of positive
definiteness, it is possible to create a situation where $\hat{\mathbf{A}}$ is
always “larger” than $\mathbf{A}$. This runs counter to the intuition behind
operator augmentation in the symmetric positive definite case [4], where such
a situation is not possible.

###  5.2 Outlier Masking and the Importance of Noise Symmetry

The second condition that one needed to prove the results in Section 4 is the
presence of symmetry in the noise distribution. In this section, we will see
how if this is not required, it is possible to construct an adversarial
example where $\beta^{*}<0$ for the augmentation Eq. 17. For this example, we
take $\mathbf{A}=\mathbf{I}$ and let $\hat{\mathbf{Z}}$ have a distribution of
three equally probable atoms, given by

(39)
$\mathbf{Z}_{1}=\mathbf{I},\qquad\mathbf{Z}_{2}=k\mathbf{I},\qquad\mathbf{Z}_{3}=-(k+1)\mathbf{I}\,,$

where $k\gg 1$. Note that this distribution is mean zero. Computing the atoms
of the distribution $\hat{\mathbf{A}}^{-1}$,

(40)
$\mathbf{A}^{-1}_{1}=\frac{1}{2}\mathbf{I},\qquad\mathbf{A}^{-1}_{2}=\frac{1}{k+1}\mathbf{I},\qquad\mathbf{A}^{-1}_{3}=-\frac{1}{k}\mathbf{I}\,.$

Now, note that in order to minimize the quantity

(41)
$\mathcal{E}_{\mathbf{A}^{T}\mathbf{A}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathbb{E}\|\tilde{\mathbf{A}}_{\beta}^{-1}-\mathbf{I}\|_{F}^{2}=\frac{2}{3}\left(\frac{1-\beta}{2}-1\right)^{2}+\frac{2}{3}\left(\frac{1-\beta}{k+1}-1\right)^{2}+\frac{2}{3}\left(-\frac{1-\beta}{k}-1\right)^{2}\,,$

one can verify by taking the derivative and setting it to zero that

(42) $\beta^{*}=-1+O(1/k)\,.$

Therefore, the optimal augmentation will grow the inverse operator instead of
shrinking it.

The message of this example is that _outliers in the matrix noise can mask
distribution imbalances in the region near $\mathbf{A}$_ that can cause
$\mathbb{E}[\hat{\mathbf{A}}^{-1}]$ to both lie in the direction of
$\mathbf{A}^{-1}$ while at the same time being dominated by $\mathbf{A}^{-1}$.
Indeed, we have that
$\mathbb{E}[\hat{\mathbf{A}}^{-1}]\approx\frac{1}{2}\mathbf{I}\preceq\mathbf{I}=\mathbf{A}^{-1}$
(it bears repeating that such a feat is impossible in the symmetric positive
definite setting where
$\mathbb{E}[\hat{\mathbf{A}}^{-1}]\succeq\mathbf{A}^{-1}$ [4]). The importance
of noise symmetry is that it forces the distribution of $\hat{\mathbf{Z}}$ to
be balanced in the region around $\mathbf{A}$, even if the distribution
contains large outliers.

###  5.3 The Importance of Conditioning and a Counter-Example for the L2 Norm

Our final counter-example concerns the use of the $L_{2}$ norm in the
objective Eq. 15 rather than the residual norm. In the symmetric positive
definite case, one can prove the positivity of the optimal augmentation factor
for a large range of different objective norms [4]. However, as we will see in
this section, there are adversarial ways to break norms other than the
residual norm in the nonsymmetric case. We will focus here on giving an
example that shows how if using the $L_{2}$ norm instead of the residual norm
makes it possible to have a negative optimal augmentation factor for the
augmentation Eq. 17.

To begin, consider the ground truth matrix

(43) $\mathbf{A}=\left[\begin{array}[]{cc}1&-\epsilon\\\
\epsilon&0\end{array}\right]\,,$

where $0<\epsilon\ll 1$. For the noise $\hat{\mathbf{Z}}$, we reuse the noise
distribution from Section 5.1 — consider a two-atom distribution with atoms

(44) $\mathbf{Z}_{1}=\left[\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right],\qquad\mathbf{Z}_{2}=\left[\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right]\,,$

where each atom has equal probability. Note that this distribution is
symmetric and mean zero. Since $\epsilon$ is extremely small, this means that
the distribution of $\hat{\mathbf{A}}$ will have two atoms whose inverses are
approximately

(45) $\mathbf{A}^{-1}_{1}\approx\left[\begin{array}[]{cc}1&1\\\
-1&0\end{array}\right]^{-1}=\left[\begin{array}[]{cc}0&1\\\
-1&1\end{array}\right],\qquad\mathbf{A}^{-1}_{2}\approx\left[\begin{array}[]{cc}1&-1\\\
1&0\end{array}\right]^{-1}=\left[\begin{array}[]{cc}0&-1\\\
1&1\end{array}\right]\,.$

Note that the inverses of these atoms have the same Frobenius norm the atoms
themselves. In sharp contrast, the ill-conditioning of $\mathbf{A}$ means that
$\mathbf{A}^{-1}$ is an order of magnitude larger than $\mathbf{A}$,

(46) $\mathbf{A}^{-1}=\left[\begin{array}[]{cc}0&-\epsilon^{-1}\\\
\epsilon^{-1}&\epsilon^{-2}\end{array}\right]\,.$

Immediately, we see that in order for $\beta^{*}$ to minimize

(47)
$\mathcal{E}_{\mathbf{I}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\frac{1}{2}\|(1-\beta)\mathbf{A}^{-1}_{1}-\mathbf{A}^{-1}\|_{F}^{2}+\frac{1}{2}\|(1-\beta)\mathbf{A}^{-1}_{2}-\mathbf{A}^{-1}\|_{F}^{2}\,,$

One must therefore have $\beta^{*}\sim-\epsilon^{-2}$, for which
$\mathcal{E}_{\mathbf{I}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta^{*}})\sim\epsilon^{-1}$.
For other growth orders of $\beta$, one has
$\mathcal{E}_{\mathbf{I}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\omega(\epsilon^{-1})$.
It is therefore clear that for $\epsilon$ small enough, $\beta^{*}$ will be
negative.

In contrast to the $L_{2}$ norm, note that the residual norm for this problem
is given by

(48)
$\mathbf{A}^{T}\mathbf{A}=\left[\begin{array}[]{cc}1+\epsilon^{2}&\epsilon\\\
\epsilon&0\end{array}\right]\,,$

Therefore, we see that the reason why Theorem 4.1 holds in the residual norm
but not in the $L_{2}$ norm is because the residual norm places significantly
less weight on the part of the matrix $\mathbf{A}^{-1}$ that contributes to
the $\mathbf{A}^{-1}$’s large increase of magnitude over $\mathbf{A}$. This
means that the residual norm $\mathbf{A}^{T}\mathbf{A}$ accounts for ill-
conditioning on $\mathbf{A}$ in a way that the $L_{2}$ norm does not.

This example therefore also demonstrates the importance of conditioning in the
ground truth matrix $\mathbf{A}$. It is perfectly possible that $\mathbf{A}$
lies close to a singular matrix, while the outcomes of $\hat{\mathbf{A}}$ are
moved away from singularity by the noise imparted by $\hat{\mathbf{Z}}$. If
this is the case, $\mathbb{E}[\hat{\mathbf{A}}^{-1}]$ will be small in
magnitude compared to $\mathbf{A}^{-1}$ and shrinking the operator
$\hat{\mathbf{A}}^{-1}$ further will not reduce average error in the Frobenius
sense.

Note that symmetric positive definite setting [4] avoids this issue, since if
$\hat{\mathbf{A}}$ is symmetric positive definite everywhere, it is impossible
for $\mathbb{E}[\hat{\mathbf{A}}]$ to be close to the origin without a
significant chunk of the probability distribution also lying close to the
origin. This ensures that $\mathbb{E}[\hat{\mathbf{A}}^{-1}]$ will always
spectrally dominate $\hat{\mathbf{A}}^{-1}$, and shrinking the operator
$\hat{\mathbf{A}}^{-1}$ will reduce error.

##  6 The Small Noise Regime Allows for Nonsymmetric Noise

As we saw in Section 5.2, the presence of large adversarial outliers can
completely mask local imbalances in the noise distribution near $\mathbf{A}$.
These local imbalances can be severe enough to invalidate Theorem 4.1.
Naturally, the example presented in Section 5.2 is quite extreme, so one might
ask if the issue inherent is not necessarily the _symmetry_ of the noise, but
rather the _magnitude_ of the noise. To answer this question, we consider the
_small noise regime_ , where deviations in $\hat{\mathbf{Y}}$ are very small
relative to $\mathbb{E}[\hat{\mathbf{Y}}]=\mathbf{I}$. It turns out, that if
one assumes that terms of order $O(\mathbb{E}[\hat{\mathbf{Y}}^{3}])$ are
negligible — what we term the _small noise regime_ — then the symmetry
assumption is unnecessary, as we will see momentarily.

For this discussion, we duplicate the setting of Theorem 4.1, except we
replace the condition that the noise distribution is symmetric with the
condition that noise terms of the order $O(\mathbb{E}[\hat{\mathbf{Y}}^{3}])$
are negligible. Recall that the statement necessary for Theorem 4.1 to be true
was

(49)
$\begin{split}\mathbb{E}\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}&=\operatorname*{tr}\mathbb{E}\left[(\mathbf{I}+\hat{\mathbf{Y}})^{-T}(\mathbf{I}+\hat{\mathbf{Y}})^{-1}-(\mathbf{I}+\hat{\mathbf{Y}})^{-T}\right]\geq
0\,.\\\ \end{split}$

We define the function $f$,

(50)
$f(\hat{\mathbf{Y}})\equiv\operatorname*{tr}\left[(\mathbf{I}+\hat{\mathbf{Y}})^{-T}(\mathbf{I}+\hat{\mathbf{Y}})^{-1}-(\mathbf{I}+\hat{\mathbf{Y}})^{-T}\right]\,.$

Taking a second order Taylor expansion of $f$ about
$\hat{\mathbf{Y}}=\mathbf{0}$, we obtain

(51) $f(\hat{\mathbf{Y}})=f(\mathbf{0})+\delta
f(\hat{\mathbf{Y}})+\delta^{2}f(\hat{\mathbf{Y}})+O(\hat{\mathbf{Y}}^{3})\,.$

Note that $f(\mathbf{0})=0$ and that $\delta f(\hat{\mathbf{Y}})$ is linear in
$\hat{\mathbf{Y}}$. Hence, in expectation, both of these terms vanish,

(52)
$\mathbb{E}[f(\hat{\mathbf{Y}})]=\mathbb{E}[\delta^{2}f(\hat{\mathbf{Y}})]+O(\mathbb{E}[\hat{\mathbf{Y}}^{3}])\,.$

A calculation of $\delta^{2}f$ gives

(53)
$\begin{split}\mathbb{E}[\delta^{2}f(\hat{\mathbf{Y}})]&=\mathbb{E}\operatorname*{tr}\left[2\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}^{T}+2\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}+2\hat{\mathbf{Y}}\hat{\mathbf{Y}}-2\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}^{T}\right]\\\
&=2\,\mathbb{E}\operatorname*{tr}\left[\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}+\hat{\mathbf{Y}}\hat{\mathbf{Y}}\right]\\\
&=\mathbb{E}\operatorname*{tr}\left[\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}^{T}+\hat{\mathbf{Y}}^{T}\hat{\mathbf{Y}}+\hat{\mathbf{Y}}\hat{\mathbf{Y}}^{T}+\hat{\mathbf{Y}}\hat{\mathbf{Y}}\right]\\\
&=\mathbb{E}\operatorname*{tr}\left[(\hat{\mathbf{Y}}+\hat{\mathbf{Y}}^{T})^{T}(\hat{\mathbf{Y}}+\hat{\mathbf{Y}}^{T})\right]\geq
0\,.\end{split}$

Indeed, if $\hat{\mathbf{Y}}$ not zero almost everywhere, then the above is in
fact a strict inequality, meaning that Eq. 49 must be true to second order.

##  7 More General Augmentations and Anisotropic Priors

In the elliptic case [4], one can prove a reduction in error for a variety of
different augmentations. In particular, the energy norm in the elliptic
augmentation setting has an extensive theory regarding the approximation of
$\beta^{*}$ via truncated Taylor expansions that is valid for augmentations of
the form

(54) $\hat{\mathbf{K}}\equiv\hat{\mathbf{A}}^{-1}\mathbf{C}\,,$

where $\mathbf{C}$ is a symmetric positive definite matrix that commutes with
the auto-correlation of $\mathbf{b}$, i.e., $[\mathbf{C},\mathbf{R}]=0$. Part
of this theory states that $\beta^{*}$ must always be positive.

In the analysis in [4], adding this $\mathbf{C}$ term effectively just changes
the autocorrelation matrix from $\mathbf{R}$ to $\mathbf{R}\mathbf{C}$ in the
analysis. Note that $\mathbf{R}\mathbf{C}$ remains symmetric positive definite
by virtue of the condition that $[\mathbf{C},\mathbf{R}]=0$. This therefore
begs the question: do augmentations of the form Eq. 54 retain their nice
properties in the nonsymmetric augmentation setting — where the residual norm
plays the role of the energy norm in the elliptic augmentation setting?
Section 5.1 provides a definitive answer to this question: it is not possible
unless the modified auto-correlation is isotropic – that is, unless
$\mathbf{R}\mathbf{C}=\mathbf{I}$.

Nonetheless, while the wide number of choices regarding norms and
augmentations do not translate to the nonsymmetric augmentation setting, the
operator augmentation framework does provide a way for handling anisotropic
right-hand-side $\mathbf{b}$ with
$\mathbb{E}[\mathbf{b}\mathbf{b}^{T}]=\mathbf{R}\neq\mathbf{I}$. Namely, one
chooses the operator augmentation

(55) $\hat{\mathbf{K}}\equiv\hat{\mathbf{A}}^{-1}\mathbf{R}^{-1}\,.$

It is immediate that the results of Theorem 4.1 hold for this choice of
augmentation — as the $\mathbf{R}^{-1}$ will cancel the $\mathbf{R}$ in the
error objective. We restate this conclusion into the main theorem of this
paper:

######  Theorem 7.1.

Let $\hat{\mathbf{A}}\sim D_{\omega^{*}}$ be a random matrix, invertible
almost everywhere, such that $\mathbb{E}[\hat{\mathbf{A}}^{-2}]$ exists.
Suppose that the distribution of $\hat{\mathbf{A}}$ is symmetric about its
mean, and that $\mathbb{E}[\hat{\mathbf{A}}]=\mathbf{A}$. Let $\mathbf{b}\sim
P$ with nonsingular auto-correlation matrix
$\mathbb{E}[\mathbf{b}\mathbf{b}^{T}]=\mathbf{R}$. Under these assumptions,
consider the semi-Bayesian residual error:

(56)
$\mathcal{E}_{\mathbf{A}^{T}\mathbf{A}}^{\text{Bayes}}(\tilde{\mathbf{A}}_{\beta})=\mathbb{E}_{P}\mathbb{E}_{D_{\omega^{*}}}\left[\|\tilde{\mathbf{A}}_{\beta}^{-1}\mathbf{b}-\mathbf{A}^{-1}\mathbf{b}\|^{2}_{\mathbf{A}^{T}\mathbf{A}}\right]=\mathbb{E}_{D_{\omega^{*}}}\left[\|\tilde{\mathbf{A}}_{\beta}^{-1}-\mathbf{A}^{-1}\|^{2}_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}\right]\,.$

Consider the operator augmentation:

(57)
$\tilde{\mathbf{A}}_{\beta}^{-1}=\hat{\mathbf{A}}^{-1}-\beta\hat{\mathbf{A}}^{-1}\mathbf{R}^{-1}\,.$

Then the optimal $\beta^{*}$ with respect to the semi-Bayesian error is always
nonnegative.

######  Proof 7.2.

As in Theorem 4.1, we must verify

(58)
$\mathbb{E}\langle\hat{\mathbf{A}}^{-1}\mathbf{R}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}\geq
0\,.$

But note, by definition in Eq. 14 and symmetry of $\mathbf{R}$, that

(59)
$\langle\hat{\mathbf{A}}^{-1}\mathbf{R}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}=\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}$

The condition Eq. 58 is therefore equivalent to

(60)
$\mathbb{E}\langle\hat{\mathbf{A}}^{-1},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{I}}\geq
0\,,$

which we proved in Theorem 4.1.

##  8 Monte Carlo Estimation

To build an algorithm from this theory, we attempt to estimate the optimal
augmentation factor with bootstrap Monte Carlo. The quantity of interest is

(61)
$\beta^{*}=\frac{\mathbb{E}_{D_{\omega^{*}}}\langle\hat{\mathbf{K}},\hat{\mathbf{A}}^{-1}-\mathbf{A}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}}{\mathbb{E}_{D_{\omega^{*}}}\|\hat{\mathbf{K}}\|^{2}_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}}\,.$

Obviously, it is not possible to compute this quantity directly, since we only
have a single sample $\mathbf{A}$ and we do not have access to the ground
truth $\mathbf{A}$ or its distribution $D_{\omega^{*}}$. We must therefore try
to approximate these quantities as best as possible with the available
information. Suppose that $\hat{\omega}\in\Omega$ is the parameter instance
that generates the matrix $\hat{\mathbf{A}}=\mathcal{M}(\hat{\omega})$. To
bootstrap Eq. 61, we replace all instances of $\mathbf{A}$ with
$\hat{\mathbf{A}}$ and draw random instances from $D_{\hat{\omega}}$ instead
of $D_{\omega^{*}}$. We get

(62)
$\tilde{\beta}=\frac{\mathbb{E}_{D_{\hat{\omega}}}\langle\hat{\mathbf{K}}_{b},\hat{\mathbf{A}}^{-1}_{b}-\hat{\mathbf{A}}^{-1}\rangle_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}}{\mathbb{E}_{D_{\hat{\omega}}}\|\hat{\mathbf{K}}_{b}\|^{2}_{\mathbf{A}^{T}\mathbf{A},\mathbf{R}}}\,,$

where $\hat{\mathbf{K}}_{b}$ and $\hat{\mathbf{A}}_{b}$ are drawn from
$D_{\hat{\omega}}$. This can then be discretized with Monte Carlo:

(63)
$\hat{\beta}=\frac{\sum_{i=1}^{m}\mathbf{b}_{i}^{T}\hat{\mathbf{K}}_{b,i}^{T}\hat{\mathbf{A}}^{T}\hat{\mathbf{A}}(\hat{\mathbf{A}}^{-1}_{b,i}-\hat{\mathbf{A}}^{-1})\mathbf{b}_{i}}{\sum_{i=1}^{m}\mathbf{b}_{i}^{T}\hat{\mathbf{K}}_{b,i}^{T}\hat{\mathbf{A}}^{T}\hat{\mathbf{A}}\hat{\mathbf{K}}_{b,i}\mathbf{b}_{i}}\,,\qquad\mathbf{b}_{i}\sim
P\;\text{i.i.d.},\quad\hat{\mathbf{A}}_{b,i}\sim
D_{\hat{\omega}}\;\text{i.i.d.}\,.$

If one opts to use the augmentation provided by Eq. 55, the expression
simplifies to

(64) $\begin{aligned}
\hat{\beta}&=1-\frac{\sum_{i=1}^{m}\mathbf{b}_{i}^{T}\hat{\mathbf{A}}_{b,i}^{-T}\hat{\mathbf{A}}^{T}\mathbf{b}_{i}}{\sum_{i=1}^{m}\mathbf{q}_{i}^{T}\hat{\mathbf{A}}_{b,i}^{-T}\hat{\mathbf{A}}^{T}\hat{\mathbf{A}}\hat{\mathbf{A}}_{b,i}^{-1}\mathbf{q}_{i}}\,,\end{aligned}\qquad\begin{aligned}
\mathbf{b}_{i}&\sim\mathcal{N}(\mathbf{0},\mathbf{I})\;\text{i.i.d.},\\\
\mathbf{q}_{i}&\sim\mathcal{N}(\mathbf{0},\mathbf{R}^{-1})\;\text{i.i.d.},\\\
\hat{\mathbf{A}}_{b,i}&\sim D_{\hat{\omega}}\;\text{i.i.d.}\,.\end{aligned}$

Note that $\mathbf{b}_{i}$ and $\mathbf{q}_{i}$ do not necessarily have to be
normal – they must just have the same auto-correlation as the normal
distributions specified. For good measure, one might also threshold the above
expression to guarantee $\hat{\beta}\geq 0$. Finally, one can now estimate the
true solution to the system Eq. 1 via

(65)
$\tilde{\mathbf{x}}_{\hat{\beta}}=\tilde{\mathbf{A}}_{\hat{\beta}}^{-1}\mathbf{b}=(\hat{\mathbf{A}}^{-1}-\hat{\beta}\hat{\mathbf{A}}^{-1}\mathbf{R}^{-1})\mathbf{b}=\hat{\mathbf{A}}^{-1}(\mathbf{b}-\hat{\beta}\mathbf{R}^{-1}\mathbf{b})\,.$

###  8.1 Approximation via Taylor Expansion

Note that every Monte Carlo sample in Eq. 64 requires inverting a matrix
system. There are times where this may be too computationally expensive to be
feasible. However, if we are in the _small noise regime_ , one may take a
Taylor expansion of $\hat{\mathbf{A}}^{-1}$ about $\mathbf{A}^{-1}$ — as this
means one only has the factorize an operator once for the whole estimation
process. The expansion to second order is given by

(66)
$\hat{\mathbf{A}}^{-1}=\mathbf{A}^{-1}-\mathbf{A}^{-1}\hat{\mathbf{Z}}\mathbf{A}^{-1}+2\mathbf{A}^{-1}\hat{\mathbf{Z}}\mathbf{A}^{-1}\hat{\mathbf{Z}}\mathbf{A}^{-1}+O(\hat{\mathbf{Z}}^{3})\,.$

Inserting this expression into Eq. 64 yields

(67) $\hat{\beta}\approx
1-\frac{\operatorname*{tr}(\mathbf{I})+\frac{2}{m}\sum_{i=1}^{m}\mathbf{b}_{i}^{T}(\hat{\mathbf{Z}}_{b,i}\hat{\mathbf{A}}^{-1})^{2}\mathbf{b}_{i}}{\operatorname*{tr}(\mathbf{R}^{-1})+\frac{1}{m}\sum_{i=1}^{m}\|\hat{\mathbf{Z}}_{b,i}\hat{\mathbf{A}}^{-1}\mathbf{q}_{i}\|_{2}^{2}+\frac{4}{m}\sum_{i=1}^{m}\mathbf{q}_{i}^{T}(\hat{\mathbf{Z}}_{b,i}\hat{\mathbf{A}}^{-1})^{2}\mathbf{q}_{i}}\,,$

note that we have omitted linear terms, because terms linear in
$\hat{\mathbf{Z}}$ will be zero in expectation, by virtue of the fact that
$\mathbb{E}[\hat{\mathbf{Z}}]=0$. It is possible that higher orders of
truncation may produce better results; however, unlike the elliptic
augmentation setting [4], it is difficult to prove guarantees about the
quality of these truncated expansions.

##  9 Numerical Experiments

To confirm our theoretical results, we will examine the asymmetric operator
augmentation algorithm applied to the problem of more accurately computing a
value function for a Markov chain whose probability transition function must
be estimated from data.

###  9.1 Background

For context, a discrete Markov chain is a stochastic process $X_{i}$ for
$i\in\mathbb{Z}_{\geq 0}$ on a finite state set $V$ is a process that
satisfies the Markov property

(68) $\mathbb{P}(X_{i}=v_{i}\mid
X_{i-1}=v_{i-1},...,X_{1}=v_{1})=\mathbb{P}(X_{i}=v_{i}\mid
X_{i-1}=v_{i-1})\,.$

A discrete Markov chain is time homogeneous if the right-hand-side of Eq. 68
does not depend on the time $i$. In this situation, the Markov chain is
completely characterized by its probability transition matrix

(69) $\mathbf{P}_{v,u}=\mathbb{P}(X_{i}=u\mid X_{i-1}=v)\,,$

as well as a distribution $\mathbb{P}(X_{0}=v)$ of the initial state $X_{0}$.
There are countless examples of such Markov chains from the fields of
probability, statistics, reinforcement learning, and physics.

One is often interested in computing or approximating functionals of the
process $X_{i}$. For example, one may think of $X_{i}$ as an agent navigating
through the set $V$ via the means of some fixed policy, where the agent
obtains a fixed reward $r(v)$ whenever it transitions from the state $v$. In
many situations in reinforcement learning, one is often interested in the
average discounted reward $Q(v)$ the agent obtains over its life-cycle when
beginning at $X_{0}=v$, namely,

(70) $Q(v)=\mathbb{E}\left[\sum_{i=0}^{\infty}\gamma^{i}r(X_{i})\mid
X_{0}=v\right]\,.$

$Q(v)$ is typically referred to as the _value function_ of $X_{i}$ with
respect to reward function $r(v)$. It is used to gauge the quality of the
agent’s policy at maximizing the discounted reward it receives over its life-
time. One notes that the function $Q$ is linear in $r$. One may use first
transition analysis to express the relationship between $Q$ and $r$ in matrix
form. Isolating the first term in the infinite sum Eq. 70 gives

(71) $Q(v)=r(v)+\gamma\sum_{u}\mathbb{P}(X_{1}=u\mid X_{0}=v)Q(u)\,.$

Alternatively,

(72) $\mathbf{q}=\mathbf{r}+\gamma\mathbf{P}\mathbf{q}\,,$

where $\mathbf{q}$ is the vector with entries $\mathbf{q}_{v}=Q(v)$ and
$\mathbf{r}$ is the vector with entries $\mathbf{r}_{v}=r(v)$. This finally
produces a linear system for the value function,

(73)
$\mathbf{A}\mathbf{q}\equiv(\mathbf{I}-\gamma\mathbf{P})\mathbf{q}=\mathbf{r}\,.$

However, there are many situations where the transition matrix $\mathbf{P}$
may not be known exactly. In reinforcement learning, for example, one does not
have access to $\mathbf{P}$ itself, but rather a number of finite realizations
of the process $X_{i}$ as it traverses the state space. Therefore, instead of
having access to the ground truth $\mathbf{P}$, one usually has access to a
noisy version $\hat{\mathbf{P}}$ thereof. A naive solve will give

(74) $(\mathbf{I}-\gamma\hat{\mathbf{P}})\hat{\mathbf{q}}=\mathbf{r}\,.$

We will see how operator augmentation can be used to reduce the average error
between our estimate $\hat{\mathbf{q}}$ and the ground truth $\mathbf{q}$ in
the residual norm.

###  9.2 Noise Model

One popular way to approximate a Markov chain’s probability transition matrix
from a sample $\hat{X}_{0},\hat{X}_{1},...,\hat{X}_{t}$ of the Markov chain is
to approximate the probability $\mathbb{P}(X_{i}=v\mid X_{i-1}=u)$ by
examining the fraction of times $X_{i}$ transitions to $v$ from $u$ out of the
total number of times $X_{i}$ transitions from $u$, i.e.,

(75) $\mathbb{P}(X_{i}=v\mid X_{i-1}=u)=\frac{\mathbb{P}(X_{i}=v\wedge
X_{i-1}=u)}{\mathbb{P}(X_{i-1}=u)}\approx\frac{\\#\\{i\mid\hat{X}_{i}=v\wedge\hat{X}_{i-1}=u\\}}{\\#\\{i\mid\hat{X}_{i}=u\wedge
i\geq 1\\}}\,.$

But while this estimator is commonly used — it is difficult to ensure coverage
of the state space in a way that is natural, as one must take care to ensure
that $\\#\\{i\mid\hat{X}_{i}=u\wedge i\geq 1\\}$ is positive for every $u$.

Therefore, to model the uncertainty in a Markov chain constructed from data,
we assume that, instead of observing the first $t$ states of the Markov chain
$X_{i}$, we instead observe the first $N$ transitions of the Markov chain
$X_{i}$ out of every state $v$. We denote the first $N$ transitions out of a
state $v$ as $\hat{Y}_{v,1},\hat{Y}_{v,2},...,\hat{Y}_{v,N}$. If the Markov
chain is irreducible, $\hat{Y}_{v,1},\hat{Y}_{v,2},...,\hat{Y}_{v,N}$ are
well-defined almost surely, and have distribution

(76) $\hat{Y}_{v,1},\hat{Y}_{v,2},...,\hat{Y}_{v,N}\sim\mathbb{P}(X_{i}\mid
X_{i-1}=v)\,,\qquad\text{i.i.d.}$

One may then easily estimate the probability transition matrix by using the
empirical distribution of $\mathbb{P}(X_{i}\mid X_{i-1}=v)$ to construct each
row of $\hat{\mathbf{P}}$,

(77)
$\hat{\mathbf{P}}_{v,u}\equiv\frac{\\#\\{\hat{Y}_{v,i}=u\\}}{N}\approx\mathbb{P}(X_{i}=u\mid
X_{i-1}=v)=\mathbf{P}_{v,u}\,.$

Note that $\hat{\mathbf{P}}$ is unbiased and hence that

(78) $\mathbb{E}[\hat{\mathbf{A}}]=\mathbf{A}\,.$

We pose the question: is it possible to use operator augmentation to increase
the accuracy of the naive estimate of the value function $Q$? Our numerical
results suggest that the answer to this question is yes.

###  9.3 Numerical Results

To test the theory, we consider 3 different Markov chains on 1D, 2D, and 3D
grids respectively. For our experiments, we use the operator augmentation
$\hat{\mathbf{K}}=\hat{\mathbf{A}}$,

(79)
$\tilde{\mathbf{A}}^{-1}=\hat{\mathbf{A}}^{-1}-\hat{\beta}\hat{\mathbf{K}}=(1-\hat{\beta})\hat{\mathbf{A}}^{-1}=(1-\hat{\beta})(\mathbf{I}-\gamma\hat{\mathbf{P}})^{-1}\,,$

where we use the Taylor expanded Monte Carlo approximation $\hat{\beta}$ from
Eq. 67. We compute our augmented solution via

(80) $\tilde{\mathbf{q}}=\tilde{\mathbf{A}}^{-1}\mathbf{r}\,.$

####  9.3.1 Random Walk with Drift on a 1D Grid

For this example, we consider the Markov chain of a lazy random walk with
drift on a 1D periodic grid graph. We let the state space $V$ be the set of
integers $\\{1,2,...,K\\}$ where $K=16$ and take the probability transition
matrix to be

(81) $\mathbb{P}^{(\text{1D})}_{\ell,r}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}\ell&u=v-1\mod K\\\ 1-\ell-r&u=v\mod K\\\ r&u=v+1\mod
K\\\ 0&\text{o.w.}\end{cases}$

We test a handful of different values for $(\ell,r)$ as well as a handful of
different values for the sample count $N$ that determines the “noise” in the
matrix $\hat{\mathbf{P}}$ (higher $N$ means less noise). For the reward
function — we consider two cases: for the first, we consider deterministic
reward function given by

(82) $r^{(\text{1D})}(v)=\sin(4\pi v/K)\,.$

For the second, we consider a random isotropic reward vector prior where
$\mathbf{r}^{(\text{1D})}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$.

In addition to the above collection of transition matrices, we also consider
random walks with the ability to skip a vertex in the grid,

(83) $\mathbb{P}^{(\text{1D})}_{\ell_{1},\ell_{2},r_{1},r_{2}}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}\ell_{1}&u=v-2\mod K\\\ \ell_{2}&u=v-1\mod K\\\
1-\ell_{1}-\ell_{2}-r_{1}-r_{2}&u=v\mod K\\\ r_{1}&u=v+1\mod K\\\
r_{2}&u=v+2\mod K\\\ 0&\text{o.w.}\end{cases}$

as well as transition matrices corresponding to a random walk on a complete
graph,

(84) $\mathbb{P}^{(\text{1D})}_{\text{complete}}(X_{i}=u\mid
X_{i-1}=v)=\frac{1}{K}\,.$

For these 1D experiments, we use discount factor

(85) $\gamma^{(1D)}=0.99\,.$

We present the results of our numerical experiments for the specified reward
vectors in Table 1 as well as for isotropic reward vectors in Table 2.

####  9.3.2 Random walk with Drift on 2D and 3D Grids

Now let us consider the Markov chain of a lazy random walk with non-uniform
drift on a 2D periodic grid graph. We let the state space $V$ be the set of
tuples $\\{1,2,...,K\\}\times\\{1,2,...,K\\}$ where $K=16$ and we consider two
probability transition matrices, the first a standard uniform random walk:

(86) $\mathbb{P}^{(\text{2D})}_{\text{unif}}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}1/4&|u-v|=1\\\ 0&\text{o.w.}\end{cases}$

and the second a random walk with a non-uniform drift:

(87) $\mathbb{P}^{(\text{2D})}_{\text{nonunif}}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}\frac{1}{4}\pm\frac{1}{8}\sin(2\pi
v_{x}/K)&u=v\mp(1,0)\\\ \frac{1}{4}\pm\frac{1}{8}\sin(2\pi
v_{y}/K)&u=v\mp(0,1)\\\ 0&\text{o.w.}\end{cases}$

For the reward function, we consider both a deterministic reward

(88) $r^{(\text{2D})}(v)=-\sin(2\pi v_{x}/K)\sin(2\pi v_{y}/K)\,,$

as well as an isotropic reward vector prior where
$\mathbf{r}^{(\text{2D})}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. For these 2D
experiments, we use discount factor

(89) $\gamma^{(2D)}=0.99\,.$

We can define similar transition matrices on 3D periodic grid graphs, where
$V=\\{1,...,K\\}^{3}$ and $K=8$:

(90) $\mathbb{P}^{(\text{3D})}_{\text{unif}}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}1/6&|u-v|=1\\\ 0&\text{o.w.}\end{cases}$

as well as an analogous 3D random walk with a non-uniform drift:

(91) $\mathbb{P}^{(\text{3D})}_{\text{nonunif}}(X_{i}=u\mid
X_{i-1}=v)=\begin{cases}\frac{1}{6}\pm\frac{1}{8}\sin(2\pi
v_{x}/K)&u=v\mp(1,0,0)\\\ \frac{1}{6}\pm\frac{1}{8}\sin(2\pi
v_{y}/K)&u=v\mp(0,1,0)\\\ \frac{1}{6}\pm\frac{1}{8}\sin(2\pi
v_{z}/K)&u=v\mp(0,0,1)\\\ 0&\text{o.w.}\end{cases}$

For the reward function, we consider both a deterministic reward

(92) $r^{(\text{3D})}(v)=-\sin(2\pi v_{x}/K)\sin(2\pi v_{y}/K)\sin(2\pi
v_{z}/K)\,,$

as well as an isotropic reward vector prior where
$\mathbf{r}^{(\text{3D})}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. For these 3D
experiments, we use discount factor

(93) $\gamma^{(3D)}=0.9\,.$

We present the results of our numerical experiments for the specified reward
vectors in Table 1 as well as for isotropic reward vectors in Table 2.

###  9.4 Discussion

Chain | Samples ($N$) | Naive Error | $\pm 2\sigma$ | Augmented Error | $\pm 2\sigma$
---|---|---|---|---|---
$\mathbb{P}^{(\text{1D})}_{1/4,1/4}$ | 16 | $166\%$ | $\pm 13.5\%$ | $53.5\%$ | $\pm 1.35\%$
| 32 | $48.8\%$ | $\pm 1.50\%$ | $30.1\%$ | $\pm 0.53\%$
| 64 | $20.1\%$ | $\pm 0.45\%$ | $16.0\%$ | $\pm 0.27\%$
$\mathbb{P}^{(\text{1D})}_{1/6,2/6}$ | 16 | $115\%$ | $\pm 10.2\%$ | $42.2\%$ | $\pm 1.57\%$
| 32 | $31.7\%$ | $\pm 1.11\%$ | $20.8\%$ | $\pm 0.47\%$
| 64 | $12.7\%$ | $\pm 0.31\%$ | $10.4\%$ | $\pm 0.20\%$
$\mathbb{P}^{(\text{1D})}_{0,1/2}$ | 16 | $11.7\%$ | $\pm 1.06\%$ | $8.70\%$ | $\pm 0.67\%$
| 32 | $4.02\%$ | $\pm 0.12\%$ | $3.49\%$ | $\pm 0.04\%$
| 64 | $1.75\%$ | $\pm 0.04\%$ | $1.63\%$ | $\pm 0.03\%$
$\mathbb{P}^{(\text{1D})}_{1/8,1/8,1/8,1/8}$ | 16 | $55.2\%$ | $\pm 2.61\%$ | $29.9\%$ | $\pm 0.59\%$
| 32 | $19.5\%$ | $\pm 0.44\%$ | $15.0\%$ | $\pm 0.24\%$
| 64 | $8.55\%$ | $\pm 0.15\%$ | $7.55\%$ | $\pm 0.12\%$
$\mathbb{P}^{(\text{1D})}_{\text{complete}}$ | 16 | $6.63\%$ | $\pm 0.18\%$ | $6.09\%$ | $\pm 0.09\%$
| 32 | $3.25\%$ | $\pm 0.05\%$ | $3.10\%$ | $\pm 0.04\%$
| 64 | $1.43\%$ | $\pm 0.02\%$ | $1.41\%$ | $\pm 0.02\%$
$\mathbb{P}^{(\text{2D})}_{\text{unif}}$ | 12 | $206\%$ | $\pm 42.4\%$ | $66.9\%$ | $\pm 3.26\%$
| 24 | $91.3\%$ | $\pm 14.3\%$ | $47.5\%$ | $\pm 3.00\%$
| 48 | $43.3\%$ | $\pm 5.63\%$ | $30.1\%$ | $\pm 2.29\%$
$\mathbb{P}^{(\text{2D})}_{\text{nonunif}}$ | 12 | $113\%$ | $\pm 20.8\%$ | $52.7\%$ | $\pm 7.91\%$
| 24 | $51.7\%$ | $\pm 7.91\%$ | $33.8\%$ | $\pm 3.08\%$
| 48 | $24.7\%$ | $\pm 3.40\%$ | $19.7\%$ | $\pm 2.05\%$
$\mathbb{P}^{(\text{3D})}_{\text{unif}}$ | 4 | $92.7\%$ | $\pm 30.9\%$ | $48.3\%$ | $\pm 7.25\%$
| 8 | $38.8\%$ | $\pm 9.09\%$ | $27.8\%$ | $\pm 4.12\%$
| 16 | $18.0\%$ | $\pm 3.66\%$ | $15.2\%$ | $\pm 2.36\%$
$\mathbb{P}^{(\text{3D})}_{\text{nonunif}}$ | 4 | $77.8\%$ | $\pm 22.6\%$ | $44.3\%$ | $\pm 7.00\%$
| 8 | $33.4\%$ | $\pm 7.86\%$ | $25.1\%$ | $\pm 4.04\%$
| 16 | $15.7\%$ | $\pm 3.20\%$ | $13.5\%$ | $\pm 2.29\%$
Table 1: Deterministic Value Function Comparison: A performance comparison between the accuracy of the operator augmented estimator $\tilde{\mathbf{x}}$ for the solution $\mathbf{x}$ versus the naive estimator $\hat{\mathbf{x}}$. The error is measured as a percentage with respect to the residual norm of the true solution. Since the error is calculated via Monte Carlo, we provide a $95\%$ confidence interval in the $\pm 2\sigma$ column.  Chain | Samples ($N$) | Naive Error | $\pm 2\sigma$ | Augmented Error | $\pm 2\sigma$
---|---|---|---|---|---
$\mathbb{P}^{(\text{1D})}_{1/4,1/4}$ | 16 | $105\%$ | $\pm 25.4\%$ | $48.4\%$ | $\pm 6.16\%$
| 32 | $30.8\%$ | $\pm 2.95\%$ | $22.8\%$ | $\pm 1.70\%$
| 64 | $12.7\%$ | $\pm 1.00\%$ | $11.1\%$ | $\pm 0.77\%$
$\mathbb{P}^{(\text{1D})}_{1/6,2/6}$ | 16 | $64.1\%$ | $\pm 16.0\%$ | $33.9\%$ | $\pm 5.52\%$
| 32 | $17.2\%$ | $\pm 1.77\%$ | $13.5\%$ | $\pm 1.14\%$
| 64 | $6.84\%$ | $\pm 0.52\%$ | $6.13\%$ | $\pm 0.41\%$
$\mathbb{P}^{(\text{1D})}_{0,1/2}$ | 16 | $11.0\%$ | $\pm 3.30\%$ | $8.49\%$ | $\pm 2.40\%$
| 32 | $3.77\%$ | $\pm 0.33\%$ | $3.34\%$ | $\pm 0.26\%$
| 64 | $1.64\%$ | $\pm 0.11\%$ | $1.54\%$ | $\pm 0.10\%$
$\mathbb{P}^{(\text{1D})}_{1/8,1/8,1/8,1/8}$ | 16 | $55.1\%$ | $\pm 2.62\%$ | $29.9\%$ | $\pm 0.59\%$
| 32 | $19.5\%$ | $\pm 0.44\%$ | $15.0\%$ | $\pm 0.24\%$
| 64 | $8.55\%$ | $\pm 0.15\%$ | $7.55\%$ | $\pm 0.12\%$
$\mathbb{P}^{(\text{1D})}_{\text{complete}}$ | 16 | $6.21\%$ | $\pm 0.29\%$ | $5.81\%$ | $\pm 0.25\%$
| 32 | $2.98\%$ | $\pm 0.13\%$ | $2.89\%$ | $\pm 0.12\%$
| 64 | $1.46\%$ | $\pm 0.06\%$ | $1.43\%$ | $\pm 0.06\%$
$\mathbb{P}^{(\text{2D})}_{\text{unif}}$ | 12 | $25.7\%$ | $\pm 7.51\%$ | $20.4\%$ | $\pm 4.89\%$
| 24 | $11.3\%$ | $\pm 2.93\%$ | $10.2\%$ | $\pm 2.38\%$
| 48 | $5.35\%$ | $\pm 1.31\%$ | $5.07\%$ | $\pm 1.19\%$
$\mathbb{P}^{(\text{2D})}_{\text{nonunif}}$ | 12 | $20.2\%$ | $\pm 6.17\%$ | $16.7\%$ | $\pm 4.89\%$
| 24 | $8.96\%$ | $\pm 2.43\%$ | $8.20\%$ | $\pm 2.04\%$
| 48 | $4.26\%$ | $\pm 1.09\%$ | $4.07\%$ | $\pm 1.00\%$
$\mathbb{P}^{(\text{3D})}_{\text{unif}}$ | 4 | $35.3\%$ | $\pm 35.0\%$ | $26.2\%$ | $\pm 20.8\%$
| 8 | $14.6\%$ | $\pm 11.8\%$ | $12.7\%$ | $\pm 9.11\%$
| 16 | $6.76\%$ | $\pm 5.03\%$ | $6.33\%$ | $\pm 4.43\%$
$\mathbb{P}^{(\text{3D})}_{\text{nonunif}}$ | 4 | $33.4\%$ | $\pm 35.4\%$ | $25.1\%$ | $\pm 21.6\%$
| 8 | $13.9\%$ | $\pm 12.0\%$ | $1.22\%$ | $\pm 9.36\%$
| 16 | $6.43\%$ | $\pm 5.12\%$ | $6.03\%$ | $\pm 4.53\%$
Table 2: Isotropic Random Value Function Comparison: A performance comparison
between the accuracy of the operator augmented estimator $\tilde{\mathbf{x}}$
for the solution $\mathbf{x}$ versus the naive estimator $\hat{\mathbf{x}}$.
The value function is sampled from the prior
$\mathcal{N}(\mathbf{0},\mathbf{I})$. The error is measured as a percentage
with respect to the average residual norm of the true solutions from the
prior. Since the error is calculated via Monte Carlo, we provide a $95\%$
confidence interval in the $\pm 2\sigma$ column.

As we see in Table 1 and Table 2, operator augmentation can provide
significant reductions in error for a variety of different Markov chain
problems, measured both in the reduction of residual norm error for
deterministic right-hand-side (Table 1). As predicted by the theory, the
method also reduces the error in isotropic residual matrix norm, as seen in
(Table 2), but these improvements seem are more marginal. This behavior is
present for all different levels of sample count $N$ we tested. Therefore,
while there are theoretical limitations in the non-symmetric case of operator
augmentation that make the theory less powerful than the symmetric positive
definite case, operator augmentation still functions quite well on the Markov
chain problems we’ve tested it on.

Note that the confidence intervals for the 2D and 3D problems in the isotropic
setting seen in Table 2 are quite large. This is despite the fact that we use
a very large number of samples (256,000 for 2D, and 25,600 for 3D) to estimate
the isotropic error. However, the theory tells us that the naive isotropic
error will always be greater than the augmented isotropic error, so this chart
is provided more-so as a way to gauge the magnitude of the error reduction,
rather than existence of an error reduction.

##  10 Conclusion

We conclude this paper by noting that we have accomplished two main goals.
First, we have investigated the extent to which the symmetric positive
definite operator augmentation theory of [4] can be applied to the general
non-symmetric matrix case. We have found that under the assumptions of noise
symmetry and right-hand-side isotropy, the optimal augmentation factor is
always positive. This means that the act of bootstrapping is relatively safe,
as operator augmentation will damp noise rather than amplify it. We have also
shown that there are pathological situations that make our assumptions
necessary in general. However, in the small noise regime, the noise symmetry
assumption can be discarded.

Second, we have shown empirically that operator augmentation can _still_
reduce error for noisy Markov chain problems, even when the aforementioned
theoretical assumptions are not satisfied. This holds true across a number of
different Markov chains and for both deterministic and isotropic right-hand-
side.

One may continue this work by attempting to apply some form of operator
augmentation to real problems – for example, in control theory or
reinforcement learning, where the underlying Markov decision process may not
be fully known and must be estimated from data. Another more theoretical
possibility would be to investigate if the operator augmentation framework can
be applied to optimization to create optimization algorithms that are less
vulnerable to noise in the objective function, as is common in many real world
applications.

## References

  * [1] G. W. Anderson, A. Guionnet, and O. Zeitouni, An introduction to random matrices, vol. 118, Cambridge university press, 2010.
  * [2] E. J. Candes and Y. Plan, Matrix completion with noise, Proceedings of the IEEE, 98 (2010), pp. 925–936.
  * [3] C. Davis and W. M. Kahan, The rotation of eigenvectors by a perturbation. iii, SIAM Journal on Numerical Analysis, 7 (1970), pp. 1–46.
  * [4] P. A. Etter and L. Ying, Operator augmentation for noisy elliptic systems, arXiv preprint arXiv:2010.09656, (2020).
  * [5] W. James and C. Stein, Estimation with quadratic loss, in Breakthroughs in statistics, Springer, 1992, pp. 443–460.
  * [6] R. Keshavan, A. Montanari, and S. Oh, Matrix completion from noisy entries, in Advances in neural information processing systems, 2009, pp. 952–960.
  * [7] Y. Marzouk, T. Moselhy, M. Parno, and A. Spantini, An introduction to sampling via measure transport, arXiv preprint arXiv:1602.05023, (2016).
  * [8] T. Palmer, G. Shutts, R. Hagedorn, F. Doblas-Reyes, T. Jung, and M. Leutbecher, Representing model uncertainty in weather and climate prediction, Annu. Rev. Earth Planet. Sci., 33 (2005), pp. 163–193.
  * [9] C. Soize, A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics, Journal of sound and vibration, 288 (2005), pp. 623–652.
  * [10] T. Tao, Topics in random matrix theory, vol. 132, American Mathematical Soc., 2012.
  * [11] D. Xiu and J. S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM Journal on Scientific Computing, 27 (2005), pp. 1118–1139.
  * [12] D. Xiu and G. E. Karniadakis, The wiener–askey polynomial chaos for stochastic differential equations, SIAM journal on scientific computing, 24 (2002), pp. 619–644.