# GNisi: A graph network for reconstructing Ising models from multivariate
binarized data

Emma Slade
GSK.ai, London,
N1C 4AG, UK
<EMAIL_ADDRESS>
Sonya Kiselgof
GSK.ai, 25 Basel St,
Petah Tikvah, Israel
<EMAIL_ADDRESS>
Lena Granovsky
GSK.ai, 25 Basel St,
Petah Tikvah, Israel
<EMAIL_ADDRESS>
Jeremy L. England
School of Physics, Georgia Institute of Technology, Atlanta, GA, USA
GSK.ai, 25 Basel St, Petah Tikvah, Israel
<EMAIL_ADDRESS>

###### Abstract

Ising models are a simple generative approach to describing interacting binary
variables. They have proven useful in a number of biological settings because
they enable one to represent observed many-body correlations as the separable
consequence of many direct, pairwise statistical interactions. The inference
of Ising models from data can be computationally very challenging and often
one must be satisfied with numerical approximations or limited precision. In
this paper we present a novel method for the determination of Ising parameters
from data, called GNisi, which uses a Graph Neural network trained on known
Ising models in order to construct the parameters for unseen data. We show
that GNisi is more accurate than the existing state of the art software, and
we illustrate our method by applying GNisi to gene expression data.

##  1 Introduction

The ability to infer graphical models which effectively model interactions in
a dataset is a powerful numerical tool across many scientific disciplines. The
Ising model, as it is known in the language of statistical physics, is a
particularly well-known approach which can be used to understand the
underlying correlations in any form of binary data (“spins”) via pairwise and
single-site interactions. Inferring the parameters of an Ising model from
data, however, is computationally demanding and often intractable, as it is an
inverse problem that typically relies on sampling from the Boltzmann
distribution the Ising parameters imply.

Existing techniques to solve the inverse Ising model approximately include
mean-field inference (Kappen and Rodríguez, 1998; Lee and Daniels, 2019),
pseudo-likelihood (Aurell and Ekeberg, 2012; Ravikumar et al., 2010), and
perturbative expansion methods (Nguyen and Berg, 2012) which can capture the
general form of the graph. More accurate methods include Monte Carlo
simulation (Sutto et al., 2015; Mann et al., 2014; Ackley et al., 1985) and
the adaptive cluster expansion method (Barton et al., 2016), which, while more
accurate may be slow to converge or numerically unstable.

Given the utility of inverse Ising models for deriving direct interactions
from observed covariations in large numbers of variables, it is natural that
there are have been many successful applications in the biological domain.
Such approaches have been used for determining protein structure (Obermayer
and Levine, 2014; Marks et al., 2011; Qin and Colwell, 2018; Shakhnovich and
Gutin, 1993), viral fitness in HIV (Mann et al., 2014; Barton et al., 2015),
correlated states in neural populations (Schneidman et al., 2006),
correlations in nucleotide sequences (Poon et al., 2007, 2008) as well as
modelling a statistical description of gene expression levels (Lezon et al.,
2006).

Graph networks are a deep learning architecture designed to learn from graph-
structured data. They were originally proposed in (Gori et al., 2005;
Scarselli et al., 2008) (see _e.g._ , (Battaglia et al., 2018; Hamilton, ; Wu
et al., 2020) for an overview). The application of graph networks to predict
the behavior of complex physical systems has succeeded in a wide range of
settings (Sanchez-Gonzalez et al., 2020; Pfaff et al., 2021), including the
modelling of spin-glass (Bapst et al., 2020). In this paper we present a novel
method for solving inverse Ising problems using graph networks, which we
denote GNisi, which alleviates the slow convergence or inaccuracy of previous
methods. In our model, the graph network is trained on known Ising models,
generated using Monte Carlo methods, and then the trained network is used to
determine the Ising parameters for unseen binary data.

We will firstly compare our method to existing methods for solving inverse
Ising problems on a generated dataset against which we can compare to ground
truth. We then apply our method to genomic data from the Cancer Cell Line
Encyclopaedia project (CCLE) (Ghandi et al., 2019; Li et al., 2019; Barretina
et al., 2012; Stransky et al., 2015) under CC BY 4.0 licence
111https://creativecommons.org/licenses/by/4.0/legalcode, in order to model
the covariation of gene expression from the anonymized (DepMap, 2020) dataset,
for two subsets of genes that are known to be highly correlated.

##  2 Theoretical outline

###  2.1 The Ising model

For a dataset of bits, $\bm{x}=\left[x_{1},\dots,x_{n}\right]$, the Ising
model for that data is given by the probability distribution

$p(\bm{x})=\frac{1}{Z}\exp\left(-\sum_{i}h_{i}x_{i}-\sum_{i<j}u_{ij}x_{i}x_{j}\right)\,,$
(1)

where the partition function is the normalisation such that probabilities sum
to identity

$Z=\sum_{i}\exp^{\beta H_{i}}\,,$ (2)

where $\beta$ is the inverse temperature, and the Hamiltonian is defined as

$H=-\sum_{i}h_{i}x_{i}-\sum_{i<j}u_{ij}x_{i}x_{j}\,.$ (3)

One can show the probability distribution, Eq. (1) maximizes the Shannon
entropy

$S=-\sum_{\bm{x}}p(\bm{x})\log p(\bm{x})\,,$ (4)

under the constraint that the first and second moments of $\bm{x}$ are
appropriately defined. In particular, we require the sample means and
covariances to match the connected first and second moments

$\langle x_{i}\rangle=\sum_{\bm{x}}p(\bm{x})x_{i}$ (5) $\langle
x_{i}x_{j}\rangle-\langle x_{i}\rangle\langle
x_{j}\rangle=\sum_{\bm{x}}p(\bm{x})x_{i}x_{j}\,.$ (6)

The problem is then to find the parameters $u_{ij},h_{i}$, such that the first
and second moments match the mean values of $x_{i}$ and $x_{i}x_{j}$ observed
in the data. A much simpler task is to compute the pairwise term, $u_{ij}$ by
inverting the covariance matrix of the dataset,
$\text{cov}^{-1}(x_{i},x_{j})$, if it is invertible for the system. One may
show analytically that, in the limit of high temperature, the two matrices are
equivalent, however, there are many situations in which the data distribution
will not be well-described by such a Gaussian form.

###  2.2 Graph neural networks

Let a graph be defined as $\mathcal{G}=(\mathcal{V},\mathcal{E})$ where
$\mathcal{V}=\\{1,\dots,N\\}$ is the set of nodes and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of edges
connecting nodes in $\mathcal{V}$. We denote $\mathcal{N}_{i}$ as the set of
neighbours of node $i$. Using as similar notation as in (Battaglia et al.,
2018), we let $\bm{v}_{i}\in\mathbb{R}^{n_{v}}$ represent the attributes of
node $i$ for all $i\in\mathcal{V}$ and $\bm{e}_{ij}\in\mathbb{R}^{n_{e}}$
represent the attributes of edge $(i,j)$ for all $(i,j)\in\mathcal{E}$. We
illustrate the architecture in Fig. 1. Then, a graph network can be
characterized in terms of edge and node updates as

$\displaystyle\bm{e}_{ij}^{+}$
$\displaystyle=\phi^{e}\big{(}\bm{e}_{ij},\bm{v}_{i},\bm{v}_{j}\big{)},$
$\displaystyle\forall(i,j)\in\mathcal{E}$ (7a) $\displaystyle\bm{v}_{i}^{+}$
$\displaystyle=\phi^{v}\big{(}\rho^{e\to
v}\big{(}\\{\bm{e}_{ij}^{+}\\}_{j\in\mathcal{N}_{i}}\big{)},\bm{v}_{i}\big{)},$
$\displaystyle\forall i\in\mathcal{V}$ (7b)

where $\phi^{e}:\mathbb{R}^{n_{e}+2n_{v}}\to\mathbb{R}^{n_{e}}$,
$\phi^{v}:\mathbb{R}^{n_{e}+n_{v}}\to\mathbb{R}^{n_{v}}$, are the update
functions to be learned and $\rho^{e\to v}$ is an aggregation function
reducing a set of elements to a single one via some input’s permutation
equivariant transformation. In the particular case of Ising models, we are
working with undirected graphs. If, as we do here, one wishes to infer the
existence of the edges (_i.e._ determine whether two bits are correlated),
then we must start with a fully-connected graph, where
$\mathcal{N}_{i}=\mathcal{V}\setminus\\{i\\}\forall i$.

(a) Edge update

(b) Node update

Figure 1: Illustration of graph network updates. In red the object being
updated, and in blue the quantities used to perform the update.

##  3 Methods

###  3.1 Metropolis Monte Carlo simulations

In this paper we aim to learn the Ising model parameters, $h_{i},u_{ij}$ of an
unseen system by training a graph network on a wide range of known Ising
models. To generate the Ising models, we use the Metropolis Monte Carlo
algorithm (Metropolis et al., 1953), which calculates valid bit-strings for an
Ising model using single-spin flip dynamics. In each Monte Carlo iteration of
the algorithm, a single spin is selected at random and the change in the
energy of the system, $\Delta E$ is calculated as the difference between the
Hamiltonian (3) for the two systems. If $\Delta E\leq 0$, the change is
favourable and the spin is flipped. Otherwise, the spin is flipped only if
$\text{e}^{-\beta\Delta E}>x$, where $x$ is a random number in the range
$[0,1]$. This process is repeated for every site in the lattice, until the
model has converged and the energy can no longer decrease.

In Fig. 2 we show representative training runs for 2 Ising models at low and
high temperatures. We observe that the energy of the high temperature model
(left) does not converge, as one would expect, whereas the low temperature
system (right) converges rapidly.

Figure 2: Convergence of the energy for Ising models at high (left) and low
(right) temperatures.

In order to allow the network to learn a wide range of Ising models, we
generate the training data for a wide range of temperatures, lattice sizes and
sparsities. The Monte Carlo data is converted to fully-connected graphs for
the purposes of training; the initial node features of each graph are the
states of the spins $(0,1)$, whilst the edge attributes are the value of the
coupling determined from the Ising model parameters.

###  3.2 Training procedure

For each Ising model we present the network with 1000 samples of converged
spins generated using the Metropolis Monte Carlo algorithm, split between 800
samples in training and 200 validation. We minimize the $\mathcal{L}_{2}$ loss
between the predicted and true Ising parameters and we use early-stopping,
selecting the model with the highest accuracy on the validation set. The test
set consists of entirely unseen models, each with 1000 samples.

The full graph network is composed of 6 graph layers, with the node and edge
encoding each composed of 1 hidden layer with dimensionality
[123,119,28,126,126,126]. The ReLU activation function is used and we use the
Adam optimizer (Kingma and Ba, 2014) with a fixed learning rate of $1e^{-3}$.
The number of nodes in the hidden layers as well as number of graph layers,
the learning rate and optimizer are determined by optimising the validation
loss of the models with $\beta=1$ and 50% sparsity using the Optuna (Akiba et
al., 2019) framework. Within each graph layer, each node is updated based on
its previous embedding and the sum of the incoming edges, whilst the edges are
updated based on their previous embeddings and the embeddings of the connected
nodes.

##  4 Experimental results

In this section we will show the results of the GNisi framework on two
separate datasets. Firstly, we will discuss how GNisi performs on a test set
composed of unseen Ising models with differing physical properties where
ground truth is known. We will additionally compare with existing codes which
solve the inverse Ising task using adaptive cluster expansion (Barton et al.,
2016; Lee and Daniels, 2019), which are the current state of the art in
solving this class of problems. Secondly we apply GNisi to data from the
Cancer Cell Line Encyclopaedia project (CCLE) (Ghandi et al., 2019; Li et al.,
2019; Barretina et al., 2012; Stransky et al., 2015), which consists of 1376
cell lines, in order to model the covariation of expression of a set of highly
correlated genes.

###  4.1 Results on unseen Ising models

Quantifying how well we are able to generate an accurate Ising model for a set
of bit-strings can be done in several ways. Matching the ground truth matrix
matrix-element by matrix-element can of course only be assessed when one has
ground truth to compare with, which is not the case in general. As discussed
in Sec. 2.1, in the maximum entropy approach, the aim is to match the first
and second moments of $\bm{x}$ such that the Shannon entropy of the system is
maximized. Here we will aim to show that with GNisi we are able to
additionally match the sampled third order moments of $\bm{x}$

$\langle x_{i}x_{j}x_{k}\rangle=\sum_{\bm{x}}p(\bm{x})x_{i}x_{j}x_{k}\,,$ (8)

which provides us with much higher level of precision of how closely the
generated Ising model matches the ground truth data. Furthermore, we will show
that we are able to closely match the underlying Boltzmann distribution of the
data, Eq. (15). We note that, in order to prevent any additional overfitting
on the test set, the bit-string samples used to construct the Ising model,
which are generated from the Metropolis Monte Carlo algorithm, are not used in
quantifying the goodness of fit. In other words, we generate a distinct set of
bit-string samples when constructing the Boltzmann distributions.

We compare our results for a representative model from the test set, namely an
Ising model of size 50, at a low temperature and 25% sparsity in the matrix
elements. We can compare the Ising model obtained by GNisi with the ones
obtained with the two current open source packages for solving inverse Ising
models; ACE v.1 (Barton et al., 2016) and coniii v.2.3.0 (Lee and Daniels,
2019), both released under MIT licence. In order to provide a fair comparison,
we run the packages using default settings when possible, in order to prevent
fine-tuning the results on a model-by-model basis as there is no equivalent to
fine-tuning in GNisi. Fine-tuning of the two codes does not, at least for the
data we present here, improve the results at the cost of vastly increased
runtime. As the coniii package contains differing sampling methods for solving
the inverse Ising method, we choose a distinct method from ACE, the Monte
Carlo histogram method (Broderick et al., 2007).

In Fig. 9 we show reconstructed Ising models for the system using the
differing methods, along with the ground truth matrix in panel (b).

(a) GNisi

(b) Ground truth

(c) ACE

(d) coniii

Figure 3: Reconstructed Ising models for a system of size 50 with $\beta=10$
and 25% sparsity in the matrix elements.

The differences in the matrices in Fig. 9 are hard to quantify by eye, however
one can clearly tell ACE and coniii do not capture the correct overall
behaviour of the system. Indeed, the mean squared error between the ground
truth and the different predictions quantifies this as shown in Table 6,
showing that GNisi is better than the other methods at estimating the
collective behavior of the system. This is an advantage of our methodology as
we do not rely on matching means and covariances; by directly using the data,
we retain much more information about the system than can be contained in the
covariance matrix over the samples.

Table 1: MSE and Pearson correlation coefficient between the Ising parameters computed using the different methods and ground truth. Method | MSE | $r$
---|---|---
GNisi | 15.38 | 0.04
ACE | 101.53 | -0.02
coniii | 179.60 | 0.02

It is useful to compute the Pearson correlation coefficient $r$ between the
ground truth and reconstructed matrices which we show in Table 6 and suggests
that GNisi is not getting the individual matrix-elements correct. In Fig. 10
we plot a scatter plot for the ground truth and predicted Boltzmann
probability distributions, Eq. (1), for a sample of possible bit-strings,
distinct from the bit-strings used to generate the Ising model for each
method, as in general, we can evaluate the Boltzmann distribution for any
possible bit-string. This can be quantified by computing the partition
function, Eq. (15) as well as the Pearson correlation coefficient, which we
present in Table 7. We observe that GNisi far outperforms the other methods
and captures well the overall distribution. We can conclude that, despite the
fact that GNisi does not reproduce individual matrix-elements, it matches the
Boltzmann distribution of the ground truth very well. Therefore, it _does_
produce a valid Ising model for the data it is given and we see that GNisi
produces the most accurate Ising model compared to the other methods.

Table 2: Partition functions and correlation coefficients computed for an Ising model with differing methods. Method | $\log(Z)$ | $r$
---|---|---
Ground truth | 55.37 | $-$
GNisi | 52.39 | $0.44$
ACE | 70.52 | $-0.35$
coniii | 271.41 | $-0.07$

(a) GNisi

(b) ACE

(c) coniii

Figure 4: Scatter plots of probability distributions, Eq. (1) computed using
the different methods for a sample of possible bit-strings.

A final point of comparison is comparing the first, second and third connected
moments, $m_{n}$, of $\bm{x}$ (Eqs. (5), (6) and (8) respectively) and
comparing to ground truth. ACE and coniii are constructed such that they
produce an Ising model which matches the first and second moments, however
GNisi does not have this requirement imposed in the architecture, and it is
therefore a non-trivial check of our methodology that we can do so. In
addition, we consider the third moment of $\bm{x}$ as a method to more
powerfully compare our methodology with the ground truth. For the second
moments, for an $n$-dimensional lattice, there are ${n\choose 2}$ possible
combinations, which is not prohibitive to compute for our system. At third
order, however, we have to consider ${n\choose 3}$, and we instead will
consider only a sample of the full distribution. In Table 8 we show the mean
squared error between the moments computed using the ground truth, and with
the various methods. In all cases we show the $n$th connected moment with
either all bit-strings being 1 ($m_{n}(1)$), or 0 ($m_{n}(0)$).
Unsurprisingly, both ACE and coniii match the moments very well, as they
construct Ising models under the condition that the maximum entropy condition
is met. It is therefore impressive that GNisi, which does not attempt to match
moments, also matches the ground truth to very high precision. We show in the
Supplementary Material the results obtained using a random matrix with the
same size couplings as GNisi, to explicitly show that GNisi is correctly
learning features of the Ising model. Interestingly, ACE and coniii match
extremely well (to the precision shown here), suggesting the maximum entropy
approach strongly constrains higher order moments.

In order to explain the behavior we observe, namely that all the methods used
in this paper can match the moments of the data distribution well but, to a
greater or lesser extent, not match the individual parameters of the ground
truth Ising model, we note that we have a sloppy model (Brown and Sethna,
2003; Brown et al., 2004). A sloppy model is one which is poorly constrained
as there exist a very large number of parameter choices which can match the
data distribution well. As we need to match $n^{2}$ parameters for an Ising
model of size $n$, this behaviour can be expected to worsen as we increase the
size of the system. It is therefore not surprising that the 3 different
methods match the moments of the ground truth matrix with 3 distinct
predictions for the Ising model. It is for this reason we argue that the
ability to correctly reconstruct the Boltzmann distribution is the fundamental
test for whether a generated model is accurate; in Fig. 10 we observe that the
predictions from GNisi correlates with the ground truth bit-string by bit-
string, which strongly suggests that GNisi has correctly learned underlying
features about the Ising model.

Table 3: Mean squared error between the $n$th moment computed using ground truth and the differing methods.  MSE |
---|---
Moment | GNisi | ACE | coniii
$m_{1}(1)$ | 0.03 | 0.11 | 0.11
$m_{2}(1)$ | $1.60\cdot 10^{-5}$ | $3.27\cdot 10^{-7}$ | $3.27\cdot 10^{-7}$
$m_{2}(0)$ | $1.60\cdot 10^{-5}$ | $3.27\cdot 10^{-7}$ | $3.27\cdot 10^{-7}$
$m_{3}(0)$ | $6.33\cdot 10^{-5}$ | $6.83\cdot 10^{-7}$ | $6.83\cdot 10^{-7}$
$m_{3}(1)$ | 0.0 | 0.0 | 0.0

###  4.2 Results on the CCLE dataset

In this section we apply the pre-trained GNisi model to two distinct datasets
of genes from the CCLE dataset, each of which are known to be highly mutually
correlated in order to provide an estimate of the ground truth. We binarize
the continuous data by identifying the $q$-level quantile, and setting the
gene expression to 0 below the $q$-th quantile, and 1 above it. In order to
prevent the binarized data from being highly skewed, we set $q=0.25$
throughout. As above, where we compare GNisi and ACE against ground truth, we
will show results run using ACE default settings. Due to the poor quality of
the results with coniii on the synthetic data, we will not include a
comparison with the package here.

####  4.2.1 Highly correlated genes

In Fig. 5 we show reconstructed Ising models for the first system using the
two differing methods. We can see that both GNisi and ACE correctly produce a
highly correlated matrix, however the couplings produced by ACE are an order
of magnitude larger than those produced by GNisi.

(a) GNisi

(b) ACE

Figure 5: Reconstructed Ising models for a group of genes from the CCLE
dataset known to be highly mutually correlated.

We argued in the previous section that the fundamental test to determine the
accuracy of our reconstructed Ising model is to compare the Boltzmann
distribution constructed from the Ising model against the data. Unlike in Sec.
4.1, where we compare samples of the partition function against ground truth
in Fig. 10, in this example (and all real-world applications) we do not know
the true Ising model which describes the data. In Fig. 6 we show a sample of
the log Boltzmann probability distributions computed using the observed data
samples and the reconstructed Ising models from GNisi and ACE. As there are
$N=1376$ samples in the observed dataset, we sample $N$ possible bit-strings
from the total number of combinations and compute the probability of each
permutation, given the underlying Ising model for GNisi and ACE. It is clear
that, while GNisi does not perfectly match the observed distribution, as the
distribution has a smaller tail, it matches very closely. On the other hand,
ACE is strongly peaked at $-\log(p)\sim 7$, which is correct, but has no other
features.

We now turn to the first, second and third moments predicted by the different
models. In Table 4 we show mean squared error between the sampled moments and
the observed moments from the data distribution. We see that on the whole,
both methods agree well with the data, with GNisi having a larger errors for
most moments. Interestingly, for this subset of genes, we find that ACE
outperforms GNisi with respect to matching moments as shown in Table 4. It
should be noted, however, that matching the low-order moments of the training
data does not necessarily capture the data distribution better by all
measures, and there was no ground truth of a known underlying Ising model to
compare to in this case of real experimental data. Indeed, in this case it is
clear that GNisi produces the more accurate Ising model when judged by how
close the predicted log Boltzmann probability distribution is to the observed
distribution.

Table 4: Mean squared error between the $n$th moment computed using ground truth and ACE for the CCLE data.  MSE |
---|---
Moment | GNisi | ACE
$m_{1}(1)$ | 0.36 | 0.43
$m_{2}(1)$ | $4.15\cdot 10^{-4}$ | $1.58\cdot 10^{-6}$
$m_{2}(0)$ | $4.17\cdot 10^{-4}$ | $1.71\cdot 10^{-6}$
$m_{3}(1)$ | $3.29\cdot 10^{-4}$ | $3.29\cdot 10^{-6}$
$m_{3}(0)$ | $1.87$ | $9.13\cdot 10^{-4}$
Figure 6: Log Boltzmann probability distribution for the different methods
using CCLE data.

####  4.2.2 Genes with known mutual relation to each other and to BRCA1

Finally, we turn to the second distinct set of genes. As before, some
biological ground truth is known about the 21 genes in this dataset; here we
pick the subset of genes with known functional relation to BRCA1 and to each
other. In Fig. 7 we show the reconstructed Ising models using GNisi and ACE;
as before we see that both methods find a high level of mutual correlation in
the matrix, however GNisi finds much stronger couplings than ACE.

(a) GNisi

(b) ACE

Figure 7: Reconstructed Ising models for a group of genes from the CCLE
dataset known to be highly mutually correlated with BRCA1.

To again determine the accuracy of the reconstructions, we plot the observed
and sampled log Boltzmann probability distributions in Fig. 8. As with the
previous dataset, GNisi matches the log Boltzmann probability distribution for
the data better than ACE, although ACE performs much better with this dataset
than the previous one. We again see that GNisi matches the shape of the
distribution somewhat better than ACE. We show the MSE between moments for
this dataset in Table 5; again both GNisi and ACE perform well, with ACE
somewhat worse than GNisi in this case.

Figure 8: Log Boltzmann probability distribution for the different methods using CCLE data for genes known to be highly mutually correlated with BRCA1. Table 5: Mean squared error between the $n$th moment computed using ground truth and ACE for the CCLE data for genes known to be highly mutually correlated with BRCA1.  MSE |
---|---
Moment | GNisi | ACE
$m_{1}(1)$ | 0.43 | 0.27
$m_{2}(1)$ | $7.51\cdot 10^{-5}$ | $5.16\cdot 10^{-4}$
$m_{2}(0)$ | $7.55\cdot 10^{-5}$ | $5.18\cdot 10^{-4}$
$m_{3}(1)$ | $2.50\cdot 10^{-4}$ | $2.50\cdot 10^{-4}$
$m_{3}(0)$ | $2.89\cdot 10^{-5}$ | $0.027$

##  5 Conclusion

In this paper we have presented a novel method for solving the inverse Ising
problem. Inverse Ising methods have been successfully applied to a wide range
of biological systems, however the computational difficulty of solving the
inverse problem often prevents their use, or necessitates the use of
approximate methods such as Gaussian approximations or pseudo-likelihood
methods. With GNisi, we have presented a fast, accurate and easy-to-use method
for solving inverse Ising models using graph neural networks.

We have presented results where we apply GNisi to both synthetic and real
data, and compared our results with the current state-of-the-art codes. We
find that across the board, GNisi outperforms other methods. Additionally,
GNisi requires no fine-tuning; once the data has been collected, it is
sufficient to provide it to the trained model and an Ising model that
accurately describes the data is rapidly produced. It is important to note
that the results presented here for GNisi are with 1000 Monte Carlo samples
per model; we expect that results will improve upon the generation of more
data.

It should be noted that the choice of what random distribution of Ising
couplings to sample in the Monte Carlo simulation data used to train GNisi
should affect the kinds of Ising models it is capable of reconstructing: in an
instance where the best Ising model of the data would be a highly rare event
in the space of models considered in the training data, GNisi might be
expected to fail at constructing a useful model, whether for protein structure
or neuronal firing.

## References

  * Ackley et al. [1985] David H. Ackley, Geoffrey E. Hinton, and Terrence J. Sejnowski.  A learning algorithm for boltzmann machines.  _Cognitive Science_ , 9(1):147–169, 1985.
  * Akiba et al. [2019] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2019.
  * Aurell and Ekeberg [2012] Erik Aurell and Magnus Ekeberg.  Inverse ising inference using all the data.  _Phys. Rev. Lett._ , 108:090201, Mar 2012.
  * Bapst et al. [2020] V. Bapst, T Keck, A Grabska-Barwińska, C Donner, E.D. Cubuk, S.S. Schoenholz, A. Obika, A.W.R. Nelson, T. Back, D Hassabis, and P. Kohli.  Unveiling the predictive power of static structure in glassy systems.  _Nature Physics_ , 16(4):448–454, 2020.
  * Barretina et al. [2012] Jordi Barretina et al.  The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity.  _Nature_ , 483(7391):603–607, Mar 2012.  ISSN 1476-4687.
  * Barton et al. [2016] J. P. Barton, E. De Leonardis, A. Coucke, and S. Cocco.  ACE: adaptive cluster expansion for maximum entropy graphical model inference.  _Bioinformatics_ , 32(20):3089–3097, 06 2016\.
  * Barton et al. [2015] John P. Barton, Mehran Kardar, and Arup K. Chakraborty.  Scaling laws describe memories of host–pathogen riposte in the hiv population.  _Proceedings of the National Academy of Sciences_ , 112(7):1965–1970, 2015.  ISSN 0027-8424.
  * Battaglia et al. [2018] Peter W. Battaglia et al.  Relational inductive biases, deep learning, and graph networks.  2018\.
  * Broderick et al. [2007] Tamara Broderick, Miroslav Dudik, Gasper Tkacik, Robert E. Schapire, and William Bialek.  Faster solutions of the inverse pairwise ising problem, 2007.
  * Brown and Sethna [2003] K. S. Brown and J. P. Sethna.  Statistical mechanical approaches to models with many poorly known parameters.  _Physical Review E_ , 68:021904, 2003.
  * Brown et al. [2004] K S Brown, C C Hill, G A Calero, C R Myers, K H Lee, J P Sethna, and R A Cerione.  The statistical mechanics of complex signaling networks: nerve growth factor signaling.  _Physical Biology_ , 1(3):184–195, Oct 2004\.  ISSN 1478-3975.
  * DepMap [2020] Broad DepMap.  Depmap 20q4 public, Nov 2020.  URL https://figshare.com/articles/dataset/DepMap_20Q4_Public/13237076/4.
  * Ghandi et al. [2019] Mahmoud Ghandi et al.  Next-generation characterization of the cancer cell line encyclopedia.  _Nature_ , 569(7757):503–508, May 2019.  ISSN 1476-4687.
  * Gori et al. [2005] Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In _Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005._ , volume 2, pages 729–734. IEEE, 2005.
  * [15] William L. Hamilton.  Graph representation learning.  _Synthesis Lectures on Artificial Intelligence and Machine Learning_ , 14(3):1–159.
  * Kappen and Rodríguez [1998] H. J. Kappen and F. B. Rodríguez.  Efficient Learning in Boltzmann Machines Using Linear Response Theory.  _Neural Computation_ , 10(5):1137–1156, 07 1998\.
  * Kingma and Ba [2014] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.  cite arxiv:1412.6980Comment: Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015.
  * Lee and Daniels [2019] E.D. Lee and B.C. Daniels.  Convenient Interface to Inverse Ising (ConIII): A Python 3 Package for Solving Ising-Type Maximum Entropy Models.  _Journal of Open Research Software_ , 7(1):3, 2019\.
  * Lezon et al. [2006] Timothy R. Lezon, Jayanth R. Banavar, Marek Cieplak, Amos Maritan, and Nina V. Fedoroff.  Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns.  _Proceedings of the National Academy of Sciences_ , 103(50):19033–19038, 2006.  ISSN 0027-8424.
  * Li et al. [2019] Haoxin Li et al.  The landscape of cancer cell line metabolism.  _Nature Medicine_ , 25(5):850–860, May 2019.  ISSN 1546-170X.
  * Mann et al. [2014] Jaclyn K. Mann, John P. Barton, Andrew L. Ferguson, Saleha Omarjee, Bruce D. Walker, Arup Chakraborty, and Thumbi Ndung’u.  The fitness landscape of hiv-1 gag: Advanced modeling approaches and validation of model predictions by in vitro testing.  _PLOS Computational Biology_ , 10(8):1–11, 08 2014.
  * Marks et al. [2011] Debora S. Marks, Lucy J. Colwell, Robert Sheridan, Thomas A. Hopf, Andrea Pagnani, Riccardo Zecchina, and Chris Sander.  Protein 3d structure computed from evolutionary sequence variation.  _PLOS ONE_ , 6(12):1–20, 12 2011.
  * Metropolis et al. [1953] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.
  * Nguyen and Berg [2012] H Chau Nguyen and Johannes Berg.  Bethe-peierls approximation and the inverse ising problem.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2012(03):P03004, mar 2012.
  * Obermayer and Levine [2014] Benedikt Obermayer and Erel Levine.  Inverse ising inference with correlated samples.  _New Journal of Physics_ , 16(12):123017, Dec 2014\.  ISSN 1367-2630.
  * Pfaff et al. [2021] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter Battaglia.  Learning mesh-based simulation with graph networks.  In _International Conference on Learning Representations_ , 2021.
  * Poon et al. [2007] Art F. Y Poon, Fraser I Lewis, Sergei L. Kosakovsky Pond, and Simon D. W Frost.  An evolutionary-network model reveals stratified interactions in the v3 loop of the hiv-1 envelope.  _PLOS Computational Biology_ , 3(11):1–12, 11 2007.
  * Poon et al. [2008] Art F. Y. Poon, Fraser I. Lewis, Simon D. W. Frost, and Sergei L. Kosakovsky Pond.  Spidermonkey: rapid detection of co-evolving sites using Bayesian graphical models.  _Bioinformatics_ , 24(17):1949–1950, 06 2008\.  ISSN 1367-4803.
  * Qin and Colwell [2018] Chongli Qin and Lucy J. Colwell.  Power law tails in phylogenetic systems.  _Proceedings of the National Academy of Sciences_ , 115(4):690–695, 2018.  ISSN 0027-8424.
  * Ravikumar et al. [2010] Pradeep Ravikumar, Martin J. Wainwright, and John D. Lafferty.  High-dimensional Ising model selection using l1-regularized logistic regression.  _The Annals of Statistics_ , 38(3):1287 – 1319, 2010.
  * Sanchez-Gonzalez et al. [2020] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia.  Learning to simulate complex physics with graph networks.  In Hal Daumé III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 8459–8468. PMLR, 13–18 Jul 2020.
  * Scarselli et al. [2008] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  _IEEE transactions on neural networks_ , 20(1):61–80, 2008.
  * Schneidman et al. [2006] Elad Schneidman, 2nd Berry, Michael J, Ronen Segev, and William Bialek.  Weak pairwise correlations imply strongly correlated network states in a neural population.  _Nature_ , 440(7087):1007–1012, Apr 2006.  ISSN 1476-4687.  16625187[pmid].
  * Shakhnovich and Gutin [1993] E I Shakhnovich and A M Gutin.  Engineering of stable and fast-folding sequences of model proteins.  _Proceedings of the National Academy of Sciences_ , 90(15):7195–7199, 1993.  ISSN 0027-8424.
  * Stransky et al. [2015] Nicolas Stransky et al.  Pharmacogenomic agreement between two cancer cell line data sets.  _Nature_ , 528(7580):84–87, Dec 2015.  ISSN 1476-4687.
  * Sutto et al. [2015] L. Sutto, S. Marsili, A. Valencia, and F.L. Gervasio.  From residue coevolution to protein conformational ensembles and functional dynamics.  _Proceedings of the National Academy of Sciences of the United States of America_ , 112(44), 2015.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 2020\.

##  Appendix A Comparison between GNisi and ACE for a high temperature system

In the main paper we compared our results for a representative model from the
test set at low temperature. Here we will present the same results, performing
the same comparisons, for an Ising model from the test set at a temperature
one of order magnitude higher, namely an Ising model of size 50 with 75%
sparsity in the matrix elements for $\beta=1$. As before, we run the ACE and
coniii packages using default settings, in order to prevent fine-tuning the
results on a model-by-model basis. The MSE between the ground truth and
differing methods is shown in Table 6.

In Fig. 9 we show reconstructed Ising models for the system using the
differing methods, along with the ground truth matrix in panel (b).

(a) GNisi

(b) Ground truth

(c) ACE

(d) coniii

Figure 9: Reconstructed Ising models for a system of size 50 with $\beta=1$ and 75% sparsity in the matrix elements. Table 6: MSE between the Ising parameters computed using the different methods and ground truth. Method | MSE | $r$
---|---|---
GNisi | 14.78 | 0.03
ACE | 69.98 | 0.03
coniii | 408.06 | 0.01

In Fig. 10 we plot the scatter of predicted versus ground truth probability
distributions, for a sample of possible bit-strings, which are again distinct
from the bit-strings used to generate the Ising model for each method. For
this model, ACE and GNisi perform very similarly as can also be observed in
Table 7 where we show also the Pearson correlation coefficient, $r$.

Table 7: Partition functions computed for an Ising model with differing methods Method | $\log(Z)$ | $r$
---|---|---
Ground-truth | 42.86 | $-$
GNisi | 46.48 | 0.11
ACE | 46.73 | 0.10
coniii | 252.78 | 0.02

(a) GNisi

(b) ACE

(c) coniii

Figure 10: Scatter plots of probability distributions computed using the
different methods for a sample of possible bit-strings.

Finally we will compare first, second and third moments. In Table 8 we show
the mean squared error between the moments computed using the ground truth,
and with the various methods. Again, GNisi does very well at matching moments
despite it not being a restriction when constructing the Ising model.

Table 8: Mean squared error between the $n$th moment computed using ground truth and the differing methods.  MSE |
---|---
Moment | GNisi | ACE | coniii
$m_{1}(1)$ | 0.07 | 0.12 | 0.13
$m_{2}(1)$ | 1.32$\cdot 10^{-4}$ | $2.98\cdot 10^{-7}$ | $2.98\cdot 10^{-7}$
$m_{2}(0)$ | $1.32\cdot 10^{-4}$ | $2.98\cdot 10^{-7}$ | $2.98\cdot 10^{-7}$
$m_{3}(1)$ | $2.6\cdot 10^{-3}$ | $1.54\cdot 10^{-6}$ | $1.54\cdot 10^{-6}$
$m_{3}(0)$ | 0.0 | 0.0 | 0.0

##  Appendix B Accuracy of inverse covariance matrix approximation

In this section we will compare the accuracy of inverting the covariance
matrix over samples compared to ground truth and GNisi. For a larger system of
size 50, there are $2^{50}\sim 10^{15}$ equally likely bit-string combinations
at high temperature, which is computationally intractable to compute, and one
therefore will not expect the inverse covariance matrix computed with 1000
samples to be accurate. We will instead here focus on a much smaller system,
of lattice size 5 as the 32 possible combinations of bit strings enable us to
perform a fair comparison between GNisi and the inverse covariance matrix. The
model shown here is at inverse temperature $\beta=1$ with 55% sparsity in the
matrix elements.

(a) GNisi

(b) Ground truth

(c) $\text{inv}(\text{cov}_{ij})$

Figure 11: Reconstructed Ising models for a system of size 5 with $\beta=1$
and 55% sparsity in the matrix elements.

We can see in Fig. 11 that GNisi matches very well the ground truth, whilst
the inverse of the covariance matrix performs very poorly. This is more
explicitly shown in Table 9 where we compare the Boltzmann distributions
computed using the two methods. We see that GNisi performs well overall,
whilst the inverse of the covariance matrix captures nothing about the true
distribution.

Table 9: Partition functions computed for an Ising model with GNisi and the inverse of the covariance matrix. Method | $\log(Z)$
---|---
Ground-truth | 19.64
GNisi | 12.84
$\text{inv}(\text{cov}_{ij})$ | 0.00

##  Appendix C Derivation of maximum entropy approach

The maximum entropy approach states that, under the condition the first and
second moments of the samples match the true connected moments, the Ising
model maximises the Shannon entropy of the system. In this section we will
derive this result.

Our constraints, in the maximum entropy approximation, are

$\displaystyle\langle x_{i}\rangle$
$\displaystyle=\sum_{\bm{x}}p(\bm{x})x_{i}$ (9) $\displaystyle\langle
x_{i}x_{j}\rangle$ $\displaystyle=\sum_{\bm{x}}p(\bm{x})x_{i}x_{j}$ (10)
$\displaystyle 1$ $\displaystyle=\sum_{\bm{x}}p(\bm{x})\,.$ (11)

The constrained optimisation problem can be formulated as extremising the
Lagrangian of the system, $\mathcal{L}$. We therefore have

$\displaystyle\mathcal{L}$ $\displaystyle=-\sum_{\bm{x}}p(\bm{x})\log
p(\bm{x})-\chi\sum_{\bm{x}}p(\bm{x})-1-\sum_{i}\left(\alpha_{i}\left(\sum_{\bm{x}}p(\bm{x})x_{i}\right)-\langle
x_{i}\rangle\right)-\sum_{i,j}\left(\gamma_{ij}\sum_{\bm{x}}p(\bm{x})x_{i}x_{j}-\langle
x_{i}x_{j}\rangle\right)\,.$

So

$\displaystyle\frac{\partial\mathcal{L}}{\partial p(\bm{x})}$
$\displaystyle=0=-\log
p(\bm{x})-1-\chi-\sum_{i}\alpha_{i}x_{i}-\sum_{i,j}\gamma_{ij}x_{i}x_{j}\,,$
(12)

and

$\displaystyle\log p(\bm{x})=-1-\chi-\alpha.\bm{x}-\bm{x}.\gamma.\bm{x}\,,$
(13) $\displaystyle p(\bm{x})=Ae^{-\alpha.\bm{x}-\bm{x}.\gamma.\bm{x}}\,,$
(14)

which is the Ising model where we identify $\alpha$ as the external field and
$\gamma$ as the pairwise terms respectively.

##  Appendix D Derivation of inverse covariance matrix approximation

In this appendix we derive the approximation where the inverse of the
covariance matrix over nodes is approximately equal to the pairwise coupling
term. We derive the approximation using field theoretic methods.

The partition function for a discrete number of bit-strings is given by the
usual

$Z=\sum_{i=1}^{n}e^{-\beta H_{i}(x_{i})}\,.$ (15)

In the infinite temperature limit, all possible $2^{n}$ combination of spins
are equally possible. We can therefore take $i\to\infty$ and then

$Z=\int\mathcal{D}xe^{-\beta H_{i}(x_{i})}\,,$ (16)

where $\mathcal{D}$ signifies we take path integrals.

The constraints on the first and second moments are now given by

$\displaystyle\langle x_{i}\rangle$
$\displaystyle=\int\mathcal{D}x\,p(\bm{x})x_{i}$ (17) $\displaystyle\langle
x_{i}x_{j}\rangle$ $\displaystyle=\int\mathcal{D}x\,p(\bm{x})x_{i}x_{j}\,,$
(18)

where

$p(\bm{x})=Ae^{-\alpha.\bm{x}-\bm{x}.\gamma.\bm{x}}\,,$ (19)

and $\alpha$ and $\gamma$ are the external field and pairwise terms
respectively. Via a change of coordinates $\bm{y}=\bm{x}+\alpha\xi^{-1}$ where
$\xi=\frac{1}{2}\gamma$ we can write (19) as

$p(\bm{x})=Be^{-\frac{1}{2}\bm{y}.\xi\bm{y}}\,,$ (20)

and so (17) becomes

$\displaystyle\langle x_{i}\rangle$
$\displaystyle=\int\mathcal{D}\bm{y}\,Be^{-\frac{1}{2}\bm{y}.\xi.\bm{y}}(y_{i}-\sum_{j}\alpha_{j}\xi^{-1}_{ij})$
(21)
$\displaystyle=-\int\mathcal{D}\bm{y}\,Be^{-\frac{1}{2}\bm{y}.\xi.\bm{y}}\sum_{j}\alpha_{j}\xi^{-1}_{ij}$
(22)
$\displaystyle=-B\alpha\xi^{-1}\frac{(2\pi)^{N/2}}{\sqrt{\text{det}\xi}}\sum_{j}\alpha_{j}\xi^{-1}_{ij}\,,$
(23)

using the fact that

$\int\mathcal{D}\bm{y}\,Be^{-\frac{1}{2}\bm{y}.\xi.\bm{y}}\bm{y}=0\,,$ (24)

and properties of multi-dimensional Gaussians of size $N$, as we assume in
this approximation our data is Gaussian. We therefore have

$\langle x_{i}\rangle=-C\sum_{j}\alpha_{j}\xi^{-1}_{ij}\,.\\\ $ (25)

Eq. (18) is then

$\displaystyle\langle x_{i}x_{j}\rangle$
$\displaystyle=\int\mathcal{D}\bm{y}Be^{-\frac{1}{2}\bm{y}.\xi\bm{y}}(y_{i}-\sum_{k}\alpha_{k}\xi^{-1}_{ik})(y_{i}-\sum_{l}\alpha_{l}\xi^{-1}_{il})$
(26)
$\displaystyle=\int\mathcal{D}\bm{y}Be^{-\frac{1}{2}\bm{y}.\xi\bm{y}}(y_{i}y_{j}-\sum_{l}y_{i}\alpha_{l}\xi_{jl}^{-1}-\sum_{k}y_{j}\alpha_{k}\xi_{ik}^{-1}+\sum_{k,l}\alpha_{l}\alpha_{k}\xi_{ik}^{-1}\xi_{jl}^{-1})$
(27)
$\displaystyle\sim\int\mathcal{D}\bm{y}Be^{-\frac{1}{2}\bm{y}.\xi\bm{y}}(y_{i}y_{j}+y_{i}\langle
x_{j}\rangle+y_{j}\langle x_{i}\rangle+\langle x_{i}\rangle\langle
x_{j}\rangle)\,,$ (28)

using Eq. (25). Again using the properties of multi-dimensional Gaussians, we
find

$\displaystyle\langle x_{i}x_{j}\rangle$ $\displaystyle=D\langle
x_{i}\rangle\langle
x_{j}\rangle+\int\mathcal{D}\bm{y}Be^{-\frac{1}{2}\bm{y}.\xi\bm{y}}y_{i}y_{j}\,.$
(29)

The remaining integral is a 2-point correlation function which we can solve
analytically. The partition function Eq. (16) is, dropping all constant
factors,

$\displaystyle Z$
$\displaystyle=\int\mathcal{D}\bm{y}e^{-\frac{1}{2}\bm{y}.\xi.\bm{y}+\alpha.\bm{y}}$
(30)
$\displaystyle=\int\mathcal{D}\bm{x}e^{(\bm{x}+\xi^{-1}\alpha).\frac{1}{2}\xi(\bm{x}+\xi^{-1}\alpha)+\alpha.\bm{x}+\alpha\xi^{-1}\alpha}$
(31)
$\displaystyle\sim\int\mathcal{D}\bm{x}e^{-\frac{1}{2}\bm{x}\xi\bm{x}+\frac{1}{2}\alpha\xi^{-1}\alpha}$
(32)
$\displaystyle=e^{\frac{1}{2}\alpha\xi^{-1}\alpha}\int\mathcal{D}\bm{x}e^{-\frac{1}{2}\bm{x}\xi\bm{x}}$
(33)
$\displaystyle=e^{\frac{1}{2}\alpha\xi^{-1}\alpha}\frac{(2\pi)^{N/2}}{\sqrt{\text{det}\xi}}\,.$
(34)

Now, using the definition of functional derivatives from statistical mechanics

$\frac{1}{Z[0]}\frac{\delta^{2}Z[\alpha]}{\delta\alpha(x_{i})\delta\alpha(x_{j})}\bigg{|}_{\alpha=0}:=\frac{1}{Z[0]}\int\mathcal{D}\varphi\varphi(x_{i})\varphi(x_{j})e^{-S[\varphi]}\,,$
(35)

where $S[\varphi]$, the action, can be identified with $H$ and so

$\frac{1}{Z[0]}\frac{\delta^{2}Z[\alpha]}{\delta\alpha(x_{i})\delta\alpha(x_{j})}\bigg{|}_{\alpha=0}=\frac{1}{Z[0]}\int\mathcal{D}\varphi\varphi(x_{i})\varphi(x_{j})e^{-H}\,,$
(36)

where the integral above is equivalent to the integral in (29). Bringing it
all together we finally have

$\displaystyle\langle x_{i}x_{j}\rangle$
$\displaystyle\sim\frac{1}{Z[0]}\frac{\delta^{2}Z[\alpha]}{\delta\alpha(x_{i})\delta\alpha(x_{j})}\bigg{|}_{\alpha=0}+\langle
x_{i}\rangle\langle x_{j}\rangle$ (37)
$\displaystyle=(2\pi)^{-N/2}\left[\text{det}\xi\right]^{\frac{1}{2}}(2\pi)^{N/2}\left[\text{det}\xi\right]^{-\frac{1}{2}}\frac{\partial}{\partial\alpha_{i}}\left(\alpha_{i}\xi_{ij}^{-1}\right)\bigg{|}_{\alpha=0}+\langle
x_{i}\rangle\langle x_{j}\rangle$ (38) $\displaystyle=\xi_{ij}^{-1}+\langle
x_{i}\rangle\langle x_{j}\rangle\,.$ (39)

Therefore, the covariance matrix can be identified as being proportional to
the inverse pairwise term in the Hamiltonian

$C_{ij}:=\langle x_{i}x_{j}\rangle-\langle x_{i}\rangle\langle
x_{j}\rangle\sim\xi_{ij}^{-1}\sim\gamma_{ij}^{-1}\,.$ (40)

##  Appendix E Definition of the 3rd connected moment

In this section we state the equation used to compute the 3rd connected
moments in the main text. We use the coskewness for $n$-random variables:

$m_{3}(i,j,k)=\frac{E[i-E[i]]E[j-E[j]]E[k-E[k]]}{\sigma_{i}\sigma_{j}\sigma_{k}}\,,$
(41)

where $E[i]$ is the expected value of $i$ and $\sigma_{i}$ its standard
deviation.

##  Appendix F Performance of a random Ising model

In the main text, we state that GNisi, for the Ising model against which we
have ground truth, out-performs ACE and coniii. It is reasonable to ask
whether a random matrix with parameters of the same size as GNisi would out-
perform it. In Fig. 12, we show the scatter plot for the Boltzmann probability
distribution for a matrix generated at random with the same size couplings as
GNisi. We observe that the probability distribution is captured very poorly;
indeed we find $\log(Z)=7.05$ compared to the ground truth value of
$\log(Z)=55.37$ with the Pearson coefficient on the Boltzmann probability
distribution of $r=-0.20$. In Table 10 we show the corresponding MSE between
the moments computed using ground truth and the random matrix, again observing
that GNisi is better. We therefore can conclude that GNisi out-performs a
random matrix and is correctly learning the Ising models.

Figure 12: Scatter plot for the Boltzmann probability distribution for a random matrix against ground truth. Table 10: Mean squared error between the $n$th moment computed using ground truth and the random matrix.  MSE |
---|---
Moment | Random matrix
$m_{1}(1)$ | 0.10
$m_{2}(1)$ | $3.5\cdot 10^{-5}$
$m_{2}(0)$ | $3.5\cdot 10^{-5}$
$m_{3}(0)$ | $1.0\cdot 10^{-3}$
$m_{3}(1)$ | 0.0


# Improved Loss Function-Based Prediction Method of Extreme Temperatures in
Greenhouses

Liao Qu  Shuaiqi Huang  Yunsong Jia  Xiang Li  College of Information and
Electrical Engineering, China Agricultural University, Beijing, China
<EMAIL_ADDRESS>
###### Abstract

The prediction of extreme greenhouse temperatures to which crops are
susceptible is essential in the field of greenhouse planting. It can help
avoid heat or freezing damage and economic losses. Therefore, it’s important
to develop models that can predict them accurately. Due to the lack of extreme
temperature data in datasets, it is challenging for models to accurately
predict it. In this paper, we propose an improved loss function, which is
suitable for a variety of machine learning models. By increasing the weight of
extreme temperature samples and reducing the possibility of misjudging extreme
temperature as normal, the proposed loss function can enhance the prediction
results in extreme situations. To verify the effectiveness of the proposed
method, we implement the improved loss function in LightGBM, long short-term
memory, and artificial neural network and conduct experiments on a real-world
greenhouse dataset. The results show that the performance of models with the
improved loss function is enhanced compared to the original models in extreme
cases. The improved models can be used to guarantee the timely judgment of
extreme temperatures in agricultural greenhouses, thereby preventing
unnecessary losses caused by incorrect predictions.

###### keywords:

Greenhouse Temperature Prediction; Improved Loss Function; Extreme Temperature

##  1 Introduction

The greenhouse is a typical scenario of modern agriculture[1]. Among the
various environmental factors, the temperature has the most significant impact
on the growth of crops [2, 3]. Temperature influences many biochemical
processes such as photosynthesis, respiration, water balance, and membrane
stability of leaves [4], which would then affect the quality and yields of
crops. Therefore, accurate temperature prediction is the key to greenhouse
management. Although fluctuations within normal temperature do not affect the
growth of crops, the extreme temperature may result in severe damage to
them[5, 6]. For example, under extremely high temperatures, crops will suffer
thermal damage, which affects the process of pollination and fruiting,
possibly leading to death [7, 8]. Similarly, under extremely low temperatures,
crops will freeze or suffer frostbite[9, 10]. Hence, accurate prediction of
extreme temperatures is conducive to early prevention of such damage. In
recent years, based on computational descriptions of greenhouse environment
[11, 12] and statistical methods [13, 14, 15, 16, 17, 18], various models have
been successfully developed to predict the greenhouse temperature. However, in
contrast to mathematical models based on the internal biology mechanism of the
crops, environment, and production process in the greenhouse, machine learning
models do not require many physical factors [19], they only need historical
environment data of greenhouses. Because using machine learning in temperature
prediction is relatively straightforward, many related studies have been
conducted in recent years. [17] used an artificial neural network(ANN) with
the Levenberg-Marquardt method [20] to predict temperature and humidity inside
a naturally ventilated greenhouse, achieving error index of 0.0025 for
temperature and 0.0024 for humidity. [18] collected several micro-climate
parameters of a greenhouse including temperature, relative humidity, wind
speed, and solar radiation located in a sub-humid sub-tropical region of
India. They then explored the optimal parameters for artificial neural network
models and made one-day-ahead predictions of air temperature and relative
humidity of a greenhouse and achieved a mean absolute error of 0.558°C and
1.976% in their dataset. [14] constructed a recurrent neural network using
substrate temperature, air humidity, illumination, and concentration for model
training and achieved results with root-mean-square errors of 0.751 and 0.781
for temperature and humidity, respectively. In recent machine learning
approaches, the evaluation of the model has often been limited to the overall
prediction accuracy, and the data from all temperature ranges are treated
equally, so the significant impact of extreme temperature on crops was not
considered. Further, the models in prior works often performed badly in
extreme conditions. The data used to train the models have often been
collected from historical greenhouse data, within which the amount of extreme
data is typically smaller than the amount of normal-temperature data. This
lowers the prediction accuracy at extreme temperatures, increasing the
probability that extreme temperatures may be misjudged as normal temperatures,
with a higher likelihood of severe economic losses. To address the above
problems, in this paper, we propose an improved loss function. In contrast to
the traditional mean squared error, the proposed improved loss function
assigns different weights to extreme and normal temperatures, which can help
models focus more on the extreme cases. The improved loss function calculates
differently for predicted values higher or lower than the actual value.
Therefore, models with this function are updated based on the types of wrong
predictions in extreme temperatures, thereby improving the performance of
models under extreme conditions.

The remainder of this paper is organized as follows. In section 2, the
greenhouse dataset we used and preprocessing steps are introduced, along with
a presentation of the theoretical basis of the work and the proposed
algorithm. In section 3, we describe the experimental setup and discuss the
experimental results. And section 4 provides conclusion for this work.

##  2 Materials and Methods

###  2.1 DATA Description

Our data were collected from the solar greenhouse of the Zhuozhou Farm of
China Agricultural University (39° 28’ 8.19740” N, 115° 51’ 9.31366” E). The
simulated greenhouse is shown in Figure 1. The greenhouse spans 60 m from east
to west and 8 m from north to south and has a cold-formed steel frame
structure. A suite of installed sensors collects environmental data including
temperature, humidity, pressure, illumination, and carbon dioxide
concentration at a frequency of one time per minute.

Figure 1: Layout of the experimental greenhouse. The environmental sensor is
installed at a height of 1m.

We collected approximately 14 months of data from June 17, 2018, to August 12,
2019, with approximately 410,000 valid data records. An example of the data
samples is shown in Table 1.

Table 1: Examples of data samples.
Time | Temperature(°C) | Humidity(%rh) |  Pressure($10^{2}$ Pa) | Illumination(lx) |  $CO_{2}$(ppm)
---|---|---|---|---|---
2018/06/17 19:35:23 | 27.90 | 60.8 | 997.9 | 657.0 | 611.0
2018/06/17 19:36:23 | 27.90 | 60.8 | 997.9 | 657.0 | 611.0
2018/06/17 19:37:22 | 27.90 | 60.8 | 997.9 | 657.0 | 611.0
2018/06/17 19:38:22 | 27.90 | 60.8 | 998.0 | 657.0 | 611.0
2018/06/17 19:39:22 | 27.90 | 60.8 | 998.0 | 657.0 | 611.0
2018/06/17 19:40:22 | 27.90 | 60.8 | 998.0 | 657.0 | 611.0
…… | …… | …… | …… | …… | ……

In Table 1, Temperature (°C) represents the temperature value inside the
greenhouse, and its minimum accuracy was 0.01°C. Humidity (%rh) denotes the
relative humidity value inside the greenhouse, with a minimum accuracy of 0.1%
Pressure ($10^{2}$ Pa), illumination (lx), and CO2 (ppm) denote the air
pressure, light intensity, and carbon dioxide concentration inside the
greenhouse.

###  2.2 Data Preprocessing

We first perform data preprocessing before feeding them into machine learning
algorithms. In this study, data preprocessing can be divided into four steps,
including data cleaning, data aggregation, data reconstruction, and
training/testing sets division. Since greenhouse sensors may be unstable
during data acquisition, anomalous or missing values may exist in raw data.
Therefore, in data cleaning process, we remove data with duration less than 30
minutes, and then remove outliers whose temperature is more than 40°C or below
than 0°C. Compared with the changing rate of the greenhouse environment, the
data collection frequency of one time per minute was relatively rapid. Thus,
it was necessary to reduce the data sampling frequency before training the
model. We average the collected greenhouse data at 10-minutes intervals. This
strategy not only served to reduce the data frequency, but also to balance
numerical fluctuations over short periods potentially caused by unstable
sensors. Since an early warning of extreme temperatures four hours in advance
is sufficient for farmers to perform timely emergency remedial measures, we
used four hours of data to predict temperature four hours in advance. We
reconstruct the data to achieve this goal, as shown in Figure 2.

Figure 2: Data reconstruction. We used 4h of data to predict temperature 4h
later.

Then, to better verify the performance of the model under extreme temperature
conditions, we selected data from August 2018 (summer) and February 2019
(winter) as testing set. The size of the testing set is 6570 and that of the
training set is 36244. After the above process, the changing curve of
greenhouse temperature from 2018 to 2019 is shown in Figure 3. When the
temperature is higher than 30 °C, greenhouse crops are susceptible to heat
damage [21, 22]. Plant chlorophyll can lose activity easily, and the dark
reaction of photosynthesis can also be blocked, causing damage to plants [23].
Therefore, we define the data with temperature higher than 30°C as extremely
high temperatures. Similarly, since crops are easily affected by cold damage
below 10°C [9, 10], we define data with temperatures less than 10°C as
extremely low temperatures.

Figure 3: Variation of solar greenhouse temperatures between 06/17/2018 and
7/17/2019.

###  2.3 MACHINE LEARNING BASICS IN EXTREME TEMPERATURES

The mathematical model of the conventional machine learning algorithm used to
predict temperature can be expressed as the following equation, where $x$ is
the greenhouse environment data, $\theta$ is the parameter of the machine
model, and $\hat{y}$ is the predicted future greenhouse temperature.

$\hat{y}=h(x,\theta)$

To maximization the accuracy of the prediction results of machine learning
algorithms, a loss function must be introduced to measure the deviation of the
predicted values from the ground truth. To obtain prediction results as close
to the ground truth as possible, it was necessary to solve the following
mathematical optimization problem on the training set.

$\operatorname*{arg\,min}_{\theta}\sum J(h(x_{i},\theta),y_{i})$

where $J$ is the loss function and $y_{i}$ is the true value to be predicted
in the training set.

The learning process of the conventional machine learning algorithms is to
solve the above formula iteratively by using specific optimization algorithms.
Commonly used optimization algorithms include Adam [24], RMSprop [25], and
Adagrad [26]. These optimization algorithms are based on the gradient descent
algorithm, so they need to obtain the gradient of the prediction result on the
training set, namely

$\operatorname{Grad}(\theta)=\sum\frac{\partial
J\left(h\left(x_{i},\theta\right),y_{i}\right)}{\partial\theta_{k}}\cdot\overrightarrow{\theta_{k}}$

according to the chain derivation rule

$\frac{\partial
J\left(h\left(x_{i},\theta\right),y_{i}\right)}{\partial\theta_{k}}=\frac{\partial
J\left(h\left(x_{i},\theta\right),y_{i}\right)}{\partial
h(x,\theta)}\cdot\frac{\partial h(x,\theta)}{\partial\theta_{k}}$

Therefore,

$\operatorname{Grad}(\theta)=\frac{\partial
J\left(h\left(x_{i},\theta\right),y_{i}\right)}{\partial
h(x,\theta)}\cdot\sum\frac{\partial
h(x,\theta)}{\partial\theta_{k}}\cdot\overrightarrow{\theta_{k}}$

Therefore, the loss function has a decisive influence on the learning process
of algorithms.

###  2.4 IMPROVED LOSS FUNCTION

Considering the temperature prediction problem, the predicted results can be
classified into four category, as shown in Table 2.

Table 2: Classification of prediction results.
Case | Predicted result | Ground truth
---|---|---
TP | Extreme | Extreme
TN | Normal | Normal
FP | Extreme | Normal
FN | Normal | Extreme

In the false positive (FP) case, the prediction accuracy has little effect on
production, whereas false negative (FN) prediction can potentially cause
significant economic losses in the greenhouse, since it mispredict the extreme
temperature as normal one. However, the most commonly used loss function for
regression problems is the mean square error,

$J_{\mathrm{RSE}}\left(h\left(x_{i},\theta\right),y_{i}\right)=\frac{1}{m}\sum_{i=1}^{m}\left(h\left(x_{i},\theta\right)-\widehat{y}_{i}\right)^{2}$

It can be observed that the weight of each data sample is the same. However,
the number of data samples of extreme cases was lower than that of normal
cases; as a result, the performance of conventional machine learning models
used in prior works may decline in extreme cases. Therefore, we propose an
improved loss function to enhance prediction accuracy for extreme cases and
reduce FN cases. We first increase the weight for extreme temperature samples
during the training phase. The value of the weight increased was calculated by
dividing the number of extreme cases and normal cases, as given below.

$\displaystyle w_{H}$
$\displaystyle=\frac{\left|U_{\text{normal}}\right|}{\left|U_{\text{high}}\right|},$
$\displaystyle w_{L}$
$\displaystyle=\frac{\left|U_{\text{normal}}\right|}{\left|U_{\text{low}}\right|}$

where $U_{low}=\\{(x_{i},y_{i})|y_{i}\leq T_{low},(x_{i},y_{i})\in U\\}$,
$U_{normal}=\\{(x_{i},y_{i})|T_{low}\textless y_{i}\textless
T_{high},(x_{i},y_{i})\in U\\}$, $U_{high}=\\{(x_{i},y_{i})|y_{i}\geq
T_{high},(x_{i},y_{i})\in U\\}$,and $U$ is the training set. In addition, to
reduce the occurrence of FN cases, we defined an extreme case importance
factor $a$. In extreme cases, we increase the weight of data whose predicted
result is within the real value and correspondingly decrease the weight of
other data. Since machine learning algorithms with backpropagation mechanisms
or leaf splitting strategies tend to pay more attention to data with larger
weight, the introduction of the factor $a$ can make the prediction results of
the model in extreme cases be more accurate. In summary, the improved loss
function definition $J(\theta)$ is as follows:

$J(\theta)=J_{\text{high}}(\theta)+J_{\text{low}}(\theta)+J_{\text{normal}}(\theta)$

where

$J_{\text{high}}(\theta)=\sum_{\left(x_{i},y_{i}\right)\in
U_{\text{high}}}w_{H}*f(x)*\left(h_{\theta}\left(x_{i}\right)-y_{i}\right)^{2}$
$J_{\text{low}}(\theta)=\sum_{\left(x_{i},y_{i}\right)\in
U_{\text{low}}}w_{L}*g(x)*\left(h_{\theta}\left(x_{i}\right)-y_{i}\right)^{2}$
$J_{\text{normal}}(\theta)=\sum_{\left(x_{i},y_{i}\right)\in
U_{\text{normal}}}\left(h_{\theta}\left(x_{i}\right)-y_{i}\right)^{2}$

Here, $f$ and $g$ are hinge functions:

$f\left(x_{i}\right)=\left\\{\begin{array}[]{r}a,y_{i}<h_{\theta}\left(x_{i}\right)\\\
1-a,y_{i}\geq
h_{\theta}\left(x_{i}\right)\end{array},\left(x_{i},y_{i}\right)\in
U_{\text{high}}\right.$
$g\left(x_{i}\right)=\left\\{\begin{array}[]{r}1-a,y_{i}<h_{\theta}\left(x_{i}\right)\\\
a,y_{i}\geq
h_{\theta}\left(x_{i}\right)\end{array},\left(x_{i},y_{i}\right)\in
U_{\text{low}}\right.$

And the calculation process can be summarized by the following pseudo code.

Algorithm 1 Improved loss function calculation process

Input: $\hat{y}=\\{\hat{y_{1}},...,\hat{y_{N}}\\}$, $a$, $T_{H}$, $T_{L}$
$\hat{y}$ is the predicted value of temperature $y$ is the true value of
temperature $a$ is the extreme temperature importance factor $T_{H}$ is the
high temperature threshold $T_{L}$ is the low temperature threshold

Output: Total loss $L_{total}$ Initialize $L_{total}=0$

$y_{H}\leftarrow\\{y_{i}|y_{i}\geq T_{H}\\}$

$y_{normal}\leftarrow\\{y_{i}|T_{L}\textless y_{i}\textless T_{H}\\}$

$y_{L}\leftarrow\\{y_{i}|y_{i}\leq T_{L}\\}$

$w_{H}=\nicefrac{{y_{normal}}}{{y_{H}}}$

$w_{L}=\nicefrac{{y_{normal}}}{{y_{L}}}$

for each pair $\hat{y_{i}}$ and $y_{i}$ do:

if $y_{i}\geq T_{H}$:

if $y_{i}\textgreater\hat{y_{i}}$:

$L\leftarrow w_{H}*a*(y_{i}-\hat{y_{i}})^{2}$

else:

$L\leftarrow w_{H}*(1-a)*(y_{i}-\hat{y_{i}})^{2}$

end if

else if $y_{i}\leq T_{L}$:

if $y_{i}\textgreater\ \hat{y_{i}}$:

$L\leftarrow w_{L}*(1-a)*(y_{i}-\hat{y_{i}})^{2}$

else:

$L\leftarrow w_{L}*a*(y_{i}-\hat{y_{i}})^{2}$

end if

else:

$L\leftarrow(y_{i}-\hat{y_{i}})^{2}$

end if

$L_{total}\leftarrow L_{total}+L$

end for

return $L_{total}$

##  3 EXPERIMENTS AND RESULTS

###  3.1 EXPERIMENTAL SETUP

We conduct experiments on our preprocessed dataset. Data from August 2018 and
February 2019, the number of which is 6,570, are used for testing, and the
remaining 36244 data are used for training. And they are conducted on a
computer workstation using the Windows 10 operating system with an Intel Core
i9-9900K CPU at 3.60 GHz, with 16 GB of memory installed as well as an Nvidia
GeForce RTX 2080 Ti graphics processing unit (11GB), and all models are
implemented using the Python3 programming language. To verify the versatility
of the improved loss function, we conduct comparative experiments on the basic
models of LightGBM (LGB) [27], long short-term memory (LSTM) [28], and
artificial neural network (ANN) [29], and evaluate the models after applying
the improved loss function. Furthermore, to further verify the optimization
effect of the improved loss function, we select a = 0.5, 0.7, and 0.9 for
comparison and compare the basic LGB model to the LGB with improved loss
function under extreme temperature conditions. When comparing models with
different extreme case important factors, we use the same hyperparameters. The
BPNN’s hyperparameters are shown in Table 3. We used PyTorch to build the BPNN
model. We used a three-layer neural network (excepting the input layer), and
the number of nodes was the same as [17].

Table 3: Settings of BPNN hyperparameters.
Algorithm | BPNN
---|---
Number of hidden layers | 3
Number of nodes in each hidden layer | 6, 3, 1
Optimizer | Adam
Learning rate | 0.001
Batch size | 2000
Epoch | 100

The hyperparameter settings used for the LSTM are shown in Table 4. We also
use PyTorch to build the LSTM model.

Table 4: Settings of LSTM hyperparameters.
Algorithm | LSTM
---|---
Number of LSTM layers | 1
Number of nodes in LSTM layer | 48
Optimizer | Adam
Learning rate | 0.02
Batch size | 2000
Epoch | 100

The key hyperparameter settings used for the LGB model are shown in Table 5.
We use python library LightGBM to implement LGB model.

Table 5: Settings of LSTM hyperparameters.
Algorithm | LGB
---|---
Boosting type | Gbdt
Learning rate | 0.05
Epoch | 100

###  3.2 MODEL PERFORMANCE EVALUATION

To effectively evaluate the predictive performance of the model, the
evaluation indicators used in this paper include mean absolute error (MAE),
mean bias and recall rate. The evaluation index calculation formula is as
follows:

$\text{MAE}=\frac{1}{m}\sum_{i=1}^{m}\left|y_{i}-\hat{y}_{l}\right|$
$\text{Mean bias}=\frac{1}{m}\sum_{i=1}^{m}\left(y_{i}-\widehat{y}_{l}\right)$
$\text{Recall}=\frac{TP}{TP+FN}$

The definitions of TP, FP, and FN are the same as those given in Table 2. The
higher the recall rate, the lower the FN, the more extreme cases were
correctly predicted. To test the effect of improved loss function, we first
examine the performance of the three basic models, then examine the
performance of a model that only increase the weight of extreme cases
($a=0.5$), and finally test the performance of a model with $a=0.9$.

Table 6: Predicted results in test set.
Model | Improved loss function | Recall | MAE (high) | MAE (low) | MAE (normal)
---|---|---|---|---|---
LGB | ✘ | 0.82 | 2.15 | 0.79 | 0.67
a=0.5 | 0.87 | 1.71 | 0.64 | 0.78
a=0.9 | 0.87 | 1.56 | 0.59 | 0.85
LSTM | ✘ | 0.68 | 2.31 | 1.90 | 1.06
a=0.5 | 0.75 | 1.61 | 1.53 | 1.57
a=0.9 | 0.78 | 1.60 | 1.33 | 1.96
ANN | ✘ | 0.78 | 1.98 | 1.11 | 0.98
a=0.5 | 0.82 | 1.72 | 0.84 | 1.62
a=0.9 | 0.83 | 1.44 | 1.11 | 1.91

As shown in Table 6, the MAE of the basic models in extreme cases is higher
than the normal case, which demonstrates that models do have performance
degradation in extreme temperature. With the introduction of the high- and
low-temperature weighting factors $w_{H}$ and $w_{L}$, all three models
($a=0.5$) obtain a lower MAE than the basic models for both the high and low
temperatures. For LSTM, when $a=0.5$, the MAE value at high temperature
decrease from 2.31 to 1.61, and the MAE value at low temperature decrease from
1.90 to 1.53. The recall rate of the models also improve, which shows that the
proposed method reduce the FN cases successfully. Although the performance of
the model under normal temperature reduce a little, the growth of the plants
is almost unaffected when temperature fluctuate in the normal range. When
$a=0.9$, the results of three models have been further improved, which shows
the function $f$ and $g$ perform well in all models, which can reduce the FN
cases and reduce the MAE under extreme temperatures. To further determine the
influence of the extreme case importance factor on the model, we take LGB
model as an example and further compare the prediction results with different
$a$. The results are shown in Table 7.

Table 7: LGB comparison results. We calculate the mean bias separately on high
temperature and low temperature. When the mean bias is less than zero, the
predicted value is less than the ground truth, and vice versa.
Improved loss function | Mean bias (high) | Mean bias (low)
---|---|---
✘ | -1.95 | 0.72
$a=0.5$ | -1.39 | 0.45
$a=0.7$ | -1.26 | 0.33
$a=0.9$ | -1.08 | 0.28

The results in Table 7 show that when the extreme case importance factor a
change from 0.5 to 0.9, the performance of the improved LGB under extreme
temperatures gradually increase. When $a=0.9$, the mean bias under high-
temperature conditions increase from ground truth-1.95°C to ground
truth-1.08°C, which means the possibility of incorrectly predicting high
temperature as normal temperature is reduced. Similarly, under low-temperature
conditions, the mean bias is reduced from ground truth+0.72°C to ground
truth+0.28°C, indicating that the improved loss function also function well on
low temperature.

Figure 4: Comparison of prediction results at high temperatures (a) and low
temperatures (b).

To show the improvement in the prediction results, we use a box plot to
compare the basic model with the improved one and visually show the
temperature difference distribution, as exhibited in Figure 4. Under high
temperature, the median value of temperature difference increases from $-0.26$
to $0.13$, the mean value increases from $-0.32$ to $0.19$, and the standard
deviation reduces from $0.9$ to $0.86$. Under low temperature, the mean value
decrease from $0.28$ to $-0.12$, the median value reduces from $0.20$ to
$-0.15$, and the standard deviation reduces from $0.57$ to $0.49$. This means
that our approach can reduce the standard deviation and can shift the
prediction results towards extreme temperatures, reducing the possibility of
missed prediction in extreme cases.

##  4 Conclusions

In order to forecast the freezing or heat damage that crops may suffer, it is
essential to predict extreme temperatures accurately and timely. In this
paper, an improved loss function that can be applied to various machine
learning algorithms was proposed. The function mainly includes two parts,
which are extreme weight and adjustment function. Then, we collected 14 months
of various greenhouse sensor data from Zhuozhou Farm’s solar greenhouse and
perform data cleaning, data aggregation, data reconstruction, and
training/testing sets division sequentially. By implementing the improved loss
function to LGB, LSTM, and ANN, the results show that compared with the basic
models, the performance of the improved model becomes better in extreme cases.
Specifically, the recall rate of the prediction results under extreme
conditions increases significantly and the MAE value decreases at the same
time. It proves that the proposed approach can better guarantee timely
judgment of extreme temperatures in greenhouses, thereby preventing
unnecessary losses caused by incorrect and missed predictions.

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

  * [1] Kang, M., Fan, X.-R., Hua, J., Wang, H., Wang, X., and Wang, F.-Y. Managing traditional solar greenhouse with cpss: A just-for-fit philosophy.  IEEE Transactions on cybernetics, 48, 3371–3380.
  * [2] Hatfield, J. L., Boote, K. J., Kimball, B. A., Ziska, L. A., Izaurralde, R. C., Ort, D. R., Thomson, A. M., and Wolfe, D. W. Climate impacts on agriculture: Implications for crop production.  Agronomy Journal, 103, 351–370.
  * [3] Luo, Q. Temperature thresholds and crop production: a review.  Climatic change, 109, 583–598.
  * [4] Kaushal, N., Bhandari, K., Siddique, K. H. M., and Nayyar, H. Food crops face rising temperatures: An overview of responses, adaptive mechanisms, and approaches to improve heat tolerance.  Cogent Food & Agriculture, 2.
  * [5] Wang, W., Vinocur, B., and Altman, A. Plant responses to drought, salinity and extreme temperatures: towards genetic engineering for stress tolerance.  Planta, 218, 1–14.
  * [6] Hatfield, J. L. and Prueger, J. H. Temperature extremes: Effect on plant growth and development.  Weather and climate extremes, 10, 4–10.
  * [7] Ferris, R., Ellis, R. H., Wheeler, T. R., and Hadley, P. Effect of high temperature stress at anthesis on grain yield and biomass of field-grown crops of wheat.  Annals of Botany, 82, 631–639.
  * [8] Siebert, S. and Ewert, F. Future crop production threatened by extreme heat.  Environmental Research Letters, 9, 041001.
  * [9] Huner, N. P., Öquist, G., Hurry, V. M., Krol, M., Falk, S., and Griffith, M. Photosynthesis, photoinhibition and low temperature acclimation in cold tolerant plants.  Photosynthesis research, 37, 19–39.
  * [10] Liu, W., Yu, K., He, T., Li, F., Zhang, D., and Liu, J. The low temperature induced physiological responses of avena nuda l., a cold-tolerant plant species.  The Scientific World Journal, 2013.
  * [11] Tong, G., Christopher, D. M., and Li, B. Numerical modelling of temperature variations in a chinese solar greenhouse.  Computers and Electronics in Agriculture, 68, 129–139.
  * [12] Vanthoor, B., Stanghellini, C., van Henten, E. J., and de Visser, P. A methodology for model-based greenhouse design: Part 1, a greenhouse climate model for a broad range of designs and climates.  Biosystems Engineering, 110, 363–377.
  * [13] Yu, H., Chen, Y., Hassan, S. G., and Li, D. Prediction of the temperature in a chinese solar greenhouse based on lssvm optimized by improved pso.  Comput. Electron. Agric., 122, 94–102.
  * [14] Hongkang, W., Li, L., Yong, W., Fanjia, M., Haihua, W., and Sigrimis, N. Recurrent neural network model for prediction of microclimate in solar greenhouse.  IFAC-PapersOnLine, 51, 790–795.
  * [15] Luo, C., Song, X., and Zheng, Y. A novel forecasting model for the long-term fluctuation of time series based on polar fuzzy information granules.  Inf. Sci., 512, 760–779.
  * [16] Lins, I. D., Araujo, M., das Chagas Moura, M., Silva, M. A., and Droguett, E. L. Prediction of sea surface temperature in the tropical atlantic by support vector machines.  Comput. Stat. Data Anal., 61, 187–198.
  * [17] Ahmad, R. B., Lazin, N. M., and Samsuri, S. F. M. Neural network modeling and identification of naturally ventilated tropical greenhouse climates.  WSEAS Transactions on Systems and Control archive, 9, 445–453.
  * [18] Singh, V. Prediction of greenhouse micro-climate using artificial neural network.  Applied Ecology and Environmental Research, 15, 767–778.
  * [19] Nayak, S. and Tiwari, G. N. Energy and exergy analysis of photovoltaic/thermal integrated with a solar greenhouse.  Energy and Buildings, 40, 2015–2021.
  * [20] Moré, J. J. The levenberg-marquardt algorithm: implementation and theory.  Numerical analysis, pp. 105–116. Springer.
  * [21] Bita, C. E. and Gerats, T. Plant tolerance to high temperature in a changing environment: scientific fundamentals and production of heat stress-tolerant crops.  Frontiers in Plant Science, 4.
  * [22] Chen, H.-H., Shen, Z.-Y., and Li, P. H. Adaptability of crop plants to high temperatures stress1.  Crop Science, 22, 719–725.
  * [23] Fahad, S., et al. Crop production under drought and heat stress: Plant responses and management options.  Frontiers in Plant Science, 8.
  * [24] Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization.  CoRR, abs/1412.6980.
  * [25] Tieleman, T., Hinton, G., et al. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude.  COURSERA: Neural networks for machine learning, 4, 26–31.
  * [26] Duchi, J. C., Hazan, E., and Singer, Y. Adaptive subgradient methods for online learning and stochastic optimization.  J. Mach. Learn. Res.
  * [27] Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., and Liu, T.-Y. Lightgbm: A highly efficient gradient boosting decision tree.  Advances in neural information processing systems, 30, 3146–3154.
  * [28] Hochreiter, S. and Schmidhuber, J. Long short-term memory.  Neural computation, 9, 1735–1780.
  * [29] Jain, A. K., Mao, J., and Mohiuddin, K. M. Artificial neural networks: A tutorial.  Computer, 29, 31–44.