1 |
\section{The line search optimisation algorithm |
2 |
\label{sectionoptim}} |
3 |
|
4 |
\subsection{General features} |
5 |
|
6 |
The line search algorithm is based on a quasi-Newton |
7 |
variable storage method which was implemented by |
8 |
\cite{gil-mar:89}. |
9 |
|
10 |
TO BE CONTINUED... |
11 |
|
12 |
\subsection{The online vs. offline version} |
13 |
|
14 |
\begin{itemize} |
15 |
% |
16 |
\item {\bf Online version} \\ |
17 |
Every call to {\it simul} refers to an execution of the |
18 |
forward and adjoint model. |
19 |
Several iterations of optimization may thus be performed within |
20 |
a single run of the main program (lsopt\_top). |
21 |
The following cases may occur: |
22 |
% |
23 |
\begin{itemize} |
24 |
\item |
25 |
cold start only (no optimization) |
26 |
\item |
27 |
cold start, followed by one or several iterations of optimization |
28 |
\item |
29 |
warm start from previous cold start with one or several iterations |
30 |
\item |
31 |
warm start from previous warm start with one or several iterations |
32 |
\end{itemize} |
33 |
% |
34 |
\item {\bf Offline version} \\ |
35 |
Every call to simul refers to a read procedure which |
36 |
reads the result of a forward and adjoint run |
37 |
Therefore, only one call to simul is allowed, |
38 |
{\tt itmax = 0}, for cold start |
39 |
{\tt itmax = 1}, for warm start |
40 |
Also, at the end, {\bf x(i+1)} needs to be computed and saved |
41 |
to be available for the offline model and adjoint run |
42 |
\end{itemize} |
43 |
|
44 |
In order to achieve minimum difference between the online and offline code |
45 |
{\bf xdiff(i+1)} is stored to file at the end of an (offline) iteration, |
46 |
but recomputed identically at the beginning of the next iteration. |
47 |
|
48 |
\subsection{Number of iterations vs. number of simulations} |
49 |
|
50 |
{\tt - itmax:} controls the max. number of iterations \\ |
51 |
{\tt - nfunc:} controls the max. number of simulations |
52 |
within one iteration |
53 |
|
54 |
\paragraph{Summary} ~ \\ |
55 |
From one iteration to the next the descent direction changes. |
56 |
Within one iteration more than one forward and adjoint |
57 |
run may be performed. |
58 |
The updated control used as input for these simulations uses the same |
59 |
descent direction, but different step sizes. |
60 |
|
61 |
\paragraph{Description} ~ \\ |
62 |
From one iteration to the next the descent direction dd changes using |
63 |
the result for the adjoint vector gg of the previous iteration. |
64 |
In lsline the updated control |
65 |
\[ |
66 |
\tt |
67 |
xdiff(i,1) = xx(i-1) + tact(i-1,1)*dd(i-1) |
68 |
\] |
69 |
serves as input for |
70 |
a forward and adjoint model run yielding a new {\tt gg(i,1)}. |
71 |
In general, the new solution passes the 1st and 2nd Wolfe tests |
72 |
so {\tt xdiff(i,1)} represents the solution sought: |
73 |
\[ |
74 |
{\tt xx(i) = xdiff(i,1)} |
75 |
\] |
76 |
If one of the two tests fails, |
77 |
an inter- or extrapolation is invoked to determine |
78 |
a new step size {\tt tact(i-1,2)}. |
79 |
If more than one function call is permitted, |
80 |
the new step size is used together |
81 |
with the "old" descent direction {\tt dd(i-1)} |
82 |
(i.e. dd is not updated using the new {\tt gg(i)}), |
83 |
to compute a new |
84 |
\[ |
85 |
{\tt xdiff(i,2) = xx(i-1) + tact(i-1,2)*dd(i-1)} |
86 |
\] |
87 |
that serves as input |
88 |
in a new forward and adjoint run, yielding {\tt gg(i,2)}. |
89 |
If now, both Wolfe tests are successful, |
90 |
the updated solution is given by |
91 |
\[ |
92 |
\tt |
93 |
xx(i) = xdiff(i,2) = xx(i-1) + tact(i-1,2)*dd(i-1) |
94 |
\] |
95 |
|
96 |
In order to save memory both the fields dd and xdiff |
97 |
have a double usage. |
98 |
% |
99 |
\begin{itemize} |
100 |
% |
101 |
\item [{\tt xdiff}] ~\\ |
102 |
- in {\it lsopt\_top}: used as {\tt x(i) - x(i-1)} for Hessian update \\ |
103 |
- in {\it lsline}: intermediate result for control update |
104 |
{\tt x = x + tact*dd} |
105 |
% |
106 |
\item [{\tt dd}] ~\\ |
107 |
- in {\it lsopt\_top, lsline}: descent vector, {\tt dd = -gg} |
108 |
and {\tt hessupd} \\ |
109 |
- in {\it dgscale}: intermediate result to compute new preconditioner |
110 |
% |
111 |
\end{itemize} |
112 |
|
113 |
\paragraph{The parameter file lsopt.par} |
114 |
|
115 |
% |
116 |
\begin{itemize} |
117 |
% |
118 |
\item {\bf NUPDATE} |
119 |
max. no. of update pairs {\tt (gg(i)-gg(i-1), xx(i)-xx(i-1))} |
120 |
to be stored in {\tt OPWARMD} to estimate Hessian |
121 |
[pair of current iter. is stored in |
122 |
{\tt (2*jmax+2, 2*jmax+3)} |
123 |
jmax must be > 0 to access these entries] |
124 |
Presently {\tt NUPDATE} must be > 0 |
125 |
(i.e. iteration without reference to previous |
126 |
iterations through {\tt OPWARMD} has not been tested) |
127 |
% |
128 |
\item {\bf EPSX} |
129 |
relative precision on xx bellow which xx should not be improved |
130 |
% |
131 |
\item {\bf EPSG} |
132 |
relative precision on gg below which optimization is |
133 |
considered successful |
134 |
% |
135 |
\item {\bf IPRINT} |
136 |
controls verbose (>=1) or non-verbose output |
137 |
% |
138 |
\item {\bf NUMITER} |
139 |
max. number of iterations of optimisation; |
140 |
NUMTER = 0: cold start only, no optimization |
141 |
% |
142 |
\item {\bf ITER\_NUM} |
143 |
index of new restart file to be created (not necessarily = NUMITER!) |
144 |
% |
145 |
\item {\bf NFUNC} |
146 |
max. no. of simulations per iteration |
147 |
(must be > 0); |
148 |
is used if step size {\tt tact} is inter-/extrapolated; |
149 |
in this case, if NFUNC > 1, a new simulation is performed with |
150 |
same gradient but "improved" step size |
151 |
% |
152 |
\item {\bf FMIN} |
153 |
first guess cost function value |
154 |
(only used as long as first iteration not completed, |
155 |
i.e. for jmax <= 0) |
156 |
% |
157 |
\end{itemize} |
158 |
|
159 |
\paragraph{OPWARMI, OPWARMD files} |
160 |
Two files retain values of previous iterations which are |
161 |
used in latest iteration to update Hessian: |
162 |
\begin{itemize} |
163 |
% |
164 |
\item {\bf OPWARMI}: contains index settings and scalar variables |
165 |
|
166 |
{\footnotesize |
167 |
\begin{tabular}{ll} |
168 |
{\tt n = nn} & no. of control variables \\ |
169 |
{\tt fc = ff} & cost value of last iteration \\ |
170 |
{\tt isize} & no. of bytes per record in OPWARMD \\ |
171 |
{\tt m = nupdate} & max. no. of updates for Hessian \\ |
172 |
{\tt jmin, jmax} & pointer indices for OPWARMD file (cf. below) \\ |
173 |
{\tt gnorm0} & norm of first (cold start) gradient gg \\ |
174 |
{\tt iabsiter} & total number of iterations with respect to cold start |
175 |
\end{tabular} |
176 |
} |
177 |
% |
178 |
\item {\bf OPWARMD}: contains vectors (control and gradient) |
179 |
|
180 |
{\scriptsize |
181 |
\begin{tabular}{cll} |
182 |
entry & name & description \\ |
183 |
\hline |
184 |
1 & {\tt xx(i)} & control vector of latest iteration \\ |
185 |
2 & {\tt gg(i)} & gradient of latest iteration \\ |
186 |
3 & {\tt xdiff(i),diag} & preconditioning vector; (1,...,1) |
187 |
for cold start \\ |
188 |
2*jmax+2 & {\tt gold=g(i)-g(i-1)} & for last update (jmax) \\ |
189 |
2*jmax+3 & {\tt xdiff=tact*d=xx(i)-xx(i-1)} & for last update (jmax) |
190 |
\end{tabular} |
191 |
} |
192 |
% |
193 |
\end{itemize} |
194 |
% |
195 |
\begin{figure}[b!] |
196 |
{\footnotesize |
197 |
\begin{verbatim} |
198 |
|
199 |
Example 1: jmin = 1, jmax = 3, mupd = 5 |
200 |
|
201 |
1 2 3 | 4 5 6 7 8 9 empty empty |
202 |
|___|___|___| | |___|___| |___|___| |___|___| |___|___| |___|___| |
203 |
0 | 1 2 3 |
204 |
|
205 |
Example 2: jmin = 3, jmax = 7, mupd = 5 ---> jmax = 2 |
206 |
|
207 |
1 2 3 | |
208 |
|___|___|___| | |___|___| |___|___| |___|___| |___|___| |___|___| |
209 |
| 6 7 3 4 5 |
210 |
|
211 |
\end{verbatim} |
212 |
} |
213 |
\caption{Examples of OPWARM file handling} |
214 |
\label{fig:opwarm} |
215 |
\end{figure} |
216 |
|
217 |
\paragraph{Error handling} |
218 |
|
219 |
|
220 |
|
221 |
\newpage |
222 |
|
223 |
\begin{figure} |
224 |
%\input{part8/lsopt_flow_1} |
225 |
{\scriptsize |
226 |
\begin{verbatim} |
227 |
lsopt_top |
228 |
| |
229 |
|---- check arguments |
230 |
|---- CALL INSTORE |
231 |
| | |
232 |
| |---- determine whether OPWARMI available: |
233 |
| * if no: cold start: create OPWARMI |
234 |
| * if yes: warm start: read from OPWARMI |
235 |
| create or open OPWARMD |
236 |
| |
237 |
|---- check consistency between OPWARMI and model parameters |
238 |
| |
239 |
|---- >>> if COLD start: <<< |
240 |
| | first simulation with f.g. xx_0; output: first ff_0, gg_0 |
241 |
| | set first preconditioner value xdiff_0 to 1 |
242 |
| | store xx(0), gg(0), xdiff(0) to OPWARMD (first 3 entries) |
243 |
| | |
244 |
| >>> else: WARM start: <<< |
245 |
| read xx(i), gg(i) from OPWARMD (first 2 entries) |
246 |
| for first warm start after cold start, i=0 |
247 |
| |
248 |
| |
249 |
| |
250 |
|---- /// if ITMAX > 0: perform optimization (increment loop index i) |
251 |
| ( |
252 |
| )---- save current values of gg(i-1) -> gold(i-1), ff -> fold(i-1) |
253 |
| (---- CALL LSUPDXX |
254 |
| ) | |
255 |
| ( |---- >>> if jmax=0 <<< |
256 |
| ) | | first optimization after cold start: |
257 |
| ( | | preconditioner estimated via ff_0 - ff_(first guess) |
258 |
| ) | | dd(i-1) = -gg(i-1)*preco |
259 |
| ( | | |
260 |
| ) | >>> if jmax > 0 <<< |
261 |
| ( | dd(i-1) = -gg(i-1) |
262 |
| ) | CALL HESSUPD |
263 |
| ( | | |
264 |
| ) | |---- dd(i-1) modified via Hessian approx. |
265 |
| ( | |
266 |
| ) |---- >>> if <dd,gg> >= 0 <<< |
267 |
| ( | ifail = 4 |
268 |
| ) | |
269 |
| ( |---- compute step size: tact(i-1) |
270 |
| ) |---- compute update: xdiff(i) = xx(i-1) + tact(i-1)*dd(i-1) |
271 |
| ( |
272 |
| )---- >>> if ifail = 4 <<< |
273 |
| ( goto 1000 |
274 |
| ) |
275 |
| (---- CALL OPTLINE / LSLINE |
276 |
| ) | |
277 |
... ... ... |
278 |
\end{verbatim} |
279 |
} |
280 |
\caption{Flow chart (part 1 of 3)} |
281 |
\label{fig:lsoptflow1} |
282 |
\end{figure} |
283 |
|
284 |
\begin{figure} |
285 |
%\input{part8/lsopt_flow_2} |
286 |
{\scriptsize |
287 |
\begin{verbatim} |
288 |
... ... |
289 |
| ) |
290 |
| (---- CALL OPTLINE / LSLINE |
291 |
| ) | |
292 |
| ( |---- /// loop over simulations |
293 |
| ) ( |
294 |
| ( )---- CALL SIMUL |
295 |
| ) ( | |
296 |
| ( ) |---- input: xdiff(i) |
297 |
| ) ( |---- output: ff(i), gg(i) |
298 |
| ( ) |---- >>> if ONLINE <<< |
299 |
| ) ( runs model and adjoint |
300 |
| ( ) >>> if OFFLINE <<< |
301 |
| ) ( reads those values from file |
302 |
| ( ) |
303 |
| ) (---- 1st Wolfe test: |
304 |
| ( ) ff(i) <= tact*xpara1*<gg(i-1),dd(i-1)> |
305 |
| ) ( |
306 |
| ( )---- 2nd Wolfe test: |
307 |
| ) ( <gg(i),dd(i-1)> >= xpara2*<gg(i-1),dd(i-1)> |
308 |
| ( ) |
309 |
| ) (---- >>> if 1st and 2nd Wolfe tests ok <<< |
310 |
| ( ) | 320: update xx: xx(i) = xdiff(i) |
311 |
| ) ( | |
312 |
| ( ) >>> else if 1st Wolfe test not ok <<< |
313 |
| ) ( | 500: INTERpolate new tact: |
314 |
| ( ) | barr*tact < tact < (1-barr)*tact |
315 |
| ) ( | CALL CUBIC |
316 |
| ( ) | |
317 |
| ) ( >>> else if 2nd Wolfe test not ok <<< |
318 |
| ( ) 350: EXTRApolate new tact: |
319 |
| ) ( (1+barmin)*tact < tact < 10*tact |
320 |
| ( ) CALL CUBIC |
321 |
| ) ( |
322 |
| ( )---- >>> if new tact > tmax <<< |
323 |
| ) ( | ifail = 7 |
324 |
| ( ) | |
325 |
| ) (---- >>> if new tact < tmin OR tact*dd < machine precision <<< |
326 |
| ( ) | ifail = 8 |
327 |
| ) ( | |
328 |
| ( )---- >>> else <<< |
329 |
| ) ( update xdiff for new simulation |
330 |
| ( ) |
331 |
| ) \\\ if nfunc > 1: use inter-/extrapolated tact and xdiff |
332 |
| ( for new simulation |
333 |
| ) N.B.: new xx is thus not based on new gg, but |
334 |
| ( rather on new step size tact |
335 |
| ) |
336 |
| (---- store new values xx(i), gg(i) to OPWARMD (first 2 entries) |
337 |
| )---- >>> if ifail = 7,8,9 <<< |
338 |
| ( goto 1000 |
339 |
| ) |
340 |
... ... |
341 |
\end{verbatim} |
342 |
} |
343 |
\caption{Flow chart (part 2 of 3)} |
344 |
\label{fig:lsoptflow2} |
345 |
\end{figure} |
346 |
|
347 |
\begin{figure} |
348 |
%\input{part8/lsopt_flow_3} |
349 |
{\scriptsize |
350 |
\begin{verbatim} |
351 |
... ... |
352 |
| ) |
353 |
| (---- store new values xx(i), gg(i) to OPWARMD (first 2 entries) |
354 |
| )---- >>> if ifail = 7,8,9 <<< |
355 |
| ( goto 1000 |
356 |
| ) |
357 |
| (---- compute new pointers jmin, jmax to include latest values |
358 |
| ) gg(i)-gg(i-1), xx(i)-xx(i-1) to Hessian matrix estimate |
359 |
| (---- store gg(i)-gg(i-1), xx(i)-xx(i-1) to OPWARMD |
360 |
| ) (entries 2*jmax+2, 2*jmax+3) |
361 |
| ( |
362 |
| )---- CALL DGSCALE |
363 |
| ( | |
364 |
| ) |---- call dostore |
365 |
| ( | | |
366 |
| ) | |---- read preconditioner of previous iteration diag(i-1) |
367 |
| ( | from OPWARMD (3rd entry) |
368 |
| ) | |
369 |
| ( |---- compute new preconditioner diag(i), based upon diag(i-1), |
370 |
| ) | gg(i)-gg(i-1), xx(i)-xx(i-1) |
371 |
| ( | |
372 |
| ) |---- call dostore |
373 |
| ( | |
374 |
| ) |---- write new preconditioner diag(i) to OPWARMD (3rd entry) |
375 |
| ( |
376 |
|---- \\\ end of optimization iteration loop |
377 |
| |
378 |
| |
379 |
| |
380 |
|---- CALL OUTSTORE |
381 |
| | |
382 |
| |---- store gnorm0, ff(i), current pointers jmin, jmax, iterabs to OPWARMI |
383 |
| |
384 |
|---- >>> if OFFLINE version <<< |
385 |
| xx(i+1) needs to be computed as input for offline optimization |
386 |
| | |
387 |
| |---- CALL LSUPDXX |
388 |
| | | |
389 |
| | |---- compute dd(i), tact(i) -> xdiff(i+1) = x(i) + tact(i)*dd(i) |
390 |
| | |
391 |
| |---- CALL WRITE_CONTROL |
392 |
| | | |
393 |
| | |---- write xdiff(i+1) to special file for offline optim. |
394 |
| |
395 |
|---- print final information |
396 |
| |
397 |
O |
398 |
\end{verbatim} |
399 |
} |
400 |
\caption{Flow chart (part 3 of 3)} |
401 |
\label{fig:lsoptflow3} |
402 |
\end{figure} |
403 |
|