aainstitutetext: Computation-based Science and Technology Research Center, The
Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus.bbinstitutetext:
University of Cyprus, 1 Panepistimiou Avenue 2109 Aglantzia, Nicosia.

# Unsupervised identification of the phase transition on the 2D-Ising model

Constantia Alexandrou a Andreas Athenodorou a Charalambos Chrysostomou a
Srijit Paul<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We investigate deep learning autoencoders for the unsupervised recognition of
phase transitions in physical systems formulated on a lattice. We use spin
configurations produced for the 2-dimensional ferromagnetic Ising model in
zero external magnetic field. We study numerically the relation between one
latent dimension extracted from the autoencoder to the critical temperature
$T_{c}$. The autoencoder reveals the two phases, one for which the spins are
ordered and the other for which spins are disordered, reflecting the
restoration of the $\mathbb{Z}_{2}$ symmetry as the temperature increases. For
the largest volume studied, the transition between the two phases occurs very
close to the theoretically extracted critical temperature. We define as a
quasi-order parameter the absolute average latent dimension ${\tilde{z}}$,
which enables us to predict the critical temperature. We show that one can
build the latent susceptibility and use it to quantify the value of the
critical temperature $T_{c}(L)$ at different lattice sizes and that these
values suffer from only small finite scaling effects. We demonstrate that
$T_{c}(L)$ extrapolates to the known theoretical value as $L\to\infty$
suggesting that the autoencoder can also be used to extract the critical
temperature of the phase transition to an adequate precision.

##  1 Introduction

Recent advances in the implementation of Artificial Intelligence (AI) for
physical systems, especially, on those which can be formulated on a lattice,
appear to be suitable for observing the corresponding underlying phase
structure Carrasquilla2017 ; vanNieuwenburg2017 ; PhysRevB.94.195105 ;
Broecker2017 ; ch2017machine ; PhysRevE.96.022140 ; PhysRevE.95.062122 ;
2017arXiv170700663B ; PhysRevB.97.205110 ; PhysRevLett.120.257204 ;
2018arXiv180402709Z ; 2018JChPh.149s4109J ; PhysRevE.98.022138 ;
2018arXiv180801731Z ; Kashiwa:2018jdi ; Giannetti:2018vif . So far methods
such as the Principal Component Analysis (PCA) PhysRevB.94.195105 ;
vanNieuwenburg2017 ; PhysRevE.95.062122 ; wetzel2017unsupervised , Supervised
Machine Learning Broecker2017 ; 2018arXiv180402709Z ; morningstar2017deep
(ML) as well as autoencoders PhysRevE.95.062122 ; PhysRevE.96.022140  appear
to successfully identify different phase regions of classical statistical
systems, such as the 2-dimensional (2D) Ising model that describes the
ferromagnetic-paramagnetic transition. These techniques were also applied on
quantum statistical systems, such as the Hubbard model ch2017machine  that
describes the transition between conducting and insulating systems. Very
recently, similar studies have been applied for simulations of quantum fields
on the lattice, such as the $SU(2)$ gauge theory Wetzel:2017ooo  with an
increased complexity in the data due to the structure of the $SU(2)$ gauge
group.

Trained neural networks can thus help distinguish phases in simple statistical
systems the structure of which is known, but more importantly in more complex
systems where the underlying phase structure is unknown. In this work, we
would like to examine whether the proposed, fully-connected (Dense), deep
learning autoencoder, which does not require supervised training, can shed
light on the phase structure of the 2D-Ising model. Deep learning autoencoders
are frequently used in cases where data hides interesting structure by
processing the raw datasets. They can, therefore, be used to discover
interesting structure in ensembles produced for a range of a parameter that
characterizes the phase space of the model in different sectors with different
physical properties. One such example is the ferromagnetic Ising model for
which at the critical temperature $T_{c}$, the system undergoes a transition
from the ordered phase to the disordered.

In this work, we investigate the action of unsupervised machine learning,
namely the deep learning autoencoder (not variational), towards the
identification of the phase transition of the 2D-Ising model. More
specifically, we produce decorrelated configurations for the 2D-Ising model
for a given range of temperatures and then we apply the autoencoder trying to
understand what information of the phase structure we can capture. Hence,
technically, this work combines the production of configurations using Monte
Carlo methods as well as the deep learning autoencoder algorithm. We observe
that the autoencoder can capture the underlying $\mathbb{Z}_{2}$ symmetry and
can indeed find out where the transition occurs by identifying as a relevant
quasi-order parameter the mean value of the absolute latent dimension.
Although this quantity is not suitable for predicting the order of the
transition, it can determine the critical temperature with small finite
scaling effects.

This article is organised as follows: In section 2 we present a brief
description of the 2D-Ising model, explaining the production of the
configurations as well as its phase structure. In section 3 we discuss the
deep learning autoencoder, explain how it works and provide the structure of
the network. Subsequently, in section 4 we provide our results, and finally,
in section 5, we conclude.

##  2  The Ferromagnetic 2-Dimensional Ising Model

One of the most interesting physical phenomena in nature is magnetism. It is
known that the ferromagnetic materials exhibit a spontaneous magnetisation in
the absence of an external magnetic field. Such magnetisation occurs only if
the temperature of the system is lower than a known critical temperature
$T_{c}$, the so called Curie temperature. If the temperature of the system is
raised so that $T>T_{c}$, then the magnetisation vanishes. In principle, the
critical temperature $T_{c}$ separates the microstates of the system from
being ordered or magnetised for $T<T_{c}$ to being randomly oriented resulting
in zero magnetisation; these two phases correspond to the ferromagnetic and
the disordered phases, respectively.

Ferro-magnetism has a quantum mechanical nature and, thus, a lot of effort is
invested towards its understanding. Albeit quantum mechanical, simple
classical models can help to gain insight into this effect. The 2D-Ising model
is a classical model that is commonly used to study magnetisation. The
2D-Ising model can be considered as a lattice with $N=N_{x}\times N_{y}$
sites, on each of which a double valued spin $s_{i}$ is located, either being
in an "up" orientation denoted by $\uparrow$ or $s_{i}=+$ or "down" denoted by
$\downarrow$ or $s_{i}=-$.

The macroscopic properties of the 2D-Ising system are determined by the nature
of the accessible micro-states. Thus, it is useful to know the dependence of
the Hamiltonian on the spin configurations. The total energy is given by

$\displaystyle H=-J\sum^{N}_{i,j=nn(i)}s_{i}s_{j}-\mu h\sum^{N}_{i=1}s_{i}\,,$
(1)

where $J$ is the self-interaction between neighbouring spins, $h$ the external
magnetic field and $\mu$ is the atomic magnetic moment. Note that in the first
sum, the notation $nn(i)$ represents nearest-neighbour pairs; the sum is taken
over all nearest-neighbouring pairs.

In the case of the canonical ensemble, in other words, when the system is
attached to a thermal reservoir and kept at a constant temperature $T$, as the
time passes the spins are left to fluctuate with rates depending on the
reservoir’s temperature. This behaviour can be captured in a Monte Carlo (MC)
simulation in the canonical ensemble. For $T=0$ the system is frozen with all
spins being at one direction either down or up. The orientation of the spins
is arbitrary, however, the dynamics enforce the system to choose one of the
two directions. This corresponds to the spontaneous symmetry breaking of the
$\mathbb{Z}_{2}$ global symmetry group. Although the Hamiltonian of the system
is invariant under $\mathbb{Z}_{2}$ transformations, the degenerate ground
states are not invariant but get interchanged under such transformations.

For small finite values of the temperature of the system, spins still form
large sectors where all spins are correlated and point to one direction. Above
the critical temperature $T_{c}$, the spins are disordered and
$\mathbb{Z}_{2}$ symmetry is restored.

The question that we address in this work is whether the behaviour described
above can be captured by a deep learning autoencoder when we pass it ensembles
for a sequence of temperatures separated by some $\delta T$. More precisely,
we seek to understand if a qualitative description of the phase structure of
the Ising model can be extracted and whether one can determine the critical
temperature $T_{c}$.

###  2.1  Swendsen-Wang algorithm

The MC simulation for the 2D-Ising model is conventionally performed using the
Metropolis algorithm. Since this algorithm is based on local updates, near the
critical temperature where the correlation length diverges, it faces the
problem of critical slowing down. In order to tackle this problem, we have
implemented the Swendsen-Wang cluster algorithm Swendsen:1987ce ; WANG1990565
, which is based on global updates of the spin configurations. This algorithm
relies on the formation of bonds between every pair of nearest
neighbours($ij$) that are aligned at a given temperature $T$, with a
probability $p_{ij}=1-\exp{(-2\beta J)}$, where $\beta=\frac{1}{k_{B}T}$
($k_{B}\equiv$ Boltzmann constant). A single cluster is defined as all the
spins, which are connected via bonds. The global update is defined as the
collective flipping with a probability of $1/2$, on all the spins in each
cluster 1969PSJJS..26…11K ; Fortuin:1971dw . This step works because of the
so-called Fortuin-Kasteleyn mapping of the Ising model on the random-cluster
model. Thus, global updates enable us to produce equilibrium configurations
close to the $T_{c}$ with a few thermalisation steps.

###  2.2  Monte-Carlo simulation setup

In this work we chose to investigate the case of zero external magnetic field
($h=0$) and for simplicity we have set $J=1$ and $k_{B}=1$. In this case, the
theoretically calculated value of the critical temperature is

$\displaystyle T_{c}=\frac{2}{\ln\left(1+\sqrt{2}\right)}=2.269185\,.$ (2)

To extract experimentally this quantity one has to investigate the order
parameter of theory, namely the magnetization. The first question that we
address is whether we can get an approximate estimate of this temperature by
using unsupervised learning. To this purpose, we choose a sequence of
different values of temperature, and for each one, we start from a frozen
configuration of spins, perform a large enough number of thermalisation sweeps
and then save the configuration. For every single temperature, we repeat the
procedure $200$ times.

###  2.3  Phase structure, observables and order parameters

The phase structure of the 2D-Ising model can be reduced to the study of the
magnetic order of the system cardy_1996 ; gould1988introduction . If we
suppose that there are $N_{\uparrow}$ spins pointing upwards and
$N_{\downarrow}$ spins pointing downwards, then the total magnetic moment
would be $N_{\uparrow}-N_{\downarrow}$ ($\mu=1$). The largest possible
magnetic moment would, therefore, be $N$. Thus, we can define the magnetic
order parameter or magnetization per spin configuration naturally as:

$\displaystyle m=(N_{\uparrow}-N_{\downarrow})/N\,,$ (3)

while the average magnetization $M=\langle m\rangle$. $M$ can get values
between $-1$ and $1$, and the average of the absolute magnetization
${\tilde{m}}=\langle|m|\rangle$ is just the magnetic order. Hence, if
${\tilde{m}}$ is close to $0$, then the system is highly disordered and, thus,
not magnetised, with approximately half of the spins pointing up and the other
half pointing down. On the other hand, if ${\tilde{m}}$ is approximately $1$,
the system is ordered and, thus, magnetised with nearly all the spins pointing
in the same direction.

The point $T=T_{c}$ is called the critical point and separates the ordered
$T<T_{c}$ phase and disordered $T>T_{c}$ phase. At $T=T_{c}$ the system is
described by a second order phase transition, i.e. à la Ehrenfest
jaeger1998ehrenfest  the first derivative of the free energy with respect to
the external field which is the order parameter is continuous while the second
derivative of the free energy is discontinuous.

##  3 Deep Learning Autoencoders

Autoencoders is a variety of artificial neural networks utilized for learning
data codings in an unsupervised manner, efficiently vincent2008extracting ;
vincent2010stacked . An autoencoder aims to define a representation (encoding)
for an assemblage of data, usually performing dimensionality reduction. An
autoencoder encodes the input data ($\\{X\\}$) from the input layer into a
latent dimension ($\\{z\\}$), and then uncompresses that latent dimension into
an approximation of the original data ($\\{\,X\\}$). This drives the
autoencoder to engage in dimensionality reduction, by learning how to ignore
the noise and recognise significant characteristics of the input data. The
first layer of an autoencoder might learn to encode simple, identifiable and
local features, and the second layer by using the output of the first layer
learns to encode more complex and less local features, until the final layer
of the encoder learns to identify and encode the most complex and global
characteristics of the input data. As Fig. 1 shows, an autoencoder consists of
two components, the encoder function $g_{\phi}$ and a decoder function
$f_{\theta}$ and the reconstructed input is $\,X=f_{\theta}(g_{\phi}(x))$.

In the training phase, the autoencoder learns the parameters $\phi$ and
$\theta$ together, where $f_{\theta}(g_{\phi}(x))$ can approximate an identity
function. Various metrics can be used to measure the error between the
original input $X$ and the reconstruction ${\tilde{X}}$, but the most simple
and most commonly used is the Mean Square Error (MSE) as this is provided in
Eq. 4, where $n_{\rm data}$ is the number of data points:

$MSE(\theta,\phi)=\frac{1}{n_{\rm data}}\sum^{n_{\rm
data}}_{i=1}(X_{i}-f_{\theta}(g_{\phi}(X_{i})))^{2}\,.$ (4)

Figure 1: Basic structure of an autoencoder

###  3.1 Proposed Autoencoder Model

For the analysis of the proposed method, an eight-layer, fully connected
(Dense), autoencoder is proposed, as Fig. 2 shows, where the encoder
compresses the configurations into a single latent dimension. Through
experimentation, we determine that the best model to detect the transition
consists of the encoder with the input layer, first, second and third hidden
layers having 625, 256, 64 and 1 neurons, respectively. The activation
function used is relu, as shown in Eq. 5, for all layers except the third
hidden layer, where tanh was used, as shown in Eq. 6. For the decoder, the
first, second and third hidden layers use 64, 256, and 625 neurons,
respectively. For the output layers, the number of neurons is set to be equal
to the number of values in the configuration under investigation. The
activation function used is relu, as given in Eq. 5, for all hidden layers,
and for the output layer, tanh is used, as per Eq. 6.

${\rm relu}:y=max(0,x)=\left\\{\begin{array}[]{ll}x,&\mathrm{if}\ x>0\\\
0&\mathrm{if}\ x\leq 0\end{array}\right\\}.\ $ (5)

${\rm tanh}:y=\frac{1-e^{-2x}}{1+e^{-2x}}\,.$ (6)

Figure 2: Proposed autoencoder model in the standard Machine Learning
nomenclature.

For the proposed autoencoder model we use the so-called dropout realization
technique hinton2012improving . The dropout regularization technique refers to
dropping out neurons from each layer, randomly, when training. Dropout is
successfully used for reducing over-fitting in neural networks by preventing
complex co-adaptations on training data. For the training of the proposed
autoencoder model the data are split into training (66.66…%) and testing
(33.33…%) sets, and the training is performed for 2000 iterations. The
implementation was performed using Keras chollet2015keras  and Tensorflow
abadi2016tensorflow .

##  4 Results

###  4.1 The latent dimension per configuration

Each configuration is re-expressed in the form of a vector, and then it is
read as an input by the autoencoder. One can think of the input as a column
with entries of 1 and -1, placed in some lexicographic order and having a
length equal to $L^{2}$. More precisely, for each different lattice size
$L=N_{x}=N_{y}$ we feed to the network all the configurations produced for all
different temperatures, and we extract the latent dimension $z_{i_{\rm
conf}}$. In other words, each configuration is assigned a number, the latent
dimension, which includes all the physically necessary information so that the
decoder re-creates the actual configuration. It should be made clear that
configurations for different lattice volumes have been fed separately into the
autoencoder, i.e. the autoencoder receives information for only one lattice
volume, and thus, it "knows" nothing about configurations produced for other
volume sizes.

In order to identify signals of the phase structure of the 2D-Ising model, as
a first step, we investigate how the latent dimension $z_{i_{\rm conf}}$
behaves as a function of the temperature $T$ for each configuration. We
produce $40000$ configurations, namely $200$ configurations for every single
temperature. The produced configurations are for $200$ different values of
temperatures within the range $T=1-4.5$ and separated by $\delta T=0.0175$. We
make sure that we cover the whole range of temperatures between the two
extreme cases of the Ising behaviour, the nearly "frozen" at $T\simeq 1$, and
the complete disordered $T\simeq 4.5$. Furthermore, we assume that we have no
prior knowledge on what is happening in between these two extremes. We note
that we could choose different temperature ranges that cover all possible
phase regions; for instance, we could choose instead $T=0.01-1000$ with
$\delta T=0.01$, but of course the computational effort would be much more
significant.

In Fig. 3 we show the latent dimension for each different configuration, as a
function of the temperature $T$, for four different lattice sizes, $L=25$,
$35$, $50$, $150$.

Figure 3: The latent dimension for each configuration as a function of the
temperature for four different lattice volumes. The dashed line represents the
analytically extracted value of the critical temperature (Eq. 2). The red
shaded area in the plot for $L=150$ is the region where (by fitting to a
constant) we expect to find the $T_{c}(L=150)$. The color on the gradient
illustrator on the right denotes the temperature $T$.

Fig. 3 has the following features:

  * •

For low temperatures we obtain two plateaus, one located at $z=1$ and one at
$z=-1$. A first simplistic explanation for this pattern would be that it
corresponds to two distinct states that are not connected through any kind of
transformation. This reflects the spontaneously broken
$\mathbb{Z}_{2}\equiv\\{-1,1\\}$ global symmetry group. One can interpret
these two plateaus as the two cases where all spins are up or down. This
interpretation is confirmed by the results presented in Fig. 4 where we show
the absolute correlation coefficient $C_{z,m}$ between the latent dimension
$z$ and the magnetisation $m$ defined as

$\displaystyle
C_{z,m}=\frac{\langle(z-{\bar{z}})(m-{\bar{m}})\rangle}{\sqrt{\langle(z-{\bar{z}})^{2}\rangle\langle(m-{\bar{m}})^{2}\rangle}}\,.$
(7)

The fact that at low temperatures the absolute correlation coefficient is $1$
demonstrates that the two different values of the latent dimension $-1$ and
$1$ correspond to the two orientations of the spins. Finally, the two plateaus
become more distinct as the lattice size increases.

  * •

At some temperature range $\Delta T_{\rm trans}$ the aforementioned behaviour
collapses to one state, which is located around $z=0$. This reflects the
restoration of $\mathbb{Z}_{2}$ symmetry. In other words, it corresponds to
the case where all the spins are disoriented.

  * •

There is a critical point where there is a change in the pattern. As the
lattice size increases the width of this transition decreases with $\Delta
T_{\rm trans}\to 0$ and this step becomes steeper and steeper. At $L\to\infty$
the transition is localised right on the critical temperature $T_{c}$
extracted analytically.

Figure 4: The absolute correlation coefficient defined in Eq. 7 as a function
of the temperature for the five different lattice sizes.

Evidently, plotting the latent dimension as a function of the temperature
demonstrates that the autoencoder predicts the two phases. Also, it provides
to a good approximation the critical temperature. In fact from Fig. 3 the
transition appears to occur right at the critical temperature for $L=150$. By
fitting the points of the latent dimension which, to a good approximation,
behave linearly to a constant as a function of $T$ we can restrict that the
collapse of the two states located at $1$ and $-1$ occurs at $T\simeq
2.28(4)$. This temporal region is denoted in Fig. 3 as the shaded area in red.
In Fig. 4 we observe that within this temporal region the value of the
absolute correlation coefficient $C_{z,m}$ starts to decrease from $1$. This
demonstrates that, although highly correlated, the latent dimension and
magnetisation are two different quantities.

Finally, we observed that for low temperatures the latent dimension is, in a
good approximation, equally distributed between the values of $z=1$ and
$z=-1$. This can be seen in Fig. 5 where we present the average latent
dimension $\langle z\rangle$ for each temperature and $L=150$ as a function of
the temperature.

Figure 5: The average latent dimension $\langle z\rangle$ as a function of the
temperature $T$ for $L=150$.

One could also investigate what happens within different "temperature windows"
introducing, however, a degree of supervision. For instance, we can use a
temperature window within the range $T=1-2$ and apply the autoencoder. The
outcome would be the behaviour presented in the left panel of Fig. 6, where
only the two ordered states are visible without the presence of a critical
point. Since there is no visible signal for a phase transition behaviour
within this range of $T$ it is reasonable to use another temperature window.
If we choose $T=3-4.5$, for instance, the corresponding latent dimension would
be the one given on the right panel of Fig. 6 where no particular pattern is
observed. A sensible next step would be to investigate what happens within a
range of temperatures located between the two previous temperature windows,
for instance $T=1-4$.

$T$$z_{i_{\rm conf}}$$T$$z_{i_{\rm conf}}$

Figure 6: The latent dimension for each configuration as a function of the
temperature for $L=100$, obtained by applying the autoencoder on
configurations produced within the range of temperatures $T=1-2$ (left) and
$T=3-4.5$ (right).

###  4.2 The absolute average latent dimension

Since the latent dimension per configuration is symmetric with respect to the
$T$ axis, it would be reasonable to define the average absolute latent
dimension as a parameter indicating the phase as

$\displaystyle{\tilde{z}}=\frac{1}{N_{\rm conf}}\sum_{i=1}^{N_{\rm
conf}}\left|z_{i_{\rm conf}}\right|\,.$ (8)

Fig. 3 shows that the latent dimension resembles the behaviour of the
magnetization per spin configuration as a function of the temperature. The
absolute average magnetization defines the order parameter of the system
distinguishing the two different phases. For the case of the autoencoder we
can define an additional quasi-order parameter as the absolute average latent
dimension.

In the left-hand-side of Fig. 7 we provide the magnetisation as a function of
the temperature while on the right-hand side we provide the absolute latent
dimension. Indeed the absolute latent dimension looks similar to the
magnetisation, albeit becoming steeper as the lattice size increases. Clearly,
the magnetization behaves as an order parameter with the characteristics of a
second order phase transition while the absolute latent dimension is
consistent with a first order phase transition. We can, therefore, conclude
that the absolute average latent dimension can be used as an order parameter
to identify the critical temperature, but cannot capture the right order of
the phase transition. The fact that ${\tilde{z}}$ as a function of the
temperature becomes steeper as the lattice size increases suggests that the
critical temperature $T_{c}(L)$ as a function of the lattice size $L$
extracted from the autoencoder data will suffer less from finite-size scaling
effects as discussed in detail in Section 4.3.

Figure 7: The average latent dimension as a function of the temperature for
five different lattice volumes.

Traditionally, $T_{c}(L)$ can be extracted by probing the peak of the magnetic
susceptibility $\chi$ at zero magnetic field $h$, where

$\displaystyle\chi=\frac{L^{2}}{T}\left(\langle m^{2}\rangle-\langle
m\rangle^{2}\right)\,.$ (9)

According to finite size scaling theory, close enough to $T_{c}$, magnetic
susceptibility $\chi$ scales as

$\displaystyle\chi\propto\left(t\right)^{-\gamma}\,,$ (10)

where $t=(T-T_{c})/T_{c}$ is the reduced temperature and $\gamma=7/4$ a
critical exponent cardy_1996 . The magnetic susceptibility measures the
ability of a spin to respond due to a change in the external magnetic field.
In the same manner we define the latent susceptibility as

$\displaystyle\chi_{\tilde{z}}=\frac{L^{2}}{T}\left(\langle{\tilde{z}}^{2}\rangle-\langle{\tilde{z}}\rangle^{2}\right)\,.$
(11)

For the extraction of $T_{c}$ one realizes, by looking at Fig. 7, that more
data points close to the critical behaviour are needed to extract the critical
temperature from the latent susceptibility. Hence, we produce configurations
for a grid of temperatures near the critical regime. More specifically, we
produce 200 configurations per temperature, for 200 different values of $T$ in
the range of $T=2-2.8$ and $\delta T=0.004$ for all the volumes considered. In
addition, for $L=100$ and $L=150$ we produce 200 configuration for each value
of $T$, in the range of $T=2.22-2.34$ with $\delta T=0.0006$. These new
configurations, however, are not used to train the autoencoder. Instead, we
use the synaptic weights extracted and predict the latent dimension for the
new configurations. Hence, this serves as a confirmation that our data do not
suffer from over-fitting.

Figure 8: The red points show the predicted latent dimension for the
configurations produced around the critical point $T_{c}$ while in black the
latent dimension extracted in the first set of configurations on which we have
trained the encoder. The blue vertical line corresponds to the analytically
extracted critical temperature $T_{c}$.

In Fig. 8 we present the results of applying the autoencoder weights on the
new configurations produced in the region close to the critical point. We
compare with results extracted using configurations produced in the range
$T=1-4.5$. Both datasets agree, and there is a nice continuation of the
behaviour of the absolute average latent dimension within the critical regime.
This demonstrates that the execution of the encoder does not suffer from any
over-fitting occurrence and at the same time more data points can be used for
the extraction of $\chi_{\tilde{z}}$. Furthermore, the plot for $L=150$,
behaves nearly as a step function with the step being right on the
theoretically extracted $T_{c}$. By fitting the second moment of the latent
dimension, as this is described in section 4.3, one sees that the transition
occurs at $T_{c}(L=150)=2.2779(3)$; this value is very close to the
theoretically extracted value $T_{c}=2.26918$.

In the following section we present the analysis of our data in order to
investigate the latent susceptibility $\chi_{\tilde{z}}$ and, to subsequently,
extract the critical temperature $T_{c}(L)$ from the corresponding peak.

Figure 9: The latent susceptibility (left) and magnetic susceptibility (right)
as a function of the temperature for five different lattice volumes. The blue
vertical line denotes the analytically extracted value for the critical
temperature (Eq. 2).

###  4.3 The Latent Susceptibility and the Critical Temperature

In the previous sections, we provided strong evidence that the latent
dimension, resulting from the proposed autoencoder, demonstrates the
underlying phase transition and that it can also be used as a rough estimate
for the critical temperature $T_{c}$. Nevertheless, as the finite lattice size
$L$ increases we need to make sure that $T_{c}(L)$ tends to the right limit
i.e. it convergences to the theoretically extracted value given in Eq. 2 as
$L\to\infty$.

To investigate the convergence of $T_{c}(L)$, we first extract $T_{c}(L)$ for
each different lattice size and then extrapolate to infinite $L$. $T_{c}(L)$
can be extracted by probing the peak of the latent susceptibility for each
$L$. The latent susceptibility as a function of the temperature for the five
different lattice sizes is presented on the left-hand side of Fig. 9. Unlike
the magnetic susceptibility, presented on the right-hand side of Fig. 9, the
latent susceptibility is much sharper with peaks being closer to the known
critical temperature $T_{c}$. This means that the critical temperature for
each $L$ is influenced by less finite-size scaling effects.

Our temporal grid is fine enough and enables an adequate extraction of the
$T_{c}(L)$ from the coordinates of Fig. 9. Hence, there is no need to use
multi-histogram reweighting Giannetti:2018vif  techniques. The latent
dimension behaves to a large extent as a step function, and thus, tends to
$\propto\delta(T-T_{c})$ as $L\to\infty$. In addition, the derivative of the
latent susceptibility appears to be continuous. So we can also use a Gaussian
fit to estimate the critical temperature.

Figure 10: The critical temperature $T_{c}(L)$ extracted from fitting the
magnetic (red) and the latent (blue) susceptibilities as a function of $1/L$
to Eq. 12.

In Fig. 10 we present $T_{c}(L)$ extracted from fitting the latent
susceptibility and the magnetic susceptibility as a function of $1/L$. Results
obtained using the latent susceptibility suffer less from finite-size scaling
effects as compared to those when using the magnetic susceptibility. Adopting,
the usual finite-size scaling behaviour

$\displaystyle T_{c}(L)-T_{c}(L=\infty)\propto L^{-1/\nu}\,,$ (12)

we fit both susceptibilities to the ansatz $T_{c}(L)=T_{c}(L=\infty)+\alpha
L^{-1/\nu}$. Our findings are listed in Table 1.

Susceptibility | $T_{c}(L=\infty)$ | $\nu$ | $\chi^{2}/{\rm dof}$
---|---|---|---
Magnetic | $2.265(8)$ | $1.08(20)$ | $0.15$
Latent | $2.266(4)$ | $1.60(14)$ | $0.41$
Table 1: The results for $T_{c}(L=\infty)$ and $\nu$ extracted by fitting the
magnetic as well as the latent susceptibilities to the ansatz
$T_{c}(L)=T_{c}(L=\infty)+\alpha L^{-1/\nu}$.

As expected, fitting the data for $T_{c}(L)$ resulting from the magnetic
susceptibility yields values of $T_{c}(L=\infty)$ and $\nu$ which are
consistent with the analytically extracted values $T_{c}=2.269184$ and
$\nu=1$. Turning now to the case of the latent dimension, it appears that the
results of $T_{c}(L)$ when fitted with a form of the known scaling behaviour
of Eq. 12, yield a value for $T_{c}(L=\infty)$, which is in accordance with
the theoretical expectation. This provides a good evidence that the deep
learning autoencoder does not only predict the phase regimes of the 2D-Ising
model as well as give an estimate for the critical temperature, but can also
lead to a precise evaluation of the critical temperature.

##  5 Conclusions and Outlook

In this work, we apply a deep learning auto-encoder on configurations produced
for the 2D Ferromagnetic Ising model for performing a classification in an
unsupervised manner. Hence, with no prior knowledge on the system, we
demonstrate that we can predict the phase structure of this system
qualitatively as well as quantitatively by determining both phase regions and
the critical temperature.

At low temperatures, by making use of the latent dimension per configuration,
the autoencoder predicts two states reflecting to the broken $\mathbb{Z}_{2}$
symmetry. As the temperature increases, these two states appear to collapse at
one state, located around zero, and the underlying symmetry is restored. This
behaviour becomes more distinct as the volume of the lattice increases and the
point where the two states collapse is getting more and more local; this
corresponds to the critical point of the phase transition.

One can define the average absolute latent dimension ${\tilde{z}}$ that
displays partially the characteristics of an order parameter; namely, it can
identify the phase but cannot capture the order of the phase. Although it
resembles the behaviour of the magnetization it becomes steeper as the size of
the volume increases, tending to a step function. The second moment of the
absolute latent dimension defines a susceptibility, named latent
susceptibility, the peak of which can determine the critical temperature
$T_{c}(L)$. By extrapolating the values of $T_{c}(L)$ to $L\to\infty$ for the
sequence of lattice sizes $L=25,\ 35,\ 50,\ 100,\ 150$, we obtain for
$T_{c}(L=\infty)=2.266(4)$ in agreement with the exact value of
$T_{c}=2.26918$ calculated analytically. This suggests that the proposed deep
learning (fully-connected) autoencoder can identify, in an unsupervised
manner, the phase structure of the 2D-Ising model but can also lead to a
precise extraction of the critical temperature at the limit of the infinite
volume. As shown in Fig. 10 the values of $T_{c}(L)$ suffer with less finite
size effects compared to those usually extracted by using the peak of the
magnetic susceptibility, and one would thus expect that the autoencoder could
give a more precise prediction for $T_{c}$. Of course to test this hypothesis
we need to extract $T_{c}(L)$ for larger volumes, for instance up to $L=1024$
similarly to Ref. Giannetti:2018vif , and obtain the extrapolated value of
$T_{c}(L=\infty)$. This requires the usage of a different autoencoder with
more layers since memory limitations make the current autoencoder insufficient
to work. This is a future extension of this work.

There are other several related directions in which this work can be extended.
Since our proposed autoencoder has been tested just in one system, it would be
important to investigate its generalisation to other physical systems with
non-trivial phase structure. An important question, which could be answered is
whether this neural network is capable of identifying the phases for cases in
which an order parameter is either not known or not existing; such an example
is the Hubbard model 1963RSPSA.276..238H  describing the transition between
conducting and insulating systems. Another relevant question is how the
autoencoder behaves in cases where the phase transition is of a different
order or an infinite order such as in the 2D $XY$ spin model where the
relevant phase transition is the Kosterlitz-Thouless which is of infinite
order Kosterlitz:1973xp . Finally, our future plans involve the testing of the
autoencoder as a tool for the unsupervised extraction of the phase structure
of physical systems with continuous symmetries. These involve quantum field
theories formulated on the lattice such as the 3D $\phi^{4}$ with $O(2)$
symmetry Arnold:2001ir  where the phase transition is of second order and
belongs to the same universality class as the 2D-Ising model, the 3D $U(1)$
gauge theory Athenodorou:2018sab  for which the phase transition is of
infinite order and belongs to the same universality class as the 2D $XY$
model, as well as the 3D $SU(N)$ gauge theory Liddle:2008kk  which has a
second-order phase transition for $N\leq 3$, a weakly first order for $N=4$
and first order for $N\geq 5$.

## Acknowledgements

We would like to thank Giannis Koutsou, Nikos Savva, Spyros Sotiriadis, Mike
Teper and Savvas Zafeiropoulos for fruitful discussions. We would also like to
express our gratitude to Barak Bringoltz, Biagio Lucini and Davide Vadacchino
for performing a critical reading of the manuscript. AA was financially
supported by the Horizon 2020 of the European Commission research and
innovation programmes VI-SEEM under grant agreement No 675121 as well as the
OpenSESAME under grand agreement No 730943. CC also acknowledges support from
the OpenSESAME project. SP is supported by the Horizon 2020 of the European
Commission research and innovation programmes HPC-LEAP under the Marie
Skłodowska Curie grant agreement No. 642069 and PRACE-5IP under grant
agreement No 730913.

## References

  * (1) J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nature Physics 13 (Feb, 2017) 431 EP –, [arXiv:1605.01735].
  * (2) E. van Nieuwenburg, Y.-H. Liu, and S. Huber, Learning phase transitions by confusion, Nature Physics 13 (Feb, 2017) 435 EP –, [arXiv:1610.02048].
  * (3) L. Wang, Discovering phase transitions with unsupervised learning, Phys. Rev. B 94 (Nov, 2016) 195105, [arXiv:1606.00318].
  * (4) P. Broecker, J. Carrasquilla, R. G. Melko, and S. Trebst, Machine learning quantum phases of matter beyond the fermion sign problem, Scientific Reports 7 (2017), no. 1 8823, [arXiv:1608.07848].
  * (5) K. Ch’ng, J. Carrasquilla, R. G. Melko, and E. Khatami, Machine learning phases of strongly correlated fermions, Physical Review X 7 (2017), no. 3 031038, [arXiv:1609.02552].
  * (6) S. J. Wetzel, Unsupervised learning of phase transitions: From principal component analysis to variational autoencoders, Phys. Rev. E 96 (Aug, 2017) 022140, [arXiv:1703.02435].
  * (7) W. Hu, R. R. P. Singh, and R. T. Scalettar, Discovering phases, phase transitions, and crossovers through unsupervised machine learning: A critical examination, Phys. Rev. E 95 (Jun, 2017) 062122, [arXiv:1704.00080].
  * (8) P. Broecker, F. F. Assaad, and S. Trebst, Quantum phase recognition via unsupervised machine learning, arXiv e-prints (July, 2017) arXiv:1707.00663, [arXiv:1707.00663].
  * (9) N. Yoshioka, Y. Akagi, and H. Katsura, Learning disordered topological phases by statistical recovery of symmetry, Phys. Rev. B 97 (May, 2018) 205110, [arXiv:1709.05790].
  * (10) J. Venderley, V. Khemani, and E.-A. Kim, Machine learning out-of-equilibrium phases of matter, Phys. Rev. Lett. 120 (Jun, 2018) 257204, [arXiv:1711.00020].
  * (11) W. Zhang, J. Liu, and T.-C. Wei, Machine learning of phase transitions in the percolation and XY models, arXiv e-prints (Apr., 2018) arXiv:1804.02709, [arXiv:1804.02709].
  * (12) R. B. Jadrich, B. A. Lindquist, and T. M. Truskett, Unsupervised machine learning for detection of phase transitions in off-lattice systems. I. Foundations, jcp 149 (Nov., 2018) 194109, [arXiv:1808.00084].
  * (13) D. Kim and D.-H. Kim, Smallest neural network to learn the ising criticality, Phys. Rev. E 98 (Aug, 2018) 022138, [arXiv:1804.02171].
  * (14) X. L. Zhao and L. B. Fu, Machine Learning Phase Transition: An Iterative Proposal, arXiv e-prints (Aug., 2018) arXiv:1808.01731, [arXiv:1808.01731].
  * (15) K. Kashiwa, Y. Kikuchi, and A. Tomiya, Phase transition encoded in neural network, arXiv:1812.01522.
  * (16) C. Giannetti, B. Lucini, and D. Vadacchino, Machine Learning as a universal tool for quantitative investigations of phase transition, arXiv:1812.06726.
  * (17) S. J. Wetzel, Unsupervised learning of phase transitions: From principal component analysis to variational autoencoders, Physical Review E 96 (2017), no. 2 022140, [arXiv:1703.02435].
  * (18) A. Morningstar and R. G. Melko, Deep learning the ising model near criticality, arXiv:1708.04622.
  * (19) S. J. Wetzel and M. Scherzer, Machine Learning of Explicit Order Parameters: From the Ising Model to SU(2) Lattice Gauge Theory, Phys. Rev. B96 (2017), no. 18 184410, [arXiv:1705.05582].
  * (20) R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58 (1987) 86–88.
  * (21) J.-S. Wang and R. H. Swendsen, Cluster monte carlo algorithms, Physica A: Statistical Mechanics and its Applications 167 (1990), no. 3 565 – 579.
  * (22) P. W. Kasteleyn and C. M. Fortuin, Phase Transitions in Lattice Systems with Random Local Properties, Physical Society of Japan Journal Supplement, Vol. 26. Proceedings of the International Conference on Statistical Mechanics held 9-14 September, 1968 in Koyto., p.11 26 (1969) 11.
  * (23) C. M. Fortuin and P. W. Kasteleyn, On the Random cluster model. 1. Introduction and relation to other models, Physica 57 (1972) 536–564.
  * (24) J. Cardy, Scaling and Renormalization in Statistical Physics.  Cambridge Lecture Notes in Physics. Cambridge University Press, 1996.
  * (25) H. Gould, J. Tobochnik, and W. Christian, An introduction to computer simulation methods, vol. 1.
  * (26) G. Jaeger, The ehrenfest classification of phase transitions: introduction and evolution, Archive for history of exact sciences 53 (1998), no. 1 51–81.
  * (27) P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol, Extracting and composing robust features with denoising autoencoders, in Proceedings of the 25th international conference on Machine learning, pp. 1096–1103, ACM, 2008.
  * (28) P. Vincent, H. Larochelle, I. Lajoie, Y. Bengio, and P.-A. Manzagol, Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion, Journal of machine learning research 11 (2010), no. Dec 3371–3408.
  * (29) G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov, Improving neural networks by preventing co-adaptation of feature detectors, arXiv:1207.0580.
  * (30) F. Chollet et al., “Keras.” https://keras.io, 2015.
  * (31) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., Tensorflow: a system for large-scale machine learning., in OSDI, vol. 16, pp. 265–283, 2016.
  * (32) J. Hubbard, Electron Correlations in Narrow Energy Bands, Proceedings of the Royal Society of London Series A 276 (Nov., 1963) 238–257.
  * (33) J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, J. Phys. C6 (1973) 1181–1203. [,349(1973)].
  * (34) P. B. Arnold and G. D. Moore, Monte Carlo simulation of O(2) $\phi^{4}$ field theory in three-dimensions, Phys. Rev. E64 (2001) 066113, [cond-mat/0103227]. [Erratum: Phys. Rev.E68,049902(2003)].
  * (35) A. Athenodorou and M. Teper, On the spectrum and string tension of U(1) lattice gauge theory in 2 + 1 dimensions, JHEP 01 (2019) 063, [arXiv:1811.06280].
  * (36) J. Liddle and M. Teper, The Deconfining phase transition in D=2+1 SU(N) gauge theories, arXiv:0803.2128.


# Fast Multiplication of Matrices with Decay

Matt Challacombe and Nicolas Bock<EMAIL_ADDRESS>Theoretical
Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545
freeon.org

###### Abstract

A fast algorithm for the approximate multiplication of matrices with decay is
introduced; the Sparse Approximate Matrix Multiply (SpAMM) reduces complexity
in the product space, a different approach from current methods that economize
within the matrix space through truncation or rank reduction. Matrix
truncation (element dropping) is compared to SpAMM for quantum chemical
matrices with approximate exponential and algebraic decay. For matched errors
in the electronic total energy, SpAMM is found to require fewer to far fewer
floating point operations relative to dropping. The challenges and
opportunities afforded by this new approach are discussed, including the
potential for high performance implementations.

Fast Matrix Multiplication, Electronic Structure, Matrix Function, Matrix
Decay, Spectral Projection, Generalized N-Body Problem, Sparse Matrix-Matrix
Multiply, Linear Scaling Quantum Chemistry

††preprint: LA-UR 10-07458

##  I Introduction

For large dense linear algebra problems, the computational advantage offered
by fast matrix-matrix multiplication can be substantial, even with seemingly
small gains in asymptotic complexity. Relative to conventional multiplication
which is $\mathcal{O}\left(n^{3}\right)$, Strassen’s algorithm achieves
$\mathcal{O}\left(n^{2.8}\right)$, while the Coppersmith and Winograd method
is $\mathcal{O}\left(n^{2.38}\right)$. For these dense methods, balancing the
trade off between cost, complexity and error is an active area of research
Demmel and Higham (1992); Demmel et al. (2007); Yuster and Zwick (2005).

On the other hand, large sparse problems are typically handled with
conventional sparse matrix techniques, with only small concessions between
multiplication algorithms. Intermediate to these regimes, a wide class of
problems exist that involve matrices with decay111A matrix $A$ is said to
decay when its matrix elements decrease exponentially, as
$|a_{i,j}|<c\lambda^{|i-j|}$ , or algebraically as
$|a_{i,j}|<\frac{c}{|i-j|^{\lambda}+1}$ with indicial separation $|i-j|$. In
non-synthetic cases, the separation $|i-j|$ typically corresponds to an
underlying physical distance $|\vec{r}_{i}-\vec{r}_{j}|$, _e.g._ of basis
functions, finite elements, _etc_. See Figure 3 as well as the excellent work
by Benzi and co-authors on this topic in References Benzi and Razouk (2007);
Benzi et al. (2010). where sparsity exists only asymptotically under an
approximate linear algebra, historically involving matrix economization
through element dropping or rank reduction. Often, problems with decay occur
in the construction of matrix functions, notably the matrix inverse Benzi and
Tuma (2000), the matrix exponential Iserles (1999), and in the case of
electronic structure theory, the Heaviside step function Challacombe (1999,
2000); Benzi and Razouk (2007); Benzi et al. (2010). The use of an approximate
matrix algebra is also an active area of interest in the solution of large
eigenproblems Simoncini and Elden (2002); Simoncini and Szyld (2003, 2005);
Challacombe (2010).

Many approaches to a sparse approximate linear algebra exist for matrices with
decay Galli (1996); Goedecker (1999); Goedecker and Scuseria (2003); Li et al.
(2005), largely predicated upon the truncation of matrix elements, with the
recent work of Benzi providing the most detailed analysis so far Benzi and
Razouk (2007); Benzi et al. (2010). In this contribution, we develop sparse
matrix multiplication as a generalized $N$-body problem Gray and Moore (2001),
and introduce a fast algorithm based on hierarchical truncation in the three-
dimensional space $i,\,j,\,k\in[1,n]$ of the product
$C_{ij}=\sum_{k}A_{ik}B_{kj}$, where $A$ and $B$ decay exponentially or
algebraically fast enough222Algebraic decay sufficient to achieve a fast
$\mathcal{O}\left(n\,\lg\,n\right)$ or better complexity is an open question
similar to that of conditional convergence in the shape dependent summation of
dipole and quadrupole components of the Lorentz field. . Viewing the product
from a length scale perspective1, if matrix elements decay as
$\mathcal{O}(1/r^{\lambda})$, then the bulk of the product interactions will
decay as $\mathcal{O}(1/r^{2\lambda})$. For small $\lambda$, the difference
between truncation in the matrix space and the product space may be
significant.

##  II A Sparse Approximate Matrix Multiply

The quadtree matrix representation,

$A^{k}=\left(\begin{array}[]{cc}A_{11}^{k+1}&A_{12}^{k+1}\\\
A_{21}^{k+1}&A_{22}^{k+1}\end{array}\right),\>k=0,\,\ldots,\,k_{{\rm\textrm{max}}},$
(1)

is the basis for recursive matrix-matrix multiplication, $C^{k}=A^{k}\cdot
B^{k}.$ For conventional recursive multiplication, the operator “$\cdot$” is
just the row-column product, while in fast multiplication, it represents an
economized sequence of operations with reduced complexity and a more
complicated error accumulation. In Reference Bini and Lotti (1980), Bini and
Lotti carried out a detailed error analysis for recursive matrix
multiplication schemes, and derived component-wise bounds of the form

$|\tilde{c}_{ij}-c_{ij}|<a\,b\,\epsilon\,n\,\lg_{2}\,n$ (2)

where $\tilde{c_{ij}}$ is a matrix element computed to within precision
$\epsilon$, $c_{ij}$ is its exact counterpart, $a=\max_{ij}|a_{ij}|$ and
$b=\max_{ij}|b_{ij}|$. While providing a sharp bound, the max norm does not
immediately lend itself to the recursive separation of interaction magnitudes.
Consider instead the framework provided by an arbitrary sub-multiplicative
matrix norm $\|\cdot\|$:

$\displaystyle\|C^{k}\|$ $\displaystyle\leq$ $\displaystyle\|A^{k}\|\|B^{k}\|$
$\displaystyle\leq$
$\displaystyle\|A_{11}^{k+1}\|\|B_{11}^{k+1}\|+\|A_{12}^{k+1}\|\|B_{21}^{k+1}\|$
$\displaystyle+$
$\displaystyle\|A_{11}^{k+1}\|\|B_{12}^{k+1}\|+\|A_{12}^{k+1}\|\|A_{22}^{k+1}\|$
$\displaystyle+$
$\displaystyle\|A_{21}^{k+1}\|\|B_{11}^{k+1}\|+\|A_{22}^{k+1}\|\|B_{21}^{k+1}\|$
$\displaystyle+$
$\displaystyle\|A_{21}^{k+1}\|\|B_{12}^{k+1}\|+\|A_{22}^{k+1}\|\|A_{22}^{k+1}\|\,.$

This structure suggests an algorithm we call the Sparse Approximate Matrix
Multiply (SpAMM), which recursively tests each of the 8 contributions in
Equation (II) for significance in accordance with a given numerical threshold
$\tau$:

${\tt
SpAMM}(A^{k},B^{k})=\begin{cases}\qquad\qquad\qquad\qquad\qquad\qquad\qquad
0&{\tt if}\quad\|A^{k}\|\|B^{k}\|<\tau\\\\[8.5359pt]
\qquad\qquad\qquad\qquad\qquad\qquad\;\;A^{k}\cdot B^{k}&{\tt elseif}\quad
k=k_{\textrm{max}}\\\\[8.5359pt] \left(\begin{array}[]{rr}{\tt
SpAMM}(A_{11}^{k+1},B_{11}^{k+1})&{\tt SpAMM}(A_{11}^{k+1},B_{12}^{k+1})\\\
+{\tt SpAMM}(A_{12}^{k+1},B_{21}^{k+1})&+{\tt
SpAMM}(A_{12}^{k+1},B_{22}^{k+1})\\\\[8.5359pt] {\tt
SpAMM}(A_{21}^{k+1},B_{11}^{k+1})&{\tt SpAMM}(A_{21}^{k+1},B_{12}^{k+1})\\\
+{\tt SpAMM}(A_{22}^{k+1},B_{21}^{k+1})&+{\tt
SpAMM}(A_{22}^{k+1},B_{22}^{k+1})\end{array}\right)&{\tt else}\end{cases}$ (4)

The truncated product space accessed by SpAMM is shown in Figure 1 for
matrices with exponential and algebraic decay.

Figure 1: Hierarchical truncation of the product space ($i,\,j,\,k$) using
$\tau=10^{-8}$ for synthetic matrices of dimension $n=512$, with
$A_{ij}=\exp\left(-\left|i-j\right|\right)$ and
$B_{ij}=\exp\left(-2\left|i-j\right|\right)$ in (A),
$A_{ij}=B_{ij}=\begin{cases}1/|i-j|^{3}&i\neq j\\\ 0&\textrm{else}\end{cases}$
in (B). Each box above the finest ($k=k_{max}$) scale represents truncation.

At each tier in SpAMM, the local truncation error is bounded by
$\|\widetilde{C}^{k}-C^{k}\|<\tau$, with an error accumulation structurally
similar to rounding in conventional recursive multiplication, except that
truncation is not guaranteed to occur at each tier in the recursion, and
truncation does not retain even the approximate magnitude of the avoided sub-
product. To see the difference between truncation and round-off, consider the
generic, norm-wise bound for recursive multiplication employed by Demmel and
co-workers in Reference Demmel et al. (2007):

$\displaystyle\|\widetilde{C}-C\|$ $\displaystyle<\mu(n)$
$\displaystyle\epsilon\,\|A\|\|B\|\,+\mathcal{O}\left(\epsilon^{2}\right),$
(5)

with $\mu(n)$ a low order polynomial $\sim n^{d}$ and $d\geq 1$. Unlike
rounding which is an error of commission, SpAMM creates errors of omission. If
$\tau<\|A\|\|B\|$, then no work is performed and we obviously
have$\|\widetilde{C}-C\|<\tau$, which is different than
$\|\widetilde{C}-C\|<\epsilon\|A\|\|B\|$ corresponding to the case of
comparable round-off and truncation parameters $\epsilon\sim\tau$. One could
certainly bound SpAMM rigorously by decreasing $\tau$ with increasing depth,
$\tau^{k+1}=\tau^{k}/8$, but that would be overly pessimistic, not taking into
account signed error accumulation or the localization and attenuation of
errors due to decay.

SpAMM is similar to the $\mathcal{H}$-matrix algebra of Hackbusch and co-
workers, where off-diagonal sub-matrices are treated as reduced rank
factorizations (truncated SVD), typically structured and grouped to reflect
properties of the underlying operators Grasedyck and Hackbusch (2003). For
problems with rapid decay, truncated SVD may behave in a similar way to simple
dropping schemes. SpAMM is different than the $\mathcal{H}$-matrix algebra as
it achieves separation uniquely in the product space and does not rely on a
reduced complexity representation of matrices. For very slow decay, the
$\mathcal{H}$-matrix algebra may certainly offer a path for intractable
problems for which SpAMM is ineffective.

Figure 2: Space filling curves map atoms in Cartesian space (A) onto an index
that is locality preserving (B), leading to clustering of matrix elements with
respect to indicial separation (C) when interactions are short ranged. The
octree generated by SpAMM in the three-dimensional product space, shown in
Figure 1, is equivalent to a second space filling curve (spatial hashing)
Warren and Salmon (1992, 1995); Samet (2006) that orders both matrices (C) and
the product space (D), features that can be exploited to achieve domain
decomposition and load balance (colors in C and D).

An understanding of the error accumulation in SpAMM must account for the decay
properties of $A$ and $B$, which in non-synthetic cases are intimately related
to the effects of ordering and structure of the underlying physical, chemical
or engineering application. Also, ordering will determine the relative
efficiencies of both matrix truncation and SpAMM, particularly under blocking.
The ordering used here is based on the Space Filling Curve (SFC) method
described in Challacombe (2000) (Hilbert atom ordering), and shown in Figure
2, involving also a second tier of ordering in the product space. The first
ordering maps atoms that are close in Cartesian-space to entries close in the
index space of the matrix. The second ordering is a natural consequence of the
SpAMM multiply, which recursively maps out a multi-level octree (see also
Figure 1) with cuboid coordinates that are equivalent to a spatial-hash;
sorting the hash produces a three-dimensional SFC in the product space Warren
and Salmon (1992, 1995); Samet (2006). This two-tiered structure is novel,
providing an encompassing scheme that parlays the locality of physical
interactions into the data locality of matrices (element clustering), and into
the irregular product space. This approach is applicable to other problems
with underlying decay properties, in which the finite-elements, radial basis
functions, _etc._ replace atoms, or graph theoretical methods (nested
dissection, Cuthill-McKee, _etc_.) replace the first-tier SFC altogether. In
either case, the correspondence between the recursive SpAMM product space and
a three-dimensional SFC provides a tool to exploit both spatial and temporal
locality in the distribution of work and data.

##  III Results and Discussion

Figure 3: Decay of normed density matrix atom-blocks
$\left\|P_{ab}\right\|_{F}$ with Cartesian separation
$\vec{|R_{a}}-\vec{R}_{b}|$ for the largest molecules in each sequence: a 450
atom water cluster (red), a 752 atom (4,3) nanotube (orange) and a 780 atom
(3,3) nanotube (blue).

In this initial contribution, we limit ourselves to exploring the numerical
and computational behavior of SpAMM applied to problems in electronic
structure theory, and compare its relative merits to the dropping of matrix
elements. In ${\cal O}\left(n\right)$ electronic structure calculations
Goedecker (1999); Goedecker and Scuseria (2003); Li et al. (2005), a primary
source of error due to sparse matrix-multiplication develops from early steps
in the iterative construction of the Heaviside matrix function $P=\theta[F-\mu
I]$ (density matrix purification), which is a projector of the effective
Hamiltonian (Fockian) $F$ Niklasson et al. (2003). Starting from the basis
spanning $F$, purification drives eigenvalues to $0$s or $1$s, whereupon error
accumulation due to an approximate linear algebra is quenched Niklasson et al.
(2003). Under a given approximate linear algebra, the electronic energy
$E_{{\rm el}}={\rm Tr}(P.F)$ is a global measure of accumulated error. In the
following Section, we carry out purification on three molecular sequences of
increasing size, water clusters, (4,3) nanotubes and (3,3) nanotubes. In each
case, a Fockian $F$ was obtained from a fully converged Self-Consistent-Field
(SCF) cycle and used as basis for the purification; our numerical experiments
probe only errors within one density matrix solve (40-60 multiplies), and do
not address error propagation throughout the SCF cycle. As the rate of decay
slows towards SCF convergence, these calculations represent the instance of
minimum decay, shown for each of the largest molecular species in Figure 3.
Within each sequence, the same number of iterative steps (40-60 matrix
multiplications) are taken using the TC2 purification algorithm Niklasson
(2002) and either: (I) matrix element dropping and exact multiplication or
(II) SpAMM.

Figure 4: Average number of 4x4 matrix multiplies (MATMULs) for a sequence of
(4,3) nanotubes at the RHF/STO-2G level of theory with dropping (red) and
SpAMM (blue). Figure 5: Average number of 4x4 matrix multiplies (MATMULs) for
a sequence of (3,3) nanotubes at the LDA/STO-2G/ level of theory with dropping
(red) and SpAMM (blue). Figure 6: Average number of 4x4 matrix multiplies
(MATMULs) for a sequence of water clusters at the RHF/6-31G** level of theory
with dropping (red) and SpAMM (blue).

In this study, we chose $k_{{\rm max}}$ to yield 4x4 blocks at the finest
level of resolution, corresponding to the most aggressive use of single
precision SSE on the x86 architecture. In all cases the matrix norm employed
is the Frobenius norm $\left\|\cdot\right\|\equiv\left\|\cdot\right\|_{F}$.
Thresholds, $\tau$, have been adjusted to roughly match relative errors in the
electronic energy,$\Delta E_{\textrm{el}}=|\tilde{E}_{{\rm el}}-E_{{\rm
el}}|/|E_{{\rm el}}|$ between the two schemes, (I) and (II), and the average
number of 4x4 MATMULs per purification step are reported. Multiplications are
only a proxy for CPU time, as neither case accounts for symbolic overheads
associated with the multiply (CSR, recursive-tree _etc_.). In the case of
dropping, SpAMM was also used in the matrix multiply but with zero threshold.
After each multiplication in the dropping scheme, a filter was applied to the
resultant, dropping blocks at the 4x4 level of resolution using the criteria
$\left\|P^{k}\right\|_{F}<\tau$. Results for the three molecular sequences are
shown in Figures 6, 4 and 5.

For comparable values of $\Delta E_{\textrm{el}}$, SpAMM was found to employ
from slightly fewer multiplies for the (4,3) nanotube, to dramatically less in
the case of the water clusters. This result is somewhat surprising, since the
spatial decay is slower in the case of the (4,3) nanotube than for the water
clusters, as shown in Figure 3. While it seems reasonable to attribute this
unexpected result to the effects of dimensionality, further study is required
to be sure. Its also worth noting that the advantage of SpAMM relative to
dropping is brought down with decreasing $\tau$; for $\tau=0$ they both revert
to the same $\mathcal{O}(n^{3})$ complexity. For both systems with exponential
decay, error appears tightly controlled albeit within an as yet unknown bound.
Note however, that unlike matrix truncation which leads to a product error
that is $\mathcal{O}\left(\tau^{2}\right)$, SpAMM leads to a truncation error
that is $\mathcal{O}\left(\tau\right)$.

Comparing the quasi one-dimensional metallic (3,3) nanotube with slow
algebraic decay to the insulating (4,3) nanotube with exponential decay, SpAMM
gains substantially over dropping in the case of slower decay. However, the
ability of SpAMM to achieve a linear scaling complexity in the case of the
metallic system remains in question, as in the case of the tightest threshold,
the SpAMM errors do not appear to be well controlled, and the cost does not
appear linear, at least not in this size regime. On the other hand, the SpAMM
result does enjoy a significant reduction in cost relative to dropping, and
the error increase is modest.

##  IV Conclusion

The Sparse Approximate Matrix Multiply (SpAMM) is a fast method for matrices
with decay, which is different from element dropping or the
$\mathcal{H}$-algebra in that it uniquely involves truncation of the product
space rather than the matrix space. For matrices with exponential or fast
algebraic decay, SpAMM can achieve stable error control comparable to element
dropping, but with a greatly reduced number of floating point operations. The
results presented here are preliminary, and have not yet explored the
interesting problems of slow decay in the asymptotic limit, ordering, error
bounds, the cost of recursion, high performance implementations or broader
gauges of accuracy and efficiency, such as the use of SpAMM in the context of
the Self-Consistent-Field cycle. Of particular concern is the relationship
between complexity, matrix decay and error control. Based on our numerical
tests, we postulate that the algorithm is at worst $\mathcal{O}(n\lg n)$ for
matrices with sufficiently fast decay.

In addition to similarities with the $\mathcal{H}$-algebra, SpAMM falls under
the rubric of the generalized $N$-body problem Gray and Moore (2001). From
this perspective, its worth noting also the connection between matrix-matrix
multiplication as $N$-body problem and matrix-matrix multiplication as spatial
join Amossen and Pagh (2009), as well as between $N$-body problems and data
base theory in general Samet (2006).

Next, we draw attention to the second tier of Space Filling Curve (SFC) shown
in Figure 2, which provides a mechanism for domain decomposition and load
balance that is proven for parallel irregular problems Warren and Salmon
(1992, 1995); Aluru and Sevilgen (1997); Devine et al. (2005). Also, the
improvement gained in Reference Buluc and Gilbert (2008) on going from a one-
dimensional to a two-dimensional matrix partitioning scheme for the parallel
SpMM suggests that partitioning the three-dimensional product space instead
may provide an even higher degree of flexibility and granularity.

The authors acknowledge support through Los Alamos LDRD award ER20110230
(computational co-design) as well as funds from the U.S. Department of Energy.
Los Alamos National Laboratory is operated by the Los Alamos National
Security, LLC, for the National Nuclear Security Administration of the U.S.
Department of Energy under Contract No. DE-AC52-06NA25396. Special
acknowledgments go to the International Ten Bar Café for tasty libations in a
scientific and collegial atmosphere, and to Michele Benzi for valuable input.

## References

  * Demmel and Higham (1992) J. W. Demmel and N. J. Higham, ACM Trans. Math. Softw. 18, 274 (1992).
  * Demmel et al. (2007) J. Demmel, I. Dumitriu, O. Holtz, and R. Kleinberg, Numer. Math. 106, 199 (2007).
  * Yuster and Zwick (2005) R. Yuster and U. Zwick, ACM Transactions on Algorithms 1, 2 (2005).
  * Benzi and Razouk (2007) M. Benzi and N. Razouk, Electronic Transactions on Numerical Analysis 28, 16 (2007).
  * Benzi et al. (2010) M. Benzi, P. Boito, and N. Razouk, Manuscript in preparation (2010).
  * Benzi and Tuma (2000) M. Benzi and M. Tuma, SIAM Journal on Scientific Computing 21, 1851 (2000).
  * Iserles (1999) A. Iserles, _How large is the exponential of a banded matrix?_ (1999).
  * Challacombe (1999) M. Challacombe, J. Chem. Phys. 110, 2332 (1999).
  * Challacombe (2000) M. Challacombe, Compututer Physics Communications 128, 93 (2000).
  * Simoncini and Elden (2002) V. Simoncini and L. Elden, Bit Numerical Mathematics 42, 159 (2002).
  * Simoncini and Szyld (2003) V. Simoncini and D. B. Szyld, SIAM Journal on Scientific Computing 25, 454 (2003).
  * Simoncini and Szyld (2005) V. Simoncini and D. B. Szyld, SIAM Review 47, 247 (2005).
  * Challacombe (2010) M. Challacombe, arXiv quant-ph, 1001.2586 (2010).
  * Galli (1996) G. Galli, Current Opinion in Solid State & Materials Science 1, 864 (1996).
  * Goedecker (1999) S. Goedecker, Reviews of Modern Physics 71, 1085 (1999).
  * Goedecker and Scuseria (2003) S. Goedecker and G. Scuseria, Computing in Science Engineering 5, 14 (2003).
  * Li et al. (2005) Z. Y. Li, W. He, and J. L. Yang, Progress in Chemistry 17, 192 (2005).
  * Gray and Moore (2001) A. G. Gray and A. W. Moore, in _Advances in Neural Information Processing Systems_ (MIT Press, 2001), vol. 4, pp. 521–527.
  * Bini and Lotti (1980) D. Bini and G. Lotti, Numerische Mathematik 36, 63 (1980).
  * Grasedyck and Hackbusch (2003) L. Grasedyck and W. Hackbusch, Computing 70, 2003 (2003).
  * Warren and Salmon (1992) M. S. Warren and J. K. Salmon, in _Supercomputing ’92_ (IEEE Comp. Soc., Los Alamitos, 1992), pp. 570–576, (1992 Gordon Bell Prize winner).
  * Warren and Salmon (1995) M. S. Warren and J. K. Salmon, _A parallel, portable and versatile treecode_ (SIAM, Philadelphia, 1995), chap. 1.
  * Samet (2006) H. Samet, _Foundations of Multidimensional and Metric Data Structures_ (Morgan Kaufmann, 2006).
  * Niklasson et al. (2003) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, J. Comp. Phys. 118, 8611 (2003).
  * Niklasson (2002) A. M. N. Niklasson, Physical Review B 66, 5 (2002).
  * Amossen and Pagh (2009) R. Amossen and R. Pagh, in _Proceedings of the 12th International Conference on Database Theory_ (ACM, 2009), pp. 121–126.
  * Aluru and Sevilgen (1997) S. Aluru and F. E. Sevilgen, in _Proceedings of the 4th IEEE Conference on High Performance Computing_ (1997), pp. 230–235.
  * Devine et al. (2005) K. D. Devine, E. G. Boman, R. T. Heaphy, B. A. Hendrickson, J. D. Teresco, J. Faik, J. E. Flaherty, and L. G. Gervasio, Applied Numerical Mathematics 52, 133 (2005), ADAPT ’03: Conference on Adaptive Methods for Partial Differential Equations and Large-Scale Computation.
  * Buluc and Gilbert (2008) A. Buluc and J. R. Gilbert, in _ICPP ’08: Proceedings of the 2008 37th International Conference on Parallel Processing_ (IEEE Computer Society, Washington, DC, USA, 2008), pp. 503–510.