1 
\chapter{Regularization}\label{Chp:ref:regularization} 
2 

3 
The general cost function $J^{total}$ to be minimized has some of the cost 
4 
function $J^f$ measuring the defect of the result from the 
5 
forward model with the data, and the cost function $J^{reg}$ introducing the 
6 
regularization into the problem and makes sure that a unique answer exists. 
7 
The regularization term is a function of, possibly vectorvalued, level set 
8 
function $m$ which represents the physical properties to be represented and is, 
9 
from a mathematical point of view, the unknown of the inversion problem. 
10 
It is the intention that the values of $m$ are between zero and one and that 
11 
actual physical values are created from a mapping before being fed into a 
12 
forward model. In general the cost function $J^{reg}$ is defined as 
13 
\begin{equation}\label{EQU:REG:1} 
14 
J^{reg}(m) = \frac{1}{2} \int_{\Omega} \left( 
15 
\sum_{k} \mu_k \cdot ( \omega^{(0)}_k \cdot m_k^2 + \omega^{(1)}_{ki}m_{k,i}^2 ) 
16 
+ \sum_{l<k} \mu^{(c)}_{lk} \cdot \omega^{(c)}_{lk} \cdot \chi(m_l,m_k) \right) \; dx 
17 
\end{equation} 
18 
where summation over $i$ is performed. The additional tradeoff factors 
19 
$\mu_k$ and $\mu^{(c)}_{lk}$ ($l<k$) are between zero and one and constant across the 
20 
domain. They are potentially modified during the inversion in order to improve the balance between the different terms 
21 
in the cost function. 
22 

23 
$\chi$ is a given symmetric, nonnegative crossgradient function\index{crossgradient 
24 
}. We use 
25 
\begin{equation}\label{EQU:REG:4} 
26 
\chi(a,b) = ( a_{,i} a_{,i}) \cdot ( b_{,j} b_{,j})  ( a_{,i} b_{,i})^2 
27 
\end{equation} 
28 
where summations over $i$ and $j$ are performed. Notice that crossgradient function 
29 
is measuring the angle between the surface normals of contours of level set functions. So 
30 
minimizing the cost function will align the surface normals of the contours. 
31 

32 
The coefficients $\omega^{(0)}_k$, $\omega^{(1)}_{ki}$ and $\omega^{(c)}_{lk}$ define weighting factors which 
33 
may depend on their location within the domain. We assume that for given level set function $k$ the 
34 
weighting factors $\omega^{(0)}_k$, $\omega^{(1)}_{ki}$ are scaled such that 
35 
\begin{equation}\label{ref:EQU:REG:5} 
36 
\int_{\Omega} ( \omega^{(0)}_k + \frac{\omega^{(1)}_{ki}}{L_i^2} ) \; dx = \alpha_k 
37 
\end{equation} 
38 
where $\alpha_k$ defines the scale which is typically set to one. $L_i$ is the width of the domain in $x_i$ direction. 
39 
Similarly we set for $l<k$ we set 
40 
\begin{equation}\label{ref:EQU:REG:6} 
41 
\int_{\Omega} \frac{\omega^{(c)}_{lk}}{L^4} \; dx = \alpha^{(c)}_{lk} 
42 
\end{equation} 
43 
where $\alpha^{(c)}_{lk}$ defines the scale which is typically set to one and 
44 
\begin{equation}\label{ref:EQU:REG:6b} 
45 
\frac{1}{L^2} = \sum_i \frac{1}{L_i^2} \;. 
46 
\end{equation} 
47 

48 
In some cases values for the level set functions are known to be zero at certain regions in the domain. Typically this is the region 
49 
above the surface of the Earths. This expressed using a 
50 
a characteristic function $q$ which varies with its location within the domain. The function $q$ is set to zero except for those 
51 
locations $x$ within the domain where the values of the level set functions is known to be zero. For these locations $x$ 
52 
$q$ takes a positive value. for a single level set function one has 
53 
\begin{equation}\label{ref:EQU:REG:7} 
54 
q(x) = \left\{ 
55 
\begin{array}{rl} 
56 
1 & \mbox{ if } m \mbox{ is set to zero at location } x \\ 
57 
0 & \mbox{ otherwise } 
58 
\end{array} 
59 
\right. 
60 
\end{equation} 
61 
For multivalued level set function the characteristic function is set componentwise: 
62 
\begin{equation}\label{ref:EQU:REG:7b} 
63 
q_k(x) = \left\{ 
64 
\begin{array}{rl} 
65 
1 & \mbox{ if component } m_k \mbox{ is set to zero at location } x \\ 
66 
0 & \mbox{ otherwise } 
67 
\end{array} 
68 
\right. 
69 
\end{equation} 
70 

71 

72 
\section{Usage} 
73 

74 
\LG{Add example} 
75 

76 
\begin{classdesc}{Regularization}{domain 
77 
\optional{, w0=\None} 
78 
\optional{, w1=\None} 
79 
\optional{, wc=\None} 
80 
\optional{, location_of_set_m=Data()} 
81 
\optional{, numLevelSets=1} 
82 
\optional{, useDiagonalHessianApproximation=\False} 
83 
\optional{, tol=1e8} 
84 
\optional{, scale=\None} 
85 
\optional{, scale_c=\None} 
86 
} 
87 

88 

89 
initializes a regularization component of the cost function for inversion. 
90 
\member{domain} defines the domain of the inversion. \member{numLevelSets} 
91 
sets the number of level set functions to be found during the inversion. 
92 
\member{w0}, \member{w1} and \member{wc} define the weighting factors 
93 
$\omega^{(0)}$, 
94 
$\omega^{(1)}$ and 
95 
$\omega^{(c)}$, respectively. A value for \member{w0} or \member{w1} or both must be given. 
96 
If more then one level set function is involved \member{wc} must be given. 
97 
\member{location_of_set_m} sets the characteristic function $q$ 
98 
to define locations where the level set function is set to zero, see equation~(\ref{ref:EQU:REG:7}). 
99 
\member{scale} and 
100 
\member{scale_c} set the scales $\alpha_k$ in equation~(\ref{ref:EQU:REG:5}) and 
101 
$\alpha^{(c)}_{lk}$ in equation~(\ref{ref:EQU:REG:6}), respectively. By default, their values are set to one. 
102 
Notice that weighting factors are rescaled to meet the scaling conditions. \member{tol} sets the 
103 
tolerance for the calculation of the Hessian approximation. \member{useDiagonalHessianApproximation} 
104 
indicates to ignore coupling in the Hessian approximation produced by the 
105 
crossgradient term. This can speedup an individual iteration step in the inversion but typically leads to more 
106 
inversion steps. 
107 
\end{classdesc} 
108 

109 
\section{Gradient Calculation} 
110 

111 

112 
The cost function kernel\index{cost function!kernel} is given as 
113 
\begin{equation}\label{ref:EQU:REG:100} 
114 
K^{reg}(m) = \frac{1}{2} 
115 
\sum_{k} \mu_k \cdot ( \omega^{(0)}_k \cdot m_k^2 + \omega^{(1)}_{ki}m_{k,i}^2 ) 
116 
+ \sum_{l<k} \mu^{(c)}_{lk} \cdot \omega^{(c)}_{lk} \cdot \chi(m_l,m_k) 
117 
\end{equation} 
118 
We need to provide the gradient of the cost function $J^{reg}$ with respect to the level set functions $m$. 
119 
The gradient is represented by two functions $Y$ and $X$ which define the 
120 
derivative of the cost function kernel with respect to $m$ and to the gradient $m_{,i}$, respectively: 
121 
\begin{equation}\label{ref:EQU:REG:101} 
122 
\begin{array}{rcl} 
123 
Y_k & = & \displaystyle{\frac{\partial K^{reg}}{\partial m_k}} \\ 
124 
X_{ki} & = & \displaystyle{\frac{\partial K^{reg}}{\partial m_{k,i}}} 
125 
\end{array} 
126 
\end{equation} 
127 
For the case of a single valued level set function $m$ we get 
128 
\begin{equation}\label{ref:EQU:REG:202} 
129 
Y = \mu \cdot \omega^{(0)} \cdot m 
130 
\end{equation} 
131 
and 
132 
\begin{equation}\label{ref:EQU:REG:203} 
133 
X_{i} = \mu \cdot \omega^{(1)}_{i} \cdot m_{,i} 
134 
\end{equation} 
135 
For a twovalued level set function $(m_0,m_1)$ we have 
136 
\begin{equation}\label{ref:EQU:REG:302} 
137 
Y_k = \mu_k \cdot \omega^{(0)}_k \cdot m_k \mbox{ for } k=0,1 
138 
\end{equation} 
139 
and for $X$ 
140 
\begin{equation}\label{ref:EQU:REG:303} 
141 
\begin{array}{rcl} 
142 
X_{0i} & = & \mu_0 \cdot \omega^{(1)}_{0i} \cdot m_{0,i} + \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot 
143 
\left( (m_{1,j}m_{1,j} ) \cdot m_{0,i}  (m_{1,j}m_{0,j} ) \cdot m_{1,i} \right) \\ 
144 
X_{1i} & = & \mu_1 \cdot \omega^{(1)}_{1i} \cdot m_{1,i} + \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot 
145 
\left( (m_{0,j}m_{0,j} ) \cdot m_{1,i}  (m_{1,j}m_{0,j} ) \cdot m_{0,i} \right) 
146 
\\ 
147 
\end{array} 
148 
\end{equation} 
149 
We also need to provide an approximation of the inverse of the Hessian operator. The operator evaluation is executes as a solution 
150 
of a linear PDE which is solved using \escript \class{LinearPDE} class. In the \escript notation we need to provide 
151 
\begin{equation}\label{ref:EQU:REG:600} 
152 
\begin{array}{rcl} 
153 
A_{kilj} & = & \displaystyle{\frac{\partial X_{ki}}{\partial m_{l,j}}} \\ 
154 
D_{kl} & = & \displaystyle{\frac{\partial Y_{k}}{\partial m_{l}}} 
155 
\end{array} 
156 
\end{equation} 
157 
For the case of a single valued level set function $m$ we get 
158 
\begin{equation}\label{ref:EQU:REG:601} 
159 
\begin{array}{rcl} 
160 
A_{ij} & =& \mu \cdot \omega^{(1)}_i \cdot \delta_{ij} \\ 
161 
D & = & \mu \cdot \omega^{(0)} 
162 
\end{array} 
163 
\end{equation} 
164 
For a twovalued level set function $(m_0,m_1)$ we have 
165 
\begin{equation}\label{ref:EQU:REG:602} 
166 
D_{kl} = \mu_k \cdot \omega^{(0)}_k \cdot \delta_{kl} 
167 
\end{equation} 
168 
and 
169 
\begin{equation}\label{ref:EQU:REG:603} 
170 
\begin{array}{rcl} 
171 
A_{0i0j} & = & \mu_0 \cdot \omega^{(1)}_{0i} \cdot \delta_{ij} + \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot 
172 
\left( (m_{1,j'}m_{1,j'} )\cdot \delta_{ij}  m_{1,i} \cdot m_{1,j} \right) \\ 
173 
A_{0i1j} & = & \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot \left( 2 \cdot m_{0,i} \cdot m_{1,j} 
174 
 m_{1,i} \cdot m_{0,j}  ( m_{1,j'} m_{0,j'} ) \cdot \delta_{ij} 
175 
\right) \\ 
176 
A_{1i0j} & = & \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot \left( 2 \cdot m_{1,i} \cdot m_{0,j} 
177 
 m_{0,i} \cdot m_{1,j}  ( m_{1,j'} m_{0,j'} ) \cdot \delta_{ij} \right) \\ 
178 
A_{1i1j} & = & \mu_1 \cdot \omega^{(1)}_{1i} \cdot \delta_{ij} + \mu^{(c)}_{01} \cdot \omega^{(c)}_{01} \cdot 
179 
\left( (m_{0,j'}m_{0,j'} ) \cdot \delta_{ij}  m_{0,i} \cdot m_{0,j} ) \right) 
180 
\end{array} 
181 
\end{equation} 
182 

183 
