FLOWER: A Flow-Matching Solver for Inverse Problems

## Abstract

Abstract We introduce Flower , a solver for inverse problems. It leverages a pre-trained flow model to produce reconstructions that are consistent with the observed measurements. Flower operates through an iterative procedure over three steps: (i) a flow-consistent destination estimation, where the velocity network predicts a denoised target; (ii) a refinement step that projects the estimated destination onto a feasible set defined by the forward operator; and (iii) a time-progression step that re-projects the refined destination along the flow trajectory. We provide a theoretical analysis that demonstrates how Flower approximates Bayesian posterior sampling, thereby unifying perspectives from plug-and-play methods and generative inverse solvers. On the practical side, Flower achieves state-of-the-art reconstruction quality while using nearly identical hyperparameters across various inverse problems.

## 1 Introduction

Inverse problems are central to computational imaging and computer vision (McCann & Unser, 2019; Zeng, 2001). Their goal is to reconstruct an underlying signal ${\bf{x}}\in\mathbb{R}^{d}$ from its observed measurements ${\bf{y}}\in\mathbb{R}^{M}$. Here, we focus on linear inverse problems, such that the acquisition of the measurements follows the model

$$ ${\bf{y}}={\bf{H}}{\bf{x}}+{\bf{n}}$ (1) $$

for some linear forward operator ${\bf{H}}\colon\mathbb{R}^{d}\to\mathbb{R}^{M}$ and additive white Gaussian noise ${\bf{n}}\sim\mathcal{N}({\bf{0}},\sigma_{n}^{2}{\bf{I}})$. From a Bayesian perspective, the simplest reconstruction approach is to obtain the maximum-likelihood estimation

$$ $\hat{{\bf{x}}}_{\mathrm{MLE}}=\operatorname*{arg\,max}_{{\bf{x}}\in\mathbb{R}^{d}}p_{{\bf{Y}}\mid{\bf{X}}={\bf{x}}}({\bf{y}})=\operatorname*{arg\,min}_{{\bf{x}}\in\mathbb{R}^{d}}\frac{1}{2\sigma_{n}^{2}}\left\rVert{\bf{H}}{\bf{x}}-{\bf{y}}\right\rVert_{2}^{2}.$ (2) $$

However, this problem is ill-posed and yields poor-quality solutions. Another approach is to obtain the maximum a posteriori estimation (MAP)

$$ $\hat{{\bf{x}}}_{\mathrm{MAP}}=\operatorname*{arg\,max}_{{\bf{x}}\in\mathbb{R}^{d}}p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}({\bf{x}})=\operatorname*{arg\,min}_{{\bf{x}}\in\mathbb{R}^{d}}\left(\frac{1}{2\sigma_{n}^{2}}\left\rVert{\bf{H}}{\bf{x}}-{\bf{y}}\right\rVert_{2}^{2}-\log p_{{\bf{X}}}({\bf{x}})\right),$ (3) $$

which requires the knowledge of the prior distribution $p_{{\bf{X}}}$ of images, a quantity that is generally unknown. The minimization problem of ([3](https://arxiv.org/html/2509.26287v1#S1.E3)) is consistent with the variational perspective of inverse problems, where the term $\left(-\log p_{{\bf{X}}}\left({\bf{x}}\right)\right)$ is replaced by a regularizer $\mathcal{R}({\bf{x}})$ that encodes some properties of the images. From classic signal processing to the advent of deep learning, the design of a good regularizer $\mathcal{R}$ has been of interest. Classic signal processing relies on the smoothness or sparsity of images to introduce wavelet- or total-variation-based regularizers (Rudin et al., 1992; Figueiredo & Nowak, 2003; Beck & Teboulle, 2009). Some methods build upon classical models and try to learn such criteria in a data-driven manner (Roth & Black, 2009; Goujon et al., 2024; Ducotterd et al., 2025; Pourya et al., 2025). Plug-and-play (PnP) approaches focus on the implicit replacement of $\mathcal{R}$ by its proximal operator, with a learned neural network that serves as a denoiser (Venkatakrishnan et al., 2013; Zhang et al., 2022; Hurault et al., 2022b; a). Although MAP estimations tend to have a good reconstruction quality, they do not necessarily provide the minimum-mean-square estimator $\hat{{\bf{x}}}_{\mathrm{MMSE}}$ that is best in terms of the peak signal-to-noise ratio (PSNR). To estimate $\hat{{\bf{x}}}_{\mathrm{MMSE}}$, one would have to compute the posterior mean $\hat{{\bf{x}}}_{\mathrm{MMSE}}=\mathbb{E}[{{\bf{X}}|{\bf{Y}}={\bf{y}}}]$. Moreover, for perceptual metrics, it is better to generate a sample from $p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}$ instead of an estimator of the distribution.

The objective of generative modeling is to sample from a target distribution $p_{{\bf{X}}}$. In practice, this distribution is unknown, and one typically only has access to a finite collection of its samples. Numerous approaches have been proposed to address this issue. Among them, diffusion models (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2020) and, more recently, flow-matching methods (Lipman et al., 2023) represent the state of the art in scalable generative modeling for images.

Flow matching, introduced by Lipman et al. (2023), constructs a continuous-time generative process by parameterizing the velocity field of an ordinary differential equation (ODE) as a neural network. It takes inspiration from optimal transport and continuous normalizing flows (Ambrosio et al., 2008; Hagemann et al., 2022) and transports an initial source distribution $p_{{\bf{X}}_{0}}$ to a target distribution $p_{{\bf{X}}_{1}}\approx p_{{\bf{X}}}$ . The choice of probability paths from $p_{{\bf{X}}_{0}}$ to $p_{{\bf{X}}_{1}}$ are numerous, with Gaussian paths recovering diffusion as a special case (Albergo & Vanden-Eijnden, 2023). However, flow matching mostly focuses on straight-line paths, which yields competitive performance and improved sampling efficiency (Liu et al., 2023; Liu, 2022).

The remarkable success of generative models in image generation motivates their extension to inverse problems, where the goal shifts from the sampling of the prior distribution $p_{\mathbf{X}1}$ to the sampling of the posterior $p_{\mathbf{X}_{1}\mid\mathbf{Y}=\mathbf{y}}$. Several inverse solvers based on diffusion models have been introduced (Chung et al., 2023; 2024; Kawar et al., 2022; Song et al., 2023; Zhu et al., 2023; Zhang et al., 2025; Mardani et al., 2024). Recent efforts also focus on flow-based solvers (Pokle et al., 2024; Martin et al., 2025). Existing approaches can be broadly grouped into two categories: (i) methods that approximate the posterior score (velocity field) with gradient corrections along the generative path; and (ii) PnP strategies that alternate between generative (diffusion or flow) updates and data-consistency steps. In this work, we introduce a novel solver based on flow matching that achieves state-of-the-art results for flow-based inverse problems. Our approach departs from existing methods by framing the problem through a Bayesian ancestral-sampling perspective, which gives rise to a simple three-step procedure with a natural plug-and-play interpretation. Our main contributions are as follows.

- 1.
Flow-matching solver for inverse problems. We introduce *Flower*, an inverse problem solver that consists of three steps: (1) a *flow-consistent destination estimation*, where the velocity network is used to predict a destination, interpretable as denoising; (2) a *measurement-aware refinement*, in which the estimated destination is projected onto the feasible set defined by the forward operator; and (3) a *time progression*, where the refined destination is re-projected along the flow path.
- 2.
Bayesian analysis and relation to PnP. We provide a Bayesian analysis in which we demonstrate how and under what considerations *Flower* generates valid posterior samples from $p_{{\bf{X}}_{1}|{\bf{Y}}={\bf{y}}}$. Specifically, we show that Step 1 computes the conditional expectation $\mathbb{E}[{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}]$, which we then use for the approximation of $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}$. Through this approximation, we show that Step 2 generates a sample $\tilde{{\bf{x}}}({\bf{x}}_{t},{\bf{y}})\sim\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$. Step 3 then updates the trajectory given the refined destination $\tilde{{\bf{x}}}({\bf{x}}_{t},{\bf{y}})$ and draws a sample ${\bf{x}}_{t+\Delta t}\sim\tilde{p}_{{\bf{X}}{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$, which by induction and ancestral sampling, produces the final sample ${\bf{x}}_{1}\sim\tilde{p}_{{\bf{X}}_{1}\mid{\bf{Y}}={\bf{y}}}$. These steps rely on three assumptions: the velocity network is optimally trained for unconditional flow matching; the acquisition of measurements follows the forward model in ([1](https://arxiv.org/html/2509.26287v1#S1.E1)); and the source and target distributions are independent. To the best of our knowledge, this Bayesian construction is novel within flow-based solvers for inverse problems. Although this construction is key to the derivation of our solver, the resulting procedure closely mirrors the PnP methods. Thus, *Flower* can be seen as a plug-and-play inverse solver with general application, with our Bayesian justification establishing a novel link between posterior sampling and the PnP frameworks in generative models.
- 3.
Numerical validation. We first examine a controlled setup with Gaussian mixture models and show that *Flower* successfully recovers posterior samples. We then evaluate our method on standard inverse problem benchmarks for flow matching. We achieve competitive performance, with nearly identical hyperparameters across all tasks.

The remainder of this paper is organized as follows. In Section [2](https://arxiv.org/html/2509.26287v1#S2), we review the fundamentals of flow matching along with the mathematical tools required for the development of our method. We then introduce *Flower* in Section [3](https://arxiv.org/html/2509.26287v1#S3) and present the associated theoretical analysis. In Section [4](https://arxiv.org/html/2509.26287v1#S4), we discuss related work and highlight their similarities and differences with our approach. Finally, we report our numerical results in Section [5](https://arxiv.org/html/2509.26287v1#S5).

## 2 Background

In this section, we first review flow matching and then proximal operators, which are required to understand our method.

### 2.1 Flow Matching

Let $p_{{\bf{X}}_{0}}$ be a source distribution that is easy to sample and let $p_{{\bf{X}}_{1}}$ be a target distribution that we want to sample from. A time-dependent flow $\psi_{t}$ transports $p_{{\bf{X}}_{0}}$ to $p_{{\bf{X}}_{1}}$ via the ODE

$$ $\frac{\mathrm{d}\psi_{t}({\bf{x}})}{\mathrm{d}t}\;=\;{\bf{v}}_{t}\!\big(\psi_{t}({\bf{x}})\big),\qquad t\in[0,1],$ (4) $$

for some velocity field ${\bf{v}}_{t}:\mathbb{R}^{d}\to\mathbb{R}^{d}$. The intermediate variables ${\bf{X}}_{t}=\psi_{t}({\bf{X}}_{0})$ follow a distribution $p_{{\bf{X}}_{t}}$. The objective of flow matching is to approximate ${\bf{v}}_{t}$ with a neural network ${\bf{v}}^{\theta}_{t}$, which will allow us to sample from $p_{{\bf{X}}_{1}}$. However, the determination of the flow-matching loss

$$ $\mathcal{L}_{\mathrm{FM}}(\theta)=\mathbb{E}_{t\sim\mathcal{U}[0,1]}\,\mathbb{E}_{{\bf{x}}_{t}\sim p_{{\bf{X}}_{t}}}\left[\left\|{\bf{v}}^{\theta}_{t}({\bf{x}}_{t})-{\bf{v}}_{t}({\bf{x}}_{t})\right\|_{2}^{2}\right]$ (5) $$

is challenging, as it requires access to the marginal velocity field ${\bf{v}}_{t}({\bf{x}}_{t})$. To address this, we focus on the conditional velocity ${\bf{v}}_{t}({\bf{x}}_{t}\mid{\bf{x}}_{1})$ and define the conditional straight-line flow and velocity

$$ ${\bf{x}}_{t}=\psi_{t}({\bf{x}}_{0}\mid{\bf{x}}_{1})=(1-t)\,{\bf{x}}_{0}+t\,{\bf{x}}_{1},\quad{\bf{v}}_{t}({\bf{x}}_{t}\mid{\bf{x}}_{1})={\bf{x}}_{1}-{\bf{x}}_{0}.$ (6) $$

This leads to the practical conditional flow-matching loss

$$ $\mathcal{L}_{\mathrm{CFM}}(\theta)=\mathbb{E}_{t\sim\mathcal{U}[0,1]}\,\mathbb{E}_{({\bf{x}}_{0},{\bf{x}}_{1})\sim\pi}\left[\left\|{\bf{v}}^{\theta}_{t}\big((1-t)\,{\bf{x}}_{0}+t\,{\bf{x}}_{1},t\big)-({\bf{x}}_{1}-{\bf{x}}_{0})\right\|_{2}^{2}\right],$ (7) $$

where $\pi\in\Pi(p_{{\bf{X}}_{0}},p_{{\bf{X}}_{1}})$ is a coupling over $({\bf{X}}_{0},{\bf{X}}_{1})$, given by joint distributions on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $p_{{\bf{X}}_{0}}$ and $p_{{\bf{X}}_{1}}$. Lipman et al. (2023) have shown that the minimization of $\mathcal{L}_{\mathrm{CFM}}$ is equivalent to the minimization of $\mathcal{L}_{\mathrm{FM}}$, since their gradients with respect to $\theta$ are equal.

Figure: Figure 1: Overview of the three steps in *Flower*. Starting from an initial sample $\mathbf{x}_{0}\sim p_{\mathbf{X}_{0}}$ at time $t$, the method: Step 1 predicts a flow-consistent destination $\hat{\mathbf{x}}_{1}(\mathbf{x}_{t})$; Step 2 refines this destination using the measurements via a proximal step and associated uncertainty sampling to obtain $\tilde{\mathbf{x}}_{1}(\mathbf{x}_{t},\mathbf{y})$; and Step 3 updates the trajectory along time by interpolating $\tilde{\mathbf{x}}_{1}(\mathbf{x}_{t},\mathbf{y})$ with new noise ${\bm{\epsilon}}\sim p_{\mathbf{X}_{0}}$. The $N$-time repetition of these steps yields the final reconstruction $\mathbf{x}_{1}$.
Refer to caption: images/flower_steps.png

The coupling $\pi$ determines how $({\bf{x}}_{0},{\bf{x}}_{1})$ are paired. With the *independent* (IND) coupling $\pi=p_{{\bf{X}}_{0}}\otimes p_{{\bf{X}}_{1}}$, the training is simple and scalable. However, the resulting interpolated paths can overlap, which may slow down convergence. At the other extreme, the *optimal transport* (OT) coupling $\pi^{\star}\in\arg\min_{\pi\in\Pi(p_{{\bf{X}}_{0}},p_{{\bf{X}}_{1}})}\mathbb{E}_{({\bf{x}}_{0},{\bf{x}}_{1})\sim\pi}\left[\|{\bf{x}}_{1}-{\bf{x}}_{0}\|_{2}^{2}\right]$ produces globally aligned pairs such that straight-line flows approximate displacement interpolation along the Wasserstein-2 geodesic. If the Monge map $T$ satisfying $T_{\#}p_{{\bf{X}}_{0}}=p_{{\bf{X}}_{1}}$ exists and is known, then no training is needed: the sampling ${\bf{x}}_{0}\sim p_{{\bf{X}}_{0}}$ and the computation of ${\bf{x}}_{1}=T({\bf{x}}_{0})$ already generate a sample from $p_{{\bf{X}}_{1}}$. In practice, $T$ is unknown and its approximation is infeasible. A practical compromise is *mini-batch OT*, which solves an entropically regularized OT problem within each batch to compute an approximate coupling $\hat{\pi}$. This improves alignment over independence with moderate computational overhead. For more details on the mathematical background of OT, such as the definition and uniqueness of the Monge map, we refer to Peyré (2025).

The choice of the source distribution $p_{{\bf{X}}_{0}}$ is crucial for effective training and sampling. In practice, $p_{{\bf{X}}_{0}}$ is often chosen as the standard normal distribution $\mathcal{N}({\bf{0}},{\bf{I}})$. With the independent coupling, this simply yields $p_{{\bf{X}}_{t}|{\bf{X}}_{1}={\bf{x}}_{1}}=\mathcal{N}(t{\bf{x}}_{1},(1-t)^{2}{\bf{I}})$. However, the computation of $p_{{\bf{X}}_{t}|{\bf{X}}_{1}={\bf{x}}_{1}}$ in the OT case is challenging due to the mini-batch approach, in which it is difficult to determine the batch a sample ${\bf{x}}_{1}$ came from, as well as its associated OT paths.

### 2.2 Proximal Operator

The *proximal operator* of a proper, lower semi-continuous convex function $f\colon\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ is defined as

$$ $\mathrm{prox}_{f}({\bf{x}})=\operatorname*{arg\,min}_{{\bf{w}}\in\mathbb{R}^{d}}\left(\frac{1}{2}\left\rVert{\bf{w}}-{\bf{x}}\right\rVert_{2}^{2}+f({\bf{w}})\right).$ (8) $$

This operator can be interpreted as a generalized projection of ${\bf{x}}$ onto a set associated with $f$, balancing proximity to ${\bf{x}}$ and regularization by $f$. Proximal operators play a central role in optimization algorithms for that solve inverse problems and are key components of proximal-gradient methods (Bubeck, 2015).

## 3 Method

Let ${\bf{v}}^{\theta}_{t}$ denote a velocity network trained to generate samples from $p_{{\bf{X}}_{1}}$ through flow matching. Therefore, starting from ${\bf{x}}_{0}\sim p_{{\bf{X}}_{0}}$, we get a sample ${\bf{x}}_{1}\sim p_{{\bf{X}}_{1}}$ if we perform $N$ iterations of the update equation

$$ ${\bf{x}}_{t+\Delta t}={\bf{x}}_{t}+\Delta t\ {\bf{v}}_{t}^{\theta}({\bf{x}}_{t}),$ (9) $$

with $\Delta t=\frac{1}{N}$.
We aim to use the pre-trained velocity network ${\bf{v}}^{\theta}_{t}$ to generate solutions ${\bf{x}}_{1}$ that are consistent with the flow and the linear forward model ${\bf{y}}={\bf{H}}{\bf{x}}_{1}+{\bf{n}}$ for ${\bf{n}}\sim\mathcal{N}({\bf{0}},\sigma_{n}^{2}{\bf{I}})$, as described in ([1](https://arxiv.org/html/2509.26287v1#S1.E1)). To achieve this goal, we introduce *Flower* which, given the measurements ${\bf{y}}$, modifies the unconditional flow path of ([9](https://arxiv.org/html/2509.26287v1#S3.E9)) and outputs ${\bf{x}}_{1}$ by iterating $N$ times over three steps. We first introduce these steps and then theoretically establish how and under what assumptions *Flower* generates a sample ${\bf{x}}_{1}$ of the conditional posterior $p_{{\bf{X}}_{1}\mid{\bf{Y}}={\bf{y}}}$. The three steps are as follows.

- 1.
Flow-consistent destination estimation
$\hat{{\bf{x}}}_{1}({\bf{x}}_{t})={\bf{x}}_{t}+(1-t){\bf{v}}_{t}^{\theta}({\bf{x}}_{t}).$
(10)
- 2.
Measurement-aware destination refinement
$\tilde{{\bf{x}}}_{1}({\bf{x}}_{t},{\bf{y}})={{\bm{\mu}}}_{t}({\bf{x}}_{t},{\bf{y}})+\gamma{\bm{\kappa}}_{t}$
(11)
for
${{\bm{\mu}}}_{t}({\bf{x}}_{t},{\bf{y}})=\mathrm{prox}_{\nu_{t}^{2}F_{{\bf{y}}}}\left(\hat{{\bf{x}}}_{1}\left({\bf{x}}_{t}\right)\right),\quad{\bm{\kappa}}_{t}\sim\mathcal{N}({\bf{0}},{\bf{\Sigma}}_{t}),$
(12)
where $F_{{\bf{y}}}({\bf{x}})=\frac{1}{2\sigma_{n}^{2}}\left\rVert{\bf{H}}{\bf{x}}-{\bf{y}}\right\rVert_{2}^{2}$ and $\mathrm{prox}_{\nu_{t}^{2}F}$ denotes the proximal operator of $\nu_{t}^{2}F_{{\bf{y}}}$, as defined in ([8](https://arxiv.org/html/2509.26287v1#S2.E8)).
We have that $\nu_{t}=\frac{(1-t)}{\sqrt{t^{2}+(1-t)^{2}}}$ and that ${\bf{\Sigma}}_{t}=\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right)^{-1}$. The hyperparameter $\gamma\in\{0,1\}$ controls the consideration of the uncertainty of the destination refinement step.
- 3.
Movement in time
${\bf{x}}_{t+\Delta t}=(1-t-\Delta t){\bm{\epsilon}}+(t+\Delta t)\tilde{{\bf{x}}}_{1}({\bf{x}}_{t},{\bf{y}}),$
(13)
where ${\bm{\epsilon}}$ is newly sampled from $p_{{\bf{X}}_{0}}$ at each iteration.

Here, $\Delta t=\tfrac{1}{N}$ and the scheme is initialized with a sample ${\bf{x}}_{0}\sim p_{{\bf{X}}_{0}}$. In Figure [1](https://arxiv.org/html/2509.26287v1#S2.F1), we present a visual illustration of these three steps. We also summarize these steps in Algorithm [1](https://arxiv.org/html/2509.26287v1#alg1) of the Appendix.

We now interpret *Flower* through a Bayesian lens. We assume that, at each iteration, the three steps collectively draw ${\bf{x}}_{t+\Delta t}$ from the transition distribution $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$. Under this assumption and with proper initialization, the procedure performs ancestral sampling along the conditional trajectory. By induction, we obtain ${\bf{x}}_{t+\Delta t}\sim p_{{\bf{X}}_{t+\Delta t}\mid{\bf{Y}}={\bf{y}}}$. We formalize this in Theorem [1](https://arxiv.org/html/2509.26287v1#Thmtheorem1), which in turn implies that the final sample ${\bf{x}}_{1}$ produced by *Flower* follows the desired posterior $p_{{\bf{X}}_{1}\mid{\bf{Y}}={\bf{y}}}$. We then detail how, in practice, the three steps realize a draw from $\tilde{p}_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ (which serves as an approximation to $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$).

###### Theorem 1 .

Let ${\bf{x}}_{0}$ be a sample from $p_{{\bf{X}}_{0}\mid{\bf{Y}}={\bf{y}}}$. If ${\bf{x}}_{t}$ is a sample from $p_{{\bf{X}}_{t}\mid{\bf{Y}}={\bf{y}}}$, then the sample ${\bf{x}}_{t+\Delta t}$ from $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ follows $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{Y}}={\bf{y}}}$.

To establish this result, we first recall ancestral sampling. It is a procedure that enables us to draw samples from a marginal distribution $p_{{\bf{Z}}_{K}}$ when the full joint distribution over a sequence of variables ${\bf{Z}}=({\bf{Z}}_{1},{\bf{Z}}_{2},\dots,{\bf{Z}}_{K})$ is defined via a chain of conditional densities and when the marginal distribution of ${\bf{Z}}_{K}$ can be written as

$$ $p_{{\bf{Z}}_{K}}({\bf{z}}_{K})=\int_{{\bf{z}}_{K-1}}\cdots\int_{{\bf{z}}_{1}}p_{{\bf{Z}}_{K}\mid{\bf{Z}}_{K-1}={\bf{z}}_{K-1}}({\bf{z}}_{t})\cdots p_{{\bf{Z}}_{2}\mid{\bf{Z}}_{1}={\bf{z}}_{1}}({\bf{z}}_{2})\,p_{{\bf{Z}}_{1}}({\bf{z}}_{1})\,\mathrm{d}{\bf{z}}_{1}\cdots\mathrm{d}{\bf{z}}_{K-1}.$ (14) $$

Ancestral sampling offers a practical way to generate samples from ${\bf{p}}_{Z_{K}}$ without direct evaluation of this integral. The process samples sequentially from the distributions

$$ ${\bf{z}}_{1}\sim p_{{\bf{Z}}_{1}},\quad{\bf{z}}_{2}\sim p_{{\bf{Z}}_{2}\mid{\bf{Z}}_{1}={\bf{z}}_{1}},\quad\dots,\quad{\bf{z}}_{K}\sim p_{{\bf{Z}}_{K}\mid{\bf{Z}}_{K-1}={\bf{z}}_{K-1}}.$ (15) $$

By following this sequence, we obtain a valid sample from the marginal ${\bf{z}}_{K}\sim p_{{\bf{Z}}_{K}}$.

###### Proof of Theorem 1 .

By using the marginal distributions and the general chain rule for joint probability, we obtain

$$ $\displaystyle p_{{\bf{X}}_{t+\Delta t}|{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})$ $\displaystyle=\int_{{\bf{x}}_{t}}p_{{\bf{X}}_{t},{\bf{X}}_{t+\Delta t}|{\bf{Y}}={\bf{y}}}({\bf{x}}_{t},{\bf{x}}_{t+\Delta t})\mathrm{d}{\bf{x}}_{t}.$ $\displaystyle=\int_{{\bf{x}}_{t}}p_{{\bf{X}}_{t}|{\bf{Y}}={\bf{y}}}({\bf{x}}_{t})p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})\mathrm{d}{\bf{x}}_{t}.$ (16) $$

It then suffices to follow the ancestral sampling procedure with $K=2$, ${\bf{Z}}_{1}={\bf{X}}_{t}|{\bf{Y}}={\bf{y}}$, and ${\bf{Z}}_{2}={\bf{X}}_{t+\Delta t}|{\bf{Y}}={\bf{y}}$ to complete the proof.
∎

###### Remark 1 .

For the inductive argument to hold, *Flower* must be initialized with a sample from the conditional distribution $p_{{\bf{X}}_{0}\mid{\bf{Y}}={\bf{y}}}$. When ${\bf{X}}_{0}$ and ${\bf{X}}_{1}$ are assumed to be independent, this reduces to a sampling from the unconditional prior $p_{{\bf{X}}_{0}}$, which is often chosen as $\mathcal{N}({\bf{0}},{\bf{I}})$.

Theorem [1](https://arxiv.org/html/2509.26287v1#Thmtheorem1) presupposes the existence of an ancestral-sampling scheme to generate samples from $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{Y}}={\bf{y}}}$. This scheme requires a sampling from the transition distribution $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$. We now describe how to realize this transition in practice. We proceed sequentially and explain the details of each step of *Flower*.

First, under the assumption that ${\bf{v}}_{t}^{\theta}$ is the optimal velocity network, we show in Proposition [1](https://arxiv.org/html/2509.26287v1#Thmproposition1) that the predicted $\hat{{\bf{x}}}_{1}({\bf{x}}_{t})$ in Step 1 equals the conditional expectation $\mathbb{E}[{\bf{X}}_{1}\mid{\bf{X}}_{t}={\bf{x}}_{t}]$. The proof is provided in Appendix [7.2.1](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS1).

###### Proposition 1 .

If ${\bf{v}}^{\theta}({\bf{x}}_{t},t)={\bf{v}}^{*}_{t}({\bf{x}})$ is a pre-trained velocity vector field that minimizes the conditional flow-matching loss, then

$$ $\hat{{\bf{x}}}_{1}({\bf{x}}_{t})=\mathbb{E}[{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}]={\bf{x}}_{t}+(1-t){\bf{v}}^{\theta}({\bf{x}}_{t},t).$ (17) $$

Since the distribution $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}$ is not directly available, we propose to approximate it with

$$ $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}=\mathcal{N}(\hat{{\bf{x}}}_{1}({\bf{x}}_{t}),\nu_{t}^{2}{\bf{I}}),$ (18) $$

an isotropic Gaussian distribution centered at $\hat{{\bf{x}}}_{1}({\bf{x}}_{t})=\mathbb{E}[{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}]$ with a time-varying covariance. As $t\to 1$, the distribution $p_{{\bf{X}}_{t}}$ approaches the target $p_{{\bf{X}}_{1}}$. For $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}$ to be consistent with this property, $\nu_{t}$ should anneal in time. We choose

$$ $\nu_{t}=\frac{(1-t)}{\sqrt{t^{2}+(1-t)^{2}}},$ (19) $$

which results in the valid covariance when $p_{{\bf{X}}_{1}}$ is a standard Gaussian distribution. Our approximation is indeed the $\Pi$GDM approximation proposed by Song et al. (2023) within diffusion solvers and later by Pokle et al. (2024) for flow matching. However, instead of having a score-based interpretation and using this approximation to obtain $\nabla_{{\bf{x}}_{t}}\log\tilde{p}_{{\bf{Y}}\mid{\bf{X}}_{t}={\bf{x}}_{t}}$, we propose to sample $\tilde{{\bf{x}}}_{1}({\bf{x}}_{t},{\bf{y}})$ from the distribution $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ that approximates ${p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$. To this end, we show in Proposition [2](https://arxiv.org/html/2509.26287v1#Thmproposition2) that $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ is indeed a Gaussian distribution, using the $\Pi$GDM approximation and the forward model of ([1](https://arxiv.org/html/2509.26287v1#S1.E1)). The proof is provided in Appendix [7.2.2](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS2).

###### Proposition 2 .

Suppose that $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}={\mathcal{N}}(\hat{{\bf{x}}}_{1}({\bf{x}}_{t}),\nu_{t}^{2}{\bf{I}})$ ($\Pi$GDM approximation) and $p_{{\bf{Y}}\mid{\bf{X}}_{1}={\bf{x}}_{1}}=\mathcal{N}({\bf{H}}{\bf{x}}_{1},\sigma_{n}^{2}{\bf{I}})$ (measurement operation). Then, $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}={\mathcal{N}}({{\bm{\mu}}}_{t}({\bf{x}}_{t},{\bf{y}}),{\bf{\Sigma}}_{t})$, where

$$ $\displaystyle{{\bm{\mu}}}_{t}({\bf{x}}_{t},{\bf{y}})$ $\displaystyle=\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right)^{-1}\left(\nu_{t}^{-2}\hat{{\bf{x}}}_{1}({\bf{x}}_{t})+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{y}}\right),$ (20) $\displaystyle{\bf{\Sigma}}_{t}$ $\displaystyle=\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right)^{-1}.$ (21) $$

Proposition [2](https://arxiv.org/html/2509.26287v1#Thmproposition2) allows us to sample from $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ as provided in Step 2 of *Flower* using the re-parameterization trick in ([11](https://arxiv.org/html/2509.26287v1#S3.E11)). However, the ${\bm{\mu}}_{t}({\bf{x}}_{t},{\bf{y}})$ in Step 2 (see ([12](https://arxiv.org/html/2509.26287v1#S3.E12))) is described using a proximal operator which differs from ([20](https://arxiv.org/html/2509.26287v1#S3.E20)). It is easy to verify the equivalence between the two, through the fact that ${\bm{\mu}}_{t}({\bf{x}}_{t},{\bf{y}})$ of ([20](https://arxiv.org/html/2509.26287v1#S3.E20)) can be written as the solution to the minimization problem

$$ $\min_{{\bf{x}}\in\mathbb{R}^{d}}\left(\frac{1}{2\sigma_{n}^{2}}\left\rVert{\bf{H}}{\bf{x}}-{\bf{y}}\right\rVert_{2}^{2}+\frac{1}{2\nu_{t}^{2}}\left\rVert{\bf{x}}-\hat{{\bf{x}}}_{1}({\bf{x}}_{t})\right\rVert_{2}^{2}\right).$ (22) $$

This directly results in ${{\bm{\mu}}}_{t}({\bf{x}}_{t},{\bf{y}})=\mathrm{prox}_{\nu_{t}^{2}F}\left(\hat{{\bf{x}}}_{1}\left({\bf{x}}_{t}\right)\right)$ for $F_{{\bf{y}}}({\bf{x}})=\frac{1}{2\sigma_{n}^{2}}\left\rVert{\bf{H}}{\bf{x}}-{\bf{y}}\right\rVert_{2}^{2}$ under the definition of the proximal operator in ([8](https://arxiv.org/html/2509.26287v1#S2.E8)). Moreover, the sampling from the anisotropic Gaussian $\mathcal{N}({\bf{0}},{\bf{\Sigma}}_{t})$ is not trivial; however, if we sample two independent ${\bm{\epsilon}}_{1}\in\mathbb{R}^{M}$ and ${\bm{\epsilon}}_{2}\in\mathbb{R}^{d}$ from standard Gaussian distributions, then we verify in Appendix [7.2.4](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS4) that

$$ ${\bm{\kappa}}_{t}={\bf{\Sigma}}_{t}({\nu_{t}}^{-1}{{\bm{\epsilon}}_{1}}+{\sigma_{n}}^{-1}{\bf{H}}^{\top}{\bm{\epsilon}}_{2})$ (23) $$

follows $\mathcal{N}({\bf{0}},{\bf{\Sigma}}_{t})$.

Step 3 of *Flower* aims to sample the distribution $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$, which is also what we required for our ancestral-sampling procedure to hold. We now show that if we have a sample $\tilde{{\bf{x}}}_{1}({\bf{x}}_{t},{\bf{y}})$ from $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$, then we could obtain a sample from the distribution $\tilde{p}_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ using ([13](https://arxiv.org/html/2509.26287v1#S3.E13)). We first compute $\tilde{p}_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ under the assumption that $p_{{\bf{X}}_{0}}$ is independent of $p_{{\bf{X}}_{1}}$ in Proposition [7.2.3](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS3) which we prove in Appendix [7.2.3](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS3).

###### Proposition 3 .

If from the pre-trained flow matching we have that $p_{{\bf{X}}_{0}}=\mathcal{N}({\bf{0}},{\bf{I}})$, if $p_{{\bf{X}}_{0}}$ is independent of $p_{{\bf{X}}_{1}}$, and if $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})=\mathcal{N}\left({\bm{\mu}}_{t}({\bf{x}}_{t},{\bf{y}}),{\bf{\Sigma}}_{t}\right)$, then it holds that

$$ $\tilde{p}_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}=\mathcal{N}\left((t+\Delta t){\bm{\mu}}_{t},(t+\Delta t)^{2}{\bf{\Sigma}}_{t}+(1-t-\Delta t)^{2}{\bf{I}}\right).$ (24) $$

From Proposition [7.2.3](https://arxiv.org/html/2509.26287v1#S7.SS2.SSS3) and by the means of the re-parametrization trick, it is easy to verify that

$$ ${\bf{x}}_{t+\Delta t}=(1-t-\Delta t){\bm{\epsilon}}+(t+\Delta t)\tilde{{\bf{x}}}_{1}({\bf{x}}_{t},{\bf{y}})$ (25) $$

follows $\tilde{p}_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ given that $\tilde{{\bf{x}}}\sim\mathcal{N}\left({\bm{\mu}}_{t}({\bf{x}}_{t},{\bf{y}}),{\bf{\Sigma}}_{t}\right)$, which happens in Step 3 of *Flower*.

###### Theorem 1 .

###### Proof of Theorem 1 .

###### Remark 1 .

###### Proposition 1 .

###### Proposition 2 .

###### Proposition 3 .

##### Application to other source–target couplings.

In our Bayesian justification, we assumed that $p_{{\bf{X}}_{0}}$ and $p_{{\bf{X}}_{1}}$ are independent. This assumption excludes, for example, the mini-batch optimal-transport coupling. Nevertheless, in practice, *Flower* can still be applied in such settings and interpreted in a PnP manner, as we discuss next. We also address this empirically in our numerical results.

##### Relation to plug-and-play.

The iterative structure of *Flower* closely resembles PnP methods: Step 1 acts as a denoising step, while Step 2 enforces data consistency. From this viewpoint, Step 3 can be seen as a re-projection onto the flow trajectory, as discussed in Martin et al. (2025). However, rather than relying solely on this interpretation, we provide a Bayesian justification of the procedure. This perspective highlights a conceptual link between PnP methods and posterior sampling.

##### The role of γ \gamma .

For exact posterior sampling to hold under our approximations, $\gamma$ should be set to one. Interestingly, in practice we find that the choice $\gamma=0$ (i.e., ignoring the uncertainty in the destination refinement step) yields a better reconstruction quality. We provide further discussion on this effect with our numerical results in Section [5](https://arxiv.org/html/2509.26287v1#S5).

## 4 Related Works

An extensive body of work adapts pre-trained diffusion or flow priors to inverse problems by modifying the dynamics to approximate the conditional posterior. We first review diffusion-based solvers, then flow-based ones. Throughout this section, we highlight how *Flower* differs from similar methods. We use the flow-matching notation with source ${\bf{X}}_{0}$ and target ${\bf{X}}_{1}$.

Among diffusion solvers, DPS (Chung et al., 2023) approximates $p_{{\bf{Y}}\mid{\bf{X}}_{t}={\bf{x}}_{t}}$ by $p_{{\bf{Y}}\mid{\bf{X}}_{1}=\hat{{\bf{x}}}_{1}({\bf{x}}_{t})}$, where $\hat{{\bf{x}}}_{1}({\bf{x}}_{t})$ is the diffusion-based denoised version of ${\bf{x}}_{t}$, which leads to a gradient correction to the diffusion dynamics. $\Pi$GDM (Song et al., 2023) approximates $\tilde{p}_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}=\mathcal{N}(\hat{{\bf{x}}}_{1}({\bf{x}}_{t}),\nu_{t}^{2}I)$ for some time-annealing $\nu_{t}$ and replaces the gradient correction of DPS with a pseudoinverse-based update. *Flower* adopts the same approximation as $\Pi$GDM, but the subsequent steps differ. DDS (Chung et al., 2024) shares the same perspective as DPS but replaces the gradient with a proximal step motivated by a manifold-preserving gradient perspective. DiffPIR (Zhu et al., 2023) arrives at a very similar structure through half-quadratic splitting, alternating proximal data updates with diffusion denoising. Both DDS and DiffPIR are structurally close to *Flower* but, unlike *Flower*, they lack the Bayesian justification that interprets the updates as posterior sampling. DAPS (Zhang et al., 2025) also uses ancestral sampling similar to *Flower*, but instead of directly computing $\tilde{p}_{{\bf{X}}_{1}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}=y}$, it applies Langevin updates for its evaluation, which is more computationally demanding but extends naturally to nonlinear inverse problems.

In the flow-matching domain, OT-ODE (Pokle et al., 2024) employs a $\Pi$GDM-based approximation, similar in spirit to ours. However, instead of adopting our ancestral sampling scheme, they, similar to the approach of $\Pi$GDM, approximate the score of the conditional distribution $\tilde{p}_{{\bf{Y}}\mid{\bf{X}}_{t}={\bf{x}}_{t}}$ in order to construct the new velocity field. Therefore, while we rely on the same approximation for $\tilde{p}_{{\bf{X}}_{1}\mid{\bf{X}}_{t}={\bf{x}}_{t}}$, our methods act in different ways.
Flow-Priors (Zhang et al., 2024) tackle inverse problems by reformulating the MAP objective as a sequence of time-dependent MAP subproblems with closed-form evaluations by taking advantage of the velocity network. However, their approach relies on the computation of $\mathrm{Tr}\,\nabla{\bf{v}}^{\theta}_{t}$, which is costly.
D-Flow (Ben-Hamu et al., 2024) adopts an implicit regularization strategy, replacing the data-fidelity objective ${\bf{x}}\mapsto\|{\bf{H}}{\bf{x}}-{\bf{y}}\|^{2}$ with a latent loss
${\bf{z}}\mapsto\|{\bf{H}}(f(1,{\bf{z}}))-{\bf{y}}\|^{2}$, where $f$ is the solution of the flow ODE.
The latent loss is non-convex with an implicit regularization effect that prevents convergence to trivial solutions. The optimization is performed by back-propagation through ODE solutions, which is computationally demanding.
PnP-Flow (Martin et al., 2025) introduces a PnP framework that employs the velocity network as a denoiser. Its update steps are similar to ours, but we replace their gradient update with a proximal operation, which leads to improved reconstruction quality. Moreover, *Flower* offers a Bayesian justification of the process while PnP-Flow is purely plug-and-play.

## 5 Numerical Results

We present our numerical results in two sections. In Section [5.1](https://arxiv.org/html/2509.26287v1#S5.SS1), we validate the Bayesian interpretation of *Flower* through a toy experiment with Gaussian mixtures, where ground-truth posterior samples are computable. In Section [5.2](https://arxiv.org/html/2509.26287v1#S5.SS2), we benchmark *Flower* against state-of-the-art flow-based inverse solvers across a range of tasks.

Figure: Figure 2: Temporal evolution of 2D *Flower* and comparison with true posterior for noise variance $\sigma_{n}=0.25$.
Refer to caption: images/flower_2d_in_time.png

### 5.1 Toy Experiment

The goal of this experiment is to validate our proposed sampling perspective in a setting where samples from the true posterior are available. To this end, we consider a Gaussian mixture model (GMM) as the target distribution of the data ${\bf{X}}\in\mathbb{R}^{2}$, with details provided in Appendix [7.3.1](https://arxiv.org/html/2509.26287v1#S7.SS3.SSS1.Px1). The forward measurement model consists of a single measurement vector ${\bf{h}}\in\mathbb{R}^{d}$ (i.e., the forward operator is ${\bf{H}}={\bf{h}}^{\top}$) corrupted by additive white Gaussian noise $n\sim{\mathcal{N}}(0,\sigma_{n}^{2})$, which results in $y={\bf{h}}^{\top}{\bf{x}}+n$. In this setup, the posterior distribution $p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}$ is itself a GMM with known parameters, whose analytical expression is given in Appendix [7.3.1](https://arxiv.org/html/2509.26287v1#S7.SS3.SSS1.Px2).
Geometrically, the noiseless measurement $y={\bf{h}}^{\top}{\bf{x}}$ defines a line in the two-dimensional plane with normal vector ${\bf{h}}$. Consequently, when sampling from the posterior $p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}$, we expect the samples to concentrate on the portions of this line that intersect regions where the prior distribution $p_{{\bf{X}}}$ has high density.

We trained the unconditional flow-matching vector field using the source $p_{{\bf{X}}_{0}}={\mathcal{N}}({\bf{0}},{\bf{I}})$ and target $p_{{\bf{X}}_{1}}=p_{{\bf{X}}}$, with further details provided in Appendix [7.3.1](https://arxiv.org/html/2509.26287v1#S7.SS3.SSS1.Px3). We ran *Flower* with $N=1000$ iterations, which corresponds to the step size $\Delta t=0.001$. The results are shown in Figure [2](https://arxiv.org/html/2509.26287v1#S5.F2), where we used ${\bf{h}}^{\top}=[1.5,1.5]$, $\sigma_{n}=0.25$, and the observation $y=1$. In the figure, we report the solution paths of *Flower* with $\gamma\in\{0,1\}$ alongside true posterior samples. We observe, that when $\gamma=1$ (the setting required by our theory), *Flower* successfully recovers the samples at $t=1$, which closely resemble the true posterior. When $\gamma=0$ (a configuration in which the uncertainty in the destination estimation step is ignored), *Flower* fails to capture samples from the tails of the distribution. Additional visual results for a different line and noise level, as well as a step-by-step illustration of *Flower* through time, are provided in Appendix [7.3.1](https://arxiv.org/html/2509.26287v1#S7.SS3.SSS1.Px4).

### 5.2 Benchmark Experiments

The goal of this section is to benchmark our method against other flow-matching-based solvers for inverse problems. For fair comparisons, we adopt the benchmark introduced by Martin et al. (2025), which also includes state-of-the-art PnP and diffusion models. We describe the datasets and experimental setup for completeness, present quantitative results in Tables [1](https://arxiv.org/html/2509.26287v1#S5.T1) and [2](https://arxiv.org/html/2509.26287v1#S5.T2), and provide qualitative examples in Figure [3](https://arxiv.org/html/2509.26287v1#S6.F3). Finally, we discuss the key observations, highlighting the performance of *Flower* and its empirical considerations.

We use two datasets for our numerical comparisons. First, we use $(128\times 128)$ human-face images from Yang et al. (2015), denoted by CelebA. Second, we use resized $(256\times 256)$ cat images from Choi et al. (2020), denoted by AFHQ-Cat. We normalized all images to the range $[-1,1]$. We train on the full training sets of both datasets. For testing, we use 100 images from the test set of each dataset. Since AFHQ-Cat does not include a validation set, we construct one using 32 images from the test set as suggested by Martin et al. (2025).

We compare *Flower* against flow-matching solvers OT-ODE (Pokle et al., 2024), D-Flow (Ben-Hamu et al., 2024), Flow-Priors (Zhang et al., 2024), and PnP-Flow (Martin et al., 2025), as well as two other baselines: PnP-GS (Hurault et al., 2022b), a state-of-the-art plug-and-play method, and DiffPIR Zhu et al. (2023), a diffusion-based inverse solver. All models (except DiffPIR) use the same U-Net backbone (Ronneberger et al., 2015) trained with Mini-Batch OT Flow Matching (Tong et al., 2024) and a Gaussian latent prior. Pre-trained weights for the flow models and PnP-GS are taken from Martin et al. (2025), trained with learning rate $10^{-4}$: on CelebA for 200 epochs (batch size 128) and on AFHQ-Cat for 400 epochs (batch size 64). For *Flower*, we additionally train a variant without latent–target coupling (Flower-IND) using the same hyperparameters, which corresponds to our theoretical setting. While Flower-IND achieves higher performance (see Appendix [7.4.1](https://arxiv.org/html/2509.26287v1#S7.SS4.SSS1)), we primarily report Flower-OT for consistency with other flow-matching baselines. Training DiffPIR with this backbone proved ineffective due to limited capacity, so following Martin et al. (2025), we adopt a pretrained model Choi et al. (2021) from the DeepInv library (Tachella et al., 2023), originally trained on FFHQ (Karras et al., 2019). This introduces some mismatch but provides the fairest diffusion-based baseline. Note that we use the latest checkpoints from Martin et al. (2025), but our averaging strategy differs from theirs. In Martin et al. (2025), results are reported by grouping four images into one batch and then averaging across 25 such batches. In contrast, we recomputed the results using 100 independent averages over the images themselves. Consequently, our reported numbers differ from those in Martin et al. (2025).

We evaluate performance on five restoration tasks: (i) denoising with Gaussian noise with $\sigma_{n}=0.2$; (ii) deblurring with a $61\times 61$ Gaussian kernel ($\sigma_{b}=1.0$ for CelebA, $\sigma_{b}=3.0$ for AFHQ-Cat) and additive noise $\sigma_{n}=0.05$; (iii) super-resolution ($2\times$ downsampling for CelebA and $4\times$ for AFHQ-Cat, with $\sigma_{n}=0.05$); (iv) random inpainting with $70\%$ of pixels removed ($\sigma_{n}=0.01$); and (v) box inpainting with a centered $40\times 40$ mask for CelebA and $80\times 80$ mask for AFHQ-Cat ($\sigma_{n}=0.05$). To report quantitative results, we use peak signal-to-noise ratio (PSNR), structural similarity index measure (SSIM), and learned perceptual image patch similarity (LPIPS). Note that for PSNR and SSIM, higher values indicate better performance, while for LPIPS, lower values are better.

To ensure fair comparisons, we adopt the optimal hyperparameters reported in Martin et al. (2025) for each method, obtained via grid search on the validation set to maximize PSNR. For PnP-Flow, we report two variants: PnP-Flow1 and PnP-Flow5, which apply one and five evaluations of the velocity network per denoising step, respectively. For *Flower*, we follow the same procedure, reporting in Tables [1](https://arxiv.org/html/2509.26287v1#S5.T1) and [2](https://arxiv.org/html/2509.26287v1#S5.T2) either the output of a single evaluation (Flower1-OT) or the average of five evaluations (Flower5-OT). A key property of *Flower* is that, apart from the number $N$ of iterations and the knowledge of the noise level $\sigma_{n}$, it uses the same hyperparameters across different inverse problems, unlike other flow models. In particular, aside from $N$, the only hyperparameter of *Flower* is $\gamma$, which controls the uncertainty of the destination refinement. Across all setups, $\gamma=0$ yields higher reconstruction quality. As discussed in Appendix [7.4.2](https://arxiv.org/html/2509.26287v1#S7.SS4.SSS2), the choice $\gamma=1$ produces samples that appear realistic but requires the averaging of multiple runs to achieve competitive PSNR, whereas $\gamma=0$ attains better quality with fewer averages. This observation is consistent with our toy experiments, where $\gamma=0$ encouraged sampling from higher-probability regions. For $N$, we always match the number of steps used by our main competitor, PnP-Flow. The full hyperparameter details are reported in Appendix [7.4.4](https://arxiv.org/html/2509.26287v1#S7.SS4.SSS4).

Key Observations.
On CelebA (Table [1](https://arxiv.org/html/2509.26287v1#S5.T1)), *Flower* achieves the best or near-best results across all tasks, with clear gains in deblurring and box inpainting. The five-step averaging further improves the results. On AFHQ-Cat (Table [2](https://arxiv.org/html/2509.26287v1#S5.T2)), *Flower* remains highly competitive and outperforms baselines in deblurring, box inpainting, and random inpainting, while PnP-GS is strongest in denoising. In these tables, bold numbers indicate the best results among single-average results of methods. Underlined numbers indicate the second best. Blue numbers highlight the overall best across all methods. We illustrate in Figure [3](https://arxiv.org/html/2509.26287v1#S6.F3) representative reconstructions across denoising, deblurring, super-resolution, and inpainting tasks. Compared to OT-ODE, D-Flow, and Flow-Priors, *Flower* consistently produces fewer artifacts, while also avoiding the over-smoothing often observed in PnP-Flow. These visual trends align with the quantitative results and highlight the robustness of *Flower* across diverse degradations. In Figure [4](https://arxiv.org/html/2509.26287v1#S6.F4), we illustrate the solution path of *Flower* for the box inpainting task shown in Figure [3](https://arxiv.org/html/2509.26287v1#S6.F3). As expected, Step 1 produces flow-based denoised images, while Step 2 enforces consistency with the measurements. In this specific box-inpainting setup, Step 2 primarily aligns the region outside the box with the measurements and preserves the result of Step 1 inside the box. Step 3 then mixes the refined destination with fresh source noise; as $t$ increases, the noise decreases and the reconstruction emerges. The injected noise is essential to prevent the velocity network from getting stuck at the previous iterate and to allow it to predict improved destinations.
As shown in Table [4](https://arxiv.org/html/2509.26287v1#S7.T4) of Appendix [7.4.3](https://arxiv.org/html/2509.26287v1#S7.SS4.SSS3), *Flower* has a runtime that is similar to PnP-Flow and OT-ODE, with only a slight overhead relative to PnP-Flow due to the proximal-projection step replacing a simple gradient update, while requiring the same minimal memory. In contrast, D-Flow and Flow-Priors are substantially slower and more memory-intensive.

**Table 1: Results on 100 test images of the dataset CelebA.**
| Method | Denoising | Deblurring | Super-resolution | Random inpainting | Box inpainting |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS |  |
| Degraded | 20.00 | 0.348 | 0.372 | 27.83 | 0.740 | 0.126 | 10.26 | 0.183 | 0.827 | 11.95 | 0.196 | 1.041 | 22.27 | 0.742 | 0.214 |
| PnP-GS | 32.64 | 0.910 | 0.035 | 34.03 | 0.924 | 0.041 | 31.31 | 0.892 | 0.064 | 29.22 | 0.875 | 0.070 | - | - | - |
| DiffPIR | 31.20 | 0.885 | 0.060 | 32.77 | 0.912 | 0.060 | 31.52 | 0.895 | 0.033 | 31.74 | 0.917 | 0.025 | - | - | - |
| OT-ODE | 30.54 | 0.859 | 0.032 | 33.01 | 0.921 | 0.029 | 31.46 | 0.907 | 0.025 | 28.68 | 0.871 | 0.051 | 29.40 | 0.920 | 0.038 |
| D-Flow | 26.04 | 0.607 | 0.092 | 31.25 | 0.854 | 0.038 | 30.47 | 0.843 | 0.026 | 33.67 | 0.943 | 0.015 | 30.70 | 0.899 | 0.026 |
| Flow-Priors | 29.34 | 0.768 | 0.134 | 31.54 | 0.858 | 0.056 | 28.35 | 0.713 | 0.102 | 32.88 | 0.871 | 0.019 | 30.07 | 0.858 | 0.048 |
| PnP-Flow1 | 31.80 | 0.905 | 0.044 | 34.48 | 0.936 | 0.040 | 31.09 | 0.902 | 0.045 | 33.05 | 0.944 | 0.018 | 30.47 | 0.933 | 0.037 |
| Flower1-OT (ours) | 32.28 | 0.914 | 0.034 | 34.98 | 0.947 | 0.026 | 32.36 | 0.923 | 0.034 | 33.08 | 0.944 | 0.018 | 31.19 | 0.945 | 0.022 |
| PnP-Flow5 | 32.30 | 0.911 | 0.056 | 34.80 | 0.940 | 0.047 | 31.49 | 0.906 | 0.056 | 33.98 | 0.953 | 0.022 | 31.09 | 0.940 | 0.043 |
| Flower5-OT (ours) | 33.14 | 0.926 | 0.038 | 35.67 | 0.954 | 0.032 | 33.09 | 0.932 | 0.040 | 33.95 | 0.953 | 0.020 | 31.87 | 0.952 | 0.023 |

**Table 2: Results on 100 test images of the dataset AFHQ-Cat.**
| Method | Denoising | Deblurring | Super-resolution | Random inpainting | Box inpainting |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS |  |
| Degraded | 20.00 | 0.314 | 0.509 | 23.94 | 0.517 | 0.444 | 11.70 | 0.208 | 0.873 | 13.36 | 0.223 | 1.081 | 21.80 | 0.740 | 0.198 |
| PnP-GS | 32.58 | 0.894 | 0.072 | 27.91 | 0.753 | 0.349 | 24.15 | 0.632 | 0.362 | 29.42 | 0.836 | 0.126 | - | - | - |
| DiffPIR | 30.58 | 0.835 | 0.189 | 27.56 | 0.728 | 0.342 | 23.65 | 0.624 | 0.402 | 31.70 | 0.881 | 0.062 | - | - | - |
| OT-ODE | 30.03 | 0.815 | 0.076 | 27.06 | 0.713 | 0.123 | 25.91 | 0.716 | 0.108 | 29.40 | 0.839 | 0.090 | 24.62 | 0.875 | 0.085 |
| D-Flow | 26.13 | 0.574 | 0.175 | 27.82 | 0.721 | 0.164 | 24.64 | 0.601 | 0.190 | 32.20 | 0.894 | 0.040 | 26.26 | 0.842 | 0.077 |
| Flow-Priors | 29.41 | 0.763 | 0.153 | 26.47 | 0.700 | 0.181 | 23.51 | 0.570 | 0.272 | 32.37 | 0.906 | 0.047 | 26.20 | 0.818 | 0.118 |
| PnP-Flow1 | 31.18 | 0.863 | 0.135 | 27.87 | 0.760 | 0.304 | 26.94 | 0.763 | 0.171 | 33.00 | 0.918 | 0.037 | 26.00 | 0.897 | 0.103 |
| Flower1-OT (ours) | 31.69 | 0.879 | 0.102 | 28.64 | 0.775 | 0.255 | 26.23 | 0.741 | 0.272 | 32.97 | 0.918 | 0.040 | 26.19 | 0.915 | 0.063 |
| PnP-Flow5 | 31.43 | 0.864 | 0.168 | 28.19 | 0.766 | 0.332 | 27.37 | 0.774 | 0.183 | 33.75 | 0.929 | 0.048 | 26.68 | 0.901 | 0.120 |
| Flower5-OT (ours) | 32.35 | 0.891 | 0.116 | 28.97 | 0.784 | 0.283 | 26.57 | 0.075 | 0.282 | 33.70 | 0.927 | 0.045 | 26.88 | 0.922 | 0.066 |

## 6 Conclusion

We introduced *Flower*, a method that leverages pre-trained flow-matching models to solve inverse problems through a simple three-step iterative procedure. By combining flow-consistent predictions, measurement-aware refinement, and time evolution, *Flower* provides a principled Bayesian interpretation while retaining the plug-and-play flexibility of existing approaches. Our analysis established the conditions under which the method recovers samples from the conditional posterior. Our experiments demonstrated both validity on toy data and state-of-the-art performance across diverse inverse problems.

Figure: Figure 3: Visual comparison for flow-matching inverse solvers.
Refer to caption: images/compare.png

Figure: Figure 4: Solution path of *Flower* for box inpainting.
Refer to caption: images/batch001_grid.png

##### Ethics Statement.

This work proposes a methodology for solving inverse problems in imaging using pre-trained generative models. Our method is designed as a general-purpose solver and does not target specific sensitive domains. Our experiments are conducted exclusively on publicly available datasets (CelebA and AFHQ-Cat) that are commonly used in the literature. No private or otherwise sensitive data were collected or used. Our method has potential positive applications in areas such as medical imaging, but, as with other generative techniques, should be applied responsibly to avoid misuse in creating misleading content.

##### Acknowledgements

The authors acknowledge support from the European Research Council (ERC Project FunLearn, Grant 101020573) and the Swiss National Science Foundation (Grant 200020_219356). We thank Anne Gagneux and Ségolène Martin for their assistance in reproducing results, and Pakshal Bohra for insightful discussions on sampling schemes.

## References

- Albergo & Vanden-Eijnden (2023)
Michael Samuel Albergo and Eric Vanden-Eijnden.
Building normalizing flows with stochastic interpolants.
In *The Eleventh International Conference on Learning Representations*, 2023.
- Ambrosio et al. (2008)
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.
*Gradient Flows in Metric Spaces and in the Space of Probability Measures*.
Lectures in Mathematics ETH Zürich. Birkhauser, 2nd edition, 2008.
- Beck & Teboulle (2009)
Amir Beck and Marc Teboulle.
A fast iterative shrinkage-thresholding algorithm for linear inverse problems.
*SIAM Journal on Imaging Sciences*, 2(1):183–202, 2009.
- Ben-Hamu et al. (2024)
Heli Ben-Hamu, Omri Puny, Itai Gat, Brian Karrer, Uriel Singer, and Yaron Lipman.
D-Flow: Differentiating through flows for controlled generation.
In *International Conference on Machine Learning (ICML)*, 2024.
- Bubeck (2015)
Sébastien Bubeck.
Convex optimization: Algorithms and complexity.
*Found. Trends Mach. Learn.*, 8(3–4):231–357, November 2015.
ISSN 1935-8237.
doi: 10.1561/2200000050.
- Choi et al. (2021)
Jooyoung Choi, Sungwon Kim, Yonghyun Jeong, Youngjune Gwon, and Sungroh Yoon.
Ilvr: Conditioning method for denoising diffusion probabilistic models.
In *CVF International Conference on Computer Vision (ICCV)*, volume 1, pp.  2, 2021.
- Choi et al. (2020)
Yunjey Choi, Youngjung Uh, Jaejun Yoo, and Jung-Woo Ha.
Stargan v2: Diverse image synthesis for multiple domains.
In *Proceedings of the IEEE/CVF conference on computer vision and pattern recognition*, pp. 8188–8197, 2020.
- Chung et al. (2023)
Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye.
Diffusion posterior sampling for general noisy inverse problems.
In *International Conference on Learning Representations*, 2023.
- Chung et al. (2024)
Hyungjin Chung, Suhyeon Lee, and Jong Chul Ye.
Decomposed diffusion sampler for accelerating large-scale inverse problems.
In *The Twelfth International Conference on Learning Representations*, 2024.
- Ducotterd et al. (2025)
Stanislas Ducotterd, Sebastian Neumayer, and Michael Unser.
Learning of patch-based smooth-plus-sparse models for image reconstruction.
In Beidi Chen, Shijia Liu, Mert Pilanci, Weijie Su, Jeremias Sulam, Yuxiang Wang, and Zhihui Zhu (eds.), *Conference on Parsimony and Learning*, volume 280 of *Proceedings of Machine Learning Research*, pp. 89–104. PMLR, 24–27 Mar 2025.
URL [https://proceedings.mlr.press/v280/ducotterd25a.html](https://proceedings.mlr.press/v280/ducotterd25a.html).
- Figueiredo & Nowak (2003)
M.A.T. Figueiredo and R.D. Nowak.
An em algorithm for wavelet-based image restoration.
*IEEE Transactions on Image Processing*, 12(8):906–916, 2003.
doi: 10.1109/TIP.2003.814255.
- Goujon et al. (2024)
Alexis Goujon, Sebastian Neumayer, and Michael Unser.
Learning weakly convex regularizers for convergent image-reconstruction algorithms.
*SIAM Journal on Imaging Sciences*, 17(1):91–115, 2024.
- Hagemann et al. (2022)
Paul Hagemann, Johannes Hertrich, and Gabriele Steidl.
Stochastic normalizing flows for inverse problems: a Markov Chains viewpoint.
*SIAM Journal on Uncertainty Quantification*, 10(3):1162–1190, 2022.
- Ho et al. (2020)
Jonathan Ho, Ajay Jain, and Pieter Abbeel.
Denoising diffusion probabilistic models.
*arXiv preprint arxiv:2006.11239*, 2020.
- Hurault et al. (2022a)
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis.
Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization.
In *International Conference on Machine Learning*, volume 162, pp. 9483–9505, 2022a.
- Hurault et al. (2022b)
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis.
Gradient step denoiser for convergent plug-and-play.
In *International Conference on Learning Representations*, 2022b.
- Karras et al. (2019)
Tero Karras, Samuli Laine, and Timo Aila.
A style-based generator architecture for generative adversarial networks.
In *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 4401–4410, 2019.
- Kawar et al. (2022)
Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song.
Denoising diffusion restoration models.
In *Advances in Neural Information Processing Systems*, 2022.
- Lipman et al. (2023)
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le.
Flow matching for generative modeling.
In *The Eleventh International Conference on Learning Representations*, 2023.
- Liu (2022)
Qiang Liu.
Rectified flow: A marginal preserving approach to optimal transport.
*arXiv preprint arXiv:2209.14577*, 2022.
- Liu et al. (2023)
Xingchao Liu, Chengyue Gong, and Qiang Liu.
Flow straight and fast: Learning to generate and transfer data with rectified flow.
In *International Conference on Learning Representations (ICLR)*, 2023.
- Mardani et al. (2024)
Morteza Mardani, Jiaming Song, Jan Kautz, and Arash Vahdat.
A variational perspective on solving inverse problems with diffusion models.
In *International Conference on Learning Representations (ICLR)*, 2024.
- Martin et al. (2025)
Ségolène Tiffany Martin, Anne Gagneux, Paul Hagemann, and Gabriele Steidl.
Pnp-flow: Plug-and-play image restoration with flow matching.
In *The Thirteenth International Conference on Learning Representations*, 2025.
- McCann & Unser (2019)
Michael T. McCann and Michael Unser.
Biomedical image reconstruction: From the foundations to deep neural networks.
*Foundations and Trends® in Signal Processing*, 13(3):283–359, 2019.
- Peyré (2025)
Gabriel Peyré.
Optimal Transport for Machine Learners, 2025.
- Pokle et al. (2024)
Ashwini Pokle, Matthew J. Muckley, Ricky T. Q. Chen, and Brian Karrer.
Training-free linear image inverses via flows.
*Transactions on Machine Learning Research*, 2024.
- Pourya et al. (2025)
Mehrsa Pourya, Erich Kobler, Michael Unser, and Sebastian Neumayer.
DEALing with image reconstruction: Deep attentive least squares.
In *Forty-second International Conference on Machine Learning*, 2025.
- Ronneberger et al. (2015)
Olaf Ronneberger, Philipp Fischer, and Thomas Brox.
U-net: Convolutional networks for biomedical image segmentation.
In *Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18*, pp. 234–241. Springer, 2015.
- Roth & Black (2009)
Stefan Roth and Michael J Black.
Fields of experts.
*International Journal of Computer Vision*, 82(2):205–229, 2009.
- Rudin et al. (1992)
Leonid I Rudin, Stanley Osher, and Emad Fatemi.
Nonlinear total variation based noise removal algorithms.
*Physica D: Nonlinear Phenomena*, 60(1-4):259–268, 1992.
- Sohl-Dickstein et al. (2015)
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.
Deep unsupervised learning using nonequilibrium thermodynamics.
In *International conference on machine learning*, pp. 2256–2265. pmlr, 2015.
- Song et al. (2023)
Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz.
Pseudoinverse-guided diffusion models for inverse problems.
In *International Conference on Learning Representations*, 2023.
- Song et al. (2020)
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.
Score-based generative modeling through stochastic differential equations.
In *International Conference on Learning Representations*, 2020.
- Tachella et al. (2023)
Julian Tachella, Dongdong Chen, Samuel Hurault, Matthieu Terris, and Andrew Wang.
DeepInverse: A deep learning framework for inverse problems in imaging, 2023.
- Tong et al. (2024)
Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector-Brooks, Guy Wolf, and Yoshua Bengio.
Improving and generalizing flow-based generative models with minibatch optimal transport.
*Transactions on Machine Learning Research*, 2024.
- Venkatakrishnan et al. (2013)
Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg.
Plug-and-play priors for model based reconstruction.
In *IEEE Global Conference on Signal and Information Processing*, pp. 945–948, 2013.
- Yang et al. (2015)
Shuo Yang, Ping Luo, Chen-Change Loy, and Xiaoou Tang.
From facial parts responses to face detection: A deep learning approach.
In *Proceedings of the IEEE international conference on computer vision*, pp. 3676–3684, 2015.
- Zeng (2001)
Gengsheng Lawrence Zeng.
Image reconstruction—A tutorial.
*Computerized Medical Imaging and Graphics*, 25(2):97–103, 2001.
- Zhang et al. (2025)
Bingliang Zhang, Wenda Chu, Julius Berner, Chenlin Meng, Anima Anandkumar, and Yang Song.
Improving diffusion inverse problem solving with decoupled noise annealing.
In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 20895–20905, June 2025.
- Zhang et al. (2022)
Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte.
Plug-and-play image restoration with deep denoiser prior.
*IEEE Transactions on Pattern Analysis and Machine Intelligence*, 44(10):6360–6376, 2022.
- Zhang et al. (2024)
Yasi Zhang, Peiyu Yu, Yaxuan Zhu, Yingshan Chang, Feng Gao, Ying Nian Wu, and Oscar Leong.
Flow priors for linear inverse problems via iterative corrupted trajectory matching.
In *The Thirty-eighth Annual Conference on Neural Information Processing Systems*, 2024.
- Zhu et al. (2023)
Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool.
Denoising diffusion models for plug-and-play image restoration.
In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pp. 1219–1229, 2023.

## 7 Appendix

### 7.1 Flower Algorithm

We outline the steps of *Flower* in Algorithm [1](https://arxiv.org/html/2509.26287v1#alg1).

Figure: Algorithm 1 FLOWER: Flow Matching Solver for Inverse Problems

### 7.2 Proofs

#### 7.2.1 Proof of Proposition 1

###### Proof.

Lipman et al. (2023) have shown that the velocity vector field that minimizes the conditional flow-matching loss is

$$ ${\bf{v}}^{*}_{t}({\bf{x}})=\mathbb{E}[{\bf{X}}_{1}-{\bf{X}}_{0}|{\bf{X}}_{t}={\bf{x}}_{t}],$ (26) $$

which then yields

$$ $\hat{{\bf{x}}}_{1}({\bf{x}}_{t})=\mathbb{E}[{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}]=\mathbb{E}[{\bf{X}}_{t}+(1-t)({\bf{X}}_{1}-{\bf{X}}_{0})|{\bf{X}}_{t}={\bf{x}}_{t}]={\bf{x}}_{t}+(1-t){\bf{v}}^{\theta}_{t}({\bf{x}}).$ (27) $$

This is the desired result under the assumption that ${\bf{v}}^{\theta}_{t}({\bf{x}})={\bf{v}}^{*}_{t}({\bf{x}})$.
∎

###### Proof.

#### 7.2.2 Proof of Proposition 2

###### Proof.

Using Bayes’ rule, we can write

$$ $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})=\frac{p_{{\bf{Y}}|{\bf{X}}_{1}={\bf{x}}_{1},{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{y}})p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{x}}_{1})}{p_{{\bf{Y}}|{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{y}})}=\frac{p_{{\bf{Y}}|{\bf{X}}_{1}={\bf{x}}_{1}}({\bf{y}})p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{x}}_{1})}{p_{{\bf{Y}}|{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{y}})},$ (28) $$

where we used the conditional independence, given ${\bf{X}}_{1}$, of ${\bf{X}}_{t}$ and the measurement ${\bf{Y}}$. We assumed that $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t}}={\mathcal{N}}({\bm{m}}_{t},\nu_{t}^{2}{\bf{I}})$, and we also have that $p_{{\bf{Y}}|{\bf{X}}_{1}={\bf{x}}_{1}}={\mathcal{N}}({\bf{H}}{\bf{x}}_{1},\sigma_{n}^{2}{\bf{I}})$ by construction. Taking the logarithm of ([28](https://arxiv.org/html/2509.26287v1#S7.E28)), we obtain

$$ $\displaystyle-2\ln\left(p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})\right)$ $\displaystyle=({\bf{y}}-{\bf{H}}{\bf{x}}_{1})^{\top}\sigma_{n}^{-2}({\bf{y}}-{\bf{H}}{\bf{x}}_{1})+({\bf{x}}_{1}-{\bm{m}}_{t})^{\top}\nu_{t}^{-2}({\bf{x}}_{1}-{\bm{m}}_{t})+C$ (29) $\displaystyle=-2{\bf{x}}_{1}^{\top}\left(\nu_{t}^{-2}{\bm{m}}_{t}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{y}}\right)+{\bf{x}}_{1}^{\top}\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right){\bf{x}}_{1}+C^{\prime}$ (30) $\displaystyle=-2{\bf{x}}_{1}^{\top}\left(\nu_{t}^{-2}{\bm{m}}_{t}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{y}}\right)+{\bf{x}}_{1}^{\top}{\bf{\Sigma}}_{t}^{-1}{\bf{x}}_{1}+C^{\prime},$ (31) $$

where $C,C^{\prime}$ are independent of ${\bf{x}}_{1}$ and considered constants, and where we defined ${\bf{\Sigma}}_{t}=\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right)^{-1}$, which is well-defined because $\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}$ is positive-definite. Completing the square with a term independent of ${\bf{x}}_{1}$, we get

$$ $-2\ln\left(p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})\right)=({\bf{x}}_{1}-{\bm{\mu}}_{t})^{\top}{\bf{\Sigma}}_{t}^{-1}({\bf{x}}_{1}-{\bm{\mu}}_{t})+C^{\prime\prime},$ (32) $$

where ${\bm{\mu}}_{t}={\bf{\Sigma}}_{t}\left(\nu_{t}^{-2}{\bm{m}}_{t}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{y}}\right)$ and $C^{\prime\prime}$ is again a constant independent of ${\bf{x}}_{1}$. This yields

$$ $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})=e^{-C^{\prime\prime}/2}\exp\left(-\frac{1}{2}({\bf{x}}_{1}-{\bm{\mu}}_{t})^{\top}{\bf{\Sigma}}_{t}^{-1}({\bf{x}}_{1}-{\bm{\mu}}_{t})\right).$ (33) $$

As $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}$ is a probability density function, $e^{-C^{\prime\prime}/2}$ corresponds to its normalization factor, which therefore proves that $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})={\mathcal{N}}({\bf{x}}_{1};{\bm{\mu}}_{t},{\bf{\Sigma}}_{t})$.
∎

###### Proof.

#### 7.2.3 Proof of Proposition 3

###### Proof.

By using the marginal distributions, the general chain rule for joint probability, and independence, we obtain that

$$ $\displaystyle p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})$ $\displaystyle=\int_{\mathbb{R}^{d}}p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})p_{{\bf{X}}_{t+\Delta t}|{\bf{X}}_{1}={\bf{x}}_{1},{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})\mathrm{d}{\bf{x}}_{1}$ (34) $\displaystyle=\int_{\mathbb{R}^{d}}p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})p_{{\bf{X}}_{t+\Delta t}|{\bf{X}}_{1}={\bf{x}}_{1},{\bf{X}}_{t}={\bf{x}}_{t}}({\bf{x}}_{t+\Delta t})\mathrm{d}{\bf{x}}_{1}$ (35) $\displaystyle=\int_{\mathbb{R}^{d}}p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1}){\mathcal{N}}\left({\bf{x}}_{t+\Delta t};(t+\Delta t){\bf{x}}_{1},(1-t-\Delta t)^{2}{\bf{I}}\right)\mathrm{d}{\bf{x}}_{1},$ (36) $$

where we used the fact that, given ${\bf{X}}_{1}={\bf{x}}_{1}$, ${\bf{X}}_{t+\Delta t}=(1-t-\Delta t){\bf{X}}_{0}+(t+\Delta t){\bf{x}}_{1}\sim{\mathcal{N}}\left((t+\Delta t){\bf{x}}_{1},(1-t-\Delta t)^{2}{\bf{I}}\right)$. By inserting the expression of $p_{{\bf{X}}_{1}|{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{1})$, we obtain that

$$ $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})=\int_{\mathbb{R}^{d}}\mathcal{N}\left({\bf{x}}_{1};{\bm{\mu}}_{t},{\bf{\Sigma}}_{t}\right)\mathcal{N}\left({\bf{x}}_{t+\Delta t};\left(t+\Delta t\right){\bf{x}}_{1},\left(1-t-\Delta t\right)^{2}{\bf{I}}\right)\mathrm{d}{\bf{x}}_{1}.$ (37) $$

This integral can be rewritten as a convolution of two Gaussian distributions, which also yields a Gaussian distribution. Explicitly, we have that

$$ $\displaystyle p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})$ $\displaystyle=\int_{\mathbb{R}^{d}}{\mathcal{N}}\left({\bf{x}}_{1};{\bm{\mu}}_{t},{\bf{\Sigma}}_{t}\right){\mathcal{N}}\left({\bf{x}}_{t+\Delta t}-\left(t+\Delta t\right){\bf{x}}_{1};{\bf{0}},\left(1-t-\Delta t\right)^{2}{\bf{I}}\right)\mathrm{d}{\bf{x}}_{1}$ (38) $\displaystyle=\int_{\mathbb{R}^{d}}{\mathcal{N}}\left(\frac{{\bf{z}}}{t+\Delta t};{\bm{\mu}}_{t},{\bf{\Sigma}}_{t}\right){\mathcal{N}}\left({\bf{x}}_{t+\Delta t}-{\bf{z}};{\bf{0}},\left(1-t-\Delta t\right)^{2}{\bf{I}}\right)\frac{\mathrm{d}{\bf{z}}}{(t+\Delta t)^{d}}$ (39) $\displaystyle=\int_{\mathbb{R}^{d}}{\mathcal{N}}\left({\bf{z}};{\bm{(}}t+\Delta t){\bm{\mu}}_{t},{\bf{(}}t+\Delta t)^{2}{\bf{\Sigma}}_{t}\right){\mathcal{N}}\left({\bf{x}}_{t+\Delta t}-{\bf{z}};{\bf{0}},\left(1-t-\Delta t\right)^{2}{\bf{I}}\right)\mathrm{d}{\bf{z}}.$ (40) $$

Using the Gaussian convolution identity

$$ $\int_{\mathbb{R}^{d}}{\mathcal{N}}({\bf{z}};{\bm{\mu}}_{1},{\bf{\Sigma}}_{1}){\mathcal{N}}({\bf{x}}-{\bf{z}};{\bm{\mu}}_{2},{\bf{\Sigma}}_{2})\,\mathrm{d}{\bf{z}}={\mathcal{N}}({\bf{x}};{\bm{\mu}}_{1}+{\bm{\mu}}_{2},{\bf{\Sigma}}_{1}+{\bf{\Sigma}}_{2}),$ (41) $$

the result simplifies to

$$ $p_{{\bf{X}}_{t+\Delta t}\mid{\bf{X}}_{t}={\bf{x}}_{t},{\bf{Y}}={\bf{y}}}({\bf{x}}_{t+\Delta t})=\mathcal{N}\left({\bf{x}}_{t+\Delta t};(t+\Delta t){\bm{\mu}}_{t},(t+\Delta t)^{2}{\bf{\Sigma}}_{t}+(1-t-\Delta t)^{2}{\bf{I}}\right),$ (42) $$

which completes the proof.
∎

###### Proof.

#### 7.2.4 Sampling from the Non-Isotropic Gaussian

We want to show that

$$ ${\bm{\kappa}}_{t}={\bf{\Sigma}}_{t}\left(\nu_{t}^{-1}{\bm{\epsilon}}_{1}+\sigma_{n}^{-1}{\bf{H}}^{\top}{\bm{\epsilon}}_{2}\right)$ (43) $$

has distribution ${\mathcal{N}}(0,{\bf{\Sigma}}_{t})$, where ${\bm{\epsilon}}_{1}\sim{\mathcal{N}}(0,{\bf{I}}_{d})$ and ${\bm{\epsilon}}_{2}\sim{\mathcal{N}}(0,{\bf{I}}_{M})$ are independent, and

$$ ${\bf{\Sigma}}_{t}=\big(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\big)^{-1}.$ (44) $$

Observe that ${\bm{\kappa}}_{t}$ is a Gaussian random vector because it is a linear transform of the independent Gaussians ${\bm{\epsilon}}_{1}$ and ${\bm{\epsilon}}_{2}$. Moreover, since ${\bm{\epsilon}}_{1}$ and ${\bm{\epsilon}}_{2}$ are zero-mean, we have that $\mathbb{E}[{\bm{\kappa}}_{t}]={\bf{0}}$. Next, we compute the covariance matrix of ${\bm{\kappa}}_{t}$ as

$$ $\mathrm{Cov}({\bm{\kappa}}_{t})={\bf{\Sigma}}_{t}\,\mathrm{Cov}\!\big(\nu_{t}^{-1}{\bm{\epsilon}}_{1}+\sigma_{n}^{-1}{\bf{H}}^{\top}{\bm{\epsilon}}_{2}\big)\,{\bf{\Sigma}}_{t}.$ (45) $$

By independence of ${\bm{\epsilon}}_{1}$ and ${\bm{\epsilon}}_{2}$, we obtain that

$$ $\mathrm{Cov}\!\big(\nu_{t}^{-1}{\bm{\epsilon}}_{1}+\sigma_{n}^{-1}{\bf{H}}^{\top}{\bm{\epsilon}}_{2}\big)=\nu_{t}^{-2}{\bf{I}}_{d}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}},$ (46) $$

which implies that

$$ $\mathrm{Cov}({\bm{\kappa}}_{t})={\bf{\Sigma}}_{t}\left(\nu_{t}^{-2}{\bf{I}}_{d}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right){\bf{\Sigma}}_{t}.$ (47) $$

But, by definition, ${\bf{\Sigma}}_{t}=\big(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\big)^{-1}$, which leads to the desired result, as

$$ $\mathrm{Cov}({\bm{\kappa}}_{t})={\bf{\Sigma}}_{t}\left(\nu_{t}^{-2}{\bf{I}}+\sigma_{n}^{-2}{\bf{H}}^{\top}{\bf{H}}\right){\bf{\Sigma}}_{t}={\bf{\Sigma}}_{t}.$ (48) $$

### 7.3 Numerical Results Extension

#### 7.3.1 Toy Experiment

##### Target Prior.

As explained in Section [5.1](https://arxiv.org/html/2509.26287v1#S5.SS1), the target distribution of our data ${\bf{X}}\in\mathbb{R}^{2}$ is a GMM with uniform mixtures given explicitly by

$$ $p_{{\bf{X}}}=\frac{1}{K}\sum_{k=1}^{K}{\mathcal{N}}({\bm{\mu}}_{k},{\bf{\Sigma}})$ (49) $$

for some $K\in\mathbb{N}$, ${\bm{\mu}}_{k}\in\mathbb{R}^{2}$, and ${\bf{\Sigma}}\in\mathbb{S}^{2}_{++}$. Specifically, we use $K=3$ with ${\bm{\mu}}_{1}=(-0.25,-0.25),{\bm{\mu}}_{2}=(-0.25,0.25),{\bm{\mu}}_{3}=(0.25,-0.25)$, and covariance matrix ${\bf{\Sigma}}=0.25^{2}{\bf{I}}_{2}$.

##### Target Posterior.

The advantage of this setup is that the posterior distribution $p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}$ can be computed exactly using Bayes’ rule. It is given by

$$ $p_{{\bf{X}}|{\bf{Y}}={\bf{y}}}=\sum_{k=1}^{K}w_{k}{\mathcal{N}}({\bm{\mu}}_{k,\text{post}},{\bf{\Sigma}}_{\text{post}})$ (50) $$

where, for all $k=1,\ldots,K$, $w_{k}\geq 0$ are some weights and

$$ $\displaystyle{\bm{\mu}}_{k,\text{post}}$ $\displaystyle=\left({\bf{\Sigma}}^{-1}+\sigma_{n}^{-2}{\bf{h}}{\bf{h}}^{\top}\right)^{-1}\left(\sigma_{n}^{-2}{\bf{h}}y+{\bf{\Sigma}}^{-1}{\bm{\mu}}_{k}\right),$ (51) $\displaystyle{\bf{\Sigma}}_{\text{post}}$ $\displaystyle=\left({\bf{\Sigma}}^{-1}+\sigma_{n}^{-2}{\bf{h}}{\bf{h}}^{\top}\right)^{-1}.$ (52) $$

##### Training Details.

The underlying unconditional velocity network is a fully connected network that takes as input a 2D vector and a scalar time, concatenated into a 3D input. It consists of two hidden layers of size 256 with SiLU activations, followed by a final linear layer that outputs a 2D vector. For training, we use a batch size of $2048$ with $20000$ training steps and a learning rate of $10^{-3}$.

##### More Results.

In Figure [7](https://arxiv.org/html/2509.26287v1#S7.F7), we present another example of *Flower* posterior sampling with ${\bf{h}}^{\top}=[1.5,-1.5]$, $\sigma_{n}=0.75$, and the observation $y=1$. Once again, we observe that *Flower* successfully generates samples that closely match the true posterior for $\gamma=1$. In contrast, for $\gamma=0$, samples from the tails of the true posterior are missing. In Figures [7](https://arxiv.org/html/2509.26287v1#S7.F7) and [7](https://arxiv.org/html/2509.26287v1#S7.F7), the solution path is illustrated across the successive steps of *Flower* for $\gamma=0$ and $\gamma=1$, respectively, which allows us to visualize the dynamics of each step directly.

Figure: Figure 5: Temporal evolution of 2D *Flower* and comparison with true posterior for noise variance $\sigma_{n}=0.75$.
Refer to caption: images/flower_2d_in_time_a.png

### 7.4 Benchmark Experiments

#### 7.4.1 Effect of Source-Target Coupling

Our goal in this section is to benchmark the effect of source–target coupling in the training of the underlying (unconditional) velocity network. To this end, we consider two variants of coupling: one based on mini-batch optimal transport (OT) coupling and the other on independent (IND) coupling. For each variant, we report two settings: a single evaluation of *Flower*; and average over five evaluations. This results in four cases, summarized in Table [3](https://arxiv.org/html/2509.26287v1#S7.T3), where we provide results for all inverse problems discussed in the main text, evaluated on 100 test images from the CelebA dataset. For this table, we fix $\gamma=0$ and $N=100$.

**Table 3: Effect of source-target coupling of the underlying velocity network of *Flower* on different inverse problems on 100 test images of the dataset CelebA.**
| Method | Denoising | Deblurring | Super-resolution | Random inpainting | Box inpainting |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS |  |
| Degraded | 20.00 | 0.348 | 0.372 | 27.83 | 0.740 | 0.126 | 10.26 | 0.183 | 0.827 | 11.95 | 0.196 | 1.041 | 22.27 | 0.742 | 0.214 |
| Flower1-OT (ours) | 32.28 | 0.914 | 0.034 | 34.98 | 0.947 | 0.026 | 32.36 | 0.923 | 0.034 | 33.08 | 0.944 | 0.018 | 31.19 | 0.945 | 0.022 |
| Flower5-OT (ours) | 33.14 | 0.926 | 0.038 | 35.67 | 0.954 | 0.032 | 33.09 | 0.932 | 0.040 | 33.95 | 0.953 | 0.020 | 31.87 | 0.952 | 0.023 |
| Flower1-IND (ours) | 32.60 | 0.918 | 0.032 | 35.22 | 0.950 | 0.026 | 32.65 | 0.927 | 0.034 | 33.23 | 0.947 | 0.017 | 31.90 | 0.950 | 0.021 |
| Flower5-IND (ours) | 33.48 | 0.930 | 0.037 | 35.90 | 0.957 | 0.031 | 33.41 | 0.935 | 0.039 | 34.24 | 0.955 | 0.020 | 32.78 | 0.958 | 0.022 |

We observe that independent coupling improves the results, which is consistent with our theoretical requirements. Nevertheless, the OT-based variant remains highly competitive, as also illustrated in the main paper. For visual comparison, in Figure [8](https://arxiv.org/html/2509.26287v1#S7.F8), we show an example from the deblurring task on an image from the CelebA dataset.

#### 7.4.2 Effect of γ \gamma

The hyperparameter $\gamma\in\{0,1\}$ in *Flower* controls whether the uncertainty of the refinement step (Step 2) is taken into account. While $\gamma=1$ is required for a Bayesian interpretation of our method, in practice we find that $\gamma=0$ yields more favorable image-reconstruction metrics. This observation is consistent with the toy experiments, where $\gamma=0$ led *Flower* to generate samples concentrated in higher-probability regions, while failing to capture the tails of the posterior.

To illustrate this effect, we provide a visual example of the deblurring task on a CelebA image, using the velocity network trained with mini-batch OT coupling for consistency with the numerical results reported in this paper. A single reconstruction with $\gamma=1$ achieves a PSNR of only $30.84$, compared to $33.01$ for a single reconstruction with $\gamma=0$ (Figure [8](https://arxiv.org/html/2509.26287v1#S7.F8)). Moreover, with $\gamma=1$, it is necessary to average over 100 reconstructions to reach a PSNR comparable to that of $\gamma=0$, where we average only over 5 reconstructions (Figure [8](https://arxiv.org/html/2509.26287v1#S7.F8)).

To further illustrate this behavior, we show in Figures [7](https://arxiv.org/html/2509.26287v1#S7.F7) and [7](https://arxiv.org/html/2509.26287v1#S7.F7) the solution paths of the three *Flower* steps over time for $\gamma=0$ and $\gamma=1$, respectively. We observe that $\gamma=1$ leads to a noisier refinement step. Consequently, we adopt $\gamma=0$ for all inverse problems, as this provides a more stable refinement step and consistently yields a better quality of reconstruction despite fewer averages.

#### 7.4.3 Computational Efficiency

We report in Table [4](https://arxiv.org/html/2509.26287v1#S7.T4) the computational time and memory usage for several methods. Each entry corresponds to the average runtime for the deblurring inverse problem, averaged over 10 CelebA test images of size $128\times 128$. All experiments were conducted on a Tesla V100-SXM2-32GB GPU.

**Table 4: Computation times and memory usage for various methods.**
| Method | OT-ODE | D-Flow | Flow Priors | PnP-Flow1 | Flower1 |
| --- | --- | --- | --- | --- | --- |
| Time (s) | 6.549 | 142.18 | 63.771 | 3.020 | 5.622 |
| Memory (GB) | 1.183 | 11.125 | 3.807 | 0.216 | 0.217 |

#### 7.4.4 Hyperparameters for All Methods

In Tables [5](https://arxiv.org/html/2509.26287v1#S7.T5) and [6](https://arxiv.org/html/2509.26287v1#S7.T6), we report the hyperparameters that we used for all methods. Most of the hyperparameters are adapted from Martin et al. (2025).

**Table 5: Hyperparameters for all methods on the CelebA dataset.**
| Method | Hyperparameters | Denoising | Deblurring | Super-resolution | Random inpainting | Box inpainting |
| --- | --- | --- | --- | --- | --- | --- |
| DiffPIR | $\zeta$ (blending) | 1.0 | 1.0 | 1.0 | 1.0 | N/A |
|  | $\lambda$ (regularization) | 1.0 | 1000.0 | 100.0 | 1.0 | N/A |
| PnP-GS | $\gamma$ (learning rate) | - | 2.0 | 2.0 | 1.0 | N/A |
|  | $\alpha$ (inertia param.) | 1.0 | 0.5 | 1.0 | 0.5 | N/A |
|  | $\sigma_{f}$ (factor for noise input) | 1.0 | 1.8 | 3.0 | 1.0 | N/A |
|  | $n_{\text{iter}}$ (number of iter.) | 1 | 35 | 20 | 23 | N/A |
| OT-ODE | $t_{0}$ (initial time) | 0.3 | 0.4 | 0.1 | 0.1 | 0.1 |
|  | $\gamma$ | time-dependent | time-dependent | constant | constant | time-dependent |
| Flow-Priors | $\lambda$ (regularization) | 100 | 1,000 | 10,000 | 10,000 | 10,000 |
|  | $\eta$ (learning rate) | 0.01 | 0.01 | 0.1 | 0.01 | 0.01 |
| D-Flow | $\lambda$ (regularization) | 0.001 | 0.001 | 0.001 | 0.01 | 0.001 |
|  | $\alpha$ (blending) | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 |
|  | $n_{\text{iter}}$ (number of iter.) | 3 | 7 | 10 | 20 | 9 |
| PnP-Flow1 | $\alpha$ (learning-rate factor) | 0.8 | 0.01 | 0.3 | 0.01 | 0.5 |
|  | $N$ (Number of time steps) | 100 | 100 | 100 | 100 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 1 | 1 | 1 | 1 | 1 |
| PnP-Flow5 | $\alpha$ (learning-rate factor) | 0.8 | 0.01 | 0.3 | 0.01 | 0.5 |
|  | $N$ (Number of time steps) | 100 | 100 | 100 | 100 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 5 | 5 | 5 | 5 | 5 |
| Flower1-OT | $\gamma$ (refinement uncertainty) | 0 | 0 | 0 | 0 | 0 |
|  | $N$ (Number of time steps) | 100 | 100 | 100 | 100 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 1 | 1 | 1 | 1 | 1 |
| Flower5-OT | $\gamma$ (refinement uncertainty) | 0 | 0 | 0 | 0 | 0 |
|  | $N$ (Number of time steps) | 100 | 100 | 100 | 100 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 5 | 5 | 5 | 5 | 5 |

**Table 6: Hyperparameters for all methods on the AFHQ-Cat dataset.**
| Method |  | Denoising | Deblurring | Super-resolution | Random inpainting | Box inpainting |
| --- | --- | --- | --- | --- | --- | --- |
| DiffPIR | $\zeta$ (blending) | 1.0 | 1.0 | 1.0 | 1.0 | N/A |
|  | $\lambda$ (regularization) | 1.0 | 1000.0 | 100.0 | 1.0 | N/A |
| PnP-GS | $\gamma$ (learning rate) | - | 2.0 | 2.0 | 1.0 | N/A |
|  | $\alpha$ (inertia param.) | 1.0 | 0.3 | 1.0 | 0.5 | N/A |
|  | $\sigma_{f}$ (factor for noise input) | 1.0 | 1.8 | 5.0 | 1.0 | N/A |
|  | $n_{\text{iter}}$ (number of iter.) | 1 | 60 | 50 | 23 | N/A |
| OT-ODE | $t_{0}$ (initial time) | 0.3 | 0.3 | 0.1 | 0.1 | 0.1 |
|  | $\gamma$ | time-dependent | time-dependent | constant | constant | time-dependent |
| Flow-Priors | $\lambda$ (regularization) | 100 | 1,000 | 10,000 | 10,000 | 10,000 |
|  | $\eta$ (learning rate) | 0.01 | 0.01 | 0.1 | 0.01 | 0.01 |
| D-Flow | $\lambda$ (regularization) | 0.001 | 0.01 | 0.001 | 0.001 | 0.01 |
|  | $\alpha$ (blending) | 0.1 | 0.5 | 0.1 | 0.1 | 0.1 |
|  | $n_{\text{iter}}$ (number of iter.) | 3 | 20 | 20 | 20 | 9 |
| PnP-Flow1 | $\alpha$ (learning-rate factor) | 0.8 | 0.01 | 0.01 | 0.01 | 0.5 |
|  | $N$ (Number of time steps) | 100 | 500 | 500 | 200 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 1 | 1 | 1 | 1 | 1 |
| PnP-Flow5 | $\alpha$ (learning-rate factor) | 0.8 | 0.01 | 0.01 | 0.01 | 0.5 |
|  | $N$ (Number of time steps) | 100 | 500 | 500 | 200 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 5 | 5 | 5 | 5 | 5 |
| Flower1-OT | $\gamma$ (refinement uncertainty) | 0 | 0 | 0 | 0 | 0 |
|  | $N$ (Number of time steps) | 100 | 100 | 100 | 100 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 1 | 1 | 1 | 1 | 1 |
| Flower5-OT | $\gamma$ (refinement uncertainty) | 0 | 0 | 0 | 0 | 0 |
|  | $N$ (Number of time steps) | 100 | 100 | 500 | 200 | 100 |
|  | $N_{\mathrm{Avg}}$ (Number of averagings) | 5 | 5 | 5 | 5 | 5 |

Figure: Figure 8: Effect of the source-target coupling for the underlying flow ($\gamma=0$).
Refer to caption: images/coupling_effect.png

Figure: Figure 9: Deblurring results with $\gamma=1$, obtained by averaging over $1$, $5$, $50$, $100$, and $200$ runs of *Flower*.
Refer to caption: images/gamma_1_ot.png

Figure: Figure 10: Solution path of *Flower* for deblurring with $\gamma=0$.
Refer to caption: images/batch002_grid.png

Figure: Figure 11: Solution path of *Flower* for deblurring with $\gamma=1$.
Refer to caption: images/batch002_grid_cov.png

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 4
presentation: 2
questions: Please refer to the weaknesses section for questions if any. Thanks.
soundness: 2
strengths: + The paper proposes a new inverse problem solver that unifies PnP intuition with a Bayesian posterior-sampling narrative.
+ Although I didn't check all the math details presented in the paper, it seems to provide a decent theoretical connection to bayesian posterior sampling using \piGDM-style gaussian approximation.
+ The practicality of this method seems good, with runtime comparable to PnP-Flow and OT-ODE; memory light (tab. 4).
summary: The paper proposes a solver for inverse problems by leveraging pre-trained flow models to produce reconstructions of observed measurements. The proposed solver use pretrained network to predict the target, which is then refined by a proximal step and noise-injection step. The paper provides theoretical analysis showing that the process aproximates bayesian posterior sampling, and provides empirical results on CelebA and AFHQ-cat datasets.
weaknesses: - The empirical contribution of this submission, though, seems relatively weak. The main numerical experiments are conducted on two small scale datasets, 100 samples from CelebA and AFHQ-Cat, respectively.
- The overall improvements upon previous methods seem not very significant (in tab. 1 & 2).
- In terms of the limitation of assumptions, i) the paper focuses only on linear problems, and ii) if I understand correctly the proximal step further assumes isotropic prior and measurement operation, which could be oversimplified.
- The most stand-out distinctions from previous flow-based inverse problem solvers are not very clear to me.
- Admittely, I am not very familiar with the current progress of this particular sub-domain. I would therefore also like to hear other colleague reviewers' opinions.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Estimating $\Sigma_t$ involves computing an inverse of a matrix. How do you efficiently solve this such that the solver itself is light weight?
2. Line 342 states that the validation set for AFHQ-Cat was constructed with 32 images. However, the caption of Table 2 suggests that 100 test images were used. One of these statements need to be updated. 
3. What is $m_t$ in line 715? Is it $\hat{x}_1(x_t)$? If that is indeed the case, it’s more readable to use $\hat{x}_1(x_t)$ instead.
soundness: 3
strengths: 1. The core method has been explained intuitively and seems easy to implement. 
2. Computational efficiency: Runtime of Flower is less compared to many other baselines as indicated in Table 4. 
3. Robustness: The method uses same hyper parameters for all the linear inverse problems included in the paper. This simplifies hyper parameter tuning.
summary: Flower is a plug-and-play solver for linear inverse problems with pre-trained flow matching models. Given a pre-trained velocity model, Flower solves a linear inverse problem by iterating over three steps: 
1. Estimation of clean image $x_1$ given the current noisy image and velocity model. 
2. Projection of the estimated image $x_1$ into a feasible set with a proximal operator 
3. Estimation of the next noisy image along the flow using the projected $x_1$.  

The correctness of the proposed solver has been explained through Bayesian analysis and the key result shows that Flower samples from an approximate posterior distribution.
weaknesses: 1. Writing should be more specific. 
    1. Contributions should indicate that this is a solver for linear inverse problems. Statements such as “Thus, Flower can be seen as a plug-and-play inverse solver with general application, with our Bayesian justification establishing a novel link between posterior sampling and the PnP frameworks in generative models.” need to be updated to clarify that this is a solver for linear inverse problems. 
    2. The paper should also indicate that the Bayesian analysis is only valid for sampling from the approximate posterior and not the true posterior. 
2. Only 100 images have been used for quantitative results which is a really small test dataset.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: - Could the authors provide an ablation plot showing the effect of varying the number of evaluations per step (e.g., Flower1-10?) to better illustrate the trade-off between reconstruction quality and computational cost?

- The current method is compatible with the standard Euler sampler used in flow matching. Could the authors discuss whether other ODE or higher-order samplers (e.g., Heun or adaptive solvers) could be integrated into Flower, and what potential advantages or challenges might arise?

- Could the authors elaborate on possible extensions of Flower to nonlinear or real-world inverse problems (e.g., MRI reconstruction, phase retrieval, or compressed sensing)? In particular, what modifications or additional components would be required to handle these more complex forward models?
soundness: 4
strengths: The paper:
-  Introduces a novel three-step flow-matching framework (Flower) for inverse problems, bridging generative flow models and plug-and-play (PnP) methods under a unified Bayesian view.

-  Provides solid theoretical grounding with a clear Bayesian justification and strong empirical results showing state-of-the-art/competitive reconstructions across diverse tasks.

- The paper is well-structured and accessible, clearly explaining algorithmic steps and visualizing the process through intuitive figures.
summary: The paper “FLOWER: A Flow-Matching Solver for Inverse Problems” introduces Flower, a general-purpose solver that leverages a pre-trained flow model to reconstruct signals consistent with given measurements. Flower operates through a three-step iterative process: (1) flow-consistent destination estimation, where the velocity network predicts a denoised target; (2) measurement-aware refinement, projecting the estimate onto the feasible set defined by the forward operator; and (3) time progression, updating the trajectory along the flow path. The authors provide a Bayesian interpretation, showing that Flower approximates posterior sampling, thereby linking plug-and-play (PnP) and generative flow-based solvers within a unified framework. Empirically, Flower achieves state-of-the-art reconstruction quality across diverse inverse problems using nearly identical hyperparameters. Its main contribution lies in bridging Bayesian posterior sampling and PnP methods through a theoretically grounded and practically efficient flow-matching formulation.
weaknesses: -  **More ablation** could be conducted on the number of evaluations performed at each step. The paper currently presents results for Flower1 and Flower5, where Flower5 shows steady improvement over Flower1. An ablation plot illustrating this trend would strengthen the analysis.

-  Although the method claims hyperparameter robustness, the validation is limited to standard imaging tasks. Demonstrating applications to **nonlinear** or **real-world** inverse problems (e.g., MRI, or compressed sensing) would further validate its generality and impact.


Decision: Accept (Poster)