Reza Sameni
email: reza.sameni@gmail.com

Spatio-Temporal Source Separation using Temporal Priors with Parameterized Uncertainties

The objective of this progressive research is to develop a spatio-temporal source separation scheme using temporal priors with parametric uncertainties. The method has potential applications in single and multi-modal data analysis. It is shown that the proposed scheme can be jointly used with classical source separation techniques based on generalized eigenvalue decomposition (GEVD) and joint approximate diagonalization of eigenmatrices (JADE). As a case study, the proposed method is evaluated over noninvasive fetal electrocardiogram signals.

I. INTRODUCTION

When dealing with multi-modal signals-such as the simultaneous recording of the phonocardiogram and electrocardiogram-a recurrent issue is how to deal with simultaneous recordings with a common (yet hidden) origin; but with totally different temporal representations. In a progressive research, we have been studying various source separation techniques, with potential application in multimodal data analysis. Herein, we present preliminary results of this progressive research, which has led to the development of spatio-temporal filters, which are applicable to signals belonging to some general class of temporal structure; but with some parametric uncertainties in their temporal priors.

II. PROBLEM DEFINITION

Consider N -channel discrete-time noisy observations x i (t) (t = 1, • • • , T ) and a set of linear (and generally time-variant) filters h i (t, τ ) (i = 1, • • • N ). Defining the filtered observed channels

y i (t) = h i (t, τ ), x i (τ ) τ = T τ =1 h i (t, τ )x i (τ ), (1)
the objective is to find linear mixtures

z m (t) = N i=1 w mi y i (t), m = 1, • • • M (2)
with maximal total channel energy under some bounding conditions on the weight vector w ki (as detailed below).

To start with, if we simplify the problem to fixed linear filter coefficient sets h i (t, τ ) = h(t, τ ) over all channels, the problem can be stated in matrix form as follows

Z = WXH (3)
where

X = S + η (4)
is the matrix form of the noisy observations

X ∈ R N ×T , S ∈ R N ×T is the desired multi-channel signal, η ∈ R N ×T is additive noise, Z ∈ R M ×T is the filtered data, H ∈ R T ×T
is the matrix form of h(t, τ ) considered as a finite impulse response (FIR) temporal filter, and W ∈ R M ×N is the spatial filter (in fact each row of W is a spatial filter). The desired source signals S may generally be mixtures of some latent variables θ (e.g., S = Aθ as in classical ICA); but we do not explicitly use this fact in the current study. With these definitions, we seek the spatial filter, such that

W * = arg max W tr(WXHH T X T W T ) tr(WXX T W T ) (5)
which is equivalent to the joint diagonalization of the following two matrices using generalized eigenvalue decomposition (GEVD):

A = XHH T X T , B = XX T (6)

III. KNOWN (PRESUMED) TEMPORAL STRUCTURE

As a linear filter, H can generally have an arbitrary structure. For time-invariant filtering, the filter takes a column-wise Toeplitz form, obtained from a prototype impulse response

vector h = [h 1 , h 2 , • • • , h L ] T (L ≤ T ).
The filter can be causal, anti-causal, or non-causal, depending on whether H is upper-triangular, lower-triangular, or non-triangular. The impulse response of the filter may have been obtained from a classical filter design procedure, or an optimal design scheme, such as a matched filter or Wiener filter (smoother) design. It can also be a linear map to a transform domain, such as the wavelet or Fourier domain. In the latter case, H takes the form of a DFT or DCT matrix.

In the more general case of nonlinear filters with different per-channel filter responses, (2) becomes

z m (t) = N i=1 w mi h i (x i (t)), (m = 1, • • • M ) (7)
(notice the difference from the classical post-nonlinear data model [START_REF] Taleb | Source separation in post-nonlinear mixtures[END_REF]). In this case, the diagonalized matrices defined in (6) can be replaced by a symmetric kernel1

K X = [ h i (x i (t)), h j (x j (t)) t ] ij (8)
Depending on the choice of the temporal filter h(•), the problem reduces to classical time-domain separation methods such as AMUSE [START_REF] Tong | Amuse: a new blind identification algorithm[END_REF], TDSEP [START_REF] Ziehe | Tdsepan efficient algorithm for blind separation using time structure[END_REF], [START_REF] Molgedey | Separation of a mixture of independent signals using time delayed correlations[END_REF], πCA [START_REF] Saul | Periodic Component Analysis: An Eigenvalue Method for Representing Periodic Structure in Speech[END_REF], [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF], and a few other time and transform domain source separation methods described in [START_REF]A Deflation Procedure for Subspace Decomposition[END_REF].

The overall problem is trivial with a known (or set of known) prior temporal filter(s)-either in the linear or kernelized form-using GEVD. But how can one solve the problem if the temporal structure is unknown? In the most general case of unknown temporal structure, the problem is ill-posed as we have more parameters than observations (H is a huge T × T unknown matrix). In the kernelized form, the problem reduces to kernel learning problem, which can be solved by prior assumptions, such as stochastic independence, on the desired sources. This approach is the basis of the well known kernel ICA [START_REF] Bach | Kernel independent component analysis[END_REF] and similar techniques.

Herein, we study the special case of temporal priors, with parametric uncertainties in their temporal structure. Smoothness priors are studied as a special case for this framework. It is shown that the problem has a low computational complexity and outperforms the fully blind scenario (especially in presence of noise), as it makes the best of the temporal priors despite their uncertainties.

IV. PARAMETRIC SPATIO-TEMPORAL FILTER DESIGN

Let us suppose that we know some general properties regarding the temporal structure of the desired sources. For instance, one might know (or assume) that the signal follows an auto-regressive model of a certain order and spectra, or that the signal is smooth up to a certain order (for example by knowledge of its spectral decay rate [9, Ch. 2]), etc. At the same time, there can be some uncertainties in the signal model, such as the exact smoothness order or level of background noise (SNR). Based on these assumptions, we propose two spatio-temporal source separation schemes to overcome the uncertainties in the temporal structure.

A. Multiple smoothness parameters using JADE

In the first approach, we assume that we have multiple temporal priors regarding the sources, or a single prior with unknown set of parameters. In case of multiple temporal priors, one can have different temporal filters represented in matrix form as

H k (k = 1, • • • , K).
The filters can be totally different (each using a different temporal prior) or may belong to the same class of filters; but with different set of parameters. In this case, the problem is to find the spatial filter that jointly (approximately) diagonalizes [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF], [START_REF] Cardoso | Perturbation of joint diagonalizers[END_REF], the matrices

A k = XH k H T k X T (k = 1, • • • , K), B = XX T (9) (the matrix B is only required in the diagonalization set, if we seek uncorrelatedness).

B. Iterative regularization factor estimation using GEVD)

In the second approach, we use the fact that the linear form assumed in [START_REF] Ziehe | Tdsepan efficient algorithm for blind separation using time structure[END_REF], implies that the spatial and temporal filters are fully decoupled. In other words, the temporal filtering does not mix the channels in space and the spatial mixture does not mix the data in time and the order of applying the filters in time and space is not important. However, due to the nonlinear (quadratic) form of the cost function in [START_REF] Saul | Periodic Component Analysis: An Eigenvalue Method for Representing Periodic Structure in Speech[END_REF], the two filters can enhance the performance of one other, in the sense that one can seek "a spatial filter that enhances some temporal structure and vice versa a temporal filter that improves the spatial structure". Based on this idea, we propose to use successive iterations of temporal and spatial filtering, each time improving the estimate of one another. The general idea is summarized in Algorithm 1. We can notice that the idea is generic and may be extended to other parametric spatio-temporal source separation schemes and for kernelized formulations. In the next section, we will study this iterative approach for the special case of signals with parameterized smoothness priors.

1: Calculate C X = XX T 2: for i = 1, • • • , I do 3: Calculate Y = XH(θ i-1 ) 4:
Calculate K X = YY T

5:

Find the spatial filter W, using GEVD on the matrix pair (K X , C X )

6: Find Z = WX = [z 1 , • • • , z N ] T 7:
Optimize the temporal filter parameter set θ i over z 1 , (or over all channels of Z i in the kernelized form)

8: return Z and W

V. CASE STUDY: PARAMETRIC TEMPORAL SMOOTHNESS

As proof of concept, we study the problem of extracting smooth signals from multichannel mixtures of signal and noise. The order of signal smoothness and the background noise level are assumed to be unknown 2 . Some related studies in this field, which have used the slow temporal variations of a signal for source separation, include [START_REF] Blaschke | Independent slow feature analysis and nonlinear blind source separation[END_REF], [START_REF] Blaschke | Independent slow feature analysis and nonlinear blind source separation[END_REF].

Numerous methods have been developed for data smoothing in the literature. In this context, smoothing regularization (also known as Tikhonov regularization) is a popular technique based on constrained least squares estimation [14, p. 307]. Consider the second-order difference operator

d 2 ∆ = [1, -2, 1] + Data model Wiener Smoother ¡ ¢ £ ¤ ¥ ¦ § © T K < ! " # T Fig.
1. An equivalent Wiener smoother scheme for Tikhonov regularization, adopted from [START_REF] Sameni | Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising[END_REF].

and its Toeplitz matrix form

D 2 =                   1 0 0 • • • 0 -2 1 0 . . . . . . 1 -2 1 . . . . . . 0
0 • • • 0 0 1                   (10)
The right multiplication of D 2 in any matrix X containing length T signals as its rows is equivalent to convolving each row with d 2 , i.e., the second order difference of the signal. This definition can be extended to order-d difference operators d d and D d [START_REF] Sameni | Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising[END_REF]. Using a constrained least squares optimization scheme, it has been shown that the following matrix acts as a smoothing operator

H(d, λ) = (I + λD d D T d ) -1 (11)
where d is the difference order and λ ≥ 0 is the regularization parameter (λ = 0 for no smoothing) considered as an unknown parameter in this work.

In a recent study [START_REF] Sameni | Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising[END_REF], it was shown that H(•) is equivalent to a linear zero-phase forward backward Wiener smoother, with the following per-channel transfer function

ψ d (z) = 1 1 + λd d (z)d d (z -1 ) = S s (z) S s (z) + S n (z) (12)
which is in the form of an optimal Wiener smoother for a signal s(t) having an auto-regressive all-pole model with a spectra denoted by S s (z), driven by input white (process) noise with variance 2 , in presence of white observation noise n(t) with a white spectra denoted by S n (z) and variance δ 2 . Fig. 1, illustrates the analogy between Tikhonov regularization and optimal Wiener smoothing in the sample-wise formulation. With this interpretation, the regularization factor λ is the ratio between the observation and process noise variances, i.e.,

λ = δ 2 2 (13)
In [START_REF] Sameni | Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising[END_REF], it was shown that this smoothing scheme can be implemented as a forward-backward smoother for instance by using the filtfilt functions in Matlab, Octave or R.

In a blind or semi-blind source separation scheme, the smoothness order d and the noise variances δ 2 and 2 (and their ratio λ) can be unknown. Therefore, one should optimize the solution over the unknown parameter pair (d, λ).

Herein, we propose two approaches for optimized source separation over the unknown set (d, λ).

A. A JADE solution

In our first approach, we choose (by intuition) different combinations of integer values for d and positive real values λ k in different orders of magnitude, to cover the possible smoothness orders and signal to noise ratios. Next, we form the following matrices:

A k = XH(d k , λ k )H(d k , λ k ) T X T (k = 1, • • • , K) B = XX T (14
) and use the joint approximate diagonalization of eigenmatrices (JADE) to simultaneously diagonalize these matrices (approximately). The joint diagonalization of these matrix set, accounts for various orders of signal smoothness and unknown noise levels.

B. An iterative GEVD solution

In our second approach, following Algorithm 1 detailed in Section IV-B, successive stages of temporal and spatial filtering are used. The only part of Algorithm 1, which requires customization for different applications is to optimize the parametric temporal filter in each iteration. For illustration we fix the temporal filter smoothness order to d = 2, which only leaves the regularization factor λ to be optimized in each iteration.

In other contexts, methods such as the L-curve have been proposed for finding the optimal regularization factor [START_REF] Boyd | Convex optimization[END_REF], [START_REF] Gander | Least squares with a quadratic constraint[END_REF]- [START_REF] Calvetti | Tikhonov regularization and the l-curve for large discrete ill-posed problems[END_REF], which also has interpretations in terms of maximum a posteriori (MAP) estimation [START_REF] Oraintara | A method for choosing the regularization parameter in generalized tikhonov regularized linear inverse problems[END_REF].

In our case, using the notations of Algorithm 1, the L-curve is the plot of X -Y versus YD d , as λ is swept from zero to infinity. The λ corresponding to the corner of the Lcurve (in our case, the point at which X -Y + YD d is minimized) is empirically chosen as the optimal regularization factor, which provides the best balance between the signal smoothness and minimum mean square error. This optimal λ is considered as the parameter for the temporal filter H(d, λ) used in the next iteration of the algorithm.

VI. RESULTS

A potential domain of interest for the proposed methods is for extracting multichannel electrocardiogram (ECG) from background noise. In [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF], a semi-blind source separation scheme was developed for extracting maternal and fetal ECG using pseudo-periodicity priors. This method has shown to be very effective in practice; however, it highly relies on the pseudo-periodicity of the ECG and its performance degrades for arrhythmic and highly irregular ECG beats. On the other hand, as noted in [START_REF] Sameni | Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising[END_REF], due to the relatively high sampling frequency of contemporary ECG recording systems, the ECG (regular or irregular) is rather smooth as compared with the background noise. Therefore, "smoothness priors" can be used for designing ECG source separation schemes. Fig. 2. The DaISy dataset consisting of five maternal abdominal and three thoracic channels [START_REF] Moor | Database for the Identification of Systems (DaISy)[END_REF]; adopted from [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF]. Fig. 3. Independent components extracted from the dataset of Fig. 2, using JADE. Notice that components 1, 2, 3, and 5 correspond to the maternal subspace and components 4 and 8 to the fetal subspace; adopted from [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF].

A. Dataset

As proof of concept, we use the well-known DaISy fetal ECG dataset for illustration [START_REF] Moor | Database for the Identification of Systems (DaISy)[END_REF]. The dataset consists of five abdominal and three thoracic channels recorded from the abdomen and chest of a pregnant woman at a sampling rate of 250Hz. The eight channels of the dataset may be seen in Fig. 2.

We use independent subspace decomposition using the JADE algorithm as benchmark [START_REF] Cardoso | Multidimensional independent component analysis[END_REF], [START_REF] Cardoso | Blind beamforming for non Gaussian signals[END_REF]. The eight independent components extracted by this algorithm are depicted in Fig. 3.

To evaluate the hereby proposed method, we first use a single temporal smoothing matrix with d = 2 and λ = 1. The covariance of the original and smoothed data are jointly diagonalized using GEVD. The results obtained from this method are shown in Fig. 4, accordingly although the extracted components are smooth, only one fetal ECG channel has been extracted (in channel 4), and the second fetal component (channel 5) is diluted by the maternal peaks. This was already anticipated as the method has only targetted smooth components and not the separation of the fetal and maternal ECGs.

In the next stage we consider a more blind scenario. We first pre-whiten the data and form the smoothing matrix using multiple regularization factors {0.008, 0.04, 0.2, 1, 5, 25, 125}. In addition, we use the pseudo-periodicity matrices corresponding to the mother and the fetus, as described in the πCA algorithm [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF]. We therefore obtain nine covariance matrices, corresponding to pseudo-periodicity and temporal smoothness priors. The joint (approximate) diagonalization of these matrices [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF], [START_REF] Cardoso | Perturbation of joint diagonalizers[END_REF], accounts for pseudo-periodicity and an unknown noise level. The resulting signals are shown in Fig. 5. It is seen that this time the fetal and maternal components have been very well extracted.

VII. CONCLUSION AND FUTURE WORK

In this work, a spatio-temporal source separation scheme was proposed for the extraction of sources with prior temporal structure having parametric uncertainties in their temporal model.

The proposed method can be extended from various aspects in future studies, including:

• The hereby developed spatio-temporal filter is based on second order temporal statistics, which was shown to be very effective for signal extraction from noise. The diagonalized matrices used in this study can be used together with higher order statistics used in the classical ICA algorithms such as JADE [START_REF] Cardoso | Source separation using higher order moments[END_REF], to obtain simultaneous Fig. 5. Smooth components extracted from the dataset of Fig. 2, using JADE on the covariance matrix and several smoothed versions of the data using d = 2 and λ = {0.008, 0.04, 0.2, 1, 5, 25, 125}, plus two periodicity matrices obtained using πCA [START_REF] Sameni | Multichannel Electrocardiogram Decomposition using Periodic Component Analysis[END_REF]. Notice that algorithm has extracted smooth components and separated the fetal and maternal ECG. signal denoising (using spectral temporal priors) and independent components.

• Combining the proposed method with the deflation algorithm proposed in [START_REF]A Deflation Procedure for Subspace Decomposition[END_REF] (for rank deficient data).

• Combining the method with generalized multi-view analysis (GMA) [START_REF] Sharma | Generalized multiview analysis: A discriminative latent space[END_REF], for multiple modalities such as simultaneous ECG and PCG databases.

• One might be able to formulate the parameterized temporal filter using kernels with the same eigenstructures (sharing the same eigen-vectors) to guarantee that all matrices are exactly diagonalized using JADE.

• The problem of finding a closed form solution for the optimal temporal filter parameter should also be studied in future studies.

Algorithm 1

 1 Iterative spatio-temporal source separation Require: Multichannel signals X ∈ R N ×T , with zero-mean (per-channel) Require: Parameterized temporal filter matrix H(θ) ∈ R T ×T , with parameter set θ Require: Initial temporal filter parameter set θ 0 Require: Number of iterations I

Fig. 4 .

 4 Fig.[START_REF] Molgedey | Separation of a mixture of independent signals using time delayed correlations[END_REF]. Smooth components extracted from the dataset of Fig.2, using GEVD on the covariance matrix and a smoothed version of the data using d = 2 and λ = 1. Notice that algorithm has extracted smooth components; but has not separated the fetal and maternal ECG.

The kernel symmetry assumption, implies translational symmetry between the correlation functions of the two filtered channels, which is generally achieved in LTI systems. The asymmetric part of the cross-correlation can be of interest for time-variant causal analysis, which is beyond the scope of the current study as it takes the eigenvalue decomposition to the complex domain (due to the asymmetry of the diagonalized matrices).

The extension of this case study to signals with arbitrary spectral models is rather straightforward using, e.g., auto-regressive modeling; but beyond the scope of this work.


Reza Sameni
email: reza.sameni@gmail.com

Towards Distributed Component Analysis (working paper)

published or not. The documents may come

I. INTRODUCTION

In the past two decades, blind source separation (BSS) has been extensively used in various applied domains including audio, video, telecommunications and biomedical signal processing. While, the classical and most common formulation of BSS is based on the assumption of linear time-invariant mixtures of punctual sources, nonlinear and time-variant models have also found great attention in recent years [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF]. However, to the author's knowledge, in all these developments the sources of interest are assumed to be (or approximated by) finite number of spatially punctual sources. On the other hand, there are many applications, in which the punctuality of the sources does not hold or is not a good approximation. For example, in biomedical applications, sensors are typically rather close to the signal sources (heart, brain, or muscles), which violates the far-field assumptions. In these cases, the sources estimated by classical algorithms, such as principal component analysis (PCA) and independent component analysis (ICA) are only approximations of the true distributed activity. In [START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF], we introduced the notion of spatially distributed sources for electrocardiogram (ECG) signals and their impact on the signals extracted by applying conventional linear time-invariant PCA and ICA [START_REF] Sameni | What ICA Provides for ECG Processing: Application to Noninvasive Fetal ECG Extraction[END_REF]. It was empirically shown that components extracted from multichannel ECG (1) have infinite number of dimensions, i.e., always appear in full-rank regardless of the number of recordings 1 ; (2) the extracted components are not necessarily invariant under the change of electrode configuration. In this work, we attempt to justify these findings and extend them using a rather rigorous electrostatic framework for a spatially distributed stochastic electrostatic density function (used to model a distributed source signal). The current study is merely O ρ(u; t)

x 1

x 2

x 3

x n

x r u du In this report, a summary of the proposed framework is presented, which will be accompanied by further details and experimental results in the future versions of the manuscript.

II. DISTRIBUTED SOURCES

We consider a distributed source with a time varying stochastic density ρ(u; t) within a finite homogeneous volume (Fig. 1). The field potentials due to this distribution can be represented as follows:

φ(r; t) = ρ(u; t) |r -u| p du (1)
where r = [x, y, z] T is the Cartesian coordinates of the observation point 2 , p is the decay exponent which depends on the nature of ρ(u; t) and the field potentials, du is the differential volume unit, and the integral is taken over V , the volume containing the distributed source. Therefore, the potential difference of each measurement (observation) point (x i ), with respect to the reference can be written as follows:

∆φ(x i ; t) . = φ(x i ; t) -φ(x r ; t) = ρ(u; t) 1 |x i -u| p - 1 |x r -u| p du . = ρ(u; t)g(x i , u)du (2)
where g(x i , u) is the Green's function for the observation point x i , which acts as a distance-dependent weighting factor for ρ(u; t)3 .

Next, suppose that we make n observations of the field potentials from different locations, each time with respect to the same reference point x r . In this case, (2) should be rewritten in vector form. Defining x

. = [x 1 , x 2 , • • • , x n ] T , φ(x; t)
. = [∆φ(x 1 ; t), ..., ∆φ(x n ; t)] T and G(x, u) . = [g(x 1 , u), ..., g(x n , u)] T , (2) can be represented in the following form:

φ(x; t) = ρ(u; t)G(x, u)du (3)
where φ(x; t) can be considered as a vector of differential observations corresponding to, e.g., the signals recorded from the body or brain surface. Without loss of generality, these observations can be considered to be zero-mean (in time). The covariance matrix of the observations is therefore found as follows:

R φ (x) = E t {φ(x; t)φ(x; t) T } = G(x, u)E t {ρ(u; t)ρ(v; t)}G(x, v) T dudv = u)r ρ (u, v)G(x, v) T dudv
(4) where r ρ (u, v) . = E t {ρ(u; t)ρ(v; t)} is the cross-correlation function of the source density for different points inside the volume V , and E t {•} represents averaging over time.

Next, using the Karhunen-Loève Transform (KLT) [Van Trees, 2004, Ch. 3], r ρ (•, •) may be expanded in terms of a unique and (possibly) infinite set of orthogonal eigenfunctions {f i (u)} ∞ i=1 , i.e.:

λ i f i (u) = r ρ (u, v)f i (v)dv, (5)
f i (v)f j (v)dv = δ i-j (6)
where {λ i } ∞ i=1 are the eigenvalues and δ i-j is the Kronecker delta function. Following Mercer's theorem [Van Trees, 2004, Ch. 3], these eigenvalues and eigenfunctions can be used to expand r ρ (•, •) as follows:

r ρ (u, v) = ∞ i=1 λ i f i (u)f i (v) (7)
Combining ( 4) and ( 7), we find:

R φ (x) = ∞ i=1 λ i G(x, u)f i (u)f i (v)G(x, v) T dudv = ∞ i=1 λ i G(x, u)f i (u)du G(x, v) T f i (v)dv = ∞ i=1 λ i α i (x)α i (x) T (8) where G(x, s) = ∞ i=1 α i (x)f i (s) α i (x) . = G(x, s)f i (s)ds = G(x, •), f i (•)
i.e., α i (x) is the coefficient vector for the orthogonal expansion of the Green's function G(x, u) in terms of the eigenfunction f i (u). Apparently, α i (x)α i (x) T ∈ R n×n are symmetric rank-1 matrices with eigenvalues equal to α i (x) 2 . Therefore, R φ (x) is an infinite weighted sum of a set of positive definite rank-1 matrices. Note that although the eigenfunctions f i (•) are orthogonal, the α i (x) are not necessarily orthogonal. Therefore, ( 8) is a redundant decomposition of R φ (x) in the n-dimensional space.

On the other hand, the matrix R φ (x), being the covariance matrix of the observed signals, can be decomposed by PCA

R φ (x) = QDQ T = n i=1 d i q i q T i = ∞ i=1 λ i α i (x)α i (x) T (9)
where

Q = [q 1 , • • • , q n ] is the orthogonal eigenmatrix, D = diag(d 1 , • • • , d n )
contains the eigenvalues on its diagonal, and the q i q T i are orthogonal rank-1 projectors onto range(Q). From (9), we can see how the eigenvectors and eigenvalues of a set of observations φ(r; t) are related to the source distribution and electrode positioning. In fact, the set of eigenvectors {q i } n i=1 can be interpreted as an orthogonal basis for the subspace spanned by the column space of (the generally) nonorthogonal and infinite set of {α i (x)} ∞ i=1 . Equation ( 9) shows how the stochastic properties of the source density and the sensor configurations, both, influence the eigenstructure of the covariance matrix of the observations. From (9) it is also seen that the covariance matrix of measurements from distributed sources is generally full-rank, no matter how many observations are made; unless if the KLT expansion of r ρ (•, •) be finite (i.e., ∃l, ∀j > l : λ j = 0) or due to a specific sensor configuration the expansion of G(x, •) in terms of f i (•) be finite (i.e., ∃l, ∀j > l : α i (x) = 0). Either of these cases, can be due to special sensor configurations or special source densities. This explains many observations made from practical distributed sources, such as cardiac, muscular, or brain signals. In these cases it has been previously reported that the decomposition of the signals recorded from the body surface (the ECG, EEG, or the EMG), using PCA or ICA can lead into multiple components corresponding to the sources of interest [START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF].

III. SPECIAL CASES Some special cases of the proposed framework are explained in the followsing.

A. White Random Field

If ρ(u; t) is a white random field, we have:

r ρ (u, v) = σ(u) 2 δ(u -v)
where δ(•) is the Dirac delta function. In this case, (4) reduces to:

R φ (x) = σ(u) 2 G(x, u)G(x, u) T du (10)
which is generally a full-rank n × n matrix. A special case that reduces (10) to a finite summation of rank-1 matrices is:

σ(u) 2 = m i=1 σ 2 i δ(u -u i ) (11)
which corresponds to m uncorrelated punctual sources inside V . In this case, we have:

R φ (x) = m i=1 σ 2 i G(x, u i )G(x, u i ) T (12)
In this case, R φ (x) has a rank of min(m, n).

B. Synchronous (Coherent) Field

If the random field ρ(u; t) is synchronous throughout the entire volume V , we have r ρ (u, v) = σ 2 . In this case the double integral in ( 4) is separable and R φ (x) may be written as follows:

R φ (x) = σG(x, u)du σG(x, v)dv T (13)
which is a rank-1 matrix (the product of a column vector by its transpose). We therefore conclude that a fully synchronous distributed source behaves like a single punctual source, which is a well known result in the electromagnetics context. This result can also be extended to isolated groups of synchronous sources.

C. Dipolar Random Field

Dipolar random fields are another interesting special case that we study. For these sources we assume:

ρ(u; t) = f (t)δ(u -u 0 ) + g(t)δ(u + u 0 ) (14)
where f (t) and g(t) are zero-mean random processes. Equation ( 14) leads to:

r ρ (u, v) = E t {ρ(u; t)ρ(v; t)} = σ 2 f δ(u -u 0 )δ(v -u 0 ) + σ 2 g δ(u + u 0 )δ(v + u 0 ) +σ f g δ(u -u 0 )δ(v + u 0 ) + σ f g δ(u + u 0 )δ(v -u 0 )
(15) where σ 2 f , σ 2 g , and σ f g , correspond to the variances of f (t), g(t), and their cross-correlations, respectively. Combining ( 15) and ( 4) we find:

R φ (x) = σ 2 f G(x, u 0 )G(x, u 0 ) T + σ 2 g G(x, -u 0 )G(x, -u 0 ) T +σ f g G(x, u 0 )G(x, -u 0 ) T + σ f g G(x, -u 0 )G(x, u 0 ) T (16)
This equation can be rearranged as the summation of three rank-1 matrices as follows:

R φ (x) = (σ 2 f + σ f g )G(x, u 0 )G(x, u 0 ) T +(σ 2 g + σ f g )G(x, -u 0 )G(x, -u 0 ) T -σ f g [G(x, u 0 ) -G(-u 0 )][G(x, u 0 ) -G(x, -u 0 )] T (17)
which has a maximum rank of three. Apparently, if the two punctual sources are uncorrelated, i.e., σ f g = 0, R φ (x) will have a maximum rank of two, which reduces to the special case studied in the previous subsection.

D. Stationary Random Field

Another special case that is rather practical is a stationary random field. For a covariance stationary random field we have [Van Trees, 2004, Ch. 3]:

r ρ (u, v) = r ρ (u -v) . = r ρ (τ ) (18)
In this case, ( 5) can be solved using classical methods of spectral decomposition [Van Trees, 2004, Ch. 3], and the λ i in ( 5) are related to the contribution of the different frequency components in the spectrum of r ρ (τ ). Examples of practical covariance stationary random fields include exponentially decaying cross-correlation functions of the form:

r ρ (τ ) = A exp(- |τ | q σ ) , q = 1, 2 (19)
which can be used to model correlation functions that decay with the increase of distance between points. The eigenvalues and eigenfunctions of such correlation functions can be found in [START_REF] Van Trees | Detection, Estimation, and Modulation Theory. Detection, Estimation, and Modulation Theory[END_REF][START_REF] Roger | Stochastic finite elements: a spectral approach[END_REF][START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF].

IV. EXTENSION TO NON-HOMOGENEOUS MEDIA

The basic equation ( 1) is the solution of Poisson's equation in homogeneous media:

∇ • [-∇σ(r)φ(r; t)] = ρ(r; t) (20)
The general problem can be studied using this equation, which enables us to consider conductivity mismatch between boundaries of brain layers. Refer to [Nunez and Srinivasan, 2006, P. 60] for the background. One can also use Green's function and related identities to convert surface integrals to volume integrals.

V. GENERALIZATION AND FUTURE WORK

In our recent research, the hereby presented framework has been extended in various aspects, including spatially stationary versus nonstationary sources and extension to nonhomogeneous media. We have developed a simulation platform for generating synthetic signals with desired distributions and a given R φ (x) for EEG signals [START_REF] Tavakkol-Shoorijeh | Distributed Component Analysis and its Applications in Biosignal Processing[END_REF]. Conventional PCA and ICA have been applied to these synthetic distributed signals and we show that the results are consistent with the ones obtained from real data (specifically ECG and EEG).

The findings of this research are expected to be helpful for better interpretation of ICA results, specifically for biomedical applications. We can show how reliable the components extracted by PCA and ICA be, when the sources are not punctual. The presented approach is hoped to be extended to higher order statistics (used in some ICA algorithms), and lead to the development of a general framework for studying spatially distributed sources, which we coin distributed component analysis (DCA).

Following this work, some of the major questions that should be answered are as follows:

1) An interesting question is what the sources extracted by conventional ICA and PCA correspond to. Our formulations show that they are projections of the distributed source onto specific basis functions; but the exact basis functions need to be derived. 2) Another important item that should be verified, is to show that for distributed sources, the components extracted by conventional PCA and ICA generally change with the location of the sensors (unlike conventional ICA). Do you know any good reference that shows this for conventional ICA (except for the variance and order of sources)? Perhaps, one approach can be to show that the mutual information of the data changes with the change of electrode positions for distributed sources. 3) Following the previous item, an interesting question is whether we can find sensor position-invariant sources corresponding to distributed sources. Note that in this approach we are not interested in independent or uncorrelated components. We just want to find similar components using any electrode configuration. This is very similar to the idea that we had in [START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF] for canonical representations of the ECG. 4) It is believed that similar, but more complicated, equations can be derived for cumulant tensors [Hyvärinen et al., 2001, Ch. 11], of φ(t). These tensors can be used to relate ICA solutions based on higher-order statistics to the aforementioned physical properties of the distributed source. 5) In the derivation of the hereby presented equations, we made no specific assumptions on the nature of the source density ρ(u; t). For the case of ECG signals, taking p = 1 (in (1) ), this density function can represent ∇ • J, the divergence of the cardiac impressed current dipole moment density defined in [Sameni, 2008, Ch. 3]. 6) In future studies, a more realistic formulation is to consider non-homogeneous propagation media, which as mentioned in [Sameni, 2008, Ch. 3], leads to some sort of 'reflection' sources due to the surfaces of conductivity discontinuity. By using Green's function and analytic identities to convert surface integrals to volume integrals it should be possible to study non-homogeneous distributions with similar formulations. 7) The hereby adopted framework was based on the integral form of Poisson's equation. In future works, one can analyse the problem from a convolution (filtering) perspective, which may lead to a filter bank formulation of the problem. 8) Overall, what is the interpretation of the components extracted by conventional ICA over distributed sources?

VI. FURTHER READING 1) An excellent book: [START_REF] Roger | Stochastic finite elements: a spectral approach[END_REF].

2) [Van Trees, 2004, Ch. 3] has some related formulations.

3) See [START_REF] Arfken | Mathematical Methods For Physicists International Student Edition[END_REF] for the nonhomogeneous case using Green's functions expansions. 4) The idea of clouds of sources: http://sccn.ucsd.edu/ ∼ arno/indexica.html 5) [START_REF] Lathauwer | Fetal electrocardiogram extraction by blind source subspace separation[END_REF] has referred to the model; but does not study components higher than 3 and has not given physiological interpretations for the extracted components. 6) See [Bronzino, 2000, Ch. 16] for a propagation model.

Fig. 1 .

 1 Fig. 1. Distributed source model and the electrode positions

For simplicity, we have normalized the equation by the conductivity of the medium, as it does not affect our calculations.

Note that g(x i , u) is also a function of the reference point xr, which we do not show for notational simplicity, as the reference is a fixed point for all the observation points.


Constantin Bacuta
email: bacuta@math.udel.edu

A UNIFIED APPROACH FOR UZAWA ALGORITHMS *

Keywords: Uzawa algorithms, saddle point system, multilevel methods, augmented Lagrangian method, Stokes problem AMS subject classifications. 74S05, 74B05, 65N22, 65N55

We present a unified approach in analyzing Uzawa iterative algorithms for saddle point problems. We study the classical Uzawa method, the augmented Lagrangian method, and two versions of inexact Uzawa algorithms. The target application is the Stokes system, but other saddle point systems, e.g., arising from mortar methods or Lagrange multipliers methods, can benefit from our study. We prove convergence of Uzawa algorithms and find optimal rates of convergence in an abstract setting on finite-or infinite-dimensional Hilbert spaces. The results can be used to design multilevel or adaptive algorithms for solving saddle point problems. The discrete spaces do not have to satisfy the LBB stability condition.

For f ∈ V * , g ∈ P * , we consider the following variational problem: Find (u, p) ∈ V × P such that

a(u, v) + b(v, p) = f , v for all v ∈ V, b(u, q) =( g, q)
for all q ∈ P. (1.3) It is known that the above variational problem has a unique solution for any f ∈ V * , g ∈ P * (see [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] or Lemma 2.1). With the forms a and b, we associate two linear operators A : V → V * and B : V → P defined by Au, v = a(u, v) for all u, v ∈ V and (Bu, q) = b(u, q) for all u ∈ V, q ∈ P.

Let B * : P → V * be the dual operator of B defined by B * q, v = (q, Bv) = (Bv, q) = b(v, q) for all v ∈ V, q ∈ P.

The problem (1.3) is equivalent to the following problem:

Find (u, p) ∈ V × P such that Au + B * p = f , Bu = g. (1.4) In this framework, we analyze Uzawa algorithms for solving the system (1.3) or (1.4). We consider that the form a gives the inner product and the norm on V. A more general case of (1.3) is considered in [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]. Our particular assumptions for the form a give rise to a simplified analysis. For the general case, we obtain sharp convergence estimates only in terms of the two constants m and M .

The rest of the paper is organized as follows. In section 2, we analyze the convergence of the classical Uzawa algorithm. The augmented Lagrangian Uzawa method is analyzed in section 3 (Fortin and Glowinski [14]). In section 4, we shall investigate the convergence of the inexact Uzawa algorithm (Bramble, Pasciak, and Vassilev [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF] and Elman and Golub [START_REF] Elman | Inexact and preconditioned Uzawa algorithms for saddle point problems[END_REF]) in the above abstract framework. Applications to discretizations on stable pairs are presented in section 5. A modified inexact Uzawa algorithm with applications in constructing multilevel methods and adaptive methods for solving (1.3) is illustrated in section 6. In section 7, we present applications of our abstract results to the Stokes system.

2. The abstract Uzawa algorithm. We begin this section with two lemmas which provide basic properties of norms and operators introduced in section 1. The proofs are based on the Riesz representation theorem (see, e.g., [START_REF]Yosida Functional Analysis[END_REF]). For completeness, we include the proofs.

Lemma 2.1. The operator A : V → V * is invertible and the Schur complement operator BA -1 B * : P → P is symmetric and a positive definite operator satisfying

(BA -1 B * p, p) = sup v∈V b(v, p) 2 |v| 2 , (2.1)
m 2 p 2 ≤ (BA -1 B * p, p) ≤ M 2 p 2 , p ∈ P. (2.2)
Consequently, the problem (1.3) (or (1.4)) has a unique solution.

Proof. From the definition of A, we get that A is a bounded injective operator. Using the Riesz representation theorem, it follows that A is also a surjective operator. Let us further note that A satisfies

Au, v = a(u, v) = a(v, u) = Av, u ,
and the changes of variable Au = u * and Av = v * lead to

u * , A -1 v * = v * , A -1 u * , u * , v * ∈ V * . (2.3) Using (2.3), we obtain (BA -1 B * p, q) = B * q, A -1 B * p = B * p, A -1 B * q = (BA -1 B * q, p) = (p, BA -1 B * q), p, q ∈ P.
To prove (2.1), we let p ∈ P be fixed and consider the following problem:

Find u ∈ V such that a(u, v) = b(v, p) for all v ∈ V. (2.4)
Since the functional v → b(v, p) is continuous on V, by the Riesz representation theorem we have that the unique solution u of (2.4) satisfies

a(u, u) = v → b(v, p) 2 V * = sup v∈V b(v, p) 2 |v| 2 . (2.5)
On the other hand, from (2.4) we have

Au = B * p or u = A -1 B * p and a(u, u) = Au, u = B * p, A -1 B * p = (p, BA -1 B * p). (2.6)
Thus, (2.1) follows from (2.5) and (2.6). The estimate (2.2) follows imediately from (2.1), (1.1), and (1.2).

To prove the existence and uniqueness of (1.3) (or (1.4)), we substitute u from the first equation of (1.4) into the second equation of (1.4). The resulting equation in p, BA -1 B * p = BA -1 fg, has a unique solution due to the fact that BA -1 B * : P → P is symmetric and a positive definite operator.

Remark 2.2. From the general theory of symmetric operators and Lemma 2.1, we have that σ(BA -1 B * ) ⊂ [m 2 , M 2 ] and m 2 , M 2 ∈ σ(BA -1 B * ). In the finitedimensional case, m 2 and M 2 are the extreme eigenvalues of the Schur complement BA -1 B * .

Lemma 2.3. The following norm estimates are valid:

φ 2 V * = a(A -1 φ, A -1 φ) = |A -1 φ| 2 , φ ∈ V * , (2.7) Au V * = |u|, u ∈ V, (2.8) B * q V * = |A -1 B * q| = (BA -1 B * q, q) 1/2 ≤ M q , q ∈ P, (2.9) B = M, hence Bu ≤ M |u|, u ∈ V. (2.10)
Proof. By the Riesz representation theorem, we have that for any φ ∈ V * the problem Find u ∈ V such that

Au, v = a(u, v) = φ, v for all v ∈ V (2.11)
has a unique solution, and the solution u satisfies

a(u, u) = sup v∈V φ, v 2 a(v, v) = φ 2 V * . (2.12)
From (2.11), we have that u = A -1 φ, which combined with (2.12) gives (2.7). The equality (2.8) is a consequence of (2.7), and (2.9) follows from (2.8) and (2.2). The last estimate follows from the definition of B and the assumption in (1.2).

Next, we present the Uzawa algorithm [START_REF] Arrow | Studies in Nonlinear Programming[END_REF] for solving the solution of the abstract problem (1.3). Given a parameter α > 0, called a relaxation parameter, the Uzawa algorithm for approximating the solution (u, p) of (1.3) can be described as follows.

Algorithm 2.4 (Uzawa method (UM)). Let p 0 be any approximation for p, and for k = 1, 2, . . . , construct (u k , p k ) by

a(u k , v) = (f , v) -b(v, p k-1 ), v ∈ V, p k = p k-1 + α(Bu k -g). (2.13)
The convergence of the UM is discussed for particular cases in, e.g., [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solutions of Boundary Value Problems[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF][START_REF] Temam | Navier-Stokes Equations[END_REF]. It shows that the UM is convergent for small enough α and that the convergence rate is the same as the convergence rate of the Richardson iterative methods for the Schur complement BA -1 B * . For completeness, we include the proof.

Theorem 2.5. Let (u, p) be the solution of (1.3) and let (u k , p k ) be the sequence of approximations built by the UM (2.13). Then, the following holds.

(i) The sequences uu k and pp k satisfy

a(u -u k , u -u k ) 1/2 ≤ M p -p k-1 , p -p k ≤ I -αBA -1 B * p -p k-1 .
(ii) For α < 2 M 2 , the UM is convergent and

I -αBA -1 B * = max{|1 -αm 2 |, |1 -αM 2 |} < 1. (iii) For α = 1 M 2 , the convergence factor is I -αBA -1 B * = 1 -m 2 M 2 . (iv)
The optimal convergence factor is achieved for

α opt = 2 M 2 + m 2 and I -α opt BA -1 B * = M 2 -m 2 M 2 + m 2 .
Proof. From the first equation of (1.3) and the first equation of (2.13), we have that

a(u -u k , v) = b(v, p k-1 -p)
for all v ∈ V. (2.14) The above relation implies

|u -u k | 2 = a(u -u k , u -u k ) = (BA -1 B * (p k-1 -p), (p k-1 -p)) ≤ M 2 p k-1 -p 2 ,
which proves the first part of (i). From the second equation of (1.4) and the second equation of (2.13), we have that

p -p k = p -p k-1 + αB(u -u k ).
Combining with (2.14), we get

p -p k = (I -αBA -1 B * )(p -p k-1 ), (2.15)
which gives the second part of (i). From Lemma 2.1, we have that (I -αBA -1 B * ) is a symmetric operator, and for any p ∈P, p = 0,

1 -αM 2 ≤ ((I -αBA -1 B * )p, p) p 2 ≤ 1 -αm 2 ,
which justifies part (ii). The rest of the proof follows from (ii).

Augmented Lagrangian Uzawa algorithm.

The main idea of the augmented Lagrangian method, introduced by Fortin and Glowinski [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solutions of Boundary Value Problems[END_REF], is to use the constraint condition for the variable p and another tuning parameter ρ > 0 in order to improve the convergence factor of the Uzawa algorithm. We will consider the approach for abstract Hilbert spaces V and P and prove sharp convergence estimates for the corresponding Uzawa algorithm.

Let (u, p) be the solution of the variational problem (1.3). Then, from the second equation of (1.4), we have that

(Bu, Bv) = (g, Bv), v ∈ V.
Thus, for any ρ > 0, (u, p) is also a solution of

a(u, v) + ρ(Bu, Bv) + b(v, p) = f , v + ρ(g, Bv), b(u, q) =( g, q). (3.1)
Using the notation a ρ (u, v) := a(u, v) + ρ(Bu, Bv) and f ρ := f + ρB * g, we have that

a ρ (u, v) + b(v, p) = f ρ , v , v ∈ V, b(u, q) =( g, q), q ∈ P. (3.2)
With the form a ρ , we associate the linear operator

A ρ : V → V * , A ρ u, v = a ρ (u, v)
for all u, v ∈ V.

Thus, an equivalent form of (3.2) is In what follows, the Uzawa algorithm applied to (3.2) will be called the augmented Lagrangian Uzawa method (ALUM).

A ρ u + B * p = f ρ , Bu = g. (3.3) Since a ρ (•,
Given a relaxation parameter α > 0, the augmented Lagrangian Uzawa algorithm for approximating the solution (u, p) of (1.3) is as follows.

Algorithm 3.1 (ALUM). Let p 0 be any approximation for p, and for

k = 1, 2, . . . , construct (u k , p k ) by a ρ (u k , v) = (f ρ , v) -b(v, p k-1 ), v ∈ V, p k = p k-1 + α(Bu k -g).
To study the convergence of (3.1), we shall calculate first

M ρ := sup p∈P sup v∈V b(v, p) p (a ρ (v, v)) 1/2 (3.4)
and

m ρ := inf p∈P sup v∈V b(v, p) p (a ρ (v, v)) 1/2 . (3.5)
Theorem 3.2. For any ρ > 0, we have

BA -1 ρ B * = ρI + (BA -1 B * ) -1 -1 , (3.6) M 2 ρ = 1 ρ + 1 M 2 and m 2 ρ = 1 ρ + 1 m 2 . (3.7)
Proof. To prove (3.6), we need two identities. First, we note that for any invertible linear operator C : P → P such that I + ρC is also invertible, we have

ρI + C -1 -1 = C -ρC(I + ρC) -1 C. (3.8)
This can be proved by checking that the proposed inverse verifies the algebraic definition of the inverse. The second identity is based on the Sherman-Morrison-Woodbury formula and can be proved again just by algebraic manipulations:

(A + ρB * B) -1 = A -1 -ρA -1 B * (I + ρBA -1 B * ) -1 BA -1 . (3.9)
From (3.9), we get

B(A + ρB * B) -1 B * = BA -1 B * -ρBA -1 B * (I + ρBA -1 B * ) -1 BA -1 B * . (3.10)
If we take C = BA -1 B * in (3.8) and combine it with (3.10), we obtain (3.6). To verify (3.7), we notice that by applying Lemma 2.1 with a ρ instead of a we have

sup v∈V b(v, p) 2 p 2 a ρ (v, v) = (BA -1 ρ B * p, p), p ∈ P. (3.11)
Thus, we get

M 2 ρ = sup p∈P (BA -1 ρ B * p, p) (p, p) = sup p∈P (ρI + (BA -1 B * ) -1 ) -1 p, p (p, p) = 1 inf q∈P ((ρI+(BA -1 B * ) -1 )q,q) (q,q) = 1 ρ + inf q∈P ((BA -1 B * ) -1 q,q) (q,q) = ⎛ ⎜ ⎝ ρ + 1 sup r∈P (BA -1 B * r,r) (r,r) ⎞ ⎟ ⎠ -1 = ρ + 1 M 2 -1
.

Here, we have used the changes of variable (ρI + (BA -1 B * ) -1 ) -1/2 p = q and (BA -1 B * ) -1/2 q = r. The proof for m ρ is similar.

The above result gives formulas for the inf-sup and sup-sup constants for the ALUM in terms of m, M , and ρ. In applications, the constant m is more difficult to obtain. The following theorem gives the convergence rate of the ALUM.

Theorem 3.3. Let (u, p) be the solution of (1.3) and let (u k , p k ) be the sequence of approximations built by Algorithm 3.1. Then, the following holds true:

(i) The sequences u -u k and p -p k satisfy a ρ (u -u k , u -u k ) 1/2 ≤ M ρ p -p k-1 , p -p k ≤ I -αBA -1 ρ B * p -p k-1 . (ii) For α < 2 M 2 ρ
, the ALUM is convergent and

I -αBA -1 ρ B * = max{|1 -αm 2 ρ |, |1 -αM 2 ρ |} < 1. (iii) For α = 1 M 2 ρ
, the convergence factor is

I -αBA -1 ρ B * = 1 - m 2 M 2 1 m 2 ρ + 1 .
(iv) The optimal convergence factor is achieved for

α opt = 2 M 2 ρ +m 2 ρ and I -α opt BA -1 ρ B * = M 2 ρ -m 2 ρ M 2 ρ + m 2 ρ = 1 - m 2 M 2 1 2m 2 ρ + 1 + m 2 /M 2 .
Proof. The result is a direct consequence of Theorem 2.5 and (3.7). A similar result for the discrete version of the Stokes system can be found in [START_REF] Xu | Multilevel Finite Element Theory[END_REF]. As it was pointed out in [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solutions of Boundary Value Problems[END_REF] and [START_REF] Xu | Multilevel Finite Element Theory[END_REF], the choice of a very large ρ improves on the rate of convergence of the ALUM, but at the same time, the operator A ρ becomes more difficult to invert. For the continuous and discrete Stokes system, estimates for the convergence factor of the ALUM were recently obtained by Nochetto and Pyo in [START_REF] Nochetto | Optimal relaxation parameter for the Uzawa method[END_REF]. The question raised in [START_REF] Nochetto | Optimal relaxation parameter for the Uzawa method[END_REF] on how much we can improve the rate of convergence of the ALUM if information about the spectral value m is available can be easily answered now by comparing part (iii) and part (iv) of Theorem 3.3 or Theorem 7.1.

Inexact Uzawa method.

Throughout the rest of the paper we will keep the notation and assumptions of section 2. In this section, following the ideas in [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF][START_REF] Elman | Inexact and preconditioned Uzawa algorithms for saddle point problems[END_REF], we shall introduce and investigate the convergence of an abstract inexact Uzawa algorithm where the exact solve of the elliptic problem (the action of A -1 ) is replaced by an approximation process, which might not be a linear operator. We describe the approximate inverse of A as a map C :

V * → V which, for φ ∈ V * , returns an approximation of ξ = A -1 φ such that |Cφ -A -1 φ| V ≤ δ φ V * for all φ ∈ V * (4.1)
for some δ ∈ (0, 1). We notice here that (4.1) is a strong condition for the infinitedimensional case. The condition can be weakened by requiring to be satisfied only for certain values φ ∈ V * . If V and P are finite-dimensional spaces, then C can be considered as a linear or nonlinear process for inverting A and (4.1) is a reasonable assumption (see [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF]). One example of nonlinear process C is the approximate inverse associated with the preconditioned conjugate gradient algorithm. A practical case would be to consider Cφ = ξ num , where ξ num is the numerical approximation of ξ defined by

a(ξ, v) = φ, v for all v ∈ V.
In any case, if Aξ = φ and Cφ is defined by Cφ = ξ ap , an approximation of ξ, then, according to (2.7), the assumption (4.1) is equivalent to

|ξ ap -ξ| V ≤ δ|ξ| V for all ξ ∈ V. (4.2)
The inexact Uzawa algorithm for approximating the solution (u, p) of (1.3) is as follows.

Algorithm 4.1 (inexact Uzawa method (IUM)). Let (u 0 , p 0 ) be any approximation for (u, p), and for k = 1, 2, . . . , construct (u k , p k ) by

u k = u k-1 + C(f -Au k-1 -B * p k-1 ), p k = p k-1 + α(Bu k -g).
Before we study the stability and convergence rate of Algorithm 4.1 we shall introduce the following notation. For k = 0, 1, . . . , let e u k = uu k , e p k = pp k , and

E k = |e u k | e p k .
Let

M := δ M(1 + δ) αM δ γ + αM 2 δ
,

where γ := I -αBA -1 B * = max{|1 -αm 2 |, |1 -αM 2 |}. On R 2 we introduce the inner product [•, •] w defined by x 1 x 2 , y 1 y 2 w = w 1 x 1 y 1 + w 2 x 2 y 2 ,
where w 1 , w 2 are any two positive numbers such that

w 1 w 2 = αδ 1 + δ ,
and δ is a positive number such that (4.1) is satisfied. We note that M is symmetric with respect to the [•, •] w inner product. We will denote the norm induced by

with δ < 1 -γ 1 -γ + 2αM 2 . (4.3)
Then, the IUM converges. If r is the spectral radius of the matrix M, then 0 < r < 1 and

E k w ≤ r k E 0 w , k = 1, 2, . . . . (4.4)
Proof. We follow the proof of a similar result in [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF] for the finite-dimensional case. From the first equation of (1.4) and the first equation of Algorithm 4.1, we have

e u k = e u k-1 -C(Ae u k-1 + B * e p k-1 ) = (A -1 -C)(Ae u k-1 + B * e p k-1 ) -A -1 B * e p k-1 . (4.5)
From the second equation of (1.4) and the second equation of Algorithm 4.1, we get

e p k = e p k-1 + αBe u k . (4.6)
If we substitute e u k from (4.5) into (4.6), then

e p k = (I -αBA -1 B * )e p k-1 + αB(A -1 -C)(Ae u k-1 + B * e p k-1 ). (4.7)
From (4.5) and (4.7), by the triangle inequality, and from the estimates (2.9) and (2.10) and the assumption (4.1), we obtain

|e u k | ≤ δ |e u k-1 | + M (1 + δ) e p k-1
and

e p k ≤ αM δ |e u k-1 | + (γ + αM 2 δ) e p k-1 .
Using the notation introduced above, we have

E k ≤ M E k-1 , (4.8)
where

x 1 x 2 ≤ y 1 y 2
means x 1 ≤ y 1 and x 2 ≤ y 2 . From (4.8), we deduce

E k ≤ M k E 0 . (4.9)
Since M is symmetric with respect to [•, •] w -inner product, we have

E k 2 w = [E k , E k ] w ≤ [M k E 0 , M k E 0 ] w = [M 2k E 0 , E 0 ] w ≤ r 2k E 0 2 w ,
which proves (4.2). To complete the proof, we have to show that r ∈ (0, 1), provided that 0 < α < 2/M 2 and (4.3) holds. The characteristic equation of the matrix M is

λ 2 -λ(δ + γ + αM 2 δ) + δ(γ -αM 2 ) = 0.
Since M has positive entries, the characteristic equation has real roots and the largest (positive) root agrees with the spectral radius of M. Consequently,

r = 1 2 δ + γ + αM 2 δ + (δ + γ + αM 2 δ) 2 -4δ(γ -αM 2 ) .
Using that γ = max{|1αm 2 |, |1 -αM 2 |} and α ∈ (0, 2/M 2 ), it is easy to verify that the function δ → r = r(δ) is an increasing function on (0, 1) and that r = 1 for

δ = δ 0 := 1 -γ 1 -γ + 2αM 2 . (4.10)
This completes the proof of the theorem.

Remark 4.3. For 0 < α ≤ 2 M 2 +m 2 we have that γ = 1αm 2 and the threshold δ 0 becomes

δ 0 = m 2 m 2 + 2M 2 , which is independent of α. For 2 M 2 +m 2 ≤ α < 2
M 2 we have that γ = αM 2 -1 and the threshold δ 0 becomes

δ 0 = 2 -αM 2 2 + αM 2 .
Nevertheless, the optimal (maximal) value of δ 0 as the function of

α ∈ [ 2 M 2 +m 2 , 2 M 2 ) is δ 0 = m 2 m 2 +2M 2 and is achieved for α = 2 M 2 +m 2 .
Thus, a good choice for α (independent of m) is α = 1/M 2 . In this case we still have δ 0 = m 2 m 2 +2M 2 . Remark 4.4. We can apply the IUM for the augmented Lagrangian formulation. The only changes in Algorithm 4.1 is that A is replaced by A ρ and f is replaced by f ρ . The convergence analysis follows from Theorem 4.2. Let us further notice that in this case |e u k | 2 = a ρ (e u k , e u k ) and for α = 1/M 2 ρ the threshold δ 0 which assures convergence for the IUM is

δ 0 (ρ) = m 2 ρ m 2 ρ + 2M 2 ρ = m 2 + ρm 2 M 2 m 2 + 2M 2 + 3ρm 2 M 2 → 1 3 as ρ → ∞.
Thus, if the IUM for the augmented Lagrangian formulation is applied with sufficiently large ρ, α = 1/M 2 ρ and with the approximation operator C satisfying

C -A -1 ρ ≤ δ 0 (ρ) < 1/3,
then the method converges. Remark 4.5. A different approach in analyzing the IUM in the finite-dimensional case is presented by Cheng in [START_REF] Cheng | On the nonlinear inexact Uzawa algorithm for saddle-point problems[END_REF]. From his analysis for α = 1 and M = 1, it follows that the IUM converges (with a different estimate for the convergence factor), under the weaker assumption that δ < δ 0 = 1/3. Cheng's result for the infinite-dimensional case seems not to have been investigated. A positive answer for this problem would be an interesting result, since in practice it is difficult to estimate the spectral value m.

Discretization with the inf-sup condition.

In this section we assume that the variational form of a PDE (or system of PDEs) leads to (1.3) and let V h and P h be two finite-dimensional spaces, V h ⊂ V, P h ⊂ P , with good approximation properties. We further assume that inf

p∈P h sup v∈V h b(v, p) p |v| = m(h) > 0 and sup p∈P h sup v∈V h b(v, p) p |v| = M (h). (5.1)
For an overview of numerical methods for solving saddle point systems, we refer the reader to the recently published review paper [START_REF] Benzi | Numerical solutions of saddle point problems[END_REF] by Benzi, Golub, and Liesen.

From Lemma 2.1 and Remark 2.2 we see that m(h), M(h) are the lowest and the largest eigenvalues of the Schur complement B h A -1 h B * h associated with the discrete spaces V h and P h . Then (see, e.g., [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]), the discrete variational problem

Find (u h , p h ) ∈ V h × P h such that a(u h , v) + b(v, p h ) = f , v , v ∈ V h , b(u h , q) =( g, q), q ∈ P h , (5.2)
has a unique solution and

|u -u h | + p -p h ≤ c inf v∈V h |u -v| + inf q∈P h p -q ,
where c is a constant depending only on m(h) and M (h). In this case, the exact or inexact Uzawa algorithms can be applied for the discrete variational problem (5.2) on V h × P h ; see, e.g., [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF]. The convergence factors depend on m(h) and M (h) and could deteriorate as h → 0 if the pair (V h , P h ) is not stable. We recall here that a pair (V h , P h ), or more precisely a family of pairs

{(V h , P h )} h , is called stable if m(h) defined in (5.1) satisfies m(h) ≥ m > 0,
with m independent of h. In the next section we use the inexact Uzawa algorithm at the continuous level to construct algorithms which avoid building stable pairs (V h , P h ).

Modified inexact Uzawa method. Eliminating the discrete inf-sup condition.

We shall apply the IUM to construct discrete approximations (u k , p k ) ∈ (V k , P k ), where V k ⊂ V and P k ⊂ P are finite-dimensional spaces such that the pairs (V k , P k ) do not have to be stable pairs.

The algorithm proposed in this section can be used for building multilevel or adaptive methods for solving the system (1.3). Adaptive methods for saddle point problems have been the subject for recent research in numerical analysis (see, e.g., [START_REF] Dahlke | Adaptive wavelet methods for saddle point problemsoptimal convergence rates[END_REF][START_REF] Bänsch | An adaptive Uzawa FEM for the Stokes problem: Convergence without the inf-sup condition[END_REF]). Our new approach, combined with standard techniques of a posteriori error estimate theory, could lead to new and efficient adaptive algorithms for solving saddle point systems. To describe our new algorithm, we assume that a sequence of nested subspaces,

V 0 ⊂ V 1 ⊂ V 2 ⊂ • • • ⊂ V,
was determined and for k = 1, 2, . . . , a linear or nonlinear process

C k : V * → V k approximating A -1 is available such that for a fixed φ ∈ V * , C k φ ∈ V k is an ap- proximation of ξ = A -1 φ. To construct a good approximate inverse C k : V * → V k
one might need to increase the space V k-1 to a space with better approximation properties using an adaptive method. Thus, the embedding assumption V k-1 ⊂ V k is needed. On the other hand, in the proposed algorithms, the variable p is updated at the continuous level and no inversion is used. Thus, the P k 's are just subsets of the space P and do not have to be nested.

The modified inexact Uzawa algorithm for approximating the solution (u, p) of (1.3) can be stated now as follows.

Algorithm 6.1 (modified inexact Uzawa method (MIUM)). Let u 0 ∈ V 0 be any approximation for u and let p 0 ∈ P be any approximation for p. For k = 1, 2, . . . , construct (u k , p k ), with u k ∈ V k , by

u k = u k-1 + C k (f -Au k-1 -B * p k-1 ), p k = p k-1 + α(Bu k -g). (6.1) Theorem 6.2. Let 0 < α < 2/M 2 , γ = max{|1 -αm 2 |, |1 -αM 2 |} = I - αBA -1 B * , and assume that for k = 1, 2, . . . , C k satisfies (C k -A -1 )(f -Au k-1 -B * p k-1 ) V ≤ δ (f -Au k-1 -B * p k-1 ) V * , (6.2) with δ < 1 -γ 1 -γ + 2αM 2 . (6.3)
Then, the MIUM converges and the convergence rate is given by (4.4).

Proof. It is similar to the proof of Theorem 4.2.

We notice here that, for a fixed α, the threshold δ which assures the convergence of the MIUM depends only on the constants m and M . In the case g = 0, we have p k ∈ P k := BV k . Nevertheless, no matter the choice of the spaces V k , P k , the pairs (V k , P k ) do not have to be stable pairs.

For the rest of this section, the first equation in (6.1) will be considered in a variational form as follows. Let d k ∈ V k be the solution of

a(d k , v) = f, v -a(u k-1 , v) -b(v, p k-1 ), v ∈ V k . (6.4) Take dk := C k (f -Au k-1 -B * p k-1
) to be an approximation of d k . For example, dk could be a numerical approximation of d k . Let us assume that D k-1 ∈ V is the solution of the continuous problem

a(D k-1 , v) = f, v -a(u k-1 , v) -b(v, p k-1 ), v ∈ V. (6.5) From the Riesz representation theorem (f -Au k-1 -B * p k-1 ) V * = |D k-1 | V .
Thus, the assumption (6.2) can be rewritten as

| dk -D k-1 | V ≤ δ |D k-1 | V . (6.6) Since d k ∈ V k is the Galerkin approximation of D k-1 ∈ V, we have that |d k | V ≤ |D k-1 | V . A sufficient condition for the assumption (6.2) is | dk -D k-1 | V ≤ δ |d k | V . (6.7)
6.1. Multilevel exact Uzawa. In this subsection we assume that the problem (6.4) can be solved exactly on V k , i.e., dk = d k . Then,

u k = u k-1 + d k and consequently, a(u k , v) = f, v -b(v, p k-1 ), v ∈ V k . If U k-1 ∈ V satisfies a(U k-1 , v) = f, v -b(v, p k-1 ), v ∈ V, then D k-1 = U k-1 -u k-1 and (6.6) is equivalent to |u k -U k-1 | V ≤ δ |u k-1 -U k-1 | V . (6.8) If η k > 0 is a computable estimator for |u k -U k-1 | V , i.e., |u k -U k-1 | V ≤ η k , (6.9)
then, using (6.7), we get that a sufficient condition for (6.8) is

η k ≤ δ |u k -u k-1 | V .
(6.10) Algorithm 6.3 (multilevel exact Uzawa). Let p 0 ∈ P be any approximation for p.

For k = 1, 2, . . . , construct (u k , p k ), with u k ∈ V k , by a(u k , v) = f, v -b(v, p k-1 ), v ∈ V k , p k = p k-1 + α(Bu k -g).
As a consequence of Theorem 6.2 we have the following. Corollary 6.4. Let 0 < α < 2/M 2 , γ = I -αBA -1 B * , and assume that (6.8) or (6.9)-(6.10) are satisfied with δ < 1-γ 1-γ+2αM 2 . Then, the multilevel exact Uzawa algorithm converges and the convergence rate is given by (4.4).

Multilevel inexact Uzawa.

In this subsection, we assume that the problem (6.4) can be solved on each V k with an absolute error k ∈ [0, δ) , i.e.,

|d k -dk | V ≤ k |d k | V . (6.11) If η k > 0 is a computable estimator for |d k -D k-1 | V , i.e., |d k -D k-1 | V ≤ η k , (6.12)
then a computable sufficient condition for (6.6) is

η k ≤ δ -k 1 + k | dk | V . (6.13)
Indeed, from (6.7) and (6.11)-(6.13) and the triangle inequality we have

| dk -D k-1 | V ≤ |d k -D k-1 | V + |d k -dk | V ≤ η k + k |d k | V ≤ δ k -k 1 + k | dk | V + k |d k | V ≤ δ k -k 1 + k (1 + k )|d k | V + k |d k | V = δ k |d k | V ≤ δ k |D k-1 | V .
We conclude this subsection with a corollary and some remarks. Corollary 6.5. Let 0 < α < 2/M 2 , γ = I -αBA -1 B * , and let d k , dk satisfy (6.11)-(6.13) with δ < 1-γ 1-γ+2αM 2 . Then, the MIUM converges and the convergence rate is given by (4.4).

6.3.

Multilevel and adaptive interpretation of the inexact Uzawa algorithm. We note that the modified inexact Uzawa algorithm can be interpreted as a multilevel algorithm. We consider that a sequence {M k } of approximating subspaces of V is constructed such that M k is strictly larger than M k-1 and that M k is built from M k-1 by a uniform refinement strategy (see, e.g., [START_REF] Bacuta | New interpolation results and applications to finite element methods for elliptic boundary value problems[END_REF][START_REF] Bramble | Computational scales of Sobolev norms with application to preconditioning[END_REF][START_REF] Bramble | The analysis of multigrid methods[END_REF][START_REF] Xu | Multilevel Finite Element Theory[END_REF]). Based on this existing sequence of nested subspaces of V, we can now build a sequence {V k } so that (6.6) holds as follows.

Take V 0 = M 0 , and for any positive integer k, assuming that V k-1 = M j is known, define V k+i := M j for i = 0, 1, . . . as long as (6.13) is satisfied for k replaced by k + i. In other words, we update u k-1 without enlarging the space V k-1 as long as (6.13) is satisfied. When (6.13) fails to hold, we solve for the u k on the next discrete level space.

The modified inexact Uzawa algorithm can be also interpreted as an adaptive method. We construct the sequence {V k } (so that (6.6) holds) by starting with a subspace V 0 of V with good approximation properties and by building the sequence {V k } k≥1 in a similar manner. If (6.13) fails to hold for V k = V k-1 , then the new discrete space V k is constructed by using an adaptive strategy which assures that (6.13) and consequently (6.6) hold.

Applications to the Stokes system. We consider the stationary Stokes equations

-Δu -∇p = f in Ω, div u = g in Ω, (7.1)
with vanishing Dirichlet boundary condition u = 0 on ∂Ω and g satisfying the constraint Ω g dx = 0.

In this section we apply the abstract Uzawa results presented in the previous sections to solve (7.1).

Let V := (H 1 0 (Ω)) d , d = 2 or = 3, and

P = L 2 0 (Ω) := h ∈ L 2 (Ω)| Ω h dx = 0 .
We assume that f ∈ (L 2 (Ω)) d and g ∈ L 2 (Ω). The variational formulation of (7.1) becomes Find u ∈ V, p ∈ P such that (∇u, ∇v) We note that, for Ω smooth enough, we have a(u, v) := (∇u, ∇v) = (curl u, curl v) + (div u, div v), u, v ∈ V. (7.3) We denote the norm induced by a with | • | V or | • |. The norm on P is the L 2standard norm and is simply denoted by • . With the above notation, the variational formulation of (7.1) becomes (1.3).

+(divv, p) = (f , v), v ∈ V. (div u, q) = ( g, q), q ∈ P, ( 7
It is known that for Ω smooth enough, the following LBB condition holds. More precisely, we have

: V → V * consists of d copies of -Δ : H 1 0 (Ω) → H -1 (Ω), Bv = div v, B * p = -∇p, and for ρ > 0, a ρ (u, v) := a(u, v) + ρ(div u, div v), u, v ∈ V.
The next two theorems are direct consequences of Theorems 2.5 and 3.3, respectively. Theorem 7.1. Let (u, p) be the solution of (7.2) and let (u k , p k ) be the sequence of approximations built by the UM (2.13). Then the statements (i)-(iv) of Theorem 2.5 hold with m = c 0 and M = 1.

Theorem 7.2. Let (u, p) be the solution of (7.2) and let (u k , p k ) be the sequence of approximations built by the ALUM (3.1). Then the statements (i)-(iv) of Theorem 3.3 hold with m = c 0 and M = 1.

According to section 5, both the UM and ALUM can be applied to any discretization of (7.2), provided that (V h , P h ), with V h ⊂ V and P h ⊂ P , is a stable pair. Let us assume that a fixed pair (V h , P h ) satisfies the discrete inf-sup and sup-sup conditions with constants m(h) = c d > 0 and M (h) = 1. If Q h : P → P h is the L 2 -orthogonal projection, then, with the new spaces, the operators associated with the forms a and b are A h and B h , respectively, where A h : V h → V * h consists of d copies of the discrete Laplacian and B h v = Q h div v. Thus, the update for the pressure becomes

p k = p k-1 + α Q h (div u k -g).
The analysis of the discrete versions of the UM and the ALUM can be carried on similarly. The only difference in describing the convergence of the two algorithms for the discrete case is that c 0 in Theorems 7.1 and 7.2 is replaced by c d .

The inexact Uzawa algorithm can be also applied for the discretization of (7.2) on (V h , P h ) (see, e.g., [START_REF] Bramble | Analysis of the inexact Uzawa algorithm for saddle point problems[END_REF]). Taking for example C h : V * h → V h to be a preconditioner for A h such that (4.1) is satisfied with δ < , we have that the IUM converges for any α ∈ (0, 2). We can also apply the inexact Uzawa algorithm for the augmented Lagrangian Uzawa formulation on (V h , P h ) (see Remark 4.4).

According to Corollary 6.5, the MIUM for solving (7.2) can be also applied for any δ < c 2 0 /(2 + c 2 0 ). The main difficulty in doing so is to find the sequence of spaces {V k } such that (6.6) or (6.13) is satisfied. Residual-type a posteriori estimators η k (see, e.g., [START_REF] Babuška | A feedback finite element method with a posteriori error estimations: Part I. The finite element method and some basic properties of the a posteriori error estimator[END_REF], [START_REF] Verfurth | A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]) could be involved in finding the right sequence {V k }. Constructing and testing multilevel or adaptive algorithms for solving the Stokes system based on the MIUM remains a challenging new problem and is a subject for future work.

Conclusion.

The paper gives a unified analysis approach of various Uzawalike algorithms for solving continuous or discrete saddle point problems. The convergence condition and the convergence factors depend upon the extreme spectral bounds of the Schur complement BA -1 B * only. To the best of our knowledge, the result concerning the optimal convergence factor of the ALUM for the infinite-dimensional case is new. The analysis of the modified inexact Uzawa algorithm at the continuous level, which was introduced in section 6, gives a general strategy for solving saddle point systems. Our inexact Uzawa algorithm is similar to the algorithm for solving the Stokes system presented in [START_REF] Bänsch | An adaptive Uzawa FEM for the Stokes problem: Convergence without the inf-sup condition[END_REF]. The differences are in the way the error bounds are imposed (see (6.6), (6.13)) and the way the pressure is updated. Our analysis, combined with standard techniques of a posteriori error estimates, could lead to new and efficient adaptive algorithms for solving saddle point systems. The main difficulty in implementing concrete algorithms based on the MIUM is finding error estimators η k such that the conditions (6.6) or (6.13) are satisfied. Finding spaces {V k } such that conditions similar to (6.6) are satisfied will be the focus of our future work.

  •) and a(•, •) give rise to equivalent norms on V, we have that (3.2) (or (3.3)) has a unique solution. Consequently, problems (1.3) and (3.2) are equivalent.

. 2 )

 2 where (•, •) represents the standard L 2 -inner product. We will denote by a(•, •) and b(•, •) the bilinear forms a(u, v) := (∇u, ∇v) = d i=1 (∇u i , ∇v i ) and b(v, p) := (div v, p), v ∈ V, p ∈ P.

5 )

 5 We notice that for the Stokes problem the operator A

  [•, •] w with • w . Theorem 4.2. Let 0 < α < 2/M 2 and assume that C satisfies (4.1)

Acknowledgment. The author would like to thank the two reviewers for the valuable suggestions towards improving this paper.

This work was partially supported by the University of Delaware Research Foundation.