Stochastic algorithms for trace estimation and inverse problems - - PowerPoint PPT Presentation

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

. .

Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions

Uri Ascher Department of Computer Science University of British Columbia May 2014

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 1 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . PDE inverse problems & many data sets

Nonlinear parameter function estimation problems involving partial differential equation (PDE) constraints Many measurements for obtaining credible reconstructions Applications

Electromagnetic data inversion in mining exploration Seismic data inversion in oil exploration Direct current (DC) resistivity Electrical impedance tomography (EIT) Diffuse optical tomography (DOT) Quantitative photoacoustic tomography (QPAT)

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 2 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . Inverse problem

After discretization and for our problems of interest: di = Fi(m) + εi i = 1, 2, . . . s Fi(m) = Piui = PiG(m)qi Calculating “G(m)qi” for each i is costly!

di ∈ Rl is the measurement obtained in the ith experiment Fi is the known forward operator for the ith experiment m ∈ Rlm is the sought-after model εi is the noise incurred in the ith experiment s is the total number of experiments ui ∈ Rlu is the ith field qi ∈ Rlu is the ith source G −1 is a square matrix discretizing the PDE with the BC Pi ia the projection matrix for the ith experiment

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 3 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . Merit function

Assume εi ∼ N(0, ˜ σI) Maximum Likelihood (ML) data misfit: φ(m) = ∥F(m) − D∥2

F

F ∈ Rl×s with Fi(m) for ith column D ∈ Rl×s with di for ith column

Maximum a Posteriori (MAP) merit function: φR,β(m) = φ(m) + βR(m)

R(m) can be implicit (e.g. dynamic regularization)

Seek to reduce merit function (not much) below the noise level.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 4 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . Numerical example: many data sets

[van den Doel, Ascher & Haber, 2013]

Forward problem: ∇ · (σ∇ui) = qi, x ∈ Ω, i = 1, . . . , s, ∂ui ∂ν |∂Ω = 0. Take Ω to be a square. Construct s data sets, choosing different current patterns: qi(x) = δx,p

iL L − δx,p iR R ,

where piL

L and piR R are located on the left and right boundaries resp.

Place each at √s different, equidistant locations, 1 ≤ iL, iR ≤ √s.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 5 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . L2G vs L1G results (3% noise)

Left: true. Center: s = 4, L2G. Right: s = 4, L1G. Left: true. Center: s = 64, L2G. Right: s = 64, L1G.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 6 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . L2G vs L1G results (1% noise)

Decrease noise to 1%, increase √s: now there is enough information for L1G to shine.

2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9

Left: s = 1024, L2G. Right: s = 1024, L1G.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 7 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Many data sets

. . Model reduction and stochastic approximation

Need to solve s PDEs just to evaluate forward operator F or misfit function! Proceed with caution. Consider framework of iterative method (e.g. Gauss-Newton for a Tikhonov regularization): at kth iteration, given current iterate mk, solve for correction µk and set mk+1 = mk + µk. Want to make do with fewer data set combinations at each iteration: sk ≪ s at iteration k. Can we apply stochastic/deterministic model reduction?

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 8 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Motivation Outline

. . Outline

Motivation

Many measurements

Implicit matrix trace estimation Stochastic/deterministic model reduction Adaptive selection of sample size Numerical experiments: EIT/DC resistivity inverse problem Conclusions

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 9 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Trace? what trace?

Let B(m) = F(m) − D ∈ Rl×s Matrix A = BTB is implicit symmetric positive semi-definite (SPSD); can only carry out matvec products A ∗ v with this s × s matrix. Let w be drawn from distribution such that E(wwT) = I. Then φ(m) = ∥B(m)∥2

F = tr(BTB) = E(wTAw)

Approximating expectation ⇔ Approximating the trace

Monte-Carlo approximation tr(A) ≈ 1 n

n

∑

i=1

wT

i Awi

Note we can obtain exact trace using n = s samples with wi scaled ith column of identity; but we want n ≪ s. Link notation with our application: n = sk

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 10 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Which probability distribution?

Can choose from several: Gaussian: normal distribution Hutchinson: Rademacher distribution – for w = (w1, . . . , ws)T, Pr(wj = 1) = Pr(wj = −1) = 1/2 Random subset (RS) of scaled columns of the identity matrix. Selling features: Hutchinson (1990) showed that Rademacher has smaller variance than normal distribution. This made the Hutchinson choice popular in industry, too. RS applies to general matrices, unlike the other two. On the other hand, for SPSD matrices including in our motivating example, RS is somewhat slower than the other two.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 11 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Probabilistic bounds

[Roosta-Khorasani & Ascher, 2013; Avron & Toledo, 2011] following [Achlioptas, 2001], want n smallest for a desired quality, look for probabilistic bounds.

Given pair (ε, δ) Probability distribution ∆ Want lower bound on n such that Pr ( |trn

∆(A) − tr(A)| ≤ ε tr(A)

) ≥ 1 − δ trn

∆(A) = 1

n

n

∑

i=1

wT

i Awi

with wi from probability distribution ∆

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 12 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Theorems – Hutchinson

Given a pair (ε, δ), denote c = c(ε, δ) = ε−2 ln(2/δ). Assume A is an s × s SPSD matrix. Using Hutchinson, the probabilistic bound holds if . .

1

n ≥ 6c(ε, δ) and also if . .

2

n > 2KHc(ε, δ) where Kj

H =

∥aj∥2 − a2

j,j

a2

j,j

= ∑

k̸=j

a2

k,j / a2 j,j,

KH = max

j

Kj

H.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 13 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Theorems – Gauss

Given a pair (ε, δ), denote c = c(ε, δ) = ε−2 ln(2/δ). Assume A is an s × s SPSD matrix. Using Gauss, the probabilistic bound holds if . .

1

n > 8KGc(ε, δ) where Kj

G =

λj tr(A), KG = max

j

Kj

G = ∥A∥

tr(A), and in particular if n ≥ 8c(ε, δ). . .

2 There is also an important necessary bound on how large N must be

for the probabilistic bound to hold. This bound kicks in when rank(A) is small.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 14 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Necessary bound – Gauss

Recall regularized Gamma functions P (α, β) := γ (α, β) /Γ (α) , Q (α, β) := Γ (α, β) /Γ (α) , where γ (α, β) , Γ (α, β) and Γ (α) are lower incomplete, upper incomplete and complete Gamma functions, resp. Also, Γ (α) = Γ (α, β) + γ (α, β). Define Φθ(x) := P (x 2, τ(1 − θ)x 2 ) + Q (x 2, τ(1 + θ)x 2 ) , where τ = ln(1 + θ) − ln(1 − θ) 2θ . Theorem: if the probabilistic bound holds then necessarily Φε(n ∗ rank(A)) ≤ δ.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 15 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Necessary bound – Gauss

Recall regularized Gamma functions P (α, β) := γ (α, β) /Γ (α) , Q (α, β) := Γ (α, β) /Γ (α) , where γ (α, β) , Γ (α, β) and Γ (α) are lower incomplete, upper incomplete and complete Gamma functions, resp. Also, Γ (α) = Γ (α, β) + γ (α, β). Define Φθ(x) := P (x 2, τ(1 − θ)x 2 ) + Q (x 2, τ(1 + θ)x 2 ) , where τ = ln(1 + θ) − ln(1 − θ) 2θ . Theorem: if the probabilistic bound holds then necessarily Φε(n ∗ rank(A)) ≤ δ.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 15 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Theorems – Random Subset

Let A be a real s × s matrix, and denote K(i,j)

U

= s tr(A) |aii − ajj| , KU = max

1≤i,j≤s i̸=j

K(i,j)

U

. Using RS, the probabilistic bound holds with replacement if n > K2

U

2 c(ε, δ) ≡ F, and without replacement if n ≥ s + 1 1 + s−1

F

.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 16 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Trace estimators

. . Example (s = 1000, ε = δ = 0.2)

Consider A = xxT/∥x∥2, where x ∈ Rs, and for some θ > 0, xj = exp(−jθ), 1 ≤ j ≤ s. Then tr(A) = 1, r = 1, KG = 1, KH = ∥x∥2x−2

s

− 1, KU = s ∥x∥2 (x2

1 − x2 s ).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 17 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Implicit matrix trace estimation Outline

. . Outline

Motivation Implicit matrix trace estimation Stochastic/deterministic model reduction

Simultaneous sources and random subset Modified Gauss-Newton for inverse problem Sampling methods Stopping criterion and uncertainty check

Adaptive selection of sample size Numerical experiments: EIT/DC resistivity inverse problem Conclusions

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 18 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Simultaneous sources and random subset

. . Unbiased estimator of misfit function

Recall our misfit function measuring deviation between predicted data F(m) and observed data D: φ(m) = ∥F(m) − D∥2

F.

Approximate this by an unbiased estimator ˆ φ(m, W ) = 1 sk ∥(F(m) − D)W ∥2

F,

where W is an s × sk matrix W = Wk = [ w1, w2, . . . , wsk ] , with wi drawn from one of Gaussian, Hutchinson or RS distributions.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 19 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Simultaneous sources and random subset

. . In the context of our inverse problems

Recall Fi(m) = PiG(m)qi, i = 1, 2, . . . , s. So, φ(m) = E ( ∥

s

∑

i=1

wi(PiG(m)qi − di)∥2 ) . Random subset method: at iteration k, use a subset of n = sk data sets to approximate φ(m). [van den Doel & Ascher, 2012] Possibly better, inexpensive options are available if the receiver locations are the same, i.e., Pi = P, i = 1, . . . , s. Then can write φ(m) = E ( ∥PG(m) (

s

∑

i=1

wiqi ) −

s

∑

i=1

widi∥2 ) . So, only one PDE solve per realization of w – comparable in RS cost. Hutchinson and Gaussian methods for choosing w lead to different methods of simultaneous sources (e.g., [Haber, Chung & Herrmann,

2012]).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 20 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Simultaneous sources and random subset

. . In the context of our inverse problems

Recall Fi(m) = PiG(m)qi, i = 1, 2, . . . , s. So, φ(m) = E ( ∥

s

∑

i=1

wi(PiG(m)qi − di)∥2 ) . Random subset method: at iteration k, use a subset of n = sk data sets to approximate φ(m). [van den Doel & Ascher, 2012] Possibly better, inexpensive options are available if the receiver locations are the same, i.e., Pi = P, i = 1, . . . , s. Then can write φ(m) = E ( ∥PG(m) (

s

∑

i=1

wiqi ) −

s

∑

i=1

widi∥2 ) . So, only one PDE solve per realization of w – comparable in RS cost. Hutchinson and Gaussian methods for choosing w lead to different methods of simultaneous sources (e.g., [Haber, Chung & Herrmann,

2012]).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 20 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Simultaneous sources and random subset

. . In the context of our inverse problems

Recall Fi(m) = PiG(m)qi, i = 1, 2, . . . , s. So, φ(m) = E ( ∥

s

∑

i=1

wi(PiG(m)qi − di)∥2 ) . Random subset method: at iteration k, use a subset of n = sk data sets to approximate φ(m). [van den Doel & Ascher, 2012] Possibly better, inexpensive options are available if the receiver locations are the same, i.e., Pi = P, i = 1, . . . , s. Then can write φ(m) = E ( ∥PG(m) (

s

∑

i=1

wiqi ) −

s

∑

i=1

widi∥2 ) . So, only one PDE solve per realization of w – comparable in RS cost. Hutchinson and Gaussian methods for choosing w lead to different methods of simultaneous sources (e.g., [Haber, Chung & Herrmann,

2012]).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 20 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Gauss-Newton and inverse problem

. . Aside – optimization process for φ(m)

In principle, want to iterate towards minimizing (assume one experiment) φ(m) = ∥F(m) − d∥2, reducing the objective (misfit) function to the noise level, but not beyond. Note also that our highly ill-posed problems require regularization. Gauss-Newton (GN) Set m = mk at the kth iteration. Obtain next iterate by approximating ( JTJ ) µk = −∇φ(m) mk+1 ← mk + µk. But the linear system is singular! Moreover, even for one experiment, F(m) = PG(m)q, the Jacobian J(m) = −PG(m)S is very expensive to calculate and store.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 21 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Gauss-Newton and inverse problem

. . Aside – optimization process for φ(m)

In principle, want to iterate towards minimizing (assume one experiment) φ(m) = ∥F(m) − d∥2, reducing the objective (misfit) function to the noise level, but not beyond. Note also that our highly ill-posed problems require regularization. Gauss-Newton (GN) Set m = mk at the kth iteration. Obtain next iterate by approximating ( JTJ ) µk = −∇φ(m) mk+1 ← mk + µk. But the linear system is singular! Moreover, even for one experiment, F(m) = PG(m)q, the Jacobian J(m) = −PG(m)S is very expensive to calculate and store.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 21 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Gauss-Newton and inverse problem

. . Aside – optimization process for φ(m)

In principle, want to iterate towards minimizing (assume one experiment) φ(m) = ∥F(m) − d∥2, reducing the objective (misfit) function to the noise level, but not beyond. Note also that our highly ill-posed problems require regularization. Gauss-Newton (GN) Set m = mk at the kth iteration. Obtain next iterate by approximating ( JTJ ) µk = −∇φ(m) mk+1 ← mk + µk. But the linear system is singular! Moreover, even for one experiment, F(m) = PG(m)q, the Jacobian J(m) = −PG(m)S is very expensive to calculate and store.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 21 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Gauss-Newton and inverse problem

. . Modified Gauss-Newton

So, use a few conjugate gradient (CG) inner iterations, preconditioned with a Laplacian. Modified (stabilized) GN: JT(Jµk) ≈ −∇φ(mk) γk ≈ argmin φ(mk + γµk) mk+1 ← mk + γkµk. Obtain inexpensive regularization by iterative regularization, dynamical regularization [Hansen, 1999; van den Doel & Ascher, 2007].

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 22 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Gauss-Newton and inverse problem

. . Modified Gauss-Newton for reduced φ

For the multi-experiment case with shared receivers, write the

• bjective as

φ(m) = ∥F(m) − D∥2

F,

F(m) = PG(m)Q, Q = [ q1, q2, . . . , qs ] At kth GN iteration, use s × sk matrix W = Wk = [ w1, w2, . . . , wsk ] to estimate φ by ˆ φ(m, W ) = 1 sk ∥(F(m) − D)W ∥2

F.

The essential GN step now reads ( sk ∑

i=1

ˆ JT

i ˆ

Ji ) µk = −∇m ˆ φ. How should we choose W ?

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 23 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Sampling methods

. . Choosing random weights for the kth iteration

Need approximate evaluations of misfit function for three different purposes:

. .

1 Determine the next iterate correction, µk, e.g. by GN.

. .

2 Decide if next iterate is acceptable.

. .

3 Check for algorithm termination.

For each of these can use Hutchinson, Gauss or random subset.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 24 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Stopping criterion and uncertainty check

. . Stopping criterion

Given tolerance ρ, the stopping criterion is φ(m) ≤ ρ. For instance, assume for all s experiments a Gaussian noise distribution with known (same) standard deviation ˜ σ. Then discrepancy principle yields ρ = η˜ σ2sl (say η = 1.5). However, calculating φ(m) requires s PDE solves! Fortunately, there is a good, unbiased estimator in ˆ φ(m, W ) with W = W e chosen randomly as before.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 25 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Stopping criterion and uncertainty check

. . Uncertainty check

For deciding on what to do with current stabilized GN iteration, note there is no direct comparison between φ(mk+1) and φ(mk). Fortunately, we can perform the inexpensive uncertainty check whether ˆ φ(m, W e) ≤ ρ. Based on this can make a more expensive check of stopping criterion

• r increase/modify data set sample.

Note: in practice we rarely know ρ accurately.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 26 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Stochastic/deterministic model reduction Stopping criterion and uncertainty check

. . Uncertainty check

For deciding on what to do with current stabilized GN iteration, note there is no direct comparison between φ(mk+1) and φ(mk). Fortunately, we can perform the inexpensive uncertainty check whether ˆ φ(m, W e) ≤ ρ. Based on this can make a more expensive check of stopping criterion

• r increase/modify data set sample.

Note: in practice we rarely know ρ accurately.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 26 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Outline

. . Outline

Motivation Implicit matrix trace estimation Stochastic/deterministic model reduction Adaptive selection of sample size Numerical experiments: EIT/DC resistivity inverse problem Conclusions

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 27 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Growing sk

. . Using only sk ≤ s data sets

How should we select sk? May not need all experiments at each iteration k, especially in early

• uter iterations (i.e., small k).

Hence, it makes sense to let sk depend on the outer iteration k and start with a small s0, increasing sample size as needed. An unchecked exponential growth of sk is sometimes desired, but not

• always. It is tied to speed of convergence and may perform poorly

especially for cases with embarrassing redundancy. May want to assess necessity of growth of sample size sk. Such a mechanism is cross validation.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 28 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Growing sk

. . Using only sk ≤ s data sets

How should we select sk? May not need all experiments at each iteration k, especially in early

• uter iterations (i.e., small k).

Hence, it makes sense to let sk depend on the outer iteration k and start with a small s0, increasing sample size as needed. An unchecked exponential growth of sk is sometimes desired, but not

• always. It is tied to speed of convergence and may perform poorly

especially for cases with embarrassing redundancy. May want to assess necessity of growth of sample size sk. Such a mechanism is cross validation.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 28 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Growing sk

. . Using only sk ≤ s data sets

How should we select sk? May not need all experiments at each iteration k, especially in early

• uter iterations (i.e., small k).

Hence, it makes sense to let sk depend on the outer iteration k and start with a small s0, increasing sample size as needed. An unchecked exponential growth of sk is sometimes desired, but not

• always. It is tied to speed of convergence and may perform poorly

especially for cases with embarrassing redundancy. May want to assess necessity of growth of sample size sk. Such a mechanism is cross validation.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 28 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Cross validation: algorithm 2

. . Algorithm

Given: Q, D, ρ, m0, sufficient decrease indicator κ ≤ 1. Initialize: m = m0 , s0 = 1. for k = 0, 1, 2, . . . until termination

. .

1 Draw W f

k ∈ Rs×sk.

. .

2 Fitting: Perform one nonlinear iteration, with W = W f

k .

. .

3 Draw W c

k , W e k , W t k ∈ Rs×sk as needed.

. .

4 Cross validation: if ˆ

φ(mk+1, W c

k ) ≤ κˆ

φ(mk, W c

k ) then

Uncertainty Check: if ˆ φ(mk+1, W e

k ) ≤ ρ, then

Stopping Criterion: terminate if φ(mk+1) ≤ ρ; alternatively, terminate if ˆ φ(mk+1, W t

k ) ≤ ρ;

• therwise set sk+1 = sk.

. .

5 else

Sample size Increase: For example set sk+1 = min(2sk, s).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 29 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Cross validation: algorithm 2

. . Notes

Using this Algorithm yields convergence (theoretical results in

[Friedlander & Schmidt, 2012; Byrd et. al., 2011]).

Estimate cost of algorithm by number of PDE solves (i.e., calculations of G(m)q). Compare cost of Algorithm for the three distributions. At each outer (stabilized GN) iteration k use the same number r of inner PCG iterations for all variants to enable fairer comparison (but see [van den Doel & Ascher, 2012] for adaptive rk).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 30 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Cross validation: algorithm 2

. . Notes

Using this Algorithm yields convergence (theoretical results in

[Friedlander & Schmidt, 2012; Byrd et. al., 2011]).

Estimate cost of algorithm by number of PDE solves (i.e., calculations of G(m)q). Compare cost of Algorithm for the three distributions. At each outer (stabilized GN) iteration k use the same number r of inner PCG iterations for all variants to enable fairer comparison (but see [van den Doel & Ascher, 2012] for adaptive rk).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 30 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Adaptive selection of sample size Cross validation: algorithm 2

. . Notes

Using this Algorithm yields convergence (theoretical results in

[Friedlander & Schmidt, 2012; Byrd et. al., 2011]).

Estimate cost of algorithm by number of PDE solves (i.e., calculations of G(m)q). Compare cost of Algorithm for the three distributions. At each outer (stabilized GN) iteration k use the same number r of inner PCG iterations for all variants to enable fairer comparison (but see [van den Doel & Ascher, 2012] for adaptive rk).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 30 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numeerical experiments Outline

. . Outline

Motivation Implicit matrix trace estimation Stochastic/deterministic model reduction Adaptive selection of sample size Numerical experiments: EIT/DC resistivity inverse problem

Problem setup Incorporating bounds Approximating piecewise constant solution Numerical experiments

Conclusions

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 31 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numeerical experiments EIT/ DC resistivity

. . EIT and DC resistivity problem setup

PDE with multiple sources ∇ · (σ(x)∇ui) = qi, i = 1, . . . , s, ∂ui ∂ν |∂Ω = 0. Conductivity σ(x) is expressed as a pointwise function of m(x). The operator G(m) is the inverse of the above PDE discretized on a staggered grid. Source qi and observation location projections Pi correspond to different selection of sources and receivers. Use any prior we may have for this very difficult problem.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 32 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numeerical experiments Incorporating a priori information

. . Incorporating bounds

Define transfer function ψ(τ) = ψ(τ; θ, α1, α2) = α tanh ( τ αθ ) + α1 + α2 2 , α = α2 − α1 2 , 0 < α1 ≤ α2. Define σ(x) = ψ(m(x)), obtaining PDE ∇ · (ψ(m)∇ui) = qi, i = 1, . . . , s. Note σ > 0. Discretize on staggered grid. Often tight bounds are available, say σmin and σmax, such that σmin ≤ σ(x) ≤ σmax. Then use σ(x) = ψ(m(x); 1, σmin, σmax),

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 33 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numeerical experiments Incorporating a priori information

. . Incorporating bounds

Define transfer function ψ(τ) = ψ(τ; θ, α1, α2) = α tanh ( τ αθ ) + α1 + α2 2 , α = α2 − α1 2 , 0 < α1 ≤ α2. Define σ(x) = ψ(m(x)), obtaining PDE ∇ · (ψ(m)∇ui) = qi, i = 1, . . . , s. Note σ > 0. Discretize on staggered grid. Often tight bounds are available, say σmin and σmax, such that σmin ≤ σ(x) ≤ σmax. Then use σ(x) = ψ(m(x); 1, σmin, σmax),

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 33 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numeerical experiments Incorporating a priori information

. . Piecewise constant solution

Can use the same, but sharper, transfer function to approximate a piecewise constant solution. This is a level set function approach. Assuming σ(x) takes only one of two values σI or σII, let σ(x; h) = ψ(m(x); h, σI, σII), and σ(x) = limh→0 σ(x; h). Note m(x) has bounded first derivatives, whereas σ(x) is discontinuous.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 34 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . 2D setup

Domain Ω is unit square. Sources are of the form qi(x) = δ(x − xi

1) − δ(x − xi 2)

with x1 positive unit point source on west boundary, x2 negative unit point source on east boundary. Vary p boundary wall locations to get s = p2 data sets. Receivers are all grid points on north and south walls. No sources or receivers at corners. Uniform 64 × 64 grid with s = 961. For bounds set σmax = 1.2 max σ(x), σmin = 1.2−1 min σ(x). W e

k from Rademacher (Hutchinson) distribution.

PCG inner iteration limit r = 20; cgtol = 1.e-3.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 35 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . Example 1: s = 961, σI = 1, σII = .1, noise 3%

Method Vanilla Random Subset Hutchinson Gaussian Work 86,490 3,190 2,279 1,618

(a) True model (b) Random subset

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 36 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . Example 2: s = 961, σI = .1, σII = 1, noise 1%

Method Vanilla Random Subset Hutchinson Gaussian Work 128,774 3,921 2,762 2,247

(c) True model (d) Hutchinson (e) Gaussian

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 37 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . Another 2D setup

Same as before except: s = 3, 969. Never evaluate φ(m). Use only Gaussian distribution

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 38 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . Example 3: s = 3, 969, σI = .1, σII = 1, noise 2%

Method Vanilla Gaussian Work 436,590 3,850 Rightmost figure below: reconstruction with s = 49 gives poorer results.

(f) True model (g) Gaussian (h) Gaussian, s=49

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 39 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 2D

. . Example 4: s = 3, 969, σI = .1, σII = 1, σIII = .01, noise 2%

Method Vanilla Gaussian Work 527,877 5,142 Rightmost figure below: reconstruction with s = 49 gives poorer results.

(i) True model (j) Gaussian (k) Gaussian, s=49

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 40 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 3D

. . 3D setup

Domain Ω is unit cube. Sources are of the form qi(x) = δ(x − xi

1) − δ(x − xi 2)

Four boreholes are at the edges of the cube. Source pairs are at

• pposing boreholes.

Receivers are all grid points at top surface. No sources at surface. Uniform 17 × 17 × 17 grid with s = 512. PCG itn limit r = 20 for bounds, r = 5 for level set method. Starting guess is m0 = 0, except for level set method, where it is a smoothed cube in middle of Ω. Note, here there is not much to compress: obtain max sk ≈ .5s.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 41 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 3D

. . Example 5: σI = 1, σII = .1, noise 2%

Method Vanilla Random Subset Hutchinson Gaussian TSVD Work 36,864 6,266 1,166 1,176 1,882

(l) True model (m) Random subset (n) TSVD

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 42 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 3D

. . A level set approach

Sharpening happens at each iteration k across the 0- level set of mk(x).

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 43 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Numerical experiments In 3D

. . Example 6: level set method

Method Vanilla Random Subset Hutchinson Gaussian TSVD Work 45,056 1,498 1,370 978 1,560

(o) Random subset (p) Gaussian (q) Hutchinson

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 44 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Conclusions Conclusions

. . Conclusions plus

We have provided 1 necessary and 6 sufficient bounds for the probabilistic convergence of Monte Carlo methods for trace estimation. Using various variants of our algorithm yields convergence (see theoretical results in [Friedlander & Schmidt, 2012; Byrd et. al., 2011]). The cross-validation algorithm shines when there is enough redundancy in the data. A rigid, deterministic convergence error tolerance is hardly ever known in practice. Relaxing requirements accordingly can lead to significant efficiency improvement. Model reduction methods for EIT/DC resistivity problems with many data sets can make a huge difference in efficiency. The Simultaneous Sources methods are roughly in the same equivalence class, and outperform the (more general) Random Subset method by a factor of 1 to 4, mean 2.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 45 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Conclusions Conclusions

. . Conclusions and future

If receivers are not shared, can directly use only Random Subset. To enable Simultaneous Sources, consider methods for data completion. There are several ways to further improve efficiency of our methods, including a cheaper stopping criterion, variable r and more specific Hessian approximations. All are under current investigation. Mais applications, pro favor.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 46 / 47

. .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. . . .. .

Conclusions References

. . References

• F. Roosta-Khorasani and U. M. Ascher, Improved bounds on sample

size for implicit matrix trace estimators.

• F. Roosta-Khorasani, K. van den Doel and U. M. Ascher, Stochastic

algorithms for inverse problems involving PDEs and many measurements, SIAM J. SISC, 2013.

• K. van den Doel and U. M. Ascher, Adaptive and stochastic

algorithms for EIT and DC resistivity problems with piecewise constant solutions and many measurements, SIAM J. SISC, 2012.

Fred Roosta & Kees van den Doel (UBC) Stochastic algorithms for trace estimation and inverse problems involving many PDE inversions May 2014 47 / 47