# Online Time-Informed Kinodynamic Motion Planning of Nonlinear Systems

Fei Meng, Jianbang Liu, Haojie Shi, Han Ma, Hongliang Ren, , and
Max Q.-H. Meng This work was supported in part by Hong Kong RGC CRF grant
C4063-18G and National Natural Science Foundation of China grant #62103181.
(Corresponding authors: Hongliang Ren and Max Q.-H. Meng)Fei Meng, Jianbang
Liu, Haojie Shi, Han Ma, and Hongliang Ren are with the Department of
Electronic Engineering, The Chinese University of Hong Kong, Hong Kong,
{feimeng, henryliu, h.shi<EMAIL_ADDRESS><EMAIL_ADDRESS>Q.-H.
Meng is with Shenzhen Key Laboratory of Robotics Perception and Intelligence
and the Department of Electronic and Electrical Engineering at Southern
University of Science and Technology in Shenzhen, China, and also a Professor
Emeritus in the Department of Electronic Engineering at The Chinese University
of Hong Kong, Hong Kong<EMAIL_ADDRESS>
###### Abstract

Sampling-based kinodynamic motion planners (SKMPs) are powerful in finding
collision-free trajectories for high-dimensional systems under differential
constraints. Time-informed set (TIS) can provide the heuristic search domain
to accelerate their convergence to the time-optimal solution. However,
existing TIS approximation methods suffer from the curse of dimensionality,
computational burden, and limited system applicable scope, e.g., linear and
polynomial nonlinear systems. To overcome these problems, we propose a method
by leveraging deep learning technology, Koopman operator theory, and random
set theory. Specifically, we propose a Deep Invertible Koopman operator with
control U model named DIKU to predict states forward and backward over a long
horizon by modifying the auxiliary network with an invertible neural network.
A sampling-based approach, ASKU, performing reachability analysis for the DIKU
is developed to approximate the TIS of nonlinear control systems online.
Furthermore, we design an online time-informed SKMP using a direct sampling
technique to draw uniform random samples in the TIS. Simulation experiment
results demonstrate that our method outperforms other existing works,
approximating TIS in near real-time and achieving superior planning
performance in several time-optimal kinodynamic motion planning problems.

##  I INTRODUCTION

Sampling-based motion planning methods are known for rapidly finding
collision-free “geometric” paths from a start to a goal point in high-
dimensional C-space. The popular RRT [1] grows a tree structure by drawing
random samples until reaching the goal, guaranteeing probabilistic
completeness. Its optimal version, RRT∗ [2], provides the property of
asymptotic optimality, whereas it costs too much time due to a uniform
sampling strategy in the whole C-space [3, 4]. To speed up the convergence to
the optimal solution, some heuristic sampling approaches have been proposed;
for example, “$L_{2}$ Informed Set” containing all the potential shorter
solutions is constructed to restrict the subsequent searches after an initial
path is found [5].

Figure 1: An illustration of online time-informed kinodynamic motion planning
of nonlinear control systems. (a) A Deep Invertible Koopman operator with
control U model, DIKU, is trained offline for the nonlinear systems to obtain
equivalent linear systems that enable forward and backward dynamics prediction
in the lifted space. (b) Our algorithm randomly samples states in the start
and goal sets. It then uses the ASKU to bidirectionally propagate the learned
linear dynamics and then recover, generating the forward reachable set
$\mathcal{X}_{k}^{f}$ (independent blue dashed convex hull) and backward
reachable tube $\mathcal{X}^{b}_{-cost:0}$ (serials of green dashed convex
hulls) in near real-time. Their intersections constitute the time-informed set
(TIS) $\Omega(cost)$ (black convex polygons). Time-informed tree growth is
achieved for the off-the-shelf SKMPs by directly sampling in the TIS. The TIS
updates according to the cost of the search tree returned after sufficient
iterations to help refine the solution (red line).

When kinodynamic constraints, i.e., kinematic and differential constraints,
are considered, the problem becomes PSPACE-hard and imposes a heavy
computational burden. Optimal trajectories under kinodynamic constraints,
instead of line segments, connecting any two intermediate tree nodes are
required to form the edges of the tree. Local trajectory optimization solvers
[6, 7] and random shooting methods [8, 1] are widely implemented to solve this
two-point boundary value problem (TPBVP). Stable Sparse RRT (SST) [8] performs
forward dynamic simulation with randomly sampled control signals for steering,
eliminating the need for TPBVP solvers. Nevertheless, it is computation-
demanding since its unbiased exploration of the whole state space.
Hierarchical rejection sampling (HRS) [9] and Markov Chain Monte-Carlo (MCMC)
[10] algorithms were proposed to generate samples within the informed set of
non-Euclidean state spaces. However, the HRS method takes a longer time as the
interested subspace shrinks or in high-dimensional space. Although the Hit-
and-Run (HNR) MCMC method [10] is efficient, it requires an initial solution
and applies to only certain systems for minimum-time problems.

Reachability-guided methods have been proposed to reduce computation cost
[11]. At each extension, the polytopic forward reachable set (FRS) of each
state is approximated to avoid selecting the nearest but unreachable state
[12]. In addition, the TIS that includes all feasible kinodynamic trajectories
is provided to restrict the sampling domain [13]. In principle, a set
operation of the FRS of the start and the backward reachable tube (BRT) of the
goal state is conducted. Tang et al. [14], which is most relevant to ours,
relaxed the Hamilton-Jacobi-Bellman (HJB) equation to compute a modified TIS,
broadening the system scope from linear [13] to polynomial nonlinear (PNL)
systems and conducting a heuristic search from the beginning. However, it
suffers from the curse of dimensionality, a heavy computational burden, and
limited system-specific applicable scope. The sampling-based reachability
analysis method that takes convex hulls for propagated states with random
control inputs from the start has great potential to solve these problems by
leveraging random set theory [15]. However, only FRSs can be approximated
since there is no generic differential equation for the backward propagation
of nonlinear systems. Considering Deep Koopman Operator (DKO) can predict the
forward evolution of nonlinear dynamics by mapping the original system into an
equivalent linear one in an observation space via deep neural network (DNN)
[16] and linear systems are invertible, we integrate it into the sampling-
based method to pave the path of constructing BRT. Furthermore, an adversarial
sampling manner can guarantee the over-approximation reachability analysis,
benefiting informed sampling. Although reachable tubes for high-dimensional
nonlinear systems can be approximated via DNN end-to-end, it takes long
computation times from hours to days [17]. By contrast, our method only needs
a short time to learn the locomotion characteristics of the nonlinear systems
to realize reachability analysis indirectly.

In this letter, we leverage the deep learning, KO and random set theories to
contribute a generic and efficient time-informed SKMP. The paradigm of the
planning procedure is shown in Fig. 1. The main contributions are as
follows111Source code: https://github.com/feimeng93/OnlineTimeInformedKinoMP.

  * •

A deep invertible Koopman operator with control U, named DIKU, that can
bidirectionally predict the long-term evolution of nonlinear control systems
is proposed by integrating the invertible neural networks (INN) [18] into our
previous DKO model [16].

  * •

TIS is approximated online for the nonlinear control systems that can be
learned. We use a sampling-based method [19] to quickly perform the
reachability analysis based on the bidirectional propagation of DIKU,
overcoming the curse of dimensionality, broadening the system’s applicable
scope and boosting the computation efficiency for the TIS significantly.

  * •

We develop an online time-informed kinodynamic motion planning algorithm where
direct sampling in the TIS is realized. We evaluate our method on six
different types of dynamic systems and achieve better planning efficiency over
the existing works.

##  II Problem Formulation and Preliminaries

###  II-A Time-optimal Kinodynamic Motion Planning

Let compact sets $\mathcal{X}\subset\mathbb{R}^{n}\,(n\geq 1)$ and
$\mathcal{U}\subset\mathbb{R}^{m}\,(m\geq 1)$ be state and admissible control
spaces, respectively. Denote $\mathcal{X}_{\text{obs}}\subseteq\mathcal{X}$ as
the obstacle space and the closure of the set difference
$\mathcal{X}_{\text{free}}=\text{cl}(\mathcal{X}\backslash\mathcal{X}_{\text{obs}})$
as the free space. Let $\bm{x}_{0}\in\mathcal{X}_{\text{free}}$ be the start
state and $\mathcal{X}^{g}\subset\mathcal{X}_{\text{free}}$ be the goal set.
Let us define a dynamical system as follows.

$\dot{\bm{x}}(t)=f(\bm{x}(t),\bm{u}(t)),$ (1)

where $t$ is the sampling time, $\bm{x}(t)\in\mathcal{X}$ is the state,
$\bm{u}(t)\in\mathcal{U}$ is the control signals, and the flow field
$f(\cdot):\mathcal{X}\times\mathcal{U}\rightarrow\mathcal{X}$ is continuously
differentiable.

Given an initial state $\bm{x}_{0}\in\mathcal{X}_{\text{free}}$ and a target
set $\mathcal{X}^{g}\subset\mathcal{X}_{\text{free}}$, we aim to determine the
optimal control signals $\bm{u}$ and minimum time-to-reach $t_{f}$ such that a
time-optimal trajectory $[\bm{x}(0),\dots,\bm{x}(t_{f})]$ connecting the start
and a goal state $\bm{x}(t_{f})\in\mathcal{X}^{g}$ is found for (1).
Mathematically, we have the time-optimal kinodynamic motion planning problem:

$\displaystyle\hskip 71.13188pt\mathop{\min}_{\bm{u}}t_{f}$ (2) s.t.
$\displaystyle\dot{\bm{x}}(t)=f(\bm{x}(t),\bm{u}(t))$
$\displaystyle\bm{x}(0)=\bm{x}_{0},\bm{x}(t_{f})\in\mathcal{X}^{g},t_{f}>0,\text{and}\;t_{f}\;\text{is
free}$
$\displaystyle\bm{x}(t)\in\mathcal{X}_{\text{free}},\,\bm{u}(t)\in\mathcal{U},\forall
t\in[0,t_{f}].$

###  II-B Time-Informed Set

The forward reachable set (FRS) is defined as the set of all possible states
that can be reached at time $t$, starting from $\bm{x}_{0}$ at time $0$, using
admissible controls,

$\displaystyle\mathcal{X}^{f}(0,t,\bm{x}_{0})\triangleq\left\\{\bm{x}_{t}\in\mathcal{X}\;|\;\forall\tau\in[0,t],\exists\bm{u}(\tau)\in\mathcal{U}\right.$
(3)
$\displaystyle\text{s.t.}\,\left.\bm{x}(0)=\bm{x}_{0},\bm{x}(t)=\bm{x}_{t},\dot{\bm{x}}(\tau)=f(\bm{x}(\tau),\bm{u}(\tau))\right\\}$

The backward reachable set (BRS) for the target set $(t_{f},\mathcal{X}^{g})$
is the set of all states at time $t$ that will reach the goal at $t_{f}$,

$\displaystyle\mathcal{X}^{b}(t,t_{f},\mathcal{X}^{g})\triangleq\left\\{\bm{x}_{t}\in\mathcal{X}\;|\;\forall\tau\in[t,t_{f}],\exists\bm{u}(\tau)\in\mathcal{U}\right.$
(4)
$\displaystyle\text{s.t.}\,\left.\bm{x}(t)=x_{t},\bm{x}(t_{f})\in\mathcal{X}^{g},\dot{\bm{x}}(\tau)=f(\bm{x}(\tau),\bm{u}(\tau))\right\\}.$

Let $\mathcal{F}(0,t,\bm{x}_{0})$ and $\mathcal{B}(t,t_{f},\mathcal{X}^{g})$
be the over-approximations of FRS and BRS, respectively. The BRT that contains
the set of all states $\bm{x}_{t}$ that can reach $\mathcal{X}^{g}$ at
$\forall\tau\in[t,t_{f}]$ is
$\mathcal{R}_{b}(t,t_{f},\mathcal{X}^{g})=\cup_{t\leq\tau\leq
t_{f}}\mathcal{B}(t,\tau,\mathcal{X}^{g})$.

The time-informed set (TIS) of time cost $t_{f}$ is [13]:

$\Omega(t_{f})\triangleq\mathop{\bigcup}\limits_{t\in[0,t_{f}]}\mathcal{F}(0,t,\bm{x}_{0})\cap\mathcal{R}_{b}(t,t_{f},\mathcal{X}^{g})\subset\mathcal{X}.$
(5)

###  II-C Koopman Invertible Autoencoder

A Koopman invertible autoencoder (KIA) is proposed by replacing the linear
auxiliary network in the classic DKO with an INN architecture [18] to model
both forward and backward dynamics [20]. INN establishes a connection between
the dual dynamics according to the bijective functions:

$\bm{g}_{t+1,1}=\bm{g}_{t,1}+f_{2}(\bm{g}_{t,2}),\;\bm{g}_{t+1,2}=\bm{g}_{t,2}+f_{1}(\bm{g}_{t+1,1});$
(6)
$\bm{g}_{t-1,2}=\bm{g}_{t,2}-f_{1}(\bm{g}_{t,1}),\;\bm{g}_{t-1,1}=\bm{g}_{t,1}-f_{2}(\bm{g}_{t-1,2});$
(7)

where $\bm{g}_{t}=[\bm{g}_{t,1};\,\bm{g}_{t,2}]$, $[;]$ is the vertical
concatenation operation, and $f_{1}(\cdot)$ and $f_{2}(\cdot)$ are translation
functions. The original lifted forward evaluation with KO $K$,
$\bm{g}(\bm{x}_{t+1})=K\bm{g}(\bm{x}_{t}),\;\bm{g}(\bm{x})=[\psi_{1}(\bm{x}),\psi_{2}(\bm{x}),...,\psi_{c}(\bm{x})]^{T}\in\mathbb{R}^{c}\subset\mathcal{C}$,
can be approximated with (6), and its inverse expression is obtained by (7) at
no cost. Although KIA can forecast and backcast temporal data, its evolution
requires an NN decoder and does not include any control input [20].

##  III Deep Koopman Reachability Analysis

###  III-A Deep Invertible Koopman Operator with Control Enabling Long-Term
Forward and Backward Propagation

It has been proven that the DNN autoencoder can discover the function
observables with linear evolution for controlled lifted linear system [16]. We
parameterize the observable $\bm{z}(\bm{x}_{t})$ as
$\bm{z}_{\theta_{z}}(\bm{x}_{t})$ by using a DNN with parameters $\theta_{z}$,
and concatenate the state and embedded feature together as the observable as
$\bm{z}_{t}=[\bm{x}_{t};\bm{z}_{\theta_{z}}(\bm{x}_{t})]\in\mathbb{R}^{n+d}$
to realize authentic recovery to nonlinear states. Furthermore, we have the
following $K$-steps feedforward prediction using NNs:

$\displaystyle\bm{z}_{t+k+1}=K_{x}\bm{z}_{t+k}$
$\displaystyle+K_{u}\bm{u}_{t+k},\;k=0,\cdots,K-1$ (8)
$\displaystyle\bm{x}_{t+k}$ $\displaystyle=C\bm{z}_{t+k},$

where $K_{x}\in\mathbb{R}^{c_{x}\times c_{x}}$ and
$K_{u}\in\mathbb{R}^{c_{x}\times c_{u}},\,c=c_{x}+c_{u},$ are fully connected
layers with parameters $\theta_{x}$ and $\theta_{u}$, respectively. Note that
nonlinear state constraints $c(\bm{x}_{t})\leq
0:\mathbb{R}^{n}\rightarrow\mathbb{R}$ can be considered by adding them to
$\bm{z}_{t}$. For instance, if there exists a constraint $\sin(x_{1})\leq 0$ ,
we can let $z_{t}^{(n+d+1)}=sin(x_{1})$ and enforce it no greater than 0 [21].
Observing the linearity in the lifted space, we could have backcasted via the
inverse matrix of $K_{x}$. However, it is computationally expensive, typically
$O(n^{3})$.

Herein, we propose the DIKU utilizing the invertible INNs to realize
bidirectional dynamics, whose architecture is shown in Fig. 2. The latent
features are divided into two parts as
$\bm{z}_{t}=[\bm{z}_{{t},1};\,\bm{z}_{{t},2}]$ and so does the control field
$K_{u}u_{t}=[B_{1};B_{2}]u_{t}=\hat{\bm{u}}_{t}=[\hat{\bm{u}}_{t,1};\,\hat{\bm{u}}_{t,2}]$.
Then, we parameterize the $f_{1}(\cdot)$ and $f_{2}(\cdot)$ by utilizing DNNs
with parameters $\theta_{1}$ and $\theta_{2}$ followed by activation
functions, and denote them as $A_{1}$ and $A_{2}$. Consequently, (8) is
rewritten as follows to improve long-term prediction ability:

$\displaystyle\bm{z}_{{t+k+1},1}$
$\displaystyle=\bm{z}_{{t+k},1}+A_{2}\bm{z}_{{t+k},2}+\hat{\bm{u}}_{{t+k},1},$
(9) $\displaystyle\bm{z}_{{t+k+1},2}$
$\displaystyle=\bm{z}_{{t+k},2}+A_{1}\bm{z}_{{t+k+1},1}+\hat{\bm{u}}_{{t+k},2}.$

The corresponding backward evolution in the lifted space is derived for free:

$\displaystyle\bm{z}_{{t-k-1},2}$
$\displaystyle=\bm{z}_{{t-k},2}-A_{1}\bm{z}_{{t-k},1}-\hat{\bm{u}}_{{t-k},2},$
(10) $\displaystyle\bm{z}_{{t-k-1},1}$
$\displaystyle=\bm{z}_{{t-k},1}-A_{2}\bm{z}_{{t-k-1},2}-\hat{\bm{u}}_{{t-k},1}.$

Figure 2: An overview of Deep Invertible Koopman operator with control U
(DIKU) neural network model for long-horizon forward and backward dynamics
prediction

We collect the training trajectory dataset as
$\mathcal{D}=\left[X_{i}\in\mathbb{R}^{N\times n},U_{i}\in\mathbb{R}^{N\times
m},i=1,\cdots,K\right]$. We define the following $K$-steps forward and
backward dynamical loss as:

$\displaystyle L(\theta_{z},\,\theta_{1},\theta_{2},\,\theta_{u})$
$\displaystyle=\frac{1-\gamma^{K}}{1-\gamma}\sum_{i=1}^{K}\gamma^{i-1}[MSE(X_{i},CZ_{i})$
(11) $\displaystyle+MSE(X_{K-i+1},CZ_{K-i+1})]$

where $MSE$ represents the mean square loss function and $\gamma$ denotes the
weight decay constant. This $K$-steps loss is conductive to fine predictions
in the long term.

###  III-B Online Approximate Time-Informed Set via Random Set and Koopman
Operator Theories

Random set theory [15] guarantees asymptotic convergence to the convex hull of
the reachable sets, whose relevant theorem is restated as follows:

###### Lemma III.1

[19] Let $\left\\{(\bm{x}_{0}^{j},\bm{u}^{j})\right\\}_{j=1}^{M}$ be i.i.d.
sampled parameters in $\mathcal{X}_{0}\times\mathcal{U}^{k-1}$, where $M$
denotes the number of states. Define $\bm{x}_{t}^{j}$ as (1). Assume that the
sampling distribution of the parameters satisfies
$\mathbb{P}(\bm{x}_{t}^{j}\in G_{t})>0$ for any open set $G_{t}$ s.t.
$\mathcal{X}_{t}\cap G_{t}\neq\emptyset$. Then, as $M\rightarrow\infty$,
$\mathcal{X}_{t}$ converges to the convex hull of the RS
$\text{Co}(\mathcal{X}_{t})$ almost surely.

A random set is defined as a map from a probability space to a family of sets.
Consider sampling i.i.d. initial tuples
$(\bm{x}_{0}^{j},\bm{u}_{0}^{j})\in\mathcal{X}^{f}_{0}\times\mathcal{U}$ and
goal tuples
$({\bm{x}}_{g}^{j},{\bm{u}}_{g}^{j})\in{\mathcal{X}}^{b}_{0}\times\mathcal{U}$
according to arbitrary probability distribution, where $j=1,\cdots,M$. Then,
we encode the sampled states, obtaining the tuples in the lifted space, i.e.,
$(\bm{z}_{0}^{j},\bm{u}_{0}^{j})\in\mathcal{Z}^{f}_{0}\times\mathcal{U}$ and
$({\bm{z}}_{g}^{j},{\bm{u}}_{g}^{j})\in{{{\mathcal{Z}}}}^{b}_{0}\times\mathcal{U}$.
Next, we define a time step size $\delta$, and forecast the initial
observables with the consistent $\bm{u}_{0}^{j}$ by (9) and backcast the goal
observables with the ${\bm{u}}_{g}^{j}$ through (10) over $T$ time period
simultaneously. We recover the predicted states from the generated
observables. The convex hulls of the resulting states forward and backward
propagation are denoted as $\text{Co}(\mathcal{X}^{f}_{0:T})$ and
$\text{Co}({\mathcal{X}}^{b}_{-T:0})$, respectively.

In most cases, the convex polytopic sets $\mathcal{X}^{f}_{0}$ and
${\mathcal{X}}^{b}_{0}$ are mapped to nonconvex sets $\mathcal{Z}^{f}_{0}$ and
${{\mathcal{Z}}}^{b}_{0}$ since $\bm{z}(\bm{x})$ are nonlinear functions of
$\bm{x}$. Our convex hulls will be over-approximations, sufficient for a
heuristic sampling. To avoid returning a subset of the true and convex
reachable set, we apply the adversarial sampling [19] to inflate the resulting
convex hulls. To over-approximate FRS, after obtaining the set of the total
$T$ period forward propagated trajectories $\mathcal{X}^{f}_{0:T}=\cup_{0\leq
t\leq T}\\{\bm{x}_{t}^{j}\\}_{j=1}^{M}$, we use projected gradient ascent to
maximize
$\mathcal{L}^{M}(\bm{p})=\frac{1}{T}\sum_{t=\delta}^{T}\lVert\bm{x}_{t}(\bm{p}_{0})-\bm{c}_{t}^{M}\rVert^{2}_{Q_{t}^{M}}$
to regenerate initial tuples $\bm{p}_{0}=(\bm{x}_{0}^{j},\bm{u}_{0}^{j})$,
where $Q_{t}^{M}$ are the positive definite matrices and $\bm{c}_{t}^{M}$
represent the geometric center of $\mathcal{X}^{f}_{t}$.
$\bm{p}_{0}\leftarrow\bm{p}_{0}+\eta\nabla_{p}\mathcal{L}^{M}(\bm{x}_{\delta:T}^{j})$,
where $\eta$ is a constant. After projecting the new $p_{0}$ onto
$\mathcal{X}^{f}_{0}\times\mathcal{U}$, we propagate the new samples again
using our DIKU and recover the trajectories to obtain
$\\{\bm{x}_{0:T}^{j}\\}_{j=1}^{M}$ that inflates
$\text{Co}(\mathcal{X}^{f}_{0:T})$. The propagated trajectories are updated as
$\mathcal{X}^{f}_{0:T}\leftarrow\mathcal{X}^{f}_{0:T}\cup\\{\bm{x}_{0:T}^{j}\\}_{j=1}^{M}$.
According to Lemma III.1, the over-approximation of FRS at time $t$ is
$\mathcal{F}(0,t,\bm{x}_{0})=\text{Co}(\mathcal{X}^{f}_{t})$.

Similarly, we resample the goal tuples $p_{g}=(\bm{x}_{g}^{j},\bm{u}_{g}^{j})$
and backcast the $p_{g}$ by utilizing our backward DIKU and retrieve
$\\{{\bm{x}}_{-T:0}^{j}\\}_{j=1}^{M}$. The adversarial sampling procedure of
backward trajectories is
${\mathcal{X}}^{b}_{-T:0}\leftarrow{\mathcal{X}}^{b}_{-T:0}\cup\\{{\bm{x}}_{-T:0}^{j}\\}_{j=1}^{M}$.
Given a terminal time instant $t_{f}$, the BRT over the time interval
$[t,t_{f}]$ is over-approximated as
$\mathcal{R}_{b}(t,t_{f},\mathcal{X}^{g})=\text{Co}({\mathcal{X}}^{b}_{t-t_{f}:0})$.

TIS is finally over-approximated with adequate number of samples $M$ as below
according to (5)

$\Omega(t_{f})=\\!\\!\\!\mathop{\bigcup}\limits_{0\leq t\leq
t_{f}}\\!\\!\\!\text{Co}(\mathcal{X}^{f}_{t})\cap\text{Co}({\mathcal{X}}^{b}_{t-t_{f}:0}).$
(12)

Although we can encode the states in the sets $\mathcal{X}^{f}_{0}$ and
${\mathcal{X}}^{b}_{0}$, it is intractable to embed the sets that have
arbitrary shape representations into $\mathcal{Z}^{f}_{0}$ and
${{\mathcal{Z}}}^{b}_{0}$ [22]. Besides, inflation errors would be added in
the process due to the extra $\bm{z}_{\theta_{z}}(\bm{x})$ even though we
predefine the shapes of $\mathcal{Z}^{f}_{0}$ and
${{\mathcal{Z}}}^{b}_{t_{f}}$ based on experience. Thus, we only propagate the
state but do not sample adversarially in the lifted space. System dynamics can
be used directly with the sampling-based technique [19] to obtain a better
accuracy while the forward DIKU, along with the technique, guarantees online
FRS approximation, even possible for some unknown systems.

##  IV Online Time-Informed SKMP of Nonlinear Systems

Input: $(\bm{x}_{0},\mathcal{X}^{g});\;Q=\\{\bm{x}_{0}\\};\;E=\emptyset$

Output: $Tree(Q,E)$

1 $\tau=0;\;\mathcal{X}^{b}_{0}=\mathcal{X}^{g};\;M=1000;\;\eta;\;Tree(Q,E);$

2 while  _$\bm{x}_{0}\notin\rm{Co}(\mathcal{X}^{\textit{b}}_{\tau})$_ do

3 $\tau=\tau-\delta$;

4$cost=-\tau$;

5 $\Omega(cost)=$ASKU$(cost,\bm{x}_{0},\mathcal{X}^{g},\eta,M)$;
$\triangleright$ Sec. III-B

6 $cost^{\prime}=\infty$;

7 while  _$|(cost-cost^{\prime})/cost^{\prime}| >eps$_ do

8 $cost^{\prime}=cost$;
$[Tree,cost_{new}]=$TreeGrowth$(Tree,cost,\Omega(cost))$; if
_$cost_{new}\neq\infty$_ then

9 $cost=cost_{new}$;

10 if  _$cost_{new} <cost$_ then

11
$\Omega(cost)=\Omega(cost^{\prime})\backslash(\text{Co}(\mathcal{X}^{f}_{cost^{\prime}})\cap\text{Co}({\mathcal{X}}^{b}_{0}))\backslash(\bigcup\limits_{0\leq
t\leq
cost+\delta}\text{Co}(\mathcal{X}^{f}_{t})\cap\text{Co}({\mathcal{X}}^{b}_{t-cost^{\prime}:t-cost-\delta}))$
; $Tree=$Prune$(\Omega(cost),Tree,cost)$;

12

13 else

14 $cost=cost+\delta_{2}$;
$\Omega(cost)=\Omega(cost^{\prime})\cup(\text{Co}(\mathcal{X}^{f}_{cost})\cap\text{Co}({\mathcal{X}}^{b}_{0}))\cup(\bigcup\limits_{0\leq
t\leq
cost-\delta}\text{Co}(\mathcal{X}^{f}_{t})\cap\text{Co}({\mathcal{X}}^{b}_{t-cost:t-cost^{\prime}-\delta}))$
;

15

Algorithm 1 Online Time-Informed Sampling-based Kinodynamic Motion Planning
Figure 3: Comparison of forward and backward dynamics prediction by our DIKU
and the DKU with consistency loss. The early predictions are more accurate
than the late ones, which can be observed from the lines with plus markers of
forward prediction from left to right and from those lines with dot markers of
backward from right to left. The lower blue lines indicate smaller
bidirectional dynamics evolution errors by DIKU, especially for nonlinear
systems.

Algorithm $1$ outlines the procedure of our planning method. We obtain the
optimal time cost ignoring obstacles by determining the index of BRS that
contains the start state $\bm{x}_{0}$. The TIS corresponding to this $cost$ is
constructed using our ASKU approach in Section III-B. Then, we grow a search
tree within $\Omega(cost)$ and compute the tree’s time cost $cost_{new}$ while
considering obstacles. $\Omega(cost)$ reduces and we eliminate the nodes
outside it once we find a better suboptimal solution. In contrast, we extend
the optimal time of the planning problem $\delta_{2}$s and expand the TIS
since a collision-free solution less than $cost$ might not exist.

Algorithm $2$ describes the procedure to grow a search tree with a direct
sampling technique, HNR sampling method [23], within the components of TIS. A
time $t$ is first uniformly sampled in $[0,cost]$ as the cost-to-come of the
sample to be generated. If a random number is generated below the threshold
$\mu$, we draw a sample within the TIS to accelerate the convergence speed;
otherwise, within the entire state space. We compute the longest time cost
from $t$ to reach the goal. A linear programming problem, subject to the
linear inequalities of the FRS at $t$ and BRS at $t^{\prime}$, is formulated
to check if there exists any intersection. We decrease the cost-to-go until
finding the minimum value. A valid cost-to-go $\bar{t}$ is then uniformly
sampled. Finally, we use the vertex inclusion algorithm [13] and an off-the-
shelf SKMP to find a solution tree and return its cost $cost_{new}$, if there
exists.

##  V Simulation Experiments

In this section, we conduct simulation experiments to 1) verify our DIKU’s
bidirectional long-term prediction performance, 2) demonstrate the high
performance of our ASKU’s reachability analysis for TIS, and 3) show the
improved planning efficiency of our online time-informed SKMP.

Input: $Tree,cost,\Omega(cost)$

Output: $Tree,cost_{new}$

1 for  _$i=1:N$_ do

2 $t=\text{UniformSample}([0,cost])$; $\triangleright$ cost-to-come

3 if  _$\rm{Rand}(0,1) <\mu$_ then

4 $t^{\prime}=t-cost$; $\triangleright$ the maximum cost-to-go while
_$\\{\bm{x}\\!\in\\!\mathcal{X}|\rm{Co}(\mathcal{X}^{\textit{f}}_{\textit{t}})\\!\cap\\!\rm{Co}(\mathcal{X}^{\textit{b}}_{\textit{t}^{\prime}}),\textit{t}^{\prime}\\!\leq\\!0\\}\\!\neq\\!\emptyset$_
do

5 $t^{\prime}=t^{\prime}+\delta$; $\triangleright$ find the minimum cost-to-go

6 $\bar{t}=\text{UniformSample}([t-cost,t^{\prime}))$;
$\bm{x}_{new}=$HNR$(\\{\bm{x}\in\mathcal{X}\,|\,\text{Co}(\mathcal{X}^{f}_{t})\cap\text{Co}(\mathcal{X}^{b}_{\bar{t}})\\})$;

7 else

8 $\bm{x}_{new}=$Sample$(\mathcal{X})$;

9 $[Tree,cost_{new}]=\text{SKMP}(Tree,\bm{x}_{new},t)$;

10if  _$\bf{no}$ node $\in\mathcal{X}^{g}$_ then

11 $cost_{new}=\infty$;

Algorithm 2 Tree Growth

###  V-A Experiment Environments

We have six dynamic systems for evaluation, including a (a) 2D linear system
in [14] (2D-L), (b) 3D polynomial nonlinear system in [14] (3D-PNL), (c)
CartPole in [16], (d) DampingPendulum in [16], (e) two-link robot in [16], and
(f) planar quadrotor in [24]. For each system, a total of 10000, 2500, and
5000 trajectories are generated by simulating the dynamics, designated for
training, validation, and testing, respectively. Our experiments were
conducted in Python and C++ using a 64-bit workstation with an Intel Core
i9-13900K processor, NVIDIA 4090 GPU, and 128GB RAM, running Ubuntu 20.04 OS.

###  V-B Long-Term Forward and Backward Dynamics Prediction

We evaluate the long horizon forward and backward dynamical prediction
performance of our DIKU model on the six dynamics and compare them with the
existing method. We extend the Deep KoopmanU with Control (DKU) in [16] to
compute the backward counterpart of (8) through additionally training an
independent linear NN of $K_{x}^{-1}$. We impose the weights of parameterized
$K_{x}$ and $K_{x}^{-1}$ in a consistency loss as in [25]. One INN module is
constructed, where $\bm{z}_{t,2}$ is randomly selected. The hidden layers of
$A_{1}$ and $A_{2}$ are $[128,64,32]$ with the ReLU [26] in the dynamics
(a)-(e), while with linear activation functions in (f). The structures and
weights of embedding and control nets, training and testing datasets, learning
rate $0.001$, batch size $1024$, Adam optimizer [27], and K step loss with
$\gamma=0.9$ of the extended DKU and DIKU are identical for each system.

The 100-step bidirectional prediction results are shown in Fig. 3. We
calculate the mean of the log10 of the maximum error. It can be seen that our
DIKU (blue) can bidirectionally predict much more accurate states in the long
horizon than the DKU (yellow). Unlike the DKU which uses different structures,
parameters, and a soft loss function to train $K_{x}$ and $K_{x}^{-1}$ for
bidirectional prediction, our model allows direct back prediction through the
INN structure (10), benefiting a better inference and convergence in training.
Another reason is that the control field’s nonlinearity can be reflected by
the INN’s translation functions. Consequently, DIKU paves the way to
efficiently perform accurate long-term bidirectional propagation in sampling-
based reachability analysis.

###  V-C Online Time-Informed Set Approximation

We compare our basic sampling-based approach (ours) and its over-approximation
version, the ASKU, (ours(AS)) with level set method [28] (ground truth),
ellipsoidal toolbox [29] (ET), and relaxed HJB equation method [14] (RHJB) for
the dynamics (a). Set identical time step $0.05$s and horizon $T=5$s. We
evaluate the volume convergence of ours for computing BRT and show the results
in Fig. 5. We set $M=1000$, $\eta=0.015$, and $\eta=0.04$ for forward and
backward prediction to let the gradient and $\bm{p}_{0}$ have the same order
of magnitude. We use the ellipsoidal parameters of initial and goal sets and
control limits in [14]. The reachable sets and tubes that constitute TIS are
shown in Fig. LABEL:fig:tis. Since the ground truth has to discretize the
spaces to calculate the numerical result of the HJB partial differential
equation, we have tried to reduce the number of girds to save time while
guaranteeing the quality of the reachable sets. Table I records the
computation costs of reachable tubes. The ground truth takes plenty of time to
return a fine FRT. ET must be calculated in 30 random directions, resulting in
a heavy computation burden and a more compact ellipsoid than RHJB. RHJB
predefines the centers of the external ellipsoids instead of treating them as
decision variables since that makes solving SOS optimization intractable.
However, the centers highly influence the convergence of RHJB, leading to no
solution sometimes. Although sometimes feasible for the system, it is time-
consuming and may cause drift, as seen in the purple BRT in Fig.
LABEL:fig:tis. These methods also suffer from the curse of dimensionality. In
contrast, ours and ASKU can be completed in near real-time. Note that ASKU can
be applied to linear and nonlinear systems that can be learned regardless of
the number of dimensions, thanks to DIKU’s NN bidirectional propagation.
Facing the true convex reachable set, our ASKU can still return over-
approximation online.

Figure 5: Volume convergence (green line) and speed (blue line) on DIKU with a
varying number of samples $M$ for our basic BRT approximation. The volume
accuracy is computed by being divided by the volume of ground truth. A high
degree of accuracy can be obtained at the cost of computation speed.

###  V-D Online Time-Informed Kinodynamic Motion Planning

TABLE I: Computation times [s] of the Forward/Backward reachable tubes for the example 2D system in [14] | Level Set [28] | ET [29] | RHJB [14] | Ours | Ours (AS)
---|---|---|---|---|---
FRT | 47.46 | 51.13 | 22.17 | 0.03 | 0.16
BRT | 5.71 | 6.07 | 23.31 | 0.03 | 0.16

We conduct the problem (2) on all dynamics between SST [8] and our algorithm
while only comparing the approach in [14] with default parameters on (a) and
(b). We set $\delta_{2}=0.5$s, and $M=1000$ and valid $\eta$ for our method.
Set $\mu=0.9$ to accelerate our search progress and $N=500$ in [14] and ours.
A library of TIS computed offline exists only for [14]. The running time of
approximating TIS is not counted in [14] but it is included in ours because
ours is an online time-informed method. The computation time is measured as
the condition (line 7, Alg. 1) is no longer satisfied.

TABLE II: Comparison results of the running time, path qualities, and the number of nodes for benchmarks and ours |  | $T_{IN}[s]$ | $C_{IN}[s]$ | $T_{OP}[s]$ | $C_{OP}[s]$ | $N_{node}(\times 10^{2})$
---|---|---|---|---|---|---
2D-L | Ours | $\bm{0.87}$ | $\bm{4.68\\!\\!\pm\\!\\!0.10}$ | $\bm{6.87}$ | $\bm{4.44\\!\\!\pm\\!\\!0.09}$ | $\bm{3.21\\!\\!\pm\\!\\!0.7}$
[14] | $1.02$ | $4.78\\!\\!\pm\\!\\!0.12$ | $7.02$ | $4.46\\!\\!\pm\\!\\!0.12$ | $3.4\\!\\!\pm\\!\\!0.81$
SST | $1.78$ | $9.11\\!\\!\pm\\!\\!2.91$ | $25.78$ | $4.87\\!\\!\pm\\!\\!0.60$ | $12.9\\!\\!\pm\\!\\!0.29$
3D-PNL | Ours | $\bm{1.45}$ | $\bm{4.22\\!\\!\pm\\!\\!0.21}$ | $\bm{4.45}$ | $\bm{4.07\\!\\!\pm\\!\\!0.08}$ | $\bm{1.67\\!\\!\pm\\!\\!0.27}$
[14] | $1.68$ | $4.31\\!\\!\pm\\!\\!0.22$ | $4.68$ | $4.08\\!\\!\pm\\!\\!0.10$ | $1.71\\!\\!\pm\\!\\!0.28$
SST | $26.37$ | $6.38\\!\\!\pm\\!\\!2.08$ | $56.37$ | $5.33\\!\\!\pm\\!\\!1.40$ | $28.51\\!\\!\pm\\!\\!7.58$
CartPole | Ours | $\bm{5.26}$ | $\bm{4.91\\!\\!\pm\\!\\!1.74}$ | $\bm{11.26}$ | $\bm{4.18\\!\\!\pm\\!\\!1.73}$ | $\bm{5.87\\!\\!\pm\\!\\!1.91}$
SST | $14.69$ | $8.21\\!\\!\pm\\!\\!4.90$ | $29.69$ | $6.15\\!\\!\pm\\!\\!4.21$ | $19.44\\!\\!\pm\\!\\!5.59$
Damping | Ours | $\bm{6.23}$ | $\bm{2.08\\!\\!\pm\\!\\!0.35}$ | $\bm{10.57}$ | $\bm{1.85\\!\\!\pm\\!\\!0.39}$ | $\bm{5.26\\!\\!\pm\\!\\!1.78}$
Pendulum | SST | $21.2$ | $2.38\\!\\!\pm\\!\\!1.31$ | $33.2$ | $1.91\\!\\!\pm\\!\\!0.96$ | $21.25\\!\\!\pm\\!\\!6.51$
Two-link | Ours | $\bm{8.34}$ | $\bm{4.72\\!\\!\pm\\!\\!1.17}$ | $\bm{17.34}$ | $\bm{4.39\\!\\!\pm\\!\\!1.56}$ | $\bm{6.52\\!\\!\pm\\!\\!2.03}$
Robot | SST | $27.12$ | $7.10\\!\\!\pm\\!\\!5.20$ | $45.12$ | $4.60\\!\\!\pm\\!\\!4.60$ | $26.35\\!\\!\pm\\!\\!12.75$
Planar | Ours | $\bm{15.64}$ | $\bm{8.50\\!\\!\pm\\!\\!1.02}$ | $\bm{24.64}$ | $\bm{7.72\\!\\!\pm\\!\\!1.12}$ | $\bm{45.54\\!\\!\pm\\!\\!13.04}$
Quadrotor | SST | $39.91$ | $10.52\\!\\!\pm\\!\\!2.83$ | $57.91$ | $8.90\\!\\!\pm\\!\\!2.62$ | $234.03\\!\\!\pm\\!\\!52.21$

Fig. LABEL:fig:problem depicts the schematics of the planning problems and
Table II reports the statistical results of 100 trials. $T_{IN}$ and $T_{OP}$
represent the computation time of discovering the initial solution and final
suboptimal solution, respectively. $C_{IN}$ and $C_{OP}$ stand for the
solutions’ costs, respectively. $N_{node}$ denotes the total number of nodes
at the conclusion. Since our method attempts to grow the search tree within
$\Omega$ from the beginning, better initial and final collision-free solutions
are returned much faster than SST in all environments. It can be observed that
our algorithm finds better initial and final solutions more quickly compared
with [14] because we directly draw samples in a tighter TIS. It can also apply
to the nonlinear systems (c)-(f) other than the PNL system (b). The sampling
acceptance rate among the six systems is $82\%$ on average. In terms of memory
burdens, our algorithm and [14] require fewer nodes than SST in the
environments of (a) and (b) due to the elimination effect of the Prune()
function. In the other environments, our method retains the tree with
$19.5\%\sim 30.2\%$ nodes since we sample within the $\Omega(cost)$ that
contains the optimal solution.

##  VI CONCLUSION

In this work, we realize online TIS approximation and propose an online time-
informed SKMP of nonlinear systems. Our proposed Deep Invertible Koopman
operator with control U (DIKU) model accurately infers forward and backward
dynamics propagation over a long horizon by designing the INN-based auxiliary
network. The developed ASKU method can over-approximate TIS in convex sets for
a variety of systems at little cost through conducting adversarial sampling
for DIKU bidirectional propagation. The time-optimal SKMP that directly
samples in the TIS takes less planning time than the [14] algorithm although
counting the TIS computation costs. Our next research direction is to improve
our DIKU and realize bidirectional search tree growth.

## References

  * [1] Steven M LaValle and James J Kuffner Jr.  Randomized kinodynamic planning.  The international journal of robotics research, 20(5):378–400, 2001.
  * [2] Sertac Karaman and Emilio Frazzoli.  Sampling-based algorithms for optimal motion planning.  The international journal of robotics research, 30(7):846–894, 2011.
  * [3] Fei Meng, Liangliang Chen, Han Ma, Jiankun Wang, and Max Q-H Meng.  Nr-rrt: Neural risk-aware near-optimal path planning in uncertain nonconvex environments.  IEEE Transactions on Automation Science and Engineering, 2022.
  * [4] Fei Meng, Liangliang Chen, Han Ma, Jiankun Wang, and Max Q-H Meng.  Learning-based risk-bounded path planning under environmental uncertainty.  IEEE Transactions on Automation Science and Engineering, 2023.
  * [5] Jonathan D Gammell, Siddhartha S Srinivasa, and Timothy D Barfoot.  Informed rrt: Optimal sampling-based path planning focused via direct sampling of an admissible ellipsoidal heuristic.  In 2014 IEEE/RSJ international conference on intelligent robots and systems, pages 2997–3004. IEEE, 2014.
  * [6] Alejandro Perez, Robert Platt, George Konidaris, Leslie Kaelbling, and Tomas Lozano-Perez.  Lqr-rrt*: Optimal sampling-based motion planning with automatically derived extension heuristics.  In 2012 IEEE International Conference on Robotics and Automation, pages 2537–2542. IEEE, 2012.
  * [7] Dustin J. Webb and Jur van den Berg.  Kinodynamic rrt*: Asymptotically optimal motion planning for robots with linear dynamics.  In 2013 IEEE International Conference on Robotics and Automation, pages 5054–5061, 2013.
  * [8] Yanbo Li, Zakary Littlefield, and Kostas E Bekris.  Asymptotically optimal sampling-based kinodynamic planning.  The International Journal of Robotics Research, 35(5):528–564, 2016.
  * [9] Tobias Kunz, Andrea Thomaz, and Henrik Christensen.  Hierarchical rejection sampling for informed kinodynamic planning in high-dimensional spaces.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 89–96. IEEE, 2016.
  * [10] Daqing Yi, Rohan Thakker, Cole Gulino, Oren Salzman, and Siddhartha Srinivasa.  Generalizing informed sampling for asymptotically-optimal sampling-based kinodynamic planning via markov chain monte carlo.  In 2018 IEEE International Conference on Robotics and Automation (ICRA), pages 7063–7070. IEEE, 2018.
  * [11] Alexander Shkolnik, Matthew Walter, and Russ Tedrake.  Reachability-guided sampling for planning under differential constraints.  In 2009 IEEE International Conference on Robotics and Automation, pages 2859–2865. IEEE, 2009.
  * [12] Albert Wu, Sadra Sadraddini, and Russ Tedrake.  R3t: Rapidly-exploring random reachable set tree for optimal kinodynamic planning of nonlinear hybrid systems.  In 2020 IEEE International Conference on Robotics and Automation (ICRA), pages 4245–4251. IEEE, 2020.
  * [13] Sagar Suhas Joshi, Seth Hutchinson, and Panagiotis Tsiotras.  Tie: Time-informed exploration for robot motion planning.  IEEE Robotics and Automation Letters, 6(2):3585–3591, 2021.
  * [14] Yongxing Tang, Zhanxia Zhu, and Hongwen Zhang.  A reachability-based spatio-temporal sampling strategy for kinodynamic motion planning.  IEEE Robotics and Automation Letters, 8(1):448–455, 2022.
  * [15] Ilya S Molchanov and Ilya S Molchanov.  Theory of random sets, volume 19.  Springer, 2005.
  * [16] Haojie Shi and Max Q-H Meng.  Deep koopman operator with control for nonlinear systems.  IEEE Robotics and Automation Letters, 7(3):7700–7707, 2022.
  * [17] Somil Bansal and Claire J Tomlin.  Deepreach: A deep learning approach to high-dimensional reachability.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 1817–1824. IEEE, 2021.
  * [18] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using real nvp.  arXiv preprint arXiv:1605.08803, 2016.
  * [19] Thomas Lew and Marco Pavone.  Sampling-based reachability analysis: A random set theory approach with adversarial sampling.  In Conference on robot learning, pages 2055–2070. PMLR, 2021.
  * [20] Kshitij Tayal, Arvind Renganathan, Rahul Ghosh, Xiaowei Jia, and Vipin Kumar.  Koopman invertible autoencoder: Leveraging forward and backward dynamics for temporal modeling.  arXiv preprint arXiv:2309.10291, 2023.
  * [21] Milan Korda and Igor Mezić.  Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control.  Automatica, 93:149–160, 2018.
  * [22] Stanley Bak, Sergiy Bogomolov, Parasara Sridhar Duggirala, Adam R Gerlach, and Kostiantyn Potomkin.  Reachability of black-box nonlinear systems after koopman operator linearization.  IFAC-PapersOnLine, 54(5):253–258, 2021.
  * [23] Seksan Kiatsupaibul, Robert L. Smith, and Zelda B. Zabinsky.  An analysis of a variation of hit-and-run for uniform sampling from general regions.  21(3), feb 2011.
  * [24] Carl Folkestad and Joel W Burdick.  Koopman nmpc: Koopman-based learning and nonlinear model predictive control of control-affine systems.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 7350–7356. IEEE, 2021.
  * [25] Omri Azencot, N Benjamin Erichson, Vanessa Lin, and Michael Mahoney.  Forecasting sequential data using consistent koopman autoencoders.  In International Conference on Machine Learning, pages 475–485. PMLR, 2020.
  * [26] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  arXiv preprint arXiv:1803.08375, 2018.
  * [27] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [28] Ian M Mitchell, Alexandre M Bayen, and Claire J Tomlin.  A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games.  IEEE Transactions on automatic control, 50(7):947–957, 2005.
  * [29] Alexander B Kurzhanski and Pravin Varaiya.  Ellipsoidal techniques for reachability analysis.  In International workshop on hybrid systems: Computation and control, pages 202–214. Springer, 2000.


capbtabboxtable[][]

# Rethinking Feature Discrimination and Polymerization for Large-scale
Recognition

Yu Liu${}^{~{}}$, Hongyang Li∗, Xiaogang Wang
The Chinese University of Hong Kong
{yuliu, yangli<EMAIL_ADDRESS>
The first two authors contribute equally.

###### Abstract

Feature matters. How to train a deep network to acquire discriminative
features across categories and polymerized features within classes has always
been at the core of many computer vision tasks, specially for large-scale
recognition systems where test identities are unseen during training and the
number of classes could be at million scale. In this paper, we address this
problem based on the simple intuition that the cosine distance of features in
high-dimensional space should be close enough within one class and far away
across categories. To this end, we proposed the congenerous cosine (COCO)
algorithm to simultaneously optimize the cosine similarity among data. It
inherits the softmax property to make inter-class features discriminative as
well as shares the idea of class centroid in metric learning. Unlike previous
work where the center is a temporal, statistical variable within one mini-
batch during training, the formulated centroid is responsible for clustering
inner-class features to enforce them polymerized around the network truncus.
COCO is bundled with discriminative training and learned end-to-end with
stable convergence. Experiments on five benchmarks have been extensively
conducted to verify the effectiveness of our approach on both small-scale
classification task and large-scale human recognition problem.

##  1 Introduction

Deep neural networks [14, 8, 9] have dramatically advanced the computer vision
community in recent years, with high performance boost in tremendous tasks,
such as image classification [8, 34, 21], object detection [4, 32, 20, 30,
27], object tracking [2], saliency detection [19], etc. The essence behind the
success of deep learning resides in both the superior expression power of non-
linear complexity in high-dimension space [9] and the large-scale datasets [3,
7] where the deep networks could, in full extent, learn complicated patterns
and representative features. The basic question still halts the core of many
vision tasks: how can we obtain discriminative features across categories and
polymerized features within one class?

The most fundamental task in computer vision is image classification [34, 14,
8]: given $N$ categories in the training set, we train a deep model that could
best discriminate features among categories in high-dimensional space. The
test set also consists of the same $N$ classes (identities). However, there
are two challenges when we apply the classification task in real-world
applications. The first is the number of categories could be exponentially
increased (e.g., MegaFace [12] contains one million test identities). This
requires the features should be discriminative to have large margin across
classes. The second is the category mismatch between training and test set.
Often we train on a set of $N$ known classes and evaluate the algorithm on $M$
new, different identities. This hastens the features should be polymerized to
share small distances within classes.

In this paper, we optimize the deeply learned features to be both
discriminative and polymerized, given the large-scale constraint: ultra
category number and open test set. Faced with the aforementioned challenges,
the community extends the classification problem into two sub tasks: data
verification and identification. Since our main focus lies in the field of
human-related tasks, we use the term face with data exchangeably thereafter.
Face verification [38, 31] aims at differentiating whether a pair of instances
belong to the same identity; whilst face identification [25, 35] strives for
predicting the class of an identity given the test gallery [12]. Person
recognition [28, 44] also belongs to the verification case, only with more
non-frontal and unstrained face (head) settings111 The term recognition has a
broad sense: detection, verification, classification, etc. Person recognition
is face verification with more unstrained samples. For not making any
confusion, we follow this terminology as did in the community. In this paper,
we refer to face verification, identification and person recognition, in a
general sense, as human recognition..

For the feature discrimination concern, the softmax loss and its variants [26,
31, 25] are proposed to provide a hyperplane with large margin in high-
dimensional space and thus effectively distinguish features across classes.
However, it ignores the intra-class similarity, making the sample distance
within one class distant also and thus amply filling in the feature space
(Fig. 1(a)); for the feature polymerization concern, one can resort to
distance metric learning [43, 41]. Previous work [36] usually introduced
neural networks, using contrastive loss [35] or triplet loss [39, 24, 33], to
learn the features, followed by a simple metric such as Euclidean or cosine
distance to verify the identity (Fig. 1(b)). Such methods are effective to
decrease the intra-class distance in some extent and yet bears an internal
drawback due to the formulation: the selection of positive and negative
samples are relative within each batch and thus could occur training
instability, especially the number of classes and data scale increase (see
Fig. 3(b)).

An alternative is to jointly combine the discriminative feature learning in a
softmax spirit and the polymerized feature learning in a metric learning
manner. The center loss [42] and its variant [45] are two typical approaches.
They linearly combine two terms as the total loss: the first being softmax and
the second a class center supervision $\|\bm{x}_{i}-\bm{c}_{l_{i}}\|$. The
hybrid solution can effectively enlarge the inter-class distance as well as
the inner-class similarity; nevertheless, the concept of class center is
updated in a statistical manner without consulting to the network parameter
update during training. Moreover, the center loss term has to be bundled with
an ad-hoc softmax (thus requires more model parameters); if not, the center
barely changes according to its formulation and thus the loss could be zero
(c.f., Fig. 1(c)-(d)).

Figure 1: Feature visualization under different loss strategies, trained on
MNIST [15] with 10 classes. (a) Softmax loss; (b) triplet loss; (c) center
loss without softmax; (d) center loss; (e) our proposed COCO. The feature
dimension is three for better illustration. The triplet loss [33] has a
constraint on feature embedding ($\|\bm{f}\|_{2}=1$) and thus the learned
features are scatted on a hypersphere.

To this end, we propose the congenerous cosine (COCO) algorithm to provide a
new perspective on jointly considering the feature discrimination and
polymerization. The intuition is that we directly optimize and compare the
cosine distance (similarity) between features: the intra-class distance should
be as minimal as possible whereas the inter-class variation should be
magnified to the most extent (see Fig. 1(e)). COCO inherits the softmax
property to make features discriminative in the high-dimensional space as well
as keeps the idea of class centroid. Note that unlike the statistical role of
$\bm{c}_{l_{i}}$ in [42], the update of our centroid is performed
simultaneously during the network training; and there is no additional loss
term in the formulation. By virtue of class centroid bundled with
discriminative training, COCO is learned end-to-end with stable convergence.
Experiments on two small-scale image classification datasets first verify the
feasibility of our proposed approach and three large-scale human recognition
benchmarks further demonstrate the effectiveness of applying COCO in an ultra
class number and open test set environment. The source code is publicly
available222 https://github.com/sciencefans/coco_loss .

##  2 Related Work

The related work of different loss methods to optimize features is compared
and discussed throughout the paper. In this section, we mainly focus on the
face and person recognition applications.

Face verification and identification have been extensively investigated [12,
31, 25, 36, 35, 42] recently due to the dramatic demand in real-life
applications. Schroff et al. [33] directly learned a mapping from face images
to a compact Euclidean space where distances directly correspond to a measure
of face similarity. DeepID2 [35] employed a contrastive loss with both the
identification and verification supervision. It increased the dimension of
hidden representations and added supervision to early convolutional layers.
Tadmor et al. [36] introduced a multibatch method to generate invariant face
signatures through training pairs. The training tended to be faster to
convergence by virtue of smaller variance of the estimator. SphereFace [25]
formulated an angular softmax loss that enabled CNN to learn discriminative
features by imposing constraints on a hypersphere manifold.

Person recognition in photo albums [1, 29, 28, 22, 44, 18] aims at recognizing
the identity of people in daily life photos, where such scenarios can be
complex with cluttered background. [1] first address the problem by proposing
a Markov random filed framework to combine all contextual cues to recognize
the identity of persons. Recently, Zhang et al. [44] introduce a large-scale
dataset called PIPA for this task. The test set on PIPA is split into two
subsets, namely test_0 and test_1 with roughly the same number of instances.
In [28], a detailed analysis of different cues is explicitly investigated and
three additional test splits are proposed for evaluation. [22] embed scene and
relation contexts in LSTM and formulate person recognition as a sequence
prediction task. Note that previous work [44, 28, 18, 22] use the training set
only for extracting features and a follow-up classifier (SVM or neural
network) is trained on test_0. The recognition system is evaluated on the
test_1 set. We argue that such a practice is infeasible and ad hoc in
realistic application since the second training on test_0 is auxiliary and
needs re-training if new samples are added.

##  3 Optimizing Congenerous Cosine Distance

###  3.1 COCO Formulation and Optimization

Let $\bm{f}^{(i)}\in\mathbb{R}^{D}$ denote the feature vector of the $i$-th
sample,where $D$ is the feature dimension. We first introduce the cosine
similarity of two features as:

$\mathcal{C}({\bm{f}^{(i)}},\bm{f}^{(j)})={\bm{f}^{(i)T}\cdot\bm{f}^{(j)}}\big{/}{\|\bm{f}^{(i)}\|\|\bm{f}^{(j)}\|}.$
(1)

The cosine similarity metric quantifies how close two samples are in the
feature space. A natural intuition to a desirable loss is to increase the
similarity of samples within a category and enlarge the centroid distance of
samples across classes. Let $l_{i},l_{j}\in\\{1,\cdots,K\\}$ be the labels of
sample $i,j$, where $K$ is the total number of categories, we have the
following loss to maximize:

$\displaystyle\mathcal{L}^{naive}=\sum_{i,j\in\mathcal{B}}\frac{\delta(l_{i},l_{j})\mathcal{C}(\bm{{f}}^{(i)},\bm{{f}}^{(j)})}{[1-\delta(l_{i},l_{j})]\mathcal{C}(\bm{{f}}^{(i)},\bm{{f}}^{(j)})+\epsilon},$
(2)

where $\mathcal{B}$ denotes the mini-batch, $\delta(\cdot,\cdot)$ is an
indicator function and $\epsilon$ is a trivial number for computation
stability. The naive design in (2) is reasonable in theory and yet suffers
from computational inefficiency. Since the complexity of the loss above is
$\mathcal{O}(C_{M}^{2})=\mathcal{O}(M^{2})$, it increases quadratically as
batch size $M$ goes bigger. Also the network suffers from unstable parameter
update and is hard to converge if we directly compute loss from two arbitrary
samples from a mini-batch: a similar drawback as does in the triplet loss
[39].

As discussed previously, an effective solution to polymerize the inner-class
feature distance is to provide a series of truncus in the network: serving as
the centroid for each class and thus enforcing features to be learned around
these hubs. To this end, we define the centroid of class $k$ as the average of
features over a mini-batch:
$\bm{c}_{k}=\frac{1}{N_{k}}\sum_{i\in\mathcal{B}}\delta(l_{i},k)\bm{{f}}^{(i)}\in\mathbb{R}^{D}$,
where $N_{k}$ is the number of samples that belong to class $k$ within the
batch. Incorporating the spirit of class centroid into Eqn. (2), one can
derive the following revised loss to maximize:

$\mathcal{L}^{revise}=\sum_{i\in\mathcal{B}}\frac{\exp\mathcal{C}(\bm{{f}}^{(i)},\bm{c}_{l_{i}})}{\sum_{m\neq
l_{i}}\exp\mathcal{C}(\bm{{f}}^{(i)},\bm{c}_{m})},$ (3)

where $m$ indexes along the category dimension. The direct intuition behind
Eqn. (3) is to measure the distance of one sample against other samples by way
of a class centroid, instead of a direct pairwise comparison as in Eqn. (2).
The numerator ensures sample $i$ is close enough to its class center $l_{i}$
and the denominator enforces a minimal distance against samples in other
classes. The exponential operator $\exp$ is to transfer the cosine similarity
to a normalized probability output.

To optimize the aforementioned loss in practice, we first normalize the
feature and centroid by $l_{2}$ norm and then scale the feature before feeding
them into the loss layer:

$\hat{\bm{c}}_{k}=\frac{\bm{c}_{k}}{\|\bm{c}_{k}\|},~{}\hat{\bm{f}}^{(i)}=\frac{\alpha\bm{f}^{(i)}}{\|\bm{f}^{(i)}\|},~{}p_{k}^{(i)}=\frac{\exp(\hat{\bm{c}}_{k}^{T}\cdot\hat{\bm{f}}^{(i)})}{\sum_{m}\exp(\hat{\bm{c}}_{m}^{T}\cdot\hat{\bm{f}}^{(i)})},\\\
$ (4)

where $\alpha$ is the scale factor and $p_{k}^{(i)}$ is an alternative
expression333 In practice, it is more feasibly convenient to include $m=l_{i}$
in the denominator and the result barely differs. to the inner term of
summation in Eqn. (3). Therefore, the proposed congenerous cosine (COCO) loss
are formulated, in a cross-entropy manner to minimize, as follows:

$\displaystyle\mathcal{L}^{COCO}\big{(}\bm{{f}}^{(i)},\bm{c}_{k}\big{)}=-\sum_{i\in\mathcal{B},k}t_{k}^{(i)}\log
p_{k}^{(i)}=-\sum_{i\in\mathcal{B}}\log p_{l_{i}}^{(i)},$ (5)

where $k$ indexes along the class dimension in $\mathbb{R}^{K}$ and
$t_{k}^{(i)}\in\\{0,1\\}$ is the binary mapping of sample $i$ based on its
label $l_{i}$. The proposed COCO inherits the properties from (2) to (3), that
is to increase the discrimination across categories as well as to polymerize
the compactness within one class in a cooperative way. Moreover, it reduces
the computation complexity compared with the naive version and can be
implemented in a cheat way.

Note that both the features and cluster centroids are trained end-to-end. We
now derive the gradients of loss $\mathcal{L}^{COCO}$ w.r.t. the unnormalized
input feature $\bm{f}^{(i)}$ and the centroid $\bm{c}_{k}$. For brevity, we
drop the sample index $i$ and superscript notation in the loss. Denote
$\nabla_{x}\mathcal{L},x\in\\{d,j\\}$ as the element-wise top gradient of loss
w.r.t. the normalized feature, i.e.,
$\partial\mathcal{L}/\partial\hat{f}_{x}$, we have the following gradient by
applying the chain rule:

$\displaystyle\frac{\partial\mathcal{L}}{\partial f_{j}}$
$\displaystyle=\sum_{d}\nabla_{d}\mathcal{L}\cdot\frac{\partial\hat{f}_{d}}{\partial{f}_{j}}=\nabla_{j}\mathcal{L}\cdot\frac{\partial\hat{f}_{j}}{\partial{f}_{j}}+\sum_{d\neq
j}\nabla_{d}\mathcal{L}\cdot\frac{\partial\hat{f}_{d}}{\partial{f}_{j}},$
$\displaystyle=\alpha\bigg{(}\frac{\nabla_{j}\mathcal{L}}{\|\bm{f}\|}-\sum_{d}\nabla_{d}\mathcal{L}\cdot\frac{f_{i}f_{d}}{\|\bm{f}\|^{3}}\bigg{)}=\frac{\alpha\big{(}\nabla_{j}\mathcal{L}-(\nabla\mathcal{L}^{T}\cdot\hat{\bm{f}})\hat{f}_{j}\big{)}}{\|\bm{f}\|},$
(6)

where $\nabla\mathcal{L}\in\mathbb{R}^{D}$ is the vector form of
$\nabla_{x}\mathcal{L}$. Considering the specific loss defined in Eqn. (5),
the top gradient $\nabla_{x}\mathcal{L}$ can be obtained:

$\nabla_{x}\mathcal{L}=\sum_{k}\frac{\partial\mathcal{L}}{\partial
z_{k}}\cdot\frac{\partial
z_{k}}{\partial\hat{f}_{x}}=\sum_{k}(p_{k}-t_{k})\cdot\hat{c}_{kx}.$ (7)

The derivation of gradient w.r.t. centroid ${\partial\mathcal{L}}/{\partial
c_{kj}}$ can be derived in a similar manner. The features are initialized from
pretrain models and the initial value of $\bm{c}_{k}$ is thereby obtained.

###  3.2 Towards an Optimal Scale Factor $\alpha$

As stated in previous subsection, we first normalize and scale the feature;
now we have the following theorem to prove that there exists an optimal value
for the factor. The derivations are provided in the supplementary.

######  Theorem 1.

Given the optimization loss $\mathcal{L}$ has an upper bound $\epsilon$, i.e,
$\mathcal{L}<\epsilon$; and the neural network has a class number of $K$, the
scale factor $\alpha$ enforced on the input feature has a lower boundary:

$\alpha>\frac{1}{2}\log\frac{K-1}{\exp\epsilon-1}.$ (8)

Figure 2: Histogram statistics of cosine distance for positive and negative
pairs on MNIST [15]. Feature discrimination and polymerization are more
evident using COCO than other losses.

###  3.3 COCO Verification on Image Classification

In this subsection, we provide results and discussions to verify the
effectiveness of COCO on small-scale image classification datasets, namely,
MNIST [15] and CIFAR-10 [13]. The dataset, network structure and
implementation details are provided in the supplementary.

An optimal scale factor matters. Take CIFAR-10 as example. When $\alpha$ is
small ($\alpha=1$, error rate 12.4), the lower bound of loss is high, leading
to a high loss in easy samples; however, the loss for hard examples still
remains at low value. This will make the gradients from easy and hard samples
stay at the same level. This is hard for the network to optimize towards hard
samples - leading to a higher error rate. When $\alpha=10$, the error rate is
7.22; if we set the factor to the optimal value defined in Eqn. (8), the error
rate is 6.25.

Feature discrimination and polymerization. Figure 2 shows the histogram of the
cosine distance among positive pairs (i.e., two arbitrary samples that belong
to the same class) and negative pairs. We can see that the distance in softmax
and triplet cases resemble each other; in the COCO case, the discrepancy
between the positive and negative is larger (center distance on x-axis between
two clusters); the intra-class similarity is also more polymerized (area under
two clusters).

Training loss. Figure 3 describes the training loss curves under different
loss strategies on CIFAR and LFW datasets. We can observe that our proposed
method undergoes a stable convergence as the epoch goes; whereas other losses
do not witness an obvious drop when the learning rate decreases. The triplet
loss is effective the number of classes is small on CIFAR; however, when the
task extends to a larger scale and includes over 1,000 categories, the triplet
term suffers from severe disturbance during training since the optimized
feature distance is altered frequently within one batch.

Quantitative results. We report the image classification error rate on MNIST
and CIFAR-10 in Table 1. Note that MNIST is quite a smalls-scale benchmark and
the performance tends to saturate. Moreover, we conduct the ablative study on
different losses: softmax, center loss, triplet loss and their combination.
For fair comparison, they all share the same CNN structure as COCO’s. We can
observe that the center loss alone bears a limited improvement (6.66 vs 6.70)
over softmax; the triplet loss has the training instability concern and thus
the error rate is much higher (12.69). Without resorting to softmax, our
method is neat in formulation and achieves the lowest error (6.25).

Figure 3:  Training loss vs epoch of different losses. (a): Image
classification on CIFAR-10 [13]; (b): Face verification on trained on MS-1M
celebrity set. Note that for the center loss there are two curves,
corresponding to the center and softmax term, respectively.

##  4 COCO for Large-scale Recognition

###  4.1 Face Verification and Identification

We follow the general convention as in [42, 33, 25] to conduct face
verification and identification: train a network using COCO loss to obtain
robust features, extract features in the test set based on the trained model,
and compare or identify the unknown classes given a large scale test setting.
The main modifications to the traditional practice are two folds: first we
employ the detection pipeline to leverage the feature selection and conduct
landmark regression, in a multi-task RPN manner [32]. One task is for
classifying the foreground or background landmarks while the other focus on
regressing them. There are five landmarks: left or right eyes, nose, and left
or right corners of mouth. The second revision is the landmark alignment via
affine transformation, in order to switch the position to a base one.

Table 1: Classification error rate (%) on MNIST and CIFAR-10 without data augmentation.  Method | MNIST | CIFAR-10
---|---|---
DropConnect [37] | 0.57 | 9.41
NIN [23] | 0.47 | 10.47
Maxout [6] | 0.45 | 11.68
DSN [17] | 0.39 | 9.69
R-CNN [5] | 0.31 | 8.69
GenPool [16] | 0.31 | 7.62
Softmax | 0.36 | 6.70
Center loss + softmax | 0.32 | 6.66
Triplet loss | 1.45 | 12.69
Triplet loss + softmax | 0.38 | 6.73
COCO | 0.30 | 6.25

###  4.2 Person Recognition

Following a general pipeline of the verification task, we first train features
using the proposed COCO algorithm on a specific body region (e.g., head) to
classify different identities; during inference, the features on both test
subsets are extracted to compute the cosine distance and to determine whether
both instances are of the same identity. We utilize the multi-region spirit
[22] during training and test to merge results from four regions, namely,
face, head, upper body and whole body. The contributions as regard to this
task are two folds: remove the second training on test_0 as did in previous
work, which credits from the property of COCO to learn discriminative and
polymerized features; align each region patch to a base location to reduce
variation among samples.

The detailed description of region detection and patch alignment are provided
in the supplementary. Given the cropped and aligned patches, we train
(finetune) models for different regions on PIPA training set using COCO loss.
At testing stage, we measure the cosine similarity between two test splits to
recognize the identity in test_1 based on the labels in test_0. The similarity
between two patches $i$ and $j$ in test_1 and test_0 is denoted as
$s(r)_{ij}$, where $r$ indicates a region. We normalize the preliminary
$s(r)_{ij}$ in order to have scores across different regions comparable:

$\displaystyle
s(r)_{ij}=\mathcal{C}(\bm{f}^{(i,r)},\bm{f}^{(j,r)}),~{}~{}~{}\hat{s}(r)_{ij}=\bigg{(}1+\exp\big{[}-\big{(}\beta_{0}+\beta_{1}s(r)_{ij}\big{)}\big{]}\bigg{)}^{-1},$
(9)

where $\beta_{0},\beta_{1}$ are parameters of the logistic regression. The
final score $S_{ij}$ is a weighted mean of the normalized scores
$\hat{s}(r)_{ij}$ of each region:
$S_{ij}=\sum_{r=1}^{R}\gamma^{r}\cdot\hat{s}(r)_{ij}$, where $R$ is the total
number of regions and $\gamma^{r}$ being the weight of each region’s score.
The identity of patch $i$ in test_1 is decided by the label corresponding to
the maximum score in the reference set: $l_{i}=\arg\max_{j*}S_{ij}$. The test
parameters of $\beta$ and $\gamma$ are determined by a validation set on PIPA.

##  5 Experiments

Dataset and evaluation metric.  Face recognition. The Labeled Face in the Wild
(LFW) dataset [10] contains 13,233 web-collected images from 5749 different
identities, with large variations in pose, expression and illuminations.
Following the standard protocol of unrestricted with labeled outside data, we
test on 6,000 face pairs. The MegaFace [12] is a very challenging dataset and
aims to evaluate the performance of face recognition algorithms at the million
scale of distractors: people who are not in the test set. It includes gallery
set and probe set. The gallery set consists of more than one million images
from 690K different individuals, as a subset of Flickr photos from Yahoo. The
probe set descends from two existing databases: Facescrub and FGNet. Person
recognition. The People In Photo Albums (PIPA) dataset [44] is divided into
train, validation, test and leftover sets, where the head of each instance is
annotated in all sets and it consists of over 60,000 instances of around 2,000
individuals. In this work, thanks to the discriminative and polymerized
features learned by COCO, we take full advantage of the training set and
remove the second training on test_0. Moreover, [28] introduced three more
challenging splits besides the original test split, namely album, time and
day. Each new split emphasizes different temporal distance (various albums,
events, days, etc.) between the two subsets of the test data. The evaluation
metric is the averaged classification accuracy over all instances on test_1.

###  5.1 Face Verification and Identification

Figure 4: The ROC curves of different methods for face verification on LFW.
Our model using COCO to train features achieves a higher AUC score.

Method | Train Data | mAcc
---|---|---
FaceNet [33] | 200M | 99.65
DeepID2 [35] * | 300K | 99.47
CenterFace [42] | 700K | 99.28
L-Softmax [26] | webface | 98.71
SphereFace [25] | webface | 99.42
Softmax | half MS-1M | 99.75
Center loss + softmax | half MS-1M | 99.78
Triplet loss | half MS-1M | 98.85
Triplet loss + softmax | half MS-1M | 99.68
COCO | half MS-1M | 99.86

Figure 5: Face verification accuracy (%) on LFW. The first three models use
private outside training data. * Result comes from an ensemble method that has
25 models.

We now apply the COCO feature learning to large-scale human recognition
system. For fair comparison of our method to other losses (softmax, center
loss, triplet loss and their combination), we use the same CNN network
structure in the following experiments. Features on LFW is available on the
github repository.

Figure 5 shows the ROC curves of different methods and Table 5 reports the
accuracy (%) for face verification on LFW [10]. It is observed that COCO model
is competitively superior (99.86) than its counterparts. Compared with FaceNet
[33] which intakes around 200M outside training data, or DeepID2 [35] which
ensembles 25 models, our approach demonstrates the effectiveness of feature
learning. Table 2 illustrates the identification accuracy on MegaFace
challenge [12], where the test identities share no interaction with those in
training and the data scale is quite large. Our model using COCO is competent
at different number of distractors. The combination of center loss with
softmax ranks second (76.57 vs 75.79 @1M), which verifies that applying center
loss alone cannot guarantee a good feature learning for large-scale task. The
triplet loss (69.13) is inferior and our model competes favorably against
those using large protocol of training data.

Method | Protocol | @1M | Method | Protocol | @1M
---|---|---|---|---|---
NTechLab | small | 58.22 | NTechLab | large | 73.30
DeepSense | small | 70.98 | DeepSense | large | 74.80
SphereFace [25] | small | 75.77 | ShanghaiTech | large | 74.05
L-Softmax [26] | small | 67.13 | Google FaceNet | large | 70.50
Method | Protocol | @1M | @100k | @10k | @1k
Softmax | small | 71.17 | 85.78 | 93.22 | 96.14
Center loss + softmax | small | 75.79 | 90.74 | 96.45 | 97.91
Triplet loss | small | 69.13 | 85.64 | 93.38 | 96.66
Triplet loss + softmax | small | 70.22 | 86.77 | 94.71 | 97.02
COCO | small | 76.57 | 91.77 | 96.72 | 98.03
Table 2: Face identification accuracy (%) on MegaFace challenge at various
number of distractors.

###  5.2 Person Recognition

Comparison with state-of-the-arts and visual results. The comparison figure
and visual results are provided in the supplementary due to page limit. We can
see that our recognition system outperforms against previous state-of-the-arts
in all four test splits. Note that the test identities share no interaction
with those in the training set and thus this task is also a challenging open
set verification problem. The visual results show some examples of the
predicted instances by our model, where complex scenes with non-frontal faces
and body occlusion can be handled properly in most scenarios. The last two
columns show failure cases where our method could not differentiate the
identities, this is probably due to the very similar appearance configuration
in these scenarios.

Table 3: Person recognition on PIPA [44]. Investigation on the merging similarity score from different body regions during inference. The top two results in each split set are marked in bold and italic, respectively. COCO loss is applied in all cases if not specified. Method | Face | Head | Upper body | Whole body | original | album | time | day
---|---|---|---|---|---|---|---|---
- | ✓ | ✓ |  |  | 84.17 | 80.78 | 74.00 | 53.75
- | ✓ |  | ✓ |  | 89.24 | 81.46 | 76.84 | 61.48
- | ✓ |  |  | ✓ | 88.40 | 82.15 | 70.90 | 57.87
- |  | ✓ | ✓ |  | 88.76 | 79.15 | 68.64 | 42.91
RNN [22] |  | ✓ | ✓ |  | 84.93 | 78.25 | 66.43 | 43.73
- |  | ✓ |  | ✓ | 87.43 | 77.54 | 67.40 | 42.30
- |  |  | ✓ | ✓ | 81.93 | 73.84 | 62.46 | 34.77
- | ✓ | ✓ | ✓ |  | 87.86 | 80.85 | 71.65 | 59.03
- | ✓ | ✓ |  | ✓ | 88.13 | 82.87 | 73.01 | 55.52
- |  | ✓ | ✓ | ✓ | 89.71 | 78.29 | 66.60 | 52.21
- | ✓ |  | ✓ | ✓ | 91.43 | 80.67 | 70.46 | 55.56
Softmax | ✓ | ✓ | ✓ | ✓ | 88.73 | 80.26 | 71.56 | 50.36
- | ✓ | ✓ | ✓ | ✓ | 92.78 | 83.53 | 77.68 | 61.73

Ablation study on different regions. Table 3 depicts the investigation on
merging the similarity score from different body regions during inference.
[22] also employs a multi-region processing step and we include it in Table 3.
For some region combination, on certain splits our model is superior (e.g.,
head plus upper body region, 88 vs 84 on original, 79 vs 78 on album, 68 vs 66
on time); whereas on other splits ours is inferior (e.g., 42 vs 43 on day).
This is probably due to distribution imbalance among splits: the upper body
region differs greatly in appearance with some instances absent of this
region, making COCO hard to learn the corresponding features. However, under
the score integration scheme, the final prediction (last line) can complement
features learned among different regions and achieves better performance
against [22].

##  6 Conclusion

In this work, we re-investigate the feature learning problem based on the
motivation that deeply learned features should be both discriminative and
polymerized. We address the problem by proposing a congenerous cosine (COCO)
loss, which optimizes the cosine distance among data features to
simultaneously enlarge inter-class variation and intra-class similarity. COCO
can be learned in a neat way with stable end-to-end training. We have
extensively conducted experiments on five benchmarks, from image
classification to face identification, to demonstrate the effectiveness of our
approach, especially applying it to the large-scale human recognition tasks.

## References

  * [1] D. Anguelov, K. chih Lee, S. B. Gokturk, and B. Sumengen.  Contextual identity recognition in personal photo albums.  In CVPR, 2007.
  * [2] Z. Chi, H. Li, H. Lu, and M.-H. Yang.  Dual deep network for visual tracking.  IEEE Trans. on Image Processing, 2017.
  * [3] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In CVPR, 2009.
  * [4] R. Girshick.  Fast R-CNN.  In ICCV, 2015.
  * [5] R. Girshick, J. Donahue, T. Darrell, and J. Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  In CVPR, 2014.
  * [6] I. J. Goodfellow, D. Warde-farley, M. Mirza, A. Courville, and Y. Bengio.  Maxout networks.  In ICML. 2013.
  * [7] Y. Guo, L. Zhang, Y. Hu, X. He, and J. Gao.  MS-Celeb-1M: a dataset and benchmark for large-scale face recognition.  arXiv preprint:1607.08221, 2016.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In CVPR, 2016.
  * [9] G. E. Hinton and R. R. Salakhutdinov.  Reducing the dimensionality of data with neural networks.  Science, 313, 2006.
  * [10] G. B. Huang, M. Ramesh, T. Berg, and E. Learned-Miller.  Labeled faces in the wild: A database for studying face recognition in unconstrained environments.  Technical Report 07-49, University of Massachusetts, Amherst, October 2007\.
  * [11] S. Ioffe and C. Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In ICML. 2015.
  * [12] I. Kemelmacher-Shlizerman, S. M. Seitz, D. Miller, and E. Brossard.  The megaface benchmark: 1 million faces for recognition at scale.  In CVPR, 2016.
  * [13] A. Krizhevsky and G. Hinton.  Learning multiple layers of features from tiny images.  In Technical Report, 2009.
  * [14] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In NIPS, pages 1106–1114, 2012.
  * [15] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  In Proceedings of the IEEE, volume 86, pages 2278–2324, 1998.
  * [16] C.-Y. Lee, P. W. Gallagher, and Z. Tu.  Generalizing pooling functions in convolutional neural networks: Mixed, gated, and tree.  arXiv preprint:1509.08985, 2015.
  * [17] C.-Y. Lee, S. Xie, P. Gallagher, Z. Zhang, and Z. Tu.  Deeply-supervised nets.  arXiv preprint: 1409.5185, 2014.
  * [18] H. Li, J. Brandt, Z. Lin, X. Shen, and G. Hua.  A multi-level contextual model for person recognition in photo albums.  In CVPR, 2016.
  * [19] H. Li, J. Chen, H. Lu, and Z. Chi.  CNN for saliency detection with low-level feature integration.  Neurocomputing, 226:212–220, 2017.
  * [20] H. Li, Y. Liu, W. Ouyang, and X. Wang.  Zoom out-and-in network with map attention decision for region proposal and object detection.  In arXiv preprint: 1709.04347, 2017.
  * [21] H. Li, W. Ouyang, and X. Wang.  Multi-bias non-linear activation in deep neural networks.  In ICML, 2016.
  * [22] Y. Li, G. Lin, B. Zhuang, L. Liu, C. Shen, and A. van den Hengel.  Sequential person recognition in photo albums with a recurrent network.  arXiv preprint:1611.09967, 2016.
  * [23] M. Lin, Q. Chen, and S. Yan.  Network in network.  In ICLR, 2014.
  * [24] H. Liu, Y. Tian, Y. Wang, L. Pang, and T. Huang.  Deep relative distance learning: Tell the difference between similar vehicles.  In CVPR, 2016.
  * [25] W. Liu, Y. Wen, Z. Yu, M. Li, B. Raj, and L. Song.  Sphereface: Deep hypersphere embedding for face recognition.  In CVPR, 2017.
  * [26] W. Liu, Y. Wen, Z. Yu, and M. Yang.  Large-margin softmax loss for convolutional neural networks.  In ICML, 2016.
  * [27] Y. Liu, H. Li, J. Yan, F. Wei, X. Wang, and X. Tang.  Recurrent scale approximation for object detection in CNN.  In IEEE International Conference on Computer Vision, 2017.
  * [28] S. J. Oh, R. Benenson, M. Fritz, and B. Schiele.  Person recognition in personal photo collections.  In ICCV, 2015.
  * [29] S. J. Oh, R. Benenson, M. Fritz, and B. Schiele.  Faceless person recognition; privacy implications in social media.  In ECCV, 2016.
  * [30] W. Ouyang, H. Li, X. Zeng, and X. Wang.  Learning deep representation with large-scale attributes.  In ICCV, 2015.
  * [31] R. Ranjan, C. D. Castillo, and R. Chellappa.  L2-constrained softmax loss for discriminative face verification.  arXiv preprint:1703.09507, 2017.
  * [32] S. Ren, K. He, R. Girshick, and J. Sun.  Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks.  In NIPS, 2015.
  * [33] F. Schroff, D. Kalenichenko, and J. Philbin.  Facenet: A unified embedding for face recognition and clustering.  In CVPR, 2015.
  * [34] K. Simonyan and A. Zisserman.  Very deep convolutional networks for large-scale image recognition.  In ICLR, 2015.
  * [35] Y. Sun, X. Wang, and X. Tang.  Deeply learned face representations are sparse, selective, and robust.  In CVPR, 2015.
  * [36] O. Tadmor, Y. Wexler, T. Rosenwein, S. Shalev-Shwartz, and A. Shashua.  Learning a metric embedding for face recognition using the multibatch method.  In NIPS, 2016.
  * [37] L. Wan, M. Zeiler, S. Zhang, Y. L. Cun, and R. Fergus.  Regularization of neural networks using dropconnect.  In ICML, 2013.
  * [38] F. Wang, X. Xiang, J. Cheng, and A. L. Yuille.  Normface: L2 hypersphere embedding for face verification.  arXiv preprint:1704.06369, 2017.
  * [39] X. Wang and A. Gupta.  Unsupervised learning of visual representations using videos.  In ICCV, 2015.
  * [40] S.-E. Wei, V. Ramakrishna, T. Kanade, and Y. Sheikh.  Convolutional pose machines.  In CVPR, 2016.
  * [41] K. Q. Weinberger and L. K. Saul.  Distance metric learning for large margin nearest neighbor classification.  Journal of Machine Learning Research, 2009.
  * [42] Y. Wen, K. Zhang, Z. Li, and Y. Qiao.  A discriminative feature learning approach for deep face recognition.  In ECCV, 2016.
  * [43] E. P. Xing, A. Y. Ng, M. I. Jordan, and S. Russell.  Distance metric learning, with application to clustering with side-information.  In NIPS, 2003.
  * [44] N. Zhang, M. Paluri, Y. Taigman, R. Fergus, and L. Bourdev.  Beyond frontal faces: Improving person recognition using multiple cues.  In CVPR, 2015.
  * [45] X. Zhang, Z. Fang, Y. Wen, Z. Li, and Y. Qiao.  Range loss for deep face recognition with long-tail.  arXiv preprint:1611.08976, 2016.

## Supplementary

###  6.1 Proof on the optimal scale factor, Theorem 8

###### Proof.

It is obvious that the infimum value of COCO loss is of the form:

$\displaystyle\inf\mathcal{L}^{COCO}=-\log\frac{\exp s^{+}}{\exp
s^{+}+(K-1)\exp s^{-}},$ (10) $\displaystyle
s^{+}=\alpha\exp(\cos\theta)\leq\alpha,$ (11) $\displaystyle
s^{-}=\alpha\exp(\cos\theta)\geq-\alpha,$ (12)

where $s$ is the probability output in the final classification layer. For
brevity, we denote $\inf s=s^{-}$ and $\sup s=s^{+}$; $\theta$ is the angle
between two features. Often we have $\mathbb{E}[\cos\theta]\geq 0$ in practice
and yet the minimal value of $\cos\theta$ can still be $-1$. Therefore,

$\displaystyle\sup\mathcal{L}^{COCO}=$
$\displaystyle-\log\frac{\exp(-\alpha)}{\exp(-\alpha)+(K-1)\exp\alpha}.$ (13)

The loss to minimize in an optimization problem usually meets some upper bound
criteria, and we denote this target value as $\epsilon$, i.e.,
$\inf\mathcal{L}^{COCO}<\epsilon$. Incorporating this prior knowledge into
previous derivations, we have:

$\displaystyle\log(\exp\alpha+(K-1)\exp(-\alpha))<\alpha+\epsilon,$ (14)
$\displaystyle\alpha>\frac{1}{2}\log\frac{K-1}{\exp\epsilon-1}.$ (15)

∎

It is empirically found that a typical value for $\epsilon$ is $10^{-4}$ and
thus the scale factor can be expressed in a deterministic closed form:

$\alpha=\frac{1}{2}\log(K-1)+3.$ (16)

###  6.2 Details on Image Classification and Face Recognition

For image classification on CIFAR, we use a simple ResNet-20 structure and
train the network from scratch. We set the momentum as 0.9 and the weight
decay to be 0.005. The base learning rate is set to be 0.1, 0.1, 0.05,
respectively. We drop the learning rate by 10way and stop to decrease the
learning rate until it reaches a minimum value (0.0001). All the convolutional
layers are initialized with Gaussian distribution with mean of zero and
standard variation of 0.05 or 0.1.

For face recognition, we use the Inception ResNet model and add a new fully-
connected layer after the final global pooling layer to generate a 128-dim
feature. The model is trained from scratch. The training data is a subset of
the Microsoft-1M celebrity dataset, which consists of 80k identities with 3M
images in total, removing all overlapping IDs. For COCO, softmax and center
loss, we train the network with SGD solver with a base learning rate 0.01 and
drop it to 10% every five epoch. The total training time is 12 epoch. The
momentum is set to 0.9. Batch normalization [11] with moving average is
employed. For the triplet loss, we have a base learning rate to 0.1 with a
bigger batch size 128*3 (it is hard to converge if the size is less than
64*3); we drop the rate to 10% every 90k iteration. The total training time is
around 300k iterations.

Evaluation metric on face recognition. For verification, the test set is LFW
[10] and the metric is under standard ‘Unrestricted, Labeled Outside Data’
setting of LFW. It verifies 6,000 pairs to determine if they belong to the
same person. We report ROC curve and the mean accuracy. For identification on
MegaFace [12], we have the test set into two subsets: the first is probe set
that consists of FaceScrub (3530 images of 80 persons); the second is gallery
set that consists of distractors (1M images with no overlap with FaceScrub).
The test metric is under standard MegaFace Challenge-1 setting. We report the
CMC curve in terms of Top-1 accuracy with 10, 100, 1000, 10k, 100k and 1M
distractors.

###  6.3 More Algorithm Details for Person Recognition

The COCO features are trained on four regions (thus four models), namely,
face, head, upper body and whole body. We take the face region as an
illustration for the alignment scheme. A face detector is first pretrained in
a RPN spirit [32]. The face detector identifies some keypoints (eye, brow,
mouth, etc.) and we align the detected face patch to a base location via an
affine transformation. Let $p,q\in\mathbb{R}^{z\times 2}$ denote $m$ keypoints
and the aligned locations, respectively. We define $P,Q$ as two affine spaces,
then the transformation $\mathcal{A}:P\mapsto Q$ is defined as: $p\mapsto
q=\bm{A}p+b$, where $\bm{A}\in\mathbb{R}^{z\times z}$ is a linear
transformation matrix in $P$ and $b\in\mathbb{R}^{z\times 2}$ being the bias
in $Q$.  Such an alignment scheme ensures samples both within and across
categories do not have large variance: if the model is learned without
alignment, it has to distinguish more patterns, e.g., different rotations
among persons, making it more prone to overfitting; with alignment, it can
better classify features of various identities despite of rotation, viewpoint,
translation, etc.  Given the cropped and aligned patches, we finetune the face
model on PIPA training set using COCO loss.

The head region is given as the ground truth for each person. To detect a
whole body, we also pre-train a detector in the RPN framework. The model is
trained on the large-scale human celebrity dataset [7], where we use the first
87021 identities in 4638717 images. The network structure is an inception
model [11] with the final pooling layer replaced by a fully connected layer.
To determine the upper body region, we conduct human pose estimation [40] to
identity keypoints of the body and the upper part is thereby located by these
points. The head, whole body and upper body models, which are used for COCO
loss training, are finetuned on PIPA training set using the pretained
inception model, following similar procedure of patch alignment stated
previously for the face region. The aligned patches of four regions are shown
in Figure 6(c).

Figure 6: Training workflow of person recognition task. (a) Each person is
labelled with a ground-truth box of the head. (b): Face, whole body and upper
body detection. (c): Each region (patch) is aligned to a base position to
alleviate the inner-class variance. (d) Each aligned patch is further fed into
a deep model to obtain representative and robust features using COCO loss.

###  6.4 Experiments on Person Recognition

Table 4: Ablation study on body regions and feature alignment using softmax as loss. We report the classification accuracy (%) where n-a denotes the non-alignment case and ali indicates the alignment case. Note that for face region, we only evaluate instances with faces. Test split | Face | Head | Upper body | Whole body
---|---|---|---|---
n-a | ali | n-a | ali | [22] | n-a | ali | [22] | n-a | ali
original | 95.47 | 97.45 | 74.23 | 82.69 | 81.75 | 76.67 | 80.75 | 79.92 | 75.04 | 79.06
album | 94.66 | 96.57 | 65.47 | 73.77 | 74.21 | 66.23 | 69.58 | 70.78 | 64.21 | 67.27
time | 91.03 | 93.36 | 55.88 | 64.31 | 63.73 | 55.24 | 57.40 | 58.80 | 55.53 | 54.62
day | 90.36 | 91.32 | 35.27 | 44.24 | 42.75 | 26.49 | 32.09 | 34.61 | 32.85 | 29.59

Ablation study on feature alignment for different regions. Table 4 reports the
performance of using feature alignment and different body regions, where
several remarks could be observed. First, the alignment case in each region
performs better by a large margin than the non-alignment case, which verifies
the motivation of patch alignment to alleviate inner-class variance stated in
the main paper. Second, for the alignment case, the most representative
features to identify a person reside in the region of face, followed by head,
upper body and whole body at last. Such a clue is not that obvious for the
non-alignment case. Third, we notice that for the whole body region, accuracy
in the non-alignment case is higher than that of the alignment case in time
and day. This is probably due to the improper definition of base points on
these two sets.

Table 5: Verification accuracy (%) comparison to state-of-the-arts on PIPA [44] for person recognition.  Methods | original | album | time | day
---|---|---|---|---
PIPER [44] | 83.05 | - | - | -
RNN [22] | 84.93 | 78.25 | 66.43 | 43.73
Naeil [28] | 86.78 | 78.72 | 69.29 | 46.61
Multi-context [18] | 88.75 | 83.33 | 77.00 | 59.35
Ours | 92.78 | 83.53 | 77.68 | 61.73

Visual results. Figure 7 shows some visualization results for person
recognition task.

Figure 7: Visualization of our method on the PIPA test_1 set. Given the image
as input to be predicted in (a), its nearest neighbour in the feature space
from text_0 is shown in (b). The identity of the input is determined by the
label of its neighbour. Our model can handle complex scenes with non-frontal
faces and body occlusion. The last two columns show failure cases.