Convergence Properties of Score-Based Models for Linear Inverse Problems Using Graduated Optimisation (2024)

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. The code is available¹¹1https://github.com/alexdenker/GradOpt-SGM.

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:endfor

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

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:endwhile

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 hom*otopy [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, givinga 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

(a)
$\|{\mathbf{d}}_{i_{j}}\|\leq\varepsilon_{1},\ \forall j\in\mathbb{N}$ and
(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.

Convergence Properties of Score-Based Models for Linear Inverse Problems Using Graduated Optimisation (1)

Convergence Properties of Score-Based Models for Linear Inverse Problems Using Graduated Optimisation (2)

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.

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 modelconsisting 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 stationary points 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.

Convergence Properties of Score-Based Models for Linear Inverse Problems Using Graduated Optimisation (7)

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 Table1 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.

			Gradient Descent	Alg.1
Ellipses	PSNR	mean	$28.96$	$32.93$
	PSNR	std	$0.43$	$0.10$
	SSIM	mean	$0.771$	$0.933$
	SSIM	std	$0.019$	$0.001$
AAPM	PSNR	mean	$35.89$	$37.33$
	PSNR	std	$0.032$	$0.028$
	SSIM	mean	$0.899$	$0.913$
	SSIM	std	$0.032$	$0.0004$

	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 Table2 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.

Convergence Properties of Score-Based Models for Linear Inverse Problems Using Graduated Optimisation (8)

5 Discussion

In the toy example in Section 4.1 it seems that the gradient-like condition in Definition3.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 learningmodels,”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 differentialequations,”in ICLR, 2021.
[5]Y.Song and S.Ermon,“Generative modeling by estimating gradients of the datadistribution,”NeurIPS, vol. 32, 2019.
[6]Y.Sun, Z.Wu, Y.Chen, B.T. Feng, and K.L. Bouman,“Provable probabilistic imaging using score-based generativepriors,”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 hom*otopy continuation and convexenvelopes,”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 unrestringierterOptimierungsaufgaben,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 contrastivedivergence,”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 etal.,“Reduce, reuse, recycle: Compositional generation with energy-baseddiffusion 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 structuralsimilarity,”IEEE Trans. Image Process., vol. 13, no. 4, pp. 600–612, 2004.
[21]R.Barbano and etal.,“An educated warm start for deep image prior-based micro ctreconstruction,”IEEE Trans. Comput. Imaging, vol. 8, pp. 1210–1222, 2022.
[22]C.H. McCollough etal.,“Low-dose ct for the detection and classification of metastaticliver 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 Theorem3.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}}}\|\eqqcolons%_{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.∎

	$\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).$