<figure><img src="image_1.png" /></figure>

HAL
open science

# A Morse-Bott Framework for Blind Inverse Problems
## Local Recovery Guarantees and the Failure of the MAP

Minh Hai Nguyen, Edouard Pauwels, Pierre Weiss

► To cite this version:

Minh Hai Nguyen, Edouard Pauwels, Pierre Weiss. A Morse-Bott Framework for Blind Inverse Problems Local Recovery Guarantees and the Failure of the MAP. 2026. hal-05533248

HAL Id: hal-05533248

https://hal.science/hal-05533248v1

Preprint submitted on 2 Mar 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Morse-Bott Framework for Blind Inverse Problems
## Local Recovery Guarantees and the Failure of the MAP *

Minh-Hai Nguyen †, Edouard Pauwels ‡, and Pierre Weiss †

**Abstract.** Maximum A Posteriori (MAP) estimation is a cornerstone framework for blind inverse problems, where an image and a forward operator are jointly estimated as the maximizers of a posterior distribution. In applications such as blind deblurring, this principle is used to recover sharp images from degraded observations. In this paper, we analyze the recovery guarantees of MAP-based methods by adopting a *Morse-Bott framework*. We model the image prior potential as a Morse-Bott function, where natural images are modeled as residing locally on a critical submanifold. This means that while the potential is locally flat along the “natural” directions of the image manifold, it is strictly convex in the directions normal to it. We demonstrate that this Morse-Bott hypothesis aligns with the structural properties of state-of-the-art learned priors, a finding we validate through an experimental analysis of the potential landscape and its Hessian spectrum. Our theoretical results show that, in a neighborhood of the ground-truth image and operator, the posterior admits local minimizers that are stable both with respect to initialization (gradient descent converge to the same minimizer) and to small perturbations of the data (solutions vary smoothly with the observations). This local stability explains the empirical success of well designed gradient-based optimization in these settings. However, we also demonstrate that this local stability is a **local** property: the “blurry trap”, well-known for sparse priors in blind deconvolution, persists even with state-of-the-art learned priors. Our findings demonstrate that the failure of MAP in blind deconvolution is not a limitation of prior quality, but an intrinsic characteristic of the landscape. We conclude that successful recovery of posterior maximization depends on strategic initialization within the basin of favorable local minima, and we validate this with numerical experiments on both synthetic and real-world data. The code for reproducing the numerical experiments in this paper is available at https://github.com/mh-nguyen712/blind_deblur_diffusion

**Key words.** blind deconvolution, maximum a posteriori, learned priors, Morse-Bott, local recovery

**MSC codes.** 68Q25, 68R10, 68U05

1. Introduction. Consider a parameterized linear forward operator $A(\theta) : \mathbb{R}^N \to \mathbb{R}^M$, where $\theta \in \mathbb{R}^P$ is a parameter vector describing the action of the operator. Blind inverse problems in imaging consist in recovering an image $\bar{x} \in \mathbb{R}^N$ and the parameter $\bar{\theta}$ from an indirect observation $y$ given by

$$ (1.1) \qquad y = A(\bar{\theta})(\bar{x}) + b $$

where $b$ represents additive noise. This problem is fundamental in areas such as biomedical imaging, astronomy, and computer vision. The parameter $\bar{\theta}$ typically describes operator properties such as point spread functions in optics, unknown projection angles in tomography, or sensitivity maps in MRI.

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** This work was funded by the ANR Micro-Blind ANR-21-CE48-0008 and the computational resources were provided by the IDRIS-CNRS under the allocation 2021-AD011012210R3.</footnote>

<footnote>†IRIT, Toulouse University, Toulouse, France (minh-hai.nguyen@irit.fr, pierre.armand.weiss@gmail.com).</footnote>

<footnote>‡Toulouse School of Economics, Université Toulouse Capitole et IUF, Toulouse, France (edouard2.pauwels@ut-capitole.fr).</footnote>
It is well-known that blind inverse problems are severely ill-posed, as there exist infinitely many pairs $(x, \theta)$ capable of explaining the observation $y$. Consequently, algorithms accompanied by theoretical recovery guarantees for $(\bar{x}, \bar{\theta})$ are scarce and often rely on restrictive or unrealistic assumptions [15, 17, 9, 21].

Taking a Bayesian viewpoint, we assume that $\bar{x}$ and $\bar{\theta}$ are realizations of two independent random vectors $x$ and $\theta$ with prior distributions $p_x$ and $p_\theta$. Using Bayes' formula, the posterior distribution is:

$$ (1.2) \qquad p(x, \theta|y) \propto p(y|x, \theta)p_x(x)p_\theta(\theta). $$

Under the assumption of white Gaussian noise $b \sim N(0, \sigma^2\text{Id})$, the negative-log-posterior reads:

$$ (1.3) \qquad L(x, \theta; y) = \frac{\|A(\theta)(x) - y\|_2^2}{2\sigma^2} - \log p_x(x) - \log p_\theta(\theta). $$

In all the paper, we assume that the parameter $\theta$ is uniformly distributed over a compact set $\Theta \subset \mathbb{R}^p$, leading to a constant prior $-\log p_\theta(\theta)$ that can be ignored in (1.3).

The Maximum A Posteriori (MAP) estimator seeks to minimize (1.3) with respect to the pair $(x, \theta)$. While handcrafted priors—promoting sparsity in spatial [27], wavelet [30], or gradient domains [6]—have historically dominated the field, recent advances in generative modeling [36, 33] have demonstrated unprecedented capability in capturing the complexity of natural images. Despite the power of these learned priors, their interaction with the MAP framework remains poorly understood. Specifically, the field lacks a rigorous explanation for why global MAP optimization leads to degenerate solutions, while local optimization – when properly initialized – remains remarkably effective.

Although the theoretical results derived hereafter extend to general blind inverse problems, we will focus on blind deconvolution in the experimental part. For this specific problem, popular approaches for describing the kernel prior include a direct modeling in the pixel domain with non-negativity and normalization constraints (simplex) [6, 26], more specific parametric models such as linear motion and out-of-focus blurs [25], diffraction-limited families [34, 10] or recent generative modeling [8]. In this setting, it has often been observed that global MAP optimization leads to degenerate solutions [20, 3], while local optimization with careful initialization may succeed.

**Contributions.** The progress in learned priors raise a fundamental question: *How does the geometry of learned prior landscapes shape the posterior in blind inverse problems, and can we formalize the conditions for local recovery versus global failure?* This paper addresses this above question through a rigorous theoretical and empirical analysis of the posterior landscape. Our contributions are as follows:

* We introduce a *Morse-Bott model* of image priors. We formalize the intuition that images lie on low-dimensional manifolds by modeling the prior potential as a Morse-Bott function, where the “natural image manifold” corresponds to a non-degenerate critical set.

* We provide *empirical evidence* showing that this hypothesis fits state-of-the-art diffusion priors, utilizing a novel analysis of the Hessian spectral gap to verify the non-degeneracy conditions required by the theory.
* We establish local recovery guarantees for general inverse and blind inverse problems. Under a Morse-Bott assumption on the prior, we prove that the ground-truth pair $(\bar{x}, \bar{\theta})$, with $\bar{x}$ a second-order critical point of the potential, acts as a stable local minimizer of the joint posterior, provided the forward operator satisfies certain injectivity conditions quantified through explicit stability constants. The stability notion is twofold: (i) with respect to initialization, in the sense that gradient descent starting in a neighborhood of $(\bar{x}, \bar{\theta})$ converges to the same minimizer, and (ii) with respect to small perturbations of the measurements $y$, in the sense that the recovered solution depends smoothly on the data.

* We provide numerical validation of these conditions for various blur kernels in the context of blind deconvolution.

* In the specific case of blind deconvolution, we identify the *blurry trap* as a fundamental limitation of global MAP. We validate empirically that for natural image distributions, learned priors consistently favor blurry images over sharp ones. This explains why the MAP coincides with a “no-blur” solution regardless of the prior’s generative quality.

* We propose an *informed initialization strategy* to navigate this landscape. By leveraging the local stability results, we demonstrate that successful recovery in blind deconvolution is a problem of basin selection, which can be achieved through a “large-to-small” kernel initialization.

**2. The Morse-Bott Landscape.** To analyze the convergence properties of blind deconvolution, we must formalize the geometry of the prior potential $q(x) = -\log p_x(x)$. Strong recovery guarantees for linear inverse problems have traditionally relied on convex, non-smooth priors—a framework instrumental to the field of compressed sensing [5, 39]. However, these handcrafted priors are increasingly recognized as insufficient for capturing the rich, high-dimensional complexity of natural images.

In this work, we pivot from convex to non-convex priors. We propose to see modern learned priors as *Morse-Bott potentials*. This framework accounts for the fact that natural images do not occupy the entire ambient space, but rather reside on a low-dimensional manifold [29].

**Definition 2.1 (Morse-Bott Potential).** Let $q: \mathbb{R}^N \to \mathbb{R}$ be a $C^2$ potential. We denote the set of its critical points as $C = \{\bar{x} \in \mathbb{R}^N | \nabla q(\bar{x}) = 0\}$. The potential $q$ is said to be Morse-Bott locally at $\bar{x} \in C$ if:

1. There exists a smooth submanifold $\mathcal{M}$ such that $C = \mathcal{M}$ locally around $\bar{x}$.

2. The Hessian $\nabla^2 q(\bar{x})$ is non-degenerate in the directions normal to $\mathcal{M}$. Equivalently, the kernel of the Hessian coincides exactly with the tangent space, $T_{\bar{x}}\mathcal{M}$, to $\mathcal{M}$ at $\bar{x}$:

$$ \ker(\nabla^2 q(\bar{x})) = T_{\bar{x}}\mathcal{M}. $$

The potential $q$ is called *Morse-Bott*, if it is *Morse-Bott* at any $\bar{x} \in C$.

Examples of Morse-Bott potentials include functions of the form $q(x) = \|Wx\|^2$ where $W$ is a rank-deficient matrix. More generally, given a smooth submanifold $\mathcal{M}$, the distance function to $\mathcal{M}$, $d_{\mathcal{M}}$ is Morse-Bott locally at each $\bar{x} \in \mathcal{M}$ [19, Problem 6-5]. Indeed, $q$ being locally Morse-Bott at a local minimizer $\bar{x} \in \mathcal{M}$ expresses roughly the fact that local minimizers of $q$ form smooth submanifold locally around $\bar{x}$, such that $x \mapsto q(x) - q(\bar{x})$ is of order $d_{\mathcal{M}}^2(x)$ up to positive constants [31].
Thus, the Morse–Bott assumption expresses a compatibility between the geometry of the minimizer set and the transverse curvature of the prior. In an imaging context, this agrees with the well established fact that natural images live on low-dimensional structures [29]. In the following sections, we explore how state-of-the-art diffusion-based priors align with the Morse-Bott property and provide additional insights on the potential’s landscape.

**2.1. Learned Priors exhibit the local Morse-Bott property.** To substantiate the Morse-Bott characterization of learned priors, we conduct an empirical analysis of the potential landscape *q*. We utilize state-of-the-art pretrained diffusion models from [36, 14] on the FFHQ, ImageNet, and AFHQ datasets at various resolutions, see Appendix A for more details. For simplicity, we refer to the potential *q* associated with these models as FFHQ-256, FFHQ-64, ImageNet-64, and AFHQ-64, respectively.

We shall verify two key properties: (i) the prior exhibits a low dimensional structure of local minimizers close to natural images, and (ii) there is a spectral gap similar to the Morse-Bott non-degeneracy condition. Taken together, these observations suggest that the Morse-Bott assumption is a relevant model for non-convex diffusion priors.

**2.1.1. Locating the local minimizers via Gradient Descent.** Following the established hypothesis that natural images reside on low-dimensional structures [29], we first investigate whether the diffusion prior potential *q* admits these structures as local minima. We perform gradient descent on the potential *q* using the score $\nabla q$ provided by the pre-trained diffusion models. The optimization is initialized from three distinct sources: natural images $x_i$ from the training set, images from unrelated datasets, and pure white Gaussian noise.

Trajectories starting from $x_i$ reach near-critical points $\bar{x}_i$ satisfying $\nabla q(\bar{x}_i) \approx 0$. Under mild smoothness assumptions, gradient descent avoids strict saddle points almost surely [18] and converges to second-order critical points where $\nabla^2 q(\bar{x}_i) \ge 0$. Figure 1 illustrates that these critical points $\bar{x}_i$ can be interpreted as “denoised” or “idealized” versions of the initial points. The monotonic decrease and stabilization of *q* in Figure 2 suggest that the optimization successfully reaches local minimizers.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Critical points obtained by a gradient descent starting from real images in FFHQ-256 using the FFHQ-256 potential. The critical points behave as denoised versions of the initial image. The potential $q = -\log p_x$ displayed on the top of each image is about twice lower for critical points.</figcaption></figure>

**2.1.2. Landscape Curvature.** In Figure 3, we report the relative distance $\|x_i - \bar{x}_i\|/\|\bar{x}_i\|$ and the corresponding variation in potential. We observe that a wide variety of initialization, including pure noise, reside in proximity (typically within 10%) to a critical point. However, the potential $q(x_i)$ for noise is orders of magnitude higher than for natural images.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the convergence of a gradient descent on the potential $q$ using different models. The potential stabilizes after a few thousands iterations.</figcaption></figure>

To characterize this behavior locally, let $d = \frac{x-\bar{x}}{\|x-\bar{x}\|}$ be the unit direction connecting an initialization to its associated critical point, and define the one-dimensional restriction

$$q^d(s) = q(\bar{x} + sd).$$

Since $\nabla q(\bar{x}) \approx 0$, a second-order Taylor expansion yields

$$q(x) - q(\bar{x}) \approx \frac{1}{2}\|x - \bar{x}\|^2 \langle \nabla^2 q(\bar{x}) d, d \rangle.$$

Equivalently, the curvature in direction $d$ is given by the Rayleigh quotient

$$(q^d)''(0) = \langle \nabla^2 q(\bar{x}) d, d \rangle \approx \frac{2(q(x) - q(\bar{x}))}{\|x - \bar{x}\|^2}.$$

Hence, the large increase observed in potential relative to the small distance $\|x - \bar{x}\|$ indicates a large curvature in directions transverse to the set of critical points. This behavior is consistent with a Morse-Bott structure, in which the potential is sharply curved in normal directions while remaining comparatively flat along the underlying low-dimensional manifold.

**2.1.3. Hessian Spectral Analysis.** The crucial part of the Morse-Bott property is the existence of a spectral gap in the Hessian $\nabla^2 q(\bar{x})$ whose kernel identifies the tangent space $T_{\bar{x}}\mathcal{M}$. We compute the full spectrum of the Hessian at several critical points $\bar{x}_i$. Since $\nabla q$ is parameterized by a neural network, we evaluate the Hessian-vector products $\nabla^2 q(\bar{x}) \cdot e_k$ via automatic differentiation. For a $3 \times 64 \times 64$ resolution image, the ambient dimension is $N = 12288$. We construct the dense $N \times N$ Hessian matrix as

$$\nabla^2 q(\bar{x}) = \begin{pmatrix} \nabla^2 q(\bar{x}) \cdot e_1 \\ \vdots \\ \nabla^2 q(\bar{x}) \cdot e_N \end{pmatrix} \in \mathbb{R}^{N \times N},$$

and compute its singular values in double precision.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distance and difference of potential between natural images and the corresponding critical points. The dataset used to train the model is highlighted in blue. Observe that the relative distance between the critical points and the initial images is less than 10% even when starting from pure noise. Yet, the potential is much lower.</figcaption></figure>

The resulting spectra, shown in Figure 4, reveal a clear separation of scales. Only a small portion of the eigenvalues are near-zero, while the vast majority are large and positive. This concentration of small eigenvalues is compatible with the hypothesis of a local image manifold of small intrinsic dimension $\mathcal{M}$ [29] and the non-degeneracy condition $\ker(\nabla^2 q(\bar{x})) = T_{\bar{x}}\mathcal{M}$. We shall see that this property ensures that the prior acts as a stable anchor for local recovery, despite the global complexities of the landscape.

**Remark 2.2.** Our approach to intrinsic dimensionality via Hessian spectra was carried out independently, and aligns with recent observations by Pidstrigach [28] and Stanczuk et al. [37], who similarly utilized score-based models to explore manifold structures. Consequently, the Morse-Bott property provides a formal theoretical framework which is aligned with these recent empirical observations.

**Remark 2.3 (On the choice of diffusion priors).** While our empirical analysis focuses on diffusion-based priors, the Morse-Bott framework and the recovery guarantees established in our results apply to any $\mathcal{C}^2$ potential. We choose diffusion models as our primary study case for two reasons: (i) they represent the current state-of-the-art in generative modeling of natural images; (ii) they provide direct access to the score function $\nabla q(x)$, enabling a precise numerical mapping of the Hessian spectrum. We expect the geometry of the posterior landscape to be a universal feature of high-performance priors trained on natural image manifolds.

**Remark 2.4 (Conservative fields).** In our experiments, the score $\nabla q$ is approximated using a pre-trained denoiser via Tweedie's formula [11]. The theoretical score $\nabla q(x)$ is a conservative vector field (i.e., its Jacobian is symmetric, and it integrates to a scalar potential). However, as noted in [32, 12], a generic neural network is not guaranteed to satisfy this integrability
condition perfectly, meaning it may not correspond exactly to the gradient of any function.

While some recent works propose architectures that are conservative by construction [12], we follow the standard approach in Score-Based Generative Modeling [36]: we rely on the fact that the optimal MMSE denoiser is inherently conservative. Consequently, we treat the non-integrability of the learned network as a negligible approximation error, assuming that the trained denoiser is sufficiently close to the true score of the distribution. To substantiate this claim, we conducted an experiment over 20 critical points $\bar{x}$, for which we evaluated the approximated score $\nabla q(\bar{x})$ and its Jacobian $\nabla^2 q(\bar{x})$. The lack of symmetry of the Jacobian can be measured with $\frac{\|\nabla^2 q(\bar{x}) - \nabla^2 q(\bar{x})^T\|_F}{\|\nabla^2 q(\bar{x})\|_F}$. The interquartile range of this quantity across the 20 critical points is [0.041, 0.047], indicating a modest deviation from integrability that we treat as a negligible approximation error.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spectra of the Hessian of the potential $\nabla^2 q(\bar{x}_i)$ for 20 different critical points $\bar{x}_i$. By thresholding all values below $10^{-5}$, the local image manifold dimension is estimated at 36 for FFHQ-64 and 16 for AFHQ-64.</figcaption></figure>

**3. Recovery Guarantees and Stability for Local Minimizers.** In this section, we establish the core theoretical result of this paper: the existence and stability of local minimizers in the posterior landscape. We show that under the Morse-Bott assumption, second order critical points of the prior remain stable anchors for local recovery. We start with non-blind inverse problems.

**3.1. Non blind inverse problems.** We first recall the notion of second order critical point.

**Definition 3.1 (Second order critical point).** A second order critical point of $q$ is a point $\bar{x}$ such that $\nabla q(\bar{x}) = 0$ and $\nabla^2 q(\bar{x}) \succeq 0$.

These points are the generic limits of gradient descents [18].

**Theorem 3.2 (Recovery guarantees – non-blind case).** Let $L(x; y) = \frac{1}{2\sigma^2} \|Ax - y\|^2 + q(x)$ and assume that the following conditions are satisfied:

1. $\bar{x}$ is a ground-truth input with $A\bar{x} = \bar{y}$,

2. $q: \mathbb{R}^N \to \mathbb{R}$ is a $C^2$ prior, locally Morse-Bott at $\bar{x}$ with local critical submanifold $\mathcal{M}$ and $\bar{x}$ is a second order critical point of $q$,

3. the identifiability condition $T_{\bar{x}}\mathcal{M} \cap \ker A = \{0\}$ holds.
219 Then $\nabla L(\bar{x}; \bar{y}) = 0$ and the Hessian $\bar{H} \stackrel{\text{def}}{=} \frac{1}{\sigma^2} A^T A + \nabla^2 q(\bar{x})$ of $L$ at $\bar{x}$ is positive definite:

220 $\bar{H} \succeq \mu I$ for some $\mu > 0$.

221 Therefore, there exist neighborhoods $U$ of $\bar{x}$ and $V$ of $\bar{y}$ such that for every $y \in V$:

222 1. $L(\cdot; y)$ admits a unique local minimizer $x(y) \in U$,

223 2. the map $y \mapsto x(y)$ is $C^1$,

224 3. the stability estimate holds:

225 $$ \|x(y) - \bar{x}\| \le \frac{\|A\|}{\mu\sigma^2} \|y - \bar{y}\| + o(\|y - \bar{y}\|). $$

226 The identifiability condition requires the operator $A$ to be injective when restricted to 227 the tangent space of the image manifold. It means that two distinct images on the manifold 228 near $\bar{x}$ cannot produce the same observation. Because $T_{\bar{x}}M$ is possibly low-dimensional, 229 this condition is significantly easier to satisfy than the global injectivity required in classical 230 inverse problems. The minimal eigenvalue $\mu$ in the theorem above is a quantitative measure 231 of this transversality condition. The larger $\mu$, the more robust to noise is the local minimizer. 232 We recall that a $C^2$-function $f$ is called locally strongly convex at $\bar{x}$, if $H_f(\bar{x}) \succeq \mu \text{Id}$ with 233 $\mu > 0$. Since $\bar{H}$ is positive definite, $L(\cdot; y)$ is locally strongly convex near $\bar{x}$. This implies 234 that a gradient descent initialized sufficiently close to $\bar{x}$ will converge to the local minimizer 235 at exponential rate [22, Thm 1.2.4], implying stability with respect to initialization.

236 **3.2. Blind inverse problems.** Theorem 3.2 is a special case of the following result which 237 describes the posterior landscape for blind linear inverse problems.

238 Theorem 3.3 (Stable recovery conditions – blind case). Consider $\theta \in \mathbb{R}^P \mapsto A(\theta) \in \mathbb{R}^{M \times N}$ 239 a $C^2$ parameterization of the linear operator. Let $L(x, \theta; y) \stackrel{\text{def}}{=} \frac{1}{2\sigma^2} \|A(\theta)x - y\|^2 + q(x)$ and 240 $J_x(\theta) = \frac{\partial}{\partial \theta}(A(\theta)x) \in \mathbb{R}^{M \times P}$ denote the Jacobian of the measurements w.r.t. $\theta$. Assume that:

241 • $(\bar{x}, \bar{\theta})$ is a ground-truth pair with $A(\bar{\theta})\bar{x} = \bar{y}$,

242 • $q: \mathbb{R}^N \to \mathbb{R}$ is $C^2$ and locally Morse-Bott at $\bar{x}$ with critical submanifold $M$ and $\bar{x}$ is a 243 second order critical point of $q$.

244 Then, $\nabla_{x,\theta}L(\bar{x}, \bar{\theta}; \bar{y}) = 0$ and the Hessian $\bar{H}$ of $L$ at $(\bar{x}, \bar{\theta})$

245 $$ (3.1) \qquad \bar{H} \stackrel{\text{def}}{=} \begin{pmatrix} \partial_{x,x}^2 L(\bar{x}, \bar{\theta}; \bar{y}) & \partial_{x,\theta}^2 L(\bar{x}, \bar{\theta}; \bar{y}) \\ \partial_{\theta,x}^2 L(\bar{x}, \bar{\theta}; \bar{y}) & \partial_{\theta,\theta}^2 L(\bar{x}, \bar{\theta}; \bar{y}) \end{pmatrix} = \begin{pmatrix} \bar{H}_{xx} & \bar{H}_{x\theta} \\ \bar{H}_{\theta x} & \bar{H}_{\theta\theta} \end{pmatrix}. $$

246 is positive definite if and only if the following geometric compatibility conditions hold:

247 1. Identifiability on the manifold: $T_{\bar{x}}M \cap \ker A(\bar{\theta}) = \{0\}$,

248 2. Identifiability of parameters: $\ker J_{\bar{x}}(\bar{\theta}) = \{0\}$,

249 3. Manifold-operator decoupling: $(A(\bar{\theta})T_{\bar{x}}M) \cap \text{Im } J_{\bar{x}}(\bar{\theta}) = \{0\}$.

250 **Corollary 3.4 (Stability bounds).** Under the conditions in Theorem 3.3, there exist neigh- 251 borhoods $U$ of $(\bar{x}, \bar{\theta})$ and $V$ of $\bar{y}$ such that for every $y \in V$:

252 1. $L(\cdot, \cdot; y)$ admits a unique local minimizer $(x(y), \theta(y)) \in U$,

253 2. the map $y \mapsto (x(y), \theta(y))$ is $C^1$,
3. the following first-order expansions hold:

$$x(y) - \bar{x} = \mathrm{Jac}_x(\bar{y})(y - \bar{y}) + o(\|y - \bar{y}\|) \quad \text{and} \quad \theta(y) - \bar{\theta} = \mathrm{Jac}_{\theta}(\bar{y})(y - \bar{y}) + o(\|y - \bar{y}\|).$$

The analytical expressions for the Jacobians are given by

$$(3.2) \qquad \mathrm{Jac}_x(\bar{y}) = \frac{1}{\sigma^2} S_x^{-1} (A(\bar{\theta})^T - \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T),$$

$$(3.3) \qquad \mathrm{Jac}_{\theta}(\bar{y}) = -\frac{1}{\sigma^2} \bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} (A(\bar{\theta})^T - \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T) + \frac{1}{\sigma^2} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T,$$

where $S_x \stackrel{\text{def}}{=} \bar{H}_{xx} - \bar{H}_{x\theta}\bar{H}_{\theta\theta}^{-1}\bar{H}_{\theta x}$ is the Schur complement of $\bar{H}$ with respect to $\bar{H}_{\theta\theta}$.

The proof of Theorem 3.3 and Corollary 3.4, provided in Appendix C.2, leverages the Morse-Bott non-degeneracy condition, $\bar{x}$ being a second-order critical point and the implicit function theorem to show that the joint Hessian is strictly positive-definite at $(\bar{x}, \bar{\theta})$.

**Remark 3.5 (Dependence on noise).** The first-order expansions provide a precise characterization of the local stability of the recovery. In particular, for bounded noise $\|y - \bar{y}\| \le \varepsilon$, we immediately obtain

$$\|x(y) - \bar{x}\| \le \sigma_{\max}(\mathrm{Jac}_x(\bar{y})) \varepsilon + o(\varepsilon),$$

$$\|\theta(y) - \bar{\theta}\| \le \sigma_{\max}(\mathrm{Jac}_{\theta}(\bar{y})) \varepsilon + o(\varepsilon)$$

where $\sigma_{\max}(\cdot)$ denotes maximal singular value, that coincides with the spectral norm. For stochastic noise $\mathbf{b} \sim \mathcal{N}(0, \sigma^2 \mathrm{Id})$, let $\mathbf{y} = \bar{y} + \mathbf{b} = A(\bar{\theta})\bar{x} + \mathbf{b}$ denote the random observation, in which case $\mathbb{E}[\|\mathbf{y} - \bar{y}\|^2] = M\sigma^2$. On the other hand, we can estimate the expected mean squared error (MSE) as:

$$\mathbb{E}[\|x(\mathbf{y}) - \bar{x}\|^2] \approx \sigma^2 \mathrm{Tr}(\mathrm{Jac}_x(\bar{y})^T \mathrm{Jac}_x(\bar{y})) = \sigma^2 \|\mathrm{Jac}_x(\bar{y})\|_F^2,$$

$$\mathbb{E}[\|\theta(\mathbf{y}) - \bar{\theta}\|^2] \approx \sigma^2 \mathrm{Tr}(\mathrm{Jac}_{\theta}(\bar{y})^T \mathrm{Jac}_{\theta}(\bar{y})) = \sigma^2 \|\mathrm{Jac}_{\theta}(\bar{y})\|_F^2.$$

Thus, while the spectral norm times the number of measurements gives a worst-case bound, the Frobenius norm provides a more realistic measure of the expected recovery error. In Subsection 5.2, we will use these analytical expressions to numerically evaluate the stability of various blur families, providing a concrete ranking of their identifiability on the image manifold.

**Remark 3.6 (On the size of the basin).** The radius of the neighborhood $V$ provided by the implicit function theorem can in principle be quantified in terms of the local regularity of the posterior. More precisely, if the Hessian $\nabla^2_{x,\theta}L$ is locally Lipschitz with constant $L_H$ near $(\bar{x}, \bar{\theta})$, then a valid basin size scales proportionally to $\mu/L_H$, where $\mu = \lambda_{\min}(\bar{H})$. Intuitively, this ratio describes how far one can move from the ground truth before the local strong convexity (quantified by $\mu$) is destroyed by the nonlinearity of the landscape (quantified by $L_H$). Since estimating the third-order derivatives required for $L_H$ is computationally prohibitive for neural priors and our focus is on geometric identifiability and local stability mechanisms, we do not pursue explicit global bounds here.
**A Bridge to Non-Convex Recovery.** Theorem 3.3 characterizes a class of images that can be recovered as stable local minima of the posterior. This result extends the recovery theories developed in compressed sensing [1, 15] to the domain of arbitrary, learned priors. While traditional results often require global convexity or range constraints [4], our Morse-Bott framework accommodates the complex, non-convex landscapes generated by modern diffusion models.

**3.3. Interpretation of the Geometric Conditions.** The local stability of the inverse problem relies on a geometric compatibility between the image manifold $\mathcal{M}$ and the forward operator $A(\theta)$. While condition (1) guarantees that the image $\bar{x}$ is uniquely recoverable when the operator is known, the blind setting in Theorem 3.3 introduces two additional requirements to resolve the joint ambiguity between the signal $x$ and the parameters $\theta$. We provide below an intuitive interpretation of these conditions, which collectively ensure that the prior $q$ and the operator $A$ interact transversely.

**Identifiability of Parameters (2).** This is a standard requirement for the identifiability of the forward operator. It ensures that variations in the parameter vector $\theta$ lead to measurable changes in the synthesized image. Without this, there would exist a continuous family of parameters indistinguishable from $\bar{\theta}$, even if the ground-truth image $\bar{x}$ were perfectly known.

**The Decoupling Condition (3).** This condition prevents the “cross-talk” between image variations and operator variations. Intuitively, it states that a change in the image along the manifold cannot be mistaken for a change in the operator’s parameters.

To see this, consider a first-order perturbation $(dx, d\theta)$ around $(\bar{x}, \bar{\theta})$. At first order, the synthesized image varies as:

$$A(\bar{\theta} + d\theta)(\bar{x} + dx) \approx \bar{y} + \underbrace{A(\bar{\theta})dx}_{\text{image change}} + \underbrace{J_{\bar{x}}(\bar{\theta})d\theta}_{\text{operator change}}$$

If the subspaces $A(\bar{\theta})T_{\bar{x}}\mathcal{M}$ and $\text{Im } J_{\bar{x}}(\bar{\theta})$ had a non-trivial intersection, we could find a non-zero pair $(dx, d\theta)$ such that these two terms cancel out $(Adx = -Jd\theta)$. In such a scenario, the synthesized image would remain constant while the potential $q$ (for $dx \in T_{\bar{x}}\mathcal{M}$) remains flat. In such a case, the posterior would contain a ridge of degenerate solutions, and strict local minimality would be lost. Condition (3) ensures that the signal subspace and the operator subspace are transverse, guaranteeing that the ground truth is an isolated, stable point in the landscape.

**4. The MAP Paradox: Global Bias and the “Blurry Trap”.** While Theorem 3.3 guarantees that the ground-truth $(\bar{x}, \bar{\theta})$ is a stable local minimizer, it does not ensure that it is the global one. In fact, a long-standing paradox in blind deconvolution is that global MAP optimization consistently favors a degenerate “no-blur” solution, where the operator $A(\theta)$ collapses to the identity and the image $x$ remains blurry. This phenomenon, well-known for handcrafted priors [20, 3], is often attributed to the inherent simplicity of sparsity-promoting models.

In this section, we show that this failure persists even with state-of-the-art diffusion priors. We demonstrate that the “blurry trap” is not a limitation of prior quality, but a fundamental structural property of the learned potential landscape $q$. Specifically, while the prior ad-
mits local minimizers around natural images, its global landscape is “slanted” toward blurry structures, effectively decoupling the global MAP estimator from the ground truth.

**4.1. Empirical Diagnostic: Blurry Images are More Likely.** The failure of the MAP estimate is rooted in a simple but disastrous property: the prior potential $q$ tends to decrease as the image becomes blurrier. To verify this for diffusion-based priors, we consider a collection of sharp images $x$ and their blurred counterparts $h_{\theta} \star x$, using the four blur families illustrated in Figure 5.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Different 1D blur families used in the forthcoming experiments. The parameter $\theta$ roughly accounts for the diameter of the PSF with $h_0 = \delta$ and the higher $\theta$, the more blur.</figcaption></figure>

*Monotonicity of the Potential.* Starting from a sharp image $x$ from the training set, we evaluate the potential $q(h_{\theta} \star x)$ as a function of the blur scale $\theta$. As shown in Figure 6 and qualitatively in Figure 7, the potential $q$ decreases strictly as $\theta$ increases across all models (FFHQ-256, ImageNet-64, AFHQ-64). This implies that $p_x(h_{\theta} \star x) > p_x(x)$, meaning the prior consistently assigns a higher likelihood to blurred versions of a natural image than to the sharp original. This behavior aligns with concurrent findings by Karczewski et al. [13], who observed that high-density regions of diffusion models gravitate toward smooth structures. Our results confirm that shifting from handcrafted to learned priors does not resolve the preference for blur; it merely encodes this bias into a more complex, non-convex landscape.

*A Specificity of Natural Image Statistics.* To determine if this preference for blur is an artifact of neural network architectures (e.g., a bias inherent to the U-Net), we trained an identical diffusion model on a binarized MNIST dataset. Unlike natural images, binarized digits consist of piecewise constant regions where blur creates unnatural intermediate intensities. As shown in Figure 8, the potential for the MNIST model *increases* with the blur level $\theta$. This experiment suggests that the “blurry trap” is not a bias of the neural architecture, but an inherent characteristic of models trained on the statistics of natural images. In natural photographs, smooth regions and defocus blur are prevalent, leading the learned prior to “accept” blurred structures as high-probability states.

**4.2. The Global Failure of MAP.** We now formalize how this global slant in the potential landscape leads to the systematic failure of the MAP estimator. Consider the negative log-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The evolution of the potential $q(h_\theta * x)$ on various models and 100 different images. From top-left to bottom-right: FFHQ-256, ImageNet-64, FFHQ-64 and AFHQ-64. The initial images are taken from the same training dataset as the pre-trained model. The potential is consistently decreasing with the blur level $\theta$, meaning that blurry images are more likely.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Variation of the potential $q$ with Gaussian blur on the FFHQ-256 model. A smaller potential means a more likely image. Observe that the potential is decreasing strictly with the blur level $\theta$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Variation of the potential $q$ with the blur level, for a model trained on a binarized MNIST dataset. The tendency of the potential to decrease with blur is not an architectural artifact, but depends on the statistics of natural images.</figcaption></figure>
posterior:

$$ (4.1) \qquad L(x, \theta; y) = \frac{1}{2\sigma^2} \|h_{\theta} \star x - y\|^2 + q(x). $$

In the following theorem, we provide a sufficient condition for the MAP estimator to collapse to the no-blur solution.

**Lemma 4.1 (Global MAP Failure).** Let $\mathcal{H} = \{h_{\theta}, \theta \in \Theta\}$ be the set of blur kernels, and assume the Dirac delta $\delta \in \mathcal{H}$. If the prior $q$ satisfies the blurry preference condition:

$$ (4.2) \qquad q(h \star x) \le q(x) \quad \forall x \in \mathbb{R}^N, \forall h \in \mathcal{H}, $$

then for any observation $y$, the set of global minimizers in $(x, \theta)$ of $L(x, \theta; y)$ contains a degenerate pair $(\hat{x}, \delta)$, where $\hat{x}$ is the MAP-denoised version of the blurry observation $y$.

**Lemma 4.1** demonstrates that the “Blurry Trap” is not a failure of the prior’s generative quality, but a topological consequence of the prior’s global slant. While Theorem 3.3 ensures that the ground truth $(\bar{x}, \bar{\theta})$ is a stable local minimum, Lemma 4.1 proves that if the blurry preference (4.2) holds globally, the ground truth can never be the global maximizer of the posterior.

Our numerical investigations in subsection 4.1 verify that condition (4.2) is satisfied by state-of-the-art diffusion priors. Consequently, blind deconvolution with MAP is characterized by competing basins of attraction: a favorable local basin around the ground truth, and a dominant global basin around the degenerate no-blur solution.

**Remark 4.2 (Local stability vs. global slant).** The local stability and the blurry preference (global failure) may coexist in the same landscape. As shown in our 1D posterior profiles Figure 9, the ground truth sits in a stable valley, but the “floor” of the image manifold $\mathcal{M}$ is tilted downward toward blurry states. Successful recovery thus depends not only on the existence of a solution, but also on an optimization strategy capable of selecting the correct local basin.

**Remark 4.3 (Local pits vs. global trend).** Recent work [13] suggests that diffusion priors may admit sharp, piece-wise constant (or “cartoon-like”) structures as local potential pits and report evidence showing that “Blurring an image increases its likelihood”. This work was conducted independently of our own, both came out roughly at the same time, independently observing this behavior. In fact [13] can be seen as an independent confirmation of the preference for blur. The preference for blur we report is a dominant statistical trend: on average, the potential floor slants toward blurry states. However, the prior’s landscape is locally non-monotonic and erratic, containing numerous sharp local minimizers with relatively small basins of attraction. While these “pits” allow the model to generate sharp samples, they do not counteract the global downward pull described by (4.2), which ultimately drives the global MAP estimator toward the degenerate no-blur solution.

**4.3. Extension to Other Blind Inverse Problems.** While we focused our experiments on blind deconvolution, the Morse-Bott framework and the “blurry trap” are applicable to any blind inverse problem where the family of operators $A(\theta)$ exhibits varying spectral decay. For instance, in blind super-resolution, the unknown anti-aliasing kernel plays a role identical to the blur kernel here; a prior preference for smoothness will systematically bias the estimator
toward over-smoothed kernels. Similarly, in blind MRI motion correction or calibration, any parameterization that allows the forward mapping to act as a low-pass filter will likely trigger the same degenerate global MAP behavior. In all these cases, the identifiability conditions in Theorem 3.3 remain the primary tool for guaranteeing that a stable, sharp solution exists locally, even when it is not the global maximizer. This highlights that the fundamental challenge of blind inverse problems is the topological inconsistency of the MAP estimator.

**5. Numerical Experiments.** The theoretical duality established in the previous sections suggests that the posterior landscape is characterized by the coexistence of a stable recovery basin (Section 3) and a dominant, degenerate global attractor (Section 4). To visualize this competition, we explore the joint posterior landscape along the univariate blur scale $\theta \in [0, \theta_{max}]$. Specifically, we analyze the marginalized negative log-posterior profile, or the lower envelope of the landscape:

$$ (5.1) \qquad \ell_y^{\text{opt}}(\theta) \stackrel{\text{def}}{=} \inf_{x \in \mathbb{R}^N} L(x, \theta; y). $$

By optimizing over the image $x$ for a fixed kernel parameter $\theta$, we reveal the basins of attraction that determine the success or failure of joint MAP estimation.

**5.1. Noiseless Measurements: The Coexistence of Basins.** We first consider the ideal noiseless scenario, $\bar{y} = A(\bar{\theta})\bar{x}$, where $A(\theta)x = h_{\theta} \star x$. To facilitate a continuous mapping of the landscape, we ensure the invertibility of the operator for all $\theta$ by considering slightly perturbed kernels $\tilde{h}_{\theta} = h_{\theta} + \lambda\delta$, where $\lambda > 0$ is a small regularization constant. In this regime, the marginalized profile can be computed explicitly via inverse filtering: $\hat{x}(\theta) = A(\theta)^{-1}\bar{y}$, which implies $\ell_{\bar{y}}^{\text{opt}}(\theta) = q(A(\theta)^{-1}\bar{y})$. We use the implementation from the DeepInverse library [38] for efficient inverse filtering.

As shown in Figure 9, the marginalized profiles for the FFHQ-256 model provide a striking visualization of the MAP paradox:

* The Blurry Trap (global failure): Confirming Lemma 4.1, the global minimum of the profile occurs at $\theta = 0$ (the Dirac kernel). This demonstrates that the no-blur solution is the global attractor of the landscape, regardless of the prior's generative quality.

* The Morse-Bott valley (local stability): Centering the analysis on a critical point $\bar{x}$ reveals a distinct local minimum at the ground-truth parameter $\bar{\theta}$. This validates Theorem 3.3, proving that the “correct” solution sits within a stable, non-degenerate recovery basin.

**Geometric Identifiability and Fourier Decay.** Interestingly, the profiles for motion and defocus blurs exhibit much sharper local minima than those for Gaussian or Airy blurs. In the context of Theorem 3.3, this suggests that motion and defocus operators satisfy the manifold-operator decoupling condition (3) with higher sensitivity. It can be seen by smaller values of $\|Jac_{\theta}(\bar{y})\|_{F}$ and $\sigma_{\max}(Jac_{\theta}(\bar{y}))$ for the motion kernel compared to the Gaussian kernel, that we numerically validated in Table 1. Because the Fourier transforms of these kernels vanish more slowly, the operator $A(\bar{\theta})$ remains well-conditioned on the tangent space $T_{\bar{x}}\mathcal{M}$, leading to a more pronounced basin of attraction and higher stability against perturbations.

**5.2. Numerical Verification of Geometric Identifiability.** The exact recovery guarantees in Theorem 3.3 rely on the three geometric transversality conditions 1, 2 and 3, or equivalently,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Marginalized posterior profiles $l_{\bar{y}}^{\text{opt}}(\theta)$ for noiseless measurements. Left: $\bar{x}$ is a natural image from the dataset. Right: $\bar{x}$ is a critical point of the prior. The emergence of a local minimum at $\theta = 0$ for critical points validates the Morse-Bott stability theory (Theorem 3.3), while the global minimum at $\theta = 0$ validates the global failure (Lemma 4.1).</figcaption></figure>

to the invertibility of the Hessian $\bar{H}$ (3.1). By Schur's formula, the invertibility of $\bar{H}$ is equivalent to the invertibility of the Schur complement $S_x$ and of $\bar{H}_{\theta\theta}$. These conditions are verified in all the experiments below. However, the mere invertibility of this matrix is qualitative. We substantiate this invertibility by computing the stability bounds given in Remark 3.5, which provide more meaningful information. We evaluate these values using the FFHQ-64 pretrained model. The ground truth images $\bar{x}$ are taken as critical points of the potential, obtained as in Subsection 2.1.1, and we evaluate the stability bounds using different blur kernels displayed in Table 1:

1. *Isotropic Gaussian blur*: a standard isotropic Gaussian kernel ($\bar{\theta} = 1.0, P = 1$), known to be ill-conditioned.

2. *Motion blur*: a generic motion kernel (length $\bar{\theta} = 5, P = 1$), which preserves more high-frequency information.

3. *Diffraction blur*: a diffraction blur with 5 first aberrations ($P = 5$, with the same amplitude of 0.15 for all aberrations).

4. *Blur in the simplex*: we also consider a (randomly generated) blur kernel living in the simplex of dimension $15 \times 15$ ($P = 225$).

We build the Jacobian matrices $Jac_x(\bar{y})$ and $Jac_\theta(\bar{y})$ (see (3.2) and (3.3)), in double precision floating point, using PyTorch.

Based on Remark 3.5, we then estimate the following quantities:

$$ \varepsilon_{\max,x}(\bar{y}) \stackrel{\text{def}}{=} \frac{\sigma_{\max}(Jac_x(\bar{y}))\sigma\sqrt{M}}{\|\bar{x}\|} \quad \text{and} \quad \varepsilon_{\text{avg},x}(\bar{y}) \stackrel{\text{def}}{=} \frac{\sigma\|Jac_x(\bar{y})\|_F}{\|\bar{x}\|} $$

$$ \varepsilon_{\max,\theta}(\bar{y}) \stackrel{\text{def}}{=} \frac{\sigma_{\max}(Jac_{\theta}(\bar{y}))\sigma\sqrt{M}}{\|\bar{\theta}\|} \quad \text{and} \quad \varepsilon_{\text{avg},\theta}(\bar{y}) \stackrel{\text{def}}{=} \frac{\sigma\|Jac_{\theta}(\bar{y})\|_F}{\|\bar{\theta}\|}, $$

which quantify the worst-case and average relative error in parameter recovery for additive white Gaussian noise $b \sim N(\sigma^2, \text{Id})$. The interquartile ranges are displayed in Table 1 to show the typical range for these values. As observed in Table 1, the distinction made in Remark 3.5 between worst-case bounded noise and stochastic noise is critical. For instance, the blur in the simplex presents a high spectral norm bound for the image $([1.13, 3.68] \times 10^1)$, suggesting
<caption>Table 1</caption>

Numerical verification of geometric conditions. The interquartile range $[Q_1, Q_3]$ for the maximal and typical relative errors over 20 random critical images $\bar{x}$ (computed by gradient descent, see Subsection 2.1.1) are presented. The parameterization of the PSF has size $P$ and the noise level is set to $\sigma = 0.01$. The images are normalized ensuring that $\bar{x} \in [-1, 1]^N$. The dimension of the nullspace of $A(\bar{\theta})$ is estimated by thresholding the singular values of $A(\bar{\theta})$ at $10^{-3}$.

<table><thead><tr><th>Metric</th><th>Gaussian<br>(P = 1)</th><th>Motion<br>(P = 1)</th><th>Diffraction<br>(P = 5)</th><th>Blur in the simplex<br>(P = 225)</th></tr></thead><tbody><tr><td>Blur kernel</td><td></td><td></td><td></td><td></td></tr><tr><td>dim(kerA(θ))</td><td>567</td><td>0</td><td>1932</td><td>6</td></tr><tr><td>ε<sub>max,x</sub></td><td>[0.43, 1.39] × 10<sup>2</sup></td><td>[1.77, 2.93] × 10<sup>0</sup></td><td>[3.07, 9.68] × 10<sup>1</sup></td><td>[1.13, 3.68] × 10<sup>1</sup></td></tr><tr><td>ε<sub>max,θ</sub></td><td>[1.27, 2.83] × 10<sup>1</sup></td><td>[3.32, 6.16] × 10<sup>-1</sup></td><td>[1.60, 2.53] × 10<sup>1</sup></td><td>[3.52, 11.0] × 10<sup>1</sup></td></tr><tr><td>ε<sub>avg,x</sub></td><td>[2.03, 3.42] × 10<sup>0</sup></td><td>[2.47, 3.37] × 10<sup>-1</sup></td><td>[3.96, 6.18] × 10<sup>0</sup></td><td>[2.38, 4.71] × 10<sup>-1</sup></td></tr><tr><td>ε<sub>avg,θ</sub></td><td>[1.15, 2.55] × 10<sup>-1</sup></td><td>[2.99, 5.53] × 10<sup>-3</sup></td><td>[2.13, 3.37] × 10<sup>-1</sup></td><td>[0.55, 1.45] × 10<sup>0</sup></td></tr></tbody></table>

severe vulnerability to adversarial perturbations. However, its expected MSE remains low ($[2.38, 4.71] \times 10^{-1}$), indicating that stochastic noise rarely aligns with the operator's most sensitive singular vectors.

The results also highlight a stark contrast in the conditioning of different blur operators. The Gaussian kernel and diffraction blur kernels exhibit the most severe ill-posedness, yielding the highest worst-case sensitivity and expected mean squared error for the image reconstruction. Conversely, the motion blur preserves more structural high-frequency information, resulting in an expected MSE that is orders of magnitude lower. Finally, the recovery of the blur parameters $\theta$ is significantly more stable than that of the image $x$, except for the motion blur. This is likely because the parameter dimension $P$ is significantly smaller than the image dimension $N$. For motion blur in the simplex, the kernel lives in a much higher dimension, which may explain a worse identifiability. In addition, for this case, the fundamental ambiguity in blind deconvolution, where the image and blur are shifted in opposite directions, leaving the blurry image unchanged, remains. Yet, the inversibility of the Hessian shows that this ambiguity is mostly solved by the prior. These experiments overall suggest that the parameter can be recovered rather reliably, while the reconstruction uncertainty on the image $x$ remains large: the prior's geometry is not able to compensate for the loss of information.

**5.3. Noisy Measurements and Basin Collapse.** In the presence of noise ($\sigma > 0$), the marginalized profile $\ell_y^{\text{opt}}(\theta)$ is no longer explicitly available. We approximate it by minimizing $x \mapsto L(x, \theta; y)$ using a momentum-based gradient descent initialized from multiple points, retaining the endpoint with the lowest cost value.

Figure 10 illustrates the evolution of the marginalized profile across increasing noise levels. At low noise levels, the posterior maintains a clear local recovery basin near $\bar{\theta}$, particularly
when the image component is close to a critical point of the prior. However, as the noise increases, we observe a *basin collapse*: the local recovery minimum is absorbed by the global downward pull of the "blurry trap." This phenomenon suggests the existence of a *critical*

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evolution of the marginalized profile $\ell_y^{\text{opt}}(\theta)$ under noisy measurements with natural images (left) and critical points (right). As the noise level $\sigma$ increases, the local recovery basin near $\bar{\theta}$ vanishes, illustrating the stability limits of the Morse-Bott landscape.</figcaption></figure>

noise threshold beyond which the Morse-Bott structure of the recovery basin is destroyed. For noise levels above this threshold, the ground truth is no longer a local minimizer, and the MAP estimator becomes fundamentally inconsistent. These results underscore that blind deconvolution is not merely a problem of finding the global minimum, but of ensuring that a favorable local basin exists and is successfully selected by the optimization algorithm.

**6. Navigating the Posterior: Initialization and Optimization.** The theoretical and numerical mapping of the posterior in the preceding sections reveals that blind deconvolution is essentially a problem of competing attractors. While the ground-truth pair $(\bar{x}, \bar{\theta})$ resides within a stable Morse-Bott recovery basin, the global landscape is dominated by a "blurry trap" that pulls the optimizer toward a degenerate Dirac solution. Successful recovery thus hinges on a *basin selection strategy*: the optimization must be initialized within the recovery basin and shielded from the global downward pull.

To navigate this landscape and validate the existence of the recovery basin, we propose a *heuristic* alternating optimization scheme, detailed in Algorithm 6.1. This solver is not designed to guarantee global convergence, but rather to exploit the topological insights of the Morse-Bott framework using three guiding principles:

1. *Basin-aware initialization:* We initialize $x_0 = y$ and $\theta_0$ such that $h_{\theta_0}$ is the “maximal” kernel in the parameter space (e.g., maximal spatial support for simplex kernels or high aberration coefficients for Zernike models). When initialized with a small kernel (near the Dirac spike), the data fidelity term $\|h_{\theta} \star x - y\|^2$ is minimized by an image $x$ that remains close to the blurry observation $y$. In this regime, the optimization stays
<caption>**Algorithm 6.1** Blind deconvolution with alternating optimization</caption>

1: **Input:** blurry image $y$, number of alternating iterations $K$, number of iteration for image update $K_x$, step-size $\gamma_x, \gamma_\theta > 0$.
2: Initialize $\theta_0$ and $x_0 = y$
3: **for** $k = 0, 1, \dots, K-1$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;$$x_{k+1,0} = \begin{cases} y & \text{if } \mod(k, 100) = 0 \\ x_{k,K_x} & \text{otherwise} \end{cases}$$
&nbsp;&nbsp;&nbsp;&nbsp;Minimization of the posterior with respect to the image
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 0, 1, \dots, K_x - 1$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$x_{k+1,i+1} = \text{prox}_{\gamma_x f(\cdot, \theta_k)}(x_{k+1,i} - \gamma_x \nabla q(x_{k+1,i}))$$
7: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
8: &nbsp;&nbsp;&nbsp;&nbsp;$$x_{k+1} = x_{k+1,K_x}$$
&nbsp;&nbsp;&nbsp;&nbsp;A gradient descent step for the kernel
9: &nbsp;&nbsp;&nbsp;&nbsp;$$\theta_{k+1} = \theta_k - \gamma_\theta \nabla_\theta f(x_{k+1}, \theta_k)$$
10: **end for**
11: **return** $x_K, \theta_K$

in the “low-frequency” part of the image manifold $M$, where the potential $q$ is most slanted toward the blurry trap. By contrast, a large kernel initialization forces the image component $x$ to contain sharp features and high-frequency content to satisfy the data fidelity.

2. **Asymmetric alternating optimization:** We perform multiple proximal gradient steps on the image $x$ for every single gradient step on the kernel $\theta$. This ensures that the image component remains close to the “manifold floor” of the Morse-Bott potential, where the stability conditions of Theorem 3.3 are strictly satisfied.

3. **Periodic reset (Heuristic escape):** We reset $x_k = y$ periodically (e.g., every 100 iterations). While this step breaks the strict variational structure of the algorithm (disrupting the monotonic decrease of the objective), it effectively “jumps” the optimization back into the data-consistent space. This heuristic serves to reset the data-consistency of the iterate, preventing stagnation in the small, erratic local pits that characterize learned prior landscapes (see Remark 4.3), allowing it to continue migrating toward the stable basin.

In Algorithm 6.1, $f(x, \theta) = \frac{1}{2\sigma^2} \|h_\theta \star x - y\|^2$ denotes the data fidelity term. The proximal operator $\text{prox}_{\gamma_x f(\cdot, \theta)}$ corresponds to a linear system resolution, which we compute efficiently using the conjugate gradient method. The gradient $\nabla q$ is provided by the score of the pretrained diffusion model, and $\nabla_\theta f$ is computed via automatic differentiation.

**6.1. Parameterized Blur Kernels: Diffraction-Limited Systems.** We first validate this strategy on diffraction-limited systems, where the PSF $h_\theta$ is defined by the squared magnitude of the Fourier transform of the pupil function $|\mathcal{F}(\exp(-i2\pi\phi_\theta))|^2$. The phase mask $\phi_\theta$ is parameterized by the first $P = 5$ Zernike polynomials (defocus, astigmatism, and spherical aberrations) [24].

In this family, “larger” kernels correspond to higher-order aberrations (larger coefficients
in $\theta$). As shown in Figure 11, initializing with a complex, large-aberration kernel successfully navigates the landscape to recover the ground truth. Conversely, starting with a “simpler” kernel (closer to the Airy pattern) leads to an immediate collapse into the no-blur solution, confirming the global pull of the blurry trap.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Example of joint posterior minimization for a diffraction-limited blur kernel. The Algorithm 6.1 is initialized at $y$ (observed with noise level $\sigma = 0.01$) and the kernels $h_{\text{init.1}}$ or $h_{\text{init.2}}$. The results for $h_{\text{init.1}}$ in the 3rd row is satisfactory, while the result for $h_{\text{init.2}}$ in the 4-th row is close to the no-blur solution. This illustrates the importance of initialization when minimizing the posterior jointly. Additional examples are given in Figure 15.</figcaption></figure>

**6.2. General Blur Kernels in the Simplex.** We further extend the evaluation to general kernels living on the unit simplex, where the mapping $\theta \mapsto h_{\theta}$ is the orthogonal projection onto the simplex. We contrast a “small” initialization (a narrow Gaussian) with a “large” initialization (a uniform kernel).

**Synthetic validation.** Results on four distinct kernels in Figure 12 confirm our theoretical predictions: the small initialization is captured by the Dirac spike, while the uniform initialization, representing the maximal-entropy state furthest from $\theta = 0$, successfully enters the recovery basin. Further results are provided in Figure 15 for diffraction-limited kernels and in Figure 14 for simplex kernels, demonstrating the consistent efficacy of the proposed strategy
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Illustration of the Algorithm 6.1 and the influence of the kernel initialization $h(\hat{\theta}_0)$. A “large” initialization (uniform kernel, top two rows) of the kernel results in a satisfactory recovery of the image and the kernel. A “small” initialization (a Gaussian kernel with standard-deviation of 1.5, bottom two rows) of the kernel results in a convergence to the no-blur solution. We use the FFHQ-256 model from [36].</figcaption></figure>

across various blur types and image models.

**Real-world data: the Köhler dataset.** To demonstrate robustness beyond exact training distributions, we evaluate the strategy on the Köhler dataset [16] using a pretrained DRUNet denoiser [41] as a prior. Following the Plug-and-Play (PnP) framework [23], the denoiser can be related to the gradient of an implicit prior. This is justified by the score-matching perspective and Tweedie's formula, where the denoiser output approximates a step along the manifold's gradient $\nabla q$.

By initializing with a 51 × 51 uniform kernel and applying Algorithm 6.1, we recover high-fidelity images and accurate motion trajectories (Figure 13). This success on real data, achieved with a generic denoiser, highlights the universal nature of the Morse-Bott recovery basin and the practical efficacy of topologically informed initialization.

**7. Conclusion.** Our most significant positive result is the general recovery guarantee provided in Theorem 3.3. We establish that for a broad class of blind inverse problems
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Results on the Köhler dataset [16]. Top: real blurry images and true kernels. Middle: deblurred images and estimated kernels. Bottom: estimated kernel along iterations.</figcaption></figure>

$y = A(\theta)x + b$, local recovery is not only possible but stable, provided the operator and the manifold satisfy three geometric conditions: identifiability on the manifold, parameter identifiability, and manifold-operator decoupling. This framework offers a rigorous template for analyzing a wide variety of blind problems beyond deconvolution, such as blind super-resolution or blind calibration in medical imaging.

However, we also demonstrate that this local stability is often overshadowed by a global spectral bias. In the case of blind deconvolution, we demonstrate that the “blurry trap”—a phenomenon well-documented for sparse priors [20]—persists even with state-of-the-art diffusion models. We show that the learned prior’s inherent preference for low-frequency content
drives the global MAP estimator toward degenerate solutions, proving that higher generative quality does not resolve the geometric degeneracy rooted in the forward operator. This suggests that for any family of operators characterized by varying Fourier decay, the global MAP estimator remains fundamentally unreliable.

Ultimately, our findings shift the focus of blind inverse problem research from the search for global maximizers to the characterization of local recovery basins. We conclude that successful recovery in high-dimensional imaging is a problem of *topologically informed optimization*: one must design initialization strategies that place the system within the recovery basin, thereby bypassing the global traps inherent to the statistics of natural images.

**Acknowledgments.** The authors acknowledge a support from the ANR Micro-Blind (ANR-21-CE48-0008) and from the ANR CLEAR-Microscopy (ANR-25-CE45-3780). This work was performed using HPC resources from GENCI-IDRIS (Grant AD011012210).
## Appendix A. Details on the experimental setup.

*Diffusion models.* In our experiments, we focus on state-of-the-art priors based on diffusion models. In what follows, bold fonts are used to denote random vectors. Diffusion models are powerful generative models and based on the following stochastic differential equation (SDE) on $\mathbb{R}^d$ [36]:

$$ (A.1) \qquad d\mathbf{x}_t = f(\mathbf{x}_t, t) dt + g(t) dw_t, $$

where $w_t$ is the standard Brownian process, $f: \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}^d$ denotes a drift term and $g: \mathbb{R} \to \mathbb{R}$ denotes a diffusion coefficient. The choice of $f$ and $g$ define the diffusion process, with various qualitative properties [36, 14]. In this model, $\mathbf{x}_0$ follows the clean image distribution $p_{\mathbf{x}}$. For long time horizon $T$ and under technical conditions, the solution of the SDE $\mathbf{x}_T$ approximately follows a normal distribution with no information left on $p_{\mathbf{x}}$, whose parameters depend on $f$ and $g$. This process can be reversed [2] using the following SDE running backward in time, which serves as a generative process:

$$ (A.2) \qquad d\mathbf{x}_t = [f(\mathbf{x}_t, t) - g(t)^2 \nabla \log p_t(\mathbf{x}_t)] dt + g(t) dw_t. $$

Here $p_t$ denotes the density at time $t$ and the solutions to this reversed-time SDE are distributed according to $p_{\mathbf{x}}$ at time $t = 0$. To generate new samples, we start from $\mathbf{x}_T \sim \mathcal{N}(0, \text{Id})$ and solve the backward SDE (A.2) for $t \in [T, 0]$ by a discretization scheme, e.g., Euler-Maruyama. Training a diffusion model amounts to estimating the Stein's score function of each marginal distribution, $\nabla \log p_t$, by score matching [40]. The score $\nabla \log p_t$ is typically parameterized by a time-dependent neural network $s_{\theta}(\cdot, t)$. For standard diffusion processes and any $t > 0$, $p_t$ corresponds to a Gaussian mixture, whose variance increases with $t$. In particular, for any $t > 0$, $p_t$ is smooth and strictly positive. In contrast, the data distribution $p_0$ maybe singular (e.g. supported on a low-dimensional manifold), making $\log p_0$ and its gradient ill-defined.

In what follows, we let $q_t \stackrel{\text{def}}{=} -\log p_t$ denote the potential (negative-log-prior) term. We evaluate potentials with a small value of time $t_{\epsilon} = 10^{-3}$. This value was chosen to obtain a sufficiently smooth approximation of $p_0$ compatible with gradient and Hessian computation, while still providing a good approximation of the prior $p_0$. In what follows, we simply write $q$ to denote $-\log p_{t_{\epsilon}}$, and we subsequently use $\nabla q(\cdot)$ to describe the action of the neural network $-s_{\theta}(\cdot, t_{\epsilon})$. This quantity is proportional to a difference between a noisy image and its denoised version via Tweedie's formula [40].

*Likelihood computation.* There exists a deterministic process with trajectories sharing the same marginal distribution as the SDE in (A.1). It can be obtained by the Fokker-Plank equation and it is defined by the following ordinary differential equation (ODE) [36, 14]:

$$ (A.3) \qquad d\mathbf{x}_t = v(\mathbf{x}_t, t) dt, \quad \text{where } v(\mathbf{x}_t, t) \stackrel{\text{def}}{=} f(\mathbf{x}_t, t) + \frac{1}{2} g(t)^2 \nabla q_t(\mathbf{x}_t). $$

If $\mathbf{x}_0 \sim p_{\mathbf{x}}$, then the trajectory $\mathbf{x}_t$ satisfies $\mathbf{x}_t \sim p_t(\mathbf{x}_t)$ for all $t \in [0, T]$. We can compute the potential of an image $x_0$ with the instantaneous change of variables formula [7]:

$$ (A.4) \qquad -\log p_{\mathbf{x}}(x_0) = -\log p_T(x_T) - \int_0^T \nabla \cdot v(x_t, t) dt, $$
where $x_t$ is the trajectory of the ODE in (A.3). The divergence term $\nabla \cdot v(x_t, t)$ can be approximated by the Skilling-Hutchinson trace estimator [35]:

$$ \nabla \cdot v(x_t, t) = \mathbb{E}_{\mathbf{u}} [\mathbf{u}^T J_v(x_t, t) \mathbf{u}] $$

where $J_v$ is the Jacobian of $v(\cdot, t)$, and $\mathbf{u}$ is a random variable with zero mean and identity covariance. In our implementation, we use a Rademacher random variable. The vector-Jacobian product $\mathbf{u}^T J_v(x_t, t)$ is computed using reverse-mode automatic differentiation with roughly the same cost as evaluating $v(x_t, t)$. The Skilling-Hutchinson trace estimator is unbiased so by averaging multiple runs, we can approximate the divergence. The integral can be evaluated using numerical ODE solver. Finally, at time $T$, the random variable $\mathbf{x}_T$ follows a Gaussian distribution with zero mean and standard deviation $\sigma(T)$ (depending on the model and on the dataset), which allows us to compute the term $\log p_T(x_T)$ in (A.4).

In all of our experiments, we use the Runge-Kutta 45 numerical scheme, provided by `scipy.integrate.solve_ivp` with tolerance value `atol=1e-5` and `rtol=1e-5`. We use a similar implementation to [36], but we compute the exact value of the potential instead of the bit-ber-dim value in the original implementation. This process is numerically intensive, as it requires hundreds of evaluations of the neural network $\nabla q_t$.

*Datasets, pretrained models and implementation details.* We use a model<fnref n="1" /> trained on the 256 × 256 FFHQ dataset from Song et al. [36] and three models<fnref n="2" /> trained on the ImageNet, FFHQ and AFHQ datasets at a 64 × 64 resolution from the widely used EDM framework by Karras et al. [14]. The latter are more recent and produce higher quality images. For brevity, we refer to these models throughout the paper as FFHQ-256, ImageNet-64, FFHQ-64, and AFHQ-64 respectively. ImageNet-64 is a class-conditional model, throughout the paper, we use the first 10 classes of the ImageNet dataset. We summarize the pre-trained models used in this paper in Table 2.

<caption>Table 2</caption>

*Pre-trained diffusion models used in our study.* The computation time was evaluated using a single NVIDIA V100 SMX2-16GB GPU and with 10 CPU cores from an Intel Cascade Lake 6248 processor.

<table><thead><tr><th>Pre-trained model</th><th>FFHQ-256 [36]</th><th>ImageNet-64 [14]</th><th>FFHQ-64 [14]</th><th>AFHQ-64 [14]</th></tr></thead><tbody><tr><td>Train dataset</td><td>FFHQ</td><td>ImageNet</td><td>FFHQ</td><td>AFHQ</td></tr><tr><td>SDE/ODE formulation</td><td>VE SDE</td><td>EDM</td><td>EDM</td><td>EDM</td></tr><tr><td>Resolution</td><td>256 × 256 × 3</td><td>64 × 64 × 3</td><td>64 × 64 × 3</td><td>64 × 64 × 3</td></tr><tr><td>Num. parameters (million)</td><td>66</td><td>296</td><td>63</td><td>63</td></tr><tr><td>Model eval. (sec./image)</td><td>0.036</td><td>0.014</td><td>0.006</td><td>0.006</td></tr><tr><td>Potential eval. (sec/image)</td><td>8.461</td><td>10.26</td><td>5.05</td><td>5.05</td></tr></tbody></table>

641

<footnote n="1">¹Pre-trained model: https://github.com/yang-song/score_sde_pytorch</footnote>

<footnote n="2">²Pre-trained model: https://github.com/NVlabs/edm</footnote>
642

643

**Appendix B. Additional numerical result.** We provide additional numerical results to complement those in subsection 6.2.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** More numerical results with a simplex constraint using different pairs $(\bar{x}, \bar{h})$, similarly as Figure 12. The algorithm *Algorithm 6.1* is initialized by a uniform kernel. The noise level is $\sigma = 0.01$. The kernel's shape is recovered approximately and the estimated images look sharper.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** More examples on the influence of initialization for blind deconvolution, with diffraction-limited blur kernels as in Figure 11. The initialization are $h_{\text{init.1}} = h(\bar{\theta}+0.15 \cdot \text{sign}(\bar{\theta}))$ and $h_{\text{init.2}} = h(\bar{\theta}-0.075 \cdot \text{sign}(\bar{\theta}))$ and $x_0 = y$. Again, using large kernels for initialization leads to a decent kernel reconstruction, while small kernels lead to the no-blur solution.</figcaption></figure>
**Appendix C. Proofs.** We start by the following lemma which characterizes the positiveness of a block matrix.

**Lemma C.1.** Let $A \in \mathbb{R}^{M \times N}$, $J \in \mathbb{R}^{M \times P}$ and $D \in \mathbb{R}^{N \times N}$ with $D$ symmetric, positive semi-definite. Consider the block matrix:

$$H = \begin{pmatrix} A^T A + D & A^T J \\ J^T A & J^T J \end{pmatrix}.$$

Then $H$ is positive definite if and only if:

1. $\ker J = \{0\}$,

2. $A \ker D \cap \text{Im } J = \{0\}$,

3. $\ker D \cap \ker A = \{0\}$.

*Proof.* Since $D$ is symmetric, positive semi-definite, let $\sqrt{D}$ be its symmetric semi-definite square root. Define the matrix $L$ as:

$$L = \begin{pmatrix} A & J \\ \sqrt{D} & 0 \end{pmatrix}.$$

We observe that $H = L^T L$. Thus, $H$ is positive semi-definite, and it is positive definite if and only if $\ker L = \{0\}$.

A vector $(u, v) \in \mathbb{R}^N \times \mathbb{R}^P$ is in $\ker L$ if and only if:

$$(C.1) \qquad \begin{cases} Au + Jv = 0 \\ \sqrt{D}u = 0 \end{cases} \iff \begin{cases} Au + Jv = 0 \\ u \in \ker D \end{cases}$$

*Sufficiency:* Assume the three conditions hold. Let $(u, v) \in \ker L$. From (C.1), $u \in \ker D$ and $Au = -Jv$. Thus $Au \in A \ker D \cap \text{Im } J$. By condition (2), this implies $Au = 0$. Since $u \in \ker D$ and $Au = 0$, condition (3) implies $u = 0$. The system then reduces to $Jv = 0$, which by condition (1) implies $v = 0$. Thus $\ker L = \{0\}$.

*Necessity:* Assume that $H$ is positive definite, so that $\ker L = \{0\}$.

* If $\ker J \neq \{0\}$, there exists $v \neq 0$ such that $Jv = 0$. Then $(0, v) \in \ker L$.

* If $y \in A \ker D \cap \text{Im } J$ is non-zero, there exist $u \in \ker D$ and $v$ such that $Au = y = -Jv$. Then $Au + Jv = 0$ and $u \in \ker D$, so $(u, v) \in \ker L$. Note that $u \neq 0$ because $Au = y \neq 0$.

* If $u \in \ker D \cap \ker A$ with $u \neq 0$, then $(u, 0) \in \ker L$.

In all cases, violation of a condition implies $\ker L \neq \{0\}$, which is a contradiction. ■

**C.1. Proof of Theorem 3.2.** We now prove the stability of the non-blind MAP estimator, as stated in Theorem 3.2.

*Proof.* Define

$$F(x, y) \stackrel{\text{def}}{=} \nabla_x L(x; y) = \frac{1}{\sigma^2} A^T (Ax - y) + \nabla q(x).$$

Since $\bar{y} = A\bar{x}$ and $\bar{x}$ lies on the critical submanifold $M$ of $q$, we have

$$\nabla q(\bar{x}) = 0,$$
hence

$$F(\bar{x}, \bar{y}) = \frac{1}{\sigma^2} A^T (A\bar{x} - \bar{y}) + \nabla q(\bar{x}) = 0.$$

The differential of $F$ with respect to $x$ at $(\bar{x}, \bar{y})$ is

$$\bar{H} \stackrel{\text{def}}{=} \partial_x F(\bar{x}, \bar{y}) = \frac{1}{\sigma^2} A^T A + \nabla^2 q(\bar{x}).$$

Since $q$ is locally Morse-Bott at $\bar{x}$ with critical submanifold $\mathcal{M}$, we have

$$\ker \nabla^2 q(\bar{x}) = T_{\bar{x}} \mathcal{M}.$$

Moreover, since $\nabla^2 q(\bar{x}) \succeq 0$, applying Lemma C.1 with $D = \nabla^2 q(\bar{x})$ and $J = 0$, we deduce that $\bar{H}$ is positive definite:

$$\mu \stackrel{\text{def}}{=} \lambda_{\min}(\bar{H}) > 0,$$

where $\mu$ as the smallest eigenvalue of $\bar{H}$. Since $\bar{H}$ is invertible, the implicit function theorem yields neighborhoods $U$ of $\bar{x}$ and $V$ of $\bar{y}$, and a $C^1$ map $y \mapsto x(y)$ such that

$$F(x(y), y) = 0, \quad x(\bar{y}) = \bar{x}, \quad \nabla_x^2 L(x(y), y) \succ 0.$$

The last inequality is a consequence of the continuity of $\nabla_x^2 L$ and the fact that $\nabla_x^2 L(\bar{x}, \bar{y}) = \bar{H} \succ 0$.

The first equality reads $\frac{1}{\sigma^2} A^T (Ax(y) - y) + \nabla q(x(y)) = 0$. Differentiating w.r.t. $y$ yields

$$\left( \frac{1}{\sigma^2} A^T A + \nabla^2 q(x(y)) \right) \mathrm{Jac}_x(y) - \frac{1}{\sigma^2} A^T = 0, \quad \forall y \in V.$$

Taking $y = \bar{y}$, we have

$$\bar{H} \mathrm{Jac}_x(\bar{y}) - \frac{1}{\sigma^2} A^T = 0, \quad \text{hence} \quad \mathrm{Jac}_x(\bar{y}) = \frac{1}{\sigma^2} \bar{H}^{-1} A^T.$$

Therefore, first-order approximation yields

$$x(y) = x(\bar{y}) + \mathrm{Jac}_x(\bar{y})(y - \bar{y}) + o(\|y - \bar{y}\|)$$

or equivalently

$$x(y) - \bar{x} = \frac{1}{\sigma^2} \bar{H}^{-1} A^T (y - \bar{y}) + o(\|y - \bar{y}\|).$$

Since $\|\bar{H}^{-1}\| = 1/\lambda_{\min}(\bar{H}) = 1/\mu$, we obtain the stability bound as

$$\|x(y) - \bar{x}\| \le \frac{\|A\|}{\mu\sigma^2} \|y - \bar{y}\| + o(\|y - \bar{y}\|).$$

### C.2. Proof of Theorem 3.3 and Corollary 3.4.
We provide a joint proof of Theorem 3.3 and Corollary 3.4, as they are closely related.

*Proof.* The proof proceeds as follows
704 First-order optimality. Define the gradient of L as

$$
705 \qquad F(x, \theta, y) \stackrel{\text{def}}{=} \nabla_{x,\theta} L(x, \theta; y).
$$

706 A direct computation gives

$$
707 \qquad F(x, \theta, y) = \left( \frac{1}{\sigma^2} A(\theta)^T (A(\theta)x - y) + \nabla q(x) \right),
$$

708 where $J_x(\theta) = \frac{\partial}{\partial \theta}(A(\theta)x) \in \mathbb{R}^{M \times P}$ denotes the Jacobian of the measurements w.r.t the
709 parameters $\theta$. Since $A(\bar{\theta})\bar{x} = \bar{y}$ and $\nabla q(\bar{x}) = 0$, we obtain

$$
710 \qquad F(\bar{x}, \bar{\theta}, \bar{y}) = 0.
$$

711 Expression of the joint Hessian. Differentiating F with respect to (x, θ) at (x̄, θ̄, ŷ) gives the
712 joint Hessian

$$
713 \qquad \bar{H} \stackrel{\text{def}}{=} \nabla_{x,\theta}^2 L(\bar{x}, \bar{\theta}; \bar{y}) = \begin{pmatrix} \bar{H}_{xx} & \bar{H}_{x\theta} \\ \bar{H}_{\theta x} & \bar{H}_{\theta\theta} \end{pmatrix} = \frac{1}{\sigma^2} \begin{pmatrix} A(\bar{\theta})^T A(\bar{\theta}) + \sigma^2 \nabla^2 q(\bar{x}) & A(\bar{\theta})^T J_{\bar{x}}(\bar{\theta}) \\ J_{\bar{x}}(\bar{\theta})^T A(\bar{\theta}) & J_{\bar{x}}(\bar{\theta})^T J_{\bar{x}}(\bar{\theta}) \end{pmatrix}.
$$

714 Under the three geometric compatibility conditions and $\nabla^2 q(\bar{x}) \succeq 0$, Lemma C.1 implies
715 that $\bar{H}$ is positive definite:

$$
716 \qquad \bar{H} \succeq \mu I_{N+P}, \quad \mu \stackrel{\text{def}}{=} \lambda_{\min}(\bar{H}) > 0 \text{ is the smallest eigenvalue.}
$$

Application of the implicit function theorem. Since F is C¹ and D<sub>(x,θ)</sub>F(x̄, θ̄, ȳ) = H̄ is
invertible, the implicit function theorem yields neighborhoods U of (x̄, θ̄) and V of ȳ and a
unique C¹ map y ↦ (x(y), θ(y)) such that

$$
717 \qquad F(x(y), \theta(y), y) = 0, \quad (x(\bar{y}), \theta(\bar{y})) = (\bar{x}, \bar{\theta}), \quad \nabla^2_{x,\theta} L(x(y), \theta(y); y) > 0.
$$

717 The last inequality is a consequence of the continuity of $\nabla^2_{x,\theta} L$ and $\nabla^2_{x,\theta} L(\bar{x}, \bar{\theta}; \bar{y}) > 0$.

718 First-order expansion. The first equality reads

$$
719 \qquad (C.2) \qquad \frac{1}{\sigma^2} \left( A(\theta(y))^T (A(\theta(y))x(y) - y) + \sigma^2 \nabla q(x(y)) \right) = 0.
$$

720 Differentiating the above equality w.r.t $y$ then evaluating at $(\bar{x}, \bar{\theta})$, we have

$$
721 \quad (C.3) \quad \begin{pmatrix} \bar{H}_{xx} & \bar{H}_{x\theta} \\ \bar{H}_{\theta x} & \bar{H}_{\theta\theta} \end{pmatrix} \begin{pmatrix} \mathrm{Jac}_x(\bar{y}) \\ \mathrm{Jac}_\theta(\bar{y}) \end{pmatrix} = \frac{1}{\sigma^2} \begin{pmatrix} A(\bar{\theta})^T \\ J_{\bar{x}}(\bar{\theta})^T \end{pmatrix},
$$

722 where Jac<sub>x</sub> and Jac<sub>θ</sub> denote the Jacobian matrices of x(y) and θ(y) respectively. Therefore,

$$
723 \quad (C.4) \qquad \begin{pmatrix} \mathrm{Jac}_x(\bar{y}) \\ \mathrm{Jac}_\theta(\bar{y}) \end{pmatrix} = \frac{1}{\sigma^2} \bar{H}^{-1} \begin{pmatrix} A(\bar{\theta})^T \\ J_{\bar{x}}(\bar{\theta})^T \end{pmatrix}.
$$
724 Using the Schur complement formula for the inverse of a block matrix, we have

$$725 \quad (C.5) \qquad \bar{H}^{-1} = \begin{pmatrix} S_x^{-1} & -S_x^{-1} \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} \\ -\bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} & \bar{H}_{\theta\theta}^{-1} + \bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} \end{pmatrix},$$

726 where $S_x = \bar{H}_{xx} - \bar{H}_{x\theta}\bar{H}_{\theta\theta}^{-1}\bar{H}_{\theta x}$ is the Schur complement of $\bar{H}$ with respect to $\bar{H}_{\theta\theta}$. We deduce that

$$728 \quad (C.6) \qquad \mathrm{Jac}_x(\bar{y}) = \frac{1}{\sigma^2} S_x^{-1} (A(\bar{\theta})^T - \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T),$$

$$729 \quad \mathrm{Jac}_{\theta}(\bar{y}) = -\bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} \frac{1}{\sigma^2} A(\bar{\theta})^T + (\bar{H}_{\theta\theta}^{-1} + \bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1}) \frac{1}{\sigma^2} J_{\bar{x}}(\bar{\theta})^T$$

$$730 \quad (C.7) \qquad = -\frac{1}{\sigma^2} \bar{H}_{\theta\theta}^{-1} \bar{H}_{\theta x} S_x^{-1} (A(\bar{\theta})^T - \bar{H}_{x\theta} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T) + \frac{1}{\sigma^2} \bar{H}_{\theta\theta}^{-1} J_{\bar{x}}(\bar{\theta})^T \quad \blacksquare$$

### C.3. Proof of Lemma 4.1.

731 Proof. Let $(\hat{x}, \hat{h})$ be a minimizer of (4.1). If $\hat{h} \neq \delta$, we can replace $(\hat{x}, \hat{h})$ with $(\hat{h} \star \hat{x}, \delta)$ and 732 get a lower cost function. Indeed, the point $(\hat{h} \star \hat{x}, \delta)$ does not modify the value of the data 733 fidelity term in (4.1) and since we assumed that $q(\hat{x}) \ge q(\hat{h} \star \hat{x})$, $(\hat{h} \star \hat{x}_0, \delta)$ is also a solution. 734 This shows that the Dirac kernel is a solution to the joint MAP problem. If the inequality 735 were strict in (4.2), the solution would be unique and satisfy $\hat{h}_{\mathrm{MAP}} = \delta$. Note that when $\hat{h}_{\mathrm{MAP}}$ 736 is the Dirac mass, the estimated image coincides with a MAP-denoised image:

$$738 \quad (C.8) \qquad \hat{x}_{\mathrm{MAP}}(y) = \arg \min_x \frac{1}{2\sigma^2} \|x - y\|_2^2 + q(x). \quad \blacksquare$$

### REFERENCES

740 [1] A. AHMED, B. RECHT, AND J. ROMBERG, *Blind deconvolution using convex programming*, IEEE Transactions on Information Theory, 60 (2013), pp. 1711–1732.

741 [2] B. D. ANDERSON, *Reverse-time diffusion equation models*, Stochastic Processes and their Applications, 12 (1982), pp. 313–326, https://doi.org/10.1016/0304-4149(82)90051-5, https://www.sciencedirect.com/science/article/pii/0304414982900515.

742 [3] A. BENICHOUX, E. VINCENT, AND R. GRIBONVAL, *A fundamental pitfall in blind deconvolution with sparse and shift-invariant priors*, in 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE, 2013, pp. 6108–6112.

743 [4] A. BORA, A. JALAL, E. PRICE, AND A. G. DIMAKIS, *Compressed sensing using generative models*, in International conference on machine learning, PMLR, 2017, pp. 537–546.

744 [5] E. J. CANDES AND T. TAO, *Decoding by linear programming*, IEEE transactions on information theory, 51 (2005), pp. 4203–4215.

745 [6] T. F. CHAN AND C.-K. WONG, *Total variation blind deconvolution*, IEEE transactions on Image Processing, 7 (1998), pp. 370–375.

746 [7] R. T. CHEN, Y. RUBANOVA, J. BETTENCOURT, AND D. K. DUVENAUD, *Neural ordinary differential equations*, Advances in neural information processing systems, 31 (2018).

747 [8] H. CHUNG, J. KIM, S. KIM, AND J. C. YE, *Parallel diffusion models of operator and image for blind inverse problems*, in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2023, pp. 6059–6069.

748 [9] V. DEBARNOT AND P. WEISS, *Blind inverse problems with isolated spikes*, Information and Inference: A Journal of the IMA, 12 (2023), pp. 26–71.
761

762

[10] V. DEBARNOT AND P. WEISS, *Deep-blur: Blind identification and deblurring with convolutional neural networks*, Biological Imaging, (2024).

763

764

[11] B. EFRON, *Tweedie's formula and selection bias*, Journal of the American Statistical Association, 106 (2011), pp. 1602-1614.

765

766

[12] S. HURault, A. LECLAIRE, AND N. PAPADAKIS, *Gradient step denoiser for convergent plug-and-play*, in International Conference on Learning Representations, 2021.

767

768

769

[13] R. KARCZEWSKI, M. HEINONEN, AND V. GARG, *Diffusion models as cartoonists: The curious case of high density regions*, in The Thirteenth International Conference on Learning Representations, 2025, https://openreview.net/forum?id=RiS2cxpENN.

770

771

[14] T. KARRAS, M. AITTALA, T. AILA, AND S. LAINE, *Elucidating the design space of diffusion-based generative models*, in Proc. NeurIPS, 2022.

772

773

774

[15] M. KECH AND F. KRAHMER, *Optimal injectivity conditions for bilinear inverse problems with applications to identifiability of deconvolution problems*, SIAM Journal on Applied Algebra and Geometry, 1 (2017), pp. 20-37.

775

776

777

778

[16] R. KÖHLER, M. HIRSCH, B. MOHLER, B. SCHÖLKOPF, AND S. HARMELING, *Recording and playback of camera shake: Benchmarking blind deconvolution with a real-world database*, in Computer Vision-ECCV 2012: 12th European Conference on Computer Vision, Florence, Italy, October 7-13, 2012, Proceedings, Part VII 12, Springer, 2012, pp. 27-40.

779

780

[17] F. KRAHMER AND D. STÖGER, *On the convex geometry of blind deconvolution and matrix completion*, Communications on Pure and Applied Mathematics, 74 (2021), pp. 790-832.

781

782

[18] J. D. LEE, M. SIMCHOWITZ, M. I. JORDAN, AND B. RECHT, *Gradient descent only converges to minimizers*, in Conference on learning theory, PMLR, 2016, pp. 1246-1257.

783

784

[19] J. M. LEE, *Smooth manifolds*, in Introduction to smooth manifolds, Springer, 2003, pp. 1-29.

785

786

[20] A. LEVIN, Y. WEISS, F. DURAND, AND W. T. FREEMAN, *Understanding and evaluating blind deconvolution algorithms*, in 2009 IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 1964-1971, https://doi.org/10.1109/CVPR.2009.5206815.

787

788

[21] S. LING AND T. STROHMER, *Self-calibration and bilinear inverse problems via linear least squares*, SIAM Journal on Imaging Sciences, 11 (2018), pp. 252-292.

789

790

[22] Y. NESTEROV, *Introductory lectures on convex optimization: A basic course*, vol. 87, Springer Science & Business Media, 2013.

791

792

[23] M. H. NGUYEN AND P. WEISS, *Comparing Plug-and-Play and Unrolled networks*. working paper or preprint, Sept. 2024, https://hal.science/hal-04703008.

793

794

[24] R. J. NOLL, *Zernike polynomials and atmospheric turbulence*, JOSA, 66 (1976), pp. 207-211.

795

796

[25] J. P. OLIVEIRA, M. A. FIGUEIREDO, AND J. M. BIOCAS-DIAS, *Parametric blur estimation for blind restoration of natural images: Linear motion and out-of-focus*, IEEE Transactions on Image Processing, 23 (2013), pp. 466-477.

797

798

799

[26] J. PAN, Z. HU, Z. SU, AND M.-H. YANG, *Deblurring text images via l0-regularized intensity and gradient prior*, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2014, pp. 2901-2908.

800

801

[27] J. PAN, D. SUN, H. PFISTER, AND M.-H. YANG, *Blind image deblurring using dark channel prior*, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 1628-1636.

802

803

[28] J. PIDSTRIGACH, *Score-based generative models detect manifolds*, Advances in Neural Information Processing Systems, 35 (2022), pp. 35852-35865.

804

805

[29] P. POPE, C. ZHU, A. ABDELKADER, M. GOLDBLUM, AND T. GOLDSTEIN, *The intrinsic dimension of images and its impact on learning*, in International Conference on Learning Representations.

806

807

[30] N. PUSTELNIK, A. BENAZZA-BENHAYIA, Y. ZHENG, AND J.-C. PESQUET, *Wavelet-based image deconvolution and reconstruction*, Wiley encyclopedia of electrical and electronics engineering, (2016).

808

809

[31] Q. REBJOCK AND N. BOUMAL, *Fast convergence to non-isolated minima: four equivalent conditions for c 2 functions*, Mathematical Programming, (2024), pp. 1-49.

810

811

[32] E. T. REEHORST AND P. SCHNITER, *Regularization by denoising: Clarifications and new interpretations*, IEEE transactions on computational imaging, 5 (2018), pp. 52-67.

812

813

814

[33] R. ROMBACH, A. BLATTMANN, D. LORENZ, P. ESSER, AND B. OMMER, *High-resolution image synthesis with latent diffusion models*, in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2022, pp. 10684-10695.
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839

[34] A. SHAJKOFCI AND M. LIEBLING, *Spatially-Variant CNN-Based Point Spread Function Estimation for Blind Deconvolution and Depth Estimation in Optical Microscopy*, IEEE Transactions on Image Processing, 29 (2020), pp. 5848-5861.

[35] J. SKILLING, *The Eigenvalues of Mega-dimensional Matrices*, Springer Netherlands, Dordrecht, 1989, pp. 455-466, https://doi.org/10.1007/978-94-015-7860-8_48, https://doi.org/10.1007/978-94-015-7860-8_48.

[36] Y. SONG, J. SOHL-DICKSTEIN, D. P. KINGMA, A. KUMAR, S. ERMON, AND B. POOLE, *Score-based generative modeling through stochastic differential equations*, in International Conference on Learning Representations, 2020.

[37] J. P. STANCZUK, G. BATZOLIS, T. DEVENEY, AND C.-B. SCHÖNLIEB, *Diffusion models encode the intrinsic dimension of data manifolds*, in Forty-first International Conference on Machine Learning, 2024.

[38] J. TACHELLA, M. TERRIS, S. HURault, A. WANG, L. DAVY, J. SCANVIC, V. SECHAUD, R. VO, T. MOREAU, T. DAVIES, D. CHEN, N. LAURENT, B. MONROY, J. DONG, Z. HU, M.-H. NGUYEN, F. SARRON, P. WEISS, P. ESCANDE, M. MASSIAS, T. MODRZYK, B. LEVAC, T. I. LIAUDAT, M. SONG, J. HERTRICH, S. NEUMAYER, AND G. SCHRAMM, *Deepinverse: A python package for solving imaging inverse problems with deep learning*, Journal of Open Source Software, 10 (2025), p. 8923, https://doi.org/10.21105/joss.08923, http://dx.doi.org/10.21105/joss.08923.

[39] S. VAITER, G. PEYRÉ, AND J. FADILI, *Model consistency of partly smooth regularizers*, IEEE Transactions on Information Theory, 64 (2017), pp. 1725-1737.

[40] P. VINCENT, *A connection between score matching and denoising autoencoders*, Neural computation, 23 (2011), pp. 1661-1674.

[41] K. ZHANG, Y. LI, W. ZUO, L. ZHANG, L. VAN GOOL, AND R. TIMOFTE, *Plug-and-play image restoration with deep denoiser prior*, IEEE Transactions on Pattern Analysis and Machine Intelligence, 44 (2021), pp. 6360-6376.