169
views
0
8 months ago by
I'm trying to use dolfin-adjoint for Bayesian computations, but I came to the point where I'm not able to implement my functional of interest.
The problem I'm trying to solve is the Neumann-to-Dirichlet problem modelling Electric Impedance Tomography.
The system of equations that I obtain with dolfin is $A\left(\sigma\right)u\left(\sigma\right)=b$A(σ)u(σ)=b .
I define the operator   $F_g\left(\sigma\right)$Fg(σ)  mapping   $\sigma\mapsto u$σu , the trace operator $T$T  and the observation operator  $P_g\left(\sigma\right):=TF_g\left(\sigma\right)$Pg(σ):=TFg(σ) .
I assume to have some noise-perturbed data   $\phi=P_g\left(\sigma\right)+\varepsilon$ϕ=Pg(σ)+ε  , where  $\varepsilon\sim GP\left(0,C\right)$εGP(0,C) is a Gaussian process with circular covariance  $C\left(x,y\right)$C(x,y) .
My functional of interest is then the log-likelihood of  $\phi$ϕ  given  $\sigma$σ , i.e. $\log\pi\left(\phi|\sigma\right)$logπ(ϕ|σ). By the definition of the noise this is

$\log\pi\left(\phi|\sigma\right)\propto-\frac{1}{2}\left(\phi-P_g\left(\sigma\right)\right)^T\Sigma^{-1}\left(\phi-P_g\left(\sigma\right)\right)$logπ(ϕ|σ)12 (ϕPg(σ))TΣ1(ϕPg(σ))

where   $\Sigma=C\left(\mathbf{x},\mathbf{x^'}\right)$Σ=C(x,x') is the covariance evaluated at the grid points. Then my functional is   $J\left(\sigma\right):=\log\pi\left(\phi|\sigma\right)$J(σ):=logπ(ϕ|σ)  .
Now I can write out the equations for  $\nabla_{\sigma}J$σJ using the chain rule

$\nabla_{\sigma}J\left(P_g\left(\sigma\right)\right)=-\lambda^T\nabla_{\sigma}A\left(\sigma\right)F_g\left(\sigma\right)$σJ(Pg(σ))=λTσA(σ)Fg(σ)

where  $\lambda$λ  is the solution of the adjoint equation

$A\left(\sigma\right)^T\lambda=T^T\left(\nabla_yJ\left(y\right)|_{y=P_g\left(\sigma\right)}\right)^T$A(σ)Tλ=TT(yJ(y)|y=Pg(σ))T .

This formulation requires only a single solve in order to obtain the gradient.
However I don't know how to implement the functional  $J$J  in the dolfin formalism and I can see only two ways out:
1. I find a way to implement  $J$J  (most elegant solution)
2. I need to extract  $\nabla_{\sigma}A\left(\sigma\right)$σA(σ)  and (possibly) the latest factorization used of  $A\left(\sigma\right)$A(σ). Is this possible at all?
Is anybody able to point me into some useful direction? Thanks.
Community: FEniCS Project