supplement.pdf

# SUD2: Supervision by Denoising Diffusion Models for Image Reconstruction

Matthew A. Chan
University of Maryland, College Park
<EMAIL_ADDRESS>Sean I. Young
Massachusetts Institute of Technology
<EMAIL_ADDRESS>Christopher A. Metzler
University of Maryland, College Park
<EMAIL_ADDRESS>
###### Abstract

Many imaging inverse problems—such as image-dependent in-painting and
dehazing—are challenging because their forward models are unknown or depend on
unknown latent parameters. While one can solve such problems by training a
neural network with vast quantities of paired training data, such paired
training data is often unavailable. In this paper, we propose a generalized
framework for training image reconstruction networks when paired training data
is scarce. In particular, we demonstrate the ability of image denoising
algorithms and, by extension, denoising diffusion models to supervise network
training in the absence of paired training data.

##  1 Introduction

Imaging inverse problems can generally be described in terms of a forward
operator $\mathcal{F}(\cdot)$ that maps a scene $x$ to a measurement $y$
according to

$y=\mathcal{F}(x).$ (1)

The goal of an image reconstruction algorithm is to recover $x$ from $y$.

Historically, computational imaging research has focused on solving inverse
problems with known forward models. For instance, computed tomography’s
forward model can be represented as a Radon transform, and magnetic resonance
imaging’s forward model can be represented as 2D Fourier Transform. Knowledge
of these forward models allows one to reconstruct scenes $x$ from measurements
$y$ using any number of classical or learning-based algorithms [18].

Since the onset of the deep learning era, significant progress has been made
in solving inverse problems which lack explicit forward models. By leveraging
large amounts of training pairs $\\{x_{i},y_{i}\\}^{N}_{i=0}$, neural networks
learn to directly map samples from a source distribution, $y_{i}$, to images
from a target distribution, $x_{i}$. In doing so, the network implicitly
learns the inverse operator $\mathcal{F}^{-1}$ without any explicit knowledge
of the forward model $\mathcal{F}$.

The main drawback of deep learning methods is that their performance is
directly related to the size and quality of the training dataset. As a result,
these methods often struggle whenever little to no paired training data is
available.

###  1.1 Problem setup

Our goal in this work is to train a network $f_{\theta}(\cdot)$ to reconstruct
images/scenes $x$ from measurements $y$ using three sets of training data.

  * •

A small set $P$ of paired examples $(x_{p},y_{p})$ drawn from the joint
distribution $p_{x,y}$.

  * •

A large set $U_{y}$ of unpaired measurements $y_{u}$ drawn from the marginal
distribution $p_{y}$.

  * •

A large set $U_{x}$ of unpaired images $x_{u}$ drawn from the marginal
distribution $p_{x}$.

Such mixed datasets naturally occur in applications where gathering unpaired
data is easy, but gathering paired data is a challenge. For instance, it is
straightforward to capture images with fog and images without fog, but
capturing two paired images of the same scene with and without fog (with all
lighting conditions and other nuisance variations fixed) is very challenging.
Often times, the latter paired dataset is restricted to only a few images
captured in a lab.

A paired training set, $P$, allows one to optimize $f_{\theta}(\cdot)$ by
minimizing the empirical risk

$\displaystyle\mathcal{L}_{\text{paired}}=\frac{1}{|P|}\sum_{(x_{p},y_{p})\in
P}\|x_{p}-f_{\theta}(y_{p})\|^{2},$ (2)

where $|P|$ denotes the cardinality of $P$. However, as the size of $P$
decreases $\mathcal{L}_{\text{paired}}$ becomes a poor approximation of the
true risk and $f_{\theta}(\cdot)$ overfits to the training set. As an
alternative, we seek to leverage unpaired datasets $U_{x}$ and $U_{y}$ to
improve the quality of our reconstructions.

###  1.2 Our contributions

In this work we introduce a novel semi-supervised learning framework for
addressing image reconstruction problems, such as image dehazing, which lack
an explicit forward model and for which gathering paired training data is
challenging. Our central contributions are as follows;

  * •

We generalize the supervision-by-denoising (SUD) semi-supervised learning
technique from [28] to work on any image reconstruction task, not just medical
image segmentation.

  * •

We prove that SUD implicitly performs cross-entropy minimization and use this
connection to identify its various failure modes.

  * •

Based on our analysis, we identify three techniques to improve SUD: (1) sample
correlation minimization, (2) noise injection, and (3) denoising diffusion
models.

  * •

We demonstrate that the resulting algorithm, which we call SUD2, outperforms
existing semi-supervised and unsupervised learning techniques, such as
CycleGAN, on image in-painting and image dehazing.

##  2 Related work

###  2.1 Semi-supervised learning

While crowdsourcing platforms such as Amazon’s Mechanical Turk [23] can
generate abundant labeled data for human-annotatable tasks like image
segmentation, it’s all but impossible to hand-label training data for many
imaging inverse problems. Given a foggy image, how would one generate
corresponding clean image to use for training? Even in applications where
hand-labeling is possible, it is often prohibitively expensive to perform at
scale.

Semi-supervised learning (SSL) serves as a workaround for applying deep
learning techniques in situations where paired data is scarce. Typically, SSL
methods regularize network training by leveraging information extracted from
unpaired data and generally fall into two distinct categories; pseudo-
labelling methods and consistency regularization methods [27]. Pseudo-
labelling methods [13] help supervise training by generating fake labels on
unpaired data samples, which allows the training of networks in a fully-
supervised manner. In contrast, consistency regularization methods like
temporal ensembling [11] and mean teacher models [24] enforce a regularization
objective during training that makes the network more robust to perturbations
in the data.

Our method borrows ideas from both categories of SSL methods. Similar to
pseudo-labelling methods, we generate fake labels on unpaired data by
incorporating a denoising diffusion probabilistic model (DDPM) [9] into the
training pipeline. Additionally, following consistency regularization methods,
we impose a regularization objective that encourages our image reconstruction
network to produces diverse outputs.

###  2.2 Regularization-by-denoising

Image denoising algorithms (denoisers) can serve as powerful priors on the
form and distribution of natural images [7, 25, 14]. Regularization-by-
denoising (RED) [21, 19] is a regularization technique that leverages an image
denoiser to solve classical single-image imaging inverse problems. Given an
known forward model $\mathcal{F}$ and a single measurement $y$, RED (under a
white Gaussian prior on the measurement noise) tries to reconstruct the
unknown scene $x$ by minimizing

$\arg\min_{x}\|y-\mathcal{F}(x)\|^{2}+\rho(x),$ (3)

where

$\displaystyle\rho(x)=\frac{1}{2}x^{T}\left[x-D(x)\right],$ (4)

and $D(\cdot)$ represents a denoiser like BM3D [5] or DnCNN [30]. Variations
on this idea have recently been combined with diffusion models as well [4].

One simplified interpretation of the RED regularization objective (4) is that
a good denoiser $D(\cdot)$ will project $x$ onto the manifold $\mathcal{M}$ of
natural images. If $x$ is far from that manifold $\rho(x)$ will be large,
whereas if $x$ is already close to this manifold, $D(x)$ will change it very
little and $\rho(x)$ will be small.

Similar to RED, our method leverages an image denoiser. However, instead of
using the denoiser to regularize classical inverse imaging optimization
objectives, we follow [28] and use an image denoiser to regularize the
training of deep neural networks. That is, we are using denoiser to regularize
a function rather than an image. The resulting function/reconstruction
network, unlike the RED, is non-iterative and does not require explicit
knowledge of the forward model $\mathcal{F}$.

Figure 1: Overview of the training pipeline. Block diagram of the semi-
supervised training pipeline used in our experiments. The pre-trained
diffusion model supervises training by pushing outputs of the image
reconstruction network towards the desired target image distribution. Note
that only the image reconstruction network (a standard U-net in our
experiments) is used at inference time.

##  3 Methods

###  3.1 Supervision-by-denoising

Young et al. [28] recently introduced the supervision-by-denoising (SUD)
framework which takes the ideas behind regularization-by-denoising [21] and
extends them to enable semi-supervised learning. The intuition behind SUD is
that a pretrained image denoiser, $D_{\sigma}(u)$, (which can be trained using
the set $U_{x}$ of unpaired images) encodes strong priors on the distribution
$p_{x}$. SUD enforces that the network’s reconstructions $f_{\theta}(y_{u})$
on the unpaired training data are “consistent” with the priors encoded in the
denoiser.

When used in combination with an $\ell_{2}$ loss and without temporal-
ensembling/damping, SUD effectively minimizes

$\displaystyle\mathcal{L}_{\text{paired}}+\lambda_{1}\mathcal{L}_{\text{denoiser}},$
(5)

where $\lambda_{1}$ is a tuning parameter and

$\displaystyle\mathcal{L}_{\text{denoiser}}=\frac{1}{|U_{y}|}\sum_{y_{u}\in
U_{y}}\|f_{\theta}(y_{u})-z_{u}\|^{2},$ (6)

where $z_{u}=D_{\sigma}(f_{\theta}(y_{u}))$. When updating the network weights
$\theta$ to minimize (6), SUD treats $z_{u}$ as a fixed pseudo-label and does
not propagate gradients through the denoiser. That is, SUD defines the
gradient of $\mathcal{L}_{\text{denoiser}}$ with respect to a single
reconstruction $f_{\theta}(y_{u})$ as

$\displaystyle\nabla_{f_{\theta}(y_{u})}\mathcal{L}_{\text{denoiser}}=\frac{2[f_{\theta}(y_{u})-D_{\sigma}(f_{\theta}(y_{u}))]}{|U_{y}|}.$
(7)

As demonstrated in [28], SUD is a powerful and effective semi-supervised
learning technique in the context of medical segmentation, where the goal is
to map an image to a discrete-valued segmentation map. Using only a handful
paired images and segmentation maps, Young et al. were able to train a
denoiser to segment brains, kidneys, and tumors.

Unfortunately, we found that without modification, SUD, with or without
temporal ensembling, was far less effective at general image restoration
tasks. As illustrated in Figure 2, minimizing the SUD loss for CelebA face
inpainting leads to mode collapse.

###  3.2 Understanding SUD

In this section we analyze SUD in order to identify and overcome its
weaknesses.

######  Theorem 3.1.

When $D_{\sigma}$ is a minimum-mean-squared error (MMSE) Gaussian denoiser,
minimizing $\mathcal{L}_{\text{denoiser}}$ minimizes the cross entropy between
the distribution of $f_{\theta}(y_{u})$ and the smoothed version of $p_{x}$.

###### Proof.

Let $\nu$ follow an independent zero-mean white Gaussian distribution with
variance $\sigma^{2}$. We will use $\nu$ to smooth the distributions $p_{x}$
(recall $p_{x+\nu}=p_{x}*p_{\nu}$, where $*$ denotes convolution) so that we
can take advantage of Tweedie’s Formula, as described below.

The cross entropy between $p_{f_{\theta}(y)}$ and $p_{x+\nu}$ is, by
definition,

$\displaystyle H(p_{f_{\theta}(y)},p_{x+\nu})$
$\displaystyle=-\mathbb{E}_{f_{\theta}(y)}[\ln p_{x+\nu}(f_{\theta}(y))].$ (8)

We can form a Monte-Carlo approximate of the expectation in (8) by averaging
over $U_{y}$:

$\displaystyle
H(p_{f_{\theta}(y)},p_{x+\nu})\approx-\frac{1}{|U_{y}|}\sum_{y_{u}\in
U_{y}}\ln p_{x+\nu}(f_{\theta}(y_{u})).$ (9)

Then, we can express the gradient of this loss with respect to a
reconstruction $f_{\theta}(y_{u})$ as

$\displaystyle\nabla_{f_{\theta}(y_{u})}H(p_{f_{\theta}(y)},p_{x+\nu})$
$\displaystyle\approx-\frac{\nabla_{f_{\theta}(y_{u})}\ln
p_{x+\nu}(f_{\theta}(y_{u}))}{|U_{y}|}.$ (10)

To efficiently evaluate (3.2) we turn to Tweedie’s Formula. Tweedie’s Formula
[6] states that for a signal corrupted with zero-mean additive white Gaussian
noise, $r=x+\nu$ where $\nu\sim\mathcal{N}(0,\sigma^{2}\mathbf{I})$, the
output of a MMSE denoiser $D_{\sigma}(\cdot)$ (and by extension a neural
network trained to act as a MMSE denoiser) can be expressed as

$D_{\sigma}(r)=r+\sigma^{2}\nabla_{r}\ln p_{x+\nu}(r).$ (11)

In other words, denoisers perform gradient ascent on the log-likelihood of
$p_{x+\nu}$ where the step size corresponds to the noise variance.
Accordingly, we can express the gradient of the log-likelihood in terms of the
denoiser’s residual:

$\displaystyle\nabla_{r}\ln p_{x+\nu}(r)=\frac{D_{\sigma}(r)-r}{\sigma^{2}}.$
(12)

By applying Tweedie’s formula to (3.2) we arrive at

$\displaystyle\nabla_{f_{\theta}(y_{u})}H(p_{f_{\theta}(y)},p_{x+\nu})\approx\frac{[f_{\theta}(y_{u})-D_{\sigma}(f_{\theta}(y_{u}))]}{\sigma^{2}|U_{y}|}.$
(13)

Up to constants, this is the same expression as the SUD gradients defined in
(7). As such, the SUD denoiser loss minimizes the cross entropy between
$p_{f_{\theta}(y)}$ and $p_{x+\nu}$. ∎

w/ $\mathcal{L}_{\text{reg}}$ | \begin{overpic}[width=39.0242pt]{img/mode_collapse/epoch10_1.png}\end{overpic} | \begin{overpic}[width=39.0242pt]{img/mode_collapse/epoch20_1.png}\end{overpic} | \begin{overpic}[width=39.0242pt]{img/mode_collapse/epoch30_1.png}\end{overpic} | \begin{overpic}[width=39.0242pt]{img/mode_collapse/epoch40_1.png}\end{overpic} |  \begin{overpic}[width=39.0242pt]{img/mode_collapse/epoch50_1.png}\end{overpic}
---|---|---|---|---|---
| Epoch 10 | Epoch 20 | Epoch 30 | Epoch 40 | Epoch 50

Training duration

Figure 2: SUD with and without correlation minimization. Without a correlation
loss to regularize reconstructions, the SUD denoising objective causes network
outputs to converge to a mode—in this case a washed out image.

######  Corollary 3.2.

Minimizing $\mathcal{L}_{\text{denoiser}}$ encourages mode collapse.

###### Proof.

Minimizing $\mathcal{L}_{\text{denoiser}}$ minimizes the cross entropy between
$p_{f_{\theta}(y_{u})}$ and $p_{x+\nu}$. The cross entropy of $H(p,q)$ of two
distributions $p$ and $q$ is minimized with respect to $p$ when $p$ is a dirac
distribution with a non-zero support where distribution $q$ is largest, i.e.,
a mode. ∎

Examples of SUD leading to mode collapse are presented in Figure 2.

######  Corollary 3.3.

Minimizing $\mathcal{L}_{\text{denoiser}}$ can encourage blurry
reconstructions.

###### Proof.

Minimizing $\mathcal{L}_{\text{denoiser}}$ minimizes the cross entropy between
$p_{f_{\theta}(y_{u})}$ and $p_{x+\nu}$ and will result in solutions
$f_{\theta}(y_{u})$ that maximize $p_{x+\nu}(f_{\theta}(y_{u}))$. For
sufficiently large $\sigma$, $p_{x+\nu}$ is maximized not where $p_{x}$ is
large (along the manifold of natural images) but rather at some point in
between high-probability points. ∎

A toy example illustrating how introducing noise onto a random variable can
move its distribution’s maxima is presented in Figure 3.

\begin{overpic}[width=71.54768pt]{img/Toy/Dist1.png}\end{overpic} | \begin{overpic}[width=71.54768pt]{img/Toy/Dist2.png}\end{overpic} | \begin{overpic}[width=71.54768pt]{img/Toy/Dist3.png}\end{overpic}
---|---|---
$p_{x}$ | $p_{x+\nu_{1}}$ | $p_{x+\nu_{2}}$
Figure 3: PDFs of a bimodal distribution with varying amounts of noise. Adding
sufficient noise can move the distribution’s mode.

###  3.3 Improving SUD

To fight mode collapse, we introduce an additional penalty,
$\mathcal{L}_{\text{reg}}$ into the SUD loss which encourages diverse outputs.
That is, we minimize

$\displaystyle\mathcal{L}_{\text{paired}}+\lambda_{1}\mathcal{L}_{\text{denoiser}}+\lambda_{2}\mathcal{L}_{\text{reg}}$
(14)

where $\lambda_{1}$ and $\lambda_{2}$ are scalar weights on each loss term.

######  Theorem 3.4.

When no paired training data is present, SUD with
$\lambda_{1}=\frac{1}{2\sigma^{2}}$, $\lambda_{2}=1$, and
$\mathcal{L}_{\text{reg}}=-H(p_{f_{\theta}(y)})$ minimizes the KL divergence
between the distributions of $f_{\theta}(y)$ and $x+\nu$.

###### Proof.

Our results from Theorem 3.1 indicate
$\frac{1}{2\sigma^{2}}\mathcal{L}_{\text{denoiser}}$ and
$H(p_{f_{\theta}(y)},p_{x+\nu})$ have the same gradients. Thus minimizing

$\displaystyle\frac{1}{2\sigma^{2}}\mathcal{L}_{\text{denoiser}}+\mathcal{L}_{\text{reg}}$
(15)

minimizes

$\displaystyle H(p_{f_{\theta}(y)},p_{x+\nu})-H(p_{f_{\theta}(y)}).$ (16)

The latter expression is the definition of KL-divergence between the
distributions of $f_{\theta}(y)$ and $x+\nu$. ∎

####  3.3.1 Correlation minimization

Computing and maximizing the entropy $H(p_{f_{\theta}(y_{u})})$ is
computationally intractable. To get around this hurdle, we instead encourage
sample diversity as follows: (1) we use a U-net architecture as our
reconstruction network $f_{\theta}(\cdot)$, (2) we let $a_{u}$ represent the
intermediate activation of $f_{\theta}$ to inputs $y_{u}$, and (3) we then
penalize correlations between activations.

Specifically, in each mini-batch, we compute the normalized covariance matrix
on the intermediate outputs from the encoder block of our network. Diagonal
entries of the matrix—which contain the correlation coefficient of each vector
with itself—are all equal to 1 by definition. Non-diagonal entries of the
matrix contain the Pearson correlation coefficient (PCC) between the latent
vectors in a mini-batch.

PCC identifies linear correlations between the latent vectors, where a value
of 0 signifies linearly uncorrelated samples, a value 1 signifies a positive
linear correlation between samples, and a value of $-1$ signifies a negative
linear correlation between samples. By minimizing this value for each entry in
the correlation matrix, we encourage network to produce outputs which are
uncorrelated in latent space.

####  3.3.2 Noise injection

with noise | \begin{overpic}[width=62.87344pt]{img/noise_ablation/noise_1.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/noise_ablation/noise_3.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/noise_ablation/noise_4.png}\end{overpic}
---|---|---|---
without noise | \begin{overpic}[width=62.87344pt]{img/noise_ablation/no_noise_1.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/noise_ablation/no_noise_3.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/noise_ablation/no_noise_4.png}\end{overpic}
Figure 4: SUD with and without noise injection. Injecting noise onto the
reconstructed images before denoising them results in more realistic-looking
reconstructions.

SUD compares the distribution of the reconstructions $f_{\theta}(y)$ with the
distribution $p_{x+\nu}$. As noted in the previous section, this encourages
solutions for which $p_{x+\nu}$ is large but $p_{x}$ is small: i.e., it can
produce solutions off of the manifold of natural images.

To alleviate this problem, we inject noise onto the reconstructions
$f_{\theta}(y_{u})$ before passing them through the denoiser. That is, we
redefine $\nabla_{f_{\theta}(y_{u})}\mathcal{L}_{\text{denoiser}}$ as

$\displaystyle\nabla_{f_{\theta}(y_{u})}\mathcal{L}_{\text{denoiser}}=\frac{2[f_{\theta}(y_{u})-D_{\sigma}(f_{\theta}(y_{u})+\nu_{2})]}{|U_{y}|},$
(17)

where $\nu_{2}\sim N(0,\sigma_{2}^{2}\mathbf{I})$.

This simple modification allows us to compare the smoothed distribution
$p_{f_{\theta}(y)+\nu_{2}}$ with the smoothed distribution $p_{x+\nu}$. As
demonstrated in Figure 4, noise injection produces reconstructions of
considerably higher quality.

####  3.3.3 Diffusion models

DDPM | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/diffusion_1.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/diffusion_2.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/diffusion_3.png}\end{overpic}
---|---|---|---
AE | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/denoiser_1.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/denoiser_2.png}\end{overpic} | \begin{overpic}[width=62.87344pt]{img/diffusion_ablation/denoiser_3.png}\end{overpic}
Figure 5: SUD with and without denoising diffusion models. Image in-painting
samples generated by SUD with a multi-step denoising diffusion model (top) and
a one-step autoencoder image denoiser (bottom). SUD training with the
denoising diffusion model produces higher-quality reconstructions.

As mentioned in Section 2.2, an alternative and more heuristic interpretation
of denoising algorithms is that they are projecting the reconstructions onto a
manifold $\mathcal{M}$ of allowable reconstructions, e.g., faces or “natural
images”. Traditional denoising algorithms perform this projection in a single
step. However, existing theory [26] suggests that one should navigate image
manifolds gradually, in a smooth-to-rough/coarse-to-fine manner.

Loosely inspired by this observation, we propose replacing our single-step
MMSE denoising algorithm with a multi-step denoising diffusion probabilistic
model (DDPMs) [9]. That is, we replace our denoiser $D_{\sigma}(\cdot)$ used
in (17) with an iterative forward “noising” operator $F(\cdot)$ and an
iterative reverse “denoising” operator $R(\cdot)$ such that
$D_{\sigma}(f_{\theta}(y_{u})+\nu)=R(F(f_{\theta}(y_{u})))$. The definitions
for each operator are expressed as

(18)

where $\alpha$ controls the noise variance at each time step $t$,
$z\sim\mathcal{N}(0,I)$, $\epsilon$ represents a diffusion network with
weights $\tau$, and $\bar{\alpha}_{t}=\prod^{t}_{i}\alpha_{i}$.

Conceptually, the DDPM serves to first project $r$ onto the smooth manifold of
noisy images and then gradually project $r$ onto correspondingly less smooth
manifolds of less noisy images. As demonstrated in Figure 5, training with a
DDPM produces higher quality results than training with a denoising auto-
encoder.

###  3.4 SUD2

We refer to SUD, including all the aforementioned modifications, as SUD2.
Ablation studies characterizing the relative importance of correlation
minimization, noise injection, and denoising diffusion models are provided in
the supplement.

##  4 Experiments

In the following experiments, we evaluate our method both quantitatively and
qualitatively against baselines. All methods below are trained on 5 paired
images, and the semi-supervised methods are trained on additional unpaired
images.

The backbone architecture for our image reconstruction network is a U-net [22]
consisting of 4 down-sampling and up-sampling blocks—implemented using strided
convolutions—with skip connections between them. Each down/up-sampling block
contains two convolutional layers, each with batch normalization and a Leaky
ReLU activation function [15]. The denoising diffusion models used in our
experiments have a similar U-net architecture, albeit with 6 down/up-sampling
blocks and spatial self-attention. Similarly, our denoising networks use an
autoencoder backbone identical to the image reconstruction U-net, but with
skip connections removed.

We train all of the image reconstruction networks on $4\times$ Nvidia RTX
A5000 GPUs using an Adam optimizer with an initial learning rate of $1\times
10^{-3}$, a weight decay of $1\times 10^{-4}$, and a batch size of 8 for 50
epochs. Additionally, we resize images to $256\times 256$ resolution and
normalize pixel intensities between $[-1,1]$ before training.

When evaluating our method on the experiments below, we set
$\lambda_{1}=0.01,\lambda_{2}=10$ and supervise training by walking 400 steps
along the forward diffusion process before walking 400 steps along the reverse
diffusion process. Generally, we find $[300,700]$ steps to be a reasonable
range to take. Traversing more than 700 steps along the diffusion chain adds
excess amounts of noise to the image, resulting in an image which no longer
resembles the input image. In contrast, traversal of less than 300 steps
yields images which are approximately identical to the input images—making
them poor targets for supervision.

Instead of explicitly setting the gradient of the loss, we equivalently
implement (17) by disabling backpropagation of gradients through the diffusion
model $D_{\sigma}(\cdot)$ before computing the $\ell_{2}$ loss between
$f_{\theta}(y_{u})$ and $D_{\sigma}(f_{\theta}(y_{u}))$. We also add a
perceptual loss—specifically Learned Perceptual Image Patch Similarity (LPIPS)
[31]—to $\mathcal{L}_{\text{paired}}$ and $\mathcal{L}_{\text{denoiser}}$ as
we empirically find it to help produce better quality reconstructions.

###  4.1 Image in-painting

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sup_1.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_cyclegan_1.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_labels_1.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud_1.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud2_1.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sup_2.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_cyclegan_2.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_labels_2.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud_2.png}\end{overpic}

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud2_2.png}\end{overpic}

(a) Input

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sup_3.png}\end{overpic}
(b) Supervised

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_cyclegan_3.png}\end{overpic}
(c) CycleGAN

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_labels_3.png}\end{overpic}
(d) CycleGAN-SSL

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud_3.png}\end{overpic}
(e) SUD

\begin{overpic}[width=62.87344pt]{img/teaser/celeba_sud2_3.png}\end{overpic}
(f) $\text{SUD}^{2}$

(g) Reference

Figure 6: Qualitative in-painting results. The images above show inference
results from several methods on CelebAMask-HQ in-painting. All methods are
trained on 5 pairs of masks and faces. Semi-supervised methods (CycleGAN-SSL,
SUD, SUD2) are trained on an additional 1000+12,500 unpaired faces.

Image in-painting is a generative process for reconstructing missing regions
of an image such that restored image fits a desired—often natural—image
distribution. We test our method on the CelebAMask-HQ dataset [12], which
contains 30,000 images of $512\times 512$ resolution and their corresponding
segmentation maps. This dataset allows us to learn a strong prior since the
image distribution is relatively constrained—all images are one-quarter head
shots of celebrities.

In this experiment, we mask out the subject’s face from each image and train a
few shot, semi-supervised in-painting network on 5 paired images and
1000+12,500 unpaired images. The 12,500 unpaired images are used to pre-train
a blind Gaussian denoiser, a diffusion model, and a CycleGAN [33]. Note that
the denoiser and diffusion model are only trained on 12,500 faces without a
mask applied. However, we train the CycleGAN on non-corresponding pairs of
12,500 faces without a mask applied and 12,500 faces with a mask applied to
give it the best possible performance. As an additional baseline, we also
train an image reconstruction network using the pre-trained CycleGAN as a
pseudo-label generator, which we refer to as CycleGAN-SSL.

Despite requiring far fewer images than both CycleGAN baselines, SUD2 produces
faces which most closely resemble the ground truth distribution. Furthermore,
the strong prior imposed by the pre-trained diffusion model results in
reconstructions with more defined facial structure compared to baseline
methods. Notably, as described in Corollary (3.2), the SUD baseline collapses
to a mode during training with high probability, yielding highly correlated
reconstructions.

Quantitatively, SUD2 achieves the most consistent results across our test set
of $768$ images, with a $48\%$ higher average PSNR over CycleGAN-SSL and a
$43\%$ lower average FID score compared to the supervised baseline. Although
CycleGAN achieves the best FID score overall, it tends to produce
qualitatively poor faces with disproportionately sized features, which is
reflected in its poor PSNR, SSIM, and LPIPS scores. Likewise, while the
supervised baseline achieves PSNR, SSIM, and LPIPS scores comparable to SUD2,
it often generates faces with missing features (i.e. eyes, nose, mouth), which
is indicated by its high FID score.

Image in-painting
---
Method | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$ | FID$\downarrow$
Supervised | 18.44 | 0.71 | 0.29 | 0.48
CycleGAN | 8.77 | 0.23 | 0.66 | 0.17
CycleGAN-SSL | 11.38 | 0.60 | 0.39 | 1.14
SUD | 11.28 | 0.29 | 0.69 | 3.07
$\text{SUD}^{2}$ (Ours) | 18.71 | 0.71 | 0.28 | 0.31
Figure 7: Quantitative image in-painting results. This table lists the average
scores attained by each method on the CelebAMask-HQ image in-painting test
set. The best scores are bolded for readability. Image dehazing
---
Method | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$
Supervised | 17.65 | 0.66 | 0.44
CycleGAN | 12.80 | 0.22 | 0.64
$\text{SUD}^{2}$ (Ours) | 17.74 | 0.64 | 0.37
Figure 8: Quantitative image dehazing results. This table lists the median
scores attained by each method on the REVIDE image dehazing test set. The best
scores are bolded for readability.

\begin{overpic}[width=86.72267pt]{img/haze/supervise_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:16.53}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.58}}}
\end{overpic}

\begin{overpic}[width=86.72267pt]{img/haze/cyclegan_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:15.68}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.42}}}
\end{overpic}

\begin{overpic}[width=86.72267pt]{img/haze/sud2_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:17.76}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.66}}}
\end{overpic}

\begin{overpic}[width=86.72267pt]{img/haze/haze_sup_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:21.70}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.75}}}
\end{overpic}

\begin{overpic}[width=86.72267pt]{img/haze/haze_cyclegan_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:19.66}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.70}}}
\end{overpic}

\begin{overpic}[width=86.72267pt]{img/haze/haze_sud2_1.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:22.01}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.76}}}
\end{overpic}

(a) Input

\begin{overpic}[width=86.72267pt]{img/haze/haze_sup_2.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:16.00}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.54}}}
\end{overpic} (b) Supervised

\begin{overpic}[width=86.72267pt]{img/haze/haze_cyclegan_2.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:
13.09}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.35}}}
\end{overpic} (c) CycleGAN

\begin{overpic}[width=86.72267pt]{img/haze/haze_sud2_2.png}
\put(5.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{PSNR:14.14}}}
\put(55.0,5.0){\scriptsize{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{SSIM:0.62}}}
\end{overpic} (d) $\text{SUD}^{2}$

(e) Reference

Figure 9: Dehazing results. Inference results comparing our method against
baselines on dehazing REVIDE bedrooms. All methods use the same network
architecture (U-net). All methods are provided 5 hazy / clear image pairs, 200
unpaired hazy images.

###  4.2 Dehazing

####  4.2.1 Haze forward model

Image dehazing is a challenging imaging inverse problem, where the objective
is to remove degradations caused by small particulates in the air which
obstruct visibility. Typically, dehazing methods approximate the hazy image
formation as an depth-dependent attenuation process, defined by the well-known
dichromatic atmospheric scattering model [16, 17]

$I(x,y)=J(x,y)t(x,y)+A(1-t(x,y)).$ (19)

Each pixel in the hazy image $I$ is obtained by attenuating pixels in the
clear image $J$ according to a transmission map $t$ dependent on depth $d$.
Two additional scalar parameters $\beta$ and $A$ describe the magnitude of
attenuation and intensity of atmospheric illumination respectively.

Many prior works such as DehazeNet [3, 20] incorporate the atmospheric
scattering model in a deep learning framework to train a dehazing network on
synthesized hazy images. However, as shown by [2], these methods often fail
when applied to images of real haze—indicating a poor fit to the true forward
model of haze.

####  4.2.2 Removing real haze

Currently, few methods exist for learning to dehaze exclusively on images of
real haze. CycleDehaze [8] trained on real hazy images from the O-HAZE [2] and
I-HAZE [1] benchmarks; however, the low number of hazy samples—less than
100—in those datasets, compounded with the data-hungry nature of neural
networks, inhibits the dehazing network from producing high-quality
reconstructions.

In this experiment, we aim to overcome these limitations and evaluate our
method on a subset of bedroom scenes from the REVIDE dataset [32], which
contains a total of 240 real hazy image pairs across 9 scenes. This dataset
represents an ideal use case for our method, as only a restricted number of
paired images are available to train on.

Similar to the in-painting experiment, we train a U-net on 5 paired images and
200+20,000 unpaired images, where the 5 paired and 200 unpaired images are
sourced from REVIDE and the additional 20,000 unpaired images are sourced from
the LSUN bedrooms dataset [29]. Before training the image reconstruction
network, we first pre-train both a diffusion model and a CycleGAN.
Specifically, we train the diffusion model on 20,000 non-hazy bedroom images
from the LSUN dataset [29] and train CycleGAN on the same LSUN images (along
with the 5 paired and 200 unpaired images from REVIDE).

Compared to baselines methods, we find that SUD2 produces reconstructions with
far fewer artifacts and achieves a median LPIPS score $17\%$ lower than the
supervised baseline. Visually, the dehazed images generated by SUD2 have a
more natural and smooth appearance relative to both CycleGAN and the
supervised baseline as shown in Figure 9. In contrast, due to the limited
number of hazy training samples, CycleGAN fails to learn a good mapping
between hazy and non-hazy images, resulting in poor PSNR, SSIM, and LPIPS
scores.

##  5 Limitations and future work

Our experiments demonstrate that SUD2 enables the use of deep learning
techniques even when paired data is scarce. Nonetheless, there remain
limitations to our approach which can be addressed in future works.

For instance, performance improvements gained from SUD2 are directly
correlated with the strength of the prior imposed by the unconditional
denoising diffusion model. As the target image distribution grows more
diverse, the strength of this learned prior weakens. Conditional diffusion
models offer a promising avenue for maintaining a strong prior on diverse
image distributions that can be explored in future works.

Furthermore, the reconstruction networks produced by the SUD2 training
procedure are deterministic. Extending SUD2 to train stochastic algorithms
which produce diverse outputs [10] is another interesting direction for future
work.

##  6 Conclusion

We introduce SUD2, a generalized deep learning framework for solving few-shot,
semi-supervised image reconstruction problems. Inspired by the recent success
of denoising diffusion models on image generation tasks, we leverage diffusion
models to regularize network training, encouraging solutions that lie close to
the desired image distribution.

To benchmark our method, we apply SUD2 on image reconstruction tasks and
compare against CycleGAN, a well-studied generative model that utilizes
unpaired data to perform image-to-image translation. When applied to image in-
painting, we find that our method produces significantly more structured faces
than baselines, where each facial feature is proportionately sized and
reasonably positioned on the face. Similarly, when applied to image dehazing,
we observe that results from our method are far less noisy compared to
baseline methods. Based off of both qualitative and quantitative experimental
results, SUD2 succeeds in enabling the training of deep networks on datasets
where few paired data samples are available.

##  7 Acknowledgements

M.C. and C.M. were supported in part by the AFOSR Young Investigator Program
Award FA9550-22-1-0208 and a Northrop Grumman seed grant.

## References

  * [1] Codruta O. Ancuti, Cosmin Ancuti, Radu Timofte, and Christophe De Vleeschouwer.  I-haze: a dehazing benchmark with real hazy and haze-free indoor images.  In arXiv:1804.05091v1, 2018.
  * [2] Codruta O. Ancuti, Cosmin Ancuti, Radu Timofte, and Christophe De Vleeschouwer.  O-haze: a dehazing benchmark with real hazy and haze-free outdoor images.  In IEEE Conference on Computer Vision and Pattern Recognition, NTIRE Workshop, NTIRE CVPR’18, 2018.
  * [3] Bolun Cai, Xiangmin Xu, Kui Jia, Chunmei Qing, and Dacheng Tao.  Dehazenet: An end-to-end system for single image haze removal.  Trans. Img. Proc., 25(11):5187–5198, nov 2016.
  * [4] Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  In The Eleventh International Conference on Learning Representations, 2023.
  * [5] Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian.  Image denoising by sparse 3-d transform-domain collaborative filtering.  IEEE Transactions on image processing, 16(8):2080–2095, 2007.
  * [6] Bradley Efron.  Tweedie’s formula and selection bias.  Journal of the American Statistical Association, 106:1602–1614, 12 2011.
  * [7] Karen Egiazarian, Alessandro Foi, and Vladimir Katkovnik.  Compressed sensing image reconstruction via recursive spatially adaptive filtering.  In 2007 IEEE International Conference on Image Processing, volume 1, pages I–549. IEEE, 2007.
  * [8] D. Engin, A. Genc, and H. Ekenel.  Cycle-dehaze: Enhanced cyclegan for single image dehazing.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 938–9388, Los Alamitos, CA, USA, jun 2018\. IEEE Computer Society.
  * [9] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems, 33:6840–6851, 2020.
  * [10] Zahra Kadkhodaie and Eero Simoncelli.  Stochastic solutions for linear inverse problems using the prior implicit in a denoiser.  Advances in Neural Information Processing Systems, 34:13242–13254, 2021.
  * [11] Samuli Laine and Timo Aila.  Temporal ensembling for semi-supervised learning.  arXiv preprint arXiv:1610.02242, 2016.
  * [12] Cheng-Han Lee, Ziwei Liu, Lingyun Wu, and Ping Luo.  Maskgan: Towards diverse and interactive facial image manipulation.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [13] Dong-Hyun Lee.  Pseudo-label : The simple and efficient semi-supervised learning method for deep neural networks.  ICML 2013 Workshop : Challenges in Representation Learning (WREPL), 07 2013.
  * [14] Christopher A Metzler, Arian Maleki, and Richard G Baraniuk.  From denoising to compressed sensing.  IEEE Transactions on Information Theory, 62(9):5117–5144, 2016\.
  * [15] Vinod Nair and Geoffrey E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, page 807–814, Madison, WI, USA, 2010. Omnipress.
  * [16] Srinivasa G. Narasimhan and Shree K. Nayar.  Vision and the atmosphere.  International Journal of Computer Vision, 48:233–254, 2002.
  * [17] Shree K Nayar and Srinivasa G Narasimhan.  Vision in bad weather.  In Proceedings of the seventh IEEE international conference on computer vision, volume 2, pages 820–827. IEEE, 1999.
  * [18] Gregory Ongie, Ajil Jalal, Christopher A Metzler, Richard G Baraniuk, Alexandros G Dimakis, and Rebecca Willett.  Deep learning techniques for inverse problems in imaging.  IEEE Journal on Selected Areas in Information Theory, 1(1):39–56, 2020.
  * [19] Edward T Reehorst and Philip Schniter.  Regularization by denoising: Clarifications and new interpretations.  IEEE transactions on computational imaging, 5(1):52–67, 2018.
  * [20] Wenqi Ren, Si Liu, Hua Zhang, Jinshan Pan, Xiaochun Cao, and Ming-Hsuan Yang.  Single image dehazing via multi-scale convolutional neural networks.  In Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part II 14, pages 154–169. Springer, 2016.
  * [21] Yaniv Romano, Michael Elad, and Peyman Milanfar.  The little engine that could: Regularization by denoising (red).  SIAM Journal on Imaging Sciences, 10(4):1804–1844, 2017.
  * [22] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, pages 234–241, Cham, 2015. Springer International Publishing.
  * [23] Alexander Sorokin and David Forsyth.  Utility data annotation with amazon mechanical turk.  In 2008 IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops, pages 1–8, 2008.
  * [24] Antti Tarvainen and Harri Valpola.  Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 1195–1204, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [25] Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg.  Plug-and-play priors for model based reconstruction.  In 2013 IEEE Global Conference on Signal and Information Processing, pages 945–948. IEEE, 2013.
  * [26] Michael B Wakin, David L Donoho, Hyeokho Choi, and Richard G Baraniuk.  High-resolution navigation on non-differentiable image manifolds.  In Proceedings.(ICASSP’05). IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005., volume 5, pages v–1073. IEEE, 2005.
  * [27] Xiangli Yang, Zixing Song, Irwin King, and Zenglin Xu.  A survey on deep semi-supervised learning.  IEEE Transactions on Knowledge and Data Engineering, 2022.
  * [28] Sean I Young, Adrian V Dalca, Enzo Ferrante, Polina Golland, Bruce Fischl, and Juan Eugenio Iglesias.  Sud: Supervision by denoising for medical image segmentation.  arXiv preprint arXiv:2202.02952, 2022.
  * [29] Fisher Yu, Ari Seff, Yinda Zhang, Shuran Song, Thomas Funkhouser, and Jianxiong Xiao.  Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop.  arXiv preprint arXiv:1506.03365, 2015.
  * [30] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang.  Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising.  IEEE transactions on image processing, 26(7):3142–3155, 2017.
  * [31] R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 586–595, Los Alamitos, CA, USA, jun 2018. IEEE Computer Society.
  * [32] Xinyi Zhang, Hang Dong, Jinshan Pan, Chao Zhu, Ying Tai, Chengjie Wang, Jilin Li, Feiyue Huang, and Fei Wang.  Learning to restore hazy video: A new real-world dataset and a new method.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 9239–9248, June 2021.
  * [33] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In Computer Vision (ICCV), 2017 IEEE International Conference on, 2017.

See pages 1 of supplement.pdf See pages 0 of supplement.pdf


76D05

Non-intrusive reduced-order modeling method for two-way coupled FSI. Sparse,
nonlinear data-driven model inference. Coupled FSI model for Hron-Turek FSI3
benchmark testcase.

# Adjacency-based, Non-intrusive Reduced-order Modeling for Fluid-Structure
Interactions

Leonidas Gkimisis  Max Planck Institute for Dynamics of Complex Techincal
Systems, Sandtorstraße 1, 39106 Magdeburg, Germany, Computational Methods in
Systems and Control Theory (CSC).  corresponding author,  Thomas Richter
Otto-von-Guericke Otto-von-Guericke- Universität Magdeburg, Universitätsplatz
2, 39106 Magdeburg, Germany, Faculty of Mathematics, Institute for Analysis
and Numerics.  Peter Benner<EMAIL_ADDRESS>
###### Abstract

Non-intrusive model reduction is a promising solution to system dynamics
prediction, especially in cases where data are collected from experimental
campaigns or proprietary software simulations. In this work, we present a
method for non-intrusive model reduction applied to Fluid-Structure
Interaction (FSI) problems. The approach is based on the a priori known
sparsity of the full-order system operators, which is dictated by grid
adjacency information. In order to enforce this type of sparsity, we solve a
“local”, regularized least-squares problem for each degree of freedom on a
grid, considering only the training data from adjacent degrees of freedom,
thus making computation and storage of the inferred full-order operators
feasible. After constructing the non-intrusive, sparse full-order model,
Proper Orthogonal Decomposition (POD) is used for its projection to a reduced
dimension subspace and thus the construction of a reduced-order model (ROM).
The methodology is applied to the challenging Hron-Turek benchmark FSI3, for
Re = 200. A physics-informed, non-intrusive ROM is constructed to predict the
two-way coupled dynamics of a solid with a deformable, slender tail, subject
to an incompressible, laminar flow. Results considering the accuracy and
predictive capabilities of the inferred reduced models are discussed.

###### keywords:

Non-intrusive Model Reduction, Vortex-Induced Vibrations, Fluid-Structure
Interactions

##  1 Introduction

Fluid-Structure Interactions encompass coupled, strongly nonlinear fluid-solid
dynamics systems, which arise in various engineering fields. For different
regimes, materials and geometry configurations, the specific coupling
conditions and corresponding FSI dynamical response can widely differ [1]. In
this work, we focus on incompressible fluid dynamics at low Reynolds numbers,
coupled with deformable solids, at relatively low mass ratios. Under such
conditions, the timescales and inertia of the two subsystems are comparable,
leading to a strong, two-way coupling mechanism.

The development of non-intrusive, reduced-order models for FSI problems is
motivated by needs in engineering design and control. Firstly, the high
computational cost of direct FSI numerical simulation renders multi-query
tasks such as design optimization and real-time control cumbersome. Such
engineering tasks could benefit from model reduction techniques, which have
proven to allow for considerable speed-ups of FSI simulation, without
significant loss of accuracy [2]. Secondly, in cases where experimental or FSI
simulation data are available, without access to a numerical simulation source
code, non-intrusive modeling methods are developed for system dynamics
prediction [3]. Exploiting knowledge on the structure of the underlying
dynamical system has been shown to enhance the predictive capabilities of non-
intrusive ROMs for specific FSI regimes such as aeroelasticity [4].

In this work, we implement an adjacency-based approach to FSI problems. The
method was recently studied for linear advection and diffusion problems in [5]
and [6] and has been applied to the inference of quadratic-bilinear,
incompressible fluid dynamics operators in Vortex-Induced Vibrations problems
by the authors [7]. Two coupled, data-driven models are inferred for the fluid
dynamics and solid dynamics subsystems, exploiting a physics-based structure
of the governing equations under the Arbitrary Lagrangian-Eulerian (ALE)
formulation. We present results for the FSI3 Hron-Turek benchmark [8], which
has been previously studied for intrusive model reduction in [2]. Through this
testcase, we examine the predictive capabilities of the obtained, non-
intrusive model.

##  2 Theoretical Background

###  2.1 Governing equations

For the fluid flow, we focus on the regime of 2D incompressible, laminar
flows. In cases of moving fluid domains $\Omega(t)$, the equations can be
expressed with respect to a reference configuration $\hat{\Omega}$ using the
ALE formulation. The map from the reference domain to the current
configuration is given as:
$\hat{\mathbf{x}}\mapsto\mathbf{x}=\hat{\mathbf{x}}+\hat{\mathbf{d}}_{f}(t)$,
where $\hat{\mathbf{d}}_{f}(t)$ is a deformation field. The incompressible ALE
Navier-Stokes equations are then given by

$\left\\{\begin{array}[]{ll}\begin{aligned}
\rho_{f}J\left(\partial_{t}\hat{\mathbf{u}}+\hat{\nabla}\hat{\mathbf{u}}\mathbf{F}^{-1}\left(\hat{\mathbf{u}}-\partial_{t}\hat{\mathbf{d}}_{f}\right)\right)-\operatorname{div}\left(J\hat{\boldsymbol{\sigma}}\left(\hat{\mathbf{u}},\hat{p}\right)\mathbf{F}^{-T}\right)&=J\rho_{f}\vec{g}\quad\\\
\operatorname{div}\left(J\mathbf{F}^{-1}\hat{\mathbf{u}}\right)&=0\quad\end{aligned}\end{array}\right.,$
(1)

where $\hat{\mathbf{u}}$ is the ALE velocity field, $\hat{p}$ is the pressure,
$\hat{\boldsymbol{\sigma}}$ is the ALE stress tensor given by

$\hat{\boldsymbol{\sigma}}\left(\hat{\mathbf{u}},\hat{p}\right)=-\hat{p}I+\rho_{f}\nu_{f}\left(\hat{\nabla}\hat{\mathbf{u}}\mathbf{F}^{-1}+\mathbf{F}^{-T}\hat{\nabla}\hat{\mathbf{u}}^{T}\right),$
(2)

$\mathbf{F}$ is the gradient of the ALE map and $J=\det{\mathbf{F}}$, with

$\mathbf{F}=I+\hat{\nabla}\hat{\mathbf{d}}_{f}.$ (3)

For solid dynamics modelling, we use the linear Navier-Lamé equations,
assuming no material or geometrical nonlinearities. The boundary of the solid
is split into $\Gamma_{1}^{s}$ where the displacement is zero and
$\Gamma_{2}^{s}$ where the fluid-structure coupling occurs. The solid
deformation $\mathbf{d}_{s}$ is thus modeled by

$\left\\{\begin{array}[]{ll}\partial_{tt}\mathbf{d}_{s}-\nabla\cdot\boldsymbol{\sigma}_{s}=0\\\
\boldsymbol{\sigma}_{s}=\mu(\nabla\mathbf{d}_{s}+\nabla\mathbf{d}_{s}^{T})+\lambda\;tr(\frac{1}{2}(\nabla\mathbf{d}_{s}+\nabla\mathbf{d}_{s}^{T}))I\end{array}\right..$
(4)

figureSchematic representation (not scaled) of the Hron-Turek FSI benchmark
testcase. A velocity profile is assigned at the inlet $\Gamma_{1}^{f}$ and FSI
coupling occurs on the surface $\Gamma_{2}^{s}$ of a slender, deformable tail.

The boundary conditions of the coupled fluid-solid problem are given below,
based on Section 2.1.

$\left\\{\begin{array}[]{ll}\begin{aligned}
\hat{\mathbf{u}}&=\mathbf{u}_{in}\;&on\;\Gamma_{1}^{f},\\\
J\left(\nu_{f}\hat{\nabla}\hat{\mathbf{u}}\mathbf{F}^{-1}-\hat{p}I\right)\mathbf{F}^{-T}\cdot\hat{\mathbf{n}}_{f}&=\mathbf{0}\;&on\;\Gamma_{2}^{f},\\\
\hat{\mathbf{u}}&=\mathbf{0}\;&on\;\Gamma_{3}^{f},\\\
\mathbf{d}_{s}&=\mathbf{0}\;&on\;\Gamma_{1}^{s},\\\
J\hat{\sigma}_{f}\mathbf{F}^{-T}\cdot\hat{\mathbf{n}}_{f}&=\boldsymbol{\sigma}_{s}\cdot\mathbf{n}_{s}\;&on\;\hat{\Omega}\cap\cal
S,\\\ \hat{\mathbf{u}}&=\partial_{t}\mathbf{d}_{s}\;&on\;\hat{\Omega}\cap\cal
S.\\\ \end{aligned}\end{array}\right.$ (5)

The above formulation considers the coupled problem on a reference fluid
domain $\hat{\Omega}$. We use a stiffened harmonic extension as an ALE map, to
transfer the solution to the current, deformed configuration $\Omega(t)$.

$-\hat{\nabla}^{2}\hat{\mathbf{d}}_{f}-\frac{\hat{\nabla}a(\hat{\mathbf{x}})}{a(\hat{\mathbf{x}})}\hat{\nabla}\hat{\mathbf{d}}_{f}=\mathbf{0},\;\text{in}\;\hat{\Omega},$
(6)

with boundary conditions on the interface $\hat{\Omega}\cap\cal S$, and on the
domain external boundaries:

$\left\\{\begin{array}[]{l}\begin{aligned}
&\hat{\mathbf{d}}_{f}=\hat{\mathbf{d}}_{s}-\hat{\mathbf{d}}_{s}(t=t_{1}),&\;\text{on}\;\hat{\Omega}\cap\cal
S,\\\
&\hat{\mathbf{d}}_{f}=\mathbf{0},&\;\text{on}\;\Gamma_{1}^{f}\cup\Gamma_{2}^{f}\cup\Gamma_{3}^{f}.\end{aligned}\end{array}\right.$
(7)

$a(\hat{\mathbf{x}})$ introduces a stiffening effect close to the FSI
interface $\hat{\Omega}\cap\cal S$, which allows for larger mesh displacements
without element degradation, compared to a purely harmonic extension (see
Section 5.3.5 of [1]). $t_{1}$ is selected as a reference time, for which
$\hat{\Omega}=\Omega(t_{1})$.

###  2.2 Approximate model structure

By discretizing eq. 1, eq. 4 and eq. 6, we aim to reveal an approximate
structure of the dynamical system, which will be exploited for non-intrusive
modeling. We assume that the gradient of the deformation field
$\hat{\nabla}\hat{\mathbf{d}}_{f}$ is sufficiently small, such that $F\approx
I$ (and thus $J\approx 1$) is a valid approximation of eq. 2. In that case,
the structure of eq. 1 is approximately quadratic-bilinear. After discretizing
in space and time, working similarly to [9], we cancel out the algebraic
contribution of the pressure and come up with the following, approximate model
structure for the fluid velocity field:

$\mathbf{u}^{k+1}\approx
A\mathbf{u}^{k}+H(\mathbf{u}^{k+1})^{2}+K\;\partial_{t}\mathbf{d}_{f}\otimes\mathbf{u}^{k+1}+B\;\partial_{t}\mathbf{d}_{s}+L\;\mathbf{u}_{in}+C,$
(8)

where $(\mathbf{u})^{2}$ denotes the vector of $2n_{f}\choose 2$ unique
elements of the Kronecker product $\mathbf{u}\otimes\mathbf{u}$. For the solid
velocity field, we get a second-order, forced linear system by discretizing
eq. 4 in space and time

$\partial_{t}\mathbf{d}_{s}^{k+1}=A_{s}\partial_{t}\mathbf{d}_{s}^{k}+K_{s}\;\mathbf{d}_{s}^{k}+\mathbf{f}.$
(9)

The forcing $\mathbf{f}$ originates from the dynamic coupling condition on the
interface (i.e. the projection of the fluid stress tensor on the normal to the
interface). This can be approximated by a quadratic model, as follows.

$\mathbf{f}\approx
B_{s}\mathbf{u}_{FS}^{k+1}+H_{s}{\mathbf{u}_{FS}^{k+1}}^{2},$ (10)

where $FS$ denotes nodes that are sufficiently close to the interface, given
the sparsity of the discretized stress tensor (see eq. 2). This means that
$\mathbf{u}_{FS}=R\mathbf{u}$, where $R$ is an ordered selection matrix.

Finally, the solid velocity field $\partial_{t}\mathbf{d}_{s}$ should be
integrated to update the solid displacement field $\mathbf{d}_{s}$. To this
end, a Crank-Nicholson scheme is used

$\mathbf{d}_{s}^{k+1}=\mathbf{d}_{s}^{k}+\frac{\Delta
t}{2}(\partial_{t}\mathbf{d}_{s}^{k+1}+\partial_{t}\mathbf{d}_{s}^{k}).$ (11)

We also need to discretize and solve eq. 6, in order to compute the fluid mesh
displacement in the current configuration, $\Omega(t)$, with respect to the
reference configuration $\hat{\Omega}$. This leads to

$\mathbf{d}_{f}^{k}=A_{L}S\mathbf{d}_{s}^{k},$ (12)

where $A_{L}$ is the inverse of the discretized ALE map and $S$ is an ordered
selection matrix, which picks only the solid DoFs on $\Gamma_{2}^{s}$.
Equations (8) to (12) comprise a discretized, coupled dynamical system, which
models ALE FSI problems under small deformations. In the following, we use
this structure to fit the corresponding operators to numerical FSI data and
thus construct a physics-based, non-intrusive model.

##  3 FSI Non-Intrusive Model Reduction

###  3.1 Physics-based adjacency

Having identified an approximate structure for the FSI dynamical system, we
aim to infer the unknown operators of the fluid dynamics subsystem (eq. 8),
the solid dynamics subsystem (eq. 9) and the solid forcing term (eq. 10), from
simulation data.

The discretized operators in eq. 9 are sparse, while the operators of eq. 8
and eq. 10 are only approximately sparse, due to canceling out the
contribution from the pressure. The corresponding sparsity pattern originates
from grid adjacency. Based on this observation, the employed method aims to
infer full-order, adjacency-based, sparse operators.

To proceed with the inference task, we should map the available data for the
flowfield over time $t=[0,T]$ to a reference configuration $\hat{\Omega}$,
since equation eq. 8 was derived there. Knowing the motion of the FSI
interface from the data, we define a reference time as
$t_{1}=\text{argmin}(\|\hat{\mathbf{d}}_{s}(t)\|_{2})$. Based on the reference
position of the FSI interface at $t_{1}$, we can construct a fluid and a solid
mesh (e.g. via Delaunay triangulation), with the constraint that the nodes of
the two meshes on the interface should match.

For the fluid, we compute the deformed mesh $\Omega(t)$ from eq. 12, given the
motion of the FSI interface (see eq. 7) and the constructed reference
$\hat{\Omega}=\Omega(t_{1})$. Then, the fluid velocity data can be
interpolated on the reference configuration $\hat{\Omega}$ at each of the
$n_{T}$ timesteps. For the solid, we interpolate the displacement field values
for all $n_{T}$ timesteps from the solver mesh configuration at $t_{1}$, to
the constructed solid mesh. The difference between handling the solid and the
fluid problem should be noted: Mesh construction is necessary only for the
fluid side due to the ALE formulation. However, grid adjacency information for
both the fluid and the solid mesh is essential for the following. Thus, the
solid mesh of the simulation data can be readily used as long as adjacency
information can be retrieved.

###  3.2 Least Squares formulation and regularization

The number of grid nodes is denoted by $n_{f}$ for the fluid and $n_{s}$ for
the solid mesh. Then, there are $2n_{f}$ and $2n_{s}$ kinematic degrees of
freedom (DoFs), respectively, for this 2D problem. For the $i$-th DoF
(streamwise or transverse component of $\mathbf{u}$ or
$\partial_{t}\mathbf{d}$), we denote by $q_{A}(i)$ the set of geometrically
adjacent DoFs in the fluid $A\equiv f$ or solid $A\equiv s$ mesh respectively.
For example, if $i$ is a fluid velocity component, $q_{f}(i)$ includes
exclusively DoFs which correspond to fluid mesh nodes, adjacent to the node of
$i$. Then, the local inference of the operators for DoF $i$ corresponds to
solving the following problem:

$\min_{\boldsymbol{\beta}_{i}}\left\|\boldsymbol{\beta}_{i}^{T}\mathcal{D}_{i}-\mathbf{y}_{i}^{m+1}\right\|_{2},$
(13)

where we make a distinction between the fluid and solid dynamics problems:

For the fluid dynamics subsystem,

$\boldsymbol{\beta}_{i}=[C_{i},A_{i,.},H_{i,.},K_{i,.}]^{T},\;\;\mathcal{D}_{i}^{T}=\begin{bmatrix}\mathbf{1}\\\
\mathbf{u}_{q_{f}(i)}^{m}\\\ {\mathbf{u}_{q_{f}(i)}^{m}}^{2}\\\
\mathbf{u}_{q_{f}(i)}^{m}\otimes\partial_{t}\hat{\mathbf{d}}_{f}^{m}\\\
\end{bmatrix},\;\;\mathbf{y}_{i}^{m+1}=[u_{i}^{t=2}\;...\;u_{i}^{t=n_{T_{1}}}]^{T},$
(14)

where $T_{1}$ corresponds to the training time and $m=[1,...,n_{T_{1}}-1]$.
For fluid DoFs on $\Gamma_{1}^{f}$ and $\hat{\Omega}\cap\cal{S}$ we can
directly impose the corresponding Dirichlet boundary conditions, by setting
row $i$ of $A,H,K,C$ to zero, and assigning entries of value $1$ in
corresponding positions of $L$ and $B$, accordingly.

For the solid dynamics subsystem, we introduce sets $q_{FS}(i)$ and
$q_{I}(i)$, to distinguish between internal DoFs and DoFs on $\Gamma_{2}^{s}$.
As a result, $q_{FS}(i)=q_{f}(i)$ for DoFs $i$ on $\Gamma_{2}^{s}$ and
$q_{FS}(i)=\emptyset$ otherwise, while $q_{I}(i)=q_{s}(i)$ everywhere, except
for $\Gamma_{2}^{s}$, where $q_{I}(i)=\emptyset$. Then the variables of eq. 13
for the local inference of the solid dynamics operators are

$\boldsymbol{\beta}_{i}=[A_{S_{i,.}},K_{S_{i,.}},B_{S_{i,.}},H_{S_{i,.}}]^{T},\;\;\mathcal{D}_{i}^{T}=\begin{bmatrix}\partial_{t}\mathbf{d}_{s_{q_{I}(i)}}^{k}\\\
\mathbf{d}_{s_{q_{I}(i)}}^{k}\\\ \mathbf{u}_{q_{FS}(i)}^{k}\\\
{\mathbf{u}_{q_{FS}(i)}^{k}}^{2}\\\
\end{bmatrix},\;\;\mathbf{y}_{i}^{m+1}=[\partial_{t}{d}_{s_{i}}^{t=2}\;...\;\partial_{t}{d}_{s_{i}}^{t=n_{T_{1}}}]^{T}.$
(15)

In practice, eq. 13 should be complemented with a regularization term, to
bypass potential numerical errors due to small singular values of
$\mathcal{D}_{i}$. To this end, we employ truncated SVD of $\mathcal{D}_{i}$,
for different truncation thresholds $\eta$ of the normalized singular values
of $\mathcal{D}_{i}$. The optimal $\eta$ is determined as the “corner” of an
L-curve, with the solution norm $\left\|\boldsymbol{\beta}_{i}\right\|_{2}$ on
the $x$ axis and the error
$\left\|\boldsymbol{\beta}_{i}^{T}\mathcal{D}_{i}-\mathbf{y}_{i}^{m+1}\right\|_{2}$
on the $y$ axis (see Section 3.2). In essence, the SVD truncation corresponds
to an $L_{2}$ regularization of eq. 13. Regularization plays a crucial role on
the inferred model properties and thus determining the range of possible
$\eta$ is not trivial. The process of computing the optimal threshold $\eta$
is performed for few, randomly sampled DoFs. An interpolation from the
sampled, optimal $\eta$ values is performed to the rest of the mesh DoFs, as
shown in Section 3.2, such that only one least-squares problem per $i$ is
solved for those, with a predefined truncation threshold.

figureLeft: Truncation of the singular values of $\mathcal{D}$ to threshold
$\eta$, for a random DoF $i$. Right: Corresponding L-curve: The optimal $\eta$
achieves a compromise between the least-squares error
$\left\|\boldsymbol{\beta}_{i}^{T}\mathcal{D}_{i}-\mathbf{y}_{i}^{m+1}\right\|_{2}$
and the solution norm $\left\|\boldsymbol{\beta}_{i}\right\|_{2}$.

figureInterpolated, sub-optimal $\eta$ value for the streamwise flowfield
component, using the optimal $\eta$ for 30 % of the DoFs. Regularization
values are higher along the wake of the flow and close to the deformable solid
tail.

###  3.3 Model reduction and coupling

After solving eq. 13 for both the solid and the fluid nodes, we have inferred
the sparse, full-order matrices of eq. 8 and eq. 9. In order to derive a
reduced-order model, we use a POD projection. This is done by taking the SVD
of the velocity flowfield $\mathbf{u}$ snapshot matrix, and of the solid
velocity field $\partial_{t}\mathbf{d}_{s}$ snapshot matrix, over some
training interval $[0,T_{1}]$. By retaining the first $r_{f}$ and $r_{s}$
singular modes accordingly, we obtain the orthonormal projection bases
$\Phi_{f}\in\mathbb{R}^{2n_{f}\times r_{f}}$ and
$\Phi_{s}\in\mathbb{R}^{2n_{s}\times r_{s}}$. Substituting the projected
states $\mathbf{u}=\Phi_{f}\tilde{\mathbf{u}}$ and
$\partial_{t}\mathbf{d}_{s}=\Phi_{f}\partial_{t}\tilde{\mathbf{d}}_{s}$ to eq.
8 and eq. 9 results to a coupled, non-intrusive, reduced-order model

$\left\\{\begin{array}[]{ll}\begin{aligned}
\tilde{\mathbf{u}}^{k+1}&=\tilde{A}\tilde{\mathbf{u}}^{k}+\tilde{H}\tilde{\mathbf{u}}^{k+1}\otimes\tilde{\mathbf{u}}^{k+1}+\tilde{K}\;\partial_{t}\tilde{\mathbf{d}}_{s}\otimes\tilde{\mathbf{u}}^{k+1}+\tilde{B}\;\partial_{t}\mathbf{d}_{s}+\tilde{L}\;\mathbf{u}_{in}+\tilde{C}\\\
\partial_{t}\tilde{\mathbf{d}}_{s}^{k+1}&=\tilde{A}_{s}\partial_{t}\tilde{\mathbf{d}}_{s}^{k}+\tilde{K}_{s}\;\tilde{\mathbf{d}}_{s}^{k}+\tilde{B}_{s}\tilde{\mathbf{u}}^{k+1}+\tilde{H}_{s}\tilde{\mathbf{u}}^{k+1}\otimes\tilde{\mathbf{u}}^{k+1}\\\
\tilde{\mathbf{d}}_{s}^{k+1}&=\tilde{\mathbf{d}}_{s}^{k}+\frac{\Delta
t}{2}(\partial_{t}\tilde{\mathbf{d}}_{s}^{k+1}+\partial_{t}\tilde{\mathbf{d}}_{s}^{k}).\end{aligned}\end{array}\right.$
(16)

The corresponding ROM matrices of eq. 16 can be efficiently computed since the
sparsity patterns of the full-order matrices are known. We then have for the
fluid subsystem

$\tilde{A}=\Phi_{f}^{T}A\Phi_{f},\;\tilde{H}=\Phi_{f}^{T}H\Phi_{f}\otimes\Phi_{f},\;\tilde{K}=\Phi_{f}^{T}K(A_{L}S\Phi_{s})\otimes\Phi_{f},\;\tilde{B}=\Phi_{f}^{T}B\Phi_{s},\;\tilde{L}=\Phi_{f}^{T}L,\;\tilde{C}=\Phi_{f}^{T}C,$
(17)

and for the solid

$\tilde{A}_{s}=\Phi_{s}^{T}A_{s}\Phi_{s},\;\tilde{K}_{s}=\Phi_{s}^{T}K_{s}\Phi_{s},\;\tilde{B}_{s}=\Phi_{s}^{T}B_{s}R\Phi_{f},\;\tilde{H}_{s}=\Phi_{s}^{T}H_{s}(R\Phi_{f})\otimes(R\Phi_{f})\;.$
(18)

It was observed that the stiff coupling exhibited in FSI numerical simulation
[8] transfers also to the data-driven, reduced-order model. The implicit
formulation of eq. 8 is necessary for stability. The combined convergence
criterion for each timestep
${\left\|\tilde{\mathbf{u}}_{j+1}^{k}-\tilde{\mathbf{u}}_{j}^{k}\right\|}_{2}<res_{1}$
and
${\left\|\tilde{\mathbf{d_{s}}}_{j+1}^{k}-\tilde{\mathbf{d_{s}}}_{j}^{k}\right\|}_{2}<res_{2}$
was used during simulation, where $j$ denotes the index of iterations within
timestep $k$. For each implicit iteration, a successive under-relaxation
method was used.

##  4 Numerical Results

The analyzed methodology was tested for the Hron-Turek FSI3 benchmark [8].
Simulation data were obtained using the Gascoigne 3D open-source, Finite
Element solver [10]. A monolithic approach was used, with a moderate mesh of
$4128$ elements and a timestep of $\Delta t=5\times 10^{-3}\;s$, for a total
time of $5\;s$. We discard the first $1.5\;s$ of data to avoid the influence
of numerical solution initialization. The clock time required for the
numerical simulation is approximately $1.5$ hours on a laptop.

As presented in Section 3.1, we interpolate the data to a new mesh with $4060$
nodes, for which we construct the ALE map eq. 12. We use $65\%$ of the overall
available time series ($T_{1}=2.3\;s$) to solve eq. 13, computing the optimal
SVD truncation value for $10\%$ of the mesh nodes, using $10$, logarithmically
sorted values of $\eta\in[10^{-5},10^{-1}]$. We then construct the ROM through
eq. 17, eq. 18, by selecting $r_{f}=40$, $r_{s}=5$, and solve the implicitly
coupled, reduced-order system (16) by successive under-relaxation and
$res_{1}=res_{2}=10^{-6}$. The “offline” step of mesh construction,
interpolation and FOM inference requires approximately $5$ minutes, while the
simulation time of the ROM requires less than a minute.

The predicted, reconstructed ROM flowfield at $T=3.5\;s$, well beyond the
training time $T_{1}$, is given in Figure 1, in comparison to the CFD results
for the same instance, for both velocity components. It is observed that the
main dynamical features of vortex shedding as well as the motion of the
deformable tail are well approximated, using only $45$ reduced variables, from
the initial $8120$ DoFs. The vertical velocity over time for the upper right
corner of the deformable tail is also presented in Section 4. The predicted
ROM solid motion is well in line with the numerical results and captures the
limit cycle behaviour beyond training time $T_{1}$. A minor phase shift can be
attributed to the considerable presence of transient dynamics in the training
data.

Figure 1: (a) CFD results at $t=3.5\;s$ for streamwise velocity component (b)
ROM results ($r_{f}=40$, $r_{s}=5$) at $t=3.5\;s$ for streamwise velocity
component (c) CFD results at $t=3.5\;s$ for transverse velocity component (d)
ROM results ($r_{f}=40$, $r_{s}=5$) at $t=3.5\;s$ for transverse velocity
component

Finally, it is interesting to investigate the stability of the derived FSI
model with respect to the ROM dimensions $r_{s}$ and $r_{f}$ of the solid and
fluid subsystems, respectively. In Section 4, we illustrate the ROM streamwise
velocity error with respect to the numerical simulation on the reference
configuration $\hat{\Omega}$, averaged over the domain $\hat{\Omega}$ and
testing time $t=[T_{1},\;T]$. This is denoted with $\overline{e_{x}(t)}$.
Based on $\overline{e_{x}(t)}$, we observe that $r_{s}\geq 3$ modes are
required for the solid subsystem for an accurate FSI ROM, as anticipated by
the deformable tail motion. For values $r_{s}\geq 8$, a coupling instability
was exhibited, irrespective of the $r_{f}$ value. Similarly $r_{f}\geq 20$ SVD
modes are necessary on the fluid side. However, we also observe that the
combination of $r_{f}$ and $r_{s}$ can significantly affect the ROM predictive
accuracy, due to the strong coupling on the FSI interface. For a significant
range of ROM dimensions, the average error $\overline{e_{x}(t)}$ is below
$5\%$.

figureVelocity ROM prediction (streamwise, transverse) for the upper right
corner of the deformable solid. Training time $T_{1}=2.3\;s$ is denoted with a
dashed vertical line.

figureAverage error for streamwise flowfield over testing time
$\overline{e_{x}(t)}$, for different combinations of ROM dimensions
($r_{f},\;r_{s}$).

##  5 Conclusions

In this work, we employed an adjacency-based method for non-intrusive,
reduced-order modelling of FSI problems. Accurate results were obtained for
the Hron-Turek FSI3 benchmark testcase, with an average ROM prediction error
of less than 5 %. Based on the potential of this framework, aspects of
computational efficiency, parametric model reduction, as well as as a
theoretical analysis of local non-intrusive inference comprise potential,
future research directions.

## References

  * [1] Richter, T. (2017). Fluid-structure Interactions Models, Analysis and Finite Elements, Springer, Lecture Notes in Computational Science and Engineering, 118.
  * [2] Nonino, M., Ballarin, F. & Rozza, G. (2021). A monolithic and a partitioned, reduced basis method for fluid–structure interaction problems. Fluids, 6(6). https://doi.org/10.3390/fluids6060229
  * [3] Xiao, D., Yang, P., Fang, F., Xiang, J., Pain, C. C. & Navon, I. M. (2016). Non-intrusive reduced order modelling of fluid-structure interactions. Computer Methods in Applied Mechanics and Engineering. https://doi.org/10.1016/j.cma.2015.12.029
  * [4] Zastrow, B. G., Chaudhuri, A., Willcox, K. E., Ashley, A. S. & Henson, M. C. (2023). Data-driven Model Reduction via Operator Inference for Coupled Aeroelastic Flutter, American Institute of Aeronautics and Astronautics SCITECH 2023 Forum. https://doi.org/10.2514/6.2023-0330
  * [5] Schumann, Y. & Neumann, P. (2023). On linear models for discrete operator inference in time dependent problems. Journal of Computational and Applied Mathematics, 425. https://doi.org/10.1016/j.cam.2022.115022
  * [6] Baddoo, P. J., Herrmann, B., McKeon, B. J., Nathan Kutz, J. & Brunton, S. L. (2023). Physics-informed dynamic mode decomposition. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479(2271). https://doi.org/10.1098/rspa.2022.0576
  * [7] Gkimisis, L., Richter, T. & Benner, P. (2023). Adjacency-based, non-intrusive model reduction for Vortex-Induced Vibrations, arXiv Repository. https://doi.org/10.48550/arXiv.2302.09981
  * [8] Turek, S. & Hron, J. (2006). Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. Lecture Notes in Computational Science and Engineering, 53. https://doi.org/10.1007/3-540-34596-5_15
  * [9] Benner, P., Goyal, P., Heiland, J. & Duff, I. P. (2021). Operator inference and physics-informed learning of low-dimensional models for incompressible flows. Electronic Transactions on Numerical Analysis, 56. https://doi.org/10.1553/ETNA_VOL56S28
  * [10] Becker, R., Braack, M., Meidner, D., Richter, T. & Vexler, B. (2021). The finite element toolkit gascoigne (v1.01). Zenodo, doi:10.5281/zenodo.5574969