plain##1 ##2 ##1 ##2 (##3)

# Numerical Discretization Methods
for the Discounted Linear Quadratic Control Problem

Zhanhao Zhang, Steen Hørsholt, John Bagterp Jørgensen

###### Abstract

This study focuses on the numerical discretization methods for the continuous-
time discounted linear-quadratic optimal control problem (LQ-OCP) with time
delays. By assuming piecewise constant inputs, we formulate the discrete
system matrices of the discounted LQ-OCPs into systems of differential
equations. Subsequently, we derive the discrete-time equivalent of the
discounted LQ-OCP by solving these systems. This paper presents three
numerical methods for solving the proposed differential equations systems: the
fixed-time-step ordinary differential equation (ODE) method, the step-doubling
method, and the matrix exponential method. Our numerical experiment
demonstrates that all three methods accurately solve the differential equation
systems. Interestingly, the step-doubling method emerges as the fastest among
them while maintaining the same level of accuracy as the fixed-time-step ODE
method.

##  I Introduction

As one of the fundamental problems within optimal control theory, the linear-
quadratic optimal control problem (LQ-OCP) is attractive due to its
mathematical tractability and simplicity. This class of problems has extensive
applications in engineering, economics, and other fields. However, real-world
scenarios often present time delays, significantly influencing system
performance and robustness [9]. Integrating time delays into LQ-OCPs adds
complexity as control inputs depend on both historical and current system
states. Considering potential risks and uncertainties in the future,
incorporating discounted cost functions becomes valuable for balancing
immediate gains against long-term losses when devising control strategies [5].
The continuous nature of these problems, coupled with the complexities
introduced by time delays and discounted cost functions, can render their
implementation infeasible in real-world scenarios. Therefore, there arises a
necessity for discretization techniques to facilitate the practical
implementation of discounted LQ-OCPs with time delays.

There is rich research on the solution and discretization methods of
undiscounted optimal control problems [12, 6, 15, 4, 20]. The discretization
of continuous-time undiscounted LQ-OCPs without time delays has been explored
extensively in the literature [3, 4, 10, 20]. Furthermore, recent studies have
introduced numerical discretization methods for deterministic and stochastic
LQ-OCPs with time delays [23, 22].

On the other hand, the discounted cost approach is popular when considering
trade-offs between the present and future costs or rewards in control
strategies, e.g., reinforced learning [21] and risk-sensitive optimal control
problems [14]. The discount factor is a key component in these problems. In
[21, 7], authors have investigated the influence of various discount factors
on the stability of the discounted optimal control problems for both linear
and nonlinear systems. Besides, the discounted LQ-OCPs are related to other
control algorithms, such as model predictive control (MPC) [8] and linear-
quadratic Gaussian (LQG) [11, 16]. However, as far as we know, the existing
literature has not explored the optimal control problems incorporating with
both the discounted cost function and linear time-delay systems. Therefore, in
this paper, we would like to investigate the discretization of the continuous-
time discounted LQ-OCPs subject to time-delay systems. The key problems that
we address in this paper:

  * 1.

Formulation of differential equations systems for the discretization of
discounted LQ-OCPs with and without time delays

  * 2.

Numerical methods for solving the resulting systems of differential equations.

This paper is organized as follows. Section II introduces the formulation of
discounted LQ-OCPs with and without time delays and describes the differential
equation systems for LQ discretization. In Section III, we present three
numerical methods for solving proposed systems of differential equations.
Section IV illustrates a numerical experiment that tests and compares the
three numerical methods, and the conclusions are given in Section V.

##  II Formulation of the Problems

###  II-A Discounted linear-quadratic optimal control problem

Consider the following LQ-OCP with a discounted cost function

$\displaystyle\min_{x,u,z,\tilde{z}}\quad$
$\displaystyle\phi=\int_{t_{0}}^{t_{0}+T}l_{c}(\tilde{z}(t))dt$ (1a)
$\displaystyle s.t.$ $\displaystyle x(t_{0})=\hat{x}_{0},$ (1b) $\displaystyle
u(t)=u_{k},\quad$ $\displaystyle t_{k}\leq t<t_{k+1},\,k\in\mathcal{N},$ (1c)
$\displaystyle\dot{x}(t)=A_{c}x(t)+B_{c}u(t),\quad$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (1d) $\displaystyle z(t)=C_{c}x(t)+D_{c}u(t),\,$ $\displaystyle
t_{0}\leq t<t_{0}+T,$ (1e) $\displaystyle\bar{z}(t)=\bar{z}_{k},$
$\displaystyle t_{k}\leq t<t_{k+1},\,k\in\mathcal{N},$ (1f)
$\displaystyle\tilde{z}(t)=z(t)-\bar{z}(t),$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (1g)

with the stage cost function

$\begin{split}l_{c}(\tilde{z}(t))&=\frac{1}{2}e^{-\mu
t}\left\|W_{z}\tilde{z}(t)\right\|_{2}^{2}=\frac{1}{2}e^{-\mu
t}\left(\tilde{z}(t)^{\prime}Q_{c}\tilde{z}(t)\right),\end{split}$ (2)

where $x\in\mathbb{R}^{n_{x}\times 1}$ and $u\in\mathbb{R}^{n_{u}\times 1}$
are the state and input. $z\in\mathbb{R}^{n_{z}\times 1}$ and
$\bar{z}\in\mathbb{R}^{n_{z}\times 1}$ are the output and its reference. The
discount factor $\mu\in\mathbb{R}^{+}$ is a positive constant and the weight
matrix $Q_{c}=W_{z}^{\prime}W_{z}$ is semi-positive definite.

Note that the above problem is in continuous-time with decision variables
$x(t),u(t),z(t)$, and $\tilde{z}(t)$. The control horizon $T=NT_{s}$ with
sampling time $T_{s}$ and $N\in\mathbb{Z}^{+}$, and
$\mathcal{N}=0,1,\ldots,N-1$. We assume piecewise constant input and output
target, $u(t)=u_{k}$ and $\bar{z}(t)=\bar{z}_{k}$ for $t_{k}\leq t<t_{k+1}$.
The corresponding discrete-time equivalent of (1) is

$\displaystyle\min_{x,u}\quad$
$\displaystyle\phi=\sum_{k\in\mathcal{N}}l_{k}(x_{k},u_{k})$ (3a)
$\displaystyle s.t.$ $\displaystyle x_{0}=\hat{x}_{0},$ (3b) $\displaystyle
x_{k+1}=Ax_{k}+Bu_{k},\quad$ $\displaystyle k\in\mathcal{N},$ (3c)

with the stage costs

$l_{k}(x_{k},u_{k})=\frac{1}{2}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}^{\prime}Q_{k}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+(\overbrace{M_{k}\bar{z}_{k}}^{q_{k}})^{\prime}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+\rho_{k},\,k\in\mathcal{N},$ (4)

where the quadratic and linear term matrices are $Q_{k}=e^{-\mu t_{k}}Q$ and
$M_{k}=e^{-\mu t_{k}}M$. The constant term $\rho_{k}$ is

$\rho_{k}=\frac{1}{2}\int_{t_{k}}^{t_{k+1}}e^{-\mu
t}\bar{z}_{k}^{\prime}Q_{c}\bar{z}_{k}dt=\frac{e^{-\mu t_{k}}(1-e^{-\mu
T_{s}})}{2\mu}\bar{z}_{k}^{\prime}Q_{c}\bar{z}_{k}.$ (5)

######  Proposition 1 (Discretization of discounted LQ-OCPs).

The system of differential equations

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (6a) $\displaystyle\dot{B}(t)$ $\displaystyle=A(t)B_{c},\qquad$
$\displaystyle B(0)=0,$ (6b) $\displaystyle\dot{H}_{q}(t)$
$\displaystyle=H_{cq}H_{q}(t),$ $\displaystyle H_{q}(0)=I_{xu},$ (6c)
$\displaystyle\dot{H}_{m}(t)$ $\displaystyle=H_{cm}H_{m}(t),$ $\displaystyle
H_{m}(0)=I_{xu},$ (6d) $\displaystyle\dot{Q}(t)$
$\displaystyle=H_{q}(t)^{\prime}\bar{Q}_{c}H_{q}(t),\qquad$ $\displaystyle
Q(0)=0,$ (6e) $\displaystyle\dot{M}(t)$
$\displaystyle=H_{m}(t)^{\prime}\bar{M}_{c},\qquad$ $\displaystyle M(0)=0,$
(6f) where $I_{xu}$ is an identity matrix with size $n_{xu}=n_{x}+n_{u}$ and
$\displaystyle H_{c}=\begin{bmatrix}A_{c}&B_{c}\\\
0&0\end{bmatrix},\>\>H_{cq}=H_{c}-\frac{\mu}{2}I_{xu},\>\>H_{cm}=H_{c}-\mu
I_{xu},$ (6g)
$\displaystyle\qquad\bar{M}_{c}=-\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}^{\prime}Q_{c},\qquad\bar{Q}_{c}=-\bar{M}_{c}\begin{bmatrix}C_{c}&D_{c}\end{bmatrix},$
(6h)

may be used to compute the discrete-time system matrices ($A=A(T_{s})$,
$B=B(T_{s})$, $Q=Q(T_{s})$, $M=M(T_{s})$) of discounted LQ-OCPs without time
delays.

###### Proof.

We can compute discrete-time state space matrices by $A(t)=e^{A_{c}t}$,
$B(t)=\int_{0}^{t}A(s)B_{c}ds$. The solution of the state space model thus can
be defined as

$\Gamma(t)=\begin{bmatrix}A(t)&B(t)\\\
0&I\end{bmatrix}=e^{H_{c}t},\;\;z(t)=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}\Gamma(t)\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}.$ (7)

Moreover, as $e^{-\mu t}$ is a scalar and $e^{-\mu
t}={e^{-\frac{\mu}{2}t}}^{\prime}e^{-\frac{\mu}{2}t}$, we obtain the following
expressing by replacing $z(t)$ with (7) in $l_{c}(\tilde{z}(t))$

$\displaystyle
H_{q}(t)=\Gamma(t)e^{-\frac{\mu}{2}t}=e^{H_{cq}t},\;\;\;H_{m}(t)=\Gamma(t)e^{-\mu
t}=e^{H_{cm}t},$ (8a) $\displaystyle Q(t)=\int_{0}^{t}e^{-\mu
s}{e^{H_{c}s}}^{\prime}\bar{Q}_{c}e^{H_{c}s}ds=\int_{0}^{t}{e^{H_{cq}s}}^{\prime}\bar{Q}_{c}e^{H_{cq}s}ds,$
(8b) $\displaystyle M(t)=\int_{0}^{t}e^{-\mu
s}{e^{H_{c}s}}^{\prime}\bar{M}_{c}ds=\int_{0}^{t}{e^{H_{cm}s}}^{\prime}\bar{M}_{c}ds.$
(8c)

######  Remark 1.

Note that when the discount factor $\mu$ becomes a diagonal matrix
$\mathcal{M}=\text{diag}(\mu_{1},\mu_{2},\ldots,\mu_{n_{z}})$ for
$l_{c}(\tilde{z}(t))=0.5\tilde{z}(t)^{\prime}(e^{-\mathcal{M}t}Q_{c})\tilde{z}(t)$,
the corresponding system of differential equations (6) becomes

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (9a) $\displaystyle\dot{B}(t)$ $\displaystyle=A(t)B_{c},\qquad$
$\displaystyle B(0)=0,$ (9b) $\displaystyle\dot{H}_{q}(t)$
$\displaystyle=H_{q}(t)H_{c}-0.5\mathcal{M}H_{q}(t),\quad$ $\displaystyle
H_{q}(0)=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix},$ (9c)
$\displaystyle\dot{H}_{m}(t)$
$\displaystyle=H_{m}(t)H_{c}-\mathcal{M}H_{m}(t),\quad$ $\displaystyle
H_{m}(0)=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix},$ (9d)
$\displaystyle\dot{Q}(t)$
$\displaystyle=H_{q}(t)^{\prime}Q_{c}H_{q}(t),\qquad$ $\displaystyle Q(0)=0,$
(9e) $\displaystyle\dot{M}(t)$ $\displaystyle=-H_{m}(t)^{\prime}Q_{c},\qquad$
$\displaystyle M(0)=0,$ (9f)

where $H_{q}=e^{-0.5\mathcal{M}t}[C_{c}\;\;D_{c}]e^{H_{c}t}$ and
$H_{m}=e^{-0.5\mathcal{M}t}H_{q}$.

The differential equations described in Remark 1 allow us to solve the
discrete system matrices using standard ODE methods such as Euler and Runge-
Kutta methods. In this paper, we only consider the case $e^{-\mu t}$ as a
scalar.

###  II-B Discounted linear-quadratic optimal control problem with time delays

Consider the following discounted LQ-OCP subject to a time-delay system

$\displaystyle\min_{x,u,z,\tilde{z}}\>\>$
$\displaystyle\phi=\int_{t_{0}}^{t_{0}+T}l_{c}(\tilde{z}(t))dt$ (10a)
$\displaystyle s.t.$ $\displaystyle x(t_{0})=\hat{x}_{0},$ (10b)
$\displaystyle u(t)=u_{k},\>$ $\displaystyle t_{k}\leq
t<t_{k+1},\,k\in\mathcal{N},$ (10c)
$\displaystyle\dot{x}(t)=A_{c}x(t)+B_{c}u(t-\tau),$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (10d) $\displaystyle z(t)=C_{c}x(t)+D_{c}u(t-\tau),$
$\displaystyle t_{0}\leq t<t_{0}+T,$ (10e)
$\displaystyle\bar{z}(t)=\bar{z}_{k},$ $\displaystyle t_{k}\leq
t<t_{k+1},\,k\in\mathcal{N},$ (10f)
$\displaystyle\tilde{z}(t)=z(t)-\bar{z}(t),$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (10g)

where $\tau\in\mathbb{R}^{+}_{0}$ is the time delay. The stage cost function
$l_{c}(\tilde{z}(t))$ is identical to (2). Note that when considering the MIMO
system with multiple time delays and the system dynamic equations (10d) and
(10e) become

$\displaystyle\dot{x}_{ij}(t)$
$\displaystyle=A_{c,ij}x_{ij}(t)+B_{c,ij}u_{j}(t-\tau_{ij}),$ (11a)
$\displaystyle z_{ij}(t)$
$\displaystyle=C_{c,ij}x_{ij}(t)+D_{c,ij}u_{j}(t-\tau_{ij}).$ (11b)
$\displaystyle\quad u=[u_{1};u_{2};\ldots,u_{n_{u}}],$ $\displaystyle
x=[x_{11};x_{21};\ldots;x_{n_{z}n_{u}}],$ (11c) $\displaystyle\quad
z=[z_{1};z_{2};\ldots;z_{n_{z}}],$ $\displaystyle
z_{i}=\sum_{j=1}^{n_{u}}z_{ij},$ (11d)

where $A_{c,ij}$, $B_{c,ij}$, $C_{c,ij}$, $D_{c,ij}$ and $\tau_{ij}$ for
$i=1,2,\ldots,n_{z}$ and $j=1,2,\ldots,n_{u}$ are parameters of the $[i,j]$
SISO system describing the dynamics from the $j^{th}$ input to the $i^{th}$
output.

Based on [22], the discrete-time equivalent of (11) is

$\displaystyle\overbrace{\left[\begin{array}[]{c}x_{k+1}\\\
\hdashline[2pt/2pt]u_{o,k+1}\end{array}\right]}^{=\tilde{x}_{k+1}}$
$\displaystyle=\overbrace{\left[\begin{array}[]{c;{2pt/2pt}c}A&B_{o,1}\\\
\hdashline[2pt/2pt]0&I_A\end{array}\right]}^{=\tilde{A}}\overbrace{\left[\begin{array}[]{c}x_{k}\\\
\hdashline[2pt/2pt]u_{o,k}\end{array}\right]}^{=\tilde{x}_{k}}+\overbrace{\left[\begin{array}[]{c}B_{o,2}\\\
\hdashline[2pt/2pt]I_{B}\end{array}\right]}^{=\tilde{B}}u_{k},$ (12i)
$\displaystyle z_{k}$
$\displaystyle=\overbrace{\left[\begin{array}[]{c;{2pt/2pt}c}C_{c}&D_{o,1}\end{array}\right]}^{=\tilde{C}}\tilde{x}_{k}+\overbrace{D_{o,2}}^{=\tilde{D}}u_{k},$
(12k)

where $u_{o,k}=[u_{k-\bar{m}},u_{k-\bar{m}+1},\ldots,u_{k-1}]$ is the
historical input vector. $\bar{m}=\max\left\\{m_{ij}\right\\}$ indicates the
maximum integer time delay constant with
$\frac{\tau_{ij}}{T_{s}}=m_{ij}-v_{ij}$, where $m_{ij}\in\mathbb{Z}_{0}^{+}$
and $0\leq v_{ij}<1$ are integer and fractional time delay constants.

The discrete-time system matrices of (12) are

$\displaystyle\bar{C}_{j}=\text{diag}(C_{c,1j},C_{c,2j},\ldots,C_{c,n_{z}j}),\>\>\text{for}\>\>j=1,2,\ldots,n_{u},$
(13a)
$\displaystyle\bar{D}_{c,i}=\sum_{j=1}^{n_{u}}D_{c,ij}e_{j}E_{\bar{m}+1}^{m_{ij}},\qquad\qquad\>\,\>\>\text{for}\>\>i=1,2,\ldots,n_{z}.$
(13b) $\displaystyle
C_{c}=[\bar{C}_{1},\bar{C}_{2},\ldots,\bar{C}_{n_{u}}],\qquad\quad\;\;\>D_{o}=[\bar{D}_{1};\bar{D}_{2};\cdots;\bar{D}_{n_{z}}],$
(13c) $\displaystyle B_{o,1}=B_{o}(:,1:end-n_{u}),\qquad
B_{o,2}=B_{o}(:,\bar{m}n_{u}:end),$ (13d) $\displaystyle
D_{o,1}=D_{o}(:,1:end-n_{u}),\quad\;\;\>D_{o,2}=D_{o}(:,\bar{m}n_{u}:end),$
(13e) $\displaystyle I_{A}=\begin{bmatrix}0&I&\ldots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\ldots&I\\\
0&0&\ldots&0\end{bmatrix},\qquad\quad I_{B}=\begin{bmatrix}0\\\ \vdots\\\ 0\\\
I\end{bmatrix},$ (13f)

and the matrices $A$ and $B_{o}=B_{1}+B_{2}$ are

$\displaystyle
A_{c}=\text{diag}(A_{c,11},A_{c,21},\ldots,A_{c,n_{z}n_{u}}),\quad
A(t)=e^{A_{c}t},$ (14a) $\displaystyle
V=\text{diag}(V_{11},V_{21},\ldots,V_{n_{z}n_{u}}),\qquad\quad\;\;\;V_{ij}=Iv_{ij},$
(14b) $\displaystyle
B_{1c}=\left[B_{1c,11};B_{1c,21};\ldots;B_{1c,n_{z}n_{u}}\right],B_{1c,ij}=B_{c}e_{j}E_{\bar{m}+1}^{m_{ij}},$
(14c) $\displaystyle
B_{2c}=\left[B_{2c,11};B_{2c,21};\ldots;B_{2c,n_{z}n_{u}}\right],B_{2c,ij}=B_{c}e_{j}E_{\bar{m}+1}^{m_{ij}+1},$
(14d) $\displaystyle
B_{1}=\int_{0}^{T_{s}}e^{A_{c}t}B_{1c}dt,\;\;\;\;B_{2}=V\int_{0}^{T_{s}}e^{VA_{c}t}(B_{2c}-B_{1c})dt,$
(14e)

where $e_{j}=[0,\ldots,1,\ldots,0]$ and
$E_{\bar{m}+1}^{p}=[0,\ldots,I,\ldots,0]$ for $p=1,2,\ldots,\bar{m}+1$ are the
unit vector and matrix for selecting the $j^{th}$ historical input from the
augmented input vector $\tilde{u}_{k}=[u_{o,k};u_{k}]$ such that
$u_{j,k-(\bar{m}+1)+p}=e_{j}E^{p}_{\bar{m}+1}\tilde{u}_{k}$ .

The discrete equivalent of the discounted LQ-OCP with time delays (10) has the
same expressions as the non-delay case introduced in (3) and (4), except that
the system parameters become $\tilde{x}$, $\tilde{A}$, $\tilde{B}$,
$\tilde{C}$, $\tilde{D}$.

######  Proposition 2 (Discretization of discounted LQ-OCPs with time delays).

The system of differential equations

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (15a) $\displaystyle\dot{A}_{v}(t)$
$\displaystyle=VA_{c}A_{v}(t),\qquad$ $\displaystyle A_{v}(0)=I,$ (15b)
$\displaystyle\dot{B}_{1}(t)$ $\displaystyle=A(t)B_{1c},\qquad$ $\displaystyle
B_{1}(0)=0,$ (15c) $\displaystyle\dot{B}_{2}(t)$
$\displaystyle=A_{v}(t)\bar{B}_{2c},\qquad$ $\displaystyle B_{2}(0)=0,$ (15d)
$\displaystyle\dot{H}_{q}(t)$ $\displaystyle=H_{cq}H_{q}(t),$ $\displaystyle
H_{q}(0)=I_{h},$ (15e) $\displaystyle\dot{H}_{m}(t)$
$\displaystyle=H_{cm}H_{m}(t),$ $\displaystyle H_{m}(0)=I_{h},$ (15f)
$\displaystyle\dot{Q}(t)$
$\displaystyle=\Gamma_{q}(t)^{\prime}\bar{Q}_{c}\Gamma_{q}(t),\qquad$
$\displaystyle Q(0)=0,$ (15g) $\displaystyle\dot{M}(t)$
$\displaystyle=\Gamma_{m}(t)^{\prime}\bar{M}_{c},\qquad$ $\displaystyle
M(0)=0,$ (15h)

where $I_{h}=\text{diag}(I_{xu},I_{xu},I_{xu})$ is an identity matrix and

$\displaystyle H_{cq}$
$\displaystyle=H_{c}-\frac{\mathcal{\mu}}{2}I_{h},\,H_{cm}=H_{c}-\mathcal{\mu}I_{h},\,\bar{B}_{2c}=V(B_{2c}-B_{1c}),$
(16a) $\displaystyle E_{1}=\left[I_{xu},I_{xu},-I_{xu}\right],\qquad
E_{2}=\left[I_{xu};I_{xu};I_{xu}\right],$ (16b)
$\displaystyle\Gamma_{m}(t)=E_{1}H_{m}(t)E_{2},\qquad\bar{M}_{c}=-\begin{bmatrix}C_{c}&D_{o}\end{bmatrix}Q_{c},$
(16c)
$\displaystyle\Gamma_{q}(t)=E_{1}H_{q}(t)E_{2},\qquad\;\;\,\bar{Q}_{c}=-\bar{M}_{c}\begin{bmatrix}C_{c}&D_{o}\end{bmatrix},$
(16d)

may be used to compute ($A=A(T_{s})$, $B_{o}=B_{1}(T_{s})+B_{2}(T_{s})$,
$Q=Q(T_{s})$, $M=M(T_{s})$) of discounted LQ-OCPs with time delays.

###### Proof.

In [22], the solution of the time-delay state space (11) are defined as

$z(t)=\begin{bmatrix}C_{c}&D_{o}\end{bmatrix}\Gamma(t)\begin{bmatrix}x_{k}\\\
\tilde{u}_{k}\end{bmatrix},\quad\Gamma(t)=\begin{bmatrix}A(t)&B_{o}(t)\\\
0&I\end{bmatrix},$ (17)

where $A(t)=e^{A_{c}t}$ and $B_{o}(t)=B_{1}(t)+B_{2}(t)$ are described in
(14). In this case, we cannot directly express $\Gamma(t)$ as a matrix
exponential like the non-delay case introduced in Proposition 1 since $B_{o}$
consists of $B_{1}(t)$ and $B_{2}(t)$. However, we can decompose it into the
linear combination of $A$, $A_{v}$, $B_{1}$ and $B_{2}$ as

$\begin{split}\Gamma(t)&=\overbrace{\begin{bmatrix}A(t)&B_{1}(t)\\\
0&I\end{bmatrix}}^{H_{1}(t)}+\overbrace{\begin{bmatrix}A_{v}(t)&B_{2}(t)\\\
0&I\end{bmatrix}}^{H_{2}(t)}-\overbrace{\begin{bmatrix}A_{v}(t)&0\\\
0&I\end{bmatrix}}^{H_{3}(t)}\\\ &=e^{H_{1c}t}+e^{H_{2c}t}-e^{H_{3c}t}\\\
&=E_{1}H(t)E_{2},\end{split}$ (18)

with $H_{c}=\text{diag}\left(H_{1c},H_{2c},H_{3c}\right)$ and

$H_{1c}=\begin{bmatrix}A_{c}&B_{1c}\\\
0&0\end{bmatrix},H_{2c}=\begin{bmatrix}VA_{c}&B_{2c}\\\
0&0\end{bmatrix},H_{3c}=\begin{bmatrix}VA_{c}&0\\\ 0&0\end{bmatrix},$ (19)

where $H(t)=e^{H_{c}t}=\text{diag}\left(H_{1}(t),H_{2}(t),H_{3}(t)\right)$ and
$H_{k}(t)=e^{H_{kc}t}$ for $k=1,2,3$.

Consequently, as $e^{-\mu t}={e^{-\frac{\mu}{2}}}^{\prime}e^{-\frac{\mu}{2}}$
is a scalar, we will get

$\displaystyle
H_{q}(t)=e^{H_{cq}t},\qquad\;\;\Gamma_{q}(t)=\Gamma(t)e^{-\frac{\mathcal{\mu}}{2}t}=E_{1}H_{q}(t)E_{2},$
(20a) $\displaystyle
H_{m}(t)=e^{H_{cm}t},\qquad\Gamma_{m}(t)=\Gamma(t)e^{-\mathcal{\mu}t}=E_{1}H_{m}(t)E_{2},$
(20b) $\displaystyle Q(t)=\int_{0}^{t}e^{-\mu
s}\Gamma(s)^{\prime}\bar{Q}_{c}\Gamma(s)ds=\int_{0}^{t}\Gamma_{q}(s)^{\prime}\bar{Q}_{c}\Gamma_{q}(s)ds,$
(20c) $\displaystyle M(t)=\int_{0}^{t}e^{-\mu
s}\Gamma(s)^{\prime}\bar{M}_{c}ds=\int_{0}^{t}\Gamma_{m}(s)^{\prime}\bar{M}_{c}ds.$
(20d)

######  Remark 2.

Note that the differential equations introduced in (6) and (15) are identical
when time delays are $\left\\{\tau_{ij}\right\\}=0$, and the time delay
constants become $\bar{m}=0$ and $V=0$. This corresponding to $B_{1c}=B_{c}$,
$\bar{B}_{2c}=0$ and $H_{2c}=H_{3c}=0$. The expressions introduced in
Proposition 2 become

$\displaystyle A_{v}(t)=I,$ $\displaystyle B_{o}(t)=\int_{0}^{t}A(s)B_{c}ds,$
(21a) $\displaystyle H(t)=e^{H_{c}t},$
$\displaystyle\Gamma(t)=E_{1}H(t)E_{2}=e^{H_{1c}t},$ (21b) $\displaystyle
H_{q}(t)=H(t)e^{-\frac{\mu}{2}t},$
$\displaystyle\Gamma_{q}(t)=E_{1}H_{q}(t)E_{2}=e^{(H_{1c}-\frac{\mu}{2})t},$
(21c) $\displaystyle H_{m}(t)=H(t)e^{-\mu t},$
$\displaystyle\Gamma_{m}(t)=E_{1}H_{m}(t)E_{2}=e^{(H_{1c}-\mu)t},$ (21d) where
$B_{o}$, $\Gamma_{q}$ and $\Gamma_{m}$ have the same expressions as the non-
delay discrete system matrices $B$, $H_{q}$ and $H_{m}$ introduced in (6b) and
(8a).

##  III Numerical Methods of LQ Discretization

This section introduces numerical methods for solving proposed systems of
differential equations.

###  III-A Fixed-time-step ordinary differential equation method

Consider an s-stage ODE method with the number of integration steps
$N\in\mathbb{Z}^{+}$ and the time step $\delta t=\frac{T_{s}}{N}$. Note that
$a_{i,j}$ and $b_{i}$ for $i=1,2,\ldots,s$ and $j=1,2,\ldots,s$ are the
Butcher tableau’s parameters of the ODE method. We have

$\displaystyle A_{k+1}=\Lambda A_{k},$ $\displaystyle k\in\mathcal{N},$ (22a)
$\displaystyle B_{1,k+1}=B_{1,k}+\Theta_{1}A_{k}\tilde{B}_{1c},$
$\displaystyle k\in\mathcal{N},$ (22b) $\displaystyle
A_{v,k+1}=\Lambda_{v}A_{v,k},$ $\displaystyle k\in\mathcal{N},$ (22c)
$\displaystyle B_{2,k+1}=B_{2,k}+\Theta_{2}A_{v,k}\tilde{B}_{2c},$
$\displaystyle k\in\mathcal{N},$ (22d) $\displaystyle
H_{q,k+1}=\Omega_{q}H_{q,k},$ $\displaystyle k\in\mathcal{N},$ (22e)
$\displaystyle H_{m,k+1}=\Omega_{m}H_{m,k},$ $\displaystyle k\in\mathcal{N},$
(22f) $\displaystyle
M_{k+1}=M_{k}+E_{2}^{\prime}H_{m,k}^{\prime}\tilde{M}_{c},$ $\displaystyle
k\in\mathcal{N},$ (22g) $\displaystyle
Q_{k+1}=Q_{k}+E_{2}^{\prime}H_{q,k}^{\prime}\tilde{Q}_{c}H_{q,k}E_{2},$
$\displaystyle k\in\mathcal{N},$ (22h) where
$\displaystyle\tilde{B}_{1c}=\delta tB_{1c},$
$\displaystyle\tilde{M}_{c}=\delta
t\sum_{i=1}^{s}b_{i}\Omega_{m,i}^{\prime}E_{1}^{\prime}\bar{M}_{c},$ (22i)
$\displaystyle\tilde{B}_{2c}=\delta t\bar{B}_{2c},$
$\displaystyle\tilde{Q}_{c}=\delta
t\sum_{i=1}^{s}b_{i}\Omega_{q,i}^{\prime}E_{1}^{\prime}\bar{Q}_{c}E_{1}\Omega_{q,i}.$
(22j)

The coefficients $\Lambda$, $\Lambda_{v}$, $\Theta_{1}$, $\Theta_{2}$,
$\Omega_{m}$ and $\Omega_{q}$ are functions of the Butcher tableau’s
parameters. They can be computed as

$\displaystyle A_{k,i}=A_{k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{A}_{k,j}=\Lambda_{i}A_{k},$ (23a) $\displaystyle
B_{1,k,i}=B_{1,k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{B}_{1,k,j}=B_{1,k}+\Theta_{1,i}A_{k}\tilde{B}_{1c},$
(23b) $\displaystyle A_{v,k,i}=A_{v,k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{A}_{v,k,j}=\Lambda_{v,i}A_{v,k},$ (23c)
$\displaystyle B_{2,k,i}=B_{2,k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{B}_{2,k,j}=B_{2,k}+\Theta_{2,i}A_{v,k}\tilde{B}_{2c},$
(23d) $\displaystyle H_{m,k,i}=H_{m,k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{H}_{m,k,j}=\Omega_{m,i}H_{m,k},$ (23e)
$\displaystyle H_{q,k,i}=H_{q,k}+\delta
t\sum_{j=1}^{s}a_{i,j}\dot{H}_{q,k,j}=\Omega_{q,i}H_{q,k},$ (23f)

Algorithm 1 Fixed-time-step ODE method

Input: $(A_{c},B_{c},C_{c},D_{c},\tau,Q_{c},T_{s},N)$

Output: $(A,B_{o},Q,M)$

Compute the step size $\delta t=\frac{T_{s}}{N}$

Compute system matrices
($B_{1c}$,$B_{2c}$,$\bar{M}_{c}$,$\bar{Q}_{c}$,$H_{cq}$,$H_{cm}$) using (16)

Set initial states ($k=0$, $A_{k}=I$, $A_{v,k}=I$, $B_{1,k}=0$, $B_{2,k}=0$,
$H_{m,k}=I_{h}$, $H_{q,k}=I_{h}$, $Q_{k}=0$, $M_{k}=0$)

Compute stage variable coefficients ($\Lambda_{i}$, $\Lambda_{v,i}$,
$\Theta_{1,i}$, $\Theta_{2,i}$, $\Omega_{m,i}$ and $\Omega_{q,i}$) using (23)

Compute constant coefficients ($\Lambda$, $\Lambda_{v}$, $\Theta_{1}$,
$\Theta_{2}$, $\Omega_{m}$ and $\Omega_{q}$) using (24)

while $k<N$ do

Use (22) to update
$(A_{k},B_{1,k},A_{v,k},B_{2,k},H_{m,k},H_{q,k},Q_{k},M_{k})$

Set $k\leftarrow k+1$

end while

Set ($A=A_{k}$, $B_{o}=B_{1,k}+B_{2,k}$, $Q=Q_{k}$, $M=M_{k}$)

and

$\displaystyle\Lambda=I+\delta t\sum_{i=1}^{s}b_{i}A_{c}\Lambda_{i},\quad$
$\displaystyle\Lambda_{v}=I+\delta t\sum_{i=1}^{s}b_{i}VA_{c}\Lambda_{v,i},$
(24a) $\displaystyle\Theta_{1}=\sum_{i=1}^{s}b_{i}\Lambda_{i},$
$\displaystyle\Omega_{m}=I_{h}+\delta t\sum_{i=1}^{s}b_{i}H_{cm}\Omega_{m,i},$
(24b) $\displaystyle\Theta_{2}=\sum_{i=1}^{s}b_{i}\Lambda_{v,i},$
$\displaystyle\Omega_{q}=I_{h}+\delta t\sum_{i=1}^{s}b_{i}H_{cq}\Omega_{q,i},$
(24c)

where $\Lambda_{i}$, $\Lambda_{v,i}$, $\Theta_{1,i}$, $\Theta_{2,i}$,
$\Omega_{m,i}$ and $\Omega_{q,i}$ are coefficients of stage variables
$A_{k,i}$, $A_{v,k,i}$, $B_{1,k,i}$, $B_{2,k,i}$, $H_{m,k,i}$ and $H_{q,k,i}$.

Consequently, we can obtain $A(T_{s})=A_{N}$, $B_{o}(T_{s})=B_{1,N}+B_{2,N}$,
$Q(T_{s})=Q_{N}$, $M(T_{s})=M_{N}$ with constant coefficients $\Lambda$,
$\Lambda_{v}$, $\Theta_{1}$, $\Theta_{2}$, $\Omega_{m}$ and $\Omega_{q}$ when
using fixed-time-step ODE methods. Algorithm 1 presents the fixed-time-step
ODE method for the discretization of the discounted LQ-OCPs with time delays.

###  III-B Step-doubling method

Table I: Numerical expressions of the step-doubling method |  Numerical expression |  Step-doubling function
---|---|---
$\tilde{A}(N)$ |  $\bar{\Lambda}^{N}$ |  $\tilde{A}(\frac{N}{2})\tilde{A}(\frac{N}{2})$
$\tilde{B}_{o}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\bar{\Lambda}^{i}$ |  $\tilde{B}_{o}(\frac{N}{2})\left(I+\tilde{A}(\frac{N}{2})\right)$
$\tilde{H}_{m}(N)$ |  $\Omega_{m}^{N}$ |  $\tilde{H}_{m}(\frac{N}{2})\tilde{H}_{m}(\frac{N}{2})$
$\tilde{H}_{q}(N)$ |  $\Omega_{q}^{N}$ |  $\tilde{H}_{q}(\frac{N}{2})\tilde{H}_{q}(\frac{N}{2})$
$\tilde{M}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\left(\Omega_{m}^{i}\right)^{\prime}$ |  $\tilde{M}(\frac{N}{2})\left(I_{h}+\tilde{H}_{m}(\frac{N}{2})^{\prime}\right)$
$\tilde{Q}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\left(\Omega_{q}^{i}\right)^{\prime}\tilde{Q}_{c}\left(\Omega_{q}^{i}\right)$ |  $\tilde{Q}(\frac{N}{2})+\tilde{H}_{q}(\frac{N}{2})^{\prime}\tilde{Q}(\frac{N}{2})\tilde{H}_{q}(\frac{N}{2})$
Algorithm 2 Step-doubling method

Input: $(A_{c},B_{c},C_{c},D_{c},\tau,Q_{c},T_{s},j)$

Output: $(A,B_{o},Q,M)$

Compute the number of the integration step $N=2^{j}$

Compute the step size $\delta t=\frac{T_{s}}{N}$

Compute system matrices
($B_{1c}$,$B_{2c}$,$\bar{M}_{c}$,$\bar{Q}_{c}$,$H_{cq}$,$H_{cm}$) using (16)

Compute stage variable coefficients ($\Lambda_{i}$, $\Lambda_{v,i}$,
$\Theta_{1,i}$, $\Theta_{2,i}$, $\Omega_{m,i}$ and $\Omega_{q,i}$) using (23)

Compute constant coefficients ($\Lambda$, $\Lambda_{v}$, $\Theta_{1}$,
$\Theta_{2}$, $\Omega_{m}$ and $\Omega_{q}$) using (24)

Set initial states ($i=1,\tilde{A}(i)=\bar{\Lambda}$, $\tilde{B}_{o}(i)=I$,
$\tilde{H}_{m}(i)=\Omega_{m}$, $\tilde{H}_{q}(i)=\Omega_{q}$,
$\tilde{M}(i)=I_{xu}$, $\tilde{Q}(i)=\tilde{Q}_{c}$)

while $i\leq j$ do

Update ($\tilde{M}(i)$, $\tilde{Q}(i)$) using equations from Table I

Update ($\tilde{A}(i)$, $\tilde{B}_{o}(i)$,$\tilde{H}_{m}(i)$,
$\tilde{H}_{q}(i)$) using equations from Table I

Set $i\leftarrow i+1$

end while

Use (26) to compute $(A,B_{o},Q,M)$

Consider the fixed-time-step ODE method with the integration step $N=2^{j}$
for $j\in\mathbb{Z}^{+}$ and the step size $\delta t=\frac{T_{s}}{N}$, the
matrices

$\displaystyle\tilde{A}(N)=\bar{\Lambda}^{N},$
$\displaystyle\tilde{A}(1)=\bar{\Lambda},$ (25a)
$\displaystyle\tilde{B}_{o}(N)=\displaystyle\sum_{i=0}^{N-1}\bar{\Lambda}^{i},$
$\displaystyle\tilde{B}(1)=I,$ (25b)
$\displaystyle\tilde{H}_{m}(N)=\Omega_{m}^{N},$
$\displaystyle\tilde{H}_{m}(1)=\Omega_{m},$ (25c)
$\displaystyle\tilde{H}_{q}(N)=\Omega_{q}^{N},$
$\displaystyle\tilde{H}_{q}(1)=\Omega_{q},$ (25d)
$\displaystyle\tilde{M}(N)=\displaystyle\sum_{i=0}^{N-1}\left(\Omega_{m}^{i}\right)^{\prime},$
$\displaystyle\tilde{M}(1)=I_{xu},$ (25e)
$\displaystyle\tilde{Q}(N)=\displaystyle\sum_{i=0}^{N-1}\left(\Omega_{q}^{i}\right)^{\prime}\tilde{Q}_{c}\left(\Omega_{q}^{i}\right),$
$\displaystyle\tilde{Q}(1)=\tilde{Q}_{c},$ (25f)

can be used to compute $(A,B_{o},M,Q)$

$\displaystyle A(T_{s})=\tilde{A}(N)(1:n_{x},1:n_{x}),$ $\displaystyle
B_{o}(T_{s})=\Theta_{o}\tilde{B}_{o}(N)\tilde{B}_{oc},$ (26a) $\displaystyle
M(T_{s})=E_{2}^{\prime}\tilde{M}(N)\tilde{M}_{c},$ $\displaystyle
Q(T_{s})=E_{2}^{\prime}\tilde{Q}(N)E_{2},$ (26b)

where $\bar{\Lambda}=\text{diag}(\Lambda,\Lambda_{v})$,
$\Theta_{o}=[\Theta_{1},\Theta_{2}]$,
$\tilde{B}_{oc}=[\tilde{B}_{1c};\tilde{B}_{2c}]$ are constant. The
coefficients $\Lambda$, $\Theta_{1}$, $\Lambda_{v}$, $\Theta_{2}$,
$\Omega_{m}$, $\Omega_{q}$, $\tilde{M}_{c}$ and $\tilde{Q}_{c}$ are introduced
in (22) and (24).

In [1, 2, 13], a scaling and squaring algorithm is introduced for solving the
matrix exponential problem. For $A(t)=e^{At}$, they use the $\frac{n}{2}^{th}$
step’s result $A(\frac{n}{2}\delta t)$ as the initial state to compute the
double step’s result $A(n\delta t)=A(\frac{n}{2}\delta t)A(\frac{n}{2}\delta
t)$ and repeat it until $n\delta t=T_{s}$. We can use the same idea to compute
matrices introduced in (25), and it leads to the step-doubling method. Define
$f(n)$ for
$f\in\left[\tilde{A},\tilde{B}_{o},\tilde{H}_{m},\tilde{H}_{q},\tilde{M},\tilde{Q}\right]$,
the step-doubling expression for computing $f(n)$ can be written as

$f(1)\rightarrow f(2)\rightarrow f(4)\rightarrow\ldots\rightarrow
f(\frac{N}{4})\rightarrow f(\frac{N}{2})\rightarrow f(N),$ (27a) where
$f(n)=f(\frac{n}{2})f(\frac{n}{2}),\qquad\qquad
n\in[2,4,\ldots,\frac{N}{2},N].$ (27b)

Table I describes the step-doubling expressions for ($\tilde{A}$,
$\tilde{B}_{o}$, $\tilde{H}_{m}$, $\tilde{H}_{m}$, $\tilde{M}$, $\tilde{Q}$).
The step-doubling method takes only $j$ steps to get the same result as the
fixed-time-step ODE method with $N=2^{j}$ integration steps. Algorithm 2
describes the step-doubling method for the discretization of the discounted
LQ-OCPs.

###  III-C Matrix exponential method

Based on the formulas described in [17, 18, 19], we can obtain ($A$, $B_{o}$,
$M$, $Q$) by solving the following matrix exponential problems

$\displaystyle\begin{bmatrix}\Phi_{1,11}&\Phi_{1,12}\\\
0&\Phi_{1,22}\end{bmatrix}$
$\displaystyle=\text{exp}\left(\begin{bmatrix}-H_{cq}^{\prime}&E_{1}^{\prime}\bar{Q}_{c}E_{1}\\\
0&H_{cq}\end{bmatrix}t\right),$ (28a)
$\displaystyle\begin{bmatrix}I&\Phi_{2,12}\\\ 0&\Phi_{2,22}\end{bmatrix}$
$\displaystyle=\text{exp}\left(\begin{bmatrix}0&I\\\
0&H_{cm}^{\prime}\end{bmatrix}t\right),$ (28b) $\displaystyle\Phi_{3}$
$\displaystyle=\text{exp}\left(H_{c}t\right),$ (28c) and the elements are
$\displaystyle\Phi_{1,12}=H_{q}(-t)^{\prime}\int_{0}^{t}H_{q}(s)^{\prime}E_{1}^{\prime}\bar{Q}_{c}E_{1}H_{q}(s)ds,$
(28d) $\displaystyle\Phi_{1,22}=H_{q}(t)=e^{H_{cq}t},$ (28e)
$\displaystyle\Phi_{2,12}=\int_{0}^{t}H_{m}(s)ds,$ (28f)
$\displaystyle\Phi_{3}=H(t)=\text{diag}(H_{1}(t),H_{2}(t),H_{3}(t)),$ (28g)

where the matrices $H_{c}$, $H_{cq}$, $H_{cm}$, $\bar{Q}_{c}$, $H_{q}$,
$H_{m}$, $H_{1}$, $H_{2}$ and $H_{3}$ are introduced in Proposition 2.

Set $t=T_{s}$, we can compute ($A$, $B_{o}$, $M$, $Q$) as

$\displaystyle A=\Gamma(1:n_{x},1:n_{x}),$ $\displaystyle
B_{o}=\Gamma(1:n_{x},n_{x}+1:end),$ (29a) $\displaystyle
M=E_{2}^{\prime}\Phi_{2,12}\bar{M}_{c},$ $\displaystyle
Q=E_{2}^{\prime}\Phi_{1,22}^{\prime}\Phi_{1,12}E_{2},$ (29b) where
$\Gamma=E_{1}\Phi_{3}E_{2}$.

##  IV Numerical Experiments

In this section, we will test and compare the proposed three numerical
methods.

Consider a MIMO input-output model $Z(s)=G(s)U(s)$ with the following transfer
functions

$G(s)=\begin{bmatrix}\frac{1}{(1.5s+1)(3s+1)}e^{-0.1s}&\frac{-2(2s+1)}{3.4s+1}e^{-1.6s}\\\
\frac{-0.5}{2.3s+1}e^{-2.0s}&\frac{2.4}{(1.7s+1)(0.9s+1)}e^{-0.9s}\end{bmatrix}.$
(30)

We can convert the above transfer functions into the state space models
introduced in (11). The state space matrices are in observable canonical form.

Consider the discounted LQ-OCP described in (10). The weight matrix
$Q_{c}=\text{diag}(1.0,2.0)$, the discount factor $\mu=0.2$, the control
horizon $T=20$ [s] and the sampling time $T_{s}=1$ [s]. We use the symbolic
toolbox in Matlab to compute the discrete system matrices ($A$, $B_{o}$, $Q$,
$M$) with their analytic expressions described in Proposition 2. The results
are used as the true solution for comparing the results from proposed
numerical methods.

Fig. 1 describes the error and CPU time of the fixed-time-step ODE and the
step-doubling method. We notice that the step-doubling method (dot plots) has
the same error as the fixed-time-step ODE method (line plots) with the same
discretization scheme. All tested methods have the correct convergence order
(indicated by dashed lines). The CPU time of the two methods is indicated in
the bar plots. We observe that the CPU time of the fixed-time-step ODE method
increases as the integration step or the stage number increases. The step-
doubling method’s CPU time is stable at around 1.8 [ms].

Figure 1: The error and CPU time of the fixed-time-step ODE method and the step-doubling method with different discretization schemes. The error is $e(i)=\left\|i(T_{s})-i(N)\right\|_{\infty}$ for $i\in[A,B_{o},M,Q]$, where $i(T_{s})$ is the result from the Matlab’s symbolic toolbox. Table II: CPU time and error of the scenario using classic RK4 with $N=2^{10}$ | Unit | Matrix Exp. | ODE Method | Step-doubling
---|---|---|---|---
$e(A)$ | [-] | $3.34\cdot 10^{-16}$ | $9.12\cdot 10^{-14}$ | $9.12\cdot 10^{-14}$
$e(B_{o})$ | [-] | $5.56\cdot 10^{-17}$ | $8.33\cdot 10^{-12}$ | $8.33\cdot 10^{-12}$
$e(M)$ | [-] | $5.10\cdot 10^{-16}$ | $6.88\cdot 10^{-14}$ | $6.88\cdot 10^{-14}$
$e(Q)$ | [-] | $8.10\cdot 10^{-16}$ | $2.55\cdot 10^{-13}$ | $2.55\cdot 10^{-13}$
CPU Time | [ms] | 44.1 | 60.3 | 1.96

Table II describes the error and CPU time of all three methods when using the
classic RK4 with $j=10$ and $N=2^{j}$ (not for the matrix exponential). The
matrix exponential method obtains the most precise result among all methods,
while the other two methods have the same error. The fixed-time-step ODE
method is the slowest, followed by the matrix exponential method. They spend
60.3 [ms] and 44.1 [ms], respectively. The step-doubling method is the fastest
and only takes 1.96 [ms].

##  V Conclusions

This paper discussed the discretization of discounted LQ-OCPs with and without
time delays. In the Propositions, the discrete system matrices of the
discounted LQ-OCPs are described as systems of differential equations. Then,
we introduced three numerical methods for solving the proposed systems of
differential equations. All three numerical methods are tested and compared in
the numerical experiment. Our results indicate that

  * 1.

All three methods can solve the proposed differential equations, and the
matrix exponential method is the most precise.

  * 2.

The step-doubling method is significantly faster than the other two methods,
keeping the same accuracy level as the fixed-time-step ODE method.

## References

  * [1] A. H. Al-Mohy and N. J. Higham, “A new scaling and squaring algorithm for the matrix exponential,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 31, no. 3, pp. 970–989, 2010.
  * [2] ——, “Computing the action of the matrix exponential, with an application to exponential integrators,” _SIAM Journal on Scientific Computing_ , vol. 33, no. 2, pp. 488–511, 2011.
  * [3] J. Cullum, “Discrete approximations to continuous optimal control problems,” _SIAM Journal on Control_ , vol. 7, no. 1, pp. 32–49, 1969.
  * [4] A. L. Dontchev and W. W. Hager, “The euler approximation in state constrained optimal control,” _Mathematics of Computation_ , vol. 70, no. 233, pp. 173–203, 2001.
  * [5] M. Fu _et al._ , “Risk-sensitive reinforcement learning via policy gradient search,” _arXiv preprint arXiv:1810.09126_ , 2018.
  * [6] R. Goebel and M. Subbotin, “Continuous time linear quadratic regulator with control constraints via convex duality,” _IEEE Transactions on Automatic Control_ , vol. 52, no. 5, pp. 886–892, 2007.
  * [7] M. Granzotto, R. Postoyan, L. Buşoniu, D. Nešić, and J. Daafouz, “Finite-horizon discounted optimal control: stability and performance,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 2, pp. 550–565, 2020.
  * [8] L. Grüne, W. Semmler, and M. Stieler, “Using nonlinear model predictive control for dynamic decision problems in economics,” _Journal of Economic Dynamics and Control_ , vol. 60, pp. 112–133, 2015.
  * [9] K. Gu, V. L. Kharitonov, and J. Chen, _Stability of Time-Delay Systems_. Birkhäuser, 2003.
  * [10] L. Han, M. Camlibel, J. Pang, and W. Heemels, “Convergence of discrete-time approximations of constrained linear-quadratic optimal control problems,” in _49th IEEE Conference on Decision and Control (CDC)_ , 2010, pp. 5210–5215.
  * [11] L. Hansen and T. Sargent, “Discounted linear exponential quadratic gaussian control,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 5, pp. 968–971, 1995.
  * [12] E. Hendricks, P. H. Sørensen, and O. Jannerup, _Linear Systems Control: Deterministic and Stochastic Methods_. Berlin, German: Springer, 2008.
  * [13] N. J. Higham, “The scaling and squaring method for the matrix exponential revisited,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 26, no. 4, pp. 1179–1193, 2005.
  * [14] D. Jacobson, “Optimal stochastic linear systems with exponential performance criteria and their relation to deterministic differential games,” _IEEE Transactions on Automatic control_ , vol. 18, no. 2, pp. 124–131, 1973.
  * [15] J. B. Jørgensen, G. Frison, N. F. Gade-Nielsen, and B. Damman, “Numerical methods for solution of the extended linear quadratic control problem,” _IFAC Proceedings Volumes_ , vol. 45, no. 17, pp. 187–193, 2012, 4th IFAC Conference on Nonlinear Model Predictive Control.
  * [16] H. Mena, L.-M. Pfurtscheller, and M. Voigt, “Discounted cost linear quadratic gaussian control for descriptor systems,” _International Journal of Control_ , vol. 95, no. 5, pp. 1349–1362, 2022.
  * [17] C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix,” _SIAM Rev._ , vol. 20, no. 4, p. 801–836, oct 1978.
  * [18] ——, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” _SIAM Review_ , vol. 45, no. 1, pp. 3–49, 2003.
  * [19] C. Van Loan, “Computing integrals involving the matrix exponential,” _IEEE Transactions on Automatic Control_ , vol. 23, no. 3, pp. 395–404, 1978.
  * [20] F. L. Walter Alt, Robert Baier and M. Gerdts, “Approximations of linear control problems with bang-bang solutions,” _Optimization_ , vol. 62, no. 1, pp. 9–32, 2013.
  * [21] D. Wang, J. Ren, M. Ha, and J. Qiao, “System stability of learning-based linear optimal control with general discounted value iteration,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2022.
  * [22] Z. Zhang, S. Hørsholt, and J. B. Jørgensen, “Numerical discretization methods for linear quadratic control problems with time delays,” in _12th IFAC Symposium on Advanced Control of Chemical Proccesses (ADCHEM 2024) (accepted)_ , 2023.
  * [23] Z. Zhang, J. L. Svensen, M. W. Kaysfeld, A. H. D. Andersen, S. Hørsholt, and J. B. Jørgensen, “Numerical discretization methods for the extended linear quadratic control problem,” in _European Control Conference (ECC) 2024 (accepted)_ , 2023.


# Modeling Uncertainty by Learning a Hierarchy of Deep Neural Connections

Raanan Y. Rohekar
Intel AI Lab
<EMAIL_ADDRESS>
&Yaniv Gurwicz
Intel AI Lab
<EMAIL_ADDRESS>
&Shami Nisimov
Intel AI Lab
<EMAIL_ADDRESS>
&Gal Novik
Intel AI Lab
<EMAIL_ADDRESS>
All authors contributed equally.

###### Abstract

Modeling uncertainty in deep neural networks, despite recent important
advances, is still an open problem. Bayesian neural networks are a powerful
solution, where the prior over network weights is a design choice, often a
normal distribution or other distribution encouraging sparsity. However, this
prior is agnostic to the generative process of the input data, which might
lead to unwarranted generalization for out-of-distribution tested data. We
suggest the presence of a confounder for the relation between the input data
and the discriminative function given the target label. We propose an approach
for modeling this confounder by sharing neural connectivity patterns between
the generative and discriminative networks. This approach leads to a new deep
architecture, where networks are sampled from the posterior of local causal
structures, and coupled into a compact hierarchy. We demonstrate that sampling
networks from this hierarchy, proportionally to their posterior, is efficient
and enables estimating various types of uncertainties. Empirical evaluations
of our method demonstrate significant improvement compared to state-of-the-art
calibration and out-of-distribution detection methods.

##  1 Introduction

Deep neural networks have become an important tool in applied machine
learning, achieving state-of-the-art regression and classification accuracy in
many domains. However, quantifying and measuring uncertainty in these
discriminative models, despite recent important advances, is still an open
problem. Representation of uncertainty is crucial for many domains, such as
safety-critical applications, personalized medicine, and recommendation
systems [7]. Common deep neural networks are not designed to capture model
uncertainty, hence estimating it implicitly from the prediction is often
inaccurate. Several types of uncertainties are commonly discussed [6, 20],
where the two main types are: 1) epistemic uncertainty and 2) aleatoric
uncertainty. Epistemic uncertainty is caused by the lack of knowledge,
typically in cases where only a small training data exists, or for out-of-
distribution inputs. Aleatoric uncertainty is caused by noisy data, and
contrary to epistemic uncertainty, does not vanish in the large sample limit.
In this work, we focus on two aspects of uncertainty, often required
addressing in practical uses of neural networks: _calibration_ [2] and _out-
of-distribution_ (OOD) detection.

Calibration is a notion that describes the relation between a predicted
probability of an event and the actual proportion of occurrences of that
event. It is generally measured using (strictly) _proper scoring rules_ [8],
such as negative log-likelihood (NLL) and the Brier score, both are minimized
for calibrated models. Guo et al. [9] examined the calibration of recent deep
architectures by computing the expected calibration error (ECE)—the difference
between an approximation of the empirical reliability curve and the optimal
reliability curve [22]. Miscalibration is often addressed by post processing
the outputs [25, 9], approximating the posterior over the weights of a pre-
trained network [26], or by learning an ensemble of networks [14, 19, 6, 17].

OOD detection is often addressed by designing a loss function to optimize
during parameter learning of a given network architecture [20, 3]. Many of
these methods are specifically tailored for detecting OOD, requiring some
information about OOD samples, which is often impractical to obtain in real-
world cases. In addition, these methods often are not capable of modeling
different types of uncertainty. Ensemble methods, on the other hand, learn
multiple sets of network parameters for a given structure [17], or approximate
the posterior distribution over the weights from which they sample at test
time [1, 19, 6].

In this paper, we make the distinction between _structure-based_ methods,
which include ensemble methods that replicate the structure or sample subsets
from it, and _parameter-based_ methods, which specify a loss function to be
used for a given structure. Ensemble methods, in general, do not specify the
loss function to be used for parameter learning, nor do they restrict post-
processing of their output. It is interesting to note that while the majority
of ensemble methods use distinct sets of parameters for each network [17], in
the MC-dropout method [6] the parameters are shared across multiple networks.
Common to all these methods is that they use a single network architecture
(structure), as it is generally unclear how to fuse outputs from different
structures.

We propose a method that samples network structures, where parts of one
sampled structure may be similar to parts of another sampled structure but
having different weights values. In addition, these structures may share some
parts, along with their weights, with other structures (weight sharing),
specifically in the deeper layers. All these properties are learned from the
input data.

##  2 Background

We focus on two approaches that are commonly used for modeling uncertainty: 1)
Bayesian neural networks [6, 1], and 2) ensembles [17, 19]. Both approaches
employ multiple networks to model uncertainty, where the main difference is
the use/lack-of-use of shared parameters across networks during training and
inference.

In Bayesian neural networks the weights, $\phi$, are treated as random
variables, and the posterior distribution is learned from the training data
$p(\phi|\mathbf{x},\mathbf{y})$. Then, the probability of label $y^{*}$ for a
test sample $x^{*}$ is evaluated by

$p(y^{*}|\boldsymbol{x}^{*},\mathbf{x},\mathbf{y})=\int
p(y^{*}|\boldsymbol{x}^{*},\phi)~{}p(\phi|\mathbf{x},\mathbf{y})~{}d\phi.$ (1)

However, since learning the posterior over $\phi$ and estimating Equation 1
are usually intractable, variational methods are often used, where an
approximating variational distribution is defined $q(\phi)$, and the KL-
divergence between the true posterior and the variational distribution,
$\mathrm{KL}(q(\phi)\left|\right|p(\phi|\boldsymbol{X},Y))$, is minimized.

A common practice is to set a prior $p(\phi)$ and use the Bayes rule,

$p(\phi|\boldsymbol{X},{Y})=\frac{p({Y}|\boldsymbol{X},\phi)~{}p(\phi)}{p({Y}|\boldsymbol{X})}.$
(2)

Thus, minimizing the KL-divergence is equivalent to maximizing a variational
lower bound,

$\mathcal{L}_{\mathrm{VI}}=\int
q(\phi)\log\left[p(Y|\boldsymbol{X},\phi)\right]d\phi-\mathrm{KL}\Big{(}q(\phi)\left|\right|p(\phi)\Big{)}d\phi.$
(3)

Gal & Ghahramani [6] showed that the dropout objective, when applied before
every layer, maximizes Equation 3. However, the prior $p(\phi)$ is agnostic to
the unlabeled data distribution $p(\mathbf{x})$, which may lead to unwarranted
generalization for out-of-distribution test samples. As a remedy, in this work
we propose to condition the parameters of the discriminative model, $\phi$, on
the unlabeled training data, $\boldsymbol{X}$. That is, to replace the prior
$p(\phi)$ in Equation 2 with $p(\phi|\mathbf{x})$, thereby letting unlabeled
training data to guide the posterior distribution rather than relying on some
prior assumption over the prior.

##  3 A Hierarchy of Deep Neural Networks

We first describe the key idea, then introduce a new neural architecture, and
finally, describe a stochastic inference algorithm for estimating different
types of uncertainties.

###  3.1 Key Idea

We approximate the prediction in Equation 1 by sampling from the posterior,
$\phi_{i}\sim p(\phi|\mathbf{x},\mathbf{y})$,

$P(y^{*}|\boldsymbol{x}^{*},\mathbf{x},\mathbf{y})\approx\frac{1}{m}\sum_{i=1}^{m}P(y^{*}|\boldsymbol{x}^{*},\phi_{i}).$
(4)

Since sampling from the posterior is intractable for deep neural networks, we
follow a Bayesian approach and propose a prior distribution for the
parameters. However, in contrast to the common practice of assuming a Gaussian
distribution (or some other prior independent of the data), our prior depends
on the unlabeled training data, $p(\phi|\mathbf{x})$. We define
$p(\phi|\mathbf{x})$ by first, considering a generative model for
$\mathbf{x}$, with parameters $\theta$. Next, we assume a dependency relation
between $\theta$ and $\phi$, such that the joint distribution factorizes as

$p(X,Y,\phi,\theta)=p(Y|X,\phi)~{}p(X|\theta)~{}p(\phi|\theta)~{}p(\theta).$
(5)

In essence, we assume a generative process of $\boldsymbol{X}$, to confound
the relation $\boldsymbol{X}\rightarrow\phi{}$ conditioned on $Y$. That is, in
contrast to the common practice where a “v-structure” is assumed,
$\boldsymbol{X}\rightarrow Y\leftarrow\phi$, we assume a generative function,
parametrized by $\theta$, to be a parent of $\boldsymbol{X}$ and $\phi$, as
illustrated in Figure 1. Given a training set $\\{\mathbf{x},\mathbf{y}\\}$
the posterior is

$p(\phi|\mathbf{x},\mathbf{y})\propto\int
p(\mathbf{y}|\mathbf{x},\phi)~{}p(\phi|\theta)~{}p(\theta|\mathbf{x})~{}d\theta.$
(6)

Figure 1: A causal diagram describing our assumptions.

The prior $p(\phi|\theta)$ is expected to diminish unwarranted generalization
for out-of-distribution samples—samples with low $p(\boldsymbol{X}|\theta)$.
As these samples are non-present or scarce in the training set, they have non
or negligible influence on $\theta$. Hence, for in-distribution,
$X_{\mathrm{in}}\sim p(\mathbf{x})$,
$p(\phi|\boldsymbol{X}_{\mathrm{in}},\theta)=p(\phi|\theta)$, whereas for out-
of-distribution this relation does not hold. During inference we first sample
$\theta$ from $p(\theta|\mathbf{x})$, which for an arbitrary out-of-
distribution sample is a uniform distribution. Thus, $p(\phi|\theta)$ is
expected to spread probability mass across $\phi$ for out-of-distribution
data.

Two main questions arise: 1) what generative model should be used, and 2) how
can we define the conditional relation $p(\phi|\theta)$. It is desired to
define a generative model such that the conditional distribution
$p(\phi|\theta)$ can be computed efficiently.

Our solution includes a new deep model, called BRAINet, which consists of
multiple networks coupled into a single hierarchical structure. This structure
(inter-layer neural connectivity pattern) is learned and scored from unlabeled
training data, $\mathbf{x}$, such that multiple generative structures,
$\\{\mathcal{\tilde{G}}_{1},\mathcal{\tilde{G}}_{2},\mathcal{\tilde{G}}_{3,}\ldots\\}$,
can be sampled from their posterior $p(\mathcal{\tilde{G}}|\mathbf{x})$ during
training and inference. We define $p(\theta_{i}|\mathbf{x})\equiv
p(\mathcal{\tilde{G}}_{i}|\mathbf{x})$, where for each $\theta_{i}$, the
corresponding discriminative network parameters are estimated
$\phi_{i}=\mathop{\mathrm{arg\,max}}{p(\phi|\theta_{i})}$. This can be
described as using multiple networks during inference and training, similarly
to MC-dropout [6] and Deep Ensembles [17]. However, as these structures are
sampled from a single network, they share some of their parameters,
specifically in deeper layers. This is different from MC-dropout where all the
parameters are shared across networks, and Deep Ensembles where none of the
parameters are shared.

###  3.2 BRAINet: A Hierarchy of Deep Networks

Recently, Rohekar et al. [28] introduced an algorithm, called B2N, for
learning the structure, $\mathcal{G}$, of discriminative deep neural networks
in an unsupervised manner. The B2N algorithm, learns an inter-layer
connectivity pattern, where neurons in a layer may connect to other neurons in
any deeper layer, not just to the ones in the next layer. Initially, B2N
learns a deep generative graph $\mathcal{\tilde{G}}$ with latent nodes
$\boldsymbol{H}$. This graph is constructed by unfolding the recursive calls
in the RAI structure learning algorithm [31]. RAI learns a causal
structure000A CPDAG (complete partially directed acyclic graph), a Markov
equivalence class encoding causal relations among $\boldsymbol{X}$, is learned
from nonexperimental observed data. $\mathcal{B}$ [24], and the B2N learns a
deep generative graph such that,

$p_{\mathcal{B}}(\boldsymbol{X})=\int
p_{\mathcal{\tilde{G}}}(\boldsymbol{X},\boldsymbol{H})d\boldsymbol{H}.$ (7)

Interestingly, a discriminative network structure $\mathcal{G}$ is proved to
mimic [24] a generative structure $\mathcal{\tilde{G}}$, having the exact same
structure for $(\boldsymbol{X},\boldsymbol{H})$. That is, for any
$\theta^{\prime}$, there exists $\phi^{\prime}$, which can produce the
posterior distribution over the latent variables in the generative model,
$p_{\theta^{\prime}}(\boldsymbol{H}|\boldsymbol{X})=p_{\phi^{\prime}}(\boldsymbol{H}|\boldsymbol{X},Y)$111As
the structures of $\mathcal{\tilde{G}}$ and $\mathcal{G}$ are identical,
differing only in edge direction, we assume a one-to-one mapping from each
latent node in $\mathcal{\tilde{G}}$ to its corresponding latent node in
$\mathcal{G}$..

Recently, an extension for the RAI algorithm [31] was proposed, called B-RAI
[27]. B-RAI is a Bayesian approach that learns multiple causal structures,
scores them and couples them into a hierarchy. This hierarchy is represented
by a tree, which they call GGT. The Bayesian scores of each structure is
encoded efficiently in the GGT, which allows structures
$\\{\mathcal{B}_{1},\mathcal{B}_{2},\ldots\\}$ to be sampled from this GGT
proportionally to their posterior distribution,
$P(\mathcal{B}|\boldsymbol{X})$. Based on the principles of the B2N algorithm,
we propose converting this GGT (generated by B-RAI) into a deep neural network
hierarchy. We call this hierarchy B-RAI neural network, abbreviated to
BRAINet. Then, a neural network structure, $\mathcal{G}$, can be sampled from
the BRAINet model proportionally to $P(\mathcal{B}|\boldsymbol{X})$, where
$\mathcal{G}$ has the same connectivity as a generative structure
$\mathcal{\tilde{G}}$, and where the relation in Equation 7 holds. This yields
a dependency between the generative $P(\boldsymbol{X})$ and discriminative
$P(Y|\boldsymbol{X})$ models.

####  3.2.1 BRAINet Structure Learning

Before describing the BRAINet structure learning algorithm, we provide
definitions of relevant concepts introduced by Yehezkel & Lerner [31].

######  Definition 1 (Autonomous set of nodes).

In a graph defined over $\boldsymbol{X}$, a set of nodes
$\boldsymbol{X}^{\prime}\subseteq\boldsymbol{X}$ is called autonomous given
$\boldsymbol{X}_{\mathrm{ex}}\subset\boldsymbol{X}$ if the parents’ set,
$\boldsymbol{Pa}(X)$, $\forall X\in\boldsymbol{X}^{\prime}$ is
$\boldsymbol{Pa}(X)\subset\boldsymbol{X}^{\prime}\cup\boldsymbol{X}_{\mathrm{ex}}$.

######  Definition 2 (d-separation resolution).

The resolution of a d-separation relation between a pair of non-adjacent nodes
in a graph is the size of the smallest condition set that d-separates the two
nodes.

######  Definition 3 (d-separation resolution of a graph).

The d-separation resolution of a graph is the highest d-separation resolution
in the graph.

We present a recursive algorithm, Algorithm 3.2.1, for learning the structure
of a BRAINet model. Each recursive call receives a causal structure
$\mathcal{B}$ (a CPDAG), a set of endogenous $\boldsymbol{X}$ and exogenous
$\boldsymbol{X}_{\mathrm{ex}}$ nodes, and a target conditional independence
order $n$. The CPDAG encodes $P(\boldsymbol{X}|\boldsymbol{X}_{\mathrm{ex}})$,
providing an efficient factorization of this distribution. The d-separation
resolution of $\mathcal{B}$ is assumed $n-1$.

[htbp]  BRAINet structure learning BRAINet_SL
_$(\mathcal{B},\boldsymbol{X},\boldsymbol{X}_{\mathrm{ex}},n$)_

2 Input: an initial CPDAG $\mathcal{B}$ over endogenous $\boldsymbol{X}$ &
exogenous $\boldsymbol{X}_{\mathrm{ex}}$ observed variables, and a desired
resolution $n$. Output: $L$, the deepest layer in a learned structure if  _the
maximal indegree of $\boldsymbol{X}$ in $\mathcal{B}$ is lower than $n+1$_
then $\triangleright$ exit condition

$r:=\texttt{Score}(\boldsymbol{X}|\mathcal{B})$

3 $\triangleright$ a Bayesian score (e.g., BDeu) $L:=$a gather layer for
$\boldsymbol{X}$ with score $r$

4 return _$L$_

5 for  _$t=1:s$_ do

$\mathbf{x}^{*}:=$ sample with replacement from training data $\mathbf{x}$

$\triangleright$ bootstrap sample
$\mathcal{B}^{*}:=$IncSeparation(_$\mathcal{B},n,\mathbf{x}^{*}$_)

$\triangleright$ increase d-separation resolution to $n$
$\\{\boldsymbol{X}_{\mathrm{D}},{\boldsymbol{X}_{\mathrm{A}}}_{1},\ldots,{\boldsymbol{X}_{\mathrm{A}}}_{k}\\}:=$FindAutonomous(_$\boldsymbol{X}|\mathcal{B}^{*}$_)

6 $\triangleright$ decompose for  _$i=1:k$_ do

${L_{\mathrm{A}}}_{i}:=$BRAI_NN_SL(_$\mathcal{B}^{*},{\boldsymbol{X}_{\mathrm{A}}}_{i},\boldsymbol{X}_{\mathrm{ex}},n+1$_)

7 $\triangleright$ recursively call for ancestors

$L_{\mathrm{D}}:=$BRAI_NN_SL(_$\mathcal{B}^{*},\boldsymbol{X}_{\mathrm{D}},\boldsymbol{X}_{\mathrm{ex}}\cup\\{{\boldsymbol{X}_{\mathrm{A}}}_{i}\\}_{i=1}^{k},n+1$_)

8 $\triangleright$ recursively call for descendant Create an empty layer
container $\boldsymbol{L}^{t}$ (tagged with index $t$)

9 In $\boldsymbol{L}^{t}$ create $k$ independent layers:
$L_{1}^{t},\ldots,L_{k}^{t}$

$\forall i\in\\{1,\ldots,k\\}$, connect: ${L_{\mathrm{A}}}_{i}\rightarrow
L_{i}^{t}\leftarrow L_{\mathrm{D}}$

10 $\triangleright$ connect

11 return _$L=\\{\boldsymbol{L}^{t}\\}_{t=1}^{s}$_

At the beginning of each recursive call, an exit condition is tested (line 2).
This condition is satisfied if conditional independence of order $n$ cannot be
tested (a conditional independence order is defined to be the size of the
condition set). In this case, the maximal depth is reached and an empty graph
is returned (a gather layer composed of observed nodes). From this point, the
recursive procedure will trace back, adding latent parent layers.

Each recursive call consists of three stages:

  * a)

Increase the d-separation resolution of $\mathcal{B}$ to $n$ and decompose the
input features, $\boldsymbol{X}$, into autonomous sets of nodes (lines 7–9):
one descendant set, $\boldsymbol{X}_{\mathrm{D}}$, and $k$ ancestor sets,
$\\{{\boldsymbol{X}_{\mathrm{A}}}_{i}\\}_{i=1}^{k}$.

  * b)

Call recursively to learn BRAINet structures for each autonomous set (lines
10–12).

  * c)

Merge the returned BRAINet structures into a single structure (lines 13–15).

Figure 2: An example of BRAINet structure learning with $s=2$. Network inputs,
$\boldsymbol{X}=\\{A,\ldots,E\\}$, are depicted on a gray plain. Red arrows
indicate recursive calls for learning local causal structures. Initially, no
prior information exists and a fully-connected graph over $\boldsymbol{X}$ is
assumed (a). Then, 0-order (condition set size of 0) independence is tested
between pairs of connected nodes, using two ($s$) different bootstrap samples
of the training data. This results in two CPDAGs (b, c). These CPDAGs are
further refined by recursive calls with higher order independence tests, until
no more edges can be removed (d-g). Tracing back, a deep NN is constructed by
adding NN layers, where each CPDAG leads to its corresponding neural
connectivity pattern (h-n). Each rectangle represents a layer of neurons
(denoted $L_{i}^{t}$ in Algorithm 3.2.1-line 14). A Venn diagram, above a set
of rectangles, indicates the field-of-view each layer has over
$\boldsymbol{X}$. Finally, a discriminative structure is created (o, p).

These stages are executed $s$ times (line 6), resulting in an ensemble of $s$
BRAINet structures. The deepest layers,
$\\{\boldsymbol{L}^{1},\ldots,\boldsymbol{L}^{s}\\}$, of these structures are
grouped together (line 16), while maintaining their index $1,\ldots,s$ in $L$,
and returned. Note that the caller function treats this group as a single
layer. A detailed description of Algorithm 3.2.1 can be found in Appendix A,
and complexity analysis in Appendix B. An example for learning a BRAINet
structure with $s=2$ (each recursive call returns an ensemble of two BRAINet
models) is given in Figure 2.

####  3.2.2 BRAINet Training and Inference

Figure 3: Stochastic training/inference in a BRAINet model. In a single
stochastic/training step, a subset of the network is selected. In this example
there are four possible sub-network selections: a,b,c, and d, encoded in a GGT
(e). In the GGT, one of $s$ (here $s=2$) branches, having scores
$r_{1},\ldots,r_{s}$, is sampled according to Equation 8. As an example, the
highlighted branches were sampled resulting in sub-network d. Note that, in
every possible stochastic step (a,b,c,d), all the inputs,
$\boldsymbol{X}=A,\ldots,E$, are selected, and each input is selected only
once.

The BRAINet model allows us to sample sub-networks with respect to their
relative posterior probability. The scores calculation and sub-network
selection are performed recursively, from the leaves to the root. For each
autonomous set, given $s$ sampled sub-networks and their scores,
$r_{1},\ldots,r_{s}$, returned from $s$ recursive calls, one of the $s$
results is sampled. We use the Boltzmann distribution,

$P(t;\\{r_{t^{\prime}}\\}_{t^{\prime}=1}^{s})=\frac{\exp[r_{t}/\gamma]}{\sum_{t^{\prime}=1}^{s}\exp[r_{t^{\prime}}/\gamma]},$
(8)

where $\gamma$ is a “temperature” term. When $\gamma\to\infty$, results are
sampled from a uniform distribution, and when $\gamma\to 0$ the index of the
maximal value is selected ($\arg\max$). We use $\gamma=1$ and the Bayesian
score, BDeu [11]. Finally, the sampled sub-networks, each corresponding to an
autonomous set, are merged. The score of the resulting network is the sum of
scores222Using a decomposable score, such as BDeu, the score of a CPDAG is
$r=\sum_{i=1}^{n}r(X_{i}|Pa_{i}$). of all autonomous sets merged into the
network. When training the parameters, at each step, a single sub-network is
sampled using a uniform distribution and its weights are updated. At
inference, however, there are two options. In the the first, which we call
“stochastic”, we run $T$ forward passes, each time sampling a single network
with respect to Equation 8. Then, the outputs of the sampled networks are
averaged. Figure 3 illustrates several forward passes sampled from the BRAINet
model. This is similar to dropout at inference, except that a sub-network is
sampled with respect to its posterior. Note that, in BRAINet, in contrast to
MC-dropout, weights are not sampled independently, and there is an implicit
dependency between sampling variables. In the second inference option, which
we call “simultaneous”, we run a single forward pass through the BRAINet and
recursively perform a weighted average of the $s$ activations for each
autonomous set.

####  3.2.3 BRAINet Uncertainty Estimation

Several measures of uncertainty can be computed using the BRAINet model.
Firstly, the max-softmax [29] and entropy [6] can be computed on the outputs
of a single “simultaneous” forward pass or on the average of outputs from
multiple stochastic forward passes. Secondly, using the distinct outputs of
multiple forward passes, we can compute the expected entropy,
$\mathbb{E}_{p(\phi|\mathbf{x},\mathbf{y})}\mathcal{H}[p(y^{*}|x^{*},\phi)]$,
or the mutual information,
$\mathrm{MI(y^{*},\phi|x^{*},\mathbf{x},\mathbf{y})}$, [30]. In Appendix
C-Figure 2, we qualitatively show epistemic uncertainty estimation using MI,
calculated from the BRAINet model, for images generated by VAE [16] trained on
MNIST. In addition, using the outputs of multiple stochastic passes, we can
estimate the distribution over the network output. Finally, we demonstrate an
interesting property of our method that learns a broader prior over $\phi$ as
the training set size decreases. That is, a relation between the predictive
uncertainty and the number of unique structures (connectivity patterns)
encoded in a BRAINet model (exemplified in Appendix C-Figure 1).

##  4 Empirical Evaluation

BRAINet structure learning algorithm is implemented using BNT [21] and runs
efficiently on a standard desktop CPU. In all experiments, we used MLP-layers
(dense), ReLU activations, ADAM optimization [15], a fixed learning rate, and
batch normalization [13]. Unless otherwise stated, each experiment was
repeated 10 times. Mean and STDev on test-set are reported. BRAINet structure
was learned directly for pixels in the ablation study (Section 4.1) and
calibration experiments (Section 4.2). In out-of-distribbution detection
experiments, higher level features were extracted using the first
convolutional layers of common pre-trained NN, and BRAINet structure was
learned for these features.

BRAINet parameters were learned by sampling sub-networks uniformally, and
updating the parameters by SGD. At inference, BRAINet performs Bayesian model
averaging (each sampled sub-network is weighted) using one of the two
strategies described in Section 3.2.2.

###  4.1 An Ablation Study for Evaluating the Effect of Confounding with a
Generative Process

First, we conduct an ablation study by gradually reducing the dependence of
the discriminative parameters $\phi$ on the generative parameters $\theta$,
i.e., the strength of the link $\theta\rightarrow\phi$ in Figure 1, and
measuring the performance of the resulting model for MNIST dataset [18]. In
the extreme case of disconnecting this link, the BRAINet structure will simply
become a Deep Ensembles model [17] with $s$ independent networks. Figure 4
demonstrates that even for a small dependence between $\theta$ and $\phi$, as
restricted by BRAINet, a significant improvement is achieved in performance
(high calibration and classification accuracy). $X$-axis represents the
strength of the link $\theta\rightarrow\phi$ (see Figure 1), where the values
represent the amount of mutual information that is required for a pair of
nodes in $\boldsymbol{X}$ to be considered dependent (line 8, Algorithm
3.2.1). When this value is low, all the nodes in $\boldsymbol{X}$ are
considered dependent and no structure is learned. For a mutual information
threshold of 0, a simple ensemble of $s$ networks is obtained where each
network is composed of stacked fully connected layer.

Figure 4: Results of an ablation study. The effect of conditioning the
discriminative function on the generative process, $p(\phi|\theta)$ as
measured by the test NLL, classification error, and Brier score. It is evident
that performance worsens after weakening the dependence of $\phi$ on $\theta$.

###  4.2 Calibration

BRAINet can be perceived as a compact representation of an ensemble. We
demonstrate on MNIST that it achieves higher classification accuracy and is
better calibrated than Deep Ensembles for the same model size (Figure 5).
Here, we used the simultaneous inference mode of BRAINet. Next, we evaluate
the accuracy and calibration of BRAINet as a function of stochastic forward
passes, and find it to significantly outperform Bayes-by-Backprop [1] and MC-
dropout (Figure 6). We also find that using the BRAINet structure within
either Deep Ensembles or with MC-dropout methods, further improves these later
approaches. For that we use common UCI-repository [4] regression benchmarks.
Results are reported: Appendix D-Table 1 for Deep Ensembles, and Appendix
D-Table 2 for MC-dropout. Finally, we compare BRAINet to various state-of-the-
art methods on large networks. In all benchmarks, BRAINet achieves the lowest
expected calibration error [9] (Appendix D-Table 3).

Figure 5: Performance as a function of normalized model size. Test NLL,
classification error, and Brier score of BRAINet, compared to Deep Ensembles
[17]. $X$-axis is the model size divided by the size of a single network (240K
parameters) in the Deep Ensembles model. For BRAINet, up to model size 5,
$s=2$, and from model size 7 and above, $s=3$. Different BRAINet sizes, for a
given $s$, are obtained by varying the number of neurons in the dense-layers.

Figure 6: Performance as a function of the number of forward passes (number of
sampled networks). Test NLL, classification error, and Brier score of BRAINet,
compared to MC-dropout [6] and Bayes-by-Backprop [1]. We used a BRAINet with
$s=2$. Model size is 240K parameters for MC-dropout and BRAINet, and double
for Bayes-by-Backprop.

###  4.3 Out-of-Distribution Detection

Next, we evaluate the the performance of detecting OOD samples by applying the
BRAINet structure after feature-extracting layers of common NN typologies.
First, we compare it to a baseline network and MC-dropout (We found MC-dropout
to significantly outperform Bayes-by-Backprop and Deep Ensembles on this
task). In order to evaluate the gain in performance resulting only from the
structure, we use the cross-entropy loss for parameter learning. Using other
loss functions that are suited for improving OOD detection [3, 20], may
further improve the results (see the next experiment). In this experiment we
used a ResNet-20 network [10], pre-trained on CIFAR-10 data. For MC-dropout,
we used a structure having 2 fully-connected layers, and for BRAINet, we
learned a 2-layer structure with $s=2$. Both structures have 16K parameters
and they replace the last layer of ResNet-20. SVHN dataset [23] is used as the
OOD samples. We calculated the area under the ROC and precision-recall curves,
treating OOD samples as positive classes (Table 1).

Table 1: OOD detection (SVHN dataset) by replacing the last layer of ResNet20, pre-trained on CIFAR-10, with a learned structure. Parameters are trained using cross-entropy loss. Two inference mode of BRAINet: a single simultaneous forward pass (BRAINet sm.), and multiple stochastic forward passes. MC-dropout and BRAINet use 15 forward passes (see also Appendix D-Figure 3).  |  | AUC-ROC | AUC-PR
---|---|---|---
Method | err | Max. P | Ent. | MI | E.Ent. | Max. P | Ent. | MI | E.Ent.
Baseline | 7.7 | 90.38 | 90.74 | — | — | 93.94 | 94.09 | — | —
BRAINet sm. | 8.0 | 91.62 | 92.18 | — | — | 94.91 | 94.80 | — | —
MC-dropout | 8.2 | 90.73 | 90.26 | 84.74 | 90.61 | 94.31 | 93.71 | 88.89 | 94.46
|  $\pm$0.1 |  $\pm$0.78 |  $\pm$0.81 |  $\pm$1.07 |  $\pm$0.93 |  $\pm$0.33 |  $\pm$0.34 |  $\pm$0.45 |  $\pm$0.46
BRAINet | 7.5 | 92.13 | 92.61 | 91.87 | 92.98 | 95.37 | 95.36 | 94.6 | 95.64
|  $\pm$0.1 |  $\pm$0.15 |  $\pm$0.03 |  $\pm$0.13 |  $\pm$0.10 |  $\pm$0.13 |  $\pm$0.07 |  $\pm$0.25 |  $\pm$0.05

Lastly, we demonstrate that training the parameters of BRAINet using a loss
function, specifically designed for OOD detection [3], achieves a significant
improvement over a common baseline [12] and improves state-of-the art results
(Table 2).

Table 2: OOD detection by training BRAINet parameters using a loss function designed for OOD detection. A comparison between a baseline [12], confidence-based thresholding [3], and BRAINet with the same confidence measure. Architecture: VGG-13, in-distribution: CIFAR-10, OOD: TinyImageNet. BRAINet replaces the last layer. FPR @TPR=95%: false positive rate at true positive rate of 95%. Detection error: minimum classification error over all possible thresholds. AUC-ROC and AU-PR: area under ROC and precision-recall curves. “in”/“out” indicate that in/out-of-distribution data is the positive class. An arrow indicates if lower ($\downarrow$) or higher ($\uparrow$) is better. Measure |  |  Baseline [12] |  Confidence [3] | BRAINet
---|---|---|---|---
Classification error | $\downarrow$ | 5.28 | 5.63 | 5.65
FPR @TPR=95% | $\downarrow$ | 0.438 | 0.195 | 0.124
Detection Error | $\downarrow$ | 0.120 | 0.092 | 0.076
AUC-ROC | $\uparrow$ | 0.935 | 0.970 | 0.980
AUC-PR (in) | $\uparrow$ | 0.946 | 0.974 | 0.982
AUC-PR (out) | $\uparrow$ | 0.917 | 0.965 | 0.979

##  5 Conclusions

We proposed a method for confounding the training process in deep neural
networks, where the discriminative network is conditioned on the generative
process of the input. This led to a new architecture—BRAINet: a hierarchy of
deep neural connections. From this hierarchy, local sub-networks can be
sampled proportionally to the posterior of local causal structures of the
input. Using an ablation study, we found that even a weak relation between the
generative and discriminative functions results in a significant gain in
calibration and accuracy. In addition, We found that the number of neural
connectivity patterns in BRAINet is adjusted automatically according to the
uncertainty in the input training data. We demonstrated that this enables
estimating different types of uncertainties, better than common and state-of-
the-art methods, as well as higher accuracy on both small and large datasets.
We conjecture that the resulting model can also be effective at detecting
adversarial attacks, an plan to explore this in our future work.

## References

  * Blundell et al. [2015] Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D.  Weight uncertainty in neural network.  In _International Conference on Machine Learning (ICML)_ , pp. 1613–1622, 2015.
  * Dawid [1982] Dawid, A. P.  The well-calibrated Bayesian.  _Journal of the American Statistical Association_ , 77(379):605–610, 1982.
  * DeVries & Taylor [2018] DeVries, T. and Taylor, G. W.  Learning confidence for out-of-distribution detection in neural networks.  _arXiv preprint arXiv:1802.04865_ , 2018.
  * Dua & Graff [2017] Dua, D. and Graff, C.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Efron & Tibshirani [1994] Efron, B. and Tibshirani, R. J.  _An introduction to the bootstrap_.  CRC press, 1994.
  * Gal & Ghahramani [2016] Gal, Y. and Ghahramani, Z.  Dropout as a Bayesian approximation: Representing model uncertainty in deep learning.  In _international conference on machine learning (ICML)_ , pp. 1050–1059, 2016.
  * Ghahramani [2015] Ghahramani, Z.  Probabilistic machine learning and artificial intelligence.  _Nature_ , 521(7553):452, 2015.
  * Gneiting & Raftery [2007] Gneiting, T. and Raftery, A. E.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American Statistical Association_ , 102(477):359–378, 2007.
  * Guo et al. [2017] Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q.  On calibration of modern neural networks.  In _International Conference on Machine Learning (ICML)_ , pp. 1321–1330, 2017.
  * He et al. [2016] He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Heckerman et al. [1995] Heckerman, D., Geiger, D., and Chickering, D. M.  Learning Bayesian networks: The combination of knowledge and statistical data.  _Machine learning_ , 20(3):197–243, 1995.
  * Hendrycks & Gimpel [2017] Hendrycks, D. and Gimpel, K.  A baseline for detecting misclassified and out-of-distribution examples in neural networks.  _International Conference on Learning Representations (ICLR)_ , 2017\.
  * Ioffe & Szegedy [2015] Ioffe, S. and Szegedy, C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International Conference on Machine Learning_ , pp. 448–456, 2015.
  * Izmailov et al. [2018] Izmailov, P., Podoprikhin, D., Garipov, T., Vetrov, D., and Wilson, A. G.  Averaging weights leads to wider optima and better generalization.  _arXiv preprint arXiv:1803.05407_ , 2018.
  * Kingma & Ba [2015] Kingma, D. and Ba, J.  Adam: A method for stochastic optimization.  In _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2015.
  * Kingma & Welling [2014] Kingma, D. P. and Welling, M.  Auto-encoding variational Bayes.  In _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2014.
  * Lakshminarayanan et al. [2017] Lakshminarayanan, B., Pritzel, A., and Blundell, C.  Simple and scalable predictive uncertainty estimation using deep ensembles.  In _Advances in Neural Information Processing Systems (NIPS)_ , pp. 6402–6413, 2017.
  * LeCun et al. [1998] LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Maddox et al. [2019] Maddox, W., Garipov, T., Izmailov, P., Vetrov, D., and Wilson, A. G.  A simple baseline for Bayesian uncertainty in deep learning.  _arXiv preprint arXiv:1902.02476_ , 2019.
  * Malinin & Gales [2018] Malinin, A. and Gales, M.  Predictive uncertainty estimation via prior networks.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2018.
  * Murphy [2001] Murphy, K.  The Bayes net toolbox for Matlab.  _Computing Science and Statistics_ , 33:331–350, 2001.
  * Naeini et al. [2015] Naeini, M. P., Cooper, G., and Hauskrecht, M.  Obtaining well calibrated probabilities using Bayesian binning.  In _Twenty-Ninth AAAI Conference on Artificial Intelligence_ , 2015\.
  * Netzer et al. [2011] Netzer, Y., Wang, T., Coates, A., Bissacco, A., Wu, B., and Ng, A. Y.  Reading digits in natural images with unsupervised feature learning.  In _NIPS workshop on deep learning and unsupervised feature learning_ , 2011.
  * Pearl [2009] Pearl, J.  _Causality: Models, Reasoning, and Inference_.  Cambridge university press, second edition, 2009.
  * Platt et al. [1999] Platt, J. et al.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  _Advances in large margin classifiers_ , 10(3):61–74, 1999.
  * Ritter et al. [2018] Ritter, H., Botev, A., and Barber, D.  A scalable laplace approximation for neural networks.  In _International Conference on Learning Representations (ICLR)_ , 2018.
  * Rohekar et al. [2018a] Rohekar, R. Y., Gurwicz, Y., Nisimov, S., Koren, G., and Novik, G.  Bayesian structure learning by recursive bootstrap.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2018a.
  * Rohekar et al. [2018b] Rohekar, R. Y. Y., Nisimov, S., Gurwicz, Y., Koren, G., and Novik, G.  Constructing deep neural networks by Bayesian network structure learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2018b.
  * Shiyu et al. [2018] Shiyu, L., Yixuan, L., and R., S.  Enhancing the reliability of out-of-distribution image detection in neural networks.  In _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2018.
  * Smith & Gal [2018] Smith, L. and Gal, Y.  Understanding Measures of Uncertainty for Adversarial Example Detection.  In _Uncertainty in Artificial Intelligence (UAI)_ , 2018.
  * Yehezkel & Lerner [2009] Yehezkel, R. and Lerner, B.  Bayesian network structure learning by recursive autonomy identification.  _Journal of Machine Learning Research (JMLR)_ , 10(Jul):1527–1570, 2009.

## Appendix

##  Appendix A BRAINet Structure Learning: A Detailed Description

In this section we provide a detailed description of Algorithm 1 and its three
main stages. For more details, specific to Bayesian networks (not to be
confused with Bayesian neural networks), refer to Pearl [24].

Stage a. First in line 7, a bootstrap sample $\mathbf{x}^{*}$ is created by
sampling-with-replacement from the training data, $\mathbf{x}$ (non-parametric
bootstrap). The bootstrap principle is a common approach used to approximate a
population distribution by a sample distribution [5].

Next in line 8, the bootstrap sample $\mathbf{x}^{*}$ is used for learning a
Bayesian network structure $\mathcal{B}^{*}$, with sparser connectivity than
$\mathcal{B}$, such that $\mathcal{B}$ can mimic $\mathcal{B}^{*}$,
$\mathcal{B}^{*}\preceq\mathcal{B}$ [24]. That is, for every set of parameters
$\nu$, quantifying $\mathcal{B}$, there exists a set of parameters $\nu^{*}$,
quantifying $\mathcal{B}^{*}$, such that
$p_{\mathcal{B},\nu}(X)=p_{\mathcal{B}^{*},\nu^{*}}(X)$. The $\mathcal{B}^{*}$
structure is learned by testing conditional independence of order $n$ between
pairs of nodes connected in $\mathcal{B}$. That is,
$X\mathbin{\mathchoice{\hbox to0.0pt{\hbox{$\displaystyle\perp$}\hss}\kern
3.46875pt{}\kern 3.46875pt\hbox{$\displaystyle\perp$}}{\hbox
to0.0pt{\hbox{$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{$\textstyle\perp$}}{\hbox
to0.0pt{\hbox{$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern
2.36812pt\hbox{$\scriptstyle\perp$}}{\hbox
to0.0pt{\hbox{$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern
1.63437pt\hbox{$\scriptscriptstyle\perp$}}}X^{\prime}|\boldsymbol{S}$ for
every connected pair $X\in\boldsymbol{X}$ and
$X^{\prime}\in\boldsymbol{X}_{\mathrm{ex}}$ given a condition set
$\boldsymbol{S}\subset\\{\boldsymbol{X}_{\mathrm{ex}}\cup\boldsymbol{X}\\}$ of
size $n$. Edges between conditionally independent nodes are then removed, and
the remaining edges are directed by applying two rules. First, v-structures
are identified and directed. Then, edges are continually directed, by avoiding
the creation of new v-structures and directed cycles, until no more edges can
be directed. Following the definition of Yehezkel & Lerner [31] for
d-separation resolution, we say that this function increases the graph
d-separation resolution from $n-1$ to $n$.

Finally in line 9, the procedure SplitAutonomous (line 7) identifies
autonomous sets: one descendant set, $\boldsymbol{X}_{\mathrm{D}}$, and $k$
ancestor sets,
${\boldsymbol{X}_{\mathrm{A}}}_{1},\ldots,{\boldsymbol{X}_{\mathrm{A}}}_{k}$.
This decomposition is achieved in two steps. First, the nodes having the
lowest topological order (nodes without outgoing directed edges) are grouped
into $\boldsymbol{X}_{\mathrm{D}}$, and then, $\boldsymbol{X}_{\mathrm{D}}$ is
removed (temporarily) from $\mathcal{B}$ revealing unconnected sub-structures.
The number of unconnected sub-structures is denoted by $k$ and the nodes set
of each sub-structure is denoted by ${\boldsymbol{X}_{\mathrm{A}}}_{i}$
($i\in\\{1\ldots k\\}$).

Stage b. The BRAINet structure learning algorithm is recursively called for
each autonomous (line 11) and descendant (line 12) set. An autonomous set in
$\mathcal{B}$ includes all its nodes’ parents (complying with the Markov
property) and thus, the algorithm can be called recursively and independently.
Each recursive call returns a BRAINet model: ${L_{\mathrm{A}}}_{i}$ for each
ancestor set ${\boldsymbol{X}_{\mathrm{A}}}_{i}$, and $L_{\mathrm{D}}$ for the
descendant set $\boldsymbol{X}_{\mathrm{D}}$. Note that, each recursive call,
has a smaller field-of-view (FOV) over the input $\boldsymbol{X}$ compared to
its caller function
(${\boldsymbol{X}_{\mathrm{A}}}_{i}\subset\boldsymbol{X}$).

Stage c. BRAINet models returned from the recursive calls are merged by
connecting the deepest layer of each of them in the following manner. First in
line 13, a layer container is created, denoted by $\boldsymbol{L}^{t}$, where
$t$ is an index created in line 6 and represents one of $s$ bootstrap-created
splits. Next in line 14, for the $i$-th BRAINet model, created for ancestor
set ${\boldsymbol{X}_{\mathrm{A}}}_{i}$, a layer of neurons, $L_{i}^{t}$ is
created. Finally in line 15, the BRAINet models returned from the recursive
calls are connected. The deepest layer of $L_{\mathrm{D}}$ is connected to all
the newly created layers $L_{1}^{t},\ldots,L_{k}^{t}$ and the deepest layer of
the $i$-th BRAINet model created for ${\boldsymbol{X}_{\mathrm{A}}}_{i}$ is
connected to layer $L_{i}^{t}$.

As an example, in the diagram below, we show the results of lines 7-9 of
Algorithm 1 in a tree-form (we use $s=2$). Each rectangle encapsulates the
result of multiple possible decomposition, obtained from multiple bootstrap
samples. For simplifying the explanation, we omit the distinction between
ancestor and descendant sets. For example, in the top rectangle,
$\boldsymbol{X}$ is decomposed, using statistical independence tests of order
$n$, in two ($s=2$) manners: 1)
$\\{\boldsymbol{X}_{1},\boldsymbol{X}_{2},\boldsymbol{X}_{3},\boldsymbol{X}_{4}\\}$,
2) $\\{\boldsymbol{X}_{5},\boldsymbol{X}_{6},\boldsymbol{X}_{7},\\}$. Each set
$\boldsymbol{X}_{i}$ is further recursively decomposed using statistical
independence tests of order $n+1$. For example, $\boldsymbol{X}_{7}$ is
decomposed in two manners (using different bootstrap samples).

##  Appendix B Complexity of BRAINet Structure Learning

The complexity of the proposed algorithm is essentially identical to that of
the B-RAI algorithm. That is, $\mathcal{O}(n^{k}s^{k+1})$ conditional
independence tests, and $\mathcal{O}(ns^{k})$ Bayesian scoring function
calls., where $s$ is the number of splits, $n$ is the number of input
variables ($|\boldsymbol{X}|$), and $k$ is the maximal order of conditional
independence in the data.

The running-trace of RAI in the worst-case scenario is a single path in the
GGT, and has a CI-test complexity of $\mathcal{O}(n^{k})$, and the ratio
between B-RAI and RAI is $\sum_{i=0}^{k}s^{i}$.

For the Bayesian scoring function, the complexity is $\mathcal{O}(ns^{k})$, as
only the leaves (when the exit condition is satisfied) are scored. Note that
the worst-case scenario is the case where the true underlying graph is a
complete graph, which is not typical in real-world cases. In practice,
significantly fewer CI tests and scoring operations are performed.

##  Appendix C BRAINet Uncertainty Estimation

###  C.1 A Relation between Generative Uncertainty and Predictive Uncertainty

During the construction of a BRAINet model, multiple connectivity patterns at
each point in the network are learned using different bootstrap samples. Thus,
when the epistemic uncertainty of $\theta$ is high, connectivity patterns
learned from different bootstrap samples are likely to be dissimilar. However,
when the epistemic uncertainty of $\theta$ lowers (e.g., for a larger training
set), the connectivity patterns for any two bootstrap samples are more likely
to be similar, at which point the aleatoric uncertainty dominates, and further
reducing epistemic uncertainty (e.g., by adding more training data) will not
result in a reduction of the number of unique connectivity patterns.

Figure 7: Number of unique structures (neural connectivity patterns) embedded
in a single BRAINet ($s=3$) for MNIST, as a function of the training set size.
As the epistemic uncertainty increases (small training sets), more unique
structures are automatically encoded in BRAINet, resulting in a broader prior
over the network parameters. At the other end, as the training set size
increases, the number of unique structures decreases and converges to a number
greater than one, indicating the existence of an aleatoric uncertainty.
Results are averaged over 5 experiments; error bars indicate standard
deviation.

###  C.2 An example for Epistemic Uncertainty Estimation using MI Criterion

In this section we repeat the experiment described by Smith & Gal [30]. Result
is in Figure 8, where we also plot 12 images generated by a VAE: 6 having low
MI score, and 6 having high MI score.

Figure 8: Epistemic uncertainty estimate using BRAINet (s = 2) as measured by
mutual information, visualized on the latent space of a VAE on MNIST. In this
experiment, as suggested by Smith & Gal [30], a VAE with latent space size 2
was trained and the axes are the values of these 2 latent variables. Brighter
pixels correspond to a higher mutual information, i.e. a higher uncertainty.
Images were then generated by setting values for each latent variable, 100
values in $\left[-10:+10\right]$. Colored pixels correspond to the training
data. As demonstrated, ambiguous images, $\\{\mathrm{A,...,F}\\}$, yield high
epistemic uncertainties, as opposed to the clearer images of
$\\{\mathrm{G,...,L}\\}$.

##  Appendix D Experiments

This section contains tables and figures referenced in the paper. See captions
for a detailed explanation.

Firstly, we wish to emphasize an important property of the BRAINet structure,
which allows us to sample networks that are significantly smaller than the
overall structure. This leads to a significantly smaller computational cost,
during training and inference, than other methods. For example, in a BRAINet
for MNIST, the average number of operations (multiply and add) during
inference is $\sim\mathbf{5.5\times}$ smaller than other methods (and has
$\sim\mathbf{5.5\times}$ fewer parameters to train each epoch). Moreover, in
contrast to MC-dropout, which samples neurons from each layer, BRAINet samples
full layers, making it computationally efficient using common hardware. In all
our experiments, we report only the size of the _full_ BRAINet structure and
ignore the sizes of the sampled networks.

Figure 9: Area under ROC and precision-recall curves as a function of the
number of stochastic forward passes. It is evident that BRAINet achieves high
AUC even for a samll number of forward passes, compared to MC-dropout. This
result corresponds to the OOD-detection experiment described in Table 1 of the
main paper. Architecture: ResNet20, in-distribution: CIFAR-10, OOD: SVHN.

Table 3: The effect of conditioning the discriminative function on a MAP
estimation of a generative network (BRAINet) for Deep Ensembles [17].
Regression Benchmark datasets comparing RMSE and NLL.
| RMSE | NLL
---|---|---
Dataset | ensemble | BRAINet+ensemble | ensemble | BRAINet+ensemble
Concrete |  5.46$\pm$ 0.52 | 4.17$\pm$ 0.59 |  2.92$\pm$ 0.16 | 2.56$\pm$ 0.13
Boston Housing |  2.90$\pm$ 0.71 |  2.17$\pm$ 0.56 |  2.37$\pm$ 0.25 | 1.88$\pm$ 0.21
Power Plant |  4.11$\pm$ 0.16 |  4.10$\pm$ 0.14 |  2.82$\pm$ 0.03 |  2.8$\pm$ 0.03
Yacht |  3.09$\pm$ 1.03 | 1.10$\pm$ 0.36 |  2.25$\pm$ 0.39 | 1.09$\pm$ 0.14
Kin8nm |  0.09$\pm$ 0.00 | 0.08$\pm$ 0.00 |  -1.18$\pm$ 0.02 | -1.23$\pm$ 0.03
Energy |  0.95$\pm$ 0.26 |  0.71$\pm$ 0.10 |  1.16$\pm$ 0.25 |  1.02$\pm$ 0.13
Naval Propulsion Plant |  0.00$\pm$ 0.00 |  0.00$\pm$ 0.00 |  -3.55$\pm$ 0.06 | -3.68$\pm$ 0.05
Wine |  0.64$\pm$ 0.04 |  0.59$\pm$ 0.05 |  0.92$\pm$ 0.08 | 0.76$\pm$ 0.05
Protein |  4.62$\pm$ 0.17 | 4.16$\pm$ 0.15 |  2.79$\pm$ 0.03 | 2.60$\pm$ 0.09
Table 4: The effect of conditioning the discriminative function on a MAP
estimation of a generative network (BRAINet) for MC-dropout [6]. Regression
Benchmark datasets comparing RMSE and NLL.
| RMSE | NLL
---|---|---
Dataset | MC-dropout | BRAINet+dropout | MC-dropout | BRAINet+dropout
Concrete |  4.83$\pm$ 0.52 | 3.09$\pm$ 0.78 |  2.92$\pm$ 0.09 | 2.64$\pm$ 0.12
Boston Housing |  2.80$\pm$ 0.52 | 1.90$\pm$ 0.38 |  2.39$\pm$ 0.15 | 2.16$\pm$ 0.08
Power Plant |  4.02$\pm$ 0.18 |  4.00$\pm$ 0.16 |  2.80$\pm$ 0.05 |  2.79$\pm$ 0.04
Yacht |  1.42$\pm$ 0.48 | 0.61$\pm$ 0.25 |  1.60$\pm$ 0.13 | 1.38$\pm$ 0.06
Kin8nm |  0.10$\pm$ 0.00 | 0.08$\pm$ 0.00 |  -0.95$\pm$ 0.03 | -1.10$\pm$ 0.03
Energy |  0.88$\pm$ 0.13 | 0.56$\pm$ 0.08 |  1.62$\pm$ 0.03 | 1.44$\pm$ 0.02
Naval Propulsion Plant |  (1.8$\pm$ 0.2)$e^{-3}$ | (1.2$\pm$ 0.3)$e^{-3}$ |  -4.22$\pm$ 0.01 | -4.35$\pm$ 0.02
Wine |  0.62$\pm$ 0.04 | 0.54$\pm$ 0.04 |  0.93$\pm$ 0.06 | 0.82$\pm$ 0.05
Protein |  3.60$\pm$ 0.03 | 3.46$\pm$ 0.19 |  2.69$\pm$ 0.01 |  2.67$\pm$ 0.03
Table 5: Comparison between BRAINet and various state-of-the-art methods on
large networks. In all benchmarks, BRAINet achieves the lowest expected
calibration error [9].
Dataset | Model | SGD | SWA | SWAG- | SWAG | KFAC- | SWA- | SWA- | BRAINet
---|---|---|---|---|---|---|---|---|---
|  |  |  | Diag |  | Laplace | Dropout | Temp |
CIFAR-10 | VGG-16 | 0.0483 | 0.0408 | 0.0267 | 0.0158 | 0.0094 | 0.0284 | 0.0366 | 0.0090
CIFAR-10 | PreResNet-164 | 0.0255 | 0.0203 | 0.0082 | 0.0053 | 0.0092 | 0.0162 | 0.0172 | 0.0036
CIFAR-10 | WRN28x10 | 0.0166 | 0.0087 | 0.0047 | 0.0088 | 0.0060 | 0.0094 | 0.0080 | 0.0040
CIFAR-100 | VGG-16 | 0.1870 | 0.1514 | 0.0819 | 0.0395 | 0.0778 | 0.1108 | 0.0291 | 0.0247