# Constrained optimisation on a hypersphere

I am currently trying to implement a GS2 (Gonzalez-Schlegel second order) IRC algorithm in a python code. I am following the original paper ref(1).

The main problem is in the constrained optimisation step of the GS2 algorithm. In this step, the geometry optimisation is done by holding the distance of the coordinates (euclidean distance) constant from the coordinates of a pivot point (fixed throughout optimisation).

I will reproduce the main equations for the GS2 constrained optimisation step here: $$p' = x' - x^* \tag{1}$$ $$\Delta x = x - x'\tag{2}$$ $$p = p' + \Delta x = x - x^* \tag{3}$$

where $$x'$$ and $$x$$ are the old and new coordinates for the constrained optimisation; $$x^*$$ is the pivot point that is constant. The $$p'$$ and $$p$$ similarly are the vectors from the old and new coordinates to the pivot point, respectively.

The energy is minimised under the constraint $$p^T\cdot p = k$$ (where k is some fixed scalar number). This is equivalent to optimization on the surface of a hypersphere. They express the optimisation problem with a Lagrange multiplier form: $$L = E' + g^T \cdot \Delta x + \frac{1}{2} \Delta x^T\cdot H \cdot \Delta x - \frac{1}{2} \lambda [p^T\cdot p - k]\tag{4}$$

Where the potential energy surface is expressed as a truncated Taylor series (as normal for QM stuff) and the $$g$$ and $$H$$ are the gradient vector and Hessian matrix as usual.

What I don't understand here is the half sign in front of the Lagrange muliplier $$\lambda$$.

Also, after this, they mention that at convergence, that $$\Delta x = - (H - \lambda I)^{-1} \cdot (g- \lambda p')\tag{5}$$

They then put that expression into the expression for $$p$$ (eqn. 3) and find the $$\lambda$$ (by doing a bracketing 1D scalar root search) for which the $$\Delta x$$ satisfies the constraint.

How does this expression follow from the above expression? I don't quite understand. Could anyone explain how this derivation works mathematically?

Also, another problem with these equations is that the step size is too high in most iterations so there is unpredictable behaviour. Is there any way to restrict the step size for the constrained optimisation step? (In the paper, they have a linear interpolation formula to stabilize the quasi-NR step, but it does not help too much)

References: (1) C. Gonzalez, H. B. Schlegel, J Chem. Phys., 1989, 90(4), 2155

• Perhaps you forgot to add the reference paper (1) Mar 18 at 19:27
• @VandanRevanur Thanks for catching that! fixed it. Mar 19 at 9:54

I'll propose a slightly different Lagrangian expression:

$$L = E' + g^T \Delta x + \frac{1}{2} \Delta x^T H \Delta x + \frac{1}{2} \gamma \Delta x^T \Delta x - \frac{1}{2} \lambda (p^T p - k). \tag{1}$$

Here $$\gamma$$ is really a spring constant (with the corresponding units) that penalizes large norms for $$\Delta x$$. The sign of $$\lambda$$ doesn't really matter since it's enforcing a constraint, so we keep it as per the expression given above in equation 4.

Now, we would like to minimize the Lagrangian with respect to $$\Delta x$$, and we will do so by completing squares. First we substitute in $$p = p' + \Delta x$$ (equation 3), and collect terms in like powers of $$\Delta x$$:

$$L = \frac{1}{2}(\Delta x^T[H + \gamma I - \lambda I]\Delta x) - (\lambda p' - g)^T \Delta x + \cdots \tag{2}$$

where the dotted terms are constant with respect to $$\Delta x$$. The term linear in $$\Delta x$$ combines a contribution from the forces $$g$$ and a contribution from the $$p^Tp$$ term of equation 4.

Now hark back to your high school algebra: a quadratic expression $$ax^2 - bx + c$$ is minimized at $$x^* = (2a)^{-1} b$$ (by completing the square), and this is no different:

$$\Delta x^* \equiv \mathrm{arg min}\,L(\Delta x) = (H + \gamma I - \lambda I)^{-1}(\lambda p' - g) \tag{3}$$

Other than my proposed modification (equivalently, when $$\gamma = 0$$) this expression is exactly equation 5.

1. Why half $$\lambda$$? Scaling your quadratic term by half leaves you with just $$\lambda$$ in the final answer, rather than $$2 \lambda$$. It's purely aesthetic -- but it also does leave you with one less factor of two to muck up when you're programming.
2. How to get to equation 5? Expand in terms of $$\Delta x$$ and complete the square.
3. How to constrain the step size? To me the most natural way is to put in an explicit norm restraint. From my final equation you can see that when $$\gamma \gg \lambda$$ and $$\gamma \gg ||H||$$, $$\Delta x^*$$ will just be $$(\lambda p' - g)/\gamma$$ -- that is, you'll simply be following the forces downwards, in ever decreasing step size as $$\gamma$$ increases.
• @SRMaiti (1) have you tried deriving it yourself? You need to use a combination of expanding squares (which works just like expanding $(a+b)^2$ in normal algebra), recognising that for any vector $v$, $v^T v = v^T I v$, and grouping like matrix multiplications. Gamma is a constant you choose. If it is large enough the effective Hessian $H + \gamma I$ is always positive definite (answering your last question). It looks like I reinvented Hessian modification (nmr-relax.com/manual/Hessian_modifications.html). Mar 21 at 9:32
• Lambda is calculated by re-imposing the constraint $p^T p = k$ on the solution for $\Delta x^*$ and finding the value which permits that. You will usually need some sort of iteration to accomplish this. Mar 21 at 9:35
• Once again, $\gamma$ is a constant for each minimisation step. It is a number you specify. So if you use equation (3) of my answer you have one (modified) constraint in one variable, $\lambda$, not two constraints in two variables. You may find it easier to think of this as a Hessian modification $H \rightarrow H+\gamma I$, and if $\gamma$ is larger than the most negative (or zero) eigenvalue of $H$, then $H+\gamma I$ is guaranteed positive definite (and clearly still symmetric). Mar 21 at 11:30