Integrability of Nonholonomic Heisenberg Type Systems

Integrability of Nonholonomic
Heisenberg Type Systems

Yury A. GRIGORYEV †, Alexey P. SOZONOV † and Andrey V. TSIGANOV †‡

Yu.A. Grigoryev, A.P. Sozonov and A.V. Tsiganov

† St. Petersburg State University, St. Petersburg, Russia
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
‡ Udmurt State University, Izhevsk, Russia

Received March 17, 2016, in final form November 22, 2016; Published online
November 25, 2016

We show that some modern geometric methods of Hamiltonian dynamics can be
directly applied to the nonholonomic Heisenberg type systems. As an example we
present characteristic Killing tensors, compatible Poisson brackets, Lax
matrices and classical $r$-matrices for the conformally Hamiltonian vector
fields obtained in a process of reduction of Hamiltonian vector fields by a
nonholonomic constraint associated with the Heisenberg system.

Hamiltonian dynamics; nonholonomic systems

37J60; 70G45; 70H45

##  1 Introduction

Hamiltonian mechanics emerged in 1833 as a convenient reformulation of
classical Newtonian mechanics and analytical Lagrangian mechanics. Over the
next almost two centuries, it was gradually realized that Hamiltonian
formulation of a system of differential equations

$\displaystyle\frac{d}{dt}x_{k}=X_{k},\qquad X=PdH$

has several advantages based on the main postulate that evolution of a
physical system over time is governed by a single Hamiltonian $H$ of that
system and a Poisson bivector $P$ describing geometry and topology of the
phase space. Many different mathematical methods and concepts are widely used
in the Hamiltonian formalism. For instance, we can construct first integrals
of the Hamiltonian vector field $X$ using Killing tensors (Riemannian
geometry), compatible Poisson brackets (bi-Hamiltonian geometry) or Lax
matrices (algebraic geometry, Lie algebras theory, classical $r$-matrix theory
etc.).

It is natural to ask what happens to these mathematical methods when we impose
nonholonomic constraints on the Hamiltonian system and make a suitable
reduction. Usually presence of the constraints drastically modifies or
destroys these geometric constructions even when one gets a conformally
Hamiltonian vector field

$\displaystyle\hat{X}=\mu(x)\hat{P}d\hat{H}$

after the reduction, see [2, 5, 7, 8, 9, 10, 11, 15, 22, 24] and references
within. Here $\hat{H}$ and $\hat{P}$ are Hamiltonian and Poisson bivector on
reduced phase space, whereas function $\mu(x)$ is the so-called conformal
factor.

The main aim of this note is to present a family of conformally Hamiltonian
dynamical systems for which all the geometric methods listed above can be
applied without any additional modifications. The corresponding nonholonomic
constraint is associated with the so-called Heisenberg system or nonholonomic
integrator, which plays an important role in both nonlinear control and
nonholonomic dynamics [7]. We have to underline that the integrable systems
discussed in this note look quite artificial from the viewpoint of control
theory and mechanics but they provide an example of the standard geometric
methods applicability to the conformally Hamiltonian systems.

This paper is organized as follows. Section 1.1 recalls a brief description of
the constrained motion in three-dimensional Euclidean space. Section 2
contains the main results on reduced motion of the nonholonomic Heisenberg
type systems on the plane. We will show that integrable potentials for this
non Hamiltonian vector field satisfy to the Bertrand–Darboux type equation.
Solutions of this equation and the corresponding characteristic coordinates
will be explicitly determined. Section 3 deals with application of the
standard Stäckel theory to the conformally Hamiltonian vector field. We will
discuss a construction of the Stäckel matrices, compatible Poisson structures,
Lax matrices and classical $r$-matrices.

###  1.1 Main definitions

According to [7, 18, 24] we take the standard Hamiltonian equations of motion
in Euclidean space $\mathbb{R}^{3}$

$\displaystyle\dot{q_{i}}=\frac{\partial H}{\partial
p_{i}},\qquad\dot{p}_{i}=-\frac{\partial H}{\partial q_{i}},\qquad i=1,2,3,$
(1.1)

where

$\displaystyle
H=\frac{1}{2}\big{(}p_{1}^{2}+p_{2}^{2}+p_{3}^{2}\big{)}+V(q_{1},q_{2},q_{3}).$
(1.2)

On the phase space $M=T^{*}\mathbb{R}^{3}$ we can introduce coordinates
$x=(q,p)$ in which equations (1.1) have form $\dot{x}_{i}=X_{i}$ and to
determine Hamiltonian vector field

$\displaystyle X=\sum_{i=1}^{6}X_{i}\frac{\partial}{\partial x_{i}},$

which is a linear operator on a space of the smooth functions on $M$ that
encodes the evolution of any quantity

$\displaystyle\dot{F}=X(F)=\sum_{i=1}^{6}X_{i}\frac{\partial F}{\partial
x_{i}}.$

Let us impose the constraint of first order in momenta (velocities)

$\displaystyle f=(b,p)=0,$ (1.3)

where $b=(b_{1},b_{2},b_{3})$ is a vector depending on coordinates $q$ and
$(x,y)$ means an inner product in $\mathbb{R}^{3}$. In this case equations of
motion are written in the following form

$\displaystyle\dot{q_{i}}=\frac{\partial H}{\partial
p_{i}},\qquad\dot{p}_{i}=-\frac{\partial H}{\partial q_{i}}+\lambda
b_{i},\qquad i=1,2,3,$ (1.4)

together with the constraint equation (1.3). The corresponding vector field
looks like an additive perturbation of the initial Hamiltonian vector field

$\displaystyle\hat{X}=X\big{|}_{b=0}+\lambda\left(b_{1}\frac{\partial}{\partial
p_{1}}+b_{2}\frac{\partial}{\partial p_{2}}+b_{3}\frac{\partial}{\partial
p_{3}}\right),$

where unknown Lagrange multiplier $\lambda$ has to be computed from the
condition

$\displaystyle\dot{f}=\hat{X}(f)=X(f)\big{|}_{b=0}+\lambda(b,b)=0,$

such that

$\displaystyle\lambda=\frac{X(f)|_{b=0}}{(b,b)}.$ (1.5)

Here $X(f)|_{b=0}$ denotes the vector field in the absence of constraint. Such
equations and an equivalence of Hamiltonian and Lagrangian reductions are
carefully discussed in the book [7].

For completeness and self-sufficiency of presentation we consider integrable
constraint

$\displaystyle f=p_{3}=0,\qquad b=(0,0,1),$ (1.6)

when a third component of momenta is equal to zero, and non integrable
constraint

$\displaystyle f=p_{3}-(q_{2}p_{1}-q_{1}p_{2})=0,\qquad b=(-q_{2},q_{1},1),$
(1.7)

when third components of momenta and angular momenta coincide with each other.

If we assume that the potential $V$ in $\mathbb{R}^{3}$ does not depend on
$q_{3}$, then Lagrange multipliers (1.5) are equal to

$\displaystyle\lambda=0\qquad\mbox{and}\qquad\lambda=\frac{q_{1}\partial_{2}V-q_{2}\partial_{1}V}{1+q_{1}^{2}+q_{2}^{2}},\qquad\partial_{j}=\frac{\partial}{\partial
q_{j}},$

respectively. Thus, the Hamiltonian system (1.2) with integrable constraint
(1.6) represent a Hamiltonian system on the plane $\mathbb{R}^{2}$ embedded in
a three dimensional Euclidean space $\mathbb{R}^{3}$ in which there is no
force acting on the third component. If we impose non integrable constraint
(1.7), the equations for $q_{3}$ also decouples from the rest of the system
(1.4) and we obtain a two-degrees of freedom non-Hamiltonian system on the
plane. Following to [7, 18] we will call dynamical system associated with
$\hat{X}$ as a Heisenberg type system.

Integrability of the Hamiltonian systems on the plane is the well-known
problem of classical mechanics considered by Bertrand [4] and Darboux [12].
Later on results of this investigation due to Darboux were included in the
textbook of Whittaker on analytical mechanics [25] almost verbatim. In the
next Section we compare classical Bertrand–Darboux theorem with its non-
Hamiltonian counterpart, which appears in the nonholonomic case.

##  2 Reduced systems on the plane

When $V=V(q_{1},q_{2})$, substituting integrable constraint $p_{3}=0$ (1.6)
into (1.2) one gets

$\displaystyle
H_{1}=\frac{1}{2}\bigl{(}p_{1}^{2}+p_{2}^{2}\bigr{)}+V(q_{1},q_{2})=\frac{1}{2}\mathrm{g}^{ij}p_{i}p_{j}+V(q_{1},q_{2}),\qquad\mathrm{g}=\left(\begin{matrix}1&0\\\
0&1\end{matrix}\right)$ (2.1)

with the summation convention in force; we regard the $\mathrm{g}^{ij}$ as
components of the covariant form of the metric tensor $\mathrm{g}$ on the
plane, and use this metric freely to raise and lower indices.

After standard Hamiltonian reduction by cyclic third coordinate, the original
Hamiltonian vector field (1.1) in $T^{*}\mathbb{R}^{3}$ becomes the
Hamiltonian vector field in $T^{*}\mathbb{R}^{2}$

$\displaystyle\dot{x}=X,\qquad X=PdH_{1},\qquad P=\left(\begin{matrix}0&I\\\
-I&0\end{matrix}\right),$ (2.2)

where $x=(q_{1},q_{2},p_{1},p_{2})$. The Poisson bivector $P$ in (2.2) defines
canonical Poisson bracket

$\displaystyle\\{q_{i},p_{i}\\}=1,\qquad\\{q_{i},p_{k}\\}=\\{p_{1},p_{2}\\}=0,\qquad
i\neq k.$ (2.3)

Of course, this vector field $X$ preserves Poisson structure $P$, energy
$H_{1}$ and the standard volume 2-form $\Omega=dq\wedge dp$.

Substituting non integrable constraint $p_{3}=q_{2}p_{1}-q_{1}p_{2}$ (1.7)
into (1.2) one gets

$\displaystyle\hat{H}_{1}=\frac{1}{2}\bigl{(}p_{1}^{2}+p_{2}^{2}+(q_{2}p_{1}-q_{1}p_{2})^{2}\bigr{)}+V(q_{1},q_{2})=\frac{1}{2}\hat{\mathrm{g}}^{ij}p_{i}p_{j}+V(q_{1},q_{2}),$
(2.4)
$\displaystyle\hat{\mathrm{g}}=\left(\begin{matrix}1+q_{2}^{2}&-q_{1}q_{2}\\\
-q_{1}q_{2}&1+q_{1}^{2}\end{matrix}\right).$

The second order contravariant symmetric tensor field $\hat{\mathrm{g}}$ is a
Killing tensor with respect to the standard metric tensor $\mathrm{g}$, i.e.,
it satisfies to the Killing equation

$\displaystyle[\\![\mathrm{g},\hat{\mathrm{g}}]\\!]=0,$

where $[\\![\cdot,\cdot]\\!]$ is a Schouten bracket. It allows us to construct
conformal Killing tensor of second order

$\displaystyle
G=\hat{\mathrm{g}}-\operatorname{tr}\left(\hat{\mathrm{g}}\right)\mathrm{g}=-\left(\begin{matrix}1+q_{1}^{2}&q_{1}q_{2}\\\
q_{1}q_{2}&1+q_{2}^{2}\end{matrix}\right)$

with vanishing Nijenhuis torsion and to define the standard Turiel deformation
of the standard symplectic form

$\displaystyle\hat{\omega}=d\big{(}G^{ij}p_{i}dq_{j}\big{)},$

see [21, 22] and references within. The corresponding Poisson bracket

$\displaystyle\\{q_{i},p_{i}\\}^{\wedge}=1+q_{i}^{2},\qquad\\{q_{i},p_{j}\\}^{\wedge}=q_{1}q_{2},\qquad\\{p_{1},p_{2}\\}^{\wedge}=q_{1}p_{2}-q_{2}p_{1},\qquad
i\neq j$ (2.5)

is compatible with canonical bracket (2.3). Because the Turiel deformation is
trivial deformation in the Lichnerowicz–Poisson sense, there is a change of
variables

$\displaystyle
p_{1}\to\pi_{1}=\frac{q_{1}q_{2}p_{2}-\big{(}1+q_{2}^{2}\big{)}p_{1}}{1+x_{1}^{2}+x_{2}^{2}},\qquad
p_{2}\to\pi_{2}=\frac{q_{1}q_{2}p_{1}-\big{(}1+q_{1}^{2}\big{)}p_{2}}{1+x_{1}^{2}+x_{2}^{2}},$
(2.6)

which transforms this bracket to the canonical one (2.3), similar to other
nonholonomic systems with constraints of first order in momenta [6, 10, 20].

After the nonholonomic reduction, original Hamiltonian vector field (1.1) in
$T^{*}\mathbb{R}^{3}$ becomes a conformally Hamiltonian vector field in
$T^{*}\mathbb{R}^{2}$

$\displaystyle\dot{x}=\hat{X},\qquad\hat{X}=\mu\hat{P}d\hat{H}_{1},\qquad\hat{P}=\hat{\omega}^{-1},$
(2.7)

where conformal factor

$\displaystyle\mu=\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}^{-1}$

is a nowhere vanishing smooth function on the plane. This vector field
possesses energy $\hat{H}_{1}$ and the volume 2-form $\hat{\Omega}=\mu
dq\wedge dp$, but it does not preserve the Poisson bivector $\hat{P}$.

###  2.1 Bertrand–Darboux type equation

An existence of the integrals of second order in velocities for Hamiltonian
vector field (2.2) is described by a classical Bertrand–Darboux theorem [4,
12, 25].

######  Proposition 2.1.

The function (2.1)

$\displaystyle
H_{1}=\frac{1}{2}\mathrm{g}^{ij}p_{i}p_{j}+V(q_{1},q_{2}),\qquad\mathrm{g}=\left(\begin{matrix}1&0\\\
0&1\end{matrix}\right)$ (2.8)

defines integrable Hamiltonian vector field $X$ (2.2), which has an
independent quadratic first integral

$\displaystyle H_{2}=K^{ij}p_{i}p_{j}+U(q_{1},q_{2}),$ (2.9)

if the second order contravariant symmetric tensor field $K$ obeys the Killing
tensor equation

$\displaystyle[\\![\mathrm{g},K]\\!]=0,$ (2.10)

and potential $V$ satisfies the compatibility condition

$\displaystyle d(\mathbf{K}dV)=0,\qquad\mathbf{K}=K\mathrm{g}^{-1},$ (2.11)

where $\mathbf{K}$ is the tensor field of $(1,1)$ type.

A characteristic coordinate system for (2.11) provides separation for the
potential $V$ and can be taken as one of the following four orthogonal
coordinate systems on the plane: elliptic, parabolic, polar or Cartesian.

If a non-trivial solution exists, the Bertrand–Darboux equation (2.11) can be
reduced to canonical form by transforming to characteristic coordinates, which
appear to be separation coordinates for the Hamilton–Jacobi equation related
to the natural Hamiltonian (see Bertrand [4], Darboux [12], or Whittaker [25]
and Ankiewicz and Pask [1] for a full proof of the BD theorem). The modern
proof of this classic statement may be found in [19].

In nonholonomic case we can also substitute the same ansatz (2.9) into the
equation $\dot{H}_{2}=0$ and get the following generalisation of the
Bertrand–Darboux result.

######  Proposition 2.2.

The function (2.4)

$\displaystyle\hat{H}_{1}=\frac{1}{2}\hat{\mathrm{g}}^{ij}p_{i}p_{j}+V(q_{1},q_{2}),\qquad\hat{\mathrm{g}}=\left(\begin{matrix}1+q_{2}^{2}&-q_{1}q_{2}\\\
-q_{1}q_{2}&1+q_{1}^{2}\end{matrix}\right)$ (2.12)

defines integrable conformally Hamiltonian vector field $\hat{X}$ (2.7), which
has an independent quadratic first integral

$\displaystyle\hat{H}_{2}=K^{ij}p_{i}p_{j}+U(q_{1},q_{2}),$ (2.13)

if the second order contravariant symmetric tensor field $K$ is a Killing
tensor

$\displaystyle[\\![\mathrm{g},K]\\!]=0,\qquad\mathrm{g}=\left(\begin{matrix}1&0\\\
0&1\end{matrix}\right)$ (2.14)

and potential $V$ satisfies the compatibility condition

$\displaystyle
d\big{(}\hat{\mathbf{K}}dV\big{)}=0,\qquad\hat{\mathbf{K}}=K\hat{\mathrm{g}}^{-1}.$
(2.15)

A characteristic coordinate system for (2.15) provides separation for
potential $V(q_{1},q_{2})$.

In this case existence of the second integral of motion guarantees
integrability of the given non Hamiltonian vector field $\hat{X}$ according to
the Euler–Jacobi theorem [17]. Equation (2.15) was briefly discussed in [23]
using coordinates $q_{1,2}$ and momenta $\pi_{1,2}$ (2.6).

######  Remark 2.3.

In contrast with the standard Bertrand–Darboux theorem the Killing equation
(2.14) and compatibility condition (2.15) include two different metrics
$\mathrm{g}$ and $\hat{\mathrm{g}}$. It happens because the reaction force is
entered in the second part of the vector field (1.4) only. Sequentially, $K$
is a Killing tensor with respect to the Euclidean metric $\mathrm{g}$ to the
plane, which is an induced metric associated with standard embedding
$\mathbb{R}^{2}\subset\mathbb{R}^{3}$ in Cartesian coordinates. However, $K$
is not a Killing tensor with respect to metric $\hat{\mathrm{g}}$, which
appears in the reduced Hamiltonian (2.4) after nonholonomic reduction and then
defines raising and lowering indices in the compatibility condition (2.15).

###  2.2 Solutions of the Bertrand–Darboux type equation

In [12] Darboux proceeds to find the unknown potential $V(q_{1},q_{2})$ by
solving the compatibility condition (2.11) using method of characteristics.
Solving the second compatibility condition (2.15) we can apply the same
method.

In both cases characteristic coordinates for (2.11) and (2.15) consist of
eigenvalues of the corresponding $(1,1)$ tensor field ${\mathbf{K}}$ or
$\hat{\mathbf{K}}$. In order to describe these coordinates we can start with a
well-known generic solution of the common Killing tensor equations (2.10) and
(2.14)

$\displaystyle
K=\left(\begin{matrix}c_{1}q_{2}^{2}+2c_{2}q_{2}+c_{3}&-c_{1}q_{1}q_{2}-c_{2}q_{1}-c_{4}q_{2}+c_{6}\\\
-c_{1}q_{1}q_{2}-c_{2}q_{1}-c_{4}q_{2}+c_{6}&c_{1}q_{1}^{2}+2c_{4}q_{1}+c_{5}\end{matrix}\right),$
(2.16)

which depends on six constants of integration $c_{1},\ldots,c_{6}$.

Substituting this generic solution into the standard compatibility condition
(2.11) one gets

$\displaystyle\bigl{(}A\partial_{11}+B\partial_{22}+C\partial_{12}+a\partial_{1}+b\partial_{2}\bigr{)}V(q_{1},q_{2})=0,$
(2.17)

where $\partial_{i}={\partial}/{\partial q_{i}}$,
$\partial_{ik}={\partial^{2}}/{\partial q_{i}\partial q_{k}}$ and polynomials
of second order read as

$\displaystyle A=c_{1}q_{1}q_{2}+c_{2}q_{1}+c_{4}q_{2}-c_{6},\\!\\!\qquad
B=-A,\\!\\!\qquad
C=c_{1}q_{2}^{2}-c_{1}q_{1}^{2}+2c_{2}q_{2}-2c_{4}q_{1}+c_{3}-c_{5},$

whereas polynomials of the first order are

$\displaystyle a=3(c_{1}q_{2}+c_{2}),\qquad b=-3(c_{1}q_{1}+c_{4}).$

The linear PDE (2.17) of second order was obtained by Bertrand [4] and studied
by Darboux [12]. Thus, it was later called the Bertrand–Darboux equation [19,
25].

Solving (2.17) for $V(q_{1},q_{2})$ amounts to finding admissible potentials
of the Hamiltonian systems defined by $H$ (2.1), which integrability is
afforded by the existence of first integrals (2.9) which are quadratic in the
momenta. Second potential $U(q_{1},q_{2})$ in (2.9) is a solution of the
equation

$\displaystyle d(\mathbf{K}dU)=V.$

In nonholonomic case, substituting the same generic solution (2.16) into the
compatibility condition (2.15) one gets similar linear PDE

$\displaystyle\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}\big{(}\hat{A}\partial_{11}+\hat{B}\partial_{22}+\hat{C}\partial_{12}\big{)}V(q_{1},q_{2})+\big{(}\hat{a}\partial_{1}+\hat{b}\partial_{2}\big{)}V(q_{1},q_{2})=0,$
(2.18)

with coefficients

$\displaystyle\hat{A}=\big{(}1+q_{1}^{2}\big{)}(q_{1}c_{2}-c_{6})+q_{2}\big{(}1-q_{1}^{2}\big{)}c_{4}+q_{1}q_{2}(c_{1}-c_{5}),$
$\displaystyle\hat{B}=-\big{(}1+q_{2}^{2}\big{)}(q_{2}c_{4}-c_{6})-q_{1}\big{(}1-q_{2}^{2}\big{)}c_{2}-q_{1}q_{2}(c_{1}-c_{3}),$
$\displaystyle\hat{C}=-\big{(}q_{1}^{2}-q_{2}^{2}\big{)}c_{1}+\big{(}1+q_{1}^{2}\big{)}(2q_{2}c_{2}+c_{3})-\big{(}1+q_{2}^{2}\big{)}(2q_{1}c_{4}+c_{5}),$

and

$\displaystyle\hat{a}=\big{(}q_{1}^{2}+q_{2}^{2}+3\big{)}\big{(}c_{1}q_{2}+2c_{2}q_{1}^{2}-2c_{4}q_{1}q_{2}\big{)}-q_{2}\big{(}2\big{(}q_{1}^{2}+1\big{)}c_{3}-\big{(}q_{1}^{2}-q_{2}^{2}-1\big{)}c_{5}\big{)}$
$\displaystyle\hphantom{\hat{a}=}{}+q_{1}\big{(}q_{1}^{2}-3q_{2}^{2}+1\big{)}c_{6}-\big{(}q_{1}^{2}+q_{2}^{2}-3\big{)}c_{2},$
$\displaystyle\hat{b}=-\big{(}q_{1}^{2}+q_{2}^{2}+3\big{)}\big{(}c_{1}q_{1}+2c_{4}q_{2}^{2}-2c_{2}q_{1}q_{2}\big{)}+q_{1}\big{(}2\big{(}q_{2}^{2}+1\big{)}c_{5}+\big{(}q_{1}^{2}-q_{2}^{2}+1\big{)}c_{3}\big{)}$
$\displaystyle\hphantom{\hat{b}=}{}+q_{2}\big{(}3q_{1}^{2}-q_{2}^{2}-1\big{)}c_{6}+\big{(}q_{1}^{2}+q_{2}^{2}-3\big{)}c_{4}.$

Solutions $V(q_{1}q_{2})$ of this Bertrand–Darboux type equation (2.18)
determine all the admissible Jacobi integrals (2.4), which define integrable
conformally Hamiltonian vector fields (2.7) with first integrals (2.13) of
second order in momenta. In this case potential $U(q_{1},q_{2})$ in (2.13) is
a solution of the equation

$\displaystyle d\big{(}\hat{\mathbf{K}}dU\big{)}=V.$

######  Remark 2.4.

Before solving the PDE (2.17), Darboux ingeniously observes that it can be
simplified without loss of generality [12]. Indeed, by rotating and
translating the axes, one can simplify the general solution of the Killing
tensor equation, thus bringing it to a certain canonical form. In modern
language solving the equivalence and canonical forms problem for the Killing
tensor equations (2.10), (2.14) is equivalent to analysing the orbits of the
six-dimensional vector space of solutions under the action of the Lie group of
orientation-preserving isometries.

It is well known that there are four types of orbits generated by the
following canonical Killing tensors:

Cartesian: $\displaystyle K^{(1)}=\left(\begin{matrix}1&0\\\
0&0\end{matrix}\right),$ polar: $\displaystyle
K^{(2)}=\left(\begin{matrix}q_{2}^{2}&-q_{1}q_{2}\\\
-q_{1}q_{2}&q_{1}^{2}\end{matrix}\right),$ parabolic: $\displaystyle
K^{(3)}=\left(\begin{matrix}0&-q_{2}\\\ -q_{2}&2q_{1}\\\ \end{matrix}\right),$
elliptic: $\displaystyle
K^{(4)}=\left(\begin{matrix}\kappa^{2}+q_{2}^{2}&-q_{1}q_{2}\\\
-q_{1}q_{2}&q_{1}^{2}\end{matrix}\right),\qquad\kappa\in\mathbb{R},$ (2.19)

associated with Cartesian, polar, parabolic and elliptic coordinate systems on
the plane, respectively.

Let us consider these characteristic coordinate systems associated with the
Killing tensors (2.19).

Elliptic coordinates. Substituting Killing tensor $K^{(4)}$ into the (1,1)
tensor $\mathbf{K}=K^{(4)}\mathrm{g}^{-1}$ and calculating its eigenvalues one
gets elliptic coordinates

$\displaystyle\xi_{1,2}=\kappa+q_{1}^{2}+q_{2}^{2}\pm\sqrt{\kappa^{2}-2\big{(}q_{1}^{2}-q_{2}^{2}\big{)}\kappa+\big{(}q_{1}^{2}+q_{2}^{2}\big{)}^{2}}$

up to the constant factor. In nonholonomic case, substituting the same tensor
$K^{(4)}$ into the characteristic tensor
$\hat{\mathbf{K}}=K^{(4)}\hat{\mathrm{g}}^{-1}$, one gets the following
eigenvalues

$\displaystyle\hat{\xi}_{1,2}=\frac{\big{(}1+q_{1}^{2}\big{)}\kappa+q_{1}^{2}+q_{2}^{2}\pm\sqrt{\big{(}1+q_{1}^{2}\big{)}^{2}\kappa^{2}-2\big{(}q_{1}^{2}\big{(}q_{1}^{2}+q_{2}^{2}\big{)}+q_{1}^{2}-q_{2}^{2}\big{)}\kappa+\big{(}q_{1}^{2}+q_{2}^{2}\big{)}^{2}}}{1+q_{1}^{2}+q_{2}^{2}}.$

Solutions of the compatibility equations (2.11) and (2.15) are labelled by two
arbitrary functions $F_{1}$ and $F_{2}$ on these coordinates:

$\displaystyle
V(q_{1},q_{2})=\frac{F_{1}(u_{1})-F_{2}(u_{2})}{u_{1}-u_{2}},\qquad
u_{1,2}=\xi_{1,2}\qquad\mbox{or}\qquad u_{1,2}=\hat{\xi}_{1,2}.$

Parabolic coordinates. In Hamiltonian case eigenvalues of the characteristic
tensor $\mathbf{K}=K^{(3)}\mathrm{g}^{-1}$ are standard parabolic coordinates

$\displaystyle\zeta_{1,2}=q_{1}\pm\sqrt{q_{1}^{2}+q_{2}^{2}},$ (2.20)

whereas in nonholonomic case eigenvalues of the characteristic tensor
$\hat{\mathbf{K}}=K^{(3)}\hat{\mathrm{g}}^{-1}$ read as

$\displaystyle\hat{\zeta}_{1,2}=\frac{q_{1}\pm\sqrt{q_{2}^{2}+1}\sqrt{q_{1}^{2}+q_{2}^{2}}}{1+q_{1}^{2}+q_{2}^{2}}.$
(2.21)

It is easy to express new characteristic coordinates via standard ones

$\displaystyle\hat{\zeta}_{1,2}=\frac{2\big{(}\zeta_{1}+\zeta_{2}\pm(\zeta_{1}-\zeta_{2})\sqrt{1-\zeta_{1}\zeta_{2}}\big{)}}{(\zeta_{1}-\zeta_{2})^{2}+4}.$

In both cases the desired separable potentials are equal to

$\displaystyle
V(q_{1},q_{2})=\frac{F_{1}(u_{1})-F_{2}(u_{2})}{u_{1}-u_{2}},\qquad
u_{1,2}=\zeta_{1,2}\qquad\mbox{or}\qquad u_{1,2}=\hat{\zeta}_{1,2}.$

Polar coordinates. Using one nontrivial eigenvalue of characteristic tensors
$\mathbf{K}=K^{(2)}\mathrm{g}^{-1}$ and
$\hat{\mathbf{K}}=K^{(2)}\hat{\mathrm{g}}^{-1}$ we can introduce only one
coordinate

$\displaystyle
r=\sqrt{q_{1}^{2}+q_{2}^{2}},\qquad\mbox{or}\qquad\hat{r}=\sqrt{\frac{q_{1}^{2}+q_{2}^{2}}{1+q_{1}^{2}+q_{2}^{2}}}.$

The second coordinate is a function on $q_{1}/q_{2}$, for instance, it could
be an angle $\varphi=\arctan q_{1}/q_{2}$, because solutions of the
compatibility equations (2.11) and (2.15) have the following form

$\displaystyle
V(q_{1},q_{2})=F_{1}(\rho)+\frac{F_{2}(\varphi)}{\rho^{2}},\qquad\rho=r\qquad\mbox{or}\qquad\rho=\hat{r}.$

Cartesian coordinates. Separable in the Cartesian coordinates solution of the
compatibility condition (2.11) reads as

$\displaystyle V(q_{1},q_{2})=F_{1}(q_{1})+F_{2}(q_{2}).$ (2.22)

In nonholonomic case characteristic tensor
$\hat{\mathbf{K}}=K^{(1)}\hat{\mathrm{g}}^{-1}$ yields only one coordinate

$\displaystyle\hat{\varrho}=\sqrt{\frac{1+q_{1}^{2}}{1+q_{1}^{2}+q_{2}^{2}}}$

and the separable solution of the compatibility condition (2.15) looks like

$\displaystyle
V(q_{1},q_{2})=F_{1}(\hat{\varrho})+\frac{F_{2}(\hat{\phi})}{\hat{\varrho}^{2}},\qquad\hat{\phi}=\arctan
q_{1}.$

Thus, after nonholonomic reduction one gets separable potential, which is
different, even in the form of standard potential (2.22) separable in
Cartesian coordinates.

Summing up, we have found four characteristic coordinate systems for the
Bertrand–Darboux type equation (2.18) associated with the integrable non-
Hamiltonian vector field $\hat{X}$ (2.7). Additional first integral of this
field $\hat{X}$ is a polynomial of the second order in momenta defined by the
standard Killing tensor on the plane. Consequently, we can directly apply the
standard Stäckel theory to the nonholonomic Heisenberg type systems.

##  3 Stäckel systems

In Hamiltonian case there is a one-to-one correspondence between the so-called
Stäckel systems, integrable Killing tensors which mutually commute in the
algebraic sense and separation of variables [3, 14].

The nondegenerate $n\times n$ Stäckel matrix $S$, which $j$ column depends
only on variable $u_{j}$

$\displaystyle\det S\neq 0,\qquad\frac{\partial S_{kj}}{\partial
u_{m}}=0,\qquad j\neq m$

defines $n$ functionally independent integrals of motion

$\displaystyle
H_{k}=\sum_{j=1}^{n}C_{jk}\big{(}p_{u_{j}}^{2}+U_{j}(u_{j})\big{)},\qquad
C=S^{-1},$ (3.1)

which are in involution with respect to canonical Poisson brackets

$\displaystyle\\{u_{i},p_{u_{j}}\\}=\delta_{ij},\qquad\\{u_{i},u_{j}\\}=\\{p_{u_{i}},p_{u_{j}}\\}=0.$

The common level surface of the first integrals
$H_{1}=\alpha_{1},\ldots,H_{2}=\alpha_{n}$ is diffeomorphic to the
$n$-dimensional real torus and one immediately gets

$\displaystyle
p_{u_{j}}^{2}=\sum_{k=1}^{n}\alpha_{k}S_{kj}(u_{j})-U_{j}(u_{j}).$

It allows us to calculate quadratures for the corresponding Hamiltonian vector
field $X=PdH_{1}$:

$\displaystyle\sum_{j=1}^{n}\int_{\gamma_{0}(p_{0},q_{0})}^{\gamma_{j}(p_{u_{j}},u_{j})}\frac{S_{kj}(u)du}{\sqrt{\sum\limits_{k=1}^{n}\alpha_{k}S_{kj}(u)-U_{j}(u)}}=\beta_{k},\qquad
k=1,\ldots,n,$ (3.2)

where $\beta_{1}=t$ and $\beta_{2},\ldots,\beta_{n}$ are constants of
integration. Solution of the problem is thus reduced to solving a sequence of
one-dimensional problems, which is the essence of the method of separation of
variables [3].

###  3.1 Stäckel matrices for reduced systems

In the definition of the Stäckel integrals of motion (3.1) momenta $p_{u}$
have to be canonically conjugated to eigenvalues $u$ of the characteristic
$(1,1)$ tensor $\mathbf{K}$. For instance, if we take tensor $K^{(3)}$ (2.19)
the eigenvalues of $\mathbf{K}$ are parabolic coordinates $\zeta_{1,2}$ (2.20)
and the corresponding momenta read as

$\displaystyle
p_{{\zeta}_{1,2}}=\frac{p_{1}}{2}\pm\frac{\sqrt{q_{1}^{2}+q_{2}^{2}}\mp
q_{1}}{2q_{2}}p_{2}.$

In nonholonomic case the eigenvalues of $\hat{\mathbf{K}}$ are given by
$\hat{\zeta}_{1,2}$ (2.21) and the conjugated momenta have the following form

$\displaystyle p_{\hat{\zeta}_{1,2}}=\frac{q_{2}^{2}+1\pm
q_{1}\sqrt{1+q_{2}^{2}}\sqrt{q_{1}^{2}+q_{2}^{2}}}{2\big{(}1-q_{1}^{2}\big{)}}p_{1}-\frac{q_{1}\mp\sqrt{1+q_{2}^{2}}\sqrt{q_{1}^{2}+q_{2}^{2}}}{2q_{2}}p_{2}.$

In similar manner we can calculate all the canonical variables associated with
tensors $K^{(j)}$ (2.19) and to obtain standard Stäckel matrices in
Hamiltonian case

$\displaystyle S^{(1)}=\left(\begin{matrix}1&1\\\ 1&-1\\\
\end{matrix}\right),\qquad S^{(2)}=\left(\begin{matrix}0&1\\\ 1&-r^{-2}\\\
\end{matrix}\right),\qquad S^{(3)}=\left(\begin{matrix}1&1\\\
-\zeta_{1}^{-1}&-\zeta_{2}^{-1}\end{matrix}\right)$

and

$\displaystyle
S^{(4)}=\left(\begin{matrix}\dfrac{1}{\xi_{1}-\kappa^{2}}&\dfrac{1}{\xi_{2}-\kappa^{2}}\vspace{1mm}\\\
\dfrac{-1}{\xi_{1}(\xi_{1}-\kappa^{2})}&\dfrac{-1}{\xi_{2}(\xi_{2}-\kappa^{2})}\end{matrix}\right).$

For the nonholonomic Heisenberg type systems associated with tensors $K^{(j)}$
(2.19) Stäckel matrices are equal to

$\displaystyle\hat{S}^{(1)}=\left(\begin{matrix}0&\dfrac{1}{\hat{\varrho}^{2}(1-\hat{\varrho}^{2})}\vspace{1mm}\\\
1&\dfrac{-1}{\hat{\varrho}^{4}(1-\hat{\varrho}^{2})}\end{matrix}\right),\qquad\hat{S}^{(2)}=\left(\begin{matrix}0&\dfrac{1}{1-\hat{r}^{2}}\vspace{1mm}\\\
1&\dfrac{-1}{\hat{r}^{2}(1-\hat{r}^{2})}\end{matrix}\right)$

and

$\displaystyle\hat{S}^{(3)}=\left(\begin{matrix}\dfrac{1}{1-\hat{\zeta}_{1}}&\dfrac{1}{1-\hat{\zeta}_{2}}\vspace{1mm}\\\
\dfrac{-1}{\hat{\zeta}_{1}(1-\hat{\zeta}_{1})}&\dfrac{-1}{\hat{\zeta}_{2}(1-\hat{\zeta}_{2})}\end{matrix}\right),$
$\displaystyle\hat{S}^{(4)}=\left(\begin{matrix}\dfrac{1}{(\hat{\xi}_{1}-\kappa^{2})(1-\hat{\xi}_{1})}&\dfrac{1}{(\hat{\xi}_{2}-\kappa^{2})(1-\hat{\xi}_{2})}\vspace{1mm}\\\
\dfrac{-1}{\hat{\xi}_{1}(\hat{\xi}_{1}-\kappa^{2})(1-\hat{\xi}_{1})}&\dfrac{-1}{\hat{\xi}_{2}(\hat{\xi}_{2}-\kappa^{2})(1-\hat{\xi}_{2})}\end{matrix}\right).$

It is easy to see, that imposing linear, non integrable constraint (1.7) in
polar, parabolic and elliptic cases we have to multiply $j$ column of the
standard Stäckel matrices on the function $(1-u_{j})$, where $u_{j}$ is the
eigenvalue of $\hat{\mathbf{K}}$. Similar transformations of the Stäckel
matrices for nonholomic systems on a two-dimensional sphere are discussed in
[22].

Using a new time variable $\tau$ defined by

$\displaystyle dt=\mu^{-1}d\tau\equiv\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}d\tau$

we can rewrite the conformally Hamiltonian vector field $\hat{X}$ (2.7) in the
Hamiltonian form and obtain quadratures similar to (3.2).

Thus, in nonholonomic case solution of the problem is reduced to solving a
sequence of one-dimensional problems after suitable change of time.

###  3.2 Compatible Poisson brackets

Riemannian geometry is not, a priori, concerned with symplectic or Poisson
structures. Nevertheless, it is known that any tensor field $L$ with vanishing
Nijenhuis torsion on the Riemannian manifold $Q$ yields trivial deformation of
the canonical Poisson bracket $\\{\cdot,\cdot\\}$ on its cotangent bundle
$T^{*}Q$

$\displaystyle\\{q_{i},q_{j}\\}_{L}=0,\qquad\\{q_{i},p_{j}\\}_{L}=-L^{ij},\qquad\\{p_{i},p_{j}\\}_{L}=\left(\frac{\partial
L^{kj}}{\partial q_{i}}-\frac{\partial L^{ki}}{\partial q_{j}}\right)p_{k},$
(3.3)

where $(p_{i},q_{i})$ are fibered coordinates. This Poisson bracket is
compatible with the canonical one, so that there is a recursion operator
$N=P_{L}P^{-1}$, which allows us to construct a whole family of compatible
brackets associated with Poisson bivectors $P_{L}^{(m)}=N^{m}P$ on $T^{*}Q$
[21, 22].

According to [3] generic solution of the Killing equation $K$ (2.16)
determines a conformal Killing tensor with vanishing Nijenhuis torsion

$\displaystyle L=K-\operatorname{tr}(K)\mathrm{g}$

and Turiel’s deformation (3.3) of the canonical Poisson bracket [21, 22]. The
eigenvalues of the corresponding recursion operator $N=P_{L}P^{-1}$ are
characteristic coordinates. Integrals of motion $H_{1}$ (2.8) and $H_{2}$
(2.9) are in involution with respect to the compatible Poisson brackets

$\displaystyle\\{H_{1},H_{2}\\}=\\{H_{1},H_{2}\\}_{L}=0$

and to other polynomial Poisson brackets associated with the Poisson bivectors
$P_{L}^{(m)}$. Indeed, we can find such $L$-tensors directly from the Hamilton
function $H_{1}$ using modern software [16].

For the Heisenberg type systems, we have to solve the same Killing equation
(2.14) and second equation on potential with $(1,1)$ tensor
$\hat{\mathbf{K}}=K\hat{\mathrm{g}}$. This tensor does not appear in the
Turiel construction of the Poisson bracket and we can directly prove the
following proposition.

######  Proposition 3.1.

Integrals of motion $\hat{H}_{1}$ (2.12) and $\hat{H}_{2}$ (2.13) are in
involution

$\displaystyle\big{\\{}\hat{H}_{1},\hat{H}_{2}\big{\\}}^{\wedge}=\big{\\{}\hat{H}_{1},\hat{H}_{2}\big{\\}}_{L}=0$

with respect to the compatible Poisson brackets $\\{\cdot,\cdot\\}^{\wedge}$
(2.5) and $\\{\cdot,\cdot\\}_{L}$ (3.3). The eigenvalues of the recursion
operator

$\displaystyle\hat{N}=P_{L}\hat{P}^{-1},$

where $\hat{P}$ is given by (2.7), are characteristic coordinates discussed in
the previous Section.

The proof consists of the straightforward calculations.

Thus, starting with a common Poisson bivector $P_{L}$ we can construct two
recursion operators $N$ and $\hat{N}$ and two families of the Poisson brackets
on $T^{*}\mathbb{R}^{2}$.

###  3.3 Lax matrices

We can construct Lax representations for the Stäckel systems with uniform
rational and polynomial potentials $U_{j}=U$ in (3.1), see [13] and references
within. For the sake of brevity we consider only systems associated with
tensor $K^{(3)}$ (2.19).

If $V(q_{1},q_{2})=0$, we can introduce a well-known $2\times 2$ Lax matrix
for the geodesic motion on the plane

$\displaystyle\mathcal{L}(u)=\left(\begin{matrix}h&e\\\
f&-h\end{matrix}\right)(u)\equiv\left(\begin{matrix}\dfrac{1}{2}\dfrac{de(u)}{dt}&e(u)\vspace{1mm}\\\
-\dfrac{1}{2}\dfrac{d^{2}e(u)}{dt^{2}}&-\dfrac{1}{2}\dfrac{de(u)}{dt}\end{matrix}\right),$
(3.4)

which satisfies to the Lax equation

$\displaystyle\frac{d}{dt}\mathcal{L}=[\mathcal{L},\mathcal{A}]\equiv\mathcal{L}\mathcal{A}-\mathcal{A}\mathcal{L},\qquad\text{with}\qquad\mathcal{A}=\left(\begin{matrix}0&1\\\
0&0\end{matrix}\right).$ (3.5)

This equation guarantees that the eigenvalues of $\mathcal{L}$ are conserved
quantities in the involution.

In our case function $e(u)$ depends on the so-called spectral parameter $u$

$\displaystyle
e(u)=\frac{(u-\zeta_{1})(u-\zeta_{2})}{u}=u-2q_{1}-\frac{q_{2}^{2}}{u},$

which is completely defined by the eigenvalues of the characteristic tensor
$\mathbf{K}=K^{(3)}\mathrm{g}^{-1}$ or parabolic coordinates $\zeta_{1,2}$ on
the plane.

The involution property of eigenvalues of the Lax matrix $\mathcal{L}$ is
equivalent to existence of a classical $r$-matrix $r_{12}$, that as

$\displaystyle\\{\mathcal{L}_{1}(u),\mathcal{L}_{2}(v)\\}=[r_{12}(u,v),\mathcal{L}_{1}(u)]-[r_{21}(u,v),\mathcal{L}_{2}(v)].$

Here we use the familiar notation for tensor product of $\mathcal{L}$ and unit
matrix $I$

$\displaystyle\mathcal{L}_{1}(u)=\mathcal{L}(u)\otimes\mathrm{I},\qquad\mathcal{L}_{2}(v)=\mathrm{I}\otimes\mathcal{L}(v),\qquad
r_{21}(u,v)=\Pi r_{12}(v,u)\Pi,$

whereas $\Pi$ is the permutation operator: ${\Pi}x\otimes y=y\otimes x$,
$\forall\,x,y$. Evaluating canonical Poisson brackets (2.3) between entries of
the Lax matrix $\mathcal{L}(u)$ (3.4) we obtain a constant $r$-matrix

$\displaystyle
r_{12}(u,v)=\frac{2}{v-u}\Pi\equiv\frac{2}{v-u}\left(\begin{matrix}1&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\end{matrix}\right).$ (3.6)

If we want to consider separable in parabolic coordinates potentials
$V(q_{1},q_{2})$, we have to add some items to the initial Lax matrix
$\mathcal{L}(u)$ (3.4). For instance, we can take the following additive
perturbation

$\displaystyle\mathcal{L}_{V}(u)=\mathcal{L}(u)+\left(\begin{matrix}0&0\\\
\Delta f&0\end{matrix}\right),\qquad\Delta f=\big{[}\phi(u)e^{-1}(u)\big{]},$

where $\phi(u)$ is an arbitrary constant function and $[g(u)]$ is a truncated
Laurent series expansion of $g(u)$ with respect to the variable $u$. For
instance, when $\phi(u)=-au^{3}$ the Taylor part of expansion about $u=\infty$
looks like

$\displaystyle\big{[}{-}au^{3}e^{-1}(u)\big{]}=-a\big{(}u^{2}+2uq_{1}+q_{2}^{2}+4q_{1}^{2}\big{)}.$

Substituting this expression into $\mathcal{L}_{V}(u)$ one gets a Lax matrix
for the Hénon–Heiles system on the plane with the Hamiltonian

$\displaystyle
H_{1}=\frac{1}{2}\big{(}p_{1}^{2}+p_{2}^{2}\big{)}+8aq_{1}\big{(}2q_{1}^{2}+q_{2}^{2}\big{)}.$

In nonholonomic case we can also find similar Lax matrices, classical
$r$-matrix and integrable potentials.

######  Proposition 3.2.

For the Heisenberg type systems equations of motion (2.7) when
$V(q_{1},q_{2})=0$ can be rewritten in the Lax form (3.5) if

$\displaystyle\hat{\mathcal{L}}(u)=\left(\begin{matrix}\hat{h}&\hat{e}\\\
\hat{f}&-\hat{h}\end{matrix}\right)(u)\equiv\left(\begin{matrix}\dfrac{1}{2\mu}\dfrac{d\hat{e}(u)}{dt}&\hat{e}(u)\vspace{1mm}\\\
-\dfrac{d}{\mu
dt}\left(\dfrac{1}{2\mu}\dfrac{d\hat{e}(u)}{dt}\right)-2\hat{H}_{1}\hat{e}(u)&-\dfrac{1}{2\mu}\frac{d\hat{e}(u)}{dt}\end{matrix}\right)$
(3.7)

and

$\displaystyle\hat{\mathcal{A}}=\left(\begin{matrix}0&\mu\\\
-2\mu\hat{H}_{1}&0\end{matrix}\right).$

Here $\mu=\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}^{-1}$ is a conformal factor
(2.7) and function

$\displaystyle\hat{e}(u)=\frac{\big{(}u-\hat{\zeta}_{1}\big{)}\big{(}u-\hat{\zeta}_{2}\big{)}}{u}=u-\frac{2q_{1}}{1+q_{1}^{2}+q_{2}^{2}}-\frac{q_{2}^{2}}{u\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}}$

is completely defined by characteristic coordinates $\hat{\zeta}_{1,2}$
(2.21).

The Lax matrix (3.7) was constructed by using the generic construction of the
Lax matrices for the Stäckel systems [13]. The modification consists only of
application of the nontrivial conformal factor $\mu$.

Evaluating Poisson brackets $\\{\cdot,\cdot\\}^{\wedge}$ (2.5) between entries
of the Lax matrix $\hat{\mathcal{L}}(u)$ (3.7) we find the corresponding
classical $r$-matrix

$\displaystyle\hat{r}_{12}(u,v)=r_{12}(u,v)+2\left(\begin{matrix}0&0&0&0\\\
0&0&\hat{e}(v)&0\\\ 0&0&0&0\\\ -\hat{f}(v)&0&0&0\end{matrix}\right),$ (3.8)

which explicitly depends on dynamical variables via entries of
$\hat{\mathcal{L}}(u)$ (3.7).

As above, we can consider additive perturbation of this Lax matrix

$\displaystyle\hat{\mathcal{L}}_{V}(u)=\hat{\mathcal{L}}(u)+\left(\begin{matrix}0&0\\\
\Delta\hat{f}&0\end{matrix}\right),\qquad\Delta\hat{f}=\big{(}u^{2}-1\big{)}\big{[}\phi(u)\hat{e}^{-1}(u)\big{]}$

in order to get Lax matrices for the nonholonomic Heisenberg type systems with
$V\neq 0$. For instance, when $\phi(u)=-au^{3}$ the Teylor part of expansion
about $u=\infty$ looks like

$\displaystyle\big{[}{-}au^{3}\hat{e}^{-1}(u)\big{]}=-a\left(u^{2}+\frac{2uq_{1}}{1+q_{1}^{2}+q_{2}^{2}}+\frac{q_{2}^{2}}{{1+q_{1}^{2}+q_{2}^{2}}}+\frac{4q_{1}^{2}}{\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}^{2}}\right).$

Substituting this expression into $\hat{\mathcal{L}}_{V}(u)$ one gets a Lax
matrix for the nonholonomic counterpart of the Hénon–Hieles system on the
plane with the Jacobi integral

$\displaystyle\hat{H}_{1}=\frac{1}{2}\big{(}p_{1}^{2}+p_{2}^{2}+(q_{2}p_{1}-q_{1}p_{2})^{2}\big{)}+\frac{8aq_{1}\big{(}q_{1}^{2}q_{2}^{2}+q_{2}^{4}+2q_{1}^{2}+q_{2}^{2}\big{)}}{\big{(}1+q_{1}^{2}+q_{2}^{2}\big{)}^{3}}.$

##  4 Conclusion

Imposing nonholonomic constraints to Hamiltonian systems and making a suitable
reduction one gets some special class of non Hamiltonian systems on the
reduced phase space, see book [7]. Usually we can not investigate the reduced
system using standard mathematical methods of Hamiltonian dynamics. For
instance, we do not know how to get Lax matrices, classical $r$-matrices or
compatible Poisson brackets for the conformally Hamiltonian systems associated
with the nonholonomic Chaplygin ball, nonholonomic Suslov or Veselova systems.

In this note we find Killing tensors and compatible Poisson brackets, describe
integrable potentials and characteristic coordinates, evaluate Stäckel
matrices and Stäckel quadratures, to show Lax matrices and classical
$r$-matrices for the nonholonomic Heisenberg type systems. Indeed, we prove
that some modern geometric methods of Hamiltonian mechanics can be directly
applied to the reduced conformally Hamiltonian systems. Other similar examples
can be found in [6, 8, 20, 22, 23].

Of course, we can impose other nonholonomic constraints on the original
Hamiltonian system (1.1), (1.2). It may be interesting to describe all the
constraints which lead to nontrivial integrable metrics and potentials on the
reduced phase space.

### Acknowledgements

We are very grateful to the referees for thorough analysis of the manuscript,
constructive suggestions and proposed corrections, which certainly lead to a
more profound discussion of the results. We are also deeply grateful A.V.
Borisov and I.A. Bizayev for the relevant discussion.

Section 2 was written by A.V. Tsiganov and supported by the Russian Science
Foundation (project 15-12-20035). Section 3 was written by Yu.A. Grigoryev and
A.P. Sozonov within the framework of the Russian Science Foundation (project
15-11-30007).

## References

  * [1] Ankiewicz A., Pask C., The complete Whittaker theorem for two-dimensional integrable systems and its application, J. Phys. A: Math. Gen. 16 (1983), 4203–4208.
  * [2] Bates L., Śniatycki J., Nonholonomic reduction, Rep. Math. Phys. 32 (1993), 99–115.
  * [3] Benenti S., Orthogonal separable dynamical systems, in Differential Geometry and its Applications (Opava, 1992), Math. Publ., Vol. 1, Editors O. Kowalsky, D. Krupka, Silesian University Opava, Opava, 1993, 163–184.
  * [4] Bertrand J.M., Mémoire sur quelques-unes des forms les plus simples que puissent présenter les intégrales des équations différentielles du mouvement d’un point matériel, J. Math. Pures Appl. 2 (1857), 113–140.
  * [5] Bizyaev I.A., Borisov A.V., Mamaev I.S., Hamiltonization of elementary nonholonomic systems, Russ. J. Math. Phys. 22 (2015), 444–453, arXiv:1601.00884.
  * [6] Bizyaev I.A., Tsiganov A.V., On the Routh sphere problem, J. Phys. A: Math. Theor. 46 (2013), 085202, 11 pages, arXiv:1210.7903.
  * [7] Bloch A.M., Nonholonomic mechanics and control, Interdisciplinary Applied Mathematics, Vol. 24, Springer-Verlag, New York, 2003.
  * [8] Borisov A.V., Mamaev I.S., Symmetries and reduction in nonholonomic mechanics, Regul. Chaotic Dyn. 20 (2015), 553–604.
  * [9] Borisov A.V., Mamaev I.S., Bizyaev I.A., The hierarchy of dynamics of a rigid body rolling without slipping and spinning on a plane and a sphere, Regul. Chaotic Dyn. 18 (2013), 277–328.
  * [10] Borisov A.V., Mamaev I.S., Tsiganov A.V., Non-holonomic dynamics and Poisson geometry, Russ. Math. Surv. 69 (2014), 481–538.
  * [11] Cushman R., Duistermaat H., Śniatycki J., Geometry of nonholonomically constrained systems, Advanced Series in Nonlinear Dynamics, Vol. 26, World Sci. Publ. Co. Pte. Ltd., Hackensack, NJ, 2010.
  * [12] Darboux G., Sur un probléme de mécanique, Arch. Néerl. 6 (1901), 371–376.
  * [13] Eilbeck J.C., Enol’skii V.Z., Kuznetsov V.B., Tsiganov A.V., Linear $r$-matrix algebra for classical separable systems, J. Phys. A: Math. Gen. 27 (1994), 567–578, hep-th/9306155.
  * [14] Eisenhart L.P., Separable systems of Stäckel, Ann. of Math. 35 (1934), 284–305.
  * [15] Fassò F., Sansonetto N., Conservation of energy and momenta in nonholonomic systems with affine constraints, Regul. Chaotic Dyn. 20 (2015), 449–462, arXiv:1505.01172.
  * [16] Grigoryev Yu.A., Tsiganov A.V., Symbolic software for separation of variables in the Hamilton–Jacobi equation for the $L$-systems, Regul. Chaotic Dyn. 10 (2005), 413–422, nlin.SI/0505047.
  * [17] Kozlov V.V., The Euler–Jacobi–Lie integrability theorem, Regul. Chaotic Dyn. 18 (2013), 329–343.
  * [18] Molina-Becerra M., Galán-Vioque J., Freire E., Dynamics and bifurcations of a nonholonomic Heisenberg system, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 22 (2012), 1250040, 14 pages.
  * [19] Smirnov R.G., On the classical Bertrand–Darboux problem, J. Math. Sci. 151 (2008), 3230–3244, math-ph/0604038.
  * [20] Tsiganov A., Integrable Euler top and nonholonomic Chaplygin ball, J. Geom. Mech. 3 (2011), 337–362, arXiv:1002.1123.
  * [21] Tsiganov A.V., On bi-integrable natural Hamiltonian systems on Riemannian manifolds, J. Nonlinear Math. Phys. 18 (2011), 245–268, arXiv:1006.3914.
  * [22] Tsiganov A.V., One family of conformally Hamiltonian systems, Theoret. Math. Phys. 173 (2012), 1481–1497, arXiv:1206.5061.
  * [23] Tsiganov A.V., On integrable perturbations of some nonholonomic systems, SIGMA 11 (2015), 085, 19 pages, arXiv:1505.01588.
  * [24] van der Schaft A.J., Maschke B.M., On the Hamiltonian formulation of nonholonomic mechanical systems, Rep. Math. Phys. 34 (1994), 225–233.
  * [25] Whittaker E.T., A treatise on the analytical dynamics of particles and rigid bodies, Cambridge Mathematical Library, Cambridge University Press, Cambridge, 1988.


# Overview and Comparison of Nonlinear Interference Modelling Approaches in
Ultra-Wideband Optical Transmission Systems

Daniel Semrau1, Robert I. Killey1, Polina Bayvel1
1Optical Networks Group, Department of Electronic and Electrical Engineering,
UCL (University College London), Torrington Place, London WC1E 7JE, United
Kingdom
e-mail<EMAIL_ADDRESS>

###### Abstract

The recent advances in modelling nonlinear interference of systems operating
beyond the C-band are discussed. Estimation accuracy as well as computational
complexity of current approaches are compared and addressed.

###### keywords:

optical fibre communication, ultra-wideband transmission, nonlinear
interference, inter-channel stimulated Raman scattering

##  1 INTRODUCTION

Lightwave transmission extending beyond the C-band, exploiting the entire
bandwidth of installed fibres, is increasingly being considered as a cost-
effective alternative to deploying new multi-core/-mode fibres. However, for
such ultra-wide bandwidths, the non-instantaneous nature of the nonlinear
fibre response becomes significant, giving rise to inter-channel stimulated
Raman scattering (ISRS). ISRS effectively transfers power from high to lower
frequencies within the same optical signal. Although the physics of the
interaction between the Kerr effect and ISRS are well understood [1, 2],
efficient low-complexity performance models for modern coherent systems were
(until recently) not available. Such models are key for efficient link design,
real-time optimisation and physical layer-aware networking. Recently, much
attention has been drawn to the development of such models, particularly by
extending the modulation format independent Gaussian noise (GN) model to
account for ISRS [3, 4, 5, 6, 7, 8, 9, 10, 11]. Extensions of the closed-form
formalism for arbitrary, non-Gaussian modulation formats have been reported in
[12]. Experimental demonstrations on the interaction between Kerr effect and
ISRS followed up, in order to validate the theoretical predictions [5, 6, 13,
14]. To date, a number of modelling approaches exist in integral as well as in
closed-form, varying in accuracy and complexity.

In this paper, an overview of the recently proposed approaches of modelling
nonlinear interference in ultra-wideband transmission is presented. The
estimation accuracies are compared to split-step simulations and their
mathematical and computational complexity are briefly discussed.

##  2 MODELLING APPROACHES IN INTEGRAL FORM

The GN model is a first-order solution with respect to the nonlinearity
coefficient of the nonlinear Schrödinger equation assuming Gaussian modulation
[15]. The nonlinear interference (NLI) power is given as $P_{\text{NLI}}=\eta
P_{i}^{3}$, with NLI coefficient $\eta$ and channel launch power $P_{i}$. To
extend the conventional GN model to account for ISRS, it is assumed that the
temporal gain dynamics of ISRS are negligible. This is motivated by the
averaging of many independently modulated channels involved in the scattering
process, smoothing the occurring ISRS gain in time. As a resulting
simplification, ISRS can be modelled as a frequency- and distance-dependent
signal power profile $\rho\left(z,f\right)$, which is obtained by solving the
continuous-wave Raman equations [2]. The validity of this approach has
recently been experimentally demonstrated in [13, 14]. In general, approaches
to include ISRS in the conventional GN model can be categorised into two
groups. The first approach, termed the effective attenuation approach, is to
approximate $\rho\left(z,f\right)$ with exponential decays, that have modified
attenuation coefficients or effective lengths [3, 5]. The second approach,
termed the ISRS GN model, is to fully rederive the conventional GN model based
on the exact signal power profile for higher accuracy but with increased
complexity [4, 5, 6, 7, 8].

###  2.1 The effective attenuation approach

The first approach to extend the conventional GN model for ISRS, was the
introduction of exponential decays with channel-dependent attenuation
coefficients or effective lengths to model the effect of ISRS [3, 5]. The
advantage of the approach is that the expressions from the conventional, and
moderately complex, GN model, can be used. In [3], the use of channel-
dependent, effective attenuation coefficients $\alpha_{\textnormal{eff},i}$
was proposed, that resemble the actual effective length present in the fibre
span, as

$\displaystyle\begin{split}L_{\textnormal{eff},i}=\displaystyle\int_{0}^{L}\rho(\zeta,f_{i})d\zeta=\frac{1-\textnormal{exp}(-\alpha_{\textnormal{eff},i}L)}{\alpha_{\textnormal{eff},i}},\end{split}$
(1)

where $f_{i}$, $L_{\text{eff},i}$ are frequency and effective length of
channel $i$ and $L$ is the span length. The approach is valid for lumped
amplification in the weak ISRS regime, where the signal power profile is well
approximated by an exponential decay. The advantage of the approach is its
minor additional complexity with respect to the conventional GN model. The
only additional complexity, to model ISRS, is to solve the Raman equations and
then perform regression operations to obtain the coefficients
$a_{\text{eff},i}$. Both can be carried out within seconds, yielding a low
complexity prediction model. However, for large ISRS power transfers, the
signal power profile is not accurately modelled by exponential decays and
approximation errors are expected.

###  2.2 The ISRS GN model in integral-form

In order, to precisely account for arbitrary ISRS power transfers, as well as
for distributed amplification, the GN model has been rederived to account for
an arbitrary signal power profile [4, 5, 6, 7] as

$\begin{split}&\eta=\frac{16\gamma^{2}B_{i}}{27P_{i}^{3}}\int df_{1}\int
df_{2}\ S_{\text{Tx}}(f_{1},f_{2},f_{i})\left|\int_{0}^{L}d\zeta\
\sqrt{\frac{\rho(\zeta,f_{1})\rho(\zeta,f_{2})\rho(\zeta,f_{1}+f_{2}-f_{i})}{\rho(\zeta,f_{i})}}e^{j\phi\left(f_{1},f_{2},f_{i},\zeta\right)}\right|^{2},\end{split}$
(2)

where
$\phi=-4\pi^{2}(f_{1}-f_{i})(f_{2}-f_{i})\left[\beta_{2}+\pi\beta_{3}(f_{1}+f_{2})\right]\zeta$,
$S_{\text{Tx}}=G_{\text{Tx}}(f_{1})G_{\text{Tx}}(f_{2})G_{\text{Tx}}(f_{1}+f_{2}-f_{i})$
with $G_{\text{Tx}}(f)$ being the power spectral density (PSD) of the signal
at the transmitter and $B_{i}$ is the channel bandwidth. Eq. (2) can be used
for multiple spans by interpreting $\rho\left(\zeta,f\right)$ as the signal
power profile of the entire link.

Eq. (2) is exact to first-order for arbitrary power profiles, enabling the
modelling of strong ISRS transfers and distributed amplification scenarios.
The computational complexity is higher than the effective attenuation approach
2.1, due to an additional integration dimension over distance $\zeta$.
Moreover, the Raman equations must be solved manually for every span along the
transmission link and substituted in (2).

For lumped amplification and optical bandwidths of up to $15$ THz, the ISRS GN
model can be written in analytical form, avoiding the necessity of solving the
Raman equations, easing implementation and providing more insight into the
underlying parameter dependencies. This is enabled by an approximate solution
of the Raman equations [2]. This dramatically simplifies the description of
ISRS which can then be parametrised by only a single parameter $C_{r}$, which
is a linear regression of the Raman gain spectrum. The NLI coefficient can
then be written as [7, 8]

$\eta=\frac{16\gamma^{2}B_{i}G_{1}(f_{i})}{27P_{i}^{3}}\int df_{1}\int
df_{2}\left|\sum_{k=1}^{n}\int_{0}^{L_{k}}d\zeta
S_{k}\left(f_{1},f_{2},f_{i}\right)\frac{\hat{P}_{k}e^{-\alpha\zeta-\hat{P}_{k}C_{\text{r}}L_{\text{eff}}(f_{1}+f_{2}-f_{i})}}{\int
G_{k}(\nu)e^{-\hat{P}_{k}C_{\text{r}}L_{\text{eff}}\nu}d\nu}e^{j\phi\left(f_{1},f_{2},f_{i},\tilde{L}_{k}+\zeta\right)}\right|^{2}$
(3)

where $G_{k}$ and $\hat{P}_{k}$ are the signal PSD and total power launched
into span $k$, $\tilde{L}_{k}$ is the accumulated distance of span $k$,
$L_{\text{eff}}=\frac{1-\exp\left(-\alpha\zeta\right)}{\alpha}$ and
$S_{k}=\sqrt[]{\frac{G_{k}(f_{1})G_{k}(f_{2})G_{k}(f_{1}+f_{2}-f_{i})}{G_{k}(f_{i})}}$.
Eq. (3) can be conveniently used to model the NLI in point-to-point
transmission with different launch power distributions per span due to non-
ideal gain equalisation as well as for complex ultra-wideband network
scenarios [8]. The computational complexity of (3) is comparable to that of
(2).

Figure 1: NLI coefficient, and its deviation, after 3 spans obtained by
integral approaches Eq. (1-3).

###  2.3 Comparison with split-step simulations

To compare the accuracy of the modelling approaches, a split-step simulation
was performed for $\text{119}\times\text{85}$ GBd Nyquist-spaced channels,
occupying the entire C+L band (10.11 THz), centered at
$\lambda_{\text{ref}}=1570$ nm. Transmission over a link comprising three
spans of standard single mode fibre was considered, with parameters
$\alpha=0.2\frac{\text{dB}}{\text{km}}$,
$D=18\frac{\text{ps}}{\text{nm}\cdot\text{km}}$,
$S=0.067\frac{\text{ps}}{\text{nm}^{2}\cdot\text{km}}$,
$\gamma=1.2\frac{\text{1}}{\text{W}\cdot\text{km}}$,
$A_{\text{eff}}=81.8\mu\text{m}^{2}$ and
$C_{r}=0.0236\frac{1}{\text{W}\cdot\text{km}\cdot\text{THz}}$. ISRS was
implemented by a frequency-dependent loss at each simulation step, according
to the power profile $\rho\left(f,\zeta\right)$ and ideal gain equalisation
was performed. A sequence length of $2^{17}$ symbols was considered and four
data realisations were averaged to increase simulation accuracy. In order to
accurately benchmark the proposed models, Gaussian symbols were used for
transmission. The NLI coefficient as a function of channel frequency and as a
function of ISRS power transfer is shown in Fig. 1. The ISRS power transfer is
defined as the sum of the ISRS gain/loss in decibel of the most outer WDM
channels. The ISRS GN model in semi-analytical form (2) and in analytical form
(3) have a negligible error with respect to the split-step simulations. The
effective attenuation approach is in very good agreement with simulations
despite the exponential decay approximation; with a maximum mismatch of 0.1 dB
for 4 dBm launch power. However, this mismatch increases for increasing ISRS
power transfers as ISRS does not strictly resemble an exponential decay.

In conclusion, the NLI estimates of all integral approaches are similar within
0.1 dB for powers up to 4 dBm per channel. However, for stronger ISRS power
transfers and higher required accuracy, the use of the ISRS GN model is
recommended at the expense of higher computational complexity.

##  3 MODELLING APPROACHES IN CLOSED-FORM

Although, the integral approaches above estimate the NLI with very good
accuracy, they are not directly suitable for real-time applications and
network optimisation where vast numbers of light paths must be evaluated. For
such scenarios, approximations in closed-form, that yield results within
picoseconds, have been proposed for the effective attenuation approach, as
well as for the ISRS GN model [3, 9, 10, 11, 12].

###  3.1 The effective attenuation approach

A semi-analytical closed-form approach was proposed in [3]. The approach is
based on a closed-form solution of the conventional GN model extended by the
effective attenuation as in Sec. 2.1. The NLI coefficient is approximated by
[3]

$\displaystyle\begin{split}\eta&\approx\frac{8}{27}\gamma^{2}\frac{n^{1+\epsilon}\alpha_{\textnormal{eff},i}L_{\textnormal{eff},i}^{2}}{\pi|\beta_{2,i}|B_{i}^{2}}\textnormal{asinh}\left(\frac{0.5\pi^{2}|\beta_{2,i}|B_{\textnormal{tot}}^{2}}{\alpha_{\textnormal{eff},i}}\right),\end{split}$
(4)

with $\beta_{2,i}$ being the GVD parameter of channel $i$. Eq. (4) is an
approximate solution for the NLI of the central channel and approximation
errors at the outer WDM channels are expected. The approach is semi-
analytically, as the Raman equations have to be solved and regression
functions have to be executed to obtain $\alpha_{\text{eff},i}$.

###  3.2 The ISRS GN model in closed-form

The first fully analytical closed-form approximation of the ISRS GN model was
derived in [9, 10]. The formula accounts for arbitrary launch power
distributions, channel configurations, ISRS and wavelength dependent
attenuation and dispersion. The approximation relies on a first-order
description of ISRS, introducing approximation errors for large ISRS power
transfers, and is given by [10]

$\begin{split}\eta&\approx\frac{4}{9}\frac{\gamma^{2}}{B^{2}_{i}}\frac{\pi
n^{1+\epsilon}}{\phi_{i}\bar{\alpha}\left(2\alpha+\bar{\alpha}\right)}\cdot\left[\frac{T_{i}-\alpha^{2}}{a}\text{asinh}\left(\frac{\phi_{i}B_{i}^{2}}{\pi
a}\right)+\frac{A^{2}-T_{i}}{A}\text{asinh}\left(\frac{\phi_{i}B_{i}^{2}}{\pi
A}\right)\right]\\\ &+\frac{32}{27}\sum_{k=1,k\neq
i}^{N_{\mathrm{ch}}}\left(\frac{P_{k}}{P_{i}}\right)^{2}\frac{\gamma^{2}}{B_{k}}\frac{n}{\phi_{i,k}\bar{\alpha}\left(2\alpha+\bar{\alpha}\right)}\left[\frac{T_{k}-\alpha^{2}}{\alpha}\mathrm{atan}\left(\frac{\phi_{i,k}B_{i}}{\alpha}\right)+\frac{A^{2}-T_{k}}{A}\
\mathrm{atan}\left(\frac{\phi_{i,k}B_{i}}{A}\right)\right],\end{split}$ (5)

with $\phi_{i}=\frac{3}{2}\pi^{2}\left(\beta_{2}+2\pi\beta_{3}f_{i}\right)$,
$T_{i}=\left(\alpha+\bar{\alpha}-P_{\text{tot}}C_{\text{r}}f_{i}\right)^{2}$
and
$\phi_{i,k}=2\pi^{2}\left(f_{k}-f_{i}\right)\left[\beta_{2}+\pi\beta_{3}\left(f_{i}+f_{k}\right)\right]$
and $A=\alpha+\bar{\alpha}$. A formula based on the same assumptions with a
similar result was presented in [11]. Eq. (5) enables real-time NLI estimation
in point-to-point links and mesh optical networks. Very recently, (5) has been
extended for arbitrary modulation formats [12].

###  3.3 Comparison with split-step simulations

The NLI coefficient as a function of channel frequency and as a function of
ISRS power transfer is shown in Fig. 2. While Eq. (4) shows some approximation
errors towards the outer WDM channels, it is reasonably accurate in
calculating the change of the NLI coefficient for ISRS power transfers of up
to 3 dB, with a maximum error of 0.2 dB compared to simulations. The ISRS GN
model in closed-form (5) shows very good agreement to simulations for all
channels, with an average deviation of 0.1 dB at 4 dBm launch power. The
maximum deviation of (5) compared to the result in [11] was negligible ($<$0.1
dB) over the entire parameter space studied.

Although, (4) offers some initial conclusions on the impact of ISRS on the
NLI, the recently derived ISRS GN model in closed-form (5) is superior in
computational complexity as well as in accuracy. Therefore, we recommend the
use of (5) for the real-time modelling of ultra-wideband transmission systems.

Figure 2: NLI coefficient, and its deviation, after 3 spans obtained by
closed-form approaches Eq. (4-5).

##  4 CONCLUSIONS

The recently proposed approaches to account for inter-channel stimulated Raman
scattering in nonlinear interference modelling are reviewed, in particular the
effective attenuation approach and the ISRS GN model; both in integral and in
closed-form. When the estimation accuracy is the priority, the ISRS GN model
in integral form or the less complex effective attenuation approach should be
considered. For time sensitive applications, the ISRS GN model in closed-form
is recommended for almost instantaneous, but yet accurate, NLI estimation.

## ACKNOWLEDGMENTS

Financial support is from UK EPSRC DTG PhD studentship to D. Semrau and EPSRC
TRANSNET Programme Grant and is gratefully acknowledged.

## References

  * [1] K. J. Blow et al.: Theoretical description of transient stimulated Raman scattering in optical fibers, J. Quantum Electron., 25, (12), pp. 2665–2673, December 1989
  * [2] D. N. Christodoulides et al.: Evolution of stimulated Raman crosstalk in wavelength division multiplexed systems, Photon. Technol. Lett., 8, (12), pp. 1722–1724, December 1996
  * [3] D. Semrau et al.: Achievable rate degradation of ultra-wideband coherent fiber communication systems due to stimulated Raman scattering, Optics Express, 25, (12), pp. 13024–13034, May 2017
  * [4] I. Roberts et al.: Channel power optimization of WDM systems following Gaussian noise nonlinearity model in presence of stimulated Raman scattering, J. Lightw. Technol., 35, (23), pp. 5237–5249, November 2017
  * [5] M. Cantono et al.: Modelling the impact of SRS on NLI generation in commercial equipment: an experimental investigation, Proc. OFC 2018, San Diego, United States, March, M1D.2
  * [6] M. Cantono et al.: On the interplay of nonlinear interference generation with stimulated Raman scattering for QoT estimation, J. Lightw. Technol., 36, (15), pp. 3131–3141, March 2018
  * [7] D. Semrau et al.: The Gaussian noise model in the presence of inter-channel stimulated Raman scattering, J. Lightw. Technol., 36, (14), pp. 3046–3055, April 2018
  * [8] D. Semrau et al.: The ISRS GN model, an efficient tool in modeling ultra-wideband transmission in point-to-point and network scenarios, Proc. ECOC 2018, Rome, Italy, September, Tu4G.6
  * [9] D. Semrau et al.: A closed-form approximation of the Gaussian noise model in the presence of inter-channel stimulated Raman scattering, arXiv, arXiv:1808.07940, August 2018
  * [10] D. Semrau et al.: A closed-form approximation of the Gaussian noise model in the presence of inter-channel stimulated Raman scattering, J. Lightw. Technol., 37, (9), pp. 1924–1936, January 2019
  * [11] P. Poggiolini: A generalized GN-model closed-form formula, arXiv, arXiv:1810.06545, September 2018
  * [12] D. Semrau et al.: A modulation format correction formula for the Gaussian noise model in the presence of inter-channel stimulated Raman scattering, arXiv, arXiv:1903.02506, March 2019
  * [13] G. Saavedra et al.: Inter-channel stimulated Raman scattering and its impact in wideband transmission systems, Proc. OFC 2018, San Diego, United States, Th1C.3
  * [14] K. Minoguchi et al.: Experiments on stimulated Raman scattering in S- and L-bands 16-QAM signals, for ultra-wideband coherent WDM systems, Proc. OFC 2018, San Diego, United States, March, Th1C.4
  * [15] P. Poggiolini et al.: The GN-Model of fiber non-linear propagation and its applications, J. Lightw. Technol., 32, (4), pp. 694–721, February 2014


# A unified Neural Network Approach to E-Commerce Relevance Learning

Yunjiang Jiang<EMAIL_ADDRESS>,  Yue Shang<EMAIL_ADDRESS>JD.com
Silicon Valley Research Center675 E Middlefield RoadMountain ViewCA94043USA ,
Rui Li<EMAIL_ADDRESS>,  Wen-Yun Yang<EMAIL_ADDRESS>JD.com
Silicon Valley Research Center675 E Middlefield RoadMountain ViewCA94043USA ,
Guoyu Tang<EMAIL_ADDRESS>,  Chaoyi Ma<EMAIL_ADDRESS>JD.comJD Building,
No. 18 Kechuang 11 Street, BDABeijing101111People’s Republic of China ,  Yun
Xiao<EMAIL_ADDRESS>and  Eric Zhao<EMAIL_ADDRESS>JD.com Silicon Valley
Research Center675 E Middlefield RoadMountain ViewCA94043USA

(2019)

###### Abstract.

Result relevance scoring is critical to e-commerce search user experience.
Traditional information retrieval methods focus on keyword matching and hand-
crafted or counting-based numeric features, with limited understanding of item
semantic relevance. We describe a highly-scalable feed-forward neural model to
provide relevance score for (query, item) pairs, using only user query and
item title as features, and both user click feedback as well as limited human
ratings as labels. Several general enhancements were applied to further
optimize eval/test metrics, including Siamese pairwise architecture, random
batch negative co-training, and point-wise fine-tuning. We found significant
improvement over GBDT baseline as well as several off-the-shelf deep-learning
baselines on an independently constructed ratings dataset. The GBDT model
relies on 10 times more features. We also present metrics for select subset
combinations of techniques mentioned above.

††journalyear: 2019††copyright: acmcopyright††conference: 1st International
Workshop on Deep Learning Practice for High-Dimensional Sparse Data; August 5,
2019; Anchorage, AK, USA††price: 15.00††doi: 10.1145/3326937.3341259††isbn:
978-1-4503-6783-7/19/08

##  1\. Introduction

Online shopping has become ubiquitous in the modern world. A good e-commerce
search engine, which helps users quickly zero in on their intended products
from billions of candidates, is thus essential to users’ online shopping
experiences.

The most fundamental problem in e-commerce search is learning relevance
between query and items. Compared with traditional web search or document
retrieval, relevance learning in e-commerce search has two unique
requirements.

  * •

Unlike web documents, e-commerce items typically have little business-side
textual content, mainly a short title of less than 50 characters; user
comments are less reliable and typically more skewed towards older products.
Thus simple keyword-based matching may not capture semantic matches
effectively. We observe that even serious users enter search queries with
words that do not necessarily conform to standard retail jargons. Thus, an
ideal model should be able to do fuzzy match between queries and product
titles.

  * •

Due to the profit-making nature of e-commerce, personalized recommendation and
other promotional logic often play a more important role than relevance in
final ranking. On the other hand, irrelevant items should clearly be filtered.
An ideal relevance model therefore should calibrate well with an established
relevance grading, such as a binary classification for relevant and irrelevant
ones or PEGFB scale ratings (Perfect, Excellent Good, Fair, Bad).

In this paper, we share our experiences of learning such an ideal relevance
model for e-commerce search. We choose a deep neural network based approach as
it can learn latent semantics based on raw text features. As a deep model
usually requires a huge amount of data to train and human labels are
expensive, we explore two sources of label supervision: 1) user clicks, an
indirect signal for relevance (i.e., users clicks are affected by many factors
such as relevance, price, and sale volumes), and 2) limited editorial
relevance ratings. Our main contribution is to demonstrate three highly
general techniques to enhance the quality of such a hybrid deep neural
network, to achieve accurate relevance prediction.

(a) Pairwise Siamese

(b) Point-wise ensemble fine-tuning

Figure 1. Model Architectures

###  1.1. Siamese Network for Pairwise Relevance

While user clicks are significantly correlated with result relevance
(huang2013learning, ; neural_ir_book, ), click derived signals, such as CTR,
are unsuitable as a point-wise learning target for relevance, especially in
e-commerce, because 1) they are not well-calibrated to absolute relevance
ratings, such as the widely used PEGFB scale, and 2) point-wise click labels
have very skewed distribution (typically concentrated at 0), making it easy
for the model to “cheat”. However, clicks in aggregate are invaluable at
distinguishing between good and bad results under the same query. For this
reason, we choose a pairwise architecture that learns preference between a
pair of items under the same query. To enable scoring single items during
serving, the model takes a Siamese form (bromley1994signature, ): the final
logit for the item pair is expressed as the difference of two individual item
scores from two towers of shared weights, and only one of the two identical
towers is used during serving. While most applications of Siamese network in
the literature are for pairwise similarity learning (chopra2005learning, ;
yin2016abcnn, ), we use it for pairwise discrimination. Experiments show that
this 2-tower setup actually beats the non-Siamese pairwise baseline on several
key metrics (Table 1).

###  1.2. Efficient Training with Batch Negatives

One challenge we face in using clicks as training labels is the paucity of
negative relevance examples. To ensure good user experience, top items in a
query session are usually at least somewhat relevant. Thus the model rarely
sees completely unrelated (query, item) pairs. This causes a serious
distributional skew during serving, where we intend to score thousands of
relevant and irrelevant items for a query. As a remedy, we design an efficient
batch negative algorithm (Algo 1), to supplement the original $n$ session
pairs in each batch with $n(n-1)$ additional random negative (query, item)
pairs. Intuitively, this gives the model a more global view of possible items
for each query. Amazingly all key metrics on the original examples improve
(Table 1). Case study also shows that unrelated items now score much lower
under the batch negative enhancement.

###  1.3. Fine-Tuning with Human Supervision

Lastly, we fine-tune the model with editorial labels, in a point-wise fashion.
This is needed for two reasons:

  * •

clicks, though abundant, can only approximate result relevance, and are often
misleading: a user can click on a result simply because of an attractive title
or being part of a malicious click farm to manipulate search ranking.

  * •

the ultimate system goal is to filter irrelevant items (those below the Good
rating), but pairwise models are not well calibrated to an absolute notion of
relevance; e.g., a score of 0.0 could mean Excellent or Fair depending on
query, because the Siamese model is shift invariant.

Point-wise architecture at this step is conveniently adapted from a single
tower of the base Siamese model. This leads to an additional 6% gain in eval
AUC.

##  2\. Related Works

In this section, we review related work on relevance learning, which aims to
learn a function to determine relevance between a query and document pair. We
categorize these models as traditional learning to rank based methods and deep
learning based methods.

###  2.1. Traditional Learning to Rank Methods

Due to the importance of search engines, using model machine learning
techniques to improve relevance ranking, often referred as learning to rank,
has been an active research topic in the past few decades. A lot of progresses
have been made, including point-wise models like RankNet (L2RPoint, ),
pairwise models like RankSVM (ranksvm, ) and GBRank (gbrank, )), and list-wise
models like AdaRank (AdaRank, ) and LambdaMart (LambdaMart, )); see
(l2rSurvey, ) for a complete survey. Here we discuss some key results that our
work builds upon. First is using pairwise relevance preference instead of
absolute relevance grade as the learning target introduced in (ranksvm, ;
gbrank, ); this is the basis of the Siamese pairwise architecture. We do not
use a list-wise approach for optimizing ranking results globally, because
unlike web search, relevance is only a signal in final ranking and is mainly
used as filtering in e-commerce search. Second is exploiting easy-to-get but
noisy click logs instead of expensive but accurate human labels to learn
relevance (ranksvm, ). Our model actually leverages both labels. All the
traditional methods rely on manually defined features, such as BM25, matching
positions and page-rank, and cannot effectively utilize raw text features.
Thus, they often fail to evaluate _semantic_ relevance if a document does not
contain exact terms in a query.

###  2.2. Neural Network Based Ranking Methods

In the past few years, a large number of embedding-based neural network models
have been successfully applied to learn semantic relevance between queries and
documents; (INR-061, ) provides an excellent survey. These works are roughly
divided into two camps: using neural network induced embeddings as features in
a 2-stage process, or directly in an end-to-end manner.

In the first category, (van2016learning, ) learns embeddings as features for a
final ranking objective. (yu2014latent, ) uses the embeddings to promote
result diversity. The resulting embeddings are amenable to fast online
retrieval and similarity analysis. (ai2017learning, ) in addition learns user
embeddings for personalized search.

In the second category, DSSM (huang2013learning, ) and their followup work
CDSSM (cdssm, ) pioneered the application of deep neural networks in end-to-
end relevance learning. Further, new models, such as DRMM (guo2016deep, ),
Duet (duet, ), DeepRank (deep_rank, ) have been proposed to improve content
based models via exploring traditional IR lexical matching signals (e.g.,
query terms importance, exact matching) in neural networks. This is the
approach that we follow. Also we keep the user dimension out of the equation
to stay focused on pure relevance learning.

DSSM uses two separate feed-forward towers, also referred to as Siamese
network, to encode query and result as embeddings, and optimizes a softmax
loss based on the cosine distance of those embeddings. Our notion of Siamese
network is different. Each of our training examples consists of a query and a
pair of items, where the model allows interaction between the query and either
item, but no interaction between the two items (see Figure 11(a)). We note
that, our model can be easily extended to include those new signals beyond
just query and title, or new architectures (e.g., RNN, Transformer) in each
tower of the Siamese network. Our main contributions are three general
techniques: 1) Siamese network for learning pairwise preference, 2) efficient
training with batch negatives, and 3) point-wise fine-tuning via human labels.

##  3\. Model Description

Much of the relevance ranking literature takes a single query and multiple
documents as input, such as DSSM and its variants. This is well suited for the
ranking problem. Our ultimate use case, however, is to filter out irrelevant
documents, based on a score and a global threshold. In other words, we learn a
relevance classification instead of a ranking model.

For simplicity and portability, we use title text as the only item side
feature. In practice, numeric features can be sparse, especially for tail
queries, whereas item titles are almost always available, and can help better
generalize from torso queries. Our model is first trained on user clicks, then
fine-tuned using human labels.

###  3.1. Siamese Pairwise Network

Due to the abundance of clicks in search log, compared to other user signals
such as purchase or placement into shopping cart, we choose clicks as weak
labels in the initial model warm-up step.

Since clicks are most informative of relative relevance, we first learn
preference between a pair of items under each query: each training example
consists of the triple ($Q$, $\rm{item}_{a}$, $\rm{item}_{b}$), where $Q$
stands for query, and $\rm{item}_{a}$ and $\rm{item}_{b}$ are two co-occurring
items under $Q$, swapped randomly. We call such triples session pairs.

The goal of the model is to predict the click ratio of $\rm{item}_{a}$, as a
proxy for its relative relevance, that is, we train a classifier $F$ according
to

(1)
$F(Q,\rm{item}_{a},\rm{item}_{b})\sim\frac{\rm{click\\_count}_{a}}{(\rm{click\\_count}_{a}+\rm{click\\_count}_{b})}$

The obvious simpler alternative is to instead fit a function G to the
pointwise click-through rate objective

(2) $G(Q,\rm{item})\sim\frac{\rm{click\\_count}}{\rm{view\\_count}}$

However this approach suffers badly from label imbalance in the e-commerce
search setting, where the click through labels are typically very sparse and
have a highly skewed asymmetric distribution. For instance, letting $G$ be
identically 0 would already achieve a reasonably good logloss or square loss,
which is completely useless for relevance scoring.

By randomly swapping $\rm{item}_{a}$ and $\rm{item}_{b}$, the function F
cannot rely on any distributional prior to beat random guessing. Since the
problem can be viewed as binary classification, logloss is a natural choice.

For the convenience of scoring a single item during serving, we further
restrict F to the factorized Siamese form:

(3) $\displaystyle
F(Q,\rm{item}_{a},\rm{item}_{b})=H(Q,\rm{item}_{a})-H(Q,\rm{item}_{b})$

Taking difference seems to work reasonable well compared to other more
elaborate contrastive functions. Since the click ratio labels are not
necessarily in $\\{0,1\\}$, hinge loss does not apply meaningfully to the
within-session pair examples. However it does make sense to apply hinge loss
to the batch negative sampled pairs. We found the hybrid session logistic-loss
and batch negative hinge-loss approach to improve eval logistic-loss and auc
significantly on the session pairs. However it did not improve key metrics
during the human-labelled fine-tuning step.

The single tower function H computes the following:

  1. (1)

For both query and item, we sum up the unigram/bigram embedding vectors (of
the same dimension), divided by sqrt(n), then concatenated to form a 2d tensor
E of shape (n, 2d), where n is the training/eval batch size, and d is the
embedding dimension.

  2. (2)

The prediction in logit space is then calculated as $\rm{H}^{\prime}(E)$,
where $\rm{H}^{\prime}$ is a feed-forward neural net, consisting of 3 ReLU
layers [1024, 256, 64, 1].

###  3.2. Efficient Training with Batch Negatives

input : Batch size: $n$

input : Feed-Forward network: $H^{\prime}$

input : Batched embedding tensors for queries, positive items, negative items:
$Q$, $I_{+}$, $I_{-}$

output : Logits for the original and batch negative pairs

$E_{+}\leftarrow\text{ColumnConcat}(Q,I_{+})$

$E_{-}\leftarrow\text{ColumnConcat}(Q,I_{-})$

$E\leftarrow\text{RowConcat}(E_{+},E_{-})$

foreach  _$1 <=i<=n$_ do

$I_{+}\leftarrow\text{RowCyclicPermute}(I_{+})$

$E_{\rm{BN}}\leftarrow\text{RowConcat}(Q,I_{+})$

$E\leftarrow\text{RowConcat}(E,E_{\rm{BN}})$

end foreach

$\vec{P}\leftarrow\text{RowEvenSplit}(H^{\prime}(E),n+1)$

original logits $\leftarrow P_{1}-P_{2}$

BN logits $\leftarrow\text{RowConcat}(P_{3}-P_{1},\ldots,P_{n+1}-P_{1})$

All logits $\leftarrow\text{RowConcat}(\text{original logits},\text{BN
logits})$

return _$\rm{All}\text{ }\rm{logits}$_

Algorithm 1 Batch negative augmented forward pass.

The Siamese pairwise model trained using session pairs only, as described in
the previous section, already exceeds DSSM on several evaluation metrics
(Table 1). More careful case study, however, reveals a fundamental problem.
Many items have little variation in their scores regardless of which query is
used. Clearly query is an important input to any relevance model, and our goal
is not to create a query-independent item quality model.

One explanation is that unlike web search, e-commerce item attractiveness
often plays a much more important role than relevance in determining users’
feedback. Truly irrelevant items (negative examples) are also rare in a
session, since our search system has optimized its relevance for years and
most of returned items are already relevant.

To address these two problems, we borrow ideas from word2vec (word2vec, ), and
augment the training data with $O(n^{2})$ soft negative pairs, constructed
within the same batch. In the simplest point-wise setting, we pair up the
query of each example in the batch with items of other examples. In the
pairwise case, a natural modification is to pair the query and positive item
of each example with the positive items of other examples. The motivation is
that since the training set is globally uniformly shuffled, the chance that
two examples within a single mini-batch have the same query, or even
synonymous queries, is extremely low. Thus we can assign with high confidence
negative labels to each new pair.

For $a,b$ two items co-occurring under a query $q$, let $\rm{click}(a;b,q)$
stand for the number of clicks item $a$ received in the context of $q$ and
$b$. Further, denote by $\Lambda$ the set of all such co-occurring triples
$(a,b,q)$ within a mini-batch, ordered by
$\rm{click}(a;b,q)\geq\rm{click}(b;a,q)$, followed by alphabetic tie-breaking
of the item titles.

Then the mini-batch training objective takes the form
$\rm{trainloss}=\text{original loss}+\text{BN loss}$, where

(4) original loss
$\displaystyle=\sum_{(a,b,q)\in\Lambda}\lambda(a,b,q;\ell(a,b,q))$ (5) BN loss
$\displaystyle=\sum_{(a,b,q),(a^{\prime},b^{\prime},q^{\prime})\in\Lambda,q^{\prime}\neq
q}\lambda(a,a^{\prime},q;0)$ (6) $\displaystyle\ell(a,b,q)$
$\displaystyle=\frac{\rm{click}(a;b,q)}{\rm{click}(a;b,q)+\rm{click}(b;a,q)}$
(7) $\displaystyle\lambda(a,b,q;\ell)$
$\displaystyle=\tau(\rm{logit}(a,q)-\rm{logit}(b,q),\ell),$

and

(8) $\displaystyle\tau(x,\ell)=-\ell\log\sigma(x)-(1-\ell)\log\sigma(-x)$

is the logloss function for logit $x\in\mathbb{R}$ and label $\ell\in[0,1]$.
Here $\sigma(x):=\frac{1}{1+e^{-x}}$ is the sigmoid function.

Thus, the batch negative (BN) loss is the sum of $n(n-1)$ terms, where $n$ is
the batch size, because each positive item can be paired with all positive
items from the other $n-1$ examples. The relative weights between the two loss
components do not seem to matter, as judged by their final eval metrics.
Algorithm 1 show the implementation of computing the forward pass with
$O(n^{2})$ batch negative examples efficiently within the tensorflow
framework.

Geometrically, with batch negatives, the query embeddings are not allowed to
cluster around a single direction, since unrelated (query, item) pairs must
score much lower than related pairs.

On the 6 month data, batch negative AUC reaches 0.99 within 100k steps of 128
batch size; in contrast, original example AUC takes 6 million steps to
converge. That is because the batch is now dominated by $O(n^{2})$ soft
negative examples, which are easy to learn. The inclusion of batch negative
examples has little effect on how quickly metrics converge on the original
examples, making weight tuning unnecessary.

###  3.3. Fine Tuning via Relevance Ratings Data

In the last stage, we use human labeled data to fine-tune the click model.
Since rating millions of items is prohibitively expensive, it is impractical
to train a model from scratch. Instead we take the click model output and its
word embeddings as a starting point, and learn a new set of hidden layers from
random initialization. In order to assign absolute relevance grade to each
(query, item) pair, we also switch from pairwise to a point-wise architecture,
with binary labels determined by whether the rating exceeds a certain
threshold: Perfect/Excellent/Good ratings are considered positive; see Figure
1 (1(b)) for the model architecture. Finally we add the new point-wise network
to the original click point-wise prediction, as a form of ensembling.

The point-wise transfer learned model thus predicts the probability whether a
(query, result) pair has a relevance score above Fair (2). More generally,
experiments (Table 2) show that the point-wise approach is better suited for
optimizing Precision/Recall AUC for both positive and negative examples.

models \ metrics | Eval AUC | Pair Accu | Neg PR AUC | PR AUC | NDCG@10 | Mean Avg Prec | Prec@3
---|---|---|---|---|---|---|---
Click-Trained Model (part) | 0.5846 | 0.6645 | 0.6065 | 0.8316 | 0.8947 | 0.8721 | 0.7351
w/o Negative Sampling | 0.5825 | 0.6476 | 0.5968 | 0.8314 | 0.8912 | 0.8680 | 0.7299
non-Siamese Pairwise | 0.5751 | 0.6071 | - | - | - | - | -
DSSM | 0.5692 | 0.5862 | 0.5288 | 0.6495 | 0.8681 | 0.8615 | 0.7114
Table 1. Click model comparison with simplified versions models \ metrics |  | Eval AUC
---
| Pair Accu
---
| Neg PR AUC
---
| PR AUC
---
NDCG@10 | Mean Avg Prec | Prec@3
GBDT | - | 0.6778 | 0.6463 | 0.7759 | 0.8995 | 0.8933 | 0.7225
Click-Trained Model (full) | 0.7572 | 0.6822 | 0.4832 | 0.7553 | 0.9053 | 0.9054 | 0.7240
| +Pointwise simple
---
0.8201 | 0.6678 | 0.6598 | 0.7656 | 0.8956 | 0.8942 | 0.7220
| +Pairwise ensemble
---
0.7744 | 0.6956 | 0.6546 | 0.7673 | 0.9077 | 0.9094 | 0.7232
| +Pointwise ensemble
---
0.8262 | 0.6921 | 0.6698 | 0.7779 | 0.9023 | 0.9003 | 0.7244
Table 2. Model compare results with GBDT and difference transfer learning configurations | No Batch Negative Model | With Batch Negative Model |
---|---|---|---
Item title | “cellphone” | “paper cup” | “cellphone” | “paper cup” |
| 荣耀10 64GB 渐变蓝 双摄像头 双卡双待 4G全面屏手机
---
(Honor 10 64GB Blue Dual Camera Dual Sim Full Display)
1.6292 | 1.2842 | 4.6584 | -8,2954 |
| Apple iPhone X 64GB 深空灰色
---
(Apple iPhone X 64GB Space Grey)
1.5081 | 1.1767 | 4.8843 | -7.1850 |
| 荣耀9i 4GB+64GB 幻夜黑 4G全面屏手机 双卡双待
---
(Honor 9i 4GB+64GB Black 4G Full Display Dual SIM)
1.5014 | 1.0814 | 4.7038 | -9.0472 |
Table 3. Siamese pairwise click model with/without batch negative for the
queries “cellphone” and “paper cup”

##  4\. Experiments

We present two sets of experiments. The first set compares models trained on
10% of 6 month session log (about 50m examples), and 1 million unigram/bigram
vocab. No human labels are used here. The second set focuses on the effect of
human label transfer learning. The common base model here is trained on the
full 6 month of user click data, with a vocab size of about 16 million.

###  4.1. Experiments Setup

####  4.1.1. Click Data Generation

Our primary training data consists of user click-through on search session
items. We perform several compression steps to strike a balance between
session coverage and data preservation:

  1. (1)

From each search session, generate (query, $\rm{item}_{a}$, $\rm{item}_{b}$,
$\rm{is\\_clicked}_{a}$, $\rm{is\\_clicked}_{b}$) 5-tuples, where at least one
item is clicked and below the other item in displayed position.

  2. (2)

Over a period of 180 days, aggregate the click counts to obtain the 5-tuple’s
(query, $\rm{item}_{a}$, $\rm{item}_{b}$, $\rm{click\\_cnt}_{a}$,
$\rm{click\\_cnt}_{b}$).

  3. (3)

For each query, we retain at most top 100 5-tuple’s of the previous step,
sorted by $\rm{click\\_cnt}_{a}+\rm{click\\_cnt}_{b}$.

The click data is randomly shuffled by query, and we use 90% for training and
10% for eval.

####  4.1.2. Rating Data Generation

We collect human labelled data totaling about 300k (query, item) pairs. On
average, each query has about 5 items. The queries are shuffled at random and
split into 65% training, 30% eval and 5% test, from which we construct the
pointwise datasets. The training and eval data sets are used for transfer
learning. All the experiments are evaluated on test set.

###  4.2. Evaluation Metrics

Besides the standard IR metrics NDCG@10, MAP, and Precision@3, we calculate
pairwise accuracy (Pair Accu) for two items under the same query and
Precision/recall AUC (PR-AUC) for both positive and negative labels. PR-AUC is
more discriminating than ROC-AUC for tasks with a large skew in the class
distribution (davis2006relationship, ). In addition, it’s a natural
generalization of classification accuracy, which we care about for the
filtering task. All experiments are evaluated on a 13k (query, item) pair test
set.

###  4.3. Experiment Results

According to the previous discussion, we proposed three general techniques for
deep relevance training, including Siamese structure for pairwise relevance
learning, efficient batch negative training, and model fine tuning through
human labeled point-wise data. To study the effect of these techniques, we
compare our model with several simpler versions.

The first set of experiments aim to evaluate the techniques proposed for
training a model based on click data only. Particularly, we evaluate the
following methods.

  * •

Base Click Model (part) is the one trained on 10% dataset, with both Siamese
and batch negative techniques.

  * •

w/o Siamese is a model without the Siamese structure, which takes a pairwise
preference example ($Q$, $\rm{item}_{a}$, $\rm{item}_{b}$) as input, and
outputs their relevance preference. As the model is trained to score a
pairwise preference example, it’s hard to evaluate point-wise relevance
metrics(e.g., PR AUC) and we just report eval AUC and pairwise accuracy.

  * •

w/o Batch Negative is our base click model without batch negative examples.

  * •

DSSM is a classical deep matching model for web search (huang2013learning, ).
DSSM could be regarded as a simplified version of our model without the
Siamese structure and batch negatives.

By comparing the results shown in Table-1, the base click model boosts the
performance on all metrics (outperforming DSSM by 13% for pairwise accuracy,
28% for precision-recall AUC, 13% for negative PR-AUC and 3% for NDCG),
showing that our proposed general techniques facilitate the learning of
relevance from click data and decrease the impact of distributional
mismatches.

In the second set of experiments, we compare the proposed transfer learning
technique with several variants and baselines.

  * •

GBDT (Gradient Boosted Decision Tree) is one of the most popular non-DL based
models used in relevance learning ; see (yin2016ranking, ) for detail. This is
our in-house ranking model baseline, implemented with Xgboost (xgboost, ),
with 60 features covering a variety of user feedback signals, hand-picked
matching features, including BM25, window size, etc. This is optimized for
years and is therefore a strong baseline.

  * •

Base Click Model (full) is trained on the whole dataset, with batch negative
and Siamese structure, but without transfer learning. We take it as another
baseline to highlight the effect of transfer learning.

  * •

+Pointwise ensemble is our propose method.

  * •

+Pairwise ensemble is similar to our proposed approach but uses pairwise human
label training target.

  * •

+Pointwise simple is similar to the proposed method but without the base click
network.

As shown in Table 2, the transfer learning models outperform GBDT and the base
click model in all the metrics.

Ensembling improves the point-wise transfer learning model on all metrics.
Pairwise ensemble on the other hand does the best on pairwise accuracy on the
validation set, as expected. However, since PR-AUC and negative PR-AUC relate
directly to the way we apply the model online, namely filtering of the top 100
results, transfer learning with the point-wise ensemble is our choice
(outperforming GBDT by 3.6% on negative PR-AUC).

###  4.4. Case Studies

To illustrate the power of batch negative co-training, we look at example
queries and items and how they score under the models trained with and without
batch negatives (vanilla model), both under the Siamese pairwise architecture.

We use the query “cellphone” and the top 3 results returned by our commercial
search engine, two of which are Huawei phones and the third one is an Apple
iPhone X. We then change the query to “paper cup”, which is something
completely unrelated and chosen randomly.

Under the session-pairs only model (first 2 columns in Table 3), the three
smart phones do receive higher scores under “cellphone” than under “paper
cup”, but the score ranges are very similar. In contrast, the batch-negative
enhanced model produces 10 times greater score separation under the two
queries, without significantly dilating the cluster score range itself. Thus
the latter model exhibits much greater confidence in scoring and penalizes
comically unrelated (query, item) pairs, in accordance with intuition.

##  5\. Conclusion

We present a complete description of a 2-stage neural net based relevance
scoring system. On an independent test set our model, using only a single text
feature (item title), outperforms traditional IR systems based on 60 numerical
features by a wide margin, and allows generalization into the long tail where
only text features are available. By exploiting several generic architectural
enhancements, most notably batch negative co-training, we obtain significant
qualitative improvement over the base network. Though we only study feed-
forward architectures here for the sake of serving scalability, these
techniques can also apply to more computationally intensive ones such as
CNN/RNN/attention.

## References

  * [1] Q. Ai, Y. Zhang, K. Bi, X. Chen, and W. B. Croft.  Learning a hierarchical embedding model for personalized product search.  In Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval, pages 645–654. ACM, 2017\.
  * [2] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah.  Signature verification using a” siamese” time delay neural network.  In Advances in neural information processing systems, pages 737–744, 1994.
  * [3] C. Burges, T. Shaked, E. Renshaw, A. Lazier, M. Deeds, N. Hamilton, and G. Hullender.  Learning to rank using gradient descent.  In Proceedings of the 22Nd International Conference on Machine Learning, ICML ’05, pages 89–96, New York, NY, USA, 2005. ACM.
  * [4] C. J. C. Burges.  From RankNet to LambdaRank to LambdaMART: An overview.  Technical report, Microsoft Research, 2010.
  * [5] T. Chen and C. Guestrin.  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.
  * [6] S. Chopra, R. Hadsell, Y. LeCun, et al.  Learning a similarity metric discriminatively, with application to face verification.  In CVPR (1), pages 539–546, 2005.
  * [7] J. Davis and M. Goadrich.  The relationship between precision-recall and roc curves.  In Proceedings of the 23rd international conference on Machine learning, pages 233–240. ACM, 2006.
  * [8] J. Guo, Y. Fan, Q. Ai, and W. B. Croft.  A deep relevance matching model for ad-hoc retrieval.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, pages 55–64, 2016.
  * [9] P.-S. Huang, X. He, J. Gao, L. Deng, A. Acero, and L. Heck.  Learning deep structured semantic models for web search using clickthrough data.  In Proceedings of the 22nd ACM international conference on Information & Knowledge Management, pages 2333–2338, 2013.
  * [10] T. Joachims.  Optimizing search engines using clickthrough data.  In Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’02, pages 133–142, 2002.
  * [11] T. Kenter, A. Borisov, C. Van Gysel, M. Dehghani, M. de Rijke, and B. Mitra.  Neural networks for information retrieval.  In Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR ’17, pages 1403–1406, 2017.
  * [12] T.-Y. Liu.  Learning to rank for information retrieval.  Found. Trends Inf. Retr., 3(3):225–331, Mar. 2009.
  * [13] T. Mikolov, K. Chen, G. Corrado, and J. Dean.  Efficient estimation of word representations in vector space, 2013.
  * [14] B. Mitra and N. Craswell.  An introduction to neural information retrieval.  Foundations and Trends® in Information Retrieval, 13(1):1–126, 2018.
  * [15] B. Mitra, F. Diaz, and N. Craswell.  Learning to match using local and distributed representations of text for web search.  In Proceedings of the 26th International Conference on World Wide Web, WWW ’17, pages 1291–1299, 2017.
  * [16] L. Pang, Y. Lan, J. Guo, J. Xu, J. Xu, and X. Cheng.  Deeprank: A new deep architecture for relevance ranking in information retrieval.  In Proceedings of the 2017 ACM on Conference on Information and Knowledge Management, CIKM ’17, pages 257–266, 2017.
  * [17] Y. Shen, X. He, J. Gao, L. Deng, and G. Mesnil.  A latent semantic model with convolutional-pooling structure for information retrieval.  In Proceedings of the 23rd ACM International Conference on Conference on Information and Knowledge Management, CIKM ’14, pages 101–110, 2014.
  * [18] C. Van Gysel, M. de Rijke, and E. Kanoulas.  Learning latent vector spaces for product search.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, pages 165–174. ACM, 2016.
  * [19] J. Xu and H. Li.  Adarank: A boosting algorithm for information retrieval.  In Proceedings of the 30th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR ’07, pages 391–398, 2007.
  * [20] D. Yin, Y. Hu, J. Tang, T. Daly, M. Zhou, H. Ouyang, J. Chen, C. Kang, H. Deng, C. Nobata, et al.  Ranking relevance in yahoo search.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 323–332, 2016.
  * [21] W. Yin, H. Schütze, B. Xiang, and B. Zhou.  Abcnn: Attention-based convolutional neural network for modeling sentence pairs.  Transactions of the Association for Computational Linguistics, 4:259–272, 2016.
  * [22] J. Yu, S. Mohan, D. P. Putthividhya, and W.-K. Wong.  Latent dirichlet allocation based diversified retrieval for e-commerce search.  In Proceedings of the 7th ACM international conference on Web search and data mining, pages 463–472. ACM, 2014.
  * [23] Z. Zheng, K. Chen, G. Sun, and H. Zha.  A regression framework for learning ranking functions using relative relevance judgments.  In Proceedings of the 30th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR ’07, pages 287–294, 2007.