# Efficient Hamiltonian reconstruction from equilibrium states

Adam Artymowicz
California Institute of Technology, Pasadena, CA 91125, USA

###### Abstract

We describe a novel algorithm that recovers the Hamiltonian and temperature of
a quantum state in thermal equilibrium using a restricted set of measurements.
The algorithm works by imposing a semidefinite constraint based on free energy
minimization, which we interpret in terms of thermodynamic stability. We
benchmark the algorithm on the problem of learning a nearest-neighbour
Hamiltonian on a 100-qubit spin chain.

##  1 Introduction

A defining feature of the Gibbs state
$\rho=e^{-h/T}/\operatorname{tr}(e^{-h/T})$ is that it minimizes the free
energy:

$\displaystyle\rho=\operatorname*{\arg\\!\min}F$ (1)

where $F(\sigma)=-TS(\sigma)+\operatorname{tr}(\sigma h)$ and $S$ is the von
Neumann entropy. Conversely, given $\rho$, the condition (1) uniquely
specifies $h$, and this fact can in principle be used to recover the
Hamiltonian of a Gibbs state [1]. Thanks to the strict convexity of $F$, (1)
is equivalent to the first-order condition

$\displaystyle dF[\delta\rho]\geq 0$ (2)

for any $\delta\rho$ in the tangent space at $\rho$. However, not only is
computing the derivative of $F$ difficult, (2) also involves imposing a number
of conditions that scales quadratically with the dimension of the physical
Hilbert space.

This work is based on a hierarchy of relaxations of the local minimality
condition (2). These relaxations recently appeared in [2] and use a matrix
version of the Araki-Sewell inequality [3]. For each relaxation we give a
classical algorithm that attempts to reconstruct the Hamiltonian $H$ and the
temperature $T$ from the expectation values of certain operators. Given a set
of variational Hamiltonian terms $h_{1},\ldots h_{s}$, the algorithm either
returns a candidate Hamiltonian $h\in\operatorname{span}(h_{1},\ldots,h_{s})$
and candidate temperature $T$, or returns a certificate that the state is not
a Gibbs state of any Hamiltonian in the span of $h_{1},\ldots,h_{s}$.

To analyze the algorithm, we introduce a condition which we call restricted
thermodynamic stability (RTS) that generalizes the Gibbs variational condition
(2). Indeed, it is equivalent to (2) where $\delta\rho$ is restricted to a
certain conic region in the tangent space of $\rho$. Physically, the Gibbs
variational condition is related to stability of the state against
perturbations. Correspondingly, we interpret the RTS condition as stability
with respect to a certain class of open system dynamics.

Although we prove correctness of the algorithm in the absence of measurement
noise (Corollary 1), we do not give any bounds on sample complexity (number of
copies of $\rho$ needed for accurate reconstruction) or classical
computational complexity, leaving this to future work. Instead, we include
numerical evidence that the algorithm achieves good reconstruction accuracy
using modest computational resources and in the presence of measurement noise.

Many other proposals for Hamiltonian reconstruction from Gibbs states [4, 5,
6, 7, 8, 9, 10] have recently been advanced. Of these, [7, 8, 9, 10] include
numerical tests. In comparison to these approaches, algorithm described in the
current work proves to be highly scalable while being robust to noise. Some
other advantages are its generality (the choice of variational Hamiltonian
terms $h_{1},\ldots,h_{s}$ is unrestricted), ease of implementation, and the
fact that it is physically motivated and involves few hyperparameters.

We begin in Section 2 with a description of the algorithm and a guide to
interpreting its output. Section 3 concerns a theoretical analysis of the
algorithm, including proof of correctness. In Section 4 we describe the
results of numerical simulations on a 100-qubit spin chain. We conclude in
Section 5 with some future directions for research.

The python implementation of the learning algorithm used in this work, as well
as all routines used for testing it, are available for use at

https://github.com/artymowicz/hamiltonian-learning

In the remainder of this section we establish some notation and recall some
basic facts about Gibbs states. Let $\mathcal{H}$ be the physical Hilbert
space of a system, and assume that $\dim\mathcal{H}<\infty$. Write
$\mathcal{A}$ for the algebra of all linear operators on $\mathcal{H}$. As a
rule, we will use lowercase letters to denote elements of $\mathcal{A}$, ie.
operators acting on the physical Hilbert space (this will be to differentiate
them from operators on the GNS Hilbert space which we introduce in later
sections). We will denote the adjoint of an operator $a\in\mathcal{A}$ by
$a^{*}$.

A state on $\mathcal{A}$ is a linear map
$\omega:\mathcal{A}\to\operatorname{\mathbb{C}}$ satisfying $\omega(1)=1$ and
$\omega(a^{*}a)\geq 0$ for all $a\in\mathcal{A}$. These are in one-to-one
correspondence with density matrices $\rho$ via
$\omega(a)=\operatorname{tr}(\rho a)$. A state $\omega$ is called faithful if
$\omega(a^{*}a)=0$ implies $a=0$, or equivalently, if its density matrix is
invertible. Given a selfadjoint operator $h\in\mathcal{A}$, we say $h$ is a
symmetry of $\omega$ if $\omega([h,a])=0$ for all $a\in\mathcal{A}$, or
equivalently if $[\rho,h]=0$. Given a positive number $T$, the Gibbs state
corresponding to $h$ at temperature $T$ is the state given by the density
matrix $\rho=e^{-h/T}/\operatorname{tr}(e^{-h/T})$. A Gibbs state is always
faithful, and any self-adjoint operator $h$ is a symmetry of its own Gibbs
state.

##  2 Description of the algorithm

We begin now with a description of the algorithm. Let $\omega$ be a faithful
state of $\mathcal{A}$. The algorithm requires the following input:

  1. 1.

A set of selfadjoint traceless operators $h_{1},\ldots,h_{s}\in\mathcal{A}$,

  2. 2.

A choice of linearly independent operators $b_{1}\ldots b_{r}\in\mathcal{A}$
such that
$\operatorname{span}(b_{1},\ldots,b_{r})=\operatorname{span}(b_{1}^{*},\ldots,b_{r}^{*})$,

  3. 3.

All expectation values of the form $\omega(b_{i}^{*}b_{j})$ for $1\leq i,j\leq
r$ and $\omega(b_{i}^{*}[h_{\alpha},b_{j}])$ for $1\leq i,j\leq r$ and
$1\leq\alpha\leq s$.

The operators $h_{1},\ldots,h_{r}$ will be the variational Hamiltonian terms.
We call the operators $b_{1}\ldots b_{r}$ the _perturbing operators_ –
physically, these determine the class of perturbations with respect to which
we will enforce stability (see section 3). A reasonable choice for spin
systems is the set of all geometrically $k$-local Pauli operators for some
$k>0$.

### Step 1

The first step is to orthonormalize the $b_{i}$ in an appropriate sense. Let
$a_{1},\ldots,a_{r}\in\operatorname{span}(b_{1},\ldots,b_{r})$ satisfy
$\omega(a_{i}^{*}a_{j})=\delta_{ij}$. Such a set can be found by diagonalizing
the quadratic form $\omega(b_{i}^{*}b_{j})$ (which is necessarily
nondegenerate because $\omega$ is faithful).

### Step 2

For each $\alpha=1,\ldots,s$ let $\boldsymbol{H}_{\alpha}$ to be the $r\times
r$ matrix

$\displaystyle(\boldsymbol{H}_{\alpha})_{ij}$
$\displaystyle:=\omega(a_{i}^{*}[h_{\alpha},a_{j}]),$ (3)

and define the $s\times s$ matrix

$\displaystyle
W_{\alpha\beta}:=\operatorname{tr}((\boldsymbol{H}_{\alpha}^{\dagger}-\boldsymbol{H}_{\alpha})(\boldsymbol{H}_{\beta}-\boldsymbol{H}_{\beta}^{\dagger})).$
(4)

If $W$ is invertible, then the algorithm terminates. Otherwise, let
$\tilde{h}_{1},\ldots,\tilde{h}_{q}$ be a basis for the kernel of $W$, and let
$(\widetilde{\boldsymbol{H}}_{\alpha})_{ij}=\omega(a_{i}^{*}[\tilde{h}_{\alpha},a_{j}])$.

### Step 3

Define the following $r\times r$ matrix:

$\displaystyle\boldsymbol{\Delta}_{ij}:=\omega(a_{j}a_{i}^{*}),$ (5)

and solve the semidefinite optimization problem:

$\displaystyle\underset{y\in\operatorname{\mathbb{R}}^{q},T\in\operatorname{\mathbb{R}}_{\geq
0}\mu\in\operatorname{\mathbb{R}}}{\text{maximize}}$ $\displaystyle\mu$ (6)
subject to $\displaystyle
T\log(\boldsymbol{\Delta})+\sum_{\alpha=1}^{q}y_{\alpha}\widetilde{\boldsymbol{H}}_{\alpha}-\mu
I\succeq 0,$ (7)
$\displaystyle\sum_{\alpha=1}^{q}y_{\alpha}\omega(\tilde{h}_{\alpha})=-1.$ (8)

Here (7) is the matrix Araki-Sewell inequality (33) with a regularization
parameter, and (8) is a normalization that removes the multiplicative gauge of
the Hamiltonian111A similar, but not equivalent, algorithm could be obtained
by removing the normalization condition (8) and setting $T=1$..

Let $y^{*},\mu^{*},T^{*}$ be the optimal values of the semidefinite program.
Then the algorithm outputs the candidate Hamiltonian
$h^{*}:=\sum_{\alpha}y_{\alpha}^{*}\tilde{h}_{\alpha}$, candidate temperature
$T^{*}$ and the parameter $\mu^{*}$ which is used in the interpretation of the
ouput.

### Interpretation

Below is a summary of the possible outputs of the algorithm and their
interpretations.

  1. 1.

In step 2, if $W$ is invertible then $\omega$ is not a stationary state of any
operator in the span of $h_{1},\ldots,h_{s}$.

  2. 2.

In step 3:

    1. (a)

If the optimization (6) terminates with $\mu^{*}<0$ then $\omega$ is not a
Gibbs state of any Hamiltonian in the span of $h_{1},\ldots,h_{s}$222However,
numerics suggest that if $\mu^{*}<0$ but the magnitude of $\mu^{*}$ is small,
then $\omega$ is close to being a Gibbs state of $\sum
y_{\alpha}^{*}\tilde{h}_{\alpha}$ (see Section 4)..

    2. (b)

Otherwise, $\sum y_{\alpha}^{*}\tilde{h}_{\alpha}$ and $T^{*}$ are candidates
for the Hamiltonian and temperature of the state $\omega$.

##  3 Restricted thermodynamic stability

The theoretical justification of the algorithm will be based on two notions.
The first is a relaxation of the notion of a symmetry of a state. For a
collection of operators $b_{1},\ldots,b_{r}\in\mathcal{A}$, we say a self-
adjoint operator $h\in\mathcal{A}$ is a quasi-symmetry of $\omega$ with
respect to $b_{1},\ldots,b_{r}$ if $\omega([b^{*}b,h])=0$ for all
$b\in\operatorname{span}\\{b_{1},\ldots,b_{r}\\}$. Any symmetry of $\omega$ is
a quasi-symmetry.

The second notion is a relaxation of the notion of a Gibbs state. It is
inspired by local thermodynamic stability, which was introduced by Araki and
Sewell in [3] and generalizes the usual notion of thermal equilibrium. Namely,
Araki and Sewell considered states in infinite volume for which any local
perturbation of the state increases the free energy. They showed that these
states are characterized by a correlation inequality which is alternatively
known as the Araki-Sewell, Roestropff-Araki-Sewell, or energy-entropy balance
inequality. Below we introduce a variant of the local thermodynamic stability
condition (Definition 1), and prove that it implies a matrix version of the
Araki-Sewell inequality (Proposition 3). This matrix version of the Araki-
Sewell inequality first appeared in [2], although a similar inequality can be
found in the 1985 article [11].

Let $b_{1},\ldots,b_{r}\in\mathcal{A}$ be any collection of operators (not
necessarily self-adjoint). Given an anti-Hermitean matrix $\boldsymbol{M}$ and
a positive semidefinite matrix $\boldsymbol{\Lambda}$, the Lindbladian
superoperator
$L_{\boldsymbol{M},\boldsymbol{\Lambda}}:\mathcal{A}\to\mathcal{A}$ is a
linear map defined as

$\displaystyle
L_{\boldsymbol{M},\boldsymbol{\Lambda}}(a):=\sum_{ij}\left\\{-\frac{1}{2}\boldsymbol{M}_{ij}[b_{i}^{*}b_{j},a]+\boldsymbol{\Lambda}_{ij}(b_{i}^{*}ab_{j}-\frac{1}{2}(b_{i}^{*}b_{j}a+ab_{i}^{*}b_{j}))\right\\},\
\ a\in\mathcal{A}.$ (9)

It generates a time-evolution
$e^{tL_{\boldsymbol{M},\boldsymbol{\Lambda}}}:\mathcal{A}\to\mathcal{A}$ which
describes open dynamics of the system weakly coupled to its environment[12,
13]. The matrics $\boldsymbol{M}_{ij}$ and $\boldsymbol{\Lambda}_{ij}$
describe the internal couplings and couplitngs to the enviromnent,
respectively. Given a self-adjoint $h\in\mathcal{A}$ and a temperature $T>0$
the free energy of a state $\omega$ is

$\displaystyle F(\omega)=-TS(\omega)+\omega(h)$ (10)

where $S(\omega)$ is the von Neumann entropy (ie.
$-\operatorname{tr}(\rho\log\rho)$ of the corresponding density matrix).

######  Definition 1.

Let $\omega$ be a faithful state and $h$ a selfadjoint operator. We say the
pair ($\omega$,$h$) satisfies _restricted thermodynamic stability (RTS)_ at
the temperature $T>0$ with respect to the operators $b_{1},\ldots,b_{r}$ if
for any Lindbladian $L_{\boldsymbol{M},\boldsymbol{\Lambda}}$ of the form (9)
we have

$\displaystyle\frac{d}{dt}\bigg{|}_{t=0}F(\omega_{t})\geq 0,$ (11)

where $\omega_{t}=\omega\circ e^{tL_{\boldsymbol{M},\boldsymbol{\Lambda}}}$.

We now give a characterization of the RTS condition in terms of the Gelfand-
Naimark-Segal construction (see Appendix A for an exposition). Let
$\mathcal{H}_{\omega}$ be the GNS Hilbert space of $\omega$, and $\pi_{\ell}$
and $\pi_{r}$ the left- and right- representations of $\mathcal{A}$ on
$\mathcal{H}_{\omega}$. For any self-adjoint $h\in\mathcal{A}$ we write $H$
for the corresponding GNS Hamiltonian:

$\displaystyle
H:=\pi_{\ell}(h)-\pi_{r}(h)\in\mathcal{B}(\mathcal{H}_{\omega}).$ (12)

Note that even though $h$ is self-adjoint, $H$ need not be. In fact one can
check that $H$ is self-adjoint if and only if $h$ is a symmetry of $\omega$.

The first-order change in the expectation value of $h$ under the Lindblad
evolution (9) has a straightforward expression in terms of the matrix elements
of $H$. With $\omega_{t}=\omega\circ
e^{tL_{\boldsymbol{M},\boldsymbol{\Lambda}}}$, we have

$\displaystyle\frac{d}{dt}\bigg{|}_{t=0}\omega_{t}(h)$
$\displaystyle=\omega(L_{\boldsymbol{M},\boldsymbol{\Lambda}}(h))$ (13)
$\displaystyle=\frac{1}{2}\sum_{ij}\boldsymbol{M}_{ij}\omega(b_{i}^{*}[h,b_{j}]-[b_{i}^{*},h]b_{j})+\boldsymbol{\Lambda}_{ij}\omega(b_{i}^{*}[h,b_{j}]+[b_{i}^{*},h]b_{j})$
(14) $\displaystyle=\frac{1}{2}\sum_{ij}\boldsymbol{M}_{ij}\langle
b_{i}|H-H^{\dagger}|b_{j}\rangle+\boldsymbol{\Lambda}_{ij}\langle
b_{i}|H+H^{\dagger}|b_{j}\rangle.$ (15)

Let $\Delta$ be the modular operator of $\omega$. Analogously to (15) above,
the matrix elements of $\log(\Delta)$ give the first-order changes of the
entropy of $\omega$ under Lindbladian evolution, as we now show. Using a prime
to indicate time-derivative at $t=0$, we have

$\displaystyle-\operatorname{tr}(\rho\log(\rho))^{\prime}$
$\displaystyle=-\operatorname{tr}(\rho^{\prime}\log(\rho))-\operatorname{tr}(\rho\log(\rho)^{\prime}).$
(16)

Using the power series of $\log$ about the identity operator and the cyclicity
of the trace, the second term can be seen to equal
$-\operatorname{tr}(\rho^{\prime})=0$. The first term, meanwhile, equals

$\displaystyle-\operatorname{tr}(\rho^{\prime}\log(\rho))$
$\displaystyle=-\operatorname{tr}(\rho
L_{\boldsymbol{M},\boldsymbol{\Lambda}}(\log(\rho)))$ (17)
$\displaystyle=-\omega(L_{\boldsymbol{M},\boldsymbol{\Lambda}}(\log\rho))$
(18) $\displaystyle=\boldsymbol{M}_{ij}\langle
b_{i}|\frac{\log(\Delta)-\log(\Delta)^{\dagger}}{2}|b_{j}\rangle+\boldsymbol{\Lambda}_{ij}\langle
b_{i}|\frac{\log(\Delta)+\log(\Delta)^{\dagger}}{2}|b_{j}\rangle,$ (19)

where in the last line we used (15) and the expression

$\displaystyle\log(\Delta)=\pi_{\ell}(\log(\rho))-\pi_{r}(\log(\rho))$ (20)

from Appendix A. Since $\rho$ and $\log(\rho)$ commute, $\log(\Delta)$ is
Hermitean and the first term in (19) vanishes, leaving us with

$\displaystyle\frac{d}{dt}\bigg{|}_{t=0}S(\omega_{t})$
$\displaystyle=-\sum_{i,j}\boldsymbol{\Lambda}_{ij}\langle
b_{i}|\log(\Delta)|b_{j}\rangle.$ (21)

Equations (15) and (21) allow us to characterize the RTS and quasi-symmetry
conditions as follows:

######  Proposition 1.

Let $\omega$ be a faithful state, $h\in\mathcal{A}$ a selfadjoint operator,
and $H$ given by (12). Let $b_{1},\ldots,b_{r}\in\mathcal{A}$ a collection of
operators and write
$P:\mathcal{H}_{\omega}\to\operatorname{span}(|b_{1}\rangle,\ldots,|b_{r}\rangle)$
for the orthogonal projection.

  1. i)

$h$ is a quasi-symmetry of $\omega$ with respect to the operators
$b_{1},\ldots,b_{r}$ if and only if $PHP^{\dagger}$ is self-adjoint.

  2. ii)

The pair $(\omega,h)$ has RTS at temperature $T>0$ with respect to the
operators $b_{1},\ldots,b_{r}$ if and only if

$\displaystyle P(T\log\Delta+H)P^{\dagger}\geq 0.$ (22)

###### Proof.

$i)$ $PHP^{\dagger}$ is self-adjoint if and only if

$\displaystyle 0$ $\displaystyle=\langle b|H-H^{\dagger}|b\rangle$ (23)
$\displaystyle=\omega(b^{*}[H,b])-\omega([b^{*},H]b)$ (24)
$\displaystyle=2\omega([b^{*}b,h])$ (25)

for any $b\in\operatorname{span}\\{b_{1},\ldots,b_{r}\\}$.

$ii)$ Suppose $P(T\log\Delta+H)P^{\dagger}\geq 0$. Since $\log\rho$ commutes
with $\rho$, $\log\Delta$ is hermitian, and so $PHP^{\dagger}$ hermitian too.
It follows that for any anti-hermitian $\boldsymbol{M}$ and any positive-
semidefinite $\boldsymbol{\Lambda}$ we have

$\displaystyle\frac{d}{dt}\bigg{|}_{t=0}F(\omega\circ
e^{tL_{\boldsymbol{M},\boldsymbol{\Lambda}}})=\sum_{i,j}\boldsymbol{\Lambda}_{ij}\langle
b_{i}|T\log(\Delta)+H|b_{j}\rangle\geq 0.$ (26)

Conversely, suppose $\omega$ satisfies RTS. Setting $\boldsymbol{M}=0$ in (9)
gives

$\displaystyle 0\leq\frac{d}{dt}\bigg{|}_{t=0}F(\omega\circ
e^{tL_{\boldsymbol{\Lambda}}})=\sum_{i,j}\boldsymbol{\Lambda}_{ij}\langle
b_{i}|T\log(\Delta)+H|b_{j}\rangle.$ (27)

Since this holds for any positive-semidefinite $r\times r$ matrix
$\boldsymbol{\Lambda}$, it follows that $P(T\log(\Delta)+H)P^{\dagger}\geq 0$.
∎

We can in short order deduce from Proposition 1 that the RTS condition depends
only on the span of $|b_{1}\rangle,\ldots,|b_{r}\rangle$, and that if
$(\omega,h)$ satisfy RTS then $h$ is necessarily a quasi-symmetry of $\omega$.
The next Proposition relates the RTS condition to the Gibbs condition.

######  Proposition 2.

Let $b_{1},\ldots,b_{r}\in\mathcal{A}$ be a collection of operators,
$h\in\mathcal{A}$ a selfadjoint operator, and $T>0$. If $\omega$ is the Gibbs
state of $h$ at temperature $T$, then the pair $(\omega,h)$ has RTS at
temperature $T$ with respect to the operators $b_{1},\ldots,b_{r}$. The
converse holds if $\operatorname{span}(b_{1},\ldots,b_{r})=\mathcal{A}$.

###### Proof.

If $\omega$ is the Gibbs state of $h$ at temperature $T$ then it is easy to
check using the explicit expression (20) that $T\log(\Delta)+H=0$, so thanks
to part $ii)$ of Proposition 1, $\omega$ automatically satisfies RTS for any
set of perturbing operators.

For the second statement, suppose that $\omega$ and $h$ satisfy RTS with
respect to $b_{1},\ldots,b_{r}$ and that
$\operatorname{span}(b_{1},\ldots,b_{r})=\mathcal{A}$. By part $ii)$ of
Proposition 1 we have $T\log\Delta+H\geq 0$. We will show that this implies
$T\log\Delta+H=0$. Let $J$ be the modular involution, which is defined as

$\displaystyle J|a\rangle:=|\rho^{1/2}a^{*}\rho^{-1/2}\rangle.$ (28)

For any $a\in\mathcal{A}$ we have

$\displaystyle\langle Ja|H|Ja\rangle$
$\displaystyle=\operatorname{tr}(\rho\,\rho^{-1/2}a\rho^{1/2}[h,\rho^{1/2}a^{*}\rho^{-1/2}])$
(29)
$\displaystyle=\operatorname{tr}(a\rho^{1/2}h\rho^{1/2}a^{*})-\operatorname{tr}(\rho^{1/2}a\rho
a^{*}\rho^{-1/2}h)$ (30) $\displaystyle=-\omega(a^{*}[h,a])$ (31)
$\displaystyle=-\langle a|H|a\rangle,$ (32)

where in the third line we used the fact that $[\rho,h]=0$. The same
calculation with $\log(\rho)$ replacing $h$ shows that $\langle
Ja|\log(\Delta)|Ja\rangle=-\langle a|\log(\Delta)|a\rangle$. It follows that
for every eigenvalue $\lambda$ of $T\log(\Delta)+H$, -$\lambda$ is also an
eigenvalue. Combined with the fact that $T\log\Delta+H\geq 0$, we conclude
that $T\log\Delta+H=0$, as claimed. From this we see that $\log(\rho)-h/T$ is
in the center of $\mathcal{A}$, which means it is a multiple of the identity,
and so $\rho=e^{-h/T}/\operatorname{tr}(e^{-h/T})$. ∎

Now we come to the matrix analog of the Araki-Sewell inequality.

######  Proposition 3.

If $(\omega,h)$ satisfies RTS at $T>0$ then

$\displaystyle T\log(\boldsymbol{\Delta})+\boldsymbol{H}\geq 0.$ (33)

The converse holds if
$\operatorname{span}(|b_{1}\rangle,\ldots,|b_{r}\rangle)$ is invariant under
the modular flow of $\omega$.

###### Proof.

The matrix $\log(\boldsymbol{\Delta})$ is related to $\log(\Delta)$ by an
operator version of Jensen’s inequality [14]:

$\displaystyle\log(\boldsymbol{\Delta})=\log(P\Delta P^{\dagger})\geq
P\log(\Delta)P^{\dagger},$ (34)

and so we have:

$\displaystyle T\log(\boldsymbol{\Delta})+\boldsymbol{H}\geq
P(T\log\Delta+H)P^{\dagger},$ (35)

which, thanks to Proposition 1, proves the first statement. For the second
statement, $\operatorname{span}(|b_{1}\rangle,\ldots,|b_{r}\rangle)$ is
invariant under the modular flow of $\omega$ precisely when $P$ commutes with
$\Delta$, in which case we have equality in (35). ∎

We conclude this section with a correctness result for the algorithm.

######  Corollary 1.

Suppose $\omega$ is a Gibbs state of a Hamiltonian
$h=\sum_{\alpha=1}^{s}y_{\alpha}h_{\alpha}$ at temperature $T$, and let
$(h^{*},T^{*},\mu^{*})$ be the output of the semidefinite program (6) with
variational Hamiltonian terms $h_{1},\ldots,h_{s}$.

  1. i)

(Feasibility) The triple $(h,T,0)$ is a feasible point of the program (6). In
particular $\mu^{*}\geq 0$.

  2. ii)

(Recoverability) Suppose $b_{1},\ldots,b_{r}$ span the entire algebra
$\mathcal{A}$. Then $(h^{*},T^{*},\mu^{*})=(h,T,0)$.

###### Proof.

$i)$: by Proposition 2, $\omega$ satisfies RTS. The statement then follows
from Proposition 3.

$ii)$: if $b_{1},\ldots,b_{r}$ span $\mathcal{A}$ then any solution of the
constraint (7) with $\mu\geq 0$ satisfies (33), which by Propositions 2 and 3
implies that $\omega$ is the Gibbs state $h^{*}$ at temperature $T^{*}$, and
so $(h,T)=(h^{*},T^{*})$ and $T\log(\boldsymbol{\Delta})+\boldsymbol{H}=0$.
This last equality implies in turn that $\mu=0$. ∎

### Theoretical justification of algorithm

Now we make explicit the connection between the results proved in this section
and the algorithm. The operators $a_{1},\ldots,a_{r}$ constructed in step 1
are chosen so that $|a_{1}\rangle,\ldots,|a_{r}\rangle$ form an orthonormal
basis of $\operatorname{span}(|b_{1}\rangle,\ldots,|b_{r}\rangle)$, and we
have

$\displaystyle\boldsymbol{\Delta}_{ij}:=\omega(a_{j}a_{i}^{*})=\langle
a_{i}|\Delta|a_{j}\rangle$ (36)

and

$\displaystyle\boldsymbol{H}_{ij}:=\omega(a_{i}^{*}[h,a_{j}])=\langle
a_{i}|H|a_{j}\rangle.$ (37)

for any selfadjoint $h\in\mathcal{A}$. In other words we have
$\boldsymbol{\Delta}=P\Delta P^{\dagger}$ and $\boldsymbol{H}=PHP^{\dagger}$.

In step 2, the kernel of $W$ computes the space of quasi-symmetries of
$\omega$ that lie in the span of $h_{1},\ldots,h_{\alpha}$. Indeed, by
Proposition 1 $i)$, a self-adjoint $h\in\mathcal{A}$ is a quasi-symmetry if
and only if $\boldsymbol{H}$ is self-adjoint, and we have have
$W=F^{\dagger}F$, where $F:\operatorname{\mathbb{R}}^{s}\to Mat_{r\times
r}(\operatorname{\mathbb{C}})$ is the map taking
$x\mapsto\sum_{\alpha=1}^{s}x_{\alpha}(\boldsymbol{H}_{\alpha}-\boldsymbol{H}_{\alpha}^{\dagger})$.

Correctness of step 3 is given by Corollary 1.

##  4 Numerical results

We now turn to the results of numerical simulations of the algorithm.

### Corrections to ideal algorithm

In the presence of noise in the expectation values of $\omega$, it is
inappropriate to use the exact kernel of the matrix $W$ constructed in step 2.
Instead, a threshold $\epsilon_{W}>0$ was chosen as a hyperparameter, and we
let $\tilde{h}_{1},\ldots,\tilde{h}_{q}$ be the eigenvectors of $W$ with
eigenvalue smaller than $\epsilon_{W}$. Since the matrices
$\widetilde{\boldsymbol{H}}_{\alpha}$ can no longer be assumed to be self-
adjoint, we symmetrize their definition:

$\displaystyle(\widetilde{\boldsymbol{H}}_{\alpha})_{ij}:=\frac{\omega(a_{i}^{*}[\tilde{h}_{\alpha},a_{j}])+\overline{\omega(a_{j}^{*}[\tilde{h}_{\alpha},a_{i}])}}{2}.$
(38)

The proof of feasibility in Corollary 1 still holds, since for
$\epsilon_{W}=0$ this reduces to the original algorithm, while increasing
$\epsilon_{W}$ (for a fixed $\omega$) increases $\mu^{*}$ (indeed, increasing
$\epsilon_{W}$ only increses the feasible region of the program (6)).

For sufficiently small values of $\sigma_{noise}$, the spectrum of $W$ was
found to have a low-lying part with a spectral gap to the rest of the
eigenvalues. $\epsilon_{W}$ was set according to the formula

$\displaystyle\epsilon_{W}=400\max(\sigma_{noise}^{2}\sqrt{m},\ 10^{-11}),$
(39)

where $m$ denotes the number of terms $a_{i}^{*}[h_{\alpha},a_{j}]$ such that
$[h_{\alpha},a_{j}]\neq 0$. This formula was found empirically to produce an
$\epsilon_{W}$ lying in the spectral gap of $W$, and is not expected to be
universal across different values of $n$ and choices of perturbing operators.
We note that in practice, while choosing $\epsilon_{W}$ to be too low caused
the output to be inaccurate, choosing $\epsilon_{W}$ to lie above the gap did
not significantly affect the accuracy of the result.

### Learning the XXZ Hamiltonian

The MPS purification technique [15] was used to prepare thermal states of the
following anisotropic Heisenberg ferromagnet:

$\displaystyle
h=-\sum_{i=1}^{n-1}(\sigma_{i}^{x}\sigma_{i+1}^{x}+\sigma_{i}^{y}\sigma_{i+1}^{y}+\frac{1}{2}\sigma_{i}^{y}\sigma_{i+1}^{y}).$
(40)

with $n=100$. Both the set of perturbing operators $b_{1},\ldots,b_{r}$ and
the set of variational Hamiltonian terms $h_{1},\ldots,h_{s}$ were chosen to
be the $1192$ geometrically $2$-local Pauli operators. Measurement error was
simulated by adding Gaussian noise with variance $\sigma_{noise}$ to the
expectation value of each Pauli operator. The learning algorithm itself was
implemented in Python, using the MOSEK solver [16] for the semidefinite
optimization333The convex modeling language CVXPy [17] and the open source
solver SCS [18] were used in prototyping but not in the final code..

Hamiltonian recovery error was quantified using the overlap as in [19]: let
$y\in\operatorname{\mathbb{R}}^{s}$ be the vector of recovered Hamiltonian
coefficients and $z\in\operatorname{\mathbb{R}}^{s}$ the vector of true
Hamiltonian coefficients. The Hamiltonian recovery error is then defined as
the relative angle of the two, which for small angles approximately equals the
reciprocal of the signal-to-noise ratio:

$\displaystyle\theta=\arccos\left(\frac{|\langle
y|z\rangle|}{\|y\|\|z\|}\right)\approx\frac{\|y-z\|}{\|z\|}.$ (41)

Note that this metric is not sensitive to the overall scaling of the
Hamiltonian. This degree of freedom of the Hamiltonian is fixed by the
normalization (8) anyway and instead appears in the temperature $T$.
Interestingly, the algorithm reconstructed the “projective” degrees of freedom
of the Hamiltonian terms much more accurately than it did its overall scale
(or equivalently, the temperature).

The Hamiltonian recovery error $\theta$ and the recovered temperature $T$ are
plotted against $\sigma_{noise}$ in Figure 1. A temperature-dependent noise
threshold is found between $\sigma_{noise}\approx 10^{-5}$ and
$\sigma_{noise}\approx 10^{-3}$ above which the matrix $\boldsymbol{\Delta}$
ceases to be positive definite. The algorithm could possibly be emended to
work for higher noise values by projecting onto the orthocomplement of the
nonpositive eigenspace of $\boldsymbol{\Delta}$, but we leave this to future
work.

Figure 1: Numerical results for the 100-qubit anisotropic Heisenberg model
(40) at several temperatures. _Left:_ Recovery error $\theta$ as a function of
noise amplitude $\sigma_{noise}$, averaged over 10 runs. Dotted line is (mean)
+ (standard deviation). _Right:_ Ratio of recovered temperature to actual
temperature, averaged over 10 runs. Shaded region is (mean) $\pm$ (standard
deviation).

As one shrinks the noise amplitude, the recovery error first decreases (for
high temperatures, this decrease is linear to a good approximation). This
persists up until, at some temperature-dependent critical value of the noise
amplitude, the recovery error plateaus. We interpret this two-stage behaviour
as follows. In the limit of zero measurement error, perfect recovery is not
guaranteed because the condition (33) is weaker than the Gibbs condition.
Instead, (33) defines a convex set of candidate Hamiltonians, and the
algorithm picks one of these by maximizing the regularization parameter $\mu$.
The recovery error is then on the order of the diameter of this convex set.
Thus for low enough levels of measurement noise the recovery error is roughly
constant.

We posit that the only way to lower the levels of these plateaux is to enlarge
the set of perturbing operators, which tightens the constraint (33). This is
relevant if one wants to prove asymptotic bounds on the number of copies of
the state and the computational resources needed to specify the Hamiltonian up
to an arbitrarily low error.

Note however that for the particular Hamiltonian under consideration, the
plateaux start at noise amplitudes $\sigma_{noise}$ of around $10^{-9}$ to
$10^{-8}$. Assuming that expectation values are estimated from independent
copies of the state, this would require on the order of $10^{16}$ to $10^{18}$
samples, far beyond what is experimentally feasible anyway. So for practical
applications it may be more important to understand the high-noise regime
rather than the locations of the plateaux.

##  5 Discussion and outlook

Let us conclude by describing two directions for future research. While
Corollary 1 establishes the correctness of the algorithm, it suffers from two
important limitations which must be overcome if one is to prove sample
complexity and computational complexity bounds. First, neither the feasibility
nor the recoverability statements of Corollary 1 take into account measurement
noise in the expectation values of the state, which is unavoidable whenever
these are estimated using finitely many copies of the state.

Second, the recoverability statement only holds when the set of perturbing
operators is grown to a complete set of operators. The utility of this
algorithm depends on approximate recoverability when the set of perturbing
operators is far smaller than a complete set. Section 4 gives numerical
evidence that this is indeed the case, but a proof is still lacking.

The second point above leads one to consider states that satisfy RTS for some
Hamiltonian $h$, but are not Gibbs states of $h$. The algorithm described here
can never tell these apart from true Gibbs states. Aside from their
implications for the analysis of the algorithm, it is an interesting physical
question to rule out or characterize such “false Gibbs states”, which are
thermodynamically stable against some perturbations but not others.

### Acknowledgements

The author would like to thank Anton Kapustin, Hsin-Yuan (Robert) Huang, and
Elliott Gesteau for helpful discussions and comments on the draft.

## References

  * [1] Brian Swingle and Isaac H. Kim  “Reconstructing Quantum States from Local Data”  In _Physical Review Letters_ 113.26 American Physical Society (APS), 2014  DOI: 10.1103/physrevlett.113.260501
  * [2] Hamza Fawzi, Omar Fawzi and Samuel O. Scalet  “Certified algorithms for equilibrium states of local quantum Hamiltonians”  DOI: 10.48550/arXiv.2311.18706
  * [3] Huzihiro Araki and Geoffrey L. Sewell  “KMS conditions and local thermodynamical stability of quantum lattice systems”  In _Communications in Mathematical Physics_ 52.2 Springer ScienceBusiness Media LLC, 1977, pp. 103–109  DOI: 10.1007/bf01625778
  * [4] Anurag Anshu, Srinivasan Arunachalam, Tomotaka Kuwahara and Mehdi Soleimanifar  “Sample-efficient learning of interacting quantum systems”  In _Nature Physics_ 17.8 Springer ScienceBusiness Media LLC, 2021, pp. 931–935  DOI: 10.1038/s41567-021-01232-0
  * [5] Jeongwan Haah, Robin Kothari and Ewin Tang  “Optimal learning of quantum Hamiltonians from high-temperature Gibbs states”  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ IEEE, 2022  DOI: 10.1109/focs54457.2022.00020
  * [6] Ainesh Bakshi, Allen Liu, Ankur Moitra and Ewin Tang  “Learning quantum Hamiltonians at any temperature in polynomial time”  arXiv, 2023  DOI: 10.48550/ARXIV.2310.02243
  * [7] Andi Gu, Lukasz Cincio and Patrick J. Coles  “Practical Hamiltonian learning with unitary dynamics and Gibbs states”  In _Nature Communications_ 15.1 Springer ScienceBusiness Media LLC, 2024  DOI: 10.1038/s41467-023-44008-1
  * [8] Tim J. Evans, Robin Harper and Steven T. Flammia  “Scalable Bayesian Hamiltonian learning”  arXiv, 2019  DOI: 10.48550/ARXIV.1912.07636
  * [9] Yotam Y. Lifshitz et al.  “Practical Quantum State Tomography for Gibbs states”  arXiv, 2021  DOI: 10.48550/ARXIV.2112.10418
  * [10] Eyal Bairey, Itai Arad and Netanel H. Lindner  “Learning a Local Hamiltonian from Local Measurements”  In _Physical Review Letters_ 122.2 American Physical Society (APS), 2019  DOI: 10.1103/physrevlett.122.020504
  * [11] Eberhard E. Muller  “Note on relative entropy and thermodynamical limit”  In _Helvetica Physica Acta_ 58 Birkhäuser Verlag Basel, 1985, pp. 622–632  URL: https://www.e-periodica.ch/cntmng?pid=hpa-001:1985:58::1152
  * [12] G. Lindblad  “On the generators of quantum dynamical semigroups”  In _Communications in Mathematical Physics_ 48.2 Springer ScienceBusiness Media LLC, 1976, pp. 119–130  DOI: 10.1007/bf01608499
  * [13] Robert Alicki and Karl Lendi  “Quantum Dynamical Semigroups and Applications”  Springer Berlin Heidelberg, 2007  DOI: 10.1007/3-540-70861-8
  * [14] Frank Hansen and Gert K. Pedersen  “Jensen’s Operator Inequality”  In _Bulletin of the London Mathematical Society_ 35.04 Wiley, 2003, pp. 553–564  DOI: 10.1112/s0024609303002200
  * [15] Adrian E. Feiguin and Steven R. White  “Finite-temperature density matrix renormalization using an enlarged Hilbert space”  In _Physical Review B_ 72.22 American Physical Society (APS), 2005  DOI: 10.1103/physrevb.72.220401
  * [16] Erling D. Andersen and Knud D. Andersen  “The Mosek Interior Point Optimizer for Linear Programming: An Implementation of the Homogeneous Algorithm”  In _High Performance Optimization_ Springer US, 2000, pp. 197–232  DOI: 10.1007/978-1-4757-3216-0˙8
  * [17] Steven Diamond and Stephen Boyd  “CVXPY: a python-embedded modeling language for convex optimization”  In _J. Mach. Learn. Res._ 17.1 JMLR.org, 2016, pp. 2909–2913
  * [18] Brendan O’Donoghue, Eric Chu, Neal Parikh and Stephen Boyd  “Conic Optimization via Operator Splitting and Homogeneous Self-Dual Embedding”  In _Journal of Optimization Theory and Applications_ 169.3 Springer ScienceBusiness Media LLC, 2016, pp. 1042–1068  DOI: 10.1007/s10957-016-0892-3
  * [19] Xiao-Liang Qi and Daniel Ranard  “Determining a local Hamiltonian from a single eigenstate”  In _Quantum_ 3 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2019, pp. 159  DOI: 10.22331/q-2019-07-08-159
  * [20] Ola Bratteli and Derek W. Robinson  “Operator Algebras and Quantum Statistical Mechanics I”  Springer Berlin Heidelberg, 1987  DOI: 10.1007/978-3-662-02520-8
  * [21] Masamichi Takesaki  “Fundamentals of Banach Algebras and C*-Algebras”  In _Theory of Operator Algebras I_ Springer New York, 1979, pp. 1–57  DOI: 10.1007/978-1-4612-6188-9˙1

##  Appendix A The Gelfand-Naimark-Segal construction

In this appendix we give a brief introduction to Gelfand-Naimark-Segal (GNS)
construction which is central to the theoretical analysis in section 3. The
GNS construction is a fundamental tool in the theory of operator algebras and
is well described in standard references [20, 21]. However, these references
work in the general setting of infinite-dimensional C∗-algebras, where the
discussion is plagued with analytic subtleties. We avoid these subtleties
because we work in finite dimensions, so it is worthwhile to introduce the
concepts we need in a streamlined and self-contained manner here.

Let $\omega$ be a faithful state of $\mathcal{A}=\mathcal{B}(\mathcal{H})$,
and $\rho\in\mathcal{A}$ its density matrix. The bilinear form
$(a,b)\mapsto\omega(a^{*}b)$ endows $\mathcal{A}$ with the structure of a
Hilbert space, which we call the GNS space and denote by the symbol
$\mathcal{H}_{\omega}$ (although $\mathcal{A}$ and $\mathcal{H}_{\omega}$ are
isomorphic as vector spaces, it is important to maintain their distinction as
mathematical objects). For an operator $a\in\mathcal{A}$ we denote the
corresponding vector in the GNS space by $|a\rangle\in\mathcal{H}_{\omega}$.
The GNS space $\mathcal{H}_{\omega}$ carries two distinguished representations
of the algebra $\mathcal{A}$. These are the left representation:

$\displaystyle\pi_{\ell}(a):|b\rangle\mapsto|ab\rangle,$ (42)

and the right representation:

$\displaystyle\pi_{r}(a):|b\rangle\mapsto|ba^{*}\rangle.$ (43)

Notice that these representations commute in the sense that
$\pi_{r}(a)\pi_{\ell}(b)|c\rangle=\pi_{\ell}(b)\pi_{r}(a)|c\rangle$ for all
$a,b,c\in\mathcal{A}$. The modular operator $\Delta$ is the complex-linear
operator on $\mathcal{H}_{\omega}$ corresponding to the sesquilinear form
$(a,b)\mapsto\omega(ba^{*})$. That is, for any $a,b\in\mathcal{A}$ we have

$\displaystyle\langle a|\Delta|b\rangle$ $\displaystyle:=\omega(ba^{*})$ (44)
$\displaystyle=\operatorname{tr}(\rho ba^{*})$ (45)
$\displaystyle=\operatorname{tr}(\rho a^{*}(\rho b\rho^{-1}))$ (46)
$\displaystyle=\omega(a^{*}(\rho b\rho^{-1}))$ (47) $\displaystyle=\langle
a|\pi_{\ell}(\rho)\pi_{r}(\rho^{-1})|b\rangle,$ (48)

and so $\Delta$ can be expressed as

$\displaystyle\Delta=\pi_{\ell}(\rho)\pi_{r}(\rho^{-1}).$ (49)

The modular flow of $\omega$ is the one-parameter group of unitaries on
$\mathcal{H}_{\omega}$ given by $t\mapsto\Delta^{it}$. It is generated (in the
sense of the Stone’s theorem) by the selfadjoint operator $\log(\Delta)$,
which can be expressed as

$\displaystyle\log(\Delta)=\pi_{\ell}(\log(\rho))-\pi_{r}(\log(\rho)).$ (50)


Hubble Tension in Power-Law $f(R)$ Gravity and Generalized Brans-Dicke Theory
${\bf Yousef~{}Bisabr}$***e-mail<EMAIL_ADDRESS>
Department of Physics, Shahid Rajaee Teacher Training University, Lavizan,
Tehran 16788, Iran

###### Abstract

We introduce a theoretical framework to alleviate the Hubble tension. This
framework is based on dynamics of a minimally coupled scalar field which
either belongs to the Brans-Dicke theory with a self-interacting potential or
is the scalar partner of $f(R)$ gravity. These two theories are dynamically
equivalent when the Brans-Dicke parameter is zero. We will use this dynamical
equivalence to interpret the Hubble tension in the same theoretical framework.
For both theories, we write one set of field equations in which the value of a
parameter distinguishes between the two theories. We will show that $H_{0}$
actually evolves with redshift so that its value is consistent with that
measured from the local distance ladder and it drops to the value measured
from CMB at high redshift. We argue that even though both theories exhibit
this behaviour, Brans-Dicke theory with an exponential potential is more
successful than power-law $f(R)$ gravity to relieve the Hubble tension.

Keywords : Cosmology, Hubble Tension, Modified Gravity.

##  1 Introduction

There are a number of unexplained features in the standard cosmology among
which the accelerating expansion of the Universe has been recently drawn much
amount of attention. There is not still a well-understood explanation on the
nature of this phenomenon. It is usually interpreted as evidence either for
existence of some exotic matter components or for modification of the
gravitational theory. In the first route of interpretation one can take a
perfect fluid with a sufficiently negative pressure, dubbed dark energy [1],
to produce the observed acceleration. In the second route, however, one
attributes the accelerating expansion to a modification of general relativity.
This puzzling picture has been enhanced by the so-called Hubble tension [2].
It consists of a discrepancy between the local value of the Hubble constant
$H_{0}$ based on Cepheids in the Large Magellanic Clouds $H_{0}=74.03\pm
1.42~{}Kms^{-1}Mpc^{-1}$ [3] and the Planck data of the CMB radiation
$H_{0}=67.4\pm 0.5~{}Kms^{-1}Mpc^{-1}$ [4]. Improved determination of the
Hubble constant from the Planck observations and the Hubble Space Telescope
measurements seem to demonstrate that this tension may not be caused by
systematics [5]. Hence, many authors prefer to believe that the tension might
be caused by new physics beyond the $\Lambda$CDM.
These observations suggest that $H_{0}$ as inferred by CMB coming from
$z\simeq 1100$ is lower than that inferred by the nearby sources. It therefore
seems that the Hubble constant measurement could be affected by redshift in
its determination. A recent interesting analysis presented in [6] uses the
Pantheon sample [7], the largest compilation of SNe Ia, to show that there is
an evolution in the Hubble constant which scales as
$g(z)=H_{0}(1+z)^{-\delta}$ with $\delta$ and $H_{0}$ being the evolution
parameter and the Hubble parameter at $z=0$, respectively. They divided the
Pantheon sample into some bins where each bin is populated with an equal
number of SNe Ia. Then they have fitted the values obtained in these bins with
$g(z)$ fitting function with $\alpha\sim 10^{-2}$. This suggests that $H_{0}$
can be considered as a monotonic decreasing function from $z=0$ up to $z\simeq
1100$.
The possibility of a slow decaying of $H_{0}$ with redshift opens theoretical
scenarios for its physical and dynamical interpretation. A possible
proposition is the modified gravity theories such as Brans-Dicke (BD) theory
[8], as a prototype of the scalar-tensor theories [9], or $f(R)$ gravity
models [10]. It is well-known that $f(R)$ gravity, at least within a classical
perspective, are dynamically equivalent to the class of BD theories with a
potential function and a null BD parameter [11]. Thus under a suitable
redefinition of the gravitational and matter fields, one can make their field
equations coincide. We will use this dynamical equivalence to consider Hubble
tension in a generalized BD theory (GBD) and $f(R)$ gravity in the Einstein
frame and under the same analysis framework. We will use an exponential
potential in this analysis. The choice of this kind of potential in $f(R)$
gravity is equivalent to considering a power-law $f(R)$ gravity [12]. We then
determine the Hubble flow and show that even though decaying of $H(z)$ with
redshift is possible for both models GBD is more consistent with the $g(z)$
model proposed in [6].

##  2 The Model

Let us begin with the following action†††We use units in which $8\pi
G=M_{p}^{-2}=1$.

$S^{BD}_{JF}=\frac{1}{2}\int
d^{4}x\sqrt{-\bar{g}}\\{\phi\bar{R}-\frac{\omega_{BD}}{\phi}\bar{g}^{\mu\nu}\bar{\nabla}_{\mu}\phi\bar{\nabla}_{\nu}\phi-V(\phi)+2L_{m}\\}$
(1)

where $\bar{R}$ is the Ricci scalar, $\phi$ is a scalar field, $V(\phi)$ is a
potential function and $L_{m}$ is the matter Lagrangian density. The above
action is the Jordan frame representation of the BD theory with the BD
parameter $\omega_{BD}$ and a self-interacting potential $V(\phi)$. The
potential generalizes the BD theory which is shown to be relevant for studying
some cosmological problems [13].
A conformal transformation

$\bar{g}_{\mu\nu}\rightarrow g_{\mu\nu}=\phi~{}\bar{g}_{\mu\nu}$ (2)

brings the above action into the Einstein frame [14]. Then by a redefinition
of the scalar field

$\varphi=\sqrt{\omega_{BD}+\frac{3}{2}}~{}\ln\phi$ (3)

the kinetic term takes a canonical form. In terms of the new variables
($g_{\mu\nu}$, $\varphi$) the action (1) takes then the form

$S^{BD}_{EF}=\int
d^{4}x\sqrt{-g}\\{\frac{1}{2}R-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\varphi\nabla_{\nu}\varphi-U(\varphi)+e^{-\sigma\varphi}L_{m}\\}$
(4)

where

$\sigma=8\sqrt{\frac{\pi}{2\omega_{BD}+3}}$ (5)

and $\nabla_{\mu}$ is the covariant derivative of the rescaled metric
$g_{\mu\nu}$. The Einstein frame potential is
$U(\varphi)=V(\phi(\varphi))~{}e^{-\sigma\varphi}$. Now variation of the
action (4) with respect to the metric $g_{\mu\nu}$ and $\varphi$ gives,
respectively,

$G_{\mu\nu}=(e^{-\sigma\varphi}T^{m}_{\mu\nu}+T^{\varphi}_{\mu\nu})$ (6)
$\Box\varphi-\frac{dU(\varphi)}{d\varphi}=\sigma e^{-\sigma\varphi}L_{m}$ (7)

where

$T^{\varphi}_{\mu\nu}=(\nabla_{\mu}\varphi\nabla_{\nu}\varphi-\frac{1}{2}g_{\mu\nu}\nabla_{\alpha}\varphi\nabla^{\alpha}\varphi)-U(\varphi)g_{\mu\nu}$
(8) $T^{m}_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}$ (9)

Applying the Bianchi identities to (6) gives

$\nabla^{\mu}(e^{-\sigma\varphi}T_{\mu\nu}^{m})=-\nabla^{\mu}T_{\mu\nu}^{\varphi}$
(10)

which means that $T_{\mu\nu}^{m}$ and $T_{\mu\nu}^{\varphi}$ do not separately
conserved. Here we consider a perfect fluid energy-momentum tensor as a matter
system

$T^{m}_{\mu\nu}=(\rho_{m}+p_{m})u_{\mu}u_{\nu}+p_{m}g_{\mu\nu}$ (11)

where $\rho_{m}$ and $p_{m}$ are energy density and pressure, respectively.
The four-velocity of the fluid is denoted by $u_{\mu}$. Details of the energy
exchange between matter and $\varphi$ depends on the explicit form of $L_{m}$.
There are different choices for $L_{m}$ which all of them leads to the same
energy-momentum tensor and field equations in the context of general
relativity [15] [16]. Here we take $L_{m}=p_{m}$ for the lagrangian density.
The $f(R)$ gravity models propose a modification of the Einstein-Hilbert
action so that the Ricci scalar is replaced by some arbitrary function $f(R)$.
There is a relation between the BD theory and $f(R)$ gravity in metric
formalism at the classical level. In fact metric $f(R)$ gravity is dynamically
equivalent to GBD with a null BD parameter $\omega_{BD}=0$ [11]. It means that
they are described by the same action and the resulting field equations
coincide. The action for a metric $f(R)$ gravity in the Jordan frame is
written as

$S^{f(R)}_{JF}=\frac{1}{2}\int d^{4}x\sqrt{-g}\\{f(R)+2L_{m}\\}$ (12)

It can be easily shown that this action has a scalar-tensor representation
[11] [17]

$S^{f(R)}_{JF}=\frac{1}{2}\int d^{4}x\sqrt{-g}\\{\chi R-V(\chi)+2L_{m}\\}$
(13)

where the scalar field is defined as $\chi=f^{\prime}(R)$. This is the BD
theory in Jordan frame with $\omega_{BD}=0$. A conformal transformation
$g_{\mu\nu}\rightarrow\Omega^{2}g_{\mu\nu}$ with $\Omega=\sqrt{\chi}=f^{\prime
1/2}(R)$ and the scalar field redefinition
$\xi=\sqrt{\frac{3}{2}}\ln\chi$‡‡‡This is the field redefinition (3) with
$\omega_{BD}=0$. bring the action (13) into the Einstein frame representation

$S^{f(R)}_{EF}=\int
d^{4}x\sqrt{-g}\\{\frac{1}{2}R-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\xi\nabla_{\nu}\xi-U(\xi)+e^{-\sqrt{\frac{2}{3}}\xi}L_{m}\\}$
(14)

where

$U(\xi)=\frac{Rf^{\prime}(R)-f(R)}{2f^{\prime 2}(R)}$ (15)

with prime being derivative with respect to $R$. Variations of (14) with
respect to $g^{\mu\nu}$ and $\xi$ gives the same field equations (6), (7) and
(8) when $\varphi$ is replaced by $\xi$ and $\sigma$ takes the value
$\sigma=8\sqrt{\frac{\pi}{3}}\simeq 8.2$.
As the last point in this section we comment that although Solar System
experiments set the constraint $\omega_{BD}>40000$ [18], the null BD parameter
does not generally mean that $f(R)$ gravity models do not satisfy
observational constraints. In fact, some authors have already argued that
since the post-Newtonian parameter satisfies $\gamma=\frac{1}{2}$ instead of
being equal to unity as required by observations, all f(R) theories should be
ruled out [19]. Later, it was noted that for scalar fields which have
sufficiently large masses it is possible for $\gamma$ to be close to unity
even for null BD parameter. In this case the scalar fields become short-ranged
and have no effect at Solar System scales [12].

##  3 Cosmological implementation

We apply the field equations (6) and (7) to a spatially flat Friedmann-
Robertson-Walker spacetime $ds^{2}=-dt^{2}+a^{2}(t)(dx^{2}+dy^{2}+dz^{2})$
with $a(t)$ being the scale factor. This gives

$3H^{2}=e^{-\sigma\varphi}\rho_{m}+\rho_{\varphi}$ (16)
$2\dot{H}+3H^{2}=-(e^{-\sigma\varphi}p_{m}+p_{\varphi})$ (17)
$\ddot{\varphi}+3H\dot{\varphi}+\frac{dU(\varphi)}{d\varphi}=-\sigma
e^{-\sigma\varphi}p_{m}$ (18)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter. The conservation
equations become

$\dot{\rho}_{\varphi}+3H(\omega_{\varphi}+1)\rho_{\varphi}=-\sigma
e^{-\sigma\varphi}\dot{\varphi}\rho_{m}$ (19)
$\dot{\rho}_{m}+3H(\omega_{m}+1)\rho_{m}=(1-2\omega_{m})\sigma\dot{\varphi}\rho_{m}$
(20)

The latter can be solved which gives the following solution

$\rho_{m}=\rho_{0m}a^{-3(\omega_{m}+1)}e^{(1-2\omega_{m})\sigma\varphi}$ (21)

where $\rho_{0m}$ is an integration constant. This solution can also be
written as [20]

$\rho_{m}=\rho_{0m}a^{-3(\omega_{m}+1)+\varepsilon}$ (22)

where we have defined

$\varepsilon\equiv\lambda(1-2\omega_{m})\frac{\sigma\varphi}{\ln a}$ (23)

with $\lambda$ being a parameter. This solution indicates that the evolution
of energy density is modified due to interaction of $\varphi$ with matter. For
$\varepsilon>0$, matter is created. In this case, energy is injecting from
$\varphi$ into the matter so that the latter dilutes more slowly compared to
the standard evolution $\rho_{m}\propto a^{-3(\omega_{m}+1)}$. For
$\varepsilon<0$, on the other hand, matter is annihilated and energy transfers
outside of the matter system. In this case, the rate of dilution of $\rho_{m}$
is faster than the standard one.
We now search for solutions for the field equations (16), (17) and (18) which
are able to reproduce the current Universe expansion and can specifically
alleviate the tensions in the Hubble constant measurements. In particular, we
are searching for those solutions for which $H(z)$ receives variations with
the redshift as a consequence of evolution of the non-minimally coupled scalar
field $\varphi(z)$. There are two independent equations among (16)-(18) which
ultimately give $H(z)$ and $\varphi(z)$. Moreover, the conservation equations
(19) and (20) characterize the energy exchange between $\varphi$ and matter
which modify evolution of $\rho_{m}$ and $\rho_{\varphi}$. In this interacting
system, $\varepsilon$ in (22) is generally an evolving function characterizing
the rate of the energy transfer. However, we consider the case that
$\varepsilon$ can be regarded as a constant parameter. Even though this is not
generally true during the whole expansion history of the Universe, it may hold
during particular eras or during short periods in the evolution of the
Universe. In this case, (23) reduces to

$\varphi=\frac{\varepsilon}{\lambda\sigma(1-2\omega_{m})}\ln a$ (24)

which implies that the rate of change of the scalar field is given by the
Hubble parameter, namely $\dot{\varphi}\propto H$. Now we study the Friedman
equation (16) for GBD theory and $f(R)$ gravity:
_1) GBD theory._ As the simplest and the best-studied generalization of
general relativity, it is natural to think about the BD scalar field as a
possible candidate for producing modifications of cosmic expansion. We first
consider the potential $V(\phi)=V_{0}e^{\alpha\phi}$ in the action (1) with
$V_{0}$ and $\alpha$ being constants. This is equivalent to
$U(\varphi)=V_{0}e^{(\alpha-\sigma)\varphi}$ in the action (4). In this case,
the Friedman equation (16) becomes

$\frac{H(z)}{H_{0}}=(3-\frac{\varepsilon^{2}}{2\lambda^{2}\sigma^{2}})^{-1}\\{3\Omega_{0m}(1+z)^{-3+\varepsilon(1-\frac{1}{\lambda})}+(1+z)^{\frac{\varepsilon}{\lambda\sigma}(\alpha-\sigma)}\\}^{\frac{1}{2}}$
(25)

Figure 1: The plots of $H(z)/H_{0}$ in BD theory with an exponential potential
(panel a) and $f(R)$ gravity (panel b). In these plots $\omega_{m}=0.3$;
$\alpha=0.04$; $\varepsilon=3.3$, 3.4, 3.5 and $\lambda=45,47,49$ from top to
bottom. The dashed line denotes the fitting function $g(z)/H_{0}$ used in [6].

where we have set $V_{0}=H_{0}^{2}$ and $\Omega_{m}\equiv\rho_{0m}/\rho_{c}$,
$\rho_{c}\equiv 3H^{2}_{0}$ and $a=(1+z)^{-1}$ are used§§§We have used
$\frac{a_{0}}{a}=1+z$ with $a_{0}=1$.. Here $H_{0}$ is taken as the late-time
Hubble constant or that one measured by the SNe Ia. In this equation,
$\sigma<0.05$ which is inferred from the observational bound
$\omega_{BD}>40000$ and (5). The function $H(z)/H_{0}$ is plotted in Fig.
1(a). The figure shows that the non-minimally coupled scalar field is actually
responsible for scaling of $H(z)$. Interestingly, this scaling behaviour is
similar to that of the fitting function $g(z)/H_{0}$ (the dashed line) used in
the analysis of the Pantheon sample in [6]. Both functions increase slowly
with decreasing of $z$ and go to $H_{0}$ at late times. This means that the
GBD theory is a good theoretical interpretation for the model
$g(z)=H_{0}(1+z)^{-\delta}$ with $\delta=0.02$.
_2) f(R) Gravity._ The Hubble parameter $H(z)$ in the $f(R)$ gravity models is
also given by (25) when $\sigma$ takes the value $\sigma\simeq 8.2$. In this
case, exponential potentials correspond to power-law $f(R)$ gravity [12]. The
function $H(z)/H_{0}$ is plotted in Fig. 2(b). The figure shows that although
$H(z)$ in power-law $f(R)$ gravity has the same scaling behaviour as GBD
theory, it does not goes to $H_{0}$ at late times. Thus GBD theory provides a
better theoretical framework for understanding discrepancies in the Hubble
constant measurements.

##  4 Conclusions

In this work, we have considered the possibility that the Hubble constant
measurements are affected by redshift so that measurements based on CMB at
$z\simeq 1100$ are smaller than those based on local objects such as Cepheids
in the Large Magellanic Clouds and SNe Ia. This possibility is supported by
analyses such as [6] [21] which fitted extracted $H_{0}$ values with a
function mimicking the redshift evolution $g(z)\propto(1+z)^{-0.02}$. A
theoretical route of interpretation of this observation is the modified
theories of gravity.
We have used the dynamical equivalence of $f(R)$ gravity and GBD with null BD
parameter to investigate Hubble tension under the same analysis framework. In
this analysis a pure exponential potential is used which corresponds to a
power-law $f(R)$ gravity. The two class of theories are described by the same
set of field equations in Einstein frame so that they are distinguished by the
parameter $\sigma$. For GBD theories, $\omega_{BD}>40000$ corresponds to the
upper bound $\sigma<0.05$ and for $f(R)$ gravity the null BD parameter is
equivalent to $\sigma\simeq 8.2$. We have shown that there is a monotonic
decreasing trend in the Hubble constant $H_{0}$ in both theories which mostly
happens for $z<2$. However, there are two weakness of the $f(R)$ gravity:
First, in order that power-law $f(R)$ gravity be consistent with local gravity
experiments the exponent of the curvature scalar hardly deviates from unity.
Second, contrary to GBD the asymptotic behaviour of $H_{0}$ at late times is
not consistent with that measured from the local distance ladder.

## References

  * [1] A. Melchiorri et al., Phy. Rev. D 68, 043509, (2003)
J. S. Alcaniz, Phy. Rev. D 69, 083521, (2004)
T. R. Choudhury and T. Padmanabhan, Astron. Astrophys. 429, 807(2004)

  * [2] E. Di Valentino ei al., Class. Quant. Grav. 38, 153001 (2021)
  * [3] A. G. Riess, Nature Reviews Physics (NatRP) 2, 10 (2020)
  * [4] N. Aghanim et al., A&A 641, A6 (2020)
  * [5] A. G. Riess et al., Astroph. J. 876, 85 (2019)
A. G. Riess et al., Astroph. J. Lett. 934, L7 (2022)

  * [6] M. G. Dainotti et al., Astrophys. J. 912, 150 (2021)
L. Kazantzidis and L. Perivolaropoulos, Phys. Rev. D 102, 023520 (2020)

  * [7] D. M. Scolnic, et al., Astrophys. J. 859, 101 (2018)
  * [8] C. Brans and R. H. Dicke, Phys. Rev. 124, 925 (1961)
  * [9] Y. Fujii and K. Maeda, _The Scalar-Tensor Theory of Gravitation_ , Cambridge University Press (2003)
  * [10] V. Faraoni and S. Capozziello, _Beyond Einstein Gravity: A Survey of Gravitational Theories_ , Springer (2011)
  * [11] T. P. Sotiriou, Rev. Mod. Phys. 82, 451 (2010)
  * [12] Y. Bisabr, Grav. Cosmol. 16, 239 (2010)
  * [13] C. Santos and R. Gregory, Annals Phys. 258, 111 (1997)
M. K. Mak and T. Harko, EPL 60, 155 (2002)
Y. Bisabr, Astrophys Space Sci 339, 87 (2012)

  * [14] V. Faraoni, E. Gunzig and P. Nardone, Fund. Cosmic Phys. 20, 121 (1999)
  * [15] S. W. Hawking and G. F. R. Ellis, ”The Large Scale Structure of Space-Time” (Cambridge 1973, Cambridge University Press)
  * [16] B. Schutz, Phys. Rev. D 2 2762 (1970)
J. D. Brown, Class. Quant. Grav. 10, 1579 (1993)

  * [17] S. Capozziello and V. Faraoni, _Beyond Einstein Gravity: A Survey of Gravitational Theories for Cosmology and Astrophysics_ , Springer (2011)
  * [18] C. M. Will, Liv. Rev. Rel. 9, 3 (2005)
  * [19] T. Chiba, Phys. Let. B 575, 1 (2003)
  * [20] Y. Bisabr, Phys. Rev. D 86 127503 (2012)
Y. Bisabr, Gen. Rel. Grav. 44, 427 (2012)

  * [21] M. G. Dainotti, et al., Galaxies 10, 24 (2022)


# Voting in Transfer Learning System for Ground-Based Cloud Classification

Mario Manzo  IT Services, University of Naples “L’Orientale”
80121, Naples, Italy<EMAIL_ADDRESS>Simone Pellino  Department of Applied
Science,
I.S. Mattei Aversa 81031, M.I.U.R. Rome, Italy<EMAIL_ADDRESS>
###### Abstract

Clouds classification is a great challenge in meteorological research. The
different types of clouds, currently known and present in our skies, can
produce radioactive effects that impact on the variation of atmospheric
conditions, with the consequent strong dominance over the earth’s climate and
weather. Therefore, identifying their main visual features becomes a crucial
aspect. In this paper, the goal is to adopt a pretrained deep neural networks
based architecture for clouds image description, and subsequently,
classification. The approach is pyramidal. Proceeding from the bottom up, it
partially extracts previous knowledge of deep neural networks related to
original task and transfers it to the new task. The updated knowledge is
integrated in a voting context to provide a classification prediction. The
framework trains the neural models on unbalanced sets, a condition that makes
the task even more complex, and combines the provided predictions through
statistical measures. Experimental phase on different cloud image datasets is
performed and results achieved show the effectiveness of the proposed approach
with respect to state of the art competitors.

###### keywords:

Cloud classification , Deep learning , Transfer learning , Voting-based
classification , Climate change

††journal: Expert Systems with Applications

##  1 Introduction

Clouds are a constant presence in our skies, combined with the role assumed by
ecosystems, are also important in determining atmospheric conditions, hours of
sunshine and temperature. Nowadays, dynamics atmospheric conditions,
attributed to climate change, have led to an increase in attention to behavior
of clouds Duda et al. (2013). Climate models allow to predict the climate
changes, but their precision degree is currently insufficient and attributable
to the alteration in the conditions determined by different phenomena.
Consequently, clouds behavior prediction is important for estimating climate
change. Furthermore, cloud changes are also a reason for influences on the
earth radiation budget and energy balance Rossow and Schiffer (1991); Chen et
al. (2000); Stephens (2005). A great deal of effort has been made by the
scientific community to acquire many datasets, but not fully managed due to a
lack of devices with adequate processing resources. In recent years, however,
the state of the art has seen the birth of numerous studies and analyzes that
can automate clouds attributes and classes detection scientifically relevant.
Due to the amount of ground cloud images, the recognition phase has been
extensively studied lately in the literature. Most of the standard algorithms
adopt hand-made features, such as brightness, texture, shape, and color, to
represent image content but without obtaining a model generalization due to
the complex distribution of data. Indeed, the visual information contained in
the image are unable to accurately describe the clouds due to the large
variations in appearance. Finally, the non-visual features, also known as
multimodal information, obtainable from the clouds process formation such as
temperature, humidity, pressure and wind speed can be of help. According to
recent studies, deep learning has proven effective for image management,
analysis, representation and classification also in cloud recognition field.
In particular, the success of deep neural networks, applied to the image
classification task, concern several interesting aspects mainly connected to
the software development and the large amount of data available. Specifically,
for cloud images analysis, deep neural networks are adopted both in the
segmentation and detection phases. However, image content and poor data
balance among classes have a decisive impact on performance, causing
uncertainty in the model generalization. With purpose to address the above
problems, we present a framework based on deep transfer and voting learning
for cloud image classification. It is built based on three steps. A first,
which performs image preprocessing operations such as resizing, essential for
neural networks training. A second, which modifies and retrains multiple deep
neural networks, exploiting previous knowledge. A third, which looks at the
different predictions provided by the deep neural networks and combines them
in order to provide the best decision in the classification phase. The main
contributions about proposed framework can be summarized in some keypoints.
First, a framework based on deep and voting learning designed to address the
imbalance between classes in the cloud recognition task. Second, a framework
built on multiple classification models based on deep transfer learning.
Third, the demonstration that several models, suitably combined, can
strengthen the decision in classification with respect to a single one.
Finally, experimental demonstrations compared to established existing methods
on the datasets recognized by field experts. The paper is structured as
follows. Section 2 provides an overview of state of art about clouds
classification approaches. Section 3 describes in detail proposed framework.
Section 4 provides a wide experimental phase, while section 5 concludes the
paper.

##  2 Related Work

In this section, we briefly describe the most important studies about clouds
images classification in literature. In this field, numerous works address the
task according to different aspects such as image characterization,
segmentation algorithm application to get new descriptors, complex mechanisms
of learning and classification and much more.

In Liu et al. (2018) authors present a layer named joint fusion (JFCNN) to
jointly learn two kinds of cloud features under one framework. After training
the proposed JFCNN, they extract the visual and multimodal features from two
subnetworks, which are based on the well known Resnet50 He et al. (2016) and
integrate them using a weighted strategy. The architecture consists of five
parts: two subnetworks, one joint fusion layer, one FC layer and the loss
function. The subnetworks are used for learning cloud visual features. Authors
work with two kind of extracted features, combined in multimodal way, which
contain some complementary information and different characteristics of the
ground-based cloud.

An approach named multi-evidence and multi-modal fusion network (MMFN) is
proposed in Liu et al. (2020b). The idea is to learn extended cloud
information by fusing heterogeneous features (global and local) in a unified
framework. MMFN takes advantage of multiple pieces of evidence using a main
network and an attentive network. In the attentive network, local visual
features are extracted from attentive maps which are obtained by refining
salient patterns from convolutional activation maps. Meanwhile, the main
network learns multi-modal features for ground-based cloud. In order to
combine the multi-modal and multi-evidence visual features, authors design two
fusion layers in MMFN to incorporate multi-modal features with global and
local visual features, respectively.

In Shi et al. (2017) authors propose to use deep convolutional activations-
based features (DCAFs). Cloud images are directly fed into a CNN model. Then,
the features from different convolutional and FC layers are extracted through
different pooling strategies. Finally, a multilabel linear support vector
machine (SVM) model is used for the classification step.

A convolutional neural network model, called CloudNet, for accurate ground-
based meteorological cloud classification is proposed in Zhang et al. (2018).
The model consists in five convolutional layers and two FC layers. In
addition, to optimize the network training, the image input is processed
through a robust strategy that subtracts the mean red-green-blue value of each
pixel over the training set to improve training speed and accuracy.
Furthermore, the authors have created a clouds dataset, called Cirrus Cumulus
Stratus Nimbus (CCSN), which consists of 11 categories under meteorological
standards.

In Liu and Li (2018) authors propose an approach named deep multimodal fusion
(DMF). In order to learn the visual features, CNN models have been applied to
capture texture information. The extracted features, from deeper layers, have
several eligible properties such as invariance and discrimination.
Subsequently authors employ a weighted strategy to integrate visual and
multimodal features. Finally, SVM algorithm to train the classification model
is adopted.

In Li et al. (2020) a deep tensor fusion network is presented in order to hold
spatial information of ground-based cloud images. It fuses cloud visual and
multimodal features at the tensor level.

In Liu et al. (2019) author propose a approach, called Hierarchical Multimodal
Fusion (HMF), which fuses deep multimodal and deep visual features in
different levels. The architecture is composed of two subnetworks, visual
subnetwork and multimodal subnetwork. The visual subnetwork is defined in
order to extract deep visual features from ground-based cloud images employing
Resnet50 He et al. (2016). The multimodal subnetwork is used to learn features
from a vector composed of six FC layers. Classification step through SVM is
managed.

In Sun et al. (2009) author propose a classification method of sky-condition
based on whole sky infrared cloud images, where the Local Binary Patterns
operator (LBP) and the contrast of local cloud image texture (VAR signal) are
combined to classify sky conditions. The correspondence relationship among
traditional cloud classes and instrument-measured cloud classes is suggested.
The approach analyzes the LBP spectra and VAR characteristics for five classes
of clouds.

An automatic cloud classification algorithm is developed in Heinle et al.
(2010), the approach uses a image-mask created by visually identifying image
regions containing discriminative information. Furthermore the approach
extracts a set of mainly statistical features describing the color as well as
the texture of an image. Classification step adopts the k-nearest neighbour
algorithm.

A modified texton-based classification approach that integrates both color and
texture information to improve classification results is proposed in Dev et
al. (2015). Color channel is adopted to generate image descriptors and filter
responses of images across all the categories aggregating them together.
K-means clustering is applied on the concatenated filter responses, producing
the different cluster centers. These clusters centers are the modified-textons
and constitute the texton dictionary. The discriminative histogram model for
each image category is generated by comparing the filter responses of the
pixels with the generated textons in the dictionary.

An ensemble learning method and resource allocation scheme for cloud
observation and classification is proposed in Zhang et al. (2020). Ensemble
methods, like Bagging, AdaBoost and Snapshot are used as a base classifier to
take the cross-semantic and structure features of cloud images.

##  3 Materials and Methods

In this section we describe the proposed framework which includes two
methodologies: deep neural networks Liu et al. (2017) and voting learning
Peteiro-Barral and Guijarro-Berdiñas (2013). The goal is to combine several
deep neural networks with purpose to classify clouds images. Specifically, a
set of competitive models are aligned and provide a range of confidential
decisions useful for making choices during classification. The framework is
composed of three blocks. A first, which performs preprocessing in terms of
image resize. A second, which learns different deep neural networks,
previously redesigned for the specific task. A third, which combines different
potential indications, through voting rules, provided by deep neural networks
for classification purpose. Finally, the framework runs a predetermined number
of iterations in a supervised learning context.

###  3.1 Image resize

One of the drawbacks of neural networks concerns the fixed dimension about the
input layer with reference to the images to be processed (details about
adopted neural networks can be found in table 1 at column 5). Size
normalization, according to the input layer dimension, is essential because it
is not possible to process different or large sized images for the network
training and classification stages. This step does not alter the content of
the image information in any way.

###  3.2 Network design and transfer learning

The transfer learning has been selected as training strategy. The basic idea
is to transfer the knowledge extracted from a source domain to a destination
one, in our case clouds classification. Generally, a pretrained network is
chosen as starting point in order to learn a new task. It turns out to be the
most convenient and forthcoming solution to adopt the representational power
of pretrained deep neural networks. Clearly, it is easy and fast to tune a
network with transfer learning than training a new network from scratch with
randomly initialized weights. For clouds recognition, deep learning
architectures are selected based on their task compliance. The goal is to
train networks on images by redesign their structures in the final layer
according to different outgoing classes. Table 1 supports the description
below about adopted neural models.

Alexnet Krizhevsky et al. (2012) consists of 5 convolutional layers and 3
fully connected layers. It includes the non-saturating ReLU activation
function, better then tanh and sigmoid during training phase.

Googlenet Szegedy et al. (2015) is composed of 22 deep layers. The network is
inspired by LeNet LeCun et al. (1989) but implemented a novel element which is
dubbed an inception module. This module is based on several very small
convolutions in order to drastically reduce the number of parameters. The
architecture reduced the number of parameters from 60 million (AlexNet) to 4
million. Furthermore, it includes batch normalization, image distortions and
Root Mean Square Propagation algorithm.

Densenet201 Huang et al. (2017) is a convolutional neural network with 201
deep layers. Unlike standard convolutional networks composed of $L$ layers
with $L$ one-to-one connections between the current layers and the nexts, it
contains $\frac{L(L+1)}{2}$ direct connections. Specifically, each layer
adopts the feature-maps of all preceding layers and its own feature-maps into
all subsequent layers as inputs.

Resnet18 and Resnet50 He et al. (2016) are inspired by pyramidal cells
contained in the cerebral cortex. They use particular skip connections or
shortcuts to jump over some layers. They are composed of 18 and 50 deep
layers, which with the help of a technique known as skip connection has paved
the way for residual networks.

Nasnetlarge Zoph et al. (2018) is designed on a search space, called NASNet
search space, which enables transferability. The model works by looking for
the best convolutional layer, or cell, and subsequently replicating this layer
in a stack, each with its own parameters to design a convolutional
architecture. Also, a regularization technique, called ScheduledDropPath, that
significantly improves generalization in the model is introduced.

Table 1: Description of adopted pretrained network. Network | Depth | Size (MB) | Parameters (Millions) | Input Size
---|---|---|---|---
Densenet201 | 201 | 77 | 20 |  224 $\times$ 224
Alexnet | 8 | 227 | 61 |  227 $\times$ 227
Googlenet | 8 | 27 | 7 |  224 $\times$ 224
Resnet18 | 18 | 44 | 11.7 |  224 $\times$ 224
Resnet50 | 50 | 96 | 25.6 |  224 $\times$ 224
Nasnetlarge | * | 332 | 88.9 |  331 $\times$ 331

Deep neural networks have been adapted to the clouds classification problem.
Originally, Imagenet dataset Deng et al. (2009), which includes one million
images divided into 1000 classes, is adopted to perform the main training
phase. Generally, a network elaborates an image and provides a prediction
about a class it might belong to with an attached probability. Indeed, a
network is structured to work on different layers. The first concerns the
input image and requires 3 color channels. Next, convolutional layers, which
work with the purpose to extract image features, are placed. The last
learnable and the final classification layers are adopted to classify the
input image. To make the pretrained network compliant to the classification of
new images, the last two layers are replaced with new layers. Often, the last
layer, with its learnable weights, is completely connected. It is removed and
replaced by a new one completely connected with the outputs related to classes
of new data (clouds types). Furthermore, the learning of the new layer,
connected with the transferred layers, can be speeded up by increasing the
learning rate factors. Optionally, the weights of the previous levels can be
left unchanged by resetting the learning rate to zero. This modification
avoids weights update during training and a consequent flattening of the
execution time as it is not necessary to calculate the gradients of the
relative layers. This improvement has a strong impact in the case of small
datasets to avoid overfitting.

###  3.3 Voting based learning

A voting based learning approach is adopted to manage the classification
phase. In particular, among all possible strategies, we selected stacking. It
works by training a single classifier and, subsequently, combines it with
further classifiers. Unlike a standard approach, where weak or strong learners
are adopted, we basically combined several equally powerful models that
predict an outcome with a certain probability. Finally, we joined all the
predictions for a classification result. The general model can be summarized
by the following matrix

$CN=\begin{bmatrix}\beta_{1}i_{1}&\dots&\beta_{1}i_{k}\\\ \vdots&\ddots&\\\
\beta_{n}i_{1}&&\beta_{n}i_{k}\end{bmatrix}$ (1)

each $i_{k}$ represent an image to be classified, taken from the set
$Imgs=\\{i_{1},i_{2},\ldots,i_{k}\\}$ with cardinality $k$, belonging to one
of $x$ classes. Furthermore, each $\beta_{n}$ represent a deep neural network,
taken from the set $C=\\{\beta_{1},\beta_{2},\ldots,\beta_{n}\\}$ with
cardinality $n$, which provides a decision $d\in I\\{1,\dots,x\\}$, with
reference to $i_{k}\in Imgs$ and $x$ membership classes. The set of decisions
can be rearranged through the following matrix $D$

$D=\begin{bmatrix}d_{\beta_{1}i_{1}}&\dots&d_{\beta_{1}i_{k}}\\\
\vdots&\ddots&\\\ d_{\beta_{n}i_{1}}&&d_{\beta_{n}i_{k}}\end{bmatrix}$ (2)

it describes the result of deep neural networks combination and images of the
matrix $CN$ in terms of position, such as $\beta_{n}i_{k}\to
d_{\beta_{n}i_{k}}$. In addition, a score value $s\in S\\{0,\dots,1\\}$ is
associated to each decision $d$ and provides the posterior probability
$P(i|x)$ that an image $i$ could belong to class $x$. Finally, all the score
values relating to the results of the possible combinations of matrix $CN$ are
collected in the matrix $S$

$S=\begin{bmatrix}P(i_{1}|x)_{d_{\beta_{1}i_{1}}}&\dots&P(i_{k}|x)_{d_{\beta_{1}i_{k}}}\\\
\vdots&\ddots&\\\
P(i_{1}|x)_{d_{\beta_{n}i_{1}}}&&P(i_{k}|x)_{d_{\beta_{n}i_{k}}}\end{bmatrix}$
(3)

each element of posterior probability in the matrix $S$ refers to element of
the matrix $CN$, such as $\beta_{n}i_{k}\to d_{\beta_{n}i_{k}}\to
P(i_{k}|x)_{d_{\beta_{n}i_{k}}}$. Moving on, each column of the matrix $D$ is
analyzed with statistical mode and stored in the vector $DM$

$DM=\\{dm_{d_{\beta_{1,\dots,n}i_{1}}},\dots,dm_{d_{\beta_{1,\dots,n}i_{k}}}\\},$
(4)

the generic value $dm$ contains the modal value of the class to which image
$i$ could belong with the average probability score $ds$. In essence, this is
the class to which an image could belong based on the votes given by different
deep neural networks. In this regard, the concept of statistical mode is
introduced. It can be defined as the value which is repeatedly occurred in a
given set

$mode=l+\left(\frac{f_{1}-f_{0}}{2f_{1}-f_{0}-f_{2}}\right)\times h$ (5)

where $l$ is the lower limit of the modal class, $h$ is the size of the class
interval, $f_{1}$ is the frequency of the modal class, $f_{0}$ is the
frequency of the class which precedes the modal class and $f_{2}$ is the
frequency of the class which successes the modal class. The columns of matrix
$D$ are analyzed in order to obtain the values of the most frequent decisions.
This step is performed in order to verify the highest voted classes from
different deep neural networks, contained in the $CN$ set. Moreover, the aim
of mode application is twofold. First, to extract the most frequent value.
Second, to extract its occurrences in terms of indices. For each most frequent
occurrence, modal value, the corresponding score from the matrix $S$ is
extracted. To this end, $DS$ vector is built

$DS=\\{ds_{P(i_{1}|x)_{d_{\beta_{1,\dots,n}i_{1}}}},\ldots,ds_{P(i_{k}|x)_{d_{\beta_{1,\dots,n}i_{k}}}}\\},$
(6)

where each element $ds$ contains the average decision scores with higher
frequency, extracted through the mode, with reference to corresponding column
of the matrix $D$.

##  4 Experimental results

This section describes the experimental phase. In order to train the neural
models, with purpose to perform classification task in a supervised context,
labeled data are need. Consequently, the issue to be addressed concerns the
quantity of data sufficient to produce experimental results. The content of a
large dataset, useful to training and testing, strongly affects the
classification performance. Therefore, the discriminating factor about the
effectiveness of neural models is the amount of data. Contextually, with
purpose to produce compliant performance, the settings reported in recent
cloud classification methods are adopted.

###  4.1 Datasets

The proposed framework on a state-of-art datasets, containing ground-based
clouds images, is tested. Datasets adopted are:

  1. 1.

Multimodal-Ground-based-Cloud-Database (MGCD) Liu et al. (2020a, b). It is
collected in China and consists in cloud images captured by a sky camera with
a fisheye lens under a variety of conditions and multimodal cloud information.
It includes a total amount of 1720 cloud data. Images are divided into seven
classes: cumulus, cirrus, altocumulus, clear sky, stratus, stratocumulus,
cumulonimbus. The number of item of each class varies from 140 to 350, and the
detailed numbers are listed in Table 2.

  2. 2.

Singapore Whole sky IMaging CATegories Database (SWIMCAT) dataset Dev et al.
(2015). It is composed of 784 sky/cloud patch images with 125 x 125 pixels
captured using wide angle high-resolution sky imaging system, a calibrated
ground-based WSI designed by Dev et al. (2014). The dataset is splitted into
five distinct categories: clear sky, patterned clouds, thick dark clouds,
thick white clouds, and veil clouds. Details are present in table 3.

  3. 3.

Cirrus Cumulus Stratus Nimbus (CCSN) dataset Zhang et al. (2018). It contains
only 2,543 unique cloud images with 256 x 256 pixels in the JPEG format and
contains 10 different forms in cloud observation. It is characterized by a
large set of images, making it the largest of the available public cloud
dataset. Details are shown in table 4.

Table 2: Details of MGCD dataset. Label | Cloud Type | Number of samples
---|---|---
1 | Cumulus | 160
2 | Cirrus | 300
3 | Altocumulus | 340
4 | Clear sky | 350
5 | Stratocumulus | 250
6 | Stratus | 140
7 | Cumulonimbus | 180
Table 3: Details of SWIMCAT dataset. Label | Cloud Type | Number of samples
---|---|---
A | Clear Sky | 224
B | Patterned clouds | 89
C | Thick dark clouds | 251
D | Thick white clouds | 135
E | Veil clouds | 85
Table 4: Details of CCSN Dataset. Label | Cloud Type | Number of samples
---|---|---
Ci | Cirrus | 139
Cs | Cirrostratus | 287
Cc | Cirrocumulus | 268
Ac | Altocumulus | 221
As | Altostratus | 188
Cu | Cumulus | 182
Cb | Cumulonimbus | 242
Ns | Nimbostratus | 274
Sc | Stratocumulus | 340
St | Strtus | 202
Ct | Contrails | 200

###  4.2 Discussion

Table 5 provides indications about the adopted neural models with respect to
datasets. As can be seen, different combinations are provided due to the
variable composition (total images, images per class, etc) of each datasets.
In fact, two additional neural models, Resnet18 and Nasnetlarge, have been
added for CCSN processing as it is the most complex to manage. Moreover, the
framework consists in different modules written in Matlab language. Neural
models were trained based on different parameters. Stochastic Gradient Descent
(SGDM) with Momentum is adopted as solver of training process. Its main
peculiarity concerns the oscillation along the steepest descent path towards
the optimum. Adding a momentum term to the parameter update is one way to
reduce this oscillation. Carrying on, MiniBatchSize value, the subset size of
the training set adopted to evaluate the gradient of the loss function and
update the weights, has been set to 10, optimal for the obtained results.
About MaxEpochs, maximum number of epochs to use for training, the right
compromise was reached with the value 6 optimizing execution time and
performance. An iteration is a step performed by SGDM to minimize the loss
function using MiniBatchSize. An epoch concerns the complete cycle of the
training process on training set. InitialLearnRate has been set to 3e-4. If it
is too low results in a high training time. Otherwise, if too high, the result
may be suboptimal or training may diverge. The right compromise has also been
found for the latter. To avoid discarding the same data every epoch, Shuffle
parameter been set to every-epoch. Finally, ValidationFrequency, number of
iterations between evaluations of validation metrics, has been set as the
ratio between training set size and MiniBatchSize.

Table 5: Deep neural networks adopted with respect to datasets. Networks Datasets | MGCD | SWIMCAT | CCSN
---|---|---|---
Densenet201 | ✔ | ✔ | ✔
Alexnet | ✔ | ✔ | ✔
Googlenet | ✔ | ✔ | ✔
Resnet18 |  |  | ✔
Resnet50 | ✔ | ✔ | ✔
Nasnetlarge |  |  | ✔

The classification accuracy on MGCD dataset is presented in table 6. In order
to produce a comparison with further methods, that work on the same ground-
based cloud classification task, the settings described in Liu et al. (2020b)
have been adopted. Looking at the results, several conclusions can be drawn.
The proposed model, composed of different pretrained networks, produces
promising recognition accuracy, giving a high representation of clouds images.
The implemented voting based architecture compared with the competitors gets
the better performance.

Table 6: Experimental results on MGCD dataset. Method | Acc
---|---
Our | 99.98
MMFN Liu et al. (2020b) | 88.63
DCAFs + MI Liu et al. (2020b) | 82.97
BOVW + MI Liu et al. (2020b) | 67.20
PBOVW + MI Liu et al. (2020b) | 67.15
LPB+ MI Liu et al. (2020b) | 50.53
CLPB+ MI Liu et al. (2020b) | 69.68
CloudNet + MI Liu et al. (2020b) | 80.37
BoVW Csurka et al. (2004) | 66.15
PBoVW Csurka et al. (2004) | 66.13
LBP Ojala et al. (2002) | 55.20
CLBP Guo et al. (2010) | 69.18
VGG-16 Simonyan and Zisserman (2014) | 77.95
DCAFs Shi et al. (2017) | 82.67
CloudNet Zhang et al. (2018) | 79.92
DMF Liu and Li (2018) | 79.05
DTFN Li et al. (2020) | 86.48
HMF Liu et al. (2019) | 87.90

The classification performance on SWIMCAT dataset are summarized in table 7.
In this phase, the settings present in Shi et al. (2017) are adopted. In
particular, a cross validation on 2,3,4,5 folds was performed first.
Subsequently, 40 images per class for training and 45 ones for testing have
been selected randomly. The average accuracy of 50 random runs is reported. In
the training phase, the neural models that did not contribute to improve both
the performance and the execution time were discarded, as can be seen in table
5. The results highlight that combining different multiple classification
predictions is useful to capture more spatial and local layout information of
clouds with purpose to outperform the compared methods. Furthermore, it is
important to underline that the proposed approach is even better than neural
models, such as VGG-16 Simonyan and Zisserman (2014) and CloudNet Zhang et al.
(2018), in which a single classification confidence value is provided compared
to a multiple voting based mechanism.

Table 7: Experimental results on SWIMCAT dataset. Methods Folds | 2 | 3 | 4 | 5 | 40/45
---|---|---|---|---|---
Our | 99.36 | 99.49 | 99.49 | 99.75 | 99.91
LPB Sun et al. (2009) | 85.26 | 81.60 | 83.51 | 85.03 | 93.47
Heinle Feature Heinle et al. (2010) | 90.26 | 91.89 | 92.91 | 93.43 | 93.09
Text-based method Dev et al. (2015) | - | - | - | - | 95.00
DCAF Shi et al. (2017) | 98.72 | 98.46 | 98.97 | 98.84 | 99.56

Table 8 shows the performance on CCSN dataset. In order to compare results
with further methods that work on same task, the settings described in Zhang
et al. (2020) have been adopted. Also in this case, the combination of neural
models lead to better performance. As shown in table 5, for this experimental
phase we have stacked all analyzed deep neural network. Once again, the
results demonstrate that the multiple base learners may lead to better
performance according to the combination with the different number of base
learners in the stacking.

Table 8: Experimental results on CCSN dataset. Method | Acc
---|---
Our | 95.08
Cloudnet Zhang et al. (2018) | 90.00
Zhang et al. (2020) | 80.00
MMI Liu et al. (2018) | 75.42
M_DF Liu et al. (2018) | 78.21
M_JFCNNLiu et al. (2018) | 84.55
V_DF Liu et al. (2018) | 85.10
V_JFCNN Liu et al. (2018) | 86.79
V_DF + MMI Liu et al. (2018) | 86.33
V_JFCNN + MMI Liu et al. (2018) | 89.40
V_DF + M_DF Liu et al. (2018) | 90.21
J_JFCNN Liu et al. (2018) | 78.82
JFCNN Liu et al. (2018) | 93.37

The presented satisfactory results are attributable to many relevant aspects.
The first regards the features extracted through convolutional layers of the
deep neural network. They provide a good image representation, although are
completely abstract and devoid of real meaning. Second regards the framework
capability to provide multiple representation models, that lead to a
significant improvement in performance. Another issue concerns the image size
normalization, tackled to many methods in the field. It is performed before
the features extraction, to avoid performance degradation. Again, we can look
at the robustness with respect to the underrepresented classes in the
datasets. In fact, the framework does not fail even though the samples are not
sufficient for a class representation in specific cases. The latter appears to
be an open problem in the literature, as ad hoc classifiers are often designed
for unbalanced classification different from the standard ones that produce
untrue results. Contrary, a weak point concerns the computational aspect.
First, the time required for training pretrained model is high but less than a
model created from scratch. Second, the classification step, that provides
multiple choices in decision making at each iteration, requires a lot of
effort. The latter works for the purpose of choosing which classifiers are
suitable for specific clouds images included in the test set. Finally, we have
shown that although the framework is more expensive from a computational point
of view and it produces better results than a single classifier.

##  5 Conclusions and Future Works

The challenge in ground-based cloud recognition is specifically interesting
and, not only, for its multiple aspects and variety of data. The complexity of
the task is linked to several factors such as the type of clouds and the
visual patterns contained in them. In support, convolutional neural networks
lend a big hand to understand the meaning of images with the consequent goal
of their classification. In this regard, we proposed a framework that combines
convolutional neural networks, adapted to the cloud recognition task through a
transfer learning approach, using voting rules. The results produced certainly
strengthen the theoretical thesis. A multiple model, based on several deep
neural networks, compared to a single one is a powerful factor. Through a
large experimental phase, it has been shown how the proposed approach is
competitive, and in some cases better, compared to the more advanced methods.
Although pretrained models have been adopted, the main weakness concerns the
computational complexity of learning phase that requires a long time,
sensitive to the growth of the data. Future work will certainly concern the
study and analysis of still unexplored convolutional neural networks for this
type of problem and the application of the proposed framework to further
datasets with the aim of taking a step forward in cloud recognition.

## Acknowledgements

Our thanking is for Alfredo Petrosino. He followed us during the first steps
towards the Computer Science, through a whirlwind of goals, ideas and,
especially, love and passion for the work. We will be forever grateful great
master.

## References

  * Chen et al. (2000) Chen, T., Rossow, W.B., Zhang, Y., 2000.  Radiative effects of cloud-type variations.  Journal of climate 13, 264–286.
  * Csurka et al. (2004) Csurka, G., Dance, C., Fan, L., Willamowski, J., Bray, C., 2004.  Visual categorization with bags of keypoints, in: Workshop on statistical learning in computer vision, ECCV, Prague. pp. 1–2.
  * Deng et al. (2009) Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L., 2009\.  Imagenet: A large-scale hierarchical image database, in: 2009 IEEE conference on computer vision and pattern recognition, Ieee. pp. 248–255.
  * Dev et al. (2015) Dev, S., Lee, Y.H., Winkler, S., 2015.  Categorization of cloud image patches using an improved texton-based approach, in: 2015 IEEE International Conference on Image Processing (ICIP), IEEE. pp. 422–426.
  * Dev et al. (2014) Dev, S., Savoy, F.M., Lee, Y.H., Winkler, S., 2014.  Wahrsis: A low-cost high-resolution whole sky imager with near-infrared capabilities, in: Infrared Imaging Systems: Design, Analysis, Modeling, and Testing XXV, International Society for Optics and Photonics. p. 90711L.
  * Duda et al. (2013) Duda, D.P., Minnis, P., Khlopenkov, K., Chee, T.L., Boeke, R., 2013.  Estimation of 2006 northern hemisphere contrail coverage using modis data.  Geophysical Research Letters 40, 612–617.
  * Guo et al. (2010) Guo, Z., Zhang, L., Zhang, D., 2010.  A completed modeling of local binary pattern operator for texture classification.  IEEE Transactions on Image Processing 19, 1657–1663.
  * He et al. (2016) He, K., Zhang, X., Ren, S., Sun, J., 2016.  Deep residual learning for image recognition, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770–778.
  * Heinle et al. (2010) Heinle, A., Macke, A., Srivastav, A., 2010.  Automatic cloud classification of whole sky images.  Atmospheric Measurement Techniques 3, 557–567.
  * Huang et al. (2017) Huang, G., Liu, Z., Van Der Maaten, L., Weinberger, K.Q., 2017.  Densely connected convolutional networks, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 4700–4708.
  * Krizhevsky et al. (2012) Krizhevsky, A., Sutskever, I., Hinton, G.E., 2012.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems 25, 1097–1105.
  * LeCun et al. (1989) LeCun, Y., Boser, B., Denker, J.S., Henderson, D., Howard, R.E., Hubbard, W., Jackel, L.D., 1989.  Backpropagation applied to handwritten zip code recognition.  Neural computation 1, 541–551.
  * Li et al. (2020) Li, M., Liu, S., Zhang, Z., 2020.  Deep tensor fusion network for multimodal ground-based cloud classification in weather station networks.  Ad Hoc Networks 96, 101991\.
  * Liu et al. (2019) Liu, S., Duan, L., Zhang, Z., Cao, X., 2019.  Hierarchical multimodal fusion for ground-based cloud classification in weather station networks.  IEEE Access 7, 85688–85695.
  * Liu and Li (2018) Liu, S., Li, M., 2018.  Deep multimodal fusion for ground-based cloud classification in weather station networks.  EURASIP Journal on Wireless Communications and Networking 2018, 1–8.
  * Liu et al. (2020a) Liu, S., Li, M., Zhang, Z., Cao, X., Durrani, T.S., 2020a.  Ground-based cloud classification using task-based graph convolutional network.  Geophysical Research Letters 47, e2020GL087338.
  * Liu et al. (2018) Liu, S., Li, M., Zhang, Z., Xiao, B., Cao, X., 2018\.  Multimodal ground-based cloud classification using joint fusion convolutional neural network.  Remote Sensing 10, 822\.  doi:10.3390/rs10060822.
  * Liu et al. (2020b) Liu, S., Li, M., Zhang, Z., Xiao, B., Durrani, T.S., 2020b.  Multi-evidence and multi-modal fusion network for ground-based cloud recognition.  Remote Sensing 12, 464\.
  * Liu et al. (2017) Liu, W., Wang, Z., Liu, X., Zeng, N., Liu, Y., Alsaadi, F.E., 2017.  A survey of deep neural network architectures and their applications.  Neurocomputing 234, 11–26.
  * Ojala et al. (2002) Ojala, T., Pietikainen, M., Maenpaa, T., 2002.  Multiresolution gray-scale and rotation invariant texture classification with local binary patterns.  IEEE Transactions on Pattern Analysis and Machine Intelligence 24, 971–987.
  * Peteiro-Barral and Guijarro-Berdiñas (2013) Peteiro-Barral, D., Guijarro-Berdiñas, B., 2013.  A survey of methods for distributed machine learning.  Progress in Artificial Intelligence 2, 1–11.
  * Rossow and Schiffer (1991) Rossow, W.B., Schiffer, R.A., 1991\.  Isccp cloud data products.  Bulletin of the American Meteorological Society 72, 2–20.
  * Shi et al. (2017) Shi, C., Wang, C., Wang, Y., Xiao, B., 2017\.  Deep convolutional activations-based features for ground-based cloud classification.  IEEE Geoscience and Remote Sensing Letters 14, 816–820.
  * Simonyan and Zisserman (2014) Simonyan, K., Zisserman, A., 2014\.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556 .
  * Stephens (2005) Stephens, G.L., 2005.  Cloud feedbacks in the climate system: A critical review.  Journal of climate 18, 237–273.
  * Sun et al. (2009) Sun, X., Liu, L., Gao, T., Zhao, S., 2009.  Classification of whole sky infrared cloud image based on the lbp operator.  Transactions of Atmospheric Sciences 32, 490–497.
  * Szegedy et al. (2015) Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., Erhan, D., Vanhoucke, V., Rabinovich, A., 2015\.  Going deeper with convolutions, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1–9.
  * Zhang et al. (2020) Zhang, J., Liu, P., Zhang, F., Iwabuchi, H., de Moura, A.A., de Albuquerque, V.H.C., 2020.  Ensemble meteorological cloud classification meets internet of dependable and controllable things.  IEEE Internet of Things Journal .
  * Zhang et al. (2018) Zhang, J., Liu, P., Zhang, F., Song, Q., 2018.  Cloudnet: Ground-based cloud classification with deep convolutional neural network.  Geophysical Research Letters 45, 8665–8672.
  * Zoph et al. (2018) Zoph, B., Vasudevan, V., Shlens, J., Le, Q.V., 2018\.  Learning transferable architectures for scalable image recognition, in: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 8697–8710.