---
abstract: 'We discuss efficient solutions to systems of shifted linear systems arising in computations for oscillatory hydraulic tomography (OHT). The reconstruction of hydrogeological parameters such as hydraulic conductivity and specific storage using limited discrete measurements of pressure (head) obtained from sequential oscillatory pumping tests, leads to a nonlinear inverse problem. We tackle this using the quasi-linear geostatistical approach [@kitanidis1995quasi]. This method requires repeated solution of the forward (and adjoint) problem for multiple frequencies, for which we use flexible preconditioned Krylov subspace solvers specifically designed for shifted systems based on ideas in [@gu2007flexible]. The solvers allow the preconditioner to change at each iteration. We analyze the convergence of the solver and perform an error analysis when an iterative solver is used for inverting the preconditioner matrices. Finally, we apply our algorithm to a challenging application taken from oscillatory hydraulic tomography to demonstrate the computational gains by using the resulting method.'
author:
- 'Arvind K. Saibaba'
- Tania Bakhos
- 'Peter K. Kitanidis'
title: A Flexible Krylov Solver for Shifted Systems with Application to Oscillatory Hydraulic Tomography
---

Introduction
============

Hydraulic tomography (HT) is a method for characterizing the subsurface that consists of applying pumping in wells while aquifer pressure (head) responses are measured. Using the data collected at various locations, important aquifer parameters (e.g., hydraulic conductivity and specific storage) are estimated. An example of such a technique is transient hydraulic tomography (reviewed in [@cardiff20113D]). Oscillatory hydraulic tomography (OHT) is an emerging technology for aquifer characterization that involves a tomographic analysis of oscillatory signals. Here we consider that a sinusoidal signal of known frequency is imposed at an injection point and the resulting change in pressure is measured at receiver wells. Consequently, these measurements are processed using a nonlinear inversion algorithm to recover estimates for the desired aquifer parameters. Oscillatory hydraulic tomography has notable advantages over transient hydraulic tomography; namely, a weak signal can be distinguished from the ambient noise and by using signals of different frequencies, we are able to extract additional information without having to drill additional wells.

Using multiple frequencies for OHT has the potential to improve the quality of the image. However, it involves considerable computational burden. Solving the inverse problem, i.e. reconstructing the hydraulic conductivity field from pressure measurements, requires several application of the forward (and adjoint) problem for multiple frequencies. As we shall show in section \[sec:application\], solving the forward (and adjoint) problem involves the solution of shifted systems for multiple frequencies. For finely discretized grids, the cost of solving the system of equations corresponding to each frequency can be high to the extent that it might prove to be computationally prohibitive when many frequencies are used, for example, on the order of $200$. The objective is to develop an approach in which the cost of solving the forward (and adjoint) problem for multiple frequencies is not significantly higher than the cost of solving the system of equations for a single frequency - in other words, the cost should depend only weakly on the number of frequencies.

Direct methods, such as sparse LU, Cholesky or LDL$^T$ factorization, are suited to linear systems in which the matrix bandwidth is small, so that the fill-in is somewhat limited. An additional difficulty that direct methods pose is that for solving a sequence of shifted systems, the matrix has to be re-factorized for each frequency, resulting in a considerable computational cost. By contrast, Krylov subspace methods for shifted systems are particularly appealing since they exploit the shift-invariant property of Krylov subspaces [@simoncini2007recent] to obtain approximate solutions for all frequencies by generating a single approximation space that is shift independent. Several algorithms have been developed for dealing with shifted systems. Some are based on Lanczos recurrences for symmetric systems [@meerbergen2003solution; @meerbergen2010lanczos]; others use the unsymmetric Lanczos [@freund1993solution], and some others use Arnoldi iteration [@datta1991arnoldi; @frommer1998restarted; @simoncini2003restarted; @darnell2008deflated; @gu2007flexible]. Shifted systems also occur in several other applications such as control theory, time dependent partial differential equations, structural dynamics, and quantum chromodynamics (see [@simoncini2003restarted] and references therein). Hence, several other communities can benefit from advances in efficient solvers for shifted systems. The Krylov subspace method that we propose is closest in spirit to [@gu2007flexible]. However, as we shall demonstrate, we have extended their solver significantly.

**Contributions**: Our major contributions can be summarized as follows:

-   We have extended the flexible Arnoldi algorithm discussed in [@gu2007flexible] for shifted systems of the form $(A + \sigma_jI)x_j =b$ to systems of the form $(K + \sigma_j M) x_j = b $ for $j = 1,\dots,n_f$ that employs multiple preconditioners of the form $(K+\tau M)$. In addition, we provide some analysis for the convergence of the solver.

-   When an iterative solver is used to apply the preconditioner, we derive an error analysis that gives us stopping tolerances for monitoring convergence without constructing the full residual.

-   Our motivation for the need for fast solvers for shifted systems comes from oscillatory hydraulic tomography. We describe the key steps involved in inversion for oscillatory hydraulic tomography, and discuss how the computation of the Jacobian can be accelerated by the use of the aforementioned fast solvers.

**Limitations**: The focus of this work has been on the computational aspects of oscillatory hydraulic tomography. Although the initial results are promising, several issues remain to be resolved for application to realistic problems of oscillatory hydraulic tomography. For example, we are inverting for the hydraulic conductivity assuming that the storage field is known. In practice, the storage is also unknown and needs to be estimated from the data as well. Moreover, simulating realistic conditions (higher variance in the log conductivity field, and adding measurement noise in a realistic manner) may significantly improve the performance with the addition of information from different frequencies. We will deal with these issues in another paper.

The paper is organized as follows. In section \[sec:krylov\], we discuss the Krylov subspace methods for solving shifted linear systems of equations based on the Arnoldi iteration using preconditioners that are also shifted systems. In section \[sec:geneigen\], we discuss the convergence of the iterative solver and its connection to the convergence of the eigenvalues of the generalized eigenvalue problem $Kx=\lambda Mx$. In section \[sec:inexact\], we discuss an error analysis when an iterative method is used to invert the preconditioner matrices. In section \[sec:application\], we discuss the basic constitutive equations in OHT, which can be expressed as shifted linear system of equations and discuss the geostatistical method for solving inverse problems. Finally, in section \[sec:numerical\] we present some numerical results on systems of shifted systems and then discuss numerical results involving the inverse problem arising from OHT. We observe significant speed-ups using our Krylov subspace solver.

Krylov subspace methods for shifted systems {#sec:krylov}
===========================================

The goal is to solve systems of equations of the form $$\label{eqn:multipleshifted}
  \left( K + \sigma_j M \right) x_j  = b \qquad j=1,\dots,n_f$$

Note that $\sigma_j$, for $j=1,\dots,n_f$ are (in general) complex shifts. We assume that none of these systems are singular. In particular, for our application, both $K$ and $M$ are stiffness and mass matrices respectively and are positive definite, but our algorithm only requires that they are invertible. By using a finite volume or lumped mass approach [@hughes2012finite], the mass matrices become diagonal but this assumption is not necessary. Later, in sections \[sec:forward\] and \[sec:sensitivity\], we will show how such equations arise in our applications.

(0,0) rectangle ( 10, 10); (A) at (5,5) [$A$]{}; (10.5, 0) rectangle (13.5, 10); (Vm) at (12, 5) [$V_m$]{}; (eq) at (14.5, 5) [$=$]{}; (15.5,0) rectangle (18.5, 10); (Vm1) at (17,5) [$V_m$]{}; (18.7,0) rectangle (19,10); (vm1) at (20,-0.5) [$v_{m+1}$]{};

(20,10) – (23, 10) – (23,7) – (22.6, 7) – (20,9.6) – (20, 10) ; (Hm) at (22,9) [$H_m$]{};

(22.8,6.7) circle \[radius = 0.2\]; (hm1m) at (24.3,5.6) [$h_{m+1,m}$]{};

As a brief introduction, we review the Krylov based iterative solvers for the system of equations $Ax = b$. In particular, we describe the variants generated by Arnoldi iteration, such as Full Orthogonalization Method (FOM) and Generalized Minimum RESidual method (GMRES). Krylov solvers typically generate a sequence of orthonormal vectors $v_1,\dots,v_m$ that are orthonormal. These vectors form a basis for the Krylov subspace: $${{\cal{K}}_{m}(A,b)} {\stackrel{\text{def}}{=}}{\text{span}\left\{b,Ab,\dots,A^{m-1}b\right\}}$$ At the end of the $m$th iteration, a typical relation is obtained of the form (see figure \[fig:arnoldi\]), $$AV_m = V_{m+1}\bar{H}_m$$ where, $V_m = [v_1,\dots,v_m]$ and $v_1 = b/\beta$, $\beta = {\lVert b \rVert_2}$ with $x_0= 0$ and $\bar{H}_m$ is an upper Hessenberg matrix. Then, an approximate solution to the system $Ax = b$ by searching for a solution of the form $x_m =  V_my_m$, where $y_m$ is chosen such that it minimizes the residual $r_m {\stackrel{\text{def}}{=}}b - Ax_m$ which, leads to GMRES subproblem, $$\label{eqn:gmressubproblem} 
\min_{y_m \in \mathbb{C}^m}{\lVert r_m \rVert_2}   \Rightarrow \min_{y_m \in \mathbb{C}^m} {\lVert \beta e_1 - \bar{H}_my_m \rVert_2}$$ or an oblique projection $r_m \perp {\text{span}\left\{V_m\right\}}$ which leads to the FOM subproblem $$\label{eqn:fomsubproblem} 
r_m \perp {\text{span}\left\{V_m\right\}}  \qquad \Rightarrow \qquad H_my_m = \beta e_1$$ As $m$ increases, the cost per iteration increases at least as ${\cal{O}}(m^2n)$ and the memory costs increase as ${\cal{O}}(mn)$[@saad2003iterative]. The standard remedies to reducing the number of iterations are 1) using an appropriate preconditioner, 2) truncating the orthogonalization in the Arnoldi algorithm and 3) restarting the Arnoldi algorithm periodically.

An interesting property of the Krylov subspaces is that they are shift-invariant. In other words, ${{\cal{K}}_{m}(A,b)} ={{\cal{K}}_{m}(A + \sigma I,b)}$. Therefore, the same Krylov basis generated for the system $Ax=b$ can be effectively used to solve shifted systems of the form $(A+\sigma I)x = b$. The strategy for solving the shifted systems is therefore, to first generate a basis that is applicable to all systems, and then use the shift-invariant property (for a detailed review, see [@simoncini2007recent section 14.1] and references therein) to solve a smaller subproblem of the form  or . The same idea can be extended to systems of the form  using a preconditioner of the form $(K+\tau M)$ [@meerbergen2003solution; @gu2007flexible], which solves for multiple shifted systems roughly at the cost of a single system. However, in practice, the number of iterations taken can often be large, especially for large matrices arising from realistic applications.

In order to minimize the number of iterations, Meerbergen [@meerbergen2003solution] proposes a left preconditioner of the form $K_\tau {\stackrel{\text{def}}{=}}K + \tau M $ that is factorized and inverted using a direct solver. The application of $K_\tau^{-1}$ to a vector is, in general, not cheap but the spectrum of $(K+\tau M)^{-1}(K + \sigma M)$ is often more favorable, which results in fast convergence of the Krylov methods in just a few iterations [@meerbergen2003solution]. This form of preconditioning has its roots in solving large-scale generalized eigenvalue problems and is known as [*Cayley transformation*]{} [@golub1996matrix]. In [@popolizio2008acceleration], the authors provide some analysis for choosing the best value of $\tau$ that optimally preconditions all the systems. However, we observed that (also, see [@gu2007flexible]) using a single preconditioner for all the systems may not yield optimal convergence for all systems. In [@gu2007flexible], the authors propose a flexible Arnoldi method for shifted systems that uses different values of $\tau$ resulting in different preconditioners at each iteration. This can potentially reduce the number of iterations for all the shifted systems. Before we describe our flexible algorithm in section \[sec:flexibleprecond\], we will derive the right preconditioned version of Krylov subspace method for shifted systems. This serves two purposes - it motivates our algorithm, while clarifying some of the notation.

Right preconditioning for shifted systems {#sec:rightprecond}
-----------------------------------------

As mentioned earlier, we will review the right preconditioned version of the Krylov subspace algorithm for shifted systems. Following the approach in [@meerbergen2003solution; @simoncini2007recent], we solve the system of equations  using a shifted right preconditioner of the form $K_\tau {\stackrel{\text{def}}{=}}K+\tau M $ $$\label{eqn:rightpreconditioned}
  (K + \sigma_j M) K_\tau^{-1} \bar{x}(\sigma_j) = b \qquad x(\sigma_j) =  K_\tau^{-1}\bar{x}(\sigma_j)$$ for $j=1,\dots,n_f$. We have the following identity that $$\label{eqn:shiftinvariant}
 (K+\sigma M) (K + \tau M)^{-1} = I + (\sigma - \tau) M (K+\tau M)^{-1}$$

Using the identity in equation , we have the following shift-invariance property ${{\cal{K}}_{m}(MK_\tau^{-1},b)} = {{\cal{K}}_{m}((K+\sigma M)(K+\tau M)^{-1},b)}$. Note that ${{\cal{K}}_{m}(MK_\tau^{-1},b)} $ is independent of $\sigma$. This shift-invariance property suggests an efficient algorithm for solving the system of equations . There is a distinct advantage in using iterative solvers for shifted systems; the expensive step of constructing the basis for the Krylov subspace is performed only once and using the shift-invariance property of the Krylov subspace, the sub-problem for each shift in algorithm \[alg:shiftedsolve\] can be computed at a relatively low cost.

$M$ and a right hand side $b$. Compute $v_1 = b /\beta$ and $\beta {\stackrel{\text{def}}{=}}{\lVert b \rVert_2}$ Choose $\tau$ and factorize $K_\tau {\stackrel{\text{def}}{=}}K + \tau M$ Define the $(m+1)\times m$ matrix $\bar{H}_m = \{ h_{i,k}\}_{1\leq i \leq m+1, 1 \leq k \leq  m}$. Set $\bar{H}_m = 0$. Compute ${z}_k = K_\tau^{-1} v_k$ $w_k := Mz_k$ $h_{ik} := w_k^*v_i$ Compute $w_k := w_k - h_{i,k}v_k$ $h_{k+1,k}:= \lVert w_k \rVert_2 $. If $h_{k+1,k} = 0$ stop $v_{k+1} = w_k/ h_{k+1,k}$

The algorithm proceeds as follows: first, we run $m$ steps of the Arnoldi algorithm on the matrix $M K_\tau^{-1}$ with the starting vector $b$ to get a basis for the Krylov subspace ${{\cal{K}}_{m}(MK_\tau^{-1},b)}$. This is summarized in algorithm \[alg:arnoldi\]. At the end of $m$ steps of the Arnoldi process, we construct two sets of vectors $V_{m+1} = [v_1,\dots,v_{m+1}]$ and $Z_m = [z_1,\dots,z_m ]$, and an upper Hessenberg matrix $\bar{H}_m$ that satisfy the following relations,

$$\begin{aligned}
\label{eqn:arnoldi}
 M  Z_m  = & \quad V_{m+1}\bar{H}_m \\
   (K + \tau M) Z_m =  & \quad V_m  \label{eqn:zm}\end{aligned}$$

where, $ V_m^* V_m = I $. Multiplying the first equation by $(\sigma_j-\tau)$ and adding it to the second equation gives us

$$\label{eqn:shiftedarnoldi}
 (K + \sigma_j M)Z_m = V_{m+1} \underbrace{\left(\begin{bmatrix} I \\ {0}\end{bmatrix} + (\sigma_j-\tau) \bar{H}_m  \right)}_{{\stackrel{\text{def}}{=}}\bar{H}_m(\sigma_j; \tau)} = V_{m+1}\bar{H}_m(\sigma_j;\tau)$$

In algorithm \[alg:arnoldi\], $V_m$ forms a basis for the Krylov subspace ${\cal{K}}_m(MK_\tau^{-1},b)$. However, we seek solutions of the form $x_m =  Z_my_m$ (with zero as the initial guess). Now $x_m \in {\text{span}\left\{Z_m\right\}}$ where $Z_m$ is the space spanned by the vectors $z_k = K_\tau^{-1}v_k$ for $k=1,\dots,m$. By minimizing the residual norm over all possible vectors in ${\text{span}\left\{Z_m\right\}}$, we obtain the generalized minimum residual (GMRES) method for shifted systems, whereas by imposing the Petrov-Galerkin condition $r_m \perp{\text{span}\left\{V_m\right\}}$, we obtain the full orthogonalized method (FOM) for shifted systems. This is summarized in algorithm \[alg:shiftedsolve\]. It should be noted that the way we have described this algorithm, we need to store the vectors $Z_m$. In practice, this is not necessary. We chose to present it this way in order to have consistent notation with the flexible algorithm we will describe in subsection \[sec:flexibleprecond\].

matrices $K$ and $M$, a right hand side $b$, $\sigma \in \{ \sigma_1,\dots,\sigma_{n_f}\}$ Choose $\tau$, build $K_\tau {\stackrel{\text{def}}{=}}K + \tau M$ and construct preconditioner. Set $m = 1$. Generate $V_{m+1},\bar{H}_m$ and $Z_m$ using algorithm \[alg:arnoldi\]. Construct $\bar{H}_m(\sigma_j;\tau) {\stackrel{\text{def}}{=}}I + (\sigma_j -\tau)\bar{H}_m$ (see equation ). FOM: $$H_m(\sigma_j; \tau)y_m^{fom}(\sigma_j) =  \beta e_1$$ GMRES: $$y_m^{gmres}(\sigma_j) {\stackrel{\text{def}}{=}}\min_{ y_m \in \mathbb{C}^m} \lVert \beta e_1 - \bar{H}_m(\sigma_j;\tau) y_m \rVert_2$$ Construct the approximate solution $x_m(\sigma_j) = Z_my_m(\sigma_j)$ $m\leftarrow m + 1$.

In [@meerbergen2003solution], the spectrum of $(K+\sigma M)(K+\tau M)^{-1}$ was analyzed and it was shown that the preconditioner $ K_\tau$ is well suited only for values of frequencies $\sigma$ near $\tau$. However, the values of $\sigma$ can be widely spread and a single preconditioner $K_\tau$ might not be a good choice for preconditioning all the systems. In section \[sec:solverexperiments\], we demonstrate an example in which a single preconditioner does not satisfactorily precondition all the systems. In [@gu2007flexible], the authors propose a flexible approach using a (possibly) different preconditioner at each iteration. We shall adopt this approach.

Flexible preconditioning {#sec:flexibleprecond}
------------------------

We now describe our flexible Krylov approach for solving shifted systems based on [@gu2007flexible] which we have extended to the case $M\neq I$. Following [@saad1993flexible] and [@gu2007flexible], we use a variant of GMRES which allows a change in the preconditioner at each iteration. In algorithm \[alg:arnoldi\], we considered a fixed preconditioner of the form $K_\tau {\stackrel{\text{def}}{=}}K + \tau M$ for a fixed $\tau$. Suppose we used a different preconditioner at each iteration of the form $K+\tau_k M$ for $k=1,\dots,m$, then instead of  we have, $$(K+\tau_k M) z_k =  v_k  \qquad k=1,\dots,m \label{eqn:zmmod}$$ The algorithm is summarized in algorithm \[alg:shiftedsolvemod\]. In this algorithm, in addition to saving $V_m$, we also save the matrix $Z_m$. If at every step in the flexible Arnoldi algorithm we use the same value of $\tau$, we are in the same position as in algorithm \[alg:arnoldi\]. We have $Z_m = [z_1,\dots,z_m]$, $\bar{H}_m = \{h_{ik}\}_{1\leq i \leq m+1,1\leq k \leq m}$ and $V_{m} = [v_1,\dots,v_{m}]$ which satisfies $V_m^*V_m = I_m$. In addition, we also have the following relations

$$\begin{aligned}
 MZ_m = & \quad V_{m+1}\bar{H}_m \label{eqn:arnoldimod} \\
 KZ_m + MZ_m T_m = & \quad V_m   \label{eqn:zmvmrelation}\end{aligned}$$

where, $T_m = \text{diag}\{ \tau_1,\dots,\tau_m\}$. Multiplying  by $\sigma_jI_m - T_m$ and adding , we obtain for $j = 1,\dots,n_f$ $$(K + \sigma_jM )Z_m = V_{m+1} \underbrace{\left( \begin{bmatrix} I \\ 0 \end{bmatrix}  + \bar{H}_m (\sigma_jI_m - T_m) \right)}_{{\stackrel{\text{def}}{=}}\bar{H}(\sigma_j; T_m)} = V_{m+1}\bar{H}_m(\sigma_j;T_m) \label{eqn:shiftedarnoldimod}$$

$M$ and $b$ the right hand side, $\tau_k, k=1,\dots,m$, $v_1 = b /\beta$ and $\beta {\stackrel{\text{def}}{=}}{\lVert b \rVert_2}$ Define the $(m+1)\times m$ matrix $\bar{H}_m = \{h_{i,k} \}_{1\leq i \leq m+1, 1 \leq k \leq  m}$. Solve $ (K+\tau_kM) {z}_k =  v_k$ $w_k := Mz_k$ $h_{i,k} := w_k^*v_i$ Compute $w_k := w_k - h_{i,k}v_k$ $h_{k+1,k}:= \lVert w_k \rVert_2 $. If $h_{k+1,k} = 0$ stop $v_{k+1} = w_k/ h_{k+1,k}$

We are now in a position to derive a FOM/GMRES algorithm for shifted systems with flexible preconditioning. We search for solutions which are approximations of the form $x_m(\sigma_j) = Z_m y_m(\sigma_j)$, which spans the columns of $Z_m$. Strictly speaking, $\text{span}\{Z_m\}$ is no longer a Krylov subspace. By minimizing the residual norm over all possible vectors in ${\text{span}\left\{Z_m\right\}}$, we obtain the flexible generalized minimum residual (FGMRES) method for shifted systems, whereas by imposing the Petrov-Galerkin condition $r_m \perp{\text{span}\left\{V_m\right\}}$, we obtain the flexible full orthogonalized method (FFOM) for shifted systems. This is summarized in algorithm \[alg:shiftedsolvemod\]. The residuals can be computed as $$\begin{aligned}
\label{eqn:fomgmresresidual}
r_m(\sigma_j) = & \quad b - (K+\sigma_jM)x_m(\sigma_j) \\ \nonumber
              = & \quad V_{m+1}\left(\beta e_1 - \bar{H}_m(\sigma_j;T_m)y_m(\sigma_j)\right) \\ \nonumber\end{aligned}$$

matrices $K$ and $M$, vector $b$, $\sigma \in \{ \sigma_1,\dots,\sigma_{n_f}\}$, set $m=1$.

Choose $T_m = \text{diag}\{\tau_1,\dots,\tau_m\}$. Generate $V_{m+1},\bar{H}_m$ and $Z_m$ using algorithm \[alg:arnoldimod\]. Construct $\bar{H}_m(\sigma_j;T_m){\stackrel{\text{def}}{=}}I + \bar{H}_m(\sigma_jI-T_m)$ (see equation ). FOM: $$H_m(\sigma_j; T_m)y_m^{fom}(\sigma_j) =  \beta e_1$$ GMRES: $$y_m^{gmres}(\sigma_j) {\stackrel{\text{def}}{=}}\min_{ y_m \in \mathbb{C}^m} \lVert \beta e_1 - \bar{H}_m(\sigma_j;T_m) y_m \rVert_2$$ Construct the approximate solution as $x_m(\sigma_j) = Z_m y_m(\sigma_j)$ $m\leftarrow m+1$

### Selecting values of $\tau_k$

In algorithm \[alg:arnoldimod\], at each iteration we solve a system of the form $(K+\tau_k M)z_k = v_k$ for $k=1,\dots,m$. This cost can be high if the dimension of the Arnoldi subspace $m$ is large and a different preconditioner $K+\tau_kM$ is used at every iteration. In practice, it is not necessary to form and factorize $m$ systems corresponding to different $\tau_k$. In applications, we only need choose a few different $\tau_k$ that cover the entire range of the parameters $\sigma_j$. This was also described in [@gu2007flexible]. The system of equations  is solved using a direct solver and since only a few values of $\tau_k$ are chosen, the systems can be formed and factorized. Thus, the computational cost will not be affected greatly even if the number of frequencies $n_f$ is large.

Let $\bar{\tau} = \{\bar{\tau}_1,\dots,\bar{\tau}_{n_p} \}$ be the set of values that $\tau_k$ can take. In other words, we take $n_p$ distinct preconditioners. Then, the first $m_1$ values of $\tau_k$ are assigned $\bar{\tau}_1$, the next $m_2$ values of $\tau_k$ are assigned $\bar{\tau}_2$ and so on. We also have $m = \sum_{k=1}^{n_p}m_k$.

### Restarting {#sec:restarting}

As the dimension of the subspace $m$ increases, the computational and memory costs increase significantly. A well known solution to this problem is restarting. The old basis is discarded and the Arnoldi algorithm is restarted on a new residual. However, for shifted systems, in order to preserve the shift-invariant property, one needs to ensure collinearity of the residuals of the shifted systems. For FOM the residuals are naturally collinear and the Arnoldi algorithm can be restarted by scaling each residual by some scalar that depends on the shift [@simoncini2003restarted; @gu2007flexible]. For GMRES, the approach used by [@frommer1998restarted] was extended to shifted systems with multiple preconditioners by [@gu2007flexible]. We did not explore this issue further, and the reader is referred to [@gu2007flexible] for further details.

Generalized eigenvalue problem and error estimates {#sec:geneigen}
==================================================

We start by computing the approximate eigenvalues and eigenvectors for the matrix $(K+\sigma M)M^{-1}$. Using estimates for approximate eigenvalues and eigenvectors we derive expressions for the convergence of the flexible algorithms. The approximate eigenvalues are called Ritz values. For convenience, we drop the subscript on the shifted frequency, i.e., use $\sigma$ instead of $\sigma_j$ where, $j=1,\dots,n_f$.

\[prop:ritz\] Let $Z_m$, $\bar{H}_m$ and $V_{m+1}$ be computed according to algorithm \[alg:arnoldimod\]. Calculate the eigenpairs of the generalized eigenvalue problem $$\label{eqn:genritz}
H_m(\sigma;T_m) f  = \theta  H_m f$$ Then, the Ritz pair $\left( \theta, u {\stackrel{\text{def}}{=}}V_{m+1}\bar{H}_mf \right)$ satisfy the Petrov-Galerkin condition [@saad1992numerical section 4.3.3] $$\label{eqn:petrovgalerkin}
(K + \sigma M)M^{-1}u - \theta u \perp {\text{span}\left\{V_m\right\}} \qquad u \in {\text{span}\left\{V_{m+1}\bar{H}_m\right\}}$$

We first begin by manipulating equations  and . Eliminating $Z_m$ from those equations and adding $\sigma V_{m+1}\bar{H}_m$ to both sides, we have $$\label{eqn:rationalarnoldi}
V_{m+1}\bar{H}_m(\sigma;T_m) =  (K+\sigma M)M^{-1}V_{m+1}\bar{H}_m$$ Now, consider the residual of the eigenvalue calculation for the $k$th eigenpair, where $k=1,\dots,m$ is

$$\begin{aligned}
\label{eqn:residualritzeigen}
 r^\text{eig}_k(\sigma) = &\quad  (K+\sigma M)M^{-1}V_{m+1}\bar{H}_m f_k - \theta_k V_{m+1}\bar{H}_mf_k  \\ \nonumber
           = & \quad  V_{m+1}\bar{H}_m(\sigma;T_m)f_k - \theta_k V_{m+1}\bar{H}_mf_k  \\ \nonumber
       = &\quad  V_m\left(H_m(\sigma;T_m)f_k - \theta_kH_mf_k \right) - h_{m+1,m}v_{m+1}(\tau_m +\theta_k - \sigma)e_m^*f_k      \\ \nonumber
       = & \quad - h_{m+1,m}v_{m+1}(\tau_m +\theta_k - \sigma)e_m^*f_k   \\ \nonumber\end{aligned}$$

From which we can claim that $u \in {\text{span}\left\{V_{m+1}\bar{H}_m\right\}}$ and $(K+\sigma M)M^{-1}u - \theta u \perp {\text{span}\left\{V_m\right\}}$. In other words, they satisfy the Petrov-Galerkin  and are an approximate eigenpair of $(K+\sigma M)M^{-1}$.

Furthermore, we define $\rho_k {\stackrel{\text{def}}{=}}{\lVert (K+\sigma M)M^{-1}u_k-\theta_k u_k \rVert_2}$, which is the residual norm of the $k$th eigenvalue calculations. When the residual of the eigenvalue calculations $\rho_k$ is small, say machine precision, the Ritz values are a good approximation to the eigenvalues. It is readily verified that the eigenvalues $\lambda$ of $KM^{-1}$ (and the generalized eigenvalue problem $Kx=\lambda Mx$) are related to the eigenvalues $\lambda(\sigma)$ of $(K+\sigma M)M^{-1}$ by the relation $\lambda(\sigma) = \lambda + \sigma$. The importance of the convergence of Ritz values to the convergence of the Krylov subspace solver using FOM can be established by the following result.

\[prop:fom\] Assume the requirements of proposition \[prop:ritz\]. Further, assume $F$ (matrix of generalized eigenvectors, see ) is invertible so that the generalized eigendecomposition $H_m(\sigma;T_m) = H_mF\Theta F^{-1}$ exists. The residual using FFOM satisfies the following inequality $$\label{eqn:fomresidualineq} {\lVert r_m(\sigma) \rVert_2} \leq \sum_{k=1}^m \rho_k\left|\frac{\sigma -\tau_m}{\theta_k +\tau_m-\sigma}\right| |\theta_k^{-1}| |s_k|$$ where, $s_k {\stackrel{\text{def}}{=}}e_k^*F^{-1}(H_m^{-1}\beta e_1)$ and $\rho_k$ is the residual norm of the eigenvalue calculation, defined above.

We start by writing the residual in equation  $$r_m(\sigma) =  V_m( \beta e_1 - H_m(\sigma;T_m)y_m) - v_{m+1}(\sigma-\tau_m)h_{m+1,m}e_m^*y_m$$ Since, for flexible FOM for shifted systems $y_m = H_m(\sigma;T_m)^{-1}\beta e_1$, the first term in the above expression is zero and we have, $$\label{eqn:fomresidualsimplified}
r_m(\sigma) =  \quad  -v_{m+1}(\sigma-\tau_m)h_{m+1,m}e_m^*H_m(\sigma;T_m)^{-1}\beta e_1  \\ \nonumber$$ Now, using the generalized eigendecomposition in  $H_m(\sigma;T_m) = H_m F\Theta F^{-1}$. Therefore, we have $H_m(\sigma;T_m)^{-1} = F\Theta^{-1} F^{-1}H_m^{-1}$. We can write this is as a sum of rank-$1$ vectors $$H_m(\sigma;T_m)^{-1} = \sum_{k=1}^m \theta_k^{-1}f_k e_k^*F^{-1}H_m^{-1}$$ where, $e_k$ is the $k$-th canonical basis vector and $f_k$ is the $k$-th column of $F$ for $k=1,\dots,m$. Using the residual of the eigenvalue calculation $r^\text{eig}_k(\sigma)$ in  and the expression derived above, $$\begin{aligned}
r_m(\sigma) \quad = & \quad -\sum_{k=1}^mv_{m+1}h_{m+1,m} (\sigma -\tau_m) e_m^*f_k\theta_k^{-1}s_k \\ \nonumber
        = & \quad \sum_{k=1}^m r^\text{eig}_k(\sigma) \frac{\sigma - \tau_m}{\tau_m +\theta_k - \sigma} \theta_k^{-1} s_k \\ \nonumber\end{aligned}$$ where, $s_k{\stackrel{\text{def}}{=}}e_k^*F^{-1}H_m^{-1}\beta e_1$. The proof follows from the properties of vector norms.

The inequality  provides insight into the importance of the accuracy of approximate eigenpairs for the convergence of flexible FOM for shifted systems. We follow the arguments in [@meerbergen2003solution]. In particular, the residual is very small if $\sigma \approx \tau_m$, $|\theta_k^{-1}|$, $|s_k|$ or $\rho_k$ are small. We shall ignore the case that $\sigma \approx \tau_m$ for further analysis, i.e. that the shifted system is almost exactly the preconditioned system. The eigenvalue residual norm $\rho_k$ being small implies that the Ritz values are a good approximation to the eigenvalues of $(K+\sigma M)M^{-1}$. This implies that all the eigenvalues in this interval have been computed fairly accurately. We now discuss when $|\theta_k^{-1}|$ is large. When all the values of $\tau_k$ are equal to $\tau$, the approximate eigenvalues $\theta_k$ of $KM^{-1}$ are related to approximate eigenvalues $\lambda_k$ of the preconditioned system $(K+\sigma M)(K+\tau M)^{-1}$ by the Cayley transformation $\frac{\lambda_k + \sigma}{\lambda_k + \tau}$. Therefore, $|\theta_k^{-1}|$ is large only if $|\lambda_k + \sigma| \ll |\lambda_k + \tau|$. The term $s_k$ can be rewritten as $s_k =  e_k^*F^{-1}H_m^{-1}V_m^*V_m\beta e_1 = e_k^*F^{-1}H_m^{-1}V_m^* b$. It is readily verified that $e_k^*F^{-1}H_m^{-1}V_m^*$ is orthonormal to all other approximate eigenvectors $V_{m+1}\bar{H}_mZ_m$ and thus, $s_k$ can be interpreted as the component of the right hand side $b$ in the direction of the approximate eigenvector. In other words, $s_k$ is small when the solution $x_m(\sigma)$ has a small component in the direction of $b$.

The analysis for the convergence of flexible FOM for shifted systems can be extended to flexible GMRES as well. The following result bounds the difference in the residuals obtained from $m$ steps using flexible FOM and flexible GMRES.

\[prop:gmres\] Let $Z_m$, $\bar{H}_m$ and $V_{m+1}$ be computed according to algorithm \[alg:arnoldimod\]. Further, from algorithm \[alg:shiftedsolvemod\] we define the flexible FOM quantities $y_m^\text{fom}(\sigma) = H_m(\sigma;T_m)^{-1}\beta e_1$, residual $r_m^\text{fom}(\sigma) = V_{m+1}(\beta e_1-\bar{H}_m(\sigma;T_m)y_m^\text{fom}(\sigma))$ and flexible GMRES quantities $y_m^\text{gmres}(\sigma) = \arg\min_{y\in \mathbb{C}^m}{\lVert \beta e_1-\bar{H}_m(\sigma;T_m)y \rVert_2}$, residual $r_m^\text{gmres}(\sigma) = V_{m+1}(\beta e_1-\bar{H}_m(\sigma;T_m)y_m^\text{gmres}(\sigma))$. Further, assume that $H_m(\sigma;T_m)$ is invertible. We have the following inequality $$\label{eqn:fomgmresresdiff}
{\lVert r_m^\text{fom}(\sigma) - r_m^\text{gmres}(\sigma) \rVert_2} \leq  \frac{\alpha(1+\alpha)}{1+\alpha^2}{\lVert r_m^\text{fom} \rVert_2}$$ where, $\eta {\stackrel{\text{def}}{=}}h_{m+1,m}(\sigma - \tau_m)$ and $\alpha {\stackrel{\text{def}}{=}}{\lVert \eta H_m^{-*}(\sigma;T_m)e_m \rVert_2}$.

We begin by the following observation from equation  $r_m^\text{fom}(\sigma) = -\eta e_m^*y_m^\text{fom} v_{m+1}$ and $$r_m^\text{fom}(\sigma) - r_m^\text{gmres}(\sigma) = -V_{m+1}\bar{H}_m(\sigma;T_m) \left( y_m^\text{fom}(\sigma) - y_m^\text{gmres} (\sigma)\right)$$ Next, we look at the solution to the GMRES least squares problem which can be written as the normal equations $$\bar{H}_m^*(\sigma;T_m)\bar{H}_m(\sigma;T_m)y_m^\text{gmres} = \bar{H}_m^*(\sigma;T_m)\beta e_1 = H_m^*(\sigma;T_m)\beta e_1$$ This can be rewritten as $$\begin{aligned}
 \left(H_m^*(\sigma;T_m)H_m(\sigma;T_m)+\eta^2e_me_m^*\right)y_m^\text{gmres}(\sigma) =& \quad  H_m^*(\sigma;T_m)\beta e_1 \\ \nonumber
\left(H_m(\sigma;T_m) + \eta^2H_m^{-*}e_me_m^*\right)y_m^\text{gmres}(\sigma) = & \quad \beta e_1 \end{aligned}$$ In other words, the solution to the GMRES subproblem is a rank-one perturbation of the FOM subproblem. Using the Sherman-Morrison identity $$y^\text{gmres}_m(\sigma) = \underbrace{H_m(\sigma;T_m)^{-1}\beta e_1}_{ =  y^\text{fom}_m(\sigma) } - H_m(\sigma;T_m)^{-1}\eta H_m^{-*}(\sigma;T_m)e_m \frac{(\eta e_m^*H_m(\sigma;T_m)^{-1}\beta e_1)}{ 1+  {\lVert \eta H_m^{-*}(\sigma;T_m)e_m \rVert_2}^2}$$ Then, the residual difference between FOM and GMRES can be bounded as $${\lVert r_m^\text{fom}(\sigma) - r_m^\text{gmres}(\sigma) \rVert_2}  \leq {\lVert \bar{H}_m(\sigma;T_m)H_m(\sigma;T_m)^{-1} \rVert_2} \frac{\alpha|\eta e_m^*y_m^\text{fom}(\sigma)|}{1+  \alpha^2}$$ The inequality  follows from the following observations ${\lVert \bar{H}_m(\sigma;T_m)H_m(\sigma;T_m)^{-1} \rVert_2} \leq 1 +  {\lVert \eta H_m^{-*}(\sigma;T_m)e_m \rVert_2}$ and from  we have ${\lVert r_m^\text{fom} \rVert_2}  =  |\eta e_m^*y_m^\text{fom}(\sigma)|$.

If ${\lVert \eta H_m^{-*}(\sigma;T_m)e_m \rVert_2}$ is large, then the difference between the two residuals can be large. This happens either when $\eta$ is large or $H_m(\sigma;T_m)$ is close to singular. In this case, flexible GMRES can stagnate and further progress may not occur. We now discuss situations in which breakdown occurs, i.e. $h_{m+1,m} =0$. If $h_{m+1,m} \neq 0$ and $H_m$ is full rank, then it can be shown from equation  that ${\text{span}\left\{MZ_m\right\}} \subseteq {\text{span}\left\{V_{m+1}\right\}} $ and from equation , it follows that ${\text{span}\left\{(K+\sigma_j M)Z_m\right\}} \subseteq {\text{span}\left\{V_{m+1}\right\}}$. Further, $h_{m+1,m} = 0$ if and only if $x_m(\sigma_j)$ is the exact solution and $H_m(\sigma_j;T_m)$ is non-singular. The argument closely follows [@saad1993flexible] and will not be repeated here.

Inexact preconditioning {#sec:inexact}
=======================

We observe that to compute vectors $z_k$ for $k =1,\dots,m$ in equation , we have to invert matrices of the form $K+\tau_kM$. When the problem sizes are large, iterative methods may be necessary to invert such matrices, resulting in a variable preconditioning procedure in which a different preconditioning operator is applied at each iteration. More precisely, for $k=1,\dots,m$, $$\label{eqn:zmmoditer}
\tilde{z}_k \approx (K + \tau_k M)^{-1}v_k \qquad p_k {\stackrel{\text{def}}{=}}v_k - (K +\tau_k M) \tilde{z}_k$$ where, $p_k$ is the residual that results after the iterative solver has been terminated. To simplify the discussion, we assume that the termination criteria for the iterative solver is such that ${\lVert p_k \rVert_2} \leq \varepsilon \underbrace{{\lVert v_k \rVert_2}}_{= 1} = \varepsilon$, for some $\varepsilon$. We closely follow the approach in [@simoncini2003theory]. The new flexible Arnoldi relationship is now,

$$\label{eqn:flexshiftedarnoldimod}
(K + \sigma_j M)\tilde{Z}_m  + P_m = V_{m+1}\bar{H}_m(\sigma_j;T_m)  \qquad j=1,\dots,n_f$$

where, $\tilde{Z}_m = [\tilde{z}_1,\dots,\tilde{z}_m]$ and $P_m = [p_1,\dots,p_m]$ and $\bar{H}_m(\sigma_j;T_m)$ is defined in equation . By using inexact applications of the preconditioner, the vectors $v_k$ for $k=1,\dots,m+1$ are no longer the same vectors generated from algorithm \[alg:arnoldimod\]. In particular, ${\text{span}\left\{V_m\right\}}$ is no longer a Krylov subspace generated by $A$. However, by construction, $V_m$ is still an orthogonal matrix.

Having constructed the matrix $\tilde{Z}_m$, we seek approximate solutions spanned by the columns of $\tilde{Z}_m$, i.e., solutions of the form $x_m(\sigma_j) = \tilde{Z}_m y_m(\sigma_j)$. The true residual corresponding to the approximation solution $x_m(\sigma_j) = \tilde{Z}_m y_m(\sigma_j)$ can be computed as follows, $$\begin{aligned}
r_m(\sigma_j)\quad   = &  \quad b - (K+\sigma_jM)\tilde{Z}_my_m(\sigma_j) \\ 
               = &  \quad b - V_{m+1}\bar{H}_m(\sigma_j;T_m)y_m(\sigma_j) + P_m y_m(\sigma_j) \\ 
           = &  \quad V_{m+1}\left(\beta e_1 - \bar{H}_m(\sigma_j;T_m)y_m(\sigma_j)\right) + P_my_m(\sigma_j) \end{aligned}$$

The columns of the matrix $P_m$ are not computed in practice because they require an additional matrix-vector product with $K+\tau_k M$. As a result, computing the true residual is expensive. However, in order to monitor the convergence of the iterative solver, we need bounds on the true residual. Using such bounds, we can derive stopping criteria for the flexible Krylov solvers for shifted systems with inexact preconditioning. To do this, we first derive bounds on the norm of inexact residual $\tilde{r}_m(\sigma_j)$ and a bound on the difference between the true and the inexact residual ${\lVert r_m(\sigma_j)-\tilde{r}_m(\sigma_j) \rVert_2}$. A simple application of the triangle inequality for vector norms, leads us to the desired bounds on the true residual

The inexact residual $ \tilde{r}_m(\sigma_j)  $ defined as $$\tilde{r}_m(\sigma_j)  {\stackrel{\text{def}}{=}}V_{m+1}\left(\beta e_1 - \bar{H}_m(\sigma_j;T_m)y_m(\sigma_j)\right)$$ The expression for$ \tilde{r}_m(\sigma_j)$ is similar to the exact residual $\tilde{r}_m(\sigma_j)$ ignoring the error due to early termination of the inner iterative solver, i.e., $P_my_m(\sigma_j)$. It is easy to verify that ${\lVert \tilde{r}_m(\sigma_j)  \rVert_2} = {\lVert \beta e_1 - \bar{H}_m(\sigma_j;T_m)y_m(\sigma_j) \rVert_2} $. We now derive an expression for the norm of the difference between the true and the inexact residuals,

$$\begin{aligned}
{\lVert r_m(\sigma_j)-\tilde{r}_m(\sigma_j) \rVert_2} & =  {\lVert P_my_m(\sigma_j) \rVert_2} =  {\lVert \sum_{k=1}^m e_k^Ty_m(\sigma_j)p_k \rVert_2}  \nonumber\\
                    & \leq     \sum_{k=1}^m|e_k^Ty_m(\sigma_j)|{\lVert p_k \rVert_2}  \nonumber \\
                    & \leq    \varepsilon\sum_{k=1}^m|e_k^Ty_m(\sigma_j)| = \varepsilon \lVert y_m(\sigma_j)\rVert_1 \label{eqn:resdiff}\end{aligned}$$

Finally, the norm of the true residual $r_m(\sigma_j)$ can be bounded using the following relation $$\begin{aligned}
{\lVert r_m(\sigma_j) \rVert_2} \leq & \quad {\lVert r_m(\sigma_j)-\tilde{r}_m(\sigma_j) \rVert_2} + {\lVert \tilde{r}_m(\sigma_j) \rVert_2}  \nonumber \\
            \leq &  \quad \varepsilon {\lVert y_m(\sigma_j) \rVert_{1}}  + {\lVert \beta e_1 - \bar{H}_m(\sigma_j;T_m)y_m(\sigma_j) \rVert_2}  \label{eqn:trueresidual}\end{aligned}$$

This bound on the true residual, gives us a convenient expression to monitor the convergence of the iterative solver for each system, corresponding to a given shift $\sigma_j$.

We can also derive specialized results for the flexible FOM/GMRES for shifted systems with inexact preconditioning. The approach used is and argument similar to [@simoncini2003theory Proposition 4.1]. Let $r^\text{fom}_m(\sigma_j) {\stackrel{\text{def}}{=}}b - (K+\sigma_j M)\tilde{Z}_my_m^\text{fom}(\sigma_j)$ and $r^\text{gmres}(\sigma_j) {\stackrel{\text{def}}{=}}b - (K+\sigma_j M)\tilde{Z}_my_m^\text{gmres}(\sigma_j)$ be the true residual, respectively resulting from the flexible FOM/GMRES for shifted systems. We have the following error bounds

$$\begin{aligned}
{\lVert V_m^*r_m^\text{fom}(\sigma_j) \rVert_2} \leq & \quad \varepsilon {\lVert y_m^\text{fom}(\sigma_j) \rVert_{1}} \\ 
{\lVert \left(V_{m+1}\bar{H}_m(\sigma_j;T_m)\right)^*r_m^\text{gmres} \rVert_2} \leq & \quad \varepsilon {\lVert \bar{H}_m(\sigma_j;T_m) \rVert_2} {\lVert y_m^\text{fom}(\sigma_j) \rVert_{1}}  \end{aligned}$$

One of main results of the paper [@simoncini2003theory] is that they provide theory for why the residual norm due inexact preconditioning can be allowed to grow at the later outer iterations. In particular, they provide computable bounds for the monitoring the outer Krylov solver residual when the termination criteria for the inner preconditioning is allowed to change at each iteration, from which efficient termination criteria can be derived. We have not pursued this issue and the reader is referred to [@simoncini2003theory] for further details.

Application to Oscillatory Hydraulic Tomography {#sec:application}
===============================================

In this section, we briefly review the application of Oscillatory Hydraulic Tomography and the Geostatistical approach for solving the resulting inverse problem.

The Forward Problem {#sec:forward}
-------------------

The equations governing ground water flow through an aquifer for a given domain $\Omega$ with boundary $\partial \Omega = \partial \Omega_D \cup \partial \Omega_N, \partial \Omega_D \cap \partial \Omega_N = \emptyset$ are given by,

$$\begin{aligned}
 \label{eqn:timedomain}
S_s({\textbf{x}}) \frac{\partial \phi({\textbf{x}},t)}{\partial t} - \nabla \cdot \left(K({\textbf{x}}) \nabla \phi({\textbf{x}},t)\right) & =   q({\textbf{x}},t), & {\textbf{x}}&\in \Omega \\ \nonumber
\phi({\textbf{x}},t) & = 0, & {\textbf{x}}& \in \partial \Omega_D  \\ \nonumber 
\nabla \phi({\textbf{x}},t) \cdot \textbf{n} & =  0, & {\textbf{x}}&\in \partial \Omega_N \\ \nonumber \end{aligned}$$

where $S_s({\textbf{x}})$ \[L$^{-1}$\] represents the specific storage and $K ({\textbf{x}})$ \[L/T\] represents the hydraulic conductivity. In the case of one source oscillating at a fixed frequency $\omega$ \[radians/T\] , $q({\textbf{x}},t)$ is given by $$\label{eq:oscillation}
 q({\textbf{x}},t) = Q_0\delta({\textbf{x}}-{\textbf{x}}_s) \cos(\omega t)$$ To model periodic simulations, we will assume the source to be a point source oscillating at a known frequency $\omega$ and peak amplitude $Q_0$ at the source location ${\textbf{x}}_s$. In the case of multiple sources oscillating at distinct frequencies, each source is modeled independently with its corresponding frequency as in , and then combined to produce the total response of the aquifer.

Since the solution is linear in time, we assume the solution (after some initial time has passed) can be represented as $$\label{eqn:measurementequation}
 \phi({\textbf{x}},t) = \Re(\Phi({\textbf{x}}) \exp(i\omega t) )$$ where $\Re(\cdot)$ is the real part and $\Phi({\textbf{x}})$ is known as the phasor, is a function of space only and contains information about the phase and amplitude of the signal. Assuming this solution, the equations  in the phasor domain are, $$\begin{aligned}
 \label{eqn:phasor1}
- \nabla \cdot \left(K({\textbf{x}}) \nabla \Phi({\textbf{x}})\right) + i\omega S_s({\textbf{x}})  \Phi({\textbf{x}})  = &  \quad Q_0\delta({\textbf{x}}-{\textbf{x}}_s), & {\textbf{x}}\in \Omega \\ \nonumber
\Phi({\textbf{x}})  = & \quad 0, &\quad  {\textbf{x}}\in \partial \Omega_D \\ \nonumber
\nabla \Phi({\textbf{x}}) \cdot \textbf{n}  = & \quad 0, &\qquad {\textbf{x}}\in \partial \Omega_N \nonumber\end{aligned}$$ The differential equation  along with the boundary conditions are discretized using FEniCS [@LoggMardalEtAl2012a; @LoggWells2010a; @LoggWellsEtAl2012a] by using standard linear finite elements. Solving it for several frequencies results in system of shifted equations of the form $$\label{eqn:genshifted}
  \left( K + \sigma_j M \right) x_j  = b \qquad j=1,\dots,n_f$$ where, $K$ and $M$ are the stiffness and mass matrices, respectively, that arise precisely from the discretization of .

The Geostatistical Approach
---------------------------

The Geostatistical approach (described in the following papers [@kitanidis1995quasi; @kitanidis2010bayesian; @kitanidis2007on]) is one of the prevalent approaches for solving stochastic inverse problems. The idea is to represent the unknown field as the sum of a few deterministic low-order polynomials and a stochastic term that models small-scale variability. Inference from the measurements is obtained by invoking the Bayes’ theorem, through the posterior probability density function which is the product of two parts - likelihood of the measurements and the prior distribution of the parameters. Let $s({\textbf{x}}) \in \mathbb{R}^{N_s}$ be the function to be estimated, here the log conductivity, and let it be modeled by a Gaussian random field. After discretization, the field can be written as $s \sim {\cal{N}}(X\beta,Q)$. Here $X$ is a matrix of low-order polynomials, $\beta$ are a set of drift coefficients to be determined and $Q$ is a covariance matrix with entries $Q_{ij} = \kappa({\textbf{x}}_i,{\textbf{x}}_j)$, and $\kappa(\cdot,\cdot)$ is a generalized covariance kernel [@christakos1984problem]. The measurement equation can be written as, $$y = h(s) + v, \qquad v \sim {\cal{N}}(0,R)$$ where $y \in \mathbb{R}^{N_y}$ represents the noisy measurements and $v$ is a random vector of observation error with mean zero and covariance matrix $R$. The matrices $R$, $Q$ and $X$ are part of a modeling choice and more details to choose them can be obtained from the following references [@kitanidis1995quasi]. The operator $h: \mathbb{R}^{N_s}\rightarrow \mathbb{R}^{N_y}$ is known as the parameter-to-observation map or [*measurement operator*]{}, with entries that are the coefficients of the oscillatory terms in the expression, $$\label{eqn:measurement}
 \int_\Omega \Re\left\{e^{i\omega t}\Phi({\textbf{x}}) \delta({\textbf{x}}-{\textbf{x}}_i)\right\} d{\textbf{x}}$$ where ${\textbf{x}}_i$, is the location of the measurement sensor and $i=1,\dots,n_y$, where $n_y$ is the number of measurement locations. At each measurement location, two coefficients are measured for every frequency. In all, we have $N_y=2n_fn_y$ measurements, where $n_f$ is the number of frequencies.

Following the geostatistical method for quasi-linear inversion [@kitanidis1995quasi], we compute $\hat{s}$ and $\hat\beta$ corresponding to the maximum-a-posteriori probability which is equivalent to computing the solution to a weighted nonlinear least squares problem. To solve the optimization problem, the Gauss-Newton algorithm is used. Starting with an initial estimate for the field $s_0$, the procedure is described in algorithm \[alg:quasi\].

Compute the $N_y \times N_s$ Jacobian $J$ as, $$J_k = {\left.\frac{\partial{h}}{\partial{s}}\right|_{s = {s}_k}}$$ Solve the system of equations,

$$\label{eq:inversion}
\left( \begin{array}{cc} 
        J_k Q J_k^T  + R & J_kX \\
     \left(J_k X\right)^T & 0 
       \end{array}
       \right) 
       \left( \begin{array}{c}{\xi_{k+1}}\\ {\beta_{k+1}} \end{array} \right) = 
\left( \begin{array}{c} y - h({s}_k) + J_k{s}_k \\ 0 \end{array} \right)$$

The update $s_{k+1}$ is computed by, $$s_{k+1} = X \beta_{k+1} + Q J_k^T \xi_{k+1}$$

Repeat steps $1-3$ until the desired tolerance has been reached. (If necessary, add a line search).

Algorithm \[alg:quasi\] requires, at each iteration, computation of the matrices $QJ_k^T$ and $J_kQJ_k^T$. Since the prior covariance matrix $Q$ is dense, a straightforward computation of $QJ_k^T$ can be performed in ${\cal{O}}(N_yN_s^2)$. However, for fine grids, i.e., when the number of unknowns $N_s$ is large, storing $Q$ can be expensive in terms of memory and computing $QJ_k^T$ can be computationally expensive. For regular equispaced grids and covariance kernels that are stationary or translation invariant, an FFT based method can be used to reduce the storage costs of the covariance matrix $Q$ to ${\cal{O}}(N_s)$ and cost of matrix-vector product to ${\cal{O}}(N_s\log N_s)$. For irregular grids, the Hierarchical matrix approach can be used to reduce the storage costs and cost of approximate matrix-vector product to ${\cal{O}}(N_s\log N_s)$ for a wide variety of covariance kernels [@saibaba2012efficient]. Thus, in either situation, the cost for computing $QJ_K^T$ can be done in ${\cal{O}}(N_yN_s\
\log N_s)$ and the cost of computing $J_kQJ_k^T$ is ${\cal{O}}(N_sN_y\log N_s + N_sN_y)$.

Sensitivity Matrix computation {#sec:sensitivity}
------------------------------

Computing the Jacobian matrix $J_k$ at each iteration is often an expensive step. Although explicit analytical expressions for the entries are nearly impossible, several approaches exist. One simple approach is to use finite differences, but this approach is expensive because it requires as many $N_s+1$ runs of the forward problem, i.e. one more than the number of parameters to be estimated. For large problems and on finely discretized grids, the number of unknowns can be quite large and so this procedure is not feasible.

To reduce the computational cost associated with calculating the sensitivity matrix we use the adjoint state method (see for example, [@sun1990coupled]). This approach is exact and is computationally advantageous when the number of measurements is far smaller than the number of unknowns. For a complete derivation of the adjoint state equations for oscillatory hydraulic tomography, refer to [@cardiff2012multi]. For the type of measurements described in , the entries of the sensitivity matrix can calculated by the following expression for $j=1,\dots,N_s$ $$\label{eqn:sensitivity}
 \frac{\partial h}{\partial s_j} = \int_\Omega \Re\left\{ e^{i\omega t}\left( \left[i \omega \frac{\partial S_s({\textbf{x}})}{\partial s_j} \Phi - \frac{\partial Q_0}{\partial s_j} \right] \Psi_{\omega} + \frac{\partial K({\textbf{x}})}{\partial s_j} \nabla \Phi \cdot \nabla \Psi_{\omega}\right)  \right\}d{\textbf{x}}$$ Since at measurement location corresponding to each frequency, two measurements are obtained from the coefficients of the oscillatory terms, so the Jacobian matrix has $N_y\times N_s$ entries where, $N_y = 2n_fn_y$. Here, $\Psi_{\omega}$ is the known as the [*adjoint solution*]{} that depends on the measurement location ${\textbf{x}}_m$ and the forcing frequency $\omega$. It satisfies the following system of equations $$\begin{aligned}
\label{eqn:adjoint} 
 - \nabla \cdot \left(K \nabla \Psi_{\omega } \right) +  i \omega S_s\Psi_{ \omega}  = & \quad- \delta({\textbf{x}}-{\textbf{x}}_m), &\quad {\textbf{x}}\in \Omega \\ \nonumber 
\Psi_{\omega}  = & \quad 0, &\quad {\textbf{x}}\in \partial \Omega_D  \\ \nonumber
\nabla \Psi_{\omega}({\textbf{x}}) \cdot \textbf{n}  = & \quad 0, &\quad {\textbf{x}}\in \partial \Omega_N \nonumber\end{aligned}$$ where, ${\textbf{x}}_m$ is the measurement location and $\omega$ is the particular frequency. The procedure for calculating the sensitivity matrix can thus be summarized as follows.

1\. For a given field ${s}$, solve the forward problem for $\Phi$. 2. For each measurement and frequency $\omega$, solve the adjoint problem for $\Psi_\omega$. 3. Compute the integral in  to calculate the sensitivity.

Since  is evaluated for all $s_j$ for each measurement, the adjoint state method requires only $N_y+1$ forward model solves to compute the sensitivity matrix. Thus, when the number of measurements is far fewer than the number of unknowns, the adjoint state method provides a much cheaper alternative for computing the entries of the Jacobian matrix. This is typically the case in hydraulic tomography, where having several measurement locations is infeasible because it requires digging new wells.

Further, we realize that equation  takes the same form as equation  for multiple frequencies. Thus, we can use the algorithms developed in section \[sec:krylov\] to solve the system of equations  for as many right hand sides as measurements. It is possible to devise algorithms for multiple right hand sides in the context of shifted systems [@meerbergen2010lanczos; @darnell2008deflated] but we will not adopt this approach.

Numerical Experiments and Results {#sec:numerical}
=================================

We present numerical results for the Krylov subspace solvers and its application to OHT. As mentioned before, we use the FEniCS software [@LoggMardalEtAl2012a; @LoggWells2010a; @LoggWellsEtAl2012a] to discretize the appropriate partial differential equations. We use the Python interface to FEniCS, with uBLASSparse as the linear algebra back-end. For the direct solvers we use SuperLU [@superlu99] package that is interfaced by Scipy whereas for the iterative solver we use an algebraic multigrid package PyAMG [@BeOlSc2011], with smoothed aggregation along with BiCGSTAB iterative solver. In the following sections, for brevity, we only most results for FOM solver but we observed similar results for the GMRES method as well. This is also suggested by the result in proposition \[prop:gmres\].

Krylov subspace solver {#sec:solverexperiments}
----------------------

In this section, we present some of the results of the algorithms that we have described in section \[sec:krylov\]. We now describe the test problem that we shall use for the rest of the section. We consider a $2$D aquifer in a rectangular domain with Dirichlet boundary conditions on the boundaries. For the log-conductivity field $\log K({\textbf{x}})$, we consider a random field generated using an exponential covariance kernel $\kappa({\textbf{x}},\textbf{y}) = 4\exp(-2{\lVert {\textbf{x}}-\textbf{y} \rVert_2}/L)$ using the algorithm described in [@dietrich1993fast]. Other parameters used for the model problem are summarized in table \[tab:parameters\]. We choose $200$ frequencies evenly spaced between the minimum and maximum frequencies, which results in $200$ systems each of size $90601$.

  Definition                 Parameters              Values
  -------------------------- ----------------------- -------------------------------------
  Aquifer length             L (m)                   500
  Specific storage           $\log S_s$ (m$^{-1}$)   $-11.52$
  Mean conductivity          $\mu(\log K)$ (m/s)     $ -11.02$
  Variance of conductivity   $\sigma^2(\log K) $     $1.42$
  Frequency range            $\omega$ ($s^{-1}$)     $[\frac{2\pi}{600},\frac{2\pi}{3}]$

  : Parameters Chosen For Test Problem[]{data-label="tab:parameters"}

First, we motivate the need for multiple preconditioners to solve the shifted system of equations . We begin by looking at the number of iterations taken by restarted GMRES without a preconditioner and using a single preconditioner. We choose a preconditioner of the form $K+\tau M$, for five different values of $\tau$. For illustration purposes, we use a direct solver to invert the preconditioned systems. These values represent the minimum frequency, the average frequency and the maximum frequency in the parameter range. The number of iterations corresponding to restarted GMRES ($30$) for each of the system is computed and displayed in figure \[fig:itercountsingleprecond\]. We observe that the number of iterations corresponding to the systems increases as the frequency of the system decreases. When we use a preconditioner $K + \tau M$, the systems with frequencies nearby $|\tau|$ converge rapidly. However, systems with frequencies further from $|\tau|$ converge more slowly, the further away they are from the frequency of the preconditioned system $|\tau|$. This is consistent with the analysis in section \[sec:geneigen\] and in particular, proposition \[prop:fom\]. Thus, no single preconditioner effectively preconditions all the systems in the given frequency range. Not surprisingly, choosing $|\tau|$ in the center of the frequency range seems to be the best choice. Thus, in order to make the iterative method competitive, we consider using multiple preconditioners to solve the shifted system .

![(left) The log conductivity field that we use for the test problem with $90601$ grid points, and (right) iteration count for restarted GMRES ($30$) for unpreconditioned case and also with single preconditioners of the form $K + \tau M$, with $|\tau| \in \{ \frac{2\pi}{600},\frac{2\pi}{6},\frac{2\pi}{3} \}$. Results indicate that for the particular choices made, $|\tau| = \frac{2\pi}{6}$ which is roughly at the center of the frequency range, performs best.[]{data-label="fig:itercountsingleprecond"}](logcond "fig:") ![(left) The log conductivity field that we use for the test problem with $90601$ grid points, and (right) iteration count for restarted GMRES ($30$) for unpreconditioned case and also with single preconditioners of the form $K + \tau M$, with $|\tau| \in \{ \frac{2\pi}{600},\frac{2\pi}{6},\frac{2\pi}{3} \}$. Results indicate that for the particular choices made, $|\tau| = \frac{2\pi}{6}$ which is roughly at the center of the frequency range, performs best.[]{data-label="fig:itercountsingleprecond"}](itercount.png "fig:")

We choose preconditioners of the form $K + \tau_k M , k = 1,\dots,m$, where $m$ is the maximum dimension of the Arnoldi iteration. From figure \[fig:itercountsingleprecond\], it is clear that the systems with smaller frequencies converge slower, so we choose the values of $|\tau_k|$ that are distributed closer to the origin. In particular, we choose the values of $\tau_k$ such that they are evenly spaced on a log scale in the domain $\omega \in [\frac{2\pi}{600},\frac{2\pi}{3}]$. For example, for $n_p = 5$, the distribution of $|\tau|$ is illustrated in figure \[fig:multipletau\]. Now, let the possible values that $\tau$ can take be labeled as $\bar{\tau} = \{ \bar{\tau}_1,\dots,\bar{\tau}_{n_p} \}$ where, $n_p$ is the number of distinct number of preconditioner frequencies. Then, the first $m_1$ values of $\tau_k$ are assigned $\bar{\tau}_1$, the next $m_2$ values of $\tau_k$ are assigned $\bar{\tau}_2$ and so on. We also have $m = \sum_{k=1}^{n_p}m_k$. We pick $m_k = m/n_p$. If the algorithm has not converged in $m$ iterations, we restart using the method in section \[sec:restarting\] if we are using a direct solver as the preconditioner. Else, we recycle the same sequence of preconditioners. We implemented both algorithms to invert the preconditioner matrices - using direct solver and using an iterative solver which is an algebraic multigrid preconditioned BiCGSTAB. Using $n_p = 5$ and $m=40$ along with the scheme to choose the preconditioner frequencies described above, we observed that the number of iterations (and hence, matrix-vector products) in both cases were less than $40$.

Finally, we present the comparison in terms of the run time of our algorithm compared to solving each system using a direct solver. The results are presented in figure \[fig:timing\]. In the plots, “Direct” implies that every system is solved individually by a direct solver. “Flexible” algorithm uses $5$ different preconditioners (see figure \[fig:multipletau\] with a direct solver for inverting the preconditioners, and solves the FOM subproblem, whereas “Inexact” uses the preconditioners and inverts the preconditioners using an iterative solver. We see that both the “Flexible” and “Inexact” algorithms outperform the “Direct” approach even for a small number of frequencies. A relative tolerance of ${\lVert r_m(\sigma_j) \rVert_2} / {\lVert r_0(\sigma_j) \rVert_2} \leq 10^{-10}$ was used as stopping criteria for all systems $j=1,\dots,n_f$ and all systems converged within $40$ iterations. Although, the “Inexact” algorithm seems to behave nearly independent of number of frequencies, its runtime is longer than the “Direct” approach. This is due to the fact that the PyAMG solver requires an additional $17$ matvecs on average per inner iteration, totaling $706$ matvecs with $K+\tau_kM$ with $k=1,\dots,m$ and $m=40$.

![Comparison of time for the different algorithms. “Direct” implies that every system is solved individually by a direct solver. “Flexible” algorithm uses $5$ different preconditioners (see figure \[fig:multipletau\] with a direct solver for inverting the preconditioners, and solves the FOM subproblem, whereas “Inexact” uses the preconditioners and inverts the preconditioners using an iterative solver. We see that both the “Flexible” and “Inexact” algorithms outperform the “Direct” approach even for a small number of frequencies. The system size is $90601$. []{data-label="fig:timing"}](timecomparison)

In figures \[fig:reserrordirect\] and \[fig:reserrorflexible\], we plot the residuals and error as a function of the frequency of the system. A direct solver was used for figure \[fig:reserrordirect\], whereas an iterative solver (algebraic multigrid preconditioned BiCGSTAB) was used in figure \[fig:reserrorflexible\] with a stopping criterion $\varepsilon = 10^{-12}$ (refer to section \[sec:inexact\]. A relative tolerance of ${\lVert r_m(\sigma_j) \rVert_2} / {\lVert r_0(\sigma_j) \rVert_2} \leq 10^{-10}$ was used for all systems $j=1,\dots,n_f$ and all systems converged within 40 iterations. The behavior of residual and the error is quite similar but this is to be expected.

As discussed in section \[sec:inexact\], when an inexact preconditioner is used, the flexible Arnoldi relation is no longer exact and the true residual $r_m(\sigma_j)$ and the inexact residual $\tilde{r}_m(\sigma_j)$ are no longer exactly equal. In fact, the error between them can be bounded by the relation . In figure \[fig:reserrorflexible\], we compare the difference between the true residual and the inexact residuals with the predicted bound $\varepsilon{\lVert y_m(\sigma_j) \rVert_{1}}$. We see in figure \[fig:resdiff\] that the bound is fairly accurate. The stopping tolerances were chosen to be $\varepsilon = 10^{-9},10^{-11}, 10^{-12}$. In fact, for larger stopping tolerances for the inner solver, the outer solver did not converge.

Application: Tomographic reconstruction
---------------------------------------

The objective is now to determine an known conductivity field $K({\textbf{x}})$ from discrete measurements of the head $\phi$ obtained from several pumping tests performed with multiple frequencies. Since the conductivity field needs to be positive, so that the forward problem is well-posed, we consider a log-transformation $s = \log K$. The “true" field is taken to be that in figure \[fig:trueloc\] which is a scaled version of Franke’s function [@franke1979critical]. We choose the covariance matrix $Q$ to have entries $Q_{ij} = \kappa({\textbf{x}}_i,{\textbf{x}}_j)$, corresponding to an exponential covariance kernel $$\kappa({\textbf{x}},{\textbf{y}}) = \exp\left(-\frac{4{\lVert {\textbf{x}}-{\textbf{y}} \rVert_2}}{L}\right)$$ where, $L$ is the length of the domain. We also choose $R = \eta^2 I$ and $X = [1,\dots,1]^T$. We did not try to optimize the choice of covariance kernels to get the best possible reconstruction. Our goal is to study the associated computational costs. The size of our problem is chosen to be $10201$ discretization points. We assume no noise in our measurements and choose $\eta =  10^{-6}$.

The measurements are obtained by taking as the true log conductivity field, the field in figure \[fig:trueloc\]. Then, the phasor is calculated by solving equations  with the source location and measurement locations given in figure  \[fig:trueloc\]. The measurements are collected for each frequency and two pieces of information are recored, the coefficients of the sine and cosine terms in equation . The inverse problem is then solved using these measurements. The frequency range that is chosen is $\omega \in [2\pi/150,2\pi/30]$. We pick $n_p = 5$ evenly spaced in log-scale in this particular frequency range and set $m=40$. All systems converged in $40$ iterations.

The time for computing the Jacobian is listed in figure \[fig:jacobian\]. For the iterative solver, we use the flexible FOM solver using direct solver as preconditioner. The relative stopping tolerance we used was $10^{-10}$. Since the cost to solve systems with multiple frequencies is nearly the same as the cost to solve a single system, the time for building the Jacobian is, more or less, independent of the number of frequencies. However, when a direct solver is used to independently solve the systems for multiple frequencies, the cost for constructing the Jacobian scales linearly with the number of frequencies. This results in significant reduction in the cost for solving the inverse problem, since constructing the Jacobian is the most expensive part of solving the inverse problem. The disparity in the computation times for the Jacobian between the direct approach and the iterative procedure is exacerbated further, with larger problem sizes resulting from finer discretizations.

Finally, we compare the error in the the reconstruction with multiple frequencies. Table \[table:frequencies\] lists the $L^2$ error in the reconstruction. We report two errors - the first being the $L^2$ error in the entire domain, the second being the $L^2$ error in the area enclosed by the measurement wells. While increasing the number of frequencies improves the error in the entire domain, as well as in the region enclosed by the measurement wells. This is the primary motivation for using multiple frequencies in the inversion. However, beyond a point, the addition of frequencies does not seem to reduce the error. This might be because there is no additional information that is obtained from the addition of measurements with these frequencies, and to further improve estimation accuracy, one would need to introduce more stimulation and observation points [@cardiff2012multi].

   $N_f$   Total error   Error within box
  ------- ------------- ------------------
    $1$     $0.3794$         $0.0511$
    $5$     $0.3379$         $0.0352$
   $10$     $0.3264$         $0.0337$
   $20$     $0.3180$         $0.0328$

  : L$^2$ error due to the reconstruction. We report two errors - the first being the $L^2$ error in the entire domain, the second being the $L^2$ error in the area enclosed by the measurement wells. []{data-label="table:frequencies"}

Conclusions
===========

We have presented a flexible Krylov subspace algorithm for shifted systems of the form  that uses multiple shifted preconditioners of the form $K+\tau M$. The values of $\tau$ are chosen in order to improve convergence of the solver for all the shifted systems. The number of preconditioners chosen varies based on the distribution of the shifts. A good rule of thumb is that the systems having shift $\sigma$ will converge faster if there a preconditioner with shift $\tau$ that is nearby $\sigma$. When the size of the linear systems is much larger, direct solvers are much more expensive. In such cases, preconditioning would be done using iterative solvers. The error analysis in section \[sec:inexact\] provides insight into monitor approximate residuals without constructing the true residuals. One can naturally extend the ideas in this paper to systems with multiple shifts and multiple right hand sides using either block or deflation techniques.

We applied the flexible Krylov solver to an application problem that benefited significantly from fast solvers for shifted systems. In particular, oscillatory hydraulic tomography is a technique for aquifer characterization. However, since drilling observation wells to obtain measurements is expensive, one of the advantages of oscillatory hydraulic tomography is obtaining more informative measurements by pumping at different frequencies using the same pumping locations and measurement wells. In future studies we aim to study more realistic conditions for tomography, including a joint inversion for storage and conductivity. This would be ultimately beneficial to the practitioners. We envision that fast solvers for shifted systems would be beneficial for rapid aquifer characterization using oscillatory hydraulic tomography.

Acknowledgments
===============

The research in this work was funded by NSF Award 0934596, “CMG Collaborative Research: Subsurface Imaging and Uncertainty Quantification” and by NSF Award 1215742, “ Collaborative Research: Fundamental Research on Oscillatory Flow in Hydrogeology.” The authors would also like to thank their collaborators Michael Cardiff and Warren Barrash for useful discussions and the two anonymous reviewers for their insightful comments.

[10]{}

W. N. Bell, L. N. Olson, and J. B. Schroder. : [A]{}lgebraic multigrid solvers in [Python]{} v2.0, 2011. Release 2.0.

M. Cardiff, T. Bakhos, P.K. Kitanidis, and W. Barrash. Multi-frequency oscillatory hydraulic tomography: The use of steady-state periodic signals for sensitivity analysis and inversion. , in review.

M. Cardiff and W. Barrash. 3-[D]{} transient hydraulic tomography in unconfined aquifers with fast drainage response. , 47(12):W12518, 2011.

G. Christakos. On the problem of permissible covariance and variogram models. , 20(2):251–265, 1984.

D. Darnell, R.B. Morgan, and W. Wilcox. Deflated [GMRES]{} for systems with multiple shifts and multiple right-hand sides. , 429(10):2415–2434, 2008.

B.N. Datta and Y. Saad. methods for large [S]{}ylvester-like observer matrix equations, and an associated algorithm for partial spectrum assignment. , 154:225–244, 1991.

James W. Demmel, Stanley C. Eisenstat, John R. Gilbert, Xiaoye S. Li, and Joseph W. H. Liu. A supernodal approach to sparse partial pivoting. , 20(3):720–755, 1999.

CR Dietrich and GN Newsam. A fast and exact method for multidimensional gaussian stochastic simulations. *Water Resources Research*, 290 (8):0 2861–2869, 1993.

R. Franke. A critical comparison of some methods for interpolation of scattered data. Technical report, DTIC Document, 1979.

R.W. Freund. Solution of shifted linear systems by quasi-minimal residual iterations. , pages 101–121, 1993.

A. Frommer and U. Gl[ä]{}ssner. Restarted [GMRES]{} for shifted linear systems. , 19(1):15–26, 1998.

G.H. Golub and C.F. Van Loan. , volume 3. Johns Hopkins University Press, 1996.

G. Gu, X. Zhou, and L. Lin. A flexible preconditioned [A]{}rnoldi method for shifted linear systems. , 25(5):522-530, 2007.

T.J.R. Hughes. The finite element method: linear static and dynamic finite element analysis. , 2012.

P. K. Kitanidis. Quasilinear geostatistical theory for inversing. , 31(10):2411–2419, 1995.

P. K. Kitanidis. , volume Geophysical Monograph 171, pages 19–30. AGU, Washington, D. C., 2007.

P. K. Kitanidis. , pages 71–85. John Wiley & Sons, Ltd, 2010.

Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. . Springer, 2012.

Anders Logg and Garth N. Wells. Dolfin: Automated finite element computing. , 37(2), 2010.

Anders Logg, Garth N. Wells, and Johan Hake. , chapter 10. Springer, 2012.

K. Meerbergen. The solution of parametrized symmetric linear systems. , 24(4):1038–1059, 2003.

K. Meerbergen and Z. Bai. The [Lanczos]{} method for parameterized symmetric linear systems with multiple right-hand sides. , 31(4):1642–1662, 2010.

M. Popolizio and V. Simoncini. Acceleration techniques for approximating the matrix exponential. , pages 657–683, 2008.

Y. Saad. , volume 158. SIAM, 1992.

Y. Saad. A flexible inner-outer preconditioned [GMRES]{} algorithm. , 14:461–461, 1993.

Y. Saad. . Society for Industrial Mathematics, 2003.

A.K. Saibaba and P.K. Kitanidis. Efficient methods for large-scale linear inversion using a geostatistical approach. , 48(5):W05522, 2012.

V. Simoncini. Restarted full orthogonalization method for shifted linear systems. , 43(2):459–466, 2003.

V. Simoncini and D.B. Szyld. Theory of inexact [Krylov]{} subspace methods and applications to scientific computing. , 25(2):454–477, 2003.

V. Simoncini and D.B. Szyld. Recent computational developments in [Krylov]{} subspace methods for linear systems. , 14(1):1–59, 2007.

N.Z. Sun, W.W.G. Yeh, et al. Coupled inverse problems in groundwater modeling 1. sensitivity analysis and parameter identification. , 26(10):2507–2525, 1990.
