Minimization of energy

References

[1] (1,2,3)

Nocedal, J. and Wright, S.J. Numerical optimization. New York, NY: Springer New York. eds., 1999.

[2] (1,2,3)

Ivanov, A.V., Uzdin, V.M. and Jónsson, H., 2021. Fast and robust algorithm for energy minimization of spin systems applied in an analysis of high temperature spin configurations in terms of skyrmion density. Computer Physics Communications, 260, p.107749.

On this page we recall the approach and main formulas from [1] and [2]. The minimization of the magnetic ground state in magnopy is implemented with the method described in [2].

Minimization of the energy function can be formulated as a problem of minimizing the function over the \(M\) vectors of the spin directions \(\boldsymbol{z}_{\alpha}, \alpha = 1, ..., M\)

\[E = F(\boldsymbol{z}_{1}, ..., \boldsymbol{z}_{M})\]

Directional vectors are unitary vectors and vary on the sphere. This fact introduces complications in the minimization procedure as the optimization space is not a vector space and the typical (BFGS, for instance) algorithms for linear optimizations can not be applied directly. This problem is elegantly solved via parametrization of those vectors with the exponents of skew-symmetric matrices [2]. Given an initial guess \(\boldsymbol{z}_{\alpha}^{(0)}\), any other set of directional vectors can be obtained by the following formulae:

\[ \boldsymbol{z}_{\alpha} = e^{\boldsymbol{A}_{\alpha}} \boldsymbol{z}_{\alpha}^{(0)} = \cos(\theta) \boldsymbol{z}_{\alpha}^{(0)} + \sin(\theta) (\boldsymbol{r}\times\boldsymbol{z}_{\alpha}^{(0)}) + (1 - \cos(\theta))\boldsymbol{r} (\boldsymbol{r}\cdot\boldsymbol{z}_{\alpha}^{(0)})\]

where \(\boldsymbol{E}_{\alpha}\) are skew-symmetric matrices and

\[\theta = \sqrt{\left(a_{\alpha}^x\right)^2 + \left(a_{\alpha}^y\right)^2 + \left(a_{\alpha}^z\right)^2}\]

\[\boldsymbol{r} = \dfrac{(a_{\alpha}^x, a_{\alpha}^y, a_{\alpha}^z)}{\theta}\]

The energy function can be rewritten as:

\[E = F( e^{\boldsymbol{A}_1} \boldsymbol{z}_{1}^{(0)}, ..., e^{\boldsymbol{A}_I} \boldsymbol{z}_{I}^{(0)} )\]

Skew-symmetric matrices are parametrized by three numbers as

\[\begin{split}\boldsymbol{A}_{\alpha} = \begin{pmatrix} 0 & -a_{\alpha}^z & a_{\alpha}^y \\ a_{\alpha}^z & 0 & -a_{\alpha}^x \\ -a_{\alpha}^y & a_{\alpha}^x & 0 \end{pmatrix}\end{split}\]

Finally the energy can be written as a function of vector \(\boldsymbol{x}\) from the vector space \(\mathbb{R}^{3I}\):

\[E = F(\boldsymbol{x}) \qquad \boldsymbol{x} =( a_{1}^x, a_{1}^y, a_{1}^z, ..., a_{I}^x, a_{I}^y, a_{I}^z )\]

Then energy of the system is minimized with the BFGS algorithm [1].

Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm

Formula for the inverse Hessian update:

\[H^{ij}_{k+1} = \sum_{u,v}(\delta_{i,u} - \rho_ks^i_ky^u_k) H^{uv}_k (\delta_{v,j} - \rho_ky^v_ks^j_k) + \rho_k s^i_ks^j_k, \qquad \rho_k = \dfrac{1}{\sum_i y^i_k s^i_k}\]

Given Initial guess \(\boldsymbol{x}_0\) and Initial approximation of the inverse hessian matrix \(\boldsymbol{H}_0\),

1. \(k \gets 0\)

2. While convergence is not achieved:

   1. Compute the gradient of the function \(\boldsymbol{\nabla} F(\boldsymbol{x}_k)\);

   2. Compute the search direction \(\boldsymbol{p}_k = -\boldsymbol{H}_k \boldsymbol{\nabla} F(\boldsymbol{x}_k)\);

   3. Compute length of the step \(\alpha_k\) via Line search;

   4. Set \(\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \alpha_k \boldsymbol{p}_k\) and compute gradient \(\boldsymbol{\nabla} F(\boldsymbol{x}_{k+1})\);

   5. Set \(\boldsymbol{s}_k = \boldsymbol{x}_{k+1} - \boldsymbol{x}_k\) and \(\boldsymbol{y}_k = \boldsymbol{\nabla} F(\boldsymbol{x}_{k+1}) - \boldsymbol{\nabla} F(\boldsymbol{x}_k)\);

   6. Update the hessian matrix \(\boldsymbol{H}_{k+1}\) by the BFGS formula;

   7. \(k \gets k + 1\).

Note

In our implementation we update the direction vectors at the end of each iteration (i.e. at step 2.g). Therefore, the vector \(\boldsymbol{x}_k\) is always equal to \(( 0, 0, 0, 0, 0, 0, ..., 0, 0, 0)\).

Initial guess

Initial guess is provided by the user or randomly generated. User provides three components for each direction vector \((z_{\alpha}^x, z_{\alpha}^y, z_{\alpha}^z)\).

Initial approximation of the inverse hessian matrix

We take an identity matrix as an initial approximation of the hessian matrix and scale it as

\[\boldsymbol{H}_0 = \dfrac{\boldsymbol{y}^T_k\boldsymbol{s}_k}{\boldsymbol{y}^T_k\boldsymbol{y}_k} \boldsymbol{I}\]

before the first update [1].

Gradient of the function F(x)

As we choose to update the direction vectors at each step of the BFGS algorithm, then the gradient with respect to these variables can be computed as

\[\dfrac{\partial F}{\partial\boldsymbol{a}_{\alpha}} = \boldsymbol{t}_{\alpha} = \boldsymbol{z}_{\alpha} \times \dfrac{\partial E^{(0)}}{\partial\boldsymbol{z}_{\alpha}}\]

where \(\boldsymbol{t}_{\alpha}\) is a torque vector and \(\boldsymbol{a}_{\alpha} = (a_{\alpha}^x, a_{\alpha}^y, a_{\alpha}^z)\).

The gradient of the energy is computed analytically

\[\begin{split} \dfrac{\partial E^{(0)}}{\partial z^t_{\alpha}} =& \,C_1 J_1^t(\boldsymbol{r}_{\alpha}) S_{\alpha} +\\&+ 2C_{2,1} \sum_{j} J_{2,1}^{tj}(\boldsymbol{r}_{\alpha}) z^j_{\alpha} (S_{\alpha})^2 +\\&+ 2C_{2,2} \sum_{\beta, \nu, j} J_{2,2}^{tj}(\boldsymbol{r}_{\nu,\alpha\beta}) z^j_{\beta} S_{\alpha} S_{\beta} +\\&+ 3C_{3, 1} \sum_{j, u} J^{tju}_{3, 1}(\boldsymbol{r}_{\alpha}) z^j_{\alpha} z^u_{\alpha} (S_{\alpha})^3 +\\&+ 3C_{3, 2} \sum_{\beta, \nu, j, u} J^{tju}_{3, 2}(\boldsymbol{r}_{\nu,\alpha\beta}) z^j_{\alpha} z^u_{\beta} (S_{\alpha})^2 S_{\beta} +\\&+ 3C_{3, 3} \sum_{\substack{\beta, \gamma, \\ \nu, \lambda, j, u}} J^{tju}_{3, 3}(\boldsymbol{r}_{\nu,\alpha\beta}, \boldsymbol{r}_{\lambda,\alpha\gamma}) z^j_{\beta} z^u_{\gamma} S_{\alpha} S_{\beta} S_{\gamma} +\\&+ 4C_{4, 1} \sum_{\alpha, j, u, v} J_{4, 1}^{tjuv}(\boldsymbol{r}_{\alpha}) z^j_{\alpha} z^u_{\alpha} z^v_{\alpha} (S_{\alpha})^4 +\\&+ 4C_{4, 2, 1} \sum_{\substack{\beta, \nu, \\ j, u, v}} J_{4, 2, 1}^{tjuv}(\boldsymbol{r}_{\nu,\alpha\beta}) z^j_{\alpha} z^u_{\alpha} z^v_{\beta} (S_{\alpha})^3 S_{\beta} +\\&+ 4C_{4, 2, 2} \sum_{\substack{\beta, \nu, \\ j, u, v}} J_{4, 2, 2}^{tjuv}(\boldsymbol{r}_{\nu,\alpha\beta}) z^j_{\alpha} z^u_{\beta} z^v_{\beta} (S_{\alpha})^2 (S_{\beta})^2 +\\&+ 4C_{4, 3} \sum_{\substack{\beta, \gamma, \\ \nu, \lambda, \\ j, u, v}} J_{4, 3}^{tjuv}(\boldsymbol{r}_{\nu,\alpha\beta}, \boldsymbol{r}_{\lambda,\alpha\gamma}) z^j_{\alpha} z^u_{\beta} z^v_{\gamma} (S_{\alpha})^2 S_{\beta} S_{\gamma} +\\&+ 4C_{4, 4} \sum_{\substack{\beta, \gamma, \varepsilon, \nu, \lambda, \rho, \\ \\ j, u, v}} J_{4, 4}^{tjuv}(\boldsymbol{r}_{\nu,\alpha\beta}, \boldsymbol{r}_{\lambda,\alpha\gamma}, \boldsymbol{r}_{\rho,\alpha\varepsilon}) z^j_{\beta} z^u_{\gamma} z^v_{\varepsilon} S_{\alpha} S_{\beta} S_{\gamma} S_{\varepsilon}\end{split}\]

Line search

Line search algorithm find an optimal step length (\(\alpha\)) for the search direction \(\boldsymbol{p}_k\). It is obtained by minimizing the function

\[f(\alpha) = F(\boldsymbol{x}_k + \alpha \boldsymbol{p}_k), \qquad \dfrac{d f(\alpha)}{d \alpha} = \boldsymbol{\nabla} F(\boldsymbol{x}_k + \alpha \boldsymbol{p}_k) \boldsymbol{p}_k\]

enough to satisfy strong Wolfe conditions:

\[\begin{split}F(\boldsymbol{x}_k + \alpha\boldsymbol{p}_k) &\le F(\boldsymbol{x}_k) + c_1 \alpha_k \boldsymbol{\nabla} F(\boldsymbol{x}_k) \boldsymbol{p}_k, \\ \vert\boldsymbol{\nabla} F(\boldsymbol{x}_k + \alpha\boldsymbol{p}_k)\boldsymbol{p}_k\vert &\le c_2\vert\boldsymbol{\nabla} F(\boldsymbol{x}_k)\boldsymbol{p}_k\vert\end{split}\]

Line search algorithm:

Given \(\boldsymbol{x}_k\) and \(\boldsymbol{p}_k\)

1. If \(\alpha = 1\) satisfies strong Wolfe condition, then return \(1\).

2. Set \(\alpha_0 = 0\), \(\alpha_{\text{max}} = 2\) and chose \(\alpha_1\) via Cubic interpolation;

3. \(i \gets 1\);

4. While maximum number of iterations is not achieved:

   1. Compute \(f(\alpha_i) = F(\boldsymbol{x}_k + \alpha_i \boldsymbol{p}_k)\);

   2. If \(f(\alpha_i) > f(0) + c_1 \alpha_i f^{\prime}(0)\) or \(f(\alpha_i) \ge f(\alpha_{i-1})\) and \(i > 1\), then return \(zoom(\alpha_{i-1}, \alpha_i)\);

   3. Compute \(f^{\prime}(\alpha_i) = \boldsymbol{\nabla} F(\boldsymbol{x}_k + \alpha_i \boldsymbol{p}_k) \boldsymbol{p}_k\);

   4. If \(\vert f^{\prime}(\alpha_i)\vert \le -c_2 f^{\prime}(0)\), then return \(\alpha_i\);

   5. If \(f^{\prime}(\alpha_i) \ge 0\), then return \(zoom(\alpha_i, \alpha_{i-1})\);

   6. Choose \(\alpha_{i+1}\) via Cubic interpolation;

   7. \(i \gets i + 1\).

\(zoom\) algorithm:

Given \(\alpha_{lo}\), \(\alpha_{hi}\)

1. Repeat

   1. Interpolate \(\alpha_j\) via Cubic interpolation;

   2. Compute \(f(\alpha_j) = F(\boldsymbol{x}_k + \alpha_j \boldsymbol{p}_k)\);

   3. Check that value of the function sufficiently decreases.

   4. If \(f(\alpha_j) > f(0) + c_1 \alpha_j f^{\prime}(0)\) or \(f(\alpha_j) \ge f(\alpha_{lo})\), then \(\alpha_{hi} \gets \alpha_j\)

   5. Else

      1. If \(\vert f^{\prime}(\alpha_j)\vert \le -c_2 f^{\prime}(0)\), then return \(\alpha_j\);

      2. If \(f^{\prime}(\alpha_j)(\alpha_{hi} - \alpha_{lo}) \ge 0\), then \(\alpha_{hi} \gets \alpha_{lo}\);

      3. \(\alpha_{lo} \gets \alpha_j\).

Cubic interpolation

Given \(\alpha_l\), \(\alpha_h\) and \(f(\alpha_l)\), \(f(\alpha_h)\) and \(f^{\prime}(\alpha_l)\), \(f^{\prime}(\alpha_h)\) compute new \(\alpha_m\) as

\[\begin{split}\alpha_{min} &= \alpha_h - (\alpha_h - \alpha_l) \dfrac{f^{\prime}(\alpha_h) + d_2 - d_1}{f^{\prime}(\alpha_h) - f^{\prime}(\alpha_l) + 2d_2} \\ d_1 &= f^{\prime}(\alpha_l) + f^{\prime}(\alpha_h) - 3 \dfrac{f(\alpha_l) - f(\alpha_h)}{\alpha_l - \alpha_h} \\ d_2 &= \text{sign}(\alpha_h - \alpha_l) \sqrt{d_1^2 - f^{\prime}(\alpha_l)f^{\prime}(\alpha_h)}\end{split}\]

If \(d_1^2 - f^{\prime}(\alpha_l)f^{\prime}(\alpha_h) < 0\), then \(\alpha_{min} = \alpha_l\) if \(f(\alpha_l) \le f(\alpha_h)\), otherwise \(\alpha_{min} = \alpha_h\).