# Learned Optimizers for Analytic Continuation

Dongchen Huang  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100049, China  Yi-feng Yang
<EMAIL_ADDRESS>Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100049, China  Songshan
Lake Materials Laboratory, Dongguan, Guangdong 523808, China

###### Abstract

Traditional maximum entropy and sparsity-based algorithms for analytic
continuation often suffer from the ill-posed kernel matrix or demand
tremendous computation time for parameter tuning. Here we propose a neural
network method by convex optimization and replace the ill-posed inverse
problem by a sequence of well-conditioned surrogate problems. After training,
the learned optimizers are able to give a solution of high quality with low
time cost and achieve higher parameter efficiency than heuristic full-
connected networks. The output can also be used as a neural default model to
improve the maximum entropy for better performance. Our methods may be easily
extended to other high-dimensional inverse problems via large-scale
pretraining.

Inverse problems appear in many perspectives of physics and machine learning,
such as learning Hamiltonian in the classical Chow1968 ; Hinton1986 ;
Bresler2015  or quantum sense Anshu2020 ; Bairey2019 ; Amin2018  and
recovering sparse signal from noise measurements Donoho2006 ; Candes2005 . In
quantum many-body problems, correlation functions are often computed in
imaginary time Gull2011 ; Gubernatis2016  so that an analytic continuation has
to be implemented to obtain the spectral function in real frequency in order
to extract meaningful information. The analytic continuation is nothing but a
linear inverse problem, which is, however, highly ill-posed and may often
yield unphysical solutions.

Many algorithms have been proposed to attack this problem. These include padé
approximation, stochastic methods Sandvik1998 ; Sandvik2016 ; Ghanem2020a ;
Ghanem2020b , maximum entropy methods JARRELL1996 ; Sim2018 ; Kraberger2017 ;
Silver1990 ; Gunnarsson2010b ; Reymbaut2015 ; Rumetshofer2019 , and the
Nevanlinna method Fei2021 . Classical methods such as the padé approximation
Deisz1996  and the singular value decomposition (SVD) Gunnarsson2010  usually
perform well in many tasks. The hyper parameter selection of maximum entropy
methods can also be adjusted in various ways Bryan1990 ; Bergeron2016 . But
none of them holds for all situations and a case-by-case tuning is often
needed. For instance, maximum entropy methods demand a highly empirical
selection of prior distributions.

From the view of representation learning, high dimensional data of real world
always have certain low-dimensional structures. The spectral function may have
a sparse structure if it is properly discretized. This has motivated a line of
work focusing on sparsity of the spectral function Ohzeki2018 ; YOSHIMI2019 ;
Otsuki2020 ; Otsuki2017 ; Shinaoka2017 . Unfortunately, the power of sparsity-
based methods is greatly limited by the ill-poseness of the Fermi kernel
matrix.

In this work, we develop a neural network architecture for analytic
continuation by transforming a highly ill-posed basis pursuit problem into a
sequence of well-conditioned surrogate problems. The problem sequence can be
implemented by introducing a learned optimizer with given data, which may be
viewed as a neural network whose structure can be directly derived from convex
optimization. This neural network avoids empirical design and shows higher
parameter efficiency compared with heuristic full-connected neural networks
used in other works Fournier2020 ; Yoon2018 . It can give a high quality
approximate solution of the linear inverse problem with much less time cost
than traditional maximum entropy methods. Moreover, the two approaches can
complement each other by taking advantage of neural network’s strengths and
treating its output as prior distributions of the maximum entropy, thus
yielding an improved solution with better precision.

## Analytic continuation.

We are dealing with the inverse problem to obtain the spectral function
$A(\omega)$ in real frequency from the Green’s function $G(\tau)$,

$G(\tau)=\int_{-\infty}^{\infty}K(\tau,\omega)A(\omega)d\omega,$ (1)

where $\tau$ is the imaginary time. For a fermionic Green’s function, the
kernel $K(\tau,\omega)$ takes the form,

$K(\tau,\omega)=\frac{e^{-\tau\omega}}{1+e^{-\beta\omega}},$ (2)

where $\beta$ is the inverse temperature. For analytic continuation, we first
discretize Eq. (1) and get the linear inverse problem (in matrix form):

$g(\tau_{i})=K(\tau_{i},\omega_{j})a(\omega_{j}),$ (3)

where $i=1,\dots,N_{\tau}$ and $j=1,\dots,N_{\omega}$ mark the discrete points
in imaginary time and real frequency, respectively, $g$ is the vectorized
Green’s function, $a$ is the vectorized spectral function, and $K(.,.)$ is a
matrix of the Fermi kernel. The problem can then be solved using the Bayesian
inference.

## Maximum entropy and sparsity-based methods.

The posterior distribution of the spectral function $a$ satisfies the Bayes’
theorem,

$P(a|g)=\frac{P(g|a)P(a)}{P(g)},$ (4)

from which $a$ can be derived by maximum likelihood. The denominator $P(g)$ is
independent of $a$ and can be omitted in optimization. The maximum entropy and
sparsity-based methods are two special forms of its implementation with
different choices of $P(g|a)$ and $P(a)$.

The traditional maximum entropy methods choose $P(g|a)=-\chi^{2}/2$ and
$P(a)\propto e^{-\alpha S}$, thus requiring

$a=\operatorname*{\arg\min}_{a}\frac{\chi^{2}}{2}+\alpha S,$ (5)

where $\chi^{2}=(g-Ka)^{T}\Sigma^{-1}(g-Ka)$ denotes the reconstruction error,
$\Sigma$ is the empirical covariance matrix, $K$ is the kernel matrix, and
$S=\sum_{i}\Delta\omega_{i}a(\omega_{i})\log\frac{a(\omega_{i})}{d(\omega_{i})}$
is the Kullback-Leibler (KL) divergence between the spectral function and a
prior default model $d(.)$.

The sparsity-based methods assume a Gaussian distribution with unit variance
such that $P(g|a)\propto\exp(-\frac{1}{2}\|g-Ka\|_{2}^{2})$. The prior
distribution $P(a)$ is typically chosen by hand to reflect the sparse
regularization. For Laplacian prior distribution, we have $P(a)\propto
e^{-\lambda\left\|a\right\|_{1}}(\lambda>0)$, where $\left\|.\right\|_{1}$ is
the $\ell_{1}$ norm defined as the summation of absolute values and $\lambda$
is the hyperparameter. This leads to a basis pursuit (BP) problem,

$a=\operatorname*{\arg\min}_{a}\frac{1}{2}\|g-Ka\|_{2}^{2}+\lambda\|a\|_{1}.$
(6)

Ideally, $\ell_{0}$ regularization should be used that counts the non-zero
entrices in a vector, but minimizing the $\ell_{0}$ norm is NP-hard
Natarajan1995 . Fortunately, the $\ell_{1}$ norm provides a good replacement
which is the largest convex function to approximate $\ell_{0}$. The sparsity
assumption and such kind of problems have achieved huge amount of successes in
machine learning and signal processing. The BP problem is well-posed and
guaranteed to recovery the exact spectral function as long as $a$ is sparse
and the kernel matrix $K$ satisfies some fine properties. The sparest solution
is favored in this method, in contrast to maximum entropy methods that look
for a solution of least deviation from the default model.

Figure 1: Architecture of the 3-layer RLISTA network with input of the Green’s
function $g$ and a zero vector $a_{1}$. The blue arrows illustrate the forward
propagation of the optimization process, and the purple arrow indicates the
backward propagation which tunes all parameters to feed the data.

## Solution of the BP problem.

Since Eq. (6) is convex, many popular methods can be applied and all of them
can converge to the only global minimum. Here we map it to the fixed-point
problem Boyd2004 ,

$a^{*}=S_{\tau\lambda}(a^{*}-\tau K^{T}(Ka^{*}-g)),$ (7)

where the soft-thresholding operator $S_{\tau\lambda}(.)$ is the proximal
mapping of the function $\tau\lambda\left\|.\right\|_{1}$ and takes the form,

$S_{\tau\lambda}(y)\equiv\left\\{\begin{array}[]{ll}y-\tau\lambda&\ \ \ \ \
y>\tau\lambda\\\ 0&\ \ \ \ \ |y|\leq\tau\lambda\\\ y+\tau\lambda&\ \ \ \ \
y<-\tau\lambda\end{array}\right..$ (8)

To see this, we note that a solution $a^{*}$ of the BP problem must satisfy
the optimality condition $0\in
K^{T}(Ka^{*}-g)+\lambda\partial\left\|a^{*}\right\|_{1}$, where $\partial(.)$
denotes the subdifferential. Thus, for any $\tau>0$, we have $a^{*}-\tau
K^{T}(Ka^{*}-g)\in a^{*}+\tau\lambda\partial\left\|a^{*}\right\|_{1}$. On the
other hand, for any convex function
$F:\mathbb{R}^{N}\rightarrow(-\infty,\infty]$ and its induced proximal mapping
$P_{F}(z)=\operatorname*{\arg\min}_{x}F(x)+\frac{1}{2}\left\|x-z\right\|_{2}^{2}$,
the identity $x=P_{F}(z)$ implies $z\in x+\partial F(x)$. Combining these two
immediately gives $F(x)=\tau\lambda\left\|x\right\|_{1}$ and the above fixed-
point equations. It is now understood that the unit variance assumption in Eq.
(6) corresponds to a convenient usage of the proximal mapping.

The fixed-point problem has a natural iteration scheme:

$\begin{split}a_{n+1}&=S_{\tau\lambda}(a_{n}-\tau
K^{T}(Ka_{n}-g)),\end{split}$ (9)

which is also called the iterative shrinkage-thresholding algorithm (ISTA)
Daubechies2004 ; Figueiredo2003 . Unfortunately, for analytic continuation, it
fails to converge to the physical fixed point because of the ill-poseness of
the Fermi kernel.

## Learned optimizers.

To overcome this issue, we note that ISTA can also be viewed as recurrent
neural networks (RNNs) with fixed weights determined by the Fermi kernel
matrix. This motivates us to design a neural network structure by unrolling
the fixed-point iteration and utilize the simple forward problem of Eq. (3) to
generate data, train the own weights for each individual layer, and learn an
adaptive optimizer. We propose the following learnable iterative soft
thresholding algorithm (LISTA) and its residual version (RLISTA) with a
$L$-layer neural network of fixed depth:

$\begin{split}a_{n+1}&=(1-\eta)a_{n}+\eta
S_{\tau\lambda}(W_{t}^{n}a_{n}+W_{e}^{n}g),\\\
a_{n+1}&=S_{\tau\lambda}(W_{t}^{n}a_{n}+W_{e}^{n}g),~{}~{}~{}n=1,2,\dots,L.\end{split}$
(10)

where $a_{1}$ and $g$ are inputs of the neural network and $0<\eta<1$ is a
relaxation factor. The parameters $W_{e}^{n}$ and $W_{t}^{n}$ represent the
weights to be learned on the $n$-th layer to replace the fixed weights
$W_{e}=\tau K^{T}$ and $W_{t}=I-\tau K^{T}K$ in the original problem. During
the training and inference processes, the neural networks are fed with the
Green’s function $g$ and a zero vector as $a_{1}$. The overall scheme is shown
in Fig. 1 containing naturally a residual connection He2016  and independent
weights on each layer. We call the neural network RLISTANet. Historically,
LISTA has been proposed for sparse coding Gregor2010 .

## Performance and parameter efficiency.

Our learned optimizers, LISTA and RLISTA, have higher performance than vanilla
ISTA, as compared in Fig. 2(a). For simple spectra with only one sharp peak
near the origin, all three optimizers can recover the solution well. But for
complex spectra containing more broad peaks, ISTA can only give a single sharp
peak, while both LISTA and RLISTA can produce the ground truth with high
accuracy.

Figure 2: (a): Comparison of RLISTA, LISTA and ISTA for simulated data
generated from the spectral function with only one sharp peak (left) and more
peaks (right). ISTA can only recover the single peak structure, while LISTA
and RLISTA work well in both cases. (b) Comparison of parameter efficiency for
four different neural network architectures: LISTA, RLISTA and fully-connected
network of two (FCN-2) or three (FCN-3) layers.

Our neural networks may be viewed as a variation of the fully-connected
network (FCN) but show higher parameter efficiency. This is seen by comparison
with the conventional FCN of one or two hidden layers. For simplicity, the
width of single FCN is set equal. The neural performance is measured by the
root square error (RSE),

$\text{\rm
RSE}(\hat{a})=\sqrt{\sum_{i=1}^{N_{\omega}}(a_{i}^{*}-\hat{a}_{i})^{2}},$ (11)

where $a^{*}$ is the ground truth spectral function in the test set and
$\hat{a}$ is the prediction of the neural network. As shown in Fig. 2(b), our
6-layer LISTANet can already achieve better accuracy than 2-layer FCN of ten
times more parameters and 3-layer FCN of seven times more parameters. However,
deeper LISTANet does not necessarily outperform shallow ones, possibly due to
the landscape of networks Li2018 . Deep FCNs are known difficult to train and
require more advanced techniques like normalization Ioffe2015 ; Ba2016 ;
Wu2018 . By contrast, RLISTA contains residual connection He2016  and allows
for the training of much deeper networks. Although shallow RLISTA cannot
outperform LISTA, we can always train a deeper (40-layer) RLISTA that beats
all other three. However, when the depth is larger than $20$, adding layers no
longer reduces the error, reflecting a possible bottleneck of RLISTA.
Relaxation methods are usually known to accelerate convergence and improve
stability. Exploring the connection between acceleration and convergence
conditions and the landscape of corresponding neural networks may be
beneficial for designing better architecture and optimization techniques and
help improve the performance and achieve better accuracy.

Figure 3: (a) Comparison of the first layer matrices $W_{t}$ and $W_{e}$ in
ISTA, LISTA and RLISTA optimizers. The color bars are normalized by their
Frobenius norm $\|W_{e}\|_{F}^{2}$ and $\|W_{t}\|_{F}^{2}$. Before
normalization, an identity matrix is subtracted from the matrix $W_{t}$, while
the learned $W_{e}$ and all matrices related to ISTA are multiplied by $-1$.
(b) Comparison of the average coherence of their learned matrices for all
layers, showing the benefits of learning in LISTA and RLISTA.

## The weight matrices.

What is the reason behind the substantial improvement in learned optimizers?
Noticing that ISTA may be viewed as an infinite-layer neural network with
fixed weight, the improvement must originate from the adaptation to the data
as reflected in the learned weight matrices. Figure 3 compares the normalized
learned matrices $W_{t}$ and $W_{e}$ on first layer for different optimizers.
While the matrices in ISTA are determined directly by the Fermi kernel, the
learned matrices in LISTA and RLISTA are heavily influenced by their different
iteration schemes. All entries of $W_{e}$ in LISTA distribute uniformly and
are of the same magnitude, while those in RLISTA differ heavily and contain
some entries of relatively larger values. Nevertheless, both types of learned
matrices can empirically solve the inverse problem and perform nicer than
their vanilla cousin ISTA.

To quantify the “niceness” of the learned matrices, we use the average
coherence Bajwa2010 ; Mixon2011  for a matrix
$A=[a_{1}|\dots|a_{n}]\in\mathbb{R}^{m\times n}$:

$\nu(A)=\frac{1}{n-1}\max_{i\in\\{1,\dots,n\\}}\Big{|}\sum_{j\neq
i,j=1}^{n}\left\langle\frac{a_{i}}{\left\|a_{i}\right\|_{2}},\frac{a_{j}}{\left\|a_{j}\right\|_{2}}\right\rangle\Big{|},$
(12)

which gives the largest average inner product of two columns and measures the
spread of column vectors of a matrix within a unit ball. The average coherence
values in $[0,1]$, becomes zero for an orthogonal matrix, and reaches its
upper bound for a matrix with repeated columns. Obviously, a lower average
coherence implies a “nicer” or less singular matrix. As shown in Fig. 3, the
average coherence of learned matrices values around the lower bound zero,
while the predefined matrices in ISTA have a large average coherence close to
0.5. The gap between them illustrates the benefits of learning.

The diversity of learned matrices in different optimizers might be understood
as the diversity of learned problems. Each layer of LISTA and RLISTA networks
works approximately as an independent iteration step for solving a given BP
problem of Eq. (6) with their individual learned kernel matrix. But a single
iteration can not solve the problem with high accuracy. With many learned
layers, the neural networks learn a sequence of optimization problems and
solve each of them imprecisely.

From the perspective of optimization, our method may be understood as a
generalization of the classical homotopy method Donoho2008  which solves the
BP problem by a sequence of problems with varying regularization parameter
$\lambda$. Here the neural network solves the problem by a series of problems
with varying $K$ and $\lambda$ in the region of physical Green’s functions,
which is much smaller than the whole space $\mathbb{R}^{N_{\tau}}$. Thus, a
trade-off comes from the fixed depth: how precisely do we need to solve a
single question (precision) and how many optimization problems do we need to
learn (diversity)? Careful studies of the precision-diversity trade-off may
reveal a closer relationship between optimization and neural network, which we
leave for future work.

Figure 4: Two examples comparing the maximum entropy method with a flat
default model (Flat) and a neural default model (RLISTA+). The dots are the
ground truth spectral function. The RLISTA results are generated by a sub-
optimal 6-layer network.

## A neural default model.

The neural networks can give an answer with low time cost after training, but
they are often criticized as being a “black-box”, where precise
interpretations of the weights remain unclear. Better neural network
architectures and optimization techniques can no doubt reduce the error, but
their design suffers from high cost trial and error. On the other hand, a
classical maximum entropy problem requires a default model $d(\omega)$ to
incorporate certain prior knowledge about the desired spectral function. A
better default model allows for easy hyperparameter tuning and gives better
accuracy. We propose that combining the two may provide a novel way to improve
the performance. The output of our neural network is also a probability
distribution and may be considered as a “neural default model” required in the
maximum entropy. This fixes the inexactness of the neural network and can
benefit from both the high speed of neural networks and the well-developed
algorithms of the maximum entropy.

To get an impression on the advantage of such combination, we first obtain an
inexact result from a sub-optimal 6-layer RLISTA. Figure 4 shows two examples
where both neural network and maximum entropy (with flat prior) cannot capture
all the details of the spectral function. While the neural network tends to
average two peaks, the maximum entropy method ignores the peak at high
frequency. By contrast, using the neural default model, the maximum entropy
solution can capture well the high frequency peak and even all the details of
the exact spectral function. This neural maximum entropy method (RLISTA+)
makes a promising improvement over the conventional one for analytic
continuation.

## Discussion and conclusions.

Motivated by sparsity-based methods, we have proposed a highly efficient
neural network scheme for analytic continuation in quantum many-body problems.
Our learned optimizers show low time costs and may be easily extended to other
high-dimensional inverse problems via large-scale pretraining where
traditional maximum entropy methods demand tremendous computation time. The
output of our method may also be used as a neural default model to improve the
performance of maximum entropy methods and make use of their both advantages.
We also find that constructing neural networks from fix-point iteration can
achieve better parameter efficiency than heuristic full-connected networks. By
viewing neural network as learnable fix-point iteration, we see that different
fix-point iteration schemes are not equivalent if their parameters can be
learned, despite that they all converge to the same solution for a given
problem. The learned iteration paths are more regular than their unlearned
counterpart. This, combined with the powerful theory about calculus of
variations, may help invent novel algorithms and lead to a better
understanding of the training dynamics and regularization of neural networks.

This work was supported by the National Key R&D Program of MOST of China
(Grant No. 2017YFA0303103), the National Natural Science Foundation of China
(NSFC Grant No. 11974397, No. 11774401), the Strategic Priority Research
Program of the Chinese Academy of Sciences (Grant No. XDB33010100), and the
Youth Innovation Promotion Association of CAS.

## References

  * (1) C. Chow and C. Liu, Approximating discrete probability distributions with dependence trees, IEEE Transactions on Information Theory, 14. 462 (1968).
  * (2) G. E. Hinton and T. J. Sejnowski, Learning and Relearning in Boltzmann Machines, p. 282–317.  Cambridge, MA, USA: MIT Press, 1986.
  * (3) G. Bresler, Efficiently learning Ising models on arbitrary graphs, in STOC’15, (New York, NY, USA), Association for Computing Machinery, 2015\.
  * (4) A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Sample-efficient learning of quantum many-body systems, in 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), (2020).
  * (5) E. Bairey, I. Arad, and N. H. Lindner, Learning a local Hamiltonian from local measurements, Phys. Rev. Lett. 122. 020504 (2019).
  * (6) M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Quantum Boltzmann machine, Phys. Rev. X 8. 021050 (2018).
  * (7) D. Donoho, M. Elad, and V. Temlyakov, Stable recovery of sparse overcomplete representations in the presence of noise, IEEE Transactions on Information Theory, 52. 6 (2006).
  * (8) E. Candes, M. Rudelson, T. Tao, and R. Vershynin, Error correction via linear programming, in 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05), (2005).
  * (9) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. 83. 349 (2011).
  * (10) J. Gubernatis, N. Kawashima, and P. Werner, Quantum Monte Carlo Methods: Algorithms for Lattice Models.  Cambridge University Press, 2016.
  * (11) A. W. Sandvik, Stochastic method for analytic continuation of quantum monte carlo data, Phys. Rev. B 57. 10287 (1998).
  * (12) A. W. Sandvik, Constrained sampling method for analytic continuation, Phys. Rev. E 94. 063308 (2016).
  * (13) K. Ghanem and E. Koch, Average spectrum method for analytic continuation: Efficient blocked-mode sampling and dependence on the discretization grid, Phys. Rev. B 101. 085111 (2020).
  * (14) K. Ghanem and E. Koch, Extending the average spectrum method: Grid point sampling and density averaging, Phys. Rev. B 102. 035114, 2020\.
  * (15) M. Jarrell and J. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo data, Physics Reports, 269. 133 (1996).
  * (16) J.-H. Sim and M. J. Han, Maximum quantum entropy method, Phys. Rev. B, 98. 205102 (2018).
  * (17) G. J. Kraberger, R. Triebl, M. Zingl, and M. Aichhorn, Maximum entropy formalism for the analytic continuation of matrix-valued Green’s functions, Phys. Rev. B 96. 155128 (2017).
  * (18) R. N. Silver, D. S. Sivia, and J. E. Gubernatis, Maximum-entropy method for analytic continuation of quantum Monte Carlo data, Phys. Rev. B 41. 2380 (1990).
  * (19) O. Gunnarsson, M. W. Haverkort, and G. Sangiovanni, Analytical continuation of imaginary axis data using maximum entropy, Phys. Rev. B 81. 155107 (2010).
  * (20) A. Reymbaut, D. Bergeron, and A.-M. S. Tremblay, Maximum entropy analytic continuation for spectral functions with nonpositive spectral weight, Phys. Rev. B 92. 060509 (2015).
  * (21) M. Rumetshofer, D. Bauernfeind, and W. von der Linden, Bayesian parametric analytic continuation of Green’s functions, Phys. Rev. B 100. 075137 (2019).
  * (22) J. Fei, C.-N. Yeh, and E. Gull, Nevanlinna analytical continuation, Phys. Rev. Lett. 126. 056402 (2021).
  * (23) J. J. Deisz, D. W. Hess, and J. W. Serene, Incipient antiferromagnetism and low-energy excitations in the half-filled two-dimensional Hubbard model, Phys. Rev. Lett. 76. 1312 (1996).
  * (24) O. Gunnarsson, M. W. Haverkort, and G. Sangiovanni, Analytical continuation of imaginary axis data for optical conductivity, Phys. Rev. B 82. 165125 (2010).
  * (25) R. K. Bryan, Solving Oversampled Data Problems By Maximum Entropy, pp. 221–232.  Dordrecht: Springer Netherlands, 1990.
  * (26) D. Bergeron and A.-M. S. Tremblay, Algorithms for optimized maximum entropy and diagnostic tools for analytic continuation, Phys. Rev. E 94. 023303 (2016).
  * (27) M. Ohzeki, Sparse modeling for quantum Monte-Carlo simulation, Journal of Physics: Conference Series, 1036. 012020 (2018).
  * (28) K. Yoshimi, J. Otsuki, Y. Motoyama, M. Ohzeki, and H. Shinaoka, “Spm: Sparse modeling tool for analytic continuation of imaginary-time Green’s function,” Computer Physics Communications, 244. 319 (2019).
  * (29) J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Sparse modeling in quantum many-body problems, Journal of the Physical Society of Japan, 89. 012001 (2020).
  * (30) J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Sparse modeling approach to analytical continuation of imaginary-time quantum Monte Carlo data, Phys. Rev. E, 95. 061302 (2017).
  * (31) H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, “Compressing Green’s function using intermediate representation between imaginary-time and real-frequency domains,” Phys. Rev. B 96. 035147 (2017).
  * (32) R. Fournier, L. Wang, O. V. Yazyev, and Q. Wu, Artificial neural network approach to the analytic continuation problem, Phys. Rev. Lett. 124. 056401 (2020).
  * (33) H. Yoon, J.-H. Sim, and M. J. Han, Analytic continuation via domain knowledge free machine learning, Phys. Rev. B 98. 245101 (2018).
  * (34) B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM J. Comput., 24. 227 (1995).
  * (35) S. Boyd and L. Vandenberghe, Convex Optimization.  Cambridge University Press, 2004.
  * (36) I. Daubechies, M. Defrise, and C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications on Pure and Applied Mathematics, 57. 1413 (2004).
  * (37) M. Figueiredo and R. Nowak, An EM algorithm for wavelet-based image restoration, IEEE Transactions on Image Processing, 12. 906 2003\.
  * (38) K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), (2016).
  * (39) K. Gregor and Y. LeCun, Learning fast approximations of sparse coding, in Proceedings of the 27th International Conference on Machine Learning (2010).
  * (40) H. Li, Z. Xu, G. Taylor, C. Studer, and T. Goldstein, “Visualizing the loss landscape of neural nets,” in Advances in Neural Information Processing Systems (S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, eds.), vol. 31, Curran Associates, Inc., 2018.
  * (41) S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in Proceedings of the 32nd International Conference on Machine Learning, (F. Bach and D. Blei, eds.) vol. 37 of Proceedings of Machine Learning Research, (Lille, France), pp. 448–456, PMLR, 2015.
  * (42) J. Ba, J. Kiros, and G. E. Hinton, Layer normalization, arXiv:1607.06450, (2016).
  * (43) Y. Wu and K. He, “Group normalization,” in Computer Vision – ECCV 2018 (V. Ferrari, M. Hebert, C. Sminchisescu, and Y. Weiss, eds.), pp. 3–19, Springer International Publishing, 2018.
  * (44) W. U. Bajwa, R. Calderbank, and S. Jafarpour, “Why gabor frames? Two fundamental measures of coherence and their role in model selection,” Journal of Communications and Networks, 12. 289 (2010).
  * (45) D. G. Mixon, W. U. Bajwa, and R. Calderbank, “Frame coherence and sparse signal processing,” in 2011 IEEE International Symposium on Information Theory Proceedings, (2011)
  * (46) D. L. Donoho and Y. Tsaig, “Fast solution of $\ell_{1}$-norm minimization problems when the solution may be sparse,” IEEE Transactions on Information Theory, 54. 4789 (2008).


# Collatz convergence is a Hydra Game

Alexander Rahn Nuremberg Institute of Technology, Keßlerpl. 12 90489
Nuremberg, Germany  Eldar Sultanow and Idriss J. Aberkane Potsdam University,
Chair of Business Informatics, Processes and Systems, Karl-Marx Straße 67,
14482, Potsdam, Germany
Capgemini, Bahnhofstraße 30, 90402, Nuremberg, GermanyUnesco-Unitwin Complex
Systems Digital Campus, Chair of Prof. Pierre Collet, Pierre Collet, ICUBE -
UMR CNRS 7357, 4 rue Kirschleger, 67000 Strasbourg, France. Corresponding
author<EMAIL_ADDRESS>
(January $23^{rd}$ 2021)

###### Abstract

The Collatz dynamic is known to generate a complex quiver of sequences over
natural numbers which inflation propensity remains so unpredictable it could
be used to generate reliable proof of work algorithms for the cryptocurrency
industry. Here we establish an ad hoc equivalent of modular arithmetic for
Collatz sequences to automatically demonstrate the convergence of infinite
quivers of numbers, based on five arithmetic rules we prove apply on the
entire Collatz dynamic and which we further simulate to gain insight on their
graph geometry and computational properties. We then formally demonstrate
these rules define an automaton that is playing a Hydra game on the graph of
undecided numbers we also prove is embedded in $24\mathbb{N}-7$, proving that
in ZFC the Collatz conjecture is true, before giving a promising direction to
also prove it in Peano arithmetic.

##  1 Introduction

The dynamical system generated by the $3n+1$ problem is known to create
complex quivers over $\mathbb{N}$, one of the most picturesque being the so-
called ”Collatz Feather” or ”Collatz Seaweed”, a name popularized by Clojure
programmer Oliver Caldwell in 2017. The inflation propensity of Collatz orbits
remains so unpredictable it can form the core of a reliable proof-of-work
algorithm for Blockchain solutions [1], with groundbreaking applications to
the field of number-theoretical cryptography as such algorithms are unrelated
to primes yet, being based on the class of congruential graphs, still allow
for a wide diversity of practical variants. If Bocart thus demonstrated that
graph-theoretical approaches to the $3n+1$ problem can be very fertile to
applied mathematics, the authors have also endeavored to demonstrate its pure
number-theoretical interest prior to this work [2], [3], [4], [5].

The definitive purpose of this article however is to establish fundamental
properties of the ”Collatz Feather” and infer provable consequences of those
properties to achieve a positive proof of the Collatz conjecture. Our
methodology consists of using the complete binary tree and the complete
ternary tree111The complete binary tree over odd numbers is defined as
$2\mathbb{N}+1$ endowed with the following two linear applications $\\{\cdot
2-1;\cdot 2+1\\}$. The complete ternary tree over odd numbers is defined as
$2\mathbb{N}+1$ endowed with operations $\\{\cdot 3-2;\cdot 3;\cdot 3+2\\}$
over $2\mathbb{N}+1$ as a general coordinate system for each node of the
Feather. We owe this strategy to earlier discussions with Feferman [6] on his
investigations on the continuum hypothesis, as it is known the complete binary
tree over natural numbers is one way of generating real numbers. The last
author’s discussion with Feferman argued morphisms, sections and origamis of
n-ary trees over $\mathbb{N}$ could be a promising strategy to define objects
of intermediate cardinalities between $\aleph_{n}$ and $\aleph_{n+1}$, in a
manner inspired from Conway’s construction of the surreal numbers [7], which
itself began by investigating the branching factor of the game of Go. The
initial interest therefore, was to investigate the branching factor of the
Collatz Feather and to define the cardinality of the set of its branches. Here
we begin by identifying and proving five arithmetical rules that apply
anywhere on the Collatz dynamic, with the purpose of demonstrating that
applying them from number $1$ allows to generate the entire Feather.

##  2 Five essential rules of Collatz dynamics

######  Note 2.1.

For all intent and purpose we will define Syr(x) or the ”Syracuse action” as
”the next odd number in the forward Collatz orbit of $x$”. Whenever two
numbers $a$ and $b$ have a common number in their orbit, we will also note
$a\equiv b$, a relation that is self-evidently transitive:

$(a\equiv b)\land(b\equiv c)\Rightarrow a\equiv c$

The choice of symbol ”$\equiv$” is a deliberate one to acknowledge a kinship
between our method and modular arithmetic.

######  Definition 2.1.

Actions G, V and S For any natural number $a$ are specified as follows:

  1. 1.

$G(a):=2a-1$

  2. 2.

$S(a):=2a+1$. The rank of $a$ is its number of consecutive end digits $1$ in
base $2$.

  3. 3.

$V(a):=4a+1=G\circ S(a)$

######  Definition 2.2.

Type A, B and C

  1. 1.

a number $a$ is of type A if its base 3 representation ends with digit 2

  2. 2.

a number $b$ is of type B if its base 3 representation ends with digit 0

  3. 3.

a number $c$ is of type C if its base 3 representation ends with digit 1

To remember which is which one need only remember the order of ABC: $a+1$ is
dividable by 3, and so is $c-1$, thus A is on the left of B and C is on the
right.

$1$$3$$5$$7$$15$$13$$11$$9$$17$$19$$21$$23$$25$$27$$29$$31$GGVSVSSGVSSSSVSSSVGGGG

Figure 1: Quiver connecting all odd numbers from 1 to 31 with the arrows of
actions S,V and G. The set $2\mathbb{N}+1$ is thus endowed with three unary
operations without a general inverse that are non commutative with $G\circ
S=V$. Whenever we will mention the inverse of these operations, it will be
assuming they exist on $\mathbb{N}$. Type A numbers are circled in teal, B in
gold and C in purple.

######  Theorem 2.2.

The following arithmetic rules apply anywhere over the system $2\mathbb{N}+1$
endowed with the Collatz dynamic.

  * •

Rule One: $\forall x$ odd, $V(x)\equiv(x)$

  * •

Rule Two: $\forall x,k$ odd, $S^{k}V(x)\equiv S^{k+1}V(x)$ and $\forall x,k$
even, $S^{k}V(x)\equiv S^{k+1}V(x)$

  * •

Rule Three:

$\forall\\{n;y\\}\in\mathbb{N}^{2}$, $\forall x$ odd non B, $3^{n}x\equiv
y\Rightarrow\bigwedge\limits_{i=1}^{n}(V(4^{i}3^{n-i}x))\wedge
S(V(4^{i}3^{n-i}x))\equiv y$

  * •

Rule Four:

$\forall\\{n;y\\}\in\mathbb{N}^{2}$ , $\forall x$ odd non B, $S(3^{n}x)\equiv
y\Rightarrow\bigwedge\limits_{i=1}^{n}(S(4^{i}3^{n-i}x)\wedge
S^{2}(4^{i}3^{n-i}x))\equiv y$

  * •

Rule Five

$\forall n\in\mathbb{N}$, $\forall y\in\mathbb{N}$, $\forall x$ odd non B
where $3^{n}x$ is of rank 1, $a\equiv y$, $a=G(3^{n}x)\\\
\Rightarrow\bigwedge\limits_{i=0}^{n}(S^{i}(G(3^{n-i}x))\wedge
S^{i+1}(G(3^{n-i}x)))\equiv y$

In the following we will prove these five rules.

######  Note 2.3.

In reference to Figure 1 we will call ”vertical even” a number that can be
written $V(e)$ where $e$ is even, and ”vertical odd” if it can be written
$V(o)$ where $o$ is odd. For example, $9$ is the first vertical even number
and $5$ is the first vertical odd.

###  2.1 Proving Rule One

If $a$ is written $4b+1$ then $3a+1=12b+4=4(3b+1)$ therefore $a\equiv b$.

###  2.2 Proving Rule Two

######  Lemma 2.4.

Let $a$ be a number of rank $1$ thus with an odd number $p$ so that $a=G(p)$
then $Syr(S(a))=G(3\cdot p)$. Let $a$ be a number of rank $n$ so that
$S^{-(n-1)}(a)=G(p)$ then $Syr^{n-1}(a)=G(3^{n-1}\cdot p)$

###### Proof.

If $a=2p-1$, $p$ is odd, then it follows:

${\begin{array}[]{l}S(a)=4p-1\\\ \frac{3\cdot
SS(a)+1}{2}=\frac{12p-2}{2}=6p-1=G(3\cdot p)\\\ Syr(S(a))=G(3\cdot
p)\end{array}}$

Let’s generalize to the $n$. If $Syr(S(a))$ can be written $G(3\cdot p)$ it is
also of rank 1, whereas $S(a)$ was of rank 2, therefore, the Syracuse action
has made it lose one rank. All we have to prove now is that $Syr(S^{2}(a))$ =
$S(Syr(S(a)))$ under those conditions:

${\begin{array}[]{l}\frac{3\cdot(S^{2}(a))+1}{2}=6a+5\\\
S(Syr(S(a)))=S(3a+2)=6a+5=Syr(S^{2}(a))\end{array}}$

If a is of rank $n>1$, $Syr(a)$ is of rank $n-1$, and $Syr(S(a))=S(Syr(a))$ ∎

######  Note 2.5.

The $3n+1$ action over an odd number, since it necessary yields an even one,
is in fine equivalent to adding $1$ to it, then the half of the result, then
$-1$. How many times one can add an half to an odd number $+1$ directly
depends on its base 2 representation, and in particular its number of
consecutive end digits 1. Let us take Mersenne numbers for example, which are
defined as $2^{n}-1$. One can transform them consecutively in this way a
number of time equal to their rank-1, indeed, $31$, which is written $11111$
in base 2 is of rank $5$, because $32=2^{5}$ so if one repeats the action ”add
to the number+1 the half of itself” this will yield an even result exactly
four consecutive times. Thus, any strictly ascending Collatz orbit concerns
only numbers $a$ of rank $n>1$, and is defined by

$(a+1)\cdot\left(\frac{3}{2}\right)^{n-1}-1$

######  Lemma 2.6.

Let $a$ be an odd number of rank $1$ that is vertical even, then $3a$ is of
rank 2 or more, and $9a$ is vertical even. Let $a$ be an odd number of rank
$1$ that is vertical odd, then $3a$ is of rank $2$ or more, and $9a$ is
vertical odd.

###### Proof.

If $a$ is vertical even it can be written $8k+1$ $\forall k:3a=24k+3$ and this
number admits an $S^{-1}$ that is $12k+1$, which is an odd number, therefore
$3a$ is at least of rank $2$.

Moreover, $9a=72k+9$ and this number admits a $V^{-1}$ that is $18k+2$, an
even number. Now if $a$ is vertical odd, it can be written $8k+5$ and $\forall
k:3a=24k+15$ and $9a=72k+45$. It follows that $3a$ admits an $S^{-1}$ and $9a$
admits a $V^{-1}$, respectively $12k+7$ and $18k+11$ and they are both odd. ∎

######  Lemma 2.7.

Let $a$ be a number that is vertical even, then $(a)\equiv S(a)$ and
$S^{k}(a)\equiv S^{k+1}(a)$ for any even k. Let $a$ be a number that is
vertical odd, then $S(a)\equiv S^{2}(a)$ and $S^{k}(a)\equiv S^{k+1}(a)$ for
any odd k.

###### Proof.

If $a$ is vertical even then it can be written as $G(p)$ where $p$ is
necessarily vertical (odd or even). We proved that $3p$ is then of rank $2$ or
more and also that we have $Syr(S(a))=G(p)$ so it is necessarily vertical odd
(since $3d$ is of rank $2$ or more) so $Syr(a)=V^{-1}(Syr(S(a))$ and therefore
$a\equiv S(a)$. This behavior we can now generalize to the $n$, because if $a$
is vertical even with $a=G(p)$, then the lemmas we used also provide that
$Syr^{n}(S^{n}(a))=G(3^{n}\cdot p)$ and therefore $Syr^{n}(S^{n}(a))$ will be
vertical even for any even $n$ because $3^{n}\cdot p$ will be vertical
something (even or odd, depending on $p$ only) for any even $n$.

Now if $a$ is vertical odd it can be written $G(p)$ and $p$ is necessarily of
rank $2$ or more because $G\circ S=V$. Thus $3p$ is vertical (even or odd) and
therefore $Syr(S(a))=G(3p)$ is vertical even. ∎

######  Note 2.8.

Observe that in the process of proving Rule Two we also demonstrated that any
number of rank $2$ or more is finitely turned into a rank $1$ number of type A
by the Collatz dynamic, and that any number $x$ of rank $2$ or more so that
$x\equiv S(x)$ under Rule Two is finitely mapped to a type A number that is
vertical even, therefore proving the convergence of such numbers is enough to
prove the Collatz Conjecture.

###  2.3 Proving Rules Three and Four

######  Lemma 2.9.

Let $a$ be a vertical even number with $a=G^{n+2}(S(b))$ where $n$ and $b$ are
odd, then $a\equiv 3^{\frac{n+1}{2}}(b)$. Let $a$ be a vertical even number
with $a=G^{m+2}(S(b))$ where $m$ is even (zero included) and $b$ is odd, then
$a\equiv S(3^{\frac{m}{2}}(b))$

###### Proof.

If $a=G^{n+2}(S(b))$ by definition $a=2^{n+3}b+1$. Then $3\cdot
a+1=3(2^{n+3}b+1)+1)=2^{n+3}\cdot(3b)+4.$ As this expression can be divided by
$2$ no more than twice, we have $Syr(a)=2^{n+1}3b+1=G^{n}(S(3b))$.

Note that if $n=1$ then
$V^{-1}(Syr(a))=V^{-1}(2^{2}\cdot(3a)+1)=2^{2}\cdot\frac{1}{4}\cdot(3b)=3b$
which is of course an odd number. Therefore $Syr(a)$ is vertical odd and
$V^{-1}(Syr(a))=3b$ thus we have proven that $a\equiv 3b$.

If $n=0$ then $a=2^{3}\cdot b+1$ so $3(a+1)=2^{3}\cdot{3b}+4$ therefore
$Syr(a)=S(3b)$ and thus $a\equiv S(3b)$. From this we can generalise the
progression of numbers that can be written $G^{n}(x)$ where $x$ is of rank $2$
or more. Let $b$ be any odd number:

  * •

All ”Variety S” numbers above $b$ are written $V(b\cdot
2^{2k-1})\;\;\text{or}\;\;S(b\cdot 2^{2k})=2^{2k+1}\cdot b+1$ and

  * •

all ”Variety V” numbers above $b$ are written $V(b\cdot 4^{k})$ or
equivalently $S(b\cdot 2^{2k+1})=4^{k+1}\cdot b+1$.

Any number $g$ that can be written $G^{n}(V(x))$ with $x$ odd and $n>0$ may
thus be finitely reduced under the Collatz dynamic to a number that can be
written either $S(3^{m}x)$ or $V(3^{m}x)$ by the repeated following
transformation:

$(g-1)\cdot\biggl{(}\frac{3}{4}\biggr{)}^{k}+1$

Therefore we have indeed that,

  * •

for Variety S numbers: $2^{2k+1}\cdot
b\cdot\left(\frac{3}{4}\right)^{k}+1=2b\cdot 3^{k}+1=S(b\cdot 3^{k})$, which
proves Rule Four.

  * •

for Variety V numbers: $4\cdot 4^{k}\cdot
b\cdot\left(\frac{3}{4}\right)^{k}+1=4b\cdot 3^{k}+1=V(b\cdot 3^{k})$ which
proves Rule Three because Rule One already provides that $V(b\cdot
3^{k})\equiv b\cdot 3^{k}$.

∎

###  2.4 Proving Rule Five

Any type A number of rank 1 can be written $a=G(b)$ where $b$ is of type B. In
proving Rule Two we showed that any number of rank $n>1$ is finitely mapped by
the Collatz dynamics to $G(3^{n-1}\cdot G^{-1}(S^{-(n-1)}(a)))$, which
combined with Rule Two itself gives Rule Five.

$1$$3$$5$$7$$15$$13$$11$$9$$17$$19$$21$$23$$25$$27$$29$$31$$33$$35$$37$$39$$41$$43$$45$$47$$49$$51$$53$$55$$57$$59$$61$$63$$65$$67$$127$$69$$71$$73$$75$$77$$79$$81$$83$$85$$87$$89$$91$$93$$95$$97$$99$$101$$103$$105$$107$$109$$111$$113$$115$$117$$119$$121$$123$$125$$1\equiv
3\Rightarrow 1\equiv 17$$1\equiv 17\Rightarrow 3\equiv 7$$3\equiv 7\Rightarrow
5\equiv 9$$1\equiv 9$$1\equiv 9$$3\equiv 9$$9\equiv 19$$3\equiv 7$$11\equiv
33$$45\equiv 89$$19\equiv 39$

Figure 2: Just a few applications of Rules Three, Four and Five starting from
$1\equiv 3\equiv 5$ are here plotted in gold. Rules One and Two are plotted in
black. Whenever a number is connected to $1$ by a finite path of black and/or
gold edges it is proven to converge to $1$.

##  3 The Golden Automaton

######  Definition 3.1.

On {$2\mathbb{N}+1$; G, S} where Rules One and Two are considered pre-computed
(the black edges on Figure 2) the systematic computation of Rules Three, Four
and Five from number $1$ onward is called the ”Golden Automaton”.

###  3.1 ”Golden Arithmetic”

Our purpose is to develop an ad hoc unary algebra that could found a
congruence arithmetic specifically made to prove the Collatz conjecture, and
which we intend as an epistemological extension of modular arithmetic, hence
our use of the symbol $\equiv$ in this article rather than the usual
$\thicksim$ which is seen more frequently in the Collatz-related literature.
This ”Golden arithmetic” involves words taken in the alphabet $\\{G;S;V;3\\}$,
which we will call in their order of application, just like in turtle
graphics. For example VGS3 means $3\cdot S\circ G\circ V$

Rules 3, 4 and 5 may now be reformulated as such, without loss of generality
as long as Rules One and Two are still assumed:

  * •

Rule Three: Let b be of type B, then $b\equiv VGS3^{-1}$ from b. We will call
this action $R_{b}(x)=16\frac{x}{3}+1$

  * •

Rule Four: Let c be of type C, then $c\equiv GS3^{-1}$ from c. We will call
this action $R_{c}(x)=\frac{4x-1}{3}$

  * •

Rule Five: Let a be of type A, then $a\equiv G3^{-1}$ from a. We will call
this action $R_{a}(x)=\frac{2x-1}{3}$

As Rules One and Two ensure that the quiver generated by the Golden Automaton
is branching, with each type B number that is vertical even providing both a
new A type and a new B type number to keep applying respectively rules 5 and
3, we may follow only the pathway of type A numbers to define a single non-
branching series of arrows, forming a single infinite branch of the quiver.
The latter, if computed from number 15, leads straight to 31 and 27, solving a
great deal of other numbers on the way:

$\begin{array}[]{rll}15&\equiv 81&\hskip 20.00003pt\text{Rule 3}\\\ 81&\equiv
1025&\hskip 20.00003pt\text{First type A reached by {Rule 3}}\\\ 1025&\equiv
303&\hskip 20.00003pt\textbf{Rule 5}\\\ 303&\equiv 607&\hskip
20.00003pt\text{Rule 2}\\\ 607&\equiv 809&\hskip 20.00003pt\text{Rule 4}\\\
809&\equiv 159&\hskip 20.00003pt\textbf{Rule 5}\\\ 159&\equiv 319&\hskip
20.00003pt\text{Rule 2}\\\ 319&\equiv 425&\hskip 20.00003pt\text{Rule 4}\\\
425&\equiv 283&\hskip 20.00003pt\textbf{Rule 5}\\\ 283&\equiv 377&\hskip
20.00003pt\text{Rule 4}\\\ 377&\equiv 111&\hskip 20.00003pt\textbf{Rule 5}\\\
111&\equiv 593&\hskip 20.00003pt\text{Rule 3}\\\ 593&\equiv 175&\hskip
20.00003pt\textbf{Rule 5}\\\ 175&\equiv 233&\hskip 20.00003pt\text{Rule 4}\\\
233&\equiv 103&\hskip 20.00003pt\textbf{Rule 5}\\\ 103&\equiv 137&\hskip
20.00003pt\text{Rule 4}\\\ 137&\equiv 91&\hskip 20.00003pt\textbf{Rule 5}\\\
91&\equiv 161&\hskip 20.00003pt\text{Rule 4}\\\ 161&\equiv 31&\hskip
20.00003pt\textbf{Rule 5}\\\ 31&\equiv 41&\hskip 20.00003pt\text{Rule 4}\\\
41&\equiv 27&\hskip 20.00003pt\textbf{Rule 5}\end{array}$

Again, it is in no way a problem, but rather a powerful property of the Golden
Automaton that this particular quiver branch already cover 19 steps (and
actually more) because each of them is branching into other solutions.

We may follow another interesting sequence to show that in the same way that
Mersenne number 15 finitely solves Mersenne number 31, Mersenne number 7
solves Mersenne number 127, this time we will follow a B branch up to
$Syr^{6}(127)$ which we know can be written $G(3^{6})$ because 127 is the
Mersenne of rank 7.

By Rule 4 we have the first equivalence $\boldsymbol{7\equiv 9}$ and
$\boldsymbol{9\equiv 25\equiv 49}$.

So by Rule 2 we also have $\boldsymbol{25\equiv 51}$.

Rule 3 gives $\boldsymbol{51\equiv 273}$ and again $\boldsymbol{273\equiv
1457=G(729)\equiv 127}$.

The cases of $15$ proving the convergence of $31$ and $27$ and of $7$ proving
the one of $127$ naturally lead us to the following conjecture:

######  Conjecture 3.1.

Suppose all odd numbers up to $2^{n}$ are proven to converge to $1$ under the
Collatz dynamic, then the Golden Automaton finitely proves the convergence of
those up to $2^{n+1}$

And indeed we already have that the Golden Automaton starting with $1$ proves
$3$ by Rule One, then $3$ proves all numbers from $5$ to $15$ which in turn
prove all numbers from $33$ to $127$. In the next subsection we render larger
quivers generated by the Golden Automaton to provide a better understanding of
their geometry and fundamental properties, and to demonstrate why it is so,
and more generally, why it can be proven they can reach any number in
$2\mathbb{N}+1$ in ZFC.

###  3.2 Computational Scale-up

The purpose of this subsection is to identify provable fundamental properties
of the Golden Automaton by scaling it up on the full binary tree over
$2\mathbb{N}+1$. To streamline its algorithmic scaling, we use the simplified
rules we defined in the previous subsection, again, without loss of
generality. Our precise purpose is to pave the way for a formal demonstration
that proving the convergence of odd numbers up to $n$ is always isomorphic to
a Hydra Game. In the next figures we color all the elements of
$24\mathbb{N}-7$ as for example $\\{17,41,65,\ldots\\}$ in red to as we
demonstrate in the next section they precisely from the ”heads” in the Hydra
Game.

Figure 3: Golden Automaton confined to numbers smaller than 32 Figure 4:
Golden Automaton confined to numbers smaller than 64 Figure 5: Golden
Automaton confined to numbers smaller than 128 Figure 6: Golden Automaton
confined to numbers smaller than 256

##  4 Proving the Collatz conjecture

In this section we prove that the Collatz conjecture holds in the
Zermelo–Fraenkel set theory with the axiom of choice included (abbreviated
ZFC). And we prove that Collatz conjecture also holds in the Peano Arithmetic.

###  4.1 ZFC proves the Collatz conjecture

######  Theorem 4.1.

In ZFC, the Collatz conjecture is true.

######  Definition 4.1.

A hydra is a rooted tree with arbitrary many and arbitrary long finite
branches. Leaf nodes are called heads. A head is short if the immediate parent
of the head is the root; long if it is neither short nor the root. The object
of the Hydra game is to cut down the hydra to its root. At each step, one can
cut off one of the heads, after which the hydra grows new heads according to
the following rules:

  * •

if the head was long, grow $n$ copies of the subtree of its parent node minus
the cut head, rooted in the grandparent node.

  * •

if the head was short, grow nothing

######  Lemma 4.2.

The Golden Automaton reaching any natural number is a Hydra Game over a finite
subtree of the complete binary tree over $24\mathbb{N}-7$.

###### Proof.

The essential questions to answer in demonstrating either a homomorphism
between a Hydra game and the Golden Automaton reaching any odd number, or that
the Golden Automaton is playing at worst a Hydra Game are:

  * •

What are the Hydra’s heads?

  * •

How do they grow?

  * •

Does the Golden Automaton cut them according to the rules (at worst)?

######  Definition 4.2.

A type A number that is vertical even is called an $A_{g}$. The set of $A_{g}$
numbers is $24\mathbb{N}-7$. Type B numbers that verify $b\equiv S(b)$ and
type C numbers that verify $c\equiv S(c)$ under Rule Two are called Bups and
Cups respectively.

What are the Hydra’s heads? $A_{g}$ numbers are the heads of the Hydra. They
are 12 points apart on $2\mathbb{N}+1$ (24 in nominal value, e.g. 17 to 41)
and any Bup or Cup of $rank>1$ they represent under Rule Five is smaller than
them since action $R_{a}$ is strictly decreasing so up to the $n^{th}$ $A_{g}$
there are $2n$ (Bups + Cups) of rank 2 or more and half of them are equivalent
to these $A_{g}$ (e.g. between 17 and 41 Bup 27 is equivalent to $A_{g}$ 41,
which is equivalent to Cup 31 by Rule Four

How do they grow? Between any two consecutive $A_{g}$ in $2\mathbb{N}+1$ there
are

  * •

8 non-A numbers

  * •

1 of them at most is mapped to the second $A_{g}$

  * •

3 at most are ”ups” (Bup or Cup) of rank 2 or more

Besides, we also have anywhere that:

  * •

Let b be of type B, there are $\frac{2b}{3}$ numbers of type $A_{g}$ that are
smaller than $V^{2}(b)$

  * •

Let c be of type C, there are $\frac{S(c)}{3}$ numbers of type $A_{g}$ that
are smaller than $V^{2}(c)$

  * •

Let 3c be a type B where c is of type C, there are $\frac{S(c)}{3}$ numbers of
type $A_{g}$ up to $R_{b}(3c)$ included

  * •

Let 3a be a type B where a is of type A, there are $\frac{G(a)}{3}$ numbers of
type $A_{g}$ smaller than $R_{b}(3a)$

Which is defining the growth of the heads. Indeed, any supposedly diverging
$A_{g}$ is forming a Hydra, as we have proven $24\mathbb{N}-7$ contains an
image of all undecided Collatz numbers and that any non-decreasing trajectory
identifies a subtree within this set.

Does the Golden Automaton play a Hydra game?

It could be demonstrated that the Golden Automaton is playing an even simpler
game as it is branching and thus cutting heads several at a time and in
particular cutting some long heads without them doubling 222In fact, the
reason the Golden Automaton dominates $24\mathbb{N}-7$ so fast is that it is
playing a much simpler game one could call ”Hecatonchire v. Hydra game” namely
a Hydra game where Herakles’ number of arms is also multiplying at each step
but as this is needless for the final proof we can now simply demonstrate that
even under the worst possible assumptions it follows at least the rules of a
regular Hydra game.

The computing of $15\equiv\ldots\equiv 27$ that we detailed in Subsection 3.1
is one case of the playing of a Hydra Game by the Golden Automaton; we
underlined each use of Rule 5 specifically so the reader can now report to it
more easily, because each time this rule is used, a head ($A_{g}$) has just
been cut.

The demonstration that $27$ and $31$ converge is the cutting of heads $41$ and
$161$ respectively. This single branch of the Automaton having first cut head
$17$, reaches to the head $1025$ via B-typed numbers $15$ and $81$. It is
therefore playing a Hydra with $\frac{1025+7}{24}=43$ heads of which one
($17$) is already cut at this point and of which at least $8$ are rooted (so
cutting them does not multiply any number of heads). This process being
independent of the targeted number, we now have that the reaching of any
number by the Golden Automaton is at least equivalent to the playing of a
Hydra with $n$ heads of which $0<m<n$ are rooted. Even without demonstrating
more precise limit theorems for each factors $n$ and $m$ (which could still be
a fascinating endeavor) the road is now open for a final resolution of the
Collatz conjecture. ∎

From there indeed, we know with Goodstein [8] and Kirby and Paris [9] that
assuming a system strong enough to prove Peano arithmetic is consistent and
that $\epsilon_{0}$ is well-ordered, no Hydra game can be lost. Since we have
that the reaching of any number $n$ is a Hydra Game for the Golden Automaton,
we have that the Golden Automaton cannot fail to finitely reach any natural
number.

###  4.2 Can Peano Arithmetic also prove the Collatz conjecture?

If it is now sure that any system strong enough to prove the convergence of
Goodstein’s series also proves the Collatz conjecture, it could very well be
possible to prove it from Peano arithmetic alone. In this final subsection, we
intend to outline a strategy towards such a demonstration by defining a
different game than the Hydra one and in particular, a zero-player game that
is significantly simpler than John Conway’s Game of Life and played on the
complete binary tree $\\{2\mathbb{N}+1;G,S\\}$.

In this cellular automaton, each cell is identified by a unique odd number and
can only adopt three states:

  * •

Black, meaning the odd number is not (yet) proven to converge under the
iterated Collatz transformation or equivalently that it is only equivalent to
another black number

  * •

Gold, meaning the odd number is proven to converge and the consequences of its
convergence have not yet been computed, ie. it can have an offspring

  * •

Blue, meaning the number is proven to converge and the consequences of its
convergence have been computed ie. its offspring has already been turned gold

In this ad hoc yet simpler game of life each gold cell yields and offspring
then turns blue, and whenever a cell is blue or gold its odd number is proven
to converge. Starting with one cell colored in gold at the positions $1$, it
applies the following algorithm to each gold cell in the natural order of odd
numbers:

  1. 1.

Rule 1: if a cell on $x$ is gold color cell on $V(x)$ gold

  2. 2.

Rule 2: if a cell on $x$ is gold, color cell on $S(x)$ gold depending on the
precise conditions of rule 2

  3. 3.

If a cell on $a$ of type A is gold, then color that on $R_{a}(x)$ in gold

  4. 4.

If a cell on $c$ of type C is gold, then color that on $R_{c}(x)$ in gold

  5. 5.

After applying the previous rules on a gold cell, turn it blue
Note that applying $R_{b}$ on a type B number being equivalent to Rule 1 then
$R_{c}$ the algorithm needs not implement a defined $R_{b}$

Whenever a complete series of odd numbers between $2^{n}+1$ and $2^{n+1}-1$
has been colored in gold, it ticks it and returns what we will call its
computational ”expense”, namely all the numbers colored blue and gold that are
higher than $2^{n+1}-1$, thus giving a clear measurement of the algorithmic
time it takes the Golden Automaton to prove the convergence of each complete
level of the binary tree over $2\mathbb{N}+1$. We then plot the evolution of
this expense on a linear and a logarithmic scale.

Figure 7: case n=6, illustrating the principles of the game we defined. On the
middle image, row {5;7} has been solved , with an ”expense” of 8 numbers also
solved above it. On the right image, row {9; 11; 13; 15} with an expense of 6.
As number $1$ is the neutral element of operation $R_{c}$ we leave it in gold
during all the simulation  Figure 8: case 12, seventh row completed Figure 9:
case 12, eighth row completed

To facilitate the observation of each row of the binary tree being covered by
the Golden Automaton we here zoom into each of them individually:

Figure 10: zoom of row 11 (going from 1025 to 2047: each line has about 100
dots)

Figure 11: Row 10 (513 to 1023)

Figure 12: zoom Row 9 (257 to 511) Figure 13: Row 8 (129 to 255)

Figure 14: Charting expenses to proof the convergence of each row of the
binary tree over $2\mathbb{N}+1$. For example, the amount of gold and blue
dots above the row going from 2049 to 4095 (row 12) is slightly below 175000.
The total expense emerges from the generated tree’s height. The same plot
against a logarithmic scale (Right) indicates a line-like shape.

From there we can thus provide two strategies to finalise a proof of the
Collatz conjecture without the axiom of choice (which is needed to demonstrate
no Hydra Game can be lost) and precisely within Peano arithmetic. The first
strategy would consist of using automated theorem proving to single out the
linear behavior we expose in Figure 14 as a provable property of our game. The
second - and we believe the most promising - would consist of analysing the
average reproductive rate of gold dots, and demonstrating at any level $n$ of
the binary tree they cannot fail to finitely take over any population of black
dots below it.

######  Definition 4.3.

(Reproductive Rate)

  1. 1.

The reproductive rate of any golden dot on coordinate $n$ is the number of
golden and blue dots it generates that are at or below coordinate $V(n)$.

  2. 2.

The average reproductive rate of all black dots converges to 3,5 new black
dots generated from $x$ under or equal to $V(x)$ . Indeed for any $x$, in the
time to reach $V(x)$ only the offspring below $G(x)$ gets to generate new
black dots under the rules of the binary tree, $S(x)$ can only reproduce once
by applying $G(S(x))$, thus generating $V(x)$, and all numbers between $S(x)$
and the next Mersenne number cannot reproduce. More precisely, for any odd
number $x$, there are $x+G^{-1}(x)$ odd numbers between itself and $V(x)$
included.

We can now count the average or limit reproductive rate of each type of golden
dot.

  * •

for type B numbers we always have V(b) in the offspring, and also S(b) one out
of two times, so the average reproductive rate is 1,5

  * •

for type C numbers we always have V(c) and R(c) in the offspring, and S(c) one
out of two times so the average reproductive rate is 2,5

  * •

for type A numbers we have a more interesting converging series defining the
average number of gold dots generated below V(a)+2. To easily count the
offspring we use the property of the golden Automaton that it can only output
$A_{g}$ through operations $R_{b}$ or $R_{c}$, and can only output non-Aup
through operations V and S.

The net offspring below $V(A_{g})+2$ of $A_{g}$ numbers is defined by a
branching process:

    * –

$R_{a}(A_{g})$ outputs either an Aup (in proportion of $\frac{1}{2}$) or
equally either a Bup or a Cup ($\frac{1}{4}$ each)

    * –

then again $R_{a}(Aup)$ outputs either an Aup ($\frac{1}{2}$) or equally
either a Bup or a Cup (again, $\frac{1}{4}$ each)

    * –

how long $R_{a}(A_{g})$ keeps rendering an Aup only depends on $n$ where
$A_{g}=G(3^{n}x)$ with $x$ non B (per Rule 5)

So the formula of the average offspring below V(a)+2 of all type $A_{g}$
number, which can be written $G(3^{n}x)$ with x non B is $3n+2,5$. Since one
out of two $A_{g}$ can be written $3^{1}x$, a quarter can be written $3^{2}x$
and so on we have the general formula that $A_{g}$ numbers up to $G(3^{n})$
would have an average reproductive rate, which limit we can now determine:

$\lim_{n\to\infty}\sum_{i=1}^{n}\frac{1}{2^{i}}\cdot(3{i}+2,5)=6,5$

And this formula does not even account for the accumulated offspring of all
the C and B numbers also colored gold in the process, and still below
$V(A_{g})$ so the average net reproductive rate of $A_{g}$ numbers is
converging to a strictly greater value than this limit.

As we know that the reproductive rate of black dots below $V(x)+2$ converges
to $3,5$, when the average reproductive rate of all A, B and C type numbers
generated by the Golden Automaton grows beyond $3,5$ we can be certain it can
always finitely finish any row. As we already have $\frac{1,5+2,5+6,5}{3}=3,5$
since we did not count the offspring of C and B type numbers solved by $A_{g}$
numbers and still below $V(A_{g})$ in the computing of their birthrate, we can
prove the average birthrate of golden numbers tends to equal $3,5+\epsilon$
with $\epsilon>0$, which finishes the Peano-arithmetical proof.

##  5 Dedication, Attribution and Acknowledgements

This work was supported by a personal grant to I. Aberkane from Mohammed VI
University, Morocco and by a collaboration between Capgemini, Potsdam
University, The Georg Simon Ohm University of Applied Sciences, and Strasbourg
University. I. Aberkane created the framework of studying the Collatz dynamic
in the coordinate system defined by the intersection of the binary and ternary
trees over $2\mathbb{N}+1$, identified and demonstrated the five rules and
predicted they would be isomorphic to a Hydra game over the set of undecided
Collatz numbers, which he defined as well, allowing for a final demonstration
of the Collatz Conjecture; he also outlined and computed the strategy of using
reproductive rates of dots to define a Peano-arithmetical proof. Contributing
equally, E. Sultanow and A. Rahn designed and coded an optimised, highly
scalable graphical implementation of the five rules and ran all the
simulations, confirming the Hydra game isomorphism and computing the first
ever dot plot of the Golden Automaton over odd numbers, which they optimised
as well. They were also the first team to ever simulate the five rules to the
level achieved in this article, and to confirm their emerging geometric
properties on such a scale, including the linearity of their logarithmic
scaling and the limit reproductive rates of single dots of the golden
automaton.

I. Aberkane wishes to thank the late Prof. Solomon Feferman and Prof. Alan T.
Waterman Jr, along with Prof. Paul Bourgine, Prof. Yves Burnod, Prof. Pierre
Collet, Dr. Françoise Cerquetti and Dr. Oleksandra Desiateryk.

The authors dedicate this work to the memory of John Horton Conway
(1937-2020), Solomon Feferman (1928-2016) and Alan T. Waterman Jr (1918-2008).

## References

  * [1] F. Bocart.  Inflation propensity of collatz orbits: a new proof-of-work for blockchain applications.  Journal of Risk and Financial Management, 11(4):83, 2018.
  * [2] I. Aberkane.  On the syracuse conjecture over the binary tree.  HAL, 2017.
  * [3] I. Aberkane.  Endomorphisms of the collatz quiver.  HAL, 2020.
  * [4] C. Koch, E. Sultanow, and S. Cox.  Divisions by two in collatz sequences: A data science approach.  International Journal of Pure Mathematical Sciences, 21, 2020.
  * [5] E. Sultanow, C. Koch, and S. Cox.  Collatz sequences in the light of graph theory.  Technical report, Potsdam University, 2020.
  * [6] S. Feferman.  Is the continuum hypothesis a definite mathematical problem., 2012.
  * [7] D. E. Knuth.  Surreal Numbers: How Two Ex-students Turned on to Pure Mathematics and Found Total Happiness.  Addison-Wesley, Boston, MA, 12 edition, 1974.
  * [8] R. Goodstein.  On the restricted ordinal theorem.  Journal of Symbolic Logic, 9(2):33–41, 1944.
  * [9] L. Kirby and J. Paris.  Accessible independence results for peano arithmetic.  Bulletin of the London Mathematical Society, 1982.