globalization methods for newton

[Christopher Marcotte]

Hello again, sick of me yet?

This is part question about how the newton update for the NLBVP solver works, partly a documentation request, and partly a feature request.

I had the brilliant idea to implement a globalization method for the NLBVP solver newton update, but I'm a little lost how the residual is evaluated in the newton code (i.e.,  when solving N(u*) = 0 for u*, what the rhs is when u is not u*).

From looking at the code, I thought A = N'(u), and b = -N(u). Is this correct?

I started with the simplest two modifications I knew of.

I tried steepest descent, where x = np.dot(A,b), which yields a bizarre type error:

np.copyto(self.data[pencil.local_index], data)

TypeError: Cannot cast scalar from dtype('O') to dtype('float64') according to the rule 'same_kind'

I tried using np.linalg.pinv to solve x = np.dot(np.linalg.pinv(A),b), which was very helpful in my own newton solver, but even a modification that simple throws an error.

numpy.linalg.linalg.LinAlgError: 0-dimensional array given. Array must be at least two-dimensional

Neither of these are actually good options (e.g., steepest descent only makes sense for nominally infinite dimensional problems as a continuous-fictitious-time method which preconditions a further globalized Newton step, like hookstep or backtracking line-search see https://doi.org/10.1017/jfm.2016.203 ) but I can't even implement the good ones without some more thorough documentation for the newton step. What would almost certainly be better is, rather than implementing a globalized newton method myself, the evaluation of the residual and jacobian be standardized behind some lambda functions, and they are passed to one of the optimization routines in scipy.optimize for the minimization of a quadratic model. Or perhaps there are more appropriate scipy solvers I'm not familiar with.

Apologies for the length. Thoughts?

CM

[Keaton Burns]

Hi Chris,

No problem, it’s great to see more interest in the NLBVP!  I think you’re understanding of how Dedalus is doing things is essentially correct, and the issues you’re seeing are due to using numpy functions on scipy sparse arrays.

To make sure we’re on the same page, the docstring on the NLBVP problem class here describes the problem setup.  We essentially start with N(X) = 0, but split the linear and nonlinear terms like N(X) = L.X - F(X) = 0 ==> L.X = F(X), so the solver only has to recompute the Jacobian of F, and not rebuild the L matrices, every iteration.  This problem is written as a BVP for a perturbation/update around some initial state as (L - F’(X0)).X1 = F(X0) - L.X0, which is equivalent to N’(X0).X1 = -N(X0), or simply A.X1 = b for the update.

 

The Newton step routine can be found here.  It uses the current values of solver.state to compute N(X0) and N’(X0) (lines 236 and 238).  Then it builds the A matrix and b vectors using this data (lines 243-244), and solves for the update X1 in line 245.  The routine then adds these perturbations back to the state vector.

The matrix A is a scipy sparse matrix, and b is a dense numpy vector.  Scipy sparse matrices are not natively supported in must regular numpy functions — in particular, np.dot(A,b) with these types produces a numpy object array, which causes the first error, and np.linalg.pinv(A) directly throws the second error.

Depending on your problem (and particularly the NCC cutoff), matrix A may still be fairly sparse, in which case it might be better to use the scipy sparse linear algebra tools and sparse matrix methods, like A.dot(_).  Otherwise, you can use the A.todense() method to convert to a dense numpy array, and then use any numpy linear algebra functions.

Everything else should still work, with A and b from lines 243-244 representing the Jacobian and residual on each iteration, no matter how you choose to determine the update x in line 245.

Hope that helps, and definitely give an update if you run into more issues or get something working with one of the scipy optimization routines!

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/13f8e983-46a9-47f9-9cfe-e0b89225a527%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Christopher Marcotte]

So to evaluate the residual of the nonlinear functions one has to first gather the state, update the state with the tentative update, scatter the state, then evaluate the rhs? Is there some convenient mechanism to evaluate a tentative state by computing a second residual vector and second jacobian? As is, the newton update seems very hard to extend. Maybe I'm just very unfamiliar with python, but even my dumb function to update the state and return A & b doesn't work with this newton function.

[Christopher Marcotte]

Attached is the eval function and the modified newton method.

[Keaton Burns]

Hi Chris,

Yeah I think there’s a couple python issues here.  

In “eval_residual_jacobian”, setting “oldstate = self.state.data” is making a local reference to that array object, not making a copy.  In the next line, that array is updated in place, so the old state data is never actually stored.  The easiest fix is probably to just do “solver.state.data -= dx” at the end, instead of trying to carry around a copy of the old state vector.

If you’re adding both of these functions as methods on the NLBVP solver, you’ll want to cross reference them like “self.eval_residual_jacobian(alpha*dx)” instead of "eval_residual_jacobian(self, alpha*dx)”.

Finally, when the “while” condition gets achieved, the function enters an infinite loop.  I think you need to update the entries of W with the new values of alpha here.

Best,

-Keaton

On February 23, 2017 at 4:37:07 AM, Christopher Marcotte (christophe...@gmail.com) wrote:

Attached is the eval function and the modified newton method.

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/fb4b8dc7-2041-4aea-8a37-55809485c593%40googlegroups.com.

[Christopher Marcotte]

Hi Keaton,

Thanks for looking at the code. Is there a standard way to get a copy of the current state from the solver, and then evaluate the residual and Jacobian at a new state vector? All I'm trying to do is implement damping for the Newton step.

CM

[Keaton Burns]

Hi Chris,

Sorry for the long delay.  The NLBVP solver is pretty bare-bones right now, just the single newton-iteration method.  But I think you have the right approach here, there were just those minor python issues getting in the way.  We can definitely add a function like this to the codebase to assist with interfacing with other optimizers, etc.

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/77301d10-a0eb-40c8-9412-324591376f3e%40googlegroups.com.

[Christopher Marcotte]

Hi Keaton,

That would be very helpful, thanks!

CM