On Wed, Aug 1, 2012 at 5:01 PM, Todd Munson <[email protected]> wrote:
On Aug 1, 2012, at 4:13 PM, Jed Brown wrote:
On Wed, Aug 1, 2012 at 11:10 AM, Todd Munson <[email protected]> wrote:
Yes, I am all in favor of actually specifying a VI externally and then internally formulating the appropriate augmented system. This is what I do in other contexts (i.e. PATH for GAMS/AMPL). Some things you probably want to do with such an interface is differentiate between linear and nonlinear constraints -- you do need second derivatives for nonlinear constraints and they should be convex.
Why not treat it as always nonlinear with the trivial second derivative in the linear case?
You can if you want. I would say that one benefit of splitting is that sometimes roundoff errors can be an issue and isolating the linear parts (which are just data) from the nonlinear parts may help. It depends on your Newton method though. Knowing that your constraints are linear also means that you do not need to worry about the Hessian of the Lagrangian. Adding and dropping linear constraints is would also seem to be easier since there is no penalty to be paid since the Jacobian does not need to be modified.
One interesting thing is that for nonlinear constraints, your Jacobian matrix in the upper right hand block is the summation of a general matrix (the Jacobian of F(x)) and a symmetric matrix (sum_i lambda_i nabla^2 c_i(x), where lambda is the Lagrange multiplier). Is there a good matrix storage scheme for this combination?
How does one precondition a system like that? Presumably, the user may at best be able to provide a preconditioner for JF. If JF responds well to multigrid, then it may be possible to use AMG-KKT to coarsen the constraint Hessian, but how does MG perform for the sum JF + Hc? There is also the problem of smoothing, although, if JF and Jc are sparse, then DD may be an option with a direct solve on subdomains. Anything more definitive on preconditioning these systems? Dmitry.
You can start getting very general (i.e. linear and nonlinear cone constraints, which can appear in contact problems with friction), but I would advocate for simplicity rather than generality.
I like simplicity, but I also like not having to change interfaces. I would rather write an interface for a more general case and even temporarily have an implementation that can't do the fully general thing, eventually perhaps doing something more general. We put a rather large amount of effort into making TS work with a convenient (solver-friendly) representation for DAEs, even though lots of methods cannot be used for DAEs. If the general case is not too complicated for the user, I'd rather do it that way and internally reformulate.
The problem with general cones is that they are not easy (you need both a definition of the cone and the polar cone) and the numerical methods are problematic (second-order cones, for example, are not differentiable at the origin). Here you require interior-point methods. So, the cone constraints algorithms themselves are more of a research issue.
The addition and deletion of nonlinear constraints in your example is fine. For an active set method, the constraints can be implicitly defined and you only need to bring in the new active nonlinear constraints and add them to the constraint set for the linear algebra. You can also do this for other methods. So the active nonlinear constraints can change from one Newton iteration to the next. The only place where the inactive constraints matter is in the evaluation of the residual for the line search and termination tests.
Cool, I'm glad you're happy with the math. Now we need a suitable interface, which I think means something a fair amount different from what we have now. We need constraint function, its Jacobian, and a Hessian, each of which can change size?
Yep. I leave it to you to figure this out though.
For the contact problems with friction, you typically do an explicit method though. In this setting, you know the subset of the possible active constraints at each time step and can keep them constant in the VI solver. They only change as you step forward through time. So you can optimize for this case and not have to worry about the constraints changing within each Newton iteration.
I'm not sure how you want to design the variable constraint functionality, but there is nothing wrong with it technically.
I can explain the semismooth method with slacks thing in a few lines. The simple problem is:
0 <= x perp f(x) >= 0
The compact semismooth formulation is:
Phi(x,f(x)) = 0
where Phi is the reformulation used.
The semismooth formulation with slacks is:
f(x) - s = 0 Phi(x,s) = 0
The linear algebra for the semismooth reformulation with slacks is very similar to that for the compact formulation, but there are subtle differences in the right hand side and the subdifferential.
Do you mean that the slacks are eliminated again in the linear algebra? What are the "subtle differences"?
If you eliminate the slack variables, in the slack formulation, then you get a system that looks like the semismooth system, but with a different right hand side. Apparently, the difference in the right hand side leads to fewer Newton iterations in general. Not a big deal.
Cheers, Todd.