# Reduced order models for spectral domain inversion: Embedding into the continuous problem and generation of internal data

We generate data-driven reduced order models (ROMs) for inversion of the one and two dimensional Schrödinger equation in the spectral domain given boundary data at a few frequencies. The ROM is the Galerkin projection of the Schrödinger operator onto the space spanned by solutions at these sample frequencies. The ROM matrix is in general full, and not good for extracting the potential. However, using an orthogonal change of basis via Lanczos iteration, we can transform the ROM to a block triadiagonal form from which it is easier to extract q. In one dimension, the tridiagonal matrix corresponds to a three-point staggered finite-difference system for the Schrödinger operator discretized on a so-called spectrally matched grid which is almost independent of the medium. In higher dimensions, the orthogonalized basis functions play the role of the grid steps. The orthogonalized basis functions are localized and also depend only very weakly on the medium, and thus by embedding into the continuous problem, the reduced order model yields highly accurate internal solutions. That is to say, we can obtain, just from boundary data, very good approximations of the solution of the Schrödinger equation in the whole domain for a spectral interval that includes the sample frequencies. We present inversion experiments based on the internal solutions in one and two dimensions.

## Authors

• 4 publications
• 6 publications
• 5 publications
• 2 publications
• 5 publications
• ### Lippmann-Schwinger-Lanczos algorithm for inverse scattering problems

Data-driven reduced order models (ROMs) are combined with the Lippmann-S...
01/28/2021 ∙ by Vladimir Druskin, et al. ∙ 0

• ### Adaptive spectral decompositions for inverse medium problems

Inverse medium problems involve the reconstruction of a spatially varyin...
06/17/2020 ∙ by Daniel H. Baffet, et al. ∙ 0

• ### Reduced Order Model Approach to Inverse Scattering

We study an inverse scattering problem for a generic hyperbolic system o...
10/29/2019 ∙ by Liliana Borcea, et al. ∙ 0

• ### Spectral Methods - Part 2: A comparative study of reduced order models for moisture transfer diffusive problems

This paper explores in details the capabilities of two model reduction t...
05/16/2017 ∙ by Suelen Gasparin, et al. ∙ 0

• ### Analytic approximation of transmutation operators for one-dimensional stationary Dirac operators and applications to solution of initial value and spectral problems

A method for approximate solution of initial value and spectral problems...
06/08/2020 ∙ by Nelson Gutiérrez Jiménez, et al. ∙ 0

• ### Randomization for the Efficient Computation of Parametric Reduced Order Models for Inversion

Nonlinear parametric inverse problems appear in many applications. Here,...
07/12/2020 ∙ by Selin Aslan, et al. ∙ 0

• ### Accurate 3D frequency-domain seismic wave modeling with the wavelength-adaptive 27-point finite-difference stencil: a tool for full waveform inversion

Efficient frequency-domain Full Waveform Inversion (FWI) of long-offset/...
08/19/2021 ∙ by Hossein S. Aghamiry, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In this work we consider the following problem on a bounded domain in dimensions

 (1) −Δu(x;λ)+q(x)u(x;λ)+λu(x;λ) = 0  for   x∈Ω⊂Rd ∂u∂ν(x;λ) = g  for  x∈∂Ω

for the outward unit normal to the boundary and . That is, is the solution to the Schrödinger equation with Neumann data corresponding to spectral parameter , which here we assume to be positive real. We apply Neumann boundary functions , with each and read the corresponding Dirichlet data on the boundary for multiple values of , and are interested in determining from this data. The reconstruction of from boundary measurements is known to be an ill-posed problem and, depending on the physical setup can result in reconstructions with very poor resolution [2]. The map from to the data is highly nonlinear, and due to ill-posedness is difficult to invert unless regularization is added. As examples of literature on this subject, see [12],[3],[17].

A relatively new class of approaches employ ideas from model reduction theory to do direct inversion using reduced order forward models [7],[15],[8],[9],[16] with impressive numerical results. These works have grown out of earlier work on spectrally matched finite difference grids [14], where special grid steps were chosen so that the numerical solution matched functionals of the solution at the boundary in the spectral domain. The spectrally matched grid steps depended very weakly on the medium, and this allowed for the direct extraction of the unknown coefficients [10]. There were some extensions of this to higher dimensions [5], [4], [6], however, the grid itself was fundamentally one-dimensional, blocking the generation of truly higher dimensional direct inversion methods, except for special geometries. The more recent works use a more general reduced order model approach [7],[8],[9],[15],[16] .

Here we show how one can obtain spectrally matched finite difference operators with well chosen localized Galerkin basis functions. Using the basis functions, higher dimensional extensions are natural. This was essentially what was used in [15],[7] in the time domain. In this work we consider the spectral domain, and we show explicitly how to view the reduced order models (ROMs) through the spectral Galerkin equivalence mentioned above. The Galerkin equivalence of spectrally matched grids was known on the boundary from [13], however the basis in which the Galerkin model is the same as the finite difference operator everywhere was not known.

A main advantage first observed in [7] is that, like the grid steps, the basis functions appear to depend only very weakly on the medium. It is this crucial observation that allows one to use the ROM for inversion. Furthermore, by viewing the reduced order model as embedded into the continuous problem and using the reference medium basis functions, good approximations of the internal solutions can be produced from the boundary data. We will demonstrate this with numerical examples in one and two dimensions.

The paper is organized as follows. In Section 2 we state the reconstruction problem in general. In Section 3 we consider a one dimensional version with one source. In Section 3.1 we follow the Loewner framework [1] and describe how one obtains a Galerkin model from the data directly, in Section 3.2 we describe the orthogonalization and prove the equivalence everywhere with spectrally matched grids. In Section 3.3 we present one dimensional numerical examples along with the internal solution generation/reconstruction algorithm. In Section 4 we describe the same approach in higher dimensions; including the generation of the ROM from the data in Section 4.1 and the orthogonalization procedure in Section 4.2. Finally in Section 4.3 we present two dimensional numerical examples and the internal solution generation/reconstruction algorithm, and Section 5 contains some concluding remarks.

## 2. Problem statement.

Consider the problem (1) where we apply Neumann boundary functions , with each and denote by the solution to (1) corresponding to spectral value and source . We read data in the form of the matrix valued transfer function as depending on the spectral parameter

 Frl(λ):=∫∂Ωur(x;λ)gl(x)dσx.

If, for example, each Neumann function is an approximate delta function corresponding to a source/receiver location, for , the transfer function corresponds to reading the solution coming from the th ”source” at the th ”receiver”. We take spectral values , and define

 uri:=ur(x;bi),

the solution to (1) at with . Now let us assume that we are in possession of boundary data of the form

 Firl:=Frl(bi)=∫∂Ωurigl

and

 DFirl:=dFrldλ(λ)|λ=bi

where the derivative here is with respect to . The inverse problem is to determine from the data

 (2) {Firl,DFirl}   for  i=1,…,m  and  r,l=1,…,K.

Such data can be obtained from the time domain via a Fourier transform.

###### Remark 2.1.

As we will see below, from this data one can find a ROM which matches the data exactly. However, it is not at all crucial that the data be of the Hermite form (2). In fact, matching point values of at spectral points leads to a ROM with even better approximation properties and still yields a stable system [14], [13]. Other forms of spectral data are also possible [1],[13]. However, we consider only the Hermite data (2) in this work for simplicity.

## 3. One dimensional problems.

### 3.1. One dimensional Galerkin model from spectral data.

To begin, we follow the Loewner framework for model reduction [1]. Consider solving the one dimensional, one source version of (1),

 (3) −d2dx2u(x;λ)+q(x)u(x;λ)+λu(x;λ) = 0  for  x∈(0,1) −ddxu(0;λ)=1 ddxu(1;λ)=0

where we read the data and for as would be natural coming from the time domain. ( ln principle could be anywhere in

away from the Neumann eigenvalues of the Schrödinger operator ). One sees easily that (

3) has the variational form: Find such that

 (4) ∫10u′ϕ′+∫10quϕ+λ∫10uϕ=ϕ(0)

for all . We note that by Sobolev embedding and the trace theorem the values of are well defined for all test functions in . For any given value of , let us define our true spectral domain medium response by

 F(λ)=u(0;λ)

where is the solution to (4).

Consider the exact solutions to (4) where

 uj(x):=u(x,bj).

These form the subspace

 Vm=span{u1,…,um},

and we can consider the Galerkin solution to the variational problem: Find such that

 (5) ∫10u′Gϕ′+∫10quGϕ+λ∫10uGϕ=ϕ(0)

for any . Searching for the unknown coefficients for the solution and by setting , one obtains the system

 (6) (S+λM)→c=→F

where

, the vector

, and are the mass and stiffness matrices given by

 Mij=∫10uiuj

and

 Sij=∫10u′iu′j+∫10quiuj.

Using the variational formulation (5) with and (recall is the exact solution for ), we obtain that

 (7) Sij+biMij=F(bj)   for all i,j=1,…,m.

Using this, if we are in possession of data , we can actually directly reconstruct and . By reversing and and subtracting, for we have

 (bi−bj)Mij=F(bj)−F(bi)

or

 (8) Mij=F(bj)−F(bi)bi−bj

where we have used the symmetry of and . Imagine now that one were to have data corresponding to solution for some spectral point close to . Then by the above one has the formula

 ∫10uiuz=F(z)−F(bi)bi−z.

Taking , we have that , so that

 (9) Mii=∫10u2i=−F′(bi).

Multiplying (7) by , reversing and and subtracting, we get

 (bj−bi)Sij=bjF(bj)−biF(bi)

or

 (10) Sij=bjF(bj)−biF(bi)bj−bi.

Again by taking some spectral point close to and letting one obtains

 (11) Sii=(λF)′(bi).

Hence the formulas (8),(9),(10),(11), well known in the model reduction community [1], provide the mass and stiffness matrices directly from the data. Similar formulas can be derived for other forms of spectral data. Together and provide a reduced order model for (3).

The solution to (5) also has the property that at , as a function of it matches the data

 (12) {F(bj),F′(bj):j=1,…,m}

exactly. We show this in the following Lemma, which was previously known [1].

###### Lemma 3.1.

Let be the Galerkin solution to (5) for subspace generated by exact solutions corresponding to . Let given by

 (13) Fm(λ)=m∑i=1y2iλ−θi

be the unique rational Hermite interpolant of form (

13) to matching the data (12) with positive residues and negative poles. We then have

 uG(0;λ)=Fm(λ)

for all .

###### Proof.

First we show that the Galerkin solution yields a response at of the form (13) for all . Here we use the continuous inner product . Decomposing the solution into an orthonormal basis of Galerkin eigenpairs ,

 uG=m∑i=1⟨uG,wi⟩wi,

we have

 (14) uG(0;λ)=m∑i=1⟨uG,wi⟩wi(0).

The fact that is an eigenpair means that

 ⟨uG,wi⟩=−1μi(∫w′iu′G+∫qwiuG)

which from the variational formulation for gives

 ⟨uG,wi⟩=1μi⟨wi,λuG⟩−1μiwi(0).

Solving for we have

 ⟨uG,wi⟩=wi(0)λ−μi

which from (14) yields

 FG(λ):=uG(0;λ)=m∑i=1w2i(0)λ−μi.

Hence has the same form as in (13), with positive residues and negative poles. Recall that is the rational Hermite interpolant to at the points . Now, since the trial space for the Galerkin solution contains the exact solutions for , will be the exact solution at those spectral points, which means that at ,

 uG(0;bi)=F(bi)=Fm(bi).

We claim also matches the derivatives at the spectral points as a function of . The Galerkin solution for some near can be used as a test function for , so that one has

 FG(z)=uz(0)=∫u′iu′z+∫quiuz+bi⟨ui,uz⟩.

We can also use the variational formulation for this same Galerkin solution for with as a test function to obtain

 FG(bi)=ui(0)=∫u′zu′i+∫quzui+z⟨uz,ui⟩,

and by subtracting them we have

 FG(z)−FG(bi)=(bi−z)⟨uz,ui⟩.

Taking the limit as one obtains

 F′G(bi)=−limz→bi⟨uz,ui⟩=−⟨ui,ui⟩=F′(bi)

from (9). Hence is the rational Hermite interpolant to

at the given spectral points, and so by uniqueness of this Padé approximation,

 uG(0;λ)=Fm(λ)

for all . ∎

Spectral Galerkin models such as the above are difficult to work with since they tend to yield full matrices. In particular for inverse problems, it is not obvious how one may extract the coefficients, since we wouldn’t know the basis functions internally from the data. It is for these reasons we follow along the lines of [7],[16] and do a Lanczos orthogonalization to construct a new orthonormal basis for in which the operator becomes tridiagonal. In the new basis, is the identity from the orthonormality.

### 3.2. Lanzcos orthogonalization and equivalence with spectrally matched finite difference grids.

The original idea of spectrally matched grids [14] was that, given receiver data in the spectral domain, choose a rational approximation for this data, then find the three point finite difference stencil whose forward solution at the receiver matches exactly this rational approximation. This yielded a very special nonuniform grid and spectral convergence at the receiver. Higher order convergence for the forward model was sacrificed elsewhere in domain, but this was irrelevant for inversion. It was also shown at the time that the receiver response of the finite difference scheme was equivalent to that of a spectral Galerkin method, for various types of spectral data [13]. Just as above we saw how to generate a Galerkin model from data, one can generate a finite difference model from the same data.

Suppose we read the data (12) as before. Forgetting for now about the Galerkin method, there is a unique rational approximation to of the form (13) with positive residues and negative poles which gives this Hermite interpolant to , matching the data (12). From these positive residues and negative poles, by using what is essentially the Lanczos algorithm, one can find a continued fraction form for ,

 Fm(λ)=1^γ1λ+1γ1+1^γ2λ+⋯1γm−1+1^γmλ.

Then from there one can extract which define a three-point staggered difference scheme with tridiagonal matrix for which the approximated solution at is exactly [14]. That is, the rational approximation uniquely determines positive , such that solving the finite difference scheme

 (15) −1^γj(Uj+1−Ujγj−Uj−Uj−1γj−1)+λUj = 0   for  j=1,…m −U1−U0γ0=1, Um+1−Umγm=0

yields

 U1=Fm(λ).

Note that and are ghost points only (and ghost grid steps) and the nonzero Neumann condition on the left yields a nonzero right hand side in the first component. That is, this yields an dimensional matrix system of the form

 (Lm+λI)→U=→e1^γ1

for the unknown , with entries of given by (15) in terms of the and .

###### Remark 3.2.

When , the are primary and dual finite difference grid steps. For nonzero , finite difference scheme (15) can be transformed to a discrete Schrödinger form with a discrete Liouville transform [11]. For simplicity, we won’t do that here, since we don’t use the finite difference framework for reconstructions.

Consider again approximating the solution to (3) by a Galerkin method with the finite dimensional subspace , generated by the exact solutions at the spectral points . Since the finite difference grid does not easily extend to higher dimensions, we would rather work with Galerkin systems. Recall the Galerkin system (6)

 (S+λM)→c=→F

where and the vector contains the coefficients for solution

 uG=∑ciui.

Unfortunately, unlike the finite difference model, the stiffness and mass matrices and are full. However, one can in fact find a new basis for in which the Galerkin system is tridiagonal, and exactly everywhere the same as the symmetrization of the finite difference system (15). To do this we will do Lanczos orthogonalization. Define to be the projection of the delta function onto , that is, the unique element of which satisfies

 ⟨δm,w⟩=w(0)   for all   w∈Vm,

identified with the coefficient vector . Note that Consider now the basis where we let and

 B={→d,A→d,A2→d,…,Am−1→d}.

and we orthogonalize using Gram-Schmidt with respect to the mass matrix inner product

 ⟨→x,→y⟩M:=⟨M→x,→y⟩.

This will yield continuous orthogonality of the basis functions. Note that is symmetric with respect to this inner product. We get new orthogonalized basis for :

 (16) Vm=span{^u1(x),^u2(x),…^um(x)}.

Due to the orthogonalization, the new basis is somewhat localized, more refined near zero and coarsening out towards one, as spectrally matched finite difference grid steps do. See an example in Figure 1 for . Furthermore, is tridiagonal in the new basis, and since becomes the identity, this means that is also tridiagonal in the new basis.

###### Theorem 3.3.

If we change the basis to the orthogonalized (16) and form the Galerkin system in this new basis

 (17) (^S+λ^M)→^c=→^F

(with ) to solve for , this is the symmetrization of the finite difference system (15) for . In particular,

 uG=m∑i=1√^γiUi^ui(x).
###### Proof.

We first symmetrize (15) by setting

 Vj=√^γjUj

and and then multiplying each row by in the resulting equations for the . The resulting system is in the form

 (L+λI)→V=1√^γ1→e1,

where is now symmetric. We recall that for the given impedance of the form (13) (now multiplied by ), the symmetric tridiagonal is uniquely determined.

Consider now the system (17) in the new basis. The mass matrix will be identity due to orthogonality. The stiffness matrix will be tridiagonal due to the fact that this is a Lanczos process with polynomials in , ( becomes tridiagonal and therefore ). Furthermore, since is the projection of the delta function, it will be nonzero at zero, but from orthogonality of the new basis functions we have that

 ⟨^u1,^ui⟩=^ui(0)=0

for all . Hence the right hand side is nonzero only in the first component and is therefore equal to . This means that the system (17) is of the form

 (^S+λI)→^c=^u1(0)→e1

for tridiagonal symmetric. We also know from Lemma 3.1 and the construction of the system (15) that for all , which means that

 ^c1^u1(0)=1√^γ1V1

for all . By decomposing into the normalized eigenpairs of as was done in the proof of Lemma 3.1 except now using the standard Euclidean inner product in , one can calculate that

 ^c1=^u1(0)m∑i=1z2iλ−θi

for the eigenvalues and

the squares of the first components of the normalized eigenvectors of

. Similarly

 V1=1√^γ1m∑i=1y2iλ−μi

where the eigenvalues and the squares of the first components of the normalized eigenvectors of . By the above impedance equality for all and the normalization of the residues, one must have that the poles and residues are equal (and that the multipliers from right hand sides are the same). Therefore and are the same by the uniqueness of the inverse eigenvalue problem, and , from which the result follows. ∎

###### Remark 3.4.

From the above we also found that

 ^u1(0)=1√^γ1,

which coincides with being a projection of the delta function.

That is, the solution components of the difference scheme (15) can be interpreted as coefficients of the spectrally converging Galerkin solution with respect to the orthogonal basis (16). See Figure 2 for an example of an orthogonalized basis with its equivalent staggered finite difference grid. Note that the masses of the basis functions are concentrated between the primary grid steps and the next dual grid step, and they spread out in the same way that the grid steps coarsen. In Figure 3, we normalize each function in the basis to be equal to at its corresponding primary grid point, and plot it against the nodal basis for the same space. (The nodal basis is the unique basis of satisfying .) Note that our orthogonal basis resembles the nodal one, but is less oscillatory away from the corresponding grid point. In what follows below, we propose an inversion method using the Lanczos orthogonalized basis, which makes no use of the finite difference grid steps, and is easily generalizable to other geometries.

### 3.3. Internal data generation and inversion in one dimension.

Consider solving

 (18) −u′′+qu+λu = 0   on    (0,1) −u′(0) = 1 u(1) = 0

as a perturbation of the corresponding reference problem

 (19) −u′′0+λu0 = 0   on    (0,1) −u′0(0) = 1 u0(1) = 0.

Here we use the reference medium , however, any known reference medium should work the same. We note also that one can have either a Dirichlet or Neumann condition on the right endpoint . The one dimensional numerical experiments below were done with a Dirichlet condition, a Neumann condition yielded similar results.

Algorithm

1. Read data (here synthetically generated) , for some positive for the perturbed problem (18), and compute all of the corresponding solutions for for the reference problem.

2. Use (8, 9, 10, 11) to generate the stiffness and mass matrices and for Galerkin system for the perturbed problem, and similarly compute , for the reference problem.

3. Compute projections of the delta function at zero onto by solving for in

 M→d=→F

where components . These are the coefficients of projected in the basis for . Similarly, we calculate by solving

 M0→d0=→F0,

where is the reference data. This means that are the coefficients of the projected delta function in the reference basis for the space .

4. Set , and note that is symmetric with respect to the inner product , which is same the inner product on the corresponding continuous functions. Perform Lanczos orthogonalization with respect to this inner product on the basis

 {→d,A→d,A2→d,…,Am−1→d}

and let , be the Galerkin stiffness and mass matrices in this new basis . Note that is tridiagonal.

5. Similarly, set and perform Lanczos orthogonalization with respect to the inner product on the basis

 {→d0,A0→d0,A20→d0,…,Am−10→d0}

to obtain the corresponding orthognalized basis for .

6. Now choose any spectral value and solve , so that is the Galerkin approximation in to the solution of (18) . Since we don’t know from the data, approximate by

 ~u=∑ci^u0i.
7. Compute ; or reserve for another use.

In summary, the spectral data gives the Galerkin system ROM for (3) corresponding to spectral snapshots. Simultaneously, we can do the same for a reference medium. Lanczos orthogonalization localizes the spectral snapshots, forming a basis which depends very weakly on the medium. In Figure 4, we show an example of the orthogonalized basis functions for a piecewise cubic varying from up to and back down to . Here they are plotted against the orthogonalized basis functions for the reference medium, and appear indistinguishable. Therefore, although we do not have the internal spectral snapshots for the true medium, in the localized basis they are very close to those of the reference medium. For any given now, by solving the Galerkin system, we can get the coefficients, and this yields a good approximation to internal data. See Figure 5 for an example of an internal solution. One possibility is to use the internal solution to do inversion- for example by calculating . Note that for positive , is never zero. We try this simple approach in Figure 6 and Figure 7.

## 4. Higher dimensional problems.

### 4.1. Galerkin model from spectral data in higher dimensions.

Again we follow the Loewner framework. Consider the full problem (1) with given data (2). Assuming that each Neumann source function and , the variational formulation of (1) is: Find such that

 ∫Ω∇uri⋅∇ϕ+∫Ωquriϕ+bi∫Ωuriϕ=∫∂Ωgrϕ

for all . Generating a Galerkin subspace with these exact solutions

 VmK=span{{uri}i=1,…,m,r=1,…,K}

and using them as test functions, we have that

 (20) ∫Ω∇uri⋅∇ulj+∫Ωquriulj+bi∫Ωuriulj = ∫∂Ωgrulj (21) = Fjlr

for all and . That is,

 Sirjl+biMirjl=Fjlr

where and are the stiffness and mass matrices respectively. From this and the same derivation as in one dimension, one obtains them directly from the data

 (22) Mirjl=Fjlr−Filrbi−bj,
 (23) Miril=−DFilr,
 (24) Sirjl=bjFjlr−biFilrbj−bi,

and

 (25) Siril=(λFrl)′(bi).

### 4.2. Block Lanczos orthogonalization and generation of the block tridiagonal reduced order model.

Note the above and can be viewed as matrices with entries that are blocks. First we consider our data functional which is defined by

 [δF(ϕ)]l:=∫∂Ωϕgl   % for  l=1,…,K.

Similarly to what was done in one dimension, we find the projection of each component of onto . That is, we look for such that

 ∫ΩδlmKurj=∫∂Ωurjgl=Fjlr.

That is, we identify with its coefficients in the basis for ;

 δlmK=∑i=1,…,m;r=1,…,Kdliruri,

which again can be obtained directly by solving the matrix system of blocks

 M→d=F.

We now use this as our starting vector for block Lanczos orthogonalization with respect to the basis generated by powers of the matrix :

 B={→d,A→d,A2→d,…,Am−1→d}.

and we orthogonalize using Gram-Schmidt with respect to the mass matrix inner product

 ⟨→x,→y⟩M:=⟨M→x,→y⟩.

This will yield again continuous orthogonality of the basis functions. Note that is symmetric with respect to this inner product and will be block tridiagonal in the new basis for :

 (26) VmK=span{{^uri(x)}i=1,…,m;r=1,…,K}.

In this new basis, will be identity and will be block tridiagonal, yielding a sparse, spectrally converging ROM generated entirely from the data.

### 4.3. Embedding into the continuous problem and generation of internal data.

We now describe the algorithm to generate the internal solutions, which is simply a generalization of that from Section 2.3. Here we assume we know or can compute solutions to the reference problem

 (27) −Δ(uri)0+q0(uri)0+bi(uri)0 = 0   in    Ω ∂(uri)0∂ν = gr   on ∂Ω

Algorithm

1. Read data (again here synthetically generated)

 Firl=∫∂Ωurigl

and

 dFirldλ

corresponding to some positive for and sources/receivers for the perturbed problem (1).

2. Compute all of the corresponding solutions for for the reference problem (27) and the corresponding reference data sets and .

3. Use (22, 23, 24, 25) to generate the systems of block stiffness and mass matrices and for Galerkin system for the perturbed problem, and similarly generate , for the reference problem.

4. Compute projection of the vector of functionals onto

 VmK=span{{uri}i=1,…,m,r=1,…,K}

and onto the corresponding reference Galerkin space

 V0mK=span{{(uri)0}i=1,…,m,r=1,…,K},

that is, we compute the coefficient vectors and by solving the systems

 M→d=F

and

 M0→d0=F0

respectively.

5. Set , and perform block Lanczos orthogonalization with respect to the inner product using

 {→d,A→d,A2→d,…,Am−1→d}

and let , be the Galerkin stiffness and mass matrices in this new basis for . Note that is block tridiagonal with blocks.

6. Similarly, set and perform Lanczos orthogonalization with respect to the inner product on the basis

 {→d0,A0→d0,A20→d0,…,Am−10→d0}

to obtain the corresponding orthogonalized basis for .

7. Now choose any spectral value and solve the Galerkin system in the new basis , so that for each ,