<EMAIL_ADDRESS>[Pavel Kroupa]P. Kroupa [Rachel Parziale]R. Parziale
[Moritz Haslbauer]M. Haslbauer

# The Phantom of RAMSES user guide for galaxy simulations using Milgromian and
Newtonian gravity

S. T. Nagesh  Argelander-Institut für Astronomie, Universität Bonn, Auf dem
Hügel 71, 53121 Bonn, Germany.  I. Banik  Helmholtz-Institut für Strahlen- und
Kernphysik, Universität Bonn, Nussallee 14-16, 53115 Bonn, Germany.  I. Thies
missing  Helmholtz-Institut für Strahlen- und Kernphysik, Universität Bonn,
Nussallee 14-16, 53115 Bonn, Germany and Astronomical Institute, Faculty of
Mathematics and Physics; Charles University in Prague, V Holešovičkách 2,
CZ-180 00 Praha, Czech Republic.  B. Famaey  Université de Strasbourg, CNRS
UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l’Université,
67000 Strasbourg, France  N. Wittenburg  missing  missing  Max-Planck-Institut
für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany and Helmholtz-
Institut für Strahlen- und Kernphysik, Universität Bonn, Nussallee 14-16,
53115 Bonn, Germany.

###### Abstract

This document describes the general process of setting up, running, and
analysing disc galaxy simulations using the freely available program phantom
of ramses (por). This implements Milgromian Dynamics (MOND) with a patch to
the ramses grid-based $N$-body and hydrodynamical code that uses adaptive mesh
refinement. We discuss the procedure of setting up isolated and interacting
disc galaxy initial conditions for por, running the simulations, and analysing
the results. This manual also concisely documents all previously developed
MOND simulation codes and the results obtained with them.

##  1 Introduction

Milgromian Dynamics (MOND) is an extension of Newtonian dynamics to encompass
the observed dynamics in the Solar System as well as in galaxies without
postulating invisible haloes around them [1]. MOND computes the gravitational
potential of galaxies using only the distribution of baryons. It has been very
successful in this regard, especially because it predicted some very tight
scaling relations which were subsequently observed [2, 3]. These are a
consequence of Milgrom’s formula

$\displaystyle g\leavevmode\nobreak\ =\leavevmode\nobreak\
\sqrt{g_{\mathrm{N}}a_{{}_{0}}}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \textrm{for}\leavevmode\nobreak\ \leavevmode\nobreak\
g_{\mathrm{N}}\ll a_{{}_{0}}=1.2\times 10^{-10}\,\textrm{m s}^{-2},$ (1)

where $a_{{}_{0}}$ is Milgrom’s constant, $g$ is the strength of the true
gravity, and $g_{\mathrm{N}}$ is that of the Newtonian gravity. To achieve a
generalization of gravity applicable in non-spherical systems, MOND requires a
generalized Poisson equation derived from a Lagrangian. Two classical variants
have been proposed, one with an aquadratic Lagrangian [AQUAL; 4], and one with
a Lagrangian making use of an auxiliary field, which is called the quasi-
linear formulation of MOND [QUMOND; 5]. MOND may be a consequence of the
quantum vacuum [6, 7, 8, 9]. Reviews of MOND can be found in [2, 10].

phantom of ramses [por; 11] is a numerical implementation of QUMOND, whose
field equation for the potential $\Phi$ is

$\displaystyle\nabla^{2}\Phi\leavevmode\nobreak\ \equiv\leavevmode\nobreak\
-\nabla\cdot\bm{g}\leavevmode\nobreak\ =\leavevmode\nobreak\
-\nabla\cdot\left(\nu\bm{g}_{\mathrm{N}}\right)\,,$ (2)

where $\nu$ is the MOND interpolating function with argument $y\equiv
g_{\mathrm{N}}/a_{{}_{0}}$, with $\bm{g}$ and $\bm{g}_{\mathrm{N}}$ being the
true and Newtonian gravitational acceleration vectors, respectively, and
$v\equiv\left|\bm{v}\right|$ for any vector $\bm{v}$. The current version of
por uses the simple form of the interpolating function (e.g. equation 5 of
[12])

$\displaystyle\nu\left(y\right)\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{1}{y}}\,.$ (3)

$\bm{g}_{\mathrm{N}}$ is found from the baryonic density $\rho_{b}$ using the
standard Poisson equation

$\displaystyle\nabla\cdot\bm{g}_{\mathrm{N}}\leavevmode\nobreak\
=\leavevmode\nobreak\ -4\mathrm{\pi}G\rho_{b}\,.$ (4)

The boundary condition for the MOND potential far from an isolated matter
distribution is

$\displaystyle\Phi\leavevmode\nobreak\ =\leavevmode\nobreak\
\sqrt{GMa_{{}_{0}}}\ln R\,,$ (5)

where $M$ is the mass within the simulation box, and $R$ is the distance from
its barycentre in the simulation unit of length.

A handful of Milgromian $N$-body codes were developed before por to handle
MOND computations, and these have been applied to various problems. The first
multi-grid, Milgromian $N$-body code was developed by [13] to investigate the
stability of disc galaxies. This was later extended to simulate how they might
warp due to the external field effect [EFE; 14]. Another $N$-body solver which
implemented the AQUAL formulation of MOND was developed and used to study the
evolution of spiral galaxies using pure stellar discs [15]. Gas dynamics was
later included using a sticky particle scheme at low resolution [16]. n-mody
was developed to solve the Milgromian Poisson equation in spherical
coordinates [17] and used to investigate dynamical friction [18], orbit
instabilities [19], and stellar kinematics [20, 21, 22]. Milgromian $N$-body
codes tailored to cosmological simulations have also been developed [23, 24,
25, 26]. Another $N$-body solver called raymond [27] was developed to
implement both the AQUAL [4] and QUMOND [5] formulations of MOND. raymond has
been applied to cosmological [28], galaxy cluster [29], and other problems.
However, not all the aforementioned $N$-body codes can be applied to generic
scenarios simultaneously involving particles, gas dynamics, and star
formation.

The fortran-based por code was developed by Fabian Lüghausen [11]. It is a
customized version of ramses [30], which exclusively uses Newtonian dynamics
to compute gravity. por can compute gravity using MOND by numerically solving
Equation 2. Since por is a patch to ramses, it inherits use of the adaptive
mesh refinement (AMR) technique. por is equipped to handle particles, gas
dynamics, and star formation, and can be applied to diverse problems. It
allows the user to compute gravity in both Milgromian and Newtonian frameworks
[11].

This document serves as a tutorial/manual for the general use of por, with
some suggestions for the specific case of setting up and simulating a disc
galaxy. Most of the steps and parameters described here are specific to por,
except the installation of ramses. For a detailed description of individual
parameters, it is always recommended to read the ramses
manual111https://bitbucket.org/rteyssie/ramses/src/master/. Most of the
parameters and files described here can be edited safely without disturbing
the core algorithms. Before changing parameters or files that are not
mentioned here, it is important to fully understand the workings and
consequences of the change.

In Section 2, we explain the installation procedure of ramses and por. In
Section 3, we explain how to set up MOND disc templates using Disk Initial
Conditions Environment (dice), and thereby generate rotating disc initial
conditions for por. In Section 4, we describe the workings of por, focusing on
particle-only and hydrodynamical runs with and without star formation. In
Section 5, random turbulence generation is briefly discussed. Section 6
discusses the extract_por tool used to analyse particle data in ramses
simulation outputs. In Section 7, we mention all publications based on por. We
conclude in Section 8.

##  2 Installation and setup of the code

The por patch by Fabian Lüghausen is rated to work with the 2015 version of
ramses, which has since been modified (the latest ramses version is available
here 1). Later versions are not compatible with por. It is therefore
recommended to use the 2015 version of ramses with por. The jointly tested
version of ramses and por is available here
222https://bitbucket.org/SrikanthTN/bonnpor/src/master/, in the PoR_hydro
folder.

The following steps describe the installation and compilation of ramses and
por. These procedures are adapted from the ramses manual, where they are
described further.

  1. 1.

The main folder needed for compilation of ramses is bin. In the bin folder,
there is a makefile. Now, do:

⬇

$ cd ~/PoR_hydro/ramses/bin/

  2. 2.

In the makefile, certain flags need to be changed. Makefile:

Compilation time parameters NVECTOR = 32 NDIM = 3 NPRE = 8 NVAR = 6 NENER = 0
SOLVER = hydro #PATCH = ../patch/phantom_units #PATCH =
../patch/phantom_staticparts (particle-only run) #PATCH =
../patch/hydro/phantom_merger PATCH = ../patch/hydro/phantom_extfield (hydro
run) EXEC = RAMSES

  3. 3.

All these flags are explained in the ramses manual1. F90 and FFLAGS should be
set carefully.

F90 = mpif90 -frecord-marker=4 -O3 -ffree-line-length-none -g -fbacktrace
FFLAGS = -x f95-cpp-input $(DEFINES)$

  4. 4.

F90 sets the Fortran compiler and FFLAGS is used to specify the required MPI
libraries, which are mainly used for parallel computing. This is important
given the likely high computational cost. The default makefile2 uses the
above-mentioned F90 and FFLAGS. If one’s computer is not compatible with these
default parameters, they can be changed in the makefile.

  5. 5.

Once all the required flags are set, compile the code:

⬇

$ make

After compilation, one can test the installation as described in section 2.3
of the ramses manual1.

  6. 6.

To make the files again, go to the bin folder and execute:

⬇

$ make clean

$ make

The ramses manual was written in 2002 and has not been updated since, so there
might be subsequent modifications to the parameter file. One must use the
phantom patch to do simulations in MOND. For Newtonian simulations, it is
recommended to use this patch and set the mond flag to .false. in the
namelist.

###  2.1 Compilation of the code with the por patch

We now describe the procedure to link the por patch and re-compile ramses. To
activate the por patch, the following parameters must be specified in the
makefile (# means a comment):

NDIM = 3 PATCH = ../patch/phantom_staticparts #PATCH = ../patch/phantom #PATCH
= ../patch/hydro/phantom_extfield

By default, ramses uses periodic boundary conditions, but the por patch in2
specifies and uses different boundary conditions appropriate to isolated
galaxy simulations. The 2015 version of ramses available in2 contains a
staticpart patch and a hydrodynamical patch with compulsory additional merger
and EFE patches whose effects can be disabled (Section 4.2). These patches are
task-specific customizations of por. For a particle run, phantom_staticparts
can be used while for a hydrodynamical run, phantom_extfield can be used. Only
one patch can be used at a time. One must change the path to the user’s
directory before making the file. After specifying the parameters, make the
file again.

##  3 Disk Initial Conditions Environment (DICE)

Any galaxy or cosmological simulation needs initial conditions. The ramses
user guide refers to two websites for these, but they are for cosmological
runs. The music333https://bitbucket.org/ohahn/music/src/master/ code also
provides initial conditions for the latter, and is recommended. The setup of
cosmological MOND simulations will be discussed elsewhere. This guide focuses
on galaxy simulations, for which we generate initial conditions with an
adapted version of Disk Initial Conditions Environment [dice; 31]. The
original version of dice can be found
here444https://bitbucket.org/vperret/dice/src/master/. It is not compatible
with MOND or the 2015 version of ramses2, so we used a modified version of
dice available here2. This has two versions, one for particle-only runs (the
dice_particle folder, hereafter p-dice) and the other for hydrodynamical runs
(in dice_gas, hereafter h-dice). These algorithms were developed by Graeme
Candlish, Roy Truelove, Indranil Banik, and Ingo Thies. Both are equipped to
initialize disc galaxies in MOND, but in principle other methods could be used
and advanced with the por patch. Before installing dice, CMake, GSL, and FFTW
must be installed. If this is not already the case, installation instructions
are provided here2.

###  3.1 Installation and setup

As mentioned above, the folders dice_particle and dice_gas contain p-dice and
h-dice, respectively. Extract them to /local in one’s home directory. In
dice_particle, the disc folder is required for disc galaxy simulations. Now,
in disc, the bin folder contains the makefile needed for compilation, while
the example folder contains the parameter files. To compile p-dice, execute:

⬇

$ cd dice_particle

$ cd disc

$ mkdir build

$ cd build

$ cmake ..

$ make

$ make install

To compile h-dice, execute:

⬇

$ cd dice_gas

$ mkdir build

$ cd build

$ cmake ..

$ make

$ make install

h-dice does not contain an additional disc folder like p-dice.

###  3.2 Running DICE

The dice_gas and disc (in dice_particle) folders contain four sub-folders:

  1. 1.

cmake should not be altered,

  2. 2.

build will contain the executable,

  3. 3.

src contains the source files which encode the physics required for
computation, and

  4. 4.

example contains files required to generate the initial conditions, with task-
specific configuration files like M31, M33, and generic scenarios like a disc
galaxy, disc with a bulge etc. Only the Milky Way (MW), M31, and M33 cases are
rated to work.

We used the test_mw.config configuration file for our disc galaxy:

Redshift 3.0 Galaxy ../../example/params_files/testMilkyWay.params Filename
dice_highz ICformat Gadget2 Nthreads 32

In the .config file, specify the path to the parameter file. The redshift is
unused in our MONDified dice. The testMilkyWay.params is the parameter file
used, though other params files exist in the /example/params_files folder.
Custom templates can be created using these parameter files, though only the
MW, M31, and M33 cases are rated to work. There are mainly three types of
parameters in testMilkyWay.params: global parameters, outer disc, and inner
disc. For both p-dice and h-dice, once the parameter file is set, go one
directory up and execute:

⬇

$ cd bin

$ ./dice ../../example/test_mw.config

After execution, 2 output files named Milky_Way_output_p2_k0.txt and
Milky_Way_rotation_curve.txt will be created in the bin folder. The rotation
curve is only required for hydrodynamical simulations.

##  4 Running por

After compilation of ramses with the required por patch, one can customize the
namelist file available in the PoR_namelist folder to meet a scientific goal.
PoR_namelist consists of all the namelist files we have used for our runs.
PoR.nml is a general template which can be customized. PoR-static.nml is the
file we used for our particle-only run, while Test_hydro_mw_NSFR.nml was used
for the hydrodynamical run without star formation. There is a general namelist
folder which consists of .nml files that can be used to test e.g. the
installation of ramses.

###  4.1 Particle-only run (staticpart patch)

Use the /patch/phantom_staticparts patch and PoR-static.nml, take care on the
boundary conditions:

&RUN_PARAMS poisson=.true. pic=.true. mond=.true. – Activates MOND poisson
solver nrestart=0 – used to restart a run from any output / &AMR_PARAMS .
levelmin=7 levelmax=12 ngridmax=2000000 boxlen=1024.0 npartmax=2000000 /
ngridmax and npartmax should be of order $10^{6}$ to avoid memory errors.
&OUTPUT_PARAMS foutput=8000 – Frequency of outputs in terms of coarse steps
noutput=100 – Number of outputs to be generated delta_tout=100. – Interval of
the output in Myr tend=10000. – Simulation end time, in this case 10 Gyr /
&INIT_PARAMS filetype=‘ascii’ initfile= ../path/to/Milky_Way_output_p2_k0.txt
as the input. / &POISSON_PARAMS a0_ms2=1.2e-10 m_threshold=1.e+30 – critical
part, set it based on usage gravity_type=0 cg_levelmin=999 / &BOUNDARY_PARAMS
nboundary=6 ibound_min=-1, 1, 0, 0, 0, 0, ibound_max=-1, 1, 0, 0, 0, 0,
jbound_min= 0, 0,-1, 1, 0, 0, jbound_max= 0, 0,-1, 1, 0, 0, kbound_min= 0, 0,
0, 0,-1, 1, kbound_max= 0, 0, 0, 0,-1, 1, bound_type= 1, 1, 1, 1, 1, 1, /

There are some parameters which are mandatory in all runs, such as
&Run_Params, &AMR_Params, &Output_Params, &Init_Params. Others vary based on
specific requirements. Most of these parameters are detailed in the ramses
manual1, so we only stress those specific to por. Parameters not shown here
should not be changed unless required.

The staticpart patch integrates particles below a certain mass m_threshold,
while more massive particles are kept static but are considered when
evaluating $\bm{g}_{\mathrm{N}}$ in Equation 4. This method is an effective
way to save computation time. If one wants to evolve all the stellar particles
in a particle-only simulation, then m_threshold should be set to a suitably
large value, e.g. $10^{30}M_{\odot}$. The units from the dice output are the
same as required by por for input (i.e. $M_{\odot}$, kpc, and km/s), while
units.f90 has the units used by por in which $G=1$.

One can modify &Output_Params and &AMR_Params, but it is not recommend to
tamper with other blocks. In the above example, the Milky_Way_output_p2_k0.txt
obtained from p-dice is given as the input file, with the rotation curve
unused. Once all parameters are set in the namelist file, the simulation can
be started by executing:

⬇

$ mpiexec -n 32 ../ramses3d ../filename.nml

This calls the simulation to run on 32 CPUs using parallel computing (the
number can be changed). Regardless of the directory of execution, one must
specify full paths to the ramses3d and namelist files. To run these
simulations without parallel computing, simply execute:

⬇

$ ../bin/ramses3d ../filename.nml

Users should check the computing capacity before running simulations without
parallel computing.

After starting the simulation, it might terminate with error message \-
“SEGSEV - invalid memory reference”. This is a memory error, which can be
solved by increasing npartmax up to $10^{7}$ and ngridmax up to $8\times
10^{6}$ (the codes are not rated for larger values). The npartmax variable
must be at least equal to the number of particles in the dice template.
Turning off the movie may help. CPU and memory errors are a bit alarming, but
are easily overcome and should not be a big concern for beginners.

The simulation will produce output folders for each snapshot. During the run,
if the memory allocated is too small, the simulation will stop and ask to
increase the number of grid cells. One must then go back to the namelist file
and increase ngridmax or npartmax based on what is asked. The restart protocol
is rated to work, so restart the run from the last output file by setting
nrestart to the desired output number (the default of 0 means to start from
scratch). If the run stops before finalising output_..45, set nrestart = 45
and resume the run by executing the above-mentioned command.

###  4.2 Hydrodynamical run without star formation

We performed this run with the /patch/hydro/phantom_extfield patch, which is a
modification to por. The EFE and merger scenarios are included using the MOND
Poisson solver, but both features can be turned off.

####  4.2.1 DICE with gas component

dice is again used to set up the initial conditions. Since the gas component
is included, we used h-dice in the dice_gas folder. To include a gas
component, the test_MilkyWay.params was slightly modified:

################## # Global parameters ################### # Virial velocity
of the galaxy [km/s] v200 200.0 # Virial mass of the galaxy [1e10 Msol] #
Overrides the v200 parameter m200 9.15#8.4 old Gas_fraction 0.2 Gas_T 50000.0

The highlighted lines are the new additions to the h-dice template. These
lines specify the gas component parameters. The gas fraction depends on the
galaxy, here 20% gas fraction was used for the MW. The gas temperature should
be set equal to another parameter called T2_ISM, which is present in the
namelist file of por. One must be careful that the gas fraction should be
greater than the mass fraction of the outer component, in this particular
case, more than 18%. This is because the distribution of gas in h-dice is done
in a particular way, so one should be cautious while setting the gas fraction
[32]. The template has a default mass fraction of 17.64% for the outer disc
component, with the remaining 82.36% for the inner disc [33]. This version is
only rated for two exponential disc components. For a beginner, it is
recommend to take advice at this point before proceeding further.

h-dice can be started the same way as p-dice. After starting the dice run, one
might notice a message in the terminal “Gas is too cold to satisfy the Toomre
condition everywhere. Increase $T$ by a factor of …”, and/or “WARNING: only
writing data for component 1”. The first message is just a warning about the
global disc stability, and has no impact on the results $-$ it can be ignored.
Temperature here is used as a measure of velocity dispersion including
turbulence, so it is not the true gas temperature. The second message can also
be ignored $-$ it indicates that the second component defined in dice is
treated as stars and used for calculating the potential, but not printed in
the output as the gas will be added in por [32]. The particle data written to
the disc template file contains only the stellar component.

The rotation curve file has columns for the gas disc scale height and its
radial gradient. This is critical as the gas component is created in por
itself (in the merger and extfield patch) by reading in the gas data from the
rotation curve file and some parameters to be set in the namelist file, e.g.
gas mass and temperature. Care is needed to ensure compatibility of the
parameters used for dice and por.

####  4.2.2 por with merger and external field patch

In the hydrodynamical case, we did not explicitly set up an isolated disc
galaxy, but instead adapted the merger template
condinit555http://www.physics.usyd.edu.au/~tepper/codes/ramses/trunk/doc/html/patch_2hydro_2merger_2condinit_8f90_source.html.
This sets up two disc galaxies in the simulation box, so we switched off the
second galaxy by setting its mass and velocity to zero and placing it outside
the simulation box. The namelist file for each run should be customized as
required, we show part of Test_hydro_mw_NSFR.nml as an example:

&RUN_PARAMS mond=.true. Activate_g_ext=.false. &INIT_PARAMS filetype=‘ascii’
initfile(1)=‘path/to/the/DICE/output/’ &MERGER_PARAMS rad_profile=‘double_exp’
z_profile=‘sech_sq’ Mgas_disc1=45.75 Mgas_disc2=0 IG_density_factor=1.0e-2
T2_ISM=40.d3 scale_a2=1. Vcirc_dat_file1=‘Milky_Way_rotation_curve.txt’
Vcirc_dat_file2=‘Milky_Way_rotation_curve.txt’
ic_part_file_gal1=‘Milky_Way_output_p2_k0.txt’
ic_part_file_gal2=‘Milky_Way_output_p2_k0.txt’ gal_center1= 0.,0.,0.
gal_center2= 2000,0.,0. Vgal1=0.,0.,0. Vgal2=0.,0.,0.

The namelist has other parameters, but we show only those critical to the
simulation. If one were to use a similar setup, the following suggestions are
helpful:

  1. 1.

In the ramses makefile, one must provide the path to the external field patch
and recompile ramses.

  2. 2.

The por patch can accommodate both MOND and Newtonian physics. The latter is
used if one sets mond = .false., allowing simulations with both gravity
theories using por.

  3. 3.

Setting Activate_g_ext to false turns the EFE off. Hydrodynamical simulations
with the EFE are discussed further in [32].

  4. 4.

For the initfile(1), one must give the path to the directory where
Milky_Way_output_p2_k0.txt and Milky_Way_rotation_curve.txt are present. These
files should be specified for ic_part_file_gal1 and Vcirc_dat_file1,
respectively. The same path and files can be given for the second galaxy,
which is unused here.

  5. 5.

The main things that need attention are the &Merger_Params. The gas mass of
the galaxy is in units of $10^{9}M_{\odot}$. If simulating interacting
galaxies, they should not start too close together. We switched the second
galaxy off by setting its gas mass and velocity to 0 and placing it outside
the box, e.g. box size = 500 kpc, gal_center2 = (2000, 0, 0) kpc. The first
galaxy was placed at the box centre.

  6. 6.

For isolated simulations, both galaxies should have zero initial velocity, i.e
Vgal1 and Vgal2 should be zero. gal_axis defines the disc’s spin axis. For
standard isolated simulations, use $\left(0,0,1\right)$ for counter-clockwise
rotation around the $z$-axis, or $\left(0,0,-1\right)$ for clockwise.

  7. 7.

The T2_ISM parameter in the namelist and Gas_T in the h-dice template should
be equal. The temperature floor T2_star should be set to a slightly lower
value than T2_ISM. We used T2_ISM = 40,000 K and T2_Star = 30,000 K. The
simulation is not rated to work with T2_ISM or T2_Star below 25,000 K.

After taking care of all these parameters, one can start the run, leading to
creation of the output folders in due course. The simulations are RAM and
memory intensive $-$ a hydrodynamical disc galaxy advanced for 1.5 Gyr on a
4-core laptop could take a week or two depending on the RAM and might occupy
up to 100 GB of hard disk space. These estimates would vary depending on the
parameters and machine used $-$ see section 4.2 of the ramses manual1 for more
details.

###  4.3 Hydrodynamical run with star formation

Converting gas into stars requires careful treatment of baryons, for which
ramses is well equipped and tested [34]. Since por is just a modification to
the Poisson solver, it does not affect the baryon treatment, inheriting that
of standard ramses.

To activate star formation, one has to include the &Physics_Params in the
namelist file. One can add &Physics_Params to the Test_hydro_mw_NSFR.nml and
activate star formation. Alternatively, one can use the MW_hydro_SFR.nml
provided in2, which we used for our star formation run.

&PHYSICS_PARAMS cooling=.true. g_star=1.6666D0 n_star=0.1D0 eps_star=0.0D0
t_star=3.0d0 (star forming timescale in Gyr) T2_star=4.0d4 /

All the above parameters are described in the ramses manual. The t_star
parameter is the star formation timescale in Gyr. Setting it to a finite, non-
zero value activates star formation. One can add other parameters as per
requirements.

##  5 Random turbulence generation

To allow for initial turbulence and (optionally) density fluctuations, a
random perturbation algorithm has been included based on the square-square
subdivision [35]. It is similar to the well-known diamond-square algorithm
widely used for the generation of random terrains, but provides a higher
quality of randomness and fewer artefacts. The algorithm first applies random
perturbations on a $2\times 2\times 2$ cubic array. At each subsequent step,
the cube cells are subdivided into $2\times 2\times 2$ arrays and perturbed
again, while the magnitude of the perturbation is reduced $2\times$ (unless
the user chooses a different value). Additional factors can be applied to the
magnitude of each step, following a user-defined power spectrum. The resulting
random noise is then multiplied with the density and/or the three velocity
components to get turbulence.

The algorithm requires some additional variables to be set in the
&Merger_Params in the namelist file, and an extra parameter file qqm3d.par.
The extra lines in the namelist are:

&MERGER_PARAMS … flg_qqm3d=-1 !Master switch (-1 means off) devflat_dens=1.0
!density mean unperturbed level devscal_dens=0.1 !density deviation scale
devflat_vel=1.0 !same, for velocities devscal_vel=0.1 scale_objsize=1.0 !size
of the perturbation mask

The master switch controls the overall usage of the random perturbation
algorithm. “-1” means “off”, other modes are:

  * •

0 or 10: only density is perturbed,

  * •

1: only adds absolute perturbation to velocities,

  * •

2: combines modes 0 and 1,

  * •

11 and 12: like 1 and 2, but with velocity perturbations scaled by the
circular velocity (recommended),

  * •

21 and 22: like 1 and 2, but with velocity perturbation relative to actual
velocities (experimental).

The parameter file qqm3d.par contains:

**** Setup parameters for qqm4ramses **** 8 2.5 4. 1 nsize,fsize,scalh
ini,balance values 1 1 1. init mode, deviate (0:lin, 1:Gauss), power 0.0
initial corner master values 10 0 1.0 1.0 hr_mode, stop rnd h after n iter
($<=0$:off), hreduce iter factor+power 309562 -1 seed, seed initialization
mode —- Corner perturbation scaling —- 0.2 scalh00 —- Feature power spectrum
—- 0.1 scalh01 … —– Corner initial values —– 0\. x01 …

Only the lines most relevant for beginners are shown. Other lines are mostly
experimental and should be left as they are, unless the user looks at the
source code for more details about their purpose. The most relevant values for
hr_mode are:

  * •

4: uses the lines from the power spectrum block as weightings. The magnitude
of the first non-zero perturbation is equal to devscal and will be reduced by
hreduce (typically $1/2$) for each refinement level.

  * •

10: uses a flat power spectrum with starting level fsize. Non-integer values
are used via an interpolation scheme.

The other hr_modes should not be used for scientific runs. For details, see
the source file qqm4ramses.f90 and its subroutine init_qqm3d.

##  6 Extraction of data with extract_por

For the extraction of particle data, a tool called extract_por was developed
by Ingo Thies and used here. Now including additional features related to star
formation, extract_por_sfr is available here2. This is a user-friendly tool
that does not require much time to learn. After the tool is downloaded, it can
be extracted to /home/local. Installation is done by executing:

⬇

$ make xpordata

Inside extract_por, fmtRAMSES.par is the parameter file where the extraction
parameters can be set. It contains:

‘path/to/your/outputfiles’ 38 Output No. 32 Number of CPU threads 0 COM reset
—- RADIAL BINNING SECTION —- 10\. 500 binning radius (in system units), nbins
—- Image params —- 1 1 0 flgimage 250 250 imagecenter 200 200 image width 500
500 nbx,nby 1.5 1.5 hx,hy smoothing (pixel units)

Again, not all the parameters are detailed here, the file itself being very
well commented. Only the parameters that might be important for a beginner are
shown.

  1. 1.

The path should only specify where the output folders are located, not the
output folders themselves. Thus, ../../output_0001 will not be recognised.

  2. 2.

COM reset subtracts the center of mass position and velocity. It could be used
if the object of interest lies outside the field of view.

To just extract the particle positions and velocities, only parameters until
the Partial COM section are important. After setting the parameters, execute:

⬇

$./xpordata

Based on the number of output files selected, the corresponding number of
part.asc and sfr.dat files will be created. The part.asc files contain data in
ascii format with the following column meanings:

⬇

1-3: position, 4-6: velocity, 7: mass,

8: particle ID

The sfr.dat file contains:

⬇

1: time interval in Myr, 2: SFR in M_Sun/Myr

The extraction algorithm calculates the total stellar mass in a given
snapshot, and evaluates the difference in stellar mass between two snapshots.
One can resolve the SFR better by changing delta_tout in the namelist file, or
by extracting particle birth times.

Any tool can be used to extract and plot the results from the part.asc files.
Even extract_por can be used for plotting, in which case all the sections
below Image_Params can be helpful. These sections can be used to set the
projected density, resolution etc. To use extract_por for plotting, the
following suggestions might be helpful:

  1. 1.

In Image_params, unless one has a special case like [36], using 2:rgb or
3:rgbw is not helpful. Set it to 1:gray. This works and one can set the
required projected density.

  2. 2.

The use of binning radius might be critical for resolution. In case of poor
resolution, increase the number of bins. To increase the pixel resolution/zoom
in, reduce the image field of view in fmtRAMSES.par, i.e. reduce the bin
sizes.

  3. 3.

Users could set the box width equal to the simulation box size and locate the
galaxy manually.

  4. 4.

hx and hy smoothing smoothens the image. This can be changed based on needs.

  5. 5.

All parameters below hx, hy smoothing are not to be changed.

Run extract_por and expect two output files to be produced:

  1. 1.

part.asc

  2. 2.

image.dat

The image.dat has the data required for plotting the image (particle
positions). The simplest way is to use gnuplot:

⬇

$ gnuplot –> plot “image.dat” with image

##  7 Tests and publications using por

Since its development in 2015, por has been applied to a variety of problems.
A first implementation showed that the observed dynamics in polar ring
galaxies is explained naturally in MOND [37]. [38] compared Antennae-like
galaxy encounters in MOND and in dark matter models, studying the evolution
towards merging and the triggering of star formation in both models. The
Galactic tidal streams of Sagittarius [39] and Palomar 5 [40] were
investigated as gravitational experiments, with the latter’s asymmetry
interpreted as evidence for the EFE. [41] showed that the satellite galaxy
planes of the MW and M31 might arise from a past encounter between them. [42]
showed that exponential disc galaxies form naturally in MOND out of collapsing
post-Big Bang gas clouds. [32] simulated M33, finding that its long-term
evolution is well understood in MOND, especially its weak bar and lack of a
bulge. Their work also details some of the numerical methods, especially in
h-dice and the extfield patch.

##  8 Conclusions

por [11] is a general-purpose $N$-body and hydrodynamical solver for MOND. It
is based on adapting ramses, whose modern version is not compatible with the
por patch. It is recommended to use por from here2. This manual is a generic
outline with which one can understand the basics required to set up, run, and
analyse por simulations. The above-mentioned files like the namelist and
patches like staticpart and hydro are custom-made for a specific purpose, so
care should be taken before using them for a different application. All the
algorithms and tools mentioned in this guide are available here2, and are
rated to work.

## Acknowledgements

IB is supported by an Alexander von Humboldt Foundation postdoctoral research
fellowship. BF acknowledges funding from the Agence Nationale de la Recherche
(ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No. 834148). On a historical note, when
PK joined the HISKP at Bonn University in 2013, financial means became
available that allowed Fabian Lüghausen to be hired as a PhD student co-
supervised by PK and BF, to program the por patch for ramses, and to buy
computer servers for MOND simulations. This led to the development of the por
code [11]. The authors would like to thank Jan Pflamm-Altenburg and the
referees for comments which helped to clarify this guide.

## References

  * Milgrom [1983] M. Milgrom.  A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis.  _ApJ_ , 270:365–370, July 1983.  10.1086/161130.
  * Famaey and McGaugh [2012] B. Famaey and S. S. McGaugh.  Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions.  _Living Reviews in Relativity_ , 15:10, September 2012.  10.12942/lrr-2012-10.
  * Lelli et al. [2017] F. Lelli, S. S. McGaugh, J. M. Schombert, and M. S. Pawlowski.  One Law to Rule Them All: The Radial Acceleration Relation of Galaxies.  _ApJ_ , 836:152, February 2017.  10.3847/1538-4357/836/2/152.
  * Bekenstein and Milgrom [1984] J. Bekenstein and M. Milgrom.  Does the missing mass problem signal the breakdown of Newtonian gravity?  _ApJ_ , 286:7–14, November 1984.  10.1086/162570.
  * Milgrom [2010] M. Milgrom.  Quasi-linear formulation of MOND.  _MNRAS_ , 403:886–895, April 2010.  10.1111/j.1365-2966.2009.16184.x.
  * Milgrom [1999] M. Milgrom.  The modified dynamics as a vacuum effect.  _Phys. Lett. A_ , 253:273–279, March 1999.  10.1016/S0375-9601(99)00077-8.
  * Pazy [2013] E. Pazy.  Quantum statistical modified entropic gravity as a theoretical basis for MOND.  _Phys. Rev. D_ , 87(8):084063, April 2013.  10.1103/PhysRevD.87.084063.
  * Verlinde [2016] E. P. Verlinde.  Emergent Gravity and the Dark Universe.  _SciPost Physics_ , 2:16, November 2016.  10.21468/SciPostPhys.2.3.016.
  * Smolin [2017] L. Smolin.  MOND as a regime of quantum gravity.  _Physical Review D_ , 96(8):083523, October 2017\.  10.1103/PhysRevD.96.083523.
  * Milgrom [2015] M. Milgrom.  MOND theory.  _Canadian Journal of Physics_ , 93(2):107–118, February 2015.  10.1139/cjp-2014-0211.
  * Lüghausen et al. [2015] F. Lüghausen, B. Famaey, and P. Kroupa.  Phantom of RAMSES (POR): A new Milgromian dynamics N-body code.  _Canadian Journal of Physics_ , 93:232–241, February 2015\.  10.1139/cjp-2014-0168.
  * Banik and Zhao [2018a] I. Banik and H. Zhao.  Testing gravity with wide binary stars like $\alpha$ Centauri.  _MNRAS_ , 480:2660–2688, October 2018a.  10.1093/mnras/sty2007.
  * Brada and Milgrom [1999] R. Brada and M. Milgrom.  Stability of Disk Galaxies in the Modified Dynamics.  _ApJ_ , 519:590–598, July 1999.  10.1086/307402.
  * Brada and Milgrom [2000] R. Brada and M. Milgrom.  The Modified Dynamics is Conducive to Galactic Warp Formation.  _ApJL_ , 531(1):L21–L24, March 2000.  10.1086/312510.
  * Tiret and Combes [2007] O. Tiret and F. Combes.  Evolution of spiral galaxies in modified gravity.  _A &A_, 464:517–528, March 2007.  10.1051/0004-6361:20066446.
  * Tiret and Combes [2008] O. Tiret and F. Combes.  Evolution of spiral galaxies in modified gravity. II. Gas dynamics.  _A &A_, 483:719–726, June 2008.  10.1051/0004-6361:200809357.
  * Londrillo and Nipoti [2009] P. Londrillo and C. Nipoti.  N-MODY: a code for collisionless N-body simulations in modified Newtonian dynamics.  _Memorie della Societa Astronomica Italiana Supplementi_ , 13:89, January 2009.
  * Nipoti et al. [2008] C. Nipoti, L. Ciotti, J. Binney, and P. Londrillo.  Dynamical friction in modified Newtonian dynamics.  _MNRAS_ , 386(4):2194–2198, June 2008.  10.1111/j.1365-2966.2008.13192.x.
  * Nipoti et al. [2011] C. Nipoti, L. Ciotti, and P. Londrillo.  Radial-orbit instability in modified Newtonian dynamics.  _MNRAS_ , 414(4):3298–3306, July 2011.  10.1111/j.1365-2966.2011.18632.x.
  * Wu and Kroupa [2013] X. Wu and P. Kroupa.  The dynamical phase transitions of stellar systems and the corresponding kinematics.  _MNRAS_ , 435(1):728–742, October 2013.  10.1093/mnras/stt1332.
  * Wu and Kroupa [2018] X Wu and P Kroupa.  Gas Expulsion in MOND: The Possible Origin of Diffuse Globular Clusters and Ultra-faint Dwarf Galaxies.  _ApJ_ , 853(1):60, January 2018.  10.3847/1538-4357/aaa081.
  * Wu and Kroupa [2019] X. Wu and P. Kroupa.  The kinematics of star clusters undergoing gas expulsion in Newtonian and Milgromian dynamics.  _MNRAS_ , 487(3):4012–4024, August 2019.  10.1093/mnras/stz1519.
  * Llinares et al. [2008] C. Llinares, A. Knebe, and H. Zhao.  Cosmological structure formation under MOND: a new numerical solver for Poisson’s equation.  _MNRAS_ , 391:1778–1790, December 2008.  10.1111/j.1365-2966.2008.13961.x.
  * Llinares [2011] C Llinares.  _On the linear and non-linear evolution of dust density perturbations with MOND_.  PhD thesis, Rijksuniversiteit Groningen, January 2011.
  * Angus et al. [2011] G. W. Angus, A. Diaferio, and P. Kroupa.  Using dwarf satellite proper motions to determine their origin.  _MNRAS_ , 416:1401–1409, September 2011.  10.1111/j.1365-2966.2011.19138.x.
  * Angus et al. [2013] G. W. Angus, A. Diaferio, B. Famaey, and K. J. van der Heyden.  Cosmological simulations in MOND: the cluster scale halo mass function with light sterile neutrinos.  _MNRAS_ , 436(1):202–211, November 2013.  10.1093/mnras/stt1564.
  * Candlish et al. [2015] G. N. Candlish, R. Smith, and M. Fellhauer.  RAyMOND: an N-body and hydrodynamics code for MOND.  _MNRAS_ , 446:1060–1070, January 2015.  10.1093/mnras/stu2158.
  * Candlish [2016] G. N. Candlish.  The velocity field in MOND cosmology.  _MNRAS_ , 460:2571–2585, August 2016.  10.1093/mnras/stw1130.
  * Candlish et al. [2018] G. N. Candlish, R. Smith, Y. Jaffé, and A. Cortesi.  Consequences of the external field effect for MOND disc galaxies in galaxy clusters.  _MNRAS_ , 480:5362–5379, November 2018.  10.1093/mnras/sty2228.
  * Teyssier [2002] R. Teyssier.  Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES.  _A &A_, 385:337–364, April 2002.  10.1051/0004-6361:20011817.
  * Perret et al. [2014] V. Perret, F. Renaud, B. Epinat, P. Amram, F. Bournaud, T. Contini, R. Teyssier, and J.-C. Lambert.  Evolution of the mass, size, and star formation rate in high redshift merging galaxies. MIRAGE - A new sample of simulations with detailed stellar feedback.  _A &A_, 562:A1, February 2014.  10.1051/0004-6361/201322395.
  * Banik et al. [2020] I. Banik, I. Thies, G. Candlish, B. Famaey, R. Ibata, and P. Kroupa.  The global stability of M33 in MOND.  _ApJ_ , 905(2):135, December 2020.  10.3847/1538-4357/abc623.
  * Banik and Zhao [2018b] I. Banik and H. Zhao.  The escape velocity curve of the Milky Way in modified Newtonian dynamics.  _MNRAS_ , 473:419–430, January 2018b.  10.1093/mnras/stx2350.
  * Rasera and Teyssier [2006] Y. Rasera and R. Teyssier.  The history of the baryon budget: Cosmic logistics in a hierarchical universe.  _A &A_, 445(1):1–27, January 2006.  10.1051/0004-6361:20053116.
  * Miller [1986] G. S. P. Miller.  The definition and rendering of terrain maps.  _SIGGRAPH Comput. Graph._ , 20(4):39–48, August 1986.  ISSN 0097-8930.  10.1145/15886.15890.
  * Oehm et al. [2017] W. Oehm, I. Thies, and P. Kroupa.  Constraints on the dynamical evolution of the galaxy group M81.  _MNRAS_ , 467:273–289, May 2017.  10.1093/mnras/stw3381.
  * Lüghausen et al. [2013] F. Lüghausen, B. Famaey, P. Kroupa, G. Angus, F. Combes, G. Gentile, O. Tiret, and H. Zhao.  Polar ring galaxies as tests of gravity.  _MNRAS_ , 432(4):2846–2853, July 2013.  10.1093/mnras/stt639.
  * Renaud et al. [2016] F. Renaud, B. Famaey, and P. Kroupa.  Star formation triggered by galaxy interactions in modified gravity.  _MNRAS_ , 463:3637–3652, December 2016.  10.1093/mnras/stw2331.
  * Thomas et al. [2017] G. F. Thomas, B. Famaey, R. Ibata, F. Lüghausen, and Pavel Kroupa.  Stellar streams as gravitational experiments. I. The case of Sagittarius.  _A &A_, 603:A65, July 2017.  10.1051/0004-6361/201730531.
  * Thomas et al. [2018] G. F. Thomas, B. Famaey, R. Ibata, F. Renaud, N. F. Martin, and P. Kroupa.  Stellar streams as gravitational experiments. II. Asymmetric tails of globular cluster streams.  _A &A_, 609:A44, January 2018.  10.1051/0004-6361/201731609.
  * Bílek et al. [2018] M. Bílek, I. Thies, P. Kroupa, and B. Famaey.  MOND simulation suggests the origin of some peculiarities in the Local Group.  _A &A_, 614:A59, June 2018.  10.1051/0004-6361/201731939.
  * Wittenburg et al. [2020] N. Wittenburg, P. Kroupa, and B. Famaey.  The formation of exponential disk galaxies in MOND.  _ApJ_ , 890, February 2020.  10.3847/1538-4357/ab6d73.


# On the Size and Width of the Decoder of a Boolean Threshold Autoencoder

Tatsuya Akutsu  Bioinformatics Center, Institute for Chemical Research, Kyoto
University  Avraham A. Melkman  Department of Computer Science, Ben-Gurion
University of the Negev

###### Abstract

In this paper, we study the size and width of autoencoders consisting of
Boolean threshold functions, where an autoencoder is a layered neural network
whose structure can be viewed as consisting of an _encoder_ , which compresses
an input vector to a lower dimensional vector, and a _decoder_ which
transforms the low-dimensional vector back to the original input vector
exactly (or approximately). We focus on the decoder part, and show that
$\Omega(\sqrt{Dn/d})$ and $O(\sqrt{Dn})$ nodes are required to transform $n$
vectors in $d$-dimensional binary space to $D$-dimensional binary space. We
also show that the width can be reduced if we allow small errors, where the
error is defined as the average of the Hamming distance between each vector
input to the encoder part and the resulting vector output by the decoder.

Keywords: Neural networks, Boolean functions, threshold functions,
autoencoders.

##  1 Introduction

Extensive studies have been done on artificial neural networks not only in
artificial intelligence and machine learning but also in theoretical computer
science [2, 12, 10]. Among various models of neural networks, _autoencoders_
have recently attracted much attention due to their ability to generate new
data, and have been applied to various areas including image processing [5,
11], natural language processing [11], and drug discovery [6]. An autoencoder
is a layered neural network consisting of two parts, an _encoder_ and a
_decoder_ , where the former transforms an input vector to a low-dimensional
vector and the latter transforms the low-dimensional vector to an output
vector which should be the same as or similar to the input vector [1, 4, 7,
3]. Therefore, an autoencoder maps input data to a low-dimensional
representation space. Such a mapping is obtained via unsupervised learning
that minimizes the difference between input and output data by adjusting
weights (and some other parameters).

Although autoencoders perform dimensionality reduction, a kind of data
compression, how data are compressed via autoencoders is still unclear. In
particular, the quantitative relationship between the compressive power and
the numbers of nodes and layers in autoencoders had been unclear. In order to
study the relationship, Melkman et al. analyzed the relations between the
architecture (e.g., the numbers of nodes and layers) of networks and their
compression ratios [8], using a layered _Boolean threshold network_ , which is
a discrete model of neural networks. A Boolean threshold network is equivalent
to a _threshold circuit_ [2, 12, 10] in which each node takes on values that
are either 1 (active) or 0 (inactive) and the activation rule for each node is
given by a Boolean threshold function. In [8] several architectures of
autoencoders were presented that map $n$ $D$-dimensional binary input vectors
into $d$-dimensional binary space and then recover the original input vectors.
In particular, they showed the following architectures for $d=\lceil\log
n\rceil$ (or, almost equivalently, $d=2\lceil\log\sqrt{n}\rceil$):

  * •

a 4-layers encoder with $O(\sqrt{n}+D)$ nodes;

  * •

a 5-layers autoencoder with $D/n/d/n/D$ architecture, where each parameter
means the number of nodes in the corresponding layer;

  * •

a 7-layers autoencoder with $O(D\sqrt{n})$ nodes;

  * •

a decoder with depth $n+1$ and width $O(D)$, where the width is the maximum
number of nodes per layer (except for the input and output layers) and the
depth is the number of layers minus one.

However, it was unclear whether or not these results are optimal (or near
optimal).

In this paper, we focus on the decoder part because the decoders (in the
above-quoted results) use $O(D\sqrt{n})$, $O(n+D)$, or more nodes whereas the
encoder uses only $O(\sqrt{n}+D)$ nodes. We show an $\Omega(\sqrt{Dn/d})$
lower bound on the size of the perfect decoder, where a decoder is called
_perfect_ if it can always recover the original input vectors exactly. To our
knowledge, this is the first lower bound result on the size (i.e., the number
of nodes) of the autoencoder. As a positive result, we show that there exists
a perfect decoder with width $\max(\lceil n/B\rceil+B,BD)$ and a constant
depth, where $B$ is an arbitrary integer larger than 1. By letting
$B=\lceil\sqrt{n/D}\rceil$, we obtain an $O(\sqrt{Dn})$ upper bound on the
size and width of the decoder, which is relatively close to the lower bound
and improves the previous $O(D\sqrt{n})$ upper bound. In order to construct
decoders that have even smaller width we permit them to output vectors that
are in error, where the error is defined as the average of the Hamming
distance between each vector input to the encoder part and the resulting
vector output by the decoder. We show that there exists a decoder with width
$\max(\lceil n/B\rceil+B,(B-1)D+1)$ and a constant depth whose error is at
most $D({\frac{1}{B2^{B}}}+{\frac{1}{n}})$.

##  2 Preliminaries

A function $f$: $\\{0,1\\}^{h}\rightarrow\\{0,1\\}$ is called a _Boolean
threshold function_ if it is represented as

$\displaystyle f({\bf x})=\left\\{\begin{array}[]{ll}1,&{\bf a}\cdot{\bf
x}\geq\theta,\\\ 0,&\mbox{otherwise,}\end{array}\right.$

for some $({\bf a},\theta)$, where ${\bf a}$ is an $h$-dimensional integer
vector and $\theta$ is an integer. We will also denote the same function as
$[{\bf a}\cdot{\bf x}\geq\theta]$. An acyclic network is called a _Boolean
threshold network_ (BTN) if all activation functions in the network are
Boolean threshold functions.

In this paper, we only consider layered BTNs in which nodes are divided into
$L$-layers and each node in the $i$-th layer has inputs only from nodes in the
($i-1$)th layer ($i=1,\ldots,L-1$). Then, the states of nodes in the $i$-th
layer can be represented as a $W_{i}$-dimensional binary vector where $W_{i}$
is the number of nodes in the $i$-th layer and is called the _width_ of the
layer. A layered BTN is represented as ${\bf y}={\bf f}^{(L-1)}({\bf
f}^{(L-2)}(\cdots{\bf f}^{(1)}({\bf x})\cdots))$, where ${\bf x}$ and ${\bf
y}$ are the input and output vectors, respectively. and ${\bf f}^{(i)}$ is a
list of activation functions for the ($i+1$)th layer. The $0$-th and
($L-1$)-th layers are called the _input and output layers_ , respectively, and
the corresponding nodes are called _input and output nodes_ , respectively.
The _size_ , _depth_ , _width_ of a layered BTN are defined as the number of
nodes, the number of layers minus one, and the maximum number of nodes in a
layer (excepting the input and output layers), respectively. When we consider
autoencoders, one layer ($k$-th layer where $k\in\\{1,\ldots,L-2\\}$) is
specified as the _middle layer_ , and the nodes in this layer are called the
_middle nodes_ (see also Fig. 1). Then, the _middle vector_ ${\bf z}$,
_encoder_ ${\bf f}$, and _decoder_ ${\bf g}$ are defined by

$\displaystyle{\bf z}$ $\displaystyle=$ $\displaystyle{\bf f}^{(k-1)}({\bf
f}^{(k-2)}(\cdots{\bf f}^{(1)}({\bf x})\cdots))={\bf f}({\bf x}),$
$\displaystyle{\bf y}$ $\displaystyle=$ $\displaystyle{\bf f}^{(L-1)}({\bf
f}^{(L-2)}(\cdots{\bf f}^{(k)}({\bf z})\cdots))={\bf g}({\bf z}).$

Since the input vectors should be recovered from the middle vectors at the
output layer, we assume that ${\bf x}$ and ${\bf y}$ are $D$-dimensional
binary vectors and ${\bf z}$ is a $d$-dimensional binary vector with $d\leq
D$. A list of functions is also referred to as a _mapping_. The $i$-th element
of a vector ${\bf x}$ will be denoted by $x_{i}$, which is also used to denote
the node corresponding to this element. Similarly, for each mapping ${\bf f}$,
$f_{i}$ denotes the $i$-th function.

Figure 1: Architecture of a five-layer autoencoder, where ${\bf z}={\bf
f}({\bf x})={\bf f}_{2}({\bf f}_{1}({\bf x}))$ and ${\bf y}={\bf g}({\bf
z})={\bf f}_{4}({\bf f}_{3}({\bf z}))$.

Let $X_{n}=\\{{\bf x}^{0},\ldots,{\bf x}^{n-1}\\}$ be a set of $n$
$D$-dimensional binary input vectors that are all different. We define perfect
encoder, decoder, autoencoder as follows [8].

######  Definition 1.

A mapping ${\bf f}$: $\\{0,1\\}^{D}\rightarrow\\{0,1\\}^{d}$ is called a
_perfect encoder_ for $X_{n}$ if ${\bf f}({\bf x}^{i})\neq$ ${\bf f}({\bf
x}^{j})$ holds for all $i\neq j$.

######  Definition 2.

A pair of mappings $({\bf f},{\bf g})$ with ${\bf f}$:
$\\{0,1\\}^{D}\rightarrow\\{0,1\\}^{d}$ and ${\bf g}$:
$\\{0,1\\}^{d}\rightarrow\\{0,1\\}^{D}$ is called a _perfect autoencoder_ if
${\bf g}({\bf f}({\bf x}^{i}))={\bf x}^{i}$ holds for all ${\bf x}^{i}\in
X_{n}$. Furthermore, such ${\bf g}$ is called a _perfect decoder_.

Note that a perfect decoder exists only if there exists a perfect autoencoder.
Furthermore, it is easily seen from the definitions that if $({\bf f},{\bf
g})$ is a perfect autoencoder, ${\bf f}$ is a perfect encoder.

######  Example 3.

Let $X_{4}=\\{{\bf x}^{0},{\bf x}^{1},{\bf x}^{2},{\bf x}^{3}\\}$ where ${\bf
x}^{0}=(0,0,0)$, ${\bf x}^{1}=(1,0,0)$, ${\bf x}^{2}=(1,0,1)$, ${\bf
x}^{3}=(1,1,1)$. Let $D=3$ and $d=2$. Define ${\bf z}={\bf f}({\bf x})$ and
${\bf y}={\bf g}({\bf z})$ by

$\displaystyle f_{0}({\bf x})$ $\displaystyle=$
$\displaystyle[x_{0}+x_{1}-x_{2}\geq 1],~{}~{}~{}f_{1}({\bf
x})~{}=~{}[x_{2}\geq 1],$ $\displaystyle g_{0}({\bf z})$ $\displaystyle=$
$\displaystyle[z_{0}+z_{1}\geq 1],~{}~{}~{}g_{1}({\bf
z})~{}=~{}[z_{0}+z_{1}\geq 2],~{}~{}~{}g_{2}({\bf z})~{}=~{}[z_{1}\geq 1].$

This pair of mappings has the following truth table, which shows it to be a
perfect autoencoder.

$x_{0}$ | $x_{1}$ | $x_{2}$ | $z_{0}$ | $z_{1}$ | $y_{0}$ | $y_{1}$ | $y_{2}$
---|---|---|---|---|---|---|---
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
1 | 0 | 0 | 1 | 0 | 1 | 0 | 0
1 | 0 | 1 | 0 | 1 | 1 | 0 | 1
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1

##  3 A Lower Bound on the Size of the Decoder

In this section we derive a lower bound on the number of threshold units
(i.e., the number of nodes except those in the input layer) in the decoder BTN
of a Boolean threshold autoencoder which perfectly encodes every possible set
of $n$ $D$-dimensional vectors, and which has a middle layer whose size is
less than one third the size of the output layer. Note that $n\leq 2^{d}$,
since $n$ vectors are perfectly encoded.

Denote the number of threshold units of the decoder-BTN by $N$, with the
$i$-th unit having $d_{i}$ inputs, where $N$ includes the number of output
nodes but does not include the number of middle nodes.

######  Theorem 4.

There exists a set $X_{n}$ for which there is no perfect autoencoder with a
middle layer of size $d$ and a decoder of size less than
$\sqrt{\frac{D-1}{3d}}\sqrt{n}$.

###### Proof.

The idea is to calculate an upper bound on the total number of different sets
of $n$ vectors that can be generated at the output layer of the autoencoder
from all possible sets of $d$-dimensional vectors, input to the middle layer,
by all possible different decoder BTNs. Let’s denote that number $PD$ (the in-
Principle Decoded sets). Since the autoencoder is perfect, $PD$ is not smaller
than the number of different sets of $D$-dimensional output (=input) vectors
of the autoencoder, each of size $n$, which is $\binom{2^{D}}{n}$.

Next we derive an upper bound on $PD$. Each threshold unit computes a Boolean
function of the $d$ inputs (in the middle-layer). In the course of proving
Theorem 7.4 of [2], which provides a lower bound on the size of the universal
networks, it is shown that the number of different functions representable by
a threshold unit with in-degree $d_{i}$ in an acyclic BTN which has $d$ input
nodes is at most $2^{dd_{i}+1}$. Therefore, an upper bound on the total number
of functions from $\\{0,1\\}^{d}$ to $\\{0,1\\}^{D}$ computable by a decoder
BTN with $N$ nodes is

$\prod_{i=1}^{N}2^{dd_{i}+1}=2^{N+d\sum_{i=0}^{N-1}d_{i}}\leq 2^{N+dN^{2}}.$

Since the number of all possible ordered sets of $n$ $d$-dimensional vectors
in the middle layer does not exceed $2^{d}!/(2^{d}-n)!\leq(2^{d})^{n}$ we
conclude that

$\displaystyle{\frac{(2^{D}-n)^{n}}{n^{n}}}\leq\binom{2^{D}}{n}\leq
PD\leq\left[(2^{d})^{n}\right]\cdot\left[2^{N+dN^{2}}\right]$

Taking base 2 logarithms of both sides shows that the parameters of the
decoder satisfy $nd+N+dN^{2}\geq n\log(2^{D}-n)-n\log n$. Noting that $n\leq
2^{d}<2^{D-1}$ yields $N+dN^{2}>n(D-1)-2nd$. Applying the precondition
$d\leq{\frac{D-1}{3}}$ results in the inequality

$\displaystyle N+dN^{2}>{\frac{n}{3}}(D-1),$

from which the theorem follows. ∎

Remarks:

  * •

This lower bound is meaningful only if $\sqrt{\frac{D-1}{3d}}\sqrt{n}\geq D$
because the decoder contains $D$ output nodes.

  * •

In the proof $d_{i}$ was replaced by its upper bound $N$. Therefore, the lower
bound stated in the Theorem holds for any acyclic decoder.

  * •

The stated lower bound becomes an $\Omega(\sqrt{Dn/d})$ lower bound on the
widths of layered BTN decoders of constant depth.

##  4 Decoder with Width $\max(\lceil n/B\rceil+B,BD)$

In this section we present the architecture of a decoder BTN that is able to
decode perfectly any set of $n$ $D$-dimensional vectors that has been encoded
into a layer of size $d$ (so that the encoder and decoder together perfectly
autoencode every such set), where $d=\lceil\log n\rceil$.

We will use $\beta_{0},\ldots,\beta_{d-1}$ to refer to the middle nodes as
well as to their values, with $\beta$ denoting the vector
$(\beta_{0},\ldots,\beta_{d-1})$. $Int(\beta)$ denotes the integer encoded by
$\beta$, i.e. $Int(\beta)=\sum_{i=0}^{d-1}\beta_{i}2^{d-i-1}$.

Let $B$ be an integer greater than 1. For ease of exposition we assume that
$n$ can be divided by $B$ and use $n_{B}$ to denote $n/B$. Otherwise we can
let $n_{B}=\lceil n/B\rceil$. In this section we detail the architecture of a
$d/{\frac{n}{B}}+B/BD/D$ decoder (with $B=2^{K}$ for some integer $K$), see
Figure 2. In addition to $\beta$ it has the following nodes:

  * •

$\gamma_{0},\gamma_{1},\ldots,\gamma_{n_{B}-1},\gamma_{n_{B}},\gamma_{n_{B}+1},\ldots,\gamma_{n_{B}+B-1}$,

  * •

$y_{j,b}$ for $j=0,\ldots,D-1$ and $b=0,\ldots,B-1$,

  * •

$y_{j}$ for $j=0,\ldots,D-1$, the output nodes.

The activation function of the $\gamma_{i}$ are:

$\displaystyle\gamma_{i}$ $\displaystyle=$
$\displaystyle[\lfloor{\frac{Int(\beta)}{B}}\rfloor=i],~{}~{}~{}\mbox{for
$i=0,\ldots,n_{B}-1$},$ $\displaystyle\gamma_{n_{B}+h}$ $\displaystyle=$
$\displaystyle[Int(\beta)=h~{}(\mbox{mod $B$})],~{}~{}~{}\mbox{for
$h=0,\ldots,B-1$},$

where $[Int(\beta)=k]$ takes value 1 if $Int(\beta)=k$, and value 0 otherwise.
Note that division and modulo operations can be done by threshold circuits of
width $O(poly(\log n))$ and depth 5 [9] (recall that $\beta$ is a vector of
dimension $O(\log n)$), a negligible increase to the $O(\sqrt{Dn})$ width of
our network. Furthermore, if $B=2^{K}$ for some integer $K$, we do not need
such circuits because we can simply use the first $d-K$ bits of $\beta$ to
represent $Int(\beta)/B$ and the remaining $K$ bits to represent
$Int(\beta)\mbox{ mod $B$}$. Note also that $[Int(\beta)=k]$ can be calculated
by using a single unit with an activation function ${\bf
a}\cdot\beta\geq\theta$ such that

$\displaystyle a_{i}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1,&\mbox{if $k_{i}=1$},\\\
-1,&\mbox{otherwise},\end{array}\right.$ $\displaystyle\theta$
$\displaystyle=$ $\displaystyle\sum_{i}k_{i},$

where $k_{0}k_{1}\cdots k_{d-1}$ is the binary representation of an integer
$k$.

The activation function of $y_{j,b}$ is

$y_{j,b}=[\sum_{i=0}^{n_{B}+B-1}w_{i,j,b}\gamma_{i}\geq 2],$

where $w_{i,j,b}$ is given by

$\displaystyle w_{h,j,b}$ $\displaystyle=$ $\displaystyle
x_{j}^{Bh+b}~{}~{}~{}\mbox{for $h=0,\ldots,n_{B}-1$ and for
$b=0,\ldots,B-1$},$ $\displaystyle w_{n_{B}+h,j,b}$ $\displaystyle=$
$\displaystyle[h=b]~{}~{}~{}\mbox{for $h=0,\ldots,B-1$ and for
$b=0,\ldots,B-1$}.$

Numbering the $n$ input/output vectors so that each ${\bf x}^{k}$ is mapped to
$\beta$ with $Int(\beta)=k$, it is not difficult to verify that when ${\bf
x}^{k}$ is input to the autoencoder, with $k=B\ell+r,\ 0\leq r\leq B-1$, then
$y_{j,b}=[x^{B\ell+b}_{j}+[b=r]\geq 2]$. Hence $y_{j,b}=1$ if and only if
$b=r$ and $x^{k}_{j}=1$.

The activation function of $y_{j}$ is

$y_{j}=[\sum_{b=0}^{B-1}y_{j,b}\geq 1],$

which ensures that on input ${\bf x}^{k}$ the value of $y_{j}$ is $x^{k}_{j}$,
i.e. the decoder is perfect.

Although we assumed that ${\bf x}^{k}$ is mapped to $\beta$ with
$Int(\beta)=k$, this assumption does not pose any restriction on an encoder
because arbitrary permutations of input vectors can be considered. Therefore,
we have:

######  Theorem 5.

For any perfect encoder that maps $X_{n}$ one-to-one to $d$-dimensional binary
vectors, with $d=\lceil\log n\rceil$ and $n$ sufficiently large, there exists
a perfect decoder of a constant depth and width at most $\max(\lceil
n/B\rceil+B,BD)$ where $B$ is any integer such that $1<B\leq n$.

Note that if $B=2^{K}$ for some integer $K$, we do not need a condition of
“sufficiently large $n$” (precisely, $n\gg(\log n)^{k}$ for some constant
$k$). This condition is needed only for calculating $\lceil
Int(\beta)/B\rceil$ and ($Int(\beta)\mod B$), the details of which are not
relevant to this paper.

Figure 2: Network structure of a decoder with width $\max(\lceil
n/B\rceil+B,BD)$ for $B=2$.

Upon choosing $B=\lceil\sqrt{{\frac{n}{D}}}\rceil$ in Theorem 5 we obtain the
following corollary.

######  Corollary 6.

Suppose $n>D$. Then, for any perfect encoder that maps $X_{n}$ one-to-one to
$d$-dimensional binary vectors with $d=\lceil\log n\rceil$, there exists a
constant depth perfect decoder width $O(\sqrt{Dn})$ nodes.

Interestingly, this corollary improves the $O(D\sqrt{n})$ size of a perfect
decoder given in [8], while using a simpler architecture. Furthermore, the
order of this upper bound is close to that of the lower bound given in Theorem
4.

##  5 Approximate Decoders

In Section 4, we showed a construction of a perfect decoder of width
$\max(\lceil n/B\rceil+B,BD)$. It seems not easy to reduce the width using a
constant depth. Here we show that we can reduce the width if small errors
between the original input vectors (to the encoder) and the output vectors are
permissible. This assumption is reasonable because the input and output
vectors are not necessarily the same in practice.

In order to measure the error, we employ the Hamming distance between the
original vector ${\bf x}$ and the decoded vector ${\bf y}$, and denote it by
$dist({\bf x},{\bf y})$. The average Hamming distance $dist(X,Y)$ between the
sets of vectors $X=({\bf x}^{0},\ldots,{\bf x}^{n-1})$ and $Y=({\bf
y}^{0},\ldots,{\bf y}^{n-1})$ is defined as

$dist(X,Y)={\frac{1}{n}}\sum_{i=0}^{n-1}dist({\bf x}^{i},{\bf y}^{i}).$

###  5.1 Case of $B=3$

As in Section 4, we divide the input vectors into $B$ sets. In this
subsection, we consider the case of $B=3$ for ease of presentation. We will
extend it to a general $B$ in the next subsection.

We assume w.l.o.g. that ${\frac{n}{3}}$ is a positive integer and define
$n_{3}={\frac{n}{3}}$. We construct $\gamma_{i}$ nodes $(i=0,\ldots,n_{3}+2)$,
$y_{i,0}$ nodes $(i=0,\ldots,D-1$), and $y_{i,1}$ nodes $(i=0,\ldots,D-1)$ as
follows (see Fig. 3).

Recall that $Int(\beta)$ denotes the integer coded by a binary vector $\beta$.
Then, the activation function of $\gamma_{i}$ is defined by

$\displaystyle\gamma_{i}$ $\displaystyle=$
$\displaystyle[\lfloor{\frac{Int(\beta)}{3}}\rfloor=i],~{}~{}~{}~{}\mbox{for
$i=0,\ldots,n_{3}-1$},$ $\displaystyle\gamma_{n_{3}}$ $\displaystyle=$
$\displaystyle[Int(\beta)=0~{}(\mbox{mod}~{}3)],$
$\displaystyle\gamma_{n_{3}+1}$ $\displaystyle=$
$\displaystyle[Int(\beta)=1~{}(\mbox{mod}~{}3)],$
$\displaystyle\gamma_{n_{3}+2}$ $\displaystyle=$
$\displaystyle[Int(\beta)=2~{}(\mbox{mod}~{}3)].$

The activation function for $y_{j,h},j=0,\ldots,D-1,h=0,1$, is

$y_{j,h}=[\sum^{n_{3}+2}_{i=0}w^{j}_{i,h}\gamma_{i}\geq 0].$

The values of $w^{j}_{i,h}$s ($i=0,\ldots,n_{3}-1$, $h=0,1$) depend on
$x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}$ and are determined as in Table 1, where
each weight $w^{j}_{i,h}$ is shown by its binary representation (e.g.,
$w^{j}_{i,0}=4$ and $w^{j}_{i,1}=1$ if
$x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}=010$). Note that the corresponding values
of the $y_{j,h}$s and $z_{j}$s are also shown in the table. The values of
$w^{j}_{i,h}$s ($i=n_{3},n_{3}+1,n_{3}+2$, $h=0,1$) are defined by

$\displaystyle w^{j}_{{n_{3}},0}=-2,$ $\displaystyle
w^{j}_{{n_{3}+1},0}=-4,~{}~{}~{}w^{j}_{{n_{3}+2},0}=0,$ $\displaystyle
w^{j}_{{n_{3}},1}=-2,$ $\displaystyle
w^{j}_{{n_{3}+1},1}=0,~{}~{}~{}~{}w^{j}_{{n_{3}+2},1}=-4,$

where these $w^{j}_{i,h}$s do not depend on $j$.

Then, for each $j=0,\ldots,D-1$, we construct a node $z_{j}$ with the
activation function defined by

$z_{j}=\left[y_{j,0}+y_{j,1}\geq 2\right].$

Table 1: Definitions of weight $w^{j}_{i,h}$s and values of relevant variables. $x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}$ | $w^{j}_{i,0}$ | $w^{j}_{i,1}$ | $y^{3i}_{j,0}y^{3i+1}_{j,0}y^{3i+2}_{j,0}$ | $y^{3i}_{j,1}y^{3i+1}_{j,1}y^{3i+2}_{j,1}$ | $z^{3i}_{j}z^{3i+1}_{j}z^{3i+2}_{j}$
---|---|---|---|---|---
000 | 000 | 000 | 001 | 010 | 000
001 | 001 | 100 | 001 | 111 | 001
010 | 100 | 001 | 111 | 010 | 010
011 | 101 | 101 | 111 | 111 | 111
100 | 010 | 010 | 101 | 110 | 100
101 | 011 | 110 | 101 | 111 | 101
110 | 110 | 011 | 111 | 110 | 110
111 | 111 | 111 | 111 | 111 | 111

We use $w_{h}(\ldots)$ ($h=0,1$) to denote the functions given in the above.
For example, $w_{0}(010)=100_{(2)}=4$, $w_{1}(010)=001_{(2)}=1$,
$w_{0}(110)=110_{(2)}=6$, and $w_{1}(110)=011_{(2)}=3$, where $x_{(2)}$ means
that $x$ is a binary representation of an integer.

Figure 3: Network structure for approximate decoding with $B=3$.

It is seen from this table that an error occurs only when
$x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}=011$. Therefore, if all bit patterns are
distributed uniformly at random, the expected Hamming distance between ${\bf
x}^{i}$ and ${\bf y}^{i}$ is ${\frac{D}{3\cdot 2^{3}}}={\frac{D}{24}}$. We can
modify the network so that the average Hamming distance between ${\bf x}^{i}$
and ${\bf y}^{i}$ is always no more than ${\frac{D}{24}}$.

This can be done as follows. Let
$\\#_{b_{0}b_{1}b_{2}}=\sum_{j=0}^{D-1}|\\{i|x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}=b_{0}b_{1}b_{2}\\}|$
where $b_{0}b_{1}b_{2}\in\\{0,1\\}^{3}$. Let
$c_{0}c_{1}c_{2}=\mbox{argmin}_{b_{0}b_{1}b_{2}}\\#_{b_{0}b_{1}b_{2}}$, where
the tie can be broken in any way. Then, we let $a_{0}a_{1}a_{2}=(0\oplus
c_{0})(1\oplus c_{1})(1\oplus c_{2})$. For example. $a_{0}a_{1}a_{2}=011$ when
$c_{0}c_{1}c_{2}=000$.

Here we define $\chi_{a_{0}a_{1}a_{2}}(x_{0}x_{1}x_{2})$ by

$\chi_{a_{0}a_{1}a_{2}}(x_{0}x_{1}x_{2})=(x_{0}\oplus a_{0})(x_{1}\oplus
a_{1})(x_{2}\oplus a_{2}).$

For example, $\chi_{011}(000)=011$, $\chi_{011}(001)=010$, and
$\chi_{011}(011)=000$ when $a_{0}a_{1}a_{2}=011$. Then, we define
$w^{j}_{i,h}$ by

$\displaystyle w^{j}_{i,0}$ $\displaystyle=$ $\displaystyle
w_{0}(\chi_{a_{0}a_{1}a_{2}}(x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j})),$
$\displaystyle w^{j}_{i,1}$ $\displaystyle=$ $\displaystyle
w_{1}(\chi_{a_{0}a_{1}a_{2}}(x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j})).$

Finally, we define nodes $z_{j}^{\prime}$ in the output layer by

$z_{j}^{\prime}=z_{j}\oplus((a_{0}\gamma_{n_{3}})\lor(a_{1}\gamma_{n_{3}+1})\lor(a_{2}\gamma_{n_{3}+2})),$

which can be realized by adding nodes $\gamma^{\prime}$,
$\gamma^{\prime\prime}$, and $z_{j,h}$ ($j=0,\ldots,D-1,h=0,1$) and by
defining the activation functions as:

$\displaystyle\gamma^{\prime}$ $\displaystyle=$
$\displaystyle[a_{0}\gamma_{n_{3}}+a_{1}\gamma_{n_{3}+1}+a_{2}\gamma_{n_{3}+2}\geq
1],$ $\displaystyle\gamma^{\prime\prime}$ $\displaystyle=$
$\displaystyle[\gamma^{\prime}\geq 1],$ $\displaystyle z_{i,0}$
$\displaystyle=$ $\displaystyle[z_{i}-\gamma^{\prime\prime}\geq 1],$
$\displaystyle z_{i,1}$ $\displaystyle=$
$\displaystyle[-z_{i}+\gamma^{\prime\prime}\geq 1],$ $\displaystyle
z^{\prime}_{i}$ $\displaystyle=$ $\displaystyle[z_{i,0}+z_{i,1}\geq 1].$

For example, when $c_{0}c_{1}c_{2}=000$, we have the values of relevant
variables as in Table 2.

Table 2: Values of relevant variables for the case of $c_{0}c_{1}c_{2}=000$. $x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}$ | $\chi_{011}(x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j})$ | $z^{3i}_{j}z^{3i+1}_{j}z^{3i+2}_{j}$ | $(z^{\prime})^{3i}_{j}(z^{\prime})^{3i+1}_{j}(z^{\prime})^{3i+2}_{j}$
---|---|---|---
000 | 011 | 111 | 100
001 | 010 | 010 | 001
010 | 001 | 001 | 010
011 | 000 | 000 | 011
100 | 111 | 111 | 100
101 | 110 | 110 | 101
110 | 101 | 101 | 110
111 | 100 | 100 | 111

Indeed, if $x^{3i}_{j}x^{3i+1}_{j}x^{3i+2}_{j}=010$, we have
$z^{3i}_{j}z^{3i+1}_{j}z^{3i+2}_{j}=001$ and
$(z^{\prime})^{3i}_{j}(z^{\prime})^{3i+1}_{j}(z^{\prime})^{3i+2}_{j}=010$ by
the following:

$\displaystyle w^{j}_{i,0}$ $\displaystyle=$ $\displaystyle
w_{0}(\chi_{011}(010))=w_{0}(001)=001_{(2)}=1,$ $\displaystyle w^{j}_{i,1}$
$\displaystyle=$ $\displaystyle
w_{1}(\chi_{011}(010))=w_{1}(001)=100_{(2)}=4,$ [for $3i$:] $\displaystyle
y_{j,0}$ $\displaystyle=$ $\displaystyle[1-2\geq 0]=0,$ $\displaystyle
y_{j,1}$ $\displaystyle=$ $\displaystyle[4-2\geq 0]=1,$ $\displaystyle z_{j}$
$\displaystyle=$ $\displaystyle[0+1\geq 2]=0,$ $\displaystyle z_{j}^{\prime}$
$\displaystyle=$ $\displaystyle 0\oplus(0\cdot\gamma_{n_{3}})=0,$ [for
$3i+1$:] $\displaystyle y_{j,0}$ $\displaystyle=$ $\displaystyle[1-4\geq
0]=0,$ $\displaystyle y_{j,1}$ $\displaystyle=$ $\displaystyle[4-0\geq 0]=1,$
$\displaystyle z_{j}$ $\displaystyle=$ $\displaystyle[0+1\geq 2]=0,$
$\displaystyle z_{j}^{\prime}$ $\displaystyle=$ $\displaystyle
0\oplus(1\cdot\gamma_{n_{3}+1})=1,$ [for $3i+2$:] $\displaystyle y_{j,0}$
$\displaystyle=$ $\displaystyle[1-0\geq 0]=1,$ $\displaystyle y_{j,1}$
$\displaystyle=$ $\displaystyle[4-4\geq 0]=1,$ $\displaystyle z_{j}$
$\displaystyle=$ $\displaystyle[1+1\geq 2]=1,$ $\displaystyle z_{j}^{\prime}$
$\displaystyle=$ $\displaystyle 1\oplus(1\cdot\gamma_{n_{3}+2})=0.$

If $n$ can be divided by 3, the total error will be at most ${\frac{nD}{3\cdot
2^{3}}}={\frac{nD}{24}}$. Otherwise, the total error will be at most
${\frac{nD}{3\cdot 2^{3}}}={\frac{nD}{24}}+D$ because there might be an
additional error per $j$. Since the maximum width is given by $\max(\lceil
n/3\rceil+3,2D+1)$, we have the following theorem.

######  Theorem 7.

For any perfect encoder that maps $X_{n}$ one-to-one to $d$-dimensional binary
vectors, with $d=\lceil\log n\rceil$ and $n$ sufficiently large, there exists
a decoder of a constant depth and width $\max(\lceil n/3\rceil+3,2D+1)$ whose
average Hamming distance error is at most $D({\frac{1}{24}}+{\frac{1}{n}})$.

###  5.2 General Case

As in Section 4, we define $n_{B}=n/B$ and assume w.l.o.g. that $n$ can be
divided by $B$. First, we define $\gamma_{i}$ nodes ($i=0,\ldots,n_{B}+B-1$)
by

$\displaystyle\gamma_{i}$ $\displaystyle=$
$\displaystyle[\lfloor{\frac{Int(\beta)}{B}}\rfloor=i],~{}~{}~{}~{}\mbox{for
$i=0,\ldots,n_{B}-1$},$ $\displaystyle\gamma_{n_{B}+k}$ $\displaystyle=$
$\displaystyle[Int(\beta)=k~{}(\mbox{mod}~{}B)],~{}~{}~{}~{}\mbox{for
$k=0,\ldots,B-1$}.$

We construct $y_{j,h}$ nodes ($j=0,\ldots,D-1$, $h=0,\ldots,B-2$) and $z_{j}$
nodes ($j=0,\ldots,D-1$) with the activation functions:

$\displaystyle y_{j,h}$ $\displaystyle=$
$\displaystyle\left[\sum_{i=0}^{n_{3}+B-1}w^{j}_{i,h}\gamma_{i}\geq\theta_{h}\right],$
$\displaystyle z_{j}$ $\displaystyle=$
$\displaystyle\left[\sum_{h=0}^{B-2}y_{0,h}~{}\geq~{}B-1\right].$

We define $w_{n_{B}+i,h}$s ($i=0,\ldots.B-1$, $h=0,\ldots,B-2$) by

$\displaystyle w_{n_{B}+i,h}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}-{\frac{2^{B}}{4}},&\mbox{if $i=0$
and $h\in\\{0,\ldots,B-2\\}$},\\\ -{\frac{2^{B}}{2}},&\mbox{if $i=h+1$ and
$h\in\\{0,\ldots,B-2\\}$},\\\ 0,&\mbox{otherwise},\end{array}\right.$

For integer $k\in\\{0,\ldots,2^{B}-1\\}$, let $k_{0}k_{1}\cdots k_{B-1}$ be
the binary representation of $k$ (i.e., $k=Int(k_{0}k_{1}\cdots k_{B-1})$),
where $k_{0}$ is the most significant bit. We define $w_{h}(k)$
($h=0,\ldots,B-2$) by

$\displaystyle w_{0}(k_{0}k_{1}k_{2}k_{3}\cdots k_{B-1})$ $\displaystyle=$
$\displaystyle Int(k_{1}k_{0}k_{2}k_{3}\cdots k_{B-1}),$ $\displaystyle
w_{1}(k_{0}k_{1}k_{2}k_{3}\cdots k_{B-1})$ $\displaystyle=$ $\displaystyle
Int(k_{2}k_{0}k_{1}k_{3}\cdots k_{B-1}),$ $\displaystyle
w_{2}(k_{0}k_{1}k_{2}k_{3}\cdots k_{B-1})$ $\displaystyle=$ $\displaystyle
Int(k_{3}k_{0}k_{1}k_{2}\cdots k_{B-1}),$ $\displaystyle\cdots$ $\displaystyle
w_{B-2}(k_{0}k_{1}k_{2}k_{3}\cdots k_{B-1})$ $\displaystyle=$ $\displaystyle
Int(k_{B-1}k_{0}k_{1}\cdots k_{B-2}).$

Then, we define $w^{j}_{i,h}=w_{h}(x^{Bi}_{j}x^{Bi+1}_{j}\cdots
x^{Bi+B-1}_{j})$ for $j=0,\ldots,D-1$, $i=0,\ldots,n_{B}-1$, and
$h=0,\ldots,B-2$. The following proposition is straightforward from the
definitions of $y_{i,h}$s, $w_{n_{B}+i,h}$s, $w_{h}(\cdots)$s, and $z_{j}$s.

######  Proposition 8.

If $x^{Bi}_{j}x^{Bi+1}_{j}x^{Bi+2}_{j}\cdots x^{Bi+B-1}_{j}=011\cdots 1$, then
$z^{Bi}_{j}z^{Bi+1}_{j}z^{Bi+2}_{j}\cdots z^{Bi+B-1}_{j}=111\cdots 1$ holds.
Otherwise, $z^{Bi}_{j}z^{Bi+1}_{j}z^{Bi+2}_{j}\cdots
z^{Bi+B-1}_{j}=x^{Bi}_{j}x^{Bi+1}_{j}x^{Bi+2}_{j}\cdots x^{Bi+B-1}_{j}$ holds.

As in Section 5.1, it is seen from this proposition that if all bit patterns
are distributed uniformly at random, the expected Hamming distance between
${\bf x}^{i}$ and ${\bf y}^{i}$ is ${\frac{D}{B2^{B}}}$. We can modify the
network so that the average Hamming distance between ${\bf x}^{i}$ and ${\bf
y}^{i}$ is always no more than ${\frac{D}{B2^{B}}}$ as well. The modification
method is given below, which is almost the same as that in Section 5.1.

Let $\\#_{b_{0}b_{1}\cdots
b_{B-1}}=\sum_{j=0}^{D-1}|\\{i|x^{Bi}_{j}x^{Bi+1}_{j}\cdots
x^{Bi+B-1}_{j}=b_{0}b_{1}\cdots b_{B-1}\\}|$ where $b_{0}b_{1}\cdots
b_{B-1}\in\\{0,1\\}^{B}$. Let $c_{0}c_{1}\cdots
c_{B-1}=\mbox{argmin}_{b_{0}b_{1}\cdots b_{B-1}}\\#_{b_{0}b_{1}\cdots
b_{B-1}}$, where the tie can be broken in any way. Then, we let
$a_{0}a_{1}\cdots a_{B-1}=(0\oplus c_{0})(1\oplus c_{1})\cdots(1\oplus
c_{B-1})$. For example. $a_{0}a_{1}a_{2}a_{3}=0110$ when
$c_{0}c_{1}c_{2}c_{3}=0001$.

Here we define $\chi_{a_{0}a_{1}\cdots a_{B-1}}(x_{0}x_{1}\cdots x_{B-1})$ by

$\chi_{a_{0}a_{1}\cdots a_{B-1}}(x_{0}x_{1}\cdots x_{B-1})=(x_{0}\oplus
a_{0})(x_{1}\oplus a_{1})\cdots(x_{B-1}\oplus a_{B-1}).$

Then, we define $w^{j}_{i,h}$ by

$\displaystyle w^{j}_{i,h}$ $\displaystyle=$ $\displaystyle
w_{h}(\chi_{a_{0}a_{1}\cdots a_{B-1}}(x^{Bi}_{j}x^{Bi+1}_{j}\cdots
x^{Bi+B-1}_{j})),$

for $j=0,\ldots,D-1$, $i=0,\ldots,n_{B}-1$, and $h=0,\ldots,B-2$. Finally, we
define nodes $z_{j}^{\prime}$ in the output layer by

$z_{j}^{\prime}=z_{j}\oplus(\bigvee_{h=0}^{B-1}a_{h}\gamma_{n_{B}+h}),$

which can be realized by adding nodes $\gamma^{\prime}$,
$\gamma^{\prime\prime}$, and $z_{j,h}$ ($j=0,\ldots,D-1,h=0,1$) and defining
the activation functions by:

$\displaystyle\gamma^{\prime}$ $\displaystyle=$
$\displaystyle[\sum_{h=0}^{B-1}a_{h}\gamma_{n_{B}+h}~{}\geq~{}1],$
$\displaystyle\gamma^{\prime\prime}$ $\displaystyle=$
$\displaystyle[\gamma^{\prime}\geq 1],$ $\displaystyle z_{i,0}$
$\displaystyle=$ $\displaystyle[z_{i}-\gamma^{\prime\prime}\geq 1],$
$\displaystyle z_{i,1}$ $\displaystyle=$
$\displaystyle[-z_{i}+\gamma^{\prime\prime}\geq 1],$ $\displaystyle
z^{\prime}_{i}$ $\displaystyle=$ $\displaystyle[z_{i,0}+z_{i,1}\geq 1].$

Since the maximum width is given by $\max(\lceil n/B\rceil+B,(B-1)D+1)$, we
have the following theorem.

######  Theorem 9.

For any perfect encoder that maps $X_{n}$ one-to-one to $d$-dimensional binary
vectors, with $d=\lceil\log n\rceil$ and sufficiently large $n$, there exists
a decoder of a constant depth and width $\max(\lceil n/B\rceil+B,(B-1)D+1)$
whose average Hamming distance error is at most
$D({\frac{1}{B2^{B}}}+{\frac{1}{n}})$, where $B$ is any integer such that
$3\leq B\leq n$.

Note that this theorem is a generalized version of Theorem 7 and the resulting
width is smaller than that of Theorem 5 for all $B\geq 3$ and $\lceil
n/B\rceil+B>BD$ (i.e., $D<n^{\prime}$ for $n^{\prime}\approx n/B^{2}$)
although the corresponding decoder is an approximate one and uses more layers.
Note also that if $B=2^{K}$, the architecture of the decoder can be explicitly
described as $d/\lceil n/B\rceil+B/(B-1)D+1/D+1/2D/D$ where $d=\lceil\log
n\rceil$. Furthermore, if we can consider the average case error over binary
input vectors given uniformly at random, this architecture can be reduced to
$d/\lceil n/B\rceil+B/(B-1)D/D$.

##  6 Concluding Remarks

In this paper, we showed an improved upper bound on the size/width of the
decoder in a Boolean threshold autoencoder. We also showed a lower bound on
the size of the decoder. Although these bounds are relatively close, there
still exists a gap. Therefore, closing the gap is left as an open problem as
well as giving a lower bound on the size/width of the encoder part. We also
constructed decoders that are permitted to make small errors. However, the
widths of these decoders are not significantly smaller than those of exact
decoders. Therefore, improvement of decoders allowing errors is left as future
work.

## References

  * [1] David H. Ackley, Geoffrey E. Hinton, and Terrence J. Sejnowski.  A learning algorithm for boltzmann machines.  Cognitive Science, 9:147–169, 1985.
  * [2] Martin Anthony.  Discrete Mathematics of Neural Networks, Selected Topics.  SIAM, 2001.
  * [3] Pierre Baldi.  Autoencoders, unsupervised learning, and deep architectures.  JMLR: Workshop and Conference Proceedings, 27:37–50, 2012.
  * [4] Pierre Baldi and Kurt Hornik.  Neural networks and principal component analysis:learning from examples without local minima.  Neural Networks, 2:53–38, 1989.
  * [5] Carl Doersch.  Tutorial on variational autoencoders.  arXiv, 1606.05908, 2016.
  * [6] Rafael Gómez-Bombarelli, Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Sánchez-Lengeling, Dennis Sheberla, Jorge Aguilera-Iparraguirre, Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik.  Automatic chemical design using a data-driven continuous representation of molecules.  ACS Central Science, 4(2):268–276, 2018.
  * [7] Geoffrey E. Hinton and Ruslan R. Salakhutdinov.  Reducing the dimensionality of data with neural networks.  Science, 313:504–507, 2006.
  * [8] Avraham A. Melkman, Sini Guo, Wai-Ki Ching, Pengyu Liu, and Tatsuys Akutsu.  On the compressive power of boolean threshold autoencoders.  IEEE Trans. Neural Networks and Learning Systems, in press.
  * [9] Kai-Yeung Siu, Jehoshual Bruck, Thomas Kailath, and Thomas Hofmeister.  Depth efficient neural networks for division and related problems.  IEEE Trans. Information Theory, 39:946–956, 1993.
  * [10] Kai-Yeung Siu, Vwani Rndoychowdhury, and Thomas Kailath.  Discrete Neural Computation. A Theoretical Foundation.  Prentice-Hall, 1995.
  * [11] Michael Tschannen, Olivier Bachem, and Mario Lucic.  Recent advances in autoencoder-based representation learning.  arXiv, 1812.05069, 2018.
  * [12] Jiří Šíma and Pekka Orponen.  General-purpose computation with neural networks: a survey of complexity theoretic results.  Neural Computation, 15:2727–2778, 2003.