Minimization of energy#
References#
On this page we recall the approach and main formulas from [1] and [2]. The minimization of the magnetic ground state in Magnopy is an implementation of the method described in [2].
Minimization of the energy function (\(E^{(0)}\) or \(E^{(0)} + E^{corr}\)) can be formulated as a problem of minimizing that function over the \(M\) vectors of the spin directions \(\boldsymbol{z}_{\alpha}, \alpha = 1, ..., 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 directional 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:
where \(\boldsymbol{A}_{\alpha}\) are skew-symmetric matrices parametrized by three real numbers as
and
Then energy function can be rewritten as
In other words, as a function of vector \(\boldsymbol{x}\) from the vector space \(\mathbb{R}^{3I}\):
Then, energy of the system is minimized with the BFGS algorithm [1].
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm#
Formula for the inverse Hessian update:
Given Initial guess \(\boldsymbol{x}_0\) and Initial approximation of the inverse hessian matrix \(\boldsymbol{H}_0\),
\(k \gets 0\)
While convergence is not achieved:
Compute the gradient of the function \(\boldsymbol{\nabla} F(\boldsymbol{x}_k)\);
Compute the search direction \(\boldsymbol{p}_k = -\boldsymbol{H}_k \boldsymbol{\nabla} F(\boldsymbol{x}_k)\);
Compute length of the step \(\alpha_k\) via Line search;
Set \(\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \alpha_k \boldsymbol{p}_k\) and compute gradient \(\boldsymbol{\nabla} F(\boldsymbol{x}_{k+1})\);
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)\);
Update the hessian matrix \(\boldsymbol{H}_{k+1}\) by the BFGS formula;
\(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 of each directional 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
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 [2]
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 classical energy is computed analytically
where \(\tilde{J}_{\alpha}^i\) is a single-spin renormalized parameter defined by equation S.68 of Supplementary material of paper about Magnopy.
The gradient of the correction energy \(E^{corr}\) is computed numerically by the two point formula.
Line search#
Line search algorithm defines an optimal step length (\(\alpha\)) for the search direction \(\boldsymbol{p}_k\). It is obtained by minimizing the function
enough to satisfy strong Wolfe conditions:
Line search algorithm:
Given \(\boldsymbol{x}_k\) and \(\boldsymbol{p}_k\)
If \(\alpha = 1\) satisfies strong Wolfe condition, then return \(1\).
Set \(\alpha_0 = 0\), \(\alpha_{\text{max}} = 2\) and chose \(\alpha_1\) via Cubic interpolation;
\(i \gets 1\);
While maximum number of iterations is not achieved:
Compute \(f(\alpha_i) = F(\boldsymbol{x}_k + \alpha_i \boldsymbol{p}_k)\);
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)\);
Compute \(f^{\prime}(\alpha_i) = \boldsymbol{\nabla} F(\boldsymbol{x}_k + \alpha_i \boldsymbol{p}_k) \boldsymbol{p}_k\);
If \(\vert f^{\prime}(\alpha_i)\vert \le -c_2 f^{\prime}(0)\), then return \(\alpha_i\);
If \(f^{\prime}(\alpha_i) \ge 0\), then return \(zoom(\alpha_i, \alpha_{i-1})\);
Choose \(\alpha_{i+1}\) via Cubic interpolation;
\(i \gets i + 1\).
\(zoom\) algorithm:
Given \(\alpha_{lo}\), \(\alpha_{hi}\)
Repeat
Interpolate \(\alpha_j\) via Cubic interpolation;
Compute \(f(\alpha_j) = F(\boldsymbol{x}_k + \alpha_j \boldsymbol{p}_k)\);
Check that value of the function sufficiently decreases.
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\)
Else
If \(\vert f^{\prime}(\alpha_j)\vert \le -c_2 f^{\prime}(0)\), then return \(\alpha_j\);
If \(f^{\prime}(\alpha_j)(\alpha_{hi} - \alpha_{lo}) \ge 0\), then \(\alpha_{hi} \gets \alpha_{lo}\);
\(\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
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\).