Using Bayes theorem to infer the material parameters of human - - PowerPoint PPT Presentation

▶
using bayes theorem to infer the material parameters of
SMART_READER_LITE
LIVE PREVIEW

Using Bayes theorem to infer the material parameters of human - - PowerPoint PPT Presentation

Oct 20, 2023 119 likes •544 views

slide-1
SLIDE 1

Using Bayes’ theorem to infer the material parameters of human tissue

Jack S. Hale Collaborators: Patrick E. Farrell, Stéphane P. A. Bordas

  • 1

University of Ghent, Bayesian Afternoon Seminar

slide-2
SLIDE 2

Overview

  • Motivation.
  • Bayesian approach to inversion.
  • Relate classical optimisation techniques to the Bayesian

inversion approach.

  • Using a domain specific language for variational forms to solve

the equations.

  • Low-rank updates to deal with high-dimensional posterior

covariance.

  • Example problem: sparse surface observations of a solid block.

2

slide-3
SLIDE 3

Motivation.

3

slide-4
SLIDE 4

4

Source: Phillips

slide-5
SLIDE 5

5

slide-6
SLIDE 6

6

slide-7
SLIDE 7

7

Q: What can we infer about the parameters inside the domain, just from displacement observations on the

  • utside?

Q: Which parameters am I most uncertain about?

slide-8
SLIDE 8

Framework

8

slide-9
SLIDE 9

9

πposterior(x | y) ∝ πlikelihood(y | x)πprior(x)

G x

Parameter Map

+

E ∼ N(0, Γnoise)

y yobs

X ∼ N(¯ x, Γprior)

πposterior(x|y) ∝ exp

  • −1

2||y − G(u)||2

Γ−1

noise

− 1 2||x − ¯ x||Γ−1

prior

  • Inference
slide-10
SLIDE 10

10

G : X → Y

Inverse Problems: A Bayesian perspective. Stuart, Acta Numerica (2010). Contribution: Bayesian inverse problems in an infinite- dimensional setting. When is a Bayesian inverse problem well-posed?

slide-11
SLIDE 11

11

The displacements y for a given material parameter x are defined by a the minimum point

  • f the following Lagrangian:

L(y, x) =

  • Ω ψ(y, x) dx
  • Γ t · y ds

where the energy density functional ψ is defined through the following equations: ψ(u, x) = x 2(Ic d) x ln(J) + λ 2 ln(J)2, F = ∂φ ∂X = I + y, C = FTF, IC = tr(C), J = detF.

B0 B φ

slide-12
SLIDE 12

12

Even once discretised (Finite Element Method)

Colin27 brain atlas

20% extension test, 16 Core Xeon, 1.12 million cells, ~29 secs.

Gh : Rn → Rm n = 1, 112, 000

slide-13
SLIDE 13

Problems

  • Evaluating parameter-to-observable map is very expensive.
  • Discretised parameter space can be very large.
  • Outcome: Exploring posterior with ‘traditional sampling’ is not

going to work.

13

slide-14
SLIDE 14

Solutions

  • 1. Connect Bayesian approach to ideas from

classical optimisation. Using derivatives of posterior in parameter-space (Girolami).

  • 2. Exploiting low-rank structure of prior to posterior

covariance updates (Flath 2012, Spantini 2015).

14

slide-15
SLIDE 15

πposterior(x | y) ∝ πlikelihood(y | x)πprior(x)

xMAP = arg max

x∈Rn πposterior(x | y)

xMAP xCM cov(x | y)

15

slide-16
SLIDE 16

Linearise

16

y = Ax

− ln πposterior(x|y) = 1 2||y − Ax||2

Γ−1

noise

+ 1 2||x − x0||Γ−1

prior

slide-17
SLIDE 17

Take the derivative

17

g(xMAP) := x

  • 1

2y Ax2

Γ1

noise

+ 1 2x x02

Γ1

prior

  • x=xmap

= ATΓ1

noise(y Axmap) + Γ1 prior(xmap x0)

= 0

xMAP =

  • Γ−1

prior − ATΓ−1 noiseA

−1 (ATΓnoisey + Γpriorx0)

slide-18
SLIDE 18

18

MAP estimate. Bound-constrained Quasi-Newton BLMVM with More-Thuente line search and ‘correct’ Riesz map.

slide-19
SLIDE 19

and the second derivative…

19

H := ∇xg = Γ−1

prior − ATΓ−1 noiseA

xMAP =

  • Γ−1

prior − ATΓ−1 noiseA

−1 (ATΓnoisey + Γpriorx0)

πposterior ∼ N(xMAP, H−1)

(After a fair bit of manipulation…)

slide-20
SLIDE 20

MAP estimate

20

xMAP xCM cov(x | y)

slide-21
SLIDE 21

21

xMAP xCM cov(x | y)

H−1(xMAP)

slide-22
SLIDE 22

22

πposterior ∼ N(xMAP, H−1)

πapprox

posterior ∼ N(xMAP, H−1(xMAP))

xMAP xCM cov(x | y)

H−1(xMAP)

slide-23
SLIDE 23

Tools

  • The FEniCS Project is

collection of free software for the automated, efficient solution of differential equations using the finite element method.

  • dolfin-adjoint automatically

derives the discrete adjoint, tangent linear and higher-

  • rder adjoint models from a

high-level description of the forward model.

23

http://fenicsproject.org http://www.dolfin-adjoint.org

Wells, Logg, Rognes, Kirby and many, many others… Farrell, Funke, Ham and Rognes. 2015 Wilkinson Prize for Numerical Software.

slide-24
SLIDE 24

24

The displacements y for a given material parameter x are defined by a the minimum point

  • f the following Lagrangian:

L(y, x) =

  • Ω ψ(y, x) dx
  • Γ t · y ds

where the energy density functional ψ is defined through the following equations: ψ(u, x) = x 2(Ic d) x ln(J) + λ 2 ln(J)2, F = ∂φ ∂X = I + y, C = FTF, IC = tr(C), J = detF.

B0 B φ

slide-25
SLIDE 25

from dolfin import * mesh = UnitSquareMesh(32, 32) U = VectorFunctionSpace(mesh, "CG", 1) V = FunctionSpace(mesh, "CG", 1) # solution u = Function(U) # test functions v = TestFunction(U) # incremental solution du = TrialFunction(U) x = interpolate(Constant(1.0), V) lmbda = interpolate(Constant(100.0), V) dims = mesh.type().dim() I = Identity(dims) F = I + grad(u) C = F.T*F J = det(F) Ic = tr(C) phi = (x/2.0)*(Ic - dims) - x*ln(J) + (lmbda/ 2.0)*(ln(J))**2 Pi = phi*dx # gateux derivative with respect to u in direction v F = derivative(Pi, u, v) # and with respect to u in direction du J = derivative(F, u, du) u_h = Function(U) F_h = replace(F, {u: u_h}) J_h = replace(J, {u: u_h}) solve(F_h == 0, u_h, bcs, J=J_h)

The displacements y for a given material parameter x are defined by a the minimum point

  • f the following Lagrangian:

L(y, x) =

  • Ω ψ(y, x) dx
  • Γ t · y ds

where the energy density functional ψ is defined through the following equations: ψ(u, x) = x 2(Ic d) x ln(J) + λ 2 ln(J)2, F = ∂φ ∂X = I + y, C = FTF, IC = tr(C), J = detF.

slide-26
SLIDE 26

26

slide-27
SLIDE 27

27

slide-28
SLIDE 28

28

Q: What can we infer about the parameters inside the domain, just from displacement observations on the

  • utside?

Q: Which parameters am I most uncertain about?

slide-29
SLIDE 29

29

xmap

slide-30
SLIDE 30

Optimal low-rank updates

30

Γprior Γposterior y

Critical idea: Observations are only informative on a low-dimensional subspace of the parameter space, relative to the prior. Spantini et al. (arXiV)

d2

F(A, B) =

  • i

ln2(σi) Awi = Bσiwi

Förstner metric

slide-31
SLIDE 31

31

Matrix-free Krylov-Schur (SLEPc).

slide-32
SLIDE 32

32

slide-33
SLIDE 33

33

Matches trends from Flath et al. p424 for linear parameter to observable maps.

Γprior

Γpost

slide-34
SLIDE 34

Trailing Eigenvector

34

Direction in parameter space least constrained by the

  • bservations
slide-35
SLIDE 35

35

slide-36
SLIDE 36

36

slide-37
SLIDE 37

Leading Eigenvectors

37

Direction in parameter space most constrained by the

  • bservations
slide-38
SLIDE 38

38

slide-39
SLIDE 39

39

slide-40
SLIDE 40

40

slide-41
SLIDE 41

41

Full Hessian. 4000+ actions. Low-rank update. 292 actions. Huge savings in computational cost. Scales with model dimension because observations stay the same.

slide-42
SLIDE 42

Summary

  • We are developing methods to assess uncertainty in the

recovered parameters in soft tissue.

  • This is done within the framework of Bayesian inversion.
  • FEniCS and dolfin-adjoint makes assembling the equations

relatively easy, solving them is tougher!

  • Next steps: exploring non-Gaussian nature of posterior

using Hamiltonian MCMC. Requires derivatives.

  • Paper soon on arXiv: infinite-dimensional setting, full

derivations of equations, 3D problems.

42