(A.62a) and the value of $w_{k,l}$ from equation (A.65) we obtain equation
(2.64f). This concludes the procedure of solution of the ${}_{1}D_{4}$
equation (1.3e).

###  A.4 $\boldsymbol{{}_{2}D_{4}}$ equation (1.3f)

To find the general solution of the ${}_{2}D_{4}$ equation (1.3f) we start
from the equation itself. Applying the general transformation (1.12) the
${}_{2}D_{4}$ equation we obtain the following system

$\displaystyle
v_{k,l}w_{k,l}+\delta_{2}w_{k,l}y_{k,l}+\delta_{1}w_{k,l}z_{k,l}+y_{k,l}z_{k,l}+\delta_{3}=0,$
(A.74a) $\displaystyle
v_{k,l+1}w_{k,l+1}+\delta_{2}y_{k,l}w_{k,l+1}+\delta_{1}z_{k,l}w_{k,l+1}+y_{k,l}z_{k,l}+\delta_{3}=0,$
(A.74b) $\displaystyle
w_{k,l}v_{k+1,l}+\delta_{2}w_{k,l}y_{k+1,l}+\delta_{1}w_{k,l}z_{k,l}+z_{k,l}y_{k+1,l}+\delta_{3}=0,$
(A.74c) $\displaystyle
w_{k,l+1}v_{k+1,l+1}+\delta_{2}y_{k+1,l}w_{k,l+1}+\delta_{1}z_{k,l}w_{k,l+1}+z_{k,l}y_{k+1,l}+\delta_{3}=0.$
(A.74d)

From the system (A.74) we have four different ways to compute $z_{k,l}$. This
means that we have some compatibility conditions. Indeed from (A.74a) and
(A.74c) we obtain the following equation for $v_{k+1,l}$

$\displaystyle
v_{k+1,l}=\frac{\delta_{1}w_{k,l}+y_{k+1,l}}{\delta_{1}w_{k,l}+y_{k,l}}v_{k,l}-\frac{(y_{k+1,l}-y_{k,l})(\delta_{2}w_{k,l}^{2}\delta_{1}-\delta_{3})}{(\delta_{1}w_{k,l}+y_{k,l})w_{k,l}},$
(A.75)

while from (A.74b) and (A.74d) we obtain the following equation for
$v_{k+1,l+1}$

$\displaystyle
v_{k+1,l+1}=\frac{\delta_{1}w_{k,l+1}+y_{k+1,l}}{\delta_{1}w_{k,l+1}+y_{k,l}}v_{k,l+1}-\frac{(y_{k+1,l}-y_{k,l})(\delta_{2}w_{k,l+1}^{2}\delta_{1}-\delta_{3})}{(\delta_{1}w_{k,l+1}+y_{k,l})w_{k,l+1}}.$
(A.76)

Equations (A.75) and (A.76) give rise to a compatibility condition between
$v_{k+1,l}$ and its shift in the $l$ direction $v_{k+1,l+1}$ which is given by

$\displaystyle\left[\begin{aligned}
y_{{k+1,l+1}}y_{{k,l}}+{\delta_{1}}(y_{{k,l}}w_{{k,l+1}}+y_{{k+1,l+1}}w_{{k,l+1}})\\\
{}-y_{{k+1,l}}y_{{k,l+1}}-{\delta_{1}}(y_{{k+1,l}}w_{{k,l+1}}-y_{{k,l+1}}w_{{k,l+1}})\end{aligned}\right]\big{(}v_{{k,l+1}}w_{{k,l+1}}+{\delta_{3}}-\delta_{1}\delta_{2}w_{{k,l+1}}^{2}\big{)}=0.$

Discarding the trivial solution

$\displaystyle v_{k,l}=\delta_{1}\delta_{2}w_{k,l}-\frac{\delta_{3}}{w_{k,l}}$

we obtain the following value for the field $w_{k,l}$

$\displaystyle
w_{k,l}=\frac{1}{\delta_{1}}\frac{y_{k+1,l-1}y_{k,l}-y_{k+1,l}y_{k,l-1}}{y_{k,l-1}+y_{k+1,l}-y_{k,l}-y_{k+1,l-1}},$
(A.77)

which makes (A.75) and (A.76) compatible. This gives us the first part of our
solution, as displayed in equation (2.68f). Inserting (A.77) into (A.75) we
are left to solve the following equation for $v_{k,l}$

$\displaystyle
v_{k+1,l}=\frac{y_{k+1,l}-y_{k+1,l-1}}{y_{k,l}-y_{k,l-1}}v_{k,l}-\frac{\delta_{1}\delta_{3}(y_{k,l-1}-y_{k+1,l-1})^{2}(y_{k,l}-y_{k,l-1})}{(y_{k+1,l}y_{k,l-1}-y_{k+1,l-1}y_{k,l})y_{k,l-1}^{2}}$
$\displaystyle\hphantom{v_{k+1,l}=}{}+\frac{y_{k,l-1}^{2}y_{k+1,l}\delta_{2}-\delta_{1}\delta_{3}y_{k+1,l}+\delta_{1}\delta_{3}y_{k+1,l-1}-\delta_{2}y_{k+1,l-1}y_{k,l-1}^{2}}{(y_{k,l}-y_{k,l-1})y_{k,l-1}}$
$\displaystyle\hphantom{v_{k+1,l}=}{}-\frac{\delta_{2}y_{k+1,l-1}y_{k,l-1}^{2}+\delta_{1}\delta_{3}y_{k+1,l-1}-2\delta_{3}\delta_{1}y_{k,l-1}}{y_{k,l-1}^{2}}.$

Making the transformation

$\displaystyle
v_{k,l}=(y_{k,l}-y_{k,l-1})V_{k,l}+\frac{\delta_{1}\delta_{3}}{y_{k,l}}-\delta_{2}y_{k,l}$
(A.78)

we obtain that $V_{k,l}$ satisfies the difference equation

$\displaystyle
V_{k+1,l}=V_{k,l}-\frac{\delta_{1}\delta_{3}(y_{k,l}-y_{k+1,l})^{2}}{y_{k,l}y_{k+1,l}(y_{k+1,l-1}y_{k,l}-y_{k+1,l}y_{k,l-1})}.$
(A.79)

To go further we need to specify the form of the field $y_{k,l}$. This can be
obtained from the Darboux integrability of the ${}_{2}D_{4}$ equation (1.3f).
From [26] we know that in the case of the ${}_{2}D_{4}$ equation (1.3f) we
have the following four-point, third-order $W_{1}$ first integral

$\displaystyle
W_{1}={F^{(+)}_{n}F^{(+)}_{m}}\alpha{\frac{\left[\begin{gathered}(u_{{n,m}}-u_{{n+2,m}}-{\delta_{1}}{\delta_{2}}u_{{n-1,m}})u_{n+1,m}^{2}\\\
{}+u_{{n+1,m}}u_{{n+2,m}}u_{{n-1,m}}+{\delta_{3}}u_{{n-1,m}}\end{gathered}\right]}{({\delta_{2}}u_{n+1,m}^{2}{\delta_{1}}-{\delta_{3}}-u_{{n,m}}u_{{n+1,m}})u_{{n-1,m}}}}$
$\displaystyle\hphantom{W_{1}=}{}-{F^{(+)}_{n}F^{(-)}_{m}}\alpha{\frac{u_{{n+2,m}}u_{{n-1,m}}+(-u_{{n+2,m}}+u_{{n,m}})u_{{n+1,m}}+{\delta_{3}}}{u_{{n-1,m}}u_{{n,m}}+{\delta_{3}}}}$
$\displaystyle\hphantom{W_{1}=}{}-{F^{(-)}_{n}F^{(+)}_{m}}\beta{\frac{(u_{{n+1,m}}-u_{{n-1,m}})(u_{{n+2,m}}-u_{{n,m}})u_{{n,m}}}{u_{{n+2,m}}({\delta_{2}}{\delta_{1}}{u_{{n,m}}}^{2}-u_{{n,m}}u_{{n+1,m}}-{\delta_{3}})}}$
$\displaystyle\hphantom{W_{1}=}{}+{F^{(-)}_{n}F^{(-)}_{m}}\beta{\frac{(u_{{n+1,m}}-u_{{n-1,m}})(u_{{n+2,m}}-u_{{n,m}})}{u_{{n+1,m}}u_{{n+2,m}}+{\delta_{3}}}}.$

This first integral implies the relation $W_{1}=\xi_{n}$ which is a third-
order, four-point ordinary difference equation in the $n$ direction depending
parametrically on $m$. In particular if we choose the case when $m=2l+1$ we
have the equation

$\displaystyle-{F^{(+)}_{n}}{\frac{u_{{n+2,2l+1}}u_{{n-1,2l+1}}-(u_{{n+2,2l+1}}-u_{{n,2l+1}})u_{{n+1,2l+1}}+{\delta_{3}}}{u_{{n-1,2l+1}}u_{{n,2l+1}}+{\delta_{3}}}}$
$\displaystyle\qquad{}+{F^{(-)}_{n}}{\frac{(u_{{n+1,2l+1}}-u_{{n-1,2l+1}})(u_{{n+2,2l+1}}-u_{{n,2l+1}})}{u_{{n+1,2l+1}}u_{{n+2,2l+1}}+{\delta_{3}}}}=\xi_{n}.$
(A.80)

where we have chosen without loss of generality $\alpha=\beta=1$. Using the
transformation (1.12b) then equation (A.80) is converted into the system

$\displaystyle(y_{k+1,l}-y_{k,l})z_{k,l}-y_{k+1,l}z_{k+1,l}-\delta_{3}=\xi_{2k}(z_{k+1,l}y_{k,l}+\delta_{3}),$
(A.81a)
$\displaystyle(y_{k+1,l}-y_{k,l})(z_{k+1,l}-z_{k,l})=\xi_{2k+1}(y_{k+1,l}z_{k+1,l}+\delta_{3}).$
(A.81b)

If we solve (A.81b) with respect to $z_{k+1,l}$ and substitute into (A.81a) we
obtain a _linear_ , second-order ordinary difference equation for $y_{k,l}$

$\displaystyle\xi_{2k-1}y_{k+1,l}-(1+\xi_{2k}-\xi_{2k}\xi_{2k-1})y_{k,l}+(1+\xi_{2k})y_{k-1,l}=0.$
(A.82)

We can find the solution of equation (A.82) exploiting the arbitrariness of
the functions $\xi_{2k}$ and $\xi_{2k+1}$ as in the case of the ${}_{1}D_{4}$
equation (1.3e). Let us introduce the field
$Y_{k,l}=a_{k}y_{k,l}+b_{k}y_{k-1,l}$ and assume that $Y_{k+1,l}-Y_{k,l}$
equals the left-hand side of (A.82). Then we must have

$\displaystyle\xi_{2k}=\frac{b_{k+1}-b_{k}+-a_{k}}{a_{k+1}},\qquad\xi_{2k-1}=-\frac{b_{k+1}-b_{k}+a_{k+1}-a_{k}}{b_{k}}.$

This implies that $y_{k,l}$ will solve the equation

$\displaystyle a_{k}y_{k,l}+b_{k}y_{k-1,l}=\alpha_{l}.$ (A.83)

If we define

$\displaystyle a_{k}=\frac{1}{c_{k}}\frac{1}{d_{k}-d_{k-1}},\qquad
b_{k}=-\frac{1}{c_{k-1}}\frac{1}{d_{k}-d_{k-1}},$

equation (A.83) is solved by

$\displaystyle y_{k,l}=c_{k}(\alpha_{l}d_{k}+\beta_{l}).$ (A.84)

This gives the second part of our solution displayed in equation (2.66a).

Inserting the value of $y_{k,l}$ from (A.84) into equation (A.79) we obtain

$\displaystyle
V_{k+1,l}=V_{k,l}+\frac{\delta_{3}\delta_{1}(c_{k+1}-c_{k})^{2}}{(d_{k+1}-d_{k})(\beta_{l}\alpha_{l-1}-\beta_{l-1}\alpha_{l})c_{k}^{2}c_{k+1}^{2}}$
$\displaystyle\hphantom{V_{k+1,l}=}{}-\frac{\alpha_{l}\delta_{1}\delta_{3}}{\beta_{l}\alpha_{l-1}-\beta_{l-1}\alpha_{l}}\left[\frac{1}{(\alpha_{l}d_{k+1}+\beta_{l})c_{k+1}^{2}}-\frac{1}{(\alpha_{l}d_{k}+\beta_{l})c_{k}^{2}}\right].$
(A.85)

Therefore we can represent the solution to this equation as

$\displaystyle
V_{k,l}=\gamma_{l}-\frac{\alpha_{l}\delta_{1}\delta_{3}}{(\alpha_{l}d_{k}+\beta_{l})c_{k}^{2}(\beta_{l}\alpha_{l-1}-\beta_{l-1}\alpha_{l})}+\frac{\delta_{3}\delta_{1}e_{k}}{\beta_{l}\alpha_{l-1}-\beta_{l-1}\alpha_{l}},$
(A.86)

where $e_{k}$ is given by the discrete integration

$\displaystyle
e_{k+1}-e_{k}=\frac{(c_{k+1}-c_{k})^{2}}{(d_{k+1}-d_{k})c_{k}^{2}c_{k+1}^{2}},$

which is just equation (2.67a). Inserting the value of $V_{k,l}$ from equation
(A.86) and the value of $y_{k,l}$ from equation (A.84) into equation (A.78) we
obtain equation (2.66b) which is the third part of our solution. Using this
value alongside with the value of $w_{k,l}$ from equation (A.65) we can
compute $z_{k,l}$ from the original system (A.74). We obtain a _single_
compatibility condition given by

$\displaystyle(\beta_{l}\alpha_{l+1}-\beta_{l+1}\alpha_{l})\gamma_{l+1}-(\beta_{l-1}\alpha_{l}-\beta_{l}\alpha_{l-1})\gamma_{l}=(\beta_{l}\alpha_{l+1}-\beta_{l+1}\alpha_{l})\delta_{2},$

which is just equation (2.67b). This concludes the procedure of solution of
the ${}_{2}D_{4}$ equation (1.3f).

###  A.5 $\boldsymbol{{}_{t}H_{3}^{\varepsilon}}$ equation (1.2c)

To solve the ${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) we start from the
equation itself. We apply the transformation (2.69) to write down the
${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) as the following system of two
coupled equations
$\displaystyle\alpha_{2}(p_{n,l}q_{n+1,l}+p_{n+1,l}q_{n,l})-p_{n,l}q_{n,l}-p_{n+1,l}q_{n+1,l}$
$\displaystyle\qquad{}-\alpha_{3}\big{(}\alpha_{2}^{2}-1\big{)}\left(\delta^{2}+\varepsilon^{2}\frac{q_{n,l}q_{n+1,l}}{\alpha_{3}^{2}\alpha_{2}}\right)=0,$
(A.87a)
$\displaystyle\alpha_{2}(q_{n,l}p_{n+1,l+1}+q_{n+1,l}p_{n,l+1})-q_{n,l}p_{n,l+1}-q_{n+1,l}p_{n+1,l+1}$
$\displaystyle\qquad{}-\alpha_{3}\big{(}\alpha_{2}^{2}-1\big{)}\left(\delta^{2}+\varepsilon^{2}\frac{q_{n,l}q_{n+1,l}}{\alpha_{3}^{2}\alpha_{2}}\right)=0.$
(A.87b)

As in the case of the ${}_{t}H^{\varepsilon}_{2}$ equation (1.2b) we have that
equation (A.87a) depends on $p_{n,l}$ and $p_{n+1,l}$ and that equation
(A.87b) depends on $p_{n,l+1}$ and $p_{n+1,l+1}$. So we can apply the
translation operator $T_{l}$ to (A.87a) to obtain two equations in terms of
$p_{n,l+1}$ and $p_{n+1,l+1}$

$\displaystyle\alpha_{2}(p_{n,l+1}q_{n+1,l+1}+p_{n+1,l+1}q_{n,l+1})-p_{n,l+1}q_{n,l+1}-p_{n+1,l+1}q_{n+1,l+1}$
$\displaystyle\qquad{}-\alpha_{3}\big{(}\alpha_{2}^{2}-1\big{)}\left(\delta^{2}+\varepsilon^{2}\frac{q_{n,l+1}q_{n+1,l+1}}{\alpha_{3}^{2}\alpha_{2}}\right)=0,$
(A.88a)
$\displaystyle\alpha_{2}(q_{n,l}p_{n+1,l+1}+q_{n+1,l}p_{n,l+1})-q_{n,l}p_{n,l+1}-q_{n+1,l}p_{n+1,l+1}$
$\displaystyle\qquad{}-\alpha_{3}\big{(}\alpha_{2}^{2}-1\big{)}\left(\delta^{2}+\varepsilon^{2}\frac{q_{n,l}q_{n+1,l}}{\alpha_{3}^{2}\alpha_{2}}\right)=0.$
(A.88b)

The system (A.88) is equivalent to the original system (A.87). Then since we
can assume $\alpha_{2},\alpha_{3}\neq 0$555If $\alpha_{2}=0$ or $\alpha_{3}=0$
in (A.87) we have that the system becomes trivially equivalent to $q_{n,l}=0$
and $p_{n,l}$ is left unspecified. Therefore we can discard this trivial case.
we can solve (A.88) with respect to $p_{n,l+1}$ and $p_{n+1,l+1}$:

$\displaystyle
p_{{n,l+1}}={\frac{\left[\begin{gathered}{\alpha_{2}}(q_{{n+1,l+1}}-q_{{n+1,l}})\big{(}{\varepsilon}^{2}q_{{n,l}}q_{{n,l+1}}+{\delta}^{2}{{\alpha_{3}}}^{2}\big{)}\\\
{}+{\delta}^{2}\alpha_{2}^{2}\alpha_{3}^{2}(q_{{n,l}}-q_{{n,l+1}})+{\varepsilon}^{2}q_{{n+1,l+1}}q_{{n+1,l}}(q_{{n,l}}-q_{{n,l+1}})\end{gathered}\right]}{(q_{{n+1,l+1}}q_{{n,l}}-q_{{n+1,l}}q_{{n,l+1}}){\alpha_{3}}{\alpha_{2}}}},$
(A.89c) $\displaystyle
p_{{n+1,l+1}}={\frac{\left[\begin{gathered}{\alpha_{2}}\big{(}{\varepsilon}^{2}q_{{n+1,l+1}}q_{{n+1,l}}+{\delta}^{2}{{\alpha_{3}}}^{2}\big{)}(q_{{n,l}}-q_{{n,l+1}})\\\
{}+{\delta}^{2}\alpha_{2}^{2}\alpha_{3}^{2}(q_{{n+1,l+1}}-q_{{n+1,l}})+{\varepsilon}^{2}q_{{n,l}}q_{{n,l+1}}(q_{{n+1,l+1}}-q_{{n+1,l}})\end{gathered}\right]}{(q_{{n+1,l+1}}q_{{n,l}}-q_{{n+1,l}}q_{{n,l+1}}){\alpha_{3}}{\alpha_{2}}}}.$
(A.89f)

We see that the right-hand sides of (A.89) are functions only of $q_{n,l}$,
$q_{n+1,l}$, $q_{n,l+1}$ and $q_{n+1,l+1}$ and are well defined as long as
$q_{n,l}$ is not a solution of equation (2.108), which is therefore a
_singular case_. Therefore at this point the procedure of solution bifurcates
into two cases. We treat them separately.

Singular case: $\boldsymbol{q_{n,l}}$ solve (2.108). If $q_{n,l}$ solves
equation (2.108) we first solve this equation with respect to $q_{n,l}$ and
then use the system (A.87) to specify $p_{n,l}$. Indeed equation (2.108) is a
trivial Darboux integrable equation, since it possesses the following two-
point, first-order first integrals

$\displaystyle W_{1}=\frac{q_{n+1,l}}{q_{n,l}},$ (A.90a) $\displaystyle
W_{2}=\frac{q_{n,l+1}}{q_{n,l}}.$ (A.90b)

As remarked in the introduction the existence of a two-point, first-order
first integral means that the equation is itself a first integral. Therefore
the equation (2.108) can be alternatively written as
$\left(T_{l}-\operatorname{Id}\right)W_{1}$ or
$\left(T_{n}-\operatorname{Id}\right)W_{2}$ with $W_{1}$ and $W_{2}$ given by
(A.90). From (A.90a) we obtain

$\displaystyle q_{n+1,l}=\xi_{n}q_{n,l}$ (A.91)

being $\xi_{n}$ an arbitrary function of its argument. Equation (A.91) is
casted into total difference form by defining $\xi_{n}=a_{n+1}/a_{n}$, with
$\alpha_{n}$ a new arbitrary function of its argument. Then we obtain that the
general solution of equation (2.108) is

$\displaystyle q_{n,l}=a_{n}\zeta_{l},$ (A.92)

where $\zeta_{l}$ is an arbitrary function of its argument.

Let us remark that equation (2.108) is the _logarithmic discrete wave
equation_ , since it can be mapped into the discrete wave equation (2.70)
exponentiating (2.70) and then taking $q_{n,l}\to e^{q_{n,l}}$, and it is a
discretization of the hyperbolic partial differential equation

$\displaystyle uu_{xt}-u_{x}u_{t}=0,$

which is obtained from the wave equation $v_{xt}=0$ using the transformation
$v=\log u$. This fact is worth to note since being the transformation
connecting (2.70) and (2.108) not bi-rational, integrability properties, in
this case linearization and Darboux integrability, are not _a priori_
preserved [19].

Substituting (A.92) into (A.88) we obtain the compatibility condition

$\displaystyle\zeta_{l+1}-\zeta_{l}=0,$ (A.93)

i.e., $\zeta_{l}=\zeta_{0}=\text{const}$ and the system (A.87) is now
consistent. With this we find the first piece of the general solution in this
case given by (2.113a). Therefore we are left with one equation for $p_{n,l}$,
e.g., (A.87a). Therefore inserting (A.92) with $\zeta_{l}=\zeta_{0}$ in
(A.87a) and solving with respect to $p_{n+1,l}$ we obtain

$\displaystyle
p_{n+1,l}-\frac{\alpha_{2}a_{n+1}-a_{n}}{a_{n+1}-\alpha_{2}a_{n}}p_{n,l}=\big{(}\alpha_{2}^{2}-1\big{)}\frac{\delta^{2}\alpha_{3}^{2}\alpha_{2}+\varepsilon^{2}a_{n}\zeta_{0}^{2}a_{n+1}}{\alpha_{3}\alpha_{2}\zeta_{0}(\alpha_{2}a_{n}-a_{n+1})}.$

Defining

$\displaystyle\frac{\alpha_{2}a_{n+1}-a_{n}}{a_{n+1}-\alpha_{2}a_{n}}=\frac{b_{n+1}}{b_{n}},$
(A.94)

we have that $p_{n,l}$ solves the equation

$\displaystyle\frac{p_{n+1,l}}{b_{n+1}}-\frac{p_{n,l}}{b_{n}}={\frac{{\delta}^{2}\alpha_{3}^{2}\alpha_{2}^{2}b_{{n}}-b_{{n+1}}\big{(}{\delta}^{2}\alpha_{3}^{2}+{\varepsilon}^{2}a_{{n}}^{2}\zeta_{{0}}^{2}\big{)}{\alpha_{2}}+{\varepsilon}^{2}a_{{n}}^{2}\zeta_{{0}}^{2}b_{{n}}}{b_{{n}}a_{{n}}\zeta_{{0}}{\alpha_{2}}{\alpha_{3}}b_{{n+1}}}}.$
(A.95)

Note that $b_{n}$ in (A.94) is given in terms of $a_{n}$ and $a_{n+1}$ through
discrete integration and it is the constraint given in equation (2.114a).
Equation (A.95) is solved by

$\displaystyle p_{n,l}=b_{n}(\beta_{l}+c_{n}),$

where $c_{n}$ is given by the discrete integration

$\displaystyle
c_{n+1}=c_{n}+{\frac{{\delta}^{2}\alpha_{3}^{2}\alpha_{2}^{2}b_{{n}}-b_{{n+1}}\big{(}{\delta}^{2}\alpha_{3}^{2}+{\varepsilon}^{2}a_{{n}}^{2}\zeta_{{0}}^{2}\big{)}{\alpha_{2}}+{\varepsilon}^{2}a_{{n}}^{2}\zeta_{{0}}^{2}b_{{n}}}{b_{{n}}a_{{n}}\zeta_{{0}}{\alpha_{2}}{\alpha_{3}}b_{{n+1}}}},$

i.e., as in equation (2.114b). This yields the solution of the
${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) when $q_{n,l}$ satisfy equation
(2.108).

General case: $\boldsymbol{q_{n,l}}$ do not solve (2.108). If the field
$q_{n,l}$ do not solve (2.108) we have that we can define $p_{n,l+1}$ and
$p_{n+1,l+1}$ as in (A.89c) and (A.89f) respectively. Furthermore these two
equations must be compatible. The compatibility condition is obtained applying
$T_{l}^{-1}$ to (A.89f) and imposing to the obtained expression to be equal to
(A.89c). We then find that $q_{n,l}$ must solve the following equation

$\displaystyle\delta^{2}\alpha_{2}^{2}\alpha_{3}^{2}\big{[}q_{{n+1,l+1}}q_{{n,l}}-q_{{n,l}}q_{{n-1,l+1}}+q_{{n,l+1}}(q_{{n-1,l}}-q_{{n+1,l}})\big{]}$
$\displaystyle\qquad{}-\alpha_{2}\big{(}{\varepsilon}^{2}q_{{n,l}}q_{{n,l+1}}+\delta^{2}\alpha_{3}^{2}\big{)}(q_{{n+1,l+1}}q_{{n-1,l}}-q_{{n-1,l+1}}q_{{n+1,l}})$
$\displaystyle\qquad{}+{\varepsilon}^{2}\big{[}q_{{n,l}}q_{{n+1,l+1}}(q_{{n-1,l}}-q_{{n+1,l}})-q_{{n,l+1}}q_{{n-1,l}}q_{{n+1,l}}\big{]}q_{{n-1,l+1}}$
$\displaystyle\qquad{}+{\varepsilon}^{2}q_{{n+1,l+1}}q_{{n+1,l}}q_{{n,l+1}}q_{{n-1,l}}=0.$
(A.96)

As in the case of the ${}_{t}H^{\varepsilon}_{2}$ equation (1.2b) the partial
difference equation for $q_{n,l}$ is not defined on the square quad graph of
Fig. 1, but it is defined on the six-point lattice shown in Fig. 3. Moreover
once equation (A.96) is solved we can use indifferently (A.89c) or (A.89f) to
obtain the value of the field $p_{n,l}$ since these two merely defines
$p_{n,l+1}$ in terms of $q_{n,l}$ and its shifts. Therefore if we find the
general solution of (A.96) the value of $p_{n,l}$ will follow. E.g., if we
solve (A.96) applying $T_{l}^{-1}$ to (A.89c) we will obtain (2.109d) which is
then the first part of the general solution. To find the solution of equation
(A.96) we turn to the first integrals. Like in the case of the
${}_{t}H^{\varepsilon}_{2}$ equation (1.2b) we will find an expression for
$q_{n,l}$ using the first integrals, and then we will insert it into (A.96) to
reduce the number of arbitrary functions to the right one. From [26] we know
that the ${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) possesses a four-point,
third-order integral in the $n$ direction

$\displaystyle
W_{1}=F^{(+)}_{m}\frac{(u_{{n+1,m}}-u_{{n-1,m}})(u_{{n+2,m}}-u_{{n,m}})}{\alpha_{2}^{4}{\varepsilon}^{2}{\delta}^{2}-\alpha_{2}^{3}u_{{n+1,m}}u_{{n,m}}+\alpha_{2}^{2}\big{(}u_{{n,m}}^{2}+u_{n+1,m}^{2}-2{\varepsilon}^{2}{\delta}^{2}\big{)}-{\alpha_{2}}u_{{n,m}}u_{{n+1,m}}+{\varepsilon}^{2}{\delta}^{2}}$
$\displaystyle\hphantom{W_{1}=}{}-F^{(-)}_{m}{\frac{(u_{{n+1,m}}-u_{{n-1,m}})(u_{{n+2,m}}-u_{{n,m}})}{{\alpha_{2}}(-u_{{n-1,m}}+{\alpha_{2}}u_{{n,m}})(-u_{{n+2,m}}+u_{{n+1,m}}{\alpha_{2}})}}.$
(A.97)

We consider the equation $W_{1}=\xi_{n}/\alpha_{2}$666The extra $\alpha_{2}$
is due to the arbitrariness of $\xi_{n}$ and is inserted in order to simplify
the formulas., where $W_{1}$ is given by (A.97), with $m=2l+1$

$\displaystyle\frac{(u_{n+1,2l+1}-u_{n-1,2l+1})(u_{n+2,2l+1}-u_{n,2l+1})}{(\alpha_{2}u_{n,2l+1}-u_{n-1,2l+1})(u_{n+2,2l+1}-\alpha_{2}u_{n+1,2l+1})}=\xi_{n}.$

Using the substitutions (2.69) we have

$\displaystyle\frac{(q_{n+1,l}-q_{n-1,l})(q_{n+2,l}-q_{n,l})}{(\alpha_{2}q_{n,l}-q_{n-1,l})(q_{n+2,l}-\alpha_{2}q_{n+1,l})}=\xi_{n}.$
(A.98)

This equation contains only $q_{n,l}$ and its shifts. By the transformation

$\displaystyle
Q_{n,l}=\frac{\alpha_{2}q_{n,l}-q_{n-1,l}}{q_{n+1,l}-q_{n-1,l}}$ (A.99)

equation (A.98) becomes

$\displaystyle Q_{n+1,l}+\frac{1}{\xi_{n}Q_{n,l}}=1,$ (A.100)

which is the same discrete Riccati equation as in (2.89). This means that the
solution of (A.100) is given by (2.95) with the appropriate definitions
(2.90), (2.92), (2.94). We can substitute into (A.99) the solution (2.95)

$\displaystyle\frac{q_{n+1,l}-q_{n-1,l}}{\alpha_{2}q_{n,l}-q_{n-1,l}}=\frac{(c_{n}+\zeta_{l})(c_{n+1}-c_{n-1})}{(c_{n+1}+\zeta_{l})(c_{n}-c_{n-1})}$

and we obtain an equation for $q_{n,l}$. Introducing

$\displaystyle P_{n,l}=(c_{n}+\zeta_{l})q_{n,l}$ (A.101)

we obtain that $P_{n,l}$ solves the equation

$\displaystyle
P_{n+1,l}-\alpha_{2}\frac{c_{n+1}-c_{n-1}}{c_{n}-c_{n-1}}P_{n,l}+\frac{c_{n+1}-c_{n}}{c_{n}-c_{n-1}}P_{n-1,l}=0.$
(A.102)

Using the transformation

$\displaystyle P_{n,l}=\frac{R_{n,l}}{R_{n-1,l}},$ (A.103)

we can cast equation (A.102) in discrete Riccati equation form

$\displaystyle
R_{n+1,l}+\frac{c_{n+1}-c_{n}}{c_{n}-c_{n-1}}\frac{1}{R_{n,l}}=\alpha_{2}\frac{c_{n+1}-c_{n-1}}{c_{n}-c_{n-1}}.$
(A.104)

Let $d_{n}$ be a particular solution of equation (A.104)

$\displaystyle
d_{n+1}+\frac{c_{n+1}-c_{n}}{c_{n}-c_{n-1}}\frac{1}{d_{n}}=\alpha_{2}\frac{c_{n+1}-c_{n-1}}{c_{n}-c_{n-1}}.$
(A.105)

Assuming $d_{n}$ as the new arbitrary function we can express $c_{n}$ as the
result of _two_ discrete integrations. Indeed introducing
$z_{n}=c_{n}-c_{n-1}$ in (A.105) we have

$\displaystyle\frac{z_{n+1}}{z_{n}}=\frac{(d_{n+1}-\alpha_{2})d_{n}}{\alpha_{2}d_{n}-1}.$
(A.106)

Equation (A.106) represents the first discrete integration, whereas the second
one is given by the definition

$\displaystyle c_{n}-c_{n-1}=z_{n}.$ (A.107)

Now we can linearize the discrete Riccati equation (A.104) by the
transformation

$\displaystyle R_{n,l}=d_{n}+\frac{1}{S_{n,l}}$ (A.108)

and we get the following _linear_ equation for $S_{n,l}$

$\displaystyle
S_{n+1,l}-\frac{d_{n}^{2}(c_{n}-c_{n-1})}{c_{n+1}-c_{n}}S_{n,l}=\frac{d_{n}(c_{n}-c_{n-1})}{c_{n+1}-c_{n}}.$
(A.109)

Defining

$\displaystyle d_{n}=\frac{e_{n}}{e_{n-1}},$ (A.110a) $\displaystyle
f_{n}-f_{n-1}=\frac{c_{n}-c_{n-1}}{e_{n}e_{n-1}},$ (A.110b)

the solution of (A.109) is

$\displaystyle S_{n,l}=\frac{(f_{n-1}+\beta_{l})e_{n-1}^{2}}{c_{n}-c_{n-1}}.$
(A.111)

Here we have the final form of the constraint $f_{n}$, which is the same as
the one given in (2.111a). Inserting (A.111) and (A.110) into (A.108) we
obtain

$\displaystyle
R_{n,l}=\frac{e_{n}(f_{n}+\beta_{l})}{e_{n-1}(f_{n-1}+\beta_{l})}.$ (A.112)

Inserting the definition of $R_{n,l}$ (A.103) into (A.112) we obtain

$\displaystyle\frac{P_{n,l}}{e_{n}(f_{n}+\beta_{l})}=\frac{P_{n-1,l}}{e_{n-1}(f_{n-1}+\beta_{l})},$

i.e.,

$\displaystyle P_{n,l}=\gamma_{l}e_{n}(f_{n}+\beta_{l}).$ (A.113)

Introducing (A.113) into (A.101) we obtain

$\displaystyle
q_{n,l}=\frac{\gamma_{l}e_{n}(f_{n}+\beta_{l})}{c_{n}+\zeta_{l}},$ (A.114)

where $f_{n}$ is defined by (A.110b), and $c_{n}$ is given by (A.106) and
(A.107), i.e., $c_{n}$ is the solution of the equation

$\displaystyle\frac{c_{{n+1}}-c_{{n}}}{c_{{n}}-c_{{n-1}}}={\frac{e_{{n+1}}-{\alpha_{2}}e_{{n}}}{{\alpha_{2}}e_{{n}}-e_{{n-1}}}},$
(A.115)

and $e_{n}$ is an arbitrary function, i.e., we have that $c_{n}$ must solve
equation (2.110).

Formally equation (A.114) has the form of the solution presented in formula
(2.109a), but it depends on three arbitrary functions in the $l$ direction,
namely $\zeta_{l}$, $\beta_{l}$ and $\gamma_{l}$. Therefore there must be a
constraint between these functions. This constraint can be obtained by
plugging (A.114) into (A.96), but here we have another bifurcation depending
on the value of parameter $\delta$. Indeed it is easy to see that we must
distinguish the cases when $\delta\neq 0$ and when $\delta=0$.

Case $\boldsymbol{\delta\neq 0}$. Inserting (A.114) into (A.96) if $\delta\neq
0$ factorizing the $n$ dependent part away we are left with

$\displaystyle\zeta_{l+1}-\zeta_{l}=\frac{\varepsilon^{2}}{\alpha_{2}\delta^{2}\alpha_{3}^{2}}\gamma_{l+1}\gamma_{l}(\beta_{l+1}-\beta_{l}).$

This equation tells us that the function $\zeta_{l}$ can be expressed after a
discrete integration in terms of the two arbitrary functions $\beta_{l}$ and
$\gamma_{l}$. This condition is just (2.111b). This yields the general
solution of the ${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) when $\delta\neq
0$ and the field $q_{n,l}$ do not satisfy equation (2.108).

Case $\boldsymbol{\delta=0}$. Inserting (A.114) into (A.96) if $\delta\neq 0$
factorizing the $n$ dependent part the compatibility condition is just
$\beta_{l+1}-\beta_{l}=0$, i.e., $\beta_{l}=\beta_{0}=\text{const}$. It is
easy to check that the obtained value of $q_{n,l}$ through formula (A.114) is
consistent with the substitution of $\delta=0$ in (A.87) and therefore that in
the case $\delta=0$ the value of $q_{n,l}$ is given by

$\displaystyle
q_{n,l}=\frac{\gamma_{l}e_{n}(f_{n}+\beta_{0})}{c_{n}+\zeta_{l}},$

where the functions $c_{n}$ and $f_{n}$ are defined implicitly and can be
found by discrete integration from (A.110b) and (A.115) respectively. This is
just equation (2.112a). Since equation (2.109d) is not singular if $\delta=0$
we obtain that in this case the general solution for the field $p_{n,l}$ is
given by substituting $\delta=0$ into (2.109d), i.e., just by equation
(2.112b) where $q_{n,l}$ is simply given by (2.112a). This yields the general
solution of the ${}_{t}H^{\varepsilon}_{3}$ equation (1.2c) in the case when
$\delta=0$.

### Acknowledgements

We would like to thank Professor Decio Levi for the many interesting and
fruitful discussion during the preparation of this paper. We are also grateful
to the anonymous referees whose comments greatly helped us in improving the
paper.

GG has been supported by INFN IS-CSN4 _Mathematical Methods of Nonlinear
Physics_ and by the Australian Research Council through an Australian Laureate
Fellowship grant FL120100094.

## References

  * [1] Adler V.E., Bobenko A.I., Suris Yu.B., Classification of integrable equations on quad-graphs. The consistency approach, Comm. Math. Phys. 233 (2003), 513–543, nlin.SI/0202024.
  * [2] Adler V.E., Bobenko A.I., Suris Yu.B., Discrete nonlinear hyperbolic equations: classification of integrable cases, Funct. Anal. Appl. 43 (2009), 3–17, arXiv:0705.1663.
  * [3] Adler V.E., Startsev S.Ya., On discrete analogues of the Liouville equation, Theoret. and Math. Phys. 121 (1999), 1484–1495, solv-int/9902016.
  * [4] Bellon M.P., Viallet C.M., Algebraic entropy, Comm. Math. Phys. 204 (1999), 425–437, chao-dyn/9805006.
  * [5] Bobenko A.I., Suris Yu.B., Integrable systems on quad-graphs, Int. Math. Res. Not. 2002 (2002), 573–611, nlin.SI/0110004.
  * [6] Bobenko A.I., Suris Yu.B., Discrete differential geometry. Integrable structure, Graduate Studies in Mathematics, Vol. 98, Amer. Math. Soc., Providence, RI, 2008.
  * [7] Boll R., Classification of 3D consistent quad-equations, J. Nonlinear Math. Phys. 18 (2011), 337–365, arXiv:1009.4007.
  * [8] Boll R., Corrigendum: Classification of 3D consistent quad-equations, J. Nonlinear Math. Phys. 19 (2012), 1292001, 3 pages.
  * [9] Boll R., Classification and Lagrangian structure of 3D consistent quad-equations, Ph.D. Thesis, Technische Universität Berlin, 2012.
  * [10] Bridgman T., Hereman W., Quispel G.R.W., van der Kamp P.H., Symbolic computation of Lax pairs of partial difference equations using consistency around the cube, Found. Comput. Math. 13 (2013), 517–544, arXiv:1308.5473.
  * [11] Butler S., Hay M., Simple identification of fake Lax pairs, arXiv:1311.2406.
  * [12] Butler S., Hay M., Two definitions of fake Lax pairs, AIP Conf. Proc. 1648 (2015), 180006, 5 pages.
  * [13] Calogero F., Degasperis A., Spectral transform and solitons. Vol. I. Tools to solve and investigate nonlinear evolution equations, Studies in Mathematics and its Applications, Vol. 13, North-Holland Publishing Co., Amsterdam – New York, 1982.
  * [14] Calogero F., Nucci M.C., Lax pairs galore, J. Math. Phys. 32 (1991), 72–74.
  * [15] Doliwa A., Santini P.M., Multidimensional quadrilateral lattices are integrable, Phys. Lett. A 233 (1997), 365–372, solv-int/9612007.
  * [16] Falqui G., Viallet C.M., Singularity, complexity, and quasi-integrability of rational mappings, Comm. Math. Phys. 154 (1993), 111–125, hep-th/9212105.
  * [17] Garifullin R.N., Yamilov R.I., Generalized symmetry classification of discrete equations of a class depending on twelve parameters, J. Phys. A: Math. Theor. 45 (2012), 345205, 23 pages, arXiv:1203.4369.
  * [18] Garifullin R.N., Yamilov R.I., Integrable discrete nonautonomous quad-equations as Bäcklund auto-transformations for known Volterra and Toda type semidiscrete equations, J. Phys. Conf. Ser. 621 (2015), 012005, 18 pages, arXiv:1405.1835.
  * [19] Grammaticos B., Ramani A., Viallet C.M., Solvable chaos, Phys. Lett. A 336 (2005), 152–158, math-ph/0409081.
  * [20] Gubbiotti G., Scimiterna C., Reconstructing a lattice equation: a non-autonomous approach to the Hietarinta equation, SIGMA 14 (2018), 004, 21 pages, arXiv:1705.00298.
  * [21] Gubbiotti G., Scimiterna C., Levi D., Algebraic entropy, symmetries and linearization of quad equations consistent on the cube, J. Nonlinear Math. Phys. 23 (2016), 507–543, arXiv:1603.07930.
  * [22] Gubbiotti G., Scimiterna C., Levi D., Linearizability and a fake Lax pair for a nonlinear nonautonomous quad-graph equation consistent around the cube, Theoret. and Math. Phys. 189 (2016), 1459––1471.
  * [23] Gubbiotti G., Scimiterna C., Levi D., On partial differential and difference equations with symmetries depending on arbitrary functions, Acta Polytechnica 56 (2016), 193–201, arXiv:1512.01967.
  * [24] Gubbiotti G., Scimiterna C., Levi D., The non-autonomous YdKN equation and generalized symmetries of Boll equations, J. Math. Phys. 58 (2017), 053507, 18 pages, arXiv:1510.07175.
  * [25] Gubbiotti G., Scimiterna C., Levi D., A two-periodic generalization of the $Q_{\rm V}$ equation, J. Integrable Syst. 2 (2017), xyx004, 13 pages.
  * [26] Gubbiotti G., Yamilov R.I., Darboux integrability of trapezoidal $H^{4}$ and $H^{4}$ families of lattice equations I: First integrals, J. Phys. A: Math. Theor. 50 (2017), 345205, 26 pages, arXiv:1608.03506.
  * [27] Habibullin I.T., Characteristic algebras of fully discrete hyperbolic type equations, SIGMA 1 (2005), 023, 9 pages, nlin.SI/0506027.
  * [28] Hay M., A completeness study on discrete, $2\times 2$ Lax pairs, J. Math. Phys. 50 (2009), 103516, 29 pages, arXiv:0806.3940.
  * [29] Hay M., A completeness study on certain $2\times 2$ Lax pairs including zero terms, SIGMA 7 (2011), 089, 12 pages, arXiv:1104.0084.
  * [30] Hietarinta J., A new two-dimensional lattice model that is ‘consistent around a cube’, J. Phys. A: Math. Gen. 37 (2004), L67–L73, nlin.SI/0311034.
  * [31] Hietarinta J., Searching for CAC-maps, J. Nonlinear Math. Phys. 12 (2005), suppl. 2, 223–230.
  * [32] Hietarinta J., Definitions and predictions of integrability for difference equations, in Symmetries and Integrability of Difference Equations (Beijing, 2009), London Math. Soc. Lecture Note Ser., Vol. 381, Editors D. Levi, P. Olver, Z. Thomova, P. Winternitz, Cambridge University Press, Cambridge, 2011, 83–114.
  * [33] Hietarinta J., Joshi N., Nijhoff F.W., Discrete systems and integrability, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2016.
  * [34] Hietarinta J., Viallet C.M., Searching for integrable lattice maps using factorization, J. Phys. A: Math. Theor. 40 (2007), 12629–12643, arXiv:0705.1903.
  * [35] Hietarinta J., Viallet C.M., Weak Lax pairs for lattice equations, Nonlinearity 25 (2012), 1955–1966, arXiv:1105.3329.
  * [36] Liouville J., Sur l’equation aux différences partielles $\partial^{2}\log\lambda/\partial u\partial v\pm\lambda/(aa^{2})=0$, J. Math. Pures Appl. 18 (1853), 71–72.
  * [37] Nijhoff F.W., Lax pair for the Adler (lattice Krichever–Novikov) system, Phys. Lett. A 297 (2002), 49–58, nlin.SI/0110027.
  * [38] Nijhoff F.W., Walker A.J., The discrete and continuous Painlevé VI hierarchy and the Garnier systems, Glasg. Math. J. 43A (2001), 109–123, nlin.SI/0001054.
  * [39] Papageorgiou V.G., Nijhoff F.W., Capel H.W., Integrable mappings and nonlinear integrable lattice equations, Phys. Lett. A 147 (1990), 106–114.
  * [40] Quispel G.R.W., Capel H.W., Papageorgiou V.G., Nijhoff F.W., Integrable mappings derived from soliton equations, Phys. A 173 (1991), 243–266.
  * [41] Roberts J.A.G., Tran D.T., Algebraic entropy of (integrable) lattice equations and their reductions, arXiv:1703.01069.
  * [42] Viallet C.M., Algebraic entropy for lattice equations, math-ph/0609043.
  * [43] Viallet C.M., Integrable lattice maps: $Q_{\rm V}$, a rational version of $Q_{4}$, Glasg. Math. J. 51 (2009), 157–163, arXiv:0802.0294.
  * [44] Xenitidis P.D., Papageorgiou V.G., Symmetries and integrability of discrete equations defined on a black-white lattice, J. Phys. A: Math. Theor. 42 (2009), 454025, 13 pages, arXiv:0903.3152.
  * [45] Yamilov R., Symmetries as integrability criteria for differential difference equations, J. Phys. A: Math. Gen. 39 (2006), R541–R623.
  * [46] Zhiber A.V., Sokolov V.V., Exactly integrable hyperbolic equations of Liouville type, Russian Math. Surveys 56 (2001), no. 1, 61–101.


ViT-Base that operates on $14\times 14$ tokens, same as in Appendix A.3.5. For
Peptides-func and Peptides-struct, we use the following as the base model
$f_{\theta}:(\mathbf{A},\mathbf{X},\mathbf{E})\mapsto\mathbf{Y}$:

$\displaystyle
f_{\theta}(\mathbf{A},\mathbf{X},\mathbf{E})=\textnormal{detokenize}_{\texttt{[cls]}}\left(\textnormal{transformer}\left(\textnormal{tokenize}_{\textnormal{2D}}(\mathbf{A},\mathbf{E})+\textnormal{tokenize}_{\textnormal{1D}}(\mathbf{X})\right)\right),$
(78)

where $\textnormal{tokenize}_{\textnormal{2D}}:\mathbb{R}^{444\times
444\times(1+7)}\to\mathbb{R}^{(14\times 14)\times 768}$ is 2D convolution with
kernel size and stride 32,
$\textnormal{tokenize}_{\textnormal{1D}}:\mathbb{R}^{444\times
64}\to\mathbb{R}^{196\times 768}$ is 1D convolution with kernel size and
stride 3, and $\textnormal{detokenize}_{\texttt{[cls]}}$ performs linear
projection of the global [cls] token [26] to the target dimensionality. For
PCQM-Contact, we use the following as base model
$f_{\theta}:(\mathbf{A},\mathbf{X},\mathbf{E})\mapsto\mathbf{Y}$:

$\displaystyle
f_{\theta}(\mathbf{A},\mathbf{X},\mathbf{E})=\textnormal{detokenize}_{\textnormal{2D}}\left(\textnormal{transformer}\left(\textnormal{tokenize}_{\textnormal{2D}}(\mathbf{A},\mathbf{E})+\textnormal{tokenize}_{\textnormal{1D}}(\mathbf{X})\right)\right),$
(79)

where $\textnormal{tokenize}_{\textnormal{2D}}:\mathbb{R}^{53\times
53\times(1+6)}\to\mathbb{R}^{(14\times 14)\times 768}$ is 2D convolution with
kernel size and stride 4,
$\textnormal{tokenize}_{\textnormal{1D}}:\mathbb{R}^{53\times
64}\to\mathbb{R}^{196\times 768}$ is 1D convolution with kernel size and
stride 1, and $\textnormal{detokenize}_{\textnormal{2D}}:\mathbb{R}^{(14\times
14)\times 768}\to\mathbb{R}^{53\times 53\times 1}$ is transposed 2D
convolution with kernel size and stride 4.

##### Training

We train our models with cross-entropy for classification and L1 loss for
regression using AdamW [55] optimizer with batch size 128, learning rate 1e-5
except for PCQM-Contact where we use 5e-5, and weight decay 0.1 except for
PCQM-Contact where we use 0.01. We train the models for 50k steps under
learning rate warm-up for 5k steps then linear decay to 0 with early stopping
based on validation performance, which usually takes less than 12 hours on 8
RTX 3090 GPUs with 24GB using PyTorch Lightning [29]. For the equivariant
distribution $p_{\omega}$, we use noise scale $\eta=1$, and use dropout with
probability 0.1 except for PCQM-Contact where we do not use dropout. We
further employ patch dropout regularization [54, 53] for Peptides-func/struct,
with fixed probability 0.5 for Peptides-func and probability sampled uniformly
from $(0,0.5)$ at each step for Peptides-struct. For the sampling-based
estimation, we use 1 sample for training, 10 samples for validation, and 100
samples for testing, except for PCQM-Contact where we always use 10 samples.

###  A.4 Supplementary Experiments (Continued from Section 3)

In this section, we present additional experimental results that supplement
the experiments in Section 3 but could not be included in the main text due to
space constraints.

####  A.4.1 Graph Isomorphism Learning (Section 3.1)

Table 7:  Supplementary results for $\textnormal{S}_{n}$ invariant graph
separation with $\textnormal{S}_{n}$ symmetrized GIN-ID base function.
Baseline scores for GIN-ID-GA and GIN-ID-FA are taken from [76].

method arch. sym. GRAPH8c $\downarrow$ EXP $\downarrow$ EXP-classify
$\uparrow$ GIN-ID-GA - $\textnormal{S}_{n}$ 0 0 50% GIN-ID-FA -
$\textnormal{S}_{n}$ 0 0 100% GIN-ID-Canonical. - $\textnormal{S}_{n}$ 0 0 84%
GIN-ID-PS (Ours) - $\textnormal{S}_{n}$ 0 0 100%

Table 8:  Supplementary results for
$\textnormal{S}_{n}\times\textnormal{E}(3)$ equivariant $n$-body with
$\textnormal{E}(3)$ symmetrized GNN base function. Baseline scores for GNN-FA
and GNN-Canonical. are from [76] and [41], respectively.

method arch. sym. Position MSE $\downarrow$ GNN-FA $\textnormal{S}_{n}$
$\textnormal{E}(3)$ 0.0057 GNN-Canonical. $\textnormal{S}_{n}$
$\textnormal{E}(3)$ 0.0043 GNN-Canonical. (Reproduced) $\textnormal{S}_{n}$
$\textnormal{E}(3)$ 0.00457 GNN-GA $\textnormal{S}_{n}$ $\textnormal{E}(3)$
0.00408 $\pm$ 0.00002 GNN-PS (Ours) $\textnormal{S}_{n}$ $\textnormal{E}(3)$
0.00386 $\pm$ 0.00001

In our experiments on graph isomorphism learning in Section 3.1, we mainly
experimented for $\textnormal{S}_{n}$ symmetrization of an MLP. Here, we
provide supplementary results on $\textnormal{S}_{n}$ symmetrization of a GIN
base model with node identifiers, following [76]. The results can be found in
Table LABEL:table:exp_more. In accordance with Section 3.1, our approach
successfully performs $\textnormal{S}_{n}$ symmetrization of GIN-ID.

####  A.4.2 Particle Dynamics Learning (Section 3.2)

In our experiments on $n$-body dataset in Section 3.2, we experimented for
$\textnormal{S}_{n}\times\textnormal{E}(3)$ symmetrization using a 1D sequence
transformer architecture which has $2.3\times$ parameters compared to
baselines. To provide parameter-matched comparison against baselines in
literature, we apply our approach for $\textnormal{E}(3)$ symmetrization of
$\textnormal{S}_{n}$ equivariant GNN base model that is widely used in
literature [76, 41]. We faithfully follow [76, 41] on the experimental setups
including training hyperparameters and the configuration of GNN base model,
and only add $\textnormal{E}(3)$ equivariant distribution
$p_{\omega}(g|\mathbf{x}-\bar{\mathbf{x}}\boldsymbol{1}^{\top})$, _i.e._ ,
$q_{\omega}:(\mathbf{x}-\bar{\mathbf{x}}\boldsymbol{1}^{\top},\boldsymbol{\epsilon})\mapsto\mathbf{Q}_{g}$
by utilizing the 2-layer Vector Neurons architecture described in Appendix
A.3.2 using only its $\textnormal{O}(3)$ prediction head. We use 20 samples
for training and testing. The results can be found in Table
LABEL:table:nbody_more. In accordance with the results in Section 3.2, our
approach outperforms other symmetrization approaches and achieves a new state-
of-the-art of 0.00386 MSE.

####  A.4.3 Effect of Sample Size on Training and Inference

Figure 4: Test accuracy of MLP $f_{\theta}$ symmetrized by equivariant
distribution $p_{\omega}(g|\mathbf{x})$ trained on EXP-classify dataset across
a range of training sample sizes. Inference sample size is set to 10.

Figure 5: Row- and column-wise entropy of aggregated permutation matrices
$\mathbf{P}_{g}\sim p_{\omega}(g|\mathbf{x})$ after trained on EXP-classify
dataset across a range of training sample sizes. Dashed line indicates entropy
measured with random permutation matrices from $\textnormal{Unif}(G)$. Figure
6: Variance of estimation of MLP $f_{\theta}$ symmetrized by equivariant
distribution $p_{\omega}(g|\mathbf{x})$ and trained on EXP-classify dataset
for a range of training and inference sample sizes.

In this section, we provide additional analysis on how the sample size for
estimation of symmetrized function (Eq. (4)) affects training and inference.
We use the experimental setup of EXP-classify (Section 3.1;
$\textnormal{S}_{n}$ invariance) and analyze the behavior of MLP-PS with
identical initialization and hyperparameters, only controlling sample sizes
$\in\\{1,2,5,10,20,50\\}$ for training. Specifically, we analyze (1) variance
of permutation matrices $\mathbf{P}_{g}\sim p_{\omega}(g|\mathbf{x})$ measured
indirectly by the entropy of their aggregation
$\bar{\mathbf{P}}=\sum{\mathbf{P}_{g}}/N$ as in Section 3.1, (2) sample
variance of the unbiased estimator $g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x})$
of the symmetrized function $\phi_{\theta,\omega}(\mathbf{x})$ as in Eq. (4),
and (3) sample mean and variance of the estimated task loss
$\mathcal{L}(\mathbf{y},g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x}))$ where
$\mathcal{L}$ is binary cross entropy. All measurements are repeated 100 times
and averaged over the inputs and labels $(\mathbf{x},\mathbf{y})$ of the
validation dataset.

Observations are as follows. First, models trained with smaller sample sizes
need more iterations to converge, but after sufficient training (2000 epochs),
all achieve $>95\%$ test accuracy when evaluated with 10 samples (Figure 5).
Second, models trained with smaller sample sizes tend to be more sample
efficient, _i.e._ , tend to perform a lower variance estimation. Their
distribution $p_{\omega}(g|\mathbf{x})$ tend to learn more low-variance
permutations (Figure 5), and the models tend to learn low-variance estimation
of output and loss (left and center panels of Figure 6). This indicates that
small sample size may serve as a regularizer that encourages lower variance of
the estimator. However, this regularization effect is not always beneficial in
terms of task loss (right panel of Figure 6), as training sample size 1
achieves a poor task loss for all sample sizes presumably due to the
optimization challenge caused by over-regularization. In other words, the
sample size for training introduces a tradeoff; a small sample size takes more
training iterations to converge, but serves as a regularizer that encourages
lower variance of the estimator and thus a better inference time sample
efficiency. On the other hand, larger sample sizes for inference consistently
benefits all models (Figure 6).

Interestingly, this observed tendency is consistent with the theoretical
claims in literature [69, 70] on the sampling based training of symmetrized
models, which we reprise here. When training the symmetrized model
$\phi_{\theta,\omega}(\mathbf{x})$ in Eq. (4), we cannot directly observe
$\phi_{\theta,\omega}(\mathbf{x})$, but observe samples of its unbiased
estimator $g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x})$. Thus, it can be
questionable what objective we are actually optimizing during the sampling-
based training. Based on [69, 70], it turns out that minimizing a convex loss
function $\mathcal{L}$ on the estimated output $g\cdot
f_{\theta}(g^{-1}\cdot\mathbf{x})$ is equivalent to minimizing an upper bound
to the true objective on the symmetrized output
$\phi_{\theta,\omega}(\mathbf{x})$. This is because our estimation is no
longer unbiased when computing loss, as we have the following from Jensen’s
inequality:

$\displaystyle\mathbb{E}_{p_{\omega}(g|\mathbf{x})}[\mathcal{L}(\mathbf{y},g\cdot
f_{\theta}(g^{-1}\cdot\mathbf{x}))]\geq\mathcal{L}(\mathbf{y},\mathbb{E}_{p_{\omega}(g|\mathbf{x})}[g\cdot
f_{\theta}(g^{-1}\cdot\mathbf{x})])=\mathcal{L}(\mathbf{y},\phi_{\theta,\omega}(\mathbf{x})).$
(80)

That is, minimizing the sampling-based loss is minimizing an upper-bound
surrogate to the true objective. It has been claimed that optimizing this
upper bound has an implicit low-variance regularization effect [69, 70], which
is consistent with our observations. This also roughly explains why our
distribution $p_{\omega}(g|\mathbf{x})$ does not collapse to uniform
distribution although we do not impose any low-variance regularization
explicitly; training to directly minimize the task loss with samples
implicitly nudges the distribution towards low-variance solutions.

####  A.4.4 Additional Comparison to Group Averaging

Figure 7: Sample variance of output $g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x})$
(left) and loss $\mathcal{L}(\mathbf{y},g\cdot
f_{\theta}(g^{-1}\cdot\mathbf{x}))$ (right) of an identical MLP $f_{\theta}$
symmetrized by equivariant distribution (PS) and uniform distribution (GA).

Figure 8: Norm of full gradient over training epochs with respect to the
parameters of an identically initialized MLP symmetrized by equivariant
distribution (PS) and uniform distribution (GA).

In this section, we provide additional analysis of our approach in comparison
to sampling-based group averaging in terms of sample variance and convergence.
We use the experimental setup of EXP-classify (Section 3.1;
$\textnormal{S}_{n}$ invariance) and experiment with MLP-PS and MLP-GA.

We first analyze whether using the equivariant distribution
$p_{\omega}(g|\mathbf{x})$ for symmetrization offers a lower variance
estimation, _i.e._ , a better sample efficiency, compared to group averaging
with $\textnormal{Unif}(G)$. This supplements the results in Figure 2 that
$p_{\omega}(g|\mathbf{x})$ learns to produce low-variance permutations
compared to $\textnormal{Unif}(G)$. Specifically, we fix a randomly
initialized MLP $f_{\theta}$ and symmetrize it using our approach and group
averaging. We then measure (1) the sample variance of the unbiased estimator
$g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x})$ of the symmetrized function (Eq.
(1) and Eq. (4)), and (2) the sample variance of the estimated task loss
$\mathcal{L}(\mathbf{y},g\cdot f_{\theta}(g^{-1}\cdot\mathbf{x}))$ where
$\mathcal{L}$ is binary cross entropy. All measurements are repeated 100 times
and averaged over the inputs and labels $(\mathbf{x},\mathbf{y})$ of the
validation dataset. The results are in Figure 8, showing that symmetrization
with equivariant distribution $p_{\omega}(g|\mathbf{x})$ consistently offers a
lower variance estimation than group averaging across inference sample sizes.

In addition, we analyze whether the equivariant distribution
$p_{\omega}(g|\mathbf{x})$ for symmetrization offers more stable gradients for
the base function $f_{\theta}$ during training compared to group averaging, as
conjectured in Section 2. For this, we fix a randomly initialized MLP
$f_{\theta}$ and symmetrize it using our approach and group averaging. For
every few training epochs, we measure the full gradient of the task loss over
the entire training dataset with respect to the parameters of the base MLP
$f_{\theta}$. This averages out the variance from individual data points and
provides the net direction of the gradient on the base function offered by
$p_{\omega}(g|\mathbf{x})$ or $\textnormal{Unif}(G)$. The results are in
Figure 8, showing that that symmetrization with equivariant distribution
$p_{\omega}(g|\mathbf{x})$ offers a consistently larger magnitude of the net
gradient, while group averaging with $\textnormal{Unif}(G)$ leads to near-zero
net gradients. This indicates, for $\textnormal{Unif}(G)$, the gradients from
each training data instances are oriented in a largely divergent manner and
therefore the training signal is collectively not very informative, while
using $p_{\omega}(g|\mathbf{x})$ for symmetrization leads to more consistent
gradient across training data instances, _i.e._ , it offers a more stable
training signal.

####  A.4.5 Additional Comparison to Canonicalization

(a) Inner-symmetric nodes

(b) MLP-PS (Ours)

(c) MLP-Canonical.

Figure 9: A graph $\mathbf{x}$ with stabilizer subgroup
$G_{\mathbf{x}}\cong\textnormal{S}_{6}$.

(d) Inner-symmetric nodes

(e) MLP-PS (Ours)

(f) MLP-Canonical.

Figure 10: A graph $\mathbf{x}$ with stabilizer subgroup
$G_{\mathbf{x}}\cong\textnormal{S}_{3}\times\textnormal{S}_{2}$.

(a) Inner-symmetric nodes

(b) MLP-PS (Ours)

(c) MLP-Canonical.

Figure 11: A graph $\mathbf{x}$ with stabilizer subgroup
$G_{\mathbf{x}}\cong\textnormal{D}_{4}$.

In this section, we provide additional analysis of our approach in comparison
to canonicalization [41] that uses a single group element $g$ from an
equivariant canonicalizer $C_{\omega}:\mathbf{x}\mapsto\rho(g)$. The main
claim is that there always exist certain inputs that canonicalization fails to
guarantee exact $G$ equivariance, while our approach guarantees equivariance
for all inputs in expectation as in Theorem 1.

More specifically, let us recall the definition of $G$ equivariant
canonicalizer from [41]. A canonicalizer $C_{\omega}$ is $G$ equivariant if
$C_{\omega}(g\cdot\mathbf{x})=\rho(g)C_{\omega}(\mathbf{x})$ for all $g\in G$
and $\mathbf{x}\in\mathcal{X}$. Consider an input $\mathbf{x}$ which has a
non-trivial stabilizer $G_{\mathbf{x}}={h\in G|h\cdot\mathbf{x}=\mathbf{x}}$,
_i.e._ , has inner symmetries. It can be shown that equivariant canonicalizers
are ill-defined for these inputs. To see this, let $g_{1}=gh_{1}$ and
$g_{2}=gh_{2}$ for some $g\in G$ and any $h_{1},h_{2}\in G_{\mathbf{x}}$ where
$h_{1}\neq h_{2}$. Then we have
$C_{\omega}(g_{1}\cdot\mathbf{x})=C_{\omega}(gh_{1}\cdot\mathbf{x})=C_{\omega}(g\cdot\mathbf{x})=C_{\omega}(gh_{2}\cdot\mathbf{x})=C_{\omega}(g_{2}\cdot\mathbf{x})$,
implying that
$\rho(g_{1})C_{\omega}(\mathbf{x})=\rho(g_{2})C_{\omega}(\mathbf{x})$. Since
$g_{1}\neq g_{2}$, this contradicts the group axiom, and thus an equivariant
canonicalizer cannot exist for inputs with non-trivial inner-symmetries. To
handle all inputs, canonicalization [41] adopts relaxed equivariance: a
canonicalizer $C_{\omega}$ satisfies relaxed equivariance if
$C_{\omega}(g\cdot\mathbf{x})=\rho(gh)C_{\omega}(\mathbf{x})$ up to arbitrary
action from the stabilizer $h\in G_{\mathbf{x}}$. As a result, the
symmetrization $\phi_{\theta,\omega}(\mathbf{x})=g\cdot
f_{\theta}(g^{-1}\cdot\mathbf{x})$ performed using a relaxed canonicalizer
$C_{\omega}$ only guarantees relaxed equivariance
$\phi_{\theta,\omega}(g\cdot\mathbf{x})=gh\cdot\phi_{\theta,\omega}(\mathbf{x})$
up to arbitrary action from the stabilizer $h\in G_{\mathbf{x}}$ (Theorem A.2
of [41]). In other words, canonicalization does not guarantee equivariance for
inputs with inner symmetries.

To visually demonstrate this, we perform a minimal experiment using several
graphs $\mathbf{x}$ with non-trivial stabilizers $G_{\mathbf{x}}$, _i.e._ ,
inner symmetries, taken from [92]. We fix a randomly initialized MLP
$f_{\theta}:\mathbb{R}^{n\times n}\to\mathbb{R}^{n}$ and symmetrize it using
our approach and canonicalization. When symmetrized, the MLP is expected to
provide scalar embedding of each node, which we color-code for visualization.
The results are in Figures 11, 11, and 11. For each graph, we illustrate three
panels: the leftmost one illustrates the color-coding of the inner symmetry of
nodes (automorphism), the middle one illustrates node embedding from MLP-PS,
and the rightmost one illustrates embedding from MLP-Canonical. If a method is
$G$ equivariant, it is expected to give identical embeddings for automorphic
nodes, since an equivariant model cannot distinguish them in principle [90].
As in the Figures 11, 11, and 11, in the presence of inner symmetry (left
panels), MLP with probabilistic symmetrization (middle panels) achieves $G$
equivariance and produces close embeddings for automorphic nodes. However, the
same MLP with canonicalization produces relatively unstructured embeddings
(right panels). The result illustrates a potential advantage of probabilistic
symmetrization over canonicalization when learning data with inner symmetries,
which is often found in applications such as molecular graphs [62].

###  A.5 Limitations and Broader Impacts (Continued from Section 4)

While the equivariance, universality, simplicity, and scalability of our
approach offers a potential for positive impact for deep learning for
chemistry, biology, physics, and mathematics, it also has limitations and
potential negative impacts. The main limitation of our work is that it trades
off certain desirable traits in equivariant deep learning in favor of
achieving architecture agnostic equivariance. For example, (1) our approach is
less interpretable compared to equivariant architectures due to less
structured computations in the base model, (2) our approach is presumably less
parameter and data efficient compared to equivariant architectures due to less
imposed prior knowledge on parameterization, and (3) our approach is expected
to be challenged when input size generalization is required, partially because
the maximum input size has to be specified in advance. Another genuine
weakness compared to canonicalization is that, our method is stochastic and
therefore incurs $\mathcal{O}(N)$ cost when using $N$ samples for estimation.
These limitations might lead to potential negative environmental impacts,
since less interpretability and lower efficiency implies higher reliance on
larger models with more computation cost. We acknowledge the aforementioned
limitations and impacts of our work, and will make effort to address them in
follow-up research. For example, we believe data efficiency of our approach
could improve with pretrained knowledge transfer from other domains, as it
would impose a strong prior on the hypothesis space and may work similarly to
architectural priors that benefit data efficiency. Also, for the sampling
cost, since the sampling is completely parallelizable and analogous to using a
larger batch size, we believe it can be overcome to some degree by leveraging
parallel computing techniques developed for scaling batch size.