# Proportionate Adaptive Graph Signal Recovery

Razieh Torkamani, Hadi Zayyani, and Mehdi Korki This work was supported by the
Iran National Science Foundation (INSF) (grant number 4005022).R. Torkamani
and H. Zayyani are with the Department of Electrical and Computer Engineering,
Qom University of Technology (QUT), Qom, Iran (e-mail:
<EMAIL_ADDRESS>zayyani@qut.ac.ir).M. Korki is with School of
Science, Computing and Engineering Technologies, Swinburne University of
Technology, Melbourne, Australia (e-mail: mkorki@swin.edu.au).

###### Abstract

This paper generalizes the proportionate-type adaptive algorithm to the graph
signal processing and proposes two proportionate-type adaptive graph signal
recovery algorithms. The gain matrix of the proportionate algorithm leads to
faster convergence than least mean squares (LMS) algorithm. In this paper, the
gain matrix is obtained in a closed-form by minimizing the gradient of the
mean-square deviation (GMSD). The first algorithm is the Proportionate-type
Graph LMS (Pt-GLMS) algorithm which simply uses a gain matrix in the recursion
process of the LMS algorithm and accelerates the convergence of the Pt-GLMS
algorithm compared to the LMS algorithm. The second algorithm is the
Proportionate-type Graph Extended LMS (Pt-GELMS) algorithm, which uses the
previous signal vectors alongside the signal of the current iteration. The Pt-
GELMS algorithm utilizes two gain matrices to control the effect of the signal
of the previous iterations. The stability analyses of the algorithms are also
provided. Simulation results demonstrate the efficacy of the two proposed
proportionate-type LMS algorithms.

###### Index Terms:

Graph signal recovery, Adaptive, Laplacian matrix, least mean-squares,
Proportionate.

##  I Introduction

Graph Signal Processing (GSP) [1]-[3] is a new research paradigm in signal
processing. In GSP, the signal is defined over a graph with irregular domains.
It can better represent the inherent structure of the signals defined over
nodes. The potential application of GSP includes wireless sensor networks,
biological network, social networks, financial networks, and vehicular
networks [3], to name a few. Because of the shift of signal processing to GSP,
the available tools in classical signal processing such as shifting, sampling,
Fourier transform, filters, etc. is generalized to the graph domain. We can
refer to graph sampling, graph signal recovery, graph topology learning, and
graph spectral representation as problems within GSP framework. Graph Signal
Recovery (GSR) is a basic problem in GSP which aims to recover the whole graph
signal by observing signal over only a subset of graph nodes [4]. GSR uses the
inherent relationship between the signal values defined over connecting nodes
in the graph. This relationship is defined by matrices called weighted
adjacency and Laplacian matrix of the graph which will be introduced in
sequel.

There are two types of GSR algorithms. The first type is the non-adaptive GSR
[4]-[11] which often uses some optimization problems to solve GSR problems in
a non-adaptive manner and in a batch-based framework. There is usually a great
deal of complexity involved in these algorithms, especially in the large scale
graph networks. The second type of GSR algorithms are adaptive GSR algorithms
[12]-[14], which similar to adaptive filter counterparts in classical signal
processing, require lower computational complexity and have the potential to
perform well in the time-varying nature of the graph signal and noise. Hence,
in this paper, we focus on the adaptive GSR algorithms.

Adaptive GSR algorithms first appeared in the literature in 2016. As a
pioneering work, a graph Least Mean Square (LMS) algorithm has been developed
in [12], which generalizes the well-known adaptive LMS algorithm to the graph
domain. In sequel, [13] suggested a distributed LMS algorithm for learning the
graph signal over a network. Moreover, an LMS and a Recursive Least Square
(RLS) algorithm have been developed to recover the graph signals from randomly
time-varying subset of nodes [14]. Also, in [15], the LMS algorithm has been
developed for the dynamic graphs in which there is a small perturbation in the
Laplacian matrix. In addition, a joint graph weighted adjacency matrix
learning and graph signal recovery has been suggested using a Kalman filter
for auto-regressive graph signals [16]. To provide scalability and privacy in
networks, an online kernel-based graph-adaptive learning algorithm has been
suggested in [17]. Besides, [18] proposed a distributed adaptive learning of
graph signals using in-network subspace projections. Moreover, a Normalized
LMS (NLMS) graph signal estimation algorithm has been proposed in [19], which
has faster convergence than LMS algorithm and has less computational
complexity than RLS algorithm. Also, [20] proposed two adaptive GSR algorithms
which are Extended LMS (ELMS) algorithm and Fast ELMS (FELMS) algorithm, in
which the signal vectors of previous iterations are reused alongside the
signal available at the current iteration. In addition, a Single-Kernel
Gradraker (SKG) algorithm has been suggested for GSR which uses a Gaussian
kernel and it specifies how to find a suitable variance for the kernel [21].
For adaptive estimation of the graph filter of a graph signal, a graph kernel
RLS algorithm has been developed for a different but related problem [22].

In this paper, we focus on adaptive graph LMS algorithms for GSR. We
generalize the proportionate adaptive filter concept used in adaptive
filtering application [23]-[24], and used in distributed estimation framework
[25]-[26] to the graph domain. Hence, we propose adaptive proportionate GSR
algorithm in which a gain matrix is used in the update of the adaptive
algorithm. When we have a sparse representation of the graph signal in the
domain of Graph Fourier Transforms (GFT), the proposed algorithm can be faster
to converge. For determining the gain matrix, the GMSD criterion is used and a
closed-form formula is obtained for the gain matrix. The stability analyses of
the Graph Proportionate LMS (GPLMS) algorithm are also provided. Finally,
simulation results corroborate the efficacy of the proposed algorithm in
comparison to some state-of-the-art algorithms in the literature.

This paper is organized as follows. Section II presents the essential
background on GSP and problem formulation. Section III presents the proposed
proportionate-type LMS algorithm for GSR. The extended proportionate-type LMS
algorithm is presented in section IV. Section V provides some theoretical
aspects of the proposed algorithms. In Section VI, the simulation results are
discussed. We conclude the paper in in Section VII.

We denote vectors by boldface lowercase letters and matrices by boldface
uppercase letters. The operators $(\cdot)^{T}$ and $(\cdot)^{H}$ represent the
transpose and Hermition transpose, respectively. $E[\cdot]$ represents
expectation. $\text{diag}(\cdot)$ represents a diagonal matrix with its
arguments. The $m$th element of the vector $\mathbf{a}$ is written as $a_{m}$.
The $(m,q)$th element of the matrix $\mathbf{A}$ is written as $a_{mq}$. The
identity matrix of dimension $N$ is written as $\mathbf{I}_{N}$ and the vector
$\mathbf{0}_{N}$ is a length $N$ vector of all zeros. $\mathbf{0}_{M\times N}$
is an $M\times N$ matrix with all entries being zero. The spectral radius of
the square matrix $\mathbf{A}$ is denoted by
$\rho(\mathbf{A})\triangleq\max\\{|\lambda(\mathbf{A})|\\}$.

##  II Problem Formulation and Background

Consider an undirected, connected, weighted graph $G=(\cal V,\cal E)$ with $N$
vertices indexed by $\cal V$$=\\{v_{1},v_{2},...,v_{N}\\}$ and connected
together according to the set of edges $\cal E$. The weighted adjacency matrix
$\mathbf{W}\in\mathbb{R}^{N\times N}$ is the collection of all edge weights
such that $W_{ij}>0$ if $(i,j)\in\cal E$, and $W_{ij}=0$ otherwise. Let
$d_{i}$ be the degree of node $i$ which is defined as
$d_{i}=\sum_{j=1}^{N}{W_{ij}}$. The degree matrix D is a diagonal matrix with
the node degrees as its diagonals. Hence, graph Laplacian matrix is defined as
$\mathbf{L}={\textbf{D}}-\mathbf{W}$.

For undirected graphs, the graph Laplacian is a symmetric and positive semi-
deﬁnite matrix. The eigendecomposition for Laplacian matrix is

$\mathbf{L}=\mathbf{U}\Lambda\mathbf{U}^{H},$ (1)

where $\mathbf{U}\in\mathbb{R}^{N\times N}$ is the matrix containing all the
eigenvectors of $\mathbf{L}$ as its columns, and
$\Lambda\in\mathbb{R}^{N\times N}$ is the diagonal matrix of eigenvalues of
$\mathbf{L}$.

Analogous to the classical signal processing, the graph Fourier transform
(GFT) of a graph signal ${\textbf{x}}\in\mathbb{R}^{N}$ is defined as its
projection onto an orthogonal set of vectors
$\\{{\textbf{u}}_{i}\\}_{i=1,...,N}$, i.e.

${\textbf{s}}=\mathbf{U}^{H}{\textbf{x}}.$ (2)

The basis vector $\\{{\textbf{u}}_{i}\\}$ are usually assumed to be the
eigenvector set of Laplacian matrix [12],[20], or the adjacency matrix
[27]-[28]. Thus, the frequncy-domain representation of the graph signal
conveys the intrinsic information of the graph topology. In this paper, we
follow the definition of GFT basis vector based on the Laplacian matrix, but
the result can be extended to the approach based on the adjacency matrix. The
inverse graph Fourier transform (IGFT) can be defined for reconstruction of
the graph signal x from its frequency domain s as

${\textbf{x}}=\mathbf{U}{\textbf{s}}.$ (3)

In this paper, we utilize the intrinsic sparsity of bandlimited signals and we
assume that the frequency domain representation of the graph signal, i.e. s,
is sparse. Thus, we can use the compressive sensing (CS) theory and assume the
under-sampling of the s as observations. Moreover, we assume that the
observation are sampled noisy versions of the graph signal. Thus, using the
notations ${\textbf{x}}^{o}$ and ${\textbf{s}}^{o}$ for the original graph
signal and its Fourier transform, respectively, we can write the observed
signal at time $n$ as

${\textbf{y}}[n]={\textbf{B}}[n]{\textbf{D}}[n]{\textbf{x}}^{o}+{\textbf{e}}[n]={\textbf{B}}[n]{\textbf{D}}[n]\mathbf{U}{\textbf{s}}^{o}+{\textbf{e}}[n],$
(4)

where B is the CS sensing matrix, and D is the sampling matrix. The above
equation can be rewritten as

${\textbf{y}}[n]={\textbf{A}}[n]{\textbf{s}}^{o}+{\textbf{e}}[n],$ (5)

where ${\textbf{A}}[n]={\textbf{B}}[n]{\textbf{D}}[n]\mathbf{U}$.

##  III Proportionate-type Graph LMS algorithm

The LMS algorithm is one of the most popular algorithms in adaptive filtering,
and is employed in graph signal recovery [12]-[15]. Based on the graph-LMS
approach, the optimal estimation for the original graph signal can be found by
the following recursive procedure

$\displaystyle{\textbf{s}}[n+1]={\textbf{s}}[n]+\mu{\textbf{A}}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(6)

where $\mu$ is the step-size parameter. The proportionate-type LMS (Pt-LMS)
algorithm [23] was introduced as an alternative to the conventional LMS
algorithm, which is proved to converge faster than the LMS algorithm by
assigning a different gain to each coefficient. This gain is proportional to
the magnitude of the coefficient at the current iteration. In this paper, we
propose to use this Pt-LMS algorithm in graph signal recovery problem, which
results in proportionate-type graph LMS (Pt-GLMS) algorithm. The update
equation for the proposed Pt-GLMS algorithm is

${\textbf{s}}[n+1]={\textbf{s}}[n]+\mu\mathbf{G}[n]{\textbf{A}}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(7)

where $\mathbf{G}[n]$ is the gain matrix. This gain matrix, distinguishes the
Pt-LMS from conventional LMS algorithm, and is diagonal, i.e.,
$\mathbf{G}[n]=\mathrm{diag}(g_{1}[n],...,g_{N}[n])$, where $g_{i}[n]$ is the
gain factor of the signal of the $i$’th node at time $n$, and is proportional
to the magnitude of the graph signal at node $i$, i.e.,
$g_{i}[n]\propto|s_{i}[n]|$. If $g_{i}[n]=1,\forall i=1,...,N$, i.e.,
$G[n]={\textbf{I}}_{N\times N}$, the Pt-GLMS and GLMS algorithms are
equivalent. In the literature [29], the weights $g_{i}[n]$ are calculated as

$g_{i}[n]=\frac{\gamma_{i}[n]}{\frac{1}{N}\sum_{j=1}^{N}{\gamma_{j}[N]}},$ (8)

where

$\displaystyle\gamma_{i}[n]=\mathrm{max}\\{\rho\gamma_{min}[n],F[|s_{i}[n]|]\\},$
(9)
$\displaystyle\gamma_{min}[n]=\mathrm{max}\\{\delta,F[|s_{1}[n]|],...,F[|s_{N}[n]|]\\},$
(10)

where $\delta$ is the initialization parameter, and the parameter $\rho$
prevents the inactive coefﬁcients from stalling. For $F[|s_{i}[n]|]$,
different functions have been used in the literature [29].

From (7), the factor $\mu g_{i}[n]$ implies the effective step-size at node
$i$, and indicates that in the case of large magnitude of current graph
signal, i.e., $|s_{i}[n]|$, the value of $g_{i}[n]$ is also large, and, thus,
the effective step-size $\mu g_{i}[n]$ is large, which speed up the
convergence of large coefficients. Conversely, for small magnitude of current
coefficients, the effective step-size is also small. The resulting
proportionate-type graph LMS algorithm is summarized in Algorithm 1.

Algorithm 1 Proposed Proportionate-type Graph LMS (Pt-GLMS) Algorithm

Input Observations y; Fourier basis functions $\mathbf{U}$; Sensing matrix B;
Sampling matrix D; $\mu$; $\rho$; $\delta$; number of time instances $T$ .
Initialize ${\textbf{s}}=0$, $\mathbf{G}=0$, $n=1$.

repeat  • ${\textbf{A}}[n]={\textbf{B}}[n]{\textbf{D}}[n]\mathbf{U}$ • Update
$\\{g_{i}[n]\\}_{n=1}^{N}$ using (8) or (13) •
$\mathbf{G}[n]=\mathrm{diag}(g_{1}[n],...,g_{N}[n])$ •
${\textbf{s}}[n+1]={\textbf{s}}[n]+\mu\mathbf{G}[n]{\textbf{A}}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])$
• $n\longleftarrow n+1$

until $n\leq T$

###  III-A Gain Matrix Calculation for Pt-GLMS

In this subsection, we calculate an optimal gain matrix in order to further
speed up the convergence of the proposed Pt-GLMS algorithm. To this end, we
compute the gradient mean-square deviation (GMSD), and find the optimal gains
by minimizing the GMSD at time $n$. The GMSD is defined as [20]

$\Delta[n]=\mathrm{E}||\tilde{{\textbf{s}}}[n+1]||^{2}_{Q}-\mathrm{E}||\tilde{{\textbf{s}}}[n]||^{2}_{Q},$
(11)

where $\tilde{{\textbf{s}}}[n]={\textbf{s}}^{o}-{\textbf{s}}[n]$ is the error
signal at time $n$, $Q={\textbf{A}}^{H}{\textbf{A}}$, and
$||{\textbf{t}}||^{2}_{P}={\textbf{t}}^{H}{\textbf{P}}{\textbf{t}}$. The GMSD
of the $i^{\prime}$th node at time $n$ can be written as

$\Delta_{i}[n]=\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}_{i}[n+1]||^{2}_{Q}\Big{]}-\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}_{i}[n]||^{2}_{Q}\Big{]}.$
(12)

The optimum gain for node $i$ at time instant $n$ (i.e., $g_{i}[n]$), by
setting the derivative of $\Delta_{i}[n]$ with respect to $g_{i}[n]$ equal to
zero, is written as

$g_{i}[n]=\frac{\mu\Big{[}\hat{{\textbf{e}}}^{T}[n]{\textbf{A}}_{i}[n]{\textbf{A}}^{T}_{i}[n]\hat{{\textbf{e}}}[n]-{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}^{T}_{i}[n]\Big{]}}{m^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]},$
(13)

where

$\hat{{\textbf{e}}}[n]={\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n],$ (14)

and

$m_{i}[n]=\mu{\textbf{A}}_{i}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]).$
(15)

###### Proof.

See Appendix A. ∎

Computing the values of $g_{i}[n]$ using (13)-(15), and substituting
$\mathbf{G}[n]=\mathrm{diag}(g_{1}[n],...,g_{N}[n])$, the graph signal can be
recovered by (7), which yields the minimum GMSD.

As an example, consider a graph with $N=50$ nodes. The spectral content of the
graph signal is limited to the first $15$ eigenvectors of the graph Laplacian
matrix. The observation noise is drawn from a zero-mean Gaussian distribution
with a diagonal covariance matrix. Fig. 1 shows the normalized MSD (NMSD),
which is calculated as follows

$\mathrm{NMSD}[n]=\frac{||{\textbf{s}}^{o}-{\textbf{s}}[n]||^{2}}{||{\textbf{s}}^{o}||^{2}}.$
(16)

As can be seen, using proportionate-type algorithm for the reconstruction of
graph signals increases the convergence rate of the algorithm. Moreover, the
GMSD of the Pt-GLMS using the proposed gain matrix is smaller than that of the
others.

Figure 1: Transient NMSD versus iteration index when estimating a synthetic
data with $N=50$ nodes.

##  IV Proportionate-type Graph Extended LMS Algorithm

In this subsection, we propose an algorithm to enhance the estimation
performance of the Pt-GLMS algorithm proposed in previous section. The
proposed algorithm is based on the extended LMS (ELMS) algorithm proposed in
[20] and speeds up the convergence of the Pt-GLMS algorithm by using the
observations of previous times as well as the current observation. In this
work, we propose the proportionate-type of the ELMS algorithm in GSR
framework, which assigns the gain matrix to the graph signal estimation
process. The resulting algorithm is proportionate-type graph ELMS (Pt-GELMS)
algorithm. The update equation for the proposed Pt-GELMS algorithm is

$\displaystyle{\textbf{s}}[n+1]=$
$\displaystyle{\textbf{s}}[n]+\mu\mathbf{G}[n]{\textbf{A}}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])$
$\displaystyle+\mu\mathbf{H}[n]\sum_{j=1}^{K-1}{\textbf{A}}^{T}[n-j]\Big{(}{\textbf{y}}[n-j]-{\textbf{A}}[n-j]{\textbf{s}}[n]\Big{)},$
(17)

where $K$ is the number of previous time instants used for estimating the
graph signal at the current iteration, $\mathbf{G}[n]$ and $\mathbf{H}[n]$ are
diagonal gain matrices assigned to the current and previous observations,
respectively. Similar to the previous sections, we have
$\mathbf{G}[n]=\mathrm{diag}(g_{1}[n],...,g_{N}[n])$ and
$\mathbf{H}[n]=\mathrm{diag}(h_{1}[n],...,h_{N}[n])$. Obviously, when
$\mathbf{H}[n]={\bf{0}}_{N\times N}$, the Pt-GELMS algorithm is equivalent to
the Pt-GLMS algorithm, and when $\mathbf{G}[n]=\mathbf{H}[n]={\bf{0}}_{N\times
N}$, the Pt-GELMS algorithm reduces to the GLMS algorithm. The resulting
proportionate-type graph extended LMS algorithm is summarized in Algorithm 2.

###  IV-A Gain Matrix Calculation for Pt-GELMS

In this subsection, similar to the process in section III-A, we find an
optimal value for the gain matrices by computing and minimizing the GMSD at
time $n$. To calculate the GMSD of the $i^{\prime}$th node at time $n$, the
update equation (IV) for the $i^{\prime}$th node can be written as

$s_{i}[n+1]=s_{i}[n]+g_{i}[n]m_{1}[n]+h_{i}[n]m_{2}[n],$ (18)

where

$\displaystyle
m_{1}[n]=\mu{\textbf{A}}^{T}_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(19) $\displaystyle
m_{2}[n]=\mu\sum_{j=1}^{K-1}{\textbf{A}}^{T}_{i}[n-j]({\textbf{y}}[n-j]-{\textbf{A}}[n-j]{\textbf{s}}[n]).$
(20)

Similar to section III-A, computing the GMSD of the $i^{\prime}$th node at
time $n$, i.e. $\Delta_{i}[n]$, and setting
$\frac{\partial\Delta_{i}[n]}{\partial g_{i}[n]}=0$ and
$\frac{\partial\Delta_{i}[n]}{\partial h_{i}[n]}=0$ yields

$g_{i}[n]=\frac{r_{3}[n]}{r_{4}[n]},$ (21)

and

$h_{i}[n]=\frac{r_{5}[n]}{r_{6}[n]},$ (22)

where

$\displaystyle r_{3}[n]=$ $\displaystyle-
h_{i}[n]m_{1}[n]m_{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}-\mu{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}_{i}^{T}[n]$
$\displaystyle+\mu({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]{\textbf{A}}_{i}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(23) $\displaystyle r_{4}[n]=$ $\displaystyle
m_{1}^{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]},$ (24) $\displaystyle r_{5}[n]=$
$\displaystyle-g_{i}[n]m_{1}[n]m_{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}$
$\displaystyle-\mu{\textbf{A}}_{i}[n]{\textbf{C}}_{e}\sum_{j=1}^{K-1}{{\textbf{A}}_{i}^{T}[n-j]}+\mu({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}$
$\displaystyle\times{\textbf{A}}_{i}[n]{\textbf{A}}_{i}^{T}[n]\sum_{j=1}^{K-1}{({\textbf{y}}[n-j]-{\textbf{A}}[n-j]{\textbf{s}}[n])},$
(25) $\displaystyle r_{6}[n]=$ $\displaystyle
m_{2}^{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}.$ (26)

###### Proof.

See Appendix B. ∎

Obviously, $g_{i}[n]$ and $h_{i}[n]$ are interdependent, and, thus, should be
updated in a repetitive manner.

Algorithm 2 Proposed Proportionate-type Graph extended LMS (Pt-GELMS)
Algorithm

Input Observations y; Fourier basis functions $\mathbf{U}$; Sensing matrix B;
Sampling matrix D; $\mu$; $K$; number of time instances $T$ .
Initialize ${\textbf{s}}=0$, $\mathbf{G}=0$, $\mathbf{H}=0$, $n=1$.

repeat  • ${\textbf{A}}[n]={\textbf{B}}[n]{\textbf{D}}[n]\mathbf{U}$ • Update
$\\{g_{i}[n]\\}_{n=1}^{N}$ and $\\{h_{i}[n]\\}_{n=1}^{N}$ using (21) and (22),
respectively • $\mathbf{G}[n]=\mathrm{diag}(g_{1}[n],...,g_{N}[n])$ and
$\mathbf{H}[n]=\mathrm{diag}(h_{1}[n],...,h_{N}[n])$ •
${\textbf{s}}[n+1]={\textbf{s}}[n]+\mu\mathbf{G}[n]{\textbf{A}}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])$
$+\mu\mathbf{H}[n]\sum_{j=1}^{K-1}{\textbf{A}}^{T}[n-j]\Big{(}{\textbf{y}}[n-j]-{\textbf{A}}[n-j]{\textbf{s}}[n]\Big{)}$
• $n\longleftarrow n+1$

until $n\leq T$

##  V Theoretical Analysis

In this section, some theoretical analysis of the proposed algorithms are
provided. Since the Pt-GLMS is a special case of Pt-GELMS, we only discuss the
theoretical analysis of Pt-GELMS, and the analysis for the Pt-GLMS can be
achieved by setting $\mathbf{H}[n]={\bf{0}}_{N\times N}$.

###  V-A Mean-Square Analysis

In this subsection, we study the mean-square behavior of the proposed
algorithms. To compute MSD, note that using the definition
$\tilde{{\textbf{s}}}[n]={\textbf{s}}^{o}-{\textbf{s}}[n]$, and subtracting
both sides of (IV) from ${\textbf{s}}^{o}$, we obtain

$\tilde{{\textbf{s}}}[n+1]=\Big{[}{\textbf{I}}-\mu{\textbf{B}}_{1}[n]\Big{]}\tilde{{\textbf{s}}}[n]+\mu{\textbf{B}}_{2}[n]+\mu{\textbf{B}}_{3}[n],$
(27)

where

$\displaystyle{\textbf{B}}_{1}[n]=\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{A}}[n]+\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{A}}[n-j]},$
(28)
$\displaystyle{\textbf{B}}_{2}[n]=\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{e}}[n],$
(29)
$\displaystyle{\textbf{B}}_{3}[n]=\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{e}}[n-j]}.$
(30)

######  Proposition 1.

Assume that the noise samples at different times are independent. Then, for
any bounded initial condition, the proposed Pt-GELMS algorithm asymptotically
converges in the mean-square error sense if

$0<\mu<\frac{2}{\lambda_{max}\left({\textbf{B}}_{1}\right)},$ (31)

where $\lambda_{max}(.)$ denotes the maximum eigenvalue of the matrix therein.

###### Proof.

See Appendix C. ∎

###  V-B Mean Performance

In this subsection, we study the steady-state mean performance of the proposed
Pt-GELMS algorithm. Taking the expectation of both sides of (27) yields

$\mathrm{E}\Big{[}\tilde{{\textbf{s}}}[n+1]\Big{]}=\Big{[}{\textbf{I}}-\mu{\textbf{B}}_{1}[n]\Big{]}\mathrm{E}\Big{[}\tilde{{\textbf{s}}}[n]\Big{]}.$
(32)

As the algorithm reaches the steady-state, we have
$\mathrm{E}\Big{[}\tilde{{\textbf{s}}}[n+1]\Big{]}\longrightarrow\mathrm{E}\Big{[}\tilde{{\textbf{s}}}[n]\Big{]}$.
Thus, to guarantee the convergence in the mean sense, the condition
$\rho({\textbf{I}}-\mu{\textbf{B}}_{1}[n])<1$ should be satisfied, which
implies

$0<\mu<\frac{2}{\lambda_{max}\left({\textbf{B}}_{1}\right)}.$ (33)

###  V-C Steady-State Performance

Taking the limit of (C) as $n\to\infty$ yields

$\lim_{n\rightarrow\infty}\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n]||^{2}_{({\textbf{I}}-{\textbf{Q}})\boldsymbol{\phi}}\Big{]}=\mu^{2}\mathrm{vec}({\textbf{P}})^{T}\boldsymbol{\phi}.$

To evaluate the steady-state MSD of the proposed Pt-GELMS algorithm (see
(IV)), we can set
$\phi=({\textbf{I}}-{\textbf{Q}})^{-1}\mathrm{vec}({\textbf{I}})$ and obtain

$\mathrm{MSD}=\lim_{n\rightarrow\infty}\mathrm{E}||\tilde{{\textbf{x}}}[n]||^{2}=\lim_{n\rightarrow\infty}\mathrm{E}||\tilde{{\textbf{s}}}[n]||^{2}$
$=\mu^{2}\mathrm{vec}({\textbf{P}})^{T}({\textbf{I}}-{\textbf{Q}})^{-1}\mathrm{vec}({\textbf{I}}).$
(34)

##  VI Simulation Results

Figure 2: Adaptive algorithms performance versus iteration index for $K=8$,
$M=30$, and $|S|=20$. Figure 3: Transient NMSD versus iteration index for
$M=30$, $|S|=20$, and different values of $K$.

In this section, we evaluate the performance of the proposed algorithms via
some numerical results on synthetic and real data.

In the first scenario, we consider a synthetic graph with $N=50$ nodes.
Similar to [30], the edge weights are drawn randomly from a uniform
distribution $W_{ij}\sim\mathcal{U}(0,1)$, and the weighted adjacency matrix
is derived as $\mathbf{W}=(\mathbf{W}+\mathbf{W}^{T})/2$. The graph Laplacian
matrix is calculated as $\mathbf{L}={\textbf{D}}-\mathbf{W}$, where D is the
diagonal degree matrix with the diagonal elements as
$d_{i}=\sum_{j=1}^{N}{W_{ij}}$. The spectral content of the graph is limited
to the first $15$ eigenvectors of the graph Laplacian matrix. The observation
noise in (5) is drawn from a zero-mean Gaussian distribution with a diagonal
covariance matrix ${\textbf{C}}_{e}=\sigma_{e}^{2}{\textbf{I}}$ with
$\sigma_{e}^{2}=0.01$. We compare the results of our proposed algorithms with
that of the graph LMS algorithm [12], and extended proportionate-type LMS
[20]. Fig. 2 shows the transient behavior of the NMSD in (16) versus the
iteration index for the number of previous time instants $K=8$, number of CS
measurements $M=30$, number of samples $|S|=20$, and step-size $\mu=0.01$.
Each point in the curves is the result of ensemble average over 50 independent
simulations. It can be seen that the proposed proportionate-type LMS
algorithms converge faster than the conventional LMS. Moreover, using the
proposed algorithm for estimation of gain matrices lead to a further increase
in the convergence rate.

(a) $M=10$

(b) $M=30$

(c) $M=40$

Figure 4: Transient NMSD behaviour for $K=8$, $|S|=20$, and different values
of $M$.

In Fig. 3, the transient behaviour of the NMSD versus the iteration index are
shown for $M=30$, step-size $\mu=0.01$, and for different values of $K$. The
curves are averaged over $50$ independent trials. It can be seen that
increasing the value of the parameter $K$, which incorporates further previous
signals in estimating the current signal, results in higher convergence speed
of the algorithms.

Fig. 4 depicts the transient behaviour of the NMSD for LMS, FELMS, Pt-GLMS
(with the proposed gain matrix), and Pt-GELMS algorithms, considering
different number of CS measurements, i.e., different values of $M$, and $K=8$,
$|S|=20$, and the step-size $\mu=0.01$. The results are averaged over 50
independent trials. As can bee seen from the curves, increasing the number of
CS measurements, which means increasing the number of linear combination of
signal ensembles in hand, leads to an increase in convergence rate and
decrease in NMSD. Moreover, with the same number of measurements, the proposed
Pt-GLMS and Pt-GELMS algorithms outperform their peer algorithms, i.e., LMS
and FELMS, respectively.

In Fig. 5, the transient behaviour of the NMSD for LMS, FELMS, Pt-GLMS (with
the proposed gain matrix), and Pt-GELMS algorithms are presented, considering
three different values of bandwidth, $|F|$, and $K=8$, $|S|=20$, $M=30$, and
the step-size $\mu=0.01$. The results are averaged over 50 independent trials.
As can bee seen from the curves, increasing the bandwidth of the signal
spectrum, leads to an increase in convergence rate and decrease in NMSD.

(a) $|F|=10$

(b) $|F|=15$

(c) $|F|=20$

Figure 5: Transient NMSD of the adaptive algorithms for $K=8$, $|S|=20$,
$M=30$, and different values of $|F|$.

(a) $|S|=10$

(b) $|S|=20$

(c) $|S|=30$

Figure 6: Transient NMSD of the adaptive algorithms versus iteration index for
$K=8$, $|F|=15$, $M=30$, and different values of $|S|$.

In Fig. 6, the transient behaviour of the NMSD for LMS, FELMS, Pt-GLMS (with
the proposed gain matrix), and Pt-GELMS algorithms are presented, considering
three different number of samples selected, $|S|$, and $K=8$, $|F|=15$,
$M=30$, and the step-size $\mu=0.01$. The results are averaged over 50
independent trials. As can bee seen from the curves, increasing number of
nodes selected in the sampling procedure, leads to an increase in convergence
rate and decrease in NMSD. As seen, both the proposed Pt-GELMS and Pt-GLMS
algorithms outperforms their non-proportionate type peers, i.e., FELMS and
LMS. Furthermore, in all cases, the Pt-GLMS with the proposed gain matrix
perform better than the Pt-GLMS with conventional gain matrix used in the
literature.

In the second scenario, for assessing the performance of the adaptive
algorithms with real-world data, we consider the temperature graph signal. In
this case, we use the dataset downloaded from the Intel Berkeley Research lab
(refer to [31], [10], [11]), in which the temperature values from a total of
$N=54$ sensors are acquired. We aim to estimate the temperature values using
the competing adaptive methods. Similar to the first case, we add an
observation noise which is a zero-mean Gaussian noise with a diagonal
covariance matrix ${\textbf{C}}_{e}=\sigma_{e}^{2}{\textbf{I}}$ with
$\sigma_{e}^{2}=3$ (which yields $SNR\sim 25dB$). The transient NMSD versus
the iteration index is depicted in Fig. 7 for three different numbers of
selected nodes ($M=30,40,54$), where we have used $K=6$. The figures show the
superiority of the proposed algorithm in comparison to the other adaptive
algorithms.

(a) $M=30$

(b) $M=40$

(c) $M=N(=54)$

Figure 7: Adaptive algorithms performance when estimating the temperature data
from [31].

##  VII Conclusion

In this paper, we formulated the proportionate adaptive graph signal recovery
algorithm. Hence, the proportionate adaptive filtering algorithm in the
classical signal processing is generalized to proportionate adaptive GSR
algorithm in GSP. Two proportionate-type algorithms are proposed for adaptive
GSR, of which the first is the Pt-GLMS algorithm and the second is the Pt-
GELMS algorithm. In two algorithms, the gain matrix (or matrices) are obtained
optimally in a closed form via minimizing the GMSD. Some theoretical analysis
of the proposed algorithms such as mean-square convergence analysis, mean
performance analysis, and steady-state performance analysis are also
presented. Simulation results in both synthetic and real data, demonstrate the
faster convergence of the proposed proportionate-type algorithms in comparison
to non-proportionate-type counterparts.

##  Appendix A Proof of Gain Matrix for Pt-GLMS

The GMSD of the $i^{\prime}$th node at time $n$ can be written as

$\Delta_{i}[n]=\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}_{i}[n+1]||^{2}_{Q}\Big{]}-\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}_{i}[n]||^{2}_{Q}\Big{]}.$
(35)

Moreover, from (7), we have

$s_{i}[n+1]=s_{i}[n]+\mu
g_{i}[n]{\textbf{A}}_{i}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(36)

where ${\textbf{A}}_{i}$ is the $i^{\prime}$th column of matrix A. We rewrite
the above equation in terms of the signal errors as

$\tilde{s}_{i}[n+1]=\tilde{s}_{i}[n]-g_{i}[n]m_{i}[n],$ (37)

where

$m_{i}[n]=\mu{\textbf{A}}_{i}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]).$
(38)

Thus, we have

$||\tilde{s}_{i}[n+1]||^{2}_{Q}=||\tilde{s}_{i}[n]-g_{i}[n]m_{i}[n]||^{2}_{Q}$
$\displaystyle=$
$\displaystyle\Big{[}(\tilde{s}_{i}[n]-g_{i}[n]m_{i}[n])^{2}\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}$
$\displaystyle+\>2\Big{[}(\tilde{s}_{i}[n]-g_{i}[n]m_{i}[n])\sum_{j=1}^{N}\sum^{N}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}$
$\displaystyle=$
$\displaystyle\Big{[}\tilde{s}^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}+\Big{[}g^{2}_{i}[n]m^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-\>2\Big{[}\tilde{s}_{i}[n]g_{i}[n]m_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}$
$\displaystyle+\>2\Big{[}\tilde{s}_{i}[n]\sum_{j=1}^{N}\sum^{N}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}$
$\displaystyle-\>2\Big{[}g_{i}[n]m_{i}[n]\sum_{j=1}^{N}\sum^{N}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}.$

It can be shown that the sum of the first and fourth terms in the above
equation is $||\tilde{s}_{i}[n]||^{2}_{Q}$. Moreover, the sum of the third and
fifth terms is
$-2g_{i}[n]m_{i}[n]\tilde{s}^{T}[n]{\textbf{A}}^{T}[n]{\textbf{A}}_{i}[n]$.
Thus, we have

$\Delta_{i}[n]=f_{1}[n]-2f_{2}[n],$ (39)

where

$\displaystyle
f_{1}[n]=\mathrm{E}\Big{[}g^{2}_{i}[n]m^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}$
(40) $\displaystyle
f_{2}[n]=\mathrm{E}\Big{[}g_{i}[n]m_{i}[n]\tilde{s}^{T}[n]{\textbf{A}}^{T}[n]{\textbf{A}}_{i}[n]\Big{]}.$
(41)

It can be seen that we can calculate the GMSD if $f_{1}[n]$ and $f_{2}[n]$ are
given. On the other hand, since the signal error vector $\tilde{{\textbf{s}}}$
is not available, we can not compute $f_{2}[n]$ directly. To solve this
problem, we rewrite $f_{2}[n]$ as

$\displaystyle f_{2}[n]$
$\displaystyle=\mathrm{E}\Big{[}g_{i}[n]m_{i}[n]({\textbf{A}}[n]\tilde{{\textbf{s}}}[n])^{T}{\textbf{A}}_{i}[n]\Big{]}$
$\displaystyle=\mathrm{E}\Big{[}g_{i}[n]m_{i}[n]({\textbf{A}}[n]({\textbf{s}}^{o}-{\textbf{s}}[n]))^{T}{\textbf{A}}_{i}[n]\Big{]}$
$\displaystyle=\mathrm{E}\Big{[}g_{i}[n]m_{i}[n]({\textbf{y}}[n]-{\textbf{e}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]\Big{]}.$
(42)

Assuming that the noise ${\textbf{e}}[n]$ is independent of the other
processes, we have

$\displaystyle f_{2}[n]=$
$\displaystyle\mu\mathrm{E}\Big{[}g_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]{\textbf{A}}^{T}_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])\Big{]}$
$\displaystyle-\mu\mathrm{E}\Big{[}g_{i}[n]{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}^{T}_{i}[n]\Big{]}.$
(43)

Thus, the GMSD becomes as

$\displaystyle\Delta_{i}[n]=\mathrm{E}\Big{[}g^{2}_{i}[n]m^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-2\mu\mathrm{E}\Big{[}g_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]{\textbf{A}}^{T}_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])\Big{]}$
$\displaystyle+2\mu\mathrm{E}\Big{[}g_{i}[n]{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}^{T}_{i}[n]\Big{]}.$
(44)

The optimum gain for node $i$ at time instance $n$ (i.e., $g_{i}[n]$) can be
found by setting the derivative of $\Delta_{i}[n]$ with respect to $g_{i}[n]$
to zero and obtain

$g_{i}[n]=\frac{\mu\Big{[}\hat{{\textbf{e}}}^{T}[n]{\textbf{A}}_{i}[n]{\textbf{A}}^{T}_{i}[n]\hat{{\textbf{e}}}[n]-{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}^{T}_{i}[n]\Big{]}}{m^{2}_{i}[n]\sum_{j=1}^{N}a^{2}_{ji}[n]},$
(45)

where $\hat{{\textbf{e}}}[n]={\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]$.

##  Appendix B Proof of Gain Matrix for Pt-GELMS

From (18), the signal errors are

$\tilde{s}_{i}[n+1]=\tilde{s}_{i}[n]-g_{i}[n]m_{1}[n]-h_{i}[n]m_{2}[n],$ (46)

where $m_{1}[n]$ and $m_{2}[n]$ are given by (19). Thus, we have

$\displaystyle\mathrm{E}||\tilde{s}_{i}[n+1]$
$\displaystyle||^{2}_{Q}=\mathrm{E}||\tilde{s}_{i}[n]-g_{i}[n]m_{1}[n]-h_{i}[n]m_{2}[n]||^{2}_{Q}$
$\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}\tilde{s}^{2}_{i}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}+\mathrm{E}\Big{[}g^{2}_{i}[n]m^{2}_{1}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle+\mathrm{E}\Big{[}h^{2}_{i}[n]m^{2}_{2}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle+2\mathrm{E}\Big{[}g_{i}[n]h_{i}[n]m_{1}[n]m_{2}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-2\mathrm{E}\Big{[}\tilde{s}_{i}[n]g_{i}[n]m_{1}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-2\mathrm{E}\Big{[}\tilde{s}_{i}[n]h_{i}[n]m_{2}[n]\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle+2\mathrm{E}\Big{[}\tilde{s}_{i}[n]\sum_{j=1}^{K-1}\sum^{K-1}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}$
$\displaystyle-2\mathrm{E}\Big{[}g_{i}[n]m_{1}[n]\sum_{j=1}^{K-1}\sum^{K-1}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}$
$\displaystyle-2\mathrm{E}\Big{[}h_{i}[n]m_{2}[n]\sum_{j=1}^{K-1}\sum^{K-1}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}.$

Therefore, the GMSD becomes as

$\displaystyle\Delta_{i}[n]$ $\displaystyle=$
$\displaystyle\mathrm{E}||\tilde{s}_{i}[n+1]||^{2}_{Q}-\mathrm{E}||\tilde{s}_{i}[n]||^{2}_{Q}$
$\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}(g^{2}_{i}[n]m^{2}_{1}[n]+2g_{i}[n]h_{i}[n]m_{1}[n]m_{2}[n]$
$\displaystyle+\>h^{2}_{i}[n]m^{2}_{2}[n])\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-\>2\mathrm{E}\Big{[}\tilde{s}_{i}[n](g_{i}[n]m_{1}[n]+h_{i}[n]m_{2}[n])\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]}$
$\displaystyle-\>2\mathrm{E}\Big{[}(g_{i}[n]m_{1}[n]+h_{i}[n]m_{2}[n])$
$\displaystyle\times\>\sum_{j=1}^{K-1}\sum^{K-1}_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}\tilde{s}_{k}[n]a_{ji}[n]a_{jk}[n]\Big{]}$
$\displaystyle=$ $\displaystyle r_{1}[n]-2r_{2}[n],$

where

$\displaystyle r_{1}[n]=$
$\displaystyle\mathrm{E}\Big{[}(g^{2}_{i}[n]m^{2}_{1}[n]+2g_{i}[n]h_{i}[n]m_{1}[n]m_{2}[n]$
$\displaystyle+h^{2}_{i}[n]m^{2}_{2}[n])\sum_{j=1}^{K-1}a^{2}_{ji}[n]\Big{]},$
$\displaystyle r_{2}[n]=$
$\displaystyle\mathrm{E}\Big{[}(g_{i}[n]m_{1}[n]+h_{i}[n]m_{2}[n])\tilde{{\textbf{s}}}^{T}[n]{\textbf{A}}^{T}[n]{\textbf{A}}_{i}[n]\Big{]}.$

Similar to section III-A, the signal error vector $\tilde{{\textbf{s}}}[n]$ is
not available, and thus, $r_{2}[n]$ is not computable. Thus, we solve this
problem using the following relations

${\textbf{A}}\tilde{{\textbf{s}}}[n]={\textbf{A}}({\textbf{s}}^{o}-{\textbf{s}}[n])={\textbf{y}}[n]-{\textbf{e}}[n]-{\textbf{A}}[n]{\textbf{s}}[n],$
(48)

which yields

$\displaystyle r_{2}[n]=$
$\displaystyle\mathrm{E}\Big{[}({\textbf{A}}[n]\tilde{{\textbf{s}}}[n])^{T}{\textbf{A}}_{i}[n](g_{i}[n]m_{1}[n]+h_{i}[n]m_{2}[n])\Big{]}$
$\displaystyle=$
$\displaystyle\mu\mathrm{E}\Big{[}g_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]{\textbf{A}}^{T}_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])\Big{]}$
$\displaystyle+\mu\mathrm{E}\Big{[}h_{i}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]\sum_{j=1}^{K-1}{\textbf{A}}^{T}_{i}[n-j]\Big{]}$
$\displaystyle-\mu\mathrm{E}\Big{[}g_{i}[n]{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}^{T}_{i}[n]\Big{]}$
$\displaystyle-\mu\mathrm{E}\Big{[}h_{i}[n]{\textbf{A}}_{i}[n]{\textbf{C}}_{e}\sum_{j=1}^{K-1}{\textbf{A}}^{T}_{i}[n-j]\Big{]}.$
(49)

Substituting the above equation in (B) and setting
$\frac{\partial\Delta_{i}[n]}{\partial g_{i}[n]}=0$ and
$\frac{\partial\Delta_{i}[n]}{\partial h_{i}[n]}=0$ yields

$g_{i}[n]=\frac{r_{3}[n]}{r_{4}[n]},$ (50)

and

$h_{i}[n]=\frac{r_{5}[n]}{r_{6}[n]},$ (51)

where

$\displaystyle r_{3}[n]=$ $\displaystyle-
h_{i}[n]m_{1}[n]m_{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}-\mu{\textbf{A}}_{i}[n]{\textbf{C}}_{e}{\textbf{A}}_{i}^{T}[n]$
$\displaystyle+\mu({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}{\textbf{A}}_{i}[n]{\textbf{A}}_{i}^{T}[n]({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n]),$
(52) $\displaystyle r_{4}[n]=$ $\displaystyle
m_{1}^{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]},$ (53) $\displaystyle r_{5}[n]=$
$\displaystyle-g_{i}[n]m_{1}[n]m_{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}$
$\displaystyle-\mu{\textbf{A}}_{i}[n]{\textbf{C}}_{e}\sum_{j=1}^{K-1}{{\textbf{A}}_{i}^{T}[n-j]}+\mu({\textbf{y}}[n]-{\textbf{A}}[n]{\textbf{s}}[n])^{T}$
$\displaystyle\times{\textbf{A}}_{i}[n]{\textbf{A}}_{i}^{T}[n]\sum_{j=1}^{K-1}{({\textbf{y}}[n-j]-{\textbf{A}}[n-j]{\textbf{s}}[n])},$
(54) $\displaystyle r_{6}[n]=$ $\displaystyle
m_{2}^{2}[n]\sum_{j=1}^{K-1}{a^{2}_{ji}[n]}.$ (55)

##  Appendix C Proof of Proposition 1

.

Let $\boldsymbol{\Phi}\in\mathbb{R}^{N\times N}$ be an arbitrary matrix. Then,
$\boldsymbol{\Phi}-$weighted norm of both sides of (27) yields

$\displaystyle\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n+1]||$
${}^{2}_{\boldsymbol{\Phi}}\Big{]}=\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n]||^{2}_{\boldsymbol{\Phi}^{\prime}}\Big{]}$
$\displaystyle+\mu^{2}\mathrm{E}\Big{[}{\textbf{e}}^{T}[n]{\textbf{A}}[n]\mathbf{G}^{T}[n]\boldsymbol{\Phi}\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{e}}^{T}[n]\Big{]}$
$\displaystyle+\mu^{2}\mathrm{E}\Big{[}\sum_{j=1}^{K-1}{{\textbf{e}}^{T}[n-j]{\textbf{A}}[n-j]}\mathbf{H}^{T}[n]\boldsymbol{\Phi}$
$\displaystyle\times\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{e}}[n-j]}\Big{]},$
(56)

where

$\boldsymbol{\Phi}^{\prime}=({\textbf{I}}-\mu{\textbf{B}}_{1}[n])^{T}\boldsymbol{\Phi}({\textbf{I}}-\mu{\textbf{B}}_{1}[n]).$
(57)

and

${\textbf{B}}_{1}[n]=\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{A}}[n]+\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{A}}[n-j]},\\\
$

Let $\boldsymbol{\phi}=\mathrm{vec}(\boldsymbol{\Phi})$, where the operator
$\mathrm{vec}(.)$ aggregates the columns of the matrix therein on top of each
other; then, we have

$\displaystyle\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n+1]$
$\displaystyle||^{2}_{\boldsymbol{\phi}}\Big{]}=\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n]||^{2}_{{\textbf{Q}}\boldsymbol{\phi}}\Big{]}$
$\displaystyle+\mu^{2}\mathrm{Tr}\Big{[}\boldsymbol{\Phi}\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{C}}_{e}{\textbf{A}}[n]\mathbf{G}^{T}[n]\Big{]}$
$\displaystyle+\mu^{2}\mathrm{Tr}\Big{[}\boldsymbol{\Phi}\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{C}}_{e}{\textbf{A}}[n-j]}\mathbf{H}^{T}[n]\Big{]},$
(58)

where $\mathrm{Tr}(.)$ is the trace operator, and we have

${\textbf{Q}}=({\textbf{I}}-\mu{\textbf{B}}_{1}[n])\otimes({\textbf{I}}-\mu{\textbf{B}}_{1}[n]).$
(59)

In (C), we have used the independency assumption among noise samples at
different times. Exploiting the trace property
$\mathrm{Tr}(\boldsymbol{\Phi}{\textbf{X}})=\mathrm{vec}({\textbf{X}}^{H})^{T}\mathrm{vec}(\boldsymbol{\Phi})$
in the above equation, we obtain

$\displaystyle\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n+1]$
$\displaystyle||^{2}_{\boldsymbol{\phi}}\Big{]}=\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n]||^{2}_{{\textbf{Q}}\boldsymbol{\phi}}\Big{]}$
$\displaystyle+\mu^{2}\mathrm{vec}\Big{[}\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{C}}_{e}{\textbf{A}}[n]\mathbf{G}^{T}[n]\Big{]}\boldsymbol{\phi}$
$\displaystyle+\mu^{2}\mathrm{vec}\Big{[}\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{C}}_{e}{\textbf{A}}[n-j]}\mathbf{H}[n]\Big{]}\boldsymbol{\phi}$
$\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}||\tilde{{\textbf{s}}}[n]||^{2}_{{\textbf{Q}}\boldsymbol{\phi}}\Big{]}+\mu^{2}\mathrm{vec}({\textbf{P}})^{T}\boldsymbol{\phi},$
(60)

where

$\displaystyle{\textbf{P}}=$
$\displaystyle\mathbf{G}[n]{\textbf{A}}^{T}[n]{\textbf{C}}_{e}{\textbf{A}}[n]\mathbf{G}^{T}[n]$
$\displaystyle+\mathbf{H}[n]\sum_{j=1}^{K-1}{{\textbf{A}}^{T}[n-j]{\textbf{C}}_{e}{\textbf{A}}[n-j]}\mathbf{H}[n].$
(61)

It is proved in [12] that the above convergence satisfies if

$0<\mu<\frac{2}{\lambda_{max}\left({\textbf{B}}_{1}\right)}$ (62)

where $\lambda_{max}(.)$ denotes the maximum eigenvalue of the matrix therein.

## References

  * [1] A. Sandryhaila, and J. M. F. Moura,  “Discrete signal processing on Graphs,”  IEEE Transactions on Signal Processing, vol. 61, no. 7, pp. 1644–1656, Apr. 2013.
  * [2] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,  “The emerging field of signal processing on Graphs: Extending high-dimensional data analysis to networks and other irregular domains,”  IEEE Trans. Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, May 2013.
  * [3] A. Ortega, P. Frossard, J. Kovacevic, J. F. Moura, and P. Vandergheynst,  “Graph Signal Processing: Overview, Challanges, and Applications,”  Proceedings of the IEEE, vol. 106, pp. 808–828, 2018.
  * [4] S. Chen, A. Sandryhaila, J. M. F. Moura, and J. Kovacevic,  “Signal Recovery on Graphs: Variation Minimization,”  IEEE Transaction on Signal Processing., vol. 63, no. 17, pp. 4609–4624, Sep 2015.
  * [5] D. Romero, M. Ma, and G. B. Giannakis,  “Kernel-based Reconstruction of Graph Signals,”  IEEE Transaction on Signal Processing., vol. 65, no. 3, pp. 764–778, Feb 2017.
  * [6] P. Berger, G. Hannak, and G. Matz,  “Graph Signal Recovery via Primal-Dual Algorithms for Total Variation Minimization,”  IEEE Journal of Selected Topics in Signal Processing., vol. 11, no. 6, pp. 842–855, Sep 2017.
  * [7] K. Qiu, X. Mao, X. Shen, X. Wang, T. Li, and Y. Gu,  “Time-Varying Graph Signal Reconstruction,”  IEEE Journal of Selected Topics in Signal Processing., vol. 11, no. 6, pp. 870–883, Sep 2017.
  * [8] X. Wang, M. Wang, and Y. Gu,  “A Distributed Tracking Algorithm for Reconstruction of Graph Signals,”  IEEE Journal of Selected Topics in Signal Processing., vol. 9, no. 4, pp. 728–740, June 2015.
  * [9] E. Brugnoli, E. Toscano, and C. Vetro,  “Iterative reconstruction of signals on graph,”  IEEE Signal Processing Letters., vol. 27, pp. 76–80, 2020.
  * [10] R. Torkamani, and H. Zayyani,  “Statistical Graph Signal Recovery Using Variational Bayes,”  IEEE Trans. Circuit and Systems II: Express Briefs, vol. 68, no. 6, pp. 2232–2236, June 2021.
  * [11] R. Torkamani, H. Zayyani, and F. Marvasti  “Joint Topology Learning and Graph Signal Recovery Using Variational Bayes in Non-Gaussian Noise,”  IEEE Trans. Circuit and Systems II: Express Briefs, vol. 69, no. 3, pp. 1887–1891, March 2022.
  * [12] P. Di Lorenzo, S. Barbarossa, P. Banelli, and S. Sardellitti,  “Adaptive Least Mean Squares Estimation of Graph Signals,”  IEEE Trans. on Signal and Inf. Proc. over Networks, vol. 2, no. 4, pp. 555–568, Dec 2016.
  * [13] P. Di Lorenzo, P. Banelli, S. Barbarossa, and S. Sardellitti,  “Distributed Adaptive Learning of Graph Signals,”  IEEE Transaction on Signal Processing., vol. 65, no. 16, pp. 4193–4208, Aug 2017.
  * [14] P. Di Lorenzo, P. Banelli, E. Isufi, S. Barbarossa,, and G. Leus,  “Adaptive Graph Signal Processing: Algorithms and Optimal Sampling Strategies,”  IEEE Transaction on Signal Processing., vol. 66, no. 13, pp. 3584–3598, Jul 2018.
  * [15] P. Di Lorenzo, and E. Ceci,  “Online Recovery of Time- varying Signals Defined over Dynamic Graphs,”  EUSIPCO 2018., Rome, Italy, Sep 2018.
  * [16] M. Ramezani-Mayiami,  “Joint Graph Learning and Signal Recovery via Kalman Filter for Multivariate Auto-Regressive Processes,”  EUSIPCO 2018., Rome, Italy, Sep 2018.
  * [17] Y. Shen, G. Leus, and G. B Giannakis,  “Online Graph-Adaptive Learning with Scalability and Privacy,”  IEEE Transaction on Signal Processing., vol. 67, no. 9, pp. 2471–2483, May 2019.
  * [18] P. Di Lorenzo, S. Barbarossa, and S. Sardellitti,  “Distributed Adaptive Learning of Graph Processes via In-Network Subspace Projections,”  Asilomar 2019., Pacific Grove, CA, USA, Nov 2019.
  * [19] M. J. M. Spelta, and W. A Martins,  “Normalized LMS algorithm and data-selective strategies for adaptive graph signal estimation,”  Elsevier Signal Processing., vol. 167, pp. 2471–2483, Feb 2020.
  * [20] M. J. Ahmadi, R. Arablouei, and R Abdolee,  “Efficient Estimation of Graph Signals With Adaptive Sampling,”  IEEE Transaction on Signal Processing., vol. 68, no. 9, pp. 3808–3823, 2020.
  * [21] Y. Zhao, and E. Ayanoglu,  “Gaussian Kernel Variance for an Adaptive Learning Method on Signals Over Graphs,”  IEEE Trans. on Signal and Information Proc. over Networks, vol. 8, pp. 389–403, 2022.
  * [22] V. C. Gogineni, V. Naumova, S. Werner, and Y. F. Huang,  “Graph Kernel Recursive Least-Squares Algorithms,”  APSIPA ASC 2021., Tokyo, Japan, Dec 2021.
  * [23] D. Duttweiler,  “Proportionate normalized least-mean-squares adaptation in echo cancelers,”  IEEE Trans. Speech Audio Process., vol. 8, no. 5, pp. 508–518, Sep 2000.
  * [24] K. Wagner, and M. Doroslovacki,  “Proportionate-type normalized least mean square algorithms with gain allocation motivated by mean-square-error minimization for white input,”  IEEE Trans. on Signal Proc., vol. 59, no. 5, pp. 2410–2415, May 2011.
  * [25] S. H. Yim, H. S. Lee, and W. J. Song,  “A proportionate diffusion LMS algorithm for sparse distributed estimation,”  IEEE Trans. on Circuit and Systems-II: Express Briefs., vol. 62, no. 10, pp. 992–996, Oct 2015.
  * [26] H. Zayyani, et al, “A robust generalized proportionate diffusion LMS algorithm for distributed estimation,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 68, no. 4, pp. 1552–1556, 2020.
  * [27] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,  “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,”  IEEE Signal Process. Mag., vol. 30, no. 3, pp. 83–98, May 2013.
  * [28] I. Pesenson,  “Sampling in Paley-Wiener spaces on combinatorial graphs,”  Trans. Amer. Math. Soc., vol. 360, no. 10, pp. 5603–5627, 2008.
  * [29] S. Mula, V. Ch. Gogineni, and A. S. Dhar,  “Algorithm and VLSI architecture design of proportionate-type LMS adaptive filters for sparse system identiﬁcation,”  IEEE Trans. on Very Large Scale Integration (VLSI) Systems, vol. 26, no. 9, pp. 1750–1762, 2018.
  * [30] X. Pu, S. L. Chau, X. Dong, and D. Sejdinovic,  “Kernel-based graph learning from smooth signals: A fundamental viewpoint,”  IEEE Trans. on Signal and Inf. Process. oner Networks, vol. 7, pp. 192–207, 2021.
  * [31] P. Bodik, W. Hong, C. Guestrin, S. Madden, M. Paskin and R. Thibaux,  Intel lab data, [Online]. Available:http://db.csail.mit.edu/labdata/labdata.html., 2004.


# Complex In-Hand Manipulation via Compliance-Enabled Finger Gaiting and
Multi-Modal Planning

Andrew S. Morgan1∗, Kaiyu Hang2∗, Bowen Wen3, Kostas Bekris3, and Aaron M.
Dollar1

###### Abstract

Constraining contacts to remain fixed on an object during manipulation limits
the potential workspace size, as motion is subject to the hand’s kinematic
topology. Finger gaiting is one way to alleviate such restraints. It allows
contacts to be freely broken and remade so as to operate on different
manipulation manifolds. This capability, however, has traditionally been
difficult or impossible to practically realize. A finger gaiting system must
simultaneously plan for and control forces on the object while maintaining
stability during contact switching. This work alleviates the traditional
requirement by taking advantage of system compliance, allowing the hand to
more easily switch contacts while maintaining a stable grasp. Our method
achieves complete $SO(3)$ finger gaiting control of grasped objects against
gravity by developing a manipulation planner that operates via orthogonal safe
modes of a compliant, underactuated hand absent of tactile sensors or joint
encoders. During manipulation, a low-latency 6D pose object tracker provides
feedback via vision, allowing the planner to update its plan online so as to
adaptively recover from trajectory deviations. The efficacy of this method is
showcased by manipulating both convex and non-convex objects on a real robot.
Its robustness is evaluated via perturbation rejection and long trajectory
goals. To the best of the authors’ knowledge, this is the first work that has
autonomously achieved full $SO(3)$ control of objects within-hand via finger
gaiting and without a support surface, elucidating a valuable step towards
realizing true robot in-hand manipulation capabilities.

###### Index Terms:

Dexterous Manipulation, In-hand Manipulation, Manipulation Planning, Compliant
Joints and Mechanisms

Figure 1: We explore the development of a complete $SO(3)$ planner for within-
hand manipulation using finger gaits, by controlling two orthogonal extrinsic
rotation axes. Given a start configuration (cube face A), the proposed planner
finds an action sequence along the two controlled dimensions so as to reach
the desired goal configuration (cube face B). During manipulation, the pose of
the object is tracked via a low-latency, 6D pose object tracker, providing
feedback for online replanning and disturbance compensation via a recovery
phase that uses translation control. 00footnotetext: Manuscript received:
September 9, 2021; Revised: December 17, 2021; Accepted: January 6,
2022.00footnotetext: This paper was recommended for publication by Editor Hong
Liu upon evaluation of the Associate Editor and Reviewers’
comments.00footnotetext: This work was supported by the U.S. National Science
Foundation grants IIS-1752134, IIS-1734190, & IIS-1900681 to Yale University
and NRI-1734492 to Rutgers University. * Authors contributed
equally.00footnotetext: 1Andrew Morgan and Aaron Dollar are with the
Department of Mechanical Engineering & Materials Science, Yale University,
USA. (email: Andrew.Morgan, Aaron.Dollar}@yale.edu)00footnotetext: 2Kaiyu Hang
is with the Department of Computer Science, Rice University, USA. (email:
kaiyu.hang@rice.edu)00footnotetext: 3Bowen Wen and Kostas Bekris are with the
Department of Computer Science, Rutgers University, USA. (email: {bw344,
kostas.bekris}@cs.rutgers.edu)00footnotetext: Digital Object Identifier (DOI):
see top of this page.

##  I INTRODUCTION

Within-Hand Manipulation can be characterized as the ability to reorient or
reposition an object with respect to the hand frame [1]. The quest for this
capability has been studied in the robot manipulation community for
decades–from planning and control with rigid hands [2] to efforts with soft,
compliant, or underactuated hands [3]. Notably, the majority of these works
constrain contacts to remain fixed or rolling during manipulation. This
constraint limits the object’s workspace as the actuators can only operate on
a single hand-object configuration manifold. Alternatively, finger gaiting,
i.e., the process of repositioning contacts on the object during manipulation,
can help alleviate such limitations and extend the object’s available
workspace.

Finger gaiting is an inherently difficult task for a robot. Given a robot
hand, the individual serial link fingers must work in proper unison without
collision; maintaining stability while making and breaking contact with the
object [4]. The computational complexity of this problem has traditionally
been very expensive–requiring the system to calculate and modulate forces and
planned joint trajectories online during manipulation. We, conversely, are
able to alleviate many of these complexities by leveraging compliance, i.e.
safe modal transitions, in our system. The capabilities we present extend
beyond what has been illustrated previously in the literature, some purely in
simulation [5, 6] and some demonstrated on a real robot, but with support
surfaces [7, 8, 9] .

In this letter, we build off the observation that by leveraging the passive
adaptive properties of an underactuated hand, we can convert traditional
position, force, and stability control problems to a unified motion-only
control paradigm, based on the compliant properties of the hand. Concretely,
finger gaiting requires contacts to be constantly in motion. The presence and
thus activation of certain contacts constrains what forces the hand can impart
onto the object. This phenomena creates the abstraction of modes, which can be
conceptualized as different families of motion manifolds that are subject to
the system’s current set of constraints–typically in the form of gaining or
losing contacts. Multi-modal planning within and between these constraints
thus becomes a major focus of this work.

We in this work can constrain our planning approach according to the nature of
the task. Formally, any orientation in $SO(3)$ can be achieved via a three
action trajectory comprised of two orthogonal rotations. From this
formulation, we design two core modal actions for the hand and a multi-modal
planner that runs online and is constantly updated via a vision-based, low-
latency 6D pose object tracker [10, 11]. This method continually calculates a
trajectory from start to goal and when undesired scenarios arise, such as
contact slip, the planner updates and suggests new actions accordingly. We
incorporate this approach into an open-source and underactuated hand with four
fingers [12]. In the end, we showcase the efficacy of our system through
various experiments: tracking the planned and executed trajectory of an
object, evaluating the recovery potential given undesired perturbations, and
finally, the ability to extend to novel object geometries.

The contributions in this letter are threefold. First, we develop a complete
and fast planning solution for in-hand reorientation using two extrinsic
rotation axes. Secondly, we describe the utility of compliance for switching
between modes and how it “inflates” the contact switching region. And,
finally, we present a simple yet effective robot system capable of complex
finger gaiting capabilities, underscoring continued discussion in the
community on the utility of compliance for in-hand manipulation tasks.

##  II RELATED WORK

####  II-1 Modeling Manipulation

Modeling robot manipulation has historically been arduous as the dynamics
associated with contact are difficult to predict in novel scenarios. From an
in-hand manipulation perspective, various levels of modeling have been
investigated – from contact models [13] and fingerpad curvature models [14],
to hand kinematic models and whole hand-object system models [15]. While these
approaches have elucidated many powerful hand-object relationships, inaccurate
parameterizations can often lead to task failure [16]. To help alleviate such
uncertainty, several end effectors leverage designs that are soft or
underactuated [17]. This property provides passive reconfigurability to the
mechanism, which can absorb much of the “slack” traditionally required to be
fully accounted for in system modeling and can further reduce grasp planning
times [18, 19]. In this work, we leverage such mechanisms to unify our
planning approach–creating the notion of safe regions for contact switching.

####  II-2 Extending In-hand Manipulation Capabilities

Models for in-hand manipulation are typically constrained to fixed or rolling
contact scenarios. While this advantageously simplifies assumptions for
control, it also limits the object’s available workspace to be dependent on
the kinematic topology of the hand. Without relying on external contact, e.g,
[20], or task-specific hands with roller-based fingers, e.g [21], two general
approaches can help alleviate such constraint: sliding manipulation and finger
gaiting. Leveraging the former is very difficult, as detecting and controlling
its nonlinear conditions requires various levels of advanced sensing [22, 23].
The latter, alternatively, has been largely difficult due to computational
considerations, but has been made more successful in recent years [24, 6]. The
work in [25] used finger gaits and tactile sensors to maintain grasp
stability. [7] used a 24-DOF hand with a motion capture system, in addition to
“over 100 years” of simulated data to perform impressive and fast cube
manipulation in the palm. [8] leveraged the capabilities of a soft hand with
16 degrees of actuation for similar types of manipulation. While as impressive
as these aforementioned works are, we are interested in extending beyond these
capabilities to work against gravity, i.e. maintaining object stability
without a support surface, and while utilizing a simple hand with very little
onboard sensing.

Figure 2: Multi-modal planning problems can be conceptualized as operating in
different configuration manifolds. (a) Given a single manifold $\mathcal{M}$,
the planner must find a path along the constrained layer. (b) A single mode,
or foliation, can have multiple manipulation manifolds depending on the start
configuration of the hand-object system. (c) Switching modes is possible when
the system finds a configuration, $q^{\prime}$, that lies on both manifolds.
(d) Finding a path from start configuration, $q_{s}$, to goal configuration,
$q_{g}$, can require multiple jumps between modes. Note that although
represented in this figure as folds, we assume modes to be differentiable but
not necessarily euclidean.

####  II-3 Multi-Modal Planning for Manipulation

Problems in robot manipulation are frequently multi-modal, i.e. the seemingly
continuous problems have an underlying discrete structure guided by
constraints [26]. These constraints are typically imposed by the nature of
contact, and by discretizing planning in terms of these manifolds, the planner
search space is confined according to physical constraints [27]. Numerous
sampling-based multi-modal planners have been described in the literature,
which are able to generalize well, particularly in high-DOF scenarios. Though,
sampling-based methods without informed exploration are vastly inefficient and
suffer from time-complexities associated with over-exploration. Recent work
has attempted to address this issue by using informed “leads” to guide
exploration [28]. In our work, we build off these observations and constrain
our planner’s search according to physical properties of the $SO(3)$ rotation
group. That is, our developed planner constrains the number and nature of mode
switching to reduce planning time for continual updates during online
execution.

##  III PRELIMINARIES

In this section, we introduce the preliminaries associated with multi-modal
motion planning. We first discuss constrained manifolds and then develop
notation and terminology for modal switching leveraging compliance.

###  III-A Constrained Within-Mode Planning

Consider a robot with configuration space, $\mathcal{Q}\subset\mathbb{R}^{N}$,
where $N$ is the number of joints and where $\mathcal{Q}$ can completely
define the state of the robot. Now, consider an object with configuration
space, $\mathcal{O}\subset SE(3)$, which is located on a support surface, e.g.
table top. For simplicity, let’s disregard other potential collisions in the
environment. In free space, the current robot configuration,
$q\in\mathcal{Q}$, is unconstrained and thus able to freely move. Though, when
$q$ causes links of the robot to come in contact with the object in its
current configuration, $o\in\mathcal{O}$, a constraint is imposed on the
system. The robot cannot push the object through the support surface. If, for
instance, the goal of the robot is to slide the object along its support
surface, the robot’s motion is thus constrained to a contact configuration
manifold, $\mathcal{M}$ (Fig. 2(a)).

Assume there are several manifolds, indexed by $\eta$. The available motion
manifold of the robot, $\mathcal{M}_{\eta}$, is subject to its associated
constraint function
$F^{\eta}(q,o):\mathbb{R}^{N}\rightarrow\mathbb{R}^{k_{\eta}}$ where $(1\leq
k_{\eta}<N)$. Modes are thus defined as a family, or foliation,
$\mathcal{F}^{\zeta}$ , of constrained manifolds that share the same, or
similar, constraint functions (Fig. 2(b)). Inside each $\mathcal{F}^{\zeta}$,
we assume all $\mathcal{M}_{\eta}$ to be smooth and differentiable. Therefore,
planning in the $\mathcal{M}_{\eta}\in\mathbb{R}^{(N-k_{\eta})}$ manifold of
$\mathbb{R}^{N}$ is subject to,

$\mathcal{M}_{\eta}=\\{q\in\mathcal{Q}|F^{\eta}(q,o)=\textbf{0}\\}$ (1)

###  III-B Safe Mode Transitions

Multi-modal planning requires a robot to transition between at least two
constraint manifolds, which likely lie in different foliations, in order to
reach the desired goal. Consider a robot starting in configuration $q_{s}$,
located in manifold $\mathcal{M}^{\zeta_{1}}_{\eta_{1}}$. Let’s now assume the
goal configuration $q_{g}$ is located in a different mode and thus along a
different manifold, i.e. $\mathcal{M}^{\zeta_{2}}_{\eta_{2}}$ (Fig. 2(c)). In
order for a robot to transition to $q_{g}$, a planner must find a candidate
submanifold, $\mathcal{M}_{\eta_{1}\cap\eta_{2}}$, that enables a transition.
Inside of this submanifold lies candidate configuration, $q^{\prime}$, for
transitioning. We can define $\mathcal{M}_{\eta_{1}\cap\eta_{2}}$ according
to,

$\mathcal{M}_{\eta_{1}\cap\eta_{2}}=\\{q\in\mathcal{Q}|F^{\eta_{1}}(q,o)=\textbf{0},F^{\eta_{2}}(q,o)=\textbf{0}\\}$
(2)

Note that it is possible for $\mathcal{M}_{\eta_{1}\cap\eta_{2}}$ to be empty
when there exists no $q^{\prime}$ to transition
$\mathcal{M}^{\zeta_{1}}_{\eta_{1}}$ to $\mathcal{M}^{\zeta_{2}}_{\eta_{2}}$.

In many robot scenarios, if a transition is possible from
$\mathcal{M}^{\zeta_{1}}_{\eta_{1}}$ to $\mathcal{M}^{\zeta_{2}}_{\eta_{2}}$,
the transfer region typically has a volume near zero, i.e.
$\Omega(\mathcal{M}_{\eta_{1}\cap\eta_{2}})\approx 0$ where $\Omega$ is the
convex hull function. There are situations, however, where we can relax this
constraint to form safe mode transitions. Formally, the manifold’s constraint
function, $F^{\eta}(q,o)$, is dependent on the robot’s current configuration,
$q$ and object configuration, $o$. Though due to near-unmodelable DOFs in
soft, compliant, or underactuated robots or objects, the controlled
configuration of the robot may still hold for
$\mathbb{R}^{N}\rightarrow\mathbb{R}^{k_{\eta}}$, but the true configuration
of the system may be non-deterministic. Thus, there can exist a relaxation of
this traditional modal switching constraint, which in turn can simplify and
accelerate planning [18]. We define a modal transition threshold vector,
${\rho}$, that defines an inflation of the constraint manifold along different
axes, which inherently increases the potential volume of
$\mathcal{M}_{\eta_{1}\cap\eta_{2}}$ to develop a safe transfer region,

$\mathcal{M}_{\eta\cap\eta^{\prime}}^{*}=\\{q\in\mathcal{Q}|-\rho<F^{\eta}(q)<\rho,-\rho<F^{\eta^{\prime}}(q)<\rho\\}$
(3)

The success of our finger gating platform in this work largely leverages this
concept. Explicitly defining $\rho$, however, is difficult as it is system
dependent, i.e. properties of the hand and the object determine how inflated
the switching regions become. As a proof-of-concept analysis, we estimate
$\rho$ via experimentation in Sec. VII and will leave computational
investigation for future work.

Figure 3: (a) Our planning solution for $SO(3)$ begins by enumerating outward
the goal orientation manifold by step size $\sigma_{R}$ and creating a KD-
Tree. (b) After creation, a forward search outward is performed from the start
orientation, $R_{s}$, combining a $z-$axis rotation first, then an $x-$axis
rotation. (c) This search continues until a candidate orientation is within
some distance $\rho$ from the expanded goal manifold. (d) Finally, the entire
plan is enumerated from start to goal along the acquired trajectory with
appropriate timestamps and step size $\sigma_{R}$.

##  IV ORIENTATION-BASED MOTION PLANNING

We formally define the problem of orientation-based motion planning and the
minimum number of actions required. Then, we utilize this as a guarantee for
the completeness our planner’s modal search transitions.

###  IV-A Proper Rotations in $SO(3)$

Proper rotations are sets of angles that can fully define any rotation in the
group, $SO(3)$. Let’s assume we have an object with some rotation, $R$. It is
possible to transition $R$ to any rotation configuration via a combination of
rotations along two orthogonal axes. Though there are six total combinations,
we will instantiate our planner to focus on rotations along the $z-$axis
($R_{z}$) and $x-$axis ($R_{x}$) to serve as an example.

Theorem 1: Let the rotational configuration of the object be $R\in SO(3)$.
Then there exists at least one
$(\phi,\theta,\psi)\in[0,2\pi)\times[0,\pi]\times[0,2\pi)$ such that
$R=R_{z}(\phi)R_{x}(\theta)R_{z}(\psi)$ [29].

###  IV-B Planning Orientation Transitions

Following this realization and the preliminaries in Sec. III, we develop a
complete multi-modal planning algorithm for finger gaiting that consists of
three rotations along two orthogonal axes. If, for instance, all three
orthogonal axes were available, the planning problem would be mostly
trivial–simply transition the object along each axis as much as necessary to
reach the goal. For the duration of this paper, we will refer to rotations
about axes as being extrinsic, i.e. according to the hand frame of robot, and
will adopt the notation that $R_{x}(\cdot)$ and $R_{z}(\cdot)$ are rotations
of the object about the $x-$ and $z-$axes of the hand frame, or roll and yaw,
respectively.

Algorithm 1 $SO3Plan(\cdot)$

1:$R_{s},R_{g},\rho,\sigma_{R}$ $\triangleright$ object start orientation,
object goal orientation, modal transition threshold, rotational step size

2:$\mathcal{P}$ $\triangleright$ multi-modal plan

3:$\Sigma=[0,\sigma_{R},-\sigma_{R},2\sigma_{R},-2\sigma_{R},\cdots,\pi,-\pi]$
$\triangleright$ candidate steps

4:$\mathcal{M}_{z_{2}}\leftarrow[R_{z}(\phi)R_{g}$ for $\phi$ in $\Sigma$]
$\triangleright$ expanded goal manifold

5:$KD\leftarrow KDTree.build(\mathcal{M}_{z_{2}})$ $\triangleright$ KD-Tree of
goal manifold

6:for $\psi$ in  $\Sigma$ do

7: for $\theta$ in  $\Sigma$ do

8: $R^{*}=R_{x}(\theta)R_{z}(\psi)R_{s}$ $\triangleright$ outward expanding
search

9: $\phi\leftarrow KD.query\\_closest\\_config(R^{*})$

10: if $dist_{R}(R_{z}(\phi)R_{g},R^{*})<\rho$ then $\triangleright$ distance
function

11: break loops $\triangleright$ bidirectional trajectory connects

12: end if

13: end for

14:end for

15:$\mathcal{P}\leftarrow\textbf{planner}.enumerate(R_{s},\sigma_{R},\psi,\theta,\phi)$

$\triangleright$ enumerate path outward according to connection (Fig. 3(d))

16:return $\mathcal{P}$

Let’s assume the capabilities of the robot enable the formulation of just two
mode foliations, which are able to take an object in any configuration and
rotate about the $x-$ and $z-$axes, respectively. Given an object in start
configuration $o_{s}=\\{T_{s}\in\mathbb{R}^{3},R_{s}\in SO(3)\\}$, the
objective of this planner is to reach a desired goal pose,
$o_{g}=\\{T_{g}\in\mathbb{R}^{3},R_{g}\in SO(3)\\}$. For now, we will
disregard the translational positions, $T_{s}$ and $T_{g}$, and focus on
rotation.

As to accelerate computation, the proposed planning method is bi-directional
in that it builds a search trajectory outward from both, the start and goal
configurations. More formally, there exists rotations such that
$R_{g}=R_{z_{2}}(\phi)R_{x}(\theta)R_{z_{1}}(\psi)R_{s}$ where
$R_{z_{1}}(\psi)$ is the first rotation about the $z-$axis, $R_{x}(\theta)$ is
the rotation about the $x-$axis, and $R_{z_{2}}(\phi)$ is the second rotation
about the $z-$axis. Thus, the goal of this planner is to solve for
$\psi,\theta,$ and $\phi$ and enumerate a trajectory accordingly.

Alg. 1 begins with the backward pass by populating a $z-$axis rotation outward
from $R_{g}$ with step size, $\sigma_{R}$, which is user-defined (Fig. 3(a)).
For our implementation, $\sigma_{R}$ can generally range anywhere from
0.5-5.0$\degree$ and can be adaptive according to system architectures. Let’s
say that these discretized rotations outward from $R_{g}$ can sufficiently
define the entire manifold, $\mathcal{M}_{z_{2}}$, i.e. an entire rotation of
the $z-$axis about $R_{g}$. From these configurations in
$\mathcal{M}_{z_{2}}$, we can build a KD-Tree, $KD$, which is a data structure
that enables efficient querying of the rotations included in this manifold.

Now that the backwards direction is instantiated, the goal of the forward pass
is to find rotations $R_{x}(\theta)R_{z_{1}}(\psi)R_{s}$ such that they can
“connect” to the manifold $\mathcal{M}_{z_{2}}$ represented by $KD$ within
some distance threshold $\rho$ from (3) (Fig. 3(c)). As aforementioned, $\rho$
is the inflation factor for contact switching, and is generally system
dependent. Here, we utilize a user-defined rotational distance function
$dist_{R}(\cdot)$, where for our instantiation, we calculate the L2-norm of
the difference vector along the roll, pitch, and yaw dimensions of two
rotations. Though, we note that there are other metrics available, e.g. [30].
Tangibly, this forward pass sequentially searches through candidate values for
$\psi$ and $\theta$ until the trajectory reaches the goal manifold
$\mathcal{M}_{z_{2}}$, which is found by continually querying $KD$. Once this
connection is found (Fig. 3(d)), values for $\psi$, $\theta$, and $\phi$ are
recorded, and a trajectory for plan $\mathcal{P}$ is enumerated with
timestamps along step size $\sigma_{R}$ from start to goal.

##  V Object Control

Finger-gaiting is a dynamic and uncertain process; making and breaking
contacts can often lead to slip which could further result in object ejection.
To mitigate such occurrences, our solution is to continually replan and
execute modal actions. More specifically, we rely on an online servoing-based
framework so as to reach the desired goal. When the object enters a region of
potential instability, such as being distant from the workspace center, a
recovery phase comprised of translational modal actions is enacted to relocate
the object.

####  V-1 Translational Planning

Planning translations within-hand, presented in Alg. 2, performs a simple,
greedy visual servoing control approach. Generally, the algorithm calculates
the translational difference, $\gamma$, along each of the coordinate axes from
start position, $T_{s}$, and goal position, $T_{g}$. The planner then chooses
the direction of largest difference and creates a plan to take a single step
in that direction [31].

Algorithm 2 $TranslationPlan(\cdot)$

1:$T_{s},T_{g},\sigma_{T}$ $\triangleright$ object position, object goal
position, step size

2:$\mathcal{P}$ $\triangleright$ plan

3:$\gamma.x\leftarrow T_{g}.x-T_{s}.x$ $\triangleright$ translational
difference along $x-$axis

4:$\gamma.y\leftarrow T_{g}.y-T_{s}.y$ $\triangleright$ translational
difference along $y-$axis

5:$\gamma.z\leftarrow T_{g}.z-T_{s}.z$ $\triangleright$ translational
difference along $z-$axis

6:$\mathcal{P}\leftarrow\\{\\}$

7:$\pi\leftarrow max(|\gamma_{x}|,|\gamma_{y}|,|\gamma_{z}|)$ $\triangleright$
direction of maximum deviation

8:if $\pi$ is $|\gamma_{x}|$ then

9: $\mathcal{P}\leftarrow[T_{x}+$sgn$(\gamma_{x})\sigma_{T},T_{y},T_{z}]$
$\triangleright$ desired move $\pm T_{x}$

10:else if $\pi$ is $|\gamma_{y}|$ then

11: $\mathcal{P}\leftarrow[T_{x},T_{y}+$sgn$(\gamma_{y})\sigma_{T},T_{z}]$
$\triangleright$ desired move $\pm T_{y}$

12:else if $\pi$ is $|\gamma_{z}|$ then

13: $\mathcal{P}\leftarrow[T_{x},T_{y},T_{z}+$sgn$(\gamma_{z})\sigma_{T}]$
$\triangleright$ desired move $\pm T_{z}$

14:return $\mathcal{P}$

####  V-2 Orientation and Translation Control

Our control approach combines both, Alg. 1 and Alg. 2. This method requires a
goal threshold, $\tau$, a linear interpolant update parameter, $\lambda$, and
a priori knowledge of the hand’s workspace center, $T_{center}$, for object
recovery.

In executing Alg. 3, orientation control is accomplished first. A plan
$\mathcal{P}$ is found according to the object’s current object rotation, $R$.
The robot then takes a single action along $\mathcal{P}$ via the
robot$.mode\\_action(\cdot)$ method. This system-specific function enacts a
motor position step proportional to $\sigma_{R}$ along the desired mode. From
this action, an object pose displacement, $\delta$, is calculated by taking
the distance between the original orientation, $R$, and the new orientation,
$R^{\prime}$, after the action. This value $\delta$ serves as a differential
update term along the orientation transition, which allows us to adaptively
re-estimate $\sigma_{R}$ according to an interpolant update parameter
$\lambda$. Since, as aforementioned, there is inherent uncertainty in the task
of finger gaiting, this adaptive update to the transition step size
$\sigma_{R}$ serves to aid in our plan’s accuracy for the true motion of the
hand-object system. The algorithm continues while $R$ is outside of the goal
distance threshold $\tau_{R}$. Notably, a recovery phase (lines 12-16) is
developed inside of this loop and is enacted when the object position, $T$, is
distance $\tau_{T}$ outside of $T_{center}$. At this point, the system takes
greedy transition steps towards the center of the workspace for recovery.

Alg. 3 concludes by performing steps along $\mathcal{P}$ found from the
$TranslationPlan(\cdot)$ method so as to reposition the object within-hand
after object reorientation. We implement this method last in an attempt to
decouple orientational and translational planning for the final pose.

Algorithm 3 $ObjectControl(\cdot)$

1:$o_{g}$, $\rho$, $\sigma$, $\tau$, $\lambda$, $T_{center}$ $\triangleright$
object goal pose, modal transition threshold vector, pose transition step
size, goal pose reached threshold, interpolant step size, center of object
workspace

2:$o_{f}$ $\triangleright$ final object pose

3:$T_{g},R_{g}\leftarrow o_{g}.T,o_{g}.R$ $\triangleright$ translation and
rotation of goal pose

4:$o\leftarrow\textbf{tracker}.perceive()$ $\triangleright$ object pose $\in
SE(3)$

5:$T,R\leftarrow o.T,o.R$ $\triangleright$ translation and rotation of object
pose

6:while $dist_{R}(R_{g},R)>\tau_{R}$ do $\triangleright$ orientation control

7: $\mathcal{P}\leftarrow$planner$.SO3Plan(R,R_{g},\rho,\sigma_{R})$
$\triangleright$ Alg. 1

8: robot$.mode\\_action(\mathcal{P},\sigma_{R}$) $\triangleright$ execute
modal action

9: $o^{\prime}\leftarrow\textbf{tracker}.perceive()$ $\triangleright$ new
object pose

10: $T^{\prime},R^{\prime}\leftarrow o^{\prime}.T,o^{\prime}.R$
$\triangleright$ new object translation and rotation

11: $\delta\leftarrow dist_{R}(R^{\prime},R)$ $\triangleright$ size of last
action step

12: $\sigma_{R}\leftarrow\sigma_{R}+\lambda(\sigma_{R}-\delta)$
$\triangleright$ update rotational transition step

13: $T,R\leftarrow T^{\prime},R^{\prime}$

14: while $dist_{T}(T$,$T_{center})>\tau_{T}$ do $\triangleright$ object
recovery

15:
$\mathcal{P}\leftarrow\textbf{planner}.TranslationPlan(T,T_{center},\sigma_{T})$

16: robot$.mode\\_action(\mathcal{P},\sigma_{T}$) $\triangleright$ execute
modal action

17: $o\leftarrow\textbf{tracker}.perceive()$ $\triangleright$ object pose

18: $T,R\leftarrow o.T,o.R$ $\triangleright$ object translation and rotation

19: end while

20:end while

21:while $dist_{T}(T_{g},T)>\tau_{T}$ do $\triangleright$ translation planning

22:
$\mathcal{P}\leftarrow\textbf{planner}.TranslationPlan(T,T_{g},\sigma_{T})$
$\triangleright$ Alg. 2

23: robot$.mode\\_action(\mathcal{P},\sigma_{T}$) $\triangleright$ execute
modal action

24: $o\leftarrow\textbf{tracker}.perceive()$ $\triangleright$ object pose

25: $T,R\leftarrow o.T,o.R$ $\triangleright$ object translation and rotation

26:$o_{f}\leftarrow\textbf{tracker}.perceive()$ $\triangleright$ final object
pose

27:return $o_{f}$

####  V-3 Generalization and Practical Algorithm Modifications

The authors present this solution as a generalized approach to planning for
in-hand manipulation. Surely, the design and realization of the
robot$.mode\\_action(\cdot)$ method begs questions for system-specific
scenarios, but can be found in a variety of ways: analytical modeling,
dynamics learning, reinforcement learning, etc. While we do not focus on
describing the underlying components of this function in this work, we largely
leverage an energy-modeling technique outlined in [32] for determining
predicted object transitions.

Moreover, the proposed planner can be further accelerated by subtle variations
to the aforementioned algorithms. First, $KD$ only needs to be computed once
since the goal pose is static. Also, initialized step sizes for $\sigma$ can
deviate between modes according to system-specific properties, e.g. when
gravity affects the rotation for one mode more than another.

##  VI SYSTEM SETUP

Figure 4: (a) Our system is comprised of a RGBD camera fixated to the robot’s
environment. (b) During manipulation, the pose of the object, e.g. bunny, is
tracked via a 6D pose object tracker. (c) A bottom view of the Yale Model Q
illustrates the 110$\degree$ abduction capabilities of the differentially
coupled fingers.

####  VI-1 Robot Setup

We utilize the open-source Model Q, an underactuated hand from the Yale
OpenHand Project. The hand is equipped with four two-link fingers and four
total actuators. A single motor controls the actuation of an opposing set of
fingers with joints comprised of flexures, guided by a differential. This
enables passive reconfigurability between these two fingers and into/out-of
the plane of grasping. These two coupled fingers are connected to an actuated
rotary joint located within the palm, allowing $110\degree$ rotation about the
$z-$axis of the hand (Fig. 4). The final two motors serve to actuate the
remaining two fingers individually via tendons. Being underactuated, the
hand’s joint configuration cannot be determined directly, as it is not
equipped with joint encoders or tactile sensors. See [12] for additional
design information.

The hand is mounted on a 7-DOF Barrett WAM arm. External to the robot, an RGBD
Intel Realsense D415 camera is calibrated to the robot’s environment. Its
imagery is processed through the low-latency (60Hz) 6D pose object tracker
(Fig. 4). The tracker is robust to finger occlusion, has approximately
$\pm$1.5$mm$/ 2$\degree$ of noise along any axis, and is trained solely with
synthetic object data. See [10] for additional information.

Figure 5: The Yale OpenHand Model Q is capable of operating along four modes
for fingertip-based manipulation. Modes are shown with transitions from the
start configuration (center).

####  VI-2 Mode Design

As described in Sec. IV-B, a series of rotations about two orthogonal axes
allow an object to reach any orientation in $SO(3)$. Given the hand topology
of the Model Q, an extrinsic $z-$axis rotation is easily possible via the
palm’s rotary joint. Moreover, an extrinsic $x-$axis rotation is also possible
via coordinated finger motions; grasping with the differentially actuated
opposing pair and pushing the object from either individually driven fingers
as validated by [32].

In addition to the two rotational motions, the hand topology furthermore
allows for only two translational modes. The first can translate the object
along the $z-$axis by squeezing the object towards the palm and regrasping, or
by leveraging small amounts of coordinated slip, and simple $y-$axis
translation is possible via a method described in [31] (Fig. 5). Note that
translation about the $x-$axis is not possible due to the actuation coupling
of the rotary set of fingers. By combining these two translational modes, we
are able to provide a recovery phase for the object (Fig. 1, 8), repositioning
the object towards the middle of the workspace so that manipulation can safely
continue. Note that $x-$axis translation is thus omitted.

Figure 6: (a) We experimentally validate our method with 5 different objects,
namely a cube, sphere, toy truck, Stanford Bunny, and toy duck. (b)
Tessellated faces of the cube show poses of the affixed letters.

##  VII EXPERIMENTS

Figure 7: The safety of mode transitions are attributed to the
reconfigurability of the underactuated mechanism. (a) Starting rotations
around different axes of the object elicit different amounts of
reconfiguration upon regrasping. (b) This is particularly apparent in the yaw
($z-$axis) direction of the object, where the object follows the minimum
energy configuration of the mechanism, and hence allowing us to estimate
$\rho$ from Eq. (3).

####  VII-1 Safe Modes Characterization

The passive adaptive nature of compliant hands is particularly beneficial for
the task of finger gaiting. Simply, small errors associated with state
estimation, modeling, and control can likely be dismissed by the system
“absorbing” the slack. It moreover determines the “inflation” of switching
regions, denoted by $\rho$. To help quantify and validate this adaptive
nature, we evaluate how a regrasp action, i.e. transferring the grasp between
sets of opposing fingers, affects the resultant object configuration.

Evaluated with the cube (Fig. 6), we perturb the object along each axis.
Subsequently, the hand transfers the grasp to the opposite set of opposing
fingers. In doing so, we record the amount of object reconfiguration along
each axis via the object tracker. Notably, with the roll orientation of the
object, we see a small amount of reconfiguration, but in the yaw orientation,
we note the most reconfiguration (Fig. 7), due to the geometric properties of
the cube. More specifically, from the start grasp (blue circle), to the final
regrasped configuration (red square), the object tends to fall towards the
minimum energy configuration of the system (red region) [32, 24]. We note that
although there is a large variability in the object’s start orientation ($\pm
0.5$ radians), the hand is able to regrasp the object successfully in all
cases. This underscores our safe modes discussion, and generally allows us to
estimate how large $\rho$ can be. For safety, we define $\rho=0.2$ radians.

Figure 8: We execute a single planned trajectory 8 times and record its
repeatability. (a) Controlled dimensions of the object trajectory, such as
roll and yaw, follow closely in all trials, where as uncontrolled dimensions,
such as $x-$axis translation and pitch, are allowed to drift. (b) Modes are
enacted at different times in the manipulation, including recovery phases when
the object is not in the center of the graspable workspace.  Figure 9: We test
system robustness by (a) following through an extended trajectory of cube
faces F, A, B, C, D, E, and F once again, and (b) deliberately applying
different perturbations to the object along controlled dimensions so as to
require online replanning and recovery.

####  VII-2 Single Trajectory Execution

We evaluate the repeatability of executing a single trajectory of the cube.
Here, the goal of the task is to rotate the cube from face A to face B, and
when completed, translate the object 1.8cm along the negative $y-$axis. The
resultant executions, depicted in Fig. 8(a), shows the average trajectory of
the object over 8 different trials. During this task, we note how closely
trajectories follow along the roll and yaw dimensions, as these are controlled
via our safe modes. Object recovery occurs at similar times in the trials,
where the system transitions the object along the positive $z-$axis before
completing the $x-$axis rotations (Fig. 8(b)). At the end of the orientation
control sequence, the object is appropriately translated along the $y-$axis to
reach its desired goal pose. The average trajectory data, depicted in solid
black, does not extend past the 55 second mark due to one trial finishing
earlier than others. The authors have, for clarity, extended the trajectory in
the $y-$axis translation using a dotted line to illustrate the average final
pose.

####  VII-3 Continuous Goal Trajectories and Robustness

We showcase the robustness of our method by evaluating both, an extended
control trajectory and by disturbing the object pose during execution. As
depicted in Fig. 9(a), we execute a trajectory starting on cube face F,
transitioning through faces A,B,C,D,E and finally reaching the face F once
again. In all cases, the object was able to reach the goal face within
8$\degree$.

Moreover, we illustrate our ability to replan control trajectories when
unmodeled perturbations occur (Fig. 9(b)). At different points during
manipulation, we randomly perturb the object pose 13 total times along the
least constrained dimensions, the $x-$axis rotation and the $z-$axis
translation. Recovering and adapting to these occurrences over the course of a
205 second trajectory execution, the system completes the object reorientation
task within 6$\degree$ of the goal face.

Figure 10: Our manipulation planner is able to extend to objects of convex and non-convex geometries, and reach any orientation in $SO(3)$. Obj. |  $err$ (deg) | Plan time (s) | Total time (s) | Success
---|---|---|---|---
sphere | 4.6 | 13.2 | 153.1 | 5/5
truck | 4.2 | 10.8 | 112.3 | 5/5
bunny | 9.5 | 14.0 | 128.0 | 2/5
duck | 8.9 | 7.1 | 73.2 | 3/5
Figure 11: Experimental evaluation of four different objects with error in
final orientation, planning time, total time, and number of execution
successes.

####  VII-4 Generalization to Object Geometry

Finally, we test the transferability of our method by varying object
geometries (Fig. 6). In doing so, recovery properties, such as $T_{center}$
from Alg. 3, were recalculated, and the 6D pose object tracker was retrained
with new synthetic object data; other aforementioned parameters, including
$\rho$, remained the same. Depicted in Fig. 10, we evaluate simple, purely
rotational trajectories about each of the objects. Notably, the sphere easily
followed each of its guided reference trajectories, and was able to reach goal
orientations within $\pm 5\degree$ along any axis. The toy truck, further
illustrated our capabilities, transitioning to its goal configuration in
approximately 110 seconds with 4 recovery phases. And finally, the Stanford
Bunny and the toy duck were the most difficult to manipulate due to their non-
convex properties. With the bunny, fingers would often get stuck behind its
ears while making and breaking contact. Consequently, this occurrence caused
the hand-object system to get stuck in specific object orientations and
required restart. Moreover, during duck manipulation, the size of the duck’s
head affected the object’s center of mass. This made some actions, such as
$x-$axis rotation, difficult as the motion was then more dynamic in nature.
Only a simple, two-rotation trajectory could be completed with the duck,
resulting in a shorter execution time. The precision of these two
manipulations was decreased to $\pm 10\degree$ along any rotation axis and
were consequently not quite as robust as the other objects. We provide
benchmarking baselines for our experimentation in Fig. 11 according to
guidelines outlined [33]. Please refer to the supplementary video for
experimental evaluations.

##  VIII DISCUSSIONS AND FUTURE WORK

In this work, we developed a robust and complete solution to $SO(3)$ finger
gait planning. While instantiated on an underactuated hand in this work, the
methods described can generalize to other hand-object systems. We showcased
our method with numerous tasks–highlighting the safety of our modes, the
repeatability of our trajectories, the robustness of our planning and control
approach, and our ability to generalize to new object scenarios. This work
builds on a promising approach of vision-based control for compliant within-
hand manipulation. The tasks completed in this letter extend beyond what has
been possible in previous works.

Although we showcase novel capabilities, there are some limitations. First,
hand design plays a crucial role in our manipulation capabilities. We
conceptualize an altered hand that would enable rotation about the $y-$axis,
thus shortening trajectories. While our method may have relied less on data
and/or advanced sensing, we note that our manipulation actions were
significantly slower than some other related works [7, 8]. Moreover, we
discuss the “inflation” parameter $\rho$ quite extensively, which begs for
further theoretical investigation. To this end, the development of geometry-
focused modal actions and transitions is an interesting avenue. Lastly, the
current tracking schema requires an object CAD model for training. We are
interested in combining object-agnostic tracking methods [34] together with
the current object-agnostic controller for instant application to novel
objects.

## References

  * [1] A. Bicchi, “Hands for dexterous manipulation and robust grasping: A difficult road toward simplicity,” _IEEE Transactions on robotics and automation_ , vol. 16, no. 6, pp. 652–662, 2000.
  * [2] J. Kerr and B. Roth, “Analysis of multifingered hands,” _The International Journal of Robotics Research_ , vol. 4, no. 4, pp. 3–17, 1986.
  * [3] A. Sintov, A. S. Morgan, A. Kimmel, A. M. Dollar, K. E. Bekris, and A. Boularias, “Learning a state transition model of an underactuated adaptive hand,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1287–1294, 2019.
  * [4] L. Han and J. C. Trinkle, “Dextrous manipulation by rolling and finger gaiting,” in _Proceedings. 1998 IEEE International Conference on Robotics and Automation (Cat. No. 98CH36146)_ , vol. 1. IEEE, 1998, pp. 730–735.
  * [5] B. Sundaralingam and T. Hermans, “Geometric in-hand regrasp planning: Alternating optimization of finger gaits and in-grasp manipulation,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 231–238.
  * [6] G. Khandate, M. Haas-Heger, and M. Ciocarlie, “On the feasibility of learning finger-gaiting in-hand manipulation with intrinsic sensing,” 2021.
  * [7] O. M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew, J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray _et al._ , “Learning dexterous in-hand manipulation,” _The International Journal of Robotics Research_ , vol. 39, no. 1, pp. 3–20, 2020.
  * [8] A. Bhatt, A. Sieler, S. Puhlmann, and O. Brock, “Surprisingly robust in-hand manipulation: An empirical study-supplementary material,” 2021.
  * [9] S. Abondance, C. B. Teeple, and R. J. Wood, “A dexterous soft robotic hand for delicate in-hand manipulation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5502–5509, 2020.
  * [10] B. Wen, C. Mitash, B. Ren, and K. E. Bekris, “se (3)-tracknet: Data-driven 6d pose tracking by calibrating image residuals in synthetic domains,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 10 367–10 373.
  * [11] B. Wen, C. Mitash, S. Soorian, A. Kimmel, A. Sintov, and K. E. Bekris, “Robust, occlusion-aware pose estimation for objects grasped by adaptive hands,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6210–6217.
  * [12] R. R. Ma and A. M. Dollar, “An underactuated hand for efficient finger-gaiting-based dexterous manipulation,” in _2014 IEEE International Conference on Robotics and Biomimetics (ROBIO 2014)_. IEEE, 2014, pp. 2214–2219.
  * [13] J. C. Trinkle, J.-S. Pang, S. Sudarsky, and G. Lo, “On dynamic multi-rigid-body contact problems with coulomb friction,” _Journal of Applied Mathematics and Mechanics_ , vol. 77, no. 4, pp. 267–279, 1997.
  * [14] D. J. Montana, “Contact Stability for Two-Fingered Grasps,” _IEEE Transactions on Robotics and Automation_ , vol. 8, no. 4, pp. 421–430, 1992.
  * [15] R. M. Murray, S. S. Sastry, and L. Zexiang, _A Mathematical Introduction to Robotic Manipulation_. USA: CRC Press, Inc., 1994.
  * [16] K. Hang, W. G. Bircher, A. S. Morgan, and A. M. Dollar, “Manipulation for self-identification, and self-identification for better manipulation,” _Science Robotics_ , vol. 6, no. 54, 2021.
  * [17] A. M. Dollar and R. D. Howe, “The Highly Adaptive SDM Hand: Design and Performance Evaluation,” _The International Journal of Robotics Research_ , vol. 29, no. 5, pp. 585–597, apr 2010.
  * [18] K. Hang, A. S. Morgan, and A. M. Dollar, “Pre-grasp sliding manipulation of thin objects using soft, compliant, or underactuated hands,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 662–669, 2019.
  * [19] A. Morgan, B. Wen, J. Liang, A. Boularias, A. Dollar, and K. Bekris, “Vision-driven compliant manipulation for reliable; high-precision assembly tasks,” _Robotics: Science and Systems XVII_ , Jul 2021.
  * [20] R. H. N. Chavan-Dafle and A. Rodriguez, “Planar in-hand manipulation via motion cones,” _IJRR_ , 2019.
  * [21] S. Yuan, L. Shao, C. L. Yako, A. Gruebele, and J. K. Salisbury, “Design and control of roller grasper v2 for in-hand manipulation,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020, pp. 9151–9158.
  * [22] D. L. Brock, “Enhancing the dexterity of a robot hand using controlled slip,” in _Proceedings. 1988 IEEE International Conference on Robotics and Automation_. IEEE, 1988, pp. 249–251.
  * [23] Z. Su, K. Hausman, Y. Chebotar, A. Molchanov, G. E. Loeb, G. S. Sukhatme, and S. Schaal, “Force estimation and slip detection/classification for grip control using a biomimetic tactile sensor,” in _2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids)_. IEEE, 2015, pp. 297–303.
  * [24] W. G. Bircher, A. S. Morgan, and A. M. Dollar, “Complex manipulation with a simple robotic hand through contact breaking and caging,” _Science Robotics_ , vol. 6, no. 54, 2021. [Online]. Available: https://robotics.sciencemag.org/content/6/54/eabd2666
  * [25] K. Hang, M. Li, J. A. Stork, Y. Bekiroglu, F. T. Pokorny, A. Billard, and D. Kragic, “Hierarchical fingertip space: A unified framework for grasp planning and in-hand grasp adaptation,” _IEEE Transactions on robotics_ , vol. 32, no. 4, pp. 960–972, 2016.
  * [26] K. Hauser and J.-C. Latombe, “Multi-modal motion planning in non-expansive spaces,” _The International Journal of Robotics Research_ , vol. 29, no. 7, pp. 897–915, 2010.
  * [27] M. Zucker, N. Ratliff, A. D. Dragan, M. Pivtoraiko, M. Klingensmith, C. M. Dellin, J. A. Bagnell, and S. S. Srinivasa, “Chomp: Covariant hamiltonian optimization for motion planning,” _The International Journal of Robotics Research_ , vol. 32, no. 9-10, pp. 1164–1193, 2013.
  * [28] Z. Kingston, A. M. Wells, M. Moll, and L. E. Kavraki, “Informing multi-modal planning with synergistic discrete leads,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 3199–3205.
  * [29] J. J. Craig, _Introduction to robotics: mechanics and control, 3/E_. Pearson Education India, 2009.
  * [30] D. Huynh, “Metrics for 3d rotations: Comparison and analysis,” _Journal of Mathematical Imaging and Vision_ , vol. 35, pp. 155–164, 10 2009.
  * [31] B. Calli and A. M. Dollar, “Vision-based model predictive control for within-hand precision manipulation with underactuated grippers,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 2839–2845.
  * [32] A. S. Morgan, K. Hang, and A. M. Dollar, “Object-agnostic dexterous manipulation of partially constrained trajectories,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5494–5501, 2020.
  * [33] S. Cruciani, B. Sundaralingam, K. Hang, V. Kumar, T. Hermans, and D. Kragic, “Benchmarking in-hand manipulation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 588–595, 2020.
  * [34] B. Wen and K. Bekris, “Bundletrack: 6d pose tracking for novel objects without instance or category-level 3d models,” _IROS_ , 2021.