# Real-Time Limited-View CT Inpainting and Reconstruction with Dual Domain
Based on Spatial Information

# Real-Time Limited-View CT Inpainting and Reconstruction with Dual Domain
Based on Spatial Information

###### Abstract

Low-dose Computed Tomography is a common issue in reality. Current reduction,
sparse sampling and limited-view scanning can all cause it. Between them,
limited-view CT is general in the industry due to inevitable mechanical and
physical limitation. However, limited-view CT can cause serious imaging
problem on account of its massive information loss. Thus, we should
effectively utilize the scant prior information to perform completion. It is
an undeniable fact that CT imaging slices are extremely dense, which leads to
high continuity between successive images. We realized that fully exploit the
spatial correlation between consecutive frames can significantly improve
restoration results in video inpainting. Inspired by this, we propose a deep
learning-based three-stage algorithm that hoist limited-view CT imaging
quality based on spatial information. In stage one, to better utilize prior
information in the Radon domain, we design an adversarial autoencoder to
complement the Radon data. In the second stage, a model is built to perform
inpainting based on spatial continuity in the image domain. At this point, we
have roughly restored the imaging, while its texture still needs to be finely
repaired. Hence, we propose a model to accurately restore the image in stage
three, and finally achieve an ideal inpainting result. In addition, we adopt
FBP instead of SART-TV to make our algorithm more suitable for real-time use.
In the experiment, we restore and reconstruct the Radon data that has been cut
the rear one-third part, they achieve PSNR of 40.209, SSIM of 0.943, while
precisely present the texture.

Index Terms—  limited-view CT, deep learning, spatial information, domain
transformation, real-time

##  1 Introduction

Computed Tomography has been successfully applied in medicine, biology,
industry and other fields, providing huge help for industrial production,
medical research and people’s daily life [1]. Nevertheless, the radiation dose
brought by CT scanning may somehow have a negative effect on human body that
cannot be neglect. Thus, it is crucial for CT scanning to lower its radiation
dose [2] in accordance with ALARA (as low as reasonably achievable) [3]. Low-
dose Computed Tomography (LDCT) can be realized through current reduction,
sparse sampling and limited-view scanning. Among these, limited-view CT is
really general because that we often encounter mechanical and physical
restriction in the industry which makes it difficult for the machine to scan
through an object. Despite the general application of limited-view CT, its
imaging leads to some grievous problems like blur [4], artifacts [5, 6, 7, 8]
and low signal-to-noise ratio [9, 1], they undoubtedly have a great influence
on clinical diagnosis. Thus, it is crucial for researchers to fully utilize
the limited prior information to effectively complement the fragmentary data.

Traditional analytical reconstruction algorithms, such as FBP [10], have high
requirements for data integrifty. When the radiation dose is reduced,
artifacts in reconstructed images will increase rapidly [11]. In order to
upgrade the quality of reconstructed images, many researchers have proposed
various algorithms for LDCT imaging reconstruction, and we conclude them into
several paths that are presented in Fig.1 for better comprehension.

Fig. 1: The technology roadmap of prevailing CT inpainting and reconstruction
algorithms, the dash line in this figure refers to the reconstruction step
from the Radon domain to the image domain through FBP or SART-TV.

Iterative Reconstruction Algorithms are represented by the red line in Fig.1,
which can directly reconstruct damaged Radon data into target results in the
image domain. Model-based iterative reconstruction (MBIR) algorithm [12], also
known as statistical image reconstruction (SIR) method, combines the modeling
of some key parameters to perform high-quality reconstruction of LDCT. Using
image priors in MBIR can effectively improve the image reconstruction quality
of LDCT scans [13, 14], while still have the high computational complexity.

In addition to the prior information, various regularization methods have
played a crucial role in iterative algorithms of CT reconstruction. The most
typical regularization method is the total variation (TV) method [15]. In the
light of TV, researchers came up with more reconstruction methods, such as TV-
POCS [16], TGV [17] and SART-TV [18] which was proposed on the basis of SART
[19]. Those algorithms can suppress image artifacts to a certain extent so as
to improve imaging quality. In addition, dictionary learning is often used as
a regularizer in MBIR algorithms [20, 21, 22, 23], and multiple dictionaries
are beneficial to reducing artifacts caused by limited-view CT reconstruction.

With the development of computing power, deep learning-based methods [24, 25,
26, 9, 27, 28, 29] have been applied to the restoration of LDCT reconstructed
images in recent years. The methods can be roughly divided into the below
three categories.

Image Inpainting algorithms are presented by blue lines in Fig.1, they firstly
reconstruct the damaged Radon data into the damaged image with artifacts, then
reduce the artifacts and noises in the image domain. Lots of researchers are
currently using convolutional neural network (CNN) and deep learning
architecture to perform this procedure [30, 1, 31, 32, 5, 33, 34, 6, 35, 7,
36]. Zhang et al [30] proposed a data-driven learning method based on deep
CNN. RED-CNN [1] combines the autoencoder, deconvolutional network and
shortcut connections into the residual encoder-decoder CNN for LDCT imaging.
Kang et al [31] applied deep CNN to the wavelet transform coefficients of LDCT
images, used directional wavelet transform to extract the directional
component of artifacts. Wang et al [33] developed a limited-angle
translational CT (TCT) image reconstruction algorithm based on U-Net [34].
Since Goodfellow et al. proposed Generative Adversarial Nets (GAN) [35] in
2014, GAN has been widely used in various image processing tasks, including
the post-processing of CT images. Xie et al. [7] proposed an end-to-end
conditional GAN with joint loss function, which can effectively remove
artifacts.

Fig. 2: The overall architecture of our proposed three-stage restoration and
reconstruction algorithm for limited-view CT imaging.

Sinogram Inpainting algorithms are presented by green lines in Fig.1, they
firstly restore the missing part in the Radon domain, then reconstruct it into
the image domain to get the final result [37, 38, 39, 40, 41]. Li et al. [37]
proposed an effective GAN-based repairing method named patch-GAN, which trains
the network to learn the data distribution of the sinogram to restore the
missing sinogram data. In another paper [38], Li et al. proposed SI-GAN on the
basis of [32], using a joint loss function combining the Radon domain and the
image domain to repair “ultra-limited-angle” sinogram. In 2019, Dai et al.
[39] proposed a limited-view cone-beam CT reconstruction algorithm. It slices
the cone-beam projection data into the sequence of two-dimensional images,
uses an autoencoder network to estimate the missing part, then stack them in
order and finally use FDK [42] for three-dimensional reconstruction. Anirudh
et al. [40] transformed the missing sinogram into a latent space through a
fully convolutional one-dimensional CNN, then used GAN to complement the
missing part. Dai et al. [41] calculated the geometric image moment based on
the projection-geometric moment transformation of the known Radon data, then
estimated the projection-geometric moment transformation of the unknown Radon
data based on the geometric image moment.

Sinogram Inpainting and Image Refining algorithms are presented by yellow
lines in Fig.1, they firstly restore the missing part in the Radon domain,
then reconstruct the full-view Radon data into the image domain so as to
finely repair the image to obtain higher quality [43, 44, 45, 46, 8]. In 2017,
Hammernik et al. [43] proposed a two-stage deep learning architecture, they
first learn the compensation weights that account for the missing data in the
projection domain, then they formulate the image restoration problem as a
variational network to eliminate coherent streaking artifacts. Zhao et al.
[44] proposed a GAN-based sinogram inpainting network, which achieved
unsupervised training in a sinogram-image-sinogram closed loop. Zhao et al.
[45] also proposed a two-stage method, firstly they use an interpolating
convolutional network to obtain the full-view projection data, then use GAN to
output high-quality CT images. In 2019, Lee et al. [46] proposed a deep
learning model based on fully convolutional network and wavelet transform. In
the latest research, Zhang et al. [8] proposed an end-to-end hybrid domain CNN
(hdNet), which consists of a CNN operating in the sinogram domain, a domain
transformation operation, and a CNN operating in the image domain.

Inspired by the combination of the two stages, we implement Radon data
completion through our proposed adversarial autoencoder (AAE) in stage one. In
the second and third stage, after enriching the information through Radon data
completion, we construct the Radon data into the image domain and realize the
image inpainting in a ”coarse-to-fine” [47] manner.

However, all of the above algorithms merely focus on a single image slice
while neglecting the abundant spatial correlation between consecutive image
slices. Consequently, these algorithms may still have trouble to reach an
ideal level of limited-view CT inpainting and reconstruction that can
precisely presents the image texture. During our investigation of video
inpainting [48, 49], we realize the significance of making full use of spatial
correlation and continuity between consecutive image slices. Therefore, we
propose an origin cascade model in stage two called Spatial-AAE to fully
utilize the spatial continuity, thereby breaking the limitation of two-
dimensional space.

It is also worth mentioning that, unlike other current limited-view CT
inpainting and reconstruction algorithms, we use FBP [10] instead of SART-TV
[18] to speed up the reconstruction process. Besides, our models do not limit
resolution of the input data, therefore can be well generalized to various
datasets. In our experiments, we compare our algorithm with the other four
prevalent algorithms under four sorts of damaged data, exhibiting its
prominent performance and robustness.

The organization of this paper is as follows, Sec II presents the design
details of our proposed algorithm and models, Sec III shows our experimental
results, and we finally conclude our research work in Sec IV.

##  2 Methods

This paper proposes a three-stage restoration and reconstruction algorithm for
limited-view CT imaging, and its overall architecture is shown in Fig.2. In
the first stage, after the limited-view Radon data is preprocessed, we input
it into the Adversarial Autoencoder we designed for data completion to obtain
the full-view Radon data. In the second stage, the output of stage one is
first reconstructed into the image domain, and combined with two consecutive
slices before and after to form a group, then we sent this group into our
proposed Spatial-AAE model to perform image restoration based on spatial
information. It is worth noting that through the above two stages of
restoration and reconstruction, most of the texture in the image ground truth
can be restored, but the result still cannot clearly reflect the precise
details, which may pose obstacles for the practical applications. Therefore,
we built the Refine-AAE high-precision inpainting network in stage three,
utilizing the idea of ”coarse-to-fine” [47] in deep learning to refine the
image in patches. The details of our algorithm are shown below.

###  2.1 Data Preprocessing

In order to provide more prior information, we adopt the data preprocessing
method from paper [5], as shown in Fig.3. For the limited-view Radon data
$\boldsymbol{R}_{lv}$, we first transform it into the image data
$\boldsymbol{I}_{recon}$ through inverse radon transformation, and then
convert the image into the full-view Radon data $\boldsymbol{R}_{fv}$ through
Radon transformation. We crop this full-view Radon data for preliminary
completion of the missing part in the original data, so as to obtain the fused
Radon data $\boldsymbol{R}_{merge}$.

Fig. 3: Procedure of data preprocessing.

###  2.2 Algorithm Pipeline

####  2.2.1 Stage 1: Limited-view Data Completion in the Radon Domain

For the input limited-view Radon data, we need to apply it as the prior
information to perform angle completion in the first stage. Due to the fact
that U-Net [34] is widely use in medical imaging currently, we propose an
adversarial autoencoder with U-Net as the backbone. Its overall architecture
is shown in Fig.4.

Fig. 4: The overall architecture of our proposed adversarial autoencoder in
stage one.

Fig. 5: (a) is the result of reconstructing the original limited-view Radon
data into the image; (b) is the result of reconstructing the full-view Radon
data generated from stage one into the image; (c) is the ground truth of the
image.

We modified U-Net as the autoencoder in our adversarial autoencoder, which
includes an encoder that downsamples the image to extract the representative
feature and a decoder that upsamples the feature to restore the image. The
precise structure of our autoencoder can be seen from TABLE I, where (Ic, Oc)
represents the in-channel and out-channel of the convolutional layer. In its
convolutional layers, the kernel size is 3$\times$3, the stride and padding
are both 1, and the kernel size is 2$\times$2 in all of its pooling layers. In
all of its deconvolution layers, the kernel size is 2$\times$2, and the stride
is 1. In order to upgrade the model’s ability of restoration, we combine this
autencoder with a discriminator whose structure is the same as the encoder
shown in TABLE I. As can be seen from Sec IV, adding this discriminator can
effectively improve the model’s performance.

Fig. 6: The overall architecture of our proposed Spatial-AAE model in stage
two.

####  2.2.2 Stage 2: Image Restoration Based on Spatial Information

Fig.5 shows that after we reconstruct the output from stage one into the image
domain, the image texture can be partly restored, while there are still some
artifacts and blurry area that can bring severe obstacles for clinical
diagnosis. Therefore, in the second stage, we propose the Spatial-AAE model
based on the spatial correlation between consecutive image slices to
significantly improve the quality of damaged image.

According to our knowledge, in previous studies of CT imaging restoration and
reconstruction algorithms, scholars seemed to neglect the rich spatial
information between consecutive image slices, and only repaired and
reconstructed images in two-dimensional space. During the process of
investigating and comparing the fields of image inpainting and video
inpainting, we were surprised to find that the third dimension usually
contains rich data coherence and continuity, which is very beneficial for
restoring successive images. Thus, we suppose that the effective use of the
third-dimensional information may remarkably improve the quality of restored
images. Inspired by the utilization of the third-dimensional information in
FastDVDNet [46], we come up with the Spatial-AAE network, whose overall
architecture is shown in Fig.6, it can be divided into Spatial autoencoder and
discriminator.

Table 1: Details of the Autoencoder in Stage One (a) Encoder |  | (b) Decoder
---|---|---
Layer | (Ic, Oc) | Layer | (Ic, Oc)
Conv1_1 | (1, 32) | UpConv6 | (512, 256)
Conv1_2 | (32, 32) | Concat | [UpConv6, Conv4]
Pool1 | Maxpool | Conv6_1 | (512, 256)
Conv2_1 | (32, 64) | Conv6_2 | (256, 256)
Conv2_2 | (64, 64) | UpConv7 | (256, 128)
Pool2 | Maxpool | Concat | [UpConv7, Conv3]
Conv3_1 | (64, 128) | Conv7_1 | (256, 128)
Conv3_2 | (128, 128) | Conv7_2 | (128, 128)
Pool3 | Maxpool | UpConv8 | (128, 64)
Conv4_1 | (128, 256) | Concat | [UpConv8, Conv2]
Conv4_2 | (256, 256) | Conv8_1 | (128, 64)
Pool4 | Maxpool | Conv8_2 | (64, 64)
Conv5_1 | (256, 512) | UpConv9 | (64, 32)
Conv5_2 | (512, 512) | Concat | [UpConv9, Conv1]
|  |  | Conv9_1 | (64, 32)
|  |  | Conv9_2 | (32, 12)
|  |  | Conv9_3 | (12, 1)

The input of the spatial autoencoder is five consecutive image slices
$S=\\{\boldsymbol{s}_{i-2},\boldsymbol{s}_{i-1},\boldsymbol{s}_{i},\boldsymbol{s}_{i+1},\boldsymbol{s}_{i+2}\\}$,
we divide them into three sets of data
$S_{1}=\\{\boldsymbol{s}_{i-2},\boldsymbol{s}_{i-1},\boldsymbol{s}_{i}\\},S_{2}=\\{\boldsymbol{s}_{i-1},\boldsymbol{s}_{i},\boldsymbol{s}_{i+1}\\}$
and
$S_{3}=\\{\boldsymbol{s}_{i},\boldsymbol{s}_{i+1},\boldsymbol{s}_{i+2}\\}$.
Then, they are sent into the AE block respectively, and their output is
concatenated as
$S^{\prime}=\\{\boldsymbol{s}^{\prime}_{i-1},\boldsymbol{s}^{\prime}_{i},\boldsymbol{s}^{\prime}_{i+1}\\}$,
this set of data is input into the AE block again to obtain the final restored
result. The spatial autoencoder network can be expressed as (1), where $F$ is
the spatial autoencoder model and $G$ is the AE block. The specific details of
the AE block and discriminator in Fig.6 can be seen from TABLE 1, they are the
same as they are in the AAE model of stage one.

$\boldsymbol{s}^{\prime\prime}_{i}=F(S)=G\left(G(S_{1}),G(S_{2}),G(S_{3})\right)$
(1)

####  2.2.3 Stage 3: Image Refining on Patches

It can be seen from Fig.7 that after the above two stages of dual-domain
combined inpainting and reconstruction, the original limited-view Radon data
can be restored to a relatively satisfying extent. However, the overall
details are still not precise enough.

Fig. 7: (a) contains the result of reconstructing the full-view Radon data
from stage one output into images; (b) contains the result of reconstructing
the full-view Radon data from stage two output into images; (c) contains the
ground truth of images.

Therefore, in the third stage, we utilize the idea of “coarse to fine” in deep
learning to propose the Refine-AAE model, so as to further refine the texture
of repaired images. The overall structure of the Refine-AAE network can be
seen from Fig.8. Give the input image $\boldsymbol{I}_{input}$, the model
divides it into four patches and concatenate them into a set of sequence
$\\{\boldsymbol{I}_{p1},\boldsymbol{I}_{p2},\boldsymbol{I}_{p3},\boldsymbol{I}_{p4},\\}$,
We send it into the autoencoder for inpainting in patches and obtain the
output as
$\\{\boldsymbol{I}^{{}^{\prime}}_{p1},\boldsymbol{I}^{{}^{\prime}}_{p2},\boldsymbol{I}^{{}^{\prime}}_{p3},\boldsymbol{I}^{{}^{\prime}}_{p4},\\}$.
The model integrates this output into $\boldsymbol{I}_{pred}$ and combines it
with the ground truth $\boldsymbol{I}_{GT}$ into pair for discriminator’s
judgment.

The autoencoder and discriminator in the Refine-AAE model are the same as the
Spatial-AAE model, they can be seen from TABLE I.

Fig. 8: The overall architecture of our proposed Refine-AAE model in stage
three.

###  2.3 Loss Function

In all three stages, we use multi-loss function to optimize the autoencoder
model, it can be expressed as (2).

$l_{AE}=\alpha_{1}l_{MSE}+\alpha_{2}l_{adv}+\alpha_{3}l_{reg}$ (2)

$l_{MSE}$ calculates the mean square error between the restored image and the
ground truth image, it is widely used in various image inpainting tasks
because it can provide an intuitive evaluation for the model’s prediction. The
expression of $l_{MSE}$ can be seen from (3).

$l_{MSE}=\dfrac{1}{W\times
H}\sum_{x=1}^{W}\sum_{y=1}^{H}\left(\boldsymbol{I}_{x,y}^{GT}-G_{AE}(\boldsymbol{I}^{input})_{x,y}\right)^{2}$
(3)

where $G_{AE}$ is the auto-encoder, $\boldsymbol{I}^{GT}$ and
$\boldsymbol{I}^{input}$ are the ground truth image and the input image, $W$
and $H$ are the width and height of the input image respectively. $l_{adv}$
refers to the adversarial loss. The autoencoder can fool the discriminator by
making its prediction as close to the ground truth as possible, so as to
achieve the ideal image restoration outcome. Its expression can be seen from
(4).

$l_{adv}=1-D\left(G_{AE}(\boldsymbol{I}^{input})\right)$ (4)

where $D$ is the discriminator and $G_{AE}$ is the autoencoder. $l_{reg}$ is
the regularization term of our multi-loss function. Since noises may have a
huge impact on the restoration result, we add a regularization term to
maintain the smoothness of the image and also avoid the problem of
overfitting. TV Loss is commonly used in image analysis tasks, it can reduce
the difference between adjacent pixel values in the image to a certain extent.
Its expression can be seen from (5).

$l_{reg}=\dfrac{1}{W\times H}\sum_{x=1}^{W}\sum_{y=1}^{H}\left\|\nabla
G_{AE}(\boldsymbol{I}^{input}_{x,y})\right\|$ (5)

where $G_{AE}$ is the auto-encoder, $\boldsymbol{I}_{input}$ is the input
image, $W$ and $H$ are the width and height of the input image respectively.
$\nabla$ calculates the gradient, $\left\|\right\|$ calculates the norm.

For the optimization of the discriminator, the loss function should enable the
discriminator to better distinguish between real and fake inputs. The loss
function can be seen from (6).

$l_{DIS}=1-D(\boldsymbol{I}^{GT})+D\left(G_{AE}(\boldsymbol{I}^{input})\right)$
(6)

where $D$ is the discriminator, $G^{AE}$ is the auto-encoder,
$\boldsymbol{I}^{GT}$ and $\boldsymbol{I}^{input}$ are the ground truth image
and the input image respectively. The discriminator outputs a scalar between
zero and one, when the output is closer to 1, the discriminator thinks that
the input is more likely to be real. On the opposite, when the output is
closer to 0, it thinks the input is more likely to be fake. Therefore,
$1-D(\boldsymbol{I}^{GT})$ makes the output closer to one when the
discriminator inputs real images, and
$D\left(G_{AE}(\boldsymbol{I}^{input})\right)$ makes the output closer to zero
when the discriminator inputs fake images generated by the autoencoder.

##  3 Experiment

Our experiment data comes from 1000 cases in the LIDC-IDRI [50] dataset. We
divided cases 1 to 200 into the test set, cases 201 to 400 into the validation
set, and cases 401 to 1000 into the training set. The CT imaging (size
512$\times$512) is stored as DCM files in the LIDC-IDRI dataset. After
processing it as an array, we reconstruct it to the Radon domain (size
512$\times$180) through the Radon transformation, and perform post-60-degree
clipping on it as the input data of the overall model. During the training
process, we set the learning rate to 1e-4, using ADAM [51] as our model
optimizer, and Leaky ReLU [52] as the nonlinear activation. For the multi-loss
function, we refer to the method in paper [53], where $\alpha_{1}$,
$\alpha_{2}$ and $\alpha_{3}$ are set to 1, 1e-3, and 2e-8 respectively. It is
worth mentioning that there is no fully connected layer in our model, so it
can flexibly handle input images of different resolutions and be applied to
various datasets. In addition, unlike other deep learning-based algorithms,
the reconstruction part of our algorithm adopts FBP instead of SART-TV which
requires a relatively high level of computational complexity, so our method
can be better applied to practical application scenarios such as clinical
diagnosis. Although FBP takes much shorter time than SART-TV, its
reconstruction results have a certain gap with SART-TV. In order to realize
the practical application value of our algorithm, we manage to compensate the
performance of FBP through the superiority of our model design. Also, we
increase the damage degree of Radon data in 4.2 to test the robustness of our
algorithm. We create four types of damaged Radon data and use this algorithm
to repair and reconstruct them. The experimental results prove that our
algorithm can effectively restore these data, thus owns outstanding
robustness. In 4.1, we conduct ablation experiments on models of each stage to
prove the necessity and effectiveness of our structural design. In 4.2, we
compared our algorithm with other four types of algorithms, and test these
algorithms on four various degrees of damaged data.

###  3.1 Ablation Study

####  3.1.1 Stage1

We first explore the necessity of fusing data in the Radon domain (refers to
Fig.3). We input the directly cut Radon data and the fused Radon data into the
stage one model shown in Fig.4 for data completion, and compare their outputs
with the Radon ground truth. The experimental results can be seen in TABLE II,
CR stands for the directly cut Radon, MR stands for the merged Radon, RCR
stands for the restored CR from stage1, RMR stands for the restored MR from
stage1.

It can be concluded from TABLE II that the fused Radon data can obtain better
experimental results due to its richer prior information, and provide more
texture for the subsequent image restoration steps. The visualized results can
be seen in Fig.9.

Table 2: Details of Various Data Preprocessing Methods |  CR |  MR |  RCR |  RMR
---|---|---|---|---
PSNR |  8.714 |  18.196 |  38.549 |  48.181
SSIM |  0.656 |  0.936 |  0.987 |  0.995

Fig. 9: Visualized results obtained from different data preprocessing methods,
(a) is the directly cut Radon data; (b) is the restored result of (a); (c) is
the fused Radon data; (d) is the restored result of (c); (e) is the Radon
ground truth.

In addition, we also explore the architecture of stage one’s adversarial
autoencoder model, and proved that it is essential to add the discriminator
reasonably. We restore the input data with: (1) The autoencoder shown in TABLE
I (a); (2) Combination of the autoencoder and the discriminator in TABLE I,
their experimental results can be seen from TABLE III.

Table 3: Results of Using Diffrent Model Structure in Stage One |  AE |  AE + D
---|---|---
PSNR |  40.129 |  48.181
SSIM |  0.983 |  0.995

It can be summarized from the above data that adding a discriminator can
greatly improve the data completion result. It can help stage one model to
improve the sinogram data PSNR by a relatively large margin.

Fig. 10: Visualized results obtained from different model structure in stage
one, (a) is the input fused Radon data; (b) is the restored Radon data from
structure AE; (c) is the restored Radon data from structure AE+D; (d) is the
Radon ground truth.

From the visualized comparison in Fig.10, we can see that if we only use this
single autoencoder, the inpainting result would have a large blurred area, and
adding the discriminator can improve this situation.

####  3.1.2 STAGE2

For the image restoration task in this stage, we adopt the Spatial-AAE model
described in 3.2 to make full use of spatial information. In order to reflect
the prominence of this structure, we compare this model with the AAE model
from stage one, which does not contain any spatial structure. For the same
input fused Radon data, the experimental results can be seen in TABLE IV. It
can be seen from the results that, due to the fact that the Spatial-AAE model
makes full use of the third-dimensional prior information, it can effectively
improve the overall performance of stage two.

Table 4: Results of Using Different Model Structure in Stage Two |  AAE |  Spacial-AAE
---|---|---
PSNR |  37.384 |  39.646
SSIM |  0.929 |  0.940

####  3.1.3 STAGE3

In this stage, the input image is divided and concatenated, and then sent to
the Refine-AAE model for finer inpainting. We believe that the way of
intercepting patches during the training process will have a certain impact on
the experimental results, so we test the following three types of interception
methods (As shown in Fig.11): (1) Randomly crop four patches (size
256$\times$256) from the input image (size 512$\times$512); (2) Crop the four
corners out of the input image; (3) Crop the four corners out of the input
image, and then adjust them into the same pattern through different flipping
method.

All of the methods above get an array of size (4, 256, 256), we input it into
the Refine-AAE model (refers to Fig.7) to finely repair the image, and the
experimental results of these three methods are shown in TABLE V.

Table 5: Results of Using Patch-Cropping Methods in Stage Three | Random Crop | Corner Crop | Corner Crop + Flip
---|---|---|---
PSNR | 40.111 | 40.209 | 40.060
SSIM | 0.942 | 0.943 | 0.942

Fig. 11: Methods of cropping patches in stage three

We can conclude that method (2) achieves the best image restoration result,
this is different from our initial assumption. We originally assumed that
patches generated from method (3) can enable the model to learn the mapping
easier. However, the fact is that method (2) gets the better result. We
suppose this is because different patterns in method (2) play a crucial role
in data enhancement, thus prevent the model from overfitting.

###  3.2 Algorithm Comparison

In order to reflect the superiority of our algorithm, we have compared its
performance with the following four sorts of algorithms: (1) Analytical
reconstruction algorithm FBP; (2) Iterative reconstruction algorithm SART
combined with TV regularization; (3) Image inpainting, after reconstructing
the limited-view Radon data into images through FBP, apply the AAE model to
image restoration; (4) Sinogram inpainting, first use the AAE model to
complement the Radon data, and then adopt FBP to reconstruct it to images. We
also test these algorithms on two types of input data: (1) the directly cut
Radon data; (2) the fused Radon data. For Radon data with its post 60 degrees
being cut off, the performance of the above algorithms is shown in TABLE VI
and Fig.12, MR in this table means input the fused Radon data.

Table 6: Results of Different Algorithms Applied to Different Data Preprocessing Methods Algorithms |  PSNR |  SSIM
---|---|---
(1) FBP |  11.272 |  0.364
(2) FBP+MR |  12.354 |  0.452
(3) SART-TV |  14.727 |  0.635
(4) SART-TV+MR |  21.518 |  0.807
(5) Image Inpainting (II) |  35.566 |  0.916
(6) Image Inpainting + MR |  36.388 |  0.927
(7) Sinogram Inpainting (SI) |  27.345 |  0.800
(8) Sinogram Inpainting + MR |  28.960 |  0.859
(9) Ours |  40.209 |  0.943
Fig. 12: Histograms of different algorithms applied to different data preprocessing methods. Fig. 13: Visualized results of different algorithms applied to different data preprocessing methods.  Fig. 14: Error maps of different algorithms applied to different data preprocessing methods. Fig. 15: Histograms of different algorithms applied to different data preprocessing methods on different input data |  CUT-MID-60 |  CUT-MID-90 |  CUT-MID-120
---|---|---|---
Algorithms |  PSNR |  SSIM |  PSNR |  SSIM |  PSNR |  SSIM
(1) FBP |  11.131 |  0.362 |  10.350 |  0.289 |  9.636 |  0.217
(2) FBP + MR |  12.182 |  0.446 |  11.432 |  0.391 |  10.525 |  0.309
(3) SART-TV |  14.758 |  0.610 |  12.945 |  0.515 |  10.492 |  0.372
(4) SART-TV + MR |  21.036 |  0.784 |  17.523 |  0.722 |  13.166 |  0.592
(5) Image Inpainting (II) |  31.717 |  0.895 |  30.157 |  0.873 |  28.507 |  0.846
(6) Image Inpainting + MR |  32.031 |  0.895 |  30.422 |  0.876 |  28.999 |  0.849
(7) Sinogram Inpainting (SI) |  26.834 |  0.793 |  25.673 |  0.763 |  24.606 |  0.705
(8) Sinogram Inpainting + MR |  27.789 |  0.828 |  26.582 |  0.795 |  25.210 |  0.755
(9) Ours |  34.248 |  0.919 |  32.624 |  0.900 |  30.975 |  0.876
Table 7: Results of Different Algorithms Applied to Different Data
Preprocessing Methods on Different Input Data

From the results above, we can see that the idea of merging Radon data brings
additional prior information on every type of algorithm, thus improves their
performance by different margin. Besides, merging Radon data is particularly
helpful for SART-TV. Under the condition of using the same AAE model,
restoration in the image domain is more effective than restoration in the
Radon domain. Our algorithm combines the Radon domain and the image domain,
complements, reconstructs, restores and refines the input limited-view Radon
data, can finally improves the image PSNR to 40.209 and SSIM to 0.943. It
upgrades the quality of CT imaging by a large margin, realizes the accurate
restoration of its texture. Comparison of the visualized results can be seen
in Fig.13, we also present the corresponding error maps in Fig.14 to reflect
the difference in performance between different algorithms. To prove the
robustness of our algorithm, we test various degrees of damage on the input
Radon data, including (1) cut off the middle 60 degrees (1/3 of the original
data); (2) cut off the middle 90 degrees (1/2 of the original data); (3) cut
off the middle 120 degrees (2/3 of the original data). We implement the above
nine algorithms in TABLE VI on these three types of input data, and their
performance is shown in TABLE VII and Fig.15. It can be seen that losing data
in the middle can cause more damage than in the rear. With the increase of the
cropping ratio, the inpainting performance of these algorithms has also been
greatly affected. Our algorithm however, proves its outstanding robustness
under various conditions. Even when cutting the middle 120 degrees off the
original Radon data, our method can still restore the seriously damaged
imaging to PSNR of 30.975. Also, our method can exceed the other methods in
TABLE VII by a large margin under varying degrees of damaged data.

##  4 Conclusion

In order to improve the quality of the seriously damaged limited-view CT
imaging, we propose a three-stage restoration and reconstruction algorithm
based on spatial information, which combines the Radon domain and the image
domain, and utilizes the idea of “coarse to fine” to restore the image with
high definition. In the first stage, we designed an adversarial autoencoder to
complement the limited-view Radon data. In the second stage, we first
reconstruct the Radon data into images through FBP, and then send this image
into the Spatial-AAE model we built to achieve image artifact and noise
reduction based on spatial correlation between consecutive slices. In the
third stage, we propose the Refine-AAE network to finely repair the image in
patches, so as to achieve the accurate restoration of the image texture. For
Radon data with limited angle of 120 degrees (cut off one-third of the full-
view Radon data), our algorithm can increase its PSNR to 40.209, and SSIM to
0.943. At the same time, due to the fact that our model does not restrict
input resolution, can adapt to varying degrees of damage, and also can be
quickly implemented, our algorithm has generalization, robustness and
significant practical application value.

In our future work, we hope to incorporate our three-stage model into an end-
to-end network that can be simultaneously trained and tested. As we all know,
such large amount of parameters may be hard to optimize, we plan to solve this
problem by using tricks such as data augmentation and dropout, while also
lightweight model backbone like MobileNet.

## References

  * [1] H. Chen _et al._ , “Low-dose CT with a residual encoder-decoder convolutional neural network,” _IEEE Trans. Med. Imaging_ , vol. 36, no. 12, pp. 2524–2535, 2017.
  * [2] M. K. Kalra _et al._ , “Strategies for CT radiation dose optimization,” _Eur. J. Radiol._ , vol. 230, no. 3, pp. 619–628, 2004.
  * [3] T. L. Slovis, “The ALARA concept in pediatric CT: myth or reality?” _Eur. J. Radiol._ , vol. 223, no. 1, pp. 5–6, 2002.
  * [4] A. Khare and U. S. Tiwary, “A new method for deblurring and denoising of medical images using complex wavelet transform,” in _2005 IEEE Engineering in Medicine and Biology 27th Annual Conference_ , vol. 2, 2005, pp. 1897–1900.
  * [5] S. Xie _et al._ , “Artifact removal using improved GoogLeNet for sparse-view CT reconstruction,” _Sci. Rep._ , vol. 8, no. 1, pp. 1–9, 2018\.
  * [6] T. Zhang, H. Gao, Y. Xing, Z. Chen, and L. Zhang, “DualRes-UNet: Limited angle artifact reduction for computed tomography,” in _2019 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_. IEEE, 2019, pp. 1–3.
  * [7] S. Xie, H. Xu, and H. Li, “Artifact removal using gan network for limited-angle CT reconstruction,” in _2019 Ninth International Conference on Image Processing Theory, Tools and Applications (IPTA)_. IEEE, 2019, pp. 1–4.
  * [8] Q. Zhang, Z. Hu, C. Jiang, H. Zheng, Y. Ge, and D. Liang, “Artifact removal using a hybrid-domain convolutional neural network for limited-angle computed tomography imaging,” _Phys. Med. Biol._ , 2020.
  * [9] J. Xie, L. Xu, and E. Chen, “Image denoising and inpainting with deep neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 341–349.
  * [10] A. Katsevich, “Theoretically exact filtered backprojection-type inversion algorithm for spiral CT,” _SIAM J. Appl. Math._ , vol. 62, no. 6, pp. 2012–2026, 2002.
  * [11] K. Imai, M. Ikeda, Y. Enchi, and T. Niimi, “Statistical characteristics of streak artifacts on CT images: Relationship between streak artifacts and mA s values,” _Med. Phys._ , vol. 36, no. 2, pp. 492–499, 2009.
  * [12] L. Liu, “Model-based iterative reconstruction: a promising algorithm for today’s computed tomography imaging,” _Journal of Medical imaging and Radiation sciences_ , vol. 45, no. 2, pp. 131–136, 2014.
  * [13] G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” _Med. Phys._ , vol. 35, no. 2, pp. 660–663, 2008.
  * [14] G. H. Chen _et al._ , “Time-resolved interventional cardiac C-arm cone-beam CT: An application of the PICCS algorithm,” _IEEE Trans. Med. Imaging_ , vol. 31, no. 4, pp. 907–923, 2011.
  * [15] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” _Physica D_ , vol. 60, no. 1-4, pp. 259–268, 1992.
  * [16] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” _Phys. Med. Biol._ , vol. 53, no. 17, p. 4777, 2008.
  * [17] S. Niu _et al._ , “Sparse-view x-ray CT reconstruction via total generalized variation regularization,” _Phys. Med. Biol._ , vol. 59, no. 12, p. 2997, 2014.
  * [18] E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” _J. X-ray Sci. Technol._ , vol. 14, no. 2, pp. 119–139, 2006.
  * [19] A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” _Ultrasonic imaging_ , vol. 6, no. 1, pp. 81–94, 1984.
  * [20] Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang, “Low-dose X-ray CT reconstruction via dictionary learning,” _IEEE Trans. Med. Imaging_ , vol. 31, no. 9, pp. 1682–1697, 2012.
  * [21] M. Cao and Y. Xing, “Limited angle reconstruction with two dictionaries,” in _2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC)_. IEEE, 2013, pp. 1–4.
  * [22] H. Zhang, L. Zhang, Y. Sun, J. Zhang, and L. Chen, “Low dose CT image statistical iterative reconstruction algorithms based on off-line dictionary sparse representation,” _Optik_ , vol. 131, pp. 785–797, 2017.
  * [23] M. Xu, D. Hu, and W. Wu, “$\ell$0dl: Joint image gradient $\ell$0-norm with dictionary learning for limited-angle CT,” in _Proceedings of the 10th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics_ , 2019, pp. 538–538.
  * [24] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [25] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [26] R. K. Srivastava, K. Greff, and J. Schmidhuber, “Training very deep networks,” in _Advances in neural information processing systems_ , 2015, pp. 2377–2385.
  * [27] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 38, no. 2, pp. 295–307, 2015.
  * [28] T. Würfl, F. C. Ghesu, V. Christlein, and A. Maier, “Deep learning computed tomography,” in _International conference on medical image computing and computer-assisted intervention_. Springer, 2016, pp. 432–440.
  * [29] X. Mao, C. Shen, and Y.-B. Yang, “Image restoration using very deep convolutional encoder-decoder networks with symmetric skip connections,” in _Advances in neural information processing systems_ , 2016, pp. 2802–2810.
  * [30] H. Zhang _et al._ (2016) Image prediction for limited-angle tomography via deep learning with convolutional neural network. [Online]. Available: https://arxiv.org/abs/1607.08707
  * [31] E. Kang, J. Min, and J. C. Ye, “A deep convolutional neural network using directional wavelets for low-dose X-ray CT reconstruction,” _Med. Phys._ , vol. 44, no. 10, pp. e360–e375, 2017.
  * [32] Z. Zhang, X. Liang, X. Dong, Y. Xie, and G. Cao, “A sparse-view CT reconstruction method based on combination of DenseNet and deconvolution,” _IEEE Trans. Med. Imaging_.
  * [33] J. Wang, J. Liang, J. Cheng, Y. Guo, and L. Zeng, “Deep learning based image reconstruction algorithm for limited-angle translational computed tomography,” _PLoS One_ , vol. 15, no. 1, p. e0226963, 2020.
  * [34] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [35] Goodfellow _et al._ , “Generative adversarial nets,” in _Advances in neural information processing systems_ , 2014, pp. 2672–2680.
  * [36] R. Anirudh, H. Kim, J. J. Thiagarajan, A. K. Mohan, and K. M. Champley. (2019) Improving Limited Angle CT Reconstruction with a Robust GAN Prior. [Online]. Available: https://arxiv.org/abs/1910.01634
  * [37] Z. Li _et al._ , “A sinogram inpainting method based on generative adversarial network for limited-angle computed tomography,” in _15th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine_ , vol. 11072. International Society for Optics and Photonics, 2019, p. 1107220.
  * [38] Z. Li, A. Cai, L. Wang, W. Zhang, and B. Yan, “Promising generative adversarial network based sinogram inpainting method for ultra-limited-angle computed tomography imaging,” _IEEE Sensors J._ , vol. 19, no. 18, p. 3941, 2019.
  * [39] X. Dai, J. Bai, T. Liu, and L. Xie, “Limited-view cone-beam CT reconstruction based on an adversarial autoencoder network with joint loss,” _IEEE Access_ , vol. 7, pp. 7104–7116, 2018.
  * [40] R. Anirudh _et al._ , “Lose the views: Limited angle CT reconstruction via implicit sinogram completion,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 6343–6352.
  * [41] X. Dai, T. Liu, D. Hu, S. Yan, D. Shi, and H. Deng, “Limited angle cone-beam CT image reconstruction method based on geometric image moment,” 2016.
  * [42] L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” _Josa a_ , vol. 1, no. 6, pp. 612–619, 1984.
  * [43] K. Hammernik, T. Würfl, T. Pock, and A. Maier, “A deep learning architecture for limited-angle computed tomography reconstruction,” in _Bildverarbeitung für die Medizin 2017_. Springer, 2017, pp. 92–97.
  * [44] J. Zhao, Z. Chen, L. Zhang, and X. Jin. (2018) Unsupervised learnable sinogram inpainting network (SIN) for limited angle CT reconstruction. [Online]. Available: https://arxiv.org/abs/1811.03911
  * [45] Z. Zhao, Y. Sun, and P. Cong, “Sparse-view CT reconstruction via generative adversarial networks,” in _2018 IEEE Nuclear Science Symposium and Medical Imaging Conference Proceedings (NSS/MIC)_. IEEE, 2018, pp. 1–5.
  * [46] D. Lee, S. Choi, and H.-J. Kim, “High quality imaging from sparsely sampled computed tomography data with deep learning and wavelet transform in various domains,” _Med. Phys._ , vol. 46, no. 1, pp. 104–115, 2019.
  * [47] E. Zhou, H. Fan, Z. Cao, Y. Jiang, and Q. Yin, “Extensive facial landmark localization with coarse-to-fine convolutional network cascade,” in _Proceedings of the IEEE International Conference on Computer Vision (ICCV) Workshops_ , June 2013.
  * [48] M. Claus and J. van Gemert, “ViDeNN: Deep blind video denoising,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2019, pp. 0–0.
  * [49] M. Tassano, J. Delon, and T. Veit, “FastDVDnet: Towards real-time deep video denoising without flow estimation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 1354–1363.
  * [50] S. Armato _et al._ , “WE-B-201B-02: The lung image database consortium (LIDC) and image database resource initiative (IDRI): A completed public database of CT scans for lung nodule analysis,” _Med. Phys._ , vol. 37, no. 6Part6, pp. 3416–3417, 2010.
  * [51] D. P. Kingma and J. Ba. (2014) Adam: A method for stochastic optimization. [Online]. Available: https://arxiv.org/abs/1412.6980
  * [52] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” in _in ICML Workshop on Deep Learning for Audio, Speech and Language Processing_. Citeseer, 2013.
  * [53] C. Ledig _et al._ , “Photo-realistic single image super-resolution using a generative adversarial network,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 105–114.


# Optimal Control via Combined Inference and Numerical Optimization

Daniel Layeghi, Steve Tonneau and Michael Mistry All authors are with the
School of Informatics, University of Edinburgh. Email<EMAIL_ADDRESS>
###### Abstract

Derivative based optimization methods are efficient at solving optimal control
problems near local optima. However, their ability to converge halts when
derivative information vanishes. The inference approach to optimal control
does not have strict requirements on the objective landscape. However,
sampling, the primary tool for solving such problems, tends to be much slower
in computation time. We propose a new method that combines second order
methods with inference. We utilise the Kullback Leibler (KL) control framework
to formulate an inference problem that computes the optimal controls from an
adaptive distribution approximating the solution of the second order method.
Our method allows for combining simple convex and non convex cost functions.
This simplifies the process of cost function design and leverages the
strengths of both inference and second order optimization. We compare our
method to Model Predictive Path Integral (MPPI) and iterative Linear Quadratic
Regulator (iLQG), outperforming both in sample efficiency and quality on
manipulation and obstacle avoidance tasks.

##  I INTRODUCTION

The goal of minimising an objective is a general one that applies to many
areas. In the simplest cases the age-old solution is to solve for the
minimisation analytically. In many interesting cases though the analytical
method becomes intractable and we have to resort to computational methods.
Methods for computing the minimisation generally fall into two categories:
numerical and sampling based optimization. Both methods have advantages when
applied to the correct domain. Numerical methods are better suited to domains
where the objective of concern is smooth and derivative information is
available. On the other hand, if the objective function is highly
discontinuous and multi-modal, sampling is a better choice due to its inherent
exploitative capabilities and looser requirements on the objective.
Interesting optimization problems however are not only restricted to one of
the mentioned domains. As a result, it may seem desirable to combine both
methods to work in tandem. However, the abstraction that combines both methods
and formulation of a consensus between both methods is nontrivial. This paper
aims to naturally combine the use of numerical methods to refine when dense
derivative information is available and sample to explore when derivatives
vanish.

Here we focus on the optimal control setting and its application to robotic,
interesting settings where the governing dynamics are discontinuous and the
exact structure of the objective function is unknown. For example, locomotion,
is a task where the higher level objective is to transfer from point A to B.
The exact objective function that defines the intricate behaviour to achieve
the goal is unknown and the derivative of this cost with respect to the span
of controls is usually sparse. As a result making and breaking contacts to
solve the problem creates discontinuous optimization landscapes. One can
imagine the same principles in dexterous manipulation problems or conversely
problems where avoiding contacts is key, for example, obstacle avoidance.

##  II Related Work

One approach to remove the problem of sparse derivatives is to relax the
optimization problem. [1] show impressive results on locomotion and
manipulation tasks. The key insight is the relaxation of the contact
constraints. Resulting in dynamics where contact forces exist at distance. A
drawback of this approach is the problem of local minima. Once the
optimization landscape is sufficiently smoothed there may exist many local
optima where finding the feasible is hard. In addition physically unrealistic
solutions are also a problem. To mediate the problem of local optimum [2]
introduced an adaption with additional logical decision variables over action
modes such as sliding and flying. This allows for selective smoothing of
dynamics reducing the total number of local optima. However, the choice modes
create a limiting range of dynamics. In addition, the problem remains
deterministic leading to extremely intricate behaviours that are not robust.

To tackle the sparsity problem [3] introduces linear complementary constraints
to the optimization problem this allows the optimizer to decide whether bodies
are in contact. This formulation underpins many formulations around
locomotion. [4, 5, 6] use the same framework with additional task specific
parametrization. The formulation in this approach is specific to contact
dynamics and additionally the optimization generated is large and slow.

[7, 8] proposed a separate approach to reducing this problem of sparsity. They
decompose the problem of planning and control and use integer variables and
Mixed Integer Programming to select the dynamic modes of systems, for example,
in or out of contact. Subsequently a trajectory optimization problem is then
solved to follow the planned points. The major drawback of this method is the
loose relationship between the two solvers. As a result, significant work has
been done to introduce dynamic feasibility into the planning segment[9]. Such
feasibility constraints can be highly task-dependent and difficult to
formulate.

Apart from numerical methods. [10, 11] proposed an entirely different approach
inspired by Kalman’s duality [12] that treats the optimal control problem as
an inference. In this formalisation, the dynamics are treated as graphical
models and the optimal control distribution is solved through approximate
inference by using a variant of Monte Carlo methods. [13] extended this
approach by introducing Gaussian approximations to solve nonlinear problems.
[14] proposed a path integral (PI) formulation of this approach. [15] proposed
a PI inspired approach modified by insights from differential dynamic
programming (DDP) [16]. The inference approach results in desirable outcomes
such as clear formulations for exploration and the relaxation of the
requirement on the objective such as continuity and smoothness, therefore,
simpler cost function design for achieving goals. However, the major drawback
of such methods is one of sampling, i.e. inferring the unknown distribution is
computationally costly.

We aim to combine the strength of both inference and second order trajectory
optimization. We aim to formulate a natural combination that can efficiently
sample when no derivative information is available and follow optimized
trajectories when derivatives arise.

##  III The Optimal Control Setting

Let us consider the generic discrete optimal control setting. An agent at
state $\mathbf{x}\in\mathbf{X}$ is required to choose an action or control
$\mathbf{u}\in\mathbf{U}$ such that the resulting sequence of states and
actions
$\\{(\mathbf{x}_{0},\mathbf{u}_{0}),(\mathbf{x}_{1},\mathbf{u}_{1})...(\mathbf{x}_{n})\\}$
minimises the summation associated to the total cost of each state and action
pair $J(\mathbf{X},\mathbf{U})$. This problem can be solved by using the
principle of optimality. That is the optimal solution is the combined optimal
solutions to the subproblems of the original problem. This principle is
formalised in the Bellman equation [17] which underpins most of the optimal
control theory.

$V(\mathbf{x},t)=\min_{\mathbf{u}}{\left[l(\mathbf{x},\mathbf{u},t)+V\left({\text{f}(\mathbf{x},\mathbf{u},t)}\right)\right]}$

$V(\mathbf{x},t)$ represents the value function, $l(\mathbf{x},\mathbf{u})$
the running cost and $f(\mathbf{x},\mathbf{u})$ the deterministic system
dynamics all evaluated at $\mathbf{x}$ and $\mathbf{u}$. Solving this equation
in its vanilla form suffers from curse dimensionality; the computation grows
exponentially with the number of states. In the next sections, we focus on the
numerical and inference based solution to the problem.

###  III-A Numerical Optimization of The Optimal Control

The most established numerical solution to the Bellman solves the problem in
its tangent space using a second order method. This formulation is known as
Differential Dynamic Programming (DDP)[16]. Rewriting the Bellman equation as
a function of state and control and using $Q(x_{t},u_{t})$ to represent the
summation of running cost with the value function

$V(\mathbf{x},t)=\min_{\mathbf{u}}[Q(\mathbf{x}_{t},\mathbf{u}_{t})]$

Making a second order Taylor approximation leads to

$\displaystyle Q(\delta{\mathbf{x}},\delta{\mathbf{u}})$
$\displaystyle\approx\frac{1}{2}{\begin{bmatrix}1\\\ \delta{\mathbf{x}}\\\
\delta{\mathbf{u}}\end{bmatrix}}^{\intercal}\begin{bmatrix}0&Q_{\mathbf{x}}^{\intercal}&Q_{\mathbf{u}}^{\intercal}\\\
Q_{\mathbf{x}}&Q_{\mathbf{ux}}&Q_{\mathbf{ux}}\\\
Q_{\mathbf{u}}&Q_{\mathbf{ux}}&Q_{\mathbf{uu}}\end{bmatrix}\begin{bmatrix}1\\\
\delta{\mathbf{x}}\\\ \delta{\mathbf{u}}\end{bmatrix}$

Here the subscripts represents derivatives. Minimising the approximation with
respect to the perturbed control $\delta{\mathbf{u}}$ and solving for the
optimal control update

$\displaystyle
Q_{\mathbf{u}}+Q_{\mathbf{ux}}\delta{\mathbf{x}}+(\hat{\mathbf{u}}-\mathbf{u})Q_{\mathbf{uu}}=0$
$\displaystyle\hat{\mathbf{u}}=\mathbf{u}-\frac{Q_{u}Q_{\mathbf{ux}}}{Q_{\mathbf{uu}}}\delta{\mathbf{x}}$

The above is identical to Newton’s method where the
$Q_{\mathbf{u}}Q_{\mathbf{ux}}$ is the control gradient direction and
$Q_{\mathbf{uu}}$ is the Hessian with respect to control used as curvature
information. In DDP the update step is defined through a feed-forward gain
$\textbf{k}=\frac{Q_{\mathbf{u}}}{Q_{\mathbf{uu}}}$ and a feedback gain
$\textbf{K}=\frac{Q_{\mathbf{ux}}}{Q_{\mathbf{uu}}}$. [18] is a comprehensive
explanation of the Gauss-Newton version of DDP known as iLQG. In this case,
the second order derivatives with respect to dynamics are ignored and
additional regularization terms are introduced for the robustness of Hessian
approximations. From here onwards we refer to iLQG when referring to the
numerical solution to the optimal control.

###  III-B Approximate Inference of Optimal Control

The inference approach to the optimal control problem has been tackled under
many formulations. The KL control approach is grounded in the Bellman equation
therefore drawing comparisons with the numerical approach is simpler. As a
result, our focus in this paper is the KL-control approach and its dualities
with path integral control. One can view the KL control approach to be
grounded in the stochastic version of the Bellman equation shown in section
III-A

$\underset{\mathbf{u}}{\text{min}}\left[l(\mathbf{x})+D(\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})\;\|\;\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x}))+\operatorname{\mathbb{E}}_{\mathbf{x}^{{}^{\prime}}\sim\mathbf{u}}\left[V(\mathbf{x}^{\prime})\right]\right]$

Where $p(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})$ and
$p(\mathbf{x}^{\prime}|\mathbf{x})$ represent the controlled and uncontrolled
probabilities, $D$ represents a KL divergence operator. KL divergence measures
the information difference between two distributions. $V$ is the typical value
function. Under this formulation the probability for the optimally controlled
dynamics becomes

$\mathbf{p}^{*}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})=\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})\text{exp}\left(-V(\mathbf{x}^{\prime})\right)\eta^{-1}$

Therefore the optimal transition probability
$\mathbf{p}^{*}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})$ becomes the
transition probability of the passive dynamics adjusted by the exponential
values of states and normalised by $\eta$. [19, 20] show comprehensive
derivation and a Model Predictive Control (MPC) formulation that allows to
importance sample from this distribution. Our final algorithmic implementation
uses the overall MPC structure of this method shown in [20].

##  IV Combining Sampling With Second Order Method

To combine second order optimization and approximate inference, we start by
reformulating the Bellman equation and introduce a secondary divergence term
with the solution of the Bellman equation in the tangent space. That is the
optimal trajectory obtained by DDP/iLQG algorithm. A core assumption here is
that the input $\mathbf{u}$ is assumed to be the mean of an input sampled from
a normal distribution, meaning the actual input
$v\sim\mathcal{N}(\mathbf{u},\mathbf{\Sigma})$. Modifying the KL objective we
can write

$\displaystyle\underset{\mathbf{u}}{\text{min}}[l(\mathbf{x})+(1-k)D(\mathbb{Q}\;\|\;\mathbb{P})+kD(\mathbb{Q}\;\|\;\mathbb{C})+\operatorname{\mathbb{E}}_{\mathbb{Q}}[V(\mathbf{x}^{\prime})]]$
(1) $\displaystyle
D(\mathbb{Q}\;\|\;\mathbb{P})=\operatorname{\mathbb{E}}_{\mathbb{Q}}\left[\log\left[\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})}{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})}\right]\right]$
$\displaystyle
D(\mathbb{Q}\;\|\;\mathbb{C})=\operatorname{\mathbb{E}}_{\mathbb{Q}}\left[\log\left[\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})}{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})}\right]\right]$

where $\mathbb{P}$ is the distribution of the passive dynamics for which the
input $v$ has a mean of 0 with variance $\mathbf{\Sigma}$. $\mathbb{C}$
represents the distribution of the controlled dynamics under iLQG where the
input has a mean of $\mathbf{u}_{\text{iLQG}}$ with variance
$\mathbf{\Sigma}_{\text{iLQG}}$. $\mathbb{Q}$ represents our control
distribution which has a mean of $u$ and variance $\Sigma$ under the input
$v$. $k\in[0,1]$ is the importance given to minimising each KL.

Combining expectations and rearranging leads to the optimal control trajectory
(complete derivation in the Appendix A-A)

$\displaystyle\mathbf{q}^{*}(V)=\mathbf{p}(V)^{1-k}\mathbf{c}(V)^{k}\text{exp}\big{(}-\frac{1}{\lambda}S(V)\big{)}\eta^{-1}$
(2)

Where $\lambda$ is a hyperparameter discounting trajectory costs. Equation 2
is very similar to the optimal transition described in section III-B. However
in this case the optimal trajectory is a factor of the joint distributions
$\mathbb{P}$ and $\mathbb{C}$ scaled by the importance $k$ and adjusted by the
exponentiated negative cost of trajectories $\text{exp}(S(-V))$. In similar
fashion to [20] we use equation 2 to importance sample from the optimal
distribution $\mathbb{Q^{*}}$

$\mathbf{u}^{*}=\bigintsss\mathbf{q}(V)\underbrace{\frac{\mathbf{q}^{*}(V)}{\mathbf{p}(V)^{1-k}\mathbf{c}(V)^{k}}\frac{\mathbf{p}(V)^{1-k}\mathbf{c}(V)^{k}}{\mathbf{q}(V)}}_{w(V)}v_{t}\mathrm{d}V$
(3)

Figure 1: Left to right: Cartpole scenario: The objective is to swing the pole
to upwards position. Finger-Spinner: The objective is for the fully actuated 2
link manipulator to rotate the object to the desired location marked as
transparent. Push Object: The goal is for the 3-link manipulator to move the
orange sphere to the desired location (green). Obstacle Avoidance: The goal is
for the 3 link manipulator to move to the desired position (red) without
making contact with spherical obstacles. Results video at
https://youtu.be/dZgRJ0CTpLw Table I: Total trajectory cost and success rate
for each method per environment. The costs are computed across length of
trajectories until completion of tasks. Column Time represents the total
computation time until task completion in seconds. For number of samples and
hyperparameters refer to each task.

Environment % Ours MPPI Naive Combination iLQG Cost % Success Time Cost %
Success Time Cost % Success Time Cost % Success Time Cartpole $6.691003$ 100
15.05 $8.201003$ 100 10.03 $2.07104$ 100 33.13 $2.5103$ 100 1.84 Finger-
Spinner $1.45108$ 100 7.19 $1.97108$ 100 7.68 $3.52108$ 100 10.56 N/A 0 N/A
Push Object $2.121011$ 100 15.87 $4.921011$ 100 73.93 N/A 0 N/A N/A 0 N/A
Obstacle Avoidance $6.50105$ 100 12.95 $8.79105$ 25 11.376 N/A 0 N/A N/A 0 N/A

Computing the importance sampling weights $w(V)$ with respect to the
distributions (see Appendix A-A) for a trajectory in length of time T.

$\displaystyle w(V)=\frac{1}{\eta}\exp\Big{(}$
$\displaystyle-\frac{1}{2}\frac{1}{\lambda}S(V)$ (4)
$\displaystyle+\sum_{t=0}^{T}(2{\mathbf{v}_{t}}^{T}\mathbf{\Sigma}^{-1}\mathbf{u}_{t}-{\mathbf{u}_{t}}^{T}\mathbf{\Sigma}^{-1}\mathbf{u}_{t}$
$\displaystyle+k\big{(}{(\mathbf{v}_{t}-\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\
t})}^{T}\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\
t}^{-1}(\mathbf{v}_{t}-\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\ t})$
$\displaystyle-{\mathbf{v}_{t}}^{T}\mathbf{\Sigma}^{-1}\mathbf{v}_{t}\big{)}\Big{)}$

From equation 4 $D(\mathbb{Q}\;\|\;\mathbb{C})$ translates into a cost term
within $w(V)$ that penalises the deviation of the sampled control
$\mathbf{v}_{t}$ from the solution obtained by iLQG
$\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\ t}$ where the cost gain is the
covariance of the iLQG solution
$\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\ t}$. Additionally, the
parameter $k$ dictates the importance of this cost and if set to 0 the weights
become the same as the KL control case.
$\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\ t}$ is a proxy measure of
our confidence in the solution provided by the second order method. To compute
the variance of $\mathbb{C}$, the distribution that represents the solution of
the second order method, we know that the trajectory
$\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\ t}$ minimises the quadratic
approximation of the instance of the Q function $\delta
Q(\mathbf{x}_{t},\mathbf{u}_{t})$ at
$\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\ t}$. As a result, the mean and
the variance for the normal distribution $\mathbb{C}$ can be defined by
minimizing the log-likelihood:

$\displaystyle\ln P(1|\mathbf{u}_{t},\mathbf{\Sigma}_{t})=\ln\exp(-\delta
Q_{t}(\mathbf{x}_{t},\mathbf{u}_{t}))$ (5)

minimising the log-likelihood leads to the same problem shown in section
III-A. Therefore, $\frac{\delta}{\delta\mathbf{u}}\ln
P(1|\mathbf{u}_{t},\mathbf{\Sigma}_{t})$ is at a minimum if
$\mathbf{u}_{t}=\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\ t}$. We can obtain
the variance by computing the second derivative,
$\frac{\delta^{2}}{\delta\mathbf{u}^{2}}\ln
P(1|\mathbf{u}_{t},\mathbf{\Sigma}_{t})=Q_{t\leavevmode\nobreak\
\mathbf{uu}}^{-1}=\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\ t}$. As a
result input trajectories sampled form $\mathbb{C}$ take the form
$v_{t}\sim\mathcal{N}(\mathbf{u}_{\text{iLQG}\leavevmode\nobreak\
t},Q_{\mathbf{uu}}^{-1}\leavevmode\nobreak\ t)$. In our final formulation we
scale this variance approximation with a hyperparameter $\beta$.

Algorithm 1 shows the final structure of our proposed method.

1$K$: Number of samples

2 $T$: Number of time-steps

3 $q_{\text{iLQG}\leavevmode\nobreak\ r}(\textbf{x}_{t})$ and
$q_{r}(\textbf{x}_{t})$: Running costs

4 $q_{\text{iLQG}\leavevmode\nobreak\ t}(\textbf{x}_{T})$ and
$q_{t}(\textbf{x}_{T})$: Terminal costs

5 $\textbf{u}_{0:T-1}$: Initial trajectory

6 $\textbf{u}_{\text{iLQG}\leavevmode\nobreak\ 0:T-1}$: Initial trajectory

7

8$\lambda,\leavevmode\nobreak\ k,\leavevmode\nobreak\
\beta,\leavevmode\nobreak\ \mathbf{\Sigma}$: Hyper parameters for eq. 4

9 while  _not done_ do

10 $\mathbf{x}_{0}$ = EstimateState()

11 $[\textbf{u}_{\text{iLQG}\leavevmode\nobreak\ 0:T-1},\leavevmode\nobreak\
Q_{\mathbf{uu}\leavevmode\nobreak\ 0:T-1}]$ = iLQG($\textbf{x}_{0}$)

12 $\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\
0:T-1}=\beta\leavevmode\nobreak\ Q^{-1}_{\mathbf{uu}\leavevmode\nobreak\
0:T-1}$

13 $\textbf{u}_{0:T-1}$ =
KLControl($\textbf{u}_{\text{iLQG}\leavevmode\nobreak\
0:T-1},\leavevmode\nobreak\ \mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\
0:T-1},\leavevmode\nobreak\ \textbf{x}_{0}$) eq. 4

14 apply $\textbf{u}_{0}$

15 $\textbf{u}_{0:T-1}=\textbf{u}_{1:T-1}$

16 $\textbf{u}_{\text{iLQG}\leavevmode\nobreak\ 0:T-1}=\textbf{u}_{0:T-1}$

17 end while

Algorithm 1 Inference with Trajectory Optimization (Ours)

##  V Experimental Results and Discussion

The results obtained were generated from environments that were designed to
test the validity and benefits of our method. Therefore, key attributes of
these environments shown in figure 1 are

  * •

Discontinuous dynamics.

  * •

Tasks where making or avoiding contact is necessary for completion.

  * •

Tasks that allows for designing simple non convex cost functions that induce
sparsity in derivative information.

In this section, we evaluate whether our algorithm can combine trajectory
optimization and inference in a manner that allows for following the optimized
trajectory when derivative information is available and resorting to unbiased
inference when derivatives vanish. We compare our method to 3 other
algorithms: MPPI, iLQG and the naive combination of both where the output of
iLQG is used as an initial trajectory for MPPI and vice versa. The results
obtained are obtained across 20 random seeds. For each task, the lengths of
the trajectories between solvers are the same. Table I is a summary of these
results.

###  V-A Implementation details

The implementation was entirely developed in C++. The environments were
simulated using Mujoco 2.0 [21]. We have tuned the contact dynamics in Mujoco
to be realistic, additionally, our iLQG implementation uses derivative
information at contacts. Some of the environments were modified versions of
models developed by [22]. The results were obtained on a PC with an Intel Core
i9 CPU running Ubuntu 18.04.

###  V-B Cartpole task

The Cartpole is a task we use to evaluate and compare the performance of our
method in its simplest case. In this task the state
$\mathbf{x}\in\mathbb{R}^{4}$ consists of the cart and pole position and
velocities. The desired state is $\mathbf{x}_{\text{goal}}=(0,0,0,0)$. The
running cost for this task is $(\textbf{x}_{\text{error}\leavevmode\nobreak\
t})\mathbf{Q}_{t}(\textbf{x}_{\text{error}\leavevmode\nobreak\ t})$ and a
terminal cost of $(\textbf{x}_{\text{error}\leavevmode\nobreak\
T})\mathbf{Q}_{T}(\textbf{x}_{\text{error}\leavevmode\nobreak\ T})$ where
$\mathbf{Q}_{t}=\text{diag}($1103$,$5102$,0,0)$ and
$\mathbf{Q}_{T}=\text{diag}($1105$,$5104$,$5102$,$5102$)$ with horizon
$T=0.75s$ and number of samples $K=3$. The hyperparameters $\beta$, $\lambda$
and $k$ were set to $110-3$, 0.1 and 1 respectively.

The results show that iLQG outperforms all other methods, beating ours by a
factor of 2.7. This is not surprising as the cost function is not sparse and
dynamics are continuous leading to a very smooth optimization landscape where
any derivative based method will do well. Our method however outperforms MPPI
by a factor of 1.2. This is achieved with a very low number of samples and
biasing the samples towards the iLQG trajectory by setting a low value of
$\beta$.

###  V-C Finger-Spinner task

This is a manipulation task where the states $\mathbf{x}\in\mathbb{R}^{6}$ are
the joint positions of the manipulator and the passive spinner followed by
their velocities. The desired goal state is
$\mathbf{x}_{\text{goal}}=(0,0,0,0,0,0)$. The cost function is quadratic with
respect to the states where the gain terms act on the state of the spinner.
$\mathbf{Q}_{t}=\text{diag}(0,0,$1.5104$,0,0,50)$ and
$\mathbf{Q}_{T}=\text{diag}(0,0,$1.5104$,0,0,50)$. An additional binary
scaling term
$\mathds{1}(\textbf{x}_{t})\textbf{x}_{\text{error}\leavevmode\nobreak\ t}$ is
also added to the running cost, where

$\mathds{1}(\textbf{x}_{t})=\left\\{\begin{array}[]{@{}ll@{}}0,&\text{if}\leavevmode\nobreak\
\text{contact between spinner and finger}\\\
1,&\text{otherwise.}\end{array}\right.$

This term tracks the global positions of bodies and encourages contact until
low state errors. Since this term is non differentiable, it is only added to
the inference cost in equation 4. The choice of the hyperparameters are
$[\beta,\lambda,k]=[0.25,0.5,1]$ and $T=0.75$ and $K_{\text{MPPI}}=30$ and
$K_{\text{Ours}}=10$

Our method outperforms all other methods, beating the second best (MPPI) by a
factor of 1.3 with a third of the number of samples and faster compute cycle.
The derivative based method (iLQG) entirely fails to solve this task. This is
because the cost function is sparse and the dynamics are highly discontinuous.
As a result, very few control inputs will create the contacts and the
derivative direction required to solve the task. The naive combination is
capable of solving this problem in effect by randomising the control input
until contact is made. This results in a suboptimal trajectory, which is due
to the very loose coupling between inference and numerical optimization.

###  V-D Push Object task

This is another manipulation task similar to the previous. In this task the
object to be manipulated has higher degrees of freedom and the exploration
domain for finding contact is much larger. The states
$\mathbf{x}\in\mathbb{R}^{10}$ are the joint positions of the planar
manipulator in addition to the global Cartesian position of the object that is
manipulated followed by their derivatives. The desired state is
$\mathbf{x}_{\text{goal}}=(0,0,0,0.25,-0.22,0.022,0,0,0,0)$. The cost function
is quadratic with respect to the states and only defined by a term on the
state of the object where
$\mathbf{Q}_{t}=\text{diag}(0,0,0,$1105$,$1105$,0,0,0,$5102$,$5102$)$ and
$\mathbf{Q}_{T}=(0,0,0,$1105$,$1105$,0,0,0,0,0)$. To encourage contact similar
to the previous case, an additional binary term
$\mathds{1}(\textbf{x}_{t})\textbf{x}_{\text{error}\leavevmode\nobreak\ t}$ is
added to the running cost in equation 4. The choice of the hyperparameters are
$[\beta,\lambda,k]=[1,0.25,1]$ and $T=0.75$, $K_{\text{MPPI}}=80$ and
$K_{\text{Ours}}=50$

Our method performs better than all other methods. Beating the second best
(MPPI) by a factor 2.3. Our result show that contact dynamics is the major
contributing factor to this. Our method is capable of pushing the object much
more intricately by staying close to the solution of the second order method.
iLQG on its own is not capable of solving this task due to the same sparsity
of derivatives and discontinuity reasons as the previous task. The naive
combination also fails 25% of the time and otherwise computes a trajectory
that is an order of magnitude worse.

###  V-E Obstacle Avoidance task

This task is essentially the opposite of the manipulation task. For successful
completion contact with environment has to be avoided. The states
$\mathbf{x}\in\mathbb{R}^{6}$ are the joint positions and velocities. The
desired state is $\mathbf{x}_{\text{goal}}=(9.72,0,0,0,0,0)$. The cost
function is quadratic with respect to the states, defined by a term on the
state of the manipulator where $\mathbf{Q}_{t}=\text{diag}(100,10,10,0,0,0)$
and $\mathbf{Q}_{T}=(100,100,100,10,1,1)$. To discourage contact, a binary
term penalising contacts
$\mathds{1}(\textbf{x}_{t})\textbf{x}_{\text{error}\leavevmode\nobreak\ t}$ is
added to the running cost in equation 4. The choice of the hyperparameters are
$[\beta,\lambda,k]=[20,0.1,1]$ and $T=0.75$ , $K_{\text{MPPI}}=100$ and
$K_{\text{Ours}}=50$

Similar to the previous case we outperform all other algorithms, on average
beating MPPI by a factor of 1.4 when compared to its best case. We are able to
achieve 100% success rate in comparison to the 25% success rate of MPPI with
half the number of samples. By choosing high values of $\lambda$ our algorithm
performs complex maneuvers to escape local optima produced by the obstacles
where as MPPI is not capable of generating such motions. iLQG and the naive
combination are not able to avoid obstacles and are not capable of solving
this task.

##  VI Discussion

The relationship between the covariance and the inverse Hessian ($\beta
Q_{t\leavevmode\nobreak\
\mathbf{uu}}^{-1}=\mathbf{\Sigma}_{\text{iLQG}\leavevmode\nobreak\ t}$) plays
a central role in our method. We use this information to define a confidence
measure. However, this measure simply reflects the rate of convergence and not
the quality of the solution. This poses problems where derivative based
optimization will pay high costs to obtain fast convergence.

The results show that our method is capable of finding derivative dense
locations and following optimized trajectories. This notion is tightly coupled
with the choice of hyperparameters, specifically $\lambda$ and $\beta$. This
choice becomes apparent for environments with flat optimization landscapes
such as the ones related to the tasks shown in section V-D. These cases
require exploration through high values of $\lambda$ to allow for the sparse
spread of the probability mass across trajectories. This however results in
discounting the KL divergence cost which can be detrimental to high confidence
optimized trajectories. Therefore choosing the hyperparameters in our method
is not trivial and can result in unsuccessful trajectories.

The ability to introduce non differentiable costs to the inference part of our
algorithm introduces a desirable formulation where we can remove the need to
define complex convex costs for the derivative based part of our algorithm.
For example, in section V-E we remove the requirement for a distance cost to
avoid obstacles. Another benefit of this formulation is the ability to
encourage contacts with similar cost structures and removing requirements such
as defining contact locations or introducing relaxations to guide the
derivative based method towards contacts.

The results that we have obtained are for cases where the mixing term $k$
between the passive and the controlled distribution is set to 1. This is
because $k$ induces a weighted joint distribution between the probabilities
$\mathbb{P}$ and $\mathbb{C}$. This can become problematic when the distance
between the two distributions is large. In that case, the joint can exist in a
domain where no useful control inputs exist. An automatic adaption of the
mixing term can prove useful when the control distribution is uninformative.

With regards to the computation time we show that we are able to outperform
MPPI for our tasks. This is because of the sample efficiency of our algorithm
through the combination with the second order method. This is especially
apparent in cases where contact and reasoning about their dynamics is
required. In such scenarios the second order method can provide more
information around completion of the task. Our push object task shows this. We
perform on similar wall clock times when compared to the naive implementation.
This is expected as complexity of both algorithms are the same.

##  VII Conclusion and Future Work

We present a method that combines approximate inference with second order
trajectory optimization. We formulate the approximate inference through a
variant of KL control that introduces a divergence term from the distribution
of optimized trajectories. Through our derivation, we show that this
distribution can be approximated by the solution of the trajectory
optimization and its inverse Hessian. This natural combination removes the
design of complex convex cost functions. We demonstrate this effect on our
obstacle avoidance and manipulation tasks where we only use binary cost
functions to encourage or discourage making contacts. We test our methods on
environments with discontinuities and a mix of simple differentiable and non
differentiable cost functions that nonetheless induce sparsity in derivative
information. We compare our results to ones obtained from vanilla iLQG, iLQG
with randomized warmstarts and MPPI. We show that in cases where derivative
information is not immediately present, our method outperforms all others by
exploring to find low cost trajectories and following the optimized
trajectories when derivative information appears. Some avenues can extend this
work. The hyperparameters for temperature, scaling the Hessian and the mixing
term between two distributions are central to efficient exploration and
refinement. A systematic approach for an adaption of these parameters can
bring efficiency gains. The KL control formulation of the approximate
inference is based on the assumption of constant input noise. Introducing
theoretically sound variance adaption strategies can provide better
convergence especially when combined with derivative based optimization.

##  Appendix A Appendix

###  A-A KL control derivation

To obtain the final optimal distribution obtained in section IV we start by

$\displaystyle\underset{\mathbf{u}}{\text{min}}[l(\mathbf{x})+(1-k)D(\mathbb{Q}\;\|\;\mathbb{P})+kD(\mathbb{Q}\;\|\;\mathbb{C})+\operatorname{\mathbb{E}}_{\mathbb{Q}}[V(\mathbf{x}^{\prime})]]$
$\displaystyle
D(\mathbb{Q}\;\|\;\mathbb{P})=\operatorname{\mathbb{E}}_{\mathbb{Q}}\left[\log\left[\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})}{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})}\right]\right]$
$\displaystyle
D(\mathbb{Q}\;\|\;\mathbb{C})=\operatorname{\mathbb{E}}_{\mathbb{Q}}\left[\log\left[\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})}{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})}\right]\right]$

Combining expectations give (only minimisable terms of interest shown)

$\displaystyle\underset{\mathbf{u}}{\text{min}}\left[\operatorname{\mathbb{E}}_{\mathbb{Q}}\left[\log\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})}{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})^{(1-k)}\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})^{k}\exp(-V(\mathbf{x}^{\prime}))}\right]\right]$

The expatiated negative values are not are not a transition probabilities
therefore introducing a normalisation constant

$\displaystyle\phi(\mathbf{x})=\sum\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})^{(1-k)}\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})^{k}\exp(-V(\mathbf{x}^{\prime}))$

Adding the normalization constant to the original problem and assuming
$z(\mathbf{x^{\prime}})=\exp(-V(\mathbf{x^{\prime}}))$we can write

$\displaystyle\underset{\mathbf{u}}{\text{min}}\left[D(\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})\;\|\;\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})^{(1-k)}\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})^{k}z(\mathbf{x^{\prime}})}{\phi(\mathbf{x})})\right]$

From the KL is minimised when the numerator and the denominator are equal:

$\displaystyle\mathbf{p}^{*}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u})=\frac{\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x})^{(1-k)}\mathbf{p}(\mathbf{x}^{\prime}|\mathbf{x},\mathbf{u}_{\text{iLQG}})^{k}\exp(-V(\mathbf{x}^{\prime}))}{\phi(\mathbf{x})}$

Writing the above transitions probabilities in terms of their density
functions over a trajectory of length $T$ with control input $V$

$\displaystyle\mathbf{p}(V)={}$ $\displaystyle
Z^{-1}\exp\left(-\frac{1}{2}\sum_{t=0}^{T_{1}}\mathbf{v}_{t}^{T}\mathbf{\Sigma}^{-1}\mathbf{v}_{t}\right)$
$\displaystyle\mathbf{q}(V)={}$ $\displaystyle
Z^{-1}\exp\left(-\frac{1}{2}\sum_{t=0}^{T_{1}}(\mathbf{v}_{t}-\mathbf{u}_{t})^{T}\mathbf{\Sigma}^{-1}(\mathbf{v}_{t}-\mathbf{u}_{t})\right)$
$\displaystyle\mathbf{c}(V)={}$ $\displaystyle Z^{-1}_{\text{iLQG}}$
$\displaystyle\exp\left(-\frac{1}{2}\sum_{t=0}^{T_{1}}(\mathbf{v}_{t}-\mathbf{u}_{\text{iLQG}_{t}})^{T}\mathbf{\Sigma}^{-1}_{\text{iLQG}_{t}}(\mathbf{v}_{t}-\mathbf{u}_{\text{iLQG}_{t}})\right)$

As a result the final optimal distribution becomes

$\displaystyle\mathbf{q}^{*}(V)$
$\displaystyle=\mathbf{p}(V)^{1-k}\mathbf{c}(V)^{k}\text{exp}\big{(}-\frac{1}{\lambda}S(V)\big{)}\eta^{-1}$

where $\eta$ is the normalisation constant over trajectories and $\lambda$ is
the temperature an adjustable hyperparameter.

## References

  * [1] I. Mordatch _et al._ , “Discovery of complex behaviors through contact-invariant optimization,” _ACM Transactions on Graphics (TOG)_ , vol. 31, no. 4, pp. 1–8, 2012.
  * [2] M. Toussaint _et al._ , “Describing physics for physical reasoning: Force-based sequential manipulation planning,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6209–6216, 2020.
  * [3] M. Posa _et al._ , “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [4] D. Pardo _et al._ , “Hybrid direct collocation and control in the constraint-consistent subspace for dynamic legged robot locomotion.” in _Robotics: Science and Systems_ , 2017.
  * [5] A. W. Winkler _et al._ , “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, 2018.
  * [6] O. Cebe _et al._ , “Online dynamic trajectory optimization and control for a quadruped robot,” _arXiv preprint arXiv:2008.12687_ , 2020.
  * [7] B. Aceituno-Cabezas _et al._ , “Simultaneous contact, gait, and motion planning for robust multilegged locomotion via mixed-integer convex optimization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2531–2538, 2017.
  * [8] S. Kuindersma _et al._ , “Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot,” _Autonomous robots_ , vol. 40, no. 3, pp. 429–455, 2016.
  * [9] S. Tonneau _et al._ , “Sl1m: Sparse l1-norm minimization for contact planning on uneven terrain,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6604–6610.
  * [10] H. J. Kappen _et al._ , “Optimal control as a graphical model inference problem,” _Machine learning_ , vol. 87, no. 2, pp. 159–182, 2012.
  * [11] E. Todorov, “Linearly-solvable markov decision problems,” in _Advances in neural information processing systems_ , 2007, pp. 1369–1376.
  * [12] R. F. Stengel, _Stochastic optimal control: theory and application_. John Wiley & Sons, Inc., 1986.
  * [13] M. Toussaint, “Robot trajectory optimization using approximate inference,” in _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 1049–1056.
  * [14] E. Theodorou _et al._ , “A generalized path integral control approach to reinforcement learning,” _The Journal of Machine Learning Research_ , vol. 11, pp. 3137–3181, 2010.
  * [15] J. Rajamäki _et al._ , “Sampled differential dynamic programming,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 1402–1409.
  * [16] D. H. Jacobson, “New second-order and first-order algorithms for determining optimal control: A differential dynamic programming approach,” _Journal of Optimization Theory and Applications_ , vol. 2, no. 6, pp. 411–440, 1968.
  * [17] R. Bellman, “Dynamic programming,” _Science_ , vol. 153, no. 3731, pp. 34–37, 1966.
  * [18] Y. Tassa _et al._ , “Synthesis and stabilization of complex behaviors through online trajectory optimization,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 4906–4913.
  * [19] E. A. Theodorou and E. Todorov, “Relative entropy and free energy dualities: Connections to path integral and kl control,” in _2012 ieee 51st ieee conference on decision and control (cdc)_. IEEE, 2012, pp. 1466–1473.
  * [20] G. Williams _et al._ , “Information theoretic mpc for model-based reinforcement learning,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 1714–1721.
  * [21] E. Todorov _et al._ , “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 5026–5033.
  * [22] Y. Tassa _et al._ , “Deepmind control suite,” _arXiv preprint arXiv:1801.00690_ , 2018.