# Transition System Representation of Boolean Control Networks

Daizhan Cheng, Xiao Zhang and Zhengping Ji This work is supported partly by
NNSF 62073315 of China, and China Postdoctoral Science Foundation 2021M703423
and 2022T150686.D. Cheng is with the Academy of Mathematics and Systems
Science, Chinese Academy of Sciences, Beijing 100190, P.R.China,
<EMAIL_ADDRESS>Zhang is with the National Center for Mathematics and
Interdisciplinary Sciences & the Key Laboratory of Systems and Control,
Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
Beijing 100190, P.R.China<EMAIL_ADDRESS>Ji is with the Key
Laboratory of Systems and Control, Academy of Mathematics and Systems Science
& School of Mathematical Sciences, University of Chinese Academy of Sciences,
Beijing 100190, P.R.China<EMAIL_ADDRESS>
###### Abstract

First, the topological structure of a transition system is studied. Then, two
types of transition system (TS) representations of Boolean networks (BNs) and
Boolean control networks (BCNs) are investigated. The first kind of
representation is state-based, which converts a BCN into a TS with either
distinct control or non-distinct control. The second representation is output-
based, which is also called the simulation of the original BCN. Some
applications are also studied.

##  I Introduction

Since Kauffman proposed the BN to formulate genetic networks [17], the study
of BNs and BCNs becomes a heat topic in the biological community as well as in
the control community. Consider a BN,

$\displaystyle\begin{cases}X_{1}(t+1)=f_{1}(X_{1}(t),X_{2}(t),\cdots,X_{n}(t)),\\\
X_{2}(t+1)=f_{2}(X_{1}(t),X_{2}(t),\cdots,X_{n}(t)),\\\ \vdots,\\\
X_{n}(t+1)=f_{n}(X_{1}(t),X_{2}(t),\cdots,X_{n}(t)),\\\ \end{cases}$ (1)

where $X_{i}\in{\mathcal{D}}:=\\{0,1\\}$,
$f_{i}:{\mathcal{D}}^{n}\rightarrow{\mathcal{D}}$, $i\in[1,n]$.

It is clear that every trajectory can converge to an attractor (either a fixed
point or a cycle) because it consists only of finite nodes and each node can
take only two values ${\mathcal{D}}=\\{0,1\\}$. Thus, the attractors, with
their attractor basins, form the entire topological structure of a BN.
Therefore, finding both fixed points and cycles of a given BN becomes a first
priority problem in the study of BN. Many early works considered this problem
by providing methods to solve a certain class of BNs [8, 13, 9, 12], to name a
few.

In the last two decades, the semi-tensor product (STP) of matrices has been
used to transform BN (or BCN) into an algebraic discrete-time dynamical
(control) system. The STP approach to BN (BCN) has proven to be very
efficient. Many theoretical results have been obtained.

The basic step for the STP approach can be described in a nutshell as follows:
construct the so-called vector form of $X_{i}$ as

$x_{i}:=\begin{bmatrix}X_{i},\\\ 1-X_{i},\end{bmatrix}\in\Delta,\quad
i\in[1,n],$

where $\Delta:=\Delta_{2}$, and $\Delta_{k}:=\operatorname{Col}(I_{k})$ is the
column set of the identity matrix $I_{k}$.

Using the vector form for $X_{i}$, the system (1) can be expressed in its
algebraic state space representation (ASSR) as [3]

$\displaystyle x(t+1)=Mx(t),$ (2)

where $x(t)=\ltimes_{i=1}^{n}x_{i}(t)\in\Delta_{2^{n}}$ and
$M=M_{1}*M_{2}*\cdots,*M_{n}\in{\mathcal{L}}_{2^{n}\times 2^{n}}$, with
$M_{i}$ being the structure matrix of $f_{i}$, $*$ is the Khatri-Rao product
of matrices [5].

The general formula for calculating the number of attractors is given by [2]
as

###### Theorem I.1

[2] Consider the Boolean network (1) with its ASSR (2). Then

$\displaystyle\begin{cases}N_{1}=\operatorname{trace}(M),\\\
N_{s}=\frac{\operatorname{trace}(M^{s})-\mathop{\sum}\limits_{k\in{\mathcal{P}}(s)}kN_{k}}{s},\quad
2\leq s\leq n,\end{cases}$ (3)

where $N_{s}$ is the number of cycles of length $s$, in particular, $N_{1}$ is
the number of fixed points considered as cycles of length $1$. Note that
${\mathcal{P}}(s)$ is the set of the proper factors of $s$, including $1$ and
excluding $s$.

###### Remark I.2

The proof of Theorem I.1 is based on the following three observations:

  * (i)

If $x_{i}\in\Delta_{2^{n}}$ is on a cycle of length $s$, then $x_{i}$ is a
fixed point of $M^{s}$.

  * (ii)

If $x_{i}\in\Delta_{2^{n}}$ is on a cycle of length $t$ and $t|s$, then
$x_{i}$ is also a fixed point of $M^{s}$.

  * (iii)

If $x_{i}\in\Delta_{2^{n}}$ is a fixed point of $M^{s}$, it must be either of
type (i) or of type (ii).

Consider a $k$-valued logical network. Assume it is expressed as in (1) with
only $X_{i}\in{\mathcal{D}}_{k}=\\{1,\cdots,k\\}$, $i\in[1,n]$. Let
$j\sim\delta_{k}^{j}$, $j\in[1,k]$, and

$x_{i}=\delta_{k}^{j},\quad X_{i}=j,\quad i\in[1,n],\;j\in[1,k],$

is the vector form of $X_{i}$. Then its ASSR is still expressed as in (2) with
only $M\in{\mathcal{L}}_{k^{n}\times k^{n}}$.

Taking into account the observations of Remark I.2, the following result is
also obvious.

###### Corollary I.3

[20] Consider the $k$-valued logical network (1) with its ASSR (2). Then the
formula (3) remains true.

The cycles of BCN have also been studied in several papers, e.g. [25, 19] etc.
But there is no formula similar to (3) to calculate all control attractors.

A transition system (TS) is a more general finite-valued network. A BN can be
seen as an autonomous TS and a BCN as a TS. The TS itself is a very useful
framework for representing finite automata (FA) [10, 15, 16]. In particular,
it provides a fundamental framework for hybrid systems [22, 21]. The STP
approach to FAs has also been developed [23, 24, 7].

To our best knowledge, the topological structure of a TS, or the structure of
its attractors, has not been clearly revealed. “Can the formula (3) be applied
to them?” is still an open problem.

In this paper, we first consider the topological structure of autonomous TS.
Two types of cycles are considered: simple cycles and compound cycles. Then
the number of attractors with different lengths is calculated. Then the
transformation of BCN into autonomous TS is considered. Using transformed TS,
the attractors of BCN can be calculated. Finally, some applications of such
transformations are examined.

The remainder of this paper is structured as follows:

Section II considers the topological structure of TSs, and the formula of
fixed points and cycles for networks is extended to TSs. Section III
investigates the state-based representation of BCNs with some direct
applications. Section IV discusses the output-based representation of BCNs,
called the simulation of BCNs. Furthermore, the formula for the simulation
dynamics is obtained. The output robust controls are also investigated.
Section V is a brief conclusion.

The notations used in the text are listed below:

  * •

$\operatorname{Col}(A)$ ($\operatorname{Row}(A)$): the set of columns (rows)
of $A$.

  * •

${\mathcal{D}}_{k}:=\\{1,2,3,\cdots,k\\}$.

  * •

$\delta_{n}^{i}$: the $i$-th column of the identity matrix $I_{n}$.

  * •

$\Delta_{n}:=\Delta_{n}=\\{\delta_{n}^{i}\,|\,i=1,2,\cdots,n\\}$.

  * •

${\mathcal{L}}_{m\times n}$: the set of ${m\times n}$ logical matrices, that
is, $\operatorname{Col}({\mathcal{L}})\subset\Delta_{m}$.

  * •

${\mathcal{B}}_{m\times n}$: the set of $m\times n$ Boolean matrices, that is,
$[{\mathcal{B}}]_{i,j}\in{\mathcal{D}}:=\\{0,1\\}$.

  * •

$\delta_{m}[i_{1},i_{2},\cdots,i_{n}]:=[\delta_{m}^{i_{1}},\delta_{m}^{i_{2}},\cdots,\delta_{m}^{i_{n}}]\in{\mathcal{L}}_{m\times
n}$.

  * •

$A+_{\mathcal{B}}B$: the Boolean addition of $A,B\in{\mathcal{B}}_{m\times n}$
that is, $[A+_{\mathcal{B}}B]_{i,j}=[A]_{i,j}\vee[B]_{i,j}$,
$\mathop{\sum}\limits_{\mathcal{B}}$ is the Boolean sum.

  * •

$A\times_{\mathcal{B}}B$: the Boolean product of $A,B$. For
$A\in{\mathcal{B}}_{n\times n}$,
$A^{(s)}:=\underbrace{A\times_{\mathcal{B}}A\times_{\mathcal{B}}\cdots\times_{\mathcal{B}}A}_{s}$.

##  II Transition Systems

###### Definition II.1

[1] A TS can be described by
$T=({\mathcal{X}},{\mathcal{U}},\Sigma,{\mathcal{O}},h)$, where

  * (i)

${\mathcal{X}}:=\\{X_{1},X_{2},\cdots,X_{n}\\}$ is a finite state set and
$X_{i}\in{\mathcal{D}}_{2}$ for Boolean TS (or $X_{i}\in{\mathcal{D}}_{k}$ for
$k$-valued TS);

  * (ii)

${\mathcal{U}}:=\\{U_{1},U_{2},\cdots,U_{m}\\}$ is a finite input set and
$U_{j}\in{\mathcal{D}}_{2}$ for Boolean TS (or $U_{j}\in{\mathcal{D}}_{k}$ for
$k$-valued TS);

  * (iii)

$\Sigma:{\mathcal{X}}\times{\mathcal{U}}\rightarrow 2^{\mathcal{X}}$ is the
state transition mapping;

  * (iv)

${\mathcal{O}}:=\\{O_{1},O_{2},\cdots,O_{p}\\}$ is the observing set;

  * (v)

$h:{\mathcal{X}}\rightarrow{\mathcal{O}}$ is the observation mapping.

If $|\Sigma(X,U)|\leq 1$, for $\forall X\in{\mathcal{X}}$ and $\forall
U\in{\mathcal{U}}$, then $T$ is called a deterministic TS, otherwise, it is
called a non-deterministic TS.

Then the dynamics of a TS can be expressed as

$\displaystyle\begin{cases}X(t+1)=\Sigma(U(t),X(t)),\\\
Y(t)=h(X(t)).\end{cases}$ (4)

For a TS where $|{\mathcal{X}}|=n$, $|{\mathcal{U}}|=m$, and
$|{\mathcal{O}}|=p$, a subset $X(t)\in 2^{\mathcal{X}}$ can be expressed into
vector form as
$x(t)=(x_{1}(t),x_{2}(t),\cdots,x_{n}(t))^{\mathrm{T}}\in{\mathcal{B}}^{n}$,
which is a Boolean vector, where

$x_{i}(t)=\begin{cases}1,\quad X_{i}\in X(t),\\\ 0,\quad X_{i}\not\in
X(t),\quad i\in[1,n].\end{cases}$

Using vector form expressions, similar to BN (or $k$-valued network), the ASSR
of a transition system can be expressed as

$\displaystyle\begin{cases}x(t+1)=Lu(t)x(t),\\\ y(t)=Hx(t),\\\ \end{cases}$
(5)

where $L\in{\mathcal{B}}_{n\times nm}$ is a Boolean matrix,
$H\in{\mathcal{L}}_{p\times n}$ is a logical matrix.

It is also clear that the ASSR of an autonomous TS is

$\displaystyle\begin{cases}x(t+1)=Mx(t),\\\ y(t)=Hx(t),\\\ \end{cases}$ (6)

where, $M\in{\mathcal{B}}_{n\times n}$ is a Boolean matrix.

###### Example II.2

[1] Consider a TS as $T=({\mathcal{X}},{\mathcal{U}},\Sigma,{\mathcal{O}},h)$
(Ref. to Fig. 1), where

  * (i)

${\mathcal{X}}=\\{x_{1},x_{2},x_{3},x_{4}\\}$.

  * (ii)

${\mathcal{U}}=\\{u_{1},u_{2}\\}$.

  * (iii)

$\begin{array}[]{ll}\Sigma(x_{1},u_{1})=\\{x_{2},x_{3}\\},&\Sigma(x_{2},u_{1})=\\{x_{2},x_{3}\\},\\\
\Sigma(x_{2},u_{2})=\\{x_{4}\\},&\Sigma(x_{3},u_{2})=\\{x_{2},x_{3}\\},\\\
\Sigma(x_{4},u_{1})=\\{x_{2},x_{4}\\}.&{}\\\ \end{array}$

  * (iv)

${\mathcal{O}}=\\{O_{1},O_{2},O_{3}\\}$.

  * (v)

$h(x_{1})=O_{1},\quad h(x_{2})=h(x_{4})=O_{2},\quad h(x_{3})=O_{3}.$

$x_{1}$$x_{3}$$x_{2}$$x_{4}$$u_{1}$$u_{1}$$u_{1}$$u_{1}$$u_{2}$$u_{2}$$u_{2}$$u_{1}$$u_{1}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}O_{1}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}O_{3}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}O_{2}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}O_{2}}$
Figure 1: TS of Example II.2

Let

$\begin{array}[]{ll}x_{i}=\delta_{4}^{i},&i=1,2,3,4;\\\
u_{j}=\delta_{2}^{j},&j=1,2;\\\ y_{k}=\delta_{3}^{k},&k=1,2,3.\end{array}$

Then the ASSR of $T$ is

$\displaystyle\begin{cases}x(t+1)=Lu(t)x(t),\\\ y(t)=Hx(t),\end{cases}$ (7)

where,

$L=\begin{bmatrix}0&0&0&0&0&0&0&0\\\ 1&1&0&1&0&0&1&0\\\ 1&1&0&0&0&0&1&0\\\
0&0&0&1&0&1&0&0\\\ \end{bmatrix},$ $H=\delta_{3}[1,2,3,2].$

######  Definition II.3 ([1])

Consider the transition system (4).

  * (i)

A set $X(t,U,X(0)):=\\{X(0),X(1),\cdots,X(t)\\}$ is called a trajectory
starting from $X(0)$ and driven by $U=\\{U(0),U(1),\cdots,U(t-1)\\}$, where
$X(\tau)\in{\mathcal{X}}$, $\tau\in[0,t]$, and

$X(\tau+1)\in\Sigma(U(\tau),X(\tau)),\quad\tau\in[0,t-1].$

  * (ii)

A trajectory $\\{X(\tau),X(\tau+1),\cdots,X(\tau+\ell-1)\\}$ is called a
general cycle (GC) of length $\ell$ if $X(\tau+\ell)=X(\tau)$.

  * (iii)

A cycle of length $1$ is called a fixed point.

  * (iv)

A cycle is called a simple cycle (SC) if $X(\tau)$, $X(\tau+1)$, $\cdots$,
$X(\tau+\ell-1)$ are distinct (fixed points are also considered SCs).
Otherwise, it is called a compound cycle (CC).

  * (v)

A cycle $C_{P}$ is called a (non-trivial) power cycle (PC), if $C$ is a cycle
and $C_{P}=\underbrace{C~{}C~{}\cdots~{}C}_{k}$, $k\geq 2$. Obviously, a PC
must be a CC.

###### Remark II.4

The cycle of an autonomous TS can be defined as a simplified version of
Definition II.3. We omit it and assume that the cycle of an autonomous TS is
also properly defined.

###### Example II.5

Consider an autonomous TS, denoted by $T$, which has a transition graph as
shown in Fig. 2.

$x_{1}$$x_{2}$ Figure 2: Transition System $T$

Its ASSR is easily obtained as

$\displaystyle x(t+1)=\begin{bmatrix}1&1\\\ 1&0\end{bmatrix}x(t).$ (8)

Then

  * (i)

$1$ is a fixed point.

  * (ii)

$(1,2)$ is an SC.

  * (iii)

$(1,2,1,2,1,2)$ is a PC.

  * (iv)

$(1,2,1,1,2,1,1,1,2,\cdots,\underbrace{1,\cdots,1}_{s},2),\quad s=1,2,\cdots$
are CCs of length $\frac{s(s+3)}{2}$, which can be arbitrarily long.

Consider a partition as:

$\mbox{GC}~{}=~{}\mbox{SC}~{}\bigcup~{}\mbox{CC}.$

Following the argument in Remark I.2, it is easy to obtain the following
result:

###### Proposition II.6

Consider an autonomous TS with its ASSR (6). Then formula (3) is applicable to
calculate its numbers of CCs of arbitrary length $s>0$.

As shown in Example II.5, the length $s>0$ of a CC could be arbitrarily large,
so applying formula (3) to calculate the numbers of all its CCs with arbitrary
length $s>0$ is impossible. The problem can be solved by observing the
following facts:

  * •

An SC is a cycle of length $s\leqslant n$;

  * •

Any CC can be viewed as a recursive concatenation or insertion of SCs. For
example, consider a trajectory of the form
$C:=(y,\cdots,\xi,x_{1},x_{2},\cdots,x_{\ell-1},\xi,\cdots,y)$ where
$x_{1},x_{2},\cdots,x_{\ell-1}$ are distinct; then
$(\xi,x_{1},x_{2},\cdots,x_{\ell-1})$ is an SC (of length less than $n$);
replace the subsequence $\xi,x_{1},x_{2},\cdots,x_{\ell-1},\xi$ in $C$ by
$\xi$ and keep finding subcycles that are SCs in the remaining trajectory, we
will end up partitioning the trajectory into SCs and finding it constructed by
recursively inserting SCs into an SC.

Therefore, we only have to compute all SCs using the formula (3) and the
algorithm for finding attractors of BNs (or $k$-valued networks) [2]. Then the
SCs are sufficient to describe the topological structure of the autonomous TS.

We consider a simple example.

###### Example II.7

Consider a TS, its ASSR is

$\displaystyle x(t+1)=\begin{bmatrix}0&0&1&0\\\ 1&0&0&0\\\ 1&1&0&0\\\
0&0&1&1\\\ \end{bmatrix}x(t):=Mx(t).$ (9)

Using the formula (3), it is easy to calculate that

$N_{1}=1,\quad N_{2}=1,\quad N_{3}=1,\quad N_{4}=0.$

The corresponding attractors are:

$(4),\quad(1,3),\quad(1,2,3).$

They are all SCs.

##  III TS Representation of BCNs

###  III-A Conversion of BCNs to Autonomous TSs

Consider a BCN as

$\displaystyle\begin{cases}X_{1}(t+1)=f_{1}(X_{1}(t),\cdots,X_{n}(t),U_{1}(t),\cdots,U_{m}(t)),\\\
X_{2}(t+1)=f_{2}(X_{1}(t),\cdots,X_{n}(t),U_{1}(t),\cdots,U_{m}(t)),\\\
\vdots\\\
X_{n}(t+1)=f_{n}(X_{1}(t),\cdots,X_{n}(t),U_{1}(t),\cdots,U_{m}(t)),\\\
\end{cases}$ (10)

where $X_{i},U_{j}\in{\mathcal{D}}$,
$f_{i}:{\mathcal{D}}^{n+m}\rightarrow{\mathcal{D}}$, $i\in[1,n]$, $j\in[1,m]$.

The ASSR of (10) is

$\displaystyle x(t+1)=Lu(t)x(t),$ (11)

where $x(t)=\ltimes_{i=1}^{n}x_{i}(t)\in\Delta_{2^{n}}$,
$u(t)=\ltimes_{j=1}^{m}u_{j}(t)\in\Delta_{2^{m}}$,
$L\in{\mathcal{L}}_{2^{n}\times 2^{m+n}}$.

###### Definition III.1

Consider system (10).

  * (i)

If there exists a sequence

$u(\tau),u(\tau+1),\cdots,u(\tau+s-1)$

such that

$\begin{array}[]{l}x(\tau)\xrightarrow{u(\tau)}x(\tau+1)\xrightarrow{u(\tau+1)}x(\tau+2)\cdots
x(\tau+s-1)\\\ \xrightarrow{u(\tau+s-1)}x(\tau+s).\end{array}$

Then $x(\tau),x(\tau+1),\cdots,x(\tau+s)$ is called a control trajectory with
undistinguished control of length $s$.

The sequence of state-control pairs
$(x(\tau),u(\tau)),(x(\tau+1),u(\tau+1)),\cdots,(x(\tau+s),u(\tau+s))$ is
called a control trajectory with distinguished control of length $s$.

  * (ii)

A control trajectory $x(\tau),x(\tau+1),\cdots,x(\tau+s)$ is called a control
cycle of length $s$, if $x(\tau)=x(\tau+s)$. The simple (or power, compound,
etc.) control cycle can be defined similarly.

  * (iii)

A control cycle of length $1$ is called a control fixed point.

###### Definition III.2

Consider BCN (10).

  * (i)

It is converted into an autonomous TS with undistinguished control, where the
ASSR of the TS is

$\displaystyle x(t+1)=Mx(t),$ (12)

where

$\displaystyle
M={\mathop{\sum}\limits_{{\mathcal{B}}}}_{i=1}^{2^{m}}(L\delta_{2^{m}}^{i}),$
(13)

here $\mathop{\sum}\limits_{{\mathcal{B}}}$ is the Boolean sum.

  * (ii)

It is converted into an autonomous TS with distinguished control, where the
ASSR of the TS is

$\displaystyle w(t+1)=\Xi w(t),$ (14)

where $w(t)=u(t)x(t)$,
$\Xi={\underbrace{[L^{\mathrm{T}}~{}L^{\mathrm{T}}\cdots~{}L^{\mathrm{T}}]}_{2^{m}}}^{\mathrm{T}}.$

Using Proposition II.6 to the converted autonomous TS yields the following
result, which can be used to calculate control cycles of BCNs.

###### Corollary III.3

  * (i)

Applying the formula (3) to the converted autonomous TS with undistinguished
control (12), the number of cycles of a BCN with undistinguished control of
different lengths $s$ can be calculated.

  * (ii)

Applying the formula (3) to the converted autonomous TS with distinguished
control (13), the number of cycles of a BCN with distinguished control of
different lengths $s$ can be calculated.

###### Remark III.4

With some mild revision, the above results can be naturally extended to
$k$-valued control networks. Similarly, they are also applicable to a general
TS, when it is converted into an autonomous TS. The following example shows
this.

###### Example III.5

Recall Example II.2.

  * (i)

A straightforward computation shows that its converted autonomous TS with
undistinguished control is determined by its ASSR as

$x(t+1)=M_{I}x(t),$

where the transition matrix is

$\displaystyle M_{I}=\begin{bmatrix}0&0&0&0\\\ 1&1&1&1\\\ 1&1&1&0\\\
0&1&0&1\\\ \end{bmatrix}.$ (15)

  * (ii)

The ASSR of its converted autonomous TS with distinguished control is

$z(t+1)=\Xi z(t),$

where

$\displaystyle\Xi=\begin{bmatrix}0&0&0&0&0&0&0&0\\\ 1&1&0&1&0&0&1&0\\\
1&1&0&0&0&0&1&0\\\ 0&0&0&1&0&1&0&0\\\ 0&0&0&0&0&0&0&0\\\ 1&1&0&1&0&0&1&0\\\
1&1&0&0&0&0&1&0\\\ 0&0&0&1&0&1&0&0\\\ \end{bmatrix}.$ (16)

  * (iii)

Using the formula (3) the number of CCs with undistinguished control are

    * (a)

$N_{1}=3$:

$X_{\mbox{fixed point}}=\\{2,3,4\\},$

where, $i$ stands for $x_{i}$, $i=1,2,3,4$.

    * (b)

$N_{2}=2$:

$C^{2}_{1}=(2,3);\quad C^{2}_{2}=(2,4).$

    * (c)

$N_{3}=4$:

$\begin{array}[]{ll}C^{3}_{1}=(2,2,3);&C^{3}_{2}=(2,3,3);\\\
C^{3}_{3}=(2,2,4);&C^{3}_{4}=(2,4,4).\end{array}$

    * (d)

$N_{4}=7$:

$\begin{array}[]{ll}C^{4}_{1}=(2,2,2,3);&C^{4}_{2}=(2,2,3,3);\\\
C^{4}_{3}=(2,3,3,3);&C^{4}_{4}=(2,2,2,4);\\\
C^{4}_{5}=(2,2,4,4);&C^{4}_{6}=(2,4,4,4);\\\ C^{4}_{7}=(2,3,2,4).&{}\\\
\end{array}$

    * (e)

$N_{5}=16$:

$\begin{array}[]{ll}C^{5}_{1}=(2,2,2,2,3);&C^{5}_{2}=(2,2,2,3,3);\\\
C^{5}_{3}=(2,2,3,3,3);&C^{5}_{4}=(2,3,3,3,3);\\\
C^{5}_{5}=(2,3,2,2,3);&C^{5}_{6}=(2,3,2,3,3);\\\
C^{5}_{7}=(2,2,2,2,4);&C^{5}_{8}=(2,2,2,4,4);\\\
C^{5}_{9}=(2,2,4,4,4);&C^{5}_{10}=(2,4,4,4,4);\\\
C^{5}_{11}=(2,4,2,2,4);&C^{5}_{12}=(2,4,2,4,4);\\\
C^{5}_{13}=(2,3,2,2,4);&C^{5}_{14}=(2,3,2,4,2);\\\
C^{5}_{15}=(2,3,3,2,4);&C^{5}_{16}=(2,4,4,2,3).\\\ \end{array}$ $\vdots$

Finally, the SC is

$SC:\\{\\{2\\},\\{3\\},\\{4\\},\\{2,3\\},\\{2,4\\}\\}.$

###  III-B Some Applications

  * •

Reachability of TSs:

Consider an autonomous TS (in ASSR form):

$\displaystyle x(t+1)=Mx(t),$ (17)

where $x(t)\in\Delta_{n}$, $M\in{\mathcal{B}}_{n\times n}$.

###### Definition III.6

$x_{j}$ is reachable from $x_{i}$, if the trajectory $x(t,x_{0})$, starting
from $x_{0}=x_{i}$, can reach $x_{j}$ at finite time $t_{0}$, i.e.,
$x(t_{0},x_{0})=x_{j}$.

Define the reachable matrix ${\mathcal{C}}$ as

$\displaystyle{\mathcal{C}}={\mathop{\sum}\limits_{{\mathcal{B}}}}_{s=1}^{n}M^{(s)}.$
(18)

Then the following result is well known.

######  Proposition III.7 ([3])

Assume (17) is the converted autonomous TS of a BCN $\Sigma$. Then for (17)
$x_{j}$ is reachable from $x_{i}$, if and only if, for BCN $\Sigma$: $x_{j}$
is reachable from $x_{i}$.

  * •

Decoupling:

###### Definition III.8

    * (i)

A subset $Z\subset\Delta_{n}$ is called an attractor of (17), if $x(t)\in Z$
implies $x(t+1)\in Z$ for $\forall t\in{\mathbb{N}}$.

    * (ii)

Assume $Z=\\{x_{0}\\}$ is an attractor, then $Z$ (or $x_{0}$) is called a
fixed point.

The following result is obvious:

###### Proposition III.9

    * (i)

Suppose (after a coordinate change if necessary)
$Z=\\{x_{1},x_{2},\cdots,x_{r}\\}$. Then $Z$ is an attractor, if and only if
$M$ has the following block upper triangle form:

$\displaystyle M=\begin{bmatrix}M_{1}&M_{2}\\\ 0&M_{3}\\\ \end{bmatrix}$ (19)

where $M_{1}\in{\mathcal{B}}_{r\times r}$.

    * (ii)

Suppose (after a possible coordinate change)
$Z_{i}=\\{x^{i}_{1},x^{i}_{2},\cdots,x^{i}_{r_{i}}\\}$, $i=1,2,\cdots,s$ are
sets of disjoint states. Then $Z_{i}$, $i=1,2,\cdots,s$ are attractor sets if
and only if, $M$ has the form (18), where $M_{1}$ is a block diagonal matrix
such that

$M_{1}=\begin{bmatrix}M_{1}^{1}&0&\cdots&0\\\ 0&M_{1}^{2}&\cdots&0\\\
\vdots&~{}&~{}&~{}\\\ 0&0&\cdots&M_{1}^{s}\\\ \end{bmatrix},$

where $M_{1}^{i}\in{\mathcal{B}}_{r_{i}\times r_{i}}$, $i\in[1,s]$.

##  IV Simulation of BNs and BCNs

###  IV-A Output-Based Simulation

Consider a logical control network or a deterministic TS (in ASSR form),
denoted by $\Sigma$, and defined by

$\displaystyle\begin{cases}x(t+1)=Mu(t)x(t),\\\ y(t)=Hx(t),\end{cases}$ (20)

where $x(t)\in\Delta_{n}$, $u(t)\in\Delta_{m}$, $M\in{\mathcal{L}}_{n\times
mn}$, $H\in{\mathcal{L}}_{p\times n}$.

The following definition is based on [1] with mild formation modification.

###### Definition IV.1

Consider control network $\Sigma$ (20).

  * (i)

Two states $x_{i}$ and $x_{j}$ are said to be (output) equivalent, denoted by
$x_{i}\sim x_{j}$, if $Hx_{i}=Hx_{j}$.

  * (ii)

The quotient system, denoted by $\Sigma/\sim$ is called a simulation of
$\Sigma$.

Denote by $\bar{x}$ the equivalence class of $x$;
$con(\bar{x}):=\\{x\;|\;x\sim\bar{x}\\}$; ${\mathcal{O}}(x)$ the output
trajectory for the state trajectory starting from $x$. Then we have

###### Proposition IV.2

[1]

$\displaystyle{\mathcal{O}}_{\Sigma}(con(\bar{x}))\subset{\mathcal{O}}_{\Sigma/\sim}(\bar{x}).$
(21)

From the set controllability point of view [11, 6], the following is a
straightforward result [14].

###### Proposition IV.3

The simulation $\Sigma/\sim$ of $\Sigma$ is a transition system, and its
dynamics is

$\displaystyle\begin{cases}\bar{x}(t+1)=\left[H\times_{{\mathcal{B}}}M\times_{{\mathcal{B}}}(I_{m}\otimes
H^{T})\right]u(t)\bar{x}(t),\\\ y(t)=HH^{T}\bar{x}(t).\end{cases}$ (22)

###  IV-B Output-Robust Network and Control

Consider a network (i.e., a deterministic TS):

$\displaystyle\begin{cases}x(t+1)=Mx(t),\\\ y(t)=Hx(t),\end{cases}$ (23)

where $M\in{\mathcal{L}}_{n\times n}$, $H\in{\mathcal{L}}_{p\times n}$. Assume
(23) is a system with possible disturbances, described by

$\displaystyle M=\begin{cases}M_{0},\quad\xi\in\emptyset,\\\ L\xi,\quad
L\in{\mathcal{L}}_{n\times ns},~{}\xi\in\Delta_{s}.\end{cases}$ (24)

Then we have two models: when $M=M_{0}$ it is called the nominated model; and
when $M=L\xi$ it is called the disturbed model. For the nominated model,
denoted by $\Sigma_{0}$, we can construct its simulation system
$\Sigma_{0}/\sim$. For the disturbed model, we can first get its TSR, denoted
by $\Sigma_{\xi}$, and then construct its simulation system
$\Sigma_{\xi}/\sim$.

###### Definition IV.4

System (23) is said to be output robust, if $\Sigma_{0}/\sim$ and
$\Sigma_{\xi}/\sim$ are bi-simulated, that is to say, they generate identical
output dynamics calculated from (22).

###### Example IV.5

Consider a Boolean network $\Sigma$, which has its nominated model as

$\displaystyle\begin{array}[]{l}\begin{cases}x_{1}(t+1)=\neg x_{1}(t),\\\
x_{2}(t+1)=x_{1}(t)\bar{\vee}x_{3}(t),\\\
x_{3}(t+1)=\left[x_{1}(t)\bar{\vee}x_{2}(t)\right]\vee x_{3}(t),\end{cases}\\\
~{}~{}y(t)=\left[x_{1}(t)\leftrightarrow x_{2}(t)\right]\leftrightarrow\neg
x_{3}(t);\end{array}$ (27)

and its disturbed model as

$\displaystyle\begin{array}[]{l}\begin{cases}x_{1}(t+1)=(\neg\xi(t))\wedge
x_{1}(t),\\\ x_{2}(t+1)=\left[\xi(t)\vee\neg
x_{1}(t)\right]\bar{\vee}x_{3}(t),\\\
x_{3}(t+1)=\left[x_{1}(t)\bar{\vee}x_{2}(t)\right]\vee x_{3}(t),\end{cases}\\\
~{}~{}y(t)=\left[x_{1}(t)\leftrightarrow x_{2}(t)\right]\leftrightarrow\neg
x_{3}(t).\end{array}$ (30)

It is easy to have its ASSR as in (23), where

$\begin{array}[]{l}M_{0}=\delta_{8}[7,6,7,5,1,3,1,4],\\\
L=\delta_{8}[7,6,7,5,7,6,1,4,1,3,7,5,7,6],\\\
H=\delta_{2}[2,1,1,2,1,2,2,1].\end{array}$

The transition representation of its disturbed model becomes $x(t+1)=Tx(t)$,
where

$T=\begin{bmatrix}1&0&1&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&1&0&0&0&0&0&0\\\ 0&0&0&1&0&1&0&0\\\ 0&1&0&0&0&0&0&1\\\ 1&0&1&0&1&0&1&0\\\
0&0&0&0&0&0&0&0\\\ \end{bmatrix}.$

Using Proposition IV.3, it is easy to calculate that the (output-based)
simulations of these two models are the same as

$\bar{x}(t+1)=\begin{bmatrix}0&1\\\ 1&1\\\ \end{bmatrix}\bar{x}(t).$

Hence, system $\Sigma$ is output robust.

Consider a control network (20) where $u(t)\in\Delta_{\ell}$,
$M\in{\mathcal{L}}_{n\times n\ell}$ satisfying (24). Assume there exists a
state feedback control

$\displaystyle u(t)=Gx(t),$ (31)

where $G\in{\mathcal{L}}_{\ell\times n}$, such that the closed-loop system is
output robust, then $u(t)$ is called an output robust control, which solves
the output robust control problem.

###### Example IV.6

Consider a Boolean control network $\Sigma$ with nominated model as

$\displaystyle\begin{array}[]{l}\begin{cases}x_{1}(t+1)=(\neg
x_{1}(t))\vee(\neg u(t)),\\\ x_{2}(t+1)=x_{1}(t)\bar{\vee}x_{3}(t),\\\
x_{3}(t+1)=\left[x_{1}(t)\bar{\vee}x_{2}(t)\right]\vee x_{3}(t),\end{cases}\\\
~{}~{}y(t)=\left[x_{1}(t)\leftrightarrow x_{2}(t)\right]\leftrightarrow\neg
x_{3}(t);\end{array}$ (34)

and its disturbed model as

$\displaystyle\begin{array}[]{l}\begin{cases}x_{1}(t+1)=(\neg\xi(t))\wedge
x_{1}(t)\wedge u(t),\\\ x_{2}(t+1)=\left[\xi(t)\vee\neg
x_{1}(t)\right]\bar{\vee}x_{3}(t),\\\
x_{3}(t+1)=\left[x_{1}(t)\bar{\vee}x_{2}(t)\right]\vee x_{3}(t),\end{cases}\\\
~{}~{}y(t)=\left[x_{1}(t)\leftrightarrow x_{2}(t)\right]\leftrightarrow\neg
x_{3}(t).\end{array}$ (37)

It is easy to see that if we choose

$\displaystyle u(t)=x_{1}(t),$ (38)

then the closed-loop system becomes (27)-(30), which is output robust. Thus,
the state-feedback control (38) solves the output robust problem of $\Sigma$.

Output robust control solves the disturbance decoupling problem without the
regularity assumption [4].

##  V Conclusion

This paper investigated the transition representation of BCNs. The main
contribution consists of three parts: (i) The topology of TSs was considered,
and the formula for calculating fixed points and cycles of BNs was extended to
TSs. (ii) Two types of state-based TS representations of BCNs, namely
representations with either distinct or non-distinct controls, were proposed.
(iii) An output-based TS representation, also called simulation, was studied.
Its dynamic equation was obtained. The output robust (control) was also
studied. The technique proposed in this paper is applicable to any finite
value network. In fact, some examples in this paper are not BN or BCN. But the
technique used for them is exactly the same as the one for BN or BCN.

Some related problems, such as finding the output robust controls, etc., are
left for further study.

## References

  * [1] C. Belta, B. Yordanov, E.A. Gol, Formal Methods for Discrete-Time Dynamical Systems, Springer, Switzerland, 2017.
  * [2] D. Cheng, H. Qi, A linear representation of dynamics of Boolean networks, IEEE Trans. Aut. Contr., Vol. 55, No. 10, 2251-2258, 2010.
  * [3] D. Cheng, H. Qi, Z. Li, Analysis and Control of Boolean Networks - A Semi-tensor Product Approach, Springer, London, 2011.
  * [4] D. Cheng, Disturbance decoupling of Boolean control networks, IEEE Trans. Aut. Contr., Vol. 56, No. 1, 2-10, 2011.
  * [5] D. Cheng, H. Qi, Y. Zhao, An Introduction to Semi-tensor Product of Matrices and Its Applications, World Scientific, Singapore, 2012.
  * [6] D. Cheng, C. Li, F. He, Observability of Boolean networks via set controllability approach, Sys. Contr. Lett., Vol. 155, 22-25, 2018.
  * [7] H. Deng, Y. Yan, Z. Chen, A matrix-based static approach to analysis of finite state machines, Front. Inform. Tech. Electr. Eng., 2022, Vol. 23, No. 8, 1239–1246, 2022.
  * [8] C. Farrow, J. Heidel, J. Rogers, Scalar equations for synchronous Boolean networks with biological applications, IEEE Trans. Neural Netw., Vol. 15, No. 2, 348-354, 2004.
  * [9] B. Goodwin, Temporal Organization in Cells, Academic Press, San Diego, 1963.
  * [10] P. Guan, Cellular automata public key cryptosystem, Complex Sys., Vol. 1, 51-57, 1987.
  * [11] Y. Guo, P. Wang, W. Gui, C. Yan, Set stability and set stabilization of Boolean control networks based on invariant subsets, Automatica, Vol. 61, 106-112, 2015.
  * [12] J. Heidel, J. Farrow, C. Farrow, J. Rogers, Finding cycles in synchronous Boolean networks with applications to biochenical systems, Int. J. Bifurc. Chaos, Vol. 13, No. 3, 535-552, 2003.
  * [13] S. Huang, D. Ingber, Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks. Exp. Cell Res., Vol. 261, No. 1, 91-103, 2000.
  * [14] Z. Ji, X. Zhang, D.Cheng, Aggregated (Bi-)Simulation of Finite Valued Networks, arxiv:2303.14390, 2023.
  * [15] J. Kari, Theory of cellular automata: a survey, Theoret. Comput. Sci., Vol. 334, 3-33, 2005.
  * [16] J. Kari, Automata and Formal Language, Lecture Notes, University of Turku, 2013.
  * [17] S. Kauffman, Metabolic stability and epigenesis in randomly constructed genetic nets, J. Theor. Biol., Vol. 22, No. 3, 437-467, 1969.
  * [18] K.H. Kim, Boolean Matrix Theory and Applications, Marcel Dekker, INC., New York, 1982.
  * [19] D. Laschov, M. Margaliot, Minimum-time control of Boolean networks, SIAM J. Contr. Opt., Vol. 51, No. 4, 2869-2892, 2013.
  * [20] Z. Li, D. Cheng, Algebraic approach to dynamics of multi-valued networks, Int. J. Bif. Chaos, Vol. 20, No. 3, 561-582, 2010.
  * [21] H. Lin, P.J. Antsaklis, Hybrid dynamical systems: an introduction to control and verification, Foundations and trends in Systems and Control, Vol. 1, No. 1, 1-172, 2014.
  * [22] P. Tabuada, Verification and Control of Hybrid Systems - A Symbolic Approach, Springer, New York, 2019.
  * [23] X. Xu, Y. Hong, Observability and observer design for finite automata via matrix approach, IET Control Theory Appl., 2013, Vol. 7, No. 12, 1609–1615, 2013.
  * [24] X. Xu, Y. Hong, Matrix expression to model matching of asynchronous sequential machines, IEEE Trans. Aut. Contr., Vol. 58, No. 11, 2974–2979, 2013.
  * [25] Y. Zhao, H. Qi, D. Cheng, Input-state incidence matrix of Boolean control networks and its applications, Sys. Contr. Lett., Vol. 59, No. 12, 767-774, 2010.


* (34) X. Du, M. P. Heller, S. Schlichting and V. Svensson, _Exponential approach to the hydrodynamic attractor in Yang-Mills kinetic theory_ , _Phys. Rev. D_ 106 (2022) 014016, [2203.16549].
  * (35) D. Almaalol, K. Boguslavski, A. Kurkela and M. Strickland, _Non-equilibrium attractor in high-temperature QCD plasmas_ , 2208.00513.
  * (36) Z. Chen, D. Teaney and L. Yan, _Hydrodynamic attractor of noisy plasmas_ , 2206.12778.
  * (37) H. Alalawi and M. Strickland, _Far-from-equilibrium attractors for massive kinetic theory in the relaxation time approximation_ , 2210.00658.
  * (38) A. Soloviev, _Hydrodynamic attractors in heavy ion collisions: a review_ , _Eur. Phys. J. C_ 82 (2022) 319, [2109.15081].
  * (39) J. Berges, M. P. Heller, A. Mazeliauskas and R. Venugopalan, _QCD thermalization: Ab initio approaches and interdisciplinary connections_ , _Rev. Mod. Phys._ 93 (2021) 035003, [2005.12299].
  * (40) W. Busza, K. Rajagopal and W. van der Schee, _Heavy Ion Collisions: The Big Picture, and the Big Questions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 339–376, [1802.04801].
  * (41) F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, _The Color Glass Condensate_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 463–489, [1002.0333].
  * (42) P. B. Arnold, G. D. Moore and L. G. Yaffe, _Effective kinetic theory for high temperature gauge theories_ , _JHEP_ 01 (2003) 030, [hep-ph/0209353].
  * (43) A. Kurkela and Y. Zhu, _Isotropization and hydrodynamization in weakly coupled heavy-ion collisions_ , _Phys. Rev. Lett._ 115 (2015) 182301, [1506.06647].
  * (44) T. Faulkner and J. Polchinski, _Semi-Holographic Fermi Liquids_ , _JHEP_ 06 (2011) 012, [1001.5049].
  * (45) A. Mukhopadhyay and G. Policastro, _Phenomenological Characterization of Semiholographic Non-Fermi Liquids_ , _Phys. Rev. Lett._ 111 (2013) 221602, [1306.3941].
  * (46) E. Iancu and A. Mukhopadhyay, _A semi-holographic model for heavy-ion collisions_ , _JHEP_ 06 (2015) 003, [1410.6448].
  * (47) A. Mukhopadhyay, F. Preis, A. Rebhan and S. A. Stricker, _Semi-Holography for Heavy Ion Collisions: Self-Consistency and First Numerical Tests_ , _JHEP_ 05 (2016) 141, [1512.06445].
  * (48) S. Banerjee, N. Gaddam and A. Mukhopadhyay, _Illustrated study of the semiholographic nonperturbative framework_ , _Phys. Rev. D_ 95 (2017) 066017, [1701.01229].
  * (49) A. Kurkela, A. Mukhopadhyay, F. Preis, A. Rebhan and A. Soloviev, _Hybrid Fluid Models from Mutual Effective Metric Couplings_ , _JHEP_ 08 (2018) 054, [1805.05213].
  * (50) C. Ecker, A. Mukhopadhyay, F. Preis, A. Rebhan and A. Soloviev, _Time evolution of a toy semiholographic glasma_ , _JHEP_ 08 (2018) 074, [1806.01850].
  * (51) I. Müller, _Zum Paradoxon der Wärmeleitungstheorie_ , _Z. Phys._ 198 (1967) 329–344.
  * (52) W. Israel and J. Stewart, _Transient relativistic thermodynamics and kinetic theory_ , _Annals Phys._ 118 (1979) 341–372.
  * (53) T. Mitra, S. Mondkar, A. Mukhopadhyay, A. Rebhan and A. Soloviev, _Hydrodynamic attractor of a hybrid viscous fluid in Bjorken flow_ , _Phys. Rev. Res._ 2 (2020) 043320, [2006.09383].
  * (54) R. Baier, A. H. Mueller, D. Schiff and D. T. Son, _’Bottom up’ thermalization in heavy ion collisions_ , _Phys. Lett. B_ 502 (2001) 51–58, [hep-ph/0009237].
  * (55) C. Loizides, _Experimental overview on small collision systems at the LHC_ , _Nucl. Phys. A_ 956 (2016) 200–207, [1602.09138].
  * (56) S. Schlichting and P. Tribedy, _Collectivity in Small Collision Systems: An Initial-State Perspective_ , _Adv. High Energy Phys._ 2016 (2016) 8460349, [1611.00329].
  * (57) H. Kowalski, T. Lappi and R. Venugopalan, _Nuclear enhancement of universal dynamics of high parton densities_ , _Phys. Rev. Lett._ 100 (2008) 022303, [0705.3047].
  * (58) S. Mondkar, A. Mukhopadhyay, A. Rebhan and A. Soloviev, _Quasinormal modes of a semi-holographic black brane and thermalization_ , _JHEP_ 11 (2021) 080, [2108.02788].
  * (59) J. Noronha and G. S. Denicol, _Transient Fluid Dynamics of the Quark-Gluon Plasma According to AdS/CFT_ , 1104.2415.
  * (60) M. P. Heller, R. A. Janik, M. Spaliński and P. Witaszczyk, _Coupling hydrodynamics to nonequilibrium degrees of freedom in strongly interacting quark-gluon plasma_ , _Phys. Rev. Lett._ 113 (2014) 261601, [1409.5087].
  * (61) P. Romatschke and U. Romatschke, _Relativistic Fluid Dynamics In and Out of Equilibrium_.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 5, 2019, 10.1017/9781108651998.
  * (62) I. Aniceto, R. Schiappa and M. Vonk, _The Resurgence of Instantons in String Theory_ , _Commun. Num. Theor. Phys._ 6 (2012) 339–496, [1106.5922].
  * (63) M. Coquet, X. Du, J.-Y. Ollitrault, S. Schlichting and M. Winn, _Intermediate mass dileptons as pre-equilibrium probes in heavy ion collisions_ , _Phys. Lett. B_ 821 (2021) 136626, [2104.07622].
  * (64) E. Grossi, A. Soloviev, D. Teaney and F. Yan, _Transport and hydrodynamics in the chiral limit_ , _Phys. Rev. D_ 102 (2020) 014042, [2005.02885].
  * (65) E. Grossi, A. Soloviev, D. Teaney and F. Yan, _Soft pions and transport near the chiral critical point_ , _Phys. Rev. D_ 104 (2021) 034025, [2101.10847].
  * (66) A. Florio, E. Grossi, A. Soloviev and D. Teaney, _Dynamics of the $O(4)$ critical point in QCD_, _Phys. Rev. D_ 105 (2022) 054512, [2111.03640].
  * (67) P. Kovtun, _Lectures on hydrodynamic fluctuations in relativistic theories_ , _J. Phys. A_ 45 (2012) 473001, [1205.5040].
  * (68) T. Dore, J. Noronha-Hostler and E. McLaughlin, _Far-from-equilibrium search for the QCD critical point_ , _Phys. Rev. D_ 102 (2020) 074017, [2007.15083].
  * (69) T. Mitra, A. Mukhopadhyay and A. Soloviev, _Hydrodynamic attractor and novel fixed points in superfluid Bjorken flow_ , _Phys. Rev. D_ 103 (2021) 076014, [2012.15644].
  * (70) R. Rodgers and J. G. Subils, _Boost-invariant superfluid flows_ , _JHEP_ 09 (2022) 205, [2207.02903].


# Towards Unbiased Random Features with Lower Variance For Stationary
Indefinite Kernels
††thanks: ⋆Correspondence author: Xiaolin Huang

Qin Luo, Kun Fang, Jie Yang, Xiaolin Huang⋆ Institute of Image Processing and
Pattern Recognition
Shanghai Jiao Tong University
Shanghai, China
<EMAIL_ADDRESS>
###### Abstract

Random Fourier Features (RFF) demonstrate well-appreciated performance in
kernel approximation for large-scale situations but restrict kernels to be
stationary and positive definite. And for non-stationary kernels, the
corresponding RFF could be converted to that for stationary indefinite kernels
when the inputs are restricted to the unit sphere. Numerous methods provide
accessible ways to approximate stationary but indefinite kernels. However,
they are either biased or possess large variance. In this article, we propose
the generalized orthogonal random features, an unbiased estimation with lower
variance. Experimental results on various datasets and kernels verify that our
algorithm achieves lower variance and approximation error compared with the
existing kernel approximation methods. With better approximation to the
originally selected kernels, improved classification accuracy and regression
ability is obtained with our approximation algorithm in the framework of
support vector machine and regression.

##  I Introduction

Kernel methods are extensively utilized in statistical machine learning with
numerous application, such as classification [1, 2], regression [3, 4, 5] and
dimensionality reduction [6, 7]. One prominent problem for kernel methods is
that it scales poorly with the size of datasets because of higher requirement
in computation and storage space. Given $N$ observations, the computational
complexity is $\mathcal{O}(N^{3})$ and the storage complexity is
$\mathcal{O}(N^{2})$. Random Fourier Features (RFF) [8] reduce the computation
cost and storage space to $\mathcal{O}(Ns^{2})$ and $\mathcal{O}(Ns)$ through
unbiased approximation, where $s$ is the number of Random Fourier Features and
much smaller than $N$. However, there exists relatively large variance for
initial random features, which is an obstacle for better approximation.
Different methods are considered to achieve more accurate and faster
approximation for random features, including Fastfood via Hadamard matrices
[9], Quasic Monte Carlo Sampling with a low discrepancy sequence [10] and the
usage of circulant matrices [11]. Orthogonality is the most effective property
to reduce the variance in high-dimensional approximation. Felix et al. [12]
firstly propose the orthogonality for random features and study the case for
Gaussian kernels. Choromanski et al. [13, 14] extend orthogonality theorem to
arbitrary positive definite (PD) radial basis kernel. Structural orthogonal
random features [12, 15] are further proposed to improve the efficiency for
orthogonality.

Random features demonstrate superior performance in kernel approximation,
however, two important characteristics are prerequisite for having random
features, 1) shift-invariant (also called stationary),
$k(\boldsymbol{x},\boldsymbol{y})=k(\boldsymbol{x}-\boldsymbol{y})$, 2)
positive definite, ensuring the Fourier transform of kernel as probability
measure. However, many kernels might not satisfy the stationary or PD
requirements. Typical examples include polynomial kernel[16], neural tangent
kernel (NTK) [17], a linear combination of Gaussian kernels (Delta-Gaussian
kernel) [18] and the TL1 kernel [19]. For non-stationary kernels like
polynomial kernel, they could be converted to stationary but indefinite
kernels when the inputs are restricted to the unit sphere [20]. Therefore, the
random features for non-stationary kernels or indefinite kernels could be
unified with that for stationary indefinite kernels.

Because of the wide application of non-stationary or indefinite kernels, a
series of works concentrate on random features for these kernels. For non-
stationary kernels, Random Maclaurin [21] and Tensor Sketch [22] are initially
proposed for approximating polynomial kernel. Explicit feature maps are
further put forward for addictive kernels [23] and the Gaussian kernels on the
sphere [24]. However, there exists relatively large variance for these methods
to approximate non-stationary kernels. For indefinite kernels, Pennington et
al. [20] and Liu et al. [25] firstly consider the unification for non-
stationary kernels and stationary kernels, and approximate stationary
indefinite kernels with Gaussian Mixture Model. Nevertheless, Gaussian Mixture
Model is a biased approximation. A recent work [26] utilizes generalized
measure to achieve unbiased approximation for stationary indefinite kernels
but smaller variance is expected to be obtained.

In this work, in order to provide an unbiased estimation with lower variance
for stationary indefinite kernels, we propose the generalized orthogonal
random features (GORF). Orthogonality is introduced to achieve lower variance
for unbiased random approximation towards stationary indefinite kernels. We
answer the question that how much variance reduction would be obtained through
GORF from the theoretical analysis and experimentally demonstrate that GORF
achieves variance reduction effect and lower approximation error. Notice that
different to imposing orthogonality on two random matrices separately, we
achieve orthogonality between two random matrices.

The contributions of this work are summarized as follows:

$\bullet$ We propose GORF, an unbiased approximation method with lower
variance for stationary indefinite kernels, and extend the orthogonality
method to the filed of stationary indefinite kernels.

$\bullet$ We theoretically examine the unbiasedness of GORF method, and derive
the reduced variance after applying orthogonality to random features for
indefinite kernels.

$\bullet$ We validate our GORF method on two typical indefinite kernels
(polynomial kernels on the unit sphere and delta-gaussian kernels) and several
datasets, demonstrating that our method achieves lower variance and
approximation error. In addtition, under the framework of support vector
machine and support vector regression, our method achieves better
classification accuracy and regression ability compared with the existing
kernel approximation methods.

##  II Preliminaries

###  II-A Random Fourier Features

This part briefly introduces Bochner’s theorem [27] and the corresponding
Random Fourier Features. Let $\mathcal{D}=\\{\boldsymbol{x_{i}}\\}_{i=1}^{N}$
be the dataset with $N$ training samples with
$\boldsymbol{x}_{i}\in\mathbb{R}^{d}$. Let $k(\cdot,\cdot)$ be the kernel
function and $K=[k(\boldsymbol{x}_{i},\boldsymbol{x}_{j})]_{N\times N}$ be the
kernel matrix sampled from $\mathcal{D}$. Then the Bochner’s theorem is
described as follows.

###### Theorem 1

(Bochner’s theorem [27]) A continuous and stationary function
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow{\mathbb{R}}$ is positive
definite if and only if it can be represented as

$\begin{split}k(\boldsymbol{x}-\boldsymbol{y})&=\int_{\mathbb{R}^{d}}\exp({\rm
i}\boldsymbol{w}^{T}(\boldsymbol{x}-\boldsymbol{y}))p(d\boldsymbol{w})\\\
&=\mathbb{E}_{\boldsymbol{w}\sim p(\boldsymbol{w})}[\exp({\rm
i}\boldsymbol{w}^{T}(\boldsymbol{x}-\boldsymbol{y})],\end{split}$ (1)

where $p(\boldsymbol{w})$ is the positive finite measure over
$\boldsymbol{w}$, i is the imaginary unit.

Bochner’s theorem indicates that the Fourier transform of a positive definite
function corresponds to a probability measure. Then $s$ random weights could
be sampled from the spectral distribution $p(\boldsymbol{w})$, notated as
$\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$. Monte-Carlo method is utilized to
numerically calculate the original Fourier integration.

$k(\boldsymbol{x}-\boldsymbol{y})\approx\frac{1}{s}\sum_{i=1}^{s}\exp(\rm{i}\boldsymbol{w}_{i}^{T}\boldsymbol{x})\exp(\rm{i}\boldsymbol{w}_{i}^{T}\boldsymbol{y})^{*}$
(2)

Since $k(.,.)$ is real-valued kernel function, the imaginary part of Equation
(2) is discarded:

$k(\boldsymbol{x}-\boldsymbol{y})\approx\psi(\boldsymbol{x})^{T}\psi(\boldsymbol{y})$
(3)

where
$\psi(\boldsymbol{x})=\frac{1}{\sqrt{s}}[\rm{cos}(\boldsymbol{w}_{1}^{T}\boldsymbol{x}),...,\rm{cos}(\boldsymbol{w}_{s}^{T}\boldsymbol{x}),\rm{sin}(\boldsymbol{w}_{1}^{T}\boldsymbol{x}),...,$
$\rm{sin}(\boldsymbol{w}_{s}^{T}\boldsymbol{x})]$. Then $\psi(\boldsymbol{x})$
is the Random Fourier Features for kernel function. Random Fourier Features is
unbiased estimation for positive definite kernels, and many works are aimed at
improving its approximation quality [9, 10, 12, 28].

###  II-B Signed Measure and Jordan Decomposition

Let $\mu:\mathcal{A}\rightarrow{[0,+\infty]}$ be a measure on a set $\Omega$
satisfying $\mu(\phi)=0$ and $\sigma$-additive. If $\mu(\Omega)=1$, it is a
probability measure. The Fourier transform for positive definite kernels is
probability measure. However, that for stationary indefinite kernel possesses
negative components. It is not a nonnegative Borel measure, or even not a
measure. To make Monte Carlo sampling applicable to stationary indefinite
kernels, a generalized version of measure allowing for negative values needs
to be introduced.

###### Definition 1

(Signed Measure [29]) Let $\Omega$ be some set, $\mathcal{A}$ be a
$\sigma$-algebra of subsets on $\Omega$. A signed measure is a function
$\mu:\mathcal{A}\rightarrow{[-\infty,+\infty)}$ or $(-\infty,+\infty]$
satisfying $\sigma$-additivity.

The definition of signed measure extends the traditional measure to a wider
range. And Jordan Decomposition below relates the traditionally nonnegative
measure with signed measure, which is a solid theoretical foundation for our
method.

###### Theorem 2

(Jordan Decomposition [30]) Let $\mu$ be a signed measure defined on the
$\sigma$-algebra $\mathcal{A}$ as given in Definition 1. There exist two
nonnegative measures $\mu_{+}$ and $\mu_{-}$ (one of them is finite measure)
such that $\mu=\mu_{+}-\mu_{-}$

Note that Jordan Decomposition is not unique. It can be characterized by the
Hahn decomposition theorem: the space $\Omega=\Omega_{+}\cup{\Omega_{-}}$,
where $\Omega_{+}$ is the positive set where $\mu(S)\geq 0$ for all subsets
$S\in\Omega_{+}$ and $\Omega_{-}$ is the negative set where $\mu(S)\leq 0$ for
all subsets $S\in\Omega_{-}$. The total mass is defined as:
$||\mu||=||\mu_{+}||+||\mu_{-}||$.

##  III Our Method

This section specifically illustrates our GORF method for stationary
indefinite kernels. Sec. III-A introduces the random features for stationary
indefinite kernels and verifies the unbiasedness. Sec. III-B demonstrates the
orthogonality process. Sec. III-C provides the theoretical derivation for
reduced variance after introducing orthogonality.

###  III-A Unbiased Generalized Random Features

Sec. II-B has mentioned the obstacle that Random Fourier Features are applied
to stationary indefinite kernel approximation. The spectrum distribution
$p(\boldsymbol{w})$ is not a probability measure and we could not directly
utilize Monte Carlo sampling. To make use of Monte Carlo sampling,
$p(\boldsymbol{w})$ is viewed as signed measure and can be decomposed into the
difference of two nonnegative measures via Jordan Decomposition.

$\begin{split}&k(\boldsymbol{x}-\boldsymbol{y})=k(\boldsymbol{z})\\\
&=\int_{\mathbb{R}^{d}}\exp({\rm
i}\boldsymbol{w}^{T}\boldsymbol{z})p(d\boldsymbol{w})\\\
&=\int_{\mathbb{R}^{d}}\exp({\rm
i}\boldsymbol{w}^{T}\boldsymbol{z})p_{+}(d\boldsymbol{w})-\int_{\mathbb{R}^{d}}\exp({\rm
i}\boldsymbol{\upsilon}^{T}\boldsymbol{z})p_{-}(d\boldsymbol{\upsilon})\\\
&=||p_{+}||\mathbb{E}_{\boldsymbol{w}\sim\widetilde{p}_{+}}(\exp({\rm
i}\boldsymbol{w}^{T}\boldsymbol{z}))-||p_{-}||\mathbb{E}_{\boldsymbol{\upsilon}\sim\widetilde{p}_{-}}(\exp({\rm
i}\boldsymbol{\upsilon}^{T}\boldsymbol{z})),\\\ \end{split}$ (4)

$p_{+}(\boldsymbol{w})$ and $p_{-}(\boldsymbol{\upsilon})$ are two nonnegative
finite measures, and $||p_{+}||$ and $||p_{-}||$ are total masses.
$\widetilde{p}_{+}$ and $\widetilde{p}_{-}$ are corresponding normalized
probability measures, where
$\widetilde{p}_{+}=\frac{p_{+}(\boldsymbol{w})}{||p_{+}||}$ and
$\widetilde{p}_{-}=\frac{p_{-}(\boldsymbol{\upsilon})}{||p_{-}||}$. Notice
that we suppose the stationary kernels mentioned in this paper are radial,
where $k(\boldsymbol{z})=k(||\boldsymbol{z}||)$, and accordingly their Fourier
transforms are radial, i.e., $p(\boldsymbol{w})=p(||\boldsymbol{w}||)$. Then
the calculation of total masses is shown as follows:

$||p_{+}||=\int_{0}^{\infty}|p_{+}(||\boldsymbol{w}||)|dw,\quad||p_{-}||=\int_{0}^{\infty}|p_{-}(||\boldsymbol{w}||)|dw$
(5)

According to Theorem 1, two expectations in Equation (LABEL:decompose) are
replaced by two PD function $\widetilde{k}_{+}(\boldsymbol{z})$ and
$\widetilde{k}_{-}(\boldsymbol{z})$, and then we obtain
$k_{+}(\boldsymbol{z})$ and $k_{-}(\boldsymbol{z})$ after multiplying
$||p_{+}||$ and $||p_{-}||$.

$\begin{split}k(\boldsymbol{x}-\boldsymbol{y})&=k(\boldsymbol{z})\\\
&=||p_{+}||\widetilde{k}_{+}(\boldsymbol{z})-||p_{-}||\widetilde{k}_{-}(\boldsymbol{z})\\\
&=k_{+}(\boldsymbol{z})-k_{-}(\boldsymbol{z})\end{split}$ (6)

Define
$\phi(\boldsymbol{x})=\frac{1}{\sqrt{s}}[\psi_{1}(\boldsymbol{x}),\psi_{2}(\boldsymbol{x}),...,\psi_{s}(\boldsymbol{x})]^{T}$
with $\psi_{i}(\boldsymbol{x})$ as:

$\begin{split}\psi_{i}(\boldsymbol{x})=[&\sqrt{||p_{+}||}{\rm
cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{x}),\sqrt{||p_{+}||}{\rm
sin}(\boldsymbol{w}_{i}^{T}\boldsymbol{x}),\\\ &{\rm i}\sqrt{||p_{-}||}{\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{x}),{\rm i}\sqrt{||p_{-}||}{\rm
sin}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{x})].\end{split}$ (7)

where random weights $\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$ and
$\\{\boldsymbol{\upsilon}_{i}\\}_{i=1}^{s}$ are independently sampled from
$\widetilde{p}_{+}$ and $\widetilde{p}_{-}$ (also called i.i.d sampling). Then
the original kernel function could be represented as:

$k(\boldsymbol{x}-\boldsymbol{y})\approx\frac{1}{s}\sum_{i=1}^{s}\langle\psi_{i}(\boldsymbol{x}),\psi_{i}(\boldsymbol{y})\rangle=\phi(\boldsymbol{x})^{T}\phi(\boldsymbol{y}).$
(8)

$\phi(\boldsymbol{x})$ is the generalized random features (GRFF) for
stationary indefinite kernels. Different from RFF for PD kernels, we construct
the random Fourier maps for stationary indefinite kernels in the complex space
and RFF for PD kernels is a special case for this random feature map.
Generalized random features demonstrate superiority in reducing the
computation time and storage space, which drops to $\mathcal{O}(2Ns)$ and
$\mathcal{O}(2Ns^{2})$ respectively.

Equation (LABEL:decompose) seems to show that any stationary indefinite kernel
could be expressed by the difference of two positive definite kernels, but
some stationary kernels do not possess this positive decomposition indeed.
[31] provides a RKKS judgement for the existence of positive decomposition but
this simply applies to some intuitive examples: a linear combination of PD
kernels. By virtue of measure decomposition of the signed measure, a necessary
and sufficient condition to the existence of positive decomposition is derived
for more extensive kernels.

###### Theorem 3

[26] Assume that an indefinite kernel is stationary, and its (generalized)
Fourier transform is denoted by the measure $\mu$, then $k$ possesses positive
decomposition if and only if the total mass of the measure $\mu$ except for
the origin 0 is finite, i.e., $||\mu||\leq\infty$.

Unbiasedness is the fundamental characteristics in kernel approximation. We
examine the unbiasedness of generalized random features as Theorem 4
demonstrates:

###### Theorem 4

$K_{\rm GRFF}(\boldsymbol{z})$ is unbiased for all stationary kernels, i.e.,
$\mathbb{E}(K_{{\rm GRFF}}(\boldsymbol{z}))=k(\boldsymbol{z})$. And the
variance is:
$Var(K_{GRFF}(\boldsymbol{z}))=\frac{||p_{+}||^{2}}{s}[\frac{1+\widetilde{k}_{+}(2\boldsymbol{z})}{2}-\widetilde{k}_{+}^{2}(\boldsymbol{z})]+\frac{||p_{-}||^{2}}{s}[\frac{1+\widetilde{k}_{-}(2\boldsymbol{z})}{2}-\widetilde{k}_{-}^{2}(\boldsymbol{z})]$

Proof. $\mathbb{E}(K_{{\rm
GRFF}}(\boldsymbol{z}))=\mathbb{E}(\sum_{i=1}^{s}\frac{1}{s}||p_{+}||{\rm
cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{z})-\sum_{i=1}^{s}\frac{1}{s}||p_{-}||{\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z}))=\sum_{i=1}^{s}\frac{||p_{+}||}{s}\mathbb{E}({\rm
cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{z}))-\sum_{i=1}^{s}\frac{||p_{-}||}{s}\mathbb{E}({\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z}))$.

Based on the definition of GRFF,
$\boldsymbol{w}_{i}\sim\widetilde{p}_{+}(\boldsymbol{w})$,
$\boldsymbol{\upsilon}_{i}\sim\widetilde{p}_{-}(\boldsymbol{\upsilon})$ and
Bochner’s theorem, we get $\mathbb{E}({\rm
cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{z}))=\widetilde{k}_{+}(\boldsymbol{z})$,
$\mathbb{E}({\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z}))=\widetilde{k}_{-}(\boldsymbol{z})$.

Then $\mathbb{E}(K_{{\rm
GRFF}}(\boldsymbol{z}))=\sum_{i=1}^{s}\frac{||p_{+}||}{s}\widetilde{k}_{+}(\boldsymbol{z})-\sum_{i=1}^{s}\frac{||p_{-}||}{s}\widetilde{k}_{-}(\boldsymbol{z})$
$=||p_{+}||\widetilde{k}_{+}(\boldsymbol{z})-||p_{-}||\widetilde{k}_{-}(\boldsymbol{z})=k(\boldsymbol{z})$.
This verifies that generalized random feature maps is unbiased for all
stationary kernels.

Since $\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$ and
$\\{\boldsymbol{\upsilon}_{i}\\}_{i=1}^{s}$ are i.i.d sampled, $Var(K_{\rm
GRFF}(\boldsymbol{z}))=Var[\sum_{i=1}^{s}\frac{1}{s}||p_{+}||\rm{cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{z})-\sum_{i=1}^{s}\frac{1}{s}||p_{-}||{\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z})]=\frac{||p_{+}||^{2}}{s^{2}}\sum_{i=1}^{s}Var[{\rm
cos}(\boldsymbol{w}_{i}^{T}\boldsymbol{z})]+\frac{||p_{-}||^{2}}{s^{2}}\sum_{i=1}^{s}Var[{\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z})]=\frac{||p_{+}||^{2}}{s}(\mathbb{E}[{\rm
cos}^{2}(\boldsymbol{w}_{1}^{T}\boldsymbol{z})]-\widetilde{k}_{+}^{2}(\boldsymbol{z}))+\frac{||p_{-}||^{2}}{s}(\mathbb{E}[{\rm
cos}^{2}(\boldsymbol{\upsilon}_{1}^{T}\boldsymbol{z})]-\widetilde{k}_{-}^{2}(\boldsymbol{z}))=\frac{||p_{+}||^{2}}{s}[\frac{1+\widetilde{k}_{+}(2\boldsymbol{z})}{2}-\widetilde{k}_{+}^{2}(\boldsymbol{z})]+\frac{||p_{-}||^{2}}{s}[\frac{1+\widetilde{k}_{-}(2\boldsymbol{z})}{2}-\widetilde{k}_{-}^{2}(\boldsymbol{z})]$

###  III-B Orthogonality

Although generalized random feature maps are unbiased, there exists large
variance shown in Theorem 4. In unbiased approximation, the variance
determines the overall approximation error. The sampling method is the main
cause for large variance. GORF method utilizes orthogonal sampling, aimed at
making random weight vectors not i.i.d sampled from the spectrum and obtaining
lower variance.

Let the sampled random weights $\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$ and
$\\{\boldsymbol{\upsilon}_{i}\\}_{i=1}^{s}$ constitute the columns of random
matrices $W_{\rm
pos}=[\boldsymbol{w}_{1},\boldsymbol{w}_{2},...,\boldsymbol{w}_{s}]$ and
$W_{\rm
neg}=[\boldsymbol{\upsilon}_{1},\boldsymbol{\upsilon}_{2},...,\boldsymbol{\upsilon}_{s}]$.
At the beginning of GORF method, we firstly sample $c$ $\ell_{2}$-norm of
random weights from the corresponding probability measures
$\widetilde{p}_{+}(\boldsymbol{w})$ and
$\widetilde{p}_{-}(\boldsymbol{\upsilon})$ separately,

$||\boldsymbol{w}_{i}||_{2}\sim\widetilde{p}_{+}(||\boldsymbol{w}||),\quad||\boldsymbol{\upsilon}_{i}||_{2}\sim\widetilde{p}_{-}(||\boldsymbol{\upsilon}||).$
(9)

The direction vector of each random weight is sampled from a normalized
Gaussian distribution, notated as $\\{\boldsymbol{a}_{i}\\}_{i=1}^{m}$ and
$\\{\boldsymbol{b}_{i}\\}_{i=1}^{m}$. And they construct the random matrices
$M$. Subsequently, QR decomposition is conducted to obtain orthogonal
direction vectors (Notice that the size of orthogonal matrix $2m$ takes the
maximum value between $2s$ and $2d$ ($2m={\rm max}(2s,2d)$))

$\begin{split}&\boldsymbol{a}_{j}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{2m}),\quad\boldsymbol{b}_{j}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{2m})\\\
&M=[\boldsymbol{a}_{1},...,\boldsymbol{a}_{m},\boldsymbol{b}_{1},...,\boldsymbol{b}_{m}]\\\
&M^{orth}={\rm QR}(M).\\\ \end{split}$ (10)

We retain the first $d$ rows after obtaining the orthogonal matrix $M^{orth}$,
and normalize each columns, notated as $M^{orthn}$. $\ell$2-norm and the
corresponding normalized orthogonal direction vector constitute the column
vectors for two random matrices $W_{pos}$ and $W_{neg}$,

$\boldsymbol{w}_{i}=||\boldsymbol{w}_{i}||_{2}M^{orthn}_{i},\quad\boldsymbol{\upsilon}_{i}=||\boldsymbol{\upsilon}_{i}||_{2}M^{orthn}_{s+i}.$
(11)

Algorithm 1 specifically illustrates the overall construction of GORF.

Algorithm 1 The Construction of Generalized Orthogonal Random Features

0: A stationary kernel function $k(\boldsymbol{x},\boldsymbol{y})=k(z)$,
$z=||\boldsymbol{x}-\boldsymbol{y}||$, the dimension of training samples $d$,
the number of Generalized Orthogonal Random Features $s$

0: Generalized Orthogonal Random Feature Maps $\phi(\boldsymbol{x})$ such that
$k(\boldsymbol{x}-\boldsymbol{y})\approx\phi(\boldsymbol{x})^{T}\phi(\boldsymbol{y})$

1: Obtain the measure $p(\boldsymbol{w})$ through Fourier transform, and
calculate $p_{+}(\boldsymbol{w})$, $p_{-}(\boldsymbol{w})$, $||p_{+}||$,
$||p_{-}||$, $\widetilde{p}_{+}(\boldsymbol{w})$,
$\widetilde{p}_{-}(\boldsymbol{w})$ through Jordan Decomposition and
normalization.

2: Individually sample $s$ $\ell$2-norms $||\boldsymbol{w}_{i}||_{2}$ and
$||\boldsymbol{\upsilon}_{i}||_{2}$ from $\widetilde{p}_{+}(\boldsymbol{w})$
and $\widetilde{p}_{-}(\boldsymbol{w})$

3: Sample $2m$ direction vectors $\boldsymbol{a}_{j}$ and $\boldsymbol{b}_{j}$
from $\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{2m})$. Use QR decomposition
to get orthogonal matrix $M^{orth}$ and $M^{orthn}$.

4: Use Equation (11) to obtain $\boldsymbol{w}_{i}$ and
$\boldsymbol{\upsilon}_{i}$ and form the random matrices $W_{\rm pos}$ and
$W_{\rm neg}$

5: Output the explicit feature mapping $\phi(\boldsymbol{x})$ associated
$\psi_{i}(\boldsymbol{x})$ in Equation (7)

###  III-C Theoretical Analysis

How much variance reduction could be achieved by introducing orthogonal
sampling? It is a crucial problem for GORF method. This section concentrates
on the theoretical derivation of variation reduction for orthogonal sampling.
In order to formulate the variance difference between applying the orthogonal
and i.i.d sampling, we firstly review the variance difference between Random
Fourier Features approximation and Orthogonal Random Features approximation
(ORF) for PD radial kernel.

###### Theorem 5

[13] For a PD radial kernel $k$ on $\mathbb{R}^{d}$ with Fourier measure
$p(\boldsymbol{w})$ and $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$,
writing $\boldsymbol{z}=\boldsymbol{x}-\boldsymbol{y}$, we have:

$\begin{split}G_{k}(\boldsymbol{z})&=Var(K_{{\rm
ORF}}(\boldsymbol{z}))-Var(K_{{\rm RFF}}(\boldsymbol{z}))\\\
&=\frac{s-1}{s}\mathbb{E}_{R_{1}}\left[\frac{J_{\frac{d}{2}-1}(R_{1}||\boldsymbol{z}||)\Gamma(d/2)}{(R_{1}||\boldsymbol{z}||/2)^{\frac{d}{2}-1}}\right]^{2}-\\\
&\frac{s-1}{s}\mathbb{E}_{R_{1},R_{2}}\left[\frac{J_{\frac{d}{2}-1}(\sqrt{R_{1}^{2}+R_{2}^{2}}||\boldsymbol{z}||)\Gamma(d/2)}{(\sqrt{R_{1}^{2}+R_{2}^{2}}||\boldsymbol{z}||/2)^{\frac{d}{2}-1}}\right],\end{split}$
(12)

where $R_{1},R_{2}\sim p(\boldsymbol{w})$, and $J_{\alpha}$ is the Bessel
function of the first kind of degree $\alpha$

Based on Theorem 5, we formulate the reduced variance after orthogonal
sampling for stationary indefinite kernels.

###### Theorem 6

$K_{\rm GRFF}(\boldsymbol{z})$ and $K_{\rm GORF}(\boldsymbol{z})$ are unbiased
approximation with i.i.d sampling and orthogonal sampling for stationary
indefinite kernels. For any stationary indefinite kernel $k$ on
$\mathbb{R}^{d}$ with Fourier signed measure $p(\boldsymbol{w})$, the
corresponding decomposition is $||p_{+}||\widetilde{k}_{+}$ and
$||p_{+}||\widetilde{k}_{-}$. Then we have:

$\begin{split}Var[K_{{\rm GORF}}(\boldsymbol{z})]-Var[K_{{\rm
GRFF}}(\boldsymbol{z})]=&||p_{+}||^{2}G_{\widetilde{k}_{+}}(\boldsymbol{z})+||p_{-}||^{2}\\\
&G_{\widetilde{k}_{-}}(\boldsymbol{z})+H(\boldsymbol{z}),\end{split}$ (13)

where $G_{\widetilde{k}_{+}}(\boldsymbol{z})$,
$G_{\widetilde{k}_{-}}(\boldsymbol{z})$ are calculated by Theorem 5, and
$H(\boldsymbol{z})=2||p_{+}||||p_{-}||[\mathbb{E}(a_{1})\mathbb{E}(b_{1})-\mathbb{E}(a_{1}b_{1})]$,
$a_{i}={\rm cos}(\boldsymbol{w}_{i}^{T}$ $\boldsymbol{z}),b_{i}={\rm
cos}(\boldsymbol{\upsilon}_{i}^{T}\boldsymbol{z})$

Proof. Since $\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$ and
$\\{\boldsymbol{\upsilon}_{i}\\}_{i=1}^{s}$ are not i.i.d sampled and
correlate with each other after orthogonalization,

$\begin{split}&Var(K_{\rm
GORF}(\boldsymbol{z}))=Var\left[\sum_{i=1}^{s}\frac{1}{s}||p_{+}||a_{i}-\sum_{i=1}^{s}\frac{1}{s}||p_{-}||b_{i}\right]\\\
&=\mathbb{E}\left[\left(\sum_{i=1}^{s}\frac{1}{s}||p_{+}||a_{i}-\sum_{i=1}^{s}\frac{1}{s}||p_{-}||b_{i}\right)^{2}\right]-\mathbb{E}^{2}\left[\sum_{i=1}^{s}\frac{1}{s}||p_{+}||a_{i}-\right.\\\
&\left.\sum_{i=1}^{s}\frac{1}{s}||p_{-}||b_{i}\right]\\\
&=\frac{||p_{+}||^{2}}{s^{2}}\sum_{i=1}^{s}\left[\mathbb{E}^{2}(a_{i})-\mathbb{E}(a_{i}^{2})\right]+\frac{||p_{-}||^{2}}{s^{2}}\sum_{i=1}^{s}\left[\mathbb{E}^{2}(b_{i})-\mathbb{E}(b_{i}^{2})\right]\\\
&+\frac{||p_{+}||^{2}}{s^{2}}\sum_{i,j=1,i\neq
j}^{s}[\mathbb{E}(a_{i}a_{j})-\mathbb{E}(a_{i})\mathbb{E}(a_{j})]+\frac{||p_{-}||^{2}}{s^{2}}\sum_{i,j=1,i\neq
j}^{s}\\\
&[\mathbb{E}(b_{i}b_{j})-\mathbb{E}(b_{i}^{2})]+H(\boldsymbol{z})\end{split}$
(14)

Since $\\{\boldsymbol{w}_{i}\\}_{i=1}^{s}$ and
$\\{\boldsymbol{\upsilon}_{i}\\}_{i=1}^{s}$ are sampled from
$\widetilde{p}_{+}(\boldsymbol{w})$ and $\widetilde{p}_{-}(\boldsymbol{w})$,
according to Theorem 4 and 5, then

$\begin{split}&Var(K_{\rm GORF}(\boldsymbol{z}))=Var[K_{{\rm
GRFF}}(\boldsymbol{z})]+||p_{+}||^{2}G_{\widetilde{k}_{+}}(\boldsymbol{z})+\\\
&||p_{-}||^{2}G_{\widetilde{k}_{-}}(\boldsymbol{z})+H(\boldsymbol{z})\end{split}$
(15)

For $H(\boldsymbol{z})$ in Equation (15), we could further simply this term:

$\begin{split}H(\boldsymbol{z})&=\frac{2}{s^{2}}||p_{+}||||p_{-}||[\mathbb{E}[\sum_{i=1}^{s}a_{i}]\mathbb{E}[\sum_{j=1}^{s}b_{j}]-\mathbb{E}[\sum_{i,j=1}^{s}a_{i}b_{j}]]\\\
&=2||p_{+}||||p_{-}||[\mathbb{E}(a_{1})\mathbb{E}(b_{1})-\mathbb{E}(a_{1}b_{1})]\end{split}$
(16)

Equation (13) divides the variance difference into two parts. The first part
$G_{\widetilde{k}_{+}}(\boldsymbol{z})+G_{\widetilde{k}_{-}}(\boldsymbol{z})$
represents the reduced variance when imposing orthogonality on $W_{\rm pos}$
and $W_{\rm neg}$ separately. The second term $H(\boldsymbol{z})$ introduces
the reduced variance after making the columns of two random matrices
orthogonal with each other.
$G_{\widetilde{k}_{+}}(\boldsymbol{z})+G_{\widetilde{k}_{-}}(\boldsymbol{z})+H(\boldsymbol{z})\leq
0$ corresponds to better performance of GORF over GRFF at $\boldsymbol{z}$.

##  IV Experiment

In this section the performance of our GORF method is evaluated, from the
aspects of variance reduction, approximation error and performance in specific
tasks like classification and regression. All the experiments are conducted on
a standard PC with i7-4970 CPU(3.6GHz) and 32G RAM in
Matlab111https://github.com/tomqingo/GORF.

###  IV-A Experimental Setup

The approximation to the stationary indefinite kernels is the main target in
our work. We evaluate our GORF method on two typical stationary indefinite
kernels and four datasets in our experiment.

####  IV-A1 kernel

We select two typical stationary indefinite kernels in our research,
polynomial kernel on the unit sphere and delta-gaussian kernel.

Polynomial kernel is the mostly used in kernel methods. In order to extend
Random Fourier Features to the polynomial kernel for acceleration, the
polynomial kernel is converted to a shift-invariant but indefinite kernels
when defined on $\mathcal{S}^{d-1}\times\mathcal{S}^{d-1}$
($||\boldsymbol{x}||_{2}=||\boldsymbol{y}||_{2}=1$).

$k(\boldsymbol{x},\boldsymbol{y})=(1-\frac{||\boldsymbol{x}-\boldsymbol{y}||^{2}}{a^{2}})^{m}=\alpha(q+\langle\boldsymbol{x},\boldsymbol{y}\rangle)^{m}$
(17)

where $q=a^{2}/2-1$, $\alpha=(2/a^{2})^{m}$. The spectrum distribution of
polynomial kernel on the sphere is shown as follows:

$p(\boldsymbol{w})=\sum_{i=0}^{p}\frac{m!}{(m-i)!}(1-\frac{4}{a^{2}})^{m-i}(\frac{2}{a^{2}})^{i}(\frac{2}{||\boldsymbol{w}||_{2}})^{\frac{d}{2}+i}J_{\frac{d}{2}+i}(2||\boldsymbol{w}||_{2})$
(18)

Delta-gaussian kernel is the linear combination of commonly used Gaussian
kernels, with arbitrary coefficients $\\{a_{i}\\}_{i=1}^{m}$.

$k(\boldsymbol{x}-\boldsymbol{y})=k(\boldsymbol{z})=\sum_{i=1}^{m}a_{i}e^{-\frac{||\boldsymbol{z}||^{2}}{2\sigma_{i}^{2}}}$
(19)

The corresponding spectrum distribution is shown below:

$p(\boldsymbol{w})=\sum_{i=1}^{m}a_{i}\sigma_{i}^{d}e^{-\frac{\sigma_{i}^{2}||\boldsymbol{w}||^{2}}{2}}$
(20)

In the implementation, $p_{+}(\boldsymbol{w})=\max(0,p(\boldsymbol{w}))$,
$p_{-}(\boldsymbol{w})=\max(0,-p(\boldsymbol{w}))$.

####  IV-A2 Datasets

Four representative benchmark datasets are adpoted.
$letter$222http://archive.ics.uci.edu/ml/datasets.php,
$ijcnn1$333https://www.csie.ntu.edu.tw/cjlin/libsvmtools/datasets/,
$usps$444https://www.kaggle.com/bistaumanga/usps-dataset are the
classification datasets with the data quantity and dimension shown in Table I.
$housing$555https://archive.ics.uci.edu/ml/machine-learning-databases/housing/
is the ordinary regression dataset for Boston housing prices with 506 samples
and 13 input dimensions. 405 samples are randomly selected for training set
and the left 101 samples constitute testing set. For all the datasets, all the
inputs are normalized to $[0,1]^{d}$.

TABLE I: BENCHMARK CLASSIFICATION DATASETS Datasets | d | training | testing
---|---|---|---
letter | 16 | 12000 | 6000
ijcnn1 | 22 | 49990 | 91701
usps | 256 | 7291 | 2007

###  IV-B Variance Reduction

Orthogonality is the crucial part in our GORF algorithm aimed at reducing
variance. Sec. III-C has derived the reduced variance after applying
orthogonality to generalized random features. Since the expectation in
variance reduction formula could not be calculated analytically, we then
numerically calculate the reduced variance for polynomial kernels on the unit
sphere and delta-gaussian kernel. Monte-Carlo method is utilized to calculate
the expectation using 10000 random weights sampled from
$\widetilde{p}_{+}(\boldsymbol{w})$ and
$\widetilde{p}_{-}(\boldsymbol{\upsilon})$ respectively.

(a) polynomial kernel on the unit sphere ($a=3,m=1$) (b) delta-gaussian kernel
($a_{1}=1,a_{2}=-1,\sigma_{1}=1,\sigma_{2}=10)$

Figure 1: Numerically calculated variance reduction for polynomial kernel on
the unit sphere and delta-gaussian kernel under different data dimensions.
Observe that the variance difference is negative, orthogonality achieves lower
variance for a wide range of mean square distance
($z=||\boldsymbol{x}-\boldsymbol{y}||_{2}$, $s=d$)

Fig. 1 shows the variation reduction for polynomial kernel on the unit sphere
and delta-gaussian kernel. The horizontal axis is the mean square distance
between two samples. It could be observed that orthogonality has a prominent
variance reduction effect for a wide range of mean square distance. Although
the variance reduction seems less significant as the data dimension increases
in Fig. 1, the relative value for variance reduction remains large for large
data dimension, since the variance for i.i.d sampling decreases as the data
dimension increases.

(a) $letter$

(b) $ijcnn1$

(c) $usps$

Figure 2: Comparisons of various algorithms for kernel approximation in terms of approximation error across two typical stationary indefinite kernels and three datasets with different dimensions. Top: polynomial kernel on the unit sphere. Below: delta-gaussian kernel TABLE II: COMPARISON RESULTS BETWEEN APPLYING ORTHOGONAL SAMPLING AND I.I.D SAMPLING ON STATIONARY INDEFINITE KERNELS IN TERMS OF APPROXIMATION ERROR (MEAN$\pm$STD.). THE LOWEST ERROR IS HIGHLIGHTED IN BOLDFACE Kernel | DataSet | Method | s=1/2d | s=d | s=2d | s=8d
---|---|---|---|---|---|---
polynomial | $letter$ | GRFF | 0.0859$\pm$0.0309 | 0.0547$\pm$0.0078 | 0.0469$\pm$0.0109 | 0.0261$\pm$0.0059
GORF | 0.0716$\pm$0.0175 | 0.0495$\pm$0.0139 | 0.0360$\pm$0.0110 | 0.0231$\pm$0.0078
$ijcnn1$ | GRFF | 0.1159$\pm$0.0158 | 0.0907$\pm$0.0194 | 0.0794$\pm$0.0142 | 0.0433$\pm$0.0059
GORF | 0.1072$\pm$0.0228 | 0.0775$\pm$0.0204 | 0.0487$\pm$0.0155 | 0.0397$\pm$0.0135
$usps$ | GRFF | 0.0270$\pm$0.0056 | 0.0213$\pm$0.0063 | 0.0160$\pm$0.0029 | 0.0137$\pm$0.0020
GORF | 0.0251$\pm$0.0078 | 0.0194$\pm$0.0087 | 0.0143$\pm$0.0030 | 0.0137$\pm$0.0016
delta-gaussian | $letter$ | GRFF | 0.3918$\pm$ 0.0428 | 0.2736$\pm$0.0345 | 0.1887$\pm$0.0201 | 0.1017$\pm$0.0088
GORF | 0.3154$\pm$0.0424 | 0.1133$\pm$0.0181 | 0.0760$\pm$0.0090 | 0.0376$\pm$0.0039
$ijcnn1$ | GRFF | 0.2924$\pm$0.0188 | 0.2171$\pm$0.0222 | 0.1504$\pm$0.0134 | 0.0757$\pm$0.0081
GORF | 0.2415$\pm$0.0190 | 0.1026$\pm$0.0129 | 0.0739$\pm$0.0065 | 0.0383$\pm$0.0022
$usps$ | GRFF | 0.1005$\pm$0.0061 | 0.0690$\pm$0.0050 | 0.0500$\pm$0.0024 | 0.0253$\pm$0.0023
GORF | 0.0724$\pm$0.0049 | 0.0235$\pm$0.0009 | 0.0166$\pm$0.0008 | 0.0083$\pm$0.0003

###  IV-C Approximation Error

Approximation error is the most important issue for RFF methods. Unbiasedness
and variance determine the overall approximation error. We compare our method
with other existing methods to demonstrate the superiority of reducing
approximation error in this part and investigate the effect of orthogonality
on reduction of approximation error for ablation study.

In this part, to quantify the approximation error, relative error
$\frac{||K-\hat{K}||_{F}}{||K||_{F}}$ is calculated where $K$ and $\hat{K}$
denote the exact kernel matrix and and its approximated kernel matrix on 1000
randomly selected samples. The parameters in polynomial kernel on the unit
sphere and delta-gaussian kernel are the same as that in Sec. IV-B. The
approximation error is represented by mean and standard deviation in 10
repetitive experiments.

####  IV-C1 comparison with existing methods

In this part, several kernel approximation baseline methods are utilized for
comparison in terms of approximation error, including Random Maclaurin (RM)
[21], Tensor Sketch (TS) [22], Spherical Random Features (SRF) [20] and
Double-Infinite Gaussian Mixtures Model (DIGMM) [25]. RM and TS are the
unbiased estimations for polynomial kernel. SRF and DIGMM are two biased
approximation for arbitrary stationary kernels using the Gaussian Mixture
Model (GMM).

The comparison between our GORF method and other existing method in terms of
approximation ability is demonstrated in Fig. 2. It shows that our method
achieves lowest approximation error. The lowest approximation error is
attributed to two important factors: 1) unbiasness. Our GORF method is
unbiased compared with the current stationary indefinite kernel

(a) $letter$

(b) $ijcnn1$

(c) $usps$

Figure 3: Comparisons of various algorithms for SVM classification task in
terms of accuracy across two typical stationary indefinite kernels and three
datasets with different dimensions. Top: polynomial kernel on the unit sphere.
Below: delta-gaussian kernel

approximation methods, like SRF and DIGMM. Especially for delta-gaussian
kernel whose positive and negative measure are almost the same, our GORF
method retains both of these two parts but SRF method simply retains either
the positive or the negative part, resulting a significant drop in
approximation error for our method. 2) low variance. Our method achieves lower
variance based on the unbiasedness. Our GORF method achieves faster
convergence rate than other unbiased methods like RM and TS in polynomial
kernel approximation because of low variance brought by orthogonality and RFF.
Notice that for the baseline method DIGMM, it requires $s\geq 2d$ because of
the intolerance of the non-PD matrices in Cholosky decomposition and the
complex value in polygamma function computation occurred in $s\textless 2d$.

####  IV-C2 effect of orthogonality

We evaluate orthogonality, the crucial part in GORF algorithm, in terms of
approximation error for further ablation study in this part . Table II
demonstrates that orthogonal sampling brings a further substantial
approximation error reduction which is averagely above 10$\%$ drop because of
the significant decrease in variance. For delta-gaussian kernel, the decrease
in approximation error becomes more significant, around 20$\%$ drop.

###  IV-D Performance in classification and regression

We compare our GORF method and other existing methods when considering the
specific support vector classification and regression task. In classification
experiment, datasets $letter$, $ijcnn1$ and $usps$ are utilized. and
classification accuracy is the evaluation metric. In regression experiment,
the dataset $housing$ is used and root mean square error (RMSE) between ground
truth and predictions is utilized as evaluation metric. We choose $liblinear$,
a toolkit for support vector machine and support vector regression in our
experiment. The parameter $C$ in $liblinear$ [32] is set as 1000. All the
other hyperparameters is the same as that in the experiments above.

In the classification experiment, we could observe that our GORF method
achieves better classification accuracy than the existing stationary
indefinite kernel approximation methods for almost all the number of random
features in Fig. 3. In the regression task, Table III and IV demonstrates that
our GORF method achieves lower RMSE and possesses better ability to predict
the housing prices in testing set. These results could be attributed to better
approximation to the originally selected kernels. Better approximation
contributes to better performance in classification and regression. And
unbiasedness and lower variance are the main factors of lower approximation
error in our GORF algorithm.

TABLE III: REGRESSION ERROR FOR KERNEL APPROXIMATION METHODS ON POLYNOMIAL KERNEL AND $HOUSING$ DATASET (RMSE: MEAN$\pm$STD). THE LOWEST ERROR IS HIGHLIGHTED IN BOLDFACE Methods | s=2d | s=4d | s=8d
---|---|---|---
RM | 7.153$\pm$1.772 | 5.436$\pm$0.917 | 4.491$\pm$0.008
TS | 5.414$\pm$0.879 | 4.772$\pm$0.177 | 4.657$\pm$0.316
SRF | 4.391$\pm$0.368 | 3.906$\pm$0.219 | 3.555$\pm$0.130
DIGMM | 4.897$\pm$0.368 | 4.130$\pm$0.324 | 4.000$\pm$0.475
GORF(OURS) | 4.079$\pm$0.233 | 3.817$\pm$0.204 | 3.472$\pm$0.137

##  V Conclusion

In this paper, we propose generalized orthogonal random features, an unbiased
estimation with lower variance for stationary indefinite kernel approximation.
We theoretically verify the unbiasedness and derive the variance reduction.
Experimental results show that our proposed algorithm achieves lower variance
and approximation error. In specific support vector classification and
regression problem, our method obtains better classification and regression
ability compared with the existing methods.

TABLE IV: REGRESSION ERROR FOR KERNEL APPROXIMATION METHODS ON DELTA-GAUSSIAN KERNEL AND $HOUSING$ DATASET (RMSE: MEAN$\pm$STD). THE LOWEST ERROR IS HIGHLIGHTED IN BOLDFACE Methods | s=2d | s=4d | s=8d
---|---|---|---
SRF | 5.432$\pm$0.729 | 3.845$\pm$0.379 | 3.321$\pm$0.274
GORF(OURS) | 3.739$\pm$0.360 | 3.474$\pm$0.330 | 3.164$\pm$0.452

## Acknowledgment

This work is supported by National Key Research Development Project (No.
2018AAA0100702) and National Natural Science Foundation of China (Nos.
61977046, 61876107, U1803261) and Committee of Science of Technology,
Shanghai, China (No. 19510711200).

## References

  * [1] C. Cortes and V. N. Vapnik, “Support-vector networks,” _Machine Learning_ , vol. 20, no. 3, pp. 273–297, 1995.
  * [2] B. Schlkopf, A. J. Smola, and F. Bach, _Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond_. The MIT Press, 2018.
  * [3] H. Drucker, C. J. C. Burges, L. Kaufman, J. C. Chris, B. L. Kaufman, A. Smola, and V. Vapnik, “Support vector regression machines,” _Advances in Neural Information Processing Systems_ , vol. 28, no. 7, pp. 779–784, 1997.
  * [4] A. Wilson and R. Adams, “Gaussian process kernels for pattern discovery and extrapolation,” in _Proceedings of the 30th International Conference on Machine Learning_ , Atlanta, Georgia, USA, 17–19 Jun 2013, pp. 1067–1075.
  * [5] M. Lazaro-Gredilla, J. Quinonero-Candela, C. E. Rasmussen, and A. R. Figueiras-Vidal, “Sparse spectrum gaussian process regression,” _Journal of Machine Learning Research_ , vol. 11, pp. 1865–1881, 2010.
  * [6] B. Scholkopf and A. Smola, “Kernel principal component analysis,” in _International Conference on Artificial Neural Networks_ , 1997.
  * [7] C. M. Alaíz, M. Fanuel, and J. A. K. Suykens, “Convex formulation for kernel pca and its use in semisupervised learning,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 29, no. 8, pp. 3863–3869, 2018.
  * [8] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in _Advances in Neural Information Processing Systems 20_ , 2008, pp. 1177–1184.
  * [9] Q. Le, T. Sarlos, and A. Smola, “Fastfood - computing hilbert space expansions in loglinear time,” ser. Proceedings of Machine Learning Research, vol. 28, no. 3, Atlanta, Georgia, USA, 17–19 Jun 2013, pp. 244–252.
  * [10] J. Yang, V. Sindhwani, H. Avron, and M. Mahoney, “Quasi-monte carlo feature maps for shift-invariant kernels,” ser. Proceedings of Machine Learning Research, vol. 32, no. 1, Bejing, China, 22–24 Jun 2014, pp. 485–493.
  * [11] C. Feng, Q. Hu, and S. Liao, “Random feature mapping with signed circulant matrix projection,” in _Proceedings of the 24th International Conference on Artificial Intelligence_ , ser. IJCAI’15, 2015, p. 3490–3496.
  * [12] F. X. Yu, A. T. Suresh, K. M. Choromanski, D. N. Holtmann-Rice, and S. Kumar, “Orthogonal random features,” in _Advances in Neural Information Processing Systems 29_ , 2016, pp. 1975–1983.
  * [13] K. Choromanski, M. Rowland, T. Sarlos, V. Sindhwani, R. Turner, and A. Weller, “The geometry of random features,” ser. Proceedings of Machine Learning Research, vol. 84, Playa Blanca, Lanzarote, Canary Islands, 09–11 Apr 2018, pp. 1–9.
  * [14] K. Choromanski, M. Rowland, W. Chen, and A. Weller, “Unifying orthogonal Monte Carlo methods,” ser. Proceedings of Machine Learning Research, vol. 97, Long Beach, California, USA, 09–15 Jun 2019, pp. 1203–1212.
  * [15] M. Bojarski, A. Choromanska, K. Choromanski, F. Fagan, C. Gouy-Pailler, A. Morvan, N. Sakr, T. Sarlos, and J. Atif, “Structured adaptive and random spinners for fast machine learning computations,” ser. Proceedings of Machine Learning Research, vol. 54, Fort Lauderdale, FL, USA, 20–22 Apr 2017, pp. 1020–1029.
  * [16] A. J. Smola, Z. L. Óvári, and R. C. Williamson, “Regularization with dot-product kernels,” in _Advances in Neural Information Processing Systems 13_ , 2001, pp. 308–314.
  * [17] A. Jacot, F. Gabriel, and C. Hongler, “Neural tangent kernel: Convergence and generalization in neural networks,” in _Advances in Neural Information Processing Systems 31_ , 2018, pp. 8571–8580.
  * [18] D. Oglic and T. Gaertner, “Learning in reproducing kernel krein spaces,” ser. Proceedings of Machine Learning Research, vol. 80, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018, pp. 3859–3867.
  * [19] Xiaolin, Huang, A. K. Johan, Suykens, Shuning, Wang, Joachim, Hornegger, Andreas, and Maier, “Classification with truncated $\ell_{1}$ distance kernel.” _IEEE transactions on neural networks and learning systems_ , 2017\.
  * [20] J. Pennington, F. X. X. Yu, and S. Kumar, “Spherical random features for polynomial kernels,” in _Advances in Neural Information Processing Systems 28_ , 2015, pp. 1846–1854.
  * [21] P. Kar and H. Karnick, “Random feature maps for dot product kernels,” ser. Proceedings of Machine Learning Research, vol. 22, La Palma, Canary Islands, 21–23 Apr 2012, pp. 583–591.
  * [22] N. Pham and R. Pagh, “Fast and scalable polynomial kernels via explicit feature maps,” in _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’13, 2013, p. 239–247. [Online]. Available: https://doi.org/10.1145/2487575.2487591
  * [23] A. Vedaldi and A. Zisserman, “Efficient additive kernels via explicit feature maps,” vol. 34, no. 3, 2012, pp. 480–492.
  * [24] M. Kafai and K. Eshghi, “Croification: Accurate kernel classification with the efficiency of sparse linear svm,” 2017, pp. 1–1.
  * [25] F. Liu, X. Huang, L. Shi, J. Yang, and J. A. K. Suykens, “A double-variational bayesian framework in random fourier features for indefinite kernels,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 31, no. 8, pp. 2965–2979, 2020.
  * [26] F. Liu, X. Huang, Y. Chen, and J. Suykens, “Fast learning in reproducing kernel krein spaces via signed measures,” in _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , vol. 130, 13–15 Apr 2021, pp. 388–396.
  * [27] S. Bochner, “Harmonic analysis and theory of probability,” _Physics Today_ , vol. 9, 01 1956.
  * [28] M. Munkhoeva, Y. Kapushev, E. Burnaev, and I. Oseledets, “Quadrature-based features for kernel approximation,” in _Advances in Neural Information Processing Systems_ , vol. 31, 2018, pp. 9147–9156.
  * [29] K. Athreya and S. Lahiri, “Measure theory and probability theory.” Springer, 2006.
  * [30] C. S. Kubrusly, “Essentials of measure theory.” Springer, 2015.
  * [31] C. S. Ong, X. Mary, S. Canu, and A. J. Smola, “Learning with non-positive kernels,” in _Proceedings of the Twenty-First International Conference on Machine Learning_ , ser. ICML ’04, 2004, p. 81.
  * [32] R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang, and C.-J. Lin, “Liblinear: A library for large linear classification,” _J. Mach. Learn. Res._ , vol. 9, p. 1871–1874, Jun. 2008.