# Learning Deep Nets for Gravitational Dynamics with Unknown Disturbance
through Physical Knowledge Distillation: Initial Feasibility Study

Hongbin Lin1, Qian Gao2, Xiangyu Chu1, Qi Dou3, Anton Deguet4, Peter
Kazanzides4, and K. W. Samuel Au1 Manuscript received: October, 13, 2020;
Revised January, 13, 2021; Accepted February, 8, 2021.This paper was
recommended for publication by Editor Pietro Valdastri upon evaluation of the
Associate Editor and Reviewers’ comments. This work was supported in part by
the CUHK Chow Yuk Ho Technology Centre of Innovative Medicine, in part by the
Multi-Scale Medical Robotics Centre, InnoHk, 8312051, RGC T42-409/18-R, in
part by SHIAE (BME-p1-17), and in part by the Natural Science Foundation of
China under Grant U1613202. (Corresponding author: K. W. Samuel Au)1Hongbin
Lin, Xiangyu Chu and K. W. Samuel Au are with Department of Mechanical and
Automation Engineering, The Chinese University of Hong Kong, Hong Kong.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Gao is with Shenzhen Institute of Artificial Intelligence and Robotics for
Society, School of Science and Engineering, The Chinese University of Hong
Kong, Shenzhen<EMAIL_ADDRESS>Dou is with Department of
Computer Science and Engineering, The Chinese University of Hong Kong, Hong
Kong<EMAIL_ADDRESS>Deguet and Peter Kazanzides are with
Department of Computer Science, Johns Hopkins University, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>Object Identifier (DOI): see top of
this page.

###### Abstract

Learning high-performance deep neural networks for dynamic modeling of high
Degree-Of-Freedom (DOF) robots remains challenging due to the sampling
complexity. Typical unknown system disturbance caused by unmodeled dynamics
(such as internal compliance, cables) further exacerbates the problem. In this
paper, a novel framework characterized by both high data efficiency and
disturbance-adapting capability is proposed to address the problem of modeling
gravitational dynamics using deep nets in feedforward gravity compensation
control for high-DOF master manipulators with unknown disturbance. In
particular, Feedforward Deep Neural Networks (FDNNs) are learned from both
prior knowledge of an existing analytical model and observation of the robot
system by Knowledge Distillation (KD). Through extensive experiments in high-
DOF master manipulators with significant disturbance, we show that our method
surpasses a standard Learning-from-Scratch (LfS) approach in terms of data
efficiency and disturbance adaptation. Our initial feasibility study has
demonstrated the potential of outperforming the analytical teacher model as
the training data increases.

###### Index Terms:

Medical Robots and Systems, Dynamics, Model Learning for Control.

##  I INTRODUCTION

Model-based control for modern robotics, such as compliance control, impedance
control and hybrid force control, endows robots with the ability to deal with
interaction with the environment. The performance of such control strategies
highly depends on the accuracy of dynamic models of robots, which are capable
of reasoning about joint acceleration and joint torques (forward and inverse
dynamics).

Figure 1: Comparing the number of joints, Degrees of Freedom (DOF) and
sampling complexity among traditional rigid-link robots, serpentine robots and
soft robots/continuum manipulators. Robots in the figure from left to right
are a UR5 robot arm, a CMU snake robot[1] and a Festo manipulator[2],
representing rigid-link robots, serpentine robots and soft robots/continuum
manipulators, respectively.

Conventionally, dynamic models can be derived by physics law analytically[3].
Such physics-based analytical models, parameterized by kinematic and dynamic
properties of robot system (such as mass, length of links, the center of mass,
etc.), are able to describe a global mapping over the entire state space of
the robot system. The parameters of the analytical models can be learned by
offline approaches (e.g., system identification) and online approaches (e.g.,
adaptive control)[4]. However, the analytical solutions cannot easily adapt to
high system disturbance caused by unmodeled dynamics (such as internal
compliance, cables, and so on) due to their fixed basis function[5]. In such a
case, extensive engineering efforts are required to further model the
disturbance for specific systems, which is an arduous process for a robot
system with complex disturbance [6].

To overcome the limitation of the physics-based approaches, researchers start
to turn their attention to model-free approaches leveraging machine learning
techniques as alternative solutions. Particularly, Linear Regression[7][8],
Gaussian Mixture Regression[9], Gaussian Process Regression[10], Support
Vector Machine[11], Feedforward Neural Network[12, 13, 14] and Recurrent
Neural Network[15][16] were applied to approximate the system dynamics of
robots based on measured data. These methods, known as Learning from Scratch
(LfS), did not require any prior physics knowledge. Furthermore, LfS was
potentially capable of capturing the unknown nonlinearity of system
disturbance due to its higher model capacity compared to the physics-based
approach. However, achieving high performance over the entire state space of
the robot system using the LfS approach for high Degree-Of-Freedom (DOF)
robots still remained challenging[17]. Sampling data using the LfS approach
was required to cover the majority of the state space with adequate density.
This requirement was hard to reach for high-DOF systems, since the sampling
complexity, $O(N^{n})$, exponentially increased with the dimensionality of the
state space, where $n$ is the number of DOF and $N$ is sampling points for
each DOF. Fig. 1 shows the sampling complexity for rigid-link robots,
serpentine robots[1] and soft robots[18]/continuum manipulators[2],
illustrating the challenge of sampling complexity for high-DOF robots. In
practice, the sampling problem prevents practitioners from applying the LfS
approach to high-DOF robots considering the expense of data generation (e.g.,
mechanical maintenance, long operation time).

To enjoy merits from the physics-based and the LfS approaches, researchers
attempted to explore hybrid approaches by combining these methods. One thread
was error-learning frameworks where the error of the physics-based model was
learned by Locally Weighted Projection Regression[19][20], Gaussian Process
Regression [17, 21] and Feedforward Neural Network[22]. Although the error-
learning frameworks enhanced the model capacity compared to the pure physics-
based approach for adapting to a high-disturbance system, they suffered from
low data efficiency similar to that of the standard LfS approach since the
errors over the entire state space of a system were required to be learned
from scratch without any prior. Therefore, they did not solve the problem of
sampling complexity fundamentally due to their LfS learning fashion.

Another thread of the hybrid approach was learning system dynamics with
estimators encoded by prior knowledge of physics inference. Dynamic knowledge
(e.g., Newton-Euler inference[23] and Euler-Lagrange inference[24]) was
encoded into the topology of neural networks. These approaches showed higher
data efficiency compared to the standard LfS approach. However, the ability to
adapt to the unknown disturbances remained poor since the model capacity of
physics-encoded approximators was drastically bounded by the pre-encoded
physics prior. As the result, further engineering efforts of modeling
disturbance were required to minimize the model bias caused by unmodeled
dynamics, which can hardly alleviate the engineering burden compared to those
of the physics-based approach.

To disentangle the dynamic modeling problem for high-DOF robot systems with
unknown disturbance, we propose a data-efficient, disturbance-adapting hybrid
approach using a knowledge fusion technique, Knowledge Distillation (KD)[25].
KD was originally proposed to solve the problem of model compression[26],
where the knowledge of a large cumbersome model or ensemble of models, called
Teacher Model (TM), is transferred to a small efficient model, called Student
Model. For KD, Student Model learned from TM by optimizing a weighted average
of two objective functions, the loss function of Student Model with TM and
that with training data. It is widely used in AI community due to its high
efficiency in inheriting the knowledge from the large cumbersome model (e.g.
the generalization) as well as learning from new data. In this paper, the
student Feedforward Deep Neural Network (FDNN) is learned from the knowledge
of an existing physics-based TM, called Physical Teacher Model (PTM), and the
measured data for learning one of the classical components of dynamics,
gravitational dynamics, in high-DOF rigid-link robots. The learned model is
further applied to feedforward Gravity Compensation Control (GCC) for
evaluating its performance.

To our best knowledge, we are the first to demonstrate the efficacy of KD in
the field of dynamics modeling for 6-DOF rigid-body manipulators. Moreover, it
can be potentially applied to the robot systems with higher DOF (such as soft
robots/continuum manipulators), which is our future research direction. We
have developed a fully-automatic open-source software package for the
reproduction of our method and experiment[27]. We tested our method on 3 robot
manipulators to verify the algorithm’s reliability or adaptability to handle
the variations in robots.

Overall, our contributions are:

  1. 1.

Proposal of a novel knowledge-fusion learning framework for modeling
gravitational dynamics of a robot system using FDNN, which pre-processes
input-output signals of FDNN by both trigonometric representation and
normalization, and learns from PTM and observed data using KD. Our method is
applicable for GCC in high-DOF robot systems with unknown highly nonlinear
configuration-dependent and direction-dependent disturbance given sparse
training data.

  2. 2.

Validation of the effectiveness of our method in master manipulators with
significant configuration-dependent and direction-dependent disturbance. Our
method exhibits higher data efficiency and disturbance-adapting capability
compared to the standard LfS method as well as the potential to outperform the
PTM with the increase of training data.

##  II Background

Dynamics for rigid-link robots can be formulated as

$\boldsymbol{\tau}=\boldsymbol{M}(\boldsymbol{q})\boldsymbol{\ddot{q}}+\boldsymbol{C}(\boldsymbol{q},\boldsymbol{\dot{q}})\boldsymbol{\dot{q}}+\boldsymbol{f_{f}}(\boldsymbol{\dot{q}})+\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{J}(\boldsymbol{q})^{T}\boldsymbol{h_{e}}+\boldsymbol{\epsilon},$
(1)

where $\boldsymbol{q}$,$\boldsymbol{\dot{q}}$,$\boldsymbol{\ddot{q}}$ denote
joint position, velocity and acceleration, respectively; $\boldsymbol{M}$ and
$\boldsymbol{C}$ are a inertia matrix and a matrix of the centrifugal and
Coriolis force, respectively; $\boldsymbol{f_{f}}$ is a force vector
collecting the Coulomb and viscous friction, $\boldsymbol{g}$ represents
gravity, $\boldsymbol{\epsilon}$ represents the disturbance due to unmodeled
dynamics, $\boldsymbol{J}$ is a Jacobian matrix and $\boldsymbol{h_{e}}$ is a
vector representing forces, moments exerted by the environment on the end-
effector.

To compensate the dynamic effect of gravity in a robot system, a feedforward
model-based GCC can be applied given the dynamics model of gravity. The
required compensated torques $\boldsymbol{\tau_{c}}$ in GCC is given as

$\boldsymbol{\tau_{c}}=\boldsymbol{g}(\boldsymbol{q}).$ (2)

To estimate the gravity in (2), one can utilize a physics-based, analytical
model derived by either the Euler-Lagrange approach or Newton-Euler
approach[3]. Parameters for the physics-based model can be obtained by
standard Least Square Estimation [28] based on the linearity of the analytical
model[3], if $\boldsymbol{\epsilon}$ is negligible.

However, when the disturbance in the robot system is significant, the
performance using the physics-based model is relatively poor due to its weak
disturbance-adapting capability. In our previous study, significant joint
drifts of a master manipulator, indicating poor control performance, were
observed in a Drift Test when evaluating the naive control in (2)[6]. The main
reason for the degenerated performance is that disturbance torques from the
unmodeled dynamics (such as the disturbance from electric cables in [6]) is
considerably large compared to the gravity.

To improve the naive control in (2), compensated torques for disturbance were
taken into account in GCC. In this paper, we study the configuration-dependent
direction-dependent disturbance
$\boldsymbol{\epsilon(\boldsymbol{q},\boldsymbol{\Delta q})}$ when the robot
system is static (i.e. both the velocity and acceleration of all joints are
zeros). We decompose the disturbance into three parts: configuration-dependent
disturbance $\boldsymbol{\epsilon_{c}}$, direction-dependent disturbance
w.r.t. positive joint direction $\boldsymbol{\epsilon_{d}^{+}}$ and that
w.r.t. negative joint direction $\boldsymbol{\epsilon_{d}^{-}}$.
$\boldsymbol{\epsilon}$ can be formulated as

$\boldsymbol{\epsilon}=\boldsymbol{\epsilon_{c}}(\boldsymbol{q})+\boldsymbol{\epsilon_{d}^{+}}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\epsilon_{d}^{-}}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q})),$
(3)

where $\odot$ denotes Hadamard product, joint direction $\Delta\boldsymbol{q}$
is the difference between current joint position and joint position in the
last iteration of the control loop, $u(\cdot)$ is the step function to
determine the directions. Fig. 2 shows the illustrate three types of
disturbances for each joint of a serial manipulator. In the figure, the
configuration-dependent disturbance is equal to the force from highly
nonlinear springs, while the direction-dependent disturbance is equal to the
force from highly nonlinear dampers related to joint directions instead of
joint velocity. (3) can be simplified as

$\displaystyle\boldsymbol{\epsilon}=$
$\displaystyle\boldsymbol{\epsilon_{c}}(\boldsymbol{q})(u(\Delta\boldsymbol{q})+1-u(\Delta\boldsymbol{q}))+$
(4) $\displaystyle\boldsymbol{\epsilon_{d}^{+}}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\epsilon_{d}^{-}}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q}))$
$\displaystyle=$
$\displaystyle(\boldsymbol{\epsilon_{c}}(\boldsymbol{q})+\boldsymbol{\epsilon_{d}^{+}}(\boldsymbol{q}))\odot
u(\Delta\boldsymbol{q})+$
$\displaystyle(\boldsymbol{\epsilon_{c}}(\boldsymbol{q})+\boldsymbol{\epsilon_{d}^{-}}(\boldsymbol{q}))\odot(1-u(\Delta\boldsymbol{q}))$
$\displaystyle=$ $\displaystyle\boldsymbol{\epsilon^{+}}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\epsilon^{-}}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q})),$

where
$\boldsymbol{\epsilon^{+}}=\boldsymbol{\epsilon_{c}}+\boldsymbol{\epsilon_{d}^{+}}$
and
$\boldsymbol{\epsilon^{-}}=\boldsymbol{\epsilon_{c}}+\boldsymbol{\epsilon_{d}^{-}}$
are disturbances w.r.t. positive direction $\boldsymbol{\epsilon^{+}}$ and
disturbances w.r.t. negative direction $\boldsymbol{\epsilon^{-}}$,
respectively. Additional torques to compensate disturbance in (4) are added to
(2) for improving the naive control as

$\displaystyle\boldsymbol{\tau_{c}}=$
$\displaystyle\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{\epsilon}(\boldsymbol{q},\Delta\boldsymbol{q})$
(5) $\displaystyle=$
$\displaystyle\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{\epsilon}^{+}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\epsilon}^{-}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q}))$

Figure 2: Illustration of configuration-dependent and directional-dependent
disturbance for each joint of a serial manipulator.

To achieve the control in (5), researchers proposed to estimate gravity using
a physics-based model and to approximate disturbance using either linear
models[29] or polynomial models[30, 31]. In our previous work, we showed that
polynomial models were able to provide a higher model capacity to capture the
nonlinearities of disturbance compared to linear models[31]. However, [29, 30,
31] were limited by the assumption that the configuration-dependent
disturbance for each joint was only related to its own joint position. In this
paper, we solve the limitation by using FDNN for a general configuration-
dependent disturbance which is related to all joint positions.

The rest of the paper is organized as follows. We first introduce our method
to model the gravitational dynamics with disturbance using FDNN in Section
III, followed by learning methods in Section IV. The controller for GCC is
designed in Section V. The efficiency of our method is validated by extensive
experiments of simulation and real hardware in Section VI. Discussion and
conclusion are given in Section VII and VIII, respectively.

##  III Modeling Gravitational Dynamics with Disturbance using FDNNs

The compensated torque in (5) can be rewritten as

$\displaystyle\boldsymbol{\tau_{c}}=$
$\displaystyle\boldsymbol{g}(\boldsymbol{q})(1-u(\Delta\boldsymbol{q})+u(\Delta\boldsymbol{q}))+$
(6) $\displaystyle\boldsymbol{\epsilon}^{+}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\epsilon}^{-}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q}))$
$\displaystyle=$
$\displaystyle(\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{\epsilon}^{+}(\boldsymbol{q}))\odot
u(\Delta\boldsymbol{q})+$
$\displaystyle(\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{\epsilon}^{-}(\boldsymbol{q}))\odot(1-u(\Delta\boldsymbol{q}))$
$\displaystyle=$ $\displaystyle\boldsymbol{\tau_{c}}^{+}(\boldsymbol{q})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{\tau_{c}}^{-}(\boldsymbol{q})\odot(1-u(\Delta\boldsymbol{q})),$

where $\boldsymbol{\tau_{c}}^{+}=\boldsymbol{g}+\boldsymbol{\epsilon^{+}}$ and
$\boldsymbol{\tau_{c}}^{-}=\boldsymbol{g}+\boldsymbol{\epsilon^{-}}$ are
compensated torques w.r.t. positive and negative directions, respectively.
$\boldsymbol{\tau_{c}^{+}}$ and $\boldsymbol{\tau_{c}^{-}}$ are approximated
by two sets of FDNNs as

$\boldsymbol{\hat{\tau}_{c}}=\mathcal{G}(\boldsymbol{q};\boldsymbol{W^{+}})\odot
u(\Delta\boldsymbol{q})+\mathcal{G}(\boldsymbol{q};\boldsymbol{W^{-}})\odot(1-u(\Delta\boldsymbol{q})),$
(7)

where $\hat{\cdot}$ refers to an approximation, $\mathcal{G}$ denotes
continuous functional mapping of an FDNN, $\boldsymbol{W^{+}}$ and
$\boldsymbol{W^{-}}$ are matrices of weight parameters for FDNNs w.r.t.
$\boldsymbol{\hat{\tau_{c}}^{+}}$ and $\boldsymbol{\hat{\tau_{c}}^{-}}$,
respectively.

To achieve stable and efficient control using FDNNs, input-output signals of
FDNNs require data pre-processing through a sequence of transformations. We
will discuss the transformations in the following sub-sections.

###  III-A Trigonometric Representation for Joint Position

The $\boldsymbol{q}$ in (7) is traditionally represented using a minimum
representation as
$\boldsymbol{q}=\begin{bmatrix}q_{1}&q_{2}&\dots&q_{n}\end{bmatrix}^{T}\in\mathbb{R}^{n},$
where the position of the $i^{\mathrm{th}}$ joint is
$q_{i}=\theta_{i}\in[0,2\pi)$ for a revolute joint and
$q_{i}=d_{i}\in\mathbb{R}$ for a prismatic joint without considering joint
limits. However, the minimal representation cannot guarantee the continuity of
the functional mapping $\mathcal{G}$ in joint space due to the singularities
of the representation(i.e the neighborhood of $0$ for $q_{i}$ in the case of a
revolute joint)[32], which might lead to unstable performance (Prove can be
found in [33]). To overcome the limitation, we use a non-minimum
representation, trigonometric representation, where $\boldsymbol{q}$ is
transformed as

$\displaystyle\boldsymbol{z}=$
$\displaystyle\Phi(\boldsymbol{q})=[\Phi(q_{1}),\Phi(q_{2}),\dots,\Phi(q_{n})]^{T}\in\mathbb{R}^{n+k},$
(8)

where $k$ is the number of revolute joints for a robot; The function for the
trigonometric transformation $\Phi$ is defined as

$\Phi(x)=\begin{cases}[\sin(x),\cos(x)]^{T}&\text{for a revolute joint}\\\
x&\text{for a prismatic joint}\end{cases};$ (9)

Compared to the minimal representation, the trigonometric representation can
guarantee the continuity of $\mathcal{G}$ in joint space, which was also
proved in [33].

###  III-B Input-Output Normalization

Although an FDNN can generate a mapping between raw input and output with
controlled accuracy, in practice the performance of an FDNN can be improved by
normalization to reduce the risk of getting stuck in poor local minima[34].
Furthermore, normalization contributes to preventing overfitting for joints
with larger outputs in the system.

Zero-score Normalization is utilized to scale signals into a Gaussian
distribution with zero mean and unit standard deviation. The zero-score
function is defined as

$\boldsymbol{{}^{n}a}=f_{n}(\boldsymbol{a};\boldsymbol{\mu},\boldsymbol{\sigma})=(\boldsymbol{a}-\boldsymbol{\mu})\oslash\boldsymbol{\sigma},$
(10)

where $\oslash$ is Hadamard division, $\boldsymbol{\mu}$ and
$\boldsymbol{\sigma}$ denote the mean and the standard deviation,
respectively, $\boldsymbol{{}^{n}}(\cdot)$ denotes a normalized value and
$\boldsymbol{a}$ is an input vector. Similarly, the inverse of the zero-score
function is given as

$\boldsymbol{a}=f_{n}^{-1}(\boldsymbol{{}^{n}a};\boldsymbol{\mu},\boldsymbol{\sigma})=\boldsymbol{{}^{n}a}+\boldsymbol{\mu}\odot\boldsymbol{\sigma}.$
(11)

The signals in (7) and (8) can be normalized as

$\boldsymbol{{}^{n}z}=f_{n}(\boldsymbol{z};\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}}),$
(12)
$\boldsymbol{{}^{n}\hat{\tau_{c}}}=f_{n}(\boldsymbol{\hat{\tau_{c}}};\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}}),$
(13)

where $\boldsymbol{\mu_{in}}$ and $\boldsymbol{\sigma_{in}}$ are parameter
vectors of mean and standard deviation for the input signals $\boldsymbol{z}$,
while $\boldsymbol{\mu_{out}}$ and $\boldsymbol{\sigma_{out}}$ are parameter
vectors of mean and standard deviation for the output signals
$\boldsymbol{\tau_{c}}$. Instead of mapping raw input-output signals in (7),
we map the normalized input-output signals by two sets of FDNNs as

$\displaystyle\boldsymbol{{}^{n}\hat{\tau_{c}}}=$
$\displaystyle\boldsymbol{{}^{n}\hat{\tau}_{c}}^{+}(\boldsymbol{{}^{n}z})\odot
u(\Delta\boldsymbol{q})+\boldsymbol{{}^{n}\hat{\tau}_{c}}^{-}(\boldsymbol{{}^{n}z})\odot(1-u(\Delta\boldsymbol{q}))$
(14) $\displaystyle=$
$\displaystyle\mathcal{G}(\boldsymbol{{}^{n}z};\boldsymbol{{}^{n}W_{z}^{+}})\odot
u(\Delta\boldsymbol{q})+$
$\displaystyle\mathcal{G}(\boldsymbol{{}^{n}z};\boldsymbol{{}^{n}W_{z}^{-}})\odot(1-u(\Delta\boldsymbol{q})),$

where $\boldsymbol{{}^{n}W_{z}^{+}}$ and $\boldsymbol{{}^{n}W_{z}^{-}}$ are
matrices of parameters of the FDNNs for the normalized mapping corresponding
to positive and negative directions, respectively. Combining (12), (14) and
the inverse of (13), the estimated compensated torques in (7) can be rewritten
as

$\displaystyle\boldsymbol{\hat{\tau_{c}}}=f_{n}^{-1}($
$\displaystyle\mathcal{G}(f_{n}(\boldsymbol{z};\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}});\boldsymbol{{}^{n}W_{z}^{+}})\odot
u(\Delta\boldsymbol{q})+$ (15)
$\displaystyle\mathcal{G}(f_{n}(\boldsymbol{z};\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}});\boldsymbol{{}^{n}W_{z}^{-}})\odot(1-u(\Delta\boldsymbol{q}));$
$\displaystyle\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}}).$

$\boldsymbol{{}^{n}W_{z}^{+}}$, $\boldsymbol{{}^{n}W_{z}^{-}}$,
$\boldsymbol{\mu_{in}}$, $\boldsymbol{\sigma_{in}}$, $\boldsymbol{\mu_{out}}$
are $\boldsymbol{\sigma_{out}}$ are parameters of the predicted model for GCC.
We will discuss how to learn the parameters in Section IV.

##  IV Learning FDNNs for Dynamics Modeling

In this section, we will elaborate on two forms of learning schemes, Learning
from Scratch (LfS) and learning with Physical Knowledge Distillation (PKD).

###  IV-A Learning from Scratch (LfS)

To learn gravity and the disturbance of a robot system, observed data of both
joint positions and joint torques should be sampled from the robot system. The
data set for the observed data is defined as

$\displaystyle\mathcal{D}_{system}=\\{\boldsymbol{x_{1}^{i}}\leftarrow\Phi(\boldsymbol{{q}^{i}}),\boldsymbol{x_{2}^{i}}\leftarrow\boldsymbol{{\Delta
q}^{i}},\boldsymbol{y^{i}}\leftarrow\boldsymbol{{\tau_{c}}^{i}}\\}_{i=1}^{T^{s}},$
(16)

where $\boldsymbol{{q}^{i}}$, $\boldsymbol{\Delta{q}^{i}}$ and
$\boldsymbol{{\tau_{c}}^{i}}$ represent $\boldsymbol{q}$, $\boldsymbol{\Delta
q}$ and $\boldsymbol{{\tau_{c}}}$ at the $i^{\mathrm{th}}$ configuration,
respectively, $\boldsymbol{x_{1}^{i}}$ and $\boldsymbol{x_{2}^{i}}$ are the
$i^{\mathrm{th}}$ pair of sampled input vectors while $\boldsymbol{y^{i}}$ is
the $i^{\mathrm{th}}$ sampled output vector, $T^{s}$ is the total amount of
sampled data in $\mathcal{D}_{system}$.

The parameters of the input-output normalization are identified by properties
of the distribution of
$(\boldsymbol{x_{1}},\boldsymbol{y})\in\mathcal{D}_{system}$, where
$\boldsymbol{\mu_{in}}$ and $\boldsymbol{\sigma_{in}}$ are obtained by the
mean and standard deviation of $\boldsymbol{x_{1}}$ in $\mathcal{D}_{system}$,
while $\boldsymbol{\mu_{out}}$ and $\boldsymbol{\sigma_{out}}$ are obtained by
the mean and standard deviation of $\boldsymbol{y}$ in $\mathcal{D}_{system}$.
Here, we do not normalize the $\boldsymbol{\Delta q}$ signal since it is not
included in the input or output of FDNN in (15). Therefore, the normalized
data set of $\mathcal{D}_{system}$ can be obtained by

$\displaystyle{}^{n}\mathcal{D}_{system}=\\{\boldsymbol{x_{1}}$
$\displaystyle\leftarrow(f_{n}(\boldsymbol{x_{1}};\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}}),~{}\boldsymbol{x_{2}}\leftarrow\boldsymbol{x_{2}},$
(17) $\displaystyle\boldsymbol{y}$ $\displaystyle\leftarrow
f_{n}(\boldsymbol{y};\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}})\\}_{(\boldsymbol{x_{1}},\boldsymbol{x_{2}},\boldsymbol{y})\in\mathcal{D}_{system}}.$

$\boldsymbol{{}^{n}W_{z}^{+}}$ and $\boldsymbol{{}^{n}W_{z}^{-}}$, are learned
by minimizing an loss function of the LfS, which is defined as

$\mathcal{L}^{s}=\frac{1}{T^{s}}\sum||\boldsymbol{{}^{n}\hat{\tau_{c}}}(\boldsymbol{x};\boldsymbol{{}^{n}W_{z}^{+}},\boldsymbol{{}^{n}W_{z}^{-}})-\boldsymbol{y}||_{(\boldsymbol{x},\boldsymbol{y})\in{}^{n}\mathcal{D}_{system}},$
(18)

where $\boldsymbol{x}$ is the combination of $\boldsymbol{x_{1}}$ and
$\boldsymbol{x_{2}}$. The loss function with regularization is given as

$\mathcal{L}^{s}_{r}=\mathcal{L}^{s}+\mathcal{R}(\boldsymbol{{}^{n}W_{z}^{+}},\boldsymbol{{}^{n}W_{z}^{-}}),$
(19)

where $\mathcal{R}$ is a regularization function. In our experiment, we used
$L^{2}$ Norm as the regularization function.

###  IV-B Learning with Physical Knowledge Distillation (PKD)

Since the LfS approach is limited by the problem of sampling complexity for
high-DOF robots, we propose a high data-efficient method using KD to learn
FDNN effectively. KD was originally proposed to solve the model compression
problem for classification[25]. In particular, the objective function was
designed as the weighted average of the cross entropy of Student Model with
the soft target of a TM (obtained by scaling the output of the TM using
Temperature Scaling) and that with labeled correct data. In this paper,
Temperature Scaling is not included in our method since it showed no
improvement for regression as discussed in[35].

Compared with the nerual-network-based TM in the origin version, we use the
state-of-art physics-based model in[31] as TM, called Physical Teacher Model
(PTM). (6) was modelled by the PTM analytically, denoted as a function
$\boldsymbol{{}^{p}{\tau_{c}}}(\cdot)$ that maps from $\boldsymbol{q}$ and
$\boldsymbol{\Delta q}$ to compensated torques.

To distill the knowledge of a PTM, we randomly sampled $T^{p}$ data points
from the PTM within its state space. For the $i^{\mathrm{th}}$ data point,
given the inputs $\boldsymbol{{q}^{i}}$ and $\boldsymbol{{\Delta q}^{i}}$, the
function of the PTM can predict the output torque for compensation, denoted as
$\boldsymbol{{{}^{p}\tau_{c}}(\boldsymbol{{q}^{i}},\boldsymbol{{\Delta
q}^{i}})}$. We define a data set $\mathcal{D}_{PTM}$ that collects all the
sampled data ($\boldsymbol{{q}}$, $\boldsymbol{{\Delta q}}$,
$\boldsymbol{{{}^{p}\tau_{c}}(\boldsymbol{{q}},\boldsymbol{{\Delta q}})}$)
from the PTM using the representation in (16).

The parameters of the normalization are identified by the distribution
property of the joint data set
$\mathcal{D}_{joint}=\mathcal{D}_{system}\cup\mathcal{D}_{PTM}$. In
particular, $\boldsymbol{\mu_{in}}$ and $\boldsymbol{\sigma_{in}}$ are
obtained by the mean and standard deviation of $\boldsymbol{x_{1}}$ in
$\mathcal{D}_{joint}$, while $\boldsymbol{\mu_{out}}$ and
$\boldsymbol{\sigma_{out}}$ are obtained by the mean and standard deviation of
$\boldsymbol{y}$ in $\mathcal{D}_{joint}$. Then, we normalize
$\mathcal{D}_{PTM}$ using (17) as ${}^{n}\mathcal{D}_{PTM}$.

The loss function of the PTM can be represented as

$\mathcal{L}^{p}=\frac{1}{T^{p}}\sum||\boldsymbol{{}^{n}\hat{\tau_{c}}}(\boldsymbol{x};\boldsymbol{{}^{n}W_{z}^{+}},\boldsymbol{{}^{n}W_{z}^{-}})-\boldsymbol{y}||_{(\boldsymbol{x},\boldsymbol{y})\in{}^{n}\mathcal{D}_{PTM}},$
(20)

Combining (18) and (20), the parameters of FDNN can be learned by minimizing
the loss function with regularization for PKD, which can be represented as

$\mathcal{L}_{r}^{{}^{PKD}}=(1-\lambda)\mathcal{L}^{s}+\lambda{}\mathcal{L}^{p}+\mathcal{R}(\boldsymbol{{}^{n}W_{z}^{+}},\boldsymbol{{}^{n}W_{z}^{-}}).$
(21)

where $\lambda\in[0,1]$ is a coefficient value for PKD. A high value of
$\lambda$ is related to high confidence in the PTM.

Our learning method (Algorithm 1) can be summarized as following:

  1. 1.

Sampling PTM: Sample $\boldsymbol{{}^{p}{\tau_{c}}}(\cdot)$ given $T^{p}$
randomized inputs ($\boldsymbol{q}$, $\boldsymbol{\Delta q}$) (line 3-4) and
collect sampled data to $\mathcal{D}_{PTM}$ (line 5).

  2. 2.

Identifying normalized parameters and normalizing datasets: Create joint
dataset (line 7) and calculate the parameters based on the distribution of the
joint dataset (line 8-9). Normalize $\mathcal{D}_{system}$ and
$\mathcal{D}_{PTM}$ based on the normalized parameters (line 10-11).

  3. 3.

Optimizing the objective function for PKD (line 12).

Algorithm 1 Training FDNNs for GCC using PKD

0: $\boldsymbol{{}^{n}W_{z}^{+}}$, $\boldsymbol{{}^{n}W_{z}^{-}}$,
$\boldsymbol{{}^{p}{\tau_{c}}}(\cdot)$, $\mathcal{D}_{system}$, $T^{p}$

0: $\boldsymbol{{}^{n}W_{z}^{+^{*}}}$, $\boldsymbol{{}^{n}W_{z}^{-^{*}}}$,
$\boldsymbol{\mu_{in}}$, $\boldsymbol{\sigma_{in}}$, $\boldsymbol{\mu_{out}}$,
$\boldsymbol{\sigma_{out}}$

1: $\mathcal{D}_{PTM}=\varnothing$

2: for $i=1$ to $T^{p}$ do

3: Randomize $(\boldsymbol{{q}^{i}},\boldsymbol{\Delta{q}^{i}})$ within the
state space

4:
$\boldsymbol{x_{1}^{i}}\leftarrow\Phi(\boldsymbol{{q}^{i}}),\boldsymbol{x_{2}^{i}}\leftarrow\boldsymbol{{\Delta
q}^{i}},\boldsymbol{y^{i}}\leftarrow\boldsymbol{{{}^{p}\tau_{c}}(\boldsymbol{{q}^{i}},\boldsymbol{{\Delta
q}^{i}})}$

5:
$\mathcal{D}_{PTM}\leftarrow\mathcal{D}_{PTM}\cup\\{\boldsymbol{x_{1}^{i}},\boldsymbol{x_{2}^{i}},\boldsymbol{y^{i}}\\}$

6: end for

7: $\mathcal{D}_{joint}\leftarrow\mathcal{D}_{system}\cup\mathcal{D}_{PTM}$

8: $\boldsymbol{\mu_{in}}$,
$\boldsymbol{\sigma_{in}}\leftarrow\boldsymbol{\mu}$, $\boldsymbol{\sigma}$ of
$\boldsymbol{x_{1}}$ in $\mathcal{D}_{joint}$

9: $\boldsymbol{\mu_{out}}$,
$\boldsymbol{\sigma_{out}}\leftarrow\boldsymbol{\mu}$, $\boldsymbol{\sigma}$
of $\boldsymbol{y}$ in $\mathcal{D}_{joint}$

10: ${}^{n}\mathcal{D}_{system}\leftarrow$
Norm$(\mathcal{D}_{system},\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}},\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}}$)

11: ${}^{n}\mathcal{D}_{PTM}\leftarrow$
Norm$(\mathcal{D}_{PTM},\boldsymbol{\mu_{in}},\boldsymbol{\sigma_{in}},\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}}$)

12:
$\boldsymbol{{}^{n}W_{z}^{+^{*}}},\boldsymbol{{}^{n}W_{z}^{-^{*}}}=\underset{\boldsymbol{{}^{n}W_{z}^{+}},\boldsymbol{{}^{n}W_{z}^{-}}}{\arg\min}\mathcal{L}^{PKD}_{r}$

13: return $\boldsymbol{{}^{n}W_{z}^{+^{*}}}$,
$\boldsymbol{{}^{n}W_{z}^{-^{*}}}$, $\boldsymbol{\mu_{in}}$,
$\boldsymbol{\sigma_{in}}$, $\boldsymbol{\mu_{out}}$,
$\boldsymbol{\sigma_{out}}$

##  V Gravity Compensation Control using FDNNs

We design a reliable controller for feedforward GCC using FDNNs. The designed
control is similar to our previous work[31] but in an FDNN version. The
estimated configuration-dependent and direction-dependent compensated torques
with normalization, $\boldsymbol{\hat{\tau}^{c}_{c}}$ and
$\boldsymbol{\hat{\tau}^{d}_{c}}$, are given as

$\boldsymbol{{}^{n}\hat{\tau}^{c}_{c}}=(\boldsymbol{{}^{n}\hat{\tau}_{c}^{+}}+\boldsymbol{{}^{n}\hat{\tau}_{c}^{-}})/2,$
(22)
$\boldsymbol{{}^{n}\hat{\tau}^{d}_{c}}=(\boldsymbol{{}^{n}\hat{\tau}_{c}^{+}}-\boldsymbol{{}^{n}\hat{\tau}_{c}^{-}})/2.$
(23)

The designed policy for GCC using FDNNs is defined as

$\boldsymbol{u}=f_{n}^{-1}(\boldsymbol{{}^{n}\hat{\tau}^{c}_{c}}+\boldsymbol{\xi}(\boldsymbol{\Delta
q})\odot\boldsymbol{{}^{n}\hat{\tau}^{d}_{c}};\boldsymbol{\mu_{out}},\boldsymbol{\sigma_{out}}),$
(24)

where $\boldsymbol{\xi}=[\xi_{1},\xi_{2},\dots,\xi_{n}]^{T}$ is a scale vector
of compensation ratio for the direction-dependent disturbance; the
$i^{\mathrm{th}}$ element of $\boldsymbol{\xi}$ is designed as

$\xi_{i}=\begin{cases}0&|\Delta{q_{i}}|\leq\Delta{q}_{db_{i}}\\\
\frac{|\Delta{q}_{i}|-\Delta{q}_{db_{i}}}{\Delta{q}_{s_{i}}-\Delta{q}_{db_{i}}}\mathrm{sgn}(\Delta{q}_{i})\alpha&\Delta{q}_{db_{i}}<|\Delta{q}_{i}|<\Delta{q}_{s_{i}}\\\
\mathrm{sgn}(\Delta{q}_{i})\alpha&\Delta{q}_{s_{i}}\leq|\Delta{q}_{i}|\end{cases},$
(25)

where $\Delta q_{i}$ is the $i^{\mathrm{th}}$ element of $\boldsymbol{\Delta
q}$, $\alpha$ is a compensation ratio for the directional-dependent
disturbance, $\Delta{q}_{db_{i}}$ is a dead-band value and $\Delta{q}_{s_{i}}$
is a saturated value for the $i^{\mathrm{th}}$ joint. In (25), when $|\Delta
q_{i}|$ is small, $\xi_{i}$ is designed to be zero to reject the noise from
the joint measurement. When $|\Delta q_{i}|$ is large, only a certain ratio
$\alpha$ of the estimated directional-dependent torque is compensated. When
$|\Delta q_{i}|$ is between $\Delta{q}_{db_{i}}$ and $\Delta{q}_{s_{i}}$, a
linear interpolation is ultilized to avoid discontinuity. The design for our
GCC is shown in Fig. 3.

Figure 3: Designed controller for GCC using FDNNs.

##  VI Experiments

We designed experiments to validate the effectiveness of our method based on
our previous work[31]. The evaluation was carried out on test benches
consisting of Master Tool Manipulators (MTMs) of a da Vinci Research Kit
(dVRK)[36]. Fig. 4 shows our test bench. We first elaborate on the experiment
settings, followed by showing the results of the experiments in both
simulation and real hardware.

###  VI-A Experiment Setting

The experiments were conducted on MTMs, serial manipulators with seven
motorized revolute joints and one passive pinching joint at the end. Gravity
for the last two joints can be neglected due to the symmetrical structure of
its link about its joint axis. In addition, the disturbance is trivial for the
last two joints such that the compensated torques for the last two joints were
not considered in the experiments. Thus, the MTM can be simplified as a 6-DOF
serial manipulator with the mass of Link 6 lumped with the last 2 links.

####  VI-A1 Baselines

In the experiments, we compared controllers of GCC derived from 3 types of
models: the physics-based model, FDNNs trained by LfS and PKD. The analytical
model served as a PTM to train the FDNNs of PKD. In practice, both Lfs and PKD
can train all types of approximators for regression. In our experiment, we
used a standard architecture of FDNN: 4 fully connected layers (6,30,30,30,6)
with activation function (ReLU[37]) after each layer except the last. The
number of neurons in the hidden layers of the architecture was tuned according
to the learned performance empirically.

####  VI-A2 Sampling Data

Three data sets were sampled in each experiment for training, validating and
testing. To sample evenly within the joint space for training data, systematic
sampling (i.e., sampling the joint space with fixed intervals) and random
sampling (i.e., sampling the joint space randomly) can be utilized, given the
number of sampling data is large. However, when the number of sampling is
small, the data distribution of the random sampling can hardly achieve the
uniform sparsity over the joint space. Therefore, in the simulation where data
can be infinitely sampled, we used the random sampling. In our experiments of
real robots where the sampling data was expensive and limited, the systematic
sampling was applied for the training data, while the random sampling was
applied for the validating and testing data.

####  VI-A3 Training FDNNs

For the learning process of FDNNs, Adam [38], a gradient-descent optimizer,
was utilized for learning the parameters of FDNNs. Furthermore, to prevent
overfitting, $L_{2}$ Norm was utilized as the regularized function
$\mathcal{R}$ and Early Stopping [39] was applied based on the data set for
validation. When the number of sampling data is large, these two kinds of
sampling methods have similar sampling effects. For the implementation, we
kept shared hyper-parameters of baselines using LfS and PKD the same to get a
fair comparison.

For the method of PKD, $T^{p}$ was set to $30000$ for dense sampling from the
PTM. $\lambda$ decreased linearly during the iterative steps of the
optimization. In such a manner, the guidance of PTM was dominant at the
beginning of the training process, while the guidance of sampled training data
increased during the process of the optimization.

Figure 4: Physical MTM setup.

(a) RRMSE for simulation, low-bias PTM  (b) RRMSE for simulation, high-bias
PTM

Figure 5: (a-b) Average RRMSE of 6 joints w.r.t. the amount of sampling data
from the robot system using low-bias and high-bias PTMs. Lines indicate the
mean of 10 repetitive experiments and shaded areas highlight the ranges within
1 standard deviation of the means. The physics-based model in [31] with
different estimation biases served as PTM. According to the results, PKD
outperforms LfS in terms of smaller variances and means of RRMSE in the cases
of sparse learning data, while it has similar performance compared with LfS
when the learning data is rich.

####  VI-A4 Offline Evaluation

Relative Root Mean Square Error (RRMSE) served as an indicator for evaluating
the accuracy of the predicted models in GCC based on the testing data. We
defined RRMSE as

$\epsilon_{RMS\%}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\frac{(\tau^{i}_{c}-\hat{\tau}^{i}_{c})^{2}}{\tau_{c}^{i^{2}}}}\cdot
100\%,$ (26)

where $\tau^{i}_{c}$ and $\hat{\tau}^{i}_{c}$ are the measured torque and the
estimated torque calculated by the evaluated model for the $i^{\mathrm{th}}$
testing data, respectively; $N$ is the amount of the testing data.

###  VI-B Simulation Experiments

We designed simulation experiments to quantify the data-saving effect, where a
state-of-the-art physics-based model in [31] served as a ground-truth model to
simulate gravity and disturbance of an MTM. The PTMs were obtained by adding
high noises and low noises to parameters of the physics-based model for
comparing cases of “Weak Teachers” (i.e. low-bias PTMs) and “Strong Teachers”
(i.e. high-bias PTMs), respectively. In addition, noises were added to the
measurement of joint positions and torques for the training and validating
data. The simulated experiments explored the results in the cases of $T^{s}$
from $10$ to $5000$ for learning data. For each case, $10$ repetitive
experiments were performed considering the randomness of the measurement noise
and the learning processes.

(a) RRMSE for MTM 1, comparing different baselines (b) RRMSE for MTM 1,
comparing different normalized cases.

Figure 6: Offline Test for a real MTM. (a) RRMSE for 6 joints and the average
over 6 joints. (b) Average RRMSE for all joints using PKD comparing different
normalized cases. In (a), PKD outperformed LfS for all joints, while it
slightly outperformed the PTM in terms of the average RRMSE. In (b), both
input and output normalization improved the performance.

Fig. 5a and 5b show the results of the offline evaluation in simulation. Our
approach outperformed LfS when the learning data was sparse. In the case of
the high-bias PTM, the physics-based model suffered from inaccuracy of
parameter estimation, leading to poor performance. Notably, our method
improved the performance with increasing learning data compared to the
physics-based model, while it achieved similar performance compared to the LfS
approach when the learning data was rich. Our approach significantly saved the
learning data compared to LfS (e.g., our method saved around $99\%$ training
data to reach $20\%$ accuracy for average RRMSE of all joints using the low-
bias PTMs).

###  VI-C Real Robot Experiments

Physical MTMs were controlled by the dVRK software system. To reach the
desired state ($\boldsymbol{q}$, $\boldsymbol{\Delta q}$), we first moved the
MTM to a prepared configuration $\boldsymbol{q}-\boldsymbol{\Delta q}$ and
then moved to the desired position $\boldsymbol{q}$ using a PID controller.
Joint positions, directions, and torques were recorded when the MTM reached a
steady state. To sample the training data, we evenly sampled $4$ points for
each joint for both positive and negative directions with the total amount of
sampling $8192=2\times 4^{6}$. The total amounts of validating data and
testing data were around $6000$ and $300$, respectively. We utilized the
physics-based model in [31] as the PTM.

TABLE I: Mean and standard deviation of joint drifts averaged over 6 joints for multiple MTMs in Online Drift Test | PTM | LfS | PKD
---|---|---|---
MTM 1 |  0.7 $\pm$2.4 |  6.1 $\pm$14.9 |  1.0 $\pm$5.7
MTM 2 |  2.0 $\pm$2.9 |  5.5 $\pm$12.8 |  2.1 $\pm$5.2
MTM 3 |  1.3 $\pm$0.5 |  6.9 $\pm$6.7 |  2.0 $\pm$7.2
TABLE II: Mean and standard deviation of Joint Drifts for MTM1 in Online Drift Test | Joint 1 | Joint 2 | Joint 3 | Joint 4 | Joint 5 | Joint 6
---|---|---|---|---|---|---
PTM |  0.4 $\pm$0.8 |  0.5 $\pm$0.9 |  0.5 $\pm$0.9 |  0.6 $\pm$0.6 |  2.0 $\pm$9.7 |  0.5 $\pm$2.4
LfS |  0.1 $\pm$0.1 |  0.9 $\pm$2.0 |  0.9 $\pm$1.6 |  2.5 $\pm$7.0 |  24.8 $\pm$46.2 |  7.7 $\pm$14.9
PKD |  0.4 $\pm$0.9 |  0.4 $\pm$0.8 |  0.4 $\pm$0.8 |  0.8 $\pm$1.4 |  2.3 $\pm$8.6 |  1.7 $\pm$5.7

We designed an online test, Drift Test, based on our previous work [31] to
evaluate the performance of GCC. We tested the MTM at 400 testing points with
random configurations and directions within joint limits. For each test, we
first moved the MTM to a testing point using the PID controller and
subsequently switched to the evaluated controller of GCC running at $500$ Hz
for around 2 seconds. Both joint drifts and Cartesian drifts of the tip were
recorded.

Fig. 6a shows the Offline Test for MTM 1. The average RRMSE of 6 joints for
our method ($29.4\%$) was lower than that of other methods (PTM: $31.3\%$,
LfS: $37.2\%$). In Fig. 6b, we compared four cases to evaluate the
effectiveness of the input-output normalization. Among the four cases, our
approach with the input-output normalization shows the lowest RRMSE.

Table I shows a statistical result of joint drifts averaged over 6 joints for
3 MTMs in the Online Test. Our approach reduced around 75% the mean of the
averaged drifts for 3 MTMs compared with LfS. The performance of our approach
was slightly inferior to that of the state-of-art PTM in terms of the mean and
standard deviation of the averaged joint drifts. However, their performance
was comparable since their maximum mean of drifts for 3 MTMs were in the same
scale (PTM: $2.0$; PKD: $2.1$) compared to that of LfS ($6.9$). Table II shows
the statistical result for individual joints of MTM1 in the Online Test. LfS
had the largest drifts for all joints except Joint 1. The mean of drifts for
Joint 5 with LfS was 10 times larger than that with PKD and PTM. On the other
hand, our approach achieved similar performance with PTM or even outperformed
PTM in terms of drifts of Joint 2 and Joint 3.

##  VII Discussion

Our experimental results showed that PKD can significantly boost the data
efficiency for FDNN in GCC. Besides GCC, our learning framework can be applied
to other model-based dynamic controllers that require learning gravitational
dynamics, such as PD control with gravity compensation.

While our method has successfully reduced bias of drifts of FDNN in the Online
Drift Test, our variance of drifts remained large compared to that of PTM,
which could be improved by ensemble learning methods. Another potential
improvement is incorporating the uncertainty of PTM in our learning process.

In this work, we provide a novel and efficient method to bridge the gap
between the physics-based approaches and the LfS approaches. Our method was
evaluated on 6-DOF rigid-serial-link manipulators through extensive
experiments, taking a key step toward exploiting our method on robot systems
with higher DOFs.

##  VIII Conclusion

We proposed a novel learning framework using PKD to learn from both physical
prior and observed data from the system for modeling gravitational dynamics
using FDNN, where trigonometric representation and input-output normalization
were applied for improving the performance of GCC. The result of experiments
showed that PKD possessed better data efficiency than the standard LfS
approach in GCC for high-DOF robots. In addition, it achieved competitive
performance using sparse data compared to PTM.

Future directions of this work are extending our method to learning full
dynamics with unknown disturbance as well as exploiting our learning framework
for other types of robots (e.g. soft robots/continuum manipulators).

## References

  * [1] C. Wright, A. Buchan _et al._ , “Design and architecture of the unified modular snake robot,” in _Proc. IEEE Int. Conf. Robot. and Autom._ , 2012, pp. 4347–4354.
  * [2] V. Falkenhahn _et al._ , “Model-based feedforward position control of constant curvature continuum robots using feedback linearization,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2015, pp. 762–767.
  * [3] B. Siciliano _et al._ , _Robotics: modelling, planning and control_. Springer Science & Business Media, 2010.
  * [4] J. Wu _et al._ , “An overview of dynamic parameter identification of robots,” _Robot. and computer-integrated manufacturing_ , vol. 26, no. 5, pp. 414–419, 2010.
  * [5] D. Nguyen-Tuong and J. Peters, “Using model knowledge for learning inverse dynamics,” in _Proc. IEEE Int. Conf. Robot. Autom._ , 2010, pp. 2677–2682.
  * [6] K. Lu, H. Lin, C. W. V. Hui, and K. W. S. Au, “A case study of gravity compensation for da Vinci robotic manipulator: A practical perspective,” in _IEEE Int. Conf. Robot. Autom. Workshop on ”Supervised Autonomy in Surgical Robotics”_ , 2018.
  * [7] M. Haruno _et al._ , “Mosaic model for sensorimotor learning and control,” _Neural computation_ , vol. 13, no. 10, pp. 2201–2220, 2001.
  * [8] S. Schaal _et al._ , “Scalable techniques from nonparametric statistics for real time robot learning,” _Applied Intelligence_ , vol. 17, no. 1, pp. 49–60, 2002.
  * [9] S. M. Khansari-Zadeh _et al._ , “Learning stable nonlinear dynamical systems with gaussian mixture models,” _IEEE Trans. Robot._ , vol. 27, no. 5, pp. 943–957, 2011.
  * [10] J. Kocijan _et al._ , “Gaussian process model based predictive control,” in _Proc. IEEE American Control Conf._ , vol. 3, 2004, pp. 2214–2219.
  * [11] Y. Choi _et al._ , “Local online support vector regression for learning control,” in _in Proc. IEEE Int. Symp. Comput. Intel. Robot. Autom._ , 2007, pp. 13–18.
  * [12] M. Jansen, “Learning an accurate neural model of the dynamics of a typical industrial robot,” in _Int. Conf. Artificial Neural Networks_ , 1994, pp. 1257–1260.
  * [13] W. Yu _et al._ , “PD control of robot with RBF networks compensation,” in _Proc. IEEE Int. Joint Conf. Neural Net._ , vol. 5, 2000, pp. 329–334.
  * [14] N. Yilmaz _et al._ , “Neural network based inverse dynamics identification and external force estimation on the da Vinci Research Kit,” in _Proc. IEEE Int. Conf. Robot. and Autom._ , 2020, pp. 1387–1393.
  * [15] A. S. Polydoros _et al._ , “Real-time deep learning of robotic manipulator inverse dynamics,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robots and Syst._ , 2015, pp. 3442–3448.
  * [16] E. Rueckert _et al._ , “Learning inverse dynamics models in O(n) time with LSTM networks,” in _Proc.IEEE-RAS Int. Conf. on Human. Robot._ , 2017, pp. 811–816.
  * [17] D. Nguyen-Tuong and J. Peters, “Model learning for robot control: a survey,” _Cognitive processing_ , vol. 12, no. 4, pp. 319–340, 2011.
  * [18] C. Della Santina _et al._ , “Model-based dynamic feedback control of a planar soft robot: Trajectory tracking and interaction with the environment,” _Int. J. Robot. Res._ , vol. 39, no. 4, pp. 490–513, 2020\.
  * [19] G. Petkos _et al._ , “Learning multiple models of non-linear dynamics for control under varying contexts,” in _Int. Conf. Artificial Neural Networks_. Springer, 2006, pp. 898–907.
  * [20] J. S. De La Cruz _et al._ , “On-line dynamic model learning for manipulator control,” _IFAC Proceedings Volumes_ , vol. 45, no. 22, pp. 869–874, 2012.
  * [21] R. Camoriano _et al._ , “Incremental semiparametric inverse dynamics learning,” in _Proc. IEEE Int. Conf. Robot. and Autom._ , 2016, pp. 544–550.
  * [22] D. Kappler _et al._ , “A new data source for inverse dynamics learning,” in _Proc. IEEE/RSJ Int. Conf. on Intell. Robots and Syst._ , 2017, pp. 4723–4730.
  * [23] F. D. Ledezma and S. Haddadin, “First-order-principles-based constructive network topologies: An application to robot inverse dynamics,” in _Proc. IEEE-RAS Int. Conf. Human. Robot._ , 2017, pp. 438–445.
  * [24] M. Lutter _et al._ , “Deep Lagrangian Networks: Using Physics as Model Prior for Deep Learning,” in _Proc. Int. Conf. Learn. Repr._ , May 2019\.
  * [25] G. Hinton _et al._ , “Distilling the knowledge in a neural network,” in _NIPS Deep Learning and Representation Learning Workshop_ , 2015.
  * [26] C. Buciluǎ, _et al._ , “Model compression,” in _Proc. of SIGKDD_ , 2006, pp. 535–541.
  * [27] H. Lin _et al._ , “DFNN for GCC,” https://github.com/linhongbin-ws/DFNN4GCC.git, 2020.
  * [28] D. Ma _et al._ , “Identifying mass parameters for gravity compensation and automatic torque sensor calibration,” in _Proc. IEEE Int. Conf. on Robot. and Autom._ , vol. 1, 1996, pp. 661–666.
  * [29] G. A. Fontanelli _et al._ , “Modelling and identification of the da vinci research kit robotic arms,” in _Proc. IEEE/RSJ Int. Conf. Intell. Robots and Syst._ , 2017, pp. 1464–1469.
  * [30] Y. Wang _et al._ , “A Convex Optimization-based Dynamic Model Identification Package for the da Vinci Research Kit,” _IEEE Robot. Autom. Lett._ , pp. 1–1, 2019.
  * [31] H. Lin, C.-W. V. Hui, Y. Wang, A. Deguet, P. Kazanzides, and K. S. Au, “A reliable gravity compensation control strategy for dVRK robotic arms with nonlinear disturbance forces,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 3892–3899, 2019.
  * [32] K. M. Lynch and F. C. Park, _Modern Robotics_. Cambridge University Press, 2017.
  * [33] H. Lin _et al._ , “Supplement of Physics Knowledge Distillation,” https://github.com/linhongbin-ws/DFNN4GCC/blob/controller-evaluation/resource/Supplement.pdf, 2020\.
  * [34] C. M. Bishop _et al._ , _Neural networks for pattern recognition_. Oxford university press, 1995\.
  * [35] M. R. U. Saputra _et al._ , “Distilling knowledge from a deep pose regressor network,” in _Proc. IEEE Int. Conf. Comp. Vision_ , 2019, pp. 263–272.
  * [36] P. Kazanzides _et al._ , “An Open-Source Research Kit for the da Vinci® Surgical System,” in _IEEE/RSJ Int. Conf. Robot. and Autom._ , 2014, pp. 6434–6439.
  * [37] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in _Proceedings of the 27th international conference on machine learning (ICML-10)_ , 2010, pp. 807–814.
  * [38] D. P. Kingma _et al._ , “Adam: A method for stochastic optimization,” _3rd International Conference on Learning Representations (ICLR)_ , 2015\.
  * [39] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT press, 2016.


nestedcomment

# Achieving Optimum Received Power with Elementwise Updates in the Least
Number of Steps for Discrete-Phase RISs ††thanks: The authors are with the
Center for Pervasive Communications and Computing (CPCC), Department of
Electrical Engineering and Computer Science, University of California, Irvine.
††thanks: This work is partially supported by NSF grant 2030029.

Dogan Kutay Pekcan, Student Member, IEEE Ender Ayanoglu, Fellow, IEEE

###### Abstract

The problem of optimizing discrete phases in a reconfigurable intelligent
surface (RIS) to maximize the received power at a user equipment is addressed.
Necessary and sufficient conditions to achieve this maximization are given.
These conditions are employed in an algorithm to achieve the maximization. New
versions of the algorithm are given that are proven to achieve convergence in
$N$ or fewer steps whether the direct link is completely blocked or not, where
$N$ is the number of the RIS elements, whereas previously published results
achieve this in $KN$ or $2N$ number of steps where $K$ is the number of
discrete phases, e.g., [1, 2]. Thus, for a discrete-phase RIS, the techniques
presented in this paper achieve the optimum received power in the smallest
number of steps published in the literature. In addition, in each of those $N$
steps, the techniques presented in this paper determine only one or a small
number of phase shifts with a simple elementwise update rule, which result in
a substantial reduction of computation time, as compared to the algorithms in
the literature, e.g., [3, 2].

###### Index Terms:

Intelligent reflective surface (IRS), reconfigurable intelligent surface
(RIS), discrete beamforming for IRS/RIS.

##  I Introduction

In this paper, we address the problem of finding the values
$\theta_{1},\theta_{2},\ldots,\theta_{N}$ to maximize
$|h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}|$ where $\theta_{n}\in\Phi_{K}$ and
$\Phi_{K}=\\{\omega,2\omega,\ldots,K\omega\\}$ with $\omega=\frac{2\pi}{K}$
and $j=\sqrt{-1}$. The set $\Phi_{K}$ can equivalently be described as
$\\{0,\omega,2\omega,\ldots,(K-1)\omega\\}$. The values $h_{n}\in\mathbb{C}$,
$n=1,2,\ldots,N$ are the channel coefficients and $\theta_{n}$ are the phase
values added to the corresponding $h_{n}$ by a reconfigurable intelligent
surface (RIS).

Towards achieving its goal, [1] introduced the following lemma.

Lemma 1: For an optimal solution $(\theta_{1}^{*},\ldots,\theta_{n}^{*})$ to
problem (8), each $\theta_{n}^{*}$ must satisfy

$\theta_{n}^{*}=\arg\min_{\theta_{n}\in\Phi_{K}}|(\theta_{n}+\alpha_{n}-\phase{\mu})\;{\rm
mod}\;2\pi|$ (11)

where $\phase{\mu}$ stands for the phase of $\mu$ in (10)111To prevent
confusion, we will use the same equation numbers (7)–(13) in [1]. Our own
equation numbers, not available in [1], will begin at (19) and will be
incremented from that number on. Similarly, we will introduce Lemma 2 and
Algorithm 2 in lieu of Lemma 1 and Algorithm 1 in [1]. Note that a lemma or an
algorithm with number 2 does not exist in [1].222In this paper, we define the
${\rm mod}$ function (the modulus function or the modulo operation) $x\;{\rm
mod}\;y$ as the remainder after the dividend $x>0$ is divided by the divisor
$y>0$. We write it as $x\;{\rm mod}\;y$, $x\;(mod\;y)$, or ${\rm mod}\;(x,y)$.
For $x<0$ and $y>0,$ we use the convention that the remainder should always be
the smallest such nonnegative number..

In [1], problem (8) is defined as The problem can be formally described as

$\displaystyle\underset{\mbox{\boldmath$\theta$}}{\rm max\
}\left[f({\mbox{\boldmath$\theta$}})=\frac{1}{\beta_{0}^{2}}\bigg{|}\beta_{0}e^{j\alpha_{0}}+\sum_{n=1}^{N}\beta_{n}e^{j(\alpha_{n}+\theta_{n})}\bigg{|}^{2}\right]$
(1) $\displaystyle{\rm subject\ to\ }\theta_{n}\in\Phi_{K},\ n=1,2,\ldots,N$

where $h_{n}=\beta_{n}e^{j\alpha_{n}}$ for $n=0,1,\ldots,N$, and
${\mbox{\boldmath$\theta$}}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$. Also,
$g$ is defined as

$g=h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}^{*}}$ (2)

and $\mu$ as $\mu=g/|g|$. In (1), $\beta_{0}^{2}$ is a constant and therefore
the maximization affects only the numerator.

Lemma 1 does not hold. This can be seen by numerical examples. We give one
such example in Table 1. In this table, we look at the simple case of $K=2$,
$N=2$. According to Lemma 1 in [1], the condition in (11) should satisfy (8)
for this simple case. We draw values of $h_{n}$ according to the first
paragraph of Sec. IV in [1]. We list these values in rows 2–4 of Table 1. We
define

$g_{0}(\theta_{1},\theta_{2})=h_{0}+\sum_{n=1}^{2}h_{n}e^{j\theta_{n}}$ (19)

and list the values of $g_{0}(\theta_{1},\theta_{2})$ for all possible
$\theta_{1},\theta_{2}\in\\{0,\pi\\}$. There are four such values and they are
listed in rows 5–8 of Table 1. The set of values for $\theta_{1}$ and
$\theta_{2}$ that maximize $|g_{0}|$, or equivalently, that achieve $g$ in
(9), are $\theta_{1}=\theta_{2}=\pi$ as in row 8 of Table 1. Note that this
operation results in $\phase{\mu}=2.3719$ radians as shown in column 5 of row
8 of Table 1.

| ${\rm Re}[\cdot]$ | ${\rm Im}[\cdot]$ | $|\cdot|$ |  $\phase{\cdot}\in[0,2\pi)$ (rad.)
---|---|---|---|---
$h_{0}$ | $-2.8267\times 10^{-7}$ | $2.7376\times 10^{-7}$ | $3.9350\times 10^{-7}$ | $2.3722$
$h_{1}$ | $1.0958\times 10^{-10}$ | $-1.0501\times 10^{-11}$ | $1.1008\times 10^{-10}$ | $6.1876$
$h_{2}$ | $-1.2238\times 10^{-11}$ | $-2.6605\times 10^{-11}$ | $2.6634\times 10^{-10}$ | $4.6664$
$g_{0}(\theta_{1}=0,\theta_{2}=0)$ | $-2.8257\times 10^{-7}$ | $2.7348\times 10^{-7}$ | $3.9324\times 10^{-7}$ | $2.3725$
$g_{0}(\theta_{1}=0,\theta_{2}=\pi)$ | $-2.8255\times 10^{-7}$ | $2.7401\times 10^{-7}$ | $3.9359\times 10^{-7}$ | $2.3715$
$g_{0}(\theta_{1}=\pi,\theta_{2}=0)$ | $-2.8279\times 10^{-7}$ | $2.7350\times 10^{-7}$ | $3.9341\times 10^{-7}$ | $2.3729$
$g_{0}(\theta_{1}=\pi,\theta_{2}=\pi)$ | $-2.8277\times 10^{-7}$ | $2.7403\times 10^{-7}$ | ${\bf 3.9377}\times{\bf 10^{-7}}$ | $\bf 2.3719$
Table 1: Sample calculation for attempting to find optimum $\theta_{1}^{*},\theta_{2}^{*},\ldots,\theta_{N}^{*}$ to maximize $|g_{0}|$ where $g_{0}(\theta_{1},\theta_{2},\ldots,\theta_{N})=h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}$ with $\theta_{n}\in\Phi_{K}=\\{0,\frac{2\pi}{K},\ldots,(K-1)\frac{2\pi}{K}\\}$, $n=1,2,\ldots,N$, for $K=2$ and $N=2$. Channel coefficients $h_{n}$, $n=0,1,2$ are calculated using the technique described in [1]. Rows 5–8 present all values of $g_{0}$ with all combinations of $\theta_{1},\theta_{2}\in\Phi_{2}$, showing that $|g|=\max|g_{0}(\theta_{1},\theta_{2})|$ is achieved with $\theta_{1}^{*}=\theta_{2}^{*}=\pi$.  $(\theta_{1}=0)+\alpha_{1}-\phase{\mu}$ | 3.8158
---|---
${\rm mod}((\theta_{1}=0)+\alpha_{1}-\phase{\mu},2\pi)$ | 3.8158
$(\theta_{1}=\pi)+\alpha_{1}-\phase{\mu}$ | 6.9574
${\rm mod}((\theta_{1}=\pi)+\alpha_{1}-\phase{\mu},2\pi)$ | 0.67417
$(\theta_{2}=0)+\alpha_{2}-\phase{\mu}$ | 2.2945
${\rm mod}((\theta_{2}=0)+\alpha_{2}-\phase{\mu},2\pi)$ | 2.2945
$(\theta_{2}=\pi)+\alpha_{2}-\phase{\mu}$ | 5.4361
${\rm mod}((\theta_{2}=\pi)+\alpha_{2}-\phase{\mu},2\pi)$ | 5.4361
$\cos((\theta_{1}=0)+\alpha_{1}-\phase{\mu})$ | -0.7812
$\cos((\theta_{1}=\pi)+\alpha_{1}-\phase{\mu})$ | 0.7812
$\cos((\theta_{2}=0)+\alpha_{2}-\phase{\mu})$ | -0.6672
$\cos((\theta_{2}=\pi)+\alpha_{2}-\phase{\mu})$ | 0.6672
Table 2: Continuation of the sample calculation for attempting to find optimum
$\theta_{1}^{*},\theta_{2}^{*},\ldots,\theta_{N}^{*}$ to maximize $|g_{0}|$.
Rows 1–8 present the calculation of $\min_{\theta_{n}\in\Phi_{K}}\;{\rm
mod}\;(\theta_{n}+\alpha_{n}-\phase{\mu},2\pi)$ for $n=1,2,\ldots,N$, as
specified in [1] to attempt to find the optimum values of $\theta_{n}$. This
calculation results in values $\theta_{1}=0$ and $\theta_{2}=\pi$, which are
not $\theta_{1}^{*},\theta_{2}^{*}$. Rows 9-12 present the calculation of
$\max_{\theta_{n}\in\Phi_{K}}\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$ to find
$\theta_{1}^{*},\theta_{2}^{*},\ldots,\theta_{N}^{*}$ as discussed in this
comment. This technique finds the optimum values of $\theta_{n},$
$n=1,2,\ldots,N$.

At this point, we would like to emphasize that [1] uses a particular
convention for the phases of complex numbers. They are defined to be in
$[0,2\pi)$, see the text that follows (2) in [1]. We use the same convention
in generating Table 1, see its column 5, as well as in generating Table 2.
With this convention, we list $\theta_{n}+\alpha_{n}-\phase{\mu}$ and
$(\theta_{n}+\alpha_{n}-\phase{\mu})\;{\rm mod}\;2\pi$ for possibilities of
$\theta_{n}=0$ and $\theta_{n}=\pi$ and $n=1,2$ in rows 1–8 of Table 2.333Note
that absolute value signs in (11) are not needed since the argument of the
minimum operation in (11) is in $[0,2\pi)$. It can be seen from rows 1–4 of
Table 2 that the method results in $\theta_{1}=\pi$ as the potential
$\theta_{1}^{*}$, which we know from the discussion in the previous paragraph
to be correct. When we carry out the calculation
$(\theta_{2}+\alpha_{2}-\phase{\mu})\;{\rm mod}\;2\pi$ in rows 5–8 of Table 2,
we find that the method suggests $\theta_{2}=0$ should be $\theta_{2}^{*}$.
However, we know from the exhaustive search in rows 5–8 of Table 1 that
$\theta_{2}^{*}=\pi$. Thus, Lemma 1 is not correct.

It is possible to come up with a correct lemma similar to Lemma 1. We specify
this lemma below. The solution to the problem can be achieved by making use of
the following lemma.

Lemma: For an optimal solution
$(\theta_{1}^{*},\theta_{2}^{*},\ldots,\theta_{n}^{*})$, it is necessary and
sufficient that each $\theta_{n}^{*}$ satisfy

$\theta_{n}^{*}=\arg\max_{\theta_{n}\in\Phi_{K}}\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$
(20)

where $\phase{\mu}$ stands for the phase of $\mu=g/|g|$.

Proof: We can rewrite (2) as

$\displaystyle|g|=$ $\displaystyle\
\beta_{0}e^{j(\alpha_{0}-\phase{\mu})}+\sum_{n=1}^{N}\beta_{n}e^{j(\alpha_{n}+\theta_{n}-\phase{\mu})}$
(21) $\displaystyle=$ $\displaystyle\
\beta_{0}\cos(\alpha_{0}-\phase{\mu})+j\beta_{0}\sin(\alpha_{0}-\phase{\mu})$
$\displaystyle+\sum_{n=1}^{N}\beta_{n}\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$
$\displaystyle+j\sum_{n=1}^{N}\beta_{n}\sin(\theta_{n}+\alpha_{n}-\phase{\mu}).$
(22)

Because $|g|$ is real-valued, the second and fourth terms in (22) sum to zero,
and

$|g|=\beta_{0}\cos(\alpha_{0}-\phase{\mu})+\sum_{n=1}^{N}\beta_{n}\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$
(23)

from which (20) follows as a necessary and sufficient condition for the lemma
to hold. $\blacksquare$

Reference [1] attempts to decide a range of $\mu$ for which
$\theta_{n}^{*}=k\omega$ must hold. Towards that end, it first defines a
sequence of complex numbers with respect to each $n=1,2,\ldots,N$ as

$s_{nk}=e^{j(\alpha_{n}+(k-0.5)\omega)},\ {\rm for}\ k=1,2,\ldots,K.$ (24)

Then, [1] defines, for any two points $a$ and $b$ on the unit circle $C$,
${\rm arc}(a:b)$ to be the unit circular arc with $a$ as the initial end and
$b$ as the terminal end in the counterclockwise direction; in particular, it
defines ${\rm arc}(a:b)$ as an open arc with the two endpoints $a$ and $b$
excluded. With this definition, [1] states the following proposition holds.

Proposition 1: A sufficient condition for $\theta_{n}^{*}=k\omega$ is

$\mu\in{\rm arc}(s_{nk}:s_{n,k+1}).$ (25)

Proposition 1 is compatible with the lemma given above. To see this, assume
$\mu$ satisfies (25). Then,

$\phase{\mu}\in\Big{(}\alpha_{n}+\Big{(}k-\frac{1}{2}\Big{)}\omega,\alpha_{n}+\Big{(}k+\frac{1}{2}\Big{)}\omega\Big{)}.$
(26)

Since $\omega=\frac{2\pi}{K}$,

$\alpha_{n}-\phase{\mu}\in\Big{(}(-2k-1)\frac{\pi}{K},(-2k+1)\frac{\pi}{K}\Big{)}$
(27)

considering the reversal of order due to the substraction of $\phase{\mu}$.
Now, let $\theta_{n}=k\omega=2k\frac{\pi}{K}$. Then

$\theta_{n}+\alpha_{n}-\phase{\mu}\in\Big{(}-\frac{\pi}{K},\frac{\pi}{K}\Big{)}$
(28)

and thus $\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$ is the largest among all
other possibilities for $\theta_{n}$ because the slice
$(-\frac{\pi}{K},\frac{\pi}{K})$ corresponds to the largest values of the
cosine function among all slices corresponding to different values of
$\theta_{n}\in\Phi_{K}$ for $n=1,2,\ldots,N$.

We note that $g$ in (2) is defined for optimal phases
$(\theta_{1}^{*},\theta_{2}^{*},\ldots,\theta_{N}^{*})$. The $\mu$ in (20),
$\mu=g/|g|$, which follows from (2), is the optimum one. On the other hand, in
the rest of the document, when we refer to $\mu$, it is a value we are
considering in search of the optimal $\mu$.

##  II A New Algorithm

Algorithm 2 Update for Algorithm 1 [1]

1:Initialization: Compute $s_{nk}=e^{j(\alpha_{n}+(k-0.5)\omega)}$ for
$n=1,2,\ldots,N$ and $k=1,2,\ldots,K$.

2:Eliminate duplicates among $s_{nk}$ and sort to get $e^{j\lambda_{l}}$ such
that $0\leq\lambda_{1}<\lambda_{2}<\cdots<\lambda_{L}<2\pi.$

3:Let, for $l=1,2,\ldots,L,$ ${\cal
N}(\lambda_{l})=\\{n|\phase{s_{nk}}=\lambda_{l}\\}.$

4:Set $\phase{\mu}=0$. For $n=1,2,\ldots,N$, calculate
$\theta_{n}=\arg\max_{\theta_{n}\in\Phi_{K}}\cos(\theta_{n}+\alpha_{n}-\phase{\mu})$.

5:Set $g_{0}=h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}$, ${\tt
absgmax}=|g_{0}|$.

6:for $l=1,2,\ldots,L-1$ do

7: For each $n\in{\cal N}(\lambda_{l})$, let
$(\theta_{n}+\omega\leftarrow\theta_{n})\;{\rm mod}\;\Phi_{K}$.

8: Let  $g_{l}=g_{l-1}+\sum_{n\in{\cal
N}(\lambda_{l})}h_{n}\big{(}e^{j\theta_{n}}-e^{j(\theta_{n}-\omega)\;{\rm
mod}\;\Phi_{K}}\big{)}$

9: if $|g_{l}|>{\tt absgmax}$ then

10: Let ${\tt absgmax}=|g_{l}|$

11: Store $\theta_{n}$ for $n=1,2,\ldots,N$

12: end if

13:end for

14:Read out $\theta_{n}^{*}$ as the stored $\theta_{n}$, $n=1,2,\ldots,N$.

Algorithm 1 in [1], as published, employs the criterion

$\theta_{n}^{*}=\arg\min_{\theta_{n}\in\Phi_{K}}|(\theta_{n}+\alpha_{n}-\phase{\mu})\;{\rm
mod}\;2\pi|$ (29)

instead of (20)444In this paper, we define the ${\rm mod}$ function (the
modulus function or the modulo operation) $x\;{\rm mod}\;y$ as the remainder
after the dividend $x>0$ is divided by the divisor $y>0$. We write it as
$x\;{\rm mod}\;y$, $x\;(mod\;y)$, or ${\rm mod}\;(x,y)$. For $x<0$ and $y>0,$
we use the convention that the remainder should always be the smallest such
nonnegative number.. We now specify Algorithm 2 as an alternative to Algorithm
1 in [1]555In Algorithm 2 we define $(\theta\pm\omega)\,{\rm mod}\,\Phi_{K}$
as follows. First note that the two sets
$\\{0,\omega,2\omega,\dots,(K-1)\omega\\}$ and
$\\{\omega,2\omega,3\omega,\ldots,K\omega\\}$ have the same members since
$\omega=2\pi/K$. Then, $(\theta+\omega)\,{\rm mod}\,\Phi_{K}$ can be defined
as $(\theta\pm\omega)\,{\rm mod}\,\Phi_{K}\triangleq((k\pm 1)\,{\rm
mod}\,K)\,\omega$..

Figure 1: CDF plots for SNR Boost [1] with Closest Point Projection (CPP) and
Approximation (APX) algorithms [4] and Algorithm 2, $K=2$.

We present the CDF results for SNR Boost[1] in Fig. 1 for $K=2$ and $N=16,$
$64,$ and $256,$ using the average of 10,000 realizations of the channel,
where we employed Closest Point Projection (CPP) and Approximation (APX)
algorithms from [4] and Algorithm 2, and where all algorithms ran over the
same realization in each step. Although the gains are not large, the figure
serves as a verification of the optimality of Algorithm 2, which we already
know from the analysis presented in this paper.

In Appendix, we discuss an alternative way to initialize Algorithm 2. In
addition to Algorithm 2, we will use this technique in initializing Algorithm
3 and Algorithm 4 in the sequel.

##  III Algorithm Convergence: Towards Two New Algorithms

We will now show the periodicity in the update rule in Algorithm 2, i.e., in
${\cal N}(\lambda_{l})$. With this observation, we will prove that the
algorithm can be made to take $N$ or fewer steps to converge, as opposed to
the statement in [1] that it takes $KN$ or $2N$ steps on average. This will
result in two new versions of the algorithm, i.e., Algorithm 3 and Algorithm
4, with a simple elementwise update rule, in the sequel. Towards this end, we
first make the following statement.

Claim 1: As in (24), set
$s_{nk}=e^{j(\alpha_{n}+(k-\frac{1}{2})\frac{2\pi}{K})}$, $n=1,2,\ldots,N,$
$k=1,2,\ldots,K$, $\alpha_{n}\in[0,2\pi).$ Let $\lambda_{l}=\phase{s_{nk}}$
such that $0\leq\lambda_{1}<\lambda_{2}<\cdots<\lambda_{L}<2\pi$. Let ${\cal
N}(\lambda_{l})=\\{n|\lambda_{l}=\phase{s_{nk}}\\}$. Assuming for now that
$|{\cal N}(\lambda_{l})|=1$, $l=1,2\ldots,L=NK$, which we will relax in the
sequel, we claim that ${\cal N}(\lambda_{l^{\prime}})={\cal
N}(\lambda_{l^{\prime}+N})$ for $l^{\prime}=1,2,\ldots,N(K-1).$

To prove Claim 1, we will first introduce Claim 2 and prove it.

Claim 2: Without loss of generality, we can assume that
$\alpha_{n}<\frac{2\pi}{K},$ $n=1,2,\ldots,N$.

Proof of Claim 2: Suppose that for some $n$, we have $\frac{2\pi}{K}\cdot
m\leq\alpha_{n}<\frac{2\pi}{K}\cdot(m+1),$ $m=1,2,\ldots,K-1$. Let
$\beta_{n}\triangleq\alpha_{n}-m\cdot\frac{2\pi}{K},$ so that
$\beta_{n}<\frac{2\pi}{K}$. We will write below each phase value in (24),
$\alpha_{n}+\frac{(2k-1)\pi}{K}$, for $k=1,2,\ldots,K$ (note that the $({\rm
mod}\;2\pi)$ notation below applies to both sides of the equation).

  * •

$k=1:$
$\alpha_{n}+\frac{\pi}{K}=\beta_{n}+\frac{(2m+1)\pi}{K}\ ({\rm mod}\;2\pi)$
where
$\Big{(}\frac{(2m+1)\pi}{K}\Big{)}_{m=1}^{K-1}=\Big{\\{}\frac{3\pi}{K},\frac{5\pi}{K},\ldots,\frac{(2K-1)\pi}{K}\Big{\\}},$

  * •

$k=2:$
$\alpha_{n}+\frac{3\pi}{K}=\beta_{n}+\frac{(2m+3)\pi}{K}\ ({\rm mod}\;2\pi)\\\
=\left\\{\begin{tabular}[]{ll}$\beta_{n}+\frac{(2m+3)\pi}{K}\ ({\rm
mod}\;2\pi),$&$m\leq K-2$\\\ $\beta_{n}+\frac{\pi}{K}\ ({\rm
mod}\;2\pi),$&$m=K-1$\end{tabular}\right.$
where
$\Big{(}\frac{(2m+3)\pi}{K}\Big{)}_{m=1}^{K-2}=\Big{\\{}\frac{5\pi}{K},\frac{7\pi}{K},\ldots,\frac{(2K-1)\pi}{K}\Big{\\}},$
$\vdots$

  * •

$k=K:$
$\alpha_{n}+\frac{(2K-1)\pi}{K}=\beta_{n}+\frac{(2m+2K-1)\pi}{K}\ ({\rm
mod}\;2\pi)\\\ \mbox{\hskip 59.18173pt}=\beta_{n}+\frac{(2m-1)\pi}{K}\ ({\rm
mod}\;2\pi)$
where
$\Big{(}\frac{(2m-1)\pi}{K}\Big{)}_{m=1}^{K-1}=\Big{\\{}\frac{\pi}{K},\frac{3\pi}{K},\ldots,\frac{(2K-3)\pi}{K}\Big{\\}}.$

Thus, if there is an $\alpha_{n}\geq\frac{2\pi}{K}$ to generate $K$ phase
values, there is always a $\beta_{n}$, $\beta_{n}<\frac{2\pi}{K}$ with which
one can generate the same $K$ phase values in a similar fashion. Therefore, in
order to prove Claim 1, one can work with the assumption that
$\alpha_{n}<\frac{2\pi}{K},$ for $n=1,2,\ldots,N$. $\blacksquare$

Proof of Claim 1: Assuming
$0\leq\alpha_{1}<\alpha_{2}<\cdots<\alpha_{N}<\frac{2\pi}{K},$ without loss of
generality, we will now show that ${\cal N}(\lambda_{l^{\prime}})={\cal
N}(\lambda_{l^{\prime}+N})$ for $l^{\prime}=1,2,\ldots,N(K-1).$ For this,
there are $N+1$ cases to consider.

Case 0: In this case, we assume $\alpha_{n}<\frac{\pi}{K},$ $n=1,2,\ldots,N$.
We write all possible values of $\phase{s_{nk}}$ as follows.

$n=1:$

$\phase{s_{1k}}\in\textstyle\left\\{\alpha_{1}+\frac{\pi}{K},\alpha_{1}+\frac{3\pi}{K},\ldots,\alpha_{1}+\frac{(2K-1)\pi}{K}\right\\},$
(30)

$n=2:$

$\phase{s_{2k}}\in\textstyle\left\\{\alpha_{2}+\frac{\pi}{K},\alpha_{2}+\frac{3\pi}{K},\ldots,\alpha_{2}+\frac{(2K-1)\pi}{K}\right\\},$
(31)

$\vdots$

$n=N:$

$\hskip
11.38109pt\phase{s_{Nk}}\in\textstyle\left\\{\alpha_{N}+\frac{\pi}{K},\alpha_{N}+\frac{3\pi}{K},\ldots,\alpha_{N}+\frac{(2K-1)\pi}{K}\right\\}.$
(32)

Sorting (30)–(32), we have

$\displaystyle\textstyle\alpha_{1}+\frac{\pi}{K}<\alpha_{2}+\frac{\pi}{K}<\cdots<\alpha_{N}+\frac{\pi}{K}$
$\displaystyle\\!\\!\\!\\!\\!\\!<\textstyle\alpha_{1}+\frac{3\pi}{K}<\alpha_{2}+\frac{3\pi}{K}<\cdots<\alpha_{N}+\frac{3\pi}{K}$
$\displaystyle\\!\\!\\!\\!\\!\\!\hskip 5.69054pt\vdots$
$\displaystyle\\!\\!\\!\\!\\!\\!<\textstyle\alpha_{1}+\frac{(2K-1)\pi}{K}<\alpha_{2}+\frac{(2K-1)\pi}{K}<\cdots<\alpha_{N}+\frac{(2K-1)\pi}{K}.$

Thus,

$\big{(}{\cal
N}(\lambda_{l})\big{)}_{l=1}^{L}=\Big{\\{}\underbrace{1,2,\ldots,N}_{1},\underbrace{1,2,\ldots,N}_{2},\ldots,\underbrace{1,2,\ldots,N}_{K}\Big{\\}}.$

where $L=NK.$ Therefore, for Case 0 and for $l^{\prime}=1,2,\ldots,N(K-1),$ we
have ${\cal N}(\lambda_{l^{\prime}})={\cal N}(\lambda_{l^{\prime}+N})$.

There are $N$ remaining cases. We will discuss these cases as Case $i$ where
$i=1,2,\ldots,N$.

Case $i,(i=1,2,\ldots,N)\\!\\!:$ In Case $i$, we have $i$ occurrences of
$\alpha_{n}>\frac{\pi}{K}$ as follows.

$\displaystyle 0$
$\displaystyle\leq\alpha_{1}<\alpha_{2}<\cdots<\alpha_{N-i}<\frac{\pi}{K}$
(33) $\displaystyle\leq\alpha_{N-i+1}<\cdots<\alpha_{N}<\frac{2\pi}{K}.$

We write all possible values of $\phase{s_{nk}}$ as follows.

$n=1:$

$\displaystyle\phase{s_{1k}}\in\textstyle\Big{\\{}\alpha_{1}+$
$\displaystyle\frac{\pi}{K},\alpha_{1}+\frac{3\pi}{K},\ldots,$ (34)
$\displaystyle\alpha_{1}+\frac{(2K-3)\pi}{K},\alpha_{1}+\frac{(2K-1)\pi}{K}\Big{\\}},$

$\hskip 14.22636pt\vdots$

$n=N-i:$

$\displaystyle\phase{s_{(N-i)k}}\in\textstyle\Big{\\{}$
$\displaystyle\alpha_{N-i}+\frac{\pi}{K},\alpha_{N-i}+\frac{3\pi}{K},\ldots,$
(35)
$\displaystyle\alpha_{N-i}+\frac{(2K-3)\pi}{K},\alpha_{N-i}+\frac{(2K-1)\pi}{K}\Big{\\}},$

$n=N-i+1:$

$\displaystyle\phase{s_{(N-i+1)k}}\in\textstyle\Big{\\{}$
$\displaystyle\alpha_{N-i+1}+\frac{\pi}{K},\alpha_{N-i+1}+\frac{3\pi}{K},\ldots,$
(36)
$\displaystyle\alpha_{N-i+1}+\frac{(2K-3)\pi}{K},\alpha_{N-i+1}-\frac{\pi}{K}\Big{\\}},$

$n=N-i+2:$

$\displaystyle\phase{s_{(N-i+2)k}}\in\textstyle\Big{\\{}$
$\displaystyle\alpha_{N-i+2}+\frac{\pi}{K},\alpha_{N-i+2}+\frac{3\pi}{K},\ldots,$
(37)
$\displaystyle\alpha_{N-i+2}+\frac{(2K-3)\pi}{K},\alpha_{N-i+2}-\frac{\pi}{K}\Big{\\}},$

$\hskip 14.22636pt\vdots$

$n=N:$

$\displaystyle\phase{s_{Nk}}\in\Big{\\{}\alpha_{N}+\frac{\pi}{K},$
$\displaystyle\alpha_{N}+\frac{3\pi}{K},\ldots,$ (38)
$\displaystyle\alpha_{N}+\frac{(2K-3)\pi}{K},\alpha_{N}-\frac{\pi}{K}\Big{\\}}.$

Sorting (34)–(38), we have

$\alpha_{N-i+1}-\frac{\pi}{K}<\cdots<\alpha_{N}-\frac{\pi}{K}\hskip
76.82243pt$ $\hskip
85.35826pt<\alpha_{1}+\frac{\pi}{K}<\cdots<\alpha_{N-i}+\frac{\pi}{K}$
$<\alpha_{N-i+1}+\frac{\pi}{K}<\cdots<\alpha_{N}+\frac{\pi}{K}\hskip
85.35826pt$ $\hskip
85.35826pt<\alpha_{1}+\frac{3\pi}{K}<\cdots<\alpha_{N-i}+\frac{3\pi}{K}$
$<\alpha_{N-i+1}+\frac{3\pi}{K}<\cdots<\alpha_{N}+\frac{3\pi}{K}\hskip
85.35826pt$ $\hskip
85.35826pt<\alpha_{1}+\frac{5\pi}{K}<\cdots<\alpha_{N-i}+\frac{5\pi}{K}$
$\vdots$
$<\alpha_{N-i+1}+\frac{(2K-3)\pi}{K}<\cdots<\alpha_{N}+\frac{(2K-3)\pi}{K}\hskip
2.84526pt$ $\hskip
42.67912pt<\alpha_{1}+\frac{(2K-1)\pi}{K}<\cdots<\alpha_{N-i}+\frac{(2K-1)\pi}{K}$

Thus,

$\displaystyle\big{(}{\cal N}(\lambda_{l})\big{)}$ ${}_{l=1}^{L=NK}=$
$\displaystyle\Big{\\{}$
$\displaystyle\underbrace{N-i+1,\ldots,N,1,2,\ldots,N-i}_{1},$
$\displaystyle\underbrace{N-i+1,\ldots,N,1,2,\ldots,N-i}_{2},$ (39)
$\displaystyle\hskip 69.70924pt\vdots$
$\displaystyle\underbrace{N-i+1,\ldots,N,1,2,\ldots,N-i}_{K}\Big{\\}}$

for $i=1,2,\ldots,N$. Therefore, for Case $i$, $i=1,2,\ldots,N$, and for
$l^{\prime}=1,2,\ldots,N(K-1),$ we have ${\cal N}(\lambda_{l^{\prime}})={\cal
N}(\lambda_{l^{\prime}+N})$. With this, Claim 1 is proved. $\blacksquare$

##  IV $N$ Steps Suffice When $|{\cal N}(\lambda_{l})|=1$ for All $l$

Algorithm 3 Simplified Algorithm 2 with $|{\cal N}(\lambda_{l})|=1$ for all
$l$

1:Initialization: Set $\phase{\mu}=\alpha_{0}-\frac{\pi}{K}$

2:Compute $\varphi_{n}=(\alpha_{n}-\alpha_{0})$
$\mathrm{mod}\,\frac{2\pi}{K}$, $n=1,2,\ldots,N$

3:Sort $\varphi_{n}$ such that
$0\leq\varphi_{n_{1}}<\varphi_{n_{2}}<\cdots<\varphi_{n_{N}}<\frac{2\pi}{K}$

4:Set $\theta_{n}={\rm
argmax}_{\theta_{n}^{\prime}\in\Phi_{K}}\cos(\theta_{n}^{\prime}+\alpha_{n}-\phase{\mu})$,
store $\theta_{n}$, $n=1,2,\ldots,N$

5:Set $g_{0}=h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}$, ${\tt
absgmax}=|g_{0}|$

6:for $l=1,2,\ldots,N$ do

7: Let
$(\theta_{n_{l}}+\omega\leftarrow\theta_{n_{l}})\;\mathrm{mod}\,\Phi_{K}$

8: Let
$g_{l}=g_{l-1}+h_{n_{l}}\left(e^{j\theta_{n_{l}}}-e^{j(\theta_{n_{l}}-\omega)\,\mathrm{mod}\,\Phi_{K}}\right)$

9: if $|g_{l}|>{\tt absgmax}$ then

10: Let ${\tt absgmax}=|g_{l}|$

11: Store updated $\theta_{n_{l}}$

12: end if

13:end for

14:Read out $\theta_{n}^{*}$ as the stored $\theta_{n}$, $n=1,2,\ldots,N$.

Given $|{\cal N}(\lambda_{l})|=1$ and

${\cal N}(\lambda_{l^{\prime}})={\cal N}(\lambda_{l^{\prime}+N}),\quad
l^{\prime}=1,2,\dots,N(K-1),$ (40)

we want to show that $N$ steps will suffice for convergence. Now, consider the
main problem of maximizing
$\left|h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}\right|$, where it is clear
that our discrete phase shift selections can only tune the second term in the
absolute value. Let

$g_{c}\triangleq\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}=\sum_{n=1}^{N}\beta_{n}e^{j(\alpha_{n}+\theta_{n})}.$
(41)

In each step of the Algorithm 2, we define

$g_{c,l}\triangleq g_{l}-h_{0},\quad l=1,2,\ldots,L.$ (42)

Note that, in (42), $h_{0}$, $g_{l}$, and $g_{c,l}$ are complex numbers, with
$l$ being a generation index. We know that whenever $\mu$ is anywhere in ${\rm
arc}(s_{nk}:s_{n,k+1})$, $\theta_{n}$ does not change. The angle $\theta_{n}$
only changes when $\mu$ changes from one arc to another, i.e.,

$\mu\in{\rm arc}{(e^{j\lambda_{l}}:e^{j\lambda_{l+1}})}\rightarrow\mu\in{\rm
arc}(e^{j\lambda_{l+1}}:e^{j\lambda_{l+2}})$ (43)

in which case $\theta_{n}$ must be updated as

$\theta_{n}\rightarrow\theta_{n}+\omega,\quad n\in{\cal N}(\lambda_{l+1}).$
(44)

With (43)–(44), the naive approach in Algorithm 2 gathers all possibilities
for $g_{l}$ in $NK$ steps by considering all possible arcs that $\mu$ can be
in. To show that $N$ steps will suffice, we want to point out the redundancy
in those $NK$ steps. Consider any consecutive $N$ steps in Algorithm 2. In
those steps, the phase shifts will be updated as
$\theta_{n}\rightarrow\theta_{n}+\omega$ with $n\in({\cal
N}(\lambda_{l}))_{l=l^{\prime}}^{l^{\prime}+N}$,
$l^{\prime}=1,2,\ldots,N(K-1)$. Since we have (40), the following must hold

$({\cal N}(\lambda_{l}))_{l=l^{\prime}}^{l^{\prime}+N}=\\{1,2,\ldots,N\\},$
(45)

which says that after any $N$ consecutive steps in Algorithm 2,
$\theta_{n}\rightarrow\theta_{n}+\omega,$ $n=1,2,\ldots,N$. To proceed
further, we need an intermediate result, which we discuss below.

Claim 3: Let $g_{c}=\sum_{n=1}^{N}\beta_{n}e^{j(\theta_{n}+\alpha_{n})}$ be
the cascaded channel term in (41). For any angle $\theta^{\prime}$,

$|g_{c}(\theta_{1},\theta_{2},\ldots,\theta_{n})|=|g_{c}(\theta_{1}+\theta^{\prime},\theta_{2}+\theta^{\prime},\ldots,\theta_{N}+\theta^{\prime})|.$
(46)

Proof of Claim 3: Let us write $|g_{c}|^{2}$ as follows

$\displaystyle\\!\\!\\!\\!|g_{c}|^{2}\\!\\!\\!\\!$
$\displaystyle\\!\\!\\!=\\!\\!\\!$
$\displaystyle\left|\sum_{n=1}^{N}\beta_{n}e^{j(\theta_{n}+\alpha_{n})}\right|^{2}=\sum_{n=1}^{N}\beta_{n}^{2}$
(47)
$\displaystyle+\;2\sum_{k=2}^{N}\sum_{l=1}^{k-1}\beta_{k}\beta_{l}\cos((\theta_{k}+\alpha_{k})-(\theta_{l}+\alpha_{l}))$

If $\theta_{n}\rightarrow\theta_{n}+\theta^{\prime},$ $n=1,2,\ldots,N$, then
the angles $\theta^{\prime}$ will cancel inside the cosine in (47). This
proves Claim 3. $\hfill\blacksquare$

With Claim 3,  we established that if $\theta_{n}\rightarrow\theta_{n}+\omega$
with $n=1,2,\ldots,N$,

$|g_{c,l^{\prime}}|=|g_{c,l^{\prime}+N}|,\quad l^{\prime}=1,2,\ldots,N(K-1)$
(48)

must be true. Therefore, among the $NK$ possibilities in Algorithm 2, there
are only $N$ unique values of $|g_{c}|$. Consequently, as the algorithm is
tuning $g_{c}$ to maximize $|g_{c}+h_{0}|$, it is sufficient to consider $N$
arcs that are closest to $h_{0}$.

The algorithm to implement when $|{\cal N}(\lambda_{l})|=1$ for
$l=1,2,\ldots,L=NK$ is given under Algorithm 3. The initialization technique
introduced in Appendix for Algorithm 2 is employed in Algorithm 3.

##  V Fewer Than $N$ Steps Suffice When $|{\cal N}(\lambda_{l})|>1$ for Some
$l$

With ${\cal N}(\lambda_{l})$ given in Claim 1, for $|{\cal N}(\lambda_{l})|>1$
to be true for some $l$, consider a repetition among $s_{nk}$, i.e., assume
there are $n_{1},$ $n_{2},$ $k_{1},$ and $k_{2}$ such that
$s_{n_{1},k_{1}}=s_{n_{2},k_{2}}$, i.e.,

$\displaystyle\left(\alpha_{n_{1}}+\frac{(2k_{1}-1)\pi}{K}\right)\;{\rm
mod}\;2\pi=$ (49)
$\displaystyle\quad\quad\quad\quad\quad\left(\alpha_{n_{2}}+\frac{(2k_{2}-1)\pi}{K}\right)\;{\rm
mod}\;2\pi.$

Equation (49) is possible only if $\beta_{n_{1}}=\beta_{n_{2}}$ as $\beta_{n}$
are defined in Claim 1. Therefore, all $K$ phase values represented by
$\beta_{n_{1}}$ and $\beta_{n_{2}}$ must be equal, meaning there is an
$N^{\prime}$ such that $M=(N-N^{\prime})K$. Consequently, the problem of
sorting $s_{nk}$ according to their phase values with
$0\leq\alpha_{1}<\alpha_{2}<\cdots<\alpha_{N}<\frac{2\pi}{K}$ reduces to
$0\leq\gamma_{1}<\gamma_{2}<\cdots<\gamma_{\frac{M}{K}=N-N^{\prime}}<\frac{2\pi}{K}$
for

$s_{mk}=e^{j(\gamma_{m}+(k-\frac{1}{2})\frac{2\pi}{K})},m=1,2,\ldots,\frac{M}{K},k=1,2,\ldots,K$
(50)

where
$\gamma_{m}=\min\\{\alpha_{n_{m},1},\alpha_{n_{m},2},\ldots,\alpha_{n_{m},G_{m}}\\}$
such that $\beta_{n_{m},1}=\beta_{n_{m},2}=\cdots=\beta_{n_{m},G_{m}}$. So,
this time, there are $\frac{M}{K}+1=N-N^{\prime}+1$ many cases.

For unique $s_{mk}$, let ${\cal
M}(\lambda_{l})=\\{m|\lambda_{l}=\phase{s_{mk}}\\}$. We know from (39) that
the following must hold

$\displaystyle\big{(}{\cal M}(\lambda_{l})\big{)}$
${}_{l=1}^{M=(N-N^{\prime})K}=$ $\displaystyle\Big{\\{}$
$\displaystyle\underbrace{\frac{M}{K}-i+1,\ldots,\frac{M}{K},1,2,\ldots,\frac{M}{K}-i}_{1},$
$\displaystyle\underbrace{\frac{M}{K}-i+1,\ldots,\frac{M}{K},1,2,\ldots,\frac{M}{K}-i}_{2},$
(51) $\displaystyle\hskip 75.3998pt\vdots$
$\displaystyle\underbrace{\frac{M}{K}-i+1,\ldots,\frac{M}{K},1,2,\ldots,\frac{M}{K}-i}_{K}\Big{\\}}$

for $i=1,2,\ldots,\frac{M}{K}$ where in each one of the $K$ groups there are
$\frac{M}{K}=N-N^{\prime}$ elements. To calculate ${\cal N}(\lambda_{l})$, we
define the following sets

${\cal R}_{m}=\\{n_{m,1},n_{m,2},\ldots,n_{m,G_{m}}|\hskip 142.26378pt$
$\gamma_{m}=\min\\{\alpha_{n_{m},1},\alpha_{n_{m},2},\ldots,\alpha_{n_{m},G_{m}}\\},\hskip
28.45274pt$ (52) $\hskip
56.9055pt\beta_{n_{m},1}=\beta_{n_{m},2}=\cdots=\beta_{n_{m},G_{m}}\\}$

where $G_{m}=|{\cal R}_{m}|$ and $\bigcup_{m=1}^{\frac{M}{K}}{\cal
R}_{m}=\\{1,2,\ldots,N\\}$ must hold. As a consequence, one can calculate
${\cal N}(\lambda_{l})={\cal R}_{{\cal M}(\lambda_{l})}$. Therefore, the
“update loop” in Algorithm 2 can be written as

$\displaystyle\big{(}{\cal N}(\lambda_{l})\big{)}$
${}_{l=1}^{M=(N-N^{\prime})K}=$ $\displaystyle\Big{\\{}$
$\displaystyle\underbrace{{\cal R}_{\frac{M}{K}-i+1},\ldots,{\cal
R}_{\frac{M}{K}},{\cal R}_{1},{\cal R}_{2},\ldots,{\cal
R}_{\frac{M}{K}-i}}_{1},$ $\displaystyle\underbrace{{\cal
R}_{\frac{M}{K}-i+1},\ldots,{\cal R}_{\frac{M}{K}},{\cal R}_{1},{\cal
R}_{2},\ldots,{\cal R}_{\frac{M}{K}-i}}_{2},$ (53) $\displaystyle\hskip
75.3998pt\vdots$ $\displaystyle\underbrace{{\cal
R}_{\frac{M}{K}-i+1},\ldots,{\cal R}_{\frac{M}{K}},{\cal R}_{1},{\cal
R}_{2},\ldots,{\cal R}_{\frac{M}{K}-i}}_{K}\Big{\\}}$

where the periodicity in the update rule still holds in (53), i.e., ${\cal
N}(\lambda_{l^{\prime}})={\cal N}(\lambda_{l^{\prime}+\frac{M}{K}}),$
$l^{\prime}=1,2,\ldots,\frac{M}{K}(K-1)$. With the new update rule, after any
$\frac{M}{K}$ consecutive steps in Algorithm 2, the phase shift selections
will be updated such that $\theta_{n}\rightarrow\theta_{n}+\omega,$
$n=1,2,\ldots,N$. This will result in

$|g_{c,l^{\prime}}|=\left|g_{c,l^{\prime}+\frac{M}{K}}\right|,\quad
l^{\prime}=1,2,\ldots,\frac{M}{K}(K-1).$ (54)

Therefore, the sufficiency of $\frac{M}{K}=N-N^{\prime}$ steps follows from
(48) and the text that follows it in Section IV.

Algorithm 4 implements the technique described in this section. The
initialization technique introduced in Appendix for Algorithm 2 is employed in
Algorithm 4.

Note that if the BS-UE link is completely blocked, i.e., $h_{0}=0$, the for
loop in Step 7 can end at $l=\frac{M}{K}-1=N-N^{\prime}-1$, which is one fewer
step to run Algorithm 4. This is because, we can guarantee in (54) that
$|g_{c,l^{\prime}}|=\left|g_{c,l^{\prime}+\frac{M}{K}}\right|,$ whereas we
cannot say right away that
$|g_{c,l^{\prime}}+h_{0}|=\left|g_{c,l^{\prime}+\frac{M}{K}}+h_{0}\right|$
will be satisfied. In [1], the authors reduce the number of steps from $KN$ to
$2N$ exploiting $h_{0}$. In this work, with the periodicity proof, we reduce
to $N$ steps whether the direct link is blocked or not.

Step 7 in Algorithm 3 and Step 8 in Algorithm 4 are such that the phase shifts
updates are restricted to just one or just a few. These steps ensure that the
running times of these algorithms are much less than those published in the
literature, e.g., [2].

Figure 2 shows the behaviors of $|g_{l}|^{2}$ and $|g_{c,l}|^{2}$ against $l$
for a number of channel realizations at $N=16$ and $K=4$. The periodicity of
$|g_{c,l}|^{2}$ is clearly observable in this figure. It is this periodicity
that we take advantage of in reducing the number of steps for the algorithm to
converge to $N$ with a simple elementwise update rule in this paper.

Algorithm 4 Simplified Algorithm 2 where $|{\cal N}(\lambda_{l})|>1$ for some
$l$

1:Initialization: Set $\phase{\mu}=\alpha_{0}-\frac{\pi}{K}$

2:Find $\gamma_{m}$ and ${\cal R}_{m}$ as in (50) and (52),
$m=1,2,\ldots,\frac{M}{K}$

3:Compute $\varphi_{m}=\gamma_{m}-\alpha_{0}$ $(\mathrm{mod}\frac{2\pi}{K})$,
$m=1,2,\ldots,\frac{M}{K}$

4:Sort $\varphi_{m}$ such that
$0\leq\varphi_{m_{1}}<\varphi_{m_{2}}<\cdots<\varphi_{m_{\frac{M}{K}}}<\frac{2\pi}{K}$

5:Set $\theta_{n}={\rm
argmax}_{\theta_{n}^{\prime}\in\Phi_{K}}\cos(\theta_{n}^{\prime}+\alpha_{n}-\phase{\mu})$,
store $\theta_{n}$, $n=1,2,\ldots,N$

6:Set $g_{0}=h_{0}+\sum_{n=1}^{N}h_{n}e^{j\theta_{n}}$, ${\tt
absgmax}=|g_{0}|$

7:for $l=1,2,\ldots,\frac{M}{K}=N-N^{\prime}$ do

8: Let $(\theta_{n}+\omega\leftarrow\theta_{n})\;\mathrm{mod}\,\Phi_{K}$
$n\in{\cal R}_{m_{l}}$

9: Let  $g_{l}=g_{l-1}+\sum_{n\in{\cal
R}_{m_{l}}}h_{n}\left(e^{j\theta_{n}}-e^{j(\theta_{n}-\omega)\,\mathrm{mod}\,\Phi_{K}}\right)$

10: if $|g_{l}|>{\tt absgmax}$ then

11: Let ${\tt absgmax}=|g_{l}|$

12: Store updated $\theta_{n}$ for $n\in{\cal R}_{m_{l}}$

13: end if

14:end for

15:Read out $\theta_{n}^{*}$ as the stored $\theta_{n}$, $n=1,2,\ldots,N$.

Figure 2: Variation of $|g_{l}|^{2}$ and $|g_{c,l}|^{2}$ with $l$.

##  VI Computational Complexity

We now provide computational complexity figures for our algorithms Algorithm 2
and Algorithm 3 against a number of algorithms from the literature, which are
Closest Point Projection (CPP), Approximation algorithm (APX) [4], and Divide-
and-Sort (DaS) [2]. We will tabulate these results in terms of simulation time
on the same computer (Dell XPS 15 9500 employing Intel Core i9-10885H CPU, 2.4
GHz, with 8 cores and 16 logical processors) with implementations carried out
in Matlab. In Table 4, we have the simulation time results in seconds plotted
against the number of RIS elements $N$ for $K=2$. Then, in Table 4, we have
the simulation results for $K=4$.

Method | $N=200$ | $N=500$ | $N=1000$ | $N=2000$
---|---|---|---|---
DaS | 0.6783 | 4.7491 | 14.4208 | 62.9249
Algorithm 2 | 0.3034 | 2.1840 | 5.1574 | 17.9194
APX | 0.0219 | 0.0633 | 0.0977 | 0.2339
Algorithm 3 | 0.0223 | 0.0765 | 0.1253 | 0.2952
CPP | 0.0033 | 0.0058 | 0.0063 | 0.0104

Table 3: Execution time [s] comparisons for 100 channel realizations in NLOS,
$K=2$.

Method | $N=200$ | $N=500$ | $N=1000$ | $N=2000$
---|---|---|---|---
DaS | 1.2708 | 8.0881 | 26.3001 | 122.4301
Algorithm 2 | 0.7171 | 4.8998 | 15.5991 | 51.6384
APX | 0.0261 | 0.0664 | 0.1130 | 0.2719
Algorithm 3 | 0.0256 | 0.0742 | 0.1336 | 0.3151
CPP | 0.0046 | 0.0070 | 0.0104 | 0.0162

Table 4: Execution time [s] comparisons for 100 channel realizations in NLOS,
$K=4$.

We note that among these algorithms, APX and CPP do not have optimal
performance. DaS is claimed to be optimal, but requires extremely high
execution time. Our algorithms Algorithm 2 and Algorithm 3 are optimal in
terms of performance. Note the substantial reduction in computational
complexity from Algorithm 2 to Algorithm 3.

##  VII Conclusion

In this paper, we provided necessary and sufficient conditions for
determination of optimum phase values in order to maximize the received power
at a UE which receives its transmission by means of reflections from an RIS,
when the phase values are from a discrete-valued set. Algorithms are provided
to achieve this in a number of steps equal to $N$, the number of RIS elements,
or fewer. In the literatute, the number of steps to achieve this maximum is
given as $KN$ or $2N$ on the average, e.g., [1, 2]. In conclusion, for a
discrete-phase RIS, the techniques in this paper achieve the optimum received
power in the smallest number of steps published in the literature with an
elementwise update rule. In addition, in each of those $N$ steps, the
techniques presented determine only one or a small number of phase shifts,
which result in a substantial reduction of computation time, as compared to
the algorithms in the literature, e.g., [2, 3].

Finally, we want to make the following important point. In this paper, we
addressed the ongoing problem in the literature of finding an optimal solution
to the problem (1) within the fewest number of steps, or with minimum
computational complexity. Our Algorithm 3 (or Algorithm 4 when applicable)
achieves this goal. Yet, a secondary result of this paper is that the
intuitive CPP solution, which is based on the independent uniform quantization
of $\theta_{n}$ on the unit circle with $K$ points, results in suboptimal but
very close to the optimal solution with very small complexity for all
practical cases.

##  Appendix A Element-Based Simple Update Rule

We now further simplify Algorithm 2, so that there is no need for calculating
$s_{nk}$ or $\lambda_{l}.$ What we need to have is, given an initial
$\phase{\mu}$ selection, say $e^{j{\phase{\mu}}}\in{\rm
arc}(e^{j\lambda_{i-1}}:e^{j\lambda_{i}}),$ we want to know the $N$-step
update rule ${\cal N}(\lambda_{l}),l=i,i+1,\ldots,i+N-1$ in the for loop of
Algorithm 2.

Claim: Let ${\cal U}$ be the set to define the $N$ consecutive updates in the
for loop of Algorithm 2. For an initial $\phase{\mu_{0}},$ the update rule in
the for loop of Algorithm 2 will be ${\cal
U}=\\{n_{1},n_{2},\ldots,n_{N}|0\leq\varphi_{n_{1}}<\varphi_{n_{2}}<\cdots<\varphi_{n_{N}}<\frac{2\pi}{K},\varphi_{n}=(\alpha_{n}-\phase{\mu_{0}}+\frac{\pi}{K})\;{\rm
mod}\;\frac{2\pi}{K},n=1,2,\ldots,N\\}$.

Proof: Due to space limitations, we will omit the proof from this version of
the paper. We will add it to the final version in the conference proceedings.

It is important to note that, this result relieves Algorithm 2 from the burden
to calculate $NK$ instances of both $s_{nk}$ and $\lambda_{l}$.

## References

  * [1] S. Ren, K. Shen, X. Lin, X. Chen, and Z.-Q. Luo, “A linear time algorithm for the optimal discrete IRS beamforming,” _IEEE Wireless Communications Letters_ , vol. 12, no. 3, pp. 496–500, Mar. 2023.
  * [2] R. Xiong, X. Dong, T. Mi, K. Wan, and R. C. Qiu, “Optimal discrete beamforming of RIS-aided wireless communications: An inner product maximization approach,” arXiv:2211.04167v5, Apr. 2023.
  * [3] J. Sanchez, E. Bengtsson, F. Rusek, J. Flordelis, K. Zhao, and F. Tufvesson, “Optimal, low-complexity beamforming for discrete phase reconfigurable intelligent surfaces,” in _2021 IEEE Global Communications Conference (GLOBECOM)_ , 2021, pp. 01–06.
  * [4] Y. Zhang, K. Shen, S. Ren, X. Li, X. Chen, and Z.-Q. Luo, “Configuring intelligent reflecting surface with performance guarantees: Optimal beamforming,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 5, pp. 967–979, Aug. 2022.