Download Optimization - Common Project - Fall 2008 | MATH 5750 and more Study Guides, Projects, Research Optimization Techniques in Engineering in PDF only on Docsity!
COMMON PROJECT FOR MATH 5750/6880 OPTIMIZATION
FALL 2008
Changes
- In §3, item 2, the regularization term used the wrong matrix.
- Added section on how to compute the derivative of the forward map.
- Added section on the Newton-CG method
- Added section on other regularization terms.
- New files are: cg.m, tvfun.m, tvfun_hess.m, tvfun_hessmult.m, test_tvfun.m, tvfun_init.m, fobj_GNHessian.m.
- Modified files: fobj.m
- Fixed typo in §4. In the description of the adjoint state method, the com- putation of the residual (b) should be: r(σ) = Cu − V.
- Fixed typos in §5. In the description of the Gauss-Newton Hessian ap- proximation computation (3) should be AT^ (σ)ψ = −QCT^ Cδu and (4): w = diag (|T 1 | , |T 2 | , · · · , |Tnt|)(∇ψ · ∇u). 1. Introduction In this project you will use electrical tomography problem (EIT) to image the (electrical) conductivity σ(x) for x in some domain Ω (say the human body), based only on measurements at the boundary ∂Ω (i.e. the skin). The equation describing electric conduction is
(1)
∇ · [σ∇u] = 0,for x ∈ Ω n · ∇u = I,for x ∈ ∂Ω.
This problem only makes sense if the compatibility condition
(2)
∂Ω
I(x)dS(x) = 0
holds, and the solution u is defined up to an additive constant (or grounding po- tential). We chose this constant so that,
(3)
∂Ω
u(x)dS(x) = 0.
Data is acquired in Nexp experiments, where in each experiment a different cur- rent pattern is applied to the boundary ∂Ω (the skin of the patient). Thus the data consists of several currents I(1)^... I(Nexp)^ and the resulting voltages V (1)^... V (Nexp) recorded at the boundary only (since we only have access to the skin of the patient). The EIT problem can then be formulated as a non-linear least squares problem:
(4) min σ
N ∑exp
k=
∂Ω
∣u(x; σ, I(k)) − V (k)(x)
2 dS(x) + Jreg (σ)
here u(x; σ, I(k)) solves (1) with conductivity σ and boundary current I(k), and Jreg is a regularization term that makes the problem well-posed (we discuss this later). 1
Figure 1. Triangulation of the unit disk. .
The first term in the functional measures how close the data for the conductivity σ is to the measured data. Some good references for EIT are [3] and [2] (for the adjoint state method for EIT).
- Finite element discretization Matlab code implementing a finite element discretization of this problem with piecewise linear (P1) triangular elements is provided in http://www.math.utah. edu/math5750_f08/eit.tgz. The domain we take is the unit disk (think of it as a slice of the human body) and we first need to generate a mesh:
% Generate t h e g r i d grid = g e n g r i d ( 3 2 , 0. 1 ) ; The grid is generated by the distmesh package (see [4] and http://www-math.mit. edu/~persson/mesh/) which is included in eit.tgz (addpath distmesh makes sure Matlab knows how to find the grid generator). The output of gengrid is a struct describing the grid with members
- p: a np × 2 matrix where p(i,1) and p(i,2) are the x and y coordinates of the i−th vertex in the grid. Here np is the number of vertices (a.k.a. nodes) in the grid.
- t: a nt × 3 matrix where t(i,:) contains the indices (in p) of the vertices of the i−th triangle in the discretization.
- be: a nbe × 2 matrix where be(i,:) contains the indices (in p) of the endpoints of the i−th boundary edge. For this particular grid see Figure 1. Then we take the conductivity function phantom1, which is supposed to model lungs and heart conductivity and evaluate it at the centers of the triangles of the discretization. c=t r i a n g l e c e n t e r s ( grid ) ; t r u e c o n d = phantom1 ( c ( : , 1 ) , c ( : , 2 ) ) ; The vector truecond (size: nt × 1) is the approximation to the conductivity that is needed for the finite element method. The system for the finite element is setup and solved in the following lines:
% a s s e m b l e sys tem A = s t i f f m a t ( truecond , grid ) ;
- For the regularization term, take for the moment L^2 regularization: Jreg (σ) = alphasigma’(ta.*sigma) where alpha is a regularization parameter that you have to adjust (trial and error!) to get acceptable images, and ta is the vector of areas of each triangle in the discretization that can be obtained with: ta = triangle_areas(grid). (note: multiplying element by element a vector by ta is the same as multiplying by the diagonal matrix with diagonal ta)
- note: It should be possible to avoid using limited memory BFGS, since the number of parameters is relatively small.
- Compute the Fr´echet derivative of the forward map (conductivity to measure- ments at the boundary map) and use it to compute the action of the Gauss- Newton approximation of the Hessian on some vector. We do not need to store the full Hessian to solve the Newton systems with conjugate gradient.
- Use the Newton-CG approach to solve the EIT problem.
- Try other regularization methods: H^1 and TV regularization.
- Code and instructions for the adjoint state method Please see the files fobj.m and test_fobj.m. The first file defines the objective function which is a discrete version of the objective function in (4):
J(σ) = Jmisf it(σ) + αJreg (σ)
∑^ nbe
i=
(r(i)(σ))T^ Qbdry r(i)(σ) + ασT^ Dσ,
where:
- σ is a vector of length nt, giving the value of the conductivity at the center of each triangle.
- r(i)(σ) = Cu(i)(σ) − V (i), i = 1,... , nbe is the residual for the i−th exper- iment.
- The observation matrix C restricts a vector defined on all the nodes to the nodes at the boundary (essentially taking measurements on the skin of the patient, file: obsop.m).
- The matrix D is a nt × nt diagonal matrix with D = diag (|T 1 |, |T 2 |,... , |Tnt|). Here |Ti| is the area of the i−th triangle in the triangulation. The diagonal ta of this matrix can be obtained by ta=triangle_area(grid);
- The matrix Q is a symmetric positive semi-definite matrix of size np × np with entries (file: bmassmat.m):
Qi,j =
∂Ω
φi(x)φj (x)dS(x).
Since the entries Qi,j are boundary integrals, the only nonzero entries are when i and j are boundary nodes. The restriction of Q to the boundary is the symmetric positive definite matrix Qbdry = CQCT^. If we are given a vector v with values of some function v(x) at the boundary nodes, then ∫
∂Ω
|v(x)|^2 dS(x) ≈ vT^ Qbdry v.
Follow these steps to compute the gradient with the adjoint state method. We assume for simplicity a single experiment where current I is applied to the boundary
∂Ω (skin of patient) and V is the voltage measured on ∂Ω. The matrix A(σ) is the stiffness matrix for σ where the last row has been adjusted to make it invertible (enforcing grounding condition), and for simplicity we compute the gradient of the misfit part of (5), that is ∇Jmisf it(σ). (a) Solve the Forward Problem A(σ)u = I
(b) Compute residual r(σ) = Cu − V (c) Solve the Adjoint Problem AT^ (σ)ψ = −QCT^ r(σ)
(d) Compute the gradient (using “discrete gradients” [Bx,By]=gradmat(grid);).
∇J(σ) = diag (|T 1 |,... , |Tnt|)(∇ψ · ∇u), where |Ti| is the area of the i−th triangle. The reason why we multiply by the triangle areas is that in this way:
∇J(σ)T^ δσ =
Ω
∇J(σ)(x)δσ(x)dx.
where the functions in the integrand are piecewise constant interpolations on the triangulation of the vectors δσ and ∇J(σ).
- Computing the Gauss-Newton approximation to the Hessian For simplicity we consider one single experiment where current I is applied to the boundary ∂Ω and V is measured voltage. The multiple experiment case is left to you. Let F : Rnt^ → Rnbe^ be the forward map, that is the map taking a conductivity σ (vector of size nt) and giving the voltage resulting from the experiment setting the current at the boundary to I. With the forward map the misfit part of the objective function in (5) can be written as: Jmisf it(σ) = F T^ (σ)Qbdry F (σ) The Gauss-Newton approximation to the Hessian of Jmisf it is: ∇^2 Jmisf it(σ) ≈ DF [σ]T^ Qbdry DF [σ]. All we need for the CG Newton method is to know how to compute the action of the G-N approximation on some vector v, that is w = DF [σ]T^ Qbdry DF [σ]v. Here are the steps to follow, using the same notation as in §4. (1) Solve the Forward Problem (already done in §4) A(σ)u = I (2) Solve the Linearized Problem (computes δu = DF [σ]v) A(σ)δu = −A(v)u (3) Solve the Adjoint Problem AT^ (σ)ψ = −QCT^ Cδu. (4) Compute w = DF [σ]T^ Qbdry δu w = diag (|T 1 |, |T 2 |,... , |Tnt|)(∇ψ · ∇u).
Unfortunately the TV function is not smooth because |x| is not smooth. Since we need derivatives, we replace the TV functional by a smooth approximation:
J reg(T V )(σ) =
Ω
‖∇σ(x)‖^22 + β^2 dx,
where β > 0 is yet another parameter to adjust. The files related to the (approximate) TV function are:
- tvfun.m computes smoothed TV functional, its gradient and additional information needed to computed the Hessian.
- tvfun_hess.m computes Hessian of the TV functional
- tvfun_hessmult.m applies the Hessian of TV functional to a vector
- test_tvfun.m tests gradient and Hessian computation on a randomly gen- erated unstructured grid.
- tvfun_init.m initializes data structure for TV functional
7.3. H^1 regularization (optional). This regularization favors smooth solutions:
J reg(H1) (σ) =
Ω
‖∇σ(x)‖^22 dx.
For the discrete version of this function, we use the same trick as for TV regular- ization. We define a “gradient” by using gradmat.m to build a “dual” grid out of the centers of the triangles:
c = t r i a n g l e c e n t e r s ( grid ) ; nc=s i z e ( c , 1 ) ; t = delaunay ( c ( : , 1 ) , c ( : , 2 ) ) ; ndt=s i z e ( t , 1 ) ; d g r i d. p=c ; d g r i d. t=t ; [ Gx , Gy]= gradmat ( d g r i d ) ;
Then for some vector s the H^1 functional is: s’Gx’DGxs + s’Gy’DGys, where D = spdiags(triangle_area(dgrid),0,ndt,ndt);. A simpler expression for the H^1 functional is: s’Ls where L = stiffmat(ones(ndt,1),dgrid) is the Laplacian on the “dual” grid.
References
[1] J. Alberty, C. Carstensen, and S. A. Funken. Remarks around 50 lines of Matlab: short finite element implementation. Numer. Algorithms, 20(2-3):117–137, 1999. ISSN 1017-1398. doi:10.1023/A:1019155918070. [2] L. Borcea. Electrical impedance tomography. Inverse Problems, 18:R99–R136,
- doi:10.1088/0266-5611/18/6/201. Topical Review. [3] M. Cheney, D. Isaacson, and J. C. Newell. Electrical impedance to- mography. SIAM Rev., 41(1):85–101 (electronic), 1999. ISSN 1095-7200. doi:10.1137/S0036144598333613. [4] P.-O. Persson and G. Strang. A simple mesh generator in Mat- lab. SIAM Rev., 46(2):329–345 (electronic), 2004. ISSN 0036-1445. doi:10.1137/S0036144503429121.
