ISMP 2009
The last few week have been very busy as I prepared for (and attended) the International Symposium on Math Programming. It was kind of a big deal, and I was lucky enough to meet many big names in my field, and be able to give a talk on my research. My talk focused on the theoretical aspects of Newton's method for optimization with line search regularization. It was definitely a fun talk for me, and I got lots of good comments from the people who were there (and lots of "I missed it!" from those who weren't). If you missed it (or not), and want a copy of my slides, you can find them here.KASTR - GMRES-Steihaug in PETSc
In the last few weeks, I've been working on getting the GMRES-Steihaug algorithm set up as a nonlinear solver (through SNES) in PETSc. The result is something slightly different than the original GMRES-Steihaug, but keeps the main ideas, and is called KASTR - KSP Aware Steihaug Trust Region. The main advantages of this (as I see them), are independence of any particular KSP solver or preconditioner, use of canned (known good) implementations for everything outside our algorithm, and easy manipulation of many of the parameters on the command line.Details on what's changed in the implementation, as well as preliminary numerical results, below.
Not Solving, But Doing Better
A somewhat short addendum to my last post, there is one bit of information that I was able to tease out of my last batch of data that gives a little more insight into how it's doing. In addition to reporting how many failures occurred, it was easy to also get the relative residual for the failure cases (i.e. if the method failed, how much progress did it make?).Results after the jump.
Solving or Avoiding the Problem?
As I mentioned back in my motivation post, the goal of the GMRES-Steihaug algorithm is to address the primary pitfall of the line search and dogleg strategies for nonlinear systems. In particular, we note that in the line search case, the Newton step can be arbitrarily long and be in a direction of poor decrease, eventually leading to a line search failure. On the other hand, the dogleg approximation of the trust region solution can tend to get bogged down and take short steepest descent-like steps. In this case, the iteration simply runs out of time (allowed iterations) without making significant progress.Our early numerical results seem to imply that we are solving our model problem better, but the goal of this post is to attempt to satisfactorily answer, "Are we solving these issues, or simply avoiding problem areas better?"
Good News, and Bad
Recently, I've been working on a rather critical component of the GMRES-Steihaug algorithm: When to stop iterating on J s = -F. At first, this seemed like a very straightforward issue, we simply stopped when the iterate crossed the trust region. The catch, however, is that the inexact Newton step is not monotonically increasing in length during the GMRES iteration, i.e. the step may be very long, then become very short. The reason this is a problem is that if we stop too early (when the step first crosses the trust region), we may miss the Newton step (if the Newton step is interior to the trust region).So, I sat down to do some proofs as to how bad it could be...
Flurry of Activity
Leading up to my thesis proposal, there was a flurry of activity in my work, but I didn't really get the chance to organize it outside the context of my thesis proposal. Unfortunately, this lead to my research blog lagging quite a bit behind (i.e. I didn't post at all). So, to make up for it, I'm going to try to separate the most important bits from my thesis proposal and post them here over the course of the next week or so. Towards the end, I'll be posting some of the things that didn't make it into my proposal, and finally the stuff that I'm working on now.Stay tuned!
New Numerical Results
Here are my new numerical results. This is for the same problem as before, with the Jacobian fixes I mentioned previously. Each cell is an average of four different starting points, to average out the case of a particularly good or bad starting point. In digging down into the data, the most notable cause of failure (or simply a large number of nonlinear iterations) was lousy Newton steps, i.e. GMRES(32) failed to converge after 32 cycles. This is most likely due to the fact that I am using no preconditioning, as GMRES is known to stagnate in this case. This brings up the point that in particular for the line search case, the performance with preconditioning may be substantially better as the Newton equations would be solved more accurately.Numerical Results: Fly in the Ointment
As I was going through my results in detail the other day, I found a noticeable problem, at certain points the dogleg step was failing, not because of a large Newton direction, but because the Cauchy point wasn't providing decrease for very small trust region. This jumps out because trust region theory is based on exactly this - when the trust region is small, the Cauchy point imitates the direction of the negative gradient which should provide at least some decrease.Taking a note from PETSc, I added a very short bit of code to generate a finite difference Jacobian and compared it to the analytic Jacobian I wrote by hand (element-by-element). Sure enough, there were a few entries that were incorrect (by an embarrassing amount). However, knowing which elements were off, made for a quick code fix (as it happens, I was imposing two incompatible border conditions on the same element of F, and that manifested as extraneous elements in J).
With the code updated, I'm rerunning all my simulations, updates soon.
Numerical Results
Using the algorithms provided in my previous post, as well as typical line search and dogleg codes, I ran a series of tests against our model problem. For these tests, the line search and dogleg Newton steps were computed using full GMRES with relative error 10^-5 - while somewhat expensive in this context, it ensured that there wouldn't be significant error due to an inexact Newton step. Conversely, the above algorithms used GMRES with starting restart parameter 32, and a maximum of 32 iterations. The model problem was discretized on an 8x8 mesh, giving 256 equations and unknowns. The lid velocity and Grashof number were fixed at 100 and 10000 respectively, while the Prandtl number was taken as 1, 2, 4, and 8. Additionally, the four different starting points were attempted - the zero vector, the vector of all 1s, and two random vectors with elements chosen iid on [-1,1]. For each of the 64 tests, the number of nonlinear iterations and the wall time (as given by the UNIX time command) are listed as well as the averages for each algorithm across the four different starting points.Below is a table summarizing the results.
An initial algorithm
Here's an initial algorithm:Model Problem
As a means to test any nonlinear solver, I've downloaded and installed PETSc - a fairly standard toolkit for solving all sorts of problems. One of the examples implements a nonlinear solver for a lid driven cavity problem with heat transport, namely:-Lap(U) - Grad_y(Omega) = 0 -Lap(V) + Grad_x(Omega) = 0 -Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0 -Lap(T) + PR*Div([U*T,V*T]) = 0
On the unit square.
Motivation
As a start to my research blog, I think it's worth mentioning the primary motivation of my work (at least in my own mind). In the realm of nonlinear systems, there are two primary approaches to globalization of the Newton iteration - Line search and trust region. Each approach has their strengths and weaknesses which I will attempt to summarize below.
|
|
Meta
