Preconditioning of conjugate-gradients in observation space for - - PowerPoint PPT Presentation

Preconditioning of conjugate-gradients in

• bservation space for 4D-VAR
• S. Gratton12
• S. Gurol2

Ph.L. Toint3

• J. Tshimanga1

1ENSEEIHT, Toulouse, France 2CERFACS, Toulouse, France 3FUNDP, Namur, Belgium

Large scale problems and Applications in the Earth Sciences

Gratton, Gurol, Toint,Tshimanga RICAM

Outline

1

Introduction

2

Algorithms for primal space

3

Working in the observation space Krylov subspace methods Results on realistic systems Preconditioning Convergence Properties

4

Conclusions

Gratton, Gurol, Toint,Tshimanga RICAM

Introduction

4D-Var problem: Formulation

→ Large-scale non-linear weighted least-squares problem: min

x∈Rn f(x) = 1

2||x − xb||2

B−1 + 1

2

N

• j=0

||Hj(Mj(x)) − yj||2

R−1

j

where: x ∈ Rn is the control variable The observations yj and the background xb are noisy Mj are model operators Hj are observation operators B is the covariance background error matrix Rj are covariance observation error matrices

Gratton, Gurol, Toint,Tshimanga RICAM

Introduction

4D-Var problem: Formulation

→ Large-scale nonlinear weighted least-squares problem: min

x∈Rn f(x) = 1

2||x − xb||2

B−1 + 1

2

N

• j=0

||Hj(Mj(x)) − yj||2

R−1

j

Typically solved by a standard Gauss-Newton method known as Incremental 4D-Var in data assimilation community

1

Solve linearized subproblem at iteration k (linesearch or trust-region techniques require this kind or steps) min

δx∈Rn J(δx)

= 1 2δx − [xb − x]2

B−1 + 1

2Hδx − d2

R−1

Sequence of quadratic minimization problems

2

Perform update x(k+1)(t0) = x(k)(t0) + δx(k)

Gratton, Gurol, Toint,Tshimanga RICAM

Introduction

Quick note on globalization

We want to converge to a point where the gradient is zero, from any starting point Technically, if the function is regular (gradient Lipschitz continuous, bounded from below), an ”approximate” gradient is all what is needed to have a globally convergent algorithm There are basically two techniques to obtain provable global convergence: linesearch, and trust region In linesearch, a steplength is controlled along a descent direction In trust regions, a model is minimized inside a sphere. In the case of the Taylor model truncated conjugate gradients (CG) do this In the rest of the talk, we focus on CG, and assume it is suitably truncated

Gratton, Gurol, Toint,Tshimanga RICAM

Introduction

4D-Var problem: Solution

From optimality condition (B−1 + HTR−1H)δx = B−1(xb − x) + HTR−1d The aim is to approximately solve sequences of this linear system. Solution algorithms: conjugate gradient like method Exact solution writes: xb − x0 +

• B−1 + HTR−1H

−1 HTR−1 (d − H(xb − x0))

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Preconditioned Lanczos algorithm (F 1/2 is not required!)

For i = 1, 2, ..., l

1

wi = (B−1 + HT R−1H)zi → Construction of the Krylov sequence

2

wi = wi − βivi−1

3

αi = wi, zi →Orthogonalization

4

wi+1 = wi − αivi

5

zi+1 = Fwi+1 →Apply preconditioner

6

βi+1 = zi+1, wi+11/2

7

vi+1 = wi+1/βi+1 →Normalization

8

zi+1 = zi+1/βi+1

9

V = [V, vi+1] →Orthonormal basis for Krylov subspace

10 Ti,i = αi; Ti+1,i = Ti,i+1 = βi+1 →Generate the tridiagonal matrix

Solution

1

yl = T −1

l

β0e1 →Impose the condition rk⊥Kl(A, r0)

2

δxl = FVlyl →Find the approximate solution

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Preconditioned CG algorithm (F 1/2 is not required!)

Initialization r0 = Aδx0 − b, z0 = Fr0, p0 = z0 For i = 0, 1, ...

1

qi = (B−1 + HT R−1H)pi

2

αi =< ri, zi > / < qi, pi > →Compute the step-length

3

δxi+1 = δxi + αipi →Update the iterate

4

ri+1 = ri − αiqi →Update the residual

5

ri+1 = ri+1 − RZT ri+1 →Re-orthogonalization

6

zi+1 = Fri+1 →Update the preconditioned residual

7

βi =< ri+1, zi+1 > / < ri, zi > →Ensure A-conjugate directions

8

R = [R, r/βi] →Re-orthogonalization

9

Z = [Z, z/βi] →Re-orthogonalization

10 pi+1 = zi+1 + βipi

→Update the descent direction

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Inner minimization

Minimizing J[δx0] = 1 2δx0 − [xb − x0]2

B−1 + 1

2Hδx0 − d2

R−1

amounts to solve (B−1 + HTR−1H)δx0 = B−1(xb − x0) + HTR−1d. Exact solution writes xb − x0 +

• B−1 + HTR−1H

−1 HTR−1 d − H(xb − x0)

• ,
• r equivalently (using the S-M-Woodbury formula)

xb − x0 + BHT R + HBHT−1 d − H(xb − x0)

• .

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Dual formulation : PSAS

1 Very popular when few observations compared to model

variables.

2 Relies on

xb − x0 + BHT R + HBHT−1 d − H(xb − x0)

• 3 Iteratively solve
• I + R−1HBH

T

• v = R−1(d − H(xb − x0))

for

• v

4 Set

δx0 = xb − x0 + BHT v

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Experiments

Quadratic cost with primal and dual solvers

20 40 60 80 10

2

10

3

10

4

10

5

size(B) = 1024 size(R) = 256 Iterations Cost function

RPCG PSAS

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Motivation : PSAS and CG-like algorithm

1 CG minimizes the Incremental 4D-Var function during its

• iterations. It minimizes a quadratic approximation of the non

quadratic function : Gauss-Newton in the model space.

2 PSAS does not minimize the Incremental 4D-Var function

during its iterations but works in the observation space. Our goal : put the advantages of both approaches together in a Trust-Region framework, to guarantee convergence: Keeping the variational property, to get the so-called Cauchy decrease even when iterations are truncated. Being computationally efficient whenever the number of

• bservations is significantly smaller than the size of the state

vector. Getting global convergence in the observation space !

Gratton, Gurol, Toint,Tshimanga RICAM

Algorithms for primal space

Dual CG algorithm : assumptions

1 Suppose the CG algorithm is applied to solve the Inc-4D using

a preconditioning matrix F

2 Suppose there exists Gm×m such that

FHT = BHTG For ”exact” preconditioners

• B−1 + HTR−1H

−1HT = BHT I + R−1HBHT−1

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Preconditioned CG on Incremental 4D-Var cost function

Initialization steps Loop: WHILE

1 qi−1 = Api−1 2 αi−1 = rT

i−1zi−1 / qT i−1pi−1

3 vi = vi−1 + αi−1pi−1 4 ri = ri−1 + αi−1qi−1 5 zi = Fri 6 βi = rT

i zi / rT i−1zi−1

7 pi = −zi + βipi−1

Initialization steps Loop: WHILE

1 qi−1 =

(HTR−1H + B−1)pi−1

2 αi−1 = rT

i−1zi−1 / qT i−1pi−1

3 vi = vi−1 + αi−1pi−1 4 ri = ri−1 + αi−1qi−1 5 zi = Fri 6 βi = rT

i zi / rT i−1zi−1

7 pi = −zi + βipi−1 Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

A useful observation

Suppose that

1 BHTG = FHT. 2 v0 = xb − x0.

theres coordinate vectors ri,

• pi,

vi, zi and qi such that ri = HT ri, pi = BHT pi, vi = v0 + BHT vi, zi = BHT zi, qi = HT qi

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Preconditioned CG on Incremental 4D-Var cost function (bis)

Initialization steps given v0; r0 = (HTR−1H + B−1)v0 − b, . . . Loop: WHILE

1 HT

qi−1 = HT(R−1HB−1HT + Im) pi−1

2 αi−1 = rT

i−1zi−1 /

qT

i−1

pi−1

3 BHT

vi = BHT( vi−1 + αi−1 pi−1)

4 HT

ri = HT( ri−1 + αi−1 qi−1)

5 BHT

zi = FHT ri = BHTG ri

6 βi = (rT

i zi / rT i−1zi−1)

7 BHT

pi = BHT(− zi + βi pi−1)

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Restricted PCG (version 1) : expensive

Initialization steps given v0; r0 = (HTR−1H + B−1)v0 − b, . . . Loop: WHILE

1

• qi−1 = (Im + R−1HB−1HT)

pi−1

2 αi−1 =

rT

i−1HBHT

zi−1 / qT

i−1HBHT

pi−1

3

• vi =

vi−1 + αi−1 pi−1

4

ri = ri−1 + αi−1 qi−1

5

zi = FHT ri = G ri

6 βi =

rT

i HBHT

zi / rT

i−1HBHT

zi−1

7

• pi = −

zi + βi pi−1

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

More transformations

1 Consider w and t defined by

wi = HBHT zi and ti = HBHT pi

2 From Restricted PCG (version 1)

ti = −w0 if i = 0, −wi + βiti−1 if i > 0,

3 Use these relations into Restricted PCG (version 1) 4 Transform Restricted PCG (version 1) into Restricted PCG

(version 2)

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Restricted PCG (version 2) : right inner-product!

Initialization steps Loop: WHILE

1

• qi−1 = R−1ti−1 +

pi−1

2 αi−1 = wT

i−1

ri−1 / qT

i−1ti−1

3

• vi =

vi−1 + αi−1 pi−1

4

ri = ri−1 + αi−1 qi−1

5

zi = G ri

6 wi = HBHT

zi

7 βi = wT

i

ri / wT

i−1

ri−1

8

• pi = −

zi + βi pi−1

9 ti = −wi + βiti−1 Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Dual Approach: RPCG and PSAS algorithm

Initialization λ0 = 0, r0 = R−1(d − H(xb − x)),

• z0 = G

r0, p1 = z0, k = 1 Loop on k

1

• qi =

A pi

2

αi =< ri−1, zi−1 >M / < qi, pi >M

3

• vi =

vi−1 + αi pi

4

• ri =

ri−1 − αi qi

5

βi =< ri−1, zi−1 >M / <

• ri−2,

zi−2 >M

6

• zi = G

ri

7

• pi =

zi−1 + βi pi−1

• A = R−1HBHT + Im

G is the preconditioner. M is the inner-product. PSAS Algorithm: M = R cheap matvec RPCG Algorithm: M = HBHT expensive matvec (model integration is required) G should be symmetric w.r.t. to M

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Dual Approach: Precond. RLanczos algorithm and PSAS

For i = 1, 2, ..., l

1

• wi = (I + R−1HBHT )zi

2

• wi =

wi − βi vi−1

3

αi = wi, ziM

4

• wi+1 =

wi − αi vi

5

• zi+1 = G

wi+1

6

βi+1 = zi+1, wi+11/2

M 7

• vi+1 =

wi+1/βi+1

8

• zi+1 =

zi+1/βi+1

9

• V = [

V , vi+1]

10 Ti,i = αi; Ti+1,i = Ti,i+1 = βi+1

Solution

1

yl = T −1

l

β0e1

2

δxl = BHT Vlyl M = HBHT G is the preconditioner When M = R, the iterates are mathematically equivalent to that of PSAS

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Dual approach: Primal equivalent algorithms

Assume that r0 ∈ range(HT ) FHT = BHT G where F is the preconditioner in primal space and G is the preconditioner in dual space. Rlanczos, RPCG, CONGRAD, PCG and Lanczos method in the primal space are mathematically equivalent to each

• ther.

MINRES is not equivalent, it minimizes the Euclidean norm of the gradient of the quadratic.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Krylov subspace methods

(B−1 + HT R−1H)

• A

δx = B−1(xb − x) + HT R−1d

• b

Krylov subspace methods searchs for an approximate solution δxl from a subspace δx0 + Kl(A, r0) where Kl(A, r0) = span

• r0, Ar0, A2r0, ..., Al−1r0
• , r0 = b − Aδx0

Krylov subspace methods impose the Petrov-Galerkin condition rk⊥Ll(A, r0). A is symmetric and positive definite Ll(A, r0) = Kl(A, r0) →Lanczos, Conjugate Gradient (CG) →FOM (unsymmetric case for further reference) →minimizes rkA−1 Ll(A, r0) = AKl(A, r0) →MINRES →minimizes rk2

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Results for ROMS

Observations: SST (Sea Surface Temperature) and SSH(Sea Surface Height)

• bservations from satellites. Sub-surface hydrographic observations from floats.

Number of observations (m): 105 Number of state variables (n): 106 for strong constraint and 107 for weak constraint. Computation: 64 CPUs

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Results for 3D-VAR FGAT NEMOVAR

Observations: Temperature, unbalanced salinity, unbalanced sea surface height Number of observations (m): 2 × 105 Number of state variables (n): 8 × 106 Computation: 8 processors are used

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

More on preconditioning

Use an approximation of the Hessian of the quadratic problem → Limited Memory preconditioning (Fisher (1998), Morales and Nocedal (2000), Tschimanga, Gratton, Sartenaer, Weaver (2008)) The idea is:

1 Formulate the limited memory Quasi-Newton matrix 2 Generate the preconditioner using the information from CG

iterations. For equivalence with the primal method, find G that satisfies FHT = BHT G for a given F For now, assume that H is not changing for each outer loop.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

G as a Quasi-Newton warm-start preconditioner : Gilbert, Lemarechal, Nocedal, Byrd, Zhu

Formulation of F as a Quasi-Newton Limited Memory Preconditioner Fk+1 = (I − τkpkqT

k )Fk(I − τkqkpT k ) + τkpkpT k

pk is the search direction τk = 1/(qT

k pk)

qk = (B−1 + HT R−1H)pk ∆Fk defined by ∆Fk = Fk+1 − Fk, is the solution to the problem: min

∆Fk

• W 1/2∆FkW 1/2
• F

subject to ∆Fk = ∆F T

k ,

Fk+1qk = pk Formulation for G as a Quasi-Newton Limited Memory Preconditioner Gk+1 = (I − τk pk(M qk)T )Gk(I − τk qk pT

k M)

+ τk pk pT

k M

M = HBHT , pk is the search direction,

• qk = (Im + R−1HBHT )

pk and

• τk = 1/(

qT

k HBHT

pk) ∆Gk defined by ∆Gk = Gk+1 − Gk is the solution to the problem: min

∆Gk

• (W M)1/2∆Gk(M−1W )1/2
• F

subject to M∆Gk = ∆GT

k M,

Gk+1 qk = pk Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Computationally efficient RPCG algorithm using Quasi-Newton Preconditioner

Loop: WHILE

1

• qi−1 = R−1ti−1 +

pi−1

2

αi−1 = wT

i−1

ri−1 / qT

i−1ti−1

3

• λi =

λi−1 + αi−1 pi−1

4

ri = ri−1 − αi−1 qi−1

5

li = HBHT ri

6

• zi = G

ri

7

wi = GT li

8

βi = wT

i

ri / wT

i−1

ri−1

9

• pi =

zi + βi pi−1

10 ti = wi + βiti−1 11 mqi−1 = (li−1 − li−2)/αi−1 1 Consider a new vector l is defined as li = HBHT ri 2

• zi = G

ri and wi = HBHT zi 3 HBHT G is symmetric (HF HT = HBHT G) wi = HBHT G ri = GT HBHT ri = GT li 4 Multiply the expression ri = ri−1 − αi qi with HBHT gives HBHT qi = (li − li−1)/αi Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Convergence Properties

If FA has eigenvalues µ1 ≤ µ2 ≤ ... ≤ µn, PCG algorithm satisfies the inequality: xk+1 − x∗A ≤ 2( √µn − √µ1 √µn + √µ1 )k x∗A If G A has eigenvalues ν1 ≤ ν2 ≤ ... ≤ νm, RPCG satisfies the inequality: xk+1 − x∗A ≤ 2( √νm − √ν1 √νm + √ν1 )k x∗A xk+1 − x∗A ≤ 2( √νm − √ν1 √νm + √ν1 )k x∗A ≤ 2( √µn − √µ1 √µn + √µ1 )k x∗A Same iterates, but tighter bound on convergence rate with the dual approach Improvement of tightness can be arbitrarily large on purposely chosen problems

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

When H changes!

When H changes in nonlinear iterations, FHT = BHT G is not

• satisfied. The preconditioner is not symmetric anymore wrt HBHT

and perturbed CG is in trouble. Expression for G Gk+1 = (I − τk pk(M qk)T )Gk(I − τk qk pT

k M) +

τk pk pT

k M

M = HBHT , τk = 1/( qT

k HBHT

pk), qk = (Im + R−1HBHT ) pk

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Solutions (1/2)

Straigthforward approach: re-generate G by using the recent pk and HBHT : costs one matvec per preconditioning pair Accept to handle non symmetry : use FOM algorithm

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Solutions (2/2)

Use FOM with Quasi-Newton preconditioner G where approximated M is used. Approximation with the Davidon Fletcher Powell (DFP) formula.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Towards further reduction of the cost

We have shown that RPCG allows memory and computational cost reduction whenever the number of observation is smaller than the size of the control vector Similar results are possible with other Krylov methods (GMRES, FOM, ...) The question now is: can we reduce cost further ? Possible answer: inexact (cheap) matrix-vector products (truncated B−1, R−1, simplified models, ...) (Simoncini and Szyld, van den Eshop and Sleipen, Giraud, Gratton and Langou, ...) → But, there is a need of a stable modification of RPCG.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

The Arnoldi process

Define (in the full space) A = In + BHT R−1H and set K = BHT , L = R−1H the successive nested Krylov subspaces generated by the sequence b, (γIn + KT L)b, (γIn + KT L)2b, (γIn + KT L)3b, . . . (1)

• r, equivalently, by

b, (KT L)b, (KT L)2b, (KT L)3b, . . . (2) The Arnoldi process generates an orthonormal basis of each of the these subspaces, i.e. a set of vectors {vi}k+1

i=1 with v1 = b/b such that, after k

steps, KT L Vk = Vk+1Hk, (3) where Vk ≡ [v1, . . . , vk] and Hk is a (k + 1) × k upper-Hessenberg matrix.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Related methods: GMRES, MINRES, FOM, CG

Depending on how the matrix Hk is exploited to solve the problem we have The GMRES algorithm (≡ MINRES for KT = L) yk = arg min

y

Hky − β1e1, sk = Vkyk The FOM algorithm (≡ CG for KT = L) H

k y = β1e1,

sk = Vkyk here, H

k is the leading k × k submatrix of Hk.

GMRES (FOM) use long recurrences while MINRES (CG) use short ones. Let rk = (I + KT K)Vkyk − b and fk = 1 2yT

k V T k (γI + KT K)Vkyk − bT Vkyk

→ GMRES and MINRES monotionically minimize rk while FOM and CG monotically minimize fk along the iterations.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Range-space GMRES and FOM (RSGMR and RSFOM)

As CG may be rewritten in the observation space to yield RPCG, algorithms GMRES, MINRES and FOM may be rewritten to yields similar variants. Why a range-space GMRES and FOM (RSGMR and RSFOM)? The FOM setting provides better accuracy and is much better suited to use inexact matrix-vector products. The cost of storing an orthonormal basis of the successive Krylov spaces is much lower for range-space methods than for full-space ones.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Exact and inexact products: FOM vs CG

Is CG a reasonable framework for inexact products ?

Comparing rk/(As∗) for FOM, CG with reortho and CG for exact (left) and inexact (right) products (τ = 10−9, κ ≈ 106) Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Stability and convergence with inexact product

We want to bound rk in the context of Arnoldi process under inexact matrix-vector products. Some reasons to consider this question The inexact nature of computer arithmetic implies that such such errors are inevitable To allow matrix-vector products in an inexact but cheaper form Note that the analysis is for GMRES but that in the context of FOM similar conclusions will hold. standard CG and MINRES are no longer equivalent to FOM and GMRES in the context of unsymmetric perturbations.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Two error models

Assume that each iteration i product by K, K or L is inexact, that is Li = L + EL,i, KT

i = KT + EKT ,i,

and Ki = K + EK,i for some errors EL,i, EKT ,i, and EK,i. Consider two error models to describe the inaccuracy in the matrix-vector products. Backward: EK,i+1 ≤ τK,i+1K, EKT ,i+1 ≤ τKT ,i+1K, EL,i+1 ≤ τL,i+1L, EKT ,∗ ≤ τ∗K Forward: EK,i+1 un ≤ τK,i+1Kun, EKT ,i+1 um ≤ τKT ,i+1Kum EL,i+1 un ≤ τL,i+1Lun EKT ,∗ um ≤ τ∗Kum

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Results for the backward error model

Define qk = Hkyk − βe1, G = max[K, L], ωk = max

i,...,k ˆ

vi κ(K) = condition number of K (... after some analysis ...) Theorem Assume the backward-error model. Then rk ≤

• 2(k + 1) qk + Kωk
• τ∗γ

√ kyk + 4 G2 k

i=1 |[yk]i| τi

• ≤
• 2(k + 1)
• qk + τmaxκ(K) (γ + 4 G2)yk
• .

where τmax = max[τ1, . . . , τk].

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Results for the forward error model

Theorem Assume the forward-error model. Then rk ≤

• 2(k + 1) qk +

√ 2

• τ∗γ

√ kyk + 4 G K k

i=1 |[yk]i| τi

• ≤
• 2(k + 1)
• qk + τmax (γ + 4 G K)yk
• .

Note in both sets of bounds: The first of these bounds allows for variable accuray requirements special role of τ∗.

Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Error models (1)

Is the error model important? (ǫ = 10−5, κ ≈ 102)

Backward error model Forward error model (normalized rk, normalized qk, accuracy threshold τ ) Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Error models (2)

Yes, it can definitely make the difference (ǫ = 10−5, κ ≈ 109)

Backward error model Forward error model (normalized rk, normalized qk, accuracy threshold τ ) Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Fixed vs variable accuracy threshold (1)

Can we use variable accuracy thresholds efficiently? (ǫ = 10−5, κ ≈ 102)

Fixed τ τ ≈ 1/qk (normalized rk, normalized qk, accuracy threshold τ ) Gratton, Gurol, Toint,Tshimanga RICAM

Working in the observation space

Fixed vs variable accuracy threshold (2)

Maybe..., not obvious. (ǫ = 10−5, κ ≈ 102)

Fixed τ τ ≈ 1/qk (normalized rk, normalized qk, accuracy threshold τ ) Gratton, Gurol, Toint,Tshimanga RICAM

Conclusions

Conclusions

We have dual space methods RPCG and RLanczos that generate the same iterates as PCG and Lanczos in primal space Further gains may be obtained from inexact products B1/2 operator is not required with the proposed primal solvers for B preconditioning RPCG and RLanczos were implemented in realistics systems: NEMOVAR thanks to Anthony Weaver and Andrea Piacentini, ROMS thanks to Andy Moore. Preconditioning is possible: find G such that FHT = BHT G

Gratton, Gurol, Toint,Tshimanga RICAM

Conclusions

Thank you for your attention !

Gratton, Gurol, Toint,Tshimanga RICAM