# The R-mAtrIx Net

Shailesh Lal, Suvajit Majumder, Evgeny Sobko

► To cite this version:

Shailesh Lal, Suvajit Majumder, Evgeny Sobko. The R-mAtrIx Net. Machine Learning: Science and Technology, 2024, 5 (3), pp.035003. 10.1088/2632-2153/ad56f9. hal-04092128

HAL Id: hal-04092128

https://hal.science/hal-04092128v1

Submitted on 10 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The R-mAtrIx Net

To cite this article: Shailesh Lal et al 2024 *Mach. Learn.: Sci. Technol.* **5** 035003

View the [article online](#) for updates and enhancements.

## You may also like

- [Benchmark calculations of electron impact electronic excitation of the hydrogen molecule](#)

T Meltzer, J Tennyson, Z Mašín et al.

- [Electron collisions with molecular nitrogen in its ground and electronically excited states using the R-matrix method](#)

He Su, Xinlu Cheng, Hong Zhang et al.

- [Electron collisions with molecular hydrogen from electronically excited states using the R-matrix method](#)

T Meltzer and J Tennyson
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

6 November 2023

REVISED

14 April 2024

ACCEPTED FOR PUBLICATION

11 June 2024

PUBLISHED

4 July 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# The R-mAtrlx Net

Shailesh Lal¹*, Suvajit Majumder² and Evgeny Sobko³

¹ Yanqi Lake Beijing Institute of Mathematical Sciences and Applications (BIMSA), Yanqi Island, Huairou District, Beijing 101408, People's Republic of China

² Centre for Mathematical Science, City, University of London, Northampton Square, EC1V 0HB London, United Kingdom

³ Laboratoire de Physique de l'École Normale Supérieure, PSL University, CNRS, 24 rue Lhomond, 75005 Paris, France

* Author to whom any correspondence should be addressed.

E-mail: shaileshlal@bimsa.cn

Keywords: Yang-Baxter equation, R-matrix, neural network, integrable spin chain

## Abstract

We provide a novel neural network architecture that can: i) output R-matrix for a given quantum integrable spin chain, ii) search for an integrable Hamiltonian and the corresponding R-matrix under assumptions of certain symmetries or other restrictions, iii) explore the space of Hamiltonians around already learned models and reconstruct the family of integrable spin chains which they belong to. The neural network training is done by minimizing loss functions encoding Yang-Baxter equation, regularity and other model-specific restrictions such as hermiticity. Holomorphy is implemented via the choice of activation functions. We demonstrate the work of our neural network on the spin chains of difference form with two-dimensional local space. In particular, we reconstruct the R-matrices for all 14 classes. We also demonstrate its utility as an *Explorer*, scanning a certain subspace of Hamiltonians and identifying integrable classes after clusterisation. The last strategy can be used in future to carve out the map of integrable spin chains with higher dimensional local space and in more general settings where no analytical methods are available.
## Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>1.1. Summary of neural network and training:</td><td>3</td></tr><tr><td>2. An overview of spin chains and quantum integrability</td><td>4</td></tr><tr><td>2.1. R-matrices of the difference form for spin chains with two-dimensional local Hilbert space</td><td>6</td></tr><tr><td>3. Neural networks for the R-matrix</td><td>6</td></tr><tr><td>3.1. An overview of neural networks</td><td>6</td></tr><tr><td>3.2. Machine learning the R-matrix</td><td>8</td></tr><tr><td>4. Results</td><td>10</td></tr><tr><td>4.1. Specific integrable spin chains</td><td>10</td></tr><tr><td>4.1.1. Hermitian models: XYZ spin chain and its isotropic limits</td><td>10</td></tr><tr><td>4.1.2. Two-dimensional classification</td><td>11</td></tr><tr><td>4.2. Explorer: new from existing</td><td>15</td></tr><tr><td>5. Conclusions and future directions</td><td>18</td></tr><tr><td>Data availability statement</td><td>20</td></tr><tr><td>Acknowledgments</td><td>20</td></tr><tr><td>Appendix A. Classes of 2D integrable spin chains of difference form</td><td>20</td></tr><tr><td>Appendix B. Designing the neural network</td><td>28</td></tr><tr><td>B.1. The neural network architecture and forward propagation</td><td>30</td></tr><tr><td>B.2. Comparing different activation functions</td><td>30</td></tr><tr><td>B.3. Proof of concept: training with a single hidden neuron</td><td>31</td></tr><tr><td>References</td><td>32</td></tr></table>

## 1. Introduction

Neural networks and deep learning have recently emerged as a competitive computational tool in many areas of theoretical physics and mathematics, in addition to their several impressive achievements in computer vision and natural language processing [1]. In string theory and algebraic geometry for instance, the application of these methods was initiated in [2–5]. Since then, deep learning has seen several interesting and remarkable applications in the field, both on the computational front [6–15] as well as towards the explication of foundational questions [16]. They have appeared in the context of conformal field theory, critical phenomena, spin systems and matrix models [17–26]. More generally, deep learning has found interesting applications in mathematics, ranging from the solution of partial nonlinear differential equations [27, 28], to symbolic calculations [29] and to hypothesis generation [30, 31]. Interestingly, deep learning is starting to play an increasingly important role in symbolic regression, i.e. the extraction of exact analytical expressions from numerical data [32]. While it is difficult to pin-point any one solitary reason for this confluence of several fields into deep learning, there are some important themes that do seem to play a recurring role. Firstly, deep neural networks are a highly flexible parametrized class of functions and provide us an efficient way to approximate various functional spaces and scan over them [33–36]. The same neural network, as we shall shortly see, can learn a Jacobi elliptic function as easily as it does a trigonometric or an exponential function. Such approximations train well for a variety of loss landscapes, including non-convex ones. Secondly, over the previous many years, robust frameworks for the design and optimization of neural networks have been developed, both as an explication of best practices [37–42] and the development of standardized software for implementation [43–45]. This has made it possible to reliably train increasingly deeper networks which are optimized to carry out increasingly sophisticated tasks such as the direct computation of Ricci flat metrics on Calabi Yau manifolds [12–15] and the solution of differential equations without necessarily providing the neural network data obtained from explicitly sampling the solution. Further, in recent interesting developments, deep learning has been applied to analyze various aspects of symmetry in physical systems ranging from their classification to their automated detection [19, 46–50].

The profound role played by symmetry in theoretical physics and mathematics is hard to overstate. Probably its most compelling expression in theoretical physics is found in the bootstrap program which rests on the idea that a theory may be significantly or even fully constrained just by the use of general principles and symmetries without analysis of the microscopic dynamics. For example, the S-matrix bootstrap bounds the space of allowed S-matrices relying only on unitarity, causality, crossing, analyticity and global symmetries [51–53]. This provides rigorous numerical bounds on the coupling constants and significantly restricts the space of self-consistent theories [54–56]. This line of consideration finds its ultimate realization in two dimensions once applied to integrable theories. Integrable bootstrap complements the aforementioned constraints with one extra functional Yang-Baxter (YB) equation, manifesting the scattering
factorization, which allows us to fix the S-matrix completely [57]. The same YB equation appears in the closely related context of integrable spin chains [58, 59]. Now instead of S-matrix, it restricts the R-matrix operator whose existence allows one to construct a commuting tower of higher charges and prove integrability. Practically, one has to solve the functional YB equation in a certain functional space. There is no known general method to do so, and all existing approaches are limited in the scope of application and fall into three groups. The first class of methods is algebraic in nature, exploiting the symmetry of R-matrix [60–66]. The second approach aims to directly solve the functional equation or the related differential equation [67]. The third alternative utilizes the boost operator to generate higher charges and impose their commutativity [68–70].

In this paper we shall demonstrate how neural networks and deep learning provide an efficient way to numerically solve the YB equation for integrable quantum spin chains. On an immediate front, we are motivated by recent interesting work on classical integrable systems using machine learning [46–48, 71]. The approach taken in the work [71] of learning classical Lax pairs for integrable systems by minimization of the loss functions encoding a flatness condition has a particularly close parallel to our approach. However, to the best of our knowledge, the present work is the first attempt to apply machine learning to quantum integrability, the analysis of R-matrices and the YB equation.

Our analysis utilizes neural networks to construct an approximator for the R-matrix and thereby solve functional YB equation while also allowing for the imposition of additional constraints. We look into the sub-class of all possible R-matrices, namely those that are regular and holomorphic, and incorporate the YB equation into the loss function. Upon training for the given integrable Hamiltonian, we successfully learn the corresponding R-matrix to a prescribed precision. Using spin chains with two-dimensional space as a main playground we reproduce all R-matrices of difference form which was recently classified in [68]. Moreover, this Solver can be turned into an Explorer which scans the space (or a certain subspace) of all Hamiltonians looking for integrable models, which in principle allows us to discover new integrable models inaccessible to other methods. Below we provide the summary of the neural network and its training, as well as an overview of the paper.

**1.1. Summary of neural network and training:**

The functional YB equation, see equation (2.3) below, is holomorphic in the spectral parameter $u \in \mathbb{C}$ and as such, holds over the entire complex plane. In this paper, we shall restrict our training to the interval $\Omega = (-1, 1)$ on the real line, but design our neural network so that it analytically continues to a holomorphic function over the complex plane. Each entry into the R-matrix is separately modeled by multi-layer perceptrons (MLPs) with two hidden layers of 50 neurons each, taking as input parameters the variable $u \in \Omega$. More details are available in section 3.2 and appendix B. All the neurons are *swish* activated [72], except for the output neurons which are *linear* activated. Training proceeds by optimizing the loss functions that encode the YB equations (3.7), regularity (3.12), and constraints on the form of the spin chain Hamiltonian, for instance via (3.13). Hermiticity of the Hamiltonian, if applicable, is imposed by the loss (3.15). Optimization is done using *Adam* [73] with a starting learning rate of $\eta = 10^{-3}$ which is annealed $\eta = 10^{-8}$ in steps of $10^{-1}$ by monitoring the YB loss (3.7) on validation data for saturation. *Adam's* hyperparameters $\beta_1$ and $\beta_2$ are fixed to 0.9 and 0.999 respectively. In the following, we will refer to this learning rate policy as the *standard schedule*. We apply this framework to explore the space of R-matrices using the following strategies:

1. **Exploration by attraction:** The Hamiltonian loss (3.13) is imposed by specifying target numerical values for the two-particle Hamiltonian, or some ansatz/symmetries instead (like 6-vertex, 8-vertex, etc). We also formally include here the ultimate case of *general search* when no restrictions are imposed on the Hamiltonian at all. This strategy is predominantly used in our Section 4.1.

2. **Exploration by repulsion:** We can generate new solutions by *repelling away* from an ansatz or a given spin chain. This requires us to activate the loss function (3.17) for a few epochs in order to move from the specific Hamiltonian. This strategy is employed in section 4.2.

Further, we also have two schemes for initializing training.

1. **Random initialization:** We randomly initialize the weights of the neural network using He initialization [74]. This samples the weights from either a uniform or a normal distribution centered around 0 but with a variance that scales as the inverse power of the layer width.

2. **Warm-start:** we use the weights and biases for an already learnt solution.
A brief overview of this paper is as follows. In section 2 we quickly introduce the R-matrix and other key concepts from the quantum integrability of spin chains relevant to this paper. Particularly in section 2.1, we review the classification program of 2D spin chains of difference form through the boost automorphism method [68]. Section 3 contains a review of neural networks with a view towards machine learning the R-matrix given an ansatz for the two-particle Hamiltonian. Our methodology for this computation is provided in section 3.2. We then present our results in section 4. Section 4.1 focuses on hermitian XYZ and XXZ models (section 4.1.1), and prototype examples from the 14 gauge-inequivalent classes of models in [68](section 4.1.2). The latter sub-section also contrasts training behaviour for integrable and non-integrable models. Section 4.2 presents a preliminary search strategy for new models which we illustrate within a toy-model setting: rediscovering the two integrable subclasses of 6-vertex Hamiltonians. Section 5 discusses ongoing and future research directions.

## 2. An overview of spin chains and quantum integrability

Quantum integrability, like its classical counterpart, hinges on the presence of tower of conserved charges in involution, i.e. operators that mutually commute. In this paper, we will consider quantum integrable spin chains and the goal of this section is to introduce such systems, and provide a brief overview of the R-matrix construction in their context.

The Hilbert space of the spin chain is a L-fold tensor product $V = V_1 \otimes \dots \otimes V_L$ of d-dimensional vector spaces $V_i \sim V = \mathbb{C}^d$. The Hamiltonian $H$ of a spin chain with nearest-neighbour interaction is a sum of two-site Hamiltonians $H_{i,i+1}$:

$$H = \sum_{i=1}^{L} H_{i,i+1}, \qquad (2.1)$$

where we assume periodic boundary conditions: $H_{L,L+1} \equiv H_{L,1}$. The chrestomathic example of the integrable spin-chain is spin-1/2 XYZ model:

$$H = \sum_{i=1}^{L} \sum_{\alpha} J^{\alpha} S_{i}^{\alpha} S_{i+1}^{\alpha}, \qquad (2.2)$$

where $\alpha = \{x, y, z\}$ and $S_i^\alpha$ are Pauli matrices acting in the two-dimensional space $V_i = \mathbb{C}^2$ of i-th site. In particular case when $J^x = J^y$ it reproduces XXZ model, while in the case of three equal coupling constants $J^x = J^y = J^z = J$ the Hamiltonian reduces to the XXX spin chain. These famous magnet models are just a few examples of integrable spin chains and now we turn to the general construction.

The central element for the whole construction and proof of quantum integrability is the R-matrix operator $R_{ij}(u)$ which acts in the tensor product $V_i \otimes V_j$ of two spin sites<fnref n="4" /> and satisfies the YB equation:

$$R_{ij}(u-v) R_{ik}(u) R_{jk}(v) = R_{jk}(v) R_{ik}(u) R_{ij}(u-v) \qquad (2.3)$$

where the operators on the left and right sides act in the tensor product $V_i \otimes V_j \otimes V_k$. The R-matrix is assumed to be an analytic function of the spectral parameter $u$. Further, in order to guarantee locality of the interaction in (2.1), it must reduce to the permutation operator $P_{ij}$ when evaluated at $u=0$, i.e.

$$R_{ij}(0) = P_{ij}. \qquad (2.4)$$

This condition will be referred to as *regularity* in the following sections. In principle, one can consider more general case when $R_{ij}(0)$ is a constant solution of YBE. However, such spin chains with long-range interactions are beyond the scope of this paper.

We next turn to defining the monodromy matrix $T_a(u)$. This matrix, denoted by $T_a(u) \in \text{End}(V_a \otimes \prod_{i=1}^L \otimes V_i) \times \mathbb{C}$, acts on the spin chain plus an auxiliary spin site labeled by $a$ with Hilbert space as $V_a \sim \mathbb{C}^d$. It is defined as a product of R-matrices $R_{a,i}(u)$ acting on the auxiliary site and one of the spin chain sites and is given by

$$T_a(u) = R_{a,L}(u) R_{a,L-1}(u) \dots R_{a,1}(u). \qquad (2.5)$$

<footnote n="4">⁴ In general, the R-matrix $R(u, v)$ is an analytic function of two complex arguments $u, v$, which can be viewed as momenta of two particles at the two sites. Here and throughout the paper we shall exclusively focus our analysis to a restricted class of R-matrices of *difference form* $R(u, v) = R(u - v)$ depending only on a single complex argument $w = u - v$. Many R-matrices of non-difference form $R(u, v) \neq R(u - v)$ are already known, and particularly appear in the context of AdS/CFT correspondence, but are beyond the scope of the current paper.</footnote>
The transfer matrix $T(u) \in \text{End}(\prod_{i=1}^{L} \otimes V_i) \times \mathbb{C}$ is obtained by taking a trace over the auxiliary vector space $V_a$:

$$T(u) = \text{tr}_a(\mathcal{T}_a(u)). \qquad (2.6)$$

From the YB equation one can derive the following $RTT$ relation constraining monodromy matrix entries

$$R_{12}(u-v) \mathcal{T}_1(u) \mathcal{T}_2(v) = \mathcal{T}_2(v) \mathcal{T}_1(u) R_{12}(u-v). \qquad (2.7)$$

This condition can be used to prove that the transfer matrices commute at different values of the momenta

$$[T(u), T(v)] = 0. \qquad (2.8)$$

The above condition implies that the transfer matrix $T(u)$ encodes all the commuting charges $\mathcal{Q}_i$ as series-expansion in $u$:

$$\log T(u) = \sum_{n=0}^{\infty} \mathcal{Q}_{n+1} \frac{u^n}{n!}. \qquad (2.9)$$

Hence we have<fnref n="5" />

$$\mathcal{Q}_{n+1} = \frac{d^n}{du^n} \log T(u)|_{u=0} = \frac{d^{n-1}}{du^{n-1}} \left( T^{-1}(u) \frac{d}{du} T(u) \right) \Big|_{u=0}. \qquad (2.10)$$

The Hamiltonian density $H_{i,i+1}$ introduced earlier in equation (2.1) can be generated from the R-matrix using

$$H_{i,i+1} = R_{i,i+1}^{-1}(0) \frac{d}{du} R_{i,i+1}(u)|_{u=0} = P_{i,i+1} \frac{d}{du} R_{i,i+1}(u)|_{u=0} \qquad (2.11)$$

where $P_{i,i+1}$ is the permutation operator between sites $i, i+1$. Also, we emphasize that while the charges are conventionally computed in equation (2.10) at $u=0$, this computation can equally well be done at generic values of $u$ to extract mutually commuting charges. The only difference is we no longer recover the Hamiltonian directly as one of the commuting charges.

YB equation (2.3) should be supplemented with certain analytical properties of R-matrix. For example, as was already mentioned, we assume that the R-matrix is a holomorphic function of spectral parameter $u$ and equal to the permutation matrix at $u=0$ (2.4). Furthermore, one can impose extra physical constraints like braided unitarity

$$R_{12}(u) R_{21}(-u) = g(u) \mathbf{I}, \quad g(u) = g(-u), \qquad (2.12)$$

crossing symmetry<fnref n="6" />, and possibly additional global symmetries. We shall also impose restrictions on the form of the resulting Hamiltonian. These restrictions may follow from requirements such as hermiticity and from symmetries of the spin chain. In addition, given a solution for the YB equation, one can generate a whole family of solutions by acting with the following transformations:

1. Similarity transformation: $(\Omega \otimes \Omega)R(u)(\Omega^{-1} \otimes \Omega^{-1})$ where $\Omega \in \text{Aut}(V)$ is a basis transformation. It transforms the commuting charges as $\mathcal{Q}_n \rightarrow (\otimes^L \Omega)\mathcal{Q}_n(\otimes^L \Omega^{-1})$

2. Rescaling<fnref n="7" /> of the spectral parameter: $u \rightarrow cu, \forall c \in \mathbb{C}$. This leads to a scaling in the charges as $\mathcal{Q}_n \rightarrow c^{n-1}\mathcal{Q}_n$

3. Multiplication by any scalar holomorphic function $f(u)$ preserving regularity condition: $R(u) \rightarrow f(u)R(u), f(0) = 1$. This degree of freedom can be used to set one of the entries of R-matrix to one or any other fixed function.

4. Permutation, transposition and their composition: $PR(u)P, R(u)^T, PR^T(u)P$. They transform the commuting charges as well. The Hamiltonian $\mathcal{H}$ is transformed to $PHP, PH^T P, \mathcal{H}^T$ respectively.

<footnote n="5">⁵ In practice, numerically it is more stable to work with the second formula on the right hand side than the first.</footnote>

<footnote n="6">⁶ The explicit form of the crossing symmetry varies for the different classes of models.</footnote>

<footnote n="7">⁷ For the general r-matrix of non-difference form $R(u, v)$ there is a reparametrization freedom $u \rightarrow f(u)$, however for the difference form $R(u)$ it reduces just to rescaling.</footnote>
In general, one should always be careful of these redundancies when comparing a trained solution against analytic results. Following [68], we shall fix the above symmetries when presenting our results in section 4.1.2 and appendix A. We look at gauge-equivalent solutions as well, by introducing similarity transformations in 4.1.2.

## 2.1. R-matrices of the difference form for spin chains with two-dimensional local Hilbert space

We will illustrate the work of our neural network using spin chains with two-dimensional local space as a playground. The regular difference-form integrable models in this context have recently been classified using the Boost operator in [68]. Here, we present a brief overview of the methods and results of this paper. Boost automorphism method allows one to find integrable Hamiltonians by reducing the problem to a set of algebraic equations. Let us focus on a spin chains with two-dimensional space $V = \mathbb{C}^2$ and nearest-neighbour Hamiltonian (2.1). One formally defines the boost operator $\mathcal{B}$ [75] as

$$ \mathcal{B} = \sum_{a=-\infty}^{\infty} a H_{a,a+1}, \qquad (2.13) $$

which generates higher charges $\mathbb{Q}_n, n \ge 3$, from the Hamiltonian $\mathbb{Q}_2$ via action by commutation:

$$ \mathbb{Q}_{r+1} = [\mathcal{B}, \mathbb{Q}_r]. \qquad (2.14) $$

This was used in [68] to successfully classify all 2-dimensional integrable Hamiltonians by solving the system of algebraic equations arising from imposing vanishing conditions on commutators between $\mathbb{Q}_i$, upto some finite value of $i$. Surprisingly it turns out that for the considered models, the vanishing of the first non-trivial commutator $[\mathbb{Q}_2, \mathbb{Q}_3] = 0$ is a sufficient condition to ensure the vanishing of all other commutators. Then making an ansatz for the R-matrices and solving YB equation in the small $u$ limit, the authors constructed the corresponding R-matrices and confirmed the integrability of the discovered Hamiltonians. The solutions can be organized into two classes: XYZ-type, and non-XYZ type, distinguished by the non-zero entries appearing in the Hamiltonian.

$$ H_{\text{XYZ type}} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_2 & b_2 & 0 \\ d_2 & 0 & 0 & a_2 \end{pmatrix}, \qquad H_{\text{non-XYZ type}} = \begin{pmatrix} a_1 & a_2 & a_3 & a_4 \\ 0 & b_1 & b_3 & b_3 \\ 0 & c_1 & c_2 & c_3 \\ 0 & 0 & 0 & d_1 \end{pmatrix}. \qquad (2.15) $$

Generically, all non-zero entries would be complex valued. Hermiticity, for the actual XYZ model and its XXZ and XXX limits, places additional constraints. Integrability also imposes additional algebraic constraints between the non-zero entries, none of which involve complex conjugation in contrast to hermiticity. Amongst the XYZ type models, there are 8 distinct solutions each corresponding to some set of algebraic constraints among the matrix elements of $H$. In particular, there is one 4-vertex model ($H_{4v}$) which is purely diagonal (see (A.1)). Next, there are two 6-vertex models ($H_{6v,1}$, $H_{6v,2}$) where $d_1$ and $d_2$ are constrained to vanish, among other conditions (see (A.2) and (A.4)). One of these, namely $H_{6v,1}$, is a non-hermitian generalisation of the XXZ model. There are also two 7-vertex models ($H_{7v,1}$, $H_{7v,2}$) where only $d_2$ vanishes (see (A.6) and (A.8)), and three 8-vertex models ($H_{8v,1}$, $H_{8v,2}$, $H_{8v,3}$) where all entries are non-zero(see (A.10), (A.13) and (A.17)). Here, $H_{8v,1}$ is the non-hermitian generalization of the XYZ model. Among these classes the Hamiltonians are distinguished by additional algebraic constraints on the non-zero elements which we have enumerated in appendix A. The corresponding R-matrices for these models were obtained in [67]. The non-XYZ models are similarly divided into 6-classes with Hamiltonians $H_{\text{class-1}}...H_{\text{class-6}}$ which have been explicitly enumerated in equation (A.19). Among these, the class 1 and class 2 Hamiltonians have rank less than four. For convenience, we also explicitly write down all of these R-matrices, both for the XYZ type and non-XYZ type models, in appendix A.

# 3. Neural networks for the R-matrix

This section reviews several essential facts about neural networks before presenting our own network-design for deep-learning and the associated custom loss functions. Further details regarding the network architecture and training schedule can be found in appendix B.

## 3.1. An overview of neural networks

The central computation in this paper is the utilization of neural networks to construct R-matrices that correspond to given integrable spin chain Hamiltonians. We therefore furnish a lightning overview of neural
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The schematic for a dense neural network, also known as a fully connected neural network. The four-dimensional input $(a_1^{(0)}, \dots, a_4^{(0)})$ is fed via the input layer (green) to a series of three fully connected layers (purple) containing four neurons each and finally feeds into the output layer of three neurons (orange). Every neuron in a given layer receives inputs from all neurons in the preceding layer, and in turn, its output is passed as input to all neurons in the next layer.</figcaption></figure>

networks in this section, along with the details of our implementation of the neural network solver for YB equations. We will focus on dense neural networks, also known as MLPs, schematically displayed in figure 1. These networks consist of an input layer $a^{\text{in}} \in \mathbb{R}^{n_0}$, followed by a series of fully connected layers and terminate in an output layer $a^{\text{out}} \in \mathbb{R}^{n_{l+1}}$. Data is read in to the network at the input layer and the output is collected at the output layer. There are $L$ fully connected layers in this network, where the $\ell$th layer contains $n_\ell$ neurons. Each neuron $a_m^{(\ell)}$ in a $l$th fully connected layer receives inputs from all the neurons in the previous $(l-1)$th layer and the output of the neuron is in turn fed as an input to neurons in the succeeding layer:

$$ \begin{pmatrix} a_1^{(\ell)} \\ a_2^{(\ell)} \\ \vdots \\ a_{n_\ell}^{(\ell)} \end{pmatrix} = h \left[ \begin{pmatrix} w_{1,0}^{(\ell)} & w_{1,1}^{(\ell)} & \cdots & w_{1,n_{\ell-1}}^{(\ell)} \\ w_{2,0}^{(\ell)} & w_{2,1}^{(\ell)} & \cdots & w_{2,n_{\ell-1}}^{(\ell)} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n_\ell,0}^{(\ell)} & w_{n_\ell,1}^{(\ell)} & \cdots & w_{n_\ell,n_{\ell-1}}^{(\ell)} \end{pmatrix} \begin{pmatrix} a_1^{(\ell-1)} \\ a_2^{(\ell-1)} \\ \vdots \\ a_{n_\ell-1}^{(\ell-1)} \end{pmatrix} + \begin{pmatrix} b_1^{(\ell)} \\ b_2^{(\ell)} \\ \vdots \\ b_{n_\ell}^{(\ell)} \end{pmatrix} \right], \quad \ell = 1, \dots, L+1, \quad (3.1) $$

where $w^{(l)} \in \mathcal{M}(n_l, n_{l-1}, \mathbb{R})$ is a weight matrix, $b^{(l)} \in \mathbb{R}^{n_l}$ - bias vector, $h(z)$ is in general a non-linear, non-polynomial function known as the activation function acting component-wise :

$$ h \begin{pmatrix} z_1 \\ z_2 \\ \vdots \\ z_n \end{pmatrix} = \begin{pmatrix} h(z_1) \\ h(z_2) \\ \vdots \\ h(z_n) \end{pmatrix}. \quad (3.2) $$

In (3.1) we also identify $a^{(0)} = a^{\text{in}}$ and $a^{(L+1)} = a^{\text{out}}$ with input and output layers respectively. Introducing shorthand notation for the affine transformations in equation (3.1) as $A^{(\ell)}(a^{(\ell-1)}) \equiv w^{(\ell)}a^{(\ell-1)} + b^{(\ell)}$, the neural network $a^{\text{out}}(a^{\text{in}}): \mathbb{R}^{n_0} \to \mathbb{R}^{n_{l+1}}$ can be expressed as compositions of affine transformations, and activation functions:

$$ a^{\text{out}} = h \circ A^{L+1} \circ h \circ A^L \circ \dots \circ h \circ A^1 \circ a^{\text{in}}. \quad (3.3) $$

The output function of the neural network is tuned by tuning the weights and biases. It is by now well established that such neural networks are a highly expressive framework capable of approximation of extremely complex functions, and indeed there exists a series of mathematical proofs which attest to their universal approximation property, e.g. [33–35, 76, 77]. This property, along with the feature learning capability of deep neural networks is the key driver to the automated search for R-matrices which we have implemented here.
Neural networks learn a target function $f(x)$ of the input data $x$ by optimizing a cost function $\mathcal{L}$ which provides a measure of the discrepancy between the actual and desired properties of the function $f$. The parameters $w$ and $b$ are then tuned to minimize this discrepancy. The best known and the most canonical examples of this are the supervised learning problems where the neural network is supplied with data $\mathcal{D} = \{(x,y)\}$ consisting of pairs of input vectors $x$ with their expected output values $y$. The neural network then tunes its weights and biases to minimize the cost function. Having done so, the output function $f$ thus learned by the neural network obeys

$$f(x) \approx y \quad \forall \quad (x,y) \in \mathcal{D}, \qquad (3.4)$$

while allowing for the possibility of outliers. A popular class of loss functions are

$$\mathcal{L}(\{w, b\}) = \sum_{(x,y) \in \mathcal{D}} |y - y_{\text{pred}}(x)|^q, \quad y_{\text{pred}}(x) = f(x; w, b), \qquad (3.5)$$

where $q=1$ corresponds to the mean average error and $q=2$ to the mean square error, respectively. We will shortly see that in contrast to the above classic supervised learning set-up, our loss functions impose constraints on the neural network output functions rather than train on a dataset of input/output values for the functions $\mathcal{R}(u)$ directly sample $\mathcal{R}(u)$ at various values of $u$ for training.

## 3.2. Machine learning the R-matrix

We are now ready to describe our proposed methodology for constructing R-matrices $\mathcal{R}(u)$ by optimizing a neural network using appropriate loss functions. An R-matrix has elements $R_{ij}(u)$ at least some of which are non-zero. In the following, we shall focus solely on the $R_{ij}(u)$ which are not identically zero as functions of $u$. We also restrict the training to the real values of spectral parameter $u \in \Omega = (-1, 1)$ and exclusively use holomorphic activations function in order to guarantee the holomorphy of the resulting R-matrix $\mathcal{R}(u)$. The matrix elements $R_{ij}(u)$ of this R-matrix are modeled by neural networks as

$$\mathcal{R}_{ij}(u) = a_{ij}(u) + i b_{ij}(u) : \quad \begin{array}{ccc} a_{ij;1}^{(1)} & & a_{ij;1}^{(2)} \\ & \ddots & \\ a_{ij;50}^{(1)} & & a_{ij;50}^{(2)} \end{array} ; \quad r_{ij} = \{a_{ij}, b_{ij}\}. \qquad (3.6)$$

We have decomposed the matrix element $\mathcal{R}_{ij}(u)$ into $a_{ij}(u) + i b_{ij}(u)$ in order to learn complex-valued functions $\mathcal{R}_{ij}$ while training with real MLPs on the real interval $\Omega$. In this paper, purely for uniformity, we have modeled each such $a_{ij}(u)$ and $b_{ij}(u)$ using an MLP containing two hidden layers of 50 neurons each and one linear activated output neuron. We emphasize that the identification of $a_{ij}(u)$ and $b_{ij}(u)$ to real and imaginary parts of $\mathcal{R}_{ij}(u)$ is only valid over the real line, and these functions separately continue into holomorphic functions over the complex plane whose sum $\mathcal{R}_{ij}(u)$ is holomorphic by construction. Now, $\mathcal{R}_{ij}(u)$ is required to solve the the YB equation (2.3) subject to (2.4). We may also place constraints on the corresponding two-particle $H$ given by (2.11). These criteria are encoded into loss functions which the R-matrix $\mathcal{R}_{ij}(u)$ aims to minimize by training. For example, in order to train $\mathcal{R}_{ij}(u)$ to satisfy YB equation (2.3) for all values of spectral parameter $u$ from the set $\Omega \subset \mathbb{C}$ we introduce the following loss function:

$$\mathcal{L}_{\text{YBE}} = \|\mathcal{R}_{12}(u-v)\mathcal{R}_{13}(u)\mathcal{R}_{23}(v) - \mathcal{R}_{23}(v)\mathcal{R}_{13}(u)\mathcal{R}_{12}(u-v)\|, \qquad (3.7)$$

where $||...||$ is a matrix norm defined as $||A|| = \sum_{\alpha, \beta=1}^{n} |A_{\alpha\beta}|$ for an complex-valued $n \times n$ matrix $A$. During the forward propagation we sample a mini-batch of $u$ and $v$ values, from which the corresponding $u-v$ is constructed. Along this paper, the spectral parameters $u$ and $v$ run over the discrete set of 20 000
points randomly chosen from the interval $\Omega^8$. The loss function $\mathcal{L}_{YBE}$ is positive semi-definite, vanishing only when $\mathcal{R}(u)$ solves the YB equation. In principle, one may imagine a scan across the space of all functions in which case, solutions of the YB equation would minimize the loss (3.7) to zero.

In practice of course, one cannot scan across the space of all functions and is restricted to a hypothesis class. Here the hypothesis class is implicitly defined by the design of the neural network, the choice of numbers of layers, number of neurons in each layer, as well as the activation function. Varying the weights and biases of the neural network allows us to scan across this hypothesis space. While in general the exact R-matrix may not belong to this hypothesis class, and the loss function would then be strictly positive, deep learning may allow us to approach the desired functions $\mathcal{R} \to R$ to a high degree of accuracy. In summary, if we restrict to a hypothesis class which does not include an actual solution of the YB equation then

$$ \mathcal{L}_{YBE} \ge \epsilon = \min_{\{w', b'\}} \mathcal{L}_{YBE}(\{w', b'\}) > 0, \quad (3.8) $$

where ideally $\epsilon$ would be small, indicating that we have obtained a good approximation to the true solution. We expect that scanning across wider and wider hypothesis classes would bring $\epsilon$ closer and closer to zero. Further, while the RTT equation (2.7) follows from the YB equation (2.3), it can also be imposed separately as a loss function on the network in order to improve the training:

$$ \mathcal{L}_{\text{RTT}} = \|\mathcal{R}_{12}(u-v)\mathcal{T}_1(u)\mathcal{T}_2(v) - \mathcal{T}_2(v)\mathcal{T}_1(u)\mathcal{R}_{12}(u-v)\|. \quad (3.9) $$

Next, we have constraints that must be imposed on the R-matrix at $u=0$. Following equation (2.4) and equation (2.11) in previous section, we require that<fnref n="9" />

$$ R(0) = P, \quad P \frac{d}{du} R(u) |_{u=0} = H, \quad (3.11) $$

where $H$ is the two particle Hamiltonian. They both can be encoded in the loss function as

$$ \mathcal{L}_{\text{reg}} = \|\mathcal{R}(0) - P\|, \quad (3.12) $$

$$ \mathcal{L}_H = \| P \frac{d}{du} \mathcal{R}(u) |_{u=0} - H \|. \quad (3.13) $$

Here, we should mention that we have some flexibility in the manner in which we implement the Hamiltonian constraint $\mathcal{L}_H$. Firstly, one can fix the exact numerical values for the entries of $H$ and learn corresponding R-matrix. We will also consider extensions of this loss function where we supply only algebraic constraints restricting the search space for target Hamiltonians to those with certain symmetries or belonging to certain gauge-equivalence classes. In general, such Hamiltonian constraints give us the requisite control to converge to the different classes of integrable Hamiltonians, and we will name such regime as a *exploration by attraction*.

In the same spirit, when working with the XYZ spin chain or its XXZ and XXX limits, we also require that the two-particle Hamiltonian computed from $R(u)$ is hermitian, i.e.

$$ H = H^{\dagger}. \quad (3.14) $$

We impose this condition by means of the loss function

$$ \mathcal{L}_{\dagger} = \|H - H^{\dagger}\|. \quad (3.15) $$

We shall therefore train our neural network with the loss function

$$ \mathcal{L} = \mathcal{L}_{\text{YBE}} + \mathcal{L}_{\text{reg}} + \lambda_{\text{RTT}} \mathcal{L}_{\text{RTT}} + \lambda_H \mathcal{L}_H + \lambda_{\dagger} \mathcal{L}_{\dagger}, \quad (3.16) $$

$$ R(u) = P + \tilde{R}(u), \quad (3.10) $$

where $\tilde{R}(u)$ is the target function of the neural network, which we design to identically output $\tilde{R}(0) = 0$. While this is possible in principle, in practice it turns out that since the neural network is learning a function in the vicinity of $P$, which trivially minimizes the YB equation and all other constraints imposed, it almost invariably collapses to the trivial solution and learns $\tilde{R}(u) = 0$ across all $u$. It would nonetheless be interesting to identify such architectures that successfully learn non-trivial R-matrices and this is in progress.

<footnote n="8">⁸ The number of points used during the training bounds the precision which one can reach and in our case it will be of order $10^{-4}$.</footnote>

<footnote n="9">⁹ It is tempting to consider a variation of our method which involves *residual learning a la* the ResNet family of networks [74]. As opposed to learning deviations from identity, which is typically the approach adopted in the ResNet architecture, we may define</footnote>
where putting the coefficients $\lambda_\alpha$, for $\alpha = \{RTT, H, \dagger\}$, to zero removes the corresponding loss term from being trained.

The loss function (3.16) produces a very complicated landscape and the NN should approach its minimum during the training. Usually, this search is performed with gradient based optimization methods. One might be skeptical about being stuck in some local minimum instead of finding the global minimum of such complicated loss function in a very high dimensional hypothesis space. However, recent analysis revealed that deep NNs end up having all their local minima almost at the same value as the global one [78, 79]. In other words, there are many different configurations of weights and biases resulting into a function of similar accuracy as the one corresponding to the global minimum. There are also many saddle points and some of them have big plateau and just a small fraction of descendent directions, making them practically indistinguishable from the local minima. However, most of their losses are close to the global minimum as well. Those with significantly higher losses have a bigger number of descendent directions and thus can be escaped during the learning [78, 79].

We find that the training converges to yield simultaneously low values for each of the above losses as applicable. Further, while the hyper-parameters {$\lambda$} are tunable experimentally, setting them all to 1 is a useful default. However, for fine-tuning the training it is also useful to tune these parameters to reflect the specific task at hand. We provide the requisite details in section 4 where we discuss specific training methodologies and the corresponding results. We will also discuss there a new loss function

$$ \mathcal{L}_{\text{repulsion}} = \exp(-||H - H_o||/\sigma), \qquad (3.17) $$

which is useful to fine-tune the training to access new integrable Hamiltonians $H$ in the neighbourhood of previously known integrable Hamiltonians $H_o$, we will call such regime as a *exploration by repulsion*.

As a final observation on the choice of activation functions, we note that at the level of the discussion above, any holomorphic activation function such as sigmoid, tanh, and the sinh would suffice. In practice we find that the training converges faster and more precisely using the *swish* activation [72]. This is given by

$$ \text{swish}(z) = z\sigma(z), \quad \sigma(z) = \frac{1}{1 + e^{-z}}. \qquad (3.18) $$

We have provided some comparison tests across activation functions in appendix B.2.

# 4. Results

We present our results for learning R-matrices within the restricted setting of spin chains of difference form with two dimensional local space. Our analysis will be divided into three parts. First, we will learn hermitian XYZ model and its well-known XXZ and XXX limits, comparing our deep-learning results against the analytic plots. Then we remove hermiticity and reproduce all 14 classes of solutions from [68]. The last set of experiments demonstrates how our neural network in the explorer mode can search for integrable models exploring the space of Hamiltonians.

## 4.1. Specific integrable spin chains

In this sub-section we look at specific physical models, by imposing tailored conditions on the Hamiltonian derived from the training R-matrix. This includes constraints on the Hamiltonian entries at $u = 0$, and hermiticity of the Hamiltonian.

### 4.1.1. Hermitian models: XYZ spin chain and its isotropic limits

Imposing hermiticity on the 8-vertex Hamiltonian, we learn the classic XYZ integrable spin chain and its symmetric XXZ limit. We start with the following 8-vertex model ansatz for the R-matrix

$$ R(u) = \begin{pmatrix} a & 0 & 0 & d \\ 0 & b & c & 0 \\ 0 & c & b & 0 \\ d & 0 & 0 & a \end{pmatrix} \qquad (4.1) $$

and impose the loss functions for YBE, hamiltonian constraint, regularity, and hermiticity (see equation (3.16)). The target Hamiltonians comprise of a 2-parameter family $H_{XYZ}(J_x, J_y, J_z)$ given by

$$ H_{XYZ}(J_x, J_y, J_z) = J_x S_1^x S_2^x + J_y S_1^y S_2^y + J_z S_1^z S_2^z = \begin{pmatrix} J_z & 0 & 0 & J_x - J_y \\ 0 & -J_z & 2 & 0 \\ 0 & 2 & -J_z & 0 \\ J_x - J_y & 0 & 0 & J_z \end{pmatrix}, \qquad (4.2) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. The evolution of training losses for the XYZ model, shown on the log scale. The losses tend to fall in a step-wise manner, which corresponds approximately to the learning rate schedule the network is trained with.</figcaption></figure>

where we have set $J_x + J_y$ to be equal to 2. The symmetric limit of XXZ model is realised for $J_x = J_y = 1$. A useful reparametrisation for these models is in terms of $(\eta, m)$ [80]

$$J_x = 1 + \frac{\sqrt{m} \operatorname{sn}(2\eta|m)}{2}, \quad J_y = 1 - \frac{\sqrt{m} \operatorname{sn}(2\eta|m)}{2}, \quad J_z = \operatorname{cn}(2\eta|m) \operatorname{dn}(2\eta|m). \quad (4.3)$$

The analytic solution for the XYZ R-matrix is given in terms of Jacobi elliptic functions as

$$\begin{align}
a(u) &= \frac{\operatorname{sn}(2\eta + \omega u|m)}{\operatorname{sn}(2\eta|m)} \exp\left(-\frac{\operatorname{cn}(2\eta|m)\operatorname{dn}(2\eta|m)}{2\operatorname{sn}(2\eta|m)}\omega u\right), \nonumber \\
b(u) &= \frac{\operatorname{sn}(\omega u|m)}{\operatorname{sn}(2\eta|m)} \exp\left(-\frac{\operatorname{cn}(2\eta|m)\operatorname{dn}(2\eta|m)}{2\operatorname{sn}(2\eta|m)}\omega u\right), \nonumber \\
c(u) &= \exp\left(-\frac{\operatorname{cn}(2\eta|m)\operatorname{dn}(2\eta|m)}{2\operatorname{sn}(2\eta|m)}\omega u\right), \nonumber \\
d(u) &= \sqrt{m}\operatorname{sn}(\omega u|m)\operatorname{sn}(2\eta + \omega u|m) \exp\left(-\frac{\operatorname{cn}(2\eta|m)\operatorname{dn}(2\eta|m)}{2\operatorname{sn}(2\eta|m)}\omega u\right), \tag{4.4}
\end{align}$$

where $\omega = 2\operatorname{sn}(2\eta|m)$, and $m$ is the elliptic modular parameter$^{10}$. Our model consistently learns the R-matrices for the XYZ model for generic values of the free parameters $\eta, m$. Figure 2 gives the time evolution of the different loss terms during training. Figure 3 plots the R-matrix component ratios with respect to $R_{12}$ in terms of the spectral parameter, and compares them with the corresponding analytic functions for a generic choice of deformation parameters $\eta = \frac{\pi}{3}$ and $m = 0.6$. Letting $m = 0$, we recover the XXZ models for generic values of $\eta$.

### 4.1.2. Two-dimensional classification

Here, we lift the hermiticity constraint on the Hamiltonian, thus allowing for more generic integrable models. As we shall see below, the neural network successfully learns all the 14 classes [68] of difference-form integrable (not necessarily Hermitian) spin chain models with two-dimensional space at each site. The R-matrices corresponding to each of these classes are written down explicitly in appendix A. Towards the end of this sub-section, we also present results for learning solutions in generic gauge obtained by similarity transformation of integrable Hamiltonians from the aforementioned 14 classes. We shall discuss the results in two parts: XYZ type models, and non-XYZ type models.

<footnote>$^{10}$ Usually, these expressions are written in terms of the elliptic modulus $k$ instead of the modular parameter $m = k^2$, e.g. as in [67]. We have expressed them in terms of the modular parameter following the implementation in both Python and Mathematica.</footnote>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3. (a) XYZ R-matrix as ratios with respect to the (12) component for $\eta = \pi/3$, $m = 0.6$, (b) Absolute and relative errors for R-matrix.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4. (a) Predicted vs actual R-matrix component ratio w.r.t. $R_{00}$ for XXZ-type model with $a_1 = a_2 = 0.3$, $b_1 = 0.45$, $b_2 = 0.6$, $c_1 = 0.4$, $c_2 = 0.25$, (b) Absolute error between predicted and actual R-matrix component ratios.</figcaption></figure>

The first set of Hamiltonians under consideration are generalisations of the XYZ model (discussed in the previous sub-section), with at most 8 non-zero elements in its Hamiltonian density

$$H_{8-vertex} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_2 & b_2 & 0 \\ d_2 & 0 & 0 & a_2 \end{pmatrix} \qquad (4.5)$$

where the coefficients can take generic complex values. The XYZ model corresponds to the subset with $a_1 = a_2, b_1 = b_2, c_1 = c_2, d_1 = d_2$. As discussed in section 2.1, these models can be further sub-divided into four, six, seven and eight vertex models. On the other hand, there are 6 distinct classes of non-XYZ type solutions. Here we will discuss the training results for one example each from the XYZ and non-XYZ type models, since the training behaviour is similar within these two types. Rest of the models will be presented in appendix A. Figure 4 plots the R-matrix components as ratios with respect to $R_{00}$ for a generic 6-vertex
<figure><img src="image_10.png" /><figcaption>Figure 5. (a) Predicted vs actual R-matrix components for the Hamiltonian of class-1, with coefficients $a_1 = .5, a_2 = .3, a_3 = .9, a_4 = 1.5, a_5 = .4$. Note here the R-matrices are automatically normalised since $R_{00}$ was fixed to the constant value of 1, (b) Absolute error between predicted and actual R-matrix components.</figcaption></figure>

model with $d_1 = d_2 = 0$, and $a_1 = a_2$. The figure also includes the absolute and relative errors with respect to the corresponding analytic R-matrix (see equation (A.2)).

From the non-XYZ classes, we will focus on the following 5-vertex Hamiltonians

$$H_{\text{class-1}} = \begin{pmatrix} 0 & a_1 & a_2 & 0 \\ 0 & a_5 & 0 & a_3 \\ 0 & 0 & -a_5 & a_4 \\ 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (4.6)$$

For integrability, we require the additional condition

$$a_1 a_3 = a_2 a_4. \qquad (4.7)$$

Training the Hamiltonian constraint (3.13) for generic values $a_1 = 0.5, a_2 = 0.3, a_3 = 0.9, a_4 = 1.5, a_5 = 0.4$ satisfying the above integrability condition, we get over 99.9% accuracy for training over ~100 epochs. Figure 5 plots the trained R-matrix components and absolute errors with respect to the analytic R-matrices in equation (A.20), for the above choice of target Hamiltonian.

We have also surveyed more general solutions beyond the representative solutions of the 14 classes a la [68], by changing the gauge of the R-matrix as well as the corresponding Hamiltonian. As noted earlier in section 2, we can act with a $2 \times 2$ similarity matrix $\Omega$ on the R-matrix:

$$R(u) \rightarrow R^{\Omega}(u) = (\Omega \otimes \Omega) R(u) (\Omega^{-1} \times \Omega^{-1}) \qquad (4.8)$$

$$H \rightarrow H^{\Omega} = (\Omega \otimes \Omega) H (\Omega^{-1} \times \Omega^{-1}). \qquad (4.9)$$

If $R(u)$ satisfies YB equation, so does $R^{\Omega}(u)$. A generic similarity matrix $\Omega$

$$\Omega = \begin{pmatrix} v_{11} & v_{12} \\ v_{21} & v_{22} \end{pmatrix} \qquad (4.10)$$

with non-zero off-diagonal entries $v_{12}, v_{21} \neq 0$, results in conjugated R-matrices and Hamiltonians with all 16 non-zero entries. We trained 16-vertex Hamiltonians resulting from XYZ model in the general gauge and recovered the corresponding R-matrix with a relative error of order $O(0.1\%)$. Generic XYZ type models, as well as non-XYZ type models gave similar results for different gauges. Figure 6 plots the learnt R-matrix components for XXZ model with $\eta = \pi/3$ conjugated by the matrix $\Omega = \begin{pmatrix} 0.4 & 0.5 \\ -1.2 & 1 \end{pmatrix}$. For comparison with analytic formulae, we normalised our results by taking ratios with respect to a fixed component $R_{00}$, i.e. we
<figure><img src="image_11.png" /><figcaption>Figure 6. (a) Predicted R-matrix component ratios w.r.t. R$_{00}$, for conjugated XXZ model with η = π/ 3 and similarity matrix</figcaption></figure>

Ω = ( 0 . 4 − 1 . 2 0 . 5 1 ) , (b) Absolute error between predicted and actual R-matrix ratios.

plot R$_{ij}$ R$_{00}$ . As a result of starting from the XXZ model, the R-matrix R$_{XXZt}$ in the general gauge has following highly symmetric form

$$
R_{XXZg} = \begin{pmatrix}
R_{00} & R_{01} & R_{01} & R_{03} \\
R_{10} & R_{11} & R_{12} & -R_{01} \\
R_{10} & R_{12} & R_{11} & -R_{01} \\
R_{30} & -R_{10} & -R_{10} & R_{00}
\end{pmatrix} . \qquad (4.11)
$$

Thus we only plot the entries R$_{00}$ , R$_{01}$ , R$_{03}$ , R$_{10}$ , R$_{11}$ , R$_{12}$ , R$_{30}$ . Since there exists overall normalisation ambiguity, we should only compare ratio of R-matrix entries with the analytic solution written in the same gauge.

Next we discuss the difference in the training of integrable vs non-integrable models with our neural network. We will focus on two representative examples : 6-vertex model with Hamiltonian H$_{6}$$_{v}$$_{,}$$_{1}$ , and class 4 models with Hamiltonian H$_{class-4}$ . Similar results hold across all the 14 classes.

For 6-vertex models with Hamiltonians following equation (4.5) with d$_{1}$ = d$_{2}$ = 0, generic values of the coefficients a$_{i}$ , b$_{i}$ , c$_{i}$ , d$_{i}$ for i = 1 , 2 leads to non-integrable models, unless

$$
a_1 = a_2, \quad a_1 + a_2 = b_1 + b_2. \tag{4.12}
$$

These are the models with Hamiltonian H$_{6}$$_{v}$$_{,}$$_{1}$ , H$_{6}$$_{v}$$_{,}$$_{2}$ in appendix A. Figure 7(a) compares the training for a generic Hamiltonian with coefficients satisfying none of the above conditions against the training for H$_{6}$$_{v}$$_{,}$$_{1}$ -type model. We see that while the Hamiltonian constraint (3.13) saturates to similarly low values in both cases, the YB loss saturates at approximately one order of magnitude higher.

Similar behavior holds for the non-XYZ type models as well. The training for a generic class-4 Hamiltonian with coefficients a$_{1}$ = 0 . 5 , a$_{2}$ = 0 . 3 , a$_{3}$ = 0 . 4 , a$_{4}$ = 0 . 9 (see equation (A.19)) and a non-integrable deformation is shown in figure 7(b).

One can further discriminate between integrable and non-integrable models by checking the point-wise values of the YB losses in the two cases. Let us define the metric

$$
\tilde{\mathcal{L}} = \frac{\| \mathcal{R}_{12}(u-v) \mathcal{R}_{13}(u) \mathcal{R}_{23}(v) - \mathcal{R}_{23}(v) \mathcal{R}_{13}(u) \mathcal{R}_{12}(u-v) \|}{\| \mathcal{R}_{12}(u-v) \mathcal{R}_{13}(u) \mathcal{R}_{23}(v) \|}, \quad (4.13)
$$

which measures the relative error in the approximate solution of the YB equation. This metric is evaluated for the trained R-matrix for both integrable and non-integrable models in figure 8(for H$_{6}$$_{v}$$_{,}$$_{1}$ model), and figure 9(for H$_{class-4}$ model). We see that the normalized error can be up to two orders of magnitude larger for
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Comparing the training history of the Type XYZ and non-XYZ models against corresponding non-integrable Hamiltonians. There is approximately an order of magnitude difference between the YB losses for the integrable case vs the non-integrable case after the training saturates, indicated by the gray region in the graph. The step-wise drops in the loss functions approximately correspond to the learning rate schedule. The presented Hamiltonians are the same as on figures 8 and 9.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** (a) The normalized YB error (4.13) plotted in the logarithmic scale at the end of training for the Hamiltonian $H_{6v,1}$ with $a_1 = 0.3, a_2 = 0.3, b_1 = 0.45, b_2 = 0.6, c_1 = 0.4, c_2 = 0.25$, and (b) its non-integrable deformation with $a_1 = -1.3, a_2 = 1.3$ and other parameters kept constant. In order to keep all three arguments appearing in YB equation inside the same interval $|u|, |v|, |u-v| < 0.8$ we cut out the area $|u-v| > 0.8$ with chess-pattern triangles.</figcaption></figure>

the non-integrable case. Note that irrespective of the choice of Hamiltonian, there are two lines along $u=v$ and $v=0$ on which the YB equation is trivially satisfied, due to regularity. This metric also can detect anomalous situations when the learned solution once satisfied the Hamiltonian constraint at $u=0$ quickly evolves to a true solution of YB equation producing relatively small YB loss (3.7). In this case we will see the big spike in (4.13) around zero which will indicate the fakeness of the found solution.

The above consideration shows that one can define the metrics which together indicate the closeness of the given system to the integrable Hamiltonian. However, the final conclusion in the binary form of 'integrable/nonintegrable' regarding the given spin chain can be made only asymptotically, namely increasing the number of neurons, density of points and training time one can get the normalized YB loss (4.13) uniformly decreasing to zero for integrable Hamiltonians while for nonintegrable case it will be bounded from below by some positive value. Also let us stress that such problem is specific for the solver mode once we stick to a given Hamiltonian, while in the case of relaxed Hamiltonian restrictions as we will see in the next section, the neural network moves to the true solution of the YB equation.

## 4.2. Explorer: new from existing

In this section we will present two kinds of experiments that illustrate how the neural network presented above can be used to scan the landscape of spin-chains with two-dimensional space for integrable models. The training schedule adopted in this section is visualized in figure 10 and relies essentially on two new ingredients which distinguish it from the previous solver framework. These are *warm-start* and *repulsion*. We
<figure><img src="image_14.png" /><figcaption>**Figure 9.** (a) The normalized Yang-Baxter error (4.13) plotted in the logarithmic scale at the end of training for the Hamiltonian $H_{\text{class-4}}$ from (A.19), with $a_1 = 0.5, a_2 = 0.3, a_3 = 0.4, a_4 = 0.9$. (b) Non-integrable deformation with same Hamiltonian parameters as in the integrable case, except for $H_{13} = -0.9$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Visualizing the Explorer scheme. We start with random initializations, marked by lightning symbols, and perform solver learning represented by red curve arrows. Once we reach an submanifold of integrable Hamiltonians, we explore it using repulsion to identify new integrable models.</figcaption></figure>

will illustrate each by an example. In the first case we shall simply use warm-start, and in the second, we shall combine warm-start with repulsion. Finally, we shall use unsupervised learning methods such as t-SNE and Independent Component Analysis to identify distinct classes of Hamiltonians within the set of integrable models thus discovered. Collectively, these strategies make up our *explorer* framework.

The first key new ingredient is a *warm-start initialization*. As mentioned previously, the standard solver framework of the previous section uses He initialization [74] to instantiate the weights and biases of the neural network. In warm-start initialization, we use the knowledge of integrable systems previously discovered by the neural network to find new systems in its vicinity. The idea, at least intuitively, is that it should be possible to find new integrable systems more efficiently than with the random initialization by exploring the vicinity in weight-space of previously determined solutions using an iterative procedure such as gradient-based optimization. On doing so, we find a significant acceleration in training convergence, with new solutions being discovered typically in about 5 epochs of training after warm-start initialization. For definiteness, we consider the hermitian XYZ model discussed earlier in section 4.1.1. This has a two-parameter family of solutions, corresponding to independent choices for the parameters $\eta$ and $m$ of the Jacobi elliptic function, as seen from equation (4.4). The XXZ model is embedded into this space as the $m=0$ subspace of solutions.

We now describe how the above strategy can be used to quickly generate the cluster of XYZ R-matrices starting from a particular one which we choose from XXZ subclass. We begin with pre-training our neural
<figure><img src="image_16.png" /></figure>

network using the solver mode of the previous section, but with the learning rate of the Adam optimizer set to $10^{-3}$. The pre-training is stopped when all losses saturate below $O(10^{-3})$, which typically requires about 50 epochs of training. We carried out this pre-training setting arbitrary reference values of $\eta$, but with $m$ fixed to zero. The results shown here correspond to $\eta = \frac{\pi}{4}$. The weights thus obtained correspond to our warm-start values. Then we shift the target Hamiltonian values to correspond to $\eta \to \eta + \delta\eta$, where $\delta\eta$ are randomly chosen $O(10^{-1})$ numbers, and $m$ can take on non-zero values as well. We then retrain the model with a smaller learning rate, $10^{-4}$ for a few epochs until all loss terms fall to $O(10^{-4})$, which typically takes about 5 epochs, upon which we update the target Hamiltonian by updating $\eta$ and $m$ and continue training. This strategy generates about 10 XYZ models within the same time-scale (i.e. about 100 to 200 epochs of training) as we earlier needed for a single model. For best results, while we randomly update $\eta$, we systematically anneal the modular parameter $m$ to upwards of zero. A sample of this training is visualized in figure 11(b).

Our next key new ingredient for the Explorer mode is repulsion, which is added to the previous strategy of warm-start initialization. In principle, it should allow us to rediscover all 14 classes of integrable spin chains. However, for sake of simplicity, we will illustrate it now with a toy-model example and return to the general analysis later [81]. Namely, we consider the class of 6-vertex Hamiltonians with unrestricted $a_1$ and $a_2$. It includes both integrable 6-vertex classes $H_{6v,1}$, $H_{6v,2}$ (A.2) and (A.4) as well as nonintegrable models. In order to mimic the general situation when all integrable classes intersect at zero, we begin by pre-training the neural network to a Hamiltonian belonging to the intersection of the classes $H_{6v,1}$ and $H_{6v,2}$, i.e. whose matrix element satisfy the constraints $a_1 = a_2$ and $a_1 + a_2 = b_1 + b_2$ simultaneously. The results mentioned in this paper correspond to setting

$$a_1 = a_2 = \frac{b_1 + b_2}{2}; \quad b_1 = 0.6, b_2 = 0.8, c_1 = 0.5, c_2 = 0.9. \qquad (4.14)$$

Having arrived at this model, we would like to navigate to neighboring models not by specifying target values of the Hamiltonian, but by scanning the neighborhood of the current model. To do so, we employ a two step strategy. First, we navigate to two<fnref n="11" /> new 6-vertex integrable Hamiltonians by random scanning the vicinity of the current model without giving specific target values. We shall use these new models as our warm-start points. From each of them, we navigate away by using the repulsion loss term (3.17) for 1 epoch, followed by training for another 5 epochs. Note in this step, we still train within the restricted class of 6-vertex models by fixing the corresponding entries of the R-matrix to zero. We repeat this schedule 25 times starting from either of the saved models. This way, we generate fifty 6-vertex integrable Hamiltonians with over 1% accuracy<fnref n="12" />. The training curve displaying how the YB loss evolves is shown in figure 12.

<footnote n="11">¹¹ We stop the scanning once we found a representative from each of two classes because we know that there are only two integrable families here. In general case one of course should generate sufficiently many points in order to find all classes. We will return to this subtle point later in [81].</footnote>

<footnote n="12">¹² If we further train the individual models for more epochs, we can improve the accuracy of the obtained solution to similar levels as obtained in the examples presented in section 4.1.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Time evolution of the Yang Baxter loss as the neural network explores the space of integrable Hamiltonians of 6-vertex models $H_{6v,1,6v,2}$ by repulsion. The loss evolves together until the 50th epoch after which it fragments slightly as the training converges to the two warm-start points on the 60th epoch. For the remaining epochs the losses fragment completely as the neural network seeks out different new Hamiltonians and is terminated when the loss reaches the neighborhood of $1 \times 10^{-4}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Clustering of Hamiltonians from the 2 classes of gauge-inequivalent 6-vertex models obtained by Explorer using repulsion from solution at intersection of both classes.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

The learnt models are classified into two classes using clusterisation methods as shown in figure 13. Figure 14 plots the trained models in terms of coordinates defined by the integrability conditions of the Hamiltonians $H_{6v,1}$, $H_{6v,2}$. Models lying near the two axes were classified correctly into the two classes in figure 13 with 100% accuracy.

## 5. Conclusions and future directions

In this paper we constructed a neural network for solving the YB equation (2.3) in various contexts. Firstly, it can learn the R-matrix corresponding to a given integrable Hamiltonian or search for an integrable spin chain and the corresponding R-matrix from a certain class specified by imposed symmetries or other restrictions. We refer to this as the solver mode. Next, in the explorer mode, it can search for new integrable models by scanning the space of Hamiltonians.

We demonstrated the use of our neural network on spin chains of difference form with two-dimensional local space. In the solver mode, the network successfully learns all fourteen distinct classes of R-matrices identified in [68] to accuracies of the order of 99.9-99.99%. We demonstrated the work of the Explorer mode, restricting the search to the space of spin chains containing both classes of 6-vertex models as well as nonintegrable Hamiltonians. Starting from the hamiltonian at the intersection of two classes, Explorer found 50 integrable Hamiltonians which after clusterisation clearly fall into two families corresponding to two integrable classes of 6-vertex model. Working in the explorer mode, we find that warm-starting our training from the vicinity of a previously learnt integrable model greatly speeds up convergence, allowing us
<figure><img src="image_20.png" /></figure>

to identify typically about 50 new integrable models in the same time that random initialization takes to
converge to a single model.

The main focus of this paper was creating the neural network architecture and demonstrating its
robustness in various solution generating frameworks using known integrable models as a testing ground.
However, we expect that this program can be extended to various scenarios such as the exploration and
classification of the space of integrable Hamiltonians with dimensionality of the local Hilbert space greater
than two. This would be of great interest since the general classification of models is currently limited to two
dimensions. Our experiments with exploration and clustering are a promising starting point in this regard.
In our setup the strategy is quite straightforward [81]. Because all integrable families of Hamiltonians can be
multiplied by arbitrary scalar, we should only scan the Hamiltonians on the unit sphere which is compact.
Scanning over sufficiently dense set of points on the sphere will allow us to identify integrable Hamiltonians
from various classes. Then we can use the Explorer to reconstruct the whole corresponding families and
perform clusterisation in order to identify them. We already found various integrable Hamiltonians with 15,
20, 21, 24, 27, 33, 36,... vertices for spin-chains with three-dimensional local Hilbert space and will present
them in our forthcoming paper [81]. On another footing, it would also be interesting to extend our study to
R-matrices of non-difference form as these are particularly relevant to the AdS/CFT correspondence [82–85].

While our network learns a numerical approximation to the R-matrix, it can also be useful for the
reconstruction of analytical solutions using symbolic regression [32, 86]. Alternately, one may try to use the
learnt numerical solution for the reconstruction of the symmetry algebra such as the Yangian and then arrive
at the analytical solution. Remarkably, machine learning is already proving helpful in the analysis of
symmetry in physical systems. In particular, one may verify the presence of a conjectured symmetry or even
automate its search using machine learning [19, 46–50, 87]. It would be very interesting to explicate the
interplay of our program in this broader line of investigation.

In addition, the flexibility of our approach would also allow us to implement various additional
symmetries or other restrictions, both at the level of the R-matrix and the Hamiltonian. It would therefore be
very interesting to develop an ‘R-matrix bootstrap’ in the spirit of the two-dimensional S-matrix bootstrap
and analyze the interplay between various symmetries. For example, all 14 families of R-matrices considered
in this paper satisfy the condition of braided unitarity (2.12) and it would be interesting to rediscover them
from the use of braided unitarity and other symmetries without imposing the YB condition, similar to how
integrable two-dimensional S-matrices have been identified in the S-matrix Bootstrap approach [54, 56, 88].
With mild modifications, we can adapt our architecture to the analysis of YB equation for the integrable S-matrices in two dimensions. The only new feature to implement is the analytic structure in the θ-plane. It can be naturally realized with the use of holomorphic networks.

YB equation and its generalizations appear in the context of quantum information where R-matrix plays the role of an operator that creates entanglement between qubits [89–92]. It would be very interesting to apply our technique to the construction of new universal quantum gates.

Learning solutions for different classes with the same architecture, we noticed that the number of epochs needed to reach the same precision varies for different classes while being roughly the same for the Hamiltonians from the same classes. Thus, it would be very tempting to use the training of losses to define the complexity of spin chains. Ideally, we should be able to go beyond the class of integrable models and see that they sit at the minima of complexity, matching common beliefs that the integrable models are the ‘simplest’ ones.

## Data availability statement

No new data were created or analysed in this study.

## Acknowledgments

We thank Jiakang Bao, Jim Halverson, Ed Hirst, Vladimir Kazakov, Sven Krippendorf, Praneeth Netrapalli, Hongfei Shu, and Hao Zhou for interesting discussions. We especially thank Yang-Hui He for initial collaboration and several helpful discussions. SM thanks participants at String Data 2022, where initial version of the work was presented. SM is grateful to NORDITA for hospitality, during which part of the work was completed.

## Appendix A. Classes of 2D integrable spin chains of difference form

In this appendix, we list the 14 gauge-inequivalent integrable Hamiltonians of difference form and the corresponding R-matrices. Amongst the XYZ type models, the simplest solution is a diagonal 4-vertex model with Hamiltonians and R-matrices as follows:

$$H_{4v} = \begin{pmatrix} a_1 & 0 & 0 & 0 \\ 0 & b_1 & 0 & 0 \\ 0 & 0 & b_2 & 0 \\ 0 & 0 & 0 & a_2 \end{pmatrix} \leftrightarrow R_{4v}(u) = \begin{pmatrix} e^{a_1 u} & 0 & 0 & 0 \\ 0 & 0 & e^{b_2 u} & 0 \\ 0 & e^{b_1 u} & 0 & 0 \\ 0 & 0 & 0 & e^{a_2 u} \end{pmatrix}. \quad (A.1)$$

Figure 15 plots the training curve for R-matrix components as ratios with respect to (00) component, against the analytic functions for parameters $a_1 = 0.9, b_1 = 0.4, b_2 = 0.6, a_2 = 0.75$.

In 6-vertex models, we have two distinct classes depending on whether the Hamiltonian entries $H^{00}$ and $H^{33}$ are equal or not. In the first case, the R-matrix $R_{6v,1}(u)$ and its associated Hamiltonian $H_{6v,1}$ are given by

$$H_{6v,1} = \begin{pmatrix} a_1 & 0 & 0 & 0 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_2 & b_2 & 0 \\ 0 & 0 & 0 & a_1 \end{pmatrix} \leftrightarrow R_{6v,1}(u) = \begin{pmatrix} R_{6v,1}^{00}(u) & 0 & 0 & 0 \\ 0 & R_{6v,1}^{11}(u) & R_{6v,1}^{12}(u) & 0 \\ 0 & R_{6v,1}^{21}(u) & R_{6v,1}^{22}(u) & 0 \\ 0 & 0 & 0 & R_{6v,1}^{33}(u) \end{pmatrix} \quad (A.2)$$

where

$$\begin{align*} R_{6v,1}^{00}(u) &= R_{6v,1}^{33}(u) = e^{(b_1+b_2)u/2} \left( \cosh(\omega u) + \frac{2a_1 - b_1 - b_2}{2\omega} \sinh(\omega u) \right) \\ R_{6v,1}^{11}(u) &= \frac{c_2}{\omega} e^{(b_1+b_2)u/2} \sinh(\omega u) \\ R_{6v,1}^{12}(u) &= e^{b_2 u} \\ R_{6v,1}^{21}(u) &= e^{b_1 u} \\ R_{6v,1}^{22}(u) &= \frac{c_1}{\omega} e^{(b_1+b_2)u/2} \sinh(\omega u), \quad \omega = \frac{\sqrt{(2a_1 - b_1 - b_2)^2 - 4c_1 c_2}}{2}. \quad (A.3) \end{align*}$$

Figure 4 gives a representative training vs actual plot for this class.
<figure><img src="image_21.png" /><figcaption>Figure 15. (a) 4-vertex model, with $H_2$ parameters $a_1 = 0.9, b_1 = 0.4, b_2 = 0.6, a_2 = 0.75$, (b) errors.</figcaption></figure>

For the case $H^{00} \neq H^{33}$, the R-matrix $R_{6v,2}(u)$ is given by

$$H_{6v,2} = \begin{pmatrix} a_1 & 0 & 0 & 0 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_2 & b_2 & 0 \\ 0 & 0 & 0 & a_2 \end{pmatrix} \leftrightarrow R_{6v,2}(u) = \begin{pmatrix} R_{6v,2}^{00}(u) & 0 & 0 & 0 \\ 0 & R_{6v,2}^{11}(u) & R_{6v,2}^{12}(u) & 0 \\ 0 & R_{6v,2}^{21}(u) & R_{6v,2}^{22}(u) & 0 \\ 0 & 0 & 0 & R_{6v,2}^{33}(u) \end{pmatrix} \quad (A.4)$$

where $a_2 = b_1 + b_2 - a_1$ and

$$R_{6v,2}^{00}(u) = e^{(a_1+a_2)u/2} \left( \cosh(\omega u) + \frac{a_1-a_2}{2\omega} \sinh(\omega u) \right)$$

$$R_{6v,2}^{11}(u) = \frac{c_2}{\omega} e^{(a_1+a_2)u/2} \sinh(\omega u)$$

$$R_{6v,2}^{12}(u) = e^{b_2 u}$$

$$R_{6v,2}^{21}(u) = e^{b_1 u}$$

$$R_{6v,2}^{22}(u) = \frac{c_1}{\omega} e^{(a_1+a_2)u/2} \sinh(\omega u)$$

$$R_{6v,2}^{33}(u) = e^{(a_1+a_2)u/2} \left( \cosh(\omega u) - \frac{a_1-a_2}{2\omega} \sinh(\omega u) \right), \quad \omega = \frac{\sqrt{(a_1-a_2)^2 - 4c_1c_2}}{2} \quad (A.5)$$

Figure 16 gives a representative training vs actual plot for Hamiltonian parameters $a_1 = 1., a_2 = 0.2, b_1 = 0.45, b_2 = 0.75, c_1 = 0.4, c_2 = 0.6$. Next we have the 7-vertex models, which consists of two classes of solution distinguished by the Hamiltonian entries $H^{00}, H^{33}$ being equal or not. In the first case, we have

$$H_{7v,1} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & a_1 + b_1 & c_1 & 0 \\ 0 & -c_1 & a_1 - b_1 & 0 \\ 0 & 0 & 0 & a_1 \end{pmatrix} \leftrightarrow R_{7v,1}(u) = \begin{pmatrix} R_{7v,1}^{00}(u) & 0 & 0 & R_{7v,1}^{03}(u) \\ 0 & R_{7v,1}^{11}(u) & R_{7v,1}^{12}(u) & 0 \\ 0 & R_{7v,1}^{21}(u) & R_{7v,1}^{22}(u) & 0 \\ 0 & 0 & 0 & R_{7v,1}^{33}(u) \end{pmatrix} \quad (A.6)$$

where

$$R_{7v,1}^{00}(u) = R_{7v,1}^{33}(u) = e^{a_1 u} \cosh(c_1 u)$$

$$R_{7v,1}^{11}(u) = -R_{7v,1}^{22}(u) = e^{a_1 u} \sinh(c_1 u)$$

$$R_{7v,1}^{12}(u) = e^{(a_1 - b_1)u}$$
<figure><img src="image_22.png" /><figcaption>Figure 16. (a) 6-vertex model with Hamiltonian of type $H_{6v,2}$, with parameters $a_1 = 1, a_2 = 0.2, b_1 = 0.45, b_2 = 0.75, c_1 = 0.4, c_2 = 0.6$, (b) errors.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 17. (a) 7-vertex model with Hamiltonian of type $H_{7v,1}$, with parameters $a_1 = 1, b_1 = 0.45, c_1 = 0.6, d_1 = 0.75$, (b) errors.</figcaption></figure>

$$ 
\begin{aligned}
R_{7v,1}^{21}(u) &= e^{(a_1+b_1)u} \\
R_{7v,1}^{03}(u) &= -\frac{d_1}{2b_1} \left( e^{(a_1-b_1)u} - e^{(a_1+b_1)u} \right).
\end{aligned} 
\qquad (A.7) 
$$

Figure 17 plots the predicted R-matrix components as ratios with respect to the (12) component against the above analytic results, and their differences for a generic choice of parameters

$a_1 = 1, b_1 = 0.45, c_1 = 0.6, d_1 = 0.75$.

In the second case for $H^{00} \neq H^{33}$, we have

$$ H_{7v,2} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & a_1 - c_2 & c_1 & 0 \\ 0 & c_2 & a_1 - c_1 & 0 \\ 0 & 0 & 0 & a_2 \end{pmatrix} \leftrightarrow R_{7v,2}(u) = \begin{pmatrix} R_{7v,2}^{00}(u) & 0 & 0 & R_{7v,1}^{03}(u) \\ 0 & R_{7v,2}^{11}(u) & R_{7v,1}^{12}(u) & 0 \\ 0 & R_{7v,2}^{21}(u) & R_{7v,1}^{22}(u) & 0 \\ 0 & 0 & 0 & R_{7v,1}^{33}(u) \end{pmatrix} \quad (A.8) $$
<figure><img src="image_24.png" /><figcaption>Figure 18. (a) 7-vertex model with Hamiltonian of type $H_{7v,2}$, with parameters $a_1 = 1, c_1 = 0.45, c_2 = 0.75, d_1 = 0.5$, (b) errors.</figcaption></figure>

where $a_2 = a_1 - c_1 - c_2$ and

$$
\begin{align*}
R_{7v,2}^{00}(u) &= \frac{e^{(a_1 - \frac{c_1+c_2}{2})u}}{c_1 - c_2} \left( (c_1 - c_2) \cosh\left(\frac{c_1 - c_2}{2}u\right) + (c_1 + c_2) \sinh\left(\frac{c_1 - c_2}{2}u\right) \right) \\
R_{7v,2}^{11}(u) &= \frac{2c_2}{c_1 - c_2} e^{(a_1 - \frac{c_1+c_2}{2})u} \sinh\left(\frac{c_1 - c_2}{2}u\right) \\
R_{7v,2}^{22}(u) &= \frac{2c_1}{c_1 - c_2} e^{(a_1 - \frac{c_1+c_2}{2})u} \sinh\left(\frac{c_1 - c_2}{2}u\right) \\
R_{7v,2}^{12}(u) &= e^{(a_1 - c_1)u} \\
R_{7v,2}^{21}(u) &= e^{(a_1 + c_2)u} \\
R_{7v,2}^{03}(u) &= \frac{2d_1}{c_1 - c_2} e^{(a_1 - \frac{c_1+c_2}{2})u} \sinh\left(\frac{c_1 - c_2}{2}u\right) \\
R_{7v,2}^{33}(u) &= \frac{e^{(a_1 - \frac{c_1+c_2}{2})u}}{c_1 - c_2} \left( (c_1 - c_2) \cosh\left(\frac{c_1 - c_2}{2}u\right) - (c_1 + c_2) \sinh\left(\frac{c_1 - c_2}{2}u\right) \right) \tag{A.9}
\end{align*}
$$

Figure 18 plots the predicted R-matrix components as ratios with respect to the (12) component against the above analytic results, and their differences for a generic choice of parameters

$a_1 = 1, c_1 = 0.45, c_2 = 0.75, d_1 = 0.5.$

8-vertex models have 3 gauge-inequivalent classes labelled $H_{8v,i}$, $i = 1, 2, 3$. One of these models, namely $H_{8v,1}$, is a generalisation of the XYZ model

$$
H_{8v,1} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_1 & b_1 & 0 \\ d_2 & 0 & 0 & a_1 \end{pmatrix} \leftrightarrow R_{8v,1}(u) = \begin{pmatrix} R_{8v,1}^{00}(u) & 0 & 0 & R_{8v,1}^{03}(u) \\ 0 & R_{8v,1}^{11}(u) & R_{8v,1}^{12}(u) & 0 \\ 0 & R_{8v,1}^{21}(u) & R_{8v,1}^{22}(u) & 0 \\ R_{8v,1}^{30}(u) & 0 & 0 & R_{8v,1}^{33}(u) \end{pmatrix} \quad (A.10)
$$

where

$$
\begin{align*}
R_{8v,1}^{00}(u) &= R_{8v,1}^{33}(u) = \frac{\operatorname{sn}(u + 2\eta, m)}{\operatorname{sn}(2\eta, m)} e^{b_1 u} \\
R_{8v,1}^{11}(u) &= R_{8v,1}^{22}(u) = \frac{\operatorname{sn}(u, m)}{\operatorname{sn}(2\eta, m)} e^{b_1 u} \\
R_{8v,1}^{12}(u) &= R_{8v,1}^{21}(u) = e^{b_1 u} \\
R_{8v,1}^{03}(u) &= \sqrt{\frac{d_1}{d_2}} \sqrt{m} \operatorname{sn}(u + 2\eta, m) \operatorname{sn}(u, m) e^{b_1 u} \\
R_{8v,1}^{30}(u) &= \sqrt{\frac{d_2}{d_1}} \sqrt{m} \operatorname{sn}(u + 2\eta, m) \operatorname{sn}(u, m) e^{b_1 u} \tag{A.11}
\end{align*}
$$
<figure><img src="image_25.png" /><figcaption>Figure 19. (a) 8-vertex model with Hamiltonian of type $H_{8v,1}$, with parameters $b_1 = 0.4, \eta = 0.8, m = 0.5, \delta_1 = 1.3, \delta_2 = 0.7$, (b) errors.</figcaption></figure>

with Hamiltonian coefficients given by

$$a_1 = b_1 + \frac{\mathrm{cn}(2\eta, m) \mathrm{dn}(2\eta, m)}{\mathrm{sn}(2\eta, m)}, \quad c_1 = \frac{1}{\mathrm{sn}(2\eta, m)}, \quad d_1 = \sqrt{m}\delta_1 \mathrm{sn}(2\eta, m), \quad d_2 = \sqrt{m}\delta_2 \mathrm{sn}(2\eta, m) \quad (A.12)$$

for free parameters $b_1, \eta, m, \delta_1, \delta_2$. Figure 19 plots the predicted R-matrix components as ratios with respect to the (12) component against the above analytic results, and their differences for a generic choice of parameters $b_1 = 0.4, \eta = 0.8, m = 0.5, \delta_1 = 1.3, \delta_2 = 0.7$.

The second class of 8-vertex XYZ-type solution has Hamiltonian $H_{8v,2}$ and R-matrix $R_{8v,2}(u)$ defined as follows

$$H_{8v,2} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & b_1 & c_1 & 0 \\ 0 & c_1 & b_1 & 0 \\ d_2 & 0 & 0 & 2b_1 - a_1 \end{pmatrix} \leftrightarrow R_{8v,2}(u) = \begin{pmatrix} R_{8v,2}^{00}(u) & 0 & 0 & R_{8v,2}^{03}(u) \\ 0 & R_{8v,2}^{11}(u) & R_{8v,2}^{12}(u) & 0 \\ 0 & R_{8v,2}^{21}(u) & R_{8v,2}^{22}(u) & 0 \\ R_{8v,2}^{30}(u) & 0 & 0 & R_{8v,2}^{33}(u) \end{pmatrix} \quad (A.13)$$

where

$$\begin{align*}
R_{8v,2}^{00}(u) &= \left( \frac{\mathrm{cn}(u, m)}{\mathrm{dn}(u, m)} + \frac{\mathrm{sn}(u, m)\mathrm{cn}(2\eta, m)}{\mathrm{sn}(2\eta, m)} \right) e^{b_1 u} \\
R_{8v,2}^{11}(u) &= R_{8v,1}^{22}(u) = \frac{\mathrm{sn}(u, m)}{\mathrm{sn}(2\eta, m)} e^{b_1 u} \\
R_{8v,2}^{12}(u) &= R_{8v,1}^{21}(u) = e^{b_1 u} \\
R_{8v,2}^{03}(u) &= \frac{\delta_1}{\beta_1} \frac{\mathrm{sn}(u, m)\mathrm{cn}(u, m)}{\mathrm{dn}(u, m)\mathrm{sn}(2\eta, m)} e^{b_1 u} \\
R_{8v,2}^{30}(u) &= \frac{\delta_2}{\beta_1} \frac{\mathrm{sn}(u, m)\mathrm{cn}(u, m)}{\mathrm{dn}(u, m)\mathrm{sn}(2\eta, m)} e^{b_1 u} \\
R_{8v,2}^{33}(u) &= \left( \frac{\mathrm{cn}(u, m)}{\mathrm{dn}(u, m)} - \frac{\mathrm{sn}(u, m)\mathrm{cn}(2\eta, m)}{\mathrm{sn}(2\eta, m)} \right) e^{b_1 u} \tag{A.14}
\end{align*}$$

with the Hamiltonian coefficients given by

$$a_1 = b_1 + \frac{\mathrm{cn}(2\eta, m)}{\mathrm{sn}(2\eta, m)}, \quad c_1 = \frac{1}{\mathrm{sn}(2\eta, m)}, \quad d_1 = \frac{\delta_1}{\beta_1 \mathrm{sn}(2\eta, m)}, \quad d_2 = \frac{\delta_2}{\beta_1 \mathrm{sn}(2\eta, m)} \quad (A.15)$$

$$m = \frac{\delta_1 \delta_2}{\alpha_1^2 - \beta_1^2}, \quad \mathrm{cn}(2\eta, m) = \frac{\alpha_1}{\beta_1}, \quad \mathrm{sn}(2\eta, m) = \sqrt{1 - \frac{\alpha_1^2}{\beta_1^2}} \quad (A.16)$$
<figure><img src="image_26.png" /><figcaption>Figure 20. (a) 8-vertex model with Hamiltonian of type $H_{8v,2}$, with parameters $b_1 = 0.4, \alpha_1 = 0.5, \beta_1 = 0.7, \delta_1 = 0.3, \delta_2 = 0.2$, (b) errors.</figcaption></figure>

for free parameters $b_1, \alpha_1, \beta_1, \delta_1, \delta_2$. Figure 20 plots the predicted R-matrix components as ratios with respect to the (12) component against the above analytic results, and their differences for a generic choice of parameters $b_1 = 0.4, \alpha_1 = 0.5, \beta_1 = 0.7, \delta_1 = 0.3, \delta_2 = 0.2$.

The third class of 8-vertex XYZ-type solution has Hamiltonian $H_{8v,3}$ and R-matrix $R_{8v,3}(u)$ defined as follows

$$H_{8v,3} = \begin{pmatrix} a_1 & 0 & 0 & d_1 \\ 0 & a_1 & -b_1 & 0 \\ 0 & b_1 & a_1 & 0 \\ d_2 & 0 & 0 & a_1 \end{pmatrix} \leftrightarrow R_{8v,3}(u) = \begin{pmatrix} R_{8v,3}^{00}(u) & 0 & 0 & R_{8v,3}^{03}(u) \\ 0 & R_{8v,3}^{11}(u) & R_{8v,3}^{12}(u) & 0 \\ 0 & R_{8v,3}^{21}(u) & R_{8v,3}^{22}(u) & 0 \\ R_{8v,3}^{30}(u) & 0 & 0 & R_{8v,3}^{33}(u) \end{pmatrix} \quad (A.17)$$

where

$$ \begin{align} R_{8v,3}^{00}(u) &= R_{8v,1}^{33}(u) = \frac{\cosh(b_1 u)}{\cos(\sqrt{d_1 d_2} u)} e^{a_1 u} \nonumber \\ R_{8v,3}^{11}(u) &= -R_{8v,1}^{22}(u) = \frac{\sinh(b_1 u)}{\cos(\sqrt{d_1 d_2} u)} e^{a_1 u} \nonumber \\ R_{8v,3}^{12}(u) &= R_{8v,1}^{21}(u) = e^{a_1 u} \nonumber \\ R_{8v,3}^{03}(u) &= \sqrt{\frac{d_1}{d_2}} e^{a_1 u} \tan(\sqrt{d_1 d_2} u) \nonumber \\ R_{8v,1}^{30}(u) &= \sqrt{\frac{d_2}{d_1}} e^{a_1 u} \tan(\sqrt{d_1 d_2} u) \tag{A.18} \end{align} $$

Figure 21 plots the predicted R-matrix components as ratios with respect to the (12) component against the above analytic results, and their differences for a generic choice of parameters $a_1 = 1, b_1 = -0.45, d_1 = 0.6, d_2 = 0.75$.
<figure><img src="image_27.png" /><figcaption>Figure 21. (a) 8-vertex model with Hamiltonian of type $H_{8v,3}$, with parameters $a_1 = 1$, $b_1 = -0.45$, $d_1 = 0.6$, $d_2 = 0.75$, (b) errors.</figcaption></figure>

For non-XYZ type models, the 6 gauge-inequivalent Hamiltonians are of the form

$$H_{\text{class-1}} = \begin{pmatrix} 0 & a_1 & a_2 & 0 \\ 0 & a_5 & 0 & a_3 \\ 0 & 0 & -a_5 & a_4 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad \text{with } a_1 a_3 = a_2 a_4; \quad H_{\text{class-2}} = \begin{pmatrix} 0 & a_2 & a_3 - a_2 & a_5 \\ 0 & a_1 & 0 & a_4 \\ 0 & 0 & -a_1 & a_3 - a_4 \\ 0 & 0 & 0 & 0 \end{pmatrix},$$

$$H_{\text{class-3}} = \begin{pmatrix} -a_1 & (2a_1 - a_2)a_3 & (2a_1 + a_2)a_3 & 0 \\ 0 & a_1 - a_2 & 0 & 0 \\ 0 & 0 & a_1 + a_2 & 0 \\ 0 & 0 & 0 & -a_1 \end{pmatrix}, \quad H_{\text{class-4}} = \begin{pmatrix} a_1 & a_2 & a_2 & a_3 \\ 0 & -a_1 & 0 & a_4 \\ 0 & 0 & -a_1 & a_4 \\ 0 & 0 & 0 & a_1 \end{pmatrix},$$

$$H_{\text{class-5}} = \begin{pmatrix} a_1 & a_2 & -a_2 & 0 \\ 0 & -a_1 & 2a_1 & a_3 \\ 0 & 2a_1 & -a_1 & -a_3 \\ 0 & 0 & 0 & a_1 \end{pmatrix}, \quad H_{\text{class-6}} = \begin{pmatrix} a_1 & a_2 & a_2 & 0 \\ 0 & -a_1 & 2a_1 & -a_2 \\ 0 & 2a_1 & -a_1 & -a_2 \\ 0 & 0 & 0 & a_1 \end{pmatrix} \quad (\text{A.19})$$

Corresponding R-matrices are

$$R_{\text{class-1}}(u) = \begin{pmatrix} 1 & \frac{a_1(e^{a_5 u} - 1)}{a_5} & \frac{a_2(e^{a_5 u} - 1)}{a_5} & \frac{(a_1 a_3 + a_2 a_4)}{a_5^2} (\cosh(a_5 u) - 1) \\ 0 & 0 & e^{-a_5 u} & \frac{a_4(1 - e^{-a_5 u})}{a_5} \\ 0 & e^{a_5 u} & 0 & \frac{a_3(1 - e^{-a_5 u})}{a_5} \\ 0 & 0 & 0 & 1 \end{pmatrix} \quad (\text{A.20})$$

$$R_{\text{class-2}}(u) = uP \left( \frac{a_1}{\sinh(a_1 u)} + H_{\text{class-5}} + \frac{\tanh(a_1 u)}{a_1} H_{\text{class-5}}^2 \right), \quad P = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \quad (\text{A.21})$$

$$R_{\text{class-3}}(u) = \begin{pmatrix} e^{-a_1 u} & a_3 (e^{(a_1 - a_2)u} - e^{-a_1 u}) & a_3 (e^{(a_1 + a_2)u} - e^{-a_1 u}) & 0 \\ 0 & 0 & e^{(a_1 + a_2)u} & 0 \\ 0 & e^{(a_1 - a_2)u} & 0 & 0 \\ 0 & 0 & 0 & e^{-a_1 u} \end{pmatrix} \quad (\text{A.22})$$

$$R_{\text{class-4}}(u) = \begin{pmatrix} e^{a_1 u} & \frac{a_2 \sinh(a_1 u)}{a_1} & \frac{a_2 \sinh(a_1 u)}{a_1} & \frac{e^{a_1 u} (a_2 a_4 + a_1 a_3 \coth(a_1 u)) \sinh^2(a_1 u)}{a_1^2} \\ 0 & 0 & e^{-a_1 u} & \frac{a_4 \sinh(a_1 u)}{a_1} \\ 0 & e^{-a_1 u} & 0 & \frac{a_4 \sinh(a_1 u)}{a_1} \\ 0 & 0 & 0 & e^{a_1 u} \end{pmatrix} \quad (\text{A.23})$$
<figure><img src="image_28.png" /><figcaption>Figure 22. (a), (b) class 2, with $H_2$ parameters $a_1 = -0.6$, $a_2 = 0.381 + 0.123i$, $a_3 = 0.447i$, $a_4 = 0.7i$, $a_5 = -0.3i$: real and imaginary parts, (c), (d) class 3, with $H_2$ parameters $a_1 = 1$, $a_2 = 0.5$, $a_3 = 0.7$: real and imaginary parts, (e), (f) class 4, with $H_2$ parameters $a_1 = 0.5$, $a_2 = -0.5i$, $a_3 = 0.5$, $a_4 = -0.5$: real and imaginary parts.</figcaption></figure>

$$R_{\text{class-5}}(u) = (1 - a_1 u) \begin{pmatrix} 2a_1 u + 1 & a_2 u & -a_2 u & a_2 a_3 u^2 \\ 0 & 2a_1 u & 1 & -a_3 u \\ 0 & 1 & 2a_1 u & a_3 u \\ 0 & 0 & 0 & 2a_1 u + 1 \end{pmatrix} \quad (A.24)$$

$$R_{\text{class-6}}(u) = (1 - a_1 u) \begin{pmatrix} 2a_1 u + 1 & a_2 u (2a_1 u + 1) & a_2 u (2a_1 u + 1) & -a_2^2 u^2 (2a_1 u + 1) \\ 0 & 2a_1 u & 1 & -a_2 u (2a_1 u + 1) \\ 0 & 1 & 2a_1 u & -a_2 u (2a_1 u + 1) \\ 0 & 0 & 0 & 2a_1 u + 1 \end{pmatrix}. \quad (A.25)$$
<figure><img src="image_29.png" /><figcaption>Figure 23. (a), (b) class 5, with $H_2$ parameters $a_1 = -0.5$, $a_2 = -0.5i$, $a_3 = -0.5i$: real and imaginary parts.</figcaption></figure>

In the class-2 solution above, the non-zero R-matrix components are explicitly given by

$$
\begin{align*}
R_{\text{class-2}}^{00}(u) &= R_{\text{class-2}}^{33}(u) = \frac{a_1 u}{\sinh u}, \\
R_{\text{class-2}}^{01}(u) &= a_2 u \left(1 + \tanh\left(\frac{a_1 u}{2}\right)\right), \\
R_{\text{class-2}}^{02}(u) &= (a_2 - a_3) u \left(-1 + \tanh\left(\frac{a_1 u}{2}\right)\right), \\
R_{\text{class-2}}^{03}(u) &= a_5 u \left(1 + \frac{((a_4 - a_3)(a_2 - a_3) + a_2 a_4) \tanh \frac{a_1 u}{2}}{a_1 a_5}\right), \\
R_{\text{class-2}}^{12}(u) &= a_1 u \left(-1 + \frac{1}{\sinh(a_1 u)} + \tanh \frac{a_1 u}{2}\right), \\
R_{\text{class-2}}^{13}(u) &= (a_4 - a_3) u \left(-1 + \tanh\left(\frac{a_1 u}{2}\right)\right), \\
R_{\text{class-2}}^{21}(u) &= a_1 u \left(1 + \frac{1}{\sinh(a_1 u)} + \tanh \frac{a_1 u}{2}\right), \\
R_{\text{class-2}}^{23}(u) &= a_4 u \left(1 + \tanh\left(\frac{a_1 u}{2}\right)\right). \tag{A.26}
\end{align*}
$$

Amongst the above non-XYZ type models, we have already looked into the training for Class 1 model in section 4.1. Figures 22–24 plot the training vs actual R-matrix components for classes 2,3,4, class 5, and class 6 respectively, with generic Hamiltonian parameters. Also we note that allowing for complex parameters results in generically complex R-matrices. We compare the predictions against the actual formulae by taking ratios with respect to the real part of the (00) component for classes 2–5, and (12) component for class 6.
<figure><img src="image_30.png" /><figcaption>Figure 24. (a), (b) class 6, with $H_2$ parameters $a_1 = -0.4$, $a_2 = -0.6$: real and imaginary parts.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 25. Visualizing the map $\phi$ which is learnt by the Siamese architecture. The points $x_1$ and $x_3$ are similar to each other while $x_2$ is dissimilar to both of them.</figcaption></figure>

## Appendix B. Designing the neural network

This appendix contains an extensive overview of the architecture of our neural network solver, as well as details of the hyperparameters with which the network is trained. Our starting point is the close analogy between our problem of machine learning R-matrices by imposing constraints and the design of the Siamese neural networks [93, 94]. These were designed to function in settings where the canonical supervised learning approach of (3.5) for classification becomes infeasible due to the large number of target classes {$y$} and the paucity of training examples {$x_\alpha$} corresponding to each class $y_\alpha$. In such a situation, one may instead define a similarity relation

$$x_{\alpha_1} \sim x_{\alpha_2} \iff y(x_{\alpha_1}) = y(x_{\alpha_2}), \qquad (B.1)$$

and train the neural network to learn a function $\phi(x) : \mathbb{R}^D \to \mathbb{R}^d$ such that the Euclidean distance between representatives $\phi(x)$ of two input vectors $x_1, x_2$ that are similar to each other is small, while the distance between dissimilar data is large. Schematically,

$$d(x_a, x_b) = |\phi(x_a) - \phi(x_b)|^2 \approx 0 \iff x_a \sim x_b. \qquad (B.2)$$

This is visualized in figure 25.

There are many loss functions by which such networks may be trained, see for example [93–96]. For definiteness, we mention the contrastive loss function of [93, 94], given by

$$\mathcal{L} = Yd(x_1, x_2) + (1 - Y) \max(r_o - d(x_1, x_2), 0), \qquad (B.3)$$
<figure><img src="image_32.png" /><figcaption>**Figure 26.** Visualizing the forward propagation of the neural network $\mathcal{R}(u)$. This has a very strong parallel to figure 25, with the function $\mathcal{R}(u)$ playing the role of the map $\phi$. The only difference is that $\mathcal{R}(u)$ also has additional constraints on $\mathcal{R}(0)$ and $\mathcal{R}'(0)$ which are unique to our problem.</figcaption></figure>

where $Y = 1$ if $x_1 \sim x_2$ and $Y = 0$ otherwise. Clearly this loss function causes the network to learn a function $\phi$ such that similar inputs $x$ are clustered together while dissimilar inputs are pushed at least a distance $r_o$ apart. This therefore realizes our naive criterion for $\phi$ laid out in equation (B.2). We also see very explicitly that the loss function in equation (B.3) does not directly depend on the values $y$ in contrast to equation (3.5). Instead, the network is trained to learn a function $\phi(x)$ which obeys a property which is not given point-wise for each input $x$ but instead is expressed as a non-linear constraint (B.2) on $\phi(x)$ evaluated at two points $x_1$ and $x_2$.

## B.1. The neural network architecture and forward propagation

We now provide some more details about our implementation of $\mathcal{R}(u)$ and the training done to converge to solutions of the YB equation (2.3) consistent with additional requirements such as regularity (2.4). As already mentioned in section 3.2, each matrix element $\mathcal{R}_{ij}$ is decomposed into the sum $a_{ij} + ib_{ij}$ which are individually modeled by MLPs. In principle each MLP is independent of the rest and can be individually designed. We shall however take all MLPs to contain two hidden layers of 50 neurons each, followed by a single output neuron which is linear activated$^{13}$. The possible activations for the hidden layers are compared in appendix B.2 below. To proceed further, note that our loss function involves a term (3.7) which takes arguments $\mathcal{R}(u)$, $\mathcal{R}(v)$ and $\mathcal{R}(u-v)$ where $u, v$ are valued in $\Omega$. This clearly has a very strong parallel with the Siamese Networks introduced above. At least intuitively, one may regard our problem as training a 'triplet' of identical neural networks $\mathcal{R}$ to optimize the loss function (3.7). In addition however, we also have to train the function on loss functions such as (3.12) and (3.13). These constraints, along with the Siamese schematic shown in figure 25 motivate our design visualized in figure 26.

During the forward propagation we sample a minibatch of $u$ and $v$ values, from which the corresponding $u-v$ is constructed. Next, the $\mathcal{R}$ matrix is constructed at $u, v$ and $u-v$ via equation (3.6). We also evaluate $\mathcal{R}(0)$ and $\mathcal{R}'(0)$ thus completing the forward propagation. Next, we compute the losses (3.7), (3.12) and (3.13) as well as possibly (3.15). The loss function is trained on by using the Adam optimizer [73], with an initial learning rate $\eta$ of $10^{-3}$ which is annealed to $10^{-8}$ in steps of $10^{-1}$ by monitoring the saturation in the YB loss computed for the validation set over 5–10 epochs. The effect of this annealing in the learning rate is also visible in the training histories in figures 2 and 27 where the step-wise falls in the losses correspond to the drops in the learning rate. Across the board, training converges in about 100 epochs and is terminated by early stopping.

## B.2. Comparing different activation functions

We now turn to a brief comparison of the performance of different activation functions with the above set up. Again for uniformity, we will use one activation throughout for all the MLPs $a_{ij}$ and $b_{ij}$, but for the

<footnote n="13">13 One might also construct an alternate formulation of the neural network where a single MLP of the kind shown in figure 1 accepts a real input $u$ and outputs all the requisite real scalar functions that comprise $\mathcal{R}(u)$. So far we have observed that such a network does not perform as well as our current formulation of independent neural networks for each real function. Nonetheless, it is possible that this formulation may eventually prove competitive with our current one and the question remains under investigation currently.</footnote>
<figure><img src="image_33.png" /><figcaption>**Figure 27.** The evolution of the Yang-Baxter loss (left) and the Hamiltonian loss (right) for a variety of activation functions when training for 200 epochs. The swish activation tends to outperform the others. The precise numbers are given in table 1.</figcaption></figure>

<table><caption>**Table 1.** Performance of different activation functions on learning the Hamiltonian (B.4). The saturation epoch is the approximate epoch after which the model did not train further. The final values of the Yang-Baxter and Hamiltonian losses after saturation is also mentioned. We observed that the swish activation converges sooner and to lower losses. This is stable across multiple runs. See also figure 27.</caption><thead><tr><th>Activation</th><th>Final Yang-Baxter loss</th><th>Final Hamiltonian loss</th><th>Saturation epoch</th></tr></thead><tbody><tr><td>sigmoid</td><td>2.5 × 10<sup>-3</sup></td><td>6.12 × 10<sup>-7</sup></td><td>150</td></tr><tr><td>tanh</td><td>1.90 × 10<sup>-4</sup></td><td>5.25 × 10<sup>-7</sup></td><td>125</td></tr><tr><td>swish</td><td>6.49 × 10<sup>-5</sup></td><td>1.51 × 10<sup>-7</sup></td><td>75</td></tr><tr><td>elu</td><td>2.75 × 10<sup>-4</sup></td><td>5.63 × 10<sup>-7</sup></td><td>100</td></tr><tr><td>relu</td><td>5.52 × 10<sup>-4</sup></td><td>4.63 × 10<sup>-7</sup></td><td>100</td></tr></tbody></table>

output neuron which is linearly activated. We then compared the performance of this neural network architecture while learning the Hamiltonian

$$H_{6v,1} = \begin{pmatrix} 0.3 & 0 & 0 & 0 \\ 0 & 0.45 & 0.4 & 0 \\ 0 & 0.25 & 0.6 & 0 \\ 0 & 0 & 0 & 0.3 \end{pmatrix}, \qquad (B.4)$$

which is 6-vertex Type 1 in the classification of [68], see equation (A.2) above. The neural network was trained with the loss functions (3.7), (3.13) and (3.12) and setting $\lambda_H$ and $\lambda_{reg}$ to 1 each. The training was carried out for 200 epochs on observing that the networks did not perform better on training for longer. Further, we set a batch size of 16 and optimized using Adam with a starting learning rate of $10^{-3}$ which was annealed to $10^{-8}$ using the saturation in the YB loss over the validation set as the criterion as mentioned above. We conducted this training using the activations sigmoid, tanh, swish, all of which are holomorphic, as well as elu and relu. The last two are not holomorphic but have been included for completeness. The evolution of the YB and the Hamiltonian loss for all these activations is shown in figure 27 and table 1. On the whole, we see that the swish activation tends to outperform the others quite significantly. While these are the results of a single run, we found that the result is consistent across several runs and tasks, leading us to adopt the swish activation uniformly across the board for all the analyses shown in this paper.

### B.3. Proof of concept: training with a single hidden neuron

As a final observation we present a simple proof of concept of our approach of solving the YB equation along with other constraints by optimizing suitable loss functions. Here, instead of attempting to deep learn the solution, we use a single layer of neurons for each R-matrix function and pick an activation function by the form of the known analytic solution. In effect, rather than rely on the feature learning properties of a deep MLP as we have done in the rest of our paper, we ourselves provide activation functions which should furnish a natural basis to express the known analytic solutions in.
<figure><img src="image_34.png" /><figcaption>Figure 28. The figure on the left shows the XXZ model R-matrix with $\eta = \frac{\pi}{3}$ obtained by trained on a single sin activated hidden neuron in the range $u \in (-1, 1)$ shown in gray. The solution remains valid outside the training domain as well. The figure on the right shows the corresponding training curves.</figcaption></figure>

For definiteness, consider the XXZ model at $\eta = \frac{\pi}{3}$. The non-zero entries in this R-matrix are

$$R_{00}(u) = \sin(u + \eta), \quad R_{11}(u) = \sin(u - \eta) = R_{22}(u), \quad R_{12}(u) = \sin(\eta) = R_{21}(u), \quad (B.5)$$

as may be observed by setting $m=0$ in equation (4.4). We define the networks $a_{ij}$ and $b_{ij}$ to have a single hidden layer of a solitary neuron activated by the sin function. This means that the functions learnt by the network are simply of the form

$$a_{ij} = \tilde{W} \circ \sin(W \circ u), \quad (B.6)$$

and similarly for the $b_{ij}$. The $W$ and $\tilde{W}$ are the weight and bias of the hidden and the output neuron respectively and the composition $W \circ u$ is shorthand for the affine transformation $wu + b$. Next, we train the network imposing the losses (3.7), (3.12), (3.13) and (3.15), each with weight $\lambda = 1$, and the Adam optimizer with our standard learning rate scheduling. Figure 28 plots the trained XXZ R-matrix components for $u$ lying in the range $(-10, 10)$.

Note that since we trained with an activation function that presupposed our knowledge of the exact solution—in effect, the true R-matrix lay within our hypothesis class—the model trained to losses of the order of $10^{-8}$ which is several orders of magnitude below the typical end of training losses we observed in the standardized framework. Further, we also obtain an excellent performance even out of the domain of training, which is usually not the case in machine learning.

## ORCID iD

Evgeny Sobko https://orcid.org/0000-0003-1906-1359

## References

[1] LeCun Y, Bengio Y and Hinton G 2015 Deep learning Nature 521 436

[2] He Y-H 2017 Machine-learning the string landscape Phys. Lett. B 774 564

[3] Carifio J, Halverson J, Krioukov D and Nelson B D 2017 Machine learning in the string landscape J. High Energy Phys. JHEP09(2017)157

[4] Krefl D and Seong R-K 2017 Machine learning of Calabi-Yau volumes Phys. Rev. D 96 066014

[5] Ruehle F 2017 Evolving neural networks with genetic algorithms to study the String Landscape J. High Energy Phys. JHEP08(2017)038

[6] Brodie C R, Constantin A, Deen R and Lukas A 2020 Machine learning line bundle cohomology Fortsch. Phys. 68 1900087

[7] Deen R, He Y-H, Lee S-J and Lukas A 2022 Machine learning string standard models Phys. Rev. D 105 046001

[8] He Y-H and Lukas A 2021 Machine learning Calabi-Yau four-folds Phys. Lett. B 815 136139

[9] Erbin H and Finotello R 2021 Machine learning for complete intersection Calabi-Yau manifolds: a methodological study Phys. Rev. D 103 126014
[10] Erbin H, Finotello R, Schneider R and Tamaazousti M 2022 Deep multi-task mining Calabi–Yau four-folds *Mach. Learn. Sci. Tech.* 3 015006

[11] Gao X and Zou H 2022 Applying machine learning to the Calabi-Yau orientifolds with string vacua *Phys. Rev. D* 105 046017

[12] Ashmore A, Calmon L, He Y-H and Ovrut B A 2023 Calabi-Yau metrics, energy functionals and machine-learning *Int. J. Data Sci. Math. Sci.* 1 49

[13] Anderson L B, Gerdes M, Gray J, Krippendorf S, Raghuram N and Ruehle F 2021 Moduli-dependent Calabi-Yau and SU (3)-structure metrics from Machine Learning *J. High Energy Phys.* jhep04(2021)001

[14] Douglas M, Lakshminarasimhan S and Qi Y 2022 Numerical Calabi-Yau metrics from holomorphic networks *Mathematical and Scientific Machine Learning* pp 223–52

[15] Larfors M, Lukas A, Ruehle F and Schneider R 2022 Numerical metrics for complete intersection and Kreuzer-Skarke Calabi-Yau manifolds *Mach. Learn. Sci. Tech.* 3 035014

[16] He Y-H, Lal S and Zaz M Z 2021 The world in a grain of sand: condensing the string vacuum degeneracy (arXiv:2111.04761)

[17] Morningstar A and Melko R G 2017 Deep learning the ising model near criticality (arXiv:1708.04622)

[18] Zhang Y, Melko R G and Kim E-A 2017 Machine learning Z 2 quantum spin liquids with quasiparticle statistics *Phys. Rev. B* 96 245119

[19] Chen H-Y, He Y-H, Lal S and Zaz M Z 2020 Machine learning etudes in conformal field theories (arXiv:2006.16114)

[20] Kántor G, Niarchos V and Papageorgakis C 2022 Conformal bootstrap with reinforcement learning *Phys. Rev. D* 105 025018

[21] Kántor G, Niarchos V and Papageorgakis C 2022 Solving conformal field theories with artificial intelligence *Phys. Rev. Lett.* 128 041601

[22] Kántor G, Niarchos V, Papageorgakis C and Richmond P 2023 6D (2,0) bootstrap with the soft-actor-critic algorithm *Phys. Rev. D* 107 025005

[23] Kuo E-J, Seif A, Lundgren R, Whitsitt S and Hafezi M 2022 Decoding conformal field theories: from supervised to unsupervised learning *Phys. Rev. Res.* 4 043031

[24] Basu P, Bhattacharya J, Jakka D P S, Mosomane C and Shukla V 2022 Machine learning of Ising criticality with spin-shuffling (arXiv:2203.04012)

[25] Shiina K, Mori H, Okabe Y and Lee H K 2020 Machine-learning studies on spin models *Sci. Rep.* 10 2177

[26] Han X and Hartnoll S A 2020 Deep quantum geometry of matrices *Phys. Rev. X* 10 011069

[27] Weinan E and Yu B 2017 The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems (arXiv:1710.00211)

[28] Raissi M, Perdikaris P and Karniadakis G E 2017 Physics informed deep learning (part I): data-driven solutions of nonlinear partial differential equations (arXiv:1711.10561)

[29] Lample G and Charton F 2019 Deep learning for symbolic mathematics (arXiv:1912.01412)

[30] Davies A et al 2021 Advancing mathematics by guiding human intuition with AI *Nature 600* 70

[31] He Y-H 2021 Machine-learning mathematical structures (arXiv:2101.06317)

[32] Udrescu S-M and Tegmark M 2020 AI Feynman: a physics-inspired method for symbolic regression *Sci. Adv.* 6 eaay2631

[33] Cybenko G 1989 Approximation by superpositions of a sigmoidal function *Math. Control Signals Syst.* 2 303

[34] Hornik K, Stinchcombe M and White H 1989 Multilayer feedforward networks are universal approximators *Neural Netw.* 2 359

[35] Lu Z, Pu H, Wang F, Hu Z and Wang L 2017 The expressive power of neural networks: a view from the width *Advances in Neural Information Processing Systems* p 30

[36] Telgarsky M 2015 Representation benefits of deep feedforward networks (arXiv:1509.08101)

[37] Glorot X and Bengio Y 2010 Understanding the difficulty of training deep feedforward neural networks *Proc. Thirteenth Int. Conf. on Artificial Intelligence and Statistics* p 249

[38] Hinton G E, Srivastava N, Krizhevsky A, Sutskever I and Salakhutdinov R R 2012 Improving neural networks by preventing co-adaptation of feature detectors (arXiv:1207.0580)

[39] Srivastava N, Hinton G, Krizhevsky A, Sutskever I and Salakhutdinov R 2014 Dropout: a simple way to prevent neural networks from overfitting *J. Mach. Learn. Res.* 15 1929

[40] Ioffe S and Szegedy C 2015 Batch normalization: accelerating deep network training by reducing internal covariate shift *Int. Conf. on Machine Learning* pp 448–56

[41] Smith L N 2018 A disciplined approach to neural network hyper-parameters: Part 1–learning rate, batch size, momentum, and weight decay (arXiv:1803.09820)

[42] Zhang H, Dauphin Y N and Ma T 2019 Fixup initialization: residual learning without normalization (arXiv:1901.09321)

[43] Chollet F et al 2018 Keras: the python deep learning library *Astrophysics Source Code Library ascl*

[44] Abadi M et al TensorFlow: large-scale machine learning on heterogeneous systems (available at: www.tensorflow.org/)

[45] Paszke A et al 2019 PyTorch: an imperative style, high-performance deep learning library *Advances in Neural Information Processing Systems* vol 32 (Curran Associates, Inc.) pp 8024–35

[46] Liu Z and Tegmark M 2022 Machine learning hidden symmetries *Phys. Rev. Lett.* 128 180201

[47] Liu Z and Tegmark M 2021 Machine learning conservation laws from trajectories *Phys. Rev. Lett.* 126 180604

[48] Bondesan R and Lamacraft A 2019 Learning symmetries of classical integrable systems (arXiv:1906.04645)

[49] Wetzel S J, Melko R G, Scott J, Panju M and Ganesh V 2020 Discovering symmetry invariants and conserved quantities by interpreting siamese neural networks *Phys. Rev. Res.* 2 033499

[50] Forestano R T, Matchev K T, Matcheva K, Roman A, Unlu E and Verner S 2023 Deep learning symmetries and their lie groups, algebras, and subalgebras from first principles (arXiv:2301.05638)

[51] Chew G F 1962 S-matrix theory of strong interactions without elementary particles *Rev. Mod. Phys.* 34 394

[52] Eden R J, Landshoff P V, Olive D I and Polkinghorne J C 1966 *The Analytic S-Matrix* (Cambridge University Press)

[53] Kruczenski M, Penedones J and van Rees B C 2022 Snowmass white paper: S-matrix bootstrap (arXiv:2203.02421)

[54] Paulos M F, Penedones J, Toledo J, van Rees B C and Vieira P 2017 The S-matrix bootstrap II: two dimensional amplitudes *J. High Energy Phys.* JHEP11(2017)143

[55] Paulos M F, Penedones J, Toledo J, van Rees B C and Vieira P 2019 The S-matrix bootstrap. Part III: higher dimensional amplitudes *J. High Energy Phys.* JHEP12(2019)040

[56] He Y, Irrgang A and Kruczenski M 2018 A note on the S-matrix bootstrap for the 2d O(N) bosonic model *J. High Energy Phys.* JHEP11(2018)093

[57] Zamolodchikov A B and Zamolodchikov A B 1977 Relativistic factorized S matrix in two-dimensions having O(N) isotopic symmetry *JETP Lett.* 26 457
[58] Faddeev L D 1996 How algebraic Bethe ansatz works for integrable model *Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation* pp 149–219

[59] Perk J H H and Au-Yang H 2006 Yang-Baxter equations *Encycl. Math. Phys.* 5 465

[60] Kulish P P, Reshetikhin N Y and Sklyanin E K 1981 Yang-Baxter equation and representation theory: I *Lett. Math. Phys.* 5 393

[61] Jimbo M 1986 Quantum R matrix for the generalized toda system *Commun. Math. Phys.* 102 537

[62] Bazhanov V V 1985 Trigonometric solution of triangle equations and classical lie algebras *Phys. Lett.* B 159 321

[63] Bazhanov V V 1987 Integrable quantum systems and classical lie algebras *Commun. Math. Phys.* 113 471

[64] Sogo K, Uchinami M, Akutsu Y and Wadati M 1982 Classification of exactly solvable two-component models *Prog. Theor. Phys.* 68 508

[65] Perk J H H and Schultz C L 1981 New families of commuting transfer matrices in q state vertex models *Phys. Lett. A* 84 407

[66] Jones V F R 1991 Baxterization *Int. J. Mod. Phys. A* 6 2035

[67] Vieira R 2018 Solving and classifying the solutions of the Yang-Baxter equation through a differential approach. Two-state systems *J. High Energy Phys.* JHEP10(2018)110

[68] De Leeuw M, Pribytok A and Ryan P 2019 Classifying integrable spin-1/2 chains with nearest neighbour interactions *J. Phys. A: Math. Theor.* 52 505201

[69] de Leeuw M, Paletta C, Pribytok A, Retore A L and Ryan P 2020 Classifying nearest-neighbor interactions and deformations of AdS *Phys. Rev. Lett.* 125 031604

[70] de Leeuw M, Paletta C, Pribytok A, Retore A L and Ryan P 2021 Yang-Baxter and the boost: splitting the difference *SciPost Phys.* 11 069

[71] Krippendorf S, Lust D and Syvaeri M 2021 Integrability Ex Machina *Fortsch. Phys.* 69 2100057

[72] Ramachandran P, Zoph B and Le Q V 2017 Searching for activation functions (arXiv:1710.05941)

[73] Kingma D P and Ba J 2014 Adam: a method for stochastic optimization (arXiv:1412.6980)

[74] He K, Zhang X, Ren S and Sun J 2016 Deep residual learning for image recognition *Proc. IEEE Conference on Computer Vision and Pattern Recognition* pp 770–8

[75] Tetel’man M 1982 Lorentz group for two-dimensional integrable lattice systems *Soviet Journal of Experimental and Theoretical Physics* 55 306

[76] Hoffman J, Roberts D A and Yaida S 2019 Robust learning with Jacobian regularization (arXiv:1908.02729)

[77] Park S, Yun C, Lee J and Shin J 2020 Minimum width for universal approximation (arXiv:2006.08859)

[78] Choromanska A, Henaff M, Mathieu M, Arous G B and LeCun Y 2015 The loss surfaces of multilayer networks *Artificial Intelligence and Statistics* pp 192–204

[79] Dauphin Y N, Pascanu R, Gulcehre C, Cho K, Ganguli S and Bengio Y 2014 Identifying and attacking the saddle point problem in high-dimensional non-convex optimization *Advances in Neural Information Processing Systems* p 27

[80] Takhtadzhan L and Faddeev L D 1979 The quantum method of the inverse problem and the Heisenberg XYZ model *Russ. Math. Surv.* 34 11

[81] Lal S, Majumder S and Sobko E Drawing the map of integrable spin chains in progress

[82] Beisert N and Staudacher M 2003 The $N = 4$ SYM integrable super spin chain *Nucl. Phys. B* 670 439

[83] Borsato R, Sax O O, Sfondrini A, Stefanski B and Torrielli A 2013 The all-loop integrable spin-chain for strings on AdS3 × S3 × T4: the massive sector *J. High Energy Phys.* JHEP08(2013)043

[84] Majumder S, Sax O O, Stefański B and Torrielli A 2021 Protected states in AdS3 backgrounds from integrability *J. Phys. A: Math. Theor.* 54 415401

[85] Frolov S and Sfondrini A 2022 Mirror thermodynamic Bethe ansatz for AdS3/CFT2 *J. High Energy Phys.* JHEP03(2022)138

[86] Schmidt M and Lipson H 2009 Distilling free-form natural laws from experimental data *Science* 324 81

[87] Quessard R, Barrett T and Clements W 2020 Learning disentangled representations and group structure of dynamical environments *Advances in Neural Information Processing Systems* vol 33 p 19727

[88] Paulos M F and Zheng Z 2020 Bounding scattering of charged particles in 1 + 1 dimensions *J. High Energy Phys.* JHEP05(2020)145

[89] Zhang Y, Kauffman L H and Ge M-L 2005 Yang-Baxterizations, universal quantum gates and Hamiltonians *Quantum Inf. Proc. 4* 159

[90] Rowell E C, Zhang Y, Wu Y-S and Ge M-L 2010 Extraspecial two-Groups, generalized Yang-Baxter equations and braiding quantum gates *Quantum Inf. Comput.* 10 0685

[91] Padmanabhan P, Sugino F and Trancanelli D 2020 Quantum entanglement, supersymmetry and the generalized Yang-Baxter equation *Quantum Inf. Comput.* 20 37

[92] Padmanabhan P, Sugino F and Trancanelli D 2020 Braiding quantum gates from partition algebras *Quantum* 4 311

[93] Bromley J, Bentz J W, Bottou L, Guyon I, LeCun Y, Moore C, Säckinger E and Shah R 1993 Signature verification using a “siamese” time delay neural network *Int. J. Pattern Recognit. Artif. Intell.* 7 669

[94] Hadsell R, Chopra S and LeCun Y 2006 Dimensionality reduction by learning an invariant mapping *IEEE Computer Society Conf. on Computer Vision and Pattern Recognition (CVPR’06)* pp 1735–42

[95] Chechik G, Sharma V, Shalit U and Bengio S 2010 Large scale online learning of image similarity through ranking *J. Mach. Learn. Res.* 11 1109–1135

[96] Schroff F, Kalenichenko D and Philbin J Facenet: a unified embedding for face recognition and clustering 2015 *Proc. IEEE Conf. on Computer Vision and Pattern Recognition* pp 815–23