# Verification and Synthesis of Compatible Control Lyapunov and Control
Barrier Functions

Hongkai Dai∗1, Chuanrui Jiang∗2, Hongchao Zhang2, and Andrew Clark2 1Hongkai
Dai is with the Toyota Research Institute, Los Altos, CA, USA. Email:
<EMAIL_ADDRESS>Jiang, Hongchao Zhang and Andrew Clark are
with the Electrical and Systems Engineering Department, McKelvey School of
Engineering, Washington University in St. Louis, St. Louis, MO 63130
<EMAIL_ADDRESS>hongchao<EMAIL_ADDRESS>denotes equal
contributionThis work was supported by the Air Force Office of Scientific
Research grants FA9550-22-1-0054 and FA9550-23-1-0208 and NSF grant
CNS-1941670.

###### Abstract

Safety and stability are essential properties of control systems. Control
Barrier Functions (CBFs) and Control Lyapunov Functions (CLFs) have been
proposed to ensure safety and stability respectively. However, previous
approaches typically verify and synthesize the CBFs and CLFs separately,
satisfying their respective constraints, without proving that the CBFs and
CLFs are compatible with each other, namely at every state, there exists
control actions that satisfy both the CBF and CLF constraints simultaneously.
There exists some recent works that synthesized compatible CLF and CBF, but
relying on nominal polynomial or rational controllers, which is just a
sufficient but not necessary condition for compatibility. In this work, we
investigate verification and synthesis of compatible CBF and CLF independent
from any nominal controllers. We derive exact necessary and sufficient
conditions for compatibility, and further formulate Sum-Of-Squares program for
the compatibility verification. Based on our verification framework, we also
design an alternating nominal-controller-free synthesis method. We evaluate
our method in a linear toy, a non-linear toy, and a power converter example.

##  I Introduction

Safety and stability are essential for control systems that must be ensured to
prevent catastrophic economic damage and loss of human life, while still
achieve the desired goal [1, 2, 3, 4]. Lyapunov stability guarantees that all
trajectories starting within the system’s region-of-attraction (RoA) will
converge to the goal state. Meanwhile, safety properties are normally
formulated as the positive invariance [5] of given regions in the state space.

The importance of safety has motivated numerous methods such as artificial
potential fields [6], reachability analysis[7], Hamilton–Jacobi analysis[8],
model predictive control [9], and control barrier functions [10] to synthesize
safe control strategies. Among these methods, energy-based methodologies such
as barrier certificates [11] and Control Barrier Functions (CBFs) [10] have
been proposed that construct an energy function that is positive if the system
is safe, and then demonstrate safety by proving that the energy function
remains positive for all time. Similarly, Lyapunov stability is certified by a
Lyapunov function [12] or a Control Lyapunov function (CLF) [13], which
remains non-negative and decreasing along any trajectories, and achieves the
minimal value at the goal state. One major benefit of CLFs/CBFs is that they
can be incorporated as constraints within optimization-based controllers,
resulting in CLF-CBF control policies for joint safety and stability [5, 14].
Such techniques have been widely used in safety critical scenarios such as
adaptive cruise control[5, 15], bipedal robot control [10, 16], and multi-
waypoint navigation [3], with many different approaches to synthesize the CLF-
CBF functions [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

A key challenge in CLF-CBF-based safe control policies is the _compatibility_
of the CLF and CBF constraints, namely the existence of at least one
admissible control action satisfying both CLF and CBF constraints
simultaneously at every state in the safe region. Most previous approaches
address incompatibility by relaxing the CLF constraints [5], and thus
prioritizing safety over stability. However, this may sacrifice the stability
and result in the system remaining in undesired states [32]. Recent efforts
[33, 34] have utilized Sum-Of-Squares (SOS) methods [35, 36] to verify and
synthesize compatible CLF and CBF based on nominal controllers. In such
methods, a nominal controller is parameterized as a polynomial (or rational)
function of state; these methods search for the CLF/CBF together with this
nominal controller, such that this nominal controller simultaneously satisfies
both the CLF and CBF constraints. The disadvantages of this approach is that
it is only a sufficient but not necessary condition for compatibility, as
there might exist a non-polynomial (or non-rational) compatible controller,
which cannot be captured by the explicit parameterization. A related but
orthogonal direction of research synthesizes closed-form controllers for CBF
and CLF that are known to be compatible [37].

In this paper, we study the problem of exact verification and synthesis of
compatible CLF/CBF for all states in the safe region. We consider a nonlinear
control-affine system, both with and without input limits, and derive
necessary and sufficient conditions for compatibility using Farkas’ Lemma and
the algebraic-geometric Postivstellensatz. We formulate the verification as a
SOS program and further simplify it using S-procedure. To synthesize
compatible CLF/CBF, we propose a bi-linear alternating method based on the
simplified verification program. To summarize, this paper makes the following
contributions.

  * •

We derived the necessary and sufficient condition for CLF-CBF compatibility,
and formulated a SOS program independent from any nominal controllers for
compatibility verification.

  * •

We next proposed an alternating method based on verification program to
synthesis compatible CBF and CLF, also independent from any nominal
controllers.

  * •

We demonstrated our method in a linear toy, a non-linear toy and a power
converter examples.

The rest of the paper is organized as follows. Section II introduces the
mathematical background on Positivstellensatz and Farkas’ lemma. Section III
presents the definition of the system model, CBF and CLF. The proposed
verification and synthesis framework are introduced in Section IV. Section V
gives the simulation results, and we conclude the paper in Section VI.

Notation: We use lower case letter ($x$) to denote a scalar, bold lower case
($\mathbf{x}$) to denote a vector, upper case letter ($X$) for a matrix, and
calligraphic letter ($\mathcal{X}$) to denote a set.

##  II Mathematical background

A polynomial $p(\mathbf{x})$ is sum-of-squares (sos) iff

$p(\mathbf{x})=\sum_{i=1}^{n}(p_{i}(\mathbf{x}))^{2}$

for some polynomials $p_{i}(\mathbf{x})$. A sos polynomial is non-negative. It
can be verified that a polynomial is sos through convex optimization.

The cone generated from a set of polynomials $\phi_{1},\ldots,\phi_{s}$ is
given by

$\Sigma[\phi_{1},\ldots,\phi_{s}]=\left\\{\sum_{S\subseteq\\{1,\ldots,s\\}}{\alpha_{S}(\mathbf{x})\prod_{i\in
S}{\phi_{i}(\mathbf{x})}}:\right.\\\ \left.\alpha_{S}(\mathbf{x})\in\mbox{ SOS
}\forall S\subseteq\\{1,\ldots,s\\}\right\\}.$

The ideal generated from polynomials $\psi_{1},\ldots,\psi_{r}$ is given by

$\mathbb{I}[\psi_{1},\ldots,\psi_{r}]=\left\\{\sum_{l=1}^{r}{a_{i}(\mathbf{x})\psi_{l}(\mathbf{x})}:a_{1},\ldots,a_{r}\mbox{
polynomials}\right\\}.$

######  Theorem 1 (Positivstellensatz (P-satz)[36])

Let $\left(\phi_{j}\right)_{j=1,\ldots,s}$,
$\left(\psi_{\ell}\right)_{\ell=1,\ldots,r}$ be finite families of polynomials
in $\mathbb{R}\left[\mathbf{x}\right]$. Denote the cone generated from
$\left(\phi_{j}\right)_{j=1,\ldots,s}$ by $\Sigma[\phi_{1},\ldots,\phi_{s}]$ ,
and the ideal generated from $\left(\psi_{\ell}\right)_{\ell=1,\ldots,r}$ by
$\mathbb{I}[\psi_{1},\ldots,\psi_{r}]$. Then, the following statements are
equivalent:

  1. 1.

The set

$\left\\{\begin{array}[]{l|ll}\mathbf{x}\in\mathbb{R}^{n}&\begin{array}[]{ll}\phi_{j}(\mathbf{x})\geq
0,&j=1,\ldots,s\\\
\psi_{\ell}(\mathbf{x})=0,&l=1,\ldots,r\end{array}\end{array}\right\\}$

is empty.

  2. 2.

There exist $\theta(\mathbf{x})\in\Sigma,\lambda(\mathbf{x})\in\mathbb{I}$
such that $\lambda(\mathbf{x})+\theta(\mathbf{x})+1=0$.

In Positivstellensatz, (2) is a necessary and sufficient condition for (1). On
the other hand, (2) requires computing the cone $\Sigma$, which needs
exponential number of polynomials $\alpha_{S}(\mathbf{x})$ (as $S$ permutes
over all subsets of $\\{1,\ldots,s\\}$). To remedy this computation challenge,
people often resort to S-procedure[36], which replaces the cone $\Sigma$ with
the quadratic module $\mathcal{Q}$

$\mathcal{Q}[\phi_{1},\ldots,\phi_{s}]=\\\
\\{\alpha_{0}(\mathbf{x})+\alpha_{1}(\mathbf{x})\phi_{1}(\mathbf{x})+\ldots+\alpha_{s}(\mathbf{x})\phi_{s}(\mathbf{x}),\alpha_{i}(\mathbf{x})\in\text{SOS}\\}.$

######  Lemma 1 (S-procedure)

Eq. (2) in theorem 1 is a sufficient condition for (1) if we replace the cone
$\theta\in\Sigma$ with the quadratic module
$\theta\in\mathcal{Q}[\phi_{1},\ldots,\phi_{s}]$.111Under some stronger
assumptions on $\phi(\mathbf{x}),\psi(\mathbf{x})$, such as Archimedean,
S-procedure also gives a necessary and sufficient condition [38].

Farkas’ Lemma gives conditions for non-existence of solutions to linear
inequalities and is presented as follows.

######  Lemma 2 (Farkas’ Lemma [39])

Let $\Lambda\in\mathbf{R}^{m\times n}$ be a matrix, and
$\mathbf{\xi}\in\mathbb{R}^{m}$ be a vector. Then the system of inequalities
$\Lambda{\mathbf{x}}\leq\mathbf{\xi}$ has a solution if and only if the set
$\\{\mathbf{z}\;|\;\mathbf{z}\geq
0,\mathbf{z}^{T}\Lambda=\mathbf{0}^{T},\mathbf{\xi}^{T}\mathbf{z}=-1\\}$ is
empty.

##  III Problem statement

###  III-A System Model

We consider a continuous-time nonlinear control-affine system in the form

$\displaystyle\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+g(\mathbf{x})\mathbf{u}$
(1a)
$\displaystyle\mathbf{u}\in\mathcal{U}=\\{\mathbf{u}\;|\;A\mathbf{u}\leq\mathbf{c}\\},$
(1b)

where $\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is the state and
$\mathbf{u}\in\mathcal{U}\subset\mathbb{R}^{n_{u}}$ is the control input. We
also assume that $f:\mathcal{X}\rightarrow\mathbb{R}^{n_{x}}$ and
$g:\mathcal{X}\rightarrow\mathbb{R}^{n_{x}\times n_{u}}$ are polynomial
functions of $\mathbf{x}$. The set of admissible control $\mathcal{U}$ is
represented as a polyhedron $A\mathbf{u}\leq\mathbf{c},A\in\mathbb{R}^{p\times
n_{u}},\mathbf{c}\in\mathbb{R}^{p}$. Note that this parameterization of
$\mathcal{U}$ includes the entire set $\mathcal{R}^{n_{u}}$ (unconstrained
control) as a special case, where $A$ and $\mathbf{c}$ are empty, hence we can
model systems with or without input limits. A control policy
$\mu:\mathcal{X}\rightarrow\mathcal{U}$ maps a state $\mathbf{x}$ to a control
input $\mathbf{u}$. Without loss of generality we assume the equilibrium state
is $\mathbf{x}=\mathbf{0}$. The state space $\mathcal{X}$ is divided into a
safe region $\mathcal{X}_{s}$ and an unsafe region
$\mathcal{X}_{u}=\mathcal{X}-\mathcal{X}_{s}$. Our goal is to design a control
policy $\mu$ which renders $\mathcal{X}_{s}$ positively invariant, and the
trajectories converging to $\mathbf{x}=\mathbf{0}$.

###  III-B Control Barrier Function (CBF) and Control Lyapunov Function (CLF)

We consider safety constraints that are specified as super-level sets of
differentiable functions $b:\mathcal{X}\rightarrow\mathbb{R}$. In particular,
we consider sets $\mathcal{C}$ that satisfy

$\displaystyle\mathcal{C}$ $\displaystyle=\\{\mathbf{x}\;|\;b(\mathbf{x})\geq
0\\}.$

If $\mathcal{C}\subseteq\mathcal{X}_{s}$ and $\mathcal{C}$ is positively
invariant, then any trajectory with $\mathbf{x}(0)\in\mathcal{C}$ will be safe
for all time $t$.

######  Definition 1 (CBF [10])

A function $b(\mathbf{x})$ is a CBF if

$\sup_{\mathbf{u}\in\mathcal{U}}\\{L_{f}b(\mathbf{x})+L_{g}b(\mathbf{x})\mathbf{u}+\kappa_{b}(b(\mathbf{x}))\\}\geq
0$ (2)

for all $\mathbf{x}\in\mathcal{C}$, where $L_{f}b(\mathbf{x})=\frac{\partial
b(\mathbf{x})}{\partial\mathbf{x}}\mathbf{f}(\mathbf{x})$,
$L_{g}b(x)=\frac{\partial b(\mathbf{x})}{\partial\mathbf{x}}g(\mathbf{x})$,
and $\kappa_{b}(\cdot)$ is an extended class $\kappa$ function222Recall
function $\kappa(.)$ is extended class $\kappa$ if it is monotone, increasing,
and satisfies $\kappa(0)=0$..

A special class of Definition 1 consists of Exponential Control Barrier
Functions, for which $\kappa_{b}(b(\mathbf{x}))=\kappa_{b}b(\mathbf{x})$ where
$\kappa_{b}$ is a positive constant. Constraint (2) ensures the set
$\mathcal{C}$ positive invariant for system represented by (1) [10].

A control Lyapunov function is defined as follows.

######  Definition 2 (CLF [40])

A differentiable function $V:\mathcal{X}\rightarrow\mathbb{R}$ is a Control
Lyapunov Function if $V(\mathbf{0})=0$, $V(\mathbf{x})>0\text{ if
}\mathbf{x}\neq 0$ and

$\inf_{\mathbf{u}\in\mathcal{U}}\\{L_{f}V(\mathbf{x})+L_{g}V(\mathbf{x})\mathbf{u}+\kappa_{V}(V(\mathbf{x}))\\}\leq
0$ (3)

for all $\\{\mathbf{x}\;|\;V(\mathbf{x})\leq 1\\}$, where $\kappa_{V}(.)$ is
also an extended class $\kappa$ function. This proves that the sub-level set
$\\{\mathbf{x}\;|\;V(\mathbf{x})\leq 1\\}$ is an inner approximation of the
region-of-attraction (ROA).

If a system satisfies constraint (3) for all time $t$, then the system can be
stabilized to the origin [40]. Also, in this paper, we use exponential
stability so that $\kappa_{V}(V(\mathbf{x}))=\kappa_{V}V(\mathbf{x})$, where
$\kappa_{V}$ is a positive constant.

##  IV Methodology

In this section, we first present our algorithm to verify compatibility of CLF
and CBF through convex optimization (section IV-A); then build upon the
verification formulation, we present our algorithm to synthesize compatible
CLF/CBF through solving a sequence of convex optimization problems (section
IV-B).

###  IV-A Verification of compatibility

For a candidate CBF $b(\mathbf{x})$ and a candidate CLF $V(\mathbf{x})$,
define 0-super-level set $\mathcal{C}=\\{\mathbf{x}\;|\;b(\mathbf{x})\geq
0\\}$ and the 1-sub-level set
$\mathcal{D}:=\\{\mathbf{x}\;|\;V(\mathbf{x})\leq 1\\}$ (this 1-sub-level set
will be an inner approximation of the region-of-attraction). We formally state
the compatible CLF-CBF as follows

######  Problem 1 (Compatible CLF-CBF)

Given positive constants $\kappa_{V}$, $\kappa_{b}$, a positive definite
function $V(\mathbf{x})$, and a function $b(\mathbf{x})$ such that
$\mathcal{C}\subseteq\mathcal{X}_{s}$, verify that
$\forall\mathbf{x}\in\mathcal{C}\cap\mathcal{D}$, there exists $\mathbf{u}$
such that

$\begin{bmatrix}L_{f}V(\mathbf{x})\\\ -L_{f}b(\mathbf{x})\\\
\mathbf{0}\end{bmatrix}+\begin{bmatrix}L_{g}V(\mathbf{x})\\\
-L_{g}b(\mathbf{x})\\\
A\end{bmatrix}\mathbf{u}\leq\begin{bmatrix}-\kappa_{V}V(\mathbf{x})\\\
\kappa_{b}b(\mathbf{x})\\\ \mathbf{c}\end{bmatrix}$ (4)

namely there exists $\mathbf{u}\in\mathcal{U}$ satisfying both the CBF
condition (1) and the CLF condition (2) simultaneously. The region
$\mathcal{C}\cap\mathcal{D}$ is certified to be safe and stabilizable.

We call the functions $V(\mathbf{x})$ and $b(\mathbf{x})$ compatible CLF and
CBF if they satisfy the conditions of Problem 1, with
$\mathcal{C}\cap\mathcal{D}$ as the compatible region.

####  IV-A1 Exact Verification

Verifying the condition (4) is equivalent to verifying that there exists a
solution $\mathbf{u}$ for the inequality

$\Lambda(\mathbf{x})\mathbf{u}\leq\xi(\mathbf{x}),\
\forall\mathbf{x}\in\mathcal{C}\cap\mathcal{D}$ (5)

where

$\Lambda(\mathbf{x})=\begin{bmatrix}L_{g}V(\mathbf{x})\\\
-L_{g}b(\mathbf{x})\\\
A\end{bmatrix},\mathbf{\xi}(\mathbf{x})=\begin{bmatrix}-\kappa_{V}V(\mathbf{x})-L_{f}V(\mathbf{x})\\\
\kappa_{b}b(\mathbf{x})+L_{f}b(\mathbf{x})\\\ \mathbf{c}\end{bmatrix}.$

Note that $\Lambda(\mathbf{x})\in\mathbb{R}^{(p+2)\times n_{u}}$ and
$\mathbf{\xi}(\mathbf{x})\in\mathbb{R}^{p+2}$. We let
$\Lambda_{i}(\mathbf{x})$ denote the $i$-th column of $\Lambda(\mathbf{x})$.

Based on Farkas’ Lemma, we give an equivalent condition (Lemma 3) for CLF-CBF
compatibility.

###### Lemma 3

The CBF $b(\mathbf{x})$ and CLF $V(\mathbf{x})$ are compatible if and only if
the following set is empty:

$\\{(\mathbf{x},\mathbf{z})\;|\;b(\mathbf{x})\geq 0,V(\mathbf{x})\leq 1,\\\
\mathbf{z}^{T}\Lambda(\mathbf{x})=\mathbf{0}^{T},\mathbf{\xi}(\mathbf{x})^{T}\mathbf{z}=-1,\mathbf{z}\geq\mathbf{0}\\},$
(6)

where $\mathbf{z}\in\mathbb{R}^{p+2}$ and
$(\mathbf{x},\mathbf{z})\in\mathbb{R}^{n_{x}+p+2}$ is a concatenated vector.

The conclusion of Lemma 3 can be further simplified. Define the element-wise
squared operation of vector $\mathbf{y}=[y_{1},\dots,y_{p+2}]^{T}$ as
$\mathbf{y}^{2}=[y_{1}^{2},\dots,y_{p+2}^{2}]^{T}$. Then an equivalent
formulation of Lemma 3 is as follows.

###### Lemma 4

The CBF $b(\mathbf{x})$ and CLF $V(\mathbf{x})$ are compatible if and only if
the following set is empty:

$\\{(\mathbf{x},\mathbf{y})\;|\;b(\mathbf{x})\geq 0,V(\mathbf{x})\leq 1,\\\
(\mathbf{y}^{2})^{T}\Lambda(\mathbf{x})=\mathbf{0}^{T},\mathbf{\xi}(\mathbf{x})^{T}\mathbf{y}^{2}=-1\\},$
(7)

where $\mathbf{y}^{2}\in\mathbb{R}_{\geq 0}^{p+2}$ and
$(\mathbf{x},\mathbf{y})\in\mathbb{R}^{n_{x}+p+2}$ is a concatenated vector.
Namely, we replace the vector $\mathbf{z}$ and the non-negative constraint
$\mathbf{z}\geq 0$ with a new vector $\mathbf{y}^{2}$ which is non-negative by
construction.

Applying P-satz to Lemma 4 gives the following theorem. Noted that
$(\mathbf{y}^{2})^{T}\Lambda(\mathbf{x})\in\mathbb{R}^{n_{u}}$, each element
of $(\mathbf{y}^{2})^{T}\Lambda(\mathbf{x})$ is computed by
$(\mathbf{y}^{2})^{T}\Lambda_{i}(\mathbf{x}),i=1,2\dots n_{u}$ , and also,
$\mathbf{\xi}(\mathbf{x})^{T}\mathbf{y}^{2}$ is a scalar.

###### Theorem 2

The given CLF $V(\mathbf{x})$ and CBF $b(\mathbf{x})$ are compatible if and
only if there exist polynomials
$\lambda(\mathbf{x},\mathbf{y}),\theta(\mathbf{x},\mathbf{y})$ satisfying:

  1. 1.

$\theta(\mathbf{x},\mathbf{y})\in\Sigma[b(\mathbf{x}),1-V(\mathbf{x})]$

  2. 2.

$\lambda(\mathbf{x},\mathbf{y})\in\mathbb{I}[(\mathbf{y}^{2})^{T}\Lambda_{1}(\mathbf{x}),\dots,(\mathbf{y}^{2})^{T}\Lambda_{n_{u}}(\mathbf{x}),\mathbf{\xi}(\mathbf{x})^{T}\mathbf{y}^{2}+1]$

  3. 3.

$\theta(\mathbf{x},\mathbf{y})+\lambda(\mathbf{x},\mathbf{y})+1=0$

As explained in section II, directly applying P-satz might be computationally
challenging. To remedy this, we could apply the S-procedure (Lemma 1) and
formulate the following convex optimization problem to verify compatibility

Find
$\displaystyle\;\mathbf{s}_{0}(\mathbf{x},\mathbf{y}),s_{1}(\mathbf{x},\mathbf{y}),s_{2}(\mathbf{x},\mathbf{y}),s_{3}(\mathbf{x},\mathbf{y})$
(8a) $\displaystyle\begin{split}\text{s.t
}&-1-\mathbf{s}_{0}(\mathbf{x},\mathbf{y})^{T}((\mathbf{y}^{2})^{T}\Lambda(\mathbf{x}))-s_{1}(\mathbf{x},\mathbf{y})(\mathbf{\xi}(\mathbf{x})^{T}\mathbf{y}^{2}+1)\\\
&\qquad-
s_{2}(\mathbf{x},\mathbf{y})(1-V(\mathbf{x}))-s_{3}(\mathbf{x},\mathbf{y})b(\mathbf{x})\in\text{SOS}\end{split}$
(8b) $\displaystyle
s_{2}(\mathbf{x},\mathbf{y}),s_{3}(\mathbf{x},\mathbf{y})\in\text{SOS}.$ (8c)

A given pair of CLF $V(\mathbf{x})$ and CBF $b(\mathbf{x})$ are compatible if
we can find the polynomials $s_{i}(\mathbf{x},\mathbf{y}),i=0,1,2,3$ in
program (8).

####  IV-A2 CBF Correctness Verification

In addition to verifying the CLF-CBF compatibility, we must also ensure that
the set $\\{\mathbf{x}\;|\;b(\mathbf{x})\geq 0\\}$ is contained in the desired
safe region, i.e., $\\{\mathbf{x}|b(\mathbf{x})\geq 0\\}$ implies
$\mathbf{x}\in\mathcal{X}_{s}$ Assume an unsafe region presented by:

$\mathcal{X}_{u}=\cup_{j=1}^{K}\\{\mathbf{x}|h_{j}(\mathbf{x})\leq 0\\}.$ (9)

Using the S-procedure, we can formulate the correctness verification program:

Find $\displaystyle q_{1}(\mathbf{x})\dots
q_{K}(\mathbf{x}),p_{1}(\mathbf{x})\dots p_{K}(\mathbf{x})\in\text{SOS}$ s.t.
$\displaystyle-1+q_{j}(\mathbf{x})b(\mathbf{x})+p_{j}(\mathbf{x})h_{j}(\mathbf{x})\in\text{SOS},\;j=1,\dots
K.$ (10a)

On the other hand, if the unsafe region is defined as

$\mathcal{X}_{u}=\cap_{j=1}^{K}\\{\mathbf{x}|h_{j}(\mathbf{x})\leq 0\\},$ (11)

then the correctness verification program is given by

Find $\displaystyle q(\mathbf{x}),p_{1}(\mathbf{x})\dots
p_{K}(\mathbf{x})\in\text{SOS}$ s.t.
$\displaystyle-1+q(\mathbf{x})b(\mathbf{x})+\sum_{j=1}^{K}p_{j}(\mathbf{x})h_{j}(\mathbf{x})\in\text{SOS}.$
(12a)

If the unsafe region is only characterized by single polynomial, then (10) and
(12) are the same and we would have single $q(\mathbf{x})$ and
$p(\mathbf{x})$.

Both problems (10) and (12) can be solved

###  IV-B Compatible CBF and CLF synthesis

Based on the previous verification framework, we next propose an alternating
method to synthesize compatible CBF and CLF. We first formulate the synthesis
problem, and then introduce the overall structure of the alternating algorithm
and each component of the algorithm in detail. Finally we prove the
compatibility of the synthesized CLF and CBF. For simplicity, in the following
derivations, we use (9) as the definition of unsafe region and use the
constraints in (10) as CBF correctness constraint.

###### Problem 2

Given a system (1), with equilibrium state $\mathbf{x}^{*}=\mathbf{0}$,
constants $\kappa_{V}$, $\kappa_{b}$, and polynomials $h_{1}(\mathbf{x})\dots
h_{K}(\mathbf{x})$ characterizing the unsafe region, find a positive definite
function $V(\mathbf{x})$ and a function $b(\mathbf{x})$ such that there always
exist a $\mathbf{u}\in\mathcal{U}$ satisfying (4) for all the
$\mathbf{x}\in\mathcal{D}\cap\mathcal{C}$, where
$\mathcal{C}=\\{\mathbf{x}\;|\;b(\mathbf{x})\geq 0\\}$, and
$\mathcal{D}=\\{\mathbf{x}\;|\;V(\mathbf{x})\leq 1\\}$. We call the region
$\mathcal{C}\cap\mathcal{D}$ as the “compatible region”, and our goal is to
maximize some objective function $\gamma(\mathcal{C}\cap\mathcal{D})$ on this
compatible region. One example of the objective function is the volume
$\gamma(\mathcal{C}\cap\mathcal{D})=\text{Vol}(\mathcal{C}\cap\mathcal{D})$.

Mathematically, our goal is to solve the following optimization problem

$\displaystyle\max_{\begin{subarray}{c}V(\mathbf{x}),b(\mathbf{x}),\\\
p(\mathbf{x}),q(\mathbf{x}),s(\mathbf{x})\end{subarray}}\gamma(\mathcal{C}\cap\mathcal{D})$
(13a) s.t $\displaystyle V(\mathbf{0})=0,V(\mathbf{x})\in\text{ SOS}$ (13b)
$\displaystyle\text{Constraint }\eqref{eq:verify_S-
procedure_constraint},\eqref{eq:cbf_safe_constraint}$ (13c) $\displaystyle
s_{2}(\mathbf{x},\mathbf{y}),s_{3}(\mathbf{x},\mathbf{y}),q_{i}(\mathbf{x}),p_{i}(\mathbf{x})\in\text{
SOS},i=1,\ldots,K.$ (13d)

Unlike the verification problem (8) and (10) where $V(\mathbf{x})$ and
$b(\mathbf{x})$ are given, now we need to jointly search for
$V(\mathbf{x}),b(\mathbf{x})$ together with the Lagrangian multipliers
$p(\mathbf{x}),q(\mathbf{x}),s(\mathbf{x})$ in (13). Notice that all the
constraints in problem (13) are convex, except for (8b) and (10a), which are
bi-linear in the decision variable
$V(\mathbf{x}),b(\mathbf{x}),p(\mathbf{x}),q(\mathbf{x}),s(\mathbf{x})$, and
hence cannot be directly handled by convex optimization. To remedy this, we
resort to bilinear alternation, a common technique in SOS-based controller
design [36, 41]. Where we alternate between

  1. 1.

Lagrangian searching: fix $V(\mathbf{x}),b(\mathbf{x})$, search for the
Lagrangian multipliers $p(\mathbf{x}),q(\mathbf{x}),s(\mathbf{x})$.

  2. 2.

Update CLF-CBF: fix the Lagrangian multiplier, search for
$V(\mathbf{x}),b(\mathbf{x})$ while maximizing a surrogate function of
$\gamma(\mathcal{C}\cap\mathcal{D})$.

Step 1 is the verification problem (8) and (10) discussed in the previous
subsection. In step 2, we consider two different choices of objective function
$\gamma(\mathcal{C}\cap\mathcal{D})$

  * •

Maximize the volume $\text{Vol}(\mathcal{C}\cap\mathcal{D})$.

  * •

Given a set of predefined “candidate states”
$\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(L)}$, we want the compatible region
$\mathcal{C}\cap\mathcal{D}$ to cover as many “candidate states” as possible.

We next present how to search for $V(\mathbf{x}),b(\mathbf{x})$ to optimize
these two objectives respectively, while fixing the Lagrangian multipliers.

####  IV-B1 Maximize volume

To maximize the volume of the compatible region, we need to solve the
following optimization problem

$\displaystyle\max_{V(\mathbf{x}),b(\mathbf{x})}\;\text{Vol}(\mathcal{C}\cap\mathcal{D})$
(14a) s.t $\displaystyle V(\mathbf{0})=0,V(\mathbf{x})\in\text{ SOS}$ (14b)
$\displaystyle\text{Constraint }\eqref{eq:verify_S-
procedure_constraint},\eqref{eq:cbf_safe_constraint}.$ (14c)

Directly maximing the volume of the semi-algebraic set
$\text{Vol}(\mathcal{C}\cap\mathcal{D})$ can be challenging333There exists
some work to compute the volume of a given basic semi-algebraic set through a
sequence of semidefinite programs with increasing sizes [42], but cannot be
directly applied to maximize volume of the basic semi-algebraic set while
searching for its defining polynomials., hence we seek an alternative approach
to grow the compatible region, as shown in Appendix VII-A.

####  IV-B2 Cover many candidate states

Alternatively we can grow the compatible region $\mathcal{C}\cap\mathcal{D}$
to cover the candidate states $\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(L)}$. We
know that the candidate state $\mathbf{x}^{(i)}\in\mathcal{C}\cap\mathcal{D}$
iff $V(\mathbf{x}^{(i)})\leq 1$ and $b(\mathbf{x}^{(i)})\geq 0$. Hence we can
design the objective function as

$\displaystyle\min_{V(\mathbf{x}),b(\mathbf{x})}\sum_{i=1}^{L}c_{1}\text{ReLU}(V(\mathbf{x}^{(i)})-1)+c_{2}\text{ReLU}(-b(\mathbf{x}^{(i))}),$
(15)

where $c_{1},c_{2}>0$ are given positive weights, and Eq (15) penalizes the
violation of $x^{(i)}\notin\mathcal{C}\cap\mathcal{D}$. To avoid arbitrary
scaling of CBF $b(\mathbf{x})$ (if $b(\mathbf{x})$ is a valid CBF, the scaling
it as $k*b(\mathbf{x})$ by any positive constant $k$ is still a valid CBF, and
the cost function becomes unbounded under infinite scaling), we impose an
additional constraint $b(\mathbf{0})\leq 1$. As a result, the optimization
problem is written as

$\displaystyle\text{Objective }\eqref{eq:candidate_states_cost}$ (16a) s.t
$\displaystyle V(\mathbf{0})=0,V(\mathbf{x})\in\text{ SOS}$ (16b)
$\displaystyle\text{Constraint }\eqref{eq:verify_S-
procedure_constraint},\eqref{eq:cbf_safe_constraint}$ (16c) $\displaystyle
b(\mathbf{0})\leq 1$ (16d)

The optimization problems (14) and (16) are convex optimization problems when
searching over $V(\mathbf{x}),b(\mathbf{x})$.

To summarize, we outline our synthesis algorithm in Algorithm 1

Algorithm 1 Synthesize compatible CLF-CBF

1:Given system dynamics (1), positive constants $\kappa_{V},\kappa_{b}$ and
unsafe region $\mathcal{X}_{u}$ in (9).

2:Initialize $V^{(0)}(\mathbf{x}),b^{(0)}(\mathbf{x})$. $i=0$

3:while $i<\text{MaxIter}$ do

4: Fix $V^{(i)}(\mathbf{x}),b^{(i)}(\mathbf{x})$, solve (8) and (10) to find
the Lagrangian multipliers
$p^{(i)}(\mathbf{x}),q^{(i)}(\mathbf{x}),s^{(i)}(\mathbf{x},\mathbf{y})$.

5: Fix
$p^{(i)}(\mathbf{x}),q^{(i)}(\mathbf{x}),s^{(i)}(\mathbf{x},\mathbf{y})$, find
CLF-CBF $V^{(i+1)}(\mathbf{x}),b^{(i+1)}(\mathbf{x})$ through solving either
(14) or (16).

6: $i\leftarrow i+1$.

7:end while

In line 2 of Algorithm 1, we can initialize the CLF/CBF in many ways, for
example a good candidate is to linearize the dynamics around the equilibrium
state and use the LQR cost-to-go to design the initial guess.

With a feasible initial guess $V^{(0)}$ and $b^{(0)}$, by induction, the CLF
$V^{(i)}$ and CBF $b^{(i)}$ are also compatible in each iteration.

##  V Numerical Simulation

In this section, we evaluate our proposed exact verification method and
compatible CLF/CBF synthesis approach. We demonstrate that our verification
can differentiate compatible and incompatible CLF/CBF by solving convex
optimization problems. We show that our proposed synthesis approach construct
CLF/CBF candidates by enlarging the compatible region while ensuring the
compatibility. All case studies are written in Python based on Drake [43] and
MOSEK [44], and are conducted on a laptop with an Intel Core i5 6300HQ CPU and
8GB RAM. The detailed parameter settings can be found in our code.

###  V-A Verification

We evaluate the proposed verification framework using the 2-D first order
integrator example mentioned in [32]. In order to show the framework’s
capability of differentiating incompatible and compatible CLF-CBF, we first
construct an incompatible CBF and CLF and show that they fail to pass the
verification. We then verify a compatible CBF/CLF pair.

We evaluate our verification approach on the 2D first-order integrator
$\dot{\mathbf{x}}=\mathbf{u}$, where $\mathbf{x}\in\mathbb{R}^{2}$,
$\mathbf{u}\in\mathcal{U}\subset\mathbb{R}^{2}$. The unsafe region
$\mathcal{X}_{u}$ is represented by
$\mathcal{X}_{u}=\\{\mathbf{x}\;|\;\mathbf{x}_{1}^{2}+(\mathbf{x}_{2}-3)^{2}-1.5^{2}\leq
0\\}$. We also consider the linear input constraint
$\mathcal{U}=\\{\mathbf{u}\;|\;|\mathbf{u}_{1}|\leq 3,|\mathbf{u}_{2}|\leq
3\\}$. We let parameters $\kappa_{V}=0.01$ and $\kappa_{b}=0.01$.

####  V-A1 Incompatible Example

In this example, we use the following given pair of incompatible CLF and CBF:

$\displaystyle\text{CLF:
}V(\mathbf{x})=0.05\mathbf{x}_{1}^{2}+0.0125\mathbf{x}_{2}^{2}$
$\displaystyle\text{CBF:
}b(x)=\mathbf{x}_{1}^{2}+(\mathbf{x}_{2}-3)^{2}-1.5^{2}.$ (17)

The incompatibility is checked by finding a state $\hat{\mathbf{x}}=[0,5]^{T}$
at which there does not exist $\mathbf{u}\in\mathcal{U}$ satisfying CLF and
CBF constraints. We then give this pair of CLF-CBF to the verification SOS
program. The degrees of the Lagrangian polynomials are set to be $2$. The
convex optimization (8) reports infeasibility, which is consistent with the
incompatible state $[0,5]^{T}$.

####  V-A2 Compatible Example

We use the idea of CLBF [3] and give the verification framework a CLBF pair:

$\displaystyle\text{CLF: }V(\mathbf{x})=\mathbf{x}_{1}^{2}+\mathbf{x}_{2}^{2}$
$\displaystyle\text{CBF: }b(\mathbf{x})=1-V(\mathbf{x}).$ (18)

In this example, the polynomial degrees for compatible Lagrangian polynomials
are set to be $2$. The solver solves the verification SOS program
successfully, which certifies the compatibility of the CLF and CBF in (V-A2)
are compatible, and the given CBF is correct.

The results for these two examples show that our proposed verification
approach can verify compatibility of CBF-CLF, while also differentiating
between compatible and incompatible CBF-CLF pairs.

###  V-B Synthesis

We next evaluate our proposed bilinear alternation synthesis method. We
investigate a linear example, a nonlinear example [45] and a power converter
example [33]. We set parameters $\kappa_{V}=0.01$, and $\kappa_{b}=0.01$.
Since the compatible region may infinitely enlarged if $\mathcal{X}_{s}$ is an
open set, then the synthesis loop may not terminate. Hence, we set maximum
number of iterations as a hard stopping condition (MaxIter=25) during
experiments.

####  V-B1 Linear Example

For the linear toy example, we use the following system dynamics:

$\begin{bmatrix}\dot{\mathbf{x}}_{1}\\\
\dot{\mathbf{x}}_{2}\end{bmatrix}=\begin{bmatrix}1&2\\\
-2&3\end{bmatrix}\begin{bmatrix}\mathbf{x}_{1}\\\
\mathbf{x}_{2}\end{bmatrix}+\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\begin{bmatrix}\mathbf{u}_{1}\\\ \mathbf{u}_{2}\end{bmatrix}.$
(19)

The obstacle is represented as
$\mathcal{X}_{u}=\\{\mathbf{x}\;|\;\mathbf{x}_{1}+\mathbf{x}_{2}+3\leq 0\\}$ .
We construct a compatible CBF and CLF using a heuristic method that solves a
linear quadratic regulator (LQR) problem in a neighborhood of an equilibrium
point to initialize CLF/CBF. We set the cost matrices of LQR as
$Q=I_{2},R=I_{2}$, where $I_{2}$ is $2\times 2$ identity matrix.

(a)

(b)

(c)

Figure 1: This figure presets the compatible CLF-CBF synthesis results for all
the investigated examples. In these plots, unsafe regions (UR) are marked as
red, initial compatible regions (ICR) are marked as gray, initial compatible
boundaries (ICB) are marked as dashed lines, synthesized compatible regions
(SCB) are marked as green, and synthesized compatible boundaries (SCB) are
marked as solid lines. The compatible regions are enlarged for all the
investigated examples. Figure 1(a) and 1(b) show the results for linear
example and nonlinear example, respectively. In these two examples, the
synthesis loop stops after maximum number of iterations (25 iterations).
Figure 1(c) shows the synthesis result for power converter example. The
synthesis stops after 7 iterations because the compatible region can not be
enlarged further according to optimization (16)

Figure 1(a) shows the initial compatible region and the compatible region
after synthesis. We observe that our alternating synthesis method can enlarge
the compatible region and incorporates the initial compatible region. The
synthesis loop stops after maximum number of iterations.

####  V-B2 Nonlinear Example

In this example, the following dynamics [21] is used:

$\begin{bmatrix}\dot{\mathbf{x}}_{1}\\\
\dot{\mathbf{x}}_{2}\end{bmatrix}=\begin{bmatrix}0\\\
-\mathbf{x}_{1}+\frac{\mathbf{x}_{1}^{3}}{6}\end{bmatrix}+\begin{bmatrix}1\\\
-1\end{bmatrix}\begin{bmatrix}u\end{bmatrix}.$ (20)

We pick $\mathbf{x}_{eq}=[0,0]^{T},u_{eq}=0$ as the equilibrium point and
linearize the system around this point to compute the LQR. The obstacle for
this example is the same as linear toy example. The cost matrices of LQR are
set as $Q=I_{2},R=1$.The synthesis results for this example are shown by
Figure 1(b). We observe that the proposed alternating synthesis can enlarge
the compatible region in this example.

####  V-B3 Power Converter Example

The dynamical model of the power converter is given in [33] as follows

$\begin{bmatrix}\dot{\mathbf{x}}_{1}\\\ \dot{\mathbf{x}}_{2}\\\
\dot{\mathbf{x}}_{3}\end{bmatrix}=\begin{bmatrix}-0.05\mathbf{x}_{1}-57.9\mathbf{x}_{2}+0.00919\mathbf{x}_{3}\\\
1710\mathbf{x}_{1}+314\mathbf{x}_{3}\\\
-0.271\mathbf{x}_{1}-314\mathbf{x}_{2}\end{bmatrix}+\\\
\begin{bmatrix}0.05-57.9\mathbf{x}_{2}&-57.9\mathbf{x}_{3}\\\
1710+1710\mathbf{x}_{1}&0\\\
0&1710+1710\mathbf{x}_{1}\end{bmatrix}\begin{bmatrix}\mathbf{u}_{1}\\\
\mathbf{u}_{2}\end{bmatrix}.$ (21)

The unsafe region is presented as:

$\displaystyle\mathcal{X}_{u}=\\{\mathbf{x}|$
$\displaystyle-\mathbf{x}_{1}+0.2\leq 0,\mathbf{x}_{1}+0.8\leq 0,$
$\displaystyle-(\mathbf{x}_{2}-0.001)^{2}-\mathbf{x}_{3}^{2}+1.2^{2}\leq
0\\}.$ (22)

The equilibrium point of the system is
$\mathbf{x}_{eq}=[0,0,0]^{T},\mathbf{u}_{eq}=[0,0]^{T}$. We linearize the
system near the equilibrium point and compute the LQR. We set the cost
matrices of LQR as $Q=I_{3},R=I_{2}$, where $I_{3}$ is $3\times 3$ identity
matrix. The initialized and synthesised compatible regions for power converter
example are shown by Figure 1(c). The proposed synthesis method enlarges the
compatible region for the power converter. The syntehsis loop stops after 7
number of iterations because the compatible region can not be enlarged further
according to optimization problem (16). The returned CLF-CBF are verified to
be compatible. Together with Figure 1(a) and Figure 1(b), we can conclude that
the proposed alternating synthesis method enlarges the compatible region and
returns compaitble CLF-CBF for all the investigated examples.

##  VI Conclusion

In this paper, for control affine systems both with and without input limits,
we studied the verification and synthesis of compatible CBF and CLF
independent from any nominal controllers. We presented sufficient and
necessary conditions for compatible CBF and CLF via Farkas’ Lemma and
algebraic-geometric Positivstellensatz, and demonstrated that these conditions
can be checked by solving an SOS program. In order to reduce the number of SOS
polynomials in the SOS program, we further simplified the necessary and
sufficient condition using S-procedure and formulated a more tractable SOS
program for compatibility verification. We next proposed an alternating method
to synthesize compatible CLF and CBFs based on the verification framework. In
each iteration, the algorithm first takes the CLF and CBF from the previous
iteration and computes the Lagrangian polynomials by solving the verification
SOS program. Then, it takes the Lagrangian polynomials and updates the CLF and
CBF to obtain a larger compatible region. Finally, we evaluated our proposed
approach with a linear system, a nonlinear system and nonlinear power
converter. In our future work, we will extend the verification and synthesis
framework to high-relative degree and non-convex scenarios.

## References

  * [1] Department of Homeland Security, Department of Homeland Security Cyber-Physical Systems Page, Jan 2022.  [Online] Available: https://www.dhs.gov/science-and-technology/cpssec [Accessed: Jan 2022].
  * [2] B. Li, S. Wen, Z. Yan, G. Wen, and T. Huang, “A survey on the control Lyapunov function and control barrier function for nonlinear-affine control systems,” IEEE/CAA Journal of Automatica Sinica, vol. 10, no. 3, pp. 584–602, 2023.
  * [3] R. K. Fachri, M. Z. Romdlony, and M. R. Rosa, “Multiple waypoint navigation for mobile robot using control Lyapunov-barrier function (CLBF),” in 2022 IEEE International Conference on Cybernetics and Computational Intelligence (CyberneticsCom), pp. 230–235, IEEE, 2022.
  * [4] M. H. Cohen, T. G. Molnar, and A. D. Ames, “Safety-critical control for autonomous systems: Control barrier functions via reduced-order models,” Annual Reviews in Control, vol. 57, p. 100947, 2024.
  * [5] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in 53rd IEEE Conference on Decision and Control, pp. 6271–6278, IEEE, 2014.
  * [6] X. Fan, Y. Guo, H. Liu, B. Wei, and W. Lyu, “Improved artificial potential field method applied for AUV path planning,” Mathematical Problems in Engineering, vol. 2020, pp. 1–21, 2020.
  * [7] X. Chen and S. Sankaranarayanan, “Reachability analysis for cyber-physical systems: Are we there yet?,” in NASA Formal Methods Symposium, pp. 109–130, Springer, 2022.
  * [8] D. Lee and C. J. Tomlin, “A Hopf-Lax formula in hamilton–jacobi analysis of reach-avoid problems,” IEEE Control Systems Letters, vol. 5, no. 3, pp. 1055–1060, 2020.
  * [9] M. Schwenzer, M. Ay, T. Bergs, and D. Abel, “Review on model predictive control: An engineering perspective,” The International Journal of Advanced Manufacturing Technology, vol. 117, no. 5, pp. 1327–1349, 2021.
  * [10] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in 2019 18th European control conference (ECC), pp. 3420–3431, IEEE, 2019.
  * [11] S. Prajna, A. Jadbabaie, and G. J. Pappas, “A framework for worst-case and stochastic safety verification using barrier certificates,” IEEE Transactions on Automatic Control, vol. 52, no. 8, pp. 1415–1428, 2007.
  * [12] A. M. Lyapunov, “The general problem of the stability of motion,” International journal of control, vol. 55, no. 3, pp. 531–534, 1892.
  * [13] E. D. Sontag, “A Lyapunov-like characterization of asymptotic controllability,” SIAM journal on control and optimization, vol. 21, no. 3, pp. 462–471, 1983.
  * [14] J. Zeng, B. Zhang, and K. Sreenath, “Safety-critical model predictive control with discrete-time control barrier function,” in 2021 American Control Conference (ACC), pp. 3882–3889, IEEE, 2021.
  * [15] M. H. Cohen and C. Belta, “High order robust adaptive control barrier functions and exponentially stabilizing adaptive control Lyapunov functions,” in 2022 American Control Conference (ACC), pp. 2233–2238, IEEE, 2022.
  * [16] J. Choi, F. Castaneda, C. J. Tomlin, and K. Sreenath, “Reinforcement learning for safety-critical control under model uncertainty, using control Lyapunov functions and control barrier functions,” arXiv preprint arXiv:2004.07584, 2020.
  * [17] C. Dawson, S. Gao, and C. Fan, “Safe control with learned certificates: A survey of neural Lyapunov, barrier, and contraction methods for robotics and control,” IEEE Transactions on Robotics, 2023.
  * [18] S. Liu, C. Liu, and J. Dolan, “Safe control under input limits with neural control barrier functions,” in Conference on Robot Learning, pp. 1970–1980, PMLR, 2023.
  * [19] H. Wang, K. Margellos, and A. Papachristodoulou, “Safety verification and controller synthesis for systems with input constraints,” IFAC-PapersOnLine, vol. 56, no. 2, pp. 1698–1703, 2023.
  * [20] S. Kang, Y. Chen, H. Yang, and M. Pavone, “Verification and synthesis of robust control barrier functions: Multilevel polynomial optimization and semidefinite relaxation,” arXiv preprint arXiv:2303.10081, 2023.
  * [21] H. Dai and F. Permenter, “Convex synthesis and verification of control-Lyapunov and barrier functions with input constraints,” in 2023 American Control Conference (ACC), pp. 4116–4123, IEEE, 2023.
  * [22] A. Clark, “A semi-algebraic framework for verification and synthesis of control barrier functions,” arXiv preprint arXiv:2209.00081, 2022.
  * [23] W. Zhao, T. He, T. Wei, S. Liu, and C. Liu, “Safety index synthesis via sum-of-squares programming,” in 2023 American Control Conference (ACC), pp. 732–737, IEEE, 2023.
  * [24] S. Chen and M. Fazlyab, “Learning performance-oriented control barrier functions under complex safety constraints and limited actuation,” arXiv preprint arXiv:2401.05629, 2024.
  * [25] H. Zhang, Z. Li, H. Dai, and A. Clark, “Efficient sum of squares-based verification and construction of control barrier functions by sampling on algebraic varieties,” in 2023 62nd IEEE Conference on Decision and Control (CDC), pp. 5384–5391, 2023.
  * [26] S. Yang, Y. Chen, X. Yin, and R. Mangharam, “Learning local control barrier functions for safety control of hybrid systems,” arXiv preprint arXiv:2401.14907, 2024.
  * [27] B. Dai, P. Krishnamurthy, and F. Khorrami, “Learning a better control barrier function,” in 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 945–950, IEEE, 2022.
  * [28] W. Jin, Z. Wang, Z. Yang, and S. Mou, “Neural certificates for safe control policies,” arXiv preprint arXiv:2006.08465, 2020.
  * [29] P. Zhao, R. Ghabcheloo, Y. Cheng, H. Abdi, and N. Hovakimyan, “Convex synthesis of control barrier functions under input constraints,” IEEE Control Systems Letters, 2023.
  * [30] A. Isaly, M. Ghanbarpour, R. G. Sanfelice, and W. E. Dixon, “On the feasibility and continuity of feedback controllers defined by multiple control barrier functions,” IEEE Transactions on Automatic Control, pp. 1–15, 2024.
  * [31] H. Wang, K. Margellos, and A. Papachristodoulou, “Relaxed compatibility between control barrier and lyapunov functions,” in 2024 UKACC 14th International Conference on Control (CONTROL), pp. 125–126, 2024.
  * [32] M. F. Reis, A. P. Aguiar, and P. Tabuada, “Control barrier function-based quadratic programs introduce undesirable asymptotically stable equilibria,” IEEE Control Systems Letters, vol. 5, no. 2, pp. 731–736, 2020.
  * [33] M. Schneeberger, S. Mastellone, and F. Dörfler, “Advanced safety filter based on SOS control barrier and Lyapunov functions,” arXiv preprint arXiv:2401.06901, 2024.
  * [34] M. Schneeberger, F. Dörfler, and S. Mastellone, “SOS construction of compatible control Lyapunov and barrier functions,” arXiv preprint arXiv:2305.01222, 2023.
  * [35] A. Papachristodoulou, J. Anderson, G. Valmorbida, S. Prajna, P. Seiler, P. Parrilo, M. M. Peet, and D. Jagt, “SOSTOOLS version 4.00 sum of squares optimization toolbox for matlab,” arXiv preprint arXiv:1310.4716, 2013.
  * [36] P. A. Parrilo, Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization.  California Institute of Technology, 2000.
  * [37] M. Li and Z. Sun, “A graphical interpretation and universal formula for safe stabilization,” in 2023 American Control Conference (ACC), pp. 3012–3017, IEEE, 2023.
  * [38] G. Blekherman, P. A. Parrilo, and R. R. Thomas, Semidefinite Optimization and Convex Algebraic Geometry.  SIAM, 2012.
  * [39] J. Matousek and B. Gärtner, Understanding and Using Linear Programming.  Springer Science & Business Media, 2006.
  * [40] E. D. Sontag, “A ‘universal’construction of artstein’s theorem on nonlinear stabilization,” Systems & control letters, vol. 13, no. 2, pp. 117–123, 1989.
  * [41] R. Tedrake, I. R. Manchester, M. Tobenkin, and J. W. Roberts, “LQR-trees: Feedback motion planning via sums-of-squares verification,” The International Journal of Robotics Research, vol. 29, no. 8, pp. 1038–1052, 2010.
  * [42] D. Henrion, J. B. Lasserre, and C. Savorgnan, “Approximate volume and integration for basic semialgebraic sets,” SIAM review, vol. 51, no. 4, pp. 722–743, 2009.
  * [43] R. Tedrake and the Drake Development Team, “Drake: Model-based design and verification for robotics,” 2019.
  * [44] M. ApS, MOSEK Optimization Software Python API Documentation, 2024.
  * [45] W. Tan and A. Packard, “Searching for control Lyapunov functions using sums of squares programming,” sibi, vol. 1, no. 1, 2004.

##  VII Appendix

###  VII-A Solving optimization problem (14)

The optimization problem (14) aims at maximizing the volume of the compatible
region $\mathcal{C}\cap\mathcal{D}$, which is a complicated basic semi-
algebraic set. It is challenging to maximize the volume of the basic semi-
algebraic set while searching for its defining polynomials. To remedy this, we
consider an alternative approach which doesn’t compute the volume directly.

In Algorithm 1 we iteratively search for the new CLF/CBF (together with the
new compatible region). Our intuition is to enlarge the compatible region,
such that the new one contains the old compatible region in the previous
iteration. Formally, in iteration $i$, given the CLF/CBF in the previous
iteration as $V^{(i)}(\mathbf{x}),b^{(i)}(\mathbf{x})$ (whose compatible
region is $\mathcal{C}^{(i)}\cap\mathcal{D}^{(i)}$), together with the
Lagrangian multiplier
$p^{(i)}(\mathbf{x}),q^{(i)}(\mathbf{x}),s^{(i)}(\mathbf{x},\mathbf{y})$, our
goal is to find a new CLF $V^{(i+1)}$ and CBF $b^{(i+1)}(\mathbf{x})$. The new
compatible region $\mathcal{C}^{(i+1)}\cap\mathcal{D}^{(i+1)}$ will contain
the old compatible region, namely

$\displaystyle\begin{split}\forall\mathbf{x}\in\\{\mathbf{x}\;|\;b^{(i)}(\mathbf{x})\geq
0,1-V^{(i)}(\mathbf{x})\geq 0\\}\\\ \qquad\Rightarrow
b^{(i+1)}(\mathbf{x})\geq t_{b}\text{ and }1-V^{(i+1)}(\mathbf{x})\geq
t_{V}\end{split}.$ (23)

where we use the positive scalar $t_{b}$ and $t_{V}$ to measure the “margin”
between the old and new compatible regions, and our goal is to maximize this
margin. Similar to the (16), to avoid arbitrary scaling of $b(\mathbf{x})$, we
introduce the constraint $b(\mathbf{0})\leq 1$ to avoid the infinite scaling.

We can convert the condition (23) to a SOS constraint through the S-procedure
as

$\displaystyle\begin{split}b^{(t+1)}(\mathbf{x})-t_{b}-\phi_{1}(\mathbf{x})b^{(i)}(\mathbf{x})-\\\
\;\phi_{2}(\mathbf{x})(1-V^{(i)}(\mathbf{x}))\in\text{ SOS}\end{split}$ (24a)
$\displaystyle\begin{split}1-V^{(i+1)}(\mathbf{x}-t_{V}-\psi_{1}(\mathbf{x})b^{(i)}(\mathbf{x})-\\\
\psi_{2}(\mathbf{x})(1-V^{(i)}(\mathbf{x}))\in\text{ SOS}\end{split}$ (24b)
$\displaystyle\phi_{1}(\mathbf{x}),\phi_{2}(\mathbf{x}),\psi_{1}(\mathbf{x}),\psi_{2}(\mathbf{x})\in\text{
SOS},$ (24c)

and obtain the following optimization problem

$\displaystyle\max_{\begin{subarray}{c}V^{(i+1)}(\mathbf{x}),b^{(i+1)}(\mathbf{x}),\\\
t_{b},t_{V},\\\ \phi_{1}(\mathbf{x}),\phi_{2}(\mathbf{x}),\\\
\psi_{1}(\mathbf{x}),\psi_{2}(\mathbf{x})\end{subarray}}t_{b}+t_{V}$ (25a) s.t
$\displaystyle\;V^{(i+1)}(\mathbf{0})=0,V^{(i+1)}(\mathbf{x})\in\text{ SOS}$
(25b) $\displaystyle\text{Constraint }\eqref{eq:verify_S-
procedure_constraint},\eqref{eq:cbf_safe_constraint}$ (25c)
$\displaystyle\text{Constraint }\eqref{eq:compatible_region_containment_sos}$
(25d) $\displaystyle t_{b}\geq 0,t_{V}\geq 0$ (25e) $\displaystyle
b^{(i+1)}(\mathbf{0})\leq 1,$ (25f)

where we add the constraint $b^{(i+1)}(\mathbf{0})\leq 1$ to avoid arbitrary
scaling of $b^{(i+1)}(\mathbf{x})$, which would make the optimization problem
unbounded with infinite scaling. The optimization problem (25) is a SOS
program.


# Quantum kernels with squeezed-state encoding for machine learning

Long Hin Li<EMAIL_ADDRESS>Guangdong-Hong Kong Joint Laboratory of
Quantum Matter, Department of Physics, and HKU-UCAS Joint Institute for
Theoretical and Computational Physics at Hong Kong, The University of Hong
Kong, Pokfulam Road, Hong Kong, China  Dan-Bo Zhang<EMAIL_ADDRESS>Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Frontier Research
Institute for Physics, South China Normal University, Guangzhou 510006, China
Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum
Materials, School of Physics and Telecommunication Engineering, South China
Normal University, Guangzhou 510006, China  Z. D. Wang<EMAIL_ADDRESS>Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Department of Physics,
and HKU-UCAS Joint Institute for Theoretical and Computational Physics at Hong
Kong, The University of Hong Kong, Pokfulam Road, Hong Kong, China

###### Abstract

Kernel methods are powerful for machine learning, as they can represent data
in feature spaces that similarities between samples may be faithfully
captured. Recently, it is realized that machine learning enhanced by quantum
computing is closely related to kernel methods, where the exponentially large
Hilbert space turns to be a feature space more expressive than classical ones.
In this paper, we generalize quantum kernel methods by encoding data into
continuous-variable quantum states, which can benefit from the infinite-
dimensional Hilbert space of continuous variables. Specially, we propose
squeezed-state encoding, in which data is encoded as either in the amplitude
or the phase. The kernels can be calculated on a quantum computer and then are
combined with classical machine learning, e.g. support vector machine, for
training and predicting tasks. Their comparisons with other classical kernels
are also addressed. Lastly, we discuss physical implementations of squeezed-
state encoding for machine learning in quantum platforms such as trapped ions.

##  I Introduction

Machine learning can be enhanced by quantum computing utilizing the incredible
power of quantum computers for information processing [1, 2, 3, 4, 5, 6, 7, 8,
9, 10]. The first breakthrough comes from an exponential quantum speedup for
solving linear equations or matrix inversion [11, 12], which lies at the heart
of many machine learning methods. While quantum computers are good at linear
algebra since they in nature perform unitary transformations that are linear
in nature, it is not enough for machine learning of complicated data, where
highly nonlinear transformations are often anticipated [13]. Another important
ingredient for quantum-enhanced machine learning, only recognized very
recently, is to introduce non-linearity when encoding classical data into
quantum states [14, 7, 8, 15, 16, 17, 18]. Such a procedure is called quantum
feature maps [7, 8], which enjoy the exponential large Hilbert space as a
feature space for expressing data.

The quantum feature maps are closely related to classical kernel methods in
the same spirit as they aim to use feature spaces, such that similarities
between samples can be faithfully captured [19, 20]. One striking difference
is that feature maps on a quantum computer are explicitly constructed [18].
This gives the flexibility of designing or even training powerful quantum
kernels that may provide quantum advantages for machine learning. Varies
encoding schemes for quantum kernels have been explored. Typically, data
vectors are encoded as parameters of quantum gates, such as rotation angles of
single-qubit gates. By loading the same data repeatedly into different quantum
gates, the quantum data state will be a nonlinear function of the data, which
consists of both different high-frequency and low-frequency components [17].
Those types of quantum kernels are rather different from classical ones, and
while useful applications can be expected, they would have limitations(e.g.,
similarity between two samples may not be nicely approximated as a periodic
function of data vector). On the other hand, quantum computing and quantum
machine learning with continuous variables can inherit many properties of
classical machine learning [21, 22, 8, 23, 24], e.g., the widely-applied
Gaussian kernel can be realized on a quantum computer by coherent-state
encoding. Moreover, squeezed-state encoding has been proposed in Ref. [8].
This suggests that quantum feature maps by encoding data as continuous-
variable states can be promising, by exploiting the infinite dimensionality of
continuous-variable for encoding high-density information.

In this paper, we focus on one important class of continuous-variable states,
the squeezed state, as a quantum kernel for machine learning. We first
introduce the squeezing vacuum state, and two different encoding schemes, the
amplitude encoding and the phase encoding. Those different data encoding
methods can give different kernels, which have their unique properties for
modeling data. The kernel can be calculated as an inner product between two
quantum states that encode the two data, which further serves as input for
machine learning models, such as support vector machines. We will then test
those kernels and also compare them with some classical kernels on different
datasets classification tasks. Finally, we propose a physical platform for
implementation and also discuss further possibilities to expand the kernel
libraries.

##  II Quantum kernels with squeezed-state encoding

Let us first introduce some notations. For a sample represented by a data
vector $x$, its quantum feature map is realized by encoding the data into a
quantum state $\ket{\psi_{x}}$. The quantum kernel of two samples $x$ and
$x^{\prime}$ is an inner product between their corresponding quantum states,
namely $K(x,x^{\prime})=\braket{\psi_{x}}{\psi_{x^{\prime}}}$.

A squeezed vacuum state is a minimum uncertainty state satisfying
$\Delta\hat{x}\Delta\hat{p}=\frac{1}{4}$, where $\hat{x}$ and $\hat{p}$ and
conjugate quadrature operators with either quadrature variance squeezed [25].
It can be obtained from the vacuum state $\ket{0}$,

$\displaystyle\ket{z}_{s}=\hat{S}(z)\ket{0}$ (1)

with squeezing operator defined as
$\hat{S}(z)=e^{\frac{1}{2}(z^{*}\hat{a}^{2}-z\hat{a}^{\dagger 2})}$. The
position uncertainty is ”squeezed” if the amplitude $r$ of the squeezing
parameter $z=re^{i\theta}$ is positive. The momentum uncertainty will be
enlarged by the same proportion to keep their product constant. The squeezed
state can be expanded in terms of Fock state basis by using the property
$\hat{S}(z)\hat{a}\hat{S}^{\dagger}(z)=\hat{a}\cosh{r}+\hat{a}^{\dagger}e^{i\theta}\sinh{r}$,

$\displaystyle\ket{z}_{s}=\frac{1}{\sqrt{\cosh{r}}}\sum_{n=0}^{\infty}(-1)^{n}\frac{\sqrt{(2n)!}}{2^{n}n!}e^{in\theta}\tanh^{n}{r}\ket{2n}$
(2)

Here we introduce the notation separating the amplitude and phase of the
squeezed vacuum state $\ket{z}=\ket{r,\theta}$. The inner product between
squeezed vacuum states is

$\displaystyle\braket{r,\theta}{r^{\prime},\theta^{\prime}}_{s}=\sqrt{\frac{\operatorname{sech}{r}\operatorname{sech}{r^{\prime}}}{1-e^{i(\theta^{\prime}-\theta)}\tanh{r}\tanh{r^{\prime}}}}$
(3)

The inner product will serve as a measure of similarity between data,
equivalent to the kernel in machine learning [20]. This provides a bridge from
quantum computing to machine learning. In many instances, we use basic
encoding $\ket{\psi_{x}}=\ket{x}$ with $\ket{x}$ a tensor product of binary
qubits or amplitude encoding $\ket{\psi_{x}}=\sum_{n=0}^{N}x_{n}\ket{n}$ to
encode data in quantum algorithms. Apart from the discrete variable qubit,
encoding with coherent states are also possible [24]. Different feature maps
are able to capture different data patterns. And with the greater amount of
possible feature maps, more data patterns can be approximated. For example,
the amplitude encoding gives the usual dot product and basic encoding on
coherent state gives Gaussian kernel.

###  II.1 Phase encoding

Given an $N$-dimensional dataset $\\{{\bf x}_{m}\\}$, we can encode data into
the phase parameters of $N$ squeezed vacuum modes with the amplitude as a
hyperparameter [8]. The encoding was $\phi:x\rightarrow\ket{c,x}_{s}$, with
$\ket{c,x}_{s}$ defined in Eq. (2) and is referred as squeezing feature map
with phase encoding. We call this kernel the squeezing phase kernel. The
kernel plot for different data entries is a periodic function, as shown in
Fig. 1. The similarity between data valued at multiples of $2\pi$ separations
is high. This shares similar properties to the exponential sine squared kernel
in machine learning

$\displaystyle
K_{ess}(x,x^{\prime})=\exp\left(-\frac{2}{l^{2}}\sin^{2}\left(|x-x^{\prime}|\pi/p\right)\right),$
(4)

where $l$, $p$ are hyperparameters of the length scale and the period,
respectively. These kinds of kernels are suitable for data that has features
of periodicity. For example, the weather forecasting data and the stock market
data. From Fig. 1, the squeezing phase kernel has sharper peaks compared to
the exponential sine squared kernel, and may be more suitable for modeling
functions with sudden changes.

Figure 1: Comparison between exponential sine squared kernel for different $l$
from Eq. 4 (dash line) and squeezing phase kernel for different $c$ from Eq. 3
(solid line).

In addition, for coherent-state encoding [24], data can also be encoded into
the argument of coherent state that the kernel becomes (the real part of the
inner product),

$\displaystyle\braket{c,x}{c,x^{\prime}}_{c}$
$\displaystyle=\exp\left(-c^{2}+c^{2}e^{i(x^{\prime}-x)}\right)$ (5)
$\displaystyle\Re{\braket{c,x}{c,x^{\prime}}_{c}}$
$\displaystyle=e^{-c^{2}+c^{2}\cos{(x^{\prime}-x)}}\cos{(c^{2}\sin{(x^{\prime}-x)})}$
(6)

The hyperparameter $c$ has the similar function as the length scale $l$ in the
exponential sine squared kernel. We call this kernel the coherent phase
kernel. Although there are no parameters for tuning the desired period, this
can always be done by rescaling the data with some multiplication factor.
Therefore, both phase encoding of squeezing and coherent feature maps are
periodic kernels, which can be used to model functions that are periodic,
similar to the exponential sine squared kernel.

###  II.2 Amplitude encoding

Another way to encode into squeezed states is to encode the data into the
amplitude of squeezing parameter $\phi:x\rightarrow\ket{x,0}$, where the phase
is set as zero. There is no hyperparameter in this case since the exponential
factor vanishes for every equal value in the squeezing vacuum inner product
(Eq.(3)) and we set it as zero. The resulting kernel is

$\displaystyle\braket{x,0}{x^{\prime},0}_{s}=\sqrt{\frac{\operatorname{sech}{x}\operatorname{sech}{x^{\prime}}}{1-\tanh{x}\tanh{x^{\prime}}}}$
(7)

We call this squeezing amplitude kernel. The kernel is symmetric and positive
semidefinite, and has a shape similar to the Gaussian kernel, as shown in Fig.
2. The inner product formed by the basis encoding of the coherent state into
the amplitude has been shown to have the similar functional form of Gaussian
kernel [24], also with no extra hyperparameter,

$\displaystyle\braket{x,0}{x^{\prime},0}_{c}=\exp\left(-\frac{1}{2}|x-x^{\prime}|^{2}\right)$
(8)

The squeezing kernel has a larger spread compared to the Gaussian kernel. It
can be used to model functions with larger variance, which varies more slowly
with respect to the distance from a center point.

Figure 2: Comparison between Gaussian kernel (left) and squeezing amplitude
kernel (right).

##  III Machine learning with quantum kernels

We now apply quantum kernels with squeezed-state encoding for machine
learning. The kernels are obtained through a quantum circuit of evaluating the
inner product between two quantum states, and then serve as the input for the
support vector machine for training and prediction. We further test those
quantum kernels on some standard datasets and compare the performances with
those of classical kernels.

###  III.1 Support vector machine with quantum kernels

There are two main streams of quantum algorithms where the kernel method
applies. The first kind is solely based on quantum computers. Famous examples
include linear regression [12, 22], quantum support vector machine [3],
quantum principal component analysis [2]. The other class is the classical-
quantum hybrid approach that a part of computational task is left to classical
devices. In this regime, the optimization of the machine learning algorithm is
usually outsourced to classical computers. Quantum computers are only required
to generate quantum gates that encode data, and then to calculate the inner
product between quantum states. This approach is referred to as quantum kernel
estimation [8], which is easier to achieve than the first approach as both the
circuit depth and the quantum resources required are less. We can use quantum
computers to calculate the proposed quantum kernels above, and the result can
be used by classical computers for training traditional machine learning
models like support vector machines, principal component analysis, clustering,
and Gaussian processes. We take the support vector machine as an example.

Firstly, the inner product can be calculated on a quantum computer in two
ways. The first is to use the swap test. Taking the squeezing kernel as an
example. We need to prepare the quantum state

$\displaystyle\ket{\Psi}=\frac{1}{\sqrt{2}}(\ket{0}\ket{z}_{s}+\ket{1}\ket{z^{\prime}}_{s})$
(9)

then take a $\sigma_{x}$ measurement on the ancillary qubit, to measure the
probability of getting $+1$. The inner product can be calculated by $2p-1$.
Another approach to calculating the inner product is to apply two consecutive
squeezing operators $\hat{S}(z)$ then $\hat{S}^{\dagger}(z^{\prime})$, and
measure the probability that the qumode is in vacuum state $\ket{0}$. This
avoids the construction of states that entangle qubit and qumode in Eq. (9).

Secondly, with kernels at hand, we continue the procedure of machine learning
on classical computers. To optimize an SVM classifier, we have to maximize the
Lagrangian dual function [19, 20]

$\displaystyle
L(\vec{\alpha})=\sum_{i=1}^{M}\alpha_{i}-\frac{1}{2}\sum_{i,j=1}^{M}y_{i}y_{j}\alpha_{i}\alpha_{j}K(x_{i},x_{j})$
(10)

subject to the constraint $\sum_{i=1}^{M}\alpha_{i}y_{i}=0$ and
$\alpha_{i}\geq 0$ for every $i$. Optimized values of $\alpha_{i}$ and $b$ are
derived from Karush-Kuhn-Tucker (KKT) conditions,

$\displaystyle\begin{aligned} y_{i}y(x_{i})&=1\\\ y_{i}\left(\sum_{j\in
S}\alpha_{j}y_{j}K(x_{j},x_{i})+b\right)&=1\end{aligned}$ (11) $\displaystyle
b=\frac{1}{N_{S}}\sum_{i\in S}\left(y_{i}-\sum_{j\in
S}\alpha_{j}y_{j}K(x_{j},x_{i})\right)$ (12)

where $S$ denotes the set of support vectors. Then, prediction can be
calculated by

$\displaystyle y(x)=\sum_{i\in S}\alpha_{i}y_{i}K(x_{i},x)+b$ (13)

Therefore, a support vector classifier can be trained after computing the
kernel for each pair of data $K(x_{j},x_{i})$ on the quantum computer, and
similar procedures apply to other traditional kernel-based machine learning
models.

###  III.2 Results

As an illustration purpose, we test and compare those kernels mentioned in
Sec. II for supervised learning of classification on some standard datasets
from scikit-learn [26]. Those kernels are used to calculate randomly generated
sample datasets of each $200$ data, and the SVM model was trained with the
kernel. Then, another set of $200$ data is used for validation to calculate
the performance score. The results are shown in Fig. 3 for comparing squeezing
amplitude kernel with the standard Gaussian kernel, and in Fig. 4 for
comparing the coherent phase kernel and exponentiated sine squared kernel.
Their corresponding performance scores are listed in Table 1 and Table 2,
respectively. The decision boundary of the squeezing amplitude kernel is
comparable to the boundary created by the standard Gaussian kernel. The
variance of the squeezing amplitude kernel is slightly lower than the Gaussian
kernel but still demonstrated the ability to capture the highly nonlinear
behaviour of the dataset. On the other hand, the decision boundary created for
coherent phase kernel and exponentiated sine squared kernels are different for
each case but they shared the common shape in general. These demonstrated that
the two kernels will be useful in different data patterns.

Figure 3: Comparison of SVM training results between squeezing amplitude kernel (first row) and Gaussian kernel (second row). It was trained with Python’s scikit-learn SVC classifier using customized kernel. The three standard datasets used are ”blobs” (first column), ”circles” (second column) and ”moons” (third column). $400$ data was generated randomly in each case where the training and validating set have $200$ data each. Only the validating data are shown. The accuracy of each training is measured and shown in Table 1. | Blobs | Moons | Circles
---|---|---|---
Squeezing amplitude | 0.915 | 0.95 | 0.795
Gaussian | 0.9 | 1.0 | 0.79
Table 1: The accuracy of validation corresponds to Figure 3.

Figure 4: Comparison of SVM training results between coherent phase kernel (first row) and exponential sine squared kernel (second row). It was trained with Python’s scikit-learn SVC classifier using customized kernel. The three generated datasets used are ”moon” (first column), ”spiral” (second column) and ”sine” (third column). $400$ data was generated randomly in each case where the training and validating set have $200$ data each. Only the validating data are shown. The accuracy of each training is measured and shown in Table 2. The hyperparameters are tuned by searching for optimized parameters. | Moons | Spiral | Sine
---|---|---|---
Coherent phase | 1.0 | 0.995 | 0.96
Exp sine-squared | 1.0 | 1.0 | 0.95
Table 2: The accuracy of validation corresponds to Figure 4.

###  III.3 Physical implementation

The whole machine learning procedure only requires the quantum computer at the
stage of evaluating quantum kernels, and we briefly discuss its physical
implementation on a quantum system. The inner product has to be calculated on
a physical system that allows manipulation of both continuous variables and a
discrete-variable qubit as control. Promising candidates of quantum platforms
capable of hybrid variable quantum computing include superconducting circuit
and trapped ion quantum computer [27, 28, 29, 30, 31, 32]. We take trapped ion
as an example. The motional state of the ion can be used as continuous
variables, and various motional states have been realized experimentally [30,
31, 32]. Remarkably, squeezed vacuum state $\ket{z}_{s}$ can be created by
irradiating two Raman beams with some specific frequency difference. The
squeezing parameter will grow exponentially with driving time for the Raman
transition that occurs during the process. Moreover, internal levels of an ion
can serve as a qubit and its coupling to the motional modes can implement the
quantum circuit to generate the state in Eq. 9 that is central to the
evaluation of quantum kernels.

##  IV Conclusions

We have illustrated the idea of using different continuous-variable quantum
states to generate new nonlinear kernels that can be used to learn different
data patterns in machine learning. By encoding the data into either amplitude
or argument of the parameter of squeezing, different kernels can be obtained.
We have applied those quantum kernels with classical support vector machines
for classification tasks on some standard datasets, and their performances are
comparable or better to typical classical kernels. Furthermore, we have
suggested a quantum platform, the trapped ion quantum computer, which can
implement the quantum algorithm involving both discrete-variable qubits and
continuous variables. We mention that the squeezing and displacement operator
suggested in this article belong to a larger family called Gaussian
operations, and general Gaussian-state encoding may be used as quantum kernels
for machine learning. Moreover, since Gaussian-state can be characterized by
the covariance matrix, a fully exploitation of the infinite dimensionality of
continuous variables should refer to non-Gaussian quantum states. The case of
non-Gaussian state encoding and its applications as quantum kernels with
quantum advantages will be investigated in the future.

###### Acknowledgements.

This work was supported by the Key-Area Research and Development Program of
GuangDong Province (Grant No. 2019B030330001) and the Collaborative Research
Fund of Hong Kong (No. C6005-17G).

## References

  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017).
  * Lloyd  _et al._ [2014] S. Lloyd, M. Mohseni, and P. Rebentrost, Quantum principal component analysis, Nat. Phys. 10, 631 (2014).
  * Rebentrost  _et al._ [2014] P. Rebentrost, M. Mohseni, and S. Lloyd, Quantum support vector machine for big data classification, Phys. Rev. Lett. 113, 130503 (2014).
  * Dunjko  _et al._ [2016] V. Dunjko, J. M. Taylor, and H. J. Briegel, Quantum-enhanced machine learning, Phys. Rev. Lett. 117, 130501 (2016).
  * Lloyd  _et al._ [2016] S. Lloyd, S. Garnerone, and P. Zanardi, Quantum algorithms for topological and geometric analysis of data, Nature Communications 7, 10138 (2016).
  * Lloyd and Weedbrook [2018] S. Lloyd and C. Weedbrook, Quantum generative adversarial learning, Phys. Rev. Lett. 121, 040502 (2018).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Schuld and Killoran [2019] M. Schuld and N. Killoran, Quantum machine learning in feature hilbert spaces, Phys. Rev. Lett. 122, 040504 (2019).
  * Huang  _et al._ [2021] H.-Y. R. Huang, M. B. Broughton, M. Mohseni, R. Babbush, S. Boixo, H. Neven, and J. R. McClean, Power of data in quantum machine learning, Nature Communications 12, 1 (2021).
  * Liu  _et al._ [2021] Y. Liu, S. Arunachalam, and K. Temme, A rigorous and robust quantum speed-up in supervised machine learning, Nature Physics 10.1038/s41567-021-01287-z (2021).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103, 150502 (2009).
  * Wiebe  _et al._ [2012] N. Wiebe, D. Braun, and S. Lloyd, Quantum algorithm for data fitting, Phys. Rev. Lett. 109, 050505 (2012).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Phys. Rev. A 98, 032309 (2018).
  * Lloyd  _et al._ [2020] S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran, Quantum embeddings for machine learning (2020), arXiv:2001.03622 [quant-ph] .
  * Pérez-Salinas  _et al._ [2020] A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, and J. I. Latorre, Data re-uploading for a universal quantum classifier, Quantum 4, 226 (2020).
  * Schuld  _et al._ [2021] M. Schuld, R. Sweke, and J. J. Meyer, Effect of data encoding on the expressive power of variational quantum-machine-learning models, Phys. Rev. A 103, 032430 (2021).
  * Schuld [2021] M. Schuld, Supervised quantum machine learning models are kernel methods (2021), arXiv:2101.11020 [quant-ph] .
  * C. M. Bishop [2006] C. M. Bishop,  _Pattern recognition and machine learning_ , Vol. 1 (Springer, 2006).
  * Trevor Hastie  _et al._ [2009] Trevor Hastie, Robert Tibshirani, and Jerome Friedman,  _The Elements of Statistical Learning_, 2nd ed. (Springer, 2009).
  * Lau  _et al._ [2017] H. K. Lau, R. Pooser, G. Siopsis, and C. Weedbrook, Quantum machine learning over infinite dimensions, Phys. Rev. Lett. 118, 080501 (2017).
  * Zhang  _et al._ [2019] D.-B. Zhang, Z.-Y. Xue, S.-L. Zhu, and Z. D. Wang, Realizing quantum linear regression with auxiliary qumodes, Phys. Rev. A 99, 012331 (2019).
  * Killoran  _et al._ [2019] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, Continuous-variable quantum neural networks, Phys. Rev. Research 1, 033063 (2019).
  * Zhang  _et al._ [2020] D.-B. Zhang, S.-L. Zhu, and Z. D. Wang, Protocol for implementing quantum nonparametric learning with trapped ions, Phys. Rev. Lett. 124, 010506 (2020).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys. 84, 621 (2012).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12, 2825 (2011).
  * Leibfried  _et al._ [2003] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75, 281 (2003).
  * Häffner  _et al._ [2008] H. Häffner, C. F. Roos, and R. Blatt, Quantum computing with trapped ions, Physics Reports 469, 155 (2008).
  * Monroe and Kim [2013] C. Monroe and J. Kim, Scaling the ion trap quantum processor, Science 339, 1164 (2013).
  * Ortiz-Gutiérrez  _et al._ [2017] L. Ortiz-Gutiérrez, B. Gabrielly, L. F. Muñoz, K. T. Pereira, J. G. Filgueiras, and A. S. Villar, Continuous variables quantum computation over the vibrational modes of a single trapped ion, Opt. Commun. 397, 166 (2017).
  * Zhang  _et al._ [2018] J. Zhang, M. Um, D. Lv, J.-N. Zhang, L.-M. Duan, and K. Kim, Noon states of nine quantized vibrations in two radial modes of a trapped ion, Phys. Rev. Lett. 121, 160502 (2018).
  * Flühmann  _et al._ [2019] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, Nature 566, 513 (2019).