Incremental Reading for Question Answering

Any system which performs goal-directed continual learning must not only learn incrementally but process and absorb information incrementally. Such a system also has to understand when its goals have been achieved. In this paper, we consider these issues in the context of question answering. Current state-of-the-art question answering models reason over an entire passage, not incrementally. As we will show, naive approaches to incremental reading, such as restriction to unidirectional language models in the model, perform poorly. We present extensions to the DocQA [2] model to allow incremental reading without loss of accuracy. The model also jointly learns to provide the best answer given the text that is seen so far and predict whether this best-so-far answer is sufficient.

Introduction
Humans can read and comprehend text incrementally. For instance, given a piece of text, our mental state gets updated as we read [10]. We do not necessarily wait until the end of a long document to understand its first sentence. This incremental reading mechanism allows us to avoid consuming more input if we have already reached our goal, and is analogous to the problem faced by a goal-directed continuous learning system, which must also incrementally absorb new information, determine how to use it, and determine if its goals have been achieved. Inspired by how humans read and learn from text incrementally, we introduce incremental models for text comprehension. Our primary goal is to address the problem of incremental reading in the context of language comprehension and Question Answering (QA). We formulate the problem as designing a model for question-answering that consumes text incrementally.
By design and definition, Recurrent Neural Networks (RNNs), e.g. Long Short-Term Memory networks (LSTMs) [5], process data sequentially, and update their internal states as they read the new tokens. However, on tasks like Question Answering, in all the existing well-performing models, RNNs are employed in a bidirectional way, or a self-attention mechanism is employed [11,2,4,8]. This means these models need to processes the whole input sequence to compute the final answer. This is a reasonable approach if the input sequence is as short as a sentence, but it becomes less effective and efficient as the length of the input sequence increases.
We introduce a new incremental model based on DocQA [2], which is an RNN based model proposed for QA. The incremental DocQA performs similarly to the original system but can process the input text incrementally. We propose the use of slicing to build incremental models. Slicing RNNs were introduced in [13] with the motivation of enabling parallelization and speeding up sequence processing. Here, we explore using slicing to facilitate incremental processing of the input sequence. Besides the fact that incremental reading is more plausible from the cognitive perspective, it can also provide an inductive bias for the models, which will make it easier for them to find a more generalizable solution [1]. Moreover, incrementality allows the model to be applied to tasks where we do not have the whole input in the beginning, and we need to compute some intermediate outputs. E.g. when the input is a stream, or we are in the context of a dialogue. We observed that, even if the whole input text is available, it is not always necessary to read the whole text to answer a given question. In fact, our model achieves its highest performance when it is limited to read only a few tokens after the answer, rather than when it is allowed access to the entire context. Considering this, we also augment the incremental DocQA with an "early stopping" strategy, where it stops consuming the rest of the input text as soon as it is confident that it has found the answer. Learning to stop reading or reasoning when the goal is reached is addressed for tasks like text classification [3] and question answering [9,12,6], but the main challenge is that implementing early stopping strategies is only possible when we have an incremental model.

Architecture of the models
We use DocQA as the baseline model [2]. The architecture of DocQA is illustrated in Figure 4a. The output of this model is two vectors with the length of the given context. One of the vectors indicates the probability of each token in the context to be the start token of the answer span. The other vector indicates the probability of each token in the context to be the end of the answer span. The gold labels indicate the ground truth begin and end of the answer spans. DocQA is inherently bidirectional thus making it hard to process the input incrementally. It is possible to remove the bidirectionality by replacing the bidirectional layers with single-directional layers and replacing the global attention with an attention layer that only attends to the past, but this change significantly reduces performance.
Sliced DocQA In order to enable the DocQA model to process the context sequence incrementally despite having bidirectional RNN and attention layers, we use the concept of slicing [13]. We divide the context sequence into slices and apply the model to each slice. We call this "Sliced DocQA". We explore different ways of using Sliced DocQA for incremental question answering.

Sliced Prediction Layer
In the simple sliced model, we predict the output for each slice independently. Thus, as the model processes each slice, it predicts whether the answer span exists in that slice or not. We call this sliced prediction (see Figure 1a). To aggregate, we use a softmax layer on the concatenation of the outputs of all slices. This architecture allows us to processes the input incrementally, but each slice is now processed independently: i.e., when reading the second slice, the model can make no use of what was read in the first slice. This architecture hence ignores the order of the slices. We propose two solutions to solve this issue, discussed in more detail below. The first solution is to let the model access all representations from all slices in the prediction layer: we call this strategy global answer prediction. The second solution is a novel mechanism to transfer knowledge between slices called step transfer.
Global Answer Prediction While the sliced model processes the context sequence in slices, in the last layer, we use the encoded information from all the slices of the context to predict the answer. This means the prediction layer is not sliced. This model is illustrated in Figure 1b. Global answer prediction can be made incremental by applying the prediction once for each slice, with information from the future slices masked out.
Step Transfer To connect slices incrementally, we can transfer the knowledge between slices by using the knowledge learned until the current slice to initialize the next slice, as illustrated in Figure  1c. Thus, at the global level we have a uni-directional RNN, but, locally, we can have bidirectional or self-attention layers. To do this, we use encoded information from the current slice as an input to a fully connected network to predict the initial state for the next slice.
Early Stopping These incremental models can be used to support early stopping [3,12]. We use a supervised approach to predict when the model should stop, with a classifier that is simply trained on detecting whether an answer is contained in a given context or not. We train this classifier in a multi-task framework by adding an extra objective function to the QA system. Hence, at each training step, the model not only tries to predict an answer but also predicts whether the true answer is within the currently processed input or not. The early stopping classifier is a two-layer fully connected neural network with RELU activations and a sigmoid output layer. The input to this network is the average of the representations in the last layer of the current slice, and the output is a scalar indicating the probability of stopping. At test time the overall decision as to whether to stop early is based on thresholding the cumulative sum of these predictions, i.e. the probability of stopping at slice i is the sum of the predicted probabilities of stopping at slices 1 to i. Equation 1 shows how we compute the loss for the early stopping model.
Here extra_length i is a factor for scaling the early stopping loss, based on the distance from the gold stop point, in particular: extra_length = log(max(dist_threshold,|length_read−answer_end|)).
In this equation, dist_threshold is a minimum number of tokens after which we start scaling the loss. Thus, in case the model has not yet reached the answer and decides to stop it will be punished more if it is too early in contrast to when it will reach the answer by reading a few extra tokens.
Similarly, if the model has already passed the answer span, the more it reads, the bigger the loss will be. If we are in a distance less than the dist_threshold from the gold stop point, the scaling factor is log(dist_threshold). In the end, this loss is added to the answer prediction loss to form the total loss. If the true answer span is not before the chosen stopping point, the answer prediction loss is set to zero for all possible answers, and the model is only trained to choose a better stopping point. We Evaluate our model on SQuAD v1 [7] which is a reading comprehension dataset. It contains question and context pairs, where the answer to the question is a span in the context. We will open source the code for reproducing the experiments.

Experiments
We experiment with different slice sizes for Sliced DocQA in three different modes: with a sliced prediction layer, with a sliced prediction layer and step transfer, without step transfer but with a global prediction layer. In the sliced models, the size of each slice can play a crucial role. The extreme cases are when the length of the slices are 1, which means we have no bidirectional layer, and when the slices are as long as the whole sequence, which means we have a fully bidirectional model. Notice that when slice size is 1, the sliced model with step transfer is not equivalent to the uni-directional DocQA, since the attention layer is implemented differently. Figure 2 shows the performance of the Sliced DocQA models with respect to different slice sizes. With a sliced prediction layer, as expected, increasing the slice size leads to an increase in the performance. In this case, having a slice size of 1 means to predict the answer based on single word representations and still, we can get an accuracy of 31%. We can explain this by the fact  that 30% of answers in the SQuAD are single words, e.g. when the questions are of type when, where or who. We expect to find a slice size for which the performance of the model reaches the performance of the complete, not sliced, model. This is because, only local attention might be enough to correctly understand the meaning of the input sequence at each step, and also to predict whether each token is part of the answer or not. We observe that by slice size = 32, the model already reaches the performance of the uni-directional DocQA, and by slice size = 128 it almost reaches the performance of the normal DocQA. The interesting finding here is that there is a limit to which increasing the slice size leads to an increase in the performance, e.g. for the SQuAD, as soon as the slice size reaches 128 tokens, the increase in the performance is almost not noticeable.

Effect of slice size
Step transfer is useful When we have a step transfer mechanism, we observe that for each slice size, the performance of the model is higher. This means transferring knowledge between slices is useful. Among these models, the Sliced DocQA with step transfer models incremental comprehension and we see that for a slice size of 64 its performance is comparable with normal DocQA.
Global answer prediction is effective Surprisingly, with the global prediction layer, the effect of the slice size is much less. In general, when we have a global prediction layer, the performance of the model with all slice sizes is higher than when we have a sliced prediction layer. It is interesting to see that in this model, even with slice size of 1, which means using word embeddings to predict the answer span, we still get a performance close to when we have larger slice sizes.  Effect of early stopping In our early stopping experiments, we employed Sliced DocQA with Step Transfer as a truly incremental model. In Table 1, we compare the performance of this model with different slice sizes. As it can be seen, at slice size of 128 it almost achieves the performance of the model without early stopping (%99). Next, we investigate whether the early stopping model is more efficient regarding the number of read tokens. We looked into the performance of the model without early stopping at different context length to find the earliest point at which the model achieves its highest performance. This can be assumed as an oracle for early stopping model (i.e. best possible stopping point). In Figure 3a we compare the distribution of best stopping points with the stopping points of our early stopping model. We observe that (1)For a large number of examples, we have to read the full context to achieve the best performance. (2)For some examples, the early stopping model is reading more than it should, and for some, it stops earlier than it should. (3)In general, while with the best stopping offsets we can read about %15 less text, our model reads about %8 less. We also studied if there is a correlation between the context length and the ratio of read context length. In Figure 3b, we see for context length up to 128, we read full context, which is because our slice size is 128. After that point, we see the trend of reading smaller ratio of longer contexts.

Conclusion
In this paper, we propose a model that reads and comprehends text incrementally. As a testbed for our approach, we have chosen the question answering task. We aim to build a model that can learn incrementally from text, where the learning goal is to answer a given question. In standard question answering, we do not care how the context is presented to the model, and for the models that achieve state of the art results, e.g. [11,2], they process the full context before making any decisions. We show that it is possible to modify these models to be incremental while achieving similar performance. Having an incremental model, allows us to employ an early stopping strategy where the model avoids reading the rest of the text as soon as it reaches a state where it thinks it has the answer.


Density Power Downweighting and Robust Inference: Some New Strategies

Preserving the robustness of the procedure has, at the present time, become almost a default requirement for statistical data analysis. Since efficiency at the model and robustness under misspecification of the model are often in conflict, it is important to choose such inference procedures which provide the best compromise between these two concepts. Some minimum Bregman divergence estimators and related tests of hypothesis seem to be able to do well in this respect, with the procedures based on the density power divergence providing the existing standard. In this paper we propose a new family of Bregman divergences which is a superfamily encompassing the density power divergence. This paper describes the inference procedures resulting from this new family of divergences, and makes a strong case for the utility of this divergence family in statistical inference.

Introduction
In statistical modeling, parameter estimation is an inevitable and formidable task. Accurate estimation of the model facilitates the characterization and the subsequent understanding of the mechanism that generates the observed data. Statistical distances can be useful tools for the estimation of the model parameters.
Statistical distances can be naturally applied to the case of parametric statistical inference. The most important idea in parametric minimum distance inference is the quantification of the degree of closeness between the sample data and parametric model as a function of an unknown set of parameters through a suitable distance-like measure. Thus the estimate of the parameter is obtained by minimizing this "distance" over the parameter space.
It is worthwhile to mention here that the class of distances which we will consider are not mathematical metrics in the strict sense of the term. They may not be symmetric in their arguments and may not satisfy the triangle inequality. The only properties that we require of these measures are that they should be nonnegative, and should equal zero if and only if the arguments are identically equal. However, we will, somewhat loosely, continue to call them distances, or "statistical distances". In a practical sense, the word "divergence" is a good descriptor of these measures. We will, in fact, use the "minimum distance" and the "minimum divergence" terminologies interchangeably.
Density-based divergences form a special class of statistical distances. Several minimum distance estimators in this family have high model efficiency. In particular, the maximum likelihood estimator (MLE) also belongs to the class of density-based minimum distance estimators, being the minimizer of the likelihood disparity (Lindsay, 1994), which is a version of the Kullback-Leibler divergence. But one of the major drawbacks of the MLE is that it is notoriously nonrobust and even a small proportion of outlying observations can lead to meaningless inference. In fact it is the failure of the classical methods like maximum likelihood to deal with outliers and mild deviations from the model which had led to the emergence of the field of robustness; see, for example, [HR09], [HRRS86], [MMYSB19] and [BSP11]. However, some of the other members of the class of minimum distance estimators have been observed to do much better in the sense of combining strong robustness with high model efficiency. See, for example, [Csi63], [AS66], [Lin94], [Par05] and [BSP11] for a description of the φ-divergence class of minimum distance measures.
A more modern class of minimum distance estimators is based on the family of Bregman divergences. The Bregman divergence (Bregman, 1967) is a distance like measure between points and has been used in mathematics and information theory for some time. When the points are represented by probability distributions, the corresponding Bregman divergence is a statistical distance. See, for example, [JB90], [Csi91], [BMDG05] and [SV12] for some examples of statistical and related applications of the Bregman divergence. The principal representatives of Bregman divergence estimators in the current statistical literature are the minimum density power divergence estimators (MDPDEs), based on the density power divergence (DPD) class of [BHHJ98]. Over the last two decades, this class of divergences has provided a popular and frequently used method to balance the trade-off between robustness and efficiency in parameter estimation, hypothesis testing, and related inference. The minimum divergence estimators based on the DPD have been shown to provide a high degree of stability under model misspecification, often with minimal loss in model efficiency. Our primary purpose in this paper is to refine the minimum distance procedure based on the DPD, so as to achieve even better compromise between efficiency and robustness.
The Bregman divergence between the density g and model density f θ is given by where the index function B(·) is strictly convex and B (·) represents its first derivative with respect to its argument. In practice, where f θ is the pdf of the parametric family, g is the true density, the minimization of the above divergence over the parameter space Θ will generate the corresponding minimum distance functional which can lead to meaningful inference, depending on the form of the function B(·). The DPD, defined later in this section, is a special case of the Bregman divergence for B(y) = y 1+α α , α ≥ 0. When the model is differentiable, the general estimating equation under the divergence in Eq. 1 is or equivalently where u θ (x) = ∇ log f θ (x) is the score function of the model f θ (x), ∇ represents derivative with respect to θ and B (·) represents the second derivative of B(·) with respect to its argument. Since G is unknown, we construct an empirical version of the divergence in Eq. 1, or the estimating equation given in Eq. 3, by replacing G (the true data generating distribution) by its empirical counterpart G n . This leads to a class of unbiased (under the model) estimating equations The root of the Eq. 4 is defined to be the minimum Bregman divergence estimator (MBDE).
Here the robustness of the corresponding minimum distance estimator may be at least partially understood by observing the effect of the downweighting function B (f θ (x))f θ (x) on u θ (x) for less probable values of x under f θ . For the DPD, this weight becomes (α + 1)f α θ (x). In this paper we attempt to find a refinement of the DPD downweighting scheme, and, by reconstruction, a corresponding divergence. We will show that the corresponding minimum distance procedure provides a better compromise between robustness and efficiency in many cases compared to the minimum density power divergence estimator (MDPDE).

The Density Power Divergence
As mentioned earlier, the density power divergence (DPD) is obtained by substituting B(y) = y 1+α α in Eq. 1. The general form of this divergence, as a function of a nonnegative tuning parameter α, is For simplicity we have dropped the dummy variable in the above equation. One can define the minimum DPD functional T α (G) at G through the relation

A New Divergence
Our key philosophy for constructing new divergences and estimation strategies involves manipulating the downweighting factor B (f θ (X i ))f θ (X i ) in Eq. 4. Here we are going to develop a stronger downweighting effect compared to the MDPD estimating equation. Our exploration will generate an estimation scheme with two tuning parameters and we will explore the possibility of coming up with specific candidates which might beat the MDPDEs both in terms of efficiency and robustness.

Choosing the B Function
The downweighting effect on the score u θ (x) applied by the MDPD estimating equation is f α θ (x). As we want to impose a stronger downweighting in relation to this, we wish to choose the B function (or rather, the B function) so that as x → 0 + , xB (x) converges to zero faster than x α for α > 0 fixed. (Note, from Eq. 4, the downweighting term for u θ (x) in the general Bregman divergence is f θ B (f θ )). In particular, we will assume the following conditions on B .
(P2) xB (x) is an increasing function over x in (0, ∞). Thus the less likely observations will be downweighted more.
x β = 0, i.e., the Bregman formulation attaches weights to the score function which go to zero at a rate faster than the corresponding weights in the MDPD estimating equation.
To prove that such choice of B(·) satisfying (P1)-(P4) can help us generate divergences which have the desired properties and provide superior inference compared to the DPD, let us first demonstrate the general asymptotic properties of the minimum Bregman divergence estimators. For ease of representation, we refer to the divergence generated by the B(·) function satisfying (P1) to (P4) as φ-DPD.

General Asymptotic Properties of the MBDE
We need some regularity assumptions to prove the asymptotic properties of the general MBDE, which we list below.
(A1) The pdfs f θ of X have common support, so that the set X = {x|f θ (x) > 0} is independent of θ. The distribution G is also supported on X , on which the corresponding density g is greater than zero.
(A2) There is an open subset ω of the parameter space Θ, containing the best fitting parameter θ g (D B (g, f θ g ) = inf θ∈Θ D B (g, f θ )) such that for almost all x ∈ X , and all θ ∈ ω, the density f θ (x) is three times differentiable with respect to θ and the third partial derivatives are continuous with respect to θ. (The best fitting parameter θ g depends on the index function B(·) also, but we suppress that notation for brevity).
is bounded by some universal constant then the following hold. (a) The usual DPD defined in Eq. 5 is a special limiting case of φ-DPD. (b) If g = f θ for some θ ∈ Θ and if for the DPD there exists α, β such that the asymptotic relative efficiency (ARE) of the estimator under tuning parameter β is greater than that of the estimator under tuning parameter α, then there exists γ such that φ-DPD with tuning parameter (β, γ) generates an estimator with higher ARE than the MDPDE with tuning parameter α.
Our choice for B (·) in the LφDPD case is B (x) = 1 γ x β log 1 + γ x 0 < β ≤ 1, 0 < γ ≤ 1, x > 0. The corresponding B function may be expressed in the integral form as Obviously other choices are possible, but we have found the LφDPD to be a very useful divergence for our purpose, and for the rest of the paper all our illustrations will be in relation to the LφDPD. We will refer to the corresponding minimum distance estimator as MLφDE. θ (y) log 1 + γ f θ (y) is finite which is indeed the case for most parametric models suggesting the observed robustness of the MLφDE under those parametric models.
In Figure 1 it is clearly seen that the tuning parameter β has a significant impact on the robustness of the estimator and the influence functions redescend faster for larger values of β. On the other hand, for fixed β the influence functions are somewhat closer for different γ as seen in Figure 2. It suggests that γ has a less pronounced impact on robustness than β, although the graphs in Figure 2 indicate that larger γ lead to relatively stronger downweighting.

The Breakdown Point under the Location Model
. Consider the contamination model H ,n = (1 − )G + K n , where {K n } is a sequence of contaminating distributions. Let h ,n , g and k n be the corresponding densities. We say that there is breakdown in the minimum LφDPD functional for level contamination if there exists a sequence K n such that |T β,γ (H ,n ) − T β,γ (G)| → ∞ as n → ∞. We write below θ n = T β,γ (H ,n ) and assume that the true distribution belongs to the model family, i.e., g = f θ g . We make the following assumptions.
(BP1) min{f θ (x), k n (x)}dx → 0 as n → ∞ uniformly for |θ| ≤ c for any fixed c, i.e., the contamination distribution is asymptotically singular to the true distribution and to specified models within the parametric family.

Description and Results
Here we have performed a simulation study to analyze the performance of the LφDPD and the associated minimum distance estimators under the N (µ, 1) model at a given level of contamination. In the following study data are generated from two normal mixtures, 0.9N (0, 1) + 0.1N (5, 1) and 0.8N (0, 1) + 0.2N (5, 1), where N (0, 1) represents the target distribution and the second component is the contamination. The sample size is 50. The empirical MSE for the location model has been calculated by replicating the process 1000 times, evaluating the estimate for each replication and taking average squared error loss against the target value, i.e., µ = 0. In Table 1 the theoretical asymptotic relative efficiency of minimum LφDPD estimator and MDPDE is shown for different values of (β, γ) while in Table 2 and Table 3 the simulated mean square errors are presented under contaminated normal data under two different contamination levels.

The LφDPDersus the DPD
We briefly note our observations as may be evident from Tables 1 and 2. The asymptotic efficiencies of the minimum divergence estimators decrease with increasing β and increasing γ. Note that given an α ∈ (0, 1), it may be possible to choose β ∈ (0, α) and γ ∈ (0, 1) so that, in relation to our numerical study, MLφDE β,γ beats MDPDE α both in terms of β γ = 0 γ = 0.01 γ = 0.02 γ = 0.03 γ = 0.04 γ = 0.05 γ = 0.06 γ = 0.07 γ = 0.08 0.     Tables 1 and 2 (as also in Tables  1 and 3), there exists a better MLφDE, both in terms of asymptotic model efficiency and obtained mean square error under contamination. In most of these cases there are several (β, γ) combinations which provide the domination over a given MDPDE. Tables 2 and 3 also show that the robust minimum distance estimators hold out well against the outliers at both 10 and 20 percent contamination. Simulation results not presented here indicate that the same holds for higher levels of contamination smaller than 1/2, a consequence of the high breakdown point of the method under location models.   5 Algorithm for Finding the Optimal (β, γ) The LφDPD can generate many different kinds of estimators, starting from the most efficient estimator to highly robust estimators. For example, in the limit γ → 0 and β → 0, one gets the likelihood disparity which is minimized by the classical maximum likelihood estimator.
On the other hand, relatively larger values of β and γ lead to estimators with extremely high outlier stability. In a given situation, therefore, it is imperative that one is able to choose the most suitable tuning parameters for that particular case. Here we consider a data driven algorithm for selecting the "optimal" tuning parameters (β, γ) which would provide best compromise for the given situation. For this purpose we modify an approach of Warwick (2002), pp. 78-82, and minimize an empirical version of the asymptotic summed mean square error. The optimization technique is a two stage process. Suppose that the data are generated by a contaminated version of a model distribution, and let θ 0 be the parameter for the model component. Although the data are generated by a contaminated version, the parameter θ 0 of the model component is our target parameter. The spirit of such a set up is described in Warwick and Jones (2005). Let θ β,γ = T β,γ (G) be the corresponding minimum distance functional andθ β,γ is the solution of the unbiased equation of LφDPD with tuning parameter (β, γ) based on the data. The summed mean square error of the minimum LφDPD estimator has the asymptotic formula Here θ * is the pilot estimator playing the role of θ 0 and tr{·} represents the trace of matrix. The asymptotic covariance matrix of . So the estimated asymptotic summed mean square of the MLφDE is For the multiparameter case, the above quantity is a matrix. So trace of the matrix is used to provide a global measure of the summed mean square error for minimization. Thus when there are two parameters to be estimated (say (µ, σ) for N (µ, σ) model) then the expression to be minimized is The optimal value of (β, γ) is the minimizer of Eq. 19 under certain conditions. One important note is that in the first stage of minimization our pilot estimate for θ * is taken to be a good robust estimate based on the data as suggested in [War02]. The empirical summed mean square error is then obtained by evaluating the expressions in Eq. 18 or Eq. 19 after substitutingθ β,γ for θ β,γ and the empirical distribution G n in place of the true unknown distribution G. Let us denote this empirical summed mean square error by AMSE in the following.

Algorithm:
Given a dataset X n×1 we perform the following steps to obtain the estimate of θ.
1. Apply the method suggested in [War02] to get an optimal α for MDPDE. Suppose this value is α w . This step is the 1st stage of optimization by assuming an initial pilot estimate of θ * .

Real Data Examples
Here we take some real data sets and use our algorithm to find the optimal tuning parameters to be used in estimating the parameters of the model. We worked with two data sets, Newcomb's light speed data and Short's parallax of the sun data, under normality assumptions. We have used the minimum L 2 distance estimates as our pilot estimates of (µ, σ).

Newcomb's Data (Speed of Light)
This example involves Newcomb's light speed data (Stigler, 1977, Table 5). The data size is n = 66. Under the normal model, the MLE of the mean and standard deviation for these Figure 3: Normal density fits for Newcomb's data data are found to be equal to 26.212 and 10.664, respectively. We employ our algorithm for tuning parameter selection and Table 4 reports the optimal tuning parameters for DPD and LφDPD, as well as the parameter estimates at these optimal values. The estimators are extremely close, but the estimated asympmtotic summed mean square, for whatever it is worth, is lower in case of the MLφDE.

Short's Data (Parallax of the Sun)
This example involves Short's data for the determination of the parallax of the sun, the angle subtended by the earth's radius as if viewed and measured from the surface of the sun. From this angle and available knowledge of the physical dimensions of the earth, the mean distance from earth to the sun can be easily determined. The raw observations are presented in Table 4 of Stigler (1977). The data size is n = 53. Under the normal model, the MLE of the mean and standard deviation for these data are found to be equal to 8.378 and 0.846 respectively. We perform all the steps of the aforesaid tuning parameter selection algorithm, and the results of the analysis are now listed in Table 5. Again, the empirical asymptotic MSE for the MLφDE is slightly better than that of the MDPDE.
From Figure 3 and Figure 4, it is evident that the normal fits coming from the MDPDE

The MLφE for Independent Non-homogeneous Observations
Here we generalize the above concept to the case of independent but not identically distributed observations. [GB13] explains the methodology for this problem in the case of DPD, but here we will extend it to the case of LφDPD. Let us assume that the observed data Y 1 , ..., Y n are independent but for each i, Y i ∼ g i where the densities g 1 , ..., g n may not be same. We want to model g i by the family F i,θ = {f i (·; θ)|θ ∈ Θ} for all i = 1, 2, ..., n. We want to estimate θ by minimizing the LφDPD between the data and the model. However, the model density may not be same for each Y i 's, and hence we need to calculate the divergence between data and model separately for each data point. For this purpose, we minimize the average divergence between the data points with respect to θ ∈ Θ, where d(ĝ i , f i (·; θ)) denotes the LφDPD between the density estimate corresponding to the i-th data point and the associated model density. In the presence of only one data point Y i from density g i , the best possible density estimate of g i is the (degenerate) density which puts the entire mass on Y i so that we have where K is a constant independent of θ, the parameter of interest. Thus, for the purpose of estimation it suffices to minimize the objective function where Differentiating the above with respect to θ we get the estimating equation of the minimum LφDPD estimator for non-homogeneous observations as where u i (·) is the score function for f i (·).

Asymptotic Properties
We will now derive the asymptotic distribution of the minimum LφDPD estimatorθ n defined by the relation H n (θ n ) = min θ∈Θ H n (θ) provided such a minimum exists. Let us first present the necessary set up and conditions. Let the parametric model F i,θ be as defined above. We also assume that there exists a best fitting parameter of θ which is independent of the index i of the different densities. Let us denote it by θ g . The assumptions hold if all the true densities g i belong to the model family so that g i = f i (·; θ) for some common θ, and in that case the best fitting parameter is nothing but the true parameter θ. Next, recall that the MLφDEθ n is obtained as a solution of the estimating Eq. 23. This equation is satisfied by the minimizer of H n (θ) in Eq. 21. Similarly, we also define, for i = 1, 2, · · · , Note, at the best fitting parameter θ g , we must have We also define, for each i = 1, 2, · · · the p × p matrix J (i) whose (k, l)-th entry is given by where ∇ kl represents the partial derivative with respect to the indicated components of θ. We further define the quantities A simple calculation shows that, and where ξ i = 1 γ f i (y;θ) 0 s β log 1 + γ s ds g i (y; θ)dy.
We will make the following assumptions to establish the asymptotic properties of the MLφDE: (G1) The support X = {y|f i (y; θ) > 0} is independent of i and θ for all i; the true distributions G i are also supported on X for all i.
(G2) There is an open subset ω of the parameter space Θ, containing the best fitting parameter θ g such that for almost all y ∈ X , and all θ ∈ Θ, all i = 1, 2, · · · , the density f i (y; θ) is thrice differentiable with respect to θ and the third partial derivatives are continuous with respect to θ.
(G3) For each i = 1, 2, · · · , the three integrals f i (y; θ) f i (y;θ) 0 s β log 1 + γ s ds dy, f i (y;θ) 0 t 0 s β log 1 + γ s ds dt dy, and f i (y;θ) 0 s β log 1 + γ s ds g i (y)dy can be differentiated thrice with respect to θ, and the derivatives can be taken under the integral sign (the first indefinite integral).
(ii) The asymptotic distribution of Ω − 1 2 n Ψ n [ √ n(θ n − θ g )] is p-dimensional normal with (vector) mean 0 and covariance matrix I p , the p-dimensional identity matrix.
Note that, putting f i = f for all i, we get back the corresponding asymptotic properties of the minimum LφDPD estimator for the i.i.d. case. If f i = f, i = 1, 2, · · · , we get J (i) = J for all i; thus Ψ n = J and Ω n = K. Here J and K are as defined previously. In this case assumptions (G1)-(G5) are exactly the same as the assumptions (A1)-(A5), while assumptions (G6) and (G7) are automatically satisfied by the dominated convergence theorem. Thus the result, which establishes the consistency and asymptotic normality of the minimum LφDPD estimatorθ with n 1/2 (θ − θ g ) having the asymptotic covariance matrix Ψ −1 n Ω n Ψ −1 n = J −1 KJ −1 , emerges as a special case of Theorem 4.

Normal Linear Regression
A natural situation where the theory proposed above would be immediately applicable is the case of linear regression. We consider the linear regression model where the error i 's are i.i.d. normal variables with mean zero and variance σ 2 , x T i = (x i1 , · · · , x ip ) is the vector of the independent variables corresponding to the i-th observation and β = (β 1 , · · · , β p ) T represents the regression coefficients. We will assume that x i 's are fixed. Then y i ∼ N (x T i β, σ 2 ), and hence the y i 's are independent but not identically distributed. Thus y i 's satisfy our independent but non-homogenous set-up and hence the MLφDE of the parameter θ = (β T , σ 2 ) T can be obtained by minimizing the expression in Eq. 21 with f i ≡ N (x T i β, σ 2 ).

Real Data Examples in Regression
We now consider some real data examples to illustrate the above technique in linear regression.

Hertzsprung-Russel Data
This example involves a robust regression on the Hertzsprung-Russel data. These data, associated with the Hertzsprung-Russel diagram of the star cluster CYG OB1 containing 47 stars in the direction of Cygnus has been analyzed previously by several authors including [RL87].
We fit the simple linear regression model y = η 0 + η 1 x + under homoscedastic normal errors. Here the independent variable (x) is the logarithm of the temperature of the stars, and the dependent variable (y) is the logarithm of the light intensity of the stars. The initial regression parameter values are the least median of squares (LMS) estimates. The initial scale estimate is the scaled median absolute deviation (MAD) of the LMS residuals. We perform the previously mentioned steps of optimal tuning parameter selection and obtain the estimates for the regression coefficients, which are given in Table 6. The regression lines for LS regression, LMS regression and minimum LφDPD regression are given in the Figure  5. The robust performance of the MLφDE is self evident.

Salinity Data
This example involves the Salinity data (Table 5, Chapter 3, Rousseeuw and Leroy, 1987). These data were originally presented by [RC80]. The measurements of the salt concentration of the water and the river discharge taken in North Carolina's Pamlico Sound were recorded as the data. These data represent a multiple linear model with salinity as the dependent variable (y), and salinity lagged by two weeks (x 1 ), the number of biweekly periods elapsed We fit the multiple linear regression model y = η 0 + η 1 x 1 + η 2 x 2 + η 3 x 3 + under homoscedastic normal errors. The initial regression parameter values are the least median of squares (LMS) estimates. The initial scale estimate is the scaled median absolute deviation (MAD) of the LMS residuals.
The optimal parameters obtained through our algorithm for optimal parameter selection are presented in Table 7. The residual plots for LS regression, LMS regression and minimum LφDPD regression are given in the Figure 6. Like the LMS method (and unlike the LS method) the MLφDE gives a nice outlier resistant fit.

Hypothesis Testing using LφDPD
Now we develop the tests of parametric hypothesis based on LφDPD divergence. The most common problem is that of testing a simple null hypothesis for a parametric family of densities {f θ : θ ∈ Θ ⊂ R p } under the one sample case. Here we test Category MLφDE Tuning Parameter (β, γ) = (1, 0.9) Estimate of η 0 −8.5557324 Estimate of η 1 3.0590795 Estimate of σ 0.4266284 Table 6: Regression estimates for Hertzsprung-Russel data Figure 6: Residual plots of the fitted regression models for Salinity data using LS, LMS and minimum LφDPD estimation when a random sample X 1 , X 2 , . . . , X n is available from the population of interest. We propose our test statistic as withθ =θ β,γ being the MLφDE estimate of θ and B(·) is as defined in Eq. 16. We shall find the asymptotic distribution of T under H 0 and reject the null hypothesis for large values of T . We assume the following regularity conditions of the parametric family of distributions, (B1) The support of the distribution function F θ , i.e. the set X = {x|f θ (x) > 0} is independent of θ.
(B2) There is an open subset ω of the parameter space Θ, containing the true parameter value θ 0 such that for almost all x ∈ X , and all θ ∈ ω, the density f θ (x) is three times Category MLφDE Tuning Parameter (β, γ) = (1, 0.9) Estimate of η 0 57.16780461 Estimate of η 1 0.06010002 Estimate of η 2 −0.01301208 Estimate of η 3 −2.08372562 Estimate of σ 0.56157558 Table 7: Regression estimates for Salinity data differentiable with respect to θ and the third partial derivatives are continuous with respect to θ.
(B3) The integrals B (f θ (x))f 2 θ (x)dx can be differentiated with respect to θ, and the derivatives can be taken under the integral sign.
(B4) The p × p matrix J(θ) defined by is positive definite where E θ represents the expectation under the density f θ .
(B5) There exists functions M jkl (x) with finite expectation, j, k, l = 1, . . . , p, such that Then we have the following theorem.
We can extend this theorem and hence the testing result to the general two sample problem of testing H 0 : θ 1 = θ 2 against H 1 : θ 1 = θ 2 where there is a random sample of size n from population 1 with parameter θ 1 and that of size m from population 2 with parameter θ 2 . Letθ 1 andθ 2 be MLφDEs of the parameter in populations 1 and 2, respectively. Then under the (B1)-(B5) regularity conditions on the model, we have the following results.

Equivalence with the Score Test
A score test, developed in the same spirit under the same set up as in Theorem 5, also has the same asymptotic null distribution.
Theorem 7. The score test statistic using the LφDPD for testing the simple null in Eq. 35 can be given by and A(θ 0 ) is as described in Theorem 5. Under the null hypothesis, the asymptotic distribution of this statistic is same as that of T β,γ (θ, θ 0 ).

Divergence Difference test statistic
We assume that we have a parametric model F of densities and X 1 , . . . , X n be i.i.d. from the true distribution G with the same support as the distributions in F. Consider the null hypothesis where Θ 0 is a proper subset of Θ. The likelihood ratio test (LRT) is one of the most common tests that may be employed in this situation. Define where L(θ|X 1 , . . . , X n ) is the likelihood of θ given the data. The test statistic in this case is −2 log λ. Assume that the distribution function G is discrete. In particular let its support be X = {0, 1, 2, . . .}, which is also the common support of the family F. Then the test statistic can be expressed in terms of observed relative frequencies ν n as where LD(· , ·) stands for the likelihood disparity. Hereθ andθ 0 stands for unrestricted maximum likelihood estimator and maximum likelihood estimator under null hypothesis respectively. Eq. 38 gives a motivation to construct a new test statistic based on LφDPD.
As an analog of the likelihood ratio test, we consider the divergence difference test (DDT) based on LφDPD to test the hypothesis given in Eq. 37. Note that the test statistic in Eq. 38 can be viewed as a difference of the minimized value of likelihood disparity under null and unrestricted minimum of likelihood disparity. In the same spirit one may define the following test statistic DDT β,γ (ν n ) = 2n d β,γ (ν n , fθ θ 0 andθ are MLφDE under null hypothesis and unrestricted minimum MLφDE respectively. Also note that where B(·) is defined as Eq. 16. We will show that under certain regularity conditions the asymptotic distribution of the the test statistic DDT β,γ (ν n ) coincides with the distribution of linear combination of independent chi-squared random variables. Suppose that Θ 0 is defined by a set of r ≤ p restrictions on Θ defined by R i (θ) = 0, 1 ≤ i ≤ r. We assume that the parameter space under H 0 can be described through a parameter ξ = (ξ 1 , . . . , ξ p−r ), with p − r independent components, i.e., H 0 specifies that there exists a function b : The function b is assumed to have continuous derivativeḃ(ξ) of order p × (p − r) with rank p − r. Then the constrained estimator isθ 0 = b(ξ), whereξ is the MLφDE under the ξ formulation of the model. Let G = F θ be the true distribution which belongs to the family F with parameter θ. Under H 0 , let ξ be the true value of the reduced parameter. So we have θ = b(ξ). When the null hypothesis is true under standard regularity conditions it can be easily shown thatξ andθ 0 are consistent for ξ and θ respectively in the sense thatξ where Z n (b(ξ)) is AN (0, K B (b(ξ))). Here J B (·) and K B (·) is defined as in Theorem 1. Now we will lay out some appropriate regularity conditions under which we will derive the asymptotic distribution of DDT β,γ (ν n ) under the null hypothesis.
Remark. In the above theorems the null distribution of the test statistic turns out to be same as that of a linear combination of independent chi squared random variables. In general it is hard to get hold of critical values under this distribution for actually performing the test. Also calculations regarding this distribution become numerically hard. This gives the motivation to explore another test statistic which will lead to a simpler null distribution.

Wald Type Test
Assume a similar setup of hypothesis testing as in Eq. 37. Suppose that the null space Θ 0 ⊆ Θ ⊆ R p is defined by a set of r ≤ p restrictions on Θ defined by R i (θ) = 0, 1 ≤ i ≤ r. Let G = F θ be the true distribution which belongs to the family F with parameter θ. Assumeθ to be the MLφDE of the true parameter θ. Define R(θ) = (R 1 (θ), . . . , R r (θ)) T and D(θ) = ∂R i (θ) . Under the spirit of the original Wald test statistic, we can construct the following test statistic where Σ(θ) = J B (θ) −1 K B (θ)J B (θ) −1 under the B(·) function described in Eq. 16. Under standard regularity conditions it is easy to prove that the asymptotic distribution of W (θ) is χ 2 r under the null hypothesis. The proof follows from simple application of delta method theorem on the quantity R(θ) and the fact that under the null hypothesis √ n(θ − θ) is AN (0, Σ(θ)). The main benefit of this test statistic is that its asymptotic null distribution is simpler. Hence it is easy to perform numerical computations based on these statistics. For example, the critical values of the test statistic can be computed with ease in this case.

Real Data Example
Researchers needed to evaluate the effectiveness of an insecticide (dieldrin) in killing Anopheles farauti mosquitoes. The theory was that resistance to dieldrin was due to a single dominant gene, and that in an appropriately selected sample of the mosquitoes, there should be 50% susceptibility to insecticide. The hypothesis is where p is the probability of susceptibility. The results of such experiment is given in [Osb79]. The sample contains 465 mosquitoes where 264 of them died on being exposed to the insecticide. We can perform this test with test statistic DDT β,γ (ν n ) in Eq. 39. Here β and γ are chosen to be 0.3 and 0.05 respectively. The support of the distribution is X = {0, 1}, where the digit 1 stands for the death of a mosquito. From here it is evident that ν n (1) = 264/465. The null hypothesis is rejected if the value of the test statistic is large. In this case the asymptotic null distribution of the test statistic turns out to be 0.774χ 2 1 . Under the observed data the value of the test statistic turns out to be approximately 6.62. The 95% quantile of the aforementioned scaled chi-squared distribution is 2.97. So, under 5% level of significance the null hypothesis is rejected.

Summary
In this paper, we have developed a large class of density based divergences which includes the density power divergence family as a special limiting case. The key philosophy of stronger downweighting effect to construct the new family has been discussed. For application purposes, the family gives the data analyst a larger number of choices of possible divergences for inference purposes. We have shown several asymptotic and distributional properties of the proposed estimator. We have also shown that judicial choice of the tuning parameters leads to highly robust and efficient estimators which can often dominate the MDPDE. Though one of the parameters has a smaller effect on the robustness we have shown that both of them play an important role in the context of finite sample efficiency. We have also presented a possible data driven algorithm to obtain the "optimal" estimator in a given data set. We have also considered several hypothesis testing strategies for parameteric models which may serve as robust alternatives to the classical likelihood ratio and other likelihood based tests.
Remark. Like the MDPDE, the procedures described in this paper avoid the nonparametric density estimation and associated complications specific to classical minimum distance estimation. Another approach of this type can be found in [TB11].
Remark. In creating the test statistics for parametric hypothesis tesing using the LφDPD, we have restricted ourselves to the case where the same set of tuning parameters have been used for estimation as well as the construction of the subsequent divergences. In practice, one could allow them to vary; see, for example, [BMMP13]. In the present context, while this is possible, we do not explore this issue as we feel that there are enough tuning parameters involved already, and there are no demonstrated results indicating that such differential choices will necessarily produce improved tests.
Remark. In this paper, most of our illustrations have been with respect to the continuous model. Theoretically, however, there is nothing preventing its successful use in discrete models. All the necessary theories work out satisfactorily in this case.

Proof of Theorems
Proofs of Theorem 2, 5, 6 and 7 are skipped as they can be reproduced along the existing proofs in [BSP11], [GB13] and [GBP15].
under k n and f θn , the set A n converges to a set of zero probability as n → ∞. Thus, on A n , d(h ,n ) → d((1 − )g, 0) as n → ∞ and so by DCT | An d(h ,n , f θn ) − An d((1 − )g, 0)| → 0.
We will have a contradiction to our breakdown assumption if we can show that there exists a constant value θ * in the parameter space such that for the same sequence k n , lim sup n→∞ D(h ,f θn , f θn ) < a 1 ( ) as then the sequence {θ n } above could not minimize D(h ,f θn , f θn ) for every n. We will now show that above equation is true for all < 1/2 under the model when we choose θ * = θ g . For any fixed θ, let B n = {x : k n (x) > max{g(x), f θ (x)}}. Since g belongs to the model F, from (BP1) we get Bn g → 0, Bn f θ → 0 and B c n k n → 0 as n → ∞. Thus, under k n , the set B c n converges to a set of zero probability, while under g and f θ , the set B n converges to a set of zero probability. Thus, on B n , d(h ,n , f θ ) → d( k n , 0) = B( k n ) as n → ∞. (1) f,(1− ) decreases as ↑ 1. So, a 1 ( ) decreases as ↑ 1. Similarly it can be shown a 3 ( ) is an increasing function of . But a 1 (1/2) = a 3 (1/2); thus asymptotically there is no breakdown and lim sup n→∞ |T β,γ (H ,n )| < ∞ for < 1/2. Hence the theorem follows.


extfieldsofstudy:
Mathematics
Computer Science

text:
Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelastic models

The parameters in the governing system of partial differential equations of multicompartmental poroelastic models typically vary over several orders of magnitude making its stable discretization and efficient solution a challenging task. In this paper, inspired by the approach recently presented by Hong and Kraus~[Parameter-robust stability of classical three-field formulation of Biot's consolidation model, ETNA (to appear)] for the Biot model, we prove the uniform stability, and design stable disretizations and parameter-robust preconditioners for flux-based formulations of multiple-network poroelastic systems. Novel parameter-matrix-dependent norms that provide the key for establishing uniform inf-sup stability of the continuous problem are introduced. As a result, the stability estimates presented here are uniform not only with respect to the Lam\'e parameter $\lambda$, but also with respect to all the other model parameters such as permeability coefficients $K_i$, storage coefficients $c_{p_i}$, network transfer coefficients $\beta_{ij}, i,j=1,\cdots,n$, the scale of the networks $n$ and the time step size $\tau$. Moreover, strongly mass conservative discretizations that meet the required conditions for parameter-robust stability are suggested and corresponding optimal error estimates proved. The transfer of the canonical (norm-equivalent) operator preconditioners from the continuous to the discrete level lays the foundation for optimal and fully robust iterative solution methods. The theoretical results are confirmed in numerical experiments that are motivated by practical applications.

Introduction
Multiple-network poroelastic theory (MPET) has been introduced into geomechanics [8] to describe mechanical deformation and fluid flow in porous media as a generalization of Biot's theory [9,10]. The deformable elastic matrix is assumed to be permeated by multiple fluid networks of pores and fissures with differing porosity and permeability.
During the last decade, MPET has acquired many important applications in medicine and biomechanics and therefore become an active area of scientific research. The biological MPET model captures flow across scales and networks in soft tissue and can be used as an embedding platform for more specific models, e.g. to describe water transport in the cerebral environment and to explore hypotheses defining the initiation and progression of both acute and chronic hydrocephalus [50]. In [52,51] multicompartmental poroelastic models have been proposed to study the effects of obstructing cerebrospinal fluid (CSF) transport within an anatomically accurate cerebral environment and to demonstrate the impact of aqueductal stenosis and fourth ventricle outlet obstruction (FVOO). As a consequence, the efficacy of treating such clinical conditions by surgical procedures that focus on relieving the buildup of CSF pressure in the brain's third or fourth ventricle could be explored by means of computer simulations, which can also assist in finding medical indications of oedema formation [16].
Recently, the MPET model has also been used to better understand the influence of biomechanical risk factors associated with the early stages of Alzheimer's disease (AD), the most common form of dementia [24]. Modeling transport of fluid within the brain is essential in order to discover the underlying mechanisms that are currently being investigated with regard to AD, such as the amyloid hypothesis according to which the accumulation of neurotoxic amyloid-β (Aβ) into parenchymal senile plaques or within the walls of arteries is a root cause of this disease.
Biot's and multiple-network poroelastic models are challenging from a computational point of view in that the physical parameters for different practical applications exhibit extremely large variations. For instance, permeabilities in geophysical applications typically range from 10 −9 to 10 −21 m 2 while Young's modulus is of the order of GPa and the Poisson ratio in the range 0.1 − 0.3, see [53,39,18]. Permeabilities in biological applications typically range from 10 −14 to 10 −16 m 2 . Young's modulus of soft tissues is in the order of kPa and the Poisson ratio in the range 0.3 to almost 0.5, see, e.g., [48,49]. For that reason it is important that the problem is well posed and the numerical methods for its solution are stable over the whole range of values of the physical (model) and discretization parameters.
The stability of the time discretization and space discretization by finite difference or finite volume methods have been studied in [5,23,22,43] and will not be addressed here. Instead we focus on the issue of uniform inf-sup stable finite element discretizations of the static multiple-network poroelastic problem. It is well known that the well-posedness analysis of saddle-point problems in their weak formulation, apart from the boundedness and definiteness of the underlying bilinear form, relies on a stability estimate that is often referred to as Ladyzenskaja-Babuska-Brezzi (LBB) condition [11,19]. The LBB condition, see [6,14], is also crucial in the analysis of stable discretizations and in the derivation of a priori error etsimates for mixed problems. Inf-sup stability for the Darcy problem as well as for the Stokes and linear elasticity problems have been established under rather general conditions and various stable mixed discretizations of either of these problems have been proposed over the years, see, e.g. [11] and the references therein.
Biot's model of poroelasticity combines these equations and the parameter-robust stability of its classical three-field formulation has been established only recently in [34]. Alternative formulations that can be proven to be stable include a two-field formulation for the displacement and the pore pressure [12,1] and a new three-field formulation based on introducing the total pressure as a weighted sum of fluid and solid pressure as the third unknown besides the displacement and fluid pressure [44,39]. Contrary to this new three-field formulation as analyzed in [39], the classic three-field formulation of Biot's consolidation model considered in [34] builds on Darcy's law in order to guarantee fluid mass conservation, a property that the discrete models studied in this paper maintain. Aside from two-and three-field formulations, a four-field formulation has been considered for the Biot model in which the stress tensor is kept as a variable in the system, see [38]. The error analysis in the latter work is robust with respect to the Lamé parameter λ, but not uniform with respect to the other model parameters such as K. Another formulation for Biot's model has recently been proposed and analyzed in [7]. The authors use mixed methods based on the Hellinger-Reissner variational principle for the elasticity part of the system, and impose weakly the symmetry of the stress tensor σ, resulting in a saddle point problem for σ, u, p, and a Lagrange multiplier. They prove the parameter-robust stability of the resulting four-field formulation.
The first attempt to design parameter-robust discretizations and analyze their stability for the MPET model is presented in [37]. Motivated by [44,39], the authors of [37] propose a mixed finite element formulation based on introducing an additional total pressure variable. Utilizing energy estimates for the solutions of the continuous problem and a priori error estimates for a family of compatible semidiscretizations, they show that the formulation is robust in the limits of incompressibility, vanishing storage coefficients, and vanishing transfer between networks. The robustness with respect to the permeability coefficients remains an open question in [37].
There are various discretizations for the classic three-field formulation of Biot's model that meet the conditions for the proof of full parameter-robust stability that has been presented in [34]. In general, whenever a discretization is based on a Stokes-stable pair of finite element spaces for the displacement and pressure and a Poisson-stable pair of finite element spaces for the flux and pressure unknowns, it is possible to define a parameter-dependent norm (which in general is not uniquely determined) such that the constant in the inf-sup condition for the Biot problem does not depend on any of the model or disretization parameters. For example, the triplets CR l /RT l−1 /P l−1 (l = 1, 2) together with the stabilization techniques suggested in [25,31], see also [21], or the triplets P 2 /RT 0 /P 0 (in 2D) and P stab 2 /RT 0 /P 0 (in 3D), or P 2 /RT 1 /P 1 , or the stabilized discretization that has recently been advocated in [45], or the finite element methods proposed in [36] would qualify for such parameter-robustness. However, the above-mentioned finite element methods do not have the property of strong mass conservation in the sense of satisfying the mass balance equation pointwise and therefore locally and globally on a discrete level.
A priori error estimates for the continuous-in-time scheme and the discontinuous Galerkin (DG) spatial discretization (similar to [34]) have been presented in [32] for the Biot model. Inspired by the approach proposed in [34] in context of the static Biot problem, we make use of the DG technology in the present work for solving the MPET system by introducing novel parameter-matrix-dependent norms.
The aim of this work is to establish the results regarding the parameter-robust stability of the weak formulation of the continuous problem as well as the stability of strongly mass conservative discretizations, corresponding error estimates and parameter-robust preconditioners for the multiple-network (MPET) model. The presented stability results and error estimates and preconditioners are independent of all model and discretization parameters including the Lamé parameter λ, permeability coefficients K i , arbitrary small or even vanishing storage coefficients c pi , network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, the time step size τ and mesh size h. To our knowledge, these are the first fully parameterrobust stability results for the MPET model in a flux-based formulation.
The paper is organized as follows. In Section 2 the multiple-network poroelastic model is stated in a flux-based formulation, which can be considered as an extension of the classical three-field formulation considered in [34]. The governing partial differential equations are then rescaled and the static boundaryvalue problem resulting from semi-discretization in time by the implicit Euler method is presented in its weak formulation in the beginning of Section 3. The proofs of the uniform boundedness and the parameter-robust inf-sup stability of the underlying bilinear form are the main results that follow in this section. Section 4 then discusses a class of uniformly stable and strongly mass conservative mixed finite element discretizations that are based on H(div)-conforming discontinuous Galerkin approximations of the displacement field. Uniform boundedness and inf-sup stability are proved to be independent of all model and discretization parameters and the corresponding parameter-robust preconditioners are provided. Next, in Section 5, optimal parameter-robust error estimates are proved. Finally, Section 6 is devoted to the validation and illustration of the theoretical results in this work and Section 7 provides a brief conclusion. model are as follows: where σ = 2µǫ(u) + λdiv(u)I, (2a) In equation (2a), λ and µ denote the Lamé parameters defined in terms of the modulus of elasticity (Young's modulus) E and the Poisson ratio ν ∈ [0, 1/2) by The constants α i appearing in (1a) couple n pore pressures p i with the displacement variable u and are known in the literature as Biot-Willis parameters. The corresponding right hand side f describes the body force density. Each fluid flux v i is related to a specific negative pressure gradient −∇p i via Darcy's law in (1b). The tensors K i denote the hydraulic conductivities which give an indication of the general permeability of a porous medium. In (1c)u andṗ i express the time derivatives of the displacement u and the pressure variables p i . The constants c pi are referred to as the constrained specific storage coefficients and are connnected to compressibility of each fluid, for more see e.g. [47] and the references therein. The parameters β ij are the network transfer coefficients coupling the network pressures [50], hence β ij = β ji . The source terms g i in (1c) represent forced fluid extractions or injections into the medium.
It is assumed that the effective stress tensor σ satisfies Hooke's law (2a) where the effective strain tensor ǫ(u) is given by the symmetric part of the gradient of the displacement field, see (2b). Here I is used to denote the identity tensor.
The following boundary and initial conditions guarantee the well posedness of system (1): where for i = 1, . . . , n it is fulfilled Γ pi,D ∩ Γ pi,N = ∅, Γ pi,D ∪ Γ pi,N = Γ = ∂Ω and Γ u,D ∩ Γ u,N = ∅, Γ u,D ∪ Γ u,N = Γ. Initial conditions at the time t = 0 to complement the boundary conditions (3), have to satisfy (1a), and are given by The stress variable σ is eliminated from the MPET system by substituting the constitutive equation (2a) in (1a) thus obtaining the classical flux-based formulation of the MPET model.
To solve numerically the time-dependent problem, the backward Euler method is employed for time discretization resulting in the following system of time-step equations: where The unknown time-step functions u k , v k i , p k i for i = 1, . . . , n at any given time t = t k = t k−1 + τ are defined as whereas the right hand side time-step functions are . . , n. Later, the static problem (5)-(6) is considered and, for convenience, the superscript for the time-step functions is dropped, that is, u k , v k i and p k i will be denoted by u, v i and p i , respectively. The considered function spaces are as follows: • L 2 (Ω) is the space of square Lebesgue integrable functions equipped with the standard L 2 norm · ; • H 1 (Ω) d denotes the space of vector-valued H 1 -functions equipped with the norm · 1 for which u 2 1 := u 2 + ∇u 2 ; • H(div; Ω) := {v ∈ L 2 (Ω) d : div v ∈ L 2 (Ω)} with norm · div defined by v 2 div := v 2 + div v 2 . When the case Γ u,D = Γ pi,N = Γ and u D = 0, q i,N = 0 is considered, the notations U = H 1 0 (Ω) d and V i = H 0 (div, Ω), i = 1, . . . , n are used. To guarantee the uniqueness of the solution for the pressure variables p i , we set P i = L 2 0 (Ω) := {p ∈ L 2 (Ω) : Ω p dx = 0} for i = 1, . . . , n.

Stability analysis
First the parameter µ is eliminated from the system by dividing equations (5)-(6) by 2µ and making the substitutions:

Equation (5) then becomes
Next, equation (7b) is multiplied by α i τ −1 , equation (7c) is multiplied by α −1 i so that the substitutions For convenience, the "tilde" symbol is skipped and system (8) is written as: Further, we denote and make the rather general and reasonable assumptions that . . , n, and α ij ≥ 0 for i, j = 1, . . . , n.
Making use of these substitutions, without loss of generality, system (5) becomes where is the scaled operator from (6) . . , q n ) and V = V 1 × · · · × V n , P = P 1 × · · · × P n . Taking into account the boundary conditions, system (10) has the following weak formulation: Find (u; v; p) ∈ U × V × P , such that for any (w; z; q) ∈ U × V × P there holds Following [40], we first consider the following Hilbert spaces and weighted norms System (13), however, is not uniformly stable with respect to the parameters R −1 i under these norms as shown in [34]. Therefore, proper parameter-dependent norms for the spaces U , V i , P i , i = 1, . . . , n, have to be introduced that allow to establish the parameter-robust stability of the MPET model (13) for parameters in the ranges From experience, we know that the largest of the values R −1 i , i = 1, . . . , n is important to us, and we note that the term (ǫ(u), ǫ(w)) dominates in the elasticity form when λ ≪ 1. Hence, we define Again by trial and error, we find that we have to deal with the parameters in a "matrix" format. Therefore, we define the following n × n matrices j =i β ij and β ij = β ji , it is obvious that Λ 1 is symmetric positive semidefinite (SPSD). Since α pi ≥ 0, we have that Λ 2 is SPSD. Noting that R > 0, it follows that Λ 3 is symmetric positive definite (SPD). Moreover, it is obvious that Λ 4 is a rank-one matrix with eigenvalues λ i = 0, i = 1, . . . , n − 1 and λ n = n λ .
|Ω| Ω gdx and Λ g = [Λ 1 + Λ 2 , g c ] be the matrix that is obtained by augmenting Λ 1 + Λ 2 with the column g c . In general, we assume that Ω gdx = 0. When Λ 1 + Λ 2 is the zero matrix, this assumption is a "(classical) consistency condition". If Λ 1 +Λ 2 is nonzero and Ω gdx = 0, then g has to satisfy the "general consistency condition" rank(Λ 1 + Λ 2 ) = rank(Λ g ), where rank(X) denotes the rank of a matrix X. In this case, there must be p T c = (p 1,c , · · · , p n,c ) ∈ R n such that (Λ 1 + Λ 2 )p c = g c (in many applications, |Ω| Ω gdx, and thus Ω g 0 dx = 0. Then the solution (u; v; p) can be decomposed according to (Ω) and p c is a basic solution of (Λ 1 +Λ 2 )p c = g c . Therefore we only need to consider the case when Ω gdx = 0.
As we will see, it will play an important role in the definition of proper norms and the splitting (19) in our analysis. The crucial idea is that we equip the Hilbert spaces U , V , P with parameter-matrix-dependent norms · U , · V , · P induced by the following inner products: . It is easy to show that (20a)-(20c) are indeed inner products on U , V , P respectively. It should be noted that Div v, Div z and p, q are vectors and the SPD matrix Λ is used to define the norms. These novel parameter-matrix-dependent norms play a key role in the analysis of the uniform stability for the MPET model. We further point out that for n = 1, the norms defined by (20) are slightly different, but equivalent to the norms that were used in [34] to establish the parameter-robust inf-sup stability of the three-field formulation of Biot's model of consolidation.
The main result of this section is a proof of the uniform well-posedness of problem (13) under the norms induced by (20). Firstly, directly related to problem (13), we introduce the bilinear form which, in view of the definition of the matrices Λ 1 and Λ 2 , can be written in the form Then the following theorem shows the boundedness of A((·; ·; ·), (·; ·; ·)) in the norms induced by (20).
There exists a constant C b independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n and the network scale n, such that for any Proof. From the definition of the bilinear form, by using Cauchy's inequality, we obtain Then, another application of Cauchy's inequality, in view of the definition of Λ 4 , yields Before we study the uniform inf-sup condition for the MPET equations, we recall the following well known results, see, e.g. [14,11]: Furthermore, we summarize some useful properties of the matrix Λ in the following lemma.
thenΛ is SPD and for any n-dimensional vector x, we have Also, Proof. From the definitions of Λ 3 , Λ 4 , noting that Λ 3 is SPD and Λ 4 is SPSD, it is obvious thatΛ is SPD. From the definition of Λ, noting that Λ 1 and Λ 2 are SPSD, we infer the estimates From the definitions of Λ 3 , Λ 4 andΛ, we havẽ Now, using the Sherman-Morrison-Woodbury formula, we find Further, noting that where I n×n is the n-th order identity matrix, we obtain Now we can calculateΛ −1 as follows: We are ready to prove the uniform inf-sup condition for A((·; ·; ·), (·; ·; ·)) in the norms induced by (20).
Remark 5. We want to emphasize that the parameter ranges as specified in (17) are indeed relevant since the variations of the model parameters are quite large in many applications. For that reason, Theorem 2 and Theorem 3 are very important fundamental results that provide the parameter-robust stability of the model (13a)-(13c). We also point out that the matrix technique plays an interesting role for proving the uniform stability.
where Inferring from the theory presented in [42], Theorems 2 and 3 imply that the operator B defined in (36) is a uniform norm-equivalent (canonical) block-diagonal preconditioner for the operator A in (12), robust in all model and discretization parameters, i.e., κ(BA) = O(1).

Uniformly stable and strongly mass conservative discretizations
There are various discretizations that meet the requirements for the proof of full parameter-robust stability as presented in this section. They include conforming as well as nonconforming methods. In general, if is a Stokes-stable pair and V i,h /P i,h satisfy the H(div) inf-sup condition for i = 1, . . . , n, see (55), then the norm that we have proposed in Section 3 allows for the proof of full parameter-robust stability using similar arguments as in the proof of Theorem 10. To give a few examples, the triplets U h /V i,h /P i,h = CR l /RT l−1 /P l−1 (l = 1, 2) together with stabilization [25,21] results in a parameterrobust stable discretization of the MPET model if the norms are defined as in Section 3. The same is true for the conforming discretizations based on the spaces P 2 /RT 0 /P 0 (in 2D), P stab 2 /RT 0 /P 0 (in 3D), or P 2 /RT 1 /P 1 . However, the above-mentioned finite element methods do not have the property of strong mass conservation in the sense of Proposition 8 although they result in parameter-robust inf-sup stability under the norms we proposed in Section 3.
In recent years, DG methods have been developed to solve various problems [3,15,4,17,26] and some unified analysis for finite element including DG methods has recently been presented in [29,30]. In this section, motivated by the works [46,28,27], we propose discretizations of the MPET model problem (13). These discretizations preserve the divergence condition (namely equation (10c)) pointwise, which results in a strong conservation of mass, see Proposition 8. Furthermore, they are also locking-free when the Lamé parameter λ tends to ∞.

Preliminaries and notation
By T h we denote a shape-regular triangulation of mesh-size h of the domain Ω into triangles {K}. We further denote by E I h the set of all interior edges (or faces) of T h and by E B h the set of all boundary edges (or faces); we set E h = E I h ∪ E B h . For s ≥ 1, we define H s (T h ) = {φ ∈ L 2 (Ω), such that φ| K ∈ H s (K) for all K ∈ T h }.
As we consider discontinuous Galerkin (DG) discretizations, we also define some trace operators. Let e = ∂K 1 ∩ ∂K 2 be the common boundary (interface) of two subdomains K 1 and K 2 in T h , and n 1 and n 2 be unit normal vectors to e pointing to the exterior of K 1 and K 2 , respectively. For any edge (or face) e ∈ E I h and a scalar q ∈ H 1 (T h ), vector v ∈ H 1 (T h ) d and tensor τ ∈ H 1 (T h ) d×d , we define the averages and jumps where v ⊙ n = 1 2 (vn T + nv T ) is the symmetric part of the tensor product of v and n. When e ∈ E B h then the above quantities are defined as If n K is the outward unit normal to ∂K, it is easy to show that Also, for τ ∈ H 1 (Ω) d×d and for The finite element spaces we consider are denoted by U h = {u ∈ H(div; Ω) : u| K ∈ U (K), K ∈ T h ; u · n = 0 on ∂Ω}, ; v · n = 0 on ∂Ω}, i = 1, . . . , n, The discretizations that we analyze in the present context define the local spaces U (K)/V i (K)/Q i (K) via the triplets BDM l (K)/ RT l−1 (K)/P l−1 (K), or BDF M l (K)/RT l−1 (K)/P l−1 (K) for l ≥ 1. Note that for each of these choices, the important condition div We recall the following basic approximation properties of these spaces: For all K ∈ T h and for all u ∈ H s (K) d , there exists u I ∈ U (K) such that

DG discretization
We note that according to the definition of U h , the normal component of any u ∈ U h is continuous on the internal edges and vanishes on the boundary edges. Therefore, by splitting a vector u ∈ U h into its normal and tangential components u n and u t , we have implying that A direct computation shows that Similar to the continuous problem, we denote ,h × · · · × P n,h . With this notation at hand, the discretization of the variational problem (13) is given as follows: , and η is a stabilization parameter independent of parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Remark 7. Consider the general rescaled boundary conditions
Usually, it is assumed that the measure of Γ u,D is nonzero to guarantee the discrete Korn's inequality [13].
The standard way to incorporate the boundary conditions (46) is to modify the trial spaces according to the boundary conditions, i.e., to seek the solution in the spaces and use the test spaces given by and and η is again a stabilization parameter which is independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h. If Γ u,D = Γ p,N = Γ and u D = 0, q N = 0, then (47) reduces to (44) which will be analyzed in the remainder of this paper. If the measure of Γ u,N is nonzero, then the analysis is similar. If Γ u,D = Γ and the measure of any Γ pi,D , i = 1, · · · , n, is nonzero, then one has to modify the norms according to Remark 3.1 in [34]. This part of the analysis is left as future work.
where Q i,h denotes the L 2 -projection on P i,h .
For any u ∈ U h , we introduce the mesh dependent norms: Next, for u ∈ U h , we define the "DG"-norm and, finally, the mesh-dependent norm · U h by We now summarize several results on well-posedness and approximation properties of the DG formulation, see, e.g. [28,27]: • From the discrete version of Korn's inequality we have that the norms · DG , · h , and · 1,h are equivalent on U h , namely, • The bilinear form a h (·, ·), introduced in (45) is continuous and we have • For our choice of the finite element spaces U h , V h and P h we have the following inf-sup conditions, see, e.g. [46]: where β sd and β dd are positive constant independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.
• We also have that a h (·, ·) is coercive, and the proof of this fact parallels the proofs of similar results: where α a is a positive constant independent of parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Related to the discrete problem (44a)-(44c) we introduce the bilinear form
In view of the definitions of the norms · U h , · V and · P , the boundedness of the bilinear form A h ((u h ; v h ; p h ), (w h ; z h ; q h )) is obvious, i.e., the following theorem holds.
Theorem 9. There exists a constant C bd independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h, such that for any (u We come to our second main result.
Theorem 10. There exits a constant β 0 > 0 independent of the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h, such that The proof of this theorem can be obtained by following the proof of Theorem 3 and using the technique shown in [34].
From the above theorem, we get the following stability estimate.
be the solution of (44a)-(44c), then we have the estimate q h P and C 2 is a constant independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.

Remark 12. Define
where Then due to the theory presented in [42], Theorems 9 and 10 imply that the norm-equivalent (canonical) block-diagonal preconditioner B h for the operator induced by the bilinear form (57) is uniform with respect to variation of the model and dicretization parameters. This means that the condition number κ(B h A h ) is uniformly bounded with respect to the parameters λ, R −1 i , α pi , α ij , i, j = 1, . . . , n in the ranges specified in (17), the network scale n and the mesh size h.
To apply the preconditioner B h , one has to solve an elasticity system discretized by an H(div)-conforming discontinuous Galerkin method [28] and n coupled elliptic H(div) problems discretized by RT elements. In the lowest order case and for n = 1, optimal solvers for this task have been proposed in [35].

Error estimates
In this section, we derive the error estimates that follow from the results presented in Section 4. Let Π div B : H 1 (Ω) d → U h be the canonical interpolation operator. We also denote the L 2 -projection on P i,h by Q i,h . The following Lemma, see [28], summarizes some of the properties of Π div B and Q i,h needed for our proof.
Theorem 13. Let (u; v; p) be the solution of (13) and (u h ; v h ; p h ) be the solution of (44a)-(44c). Then the error estimates and hold, where C e,u , C e,p are constants independent of λ, R −1 i , α pi , α ij , i, j = 1, . . . , n, the network scale n and the mesh size h.
Proof. Subtracting (44a)-(44c) from (13a)-(13c) and noting the consistency of a h (·, ·), we have that for Let by the stability result (58) for the discrete problem (44a)-(44c), we obtain Using the boundedness of a h (·, ·), the second inequality in Lemma 4, the triangle inequality and noting that v I is arbitrary and (Div and Remark 14. From the above theorem, we can see that the discretizations are locking-free.

Numerical Experiments
The following numerical experiments are for three widely applied MPET models, namely the one-network, two-network and four-network models. We suppose that the domain Ω is the unit square in R 2 and during the discretization it has been partitioned as bisections of 2N 2 triangles with mesh size h = 1/N . To discretize the pressure variables we use discontinuous piecewise constant elements, the fluxes are discretized employing the lowest-order Raviart-Thomas space and the displacement we approximate with the Brezzi-Douglas-Marini elements of lowest order. All the numerical tests included in this section have been carried out in FEniCS, [2,41]. The aim of these experiments is: (i) to validate the convergence of the error estimates in the derived parameter-dependent norms; (ii) to test the robustness of the proposed block-diagonal preconditioners by using it within the MinRes algorithm.
Then the exact solution of system (10) with boundary conditions u| ∂Ω = 0, v · n| ∂Ω = 0 is given by u = ∂φ1 ∂y , − ∂φ1 ∂x , p = φ 2 − 1, v = −R 1 ∇p and p ∈ L 2 0 (Ω). We performed experiments with different sets of input parameters. In Tables 1-3 we report the error of the numerical solution in the introduced parameter-dependent norms · P , · V , · U h . Additionally, we list the number of MinRes iterations n it and average residual convergence factor with the proposed blockdiagonal preconditioner where the stopping criterion is residual reduction by 10 8 in the norm induced by the preconditioner. The robustness of the method is validated with respect to variation of the parameters λ, R −1 1 , α p1 , as introduced in (10), and the discretization parameter h. As can be seen from Tables 1-3 the error in the considered parameter-dependent norms decreases by a factor 2 when decreasing the mesh size by the same factor independently of the model parameters. The results in Table 4 suggest that the number of MinRes iterations required to achieve a prescribed solution accuracy is bounded by a constant independent of λ, R −1 1 , α p1 and h while the average residual reduction factor always remains smaller than 0.70. Note that in this table the authors have tried to present the most unfavourable setting of input parameters in order to stress test the proposed method.

The two-network model
The governing partial differential equations of the Biot-Barenblatt model in which the flux-based MPET system involves two pressures and two fluxes are given by We consider here the cantilever bracket benchmark problem proposed by the National Agency for Finite Element Methods and Standards in [20] with f = 0, g 1 = 0 and g 2 = 0. The boundary of the domain Ω = [0, 1] 2 is split into Γ 1 , Γ 2 , Γ 3 and Γ 4 denoting the bottom, right, top and left boundaries respectively, and the boundary conditions u = 0 on Γ 4 , (σ − p 1 I − p 2 I)n = (0, 0) T on Γ 1 ∪ Γ 2 , (σ − p 1 I − p 2 I)n = (0, −1) T on Γ 3 , p 1 = 2 on Γ p 2 = 20 on Γ are imposed.
The base values of the model parameters are taken from [33] and are presented in Table 5. The computed numerical results in Table 6 show robust behaviour with respect to mesh refinements and variation of the parameters including high contrasts of the hydraulic conductivities. The parameter K 2 has been varied over a wider range than K 1 as it appeared to be the more interesting case.
Here, it is important to note that the authors have attempted to present again the least optimal choice of parameters for testing their implementation.

Conclusions
In this paper, motivated by the approach recently presented by Hong and Kraus [Parameter-robust stability of classical three-field formulation of Biot's consolidation model, ETNA (to appear)] for the Biot model, we establish the uniform stability, design stable disretizations and a parameter-robust preconditioners for flux-     216 Nm −2 c p1 = c p2 = c p3 = c p4 4.5 · 10 −10 m 2 N −1 α 1 = α 2 = α 3 = α 4 0.99 β 12 = β 24 1.5 · 10 −19 m 2 N −1 s −1 β 23 2.0 · 10 −19 m 2 N −1 s −1 β 34 1.0 · 10 −13 m 2 N −1 s −1 K 1 = K 2 = K 4 = K (1.0 · 10 −10 )/(2.67 · 10 −3 ) m 2 /Nsm −2 K 3 (1.4 · 10 −14 )/(8.9 · 10 −4 ) m 2 /Nsm −2 based formulations of multiple-network poroelastic systems. Novel proper parameter-matrix-dependent norms that provide the key for establishing uniform inf-sup stability of the continuous problems are introduced. The stability results that could be obtained using the presented matrix technique are uniform not only with respect to the Lamé parameter λ but also with respect to all the other model parameters such as small or large permeability coefficients K i , arbitrary small or even vanishing storage coefficients c pi , arbitrary small or even vanishing network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, and the time step size τ . Moreover, strongly mass conservative and uniformly stable discretizations are proposed and corresponding uniform and optimal error estimates proved which are also independent of the Lamé parameter λ, the permeability coefficients K i , the storage coefficients c pi , the network transfer coefficients β ij , i, j = 1, · · · , n, the scale of the networks n, the time step size τ and the mesh size h. The transfer of the canonical (norm-equivalent) operator preconditioners from the continuous to the discrete level lays the foundation for optimal and fully robust iterative solution methods. Numerical experiments that are motivated by practical applications are presented confirming both the uniform and optimal convergence of the proposed finite element methods and the uniform robustness of the norm-equivalent preconditioners.