# Convergence Properties of Score-Based Models using Graduated Optimisation
for linear inverse problems

###### Abstract

The incorporation of generative models as regularisers within variational
formulations for inverse problems has proven effective across numerous image
reconstruction tasks. However, the resulting optimisation problem is often
non-convex and challenging to solve. In this work, we show that score-based
generative models (SGMs) can be used in a graduated optimisation framework to
solve inverse problems. We show that the resulting graduated non-convexity
flow converge to stationary points of the original problem and provide a
numerical convergence analysis of a 2D toy example. We further provide
experiments on computed tomography image reconstruction, where we show that
this framework is able to recover high-quality images, independent of the
initial value. The experiments highlight the potential of using SGMs in
graduated optimisation frameworks.

Index Terms—  graduated optimisation, score-based generative models,
optimisation, inverse problems

##  1 Introduction

Many problems in image reconstruction can be formulated as linear inverse
problems, where the goal is to recover an image ${\mathbf{x}}\in\mathcal{X}$
given noisy measurements ${\mathbf{y}}^{\delta}$ that are related via

$\displaystyle{\mathbf{y}}^{\delta}={\mathbf{A}}{\mathbf{x}}+\epsilon.$ (1)

Here ${\mathbf{A}}:\mathcal{X}\to\mathcal{Y}$ is a linear forward operator and
$\epsilon\in\mathcal{Y}$ the noise. Inverse problems are often ill-posed and
require regularisation to stabilise the reconstruction. Variational
regularisation is a popular framework to address ill-posedness [1]. It
formulates the reconstruction as an optimisation problem

$\displaystyle\hat{{\mathbf{x}}}\in\operatorname*{arg\,min}_{\mathbf{x}}\\{\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})+\alpha\mathcal{R}({\mathbf{x}})\\},$
(2)

where $\mathcal{D}:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}_{\geq 0}$
quantifies fitness to the measurements,
$\mathcal{R}:\mathcal{X}\to\mathbb{R}_{\geq 0}$ is a regulariser, and
$\alpha\geq 0$ balances the two terms. For additive Gaussian noise the datafit
is typically chosen as
$\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})=\frac{1}{2}\|{\mathbf{A}}{\mathbf{x}}-{\mathbf{y}}^{\delta}\|_{2}^{2}$.
Finding a suitable regulariser is a difficult task. In recent years a wide
variety of deep learning methods have been developed that aim to learn a
regulariser from a given dataset, see [2] for a review.

The statistical perspective on inverse problems identifies the regulariser
with the negative log-likelihood of a prior distribution $\pi({\mathbf{x}})$,
i.e. $\mathcal{R}({\mathbf{x}})=-\log\pi({\mathbf{x}})$. In this context,
learning a prior can be formulated as learning a parametrised distribution
$p_{\theta}({\mathbf{x}})$ as a proxy for $\pi({\mathbf{x}})$ [3]. If we have
access to the likelihood, we can consider the optimisation problem

$\displaystyle\hat{{\mathbf{x}}}\in\operatorname*{arg\,min}_{\mathbf{x}}\\{\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})-\log
p_{\theta}({\mathbf{x}})\\},$ (3)

as an instance of variational regularisation. Generative regularisers of this
form have been extensively studied for various inverse problems [3]. However,
as the negative log-likelihood of the generative model is often highly non-
convex, (3) is a challenging optimisation problem with many local minima and
it strongly depends on the initialisation.

We propose to tackle this problem using score-based generative models (SGMs)
[4] as regularisers. SGMs learn a sequence of gradually perturbed
distributions, starting at the data distribution and terminating in pure
noise. It was recently observed that this sequence can be used in annealed
Langevin sampling [5, 6] and graduated optimisation [7]. Graduated
optimisation [8], also known as continuation methods [9], is a heuristic for
dealing with non-convex problems in which the objective function is replaced
with a convex surrogate which can be solved efficiently. The surrogate
objective is then gradually transformed to the original non-convex objective.

In this work, we exploit the connection between SGMs and graduated
optimisation to solve the underlying non-convex optimisation problem and avoid
local minima. We showcase our algorithms on computed tomography.

##  2 Background

###  2.1 Graduated Optimisation

Graduated non-convexity is a heuristic global optimisation method for solving
non-convex minimisation problems, which creates a sequence of surrogate
optimisation problems that approximate the original problem through gradual
levels of smoothing or convexification [8]. Namely, for a non-convex function
$f:\mathcal{X}\to\mathbb{R}_{\geq 0}$ let
$F:\mathcal{X}\times[{t_{\text{min}}},{t_{\text{max}}}]\to\mathbb{R}_{\geq
0}$, such that $F({\mathbf{x}},{t_{\text{min}}})=f({\mathbf{x}})$ and that
$F({\mathbf{x}},{t_{\text{max}}})$ is convex with a unique minimiser. The
optimisation protocol consists of $I\in\mathbb{N}$ iterations over
${t_{\text{max}}}=t_{1}>\dots>t_{I}={t_{\text{min}}}$. In step $i$ we find the
minimiser ${\mathbf{x}}_{i}$ of $F({\mathbf{x}},t_{i})$, using
${\mathbf{x}}_{i-1}$ as the starting value. The process is terminated by
minimising the original function $f({\mathbf{x}})$ using ${\mathbf{x}}_{I-1}$
as initialisation.

The success of graduated optimisation strongly depends on the construction of
the embedding function $F({\mathbf{x}},t)$. Under specific constraints on the
type of embedding and the class of non-convex function under consideration, it
is possible to get global convergence results [10].

1: Initialise:
$\\!\displaystyle{\mathbf{x}}_{1}\\!\\!\in\\!\\!\mathbb{R}^{n},{t_{\text{max}}}\\!\\!=\\!t_{1}\\!\\!>\\!\\!\dots\\!\\!>\\!t_{I}\\!=\\!{t_{\text{min}}},I\\!\\!\in\\!\\!\mathbb{N},c\\!\\!\in\\!\\!(\\!0,\\!1\\!)$

2: for $i=1,\dots,I-1$ do

3: ${\mathbf{d}}_{i}=-t_{i}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},t_{i})$

4: Find: $\displaystyle\lambda_{i}$
s.t.$f({\mathbf{x}}_{i}+\\!\lambda_{i}{\mathbf{d}}_{i})\\!\leq\\!f({\mathbf{x}}_{i})\\!+\\!c\lambda_{i}\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}$

5:
$\displaystyle{\mathbf{x}}_{i+1}={\mathbf{x}}_{i}+\lambda_{i}{\mathbf{d}}_{i}$

6: end for

7: Output: $\displaystyle{\mathbf{x}}_{I}$

Algorithm 1 Graduated non-convexity flow with step size rule

1: Initialize: $\\!\displaystyle\
{\mathbf{x}}_{1}\\!\in\\!\mathbb{R}^{n},I\\!\in\\!\mathbb{N},c\\!\in\\!(0,\\!1\\!),\beta\\!\in\\!(0,\\!1\\!),\varepsilon\\!>\\!0,\
i\\!=\\!1$

2: while $1\leq i\leq
I\land\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i})\|>\varepsilon$ do

3: Determine:
$\displaystyle{\mathbf{d}}_{i}\in\mathbb{R}^{n}\quad\text{s.t.}\quad\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}<0$

4: Determine:
$\displaystyle\lambda_{i}\coloneqq\max\\{\beta^{\ell}|\ell=0,1,2,3,\dots\\}$
s.t. $f({\mathbf{x}}_{i}+\lambda_{i}{\mathbf{d}}_{i})\leq
f({\mathbf{x}}_{i})+c\lambda_{i}\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}$

5: ${\mathbf{x}}_{i+1}\leftarrow{\mathbf{x}}_{i}+\lambda_{i}{\mathbf{d}}_{i}$

6: $i\leftarrow i+1$

7: end while

Algorithm 2 Gradient-like Method

###  2.2 Score-based Generative Models

Score-based generative models (SGMs) have emerged as a powerful tool for
generative modelling [4, 11]. SGMs consist of a (predefined) forward process,
during which noise is gradually added; and a learnable reverse process, which
allows transforming noise into samples. The forward process is prescribed by
an Itô stochastic differential equation (SDE)

$\displaystyle
d{\mathbf{x}_{t}}={\mathbf{f}}({\mathbf{x}_{t}},t)dt+g(t)d\mathbf{w}_{t},\quad{\mathbf{x}}_{0}\sim
p_{0}:=\pi.$ (4)

The drift function ${\mathbf{f}}:\mathcal{X}\times\mathbb{R}\to\mathcal{X}$
and the diffusion function $g:\mathbb{R}\to\mathbb{R}$ are chosen such that
the terminal distribution approximates a Gaussian, i.e.
$p_{T}\approx\mathcal{N}(0,\mathbf{I})$. Under certain conditions there exists
a reverse diffusion process

$\displaystyle
d{\mathbf{x}}_{t}=\left[{\mathbf{f}}({\mathbf{x}}_{t},t)-g(t)^{2}\nabla_{{\mathbf{x}}_{t}}\log
p_{t}({\mathbf{x}}_{t})\right]dt+g(t)d\mathbf{w}_{t},$ (5)

which runs backwards in time [12]. SGMs approximate the score function
$\nabla_{{\mathbf{x}}_{t}}\log p_{t}({\mathbf{x}}_{t})$ by a neural network
$s_{\theta}({\mathbf{x}}_{t},t)$, which can be trained using denoising score
matching [13]

$\displaystyle\min_{\theta}\mathop{\mathbb{E}}_{\begin{subarray}{c}t\sim
U(0,1),\,{\mathbf{x}}_{0}\sim\pi\\\ {\mathbf{x}}_{t}\sim
p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})\end{subarray}}\left[\|s_{\theta}({\mathbf{x}}_{t},t)-\nabla_{{\mathbf{x}}_{t}}\log
p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})\|_{2}^{2}\right].$

For an affine drift ${\mathbf{f}}$, the transition kernel $p_{t|0}$ is given
as
$p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{t};\gamma_{t}{\mathbf{x}}_{0},\nu_{t}^{2}\mathbf{I})$
with parameters $\gamma_{t},\nu_{t}>0$ [14]. The SGM learns all intermediate
distributions, i.e.

$\displaystyle p_{\theta}({\mathbf{x}}_{t},t)\approx
p_{t}({\mathbf{x}}_{t})=\int\pi({\mathbf{x}}_{0})p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})d{\mathbf{x}}_{0}.$
(6)

Intermediate distributions are similar to Gaussian homotopy [9], as SDE
without a drift gives
$p_{t}(\cdot)=(\pi*\mathcal{N}(0,\nu_{t}^{2}\mathbf{I}))(\cdot)$.

##  3 Methods

###  3.1 Graduated non-convexity Flow

Using the SGM framework we define negative log-likelihoods

$\displaystyle\mathcal{R}({\mathbf{x}_{t}},t)\coloneqq-\log
p_{t}({\mathbf{x}_{t}}).$ (7)

We expect that $\mathcal{R}(\cdot,{t_{\text{min}}})$ is highly non-convex so
that the corresponding optimisation problem

$\displaystyle\min_{\mathbf{x}}\left\\{f({\mathbf{x}})\coloneqq\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})+\alpha_{{t_{\text{min}}}}\mathcal{R}({\mathbf{x}},{t_{\text{min}}})\right\\},$
(8)

is challenging to solve. However, due to the convolution with a Gaussian in
(6), there exists a $\tilde{t}<\infty$ such that
$\mathcal{R}(\cdot,\tilde{t})$ is convex for all $t>\tilde{t}$ [7]. Since
$\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})$ is convex, the
resulting optimisation problem is also convex at ${t_{\text{max}}}$ and can be
solved using common convex optimisation methods. This motivates using SGMs as
an embedding in a graduate optimisation scheme, giving a sequence of
optimisation problems

$\displaystyle\hat{{\mathbf{x}}}_{i+1}\\!\in\operatorname*{arg\,min}_{\mathbf{x}}\left\\{F({\mathbf{x}},t_{i})\\!\coloneqq\\!\mathcal{D}({\mathbf{A}}{\mathbf{x}},{\mathbf{y}}^{\delta})\\!+\\!\alpha_{t_{i}}\mathcal{R}({\mathbf{x}},t_{i})\right\\},\\!$
(9)

where $\alpha_{t_{i}}>0$ is a regularisation parameter and the optimisation is
initialised with $\hat{{\mathbf{x}}}_{i}$. Instead of exactly solving (9) in
each step, we perform one step of gradient descent,

$\begin{split}{\mathbf{x}}_{i+1}&={\mathbf{x}}_{i}-\lambda
t_{i}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},t_{i})\\\
&={\mathbf{x}}_{i}-\lambda
t_{i}{\mathbf{A}}^{*}({\mathbf{A}}{\mathbf{x}}_{i}-{\mathbf{y}}^{\delta})+\lambda
t_{i}\alpha_{t_{i}}s_{\theta}({\mathbf{x}}_{i},t_{i}),\end{split}$ (10)

where $s_{\theta}({\mathbf{x}}_{i},t_{i})\approx\nabla_{{\mathbf{x}}_{i}}\log
p_{t_{i}}({\mathbf{x}}_{i})$ is a pre-trained score-based model. This
corresponds to Algorithm 2 in [7] with a predefined smoothing schedule.

###  3.2 Gradient-like Methods

The critical issue in surrogate optimisation methods is ensuring the
convergence to an optimum of the original problem. For graduated non-convexity
flow in Section 3.1 this boils down to update directions used in (10), which
are not ensured to point in the direction of steepest descent. However,
convergence of iterations of the form
${\mathbf{x}}_{i+1}={\mathbf{x}}_{i}+\lambda_{i}{\mathbf{d}}_{i}$ can still be
shown provided ${\mathbf{d}}_{i}$ is a descent direction, i.e.
$\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}<0$.
Iterative methods of this type are also known as gradient-like methods [15,
Section 8.3, p. 75] or sufficient descent methods [16].

We consider the case when the maximal number of iterations $I$ (i.e. the
fineness of the smoothing schedule) tends to infinity. Let
$\\{t_{i}\\}_{i\in\mathbb{N}}\subset[{t_{\text{min}}},{t_{\text{max}}}]$ be
such that

$\displaystyle t_{1}={t_{\text{max}}},$ $\displaystyle t_{i+1}\leq
t_{i}\quad\forall i\in\mathbb{N},$
$\displaystyle\displaystyle\lim_{i\to\infty}t_{i}={t_{\text{min}}}.$ (11)

To show the convergence of $\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}$ we use
the following notion of gradient-like directions.

######  Definition 3.1 (Gradient-like Directions [15]).

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be continuously differentiable and
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}\subseteq\mathbb{R}^{n}.$ A sequence
$\\{{\mathbf{d}}_{i}\\}_{i\in\mathbb{N}}\subseteq\mathbb{R}^{n}$ is called
gradient-like with respect to $f$ and
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}},$ if for every subsequence
$\\{{\mathbf{x}}_{i_{j}}\\}_{j\in\mathbb{N}},$ converging to a non-stationary
point of $f,$ there exist $\varepsilon_{1},\varepsilon_{2}>0$ and
$N^{*}\in\mathbb{N}$ such that

  1. (a)

$\|{\mathbf{d}}_{i_{j}}\|\leq\varepsilon_{1},\ \forall j\in\mathbb{N}$ and

  2. (b)

$\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}}),{\mathbf{d}}_{i_{j}}\rangle_{\mathbb{R}^{n}}\leq-\varepsilon_{2},\
\text{for}\ j\geq N^{*}.$

Under these conditions every limit point of iterates produced by a gradient-
like algorithm is a stationary point of $f$ (see [15, Theorem 8.9]). In order
to show convergence we also require
$\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}},\cdot)$ to be Lipschitz
continuous, which ensures that search directions
$\\{{\mathbf{d}}_{i}\\}_{i\in\mathbb{N}}$ defined in (12) are gradient-like
with respect to $f$ in (8) and iterates
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}.$ With this, we have the following
convergence result.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 1: Convergence properties of the graduated non-convexity flow with
constant step size (Algorithm 1, Figures 1(a)-1(c)) and the gradient-like
method with the adaptive smoothing schedule (Algorithm 2, Figures 1(d)-1(f))
for the 2D toy example. We show the rate of trajectories, for different
choices of ${t_{\text{max}}}$, converging to the global minimum depending on
the initial starting point ${\mathbf{x}}_{1}$ (left column). Furthermore, we
display the rate of trajectories which converge to the global minimum (middle
column) and to local minima (right column) depending on the initial smoothing
parameter ${t_{\text{max}}}$ and the iteration number.

######  Theorem 3.2.

Take $F$ and $f$, given in (8) and (9), with $F(\cdot,t)$ continuously
differentiable for $t\in[{t_{\text{min}}},{t_{\text{max}}}]$ with
${t_{\text{min}}}>0$. Assume $\\{t_{i}\\}_{i\in\mathbb{N}}$ satisfy (11). Let
$\\{{\mathbf{d}}_{i}\\}_{i\in\mathbb{N}},\
\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}\subseteq\mathbb{R}^{n}$ be sequences
generated with Algorithm 2, where ${\mathbf{d}}_{i}$ is determined using the
adaptive smoothing schedule. Moreover, we assume that
$\alpha_{\cdot}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}},\cdot)$ is
Lipschitz continuous with a global Lipschitz constant $L\geq 0$ for all
${\mathbf{x}}\in\mathbb{R}^{n}$. Then every limit point of the sequence
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}$ is a stationary point of $f.$

In case of the graduated non-convexity flow for the function (9) and
iterations (10), the directions ${\mathbf{d}}_{i}$ correspond to

$\displaystyle{\mathbf{d}}_{i}$
$\displaystyle=-t_{i}\left({\mathbf{A}}^{*}({\mathbf{A}}{\mathbf{x}}_{i}-{\mathbf{y}}^{\delta})+\alpha_{t_{i}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i},t_{i})\right),$
(12)

which do not necessarily satisfy the descent condition

$\displaystyle\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}<0,$
(13)

needed for convergence. We can address this issue through an adaptive
smoothing schedule by selecting, in each iteration, the largest smoothing
parameter $t_{j^{\ast}}$ such that the descent condition is satisfied. Such an
index $j^{\ast}$ always exists as long as ${\mathbf{x}}_{i}$ is not a
stationary point of $f$, since the maximal index $j^{*}=I$ in Algorithm 2
leads to $t_{j^{*}}={t_{\text{min}}}$ and

$\Big{\langle}\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),-t_{j^{*}}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},t_{j^{*}})\Big{\rangle}_{\mathbb{R}^{n}}=-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i})\|^{2}<0.$

Theorem 3.2 can be shown using this adaptive smoothing schedule and by
ensuring the properties (a) and (b) in Definition 3.1 hold. The details are in
the Appendix.

###  3.3 Energy-based Parametrisation

Adaptive step-size iterative methods require evaluating the objective function
in each iteration to ensure convergence. However, traditional SGMs approximate
only the score function $\nabla_{\mathbf{x}}\log p_{t}({\mathbf{x}})$ and do
not offer access to the likelihood. The alternative are energy-based models
(EBM) [17], where the probabilistic model is parametrised as
$p_{\theta}({\mathbf{x}})=\frac{e^{g_{\theta}({\mathbf{x}})}}{Z(\theta)}$, for
a scalar-valued neural network $g_{\theta}$ and a normalisation constant
$Z(\theta)$. EBMs can be trained using the score matching objective by
defining
$s_{\theta}({\mathbf{x}}_{t},t)=\nabla_{{\mathbf{x}}_{t}}g_{\theta}({\mathbf{x}}_{t},t)$
[18], which can be evaluated using common automatic differentiation libraries.
This comes with a higher computational cost than the score parametrisation.
Following Du et al. [19], we parametrise the energy as

$\displaystyle
g_{\theta}({\mathbf{x}}_{t},t)=-\frac{1}{2}\|h_{\theta}({\mathbf{x}}_{t},t)\|_{2}^{2},$
(14)

where $h_{\theta}$ is implemented as a time-conditional U-Net [11]. The
perturbed cost function is then given as

$\displaystyle
F({\mathbf{x}},t)=\frac{1}{2}\|\mathbf{A}{\mathbf{x}}-{\mathbf{y}}^{\delta}\|_{2}^{2}+\frac{\alpha_{t}}{2}\|h_{\theta}({\mathbf{x}},t)\|_{2}^{2}.$
(15)

######  Remark 1.

Likelihood in the EBM framework is defined by

$\displaystyle\log
p_{\theta}({\mathbf{x}}_{t},t)=g_{\theta}({\mathbf{x}}_{t},t)-\log
Z(\theta,t).$ (16)

This introduces a time step dependency into the (intractable) normalisation
constant $Z(\theta,t)$. However, to compute the step sizes we only need to
evaluate the objective up to additive constants, which means that
$Z(\theta,t)$ does not need to be computed. Energy-based parametrisation was
used in a similar fashion to compute Metropolis correction probabilities [19].

##  4 Experiments

In this section we investigate the numerical performance of Algorithm 2. We
start with a toy example, where the data distribution is given by a Gaussian
mixture model and the score can be computed analytically. For the second
experiment we investigate computed tomography reconstruction on two datasets
with a trained SGM. We use the PSNR and SSIM [20] to evaluate the
reconstruction quality.

###  4.1 2D Toy Example

To illustrate the behaviour and the convergence properties, we consider a
Gaussian mixture model consisting of five Gaussians. The density at time step
$t$ is a Gaussian mixture model with a perturbed covariance matrix
$\Sigma_{k}^{t}=\Sigma_{k}+\frac{\sigma^{2t}-1}{2\log\sigma}\mathbf{I}$. The
diffusion process is given by the forward SDE
$d{\mathbf{x}}_{t}=\sigma^{t}d\mathbf{w}_{t}$ with perturbation kernel
$p_{t|0}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{t}|{\mathbf{x}}_{0},\frac{\sigma^{2t}-1}{2\log\sigma}\mathbf{I})$.
Furthermore, we consider a simple two dimensional inverse problem with the
forward operator $\mathbf{A}=\begin{pmatrix}1&1\\\ 0&0\end{pmatrix}$ and clean
measurements $\mathbf{y}=(2,0)^{\intercal}.$ We choose a constant
regularisation parameter $\alpha_{t}=5$ and adjust one mean of the Gaussian
mixture model to the position ${\mathbf{x}}^{*}=(1,1)^{\intercal}$ in order to
ensure that the global minimum of the cost function $f$ with
${t_{\text{min}}}=10^{-3}$ is at ${\mathbf{x}}^{*}.$ We evaluate the graduated
non-convexity flow (Algorithm 1) with a constant step size $\lambda_{i}=1$ as
well as the gradient-like method (Algorithm 2) with the adaptive smoothing
schedule. The values $t_{i}$ are evenly spaced between ${t_{\text{min}}}$ and
${t_{\text{max}}}.$ The goal is to analyse the algorithms in terms of their
convergence properties with respect to the initialisations ${\mathbf{x}}_{1},$
the initial smoothing parameter ${t_{\text{max}}}$ and the iteration number.
To do so, we run the algorithms with 1300 iterations for $10^{4}$ equally
spaced initial points ${\mathbf{x}}_{1}$ on $[-10,10]^{2}$ as well as 100
different values of ${t_{\text{max}}}\in[10^{-2},10]$, which are evenly
distributed on a log scale. The resulting rate of trajectories converging to
local minima and the global minimum are shown as a pseudocolor plot in Fig. 1.
The dependence on the initialisations is shown by the left column, where
isolines display the loss landscape of $f$ and the orange paths of iterates
show exemplary trajectories for different ${\mathbf{x}}_{1}$ and
${t_{\text{max}}}.$ The dependence on the iteration number and
${t_{\text{max}}}$ is shown by the middle and right column.

###  4.2 Computed Tomography

We evaluate our approach on two datasets: Ellipses [21] and AAPM [22].
Ellipses contains $128\times 128$px images of a varying number of ellipses
with random orientation and size. The AAPM dataset consists of CT images from
$10$ patients. We use images of $9$ patients to train the SGM, resulting in
$2824$ training images. We use the remaining patient (id C027) for evaluation.
We simulate measurements using a parallel-beam Radon transform with $60$
angles and use $10\%$ and $5\%$ relative Gaussian noise, respectively. For
both datasets, we train SGMs using the VPSDE [4]. We use a backtracking line
search using the Armijo-Goldstein condition to determine a suitable step size
with Barzilai-Borwein method [23] to find a candidate step size. We choose the
parameter as $\alpha_{t}=\alpha\nu_{t}/\gamma_{t}$, where $\nu_{t}$ is the
standard deviation and $\gamma_{t}$ the mean scaling of the perturbation
kernel. The value $\alpha$ was set according to a coarse sweep over one
example image. In Fig. 2 we show the result for the gradient like method on
one AAPM example, including the PSNR, objective value, step size and gradient-
like condition from (13). Fig. 3 shows example reconstructions for both
datasets.

Fig. 2: Result of the gradient like method for one example from the AAPM
dataset. We show the PSNR, the objective value, the gradient-like condition
and the computed step size during the iterations. The gradient-like condition
is always satisfied during the iterations.

#### Initial value

The gradient-like method (Alg. 1) is deterministic, except for the choice of
initialisation ${\mathbf{x}}_{1}$. We study the effect of ${\mathbf{x}}_{1}$
for both Ellipses and AAPM, see Table 1 by comparing Alg. 1 with gradient
descent for optimising the objective function at ${t_{\text{min}}}$, see (8).
For both datasets we use $300$ steps with a logarithmic time interval for Alg.
1. We evaluate $10$ random initialisations ${\mathbf{x}}_{1}$ for $10$ images.
We compute mean PSNR and SSIM values over all reconstructions. Further, we
compute the mean standard deviation value over images. We see that for both
Ellipses and AAPM the mean PSNR and mean SSIM are higher, while the mean
standard deviation is lower. That is, the gradient-like method achieves better
reconstructions more consistently.

Table 1: Mean and mean standard deviation of PSNR and SSIM for $10$ reconstructions on the Ellipses and AAPM dataset for $10$ different starting values. |  |  | Gradient Descent | Alg. 1
---|---|---|---|---
Ellipses | PSNR | mean | $28.96$ | $32.93$
std | $0.43$ | $0.10$
SSIM | mean | $0.771$ | $0.933$
std | $0.019$ | $0.001$
AAPM | PSNR | mean | $35.89$ | $37.33$
std | $0.032$ | $0.028$
SSIM | mean | $0.899$ | $0.913$
std | $0.032$ | $0.0004$
Table 2: PSNR and SSIM for $100$ reconstructions on the AAPM dataset for a constant step size ($1200$ iterations) and adaptive step size ($300$ iterations). | AAPM | Ellipses
---|---|---
| PSNR | SSIM | PSNR | SSIM
Constant step size | $32.93$ | $0.748$ | $30.82$ | $0.742$
Adaptive step size | $37.08$ | $0.896$ | $33.21$ | $0.935$

#### Adaptive step size

The energy-based parametrisation comes with a higher computational cost.
However, it allows uing adaptive step sizes for Alg. 2. In Table 2 we show
that this leads to a better reconstruction than a fixed step size. Compare
also the results in Fig. 2, which shows that the computed step size changes by
several orders of magnitude over the iterations. Also, we need fewer iteration
($300$ vs $1200$) to get convincing reconstructions.

Fig. 3: Example Reconstruction of Ellipses and AAPM. All images are shown in a
color range $[0,1]$.

##  5 Discussion

In the toy example in Section 4.1 it seems that the gradient-like condition in
Definition 3.1 hinders the ability of the algorithm to reach the global
optimum, as the search direction ${\mathbf{d}}_{i}$ always points in the same
halfspace as $-\nabla_{\mathbf{x}}f({\mathbf{x}}_{i})$. One possible approach
to circumvent this is to only enforce this condition after a set number of
iterations in order to ensure the convergence to a stationary point of $f$.
However, for the high-dimensional imaging experiments the gradient-like
condition does not seem to be too restrictive. As we see in Fig. 2, the
condition is satisfied for all iterations and we still achieve a good
reconstruction. Further, the adaptive step size rule leads to a improvement in
terms of both PSNR and SSIM, while using fewer iterations, which warrants the
additional computational cost of the step size search.

##  6 Conclusion

In this work we show that SGMs can be used in a graduated optimisation
framework to solve inverse problems. Further, we propose to use an energy-
based parametrisation, which enables the use of line search method for finding
a suitable step size. Using the theory of gradient-like directions, we can
prove convergence to stationary points. We hypothesise that this method is
also able to better escape saddle points and converge to local minima. This
research question is left for further work.

## References

  * [1] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen,  Variational methods in imaging, vol. 167,  Springer, 2009.
  * [2] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb,  “Solving inverse problems using data-driven models,”  Acta Numerica, vol. 28, pp. 1–174, 2019.
  * [3] M. Duff, N. Campbell, and M. J. Ehrhardt,  “Regularising inverse problems with generative machine learning models,”  J. Math. Imaging Vis., vol. 66, no. 1, pp. 37–56, 2024.
  * [4] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole,  “Score-based generative modeling through stochastic differential equations,”  in ICLR, 2021.
  * [5] Y. Song and S. Ermon,  “Generative modeling by estimating gradients of the data distribution,”  NeurIPS, vol. 32, 2019.
  * [6] Y. Sun, Z. Wu, Y. Chen, B. T. Feng, and K. L. Bouman,  “Provable probabilistic imaging using score-based generative priors,”  arXiv preprint arXiv:2310.10835, 2023.
  * [7] E. Kobler and T. Pock,  “Learning gradually non-convex image priors using score matching,”  arXiv preprint arXiv:2302.10502, 2023.
  * [8] A. Blake and A. Zisserman,  Visual reconstruction,  MIT press, 1987.
  * [9] H. Mobahi and J. W. Fisher,  “On the link between gaussian homotopy continuation and convex envelopes,”  in EMMCVPR 2015. Springer, 2015, pp. 43–56.
  * [10] E. Hazan, K. Y. Levy, and S. Shalev-Shwartz,  “On graduated optimization for stochastic non-convex problems,”  in ICML. PMLR, 2016, pp. 1833–1841.
  * [11] P. Dhariwal and A. Nichol,  “Diffusion models beat gans on image synthesis,”  NeurIPS, vol. 34, pp. 8780–8794, 2021.
  * [12] B. Anderson,  “Reverse-time diffusion equation models,”  Stoc. Proc. Appl., vol. 12, no. 3, pp. 313–326, 1982.
  * [13] P. Vincent,  “A connection between score matching and denoising autoencoders,”  Neural computation, vol. 23, no. 7, pp. 1661–1674, 2011.
  * [14] S. Särkkä and A. Solin,  Applied stochastic differential equations, vol. 10,  Cambridge University Press, 2019.
  * [15] C. Geiger and C. Kanzow,  Numerische Verfahren zur Lösung unrestringierter Optimierungsaufgaben,  Springer-Verlag, 2013.
  * [16] X. An, S. Li, and Y. Xiao,  “Sufficient descent directions in unconstrained optimization,”  Comput. Optim. Appl., vol. 48, no. 3, pp. 515–532, 2011.
  * [17] G. E. Hinton,  “Training products of experts by minimizing contrastive divergence,”  Neural computation, vol. 14, no. 8, pp. 1771–1800, 2002.
  * [18] T. Salimans and J. Ho,  “Should EBMs model the energy or the score?,”  in Energy Based Models Workshop - ICLR, 2021.
  * [19] Y. Du et al.,  “Reduce, reuse, recycle: Compositional generation with energy-based diffusion models and mcmc,”  in ICML. PMLR, 2023, pp. 8489–8510.
  * [20] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli,  “Image quality assessment: from error visibility to structural similarity,”  IEEE Trans. Image Process., vol. 13, no. 4, pp. 600–612, 2004.
  * [21] R. et al. Barbano,  “An educated warm start for deep image prior-based micro ct reconstruction,”  IEEE Trans. Comput. Imaging, vol. 8, pp. 1210–1222, 2022.
  * [22] C. H. McCollough et al.,  “Low-dose ct for the detection and classification of metastatic liver lesions: results of the 2016 low dose ct grand challenge,”  Medical physics, vol. 44, no. 10, pp. e339–e352, 2017.
  * [23] J. Barzilai and J. M. Borwein,  “Two-point step size gradient methods,”  IMA J. of numerical analysis, vol. 8, no. 1, pp. 141–148, 1988\.

## Appendix

### Proof of Theorem 3.2

######  Theorem 3.2.

Let $F$ and $f$ be given as in (8) and (9) and let $F(\cdot,t)$ be
continuously differentiable for $t\in[{t_{\text{min}}},{t_{\text{max}}}]$ with
${t_{\text{max}}}>{t_{\text{min}}}>0$. Furthermore, let
$\\{t_{i}\\}_{i\in\mathbb{N}}$ be as in (11) and
$\\{{\mathbf{d}}_{i}\\}_{i\in\mathbb{N}}\subseteq\mathbb{R}^{n}$ as well as
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}\subseteq\mathbb{R}^{n}$ be the
sequences generated based on the update rules of Algorithm 2, where
${\mathbf{d}}_{i}$ is determined by using the adaptive smoothing schedule.
Moreover, we assume that
$\alpha_{\cdot}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}},\cdot)$ is
Lipschitz continuous with a global Lipschitz constant $L\geq 0$ for all
${\mathbf{x}}\in\mathbb{R}^{n}$, i.e. there exists a constant $L\geq 0$ such
that

$\|\alpha_{t}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}},t)-\alpha_{\tilde{t}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}},\tilde{t})\|\leq
L|t-\tilde{t}|$

$\forall t,\tilde{t}\in[{t_{\text{min}}},{t_{\text{max}}}],\
\forall{\mathbf{x}}\in\mathbb{R}^{n}.$ Then every limit point of the sequence
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}$ is a stationary point of $f.$

###### Proof.

As described in Section 3.2, the adaptive smoothing schedule selects in the
$i$-th iteration of Algorithm 2 the largest smoothing parameter $t_{j^{*}}$
such that the descent condition in equation (13) with
${\mathbf{d}}_{i}=-\tilde{t}_{i}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},\tilde{t}_{i})$
and $\tilde{t}_{i}\coloneqq t_{j^{*}}$ holds. In other words, we choose
$j^{*}$ to be

$\displaystyle j^{*}\coloneqq\min\Big{\\{}j\in\mathbb{N}\ \Big{|}\
\Big{\langle}\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),-t_{j}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},$
$\displaystyle t_{j})\Big{\rangle}_{\mathbb{R}^{n}}<0$ $\displaystyle\land\
t_{j}\leq\tilde{t}_{i-1}\Big{\\}}.$

Such an index always exists as long as ${\mathbf{x}}_{i}$ is not a stationary
point of $f$ (which is ensured by Algorithm 2 due to the stopping criterion in
Line 2), as we have

$\displaystyle\lim_{j\to\infty}\Big{\langle}\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),-t_{j}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},t_{j})\Big{\rangle}_{\mathbb{R}^{n}}$
$\displaystyle=-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i})\|^{2}$
$\displaystyle<0,$

where the equation follows from the continuity of
$\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},\cdot).$ Hence, there exists a
$N\in\mathbb{N},$ such that

$\Big{\langle}\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),-t_{j}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i},t_{j})\Big{\rangle}_{\mathbb{R}^{n}}<0$

for all $j\geq N.$ The above choice of the ${\mathbf{d}}_{i}$ ensures that
$\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i}),{\mathbf{d}}_{i}\rangle_{\mathbb{R}^{n}}<0$
for all $i\in\mathbb{N},$ which is a needed criterion for the gradient-like
algorithm to converge according to [15]. The sequence
$\\{\tilde{t}_{i}\\}_{i\in\mathbb{N}}$ is a subsequence of
$\\{t_{i}\\}_{i\in\mathbb{N}}$ and fulfills the properties
$\tilde{t}_{i+1}\leq\tilde{t}_{i}$ as well as
$\lim_{i\to\infty}\tilde{t}_{i}={t_{\text{min}}}.$

Proof by contradiction: Let ${\mathbf{x}}^{*}\in\mathbb{R}^{n}$ be a non-
stationary point of $f$ and let $\\{{\mathbf{x}}_{i_{j}}\\}_{j\in\mathbb{N}}$
be a subsequence, which converges to ${\mathbf{x}}^{*}.$ The aim is to show
both properties (a) and (b) in Definition 3.1. We would then obtain that
$\\{{\mathbf{d}}_{i}\\}_{i\in\mathbb{N}}$ is gradient-like with respect to $f$
and $\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}.$ Hence, by [15, Theorem 8.9], we
would obtain that every accumulation point of
$\\{{\mathbf{x}}_{i}\\}_{i\in\mathbb{N}}$ is a stationary point of $f,$ which
would give the desired contradiction to the assumption that
$\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\neq 0.$

We first show that $\|{\mathbf{d}}_{i_{j}}\|\leq\varepsilon_{1}\ \forall
j\in\mathbb{N},$ which corresponds to property (a) of Definition 3.1. It holds
that

$\displaystyle\|{\mathbf{d}}_{i_{j}}\|$
$\displaystyle=\|-\tilde{t}_{i_{j}}\nabla_{\mathbf{x}}F({\mathbf{x}}_{i_{j}},\tilde{t}_{i_{j}})\|$
$\displaystyle=\Big{\|}\tilde{t}_{i_{j}}\left({\mathbf{A}}^{*}({\mathbf{A}}{\mathbf{x}}_{i_{j}}-{\mathbf{y}}^{\delta})+\alpha_{\tilde{t}_{i_{j}}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i_{j}},\tilde{t}_{i_{j}})\right)\Big{\|},$
$\displaystyle\leq{t_{\text{max}}}\Big{(}\|\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\|$
$\displaystyle\hskip
34.44434pt+\|\alpha_{\tilde{t}_{i_{j}}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i_{j}},\tilde{t}_{i_{j}})-\alpha_{{t_{\text{min}}}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i_{j}},{t_{\text{min}}})\|\Big{)}$
$\displaystyle\leq{t_{\text{max}}}\left(\|\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\|+L|\tilde{t}_{i_{j}}-{t_{\text{min}}}|\right)$
$\displaystyle\leq{t_{\text{max}}}\left(\|\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\|+L({t_{\text{max}}}-{t_{\text{min}}})\right).$

As $f$ is continuously differentiable and with
$\lim_{j\to\infty}{\mathbf{x}}_{i_{j}}={\mathbf{x}}^{*},$ we have that
$\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\to\nabla_{{\mathbf{x}}}f({\mathbf{x}}^{*})$
for $j\to\infty$ and the sequence
$\\{\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\\}_{j\in\mathbb{N}}$ is
bounded, i.e. there exists a constant $\tilde{c}\geq 0$ such that
$\|\nabla_{{\mathbf{x}}}f({\mathbf{x}}_{i_{j}})\|\leq\tilde{c}\ \ \forall
j\in\mathbb{N}$. Therefore, we have

$\displaystyle\|{\mathbf{d}}_{i_{j}}\|$
$\displaystyle\leq{t_{\text{max}}}\left(\tilde{c}+L({t_{\text{max}}}-{t_{\text{min}}})\right),$

which shows property (a) of Defintion 3.1. To prove property (b), we compute

$\displaystyle\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}}),{\mathbf{d}}_{i_{j}}\rangle_{\mathbb{R}^{n}}$
$\displaystyle=-\tilde{t}_{i_{j}}\langle\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}}),\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}})+\alpha_{\tilde{t}_{i_{j}}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i_{j}},\tilde{t}_{i_{j}})$
$\displaystyle\hskip
116.24963pt-\alpha_{{t_{\text{min}}}}\nabla_{\mathbf{x}}\mathcal{R}({\mathbf{x}}_{i_{j}},{t_{\text{min}}})\rangle_{\mathbb{R}^{n}}$
$\displaystyle\leq-\tilde{t}_{i_{j}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}})\|^{2}+\tilde{t}_{i_{j}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}})\|\cdot
L|\tilde{t}_{i_{j}}-{t_{\text{min}}}|$
$\displaystyle\leq-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}})\|^{2}+{t_{\text{max}}}\tilde{c}L|\tilde{t}_{i_{j}}-{t_{\text{min}}}|\eqqcolon
s_{i_{j}}.$

As $\nabla_{\mathbf{x}}f$ is continuous and
${\mathbf{x}}_{i_{j}}\to{\mathbf{x}}^{*}$ for $j\to\infty$ with
${\mathbf{x}}^{*}$ being a non-stationary point of $f,$ we have that

$\lim_{j\to\infty}\|\nabla_{\mathbf{x}}f({\mathbf{x}}_{i_{j}})\|^{2}=\|\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\|^{2}>0.$

The sequence $s_{i_{j}}$ converges and it follows

$\lim_{j\to\infty}s_{i_{j}}=-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\|^{2}<0.$

By definition of the convergence, it holds $\forall\varepsilon\\!>\\!0\
\exists N(\varepsilon)\geq 0:$

$|s_{i_{j}}-(-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\|^{2})|<\varepsilon\quad\forall
j\geq N(\varepsilon).$

By choosing
$\varepsilon\coloneqq{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\|^{2}/2>0,$
we obtain in particular

$s_{i_{j}}<-{t_{\text{min}}}\|\nabla_{\mathbf{x}}f({\mathbf{x}}^{*})\|^{2}/2<0\quad\forall
j\geq N(\varepsilon),$

which shows property (b) of Definition 3.1. ∎