Constrained optimisation on a hypersphere

Question

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

Answer 1

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.

So hopefully this answers your three questions:

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.