# Investigating Numeracy Learning Ability of a Text-to-Text Transfer Model

Kuntal Kumar Pal, Chitta Baral
Department of Computer Science
Arizona State University, Tempe, Arizona, USA,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The transformer-based pre-trained language models have been tremendously
successful in most of the conventional NLP tasks. But they often struggle in
those tasks where numerical understanding is required. Some possible reasons
can be the tokenizers and pre-training objectives which are not specifically
designed to learn and preserve numeracy. Here we investigate the ability of
text-to-text transfer learning model (T5), which has outperformed its
predecessors in the conventional NLP tasks, to learn numeracy. We consider
four numeracy tasks : numeration, magnitude order prediction, finding minimum
and maximum in a series, and sorting. We find that, although T5 models perform
reasonably well in the interpolation setting, they struggle considerably in
the extrapolation setting across all four tasks.

##  1 Introduction

Recent advances in transfer learning in NLP have led to the emergence of pre-
trained models which show a much stronger contextual representation of words
than earlier static word embeddings. They have all performed extremely well in
conventional NLP tasks. Yet, they fail to capture a better understanding of
numbers. Numbers are integral part of natural language texts which can change
the meaning of a sentence. So there is a need for NLP models which can
identify numbers represented in any surface forms like words, floats or
strings (Numeration), understand its values in various context (Magnitude
Order Prediction), compare their values with others (List-MinMax) or able to
rearrange a series of numbers based on its values (Sorting).

The transfer-learned models are pre-trained on huge amount of natural language
texts with specially designed tasks and tokenizers to create stronger word-
embeddings. This causes the numbers embedded in the texts to lose their
meaning and inherent rules of numeracy guiding them Thawani et al. (2021);
Nogueira et al. (2021). This is possibly the reason they perform worse in
numerical reasoning tasks on numbers absent in training data Nogueira et al.
(2021); Wallace et al. (2019).

Figure 1: Examples of Numeracy Tests

In this paper, we test this numeracy learning ability of a text-to-text
transfer learning generative model, T5 Raffel et al. (2020) which has
outperformed its predecessors in conventional NLP tasks. The text-to-text
format of input and output helps the model to generalize all the NLP tasks as
a unified model. We use four numeracy tests both in interpolation (training
and testing on same range of data) and extrapolation settings (training on
lower and testing on higher range of data) and study how much numeracy skill
it can acquire. Figure 1 shows some examples of each of the numeracy tests.

Our contributions in this paper are: (1) Extensive study on three versions of
T5 models (small, base, large) on four numeracy tests in interpolation and
extrapolation settings. (2) Reporting interesting observations in the behavior
of each model version across multiple experimental settings through detailed
manual error analysis. The synthetically generated data and codes are publicly
available111https://github.com/kuntalkumarpal/T5Numeracy for future numeracy
analysis in similar settings.

##  2 Numeracy Tests

We perform four essential numeracy tests to explore model’s ability to
understand numerical values.

Motivation: These four elementary tasks are simple and easy for the models,
since they do not need to generate a completely new number in a different
numerical range (like in mathematical tests : multiplication, division,
exponentiation). Here we evaluate whether the models learn the numeracy tasks
or they simply learn bias from the number range seen in training data.

###  2.1 Numeration

The probability of a number represented in multiple surface forms (word,
scientific, float, integer) increases with the increase in the volume of pre-
training corpus of the language models. It is impractical for an end-to-end
NLP model to semantically parse these numbers accurately and convert into a
single representation to retain its value or reason with. This task tests the
model’s ability to understand word representation of a number and to decode
into integer form.

###  2.2 Magnitude Order Prediction

The task is to identify the order of magnitude of a missing (masked) number
which fits the context of a natural language text. This task is important in
numerical commonsense reasoning Lin et al. (2020) and prompt-based methods Liu
et al. (2021). Here, we do not expect the model to predict the exact number
that fits the context as this may vary in different domains. Instead, this
task tests the model’s ability to understand a missing number’s context and
predict its appropriate range.

###  2.3 List-MinMax

We test the model’s ability to understand numerical values and compare among
them. Given a series of $n$ positive numbers, the task is to find the minimum
and the maximum number. This is the basis of many question answering and
commonsense numerical reasoning dataset like SQuAD Rajpurkar et al. (2016),
DROP Dua et al. (2019) and NUMBERGAME Mishra et al. (2020). We simplify the
task by generating templates so that the models can concentrate on
understanding the task rather than getting confused by the language
complexities.

###  2.4 Sorting

In addition to understanding the values of each number in a series, the model
will have to rearrange them in the correct order through this task, making it
even harder than List-MinMax. Even if a model is successful in the previous
test, it is necessary to identify whether it has actually compared among all
the numbers in the series. Hence, sorting a list of $n$ numbers in ascending
and descending orders ensures that the model compares all the numbers and
rearrange them into two different sequences.

##  3 Experiments

###  3.1 Experimental Setup:

We use T5-SM (small, 60M parameters), T5-BS (base, 220M), T5-LG (large, 770M)
and positive integers for the experiments. The results are average of three
random seeds. We perform experiments in two settings: interpolation (training
and testing on same numerical range) and extrapolation (training on lower and
testing on higher numerical range). The latter helps us to analyze whether a
model has learnt the task, or it has exploited bias in the numerical range of
the training data.

###  3.2 Data Preparation:

# TRAIN $\rightarrow{}$ | 4.9K | 1.3K | 0.9K
---|---|---|---
TP | Model | IN | EX | IN | EX | IN | EX
FL | T5-SM | 45.31 | 0.08 | 1.90 | 0.01 | 0.33 | 0.00
T5-BS | 92.16 | 1.03 | 66.47 | 0.45 | 37.20 | 0.42
T5-LG | 98.06 | 1.91 | 89.49 | 1.96 | 79.48 | 1.58
SP | T5-SM | 69.67 | 39.35 | 26.89 | 1.10 | 0.23 | 0.01
T5-BS | 99.50 | 11.31 | 81.21 | 22.44 | 73.61 | 31.06
T5-LG | 100.00 | 10.05 | 99.97 | 7.35 | 91.59 | 12.92

Table 1: Numeration EM scores w/ split (SP) and w/o split (FL) representation
on 4.9K, 1.3K, 0.9K train-data in Interpolation (IN) and Extrapolation (EX)
settings.

|  | LIST MINIMUM | LIST MAXIMUM
---|---|---|---
# ELEMENTS | 3 | 5 | 10 | 3 | 5 | 10
Range | Model | IN | EX | IN | EX | IN | EX | IN | EX | IN | EX | IN | EX
< 99 | T5-SM | 90.5 | 0.6 | 86.5 | 0.1 | 65.9 | 0.0 | 80.4 | 0.5 | 71.6 | 0.3 | 74.7 | 0.1
T5-BS | 96.2 | 33.9 | 99.1 | 13.0 | 98.2 | 2.8 | 92.3 | 22.7 | 96.8 | 6.0 | 90.4 | 1.1
T5-LG | 100.0 | 22.2 | 99.4 | 2.8 | 100.0 | 0.5 | 100.0 | 29.6 | 100.0 | 13.6 | 100.0 | 2.0
< 999 | T5-SM | 72.6 | 41.8 | 55.5 | 22.2 | 49.9 | 9.7 | 65.3 | 38.4 | 54.8 | 17.5 | 40.0 | 5.2
T5-BS | 91.5 | 67.2 | 92.1 | 42.6 | 80.4 | 27.1 | 89.1 | 65.3 | 90.8 | 47.2 | 88.3 | 25.0
T5-LG | 98.3 | 70.1 | 96.1 | 49.3 | 87.4 | 34.7 | 96.1 | 61.2 | 97.8 | 58.7 | 95.2 | 35.3
< 9999 | T5-SM | 59.1 | 44.7 | 43.5 | 30.4 | 30.7 | 17.1 | 51.2 | 47.0 | 36.0 | 27.0 | 20.9 | 11.1
T5-BS | 89.6 | 68.8 | 86.9 | 53.8 | 85.4 | 38.1 | 87.1 | 58.6 | 83.1 | 43.4 | 81.6 | 29.9
T5-LG | 97.1 | 81.3 | 93.7 | 71.8 | 94.0 | 58.2 | 96.2 | 84.9 | 94.9 | 76.4 | 94.9 | 59.1

Table 2: List-MinMax (series length: 3, 5, 10) in three different number
ranges evaluated as Interpolation (IN) and Extrapolation (EX) exact-match
scores on 1K test data.

Numeration: We create a dataset keeping in mind that at least few examples of
all unique words needed to represent each number, are present in the training
data Trask et al. (2018). In Table 1, interpolation samples are from [0,10K)
and 99K extrapolation samples are from [10K,1000K). We use
num2words222https://github.com/savoirfairelinux/num2words for generating word-
form of each integer. To simulate fewer shot setting, we carefully craft two
smaller training sets taking only 20% and 10% data. We show two number
representation schemes with split-digits (SP) and without split (FL)
hypothesizing that for a generative model it would be easier to correctly
generate individual digits instead of full integer at once.

Datasets $\rightarrow$ | AT | MC
---|---|---
Models $\downarrow$ |  $\mu$F1 |  $m$F1 |  $\mu$F1 |  $m$F1
LR | 62.49 | 30.81 | 71.25 | 60.80
CNN | 69.27 | 35.96 | 77.17 | 58.49
GRU | 70.92 | 38.43 | 78.25 | 58.08
BiGRU | 71.49 | 39.94 | 80.16 | 62.74
CRNN | 69.50 | 36.15 | 78.00 | 64.62
CNN-capsule | 63.11 | 29.41 | 75.89 | 59.22
GRU-capsule | 70.73 | 33.57 | 77.36 | 64.71
BiGRU-capsule | 71.49 | 34.18 | 77.97 | 64.34
BiLSTM-DICE | 75.56 | 46.80 | - | -
T5-SM | 69.87 | 31.36 | 66.11 | 34.68
T5-BS | 78.06 | 40.04 | 72.22 | 47.44
T5-LG | 81.40 | 44.64 | 80.29 | 59.16
Table 3: Magnitude Order Prediction for Market Comments (MC) and Article
Titles (AT) datasets of numeracy600K in micro-F1 ($\mu$F1) and macro-F1
($m$F1). Best score is in bold and second-best is underlined.

Magnitude Order Prediction: For this task we work on Numeracy600K Chen et al.
(2019) dataset. We consider this as a mask prediction task. We train models to
find the exact number that fits the mask. Then, we map the predicted numbers
into its magnitude order, save the model based on best magnitude order and
calculate the evaluation metrics on test data. Since this is a generation task
we reject those answers which are not valid floating point numbers. The
baseline results in Table 3 are from Chen et al. (2019); Sundararaman et al.
(2020). We also consider extrapolation setting by showing the cross-domain
performance (train on market comments and test on article title and vice-
versa) in Table 4.

Train on $\rightarrow$ | AT | MC
---|---|---
Models $\downarrow$ |  $\mu$F1 |  $m$F1 |  $\mu$F1 |  $m$F1
BiGRU | 25.59 | 10.58 | 31.38 | 11.08
T5-SM | 28.88 | 12.04 | 37.35 | 10.81
T5-BS | 35.53 | 14.48 | 31.51 | 12.25
T5-LG | 50.18 | 21.24 | 38.43 | 12.32
Table 4: Cross Domain (Extrapolation) Tests of Order Prediction. Train on MC,
test on AT and vice-versa.

List Min-Max & Sort: We experiment on three different number ranges: [0,100),
[0,1K), [0,10K). For interpolation tests, the numbers in the test data are
from the same ranges. The extrapolation numbers are from the maximum of
respective ranges to 100K. To prevent the model’s bias on number lengths, we
bring them closer following prior work Wallace et al. (2019). We extend the
experiment on a series of 3, 5 and 10 numbers (for each range) to study how
each of the models behave with increasing series length. We consider the same
data for sorting experiments as well. The results are in Table 2 for List-
MinMax and Table 5 for List-Sort.

|  | LIST-SORT ASCENDING | LIST-SORT DESCENDING
---|---|---|---
# ELEMENTS | 3 | 5 | 10 | 3 | 5 | 10
Range | Model | IN | EX | IN | EX | IN | EX | IN | EX | IN | EX | IN | EX
< 99 | T5-SM | 54.0 | 12.4 | 7.6 | 0.0 | 0.0 | 0.0 | 56.0 | 12.6 | 5.9 | 0.4 | 0.0 | 0.0
T5-BS | 80.6 | 12.2 | 87.2 | 0.0 | 0.4 | 0.0 | 84.3 | 12.9 | 75.5 | 0.0 | 6.2 | 0.0
T5-LG | 100.0 | 5.8 | 99.9 | 0.0 | 69.7 | 0.1 | 100.0 | 13.1 | 96.6 | 0.1 | 57.6 | 0.1
< 999 | T5-SM | 32.6 | 15.1 | 1.4 | 0.6 | 0.0 | 0.0 | 38.0 | 22.3 | 3.4 | 1.3 | 0.0 | 0.0
T5-BS | 74.7 | 45.7 | 64.0 | 8.0 | 12.5 | 0.0 | 73.1 | 42.0 | 62.6 | 9.6 | 16.8 | 0.1
T5-LG | 95.1 | 64.2 | 91.8 | 16.8 | 61.9 | 1.7 | 94.7 | 63.5 | 92.5 | 25.7 | 61.2 | 1.6
< 9999 | T5-SM | 23.4 | 17.1 | 1.0 | 0.1 | 0.0 | 0.0 | 30.4 | 21.2 | 0.7 | 0.4 | 0.0 | 0.0
T5-BS | 63.1 | 45.5 | 51.1 | 12.7 | 15.0 | 0.2 | 59.8 | 43.9 | 51.4 | 12.4 | 14.3 | 0.3
T5-LG | 94.5 | 76.0 | 87.4 | 43.2 | 74.6 | 12.6 | 94.2 | 76.1 | 86.1 | 44.4 | 75.6 | 11.9

Table 5:  List-Sort (Ascending & Descending) on series lengths: 3, 5, 10 in
three different integer ranges evaluated as Interpolation (IN) and
Extrapolation (EX) exact-match scores on 1K test data.

##  4 Results and Error Analysis

Table 1 shows, all versions of T5 benefit when they are trained with split
representation. When trained with 4.9K data, T5-SM gains 24% points in
interpolation evaluation where T5-LG gains only 2%. None of the models perform
well on unseen number data ranges. In fewer shot interpolation settings
however, only the T5-LG model maintains its performance beyond 90% which is
not surprising because of its large parameter-space. We noticed that the best
model could only partially decode numbers having multiple zeros (Figure 2). In
the first example, the model predicts an extra seven and in the second
(extrapolation), it ignored the key word ‘hundred’ as it attempts to fit this
unseen data into a similar seen number range (4 digits).

In magnitude order prediction (Table 3), T5-LG’s performance improves by 5
$\mu$F1 in article title. For extrapolation (Table 4), all T5 versions beats
previous estimates (BiGRU) by at most 25%. This shows that T5 can learn robust
numeric representations based on contexts. Both the samples in Fig 2 are hard
as they need prior explicit knowledge. Yet they are able to predict numbers in
similar feasible ranges. This shows that the model is not randomly assigning
magnitude but has learnt based on the domain and context. We found that, the
best T5 model predicted an order of 1 instead of 2 for market and article data
making a maximum error of 39.07% and 33.59% respectively.

Figure 2: Two incorrect predictions for each task.

Table 2 shows List-MinMax results. Both T5-BS and T5-LG perform over 80%
across all ranges and series lengths. T5-SM however, degrades in performance
as the range increases along with the list size. As the model learns more
variations in numbers, the extrapolation performance increases to a max of 81%
(List-Min) and 84.9% (List-Max). But the performance drops as series length
increases. The best model predicted second minimum and maximum element in the
examples of Fig 2.

From the sorting results (Table 5), we see T5-SM performance drops (18-22%
from 2-3 digits, 8-9% from 3-4) as number ranges increase across series length
of 3. T5-SM fails to generate a single correct order for a series of 10
elements and achieves less than 10% success in 5-element series across all
ranges. This degrading performance can be attributed to its mere 60M parameter
space. As the number of parameters keep increasing the models performs
consistently across each of 3, 5, 10 elements in series, both for
interpolation and extrapolation settings. With the increasing range of
training data, the models become more robust to extrapolated numbers across
all series lengths with 8-30% change in ascending order and 7-20% change in
descending order. Finally, for sorting, we find a variety of incorrect
predictions: missing order of one element, omission of one and two elements or
repeating a particular element.

Overall, none of the models were able to perform well on extrapolation samples
showing the inherent rules of numeracy is difficult for these models to learn.
But, it also shows, more variations in numbers (increasing the range) help
them perform better in extrapolation setting. The smaller model’s limited
parameter-space affects its performance in all four tasks whereas larger
models are able to pick up some numeracy skills through training. We show more
predictions in Figure 3, 4, 5, 6.

Figure 3: Some predictions for Numeration task. Figure 4: Magnitude Order
Prediction Examples.

Analysis of NT5: We test with the NT5 Yang et al. (2021) model on all our
experiments and compared the results with T5-small. For the Numeration task
with the split number representation NT5 performed 73.07 (accuracy), a 4%
improvement over T5. The performance however did not improve for the MinMax
and Sorting tasks. For 3-element sorting it dropped by 10-20%. In the
Magnitude Order Prediction, we find the cross-domain (extrapolation) $\mu$F1
score increases by 5-7% while in-domain decreases by 3-6%. This might be
because NT5 has seen more variety of contexts of numbers and can generalize
well on this task.

Figure 5: Some predictions for List-MinMax task. Figure 6: Some predictions
for List-Sort task.

##  5 Related Works

Numeracy Tests: Multiple numeracy tests have been proposed to evaluate the
static word embeddings Naik et al. (2019) like GloVe, Word2Vec, FastText and
contextual embeddings Wallace et al. (2019) like BERT through probing tasks
like numeration, magnitude comparison, addition, list-maximum. Multilingual
numeration Johnson et al. (2020) tests have been performed by probing models
like DistilBERT, XLM, and BERT. CNN, BiGRU models have been shown to perform
well in magnitude order prediction Chen et al. (2019) and T5 on addition and
subtraction tasks Nogueira et al. (2021) through training on similar texts.
We, however focus on studying how much text-to-text transfer models (T5) can
learn across four fundamental numeracy tasks in samples containing both in-
domain and out-of-domain numerical ranges.

Specially Designed Models: NALU Trask et al. (2018), NAU and NMU Madsen and
Johansen (2020), numBERT Zhang et al. (2020), GenBERT Geva et al. (2020), NT5
Yang et al. (2021) have emerged in the last few years to incorporate
arithmetic skills into models through specially designed architecture or fine-
tuning tasks which improves the performance in synthetic arithmetic or crowd-
sourced numerical reasoning tasks like DROP.

Numerical Embeddings: There are limited prior works in numeracy aware
embeddings which show good performance in extrapolation setting. One approach
Jiang et al. (2019) represents numerals as a weighted average of prototype
numeral embeddings obtained using either self organizing map or Gaussian
Mixture models. DICE Sundararaman et al. (2020) is a deterministic numeral
embedding approach, independent of corpus, which preserves the relative
magnitude between two numerals and their embeddings.

##  6 Conclusion & Future Works

We show that text-to-text models are able to learn numeracy quite well in an
interpolation setting. Our extensive experiments show that T5 models struggle
to learn with numbers outside training data ranges. We believe that, to make
further progress in transfer learning, models need to achieve such elementary
numeracy skills and this gap between interpolation and extrapolation
performance needs to be reduced. We are of the opinion that, adding more data
would not bridge this gap since domain of numbers is open. However, special
pre-training objectives for digits rather than whole numbers can be designed
to teach the inherent numeracy to models. In future, we intend to explore
these objectives centered around preserving numeracy rules in transfer-learned
models to generalize between in-domain and out-of-domain numbers.

## Acknowledgement

The authors acknowledge support from DARPA grant number FA875019C0003 for this
project.

## Ethical Considerations

In this paper, we analyze performance of three publicly available T5 models on
four numeracy tasks. For Magnitude Order Prediction task we use publicly
available dataset, Numeracy600K. We synthetically create the data for rest of
the tasks.

## References

  * Chen et al. (2019) Chung-Chi Chen, Hen-Hsen Huang, Hiroya Takamura, and Hsin-Hsi Chen. 2019.  Numeracy-600k: learning numeracy for detecting exaggerated information in market comments.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 6307–6313.
  * Dua et al. (2019) Dheeru Dua, Yizhong Wang, Pradeep Dasigi, Gabriel Stanovsky, Sameer Singh, and Matt Gardner. 2019.  Drop: A reading comprehension benchmark requiring discrete reasoning over paragraphs.  _arXiv preprint arXiv:1903.00161_.
  * Geva et al. (2020) Mor Geva, Ankit Gupta, and Jonathan Berant. 2020.  Injecting numerical reasoning skills into language models.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 946–958, Online. Association for Computational Linguistics.
  * Jiang et al. (2019) Chengyue Jiang, Zhonglin Nian, Kaihao Guo, Shanbo Chu, Yinggong Zhao, Libin Shen, and Kewei Tu. 2019.  Learning numeral embeddings.  _arXiv preprint arXiv:2001.00003_.
  * Johnson et al. (2020) Devin Johnson, Denise Mak, Andrew Barker, and Lexi Loessberg-Zahl. 2020.  Probing for multilingual numerical understanding in transformer-based language models.  In _Proceedings of the Third BlackboxNLP Workshop on Analyzing and Interpreting Neural Networks for NLP_ , pages 184–192, Online. Association for Computational Linguistics.
  * Lin et al. (2020) Bill Yuchen Lin, Seyeon Lee, Rahul Khanna, and Xiang Ren. 2020.  Birds have four legs?! numersense: Probing numerical commonsense knowledge of pre-trained language models.  In _Proceedings of EMNLP_.  To appear.
  * Liu et al. (2021) Pengfei Liu, Weizhe Yuan, Jinlan Fu, Zhengbao Jiang, Hiroaki Hayashi, and Graham Neubig. 2021.  Pre-train, prompt, and predict: A systematic survey of prompting methods in natural language processing.
  * Madsen and Johansen (2020) Andreas Madsen and Alexander Rosenberg Johansen. 2020.  Neural arithmetic units.  _arXiv preprint arXiv:2001.05016_.
  * Mishra et al. (2020) Swaroop Mishra, Arindam Mitra, Neeraj Varshney, Bhavdeep Sachdeva, and Chitta Baral. 2020.  Towards question format independent numerical reasoning: A set of prerequisite tasks.  _arXiv preprint arXiv:2005.08516_.
  * Naik et al. (2019) Aakanksha Naik, Abhilasha Ravichander, Carolyn Rose, and Eduard Hovy. 2019.  Exploring numeracy in word embeddings.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3374–3380.
  * Nogueira et al. (2021) Rodrigo Nogueira, Zhiying Jiang, and Jimmy Li. 2021.  Investigating the limitations of the transformers with simple arithmetic tasks.  _arXiv preprint arXiv:2102.13019_.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ , 21(140):1–67.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang. 2016.  SQuAD: 100,000+ questions for machine comprehension of text.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 2383–2392, Austin, Texas. Association for Computational Linguistics.
  * Sundararaman et al. (2020) Dhanasekar Sundararaman, Shijing Si, Vivek Subramanian, Guoyin Wang, Devamanyu Hazarika, and Lawrence Carin. 2020.  Methods for numeracy-preserving word embeddings.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 4742–4753.
  * Thawani et al. (2021) Avijit Thawani, Jay Pujara, Pedro A Szekely, and Filip Ilievski. 2021.  Representing numbers in nlp: a survey and a vision.  _arXiv preprint arXiv:2103.13136_.
  * Trask et al. (2018) Andrew Trask, Felix Hill, Scott E Reed, Jack Rae, Chris Dyer, and Phil Blunsom. 2018\.  Neural arithmetic logic units.  In _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc.
  * Wallace et al. (2019) Eric Wallace, Yizhong Wang, Sujian Li, Sameer Singh, and Matt Gardner. 2019.  Do NLP models know numbers? probing numeracy in embeddings.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 5307–5315, Hong Kong, China. Association for Computational Linguistics.
  * Yang et al. (2021) Peng-Jian Yang, Ying Ting Chen, Yuechan Chen, and Daniel Cer. 2021.  Nt5?! training t5 to perform numerical reasoning.  _arXiv preprint arXiv:2104.07307_.
  * Zhang et al. (2020) Xikun Zhang, Deepak Ramachandran, Ian Tenney, Yanai Elazar, and Dan Roth. 2020.  Do language embeddings capture scales? In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pages 4889–4896, Online. Association for Computational Linguistics.

##  Appendix A Appendix

###  A.1 Data Statistics Experimental Setup

Numeration: We have 4906, 2097, 2997 in train, dev and test respectively. We
make sure that all numbers within 10K are present in any of train, dev or
test. For extrapolation we select 1K integers randomly from every 10K range
from [10K,1000K) making it a total of 99K.

Magnitude Order Prediction: For this data we consider 450K, 50K and 100K
samples for train, dev and test data respectively from each of market comments
and article titles data.

List-Sort: We consider both the task of arranging in ascending and descending
orders since if a series is already sorted in ascending order the model can
directly predict by copying it from the given input.

Figure 7: More predictions for Numeration task.

###  A.2 Hyperparameters

For all the experiments we use maximum sequence length of 128 and 256 for
question context. The maximum sequence length of the answers is kept as [5,
10, 20, 25] for different tasks. We ran for 20 epochs and save a model based
on validation EM performance. Our training and validation batch size varies
between [2, 4, 8, 16, 32] based on the experiment. We work on 4 Tesla V100
GPUs. We use AdamW optimizer and StepLR scheduler with step size of 2,
learning rate of 5e-5 and gamma of 0.1.

Figure 8: More Magnitude Order Prediction Examples.

###  A.3 Results and Error analysis

Magnitude Order Prediction: We also experimented with zero-shot magnitude
order predictions. We found 553 and 8783 exact-matches out of 100K test data
using T5-large which shows that the performance is very poor without proper
fine-tuning. We show some more predictions of the best performing T5 model in
Figure 7, 8, 9, 10.

Figure 9: More predictions for List-MinMax task. Figure 10: More predictions
for List-Sort task.


# Iterate Averaging, the Kalman Filter, and 3DVAR for Linear Inverse Problems

Felix G. Jones  and  Gideon Simpson<EMAIL_ADDRESS>Department of
Mathematics, Drexel University, Philadelphia, PA, 19103, USA

###### Abstract.

It has been proposed that classical filtering methods, like the Kalman filter
and 3DVAR, can be used to solve linear statistical inverse problems. In the
work of Iglesias, Lin, Lu, & Stuart (2017), [10], error estimates were
obtained for this approach. By optimally tuning a regularization parameter in
the filters, the authors were able to show that the mean squared error could
be systematically reduced.

Building on the aforementioned work of Iglesias, Lin, Lu, & Stuart, we prove
that by (i) considering the problem in a weaker norm and (ii) applying simple
iterate averaging of the filter output, 3DVAR will converge in mean square,
unconditionally on the choice of parameter. Without iterate averaging, 3DVAR
cannot converge by running additional iterations with a fixed choice of
parameter. We also establish that the Kalman filter’s performance in this
setting cannot be improved through iterate averaging. We illustrate our
results with numerical experiments that suggest our convergence rates are
sharp.

###### Key words and phrases:

Kalman filter, 3DVAR, statistical inverse problems, averaging

###### 2010 Mathematics Subject Classification:

93E11, 65J22, 47A52

##  1\. Introduction

The focus of this work is on the inverse problem

(1.1) $y=Au^{\dagger}+\eta,$

where, given the noisy observation $y$ of $Au^{\dagger}$, we wish to infer
$u^{\dagger}$. In our setting, $A:X\to Y$ is a compact operator between
separable Hilbert spaces and $\eta\sim N(0,\gamma^{2}I)$ is white noise,
modelling measurement error. This problem is well-known to be ill-posed in the
infinite dimensional setting, as $A$ has an unbounded inverse. Methods of
solution include the use of regularized Moore-Penrose inverses and, subject to
the introduction of a prior, Bayesian formulations, [19, 6, 4, 25, 24, 12, 5].

In [10], a key inspiration for the present work, Iglesias, Lin, Lu, & Stuart
considered two classical filtering algorithms, the Kalman filter and 3DVAR,
with the goal of using them to solve (1.1). The filtering methodology for
(1.1) requires the introduction, conceptually, of the artificial dynamical
system

(1.2a) $\displaystyle u_{n}$ $\displaystyle=u_{n-1},\quad u_{0}=u^{\dagger},$
(1.2b) $\displaystyle y_{n}$
$\displaystyle=Au_{n}+\eta_{n},\quad\eta_{n}\overset{\operatorname{i.i.d.}}{\sim}N(0,\gamma^{2}I).$

Here, at algorithmic time step $n$, $u_{n}$ is the quantity of interest, and
$y_{n}$ is the noisy observation. Having ascribed a notion of time to the
problem, we can then apply a filter. This provides a mechanism for estimating
$u^{\dagger}$ in (1.1) in an online setting, where a sequence of i.i.d.
observations, $\\{y_{n}\\}$, is available. This corresponds to “Data Model 1”
of [10].

Amongst the key results of [10], reviewed in detail below, is that under
sufficiently strong assumptions, the Kalman filter will recover the truth in
mean square, unconditionally on the choice of the scalar regularization
parameter. Under somewhat weaker assumptions, the error will only be bounded,
though through minimax selection of a scalar parameter, an optimal error can
be achieved for a given number of iterations, allowing the error to be driven
to zero.

3DVAR is a simplification of Kalman that is demonstrated to have, at best,
bounded error, though, again, through minimax parameter tuning, it can perform
comparably to Kalman. Kalman is more expensive than 3DVAR, as it requires
updating an entire covariance operator at each iteration. For finite
dimensional approximations, this may require costly matrix-matrix
multiplications at each iterate.

Here, by working in a weaker norm and averaging the iterates, we are able to
establish that 3DVAR will unconditionally converge in mean square for all
admissible filter parameters. Such weaker convergence was also considered in
[3], for a related problem on 4DVAR. Further, we show that this simple iterate
averaging cannot improve the performance of the Kalman filter.

###  1.1. Filtering Algorithms

The Kalman filter is a probabilistic filter that estimates a Gaussian
distribution, $N(m_{n},C_{n})$, for $u^{\dagger}$ at each iterate. Given a
starting mean and covariance, $m_{0}$ and $C_{0}$, the updates are as follows:

(1.3a) $\displaystyle m_{n}$ $\displaystyle=K_{n}y_{n}+(I-K_{n}A)m_{n-1},$
(1.3b) $\displaystyle C_{n}$ $\displaystyle=(I-K_{n}A)C_{n-1},$ (1.3c)
$\displaystyle K_{n}$
$\displaystyle=C_{n-1}A^{\ast}(AC_{n-1}A^{\ast}+\gamma^{2}I)^{-1}.$

Here, $K_{n}$ is the so-called “Kalman gain.” $m_{n}$ is a point estimate of
$u^{\dagger}$.

While Kalman is a probabilistic filter, 3DVAR is not. It is obtained by
applying Kalman with a static covariance operator
$C_{n}=\frac{\gamma^{2}}{\alpha}\Sigma$ for some predetermined operator
$\Sigma$:

(1.4a) $\displaystyle u_{n}$
$\displaystyle=\mathcal{K}y_{n}+(I-\mathcal{K}A)u_{n-1},$ (1.4b)
$\displaystyle\mathcal{K}$
$\displaystyle=(A^{\ast}A+\alpha\Sigma^{-1})^{-1}A^{\ast}.$

We refer the reader to [14, 25, 23, 7], and references therein, for a thorough
discussion and analysis of these classical filtering methods and their
extensions.

Indeed, several important extensions of these classical methods that have
appeared in the literature have also been directly applied to statistical
inverse problems like (1.1), along with its nonlinear variation,
$y=\mathcal{G}(u^{\dagger})+\eta$. In particular, the ensemble Kalman filter
(EnKF), using an ensemble of replicas of the problem, has been successfully
applied to solve such problems in [9, 8]. See, for instance, [14, 21, 22], for
additional details and analysis of EnKF. We also mention [3], which uses
similar ideas with 4DVAR.

Continuous in time analogs of these methods and problems also exist, resulting
in the Kalman-Bucy filter and continuous in time 3DVAR, [14, 18, 25]. In [15],
these were used to solve the continuous in time analog of (1.2)

(1.5a) $\displaystyle du$ $\displaystyle=0,$ (1.5b) $\displaystyle dy$
$\displaystyle=Audt+d\eta,$

where $\eta(t)$ is now a Weiner process in the appropriate function space,
[2].

###  1.2. Key Assumptions and Prior Results

In [10], the following assumptions were invoked.

######  Assumption 1.

  1. (1)

$C_{0}=\frac{\gamma^{2}}{\alpha}\Sigma$ with
$\operatorname{Ran}(\Sigma^{\frac{1}{2}})\subset\operatorname{Dom}(A)$,
$\alpha>0$, and $\Sigma$ a self-adjoint positive definite trace class operator
with $\Sigma^{-1}$ densely defined.

  2. (2)

$\Sigma$ induces a Hilbert scale, and there exist constants $C>1$, $\nu>0$
such that $A$ induces an equivalent norm:

(1.6) $C^{-1}\|x\|_{\nu}\leq\|Ax\|\leq
C\|x\|_{\nu},\quad\|\bullet\|_{\nu}=\|\Sigma^{\frac{\nu}{2}}\bullet\|.$

  3. (3)

The initial error is sufficiently “smooth,”

(1.7) $m_{0}-u^{\dagger}\in\operatorname{Dom}(\Sigma^{-\frac{s}{2}}),\quad
0\leq s\leq\nu+2,$

where we replace $m_{0}$ with $u_{0}$ in the case of 3DVAR in the above
expression.

Under this first set of assumptions, Iglesias et al. established

######  Theorem 1.1 (Theorem 4.1 of [10]).

The Kalman filter admits the mean square error bound

$\mathbb{E}[\|m_{n}-u^{\dagger}\|^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{s}{\nu+1}}+\frac{\gamma^{2}}{\alpha}\operatorname{Tr}\Sigma$

and

######  Theorem 1.2 (Theorem 5.1 of [10]).

3DVAR admits the mean square error bound

$\mathbb{E}[\|u_{n}-u^{\dagger}\|^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{s}{\nu+1}}+\frac{\gamma^{2}}{\alpha}\operatorname{Tr}\Sigma\log
n.$

At fixed values of $\alpha$, Theorems 1.1 and 1.2 preclude convergence, and,
in the case of 3DVAR, the error may even grow. However, there are two free
parameters: the number of iterations $n$ and the regularization parameter
$\alpha$. Indeed, within a Bayesian framework, $\alpha$ can be interpreted as
the strength of a prior relative to a likelihood. For a fixed number of
iterations, $n$, $\alpha$ can be tuned to minimize the error. Indeed, the
error can be made arbitrarily small by selecting a sufficiently large $n$ with
the optimal $\alpha$.

However, in both Theorems 1.1 and 1.2, there is an unknown constant. If the
error at the given, optimal choice of $\alpha$ for a given $n$ is inadequate,
one must obtain additional data, update the value of $\alpha$, and rerun the
algorithm. A benefit of the present work is that, by using iterate averaging,
the error of 3DVAR can always be reduced with additional iterates, without
necessarily altering $\alpha$ and discarding previously computed iterations.
We will revisit the minimax estimates under a simultaneous diagonalization
assumption.

Indeed, stronger results were obtained in [10] subject to the simultaneous
diagonalization assumption:

######  Assumption 2.

  1. (1)

$\Sigma$ and $A^{\ast}A$ simultaneously diagonalize against the set
$\\{\varphi_{i}\\}$ with respective eigenvalues $\sigma_{i}$ and $a_{i}^{2}$,
and these eigenvalues satisfy

(1.8) $\sigma_{i}=i^{-1-2\epsilon},\quad a_{i}\asymp
i^{-p},\quad\epsilon>0,\quad p>0.$

  2. (2)

$m_{0}=0$ (or $u_{0}$ in 3DVAR) and $u^{\dagger}$ satisfies, for $0<\beta\leq
1+2\epsilon+2p$,

(1.9) $\sum_{i=1}^{\infty}i^{2\beta}|u^{\dagger}_{i}|^{2}<\infty.$

With this, Iglesias et al. obtain

######  Theorem 1.3 (Theorem 4.2 of [10]).

Under Assumption 2, for the Kalman filter,

$\mathbb{E}[\|m_{n}-u^{\dagger}\|^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{2\beta}{1+2\epsilon+2p}}+\gamma^{2}n^{-\frac{2\epsilon}{1+2\epsilon+2p}}\alpha^{-\frac{1+2p}{1+2\epsilon+2p}}$

and

######  Theorem 1.4 (Theorem 5.2 of [10]).

Under Assumption 2, for 3DVAR,

$\mathbb{E}[\|u_{n}-u^{\dagger}\|^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{2\beta}{1+2\epsilon+2p}}+C\gamma^{2}\alpha^{-\frac{1+2p}{1+2\epsilon+2p}}.$

Now the Kalman filter will converge at any choice of parameter, while 3DVAR
has at worst a bounded error. Again, $\alpha$ can be tuned so as to obtain a
minimax convergence rate. Indeed, in the setting where one has a fixed number
of $n$ samples, at the optimal value of $\alpha$, Theorems 1.3 and 1.4 lead to
the estimates (also found in [10]):

(1.10) $\displaystyle\mathbb{E}[\|m_{n}-u^{\dagger}\|^{2}]\lesssim$
$\displaystyle n^{-\frac{2\beta}{1+2\beta+2p}},$ (1.11)
$\displaystyle\mathbb{E}[\|u_{n}-u^{\dagger}\|^{2}]\lesssim$ $\displaystyle
n^{-\frac{2\beta}{1+2\beta+2p+2\epsilon}}\log n,$

where the first expression is for Kalman and the second is for 3DVAR. Similar
expressions are also available in the general case for Theorems 1.1 and 1.2.

Thus far, we have discussed the study of problem (1.1) in a sequential
setting, where the data, $\\{y_{n}\\}$, is assimilated one sample at a time.
In some settings, a static, fixed, number of samples, $n$, may be available
together. Instead of (1.1), we might then examine

(1.12)
$\begin{split}\bar{y}_{n}=Au^{\dagger}+\bar{\eta}_{n},\quad\bar{\eta}_{n}\sim
N(0,\tfrac{\gamma^{2}}{n}I),\\\
\bar{y}_{n}=\frac{1}{n}\sum_{k=1}^{n}y_{k},\quad\bar{\eta}_{n}=\frac{1}{n}\sum_{k=1}^{n}\eta_{k}.\end{split}$

The variance of the noise has been reduced by a factor of $n$. This can be
solved using a regularized approximation of $A^{+}$ to obtain
$\bar{u}_{n,\alpha}$. Under suitable assumptions and identifying the optimal
$\alpha=\alpha_{\star}(n)$, one can obtain (see, for instance, [1, 19, 16, 17,
12, 26])

(1.13) $\mathbb{E}[\|\bar{u}_{n,\alpha_{\star}(n)}-u^{\dagger}\|^{2}]\lesssim
n^{-\frac{2\beta}{1+2p+2\beta}}$

This precisely corresponds to the minimax solution of Kalman, (1.10), while
there is a loss for 3DVAR, (1.11). Note that this is only for the $t=0$ norm.
A generalization to the $t<0$ norm is covered in [16] and for
$t(1+2\epsilon)\leq 2p$ in [17]. As we are principally interested in the
general $t\geq 0$ case, we state and prove our own version of theorem below
using a spectral cutoff regularization.

###  1.3. Main Results

The main results of this paper are contained in the following theorems.

First, we have the elementary result that 3DVAR, without averaging, cannot
converge at fixed parameter choices:

######  Theorem 1.5.

Under Assumption 1 in dimension one, if $u_{n}$ is generated by 3DVAR, then

$\mathbb{E}[|u_{n}-u^{\dagger}|^{2}]\geq\gamma^{2}\mathcal{K}^{2}.$

As the method cannot converge in dimension one, it has no hope of converging
in higher dimensions. By time averaging,

(1.14)
$\bar{u}_{n}=\frac{1}{n}\sum_{k=1}^{n}u_{k}=\frac{1}{n}u_{n}+\frac{n-1}{n}\bar{u}_{n-1},$

we can obtain convergence for all $\alpha>0$:

######  Theorem 1.6.

Under Assumption 1, fix $t\in[0,\nu]$ and $\tau_{\rm v}\in[0,1]$, and, having
set these indices, assume that $\Sigma^{t+1-\tau_{\rm v}(1+\nu)}$ is trace
class. Then

$\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{s+t}{1+\nu}}\|z_{0}\|^{2}+\frac{\gamma^{2}}{\alpha}\operatorname{Tr}(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)})\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}$

where $z_{0}$ is the solution to

(1.15)
$\Sigma^{-\frac{1}{2}}(u_{0}-u^{\dagger})=(B^{\ast}B)^{\frac{s-1}{2(1+\nu)}}z_{0}$

and $B=A\Sigma^{\frac{1}{2}}$.

We will repeatedly make use of the operator

(1.16) $B=A\Sigma^{\frac{1}{2}}$

throughout this work. The existence of $z_{0}$ in (1.15) is a consequence of
Assumption 1 on the initial error and an equivalence of spaces result
encapsulated in Proposition 2.3, given below.

The motivation for time averaging comes from two related problems. First,
formally, (1.4), has the structure of an AR(1) process, [23]. Under typical
assumptions, an AR(1) process will not converge to a fixed value, but instead,
sample an invariant distribution. Consequently, the time average will converge
to the mean, with respect to this invariant distribution. Another motivation
comes from the stochastic root finding problem and the Robbins-Monro
algorithm. In [20], Polyak & Juditsky proved that by time averaging the
sequence of of estimates generated by Robbins-Monro, the convergence rate
could be improved. See, also, [13].

As a consequence of Theorem 1.7, we will have unconditional mean squared
convergence of the iterate averaged value, $\bar{u}_{n}$, provided:

  * •

We study the problem in a sufficiently weak weighted space ($t>0$) and/or have
sufficiently smooth data ($s>0$);

  * •

$\Sigma$ has a sufficiently well behaved spectrum, allowing $\tau_{\rm v}>0$.
Note that taking $\tau_{\rm v}=t/(1+\nu)$ will not require additional
assumptions on $\Sigma$, but will require $t>0$ for convergence.

We emphasize that iterate averaging is a post-processing step, requiring no
modification of the underlying 3DVAR iteration.

We introduce a modified version of Assumption 2,

######  Assumption $2^{\prime}$.

  1. (1)

$\Sigma$ and $A^{\ast}A$ simultaneously diagonalize against the set
$\\{\varphi_{i}\\}$ with respective eigenvalues $\sigma_{i}$ and $a_{i}^{2}$,
and these eigenvalues satisfy

(1.17) $\sigma_{i}\asymp i^{-1-2\epsilon},\quad a_{i}\asymp
i^{-p},\quad\epsilon>0,\quad p>0.$

  2. (2)

For $\beta\geq 0$, the initial error, $u_{0}-u^{\dagger}$, satisfies the
condition

(1.18) $\sum_{i=1}^{\infty}i^{2\beta}|u_{0,i}-u^{\dagger}_{i}|^{2}<\infty.$

Condition (1.18) on the initial error will automatically be satisfied if
$u_{0}$ and $u^{\dagger}$ are, separately, sufficiently smooth. The
assumptions of (1.17) and (1.18) are equivalent to those of (1.6) and (1.7)
under the identifications:

$\nu(1+2\epsilon)={2p},\quad s(1+2\epsilon)={2\beta}.$

In contrast to Assumption 2, no upper bound on $\beta$ is necessary.

######  Theorem 1.7.

Under Assumption $2^{\prime}$, and having fixed a choice of
$\left\|\bullet\right\|_{t}$ norm with $t\geq 0$, assume $\tau_{\rm
b},\tau_{\rm v}\in[0,1]$ satisfy

(1.19a) $\displaystyle\tau_{\rm b}$
$\displaystyle\leq\frac{t(1+2\epsilon)+2\beta}{2(1+2\epsilon+2p)}\equiv\bar{\tau}_{\rm
b}$ (1.19b) $\displaystyle\tau_{\rm v}$
$\displaystyle<\frac{t(1+2\epsilon)+2\epsilon}{1+2\epsilon+2p}\equiv\bar{\tau}_{\rm
v}$

then,

$\begin{split}\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-2\tau_{\rm
b}}+\frac{\gamma^{2}}{\alpha}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm
v}}\end{split}$

While our results in both the general and diagonal case establish
unconditional convergence for any choice of $\alpha$ for the iterate averaged
3DVAR, in a practical setting, there may only be $n$ iterates available. One
might then ask how well iterate averaged 3DVAR behaves if, at fixed $n$, we
choose the optimal $\alpha$, and how this would compare to the minimax
solution of (1.12). Focusing on the diagonal case, for comparison, we have the
following result for the the minimax solution of (1.12):

######  Theorem 1.8.

Under Assumption $2^{\prime}$ with $u_{0}=0$ in (1.18), if (1.12) is solved
using a spectral cutoff with regularization $\alpha$ in the $t\geq 0$ norm,
then at the optimal value of $\alpha=\alpha_{\star}(n)$,

$\mathbb{E}[\|\bar{u}_{n,\alpha_{\star}(n)}-u^{\dagger}\|_{t}^{2}]\lesssim\begin{cases}n^{-\frac{t(1+2\epsilon)+2\beta}{1+2p+2\beta}}&1+2p\neq
t(1+2\epsilon)\\\ n^{-1}\log n&1+2p=t(1+2\epsilon)\end{cases}$

This is consistent with (1.13) and the results in [16, 17]. Then, looking at
the minimax solution of 3DVAR, we obtain for two particular regimes:

######  Corollary 1.9.

With the same assumptions as Theorem 1.7, first, assume $\bar{\tau}_{\rm
b},\bar{\tau}_{\rm v}\leq 1$. Taking $\tau_{\rm b}=\bar{\tau}_{\rm b}$ and
$\tau_{\rm v}=(1-\theta)\bar{\tau}_{\rm v}$ for $\theta\in(0,1]$,

$\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim\theta^{-\frac{2\beta+t(1+2\epsilon)}{1+2p+2\beta+\theta[t(1+2\epsilon)+2\epsilon]}}n^{-\frac{t(1+2\epsilon)+2\beta}{1+2p+2\beta+\theta[t(1+2\epsilon)+2\epsilon]}}.$

If, instead, $\bar{\tau}_{\rm b},\bar{\tau}_{\rm v}>1$, then, taking
$\tau_{\rm b}=\tau_{\rm v}=1$,

$\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim n^{-1}$

Consequently:

  * •

At $t=0$, in the first case,

$\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim\theta^{-\frac{2\beta}{1+2p+2\beta+2\epsilon\theta}}n^{-\frac{2\beta}{1+2p+2\beta+2\epsilon\theta}}$

This is somewhat better than (1.11), as there is no logarithmic term, and the
factor of $2\epsilon$ has been replaced by $2\epsilon\theta$, which can be
reduced by taking $\theta$ smaller. The prefactor will grow, but it is
independent of $n$.

  * •

In the first case, where $\bar{\tau}_{\rm b},\bar{\tau}_{\rm v}\leq 1$, by
taking $\theta$ sufficiently close to zero, we can get arbitrarily close to
the optimal rate in (1.13).

  * •

The first case can be realized by taking $t$ and $\beta$ sufficiently small.
The second case, where $\bar{\tau}_{\rm b},\bar{\tau}_{\rm v}>1$, is
accessible by taking $t$ large enough.

  * •

There are two other cases to consider, $\bar{\tau}_{\rm b}\leq 1$,
$\bar{\tau}_{\rm v}>1$ and vice versa, but, for brevity we do not explore them
here.

In contrast to iterate averaged 3DVAR, there is no gain to iterate averaging
for Kalman:

######  Theorem 1.10.

For the scalar Kalman filter, take $C_{0}=\frac{\gamma^{2}}{\alpha}\sigma>0$.
Then the bias and variance of the iterate-averaged mean, $\bar{m}_{n}$ satisfy
the inequalities

$\displaystyle|\mathbb{E}[\bar{m}_{n}]-u^{\dagger}|$
$\displaystyle\geq|\mathbb{E}[{m}_{n}]-u^{\dagger}|,$
$\displaystyle\operatorname{Var}(\bar{m}_{n})$
$\displaystyle\geq\operatorname{Var}({m}_{n}).$

Consequently, we do not further explore the impact of averaging upon the
Kalman filter in this setting.

###  1.4. Outline

The structure of this paper is as follows. In Section 2 we review certain
background results needed for our main results. Section 3 examines the scalar
case, and it includes proofs of Theorems 1.5 and 1.10. We prove Theorems 1.6
and 1.7 in Section 4. Numerical examples are given in Section 5. We conclude
with a brief discussion in Section 6.

Acknowledgements: The authors thank A.M. Stuart for suggesting an
investigation of this problem. This work was supported by US National Science
Foundation Grant DMS-1818716. The content of this work originally appeared in
[11] as a part of F.G. Jones’s PhD dissertation. Work reported here was run on
hardware supported by Drexel’s University Research Computing Facility.

##  2\. Preliminary Results

In this section, we establish some identities and estimates that will be
crucial to proving our main results.

Much of our analysis relies on spectral calculus involving the following
rational functions which are closely related to the Tikhonov-Phillips
regularization $(\alpha+\lambda)^{-1}$:

(2.1) $\displaystyle r_{n,\alpha}(\lambda)$
$\displaystyle=\left(\frac{\alpha}{\alpha+\lambda}\right)^{n},$ (2.2)
$\displaystyle q_{n,\alpha}(\lambda)$
$\displaystyle=\frac{1}{\lambda}\left\\{1-\left(\frac{\alpha}{\alpha+\lambda}\right)^{n}\right\\}=\lambda^{-1}(1-r_{n,\alpha}(\lambda)).$

These are related by the identity

(2.3) $\sum_{k=1}^{m}r_{k,\alpha}(\lambda)=\alpha q_{m,\alpha}(\lambda).$

The following estimates can be found in [10] and [19], particularly Section
2.2 of the latter reference:

######  Lemma 2.1.

For $\lambda\in[0,\Lambda]$ and $n\in\mathbb{N}$,

$\displaystyle 0<r_{n,\alpha}(\lambda)\leq\frac{\alpha}{\alpha+n\lambda}\leq
1,$
$\displaystyle\lambda^{p}r_{n,\alpha}(\lambda)\leq\begin{cases}\left(\frac{\alpha
p}{n}\right)^{p},&p\in[0,n],\\\ \alpha^{n}\Lambda^{p-n},&p>n.\end{cases}$

######  Lemma 2.2.

For $\lambda\in[0,\Lambda]$, $n\in\mathbb{N}$,

$\displaystyle\lambda^{p}q_{n,\alpha}(\lambda)\leq\begin{cases}\left(\frac{n}{\alpha}\right)^{1-p},&p\in[0,1],\\\
\Lambda^{p-1},&p>1,\end{cases}$
$\displaystyle\lambda^{p}q_{n,\alpha}(\lambda)\leq\lambda^{p-1}.$

Next, we recall the following result on Hilbert scales,

######  Proposition 2.3.

There exists a constant $D>1$, such that for $|\theta|\leq 1$,

$D^{-1}\|x\|_{\theta(1+\nu)}\leq\|(B^{\ast}B)^{\frac{\theta}{2}}x\|\leq
D\|x\|_{\theta(1+\nu)}$

and

$\operatorname{Ran}\left((B^{\ast}B)^{\frac{\theta}{2}}\right)=\operatorname{Dom}\left(\Sigma^{-\frac{\theta(1+\nu)}{2}}_{0}\right).$

This result, based on a duality argument, is proven in Lemma 4.1 of [10]. See,
also, Section 8.4 of [6], particularly Corollary 8.22.

We also have a few useful identities for the filters which we state without
proof.

######  Lemma 2.4.

For the Kalman filter, the mean and covariance operators and the Kalman gains
satisfy the identities

$\displaystyle m_{n}$
$\displaystyle=\left({\gamma^{2}}{n^{-1}}C_{0}^{-1}+A^{*}A\right)^{-1}\left(A^{*}\bar{y}_{n}+{\gamma^{2}}n^{-1}C_{0}^{-1}m_{0}\right)$
$\displaystyle C_{n}^{-1}$
$\displaystyle=C_{n-1}^{-1}+\gamma^{-2}A^{\ast}A=C_{0}^{-1}+\gamma^{-2}nA^{\ast}A$
$\displaystyle K_{n}$
$\displaystyle=(\gamma^{2}C_{n-1}^{-1}+A^{*}A)^{-1}A^{*}=(\gamma^{2}C_{0}^{-1}+nA^{*}A)^{-1}A^{*}=\gamma^{-2}C_{n}A^{*}.$

######  Lemma 2.5.

For 3DVAR,

$\bar{u}_{n}=\sum_{k=0}^{n-1}\frac{n-k}{n}(I-\mathcal{K}A)^{k}\mathcal{K}\bar{y}_{n-k}+\sum_{k=0}^{n-1}\frac{1}{n}(I-\mathcal{K}A)^{k}(I-\mathcal{K}A)u_{0}.$

######  Corollary 2.6.

Letting $v_{n}=u_{n}-u^{\dagger}$,
$\bar{v}_{n}=\frac{1}{n}\sum_{k=1}^{n}v_{k}$,

$\bar{v}_{n}=\sum_{k=0}^{n-1}\frac{n-k}{n}(I-\mathcal{K}A)^{k}\mathcal{K}\bar{\eta}_{n-k}+\sum_{k=0}^{n-1}\frac{1}{n}(I-\mathcal{K}A)^{k}(I-\mathcal{K}A)v_{0}.$

######  Remark 2.7.

As this is a linear problem, it will be sufficient to study the behavior of
$\bar{v}_{n}$ to infer convergence of $\bar{u}_{n}$ to $u^{\dagger}$.

For the analysis of 3DVAR, the essential decomposition into bias and variance
terms can be read off of Corollary 2.6. These can be expressed in the more
useful forms using $q_{n,\alpha}$:

######  Lemma 2.8.

(2.4) $\displaystyle\bar{I}_{n}^{\rm bias}$
$\displaystyle=\sum_{k=0}^{n-1}\frac{1}{n}(I-\mathcal{K}A)^{k}(I-\mathcal{K}A)v_{0}=\frac{\alpha}{n}\Sigma^{\frac{1}{2}}q_{n,\alpha}(B^{\ast}B)\Sigma^{\frac{1}{2}}v_{0},$
(2.5) $\displaystyle\bar{I}_{n}^{\rm var}$
$\displaystyle=\sum_{k=0}^{n-1}\frac{n-k}{n}(I-\mathcal{K}A)^{k}\mathcal{K}\bar{\eta}_{n-k}=\frac{1}{n}\sum_{j=1}^{n}{\Sigma^{\frac{1}{2}}q_{n-j+1,\alpha}(B^{\ast}B)B^{\ast}}\eta_{j}.$

###### Proof.

First, observe that

$I-\mathcal{K}A=\Sigma^{\frac{1}{2}}\alpha(\alpha
I+B^{\ast}B)\Sigma^{-\frac{1}{2}}.$

Using this in (2.4) together with spectral calculus applied to positive self-
adjoint compact operator $B^{\ast}B$, along with (2.3),

$\begin{split}\bar{I}_{n}^{\operatorname{bias}}&=\frac{1}{n}\sum_{k=0}^{n-1}\Sigma^{{1}/{2}}\alpha^{k}(\alpha
I+B^{\ast}B)^{-{k+1}}\Sigma^{-{1}/{2}}_{0}v_{0}\\\
&=\frac{1}{n}\sum_{k=1}^{n}\Sigma^{{1}/{2}}r_{k,\alpha}(B^{\ast}B)\Sigma^{-{1}/{2}}v_{0}=\frac{\alpha}{n}\Sigma^{\frac{1}{2}}q_{n,\alpha}(B^{\ast}B)\Sigma^{-\frac{1}{2}}v_{0}.\end{split}$

Applying the same computations to (2.5), we have,

$\begin{split}\bar{I}_{n}^{\operatorname{var}}&=\sum_{k=0}^{n-1}\frac{n-k}{n}\alpha^{-1}\Sigma^{\frac{1}{2}}_{0}r_{k+1,\alpha}(B^{\ast}B)B^{\ast}\bar{\eta}_{n-k}\\\
&=\frac{1}{n}\sum_{j=1}^{n}\left\\{\sum_{k=0}^{n-j}\alpha^{-1}\Sigma^{\frac{1}{2}}r_{k+1,\alpha}(B^{\ast}B)B^{\ast}\right\\}\eta_{j}=\frac{1}{n}\sum_{j=1}^{n}{\Sigma^{\frac{1}{2}}q_{n-j+1,\alpha}(B^{\ast}B)B^{\ast}}\eta_{j}.\end{split}$

∎

##  3\. Analysis of the Scalar Problem

Before studying the general, infinite-dimensional case, it is instructive to
consider the scalar problem, where $X=Y=\mathbb{R}$ and $A$, $\Sigma$, and
$\mathcal{K}$ are now scalars. This setting will also allow us to establish
the limitations of both 3DVAR and the Kalman filter.

###  3.1. 3DVAR

First, we prove Theorem 1.5 which asserts that the 3DVAR iteration cannot
converge in mean square:

###### Proof.

Since $y_{n}\sim\mathcal{N}(Au^{\dagger},\gamma^{2})$, we write
$y_{n}=Au^{\dagger}+\eta_{n}$ for $\eta_{n}\sim\mathcal{N}(0,\gamma^{2})$. By
(1.4),

$\begin{split}u_{n}-u^{\dagger}&=\mathcal{K}\eta_{n}+\mathcal{K}Au^{\dagger}+(1-\mathcal{K}A)u_{n-1}-u^{\dagger}\\\
&=\mathcal{K}\eta_{n}+(1-\mathcal{K}A)(u_{n-1}-u^{\dagger}).\end{split}$

Consequently,

$\begin{split}\mathbb{E}[|u_{n}-u^{\dagger}|^{2}]&=\mathbb{E}[|\mathcal{K}\eta_{n}|^{2}]+\mathbb{E}[|(1-\mathcal{K}A)(u_{n-1}-u^{\dagger})|^{2}]\\\
&\geq\mathbb{E}[|\mathcal{K}\eta_{n}|^{2}]=\mathcal{K}^{2}\gamma^{2}.\end{split}$

∎

Next, studying the bias and variance of the time averaged problem, given by
(2.4) and (2.5), we prove

######  Theorem 3.1.

For scalar time averaged 3DVAR, for $\tau_{\rm b},\tau_{\rm v}\in[0,1]$

$\mathbb{E}[|\bar{u}_{n}-u^{\dagger}|^{2}]\leq{(A^{2}\Sigma)^{-2\tau_{\rm
b}}}|v_{0}|^{2}\left(\frac{n}{\alpha}\right)^{-2\tau_{\rm
b}}+\frac{\Sigma\gamma^{2}}{\alpha}(A^{2}\Sigma)^{-\tau_{\rm
v}}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}.$

Thus, we have unconditional convergence for any choice for $\alpha>0$,
something that we do not have for 3DVAR without any iterate averaging. The
rate of convergence is greatest when $\tau_{\rm b}\geq 1/2$ and $\tau_{\rm
v}=1$.

To obtain the result, we make use of the bias variance decomposition and
expressions (2.4) and (2.5). In the scalar case, $B^{\ast}B=B^{2}=\Sigma
A^{2}$, so that

(3.1)
$\left|\bar{I}_{n}^{\operatorname{bias}}\right|^{2}=\left(\frac{n}{\alpha}\right)^{-2}q_{n,\alpha}(\Sigma
A^{2})^{2}|v_{0}|^{2}.$

Applying Lemma 2.2 to this expression, we immediately obtain

######  Proposition 3.2.

For $0\leq\tau_{\rm b}\leq 1$,

(3.2)
$\left|\bar{I}_{n}^{\operatorname{bias}}\right|^{2}\leq{(A^{2}\Sigma)^{-2\tau_{\rm
b}}}|v_{0}|^{2}\left(\frac{n}{\alpha}\right)^{-2\tau_{\rm b}}.$

For the variance, we have the result

######  Proposition 3.3.

Let $\tau_{\rm v}\in[0,1]$,

(3.3)
$\mathbb{E}[|\bar{I}_{n}^{\operatorname{var}}|^{2}]\leq\frac{\Sigma\gamma^{2}}{\alpha}(A^{2}\Sigma)^{-\tau_{\rm
v}}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}.$

###### Proof.

For the scalar case of (2.5),using Lemma 2.2,

$\begin{split}\mathbb{E}[|\bar{I}_{n}^{\operatorname{var}}|^{2}]&=\frac{\gamma^{2}(A\Sigma)^{2}}{n^{2}}\sum_{j=1}^{n}q_{j,\alpha}(A^{2}\Sigma)^{2}\\\
&=\frac{\gamma^{2}\Sigma}{n^{2}}(A^{2}\Sigma)^{1-({1+\tau_{\rm
v}})}\sum_{j=1}^{n}\left[(A^{2}\Sigma)^{\frac{1+\tau_{\rm
v}}{2}}q_{j,\alpha}(A^{2}\Sigma)\right]^{2}\\\
&\leq\frac{\Sigma\gamma^{2}}{n^{2}}(A^{2}\Sigma)^{-\tau_{\rm
v}}\sum_{j=1}^{n}\left(\frac{j}{\alpha}\right)^{2\left(1-{\frac{1+\tau_{\rm
v}}{2}}\right)}\\\ &\leq\frac{\Sigma\gamma^{2}(A^{2}\Sigma)^{-\tau_{\rm
v}}}{n^{2}}n\left(\frac{n}{\alpha}\right)^{1-\tau_{\rm
v}}=\frac{\Sigma\gamma^{2}}{\alpha}(A^{2}\Sigma)^{-\tau_{\rm
v}}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}.\end{split}$

∎

###### Proof of Theorem 3.1.

The result then follows immediately by combining the two preceding
propositions.

∎

###  3.2. Kalman Filter

Here, we prove Theorem 1.10, showing there is no improvement in mean squared
convergence of Kalman under iterate averaging.

###### Proof.

Using Lemma 2.4, for the $k$-the estimate of the mean,

$\begin{split}m_{k}&=\left(\frac{\alpha}{\Sigma
k}+a^{2}\right)^{-1}\left(A\bar{y}_{k}+\frac{\alpha}{\Sigma k}m_{0}\right)\\\
&=\left(\frac{\alpha}{\Sigma
k}+A^{2}\right)^{-1}\left(A^{2}u^{\dagger}+A\bar{\eta}_{k}+\frac{\alpha}{\Sigma
k}m_{0}\right)\\\ &=\left(1+\frac{\alpha}{A^{2}\Sigma
k}\right)^{-1}u^{\dagger}+\left(1+\frac{A^{2}\Sigma
k}{\alpha}\right)^{-1}m_{0}+\left(A+\frac{\alpha}{A\Sigma
k}\right)^{-1}\bar{\eta}_{k}.\end{split}$

and without averaging,

$\displaystyle\mathbb{E}[m_{n}]-u^{\dagger}$
$\displaystyle=\left(1+\frac{A^{2}\Sigma
n}{\alpha}\right)^{-1}(m_{0}-u^{\dagger}),$
$\displaystyle\operatorname{Var}(m_{n})$
$\displaystyle=\left({A+\frac{\alpha}{A\Sigma
n}}\right)^{-2}\frac{\gamma^{2}}{n}.$

Then, with averaging, for the bias,

$\mathbb{E}[\bar{m}_{n}]-u^{\dagger}=\frac{1}{n}\sum_{k=1}^{n}\left(1+\frac{A^{2}\Sigma
k}{\alpha}\right)^{-1}(m_{0}-u^{\dagger}),$

and

$\begin{split}|\mathbb{E}[\bar{m}_{n}]-u^{\dagger}|^{2}&=\left|\frac{1}{n}\sum_{k=1}^{n}\left(1+\frac{A^{2}\Sigma
k}{\alpha}\right)^{-1}\right|^{2}|m_{0}-u^{\dagger}|^{2}\\\
&\geq\left|\frac{1}{n}\sum_{k=1}^{n}\left(1+\frac{A^{2}\Sigma
n}{\alpha}\right)^{-1}\right|^{2}|m_{0}-u^{\dagger}|^{2}=|\mathbb{E}[m_{n}]-u^{\dagger}|^{2}.\end{split}$

For the variance, first note

$\begin{split}\bar{m}_{n}-\mathbb{E}[\bar{m}_{n}]=\frac{1}{n}\sum_{k=1}^{n}\left(A+\frac{\alpha}{A\Sigma
k}\right)^{-1}\bar{\eta}_{k}&=\frac{1}{n}\sum_{k=1}^{n}\left(A+\frac{\alpha}{A\Sigma
k}\right)^{-1}\left\\{\sum_{j=1}^{k}\eta_{j}\right\\}\\\
&=\frac{1}{n}\sum_{j=1}^{n}\eta_{j}\left\\{\sum_{k=j}^{n}\left(A+\frac{\alpha}{A\Sigma
k}\right)^{-1}\right\\}.\end{split}$

Then, by dropping all but the $k=n$-th term in the inner sum,

$\begin{split}\operatorname{Var}(\bar{m}_{n})=\frac{1}{n^{2}}\sum_{j=1}^{n}\gamma^{2}\left\\{\sum_{k=j}^{n}\left(A+\frac{\alpha}{A\Sigma
k}\right)^{-1}\right\\}^{2}&\geq\frac{1}{n^{2}}\sum_{j=1}^{n}\gamma^{2}\left(A+\frac{\alpha}{A\Sigma
n}\right)^{-2}\\\ &\quad=\operatorname{Var}(m_{n})\end{split}$

∎

##  4\. Analysis of the Infinite Dimensional Problem

We return to the bias and variance of 3DVAR in the general, potentially
infinite dimensional, setting and obtain estimates on the terms. We prove the
general case in Section 4.1, and then the diagonal case in 4.2. Our minimax
results are proven in Section 4.3.

###  4.1. General Case

Here, we prove Theorem 1.6 by first establishing results on the bias and
variance.

######  Proposition 4.1.

Under Assumption 1, with $t\in[0,\nu]$,

$\|\bar{I}_{n}^{\operatorname{bias}}\|_{t}^{2}\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{s+t}{1+\nu}}\|z_{0}\|^{2}$

where $z_{0}$ solves (1.15).

The fastest possible decay available for the squared bias in Proposition 4.1
is $\operatorname{O}(n^{-2})$ when $s=\nu+2$ and $t=\nu$.

###### Proof.

We make use of bias term from Lemma 2.8, allowing us to write

$\|\bar{I}_{n}^{\operatorname{bias}}\|_{t}^{2}=\left\|\frac{\alpha}{n}\Sigma^{\frac{t+1}{2}}q_{n,\alpha}(B^{\ast}B)\Sigma^{-\frac{1}{2}}v_{0}\right\|^{2}.$

Next, we make use of (1.6) and argue as in the Appendix of [10], applying
Proposition 2.3. Since, by assumption,
$v_{0}\in\operatorname{Dom}(\Sigma^{-\frac{s}{2}})$,
$\Sigma^{-\frac{1}{2}}v_{0}\in\operatorname{Dom}(\Sigma^{-\frac{s-1}{2}})$.
Then taking $\theta=(s-1)/(1+\nu)$ in the proposition,
$\Sigma^{-\frac{1}{2}}v_{0}\in\operatorname{Ran}((B^{\ast}B)^{\frac{s-1}{2(1+\nu)}})$
allows us to conclude the existence of $z_{0}$. Therefore,

$\|\bar{I}_{n}^{\operatorname{bias}}\|_{t}^{2}=\left\|\frac{\alpha}{n}\Sigma^{\frac{t+1}{2}}q_{n,\alpha}(B^{\ast}B)(B^{\ast}B)^{\frac{s-1}{2(1+\nu)}}z_{0}\right\|^{2}.$

Next, using Proposition 2.3 again, now with $\theta=(1+t)/(1+\nu)$,

$\begin{split}\|\bar{I}_{n}^{\operatorname{bias}}\|_{t}^{2}&\lesssim\left\|\frac{\alpha}{n}(B^{\ast}B)^{\frac{t+1}{2(1+\nu)}}q_{n,\alpha}(B^{\ast}B)(B^{\ast}B)^{\frac{s-1}{2(1+\nu)}}z_{0}\right\|^{2}\\\
&\quad=\left\|\frac{\alpha}{n}(B^{\ast}B)^{\frac{s+t}{2(1+\nu)}}q_{n,\alpha}(B^{\ast}B)z_{0}\right\|^{2}\\\
&\leq\left(\sup_{0\leq\lambda\leq\|B^{\ast}B\|}\left|\frac{\alpha}{n}\lambda^{\frac{s+t}{2(1+\nu)}}q_{n,\alpha}(\lambda)\right|\right)^{2}\|z_{0}\|^{2}\leq\left(\frac{n}{\alpha}\right)^{-\frac{s+t}{1+\nu}}\|z_{0}\|^{2}.\end{split}$

The last inequality holds since, $s\leq\nu+2$ and $t\leq\nu$, so that $0\leq
s+t\leq s+\nu\leq 2\nu+2$ allowing for the application of Lemma 2.2. ∎

######  Proposition 4.2.

Under Assumption 1, for $t\geq 0$, $\tau_{\rm v}\in[0,1]$, and for this choice
of $\tau_{\rm v}$ and $t$, assume $\Sigma^{(1+t)-\tau_{\rm v}(1+\nu)}$ is
trace class. Then

$\mathbb{E}[\|\bar{I}_{n}^{\operatorname{var}}\|_{t}^{2}]\lesssim\frac{\gamma^{2}}{\alpha}\operatorname{Tr}(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)})\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}.$

######  Remark 4.3.

The fastest possible decay in the variance will be $\operatorname{O}(n^{-1})$
when $\tau_{\rm v}=1$ and $t$ is sufficiently large such that $\Sigma^{t-\nu}$
is trace class. However, the bias term requires $t\leq\nu$. This requires the
identity operator to be trace class which will not hold in infinite
dimensions.

###### Proof of Proposition 4.2.

We begin with equation (2.5) and using that for any bounded operator $T$ and
positive self adjoint trace class operator $C$,
$|\operatorname{Tr}(CT)|\leq\|T\||\operatorname{Tr}C|$,

$\begin{split}\mathbb{E}[\|\bar{I}_{n}^{\operatorname{var}}\|_{t}^{2}]&=\frac{1}{n^{2}}\sum_{j=1}^{n}\mathbb{E}[\|{\Sigma^{\frac{t+1}{2}}q_{n-j+1,\alpha}(B^{\ast}B)B^{\ast}}\eta_{j}\|^{2}]\\\
&=\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\operatorname{Tr}\left(\Sigma^{\frac{t+1}{2}}q_{j,\alpha}(B^{\ast}B)(B^{\ast}B)q_{j,\alpha}(B^{\ast}B)\Sigma^{\frac{t+1}{2}}\right)\\\
&=\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\operatorname{Tr}\left(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)}\left(\Sigma^{\tau_{\rm
v}\frac{1+\nu}{2}}(B^{\ast}B)^{\frac{1}{2}}q_{j,\alpha}(B^{\ast}B)(B^{\ast}B)\right)^{2}\right)\\\
&\leq\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\|{\Sigma^{\tau_{\rm
v}\frac{1+\nu}{2}}(B^{\ast}B)^{\frac{1}{2}}q_{j,\alpha}(B^{\ast}B)(B^{\ast}B)}\|^{2}\operatorname{Tr}(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)}).\end{split}$

Using Proposition 2.3 with $\theta=\tau_{\rm v}$ and Lemma 2.2,

$\begin{split}\|{\Sigma^{\tau_{\rm
v}\frac{1+\nu}{2}}(B^{\ast}B)^{\frac{1}{2}}q_{j,\alpha}(B^{\ast}B)(B^{\ast}B)}\|&\lesssim\|(B^{\ast}B)^{\frac{1+\tau_{\rm
v}}{2}}q_{j,\alpha}(B^{\ast}B)\|\\\
&\lesssim\sup_{\lambda\in[0,\|B^{\ast}B\|]}\lambda^{\frac{1+\tau_{\rm
v}}{2}}q_{j,\alpha}(\lambda)\lesssim\left(\frac{j}{\alpha}\right)^{1-\frac{1+\tau_{\rm
v}}{2}}\end{split}$

Therefore,

$\begin{split}\mathbb{E}[\|\bar{I}_{n}^{\operatorname{var}}\|_{t}^{2}]&\lesssim\frac{\gamma^{2}}{n^{2}}\operatorname{Tr}(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)})\sum_{j=1}^{n}\left(\frac{j}{\alpha}\right)^{1-\tau_{\rm
v}}\lesssim\frac{\gamma^{2}}{\alpha}\operatorname{Tr}(\Sigma^{t+1-\tau_{\rm
v}(1+\nu)})\left(\frac{n}{\alpha}\right)^{-\tau_{\rm v}}\end{split}$

∎

###### Proof of Theorem 1.6.

The theorem immediately follows from the two preceding propositions. ∎

###  4.2. Simultaneous Diagonalization

A sharper result is available under the simultaneous diagonalization
Assumption $2^{\prime}$ . For convenience, letting

(4.1) $\omega=\frac{1+2\epsilon}{1+2\epsilon+2p},$

we have the relationship

(4.2) $\sigma_{i}\asymp(\sigma_{i}a_{i}^{2})^{\omega}.$

######  Proposition 4.4.

Under Assumption $2^{\prime}$, let $\tau_{\rm b}\in[0,1]$ satisfy condition
(1.19a),

$\left\|\bar{I}_{n}^{\operatorname{bias}}\right\|_{t}^{2}\lesssim\left(\frac{n}{\alpha}\right)^{-2\tau_{\rm
b}}.$

###### Proof.

We start with equation (2.4) and then use (4.2) and Lemma 2.2,

$\begin{split}\left\|\bar{I}_{n}^{\operatorname{bias}}\right\|_{t}^{2}&=\sum_{i=1}^{\infty}\left\langle\frac{\alpha}{n}\Sigma^{\frac{t+1}{2}}q_{n,\alpha}(B^{\ast}B)\Sigma^{-\frac{1}{2}}v_{0},\varphi_{i}\right\rangle^{2}\\\
&=\sum_{i=1}^{\infty}\left|\frac{\alpha}{n}\sigma_{i}^{\frac{t}{2}}q_{n,\alpha}(\sigma_{i}a_{i}^{2})\right|^{2}\left|v_{0,i}\right|^{2}=\left(\frac{\alpha}{n}\right)^{2}\sum_{i=1}^{\infty}\sigma_{i}^{t}q_{n,\alpha}(\sigma_{i}a_{i}^{2})^{2}\left|v_{0,i}\right|^{2}\\\
&\asymp\left(\frac{n}{\alpha}\right)^{-2}\sum_{i=1}^{\infty}(\sigma_{i}a_{i}^{2})^{t\omega-2\tau_{\rm
b}}((\sigma_{i}a_{i}^{2})^{\tau_{\rm
b}}q_{n,\alpha}(\sigma_{i}a_{i}^{2}))^{2}\left|v_{0,i}\right|^{2}\\\
&\lesssim\left(\frac{n}{\alpha}\right)^{-2}\left(\frac{n}{\alpha}\right)^{2-2\tau_{\rm
b}}\sum_{i=1}^{\infty}(\sigma_{i}a_{i}^{2})^{t\omega-2\tau_{\rm
b}}\left|v_{0,i}\right|^{2}\end{split}$

Using (1.19a),

$\begin{split}\sum_{i=1}^{\infty}(\sigma_{i}a_{i}^{2})^{t\omega-2\tau_{\rm
b}}\left|v_{0,i}\right|^{2}&\asymp\sum_{i=1}^{\infty}i^{-(1+2\epsilon+2p)(t\omega-i\tau_{\rm
b})}\left|v_{0,i}\right|^{2}\\\
&\asymp\sum_{i=1}^{\infty}i^{-(1+2\epsilon+2p)(t\omega-2\tau_{\rm
b})-2\beta}i^{2\beta}\left|v_{0,i}\right|^{2}\\\
&\lesssim\left(\sup_{i}i^{-(1+2\epsilon+2p)(t\omega-2\tau_{\rm
b})-2\beta}\right)\sum_{i=1}^{\infty}i^{2\beta}\left|v_{0,i}\right|^{2}<\infty\end{split}$

we have the result.

∎

Comparing this to the general case, we again see that if the data is
sufficiently smooth and/or we study the probelm in a sufficiently smooth space
($\beta$ and/or $t$ large), we can again obtain $\operatorname{O}(n^{-2})$
convergence of the squared bias.

######  Proposition 4.5.

Under Assumption $2^{\prime}$, having fixed $t$, for $\tau_{\rm v}\in[0,1]$
satisfying (1.19b),

$\mathbb{E}\left[\left\|\bar{I}_{n}^{\operatorname{var}}\right\|_{t}^{2}\right]\lesssim\frac{\gamma^{2}}{\alpha}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm
v}}$

###### Proof.

Using (2.5), we begin by writing

$\begin{split}\mathbb{E}\left[\left\|\bar{I}_{n}^{\operatorname{var}}\right\|_{t}^{2}\right]&=\frac{1}{n^{2}}\sum_{j=1}^{n}\mathbb{E}\left[\left\|\Sigma^{\frac{t+1}{2}}q_{n-j+1,\alpha}(B^{\ast}B)B^{\ast}\eta_{j}\right\|^{2}\right],\\\
&=\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\operatorname{Tr}\left(\Sigma^{\frac{t+1}{2}}q_{n-j+1,\alpha}B^{\ast}B(B^{\ast}B)q_{n-j+1,\alpha}(B^{\ast}B)\Sigma^{\frac{t+1}{2}}\right),\\\
&=\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\operatorname{Tr}\left(\Sigma^{t+1}B^{\ast}Bq_{j,\alpha}(B^{\ast}B)^{2}\right).\end{split}$

Using (2.2) on each term in the sum,

$\begin{split}\operatorname{Tr}\left(\Sigma^{t+1}B^{\ast}Bq_{j,\alpha}(B^{\ast}B)^{2}\right)&=\sum_{i=1}^{\infty}\sigma_{i}^{t+2}a_{i}^{2}q_{j,\alpha}(\sigma_{i}a_{i}^{2})^{2}.\end{split}$

Then, using (4.2) and Lemma 2.2

$\begin{split}\sigma_{i}^{t+2}a_{i}^{2}q_{j,\alpha}(\sigma_{i}a_{i}^{2})^{2}&\asymp\sigma_{i}^{t+1}((\sigma_{i}a_{i}^{2})^{\frac{1}{2}}q_{j,\alpha}(\sigma_{i}a_{i}^{2}))^{2}\\\
&\asymp(\sigma_{i}a_{i}^{2})^{\omega(t+1)}((\sigma_{i}a_{i}^{2})^{\frac{1}{2}}q_{j,\alpha}(\sigma_{i}a_{i}^{2}))^{2}\\\
&\asymp(\sigma_{i}a_{i}^{2})^{\omega(t+1)-\tau_{\rm
v}}((\sigma_{i}a_{i}^{2})^{(1+\tau_{\rm
v})/2}q_{j,\alpha}(\sigma_{i}a_{i}^{2}))^{2}\\\
&\lesssim(\sigma_{i}a_{i}^{2})^{\omega(t+1)-\tau_{\rm
v}}\left(\frac{j}{\alpha}\right)^{1-\tau_{\rm v}}\end{split}$

Under assumption 1.19b

(4.3) $\sum_{i=1}^{\infty}(\sigma_{i}a_{i}^{2})^{\omega(t+1)-\tau_{\rm
v}}\asymp\sum_{i=1}^{\infty}i^{-[(1+2\epsilon)(t+1)-\tau_{\rm
v}(1+2\epsilon+2p)]}<\infty$

Consequently,

$\operatorname{Tr}\left(\Sigma^{t+1}(B^{\ast}B)q_{j,\alpha}(B^{\ast}B)^{2}\right)\lesssim\left(\frac{j}{\alpha}\right)^{1-\tau_{\rm
v}},$

and

$\frac{\gamma^{2}}{n^{2}}\sum_{j=1}^{n}\operatorname{Tr}\left(\Sigma^{t+1}(B^{\ast}B)q_{j,\alpha}(B^{\ast}B)^{2}\right)\lesssim\frac{\gamma^{2}}{\alpha}\left(\frac{n}{\alpha}\right)^{-\tau_{\rm
v}}$

∎

In contrast to the non-diagonal case, if the problem is studied in a
sufficiently weak sense (large enough $t$), one obtains
$\operatorname{O}(n^{-1})$ convergence of the variance.

###### Proof of Theorem 1.7.

This result immediately follows from the previous two propositions.

∎

###  4.3. Minimax Analysis

###### Proof of Theorem 1.8.

Recall the spectral cutoff regularization

$g_{\alpha}(\lambda)=\lambda^{-1}1_{[\alpha,\infty)}(\lambda).$

For a fixed $\alpha$, the regularized solution of (1.12) is

$\bar{u}_{n,\alpha}=g_{\alpha}(A^{\ast}A)A^{\ast}\bar{y}_{n}=g_{\alpha}(A^{\ast}A)A^{\ast}Au^{\dagger}+g_{\alpha}(A^{\ast}A)A^{\ast}\bar{\eta}_{n}.$

This allows us to write the bias-variance decomposition of the error as

$\mathbb{E}[\|\bar{u}_{n,\alpha}-u^{\dagger}\|_{t}^{2}]=\|(g_{\alpha}(A^{\ast}A)A^{\ast}A-I)u^{\dagger}\|_{t}^{2}+\mathbb{E}[\|g_{\alpha}(A^{\ast}A)A^{\ast}\bar{\eta}_{n}\|_{t}^{2}]$

For the bias term,

$\begin{split}\|(g_{\alpha}(A^{\ast}A)A^{\ast}A-I)u^{\dagger}\|_{t}^{2}&=\sum_{i=1}^{\infty}\sigma^{t}_{i}(g_{\alpha}(a_{i}^{2})a_{i}^{2}-1)^{2}|u_{i}^{\dagger}|^{2}.\end{split}$

Since $a_{i}^{2}\asymp i^{-2p}$, all terms with
$i\lesssim\alpha^{-\frac{1}{2p}}$ will vanish. This leaves us with

$\begin{split}\|(g_{\alpha}(A^{\ast}A)A^{\ast}A-I)u^{\dagger}\|_{t}^{2}&\lesssim\sum_{i=\lfloor\alpha^{-\frac{1}{2p}}\rfloor}^{\infty}\sigma_{i}^{t}|u_{i}^{\dagger}|^{2}\\\
&\lesssim\sum_{i=\lfloor\alpha^{-\frac{1}{2p}}\rfloor}^{\infty}i^{-t(1+2\epsilon)-2\beta}i^{2\beta}|u_{i}^{\dagger}|^{2}\\\
&\lesssim\alpha^{\frac{t(1+2\epsilon)+2\beta}{2p}}\sum_{i}i^{2\beta}|u_{i}^{\dagger}|^{2}.\end{split}$

For the variance term, all terms with $i\gtrsim\alpha^{-\frac{1}{2p}}$ will
vanish,

$\begin{split}\mathbb{E}[\|g_{\alpha}(A^{\ast}A)A^{\ast}\eta\|_{t}^{2}]&=\frac{\gamma^{2}}{n}\sum_{i=1}^{\infty}\sigma_{i}^{t}g_{\alpha}(a_{i}^{2})^{2}a_{i}^{2}\\\
&\lesssim\frac{\gamma^{2}}{n}\sum_{i=1}^{\lceil\alpha^{-\frac{1}{2p}}\rceil}\sigma_{i}^{t}a_{i}^{-2}\lesssim\frac{\gamma^{2}}{n}\sum_{i=1}^{\lceil\alpha^{-\frac{1}{2p}}\rceil}i^{2p-t(1+2\epsilon)}\\\
&\lesssim\begin{cases}\frac{\gamma^{2}}{n}\alpha^{-\frac{1+2p-t(1+2\epsilon)}{2p}}&2p-t(1+2\epsilon)\neq-1\\\
\frac{\gamma^{2}}{n}\log\frac{1}{\alpha}&2p-t(1+2\epsilon)=-1.\end{cases}\end{split}$

Combining the two terms, we thus have,

$\mathbb{E}[\|u-u^{\dagger}\|_{t}^{2}]\lesssim\begin{cases}\alpha^{\frac{t(1+2\epsilon)+2\beta}{2p}}+\frac{\gamma^{2}}{n}\alpha^{-\frac{1+2p-t(1+2\epsilon)}{2p}}&2p-t(1+2\epsilon)\neq-1\\\
\alpha^{\frac{t(1+2\epsilon)+2\beta}{2p}}+\frac{\gamma^{2}}{n}\log\frac{1}{\alpha}&2p-t(1+2\epsilon)=-1.\end{cases}$

Optimizing over $\alpha$ yields the result. ∎

###### Proof of Corollary 1.9.

Note that in the proof of Proposition 4.5, under our assumptions, in (4.3),
with $\tau_{\rm v}=(1-\theta)\bar{\tau}_{\rm v}$

$\begin{split}\sum_{i=1}^{\infty}(\sigma_{i}a_{i}^{2})^{\omega(t+1)-\tau_{\rm
v}}&\asymp\sum_{i=1}^{\infty}i^{-[(1+2\epsilon)(t+1)-t(1+2\epsilon)(1-\theta)-2\epsilon(1-\theta)]}\\\
&\asymp\sum_{i=1}^{\infty}i^{-1-\theta[t(1+2\epsilon)+2\epsilon]}\\\
&\lesssim\int_{1}^{\infty}x^{-1-\theta[t(1+2\epsilon)+2\epsilon]}dx=\frac{1}{\theta[t(1+2\epsilon)+2\epsilon]}\end{split}$

Consequently,

$\mathbb{E}[\|\bar{u}_{n}-u^{\dagger}\|_{t}^{2}]\lesssim\left(\frac{n}{\alpha}\right)^{-\frac{t(1+2\epsilon)+2\beta}{1+2\epsilon+2p}}+\frac{1}{\theta}\frac{1}{\alpha}\left(\frac{n}{\alpha}\right)^{-(1-\theta)\frac{t(1+2\epsilon)+2\epsilon}{1+2\epsilon+2p}}$

where the implicit constants in each term are independent of $\alpha$,
$\theta$, and $n$. The optimally scaled $\alpha$ will be

$\alpha=\theta^{-\frac{1+2p+2\epsilon}{1+2p+2\beta+\theta[t(1+2\epsilon)+2\epsilon]}}n^{\frac{2\beta-2\epsilon+\theta[t(1+2\epsilon)+2\epsilon]}{1+2p+2\beta+\theta[t(1+2\epsilon)+2\epsilon]}}$

Substituting back in, we have our result. ∎

##  5\. Numerical Experiments

In this section we illustrate our results with some numerical experiments.

###  5.1. Scalar Examples

As a simple scalar example, let $A=1$, $\gamma=0.1$, and $u^{\dagger}=0.5$.
For 3DVAR, take $u_{0}=0$, $\Sigma=1$, and $\alpha=1$, while for Kalman, take
$m_{0}=0$ and $C_{0}=1$. Running $10^{2}$ independent trials of each algorithm
for $10^{4}$ iterations, we obtain the results in Figure 1. These simulations
demonstrate our predictions from Theorems 1.5, Theorem 3.1, and Theorem 1.10,
that 3DVAR can only converge with time averaging, while Kalman will not be
improved by time averaging. The confidence bounds are computed using $10^{4}$
bootstrap samples to produce 95% confidence intervals.

Figure 1. Scalar results for 3DVAR and the Kalman filter. These results are
consistent with Theorems 1.5,1.10, and 3.1; 3DVAR will not converge without
time averaging while Kalman will not improve from time averaging. Shaded
regions reflect 95% confidence intervals at each $n$.

###  5.2. Simultaneous Diagonalization Example

Next, we consdier the case of simultaneous diagonalization, working with
functions in $L^{2}(0,2\pi;\mathbb{R})$, and

(5.1) $A=(I-\tfrac{d^{2}}{dx^{2}})^{-1},\quad\Sigma=A^{2},\quad
u^{\dagger}=0.$

The $A$ operator is equipped with periodic boundary conditions, allowing us to
easily work in Fourier space. As the problem is linear, we can separately
consider the bias and the variance. In all examples below we discretize on
$N=2^{12}$ modes, and run for $10^{4}$ iterations. This corresponds to $p=2$
and $\epsilon=1.5$ in Assumption $2^{\prime}$.

For the bias, we choose, before truncation, as the initial condition

(5.2) $u_{0}=\sum_{k=1}^{\infty}k^{-\frac{1}{2}-\beta-\delta}\cos(kx),$

with $\beta=1$ and $\delta=0.01$. Consequently, this function satisfies (1.18)
from Assumption $2^{\prime}$. The perturbation $\delta$ is introduced so that
we can best see the sharpness of our rates. Running the truncated and
discretized problem, we obtain the results shown in Figure 2 for the norms
$t=0,0.5,1,2$. As the plots show, we are in good agreement with the maximal
rate predicted by Theorem 1.7.

For the variance, taking $u_{0}=0$, we run $10^{2}$ independent trials of the
problem, and then use bootstrapping to estimate 95% confidence intervals. The
results, shown in Figure 3, again show good agreement with the maximal rate
predicted by Theorem 1.7.

Figure 2. Decay of the squared bias in our simultaneously diagonalized test
problem for different $t$-norms. All are in good agreement with the rates
predicted by Theorem 1.7. The constant $\bar{\tau}_{\rm v}$ reflects the
greatest possible decay rate from (1.19a).

Figure 3. Decay of the mean squared variance term in our simultaneously
diagonalized test problem for different $t$-norms. All are in good agreement
with the rates predicted by Theorem 1.7. Shaded regions reflect 95% confidence
intervals at each $n$. The constant $\bar{\tau}_{\rm v}$ reflects the greatest
possible decay rate from (1.19b).

##  6\. Discussion

In this work we have examined the impact of iterate averaging upon the Kalman
filter and 3DVAR as tools for solving a statistical inverse problem. We have
found that this modest post-processing step ensures that the simpler
algorithm, 3DVAR, will converge, unconditionally with respect to $\alpha$, in
mean square as the number of iterations $n\to\infty$. In contrast, there is no
performance gain when this averaging is applied to the Kalman filter.

Our simulations suggest that our rates, at least in the diagonal case, may be
sharp. For the diagonal case, we should expect to see something slower than
the Monte Carlo rate of convergence, $\operatorname{O}(n^{-1})$ unless working
in a sufficiently weak space (large $t$). In the general case, it would seem
that for the infinite dimensional problem, we will never be able to achieve
$\operatorname{O}(n^{-1})$ convergence for the reasons outlined in Remark 4.3;
the operator $\Sigma^{t-\nu}$ would need to be trace class, but $t\leq\nu$ for
the bias to converge. The sharpness of the result in the non-diagonalizable
case remains to be established. There is also potential for the extension of
this work to the analogous continuous in time problem (1.5) studied in [15].

In actual applications, the problem will always be finite dimensional, making
$\operatorname{O}(n^{-1})$ achievable. In a spectral Galerkin formulation,
truncating to $N$ modes, and, $\operatorname{Tr}\Sigma^{t-\nu}_{N}$, will
always be finite, though the constant may be large. Hence, we should expect to
see $\operatorname{O}(n^{-1})$ convergence, for sufficiently large $n$ and a
sufficiently severe dimensional truncation.

## Acknowledgements

The authors thank A.M. Stuart for suggesting an investigation of this problem.
This work was supported by US National Science Foundation Grant DMS-1818716.
The content of this work originally appeared in [11] as a part of F.G. Jones’s
PhD dissertation. Work reported here was run on hardware supported by Drexel’s
University Research Computing Facility.

## References

  * [1] L. Cavalier.  Nonparametric statistical inverse problems.  Inverse Problems, 24(3):034004, 2008.
  * [2] G. Da Prato and J. Zabczyk.  Stochastic equations in infinite dimensions.  Cambridge university press, 2014.
  * [3] L. Ding, S. Lu, and J. Cheng.  Weak-norm posterior contraction rate of the 4dvar method for linear severely ill-posed problems.  Journal of Complexity, 46:1–18, 2018.
  * [4] R. Ghanem, D. Higdon, and H. Owhadi, editors.  Handbook of Uncertainty Quantification.  Springer International Publishing, Cham, 2017.
  * [5] S. Ghosal and A. van der Vaart.  Fundamentals of Nonparametric Bayesian Inference.  Cambridge University Press, 2017.
  * [6] W. E. Heinz, M. Hanke, and A. Neubauer.  Regularization of Inverse Problems.  Kluwer Academic Publishers, 2000.
  * [7] J. Humpherys, P. Redd, and J. West.  A Fresh Look at the Kalman Filter.  SIAM Review, 54(4):801–823, Jan. 2012.
  * [8] M. A. Iglesias.  A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems.  Inverse Problems, 32:025002, 2016.
  * [9] M. A. Iglesias, K. J. Law, and A. M. Stuart.  Ensemble kalman methods for inverse problems.  Inverse Problems, 29(4):045001, 2013.
  * [10] M. A. Iglesias, K. Lin, S. Lu, and A. M. Stuart.  Filter based methods for statistical linear inverse problems.  Communications in Mathematical Sciences, 15(7):1867–1896, 2017\.
  * [11] F. G. E. Jones.  High and Infinite-Dimensional Filtering Methods.  PhD thesis, Drexel University, 2020.
  * [12] B. T. Knapik, A. W. van der Vaart, J. H. van Zanten, et al.  Bayesian inverse problems with gaussian priors.  The Annals of Statistics, 39(5):2626–2657, 2011.
  * [13] H. J. Kushner and G. Yin.  Stochastic approximation and recursive algorithms and applications.  Springer, New York, 2003.
  * [14] K. Law, A. Stuart, and K. Zygalakis.  Data Assimilation: A Mathematical Introduction.  Springer International Publishing, 2015.
  * [15] S. Lu, P. Niu, and F. Werner.  On the asymptotical regularization for linear inverse problems in presence of white noise.  SIAM/ASA Journal on Uncertainty Quantification, 9(1):1–28, 2021\.
  * [16] B. A. Mair and F. H. Ruymgaart.  Statistical inverse estimation in hilbert scales.  SIAM Journal on Applied Mathematics, 56(5):1424–1444, 1996.
  * [17] P. Mathé and S. V. Pereverzev.  Optimal Discretization of Inverse Problems in Hilbert Scales. Regularization and Self-Regularization of Projection Methods.  SIAM Journal on Numerical Analysis, 38(6):1999–2021, Jan. 2001\.
  * [18] B. Øksendal.  Stochastic differential equations.  Springer, 2003.
  * [19] S. Pereverzev and S. Lu.  Regularization Theory for Ill-Posed Problems.  De Gruyter, 2013.
  * [20] B. T. Polyak and A. B. Juditsky.  Acceleration of stochastic approximation by averaging.  SIAM journal on control and optimization, 30(4):838–855, 1992.
  * [21] C. Schillings and A. Stuart.  Convergence analysis of ensemble Kalman inversion: the linear, noisy case.  Applicable Analysis, pages 1–17, 2017.
  * [22] C. Schillings and A. M. Stuart.  Analysis of the Ensemble Kalman Filter for Inverse Problems.  SIAM Journal on Numerical Analysis, 55:1264–1290, 2017.
  * [23] R. H. Shumway and D. S. Stoffer.  Time Series Analysis and Its Applications.  Springer, 2011.
  * [24] A. M. Stuart.  Inverse problems: A Bayesian perspective.  Acta Numerica, 19:451–559, May 2010.
  * [25] T. J. Sullivan.  Introduction to uncertainty quantification, volume 63.  Springer, 2015.
  * [26] A. C. van Rooij and F. H. Ruymgaart.  Asymptotic minimax rates for abstract linear estimators.  Journal of Statistical Planning and Inference, 53(3):389–402, Aug. 1996.