# Constrained Newton-Raphson root finding

Original Question

I have a set of non-linear equations and I need to find the root where a subset of my solution vector is constrained to be greater than or equal to 0. I have implemented the Newton-Raphson algorithm but I am finding that some of the quantities that I need to remain positive are going negative. I am familiar with using Lagrange multipliers to enforce constraints in an optimization problem but I'm unsure about how to do the same for root finding.

It had occurred to me that I could map my variables which must be positive to log space (i.e., solve for $$\log(x_i)$$ rather than $$x_i$$) but since zero is an important case I was worried that this might cause numerical problems.

From a previous answer it looks like an interior point method might be appropriate but it would involve modifying my solver rather than just the residual equation which I would like to avoid.

Further Description

The equations are a series of non-linear equations arising from a material model which incorporates irrecoverable deformation ( $$E^p$$ ). The way that this is typically achieved is through internal state variables ($$\xi$$) which are governed through some evolution equations which are determined by some evolution rate ( $$\dot{\gamma}$$ ).

The evolution equations take the form of:

$$\dot{E^p} = \dot{E}\left(E^p,\xi, \dot{\gamma}\right)$$

$$\dot{\xi} = \dot{\xi}\left( \xi, \dot{\gamma}\right)$$

Subject to some additional onset condition ( yield function ):

$$F=F(E - E^p)$$

where $$E$$ is the total deformation which is a function which is <= 0. The problem is subject to a Kuhn-Tucker condition which can be summarized as:

$$\dot{\gamma}F = 0$$

Implying that the evolution rate is zero if $$F < 0$$ and $$\dot{\gamma}>0$$ if $$F=0$$.

My unknown vector is something like:

$$x=\left[ E^p, \xi, \gamma \right]$$

My residual is something like:

$$R=\left[{ E }^{p,expected} - E^p, \xi^{expected} - \xi, \dot{\gamma}F + \langle F \rangle \right]$$

Where the values marked $$\left(\cdot\right)^{expected}$$ are the results of the evolution equations and $$\langle \cdot \rangle$$ are the Macaulay brackets which are defined as:

$$\langle x \rangle = \frac{1}{2}\left( x + abs( x ) \right)$$

I include the Macaulay term because, otherwise, the case of F > 0 and $$\dot{\gamma}=0$$ would be accepted which is incorrect. I tried using if statement in the residual calculation but that caused other issues and this seems to be performing better.

• Since it is a non-linear system of equations, it is not surprising that there are multiple roots with some not satisfying the constraints. Also it wouldn't be weird to jump to a solution outside of the regions due to basins of attraction phenomenon. I guess you could take the square of the objective function and solve it like a minimization problem with constraints using Lagrange multipliers. At least, that is what I would try first. – Abdullah Ali Sivas May 11 at 22:52
• @NateM Global optimization finds maxima/minima of a function over a given set; this is closely related to the root finding problem since a root minimizes abs(f). – Maxim Umansky May 13 at 0:38
• Could you tell us more about the system of nonlinear equations? – Brian Borchers May 13 at 1:25
• @MaximUmansky Okay, I think I did understand what you meant but I still don't think this will solve the problem. The negative root IS a solution to the system of equations but it is not allowable because the variables must remain positive. I want to impose that constraint in the system of equations. If I recast the problem as an optimization problem ( in the absence of constraints ) both the correct solution and the incorrect would have the same objective value. – NateM May 13 at 13:51
• Interior point methods are typically used to solve complementarity problems like these. These methods compute Newton directions toward a solution to a perturbed system of equations and then restrict the step length to prevent the variables from going negative. You should probably look at solvers for this class of problems rather than trying to implement your own solver from scratch. – Brian Borchers May 13 at 16:37

I thought I would come back to this now that I have an answer that works. The problem was multifaceted and I tried several solutions that worked to varying degrees before eventually settling on a solution.

Things I tried that didn't work on their own:

1. Homotopy solver that tried to implement barriers

$$g( x, s ) = (1 - \frac{ 1 } { a( s ) } ) R( x ) + \frac{ 1 }{ a( s ) } b( x, s )$$

where

$$b( x, s ) = e^{ s a ( b - x ) } - 1$$

$$a( s ) = e^{ A s }$$

Where $$R$$ is the original residual, $$s$$, is the pseudo time, and $$A$$ is some value large enough to make sure the barrier holds. I found that this gave a pretty smooth variation from the barrier function to the true function if A is at least 10 or so. This worked for toy problems but not for my actual problem.

2. Newton homotopy solver:

$$g( x, s ) = R( x ) + ( 1 - s ) R( x_0 )$$

I like this homotopy and ended up using it for my final non-linear equation solve. In the solve I first try s = 1 and then cutback if required.

3. Performing multiple nested Newton-Raphson solves. Instead of trying to solve for $$E^p$$, $$\xi$$, and $$\dot{\gamma}$$ simultaneously I solved for $$\xi$$ analytically and $$E^p$$ numerically for a given value of $$\dot{\gamma}$$ as each iteration to solve for $$\dot{\gamma}$$. This dramatically helped with my convergence for many problems but didn't solve the more difficult cases.

4. Worked through my code and found a bug which was likely causing unintended behavior. Again, this helped, but didn't solve my problems.

5. Tried limiting the line-search in the Newton Raphson to always be in the correct domain. This helped, but sometimes didn't allow for convergence.

Now for what actually worked:

I had unintentionally set up my residual equation such that if the error was too great it would pass over into a new basin of attraction which pulled it in the wrong direction. I needed to modify my residual equation to prevent this behavior. Note that in the original residual for $$\dot{\gamma}$$

$$R^{\dot{\gamma}} = \langle F \rangle + \dot{ \gamma } F$$

The problem with this is that if $$F$$ gets very large, eventually the Jacobian can go positive due to the second term which will try to drive $$\dot{\gamma}$$ negative. So in the end I needed to:

1. Use a Newton homotopy solver ( this isn't necessarily required right now but helps when the problem gets stiff ).
2. Break the monolithic solve into sub-problems.
3. Modify my residual to exclude the problematic behavior.

In the context of the original equation, I believe 3 is probably the crucial point.