Method for Estimating Spin-Spin Interactions from Magnetization Curves

TL;DR

This paper introduces a machine learning approach using Bayesian inference and Monte Carlo methods to accurately estimate spin-spin interactions in magnetic systems from noisy magnetization data.

Contribution

It presents a novel method combining Bayesian inference, MCMC, and $l_1$ regularization to determine relevant spin interactions from magnetization curves.

Findings

High accuracy in estimating spin-spin interactions.

Effective selection of relevant interaction terms.

Validated with synthetic data.

Abstract

We develop a method to estimate the spin-spin interactions in the Hamiltonian from the observed magnetization curve by machine learning based on Bayesian inference. In our method, plausible spin-spin interactions are determined by maximizing the posterior distribution, which is the conditional probability of the spin-spin interactions in the Hamiltonian for a given magnetization curve with observation noise. The conditional probability is obtained by the Markov-chain Monte Carlo simulations combined with an exchange Monte Carlo method. The efficiency of our method is tested using synthetic magnetization curve data, and the results show that spin-spin interactions are estimated with a high accuracy. In particular, the relevant terms of the spin-spin interactions are successfully selected from the redundant interaction candidates by the $l_1$ regularization in the prior distribution.

Tables

Table 1. Table 1: Number of each type of interaction J k subscript 𝐽 𝑘 J_{k} per tetrahedron, distance between spins connected by each interaction J k subscript 𝐽 𝑘 J_{k} within the nearest-neighbor lattice spacing, spin-spin interactions used in the input magnetization curve, and spin-spin interactions estimated using our estimation method.

	$J_1$	$J_2$	$J_3$	$J_4$	$J_5$	$J_6$	$J_7$	$b$
Number of interactions $\left(n_k\right)$	$2$	$2$	$1$	$1$	$2$	$1$	$1$	$-$
Distance between spins $\left(r_k\right)$	$1$	$1$	$1$	$1$	$\sqrt{3}$	$2$	$2$	$-$
Input spin-spin interactions (${𝐱}^{\prime }$)	$1.000$	$4.000$	$5.000$	$6.000$	$0.000$	$0.000$	$0.000$	$0.100$
Estimated spin-spin interactions (${𝐱}^{\ast }$)	$1.074$	$3.850$	$5.012$	$6.051$	$0.011$	$-0.051$	$0.002$	$0.102$

Equations

$$m(H,\mathbf{x})=\left|\frac{1}{N|\mathbf{s}|}\sum _{i=1}^N⟨{\mathbf{s}}_i{⟩}_{H,\mathbf{x}}\right|,$$

$$\displaystyle\langle\mathbf{s}_{i}\rangle_{H,\mathbf{x}}=\frac{\displaystyle{{\rm Tr}\ \mathbf{s}_{i}\exp\left[-\left(\mathcal{H}(\{\mathbf{s}_{i}\};\mathbf{x})-H\sum_{i=1}^{N}s_{i}^{z}\right)/T\right]}}{\displaystyle{{\rm Tr}\exp\left[-\left(\mathcal{H}(\{\mathbf{s}_{i}\};\mathbf{x})-H\sum_{i=1}^{N}s_{i}^{z}\right)/T\right]}},\ \ \ \$$

$$\displaystyle P(m(H,\mathbf{x})|\mathbf{x})=\delta\left(m(H,\mathbf{x})-\left|\frac{1}{N|\mathbf{s}|}\sum_{i=1}^{N}\langle\mathbf{s}_{i}\rangle_{H,\mathbf{x}}\right|\right),\ \$$

$$m^{\mathrm{ex}}(H)=m(H,\mathbf{x})+ϵ.$$

$$P(ϵ)\propto \exp \left(-\frac{{ϵ}^2}{2{\sigma }^2}\right).$$

$$P(m^{\mathrm{ex}}(H)|m(H,\mathbf{x}))$$

$$\propto \exp \left[-\frac{1}{2{\sigma }^2}{\left(m^{\mathrm{ex}}(H)-m(H,\mathbf{x})\right)}^2\right].$$

$$P(m^{\mathrm{ex}}(H)|\mathbf{x})$$

$$\propto \int dm(H,\mathbf{x})P(m^{\mathrm{ex}}(H)|m(H,\mathbf{x}))P(m(H,\mathbf{x})|\mathbf{x})$$

$$\displaystyle\ \ \ \ \ \propto\exp\left[-\frac{1}{2\sigma^{2}}\left(m^{\rm ex}(H)-\left|\frac{1}{N|\mathbf{s}|}\sum_{i=1}^{N}\langle\mathbf{s}_{i}\rangle_{H,\mathbf{x}}\right|\right)^{2}\right].\ \ \ \ \$$

$$\displaystyle P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})=\frac{P(\{m^{\rm ex}(H_{l})\}_{l\in D}|\mathbf{x})P(\mathbf{x})}{P(\{m^{\rm ex}(H_{l})\}_{l\in D})},\ \$$

$$P(\mathbf{x})\propto \exp \left(-\lambda \sum _{k=1}^K|x_k|\right),$$

$$P(\{m^{\mathrm{ex}}(H_l){\}}_{l\in D}|\mathbf{x})=\prod _{l=1}^{L}P(m^{\mathrm{ex}}(H_l)|\mathbf{x}).$$

$$P(\mathbf{x}|\{m^{\mathrm{ex}}(H_l){\}}_{l\in D})$$

$$\propto \prod _{l=1}^{L}P(m^{\mathrm{ex}}(H_l)|\mathbf{x})P(\mathbf{x})$$

$$\displaystyle\ \propto\exp\Biggl{[}-\frac{1}{2\sigma^{2}}\sum_{l=1}^{L}\left(m^{\rm ex}(H_{l})-\left|\frac{1}{N|\mathbf{s}|}\sum_{i=1}^{N}\langle\mathbf{s}_{i}\rangle_{H_{l},\mathbf{x}}\right|\right)^{2}$$

$$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -\lambda\sum_{k=1}^{K}|x_{k}|\Biggr{]}.\ \ \ \ \ \$$

$$E(\mathbf{x}|\sigma ,\lambda ,K)$$

$$:=\frac{1}{2{\sigma }^2}\sum _{l=1}^L{\left(m^{\mathrm{ex}}(H_l)-\left|\frac{1}{N|\mathbf{s}|}\sum _{i=1}^N⟨{\mathbf{s}}_i{⟩}_{H_l,\mathbf{x}}\right|\right)}^2$$

$$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +\lambda\sum_{k=1}^{K}|x_{k}|,\ \$$

$$P(\mathbf{x}|\{m^{\mathrm{ex}}(H_l){\}}_{l\in D})\propto \exp \left[-\frac{1}{T_{\mathrm{R}}}E(\mathbf{x}|\sigma ,\lambda ,K)\right],$$

$$w({\mathbf{x}}^{\prime }|\mathbf{x})=\min \left\{1,\exp \left[-\frac{1}{T_{\mathrm{R}}}\Delta E({\mathbf{x}}^{\prime },\mathbf{x})\right]\right\},$$

$$\Delta E({\mathbf{x}}^{\prime },\mathbf{x}):=E({\mathbf{x}}^{\prime }|\sigma ,\lambda ,K)-E(\mathbf{x}|\sigma ,\lambda ,K).$$

$$w_{\mathrm{ex}}({\mathbf{x}}^j|{\mathbf{x}}^i)=\min \left\{1,\exp \left[\left(\frac{1}{T_{\mathrm{R}}^j}-\frac{1}{T_{\mathrm{R}}^i}\right)\Delta E({\mathbf{x}}^j,{\mathbf{x}}^i)\right]\right\},$$

$$\Delta E({\mathbf{x}}^j,{\mathbf{x}}^i):=E({\mathbf{x}}^j|\sigma ,\lambda ,K)-E({\mathbf{x}}^i|\sigma ,\lambda ,K),$$

$${\Delta }^{(s)}(\lambda )$$

$$:=\frac{1}{L_{\mathrm{te}}}\sum _{l_{\mathrm{te}}\in D_s}{\left(m^{\mathrm{ex}}(H_{l_{\mathrm{te}}})-\left|\frac{1}{N|\mathbf{s}|}\sum _{i=1}^N⟨{\mathbf{s}}_i{⟩}_{H_{l_{\mathrm{te}}},{\mathbf{x}}_s^{\ast }}\right|\right)}^2.$$

$$\mathcal{H}=\sum _{i,j}{J}_{ij}\left[{\mathbf{s}}_i\cdot {\mathbf{s}}_j-b_{ij}({\mathbf{s}}_i\cdot {\mathbf{s}}_j{)}^2\right]-H\sum _i{s}_i^z.$$

$$\mathrm{Type}\mathrm{I}:$$

$$\mathrm{Type}\mathrm{II}:$$

$$\mathrm{Type}\mathrm{III}:$$

Full text

I Introduction

The importance of data-driven techniques using machine learningMitchell-1997 ; Bishop-2006 is recognized in both academic and industrial fields. Machine learning is generally a statistical tool used to extract the inherent structure from a finite set of observed data. In the condensed matter physics, machine learning techniques have been used such as for interpolation of density-functional theory (DFT)Snyder-2012 ; Bartok-2013 ; Zhou-2014 ; Li-2015 and dynamical-mean field theoryArsenault-2014 calculations, and for model selection of strongly correlated systemsTakenaka-2014 and the Ginzburg-Landau equationTaverniers-2015 . In this paper, we apply Bayesian inference from the field of machine learning Toussaint-2011 ; DAgostini-2003 ; Dose-2003 to estimate the spin-spin interactions in the Hamiltonian given an observed magnetization curve as an input (Fig. 1).

Magnetization is the average of the magnetic moments normalized by the amplitude of the magnetic moment, which is induced in a magnetic material by a magnetic field and is a fundamental physical quantity in both experiments and theories of magnetism. The magnetic-field dependence of the magnetization is called the magnetization curveNagamiya-1962 ; Kanamori-1966 , and has been measured in many magnetic materials using devices such as a superconducting quantum interference device (SQUID) in the laboratory. In general, the magnetization curve for a fixed temperature increases monotonically as a nonlinear function of the magnetic field. Depending on the magnetic material, various characteristic features, such as plateausHida-1994 ; Oshikawa-1997 ; Totsuka-1998 ; Kitazawa-1999 ; Onizuka-2000 ; Kikuchi-2005 ; Kikuchi-2005 ; Ishii-2011 , rampsBindilatti-1996 ; Gratens-2004 ; Nakano-2010 , and sharp increasesKatsumata-1989 , may be observed in the magnetization curve. Consequently, the magnetization curve is an important quantity for understanding the magnetic properties of individual magnetic materials.

To understand the essence of the properties in magnetic materials, an effective Hamiltonian is often derived in materials science. Many methods have been proposed to determine Hamiltonians, and they can be divided roughly into two groups: one where the model parameter values in the Hamiltonian are calculated by ab initio electronic structure calculations when basic data of the magnetic materials is inputtedMazin-2007 ; Nakamura-2008 ; Pruneda-2008 ; Hirayama-2013 , and the other where the model parameter values are evaluated by fitting the physical quantities observed in a target magnetic materialTakubo-2007 ; Tamura-2011 ; Matsumoto-2012 ; DiStasio-2013 ; Marcotte-2013 .

The aim of this paper is to establish a new Hamiltonian estimation method, which can be classified as the latter group where the model parameters in the Hamiltonian of a target material are predicted from experimental data (Fig. 1). In particular, we propose a method, which selects relevant terms for the model parameters from redundant model parameter candidates by machine learning based on the $l_{1}$ regularization. The framework of this method can be adopted for any measurement data, but herein we focus on the case where an observed magnetization curve is used as the input data to estimate the spin-spin interactions.

The rest of the paper is organized as follows. In Sec. II, we introduce a method to estimate the model parameters in the Hamiltonian from the observed magnetization curve. We first present the forward modeling for an observed magnetization curve. By using Bayesian inference, we construct the posterior distribution, that is, the conditional probability of the model parameters in the Hamiltonian given an observed magnetization curve (the Bayes modeling). In this case, the plausible model parameters maximize the posterior distribution. Next, we explain the simulation methods, which employ the Markov-chain Monte Carlo (MCMC) method and exchange Monte Carlo method, to analyze the posterior distribution. In the last part of Sec. II, we describe cross validation to avoid over-fitting. Cross validation makes it possible to calculate prediction errors and allows a hyperparameter in the prior distribution to be determined. To validate our estimation method, we determine the spin-spin interactions from a synthetic magnetization curve in Sec. III. We prepare an input observed magnetization curve using a theoretical model with fixed spin-spin interactions to consider the case where the solution of the estimation problem is known in advance. We show that the spin-spin interactions can be accurately estimated from the observed magnetization curve using our method. Section V presents the discussion and summary.

II Method for estimating model parameters

This section presents our estimation method for model parameters in the Hamiltonian determined from the observed magnetization curve. Using Bayesian inference, the conditional probability of the spin-spin interactions are obtained for a given observed magnetization curve. Furthermore, the simulation methods used to analyze the conditional probability and the cross validation to avoid over-fitting are introduced.

II.1 Bayesian inference for the magnetization curve

II.1.1 Forward modeling

Let us assume a class of the Hamiltonian $\mathcal{H}(\{\mathbf{s}_{i}\};\mathbf{x})$ for a magnetic system with a set of model parameters denoted as $\mathbf{x}=(x_{1},...,x_{K})$ where $K$ is the number of model parameters. Here, $\mathbf{s}_{i}$ is the spin variable at site $i$ with the absolute value $|\mathbf{s}|$, and $\{\mathbf{s}_{i}\}$ denotes a set of spin configurations. In statistical physics, the magnetization under a magnetic field $H$ at a temperature $T$ for a given Hamiltonian $\mathcal{H}(\{\mathbf{s}_{i}\};\mathbf{x})$ is calculated as

[TABLE]

where $s_{i}^{z}$ is the $z$-component of $\mathbf{s}_{i}$ and $N$ is the number of spins. Throughout the paper, the $g$-factor, the Bohr magneton, and the Boltzmann constant are set to unity. In the forward modeling, the conditional probability of the observed magnetization in the given model parameters is examined. Because the magnetization $m(H,\mathbf{x})$ is uniquely obtained in the thermodynamic limit by statistical physics, the conditional probability of the magnetization $m(H,\mathbf{x})$ given $\mathbf{x}$ is expressed as

[TABLE]

where $\delta(\cdots)$ is the Dirac delta function and $P(A|B)$ expresses the conditional probability of $A$ given $B$.

The measurements in an experiment always have some uncertainty. Taking this uncertainty into account as the observation noise $\epsilon$, the observed magnetization in experiments $m^{\rm ex}(H)$ is expressed as

[TABLE]

At this point, we assume that the observation noise follows a Gaussian distribution with a mean of zero and a standard deviation of $\sigma$:

[TABLE]

The conditional probability of $m^{\rm ex}(H)$ for a given $m(H,\mathbf{x})$ is then

[TABLE]

Thus, in the forward modeling, the conditional probability of the observed magnetization $m^{\rm ex}(H)$ given the model parameters $\mathbf{x}$ is written as

[TABLE]

By using this formula, the probabilistic prediction of the magnetization observed in the experiments is obtained as the conditional probability for the given model parameters. This can be compared with the experimental data. Figure 1 shows the flow of the forward modeling.

II.1.2 Bayes modeling

The Bayes modeling provides a general framework to estimate the model parameters $\mathbf{x}$ from the observed magnetizations. In this work, we consider the case where the magnetization under various magnetic fields is observed experimentally; that is, the observed magnetization curve is expressed as $\{m^{\rm ex}(H_{l})\}_{l\in D}$ for a series of magnetic fields $\{H_{l}\}_{l\in D}$ where the data are indexed by $D=(1,...,L)$ and $L$ is the number of observed magnetizations. By using Bayes’ theorem, the posterior distribution, which is the conditional probability of the model parameters $\mathbf{x}$ given the observed magnetization curve $\{m^{\rm ex}(H_{l})\}_{l\in D}$, is expressed as

[TABLE]

where $P(\mathbf{x})$ is the prior distribution of $\mathbf{x}$. On the right-hand side, the denominator is independent of the model parameters. In contrast, the prior distribution $P(\mathbf{x})$ explicitly contributes to the estimation in the posterior distribution. For example, in the case without a priori knowledge of the model parameters, a uniform distribution of $P(\mathbf{x})$ may be used. In this work, we assume that the number of important model parameters is small, indicating that the relevant parameters are sparse in the model Hamiltonian. In this situation, the prior distribution $P(\mathbf{x})$ based on the $l_{1}$ regularization is frequently represented as

[TABLE]

where $\lambda$ is called the hyperparameter given before the analysis. The inference scheme using the prior distribution is called LASSO (least absolute shrinkage and selection operator)Tibshirani-1996 , and the value of $\lambda$ controls the strength of the regularization.

In addition, we assume that $m^{\rm ex}(H_{l})$ is independently obtained from different magnetic fields. That is, the previous measurement in a magnetization process does not affect the present measurement, which gives

[TABLE]

Combined with the prior distribution of Eq. (10), the posterior distribution of $\mathbf{x}$ given $\{m^{\rm ex}(H_{l})\}_{l\in D}$ based on Bayes’ theorem is expressed as

[TABLE]

From the view point of the maximum a posterior (MAP) estimation, the plausible model parameters $\mathbf{x}^{*}$ are obtained by the maximizer of Eq. (13). In other words, we search the model parameters so that Eq. (13) is maximized in the Bayes modeling where $\sigma$, $\lambda$, and $K$ are the previously fixed hyperparameters. Although the posterior distribution is constructed using the magnetization curve as the input data, it should be emphasized that the framework of our estimation method can be adopted for any measurement data.

II.2 Simulation method

In this study, we use the MCMC method and the exchange Monte Carlo method to analyze the posterior distribution $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})$. This combination significantly contributes to finding the global maximum of the posterior distribution in systems where many local maxima exist. When the energy function $E(\mathbf{x}|\sigma,\lambda,K)$ is defined as

[TABLE]

the posterior distribution can be written as

[TABLE]

where $T_{\rm R}=1$ in the case of Eq. (13). This distribution formally has the same form as the Boltzmann distribution by regarding $E(\mathbf{x}|\sigma,\lambda,K)$ and $T_{\rm R}$ as the energy of the system and the virtual temperature, respectively. Thus, a statistical mechanical approach is a promising for sampling from the distributionLandau-2009 . It should be noted that the dynamical variables in this scheme are the model parameters in the Hamiltonian, which is unusual in the study of magnetism as well as in statistical physics.

In our MCMC simulations, we employ the Metropolis-type transition probability from $\mathbf{x}$ to $\mathbf{x}^{\prime}$:

[TABLE]

By sampling, we investigate the properties of $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})$ and search for the model parameters $\mathbf{x}^{*}$ that maximize $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})$.

In addition, we use the exchange Monte Carlo methodHukushima-1996 to enhance the sampling efficiency and to find the global maximum of $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})$ in a system where many local maxima exist. In the method, several MCMC simulations using Eq. (16) with different virtual temperatures are performed in parallel, and the model parameters $\mathbf{x}$ are exchanged between the different MCMC simulations with an exchange probability expressed as

[TABLE]

where $\mathbf{x}^{i}$ and $T_{\rm R}^{i}$ indicate the model parameters and the virtual temperature in the $i$-th MCMC simulation, respectively. By using this technique, we can find the global maximum of $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in D})$ and obtain the plausible model parameters $\mathbf{x}^{*}$ with a high probability. Without a loss in generality, the units of the parameter $\lambda$ and the virtual temperature $T_{\rm R}$ can be taken as $\sigma^{2}$. Thus, we set $\sigma=1$ throughout the paper.

II.3 Cross validation

In general, the observed magnetization curve is well fitted when the number of model parameters is sufficiently large. Often the estimated model in machine learning cannot properly predict unknown data not used in the fitting, even when the observed data are well fitted by the estimated model. This is known as the over-fitting problem.

To determine the hyperparameter $\lambda$ while preventing over-fitting, we perform cross validation in which the hyperparameter is chosen to minimize the prediction error. Consider the case where the value of $K$ is fixed. In cross validation, the dataset $D$ is divided into $S$ groups as $D_{s}=((s-1)L/S+1,(s-1)L/S+2,...,sL/S)$ with $s=1,...,S$. One of the $S$ groups is regarded as the test data, while the remaining $S-1$ groups are used as the training data. Figure 2 shows an example for $S=4$. The number of pieces of test data and that of the training data are given by $L_{\rm te}:=N/S$ and $L_{\rm tr}:=N(S-1)/S$, respectively. For each data subset $G_{s}=D\setminus D_{s}$, the plausible model parameters $\mathbf{x}^{*}_{s}$ are determined by the MAP estimation according to $\lambda$. The mean-square deviation between the test data and the estimated magnetization curve when using $\mathbf{x}^{*}_{s}$ for various $\lambda$ is evaluated as

[TABLE]

In the cross validation, $\Delta^{(s)}(\lambda)$ represents the prediction error when the test data $D_{s}$ are treated as unknown data. Using the average of $\Delta^{(s)}(\lambda)$ over different data sets, we evaluate the optimal value of $\lambda^{*}$ where the averaged prediction error takes a minimum as a function of $\lambda$. Furthermore, the plausible model parameters $\mathbf{x}^{*}$ are obtained such that $\Delta^{(s)}(\lambda)$ is the smallest among $S$ selections for the test data for $\lambda^{*}$. It should be noted that $\lambda^{*}$ depends on the input magnetization curve by performing cross validation.

III Demonstration of estimation of spin-spin interactions

In this section, we confirm that the estimation method proposed in Sec. II can determine the spin-spin interactions in the Hamiltonian. The inputted observed magnetization curve was a zero-temperature magnetization curve calculated from a theoretical model. The results confirm that our estimation method is highly accurate and the relevant terms of the spin-spin interactions are successfully selected from the redundant interaction candidates by the $l_{1}$ regularization.

III.1 Model and observed magnetization curve

We examine the efficiency of our proposed method by synthetic data for the magnetization curve. To generate the input data, we used the classical Heisenberg model with bilinear and biquadratic interactions, which has the following Hamiltonian:

[TABLE]

Here, $b_{ij}=bJ_{ij}$, and $\mathbf{s}_{i}=(s_{i}^{x},s_{i}^{y},s_{i}^{z})$ denotes the three component vector spin with a length $|\mathbf{s}|=1/2$. The value of $b$ represents the amplitude of the biquadratic interactions. Depending on the lattice and sets of spin-spin interactions, this model exhibits fruitful magnetization curvesPenc-2004 ; Motome-2006 ; Shannon-2006 . Here, the inputted observed magnetization curve is given by the model on a tetrahedral chain at a zero temperature. The lattice structure is shown in Fig. 3. For simplicity, we consider seven types of spin-spin interactions $J_{k}\ (k=1,...,7)$. Thus, the model parameter can be expressed by $\mathbf{x}=(J_{1},J_{2},J_{3},J_{4},J_{5},J_{6},J_{7},b)$ with the dimension $K=8$.

Figure 4 (a) is the zero-temperature magnetization process $\{m(H_{l},\mathbf{x}^{\prime})\}_{l\in D}$ with the number of data points $L=160$, which is obtained by $\mathbf{x}^{\prime}=(1,4,5,6,0,0,0,0.1)$. The magnetization curve exhibits a huge 2/3 magnetization plateau. Furthermore, we assume that the observed magnetization curve $\{m^{\rm ex}(H_{l})\}_{l\in D}$ can be represented by adding Gaussian noise with a mean of zero and a standard deviation of $\sigma=0.004$. Figure 4 (b) shows the synthesized magnetization curve $\{m^{\rm ex}(H_{l})\}_{l\in D}$.

For the magnetization curve shown in Fig. 4 (b), we estimate model parameters $\mathbf{x}=(J_{1},J_{2},J_{3},J_{4},J_{5},J_{6},J_{7},b)$ using our proposed method. For the cross validation, we randomly divided the observed magnetization curve into four sets of data points $(S=4)$. Each set of training data contains 120 data points $(L_{\rm te}=120)$ in the magnetization curve.

III.2 Estimation of spin-spin interactions

This section explains the details of the numerical simulations to evaluate the posterior distribution $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in G_{s}})$ with $s=1,...,4$. Our numerical method consists of double loop calculations, which include a calculation of the magnetization curve for a given set of model parameters and the sampling of the model parameters. We used the Hamiltonian defined by Eq. (21) with periodic boundary conditions where the total number of spins was six as a theoretical model to obtain $\langle\mathbf{s}_{i}\rangle_{H_{l},\mathbf{x}}$. The zero-temperature magnetization curve was obtained by the steepest descent method. This is the inner-loop calculation which was performed in a few seconds. Monte Carlo simulations were employed to sample the model parameters $\mathbf{x}$ in the outer loop. Typically, the number of the Monte Carlo steps to update the model parameters was $10^{4}$. Furthermore, we prepared $20$ replicas with virtual temperatures $T_{\rm R}$ between $0.001$ and $10$ for the exchange Monte-Carlo method.

We considered three types of prior distributions for the model parameters $\mathbf{x}$ based on the $l_{1}$ regularization:

[TABLE]

where $n_{k}$ is the number of each type of interaction $J_{k}$ per tetrahedron, and $r_{k}$ is the distance between the spins connected by interaction $J_{k}$. Table I summarizes the values used in this demonstration. Type I is the general form of the $l_{1}$ regularization. Type II considers the number of interactions per tetrahedron. In the theoretical model, the number of interactions in the lattice depends on $J_{k}$. Thus, we constructed the $l_{1}$ regularization from the viewpoint of the sum of the absolute values of the spin-spin interactions (i.e., the number of interactions per tetrahedron was introduced in the prior distribution). Type III considers both the number and the distance of the interactions. In general, the amplitude of the spin-spin interactions decreases as the distance between spins increases. Therefore, we assumed that the effect of the long-range interactions on the physical properties was smaller than that of the short-range interactions. To add a penalty to the long-range interactions, we introduced the distance between spins in the prior distribution in Type III.

Figure 5 shows the $\lambda$ dependence of the estimated spin-spin interactions that maximize $P(\mathbf{x}|\{m^{\rm ex}(H_{l})\}_{l\in G_{s}})$ with $s=1,...,4$. The results are the averages over the spin-spin interactions estimated by the four sets of training data $\mathbf{x}_{s}^{*}$ with $s=1,...,4$, and the error bars calculated from the standard deviation. Large error bars indicate that the estimated spin-spin interactions differ across the four sets of training data. The dashed lines in Fig. 5 show the values of the spin-spin interactions used in the inputted magnetization curve, that is, the solution of the estimation problem. Between $\lambda=0.002$ and $0.008$ for Type III, which is shown in the shaded area in Fig. 5 (c), the estimated spin-spin interactions are roughly the same as the values used in the inputted magnetization curve. However, these figures do not conclusively show which estimated spin-spin interactions are plausible for the inputted magnetization curve.

Using cross validation, we calculated the prediction error $\Delta_{\rm av}(\lambda)$, which is the average of $\Delta^{(s)}(\lambda)$ with $s=1,...,4$ defined in Eq. (20). Figure 6 shows the $\lambda$ dependence of $\Delta_{\rm av}(\lambda)$ for each type of prior distribution. The minimum value of $\Delta_{\rm av}(\lambda)$ occurs at $\lambda=0.004$ in Type III, which produces the highest prediction accuracy. We conclude that the estimated spin-spin interactions for $\lambda=0.004$ in Type III are the most likely spin-spin interactions for the inputted magnetization curve. Table I, which shows the estimated spin-spin interactions for $\lambda=0.004$ in Type III where $\Delta^{(s)}(\lambda)$ is the smallest among $s=1,...,4$, confirms that the estimated spin-spin interactions are almost the same as the values used in the inputted magnetization curve. This demonstrates that our proposed estimation method can estimate the spin-spin interactions with a high accuracy from the observed magnetization curve. Furthermore, the estimated values of $J_{5}$, $J_{6}$, and $J_{7}$ are very close to zero, indicating that the important spin-spin interactions can be selected among the many candidate interactions in the Hamiltonian. In addition, Fig. 6 shows that over-fitting occurs in the case of $\lambda=0$ when using non-zero values for the additional interaction terms. Consequently, over-fitting can be avoided by introducing a regularization term based on the $l_{1}$ regularization. Finally, Fig. 7 compares the estimated magnetization curve for $\lambda=0.004$ in Type III with the inputted magnetization curve.

The model parameters in the Hamiltonian can be estimated with a high accuracy from the observed magnetization curve using the estimation method introduced in Sec. II. In addition, we confirm that the model parameters are well estimated for certain input magnetization curves with the Hamiltonian defined in Eq. (21).

IV Discussion and summary

We developed a method for estimating the model parameters in the Hamiltonian from the observed magnetization curve. Bayesian inference is initially used to determine the posterior distribution (i.e., the conditional probability of model parameters in the Hamiltonian for a given observed magnetization curve). We adopted the $l_{1}$ regularization as the prior distribution of the model parameters. In this scheme, the plausible model parameters are determined by maximizing the posterior distribution. To analyze the posterior distribution, we adopted the Markov-chain Monte Carlo method and the exchange Monte Carlo method because this combination significantly contributes to finding the global maximum of a posterior distribution where many local maxima exist. Cross validation is introduced to avoid over-fitting and to determine the hyperparameter in the prior distribution. Note that the framework of our estimation method can be used not only for magnetization curves but also for any measured data as the input data.

The effectiveness of our method is demonstrated by estimating the spin-spin interactions in the Hamiltonian from the observed magnetization curve. We used a magnetization curve synthesized by adding Gaussian noise prepared by a theoretical model with fixed spin-spin interactions. Specifically, we employed a classical Heisenberg model with bilinear and biquadratic interactions on a tetrahedral chain to consider the case where the solution of the estimation problem is known in advance. Furthermore, we introduced three types of prior distributions of the spin-spin interactions based on the $l_{1}$ regularization. These prior distributions include the hyperparameter $\lambda$, which controls the strength of the regularization. The spin-spin interactions are sampled from the posterior distribution for a given $\lambda$ by numerical simulations based on the Monte Carlo method. We also determined the value of $\lambda$ and the type of the prior distribution by minimizing the prediction error by cross validation. Spin-spin interactions consistent with the input data can be estimated. It is noted that our estimation method successfully selects the relevant spin-spin interactions from the many redundant interactions in the Hamiltonian. The prior distribution taking into account both the number and the distance of the interactions plays a relevant role in this study, confirming that our method can estimate the spin-spin interactions from an observed magnetization curve with a high accuracy.

In general, it is difficult to obtain physical quantities measured by indirect measurements such as magnetic specific heat and magnetic entropy. However, our estimation method should predict these quantities once the Hamiltonian for a target magnetic material is obtained. Furthermore, we can directly obtain the predicted spin-snapshots as well as predict the magnetic structure and structure factor. Thus, the information from the estimated Hamiltonian should provide preliminary information before using large facilities such as those for neutron diffraction measurements.

Our final objective is to assist with measurements by providing additional estimation information. However, this research shows only the case of inputting a synthetic magnetization curve obtained from a classical theoretical model, including bilinear and biquadratic interactions. To achieve our broader objectives, we must address two issues. First, we need to establish an estimation method that can be used for many types of magnetic materials described by classical or quantum Hamiltonians with many model parameters such as various types of spin-spin interactions, spin-lattice interactions, and spin anisotropy. In this case, a prior distribution, which estimates the correct model parameters, will depend on the types of target model parameters and the target material. Thus, it is important to develop a method that provides a useful prior distribution. Second, we need to improve our proposed estimation method to produce a more robust Hamiltonian that can simultaneously explain various physical properties of a target material. To achieve this, the input data in the estimation method should include multiple data sets such as some magnetization curves at different temperatures and the magnetic field dependences of the magnetic specific heat. These are the future directions of our research and will be discussed elsewhereTamura-2016 .

Finally, we emphasize that our proposed estimation method is a kind of physics-based machine learning technique because it employs a Hamiltonian based on physical laws. Our proposed method should form the basis for a new applications of machine learning in physics and interdisciplinary physics. In particular, in materials informaticsRajan-2005 ; Jain-2013 , which is an interdisciplinary field between materials science based on the physical law and informatics, our proposed estimation method should be a useful tool to establish important concepts.

Acknowledgement

We thank Masashi Hase, Haruhiko Kuroe, Hideaki Kitazawa, Noriki Terada, Shu Tanaka, and Masato Okada for the useful discussions. R. T. was partially supported by a Grant-in-Aid for Scientific Research (C) (Grant No. 25420698). K. H. was partially supported by a Grant-in-Aid for Scientific Research from JSPS, Japan (Grant No. 25120010 and 25610102). The computations in the present work were performed on Numerical Materials Simulator at NIMS, and the supercomputer at Supercomputer Center, Institute for Solid State Physics, The University of Tokyo. This work was conducted as part of the “Materials Research by Information Integration” Initiative of the Support Program for Starting Up Innovation Hub, Japan Science and Technology Agency.