Constrained Posterior Sampling: Time Series Generation with Hard Constraints

## Abstract

Abstract Generating realistic time series samples is crucial for stress-testing models and protecting user privacy by using synthetic data. In engineering and safety-critical applications, these samples must meet certain hard constraints that are domain-specific or naturally imposed by physics or nature. Consider, for example, generating electricity demand patterns with constraints on peak demand times. This can be used to stress-test the functioning of power grids during adverse weather conditions. Existing approaches for generating constrained time series are either not scalable or degrade sample quality. To address these challenges, we introduce Constrained Posterior Sampling (CPS), a diffusion-based sampling algorithm that aims to project the posterior mean estimate into the constraint set after each denoising update. Notably, CPS scales to a large number of constraints ( ∼ 100 similar-to absent 100 \sim 100 ) without requiring additional training. We provide theoretical justifications highlighting the impact of our projection step on sampling. Empirically, CPS outperforms state-of-the-art methods in sample quality and similarity to real time series by around 10% and 42%, respectively, on real-world stocks, traffic, and air quality datasets.

## 1 Introduction

Synthesizing realistic time series samples can aid in “what-if” scenario analysis, stress-testing machine learning (ML) models (Rizzato et al., ; Gowal et al., ), anonymizing private user data (Yoon et al., ), etc. Current approaches for time series generation use state-of-the-art (SOTA) generative models, such as Generative Adversarial Networks (GANs) (Yoon et al., ; Donahue et al., ) and Diffusion Models (DMs) (Tashiro et al., ; Alcaraz & Strodthoff, ; Narasimhan et al., ), to generate high-fidelity time series samples.

Figure: Figure 1: Our Proposed Constrained Posterior Sampling (CPS) Approach. CPS is a novel diffusion-based sampling approach to generate time series samples that adhere to hard constraints. Here, we show an example of generating the daily stock price time series, where CPS ensures that the generated stock prices adhere to natural constraints such as the bounds on the opening and closing prices of the stock.
Refer to caption: /html/2410.12652/assets/x1.png

However, generating realistic and high-fidelity time series samples requires strict adherence to various domain-specific constraints. For example, consider generating the daily Open-high-low-close (OHLC) chart for the stock price of an S&P 500 company. The generated time series samples should have opening and closing stock prices bounded by the high and low values. Similarly, consider generating stock price time series with a user-specified measure of volatility to stress-test trading strategies. If the generated samples do not have the exact volatility, the stress testing results might not be accurate.

On a more general note, the advent of large-scale generative models for language and vision, like GPT-4 (Bubeck et al., ) and Stable Diffusion (Podell et al., ), has increased the focus on constraining the outputs from these models, owing to usefulness and privacy reasons. Note that we cannot clearly define the notion of a constraint set in these domains. For example, verifying if the image of a hand has 6 fingers is practically hard, as all deep-learned perception models for this task have associated prediction errors. However, our key insight is that we can describe a time series through statistical features computed using well-defined functions. These features can be imposed as constraints, and we can accurately verify the constraint satisfaction. Hence, the time series domain allows for the development of a new class of constrained generation algorithms. We first outline the qualities of an ideal constrained time series generator.

- 1.
Training-free approach to include multiple constraints: Training the generative model for a specific constraint, as in the case of Loss-DiffTime (Coletta et al., ), is not scalable. A model trained to generate samples with specified mean constraints cannot adapt to argmax constraints.
- 2.
Independence from external realism enforcers: Generally, prior works involve a projection step to a feasible set defined by a set of constraints, which often destroys the sample quality. To address this, prior approaches (Coletta et al., ) rely on external models to enforce realism, in addition to the generative model, resulting in additional training and complex sampling procedures.
- 3.
Hyperparameter-free approach to constrained generation: The choice of guidance weights in guidance-based approaches with DMs significantly affects the sample quality. Optimizing for guidance weights becomes combinatorially hard while dealing with hundreds of constraints.

Figure: Figure 2: CPS outperforms existing approaches on real-world datasets. Dynamic Time Warping (DTW) measures the similarity between the real and the generated time series. The Train on Synthetic and Test on Real (TSTR) evaluates a task model on real test data when the model was trained on synthetic data. Improved TSTR indicates high generated sample quality. CPS provides 42% and 10% improvements for DTW and TSTR, respectively, over SOTA methods.
Refer to caption: /html/2410.12652/assets/x2.png

Given the following requirements, we propose Constrained Posterior Sampling  (CPS), a novel sampling procedure for diffusion-based generative models (check Fig. [1](#S1.F1)). CPS introduces a projection step that aims to project the posterior mean estimate into the constraint set after each diffusion denoising update. We rely on off-the-shelf optimization routines, thereby providing a training and hyperparameter-free approach to include multiple constraints. Additionally, CPS does not require external models to enforce realism, as the key intuition in our approach is that the subsequent denoising steps rectify the adverse effects of the projection steps toward sample quality. To this end, our contributions in this paper are:

- 1.
We present Constrained Posterior Sampling  ((CPS), Fig. [1](#S1.F1)), a scalable diffusion sampling process that generates realistic time series samples that belong to a constraint set. Without any additional training, CPS can handle a large number of constraints without sacrificing sample quality (Fig. [3](#S2.F3)).
- 2.
We provide a detailed theoretical analysis of the effect of modifying the traditional diffusion sampling process with CPS. Additionally, we perform convergence analysis for well-studied settings, such as convex constraint sets and Gaussian prior data distribution, to draw useful insights for the practical implementation of CPS.
- 3.
Through extensive experiments on six diverse real-world and simulated datasets spanning finance, traffic, and environmental monitoring, we demonstrate that CPS outperforms SOTA approaches on sample quality, similarity, and constraint violation metrics (check Fig. [2](#S1.F2)).

## 2 Preliminaries

Notations: We denote a time series sample by $x\in\mathbb{R}^{K\times L}$. Here, $K$ and $L$ refer to the number of channels and the horizon, respectively. A dataset is defined as $\mathcal{D}=\{x^{1},\dots,x^{N_{\mathrm{D}}}\}$, where the superscript $i\in[1,\dots N_{\mathrm{D}}]$ refers to the sample number, and $N_{\mathrm{D}}$ is the total number of samples in the dataset. $P_{\mathrm{data}}$ denotes the real time series data distribution. $x^{i}$ is the realization of the random vector $X^{i}$, where $X^{1},\dots X^{N_{\mathrm{D}}}\sim P_{\mathrm{data}}$. The Probability Density Function (PDF) associated with $P_{\mathrm{data}}$ is represented by $p_{\mathrm{data}}:\mathbb{R}^{K\times L}\rightarrow\mathbb{R}$, where $\int p_{\mathrm{data}}(x)dx=1$. Here, $\int$ refers to the integration operator over $\mathbb{R}^{K\times L}$. The notation $\mathcal{N}(\mu,\Sigma)$ refers to the Gaussian distribution with mean $\mu$ and covariance matrix $\Sigma$. Similarly, $\mathcal{U}(a,b)$ indicates the uniform distribution with non-zero density from $a$ to $b$. $\|\cdot\|_{2}$ is overloaded and indicates the $l_{2}$ norm in the case of a vector and the spectral norm in the case of a matrix. We denote the constraint set $\mathcal{C}$ as $\mathcal{C}=\mathcal{C}_{1}\bigcap\mathcal{C}_{2},\dots,\bigcap\mathcal{C}_{N_{\mathrm{C}}}$, where $N_{\mathrm{C}}$ is the total number of constraints and $\bigcap$ denotes intersection. Here, $\mathcal{C}_{i}=\{x\mid f_{c_{i}}(x)\leq 0\}$ with $f_{c_{i}}:\mathbb{R}^{K\times L}\rightarrow\mathbb{R}\;\forall\ c_{i}\in[1,\dots,N_{\mathrm{C}}]$. $\lambda_{\mathrm{max}}(M)$ and $\lambda_{\mathrm{min}}(M)$ refer to the largest and the smallest eigen values of the square matrix $M$. The rank of the matrix $M$ is indicated by $rank(M)$.

Example: The stocks dataset has $6$ channels ($K=6$) with $96$ timestamps in each channel ($L=96$). The first $4$ channels represent the opening price ($o$), the highest price ($h$), the lowest price ($l$), and the closing price ($c$), and each timestamp represents a day. The OHLC constraint, *i.e.,* the opening and closing prices should lie between the highest and the lowest prices, is given by $o-h\leq 0$, $c-h\leq 0$, $l-o\leq 0$, and $l-c\leq 0$. Additionally, a mean equality constraint on the closing price is expressed as $\frac{1}{L}\left(\sum_{u=1}^{L}c(u)\right)-\mu_{c}\leq 0$ and $\mu_{c}-\frac{1}{L}\left(\sum_{u=1}^{L}c(u)\right)\leq 0$, where $\mu_{c}$ is the required mean.

### 2.1 Background and Related Work

Figure: Figure 3: CPS tracks the real data samples as the number of constraints increases. Increasing the number of constraints reduces the size of the constraint set, and an ideal approach should effectively generate samples that resemble the real time series samples that belong to the constraint set. Here, we show a qualitative example from the Stocks dataset. Observe that CPS accurately tracks the real sample that concurs with the specified constraints while other approaches suffer.
Refer to caption: /html/2410.12652/assets/x3.png

GANs (Goodfellow et al., ) have been the popular choice for time series generation(Yoon et al., ; Donahue et al., ; Srinivasan & Knottenbelt, ; Ni et al., ). Recently, DMs have dominated the landscape of image, video, and audio generation (Rombach et al., ; Ho et al., ; Kong et al., ). Denoising DMs (Ho et al., ; Dhariwal & Nichol, ) generate samples by learning to gradually denoise clean data, sampled from the data distribution $P_{\mathrm{data}}$, corrupted with Gaussian noise. Denoising Diffusion Probabilistic Models (DDPMs) (Ho et al., ) define a Markovian forward noising process, where the clean data sample $x$, referred to as $z_{0}$, is transformed into $z_{T}$ with iterative Gaussian corruption for $T$ noising steps, such that $z_{T}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. With abuse of notation, $\mathbf{0}$ represents zero mean, and $\mathbf{I}$ represents the identity covariance. The forward process introduces $T$ conditional Gaussian distributions with fixed covariance matrices governed by the diffusion coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, where $\bar{\alpha}_{0}=1,\bar{\alpha}_{t}\in[0,1]$, and $\bar{\alpha}_{t-1}>\bar{\alpha}_{t}\ \forall\ t\in[1,T]$. Formally, $q_{t}(z_{t}\mid z_{0})$ is the PDF of the conditional Gaussian distribution at the forward step $t$ with mean $\sqrt{\bar{\alpha}_{t}}z_{0}$ and covariance matrix $\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}$. The PDF associated with the marginal distribution at $t=0$ is given by $q_{0}=p_{\mathrm{data}}$.

The sample generation or the reverse process is also Markovian, where we autoregressively sample from $T$ Gaussian distributions with fixed covariance matrices, indicated by PDFs $p_{\theta,t}(z_{t-1}\mid z_{t})\ \forall\ t\in[1,T]$, to get from $z_{T}$ to $z_{0}$, where $z_{T}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. The means of $p_{\theta,t}(z_{t-1}\mid z_{t})$ are learned using neural networks. DDPMs are trained to maximize the log-likelihood of observing the clean data, *i.e.,* $\mathrm{log}\ p_{\theta}(z_{0})$, where $p_{\theta}(z_{0})=\int p_{\theta}(z_{0:T})dz_{1:T}$. The joint PDF $p_{\theta}(z_{0:T})$ can be factorized as $p(z_{T})\prod_{t=1}^{T}p_{\theta,t}(z_{t-1}\mid z_{t})$, due to the Markovian nature of the reverse process, with $p(z_{T})$ representing the PDF of $\mathcal{N}(\mathbf{0},\mathbf{I})$. With successive reparametrizations, the training objective can be simplified into the following denoising objective:

$$ $\mathbb{E}_{z_{0}\sim P_{\mathrm{data}},\epsilon\sim\mathcal{N}(\mathbf{0},\mathbf{I}),t\sim\mathcal{U}(1,T)}\left[\|\epsilon-\epsilon_{\theta}(z_{t},t)\right\|^{2}_{2}],$ (1) $$

where $\epsilon_{\theta}(z_{t},t)$ is trained to estimate the noise $\epsilon$ from $z_{t}$, and $z_{t}=\sqrt{\bar{\alpha}_{t}}z_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon$, with $t$ ranging from $1$ to $T$. Denoising Diffusion Implicit Models (DDIMs) Song et al. () propose a non-Markovian forward process and, accordingly, a novel mechanism for sample generation given by

$$ $z_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}-\sigma_{t}^{2}}\epsilon_{\theta}(z_{t},t)+\sigma_{t}\epsilon.$ (2) $$

Here, $\hat{z}_{0}(z_{t};\epsilon_{\theta})=\frac{z_{t}-\sqrt{1-\bar{\alpha}_{t}}\epsilon_{\theta}(z_{t},t)}{\sqrt{\bar{\alpha}_{t}}}$ is the posterior mean estimate, and $\sigma_{t}$ is a control parameter that dictates determinism in the sampling process. Song et al. () show that Eq. [2](#S2.E2) corresponds to the following reverse process:

$$ $p_{\theta,t}(z_{t-1}\mid z_{t})=\left\{\begin{matrix}p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))&\mathrm{if}\ t=1,\\ q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta}))&\mathrm{otherwise},\\ \end{matrix}\right.$ (3) $$

where $q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta}))$ represents the PDF of the Gaussian distribution with mean $\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}-\sigma_{t}^{2}}\epsilon_{\theta}(z_{t},t)$ and covariance matrix $\sigma_{t}^{2}\mathbf{I}$. Similarly, $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))$ is the PDF of the Gaussian distribution with mean $\hat{z}_{0}(z_{1};\epsilon_{\theta})$ and covariance matrix $\sigma_{1}^{2}\mathbf{I}$. This reverse sampling process can be viewed as obtaining the posterior mean estimate $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ and transforming it to the noise level for step $t-1$. CPS builds on Eq. [2](#S2.E2).

Time Series-specific tasks like forecasting (Rasul et al., ; Yan et al., ; Biloš et al., ) and imputation (Tashiro et al., ; Alcaraz & Strodthoff, ) have been addressed using conditional DMs as well as guidance-based approaches (Li et al., ). Alcaraz & Strodthoff () and Narasimhan et al. () have explored conditional time series generation for various domains, such as medical, energy, etc. These works aim to sample from a conditional distribution. However, there are limited prior works that focus on generating constrained samples. Recent works (Wang et al., ; ) focus on generating constrained counterfactual explanations for time series classification and forecasting by perturbing selected time stamps of a synthesized seed sample. However, these approaches do not provide any mechanism to induce realism during the constrained generation process other than staying near the seed sample. Recently, (Coletta et al., ) proposed three approaches - Loss-DiffTime, a training-based approach where constrained samples are generated with constraints as conditional input to the generator, Guided DiffTime, which uses guidance gradients from differentiable constraint functions to guide the sample generation towards a constraint set, and Constrained Optimization Problem (COP), which projects a seed sample to the constraint set while using the critic function from any Wasserstein GAN (Arjovsky et al., ) as a realism enforcer. Loss-DiffTime is not scalable to new constraints as the generator has to be retrained for every new constraint, while Guided DiffTime does not guarantee constraint satisfaction even for convex constraint sets. We compare our proposed approach against Guided DiffTime and COP on many real-world datasets and highlight our advantages.

Finally, we note that constrained time series generation can be viewed as controlling the outputs of a generative model. Controlled generation is a well-established class of problems in image generation and occurs in multiple formulations, such as solving inverse problems (Rout et al., ; ; ; Chung et al., ), personalization (Rout et al., ; Ruiz et al., ), text-to-image generation (Rombach et al., ; Ramesh et al., ), and text-based image editing (Kawar et al., ; Choi et al., ). Indeed, we note that CPS (Sec. [3](#S3)) can be viewed as a constraint satisfaction (through projection) approach for time series, in the same spirit as gradient-based image personalization through diffusions (Rout et al., ). However, these works do not impose hard constraints, as described in the case of OHLC charts in Sec. [1](#S1).
Formally, the constrained time series generation problem is defined as follows:

Problem Setup. Consider a dataset $\mathcal{D}=\{x^{i}\}_{i=1}^{N_{\mathrm{D}}}$, where $N_{\mathrm{D}}$ denotes the number of samples, $x^{i}\sim P_{\mathrm{data}}$ with the density function $p_{\mathrm{data}}$ and $x^{i}\in\mathbb{R}^{K\times L}$. The goal is to generate $x^{\mathrm{gen}}\sim P_{\mathrm{data}}$ such that $x^{\mathrm{gen}}$ belongs to the constraint set $\mathcal{C}=\mathcal{C}_{1}\bigcap\mathcal{C}_{2},\dots,\bigcap\mathcal{C}_{N_{\mathrm{C}}}$, where $N_{\mathrm{C}}$ denotes the number of constraints. Here, $\mathcal{C}_{i}=\{x\mid f_{c_{i}}(x)\leq 0\}$ with $f_{c_{i}}:\mathbb{R}^{K\times L}\rightarrow\mathbb{R}$. To put it more succinctly,

$$ $\displaystyle x^{\mathrm{gen}}\coloneq\operatorname*{arg\,min}_{x}-\mathrm{log}\ p_{\mathrm{data}}(x)\ \textrm{s.t.}\ f_{c_{i}}(x)\leq 0,\ \forall\ c_{i}\in[1,N_{\mathrm{C}}],$ (4) $$

where the objective is to find a maximum likelihood sample in the constraint set.

## 3 Constrained Posterior Sampling

To generate realistic samples with high likelihood, our approach assumes the availability of a pre-trained diffusion model trained on the dataset $\mathcal{D}$. Given the diffusion denoiser model $\epsilon_{\theta}$, we propose Constrained Posterior Sampling (CPS, check Fig. [4](#S3.F4)) to restrict the domain of a generated sample without sacrificing sample quality. Described in Algorithm [1](#alg1), CPS effectively guides the diffusion denoising process towards the constraint set.

We follow the typical DDIM inference procedure. Starting with a sample from the standard normal distribution $\mathcal{N}(\mathbf{0},\mathbf{I})$ (line 1), we perform sequential denoising (lines 2 to 10). Line 3 refers to the forward pass through the denoiser to obtain the noise estimate $\epsilon_{\theta}(z_{t},t)$. After every denoising step, we obtain the posterior mean estimate $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ (line 4). We then project this estimate towards the constraint set $\mathcal{C}$ to obtain the projected posterior mean estimate $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ (line 5). Later, we perform a DDIM reverse sampling step with $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ and $\epsilon_{\theta}(z_{t},t)$ to obtain $z_{t-1}$ (lines 7-9).

Figure: Figure 4: Our proposed Constrained Posterior Sampling approach. We show the graphical model for one step of denoising in CPS, as outlined in Algorithm [1](#alg1).
Refer to caption: /html/2410.12652/assets/x4.png

The projection step in line 5 solves an optimization problem with the objective function $\frac{1}{2}(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}$ + $\gamma(t)\Pi(z))$. The first term of the objective function ensures that $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ is close to $\hat{z}_{0}(z_{t};\epsilon_{\theta})$, thereby ensuring that $z_{t-1}$ is not heavily perturbed for the denoiser to perform poorly. We define the constraint violation function $\Pi:\mathbb{R}^{K\times L}\rightarrow\mathbb{R}$ as $\Pi(z)=\sum_{c_{i}=1}^{N_{\mathrm{C}}}\mathbf{max}(0,f_{c_{i}}(z))$, such that $\Pi(z)=0$ if $z\in\mathcal{C}$ and $\Pi(z)>0$ otherwise. For the denoising step $t$, the constraint violation function is scaled by a time-varying penalty coefficient $\gamma(t)$.

Figure: Algorithm 1 Constrained Posterior Sampling

Our key intuition is to design $\gamma(t)$ as a strictly decreasing function of $t$ that takes small values for the initial denoising steps ($t$ close to $T$) and tends to $\infty$ for the final denoising steps. This ensures that the constraint satisfaction is not heavily enforced during the initial denoising steps when the signal-to-noise ratio in $z_{t}$ is very low. Given the requirements for the penalty coefficient, we choose $\gamma(t)=e^{1/(1-\bar{\alpha}_{t-1})}$ such that $\gamma(t)$ is close to 0 for the initial denoising steps ($\gamma(T)\simeq e$) and $\gamma(t)\rightarrow\infty$ for $t=1$. Note that our choice of $\gamma(t)$ ensures that $\gamma(t)$ is strictly decreasing with respect to $t$ since $\bar{\alpha}_{t}$ strictly decreases with $t$.

Observe that CPS is the DDIM sampling process with one change. We replace the posterior mean estimate $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ with the projected posterior mean estimate $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$. Additionally, CPS can be viewed similarly to the penalty-based methods to solve a constrained optimization problem. With each progressing denoising update, the penalty coefficient increases, thereby pushing the posterior mean estimate towards the constraint set.

We do not add noise after the final denoising step ($\sigma_{1}=0$). This ensures that the efforts of the final projection step towards constraint satisfaction are not compromised by additional noise. For convex constraint sets with assumptions on the convexity of the constraint definition functions $f_{c_{i}}$, we note that the projection step is an unconstrained minimization of a convex function with the optimal constraint violation value being 0 if $\gamma(1)>0$. With a suitable choice of solvers (Diamond & Boyd, ), the optimal solution can be obtained for these cases, thereby ensuring constraint satisfaction ($\Pi(\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))=0$) when $\gamma(1)$ tends to $\infty$.

Note that CPS satisfies the key requirements of an ideal constrained generation approach. CPS can handle multiple constraints without any training requirements. Further, CPS does not require additional critics to enforce realism, as our key intuition is that the successive denoising steps address the adverse effects of the projection step. Finally, CPS is hyperparameter-free as off-the-shelf solvers can perform the unconstrained optimization step in line 5. Our key observation is that unlike heuristically setting the guidance weights (Coletta et al., ), we can choose the parameters of the solvers using principled approaches from the vast optimization literature (Nocedal & Wright, ).

### 3.1 Theoretical Justification

Now, we provide a detailed analysis of the effect of modifying the traditional DDIM sampling process with CPS. For ease of explanation, we consider $z\in\mathbb{R}^{n}$. We indicate the identity matrix in $\mathbb{R}^{n\times n}$ as $\mathbf{I}_{n}$. First, we describe the exact distribution from which the samples are generated. For this, we make the following assumption.

###### Assumption 1 .

Let the constraint set be $\mathcal{C}=\{z\mid f_{\mathcal{C}}(z)=0\}$, where $f_{\mathcal{C}}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and the penalty function $\Pi(z)=\|f_{\mathcal{C}}(z)\|_{2}^{2}$ has $L$-Lipschitz continuous gradients, i.e., $\|\nabla\Pi(u)-\nabla\Pi(v)\|_{2}\leq L\|u-v\|_{2}\ \forall\ u,v\in\mathbb{R}^{n}$.

###### Theorem 1 .

Suppose Assumption [1](#Thmassumption1) holds. Given a denoiser $\epsilon_{\theta}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ for a diffusion process with noise coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, if $\gamma(t)>0\ \forall\ t\in[1,T]$, the denoising step in Algorithm [1](#alg1) is equivalent to sampling from the following conditional distribution:

$$ $p_{\theta,t}(z_{t-1}\mid z_{t})=\left\{\begin{matrix}p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))&\mathrm{if}\ t=1,\\ q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))&\mathrm{otherwise}.\end{matrix}\right.$ (5) $$

Here, $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))$ indicates the PDF of $\mathcal{N}\left(\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}),\sigma_{1}^{2}\mathbf{I}_{n}\right)$, and $q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))$ indicates the PDF of $\mathcal{N}\left(\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}-\sigma_{t}^{2}}\epsilon_{\theta}(z_{t},t),\sigma_{t}^{2}\mathbf{I}_{n}\right)$. $\sigma_{1},\dots,\sigma_{T}$ denote the DDIM control parameters, and $\gamma(t)$ indicates the penalty coefficient for the denoising step $t$ in Algorithm [1](#alg1).

Intuitively, Algorithm [1](#alg1) can be viewed as replacing $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ with $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ and following the DDIM sampling process. Therefore, the reverse process PDFs are obtained by replacing $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ with $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ in Eq. [3](#S2.E3). More formally, under Assumption [1](#Thmassumption1), the projection step (line 5) can be written as a series of gradient updates that transform $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ to $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$. Having Lipschitz continuous gradients for $f_{\mathcal{C}}$ allows for fixed step sizes which can guarantee a reduction in the value of the objective function $\frac{1}{2}(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}$ + $\gamma(t)\|f_{\mathcal{C}}(z)\|_{2}^{2})$ with each gradient update. We refer the readers to Sec. [A.1](#A1.SS1) in the Appendix for the detailed proof. Now, we investigate the convergence properties for Algorithm [1](#alg1) under the following assumption.

###### Assumption 2 .

The real data distribution is $\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$, where $\mu\in\mathbb{R}^{n}$, and the constraint set $\mathcal{C}$ is defined as $\mathcal{C}=\{z\mid Az=y\}$ with $A\in\mathbb{R}^{m\times n}$ such that $rank(A)=n\leq m$. Additionally, for the real data distribution $\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$ and the constraint set $\mathcal{C}=\{z\mid Az=y\}$, there exists a unique solution to Eq. [4](#S2.E4), indicated by $x^{*}$.

We note that Assumption [2](#Thmassumption2) ensures the existence of a unique solution to the linear problem $Ax=y$.
While there exist many efficient methods to solve such problems under this assumption, the focus of this paper is not on solving this problem efficiently. Instead, we use this well-studied problem as a framework to analyze the convergence properties of Algorithm [1](#alg1), providing valuable insights for better practical performance.

###### Theorem 2 .

Suppose Assumption [2](#Thmassumption2) holds.
For a diffusion process with noise coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, where $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]\ \forall\ t\in[0,T]$, if $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}$ and $\gamma(t)=\frac{2k(T-t+1)}{\lambda_{\mathrm{min}}(A^{T}A)}$ with any design parameter $k>1$,
then in the limit as $T\rightarrow\infty$, Algorithm [1](#alg1) returns $x^{\mathrm{gen}}$ such that:

$$ $\displaystyle\|x^{\mathrm{gen}}-x^{*}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{k}\left(\|x^{*}\|_{2}+\|\mu\|_{2}\right).$ $$

We refer the readers to Sec. [A.2](#A1.SS2) in the Appendix for detailed proof. Briefly, the proof in Sec. [A.2](#A1.SS2) indicates that the terminal error $\|x^{\mathrm{gen}}-x^{*}\|_{2}$ reduces to 0 as $T,k\rightarrow\infty$, thereby ensuring that Algorithm [1](#alg1) converges to the true solution. From the proof, we observe that under Assumption [2](#Thmassumption2), the convergence can be guaranteed when the penalty coefficient is set to very large values for the final denoising step. This is in accordance with our choice of penalty coefficients which assumes very large values for the final denoising step.

###### Assumption 1 .

###### Theorem 1 .

###### Assumption 2 .

###### Theorem 2 .

## 4 Experiments

This section describes the experimental procedure, including the wide range of datasets and metrics used to evaluate CPS against the state-of-the-art constrained generation approaches.

Datasets: We use real-world datasets from different domains, such as stocks (Yoon et al., ), air quality (Chen, ), and traffic (Hogue, ). Specifically, we test the performance of CPS on both conditional and unconditional variants of these datasets. We also evaluate our approach on a simulated sinusoidal waveforms dataset to generate sinusoids with varying amplitudes, phases, and frequencies specified as constraints.

Our evaluation procedure is framed to test any approach for generating the maximum likelihood sample from a constraint set, such that the real time series samples from the constraint set were never seen during training. To achieve this, from every sample in the test dataset, we first extract an exhaustive set of features such that only one test sample exists per set of features. These features are considered constraints, which we impose on the generative model.

Constraints: We extract the following features to be used as constraints - *mean*, *mean consecutive change*, *argmax*, *argmin*, *value at argmax*, *value at argmin*, *values at timestamps 1, 24, 48, 72, & 96*. For the stocks dataset, we additionally impose the natural OHLC constraint, *i.e.,* the opening and closing prices should be bounded by the highest and the lowest prices. Similarly, for the sinusoidal waveforms dataset, we extract the locations and values of the peaks and valleys and the trend from a peak to its adjacent valley. Note that these constraints can be written in the form $Ax\leq 0$. Projection to such constraint sets is easy and can be handled by numerous off-the-shelf solvers (Diamond & Boyd, ; Virtanen et al., ). This allows us to analyze the effect of the sampling process without worrying about the off-the-shelf solvers that influence the projection step. We provide a budget of $0.01$ for constraint violation.

Baselines: We compare against the Constrained Optimization Problem (COP) approach (Coletta et al., ) and its fine-tuning variant, which is referred to as COP-FT. COP projects a random sample from the training dataset to the required set of constraints, whereas COP-FT projects a generated sample. Both these variants rely on a discriminator to enforce realism after perturbation. We also compare our approach against Guided DiffTime (Coletta et al., ), a guidance-based diffusion sampling approach. It is unfair to compare against Loss DiffTime, a training-based approach proposed by (Coletta et al., ), as the other methods were not explicitly trained for adhering to a constraint set. All baselines, except COP, utilize the same Time Weaver-CSDI denoiser backbone (Narasimhan et al., ) for fair comparison.

Figure: Figure 5: CPS outperforms existing baselines with increasing number of constraints. Note that constraints are the features extracted from real time series samples. We gradually increase the number of constraints imposed on the generative model. Observe that CPS achieves the lowest DTW score for any number of constraints while having the best sample quality, indicated by the lowest FTSD metric. This result is in accordance with the qualitative example shown in Fig. [3](#S2.F3).
Refer to caption: /html/2410.12652/assets/x5.png

Metrics: We evaluate the performance of CPS on three fronts - sample quality, ability to track the test time series, and constraint violation. For sample quality, we use the Frechet Time Series Distance (FTSD) metric (Narasimhan et al., ; Paul et al., ) for the unconditional setting and the Joint Frechet Time Series Distance (J-FTSD) metric (Narasimhan et al., ) for the conditional setting. The FTSD metric is also referred to as Context-FID (Paul et al., ). For simplicity, we indicate both these metrics by Frechet Distance or FD. Additionally, we show the Train on Synthetic and Test on Real (TSTR) metric for sample quality. For TSTR, we choose random imputation as the task with 75% masking. We train the TimesNet model (Wu et al., ) for imputation on the synthesized training data, generated with constraints, and evaluate the trained model on the real test data for imputation performance. We report the mean squared error (MSE) on the real test set as the TSTR metric. Lower MSE indicates accurate modeling of the true data distribution.

Figure: Figure 6: CPS provides high-fidelity synthetic time series samples that match real time series data. Here, we show a qualitative comparison between the baselines (Guided DiffTime and COP-FT) and CPS for six different experimental settings. As described in Sec. [4](#S4), the real test time series samples from which the constraints are extracted are shown in blue. Observe that across datasets, CPS generates high-fidelity samples that match the ground truth, while the baselines suffer to generate meaningful qualitative results.
Refer to caption: /html/2410.12652/assets/x6.png

From our evaluation procedure, note that we aim to enforce one test sample per set of constraints. Therefore, an ideal approach is expected to generate a sample that is similar to that single test sample. To estimate this, we report the Dynamic Time Warping (DTW) metric (Müller, ) and the Structural Similarity Index Measure (SSIM) metric (Nilsson & Akenine-Möller, ). Though SSIM is typically used for images, in essence, both these metrics capture the similarity between the generated sample and the real test sample that belongs to the constraint set. Similarly, for constraint violation, we report the ratio of the generated samples that do not belong to the constraint set to the total number of test samples. We also report the average constraint violation magnitude.

A detailed discussion on the baselines, metrics, etc., is provided in Sec. [D.3](#A4.SS3) and Sec. [B](#A2) of the Appendix, respectively. Across all metrics, CPS outperforms the baselines on real-world and simulated datasets in conditional and unconditional settings. We provide intuitive reasons, backed by empirical evidence, for these performance gains by answering the following key questions:

How well does CPS generate realistic samples that belong to the constraint set? We argue about the performance of CPS based on the sample quality and the constraint violation metrics in Table [1](#S4.T1). As the constraint sets used in our experiments are convex, both CPS and COP variants can almost always ensure constraint satisfaction using off-the-shelf solvers. However, Guided DiffTime struggles severely to generate samples that belong to the constraint set. This is clearly observed in the stocks dataset, where Guided DiffTime has an average constraint violation magnitude of 1128. With respect to sample quality, we observe that CPS provides the lowest FD and TSTR values.

Even though Guided DiffTime provides comparable TSTR values for some settings, we note that the generated samples are very less likely to belong to the constraint set. Therefore, guidance gradients alone are insufficient to drive the sample generation process to the constraint set. Similarly, there exists a considerable difference in performance between both COP variants, specifically for the conditional setting. Here, our key observation is that conditional generation provides a seed sample for COP-FT that lies close to the constraint set. Therefore, projection does not degrade the sample quality by a lot. However, the sample quality degradation due to projection is significant for COP, and it can be observed through very high values of J-FTSD in Table [1](#S4.T1). Therefore, our key insight is that COP is influenced by the choice of initial seed.

**Table 1: CPS outperforms existing baselines on sample quality and similarity metrics. Yellow corresponds to sample quality metrics, and orange and violet correspond to similarity and constraint violation metrics, respectively. The best approach is shown in bold for each metric. Overall, we observe that CPS maintains high sample quality (very low FD and TSTR values) and the highest similarity with real time series samples (best values for the DTW and SSIM metrics). The mean and standard deviation for the TSTR values are obtained from the results for 3 seeds. Our key intuition is that the adverse effects of the projection step are nullified by the subsequent denoising steps. Note that as the constraints are all convex, the COP variants and CPS can achieve very low constraint violation.**
| Metric | Approach | Air Quality | Air Quality<br>(Conditional) | Traffic | Traffic<br>(Conditional) | Stocks | Waveforms |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | Guided<br>DiffTime | 0.7457 | 3.1883 | 0.5351 | 0.5638 | 1.2575 | 0.3108 |
|  | COP-FT | 0.3793 | 0.9931 | 0.8156 | 0.8135 | 0.0759 | 1.8419 |
|  | COP | 0.2165 | 27.9425 | 0.9242 | 43.2472 | 0.0701 | 1.6627 |
| Frechet<br>Distance<br>($\downarrow$) | CPS (Ours) | 0.0234 | 0.6039 | 0.2077 | 0.2812 | 0.0023 | 0.0029 |
|  | Guided<br>DiffTime | 0.29$\pm$0.015 | 0.25$\pm$0.003 | 0.30$\pm$0.01 | 0.28$\pm$0.01 | 0.05$\pm$0.001 | 0.005$\pm$0.001 |
|  | COP-FT | 0.23$\pm$0.005 | 0.19$\pm$0.002 | 0.32$\pm$0.01 | 0.28$\pm$0.01 | 0.048$\pm$0.001 | 0.023$\pm$0.001 |
|  | COP | 0.22$\pm$0.002 | 0.22$\pm$0.003 | 0.33$\pm$0.01 | 0.32$\pm$0.01 | 0.048$\pm$0.001 | 0.024$\pm$0.001 |
| TSTR<br>($\downarrow$) | CPS (Ours) | 0.19$\pm$0.003 | 0.19$\pm$0.003 | 0.29$\pm$0.01 | 0.28$\pm$0.01 | 0.041$\pm$0.001 | 0.005$\pm$0.001 |
|  | Guided<br>DiffTime | 6.74$\pm$8.18 | 4.28$\pm$5.66 | 4.38$\pm$1.25 | 1.31$\pm$1.01 | 7.84$\pm$7.24 | 1.67$\pm$1.15 |
|  | COP-FT | 3.52$\pm$2.08 | 2.01$\pm$1.24 | 4.61$\pm$1.08 | 1.26$\pm$0.87 | 0.90$\pm$1.41 | 1.19$\pm$0.64 |
|  | COP | 3.72 $\pm$ 2.14 | 3.72 $\pm$ 2.12 | 5.16 $\pm$ 1.34 | 4.94 $\pm$ 1.08 | 0.88 $\pm$ 1.39 | 1.16 $\pm$ 0.65 |
| DTW<br>($\downarrow$) | CPS (Ours) | 2.35$\pm$1.48 | 1.83$\pm$1.16 | 3.41$\pm$1.47 | 0.84$\pm$0.62 | 0.20$\pm$0.71 | 0.23$\pm$0.17 |
|  | Guided<br>DiffTime | 0.18$\pm$0.13 | 0.38$\pm$0.18 | 0.16$\pm$0.16 | 0.9$\pm$0.11 | 0.09$\pm$0.09 | 0.37$\pm$0.3 |
|  | COP-FT | 0.19$\pm$0.11 | 0.48$\pm$0.16 | 0.10$\pm$0.14 | 0.89$\pm$0.14 | 0.15$\pm$0.10 | 0.35$\pm$0.11 |
|  | COP | 0.17$\pm$0.11 | 0.17$\pm$0.11 | 0.09$\pm$0.13 | 0.09$\pm$0.13 | 0.14$\pm$0.09 | 0.39$\pm$0.12 |
| SSIM<br>($\uparrow$) | CPS (Ours) | 0.38$\pm$0.15 | 0.52$\pm$0.15 | 0.31$\pm$0.20 | 0.95$\pm$0.07 | 0.73$\pm$0.26 | 0.96$\pm$0.05 |
|  | Guided<br>DiffTime | 1.0 | 1.0 | 0.99 | 0.89 | 1.0 | 0.933 |
|  | COP-FT | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.003 |
|  | COP | 0.0 | 0.0 | 0.005 | 0.0 | 0.0 | 0.008 |
| Constraint<br>Violation<br>Rate<br>($\downarrow$) | CPS (Ours) | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
|  | Guided<br>DiffTime | 23.21 | 16.35 | 0.50 | 0.15 | 1128.22 | 5.23 |
|  | COP-FT | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0002 |
|  | COP | 0.0 | 0.0 | 0.0001 | 0.0 | 0.0 | 0.0003 |
| Constraint<br>Violation<br>Magnitude<br>($\downarrow$) | CPS (Ours) | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |

Does CPS handle unnatural artifacts that typically occur due to the projection step? While imposing constraints on the generation process, we note that even though Guided DiffTime generates a realistic sample, it fails to adhere to the constraints. On the other hand, COP variants adhere to constraints but generate samples with unnatural artifacts induced by the projection step. However, our key intuition is that CPS circumvents such artifacts using the iterative projection and denoising updates, where the adverse effects of the projection step are nullified by the subsequent denoising steps. The difference between the baseline approaches and CPS is significantly pronounced specifically in the waveforms dataset (check Fig. [6](#S4.F6)). The stark contrast between the generated sinusoid from CPS and other baselines is empirically supported with a $100\times$ reduction in the FD value.

How does CPS perform in comparison with baselines for a large number of constraints? We consider the stocks dataset as the OHLC condition introduces more than 400 constraints. With a large number of constraints, the feasible set size reduces. This necessitates the requirement of accurate guidance to generate samples from such constraint sets. In such settings, Guided DiffTime performs poorly. This can be attributed to the interaction between gradients for each constraint violation. The combination of these gradients, if not scaled appropriately, leads to poor guidance. Additionally, finding the correct set of guidance weights is practically very hard for a large number of constraints. Similarly, projection to small constraint sets affects the sample quality of COP variants, specifically when the initial seed is far away from the constraint set. While the baselines suffer with an increasing number of constraints, CPS gets rid of these issues by alternating projection and denoising updates. We observe this through the qualitative example from the stocks datasets in Fig. [3](#S2.F3). Quantitatively, we observe $14.5\%$ improvement in the TSTR and $67\%$ improvement in the FD metric when compared against the best-performing baseline.

Does CPS track the real test samples that adhere to the same set of constraints?
For a large number of constraints or a small constraint set, we expect the generated samples that satisfy the constraints to have a high degree of similarity with the real test samples from which we extract the constraints. To this end, we denote tracking real test samples as the property to have better similarity scores with the real sample as the number of constraints increases. In Fig. [5](#S4.F5), we note that CPS outperforms all baselines in the DTW metric for any number of constraints, thereby showing higher similarity with the real test samples. Note that out of all approaches, CPS has the best reduction in the DTW scores as the number of constraints increases. Simultaneously, we also note that the sample quality is unaffected or even improves for CPS with increasing constraints (lower FD scores). We observe that CPS’s performance is consistent across multiple real-world datasets, with significant improvements in the DTW values of around 33% for air quality, 77% for stocks, and 22% for the traffic dataset with respect to the best performing baseline.

Limitations. Although CPS outperforms all the compared baselines in standard evaluation metrics, we note that the projection step (line 5) in Algorithm [1](#alg1) can be time-consuming for some applications. This increases the overall sampling time of CPS, the trade-off being superior performance.
In time-critical applications, the sampling time can be reduced further by leveraging higher order moments and different initialization schemes (Rout et al., ). Additionally, the projection step is not necessary after every denoising step and can be adapted to the constraint violation magnitude.

## 5 Conclusion

We proposed Constrained Posterior Sampling – a novel training-free approach for constrained time series generation. CPS is designed such that it exploits off-the-shelf optimization routines to perform a projection step towards the constraint set after every denoising step. Through an array of sample quality and constraint violation metrics, we empirically show that CPS outperforms the state-of-the-art baselines in generating realistic samples that belong to a constraint set.

Future work. We aim to apply our approach for constrained trajectory generation in the robotics domain
with dynamics constraints, typically modeled by neural networks. Additionally, constrained time series generation readily applies to style transfer applications. Hence, we plan on extending the current work to perform style transfer from one time series to another by perturbing statistical features.

Reproducibility. The pseudo-code and hyper-parameter details have been provided in the Appendix to help reproduce the results reported in the paper. The source code will be released post publication.

## Appendix

## Appendix A Proofs

In this section, we provide the detailed proof for the theorems stated in the manuscript.

### A.1 Proof of Theorem 1

We first describe the assumption on the constraint set. The constraint set is defined as $\mathcal{C}=\{z\mid f_{\mathcal{C}}(z)=0\}$, where $f_{\mathcal{C}}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, and the penalty function $\Pi(z)=\|f_{\mathcal{C}}(z)\|_{2}^{2}$ has $L$-Lipschitz continuous gradients, *i.e.,* $\|\nabla\Pi(u)-\nabla\Pi(v)\|_{2}\leq L\|u-v\|_{2}\ \forall\ u,v\in\mathbb{R}^{n}$.

Line 7 of the Algorithm [1](#alg1) modifies the traditional DDIM sampling by replacing $\hat{z}_{0}(z_{t};\epsilon_{\theta})$ with $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$. Without this modification, the DDIM sampling denotes the following reverse process when started with $x_{T}\sim\mathcal{N}(\mathbf{0}_{n},\mathbf{I}_{n})$, where $\mathbf{0}_{n}$ indicates the zero mean vector in $\mathbb{R}^{n}$ and $\mathbf{I}_{n}$ is the identity matrix in $\mathbb{R}^{n\times n}$:

$$ $p_{\theta,t}(z_{t-1}\mid z_{t})=\left\{\begin{matrix}p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))&\mathrm{if}\ t=1,\\ q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta}))&\mathrm{otherwise},\\ \end{matrix}\right.$ (6) $$

where $q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta}))$ represents the PDF of the Gaussian distribution $\mathcal{N}\left(\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}-\sigma_{t}^{2}}\epsilon_{\theta}(z_{t},t),\sigma_{t}^{2}\mathbf{I}_{n}\right)$ with $\sigma_{t}$ as the DDIM control parameter. Similarly, $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))$ is the PDF of the Gaussian distribution with mean $\hat{z}_{0}(z_{1};\epsilon_{\theta})$ and covariance matrix $\sigma_{1}^{2}\mathbf{I}_{n}$ (Song et al., ).

Note that sampling from $q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta}))$ provides the DDIM sampling step (check Eq. [2](#S2.E2)).

We reiterate that the main modification with respect to the DDIM sampling approach is the projection step in line 5 of Algorithm [1](#alg1). Therefore, we first analyze the projection step,

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=\operatorname*{arg\,min}_{z}\frac{1}{2}\left(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}+\gamma(t)\|f_{\mathcal{C}}(z)\|_{2}^{2}\right).$ (7) $$

Here, $\hat{z}_{0}(z_{t};\epsilon_{\theta})=\frac{z_{t}-\sqrt{1-\bar{\alpha}_{t}}\epsilon_{\theta}(z_{t},t)}{\sqrt{\bar{\alpha}_{t}}}$ (line 3, predicted $z_{0}$). We will denote the objective function $\frac{1}{2}\left(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}+\gamma(t)\|f_{\mathcal{C}}(z)\|_{2}^{2}\right)$ as $g(z)$. Note that we replaced the constraint violation function $\Pi(z)$ by $\|f_{\mathcal{C}}(z)\|_{2}^{2}$ for this case. Given that $f_{\mathcal{C}}$ is a differentiable and convex with $\|f_{\mathcal{C}}\|_{2}^{2}$ having $L$-Lipschitz continuous gradients, Eq. [7](#A1.E7) can be written as a series of gradient updates with a suitable step size such that the value of the objective function decreases for each gradient update.

From the statement, we observe that $\gamma(t)>0\ \forall\ t\in[1,T]$. Under this condition and Assumption [1](#Thmassumption1), note that the function $g(z)$ is convex and has $\left(\frac{2+\gamma(t)L}{2}\right)$-Lipschitz continuous gradients, as $\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}$ has 2-Lipschitz continuous gradients, $\gamma(t)\|f_{\mathcal{C}}(z)\|_{2}^{2}$ has $\left(\gamma(t)L\right)$-Lipschitz continuous gradients, and the fraction $\frac{1}{2}$ makes $g(z)$ to have $\left(\frac{2+\gamma(t)L}{2}\right)$-Lipschitz continuous gradients. Let $\eta$ be the step size of the projection step. From Nocedal & Wright (), we know that $\eta\in(0,2/(2+\gamma(t)L))$ ensures that the objective function in Eq. [7](#A1.E7) reduces after each gradient update. We denote the gradient update as

$$ $\displaystyle\left.\begin{matrix}\leftidx^{n}{\hat{z}_{0}(z_{t};\epsilon_{\theta})}=\ \leftindex^{n-1}{\hat{z}_{0}(z_{t};\epsilon_{\theta})}-\eta\nabla_{z}(g(z))\end{matrix}\right|_{\leftindex^{n-1}{\hat{z}_{0}(z_{t};\epsilon_{\theta})}},$ (8) $$

where $\leftidx^{0}{\hat{z}_{0}(z_{t};\epsilon_{\theta})}=\hat{z}_{0}(z_{t};\epsilon_{\theta})$ and $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=\ \leftindex^{N_{\mathrm{pr}}}{\hat{z}_{0}(z_{t};\epsilon_{\theta})}$. Here, $N_{\mathrm{pr}}$ is the total number of gradient update steps.

The iteration in Eq. [8](#A1.E8) always leads to $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$ deterministically. Therefore, the projection step can be considered sampling from a Dirac delta distribution centered at $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$, *i.e.,* $\delta(z-\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))$. Consequently, using the law of total probability, the reverse process corresponding to the denoising step $t\ \forall\ t\in[2,T]$ in Algorithm [1](#alg1) is given by

$$ $\displaystyle p_{\theta,t}(z_{t-1}\mid z_{t})=\int{p_{\theta,t}(z_{t-1},\hat{z}_{0}\mid z_{t})}d\hat{z}_{0},$ $$

where $\hat{z}_{0}\in\mathbb{R}^{n}$. This can be simplified using Bayes’ rule as

$$ $\displaystyle p_{\theta,t}(z_{t-1}\mid z_{t})=\int{\delta(\hat{z}_{0}-\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0})}d\hat{z}_{0}.$ $$

The above equation stems from the fact that the distribution of $z_{0}$ conditioned on $z_{t}$ is a Dirac delta distribution centered at $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$. Since $\delta(x-y)=\delta(y-x)$ and using the sifting property of a Dirac delta function $\left(\int f(z)\delta(a-z)dz=f(a)\right)$, we get

$$ $p_{\theta,t}(z_{t-1}\mid z_{t})=q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))\ \forall\ t\ \in\ [2,T].$ (9) $$

Similarly, we repeat the steps for $t=1$, which gives

$$ $\displaystyle p_{\theta,1}(z_{0}\mid z_{1})$ $\displaystyle=\int{p_{\theta,1}(z_{0},\hat{z}_{0}\mid z_{t})}d\hat{z}_{0},$ $\displaystyle p_{\theta,1}(z_{0}\mid z_{1})$ $\displaystyle=\int{\delta(\hat{z}_{0}-\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0})}d\hat{z}_{0},$ $\displaystyle p_{\theta,1}(z_{0}\mid z_{1})$ $\displaystyle=p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta})).$ $$

Combining the two, we get

$$ $\displaystyle p_{\theta,t}(z_{t-1}\mid z_{t})=\left\{\begin{matrix}p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))&\mathrm{if}\ t=1,\\ q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))&\mathrm{otherwise},\end{matrix}\right.$ (10) $$

where $q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta}))$ represents the PDF of the Gaussian distribution $\mathcal{N}(\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}-\sigma_{t}^{2}}\epsilon_{\theta}(z_{t},t),\sigma_{t}^{2}\mathbf{I}_{n})$ with $\sigma_{t}$ as the DDIM control parameter. Similarly, $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))$ is the PDF of the Gaussian distribution with mean $\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta})$ and covariance matrix $\sigma_{1}^{2}\mathbf{I}_{n}$ (Song et al., ). This is the same as Eq. [5](#S3.E5).

We note that the value of $\sigma_{1}$ is set to 0 in Algorithm [1](#alg1). However, similar to (Song et al., ), for theoretical analysis, we consider a negligible value for $\sigma_{1}$ ($\sim 10^{-12}$) to ensure that the generative process is supported everywhere. In other words, $\sigma_{1}$ is chosen to be so low such that for $\sigma_{1}\simeq 0$, $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta}))\simeq\delta(z_{0}-\hat{z}_{0,\mathrm{pr}}(z_{1};\epsilon_{\theta})).$

Now, we show that the exact DDIM reverse process (check Eq. [6](#A1.E6)) can be obtained from Eq. [10](#A1.E10) in the case where there are no constraints. Here, note that in the absence of any constraint, the projection step can be written as $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=\operatorname*{arg\,min}_{z}\frac{1}{2}\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}$, in which case $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=\hat{z}_{0}(z_{t};\epsilon_{\theta})$.

For $t\in[2,T]$, using the law of total probability, we get

$$ $\displaystyle p_{\theta,t}(z_{t-1}\mid z_{t})=\int{\delta(\hat{z}_{0}-\hat{z}_{0}(z_{t};\epsilon_{\theta}))q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0})}d\hat{z}_{0},$ (11) $$

which simplifies further to

$$ $\displaystyle p_{\theta,t}(z_{t-1}\mid z_{t})=q_{\sigma,t}(z_{t-1}\mid z_{t},\hat{z}_{0}(z_{t};\epsilon_{\theta})).$ (12) $$

The above equation stems from the same sifting property of Dirac delta functions. The same applies to $t=1$, except that after the projection step, since there is no necessity for constraint satisfaction, we sample from $p_{\theta,\mathrm{init}}(z_{0}\mid\hat{z}_{0}(z_{1};\epsilon_{\theta}))$, which is a Gaussian distribution with mean $\hat{z}_{0}(z_{1};\epsilon_{\theta})$ and covariance matrix $\sigma_{1}^{2}\mathbf{I}_{n}$.

Combining both cases, we observe that without any constraints, the exact DDIM reverse process can be recovered from Algorithm [2](#S2.E2) for all $t\in[1,T]$.

### A.2 Proof of Theorem 2

We note that the intermediate samples in a $T$-step reverse sampling process are denoted as $z_{T},\dots,z_{0}$, where $z_{0}=x^{\mathrm{gen}}$ and $z_{T}\sim\mathcal{N}(\mathbf{0}_{n},\mathbf{I}_{n})$. Once again, we reiterate the assumptions. We consider the real data distribution to be Gaussian with mean $\mu\in\mathbb{R}^{n}$ and covariance matrix $\mathbf{I}_{n}$, *i.e.,* $\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$. The constraint set $\mathcal{C}$ is defined as $\mathcal{C}=\{z\mid Az=b\}$ with $A\in\mathbb{R}^{m\times n}$ such that $rank(A)=n$, where $m\geq n$. Additionally, for the real data distribution $\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$ and the constraint set $\mathcal{C}=\{z\mid Az=y\}$, there exists a unique solution to Eq. [4](#S2.E4), indicated by $x^{*}$.

Given that $rank(A)=n$ for $A\in\mathbb{R}^{m\times n}$ with $m\geq n$, we note that $(A^{T}A)^{-1}$ exists. Consequently, $\lambda_{\mathrm{min}}(A^{T}A)>0$.

From the theorem statement, we have $\gamma(t)=\frac{2k(T-t+1)}{\lambda_{\mathrm{min}}(A^{T}A)}$, with $k>1$. Immediately, we note that for all $t\in[1,T]$, $\gamma(t)>0$. More specifically, $t\in[1,T]$, $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$.

The proof is divided into 2 parts. First, we obtain the expression for $z_{t-1}$ in terms of $z_{t}$. Then, we obtain an upper bound for $\|z_{0}-x^{*}\|_{2}$, which is the same as $\|x^{\mathrm{gen}}-x^{*}\|_{2}$, as from Algorithm [1](#alg1) we note that $z_{0}=x^{\mathrm{gen}}$.

First, we note that for deterministic sampling, we have the DDIM control parameters $\sigma_{1}\dots\sigma_{T}=0$. Therefore, the DDIM reverse sampling step from Algorithm [1](#alg1) (line7) can be written as

$$ $\displaystyle z_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})+\sqrt{1-\bar{\alpha}_{t-1}}\epsilon_{\theta}(z_{t},t).$ (13) $$

Since the true data distribution is Gaussian, the optimal denoiser $\epsilon^{*}(z_{t},t)$ can be expressed analytically for any diffusion step $t$. Therefore, the deterministic sampling step can be written as

$$ $\displaystyle z_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})+\sqrt{1-\bar{\alpha}_{t-1}}\epsilon^{*}(z_{t},t).$ $$

We can obtain an analytical expression for the optimal denoiser from Lemma [1](#A4.EGx29). Using Eq. [27](#A1.E27) from Lemma [1](#A4.EGx29), we note that the optimal denoiser at the diffusion step $t$ is

$$ $\epsilon^{*}(z_{t},t)=-\sqrt{1-\bar{\alpha}_{t}}(\sqrt{\bar{\alpha}_{t}}\mu-z_{t}).$ (14) $$

Now, we obtain the expression for $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})$. Note that the constraint violation function is defined as $\Pi(z)=\|y-Az\|_{2}^{2}$. Consequently, we note that the objective function in line 5 of Algorithm [1](#alg1), *i.e.,* $\frac{1}{2}(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}+\gamma(t)\|y-Az\|_{2}^{2})$, is convex with respect to $z$ for $\gamma(t)>0$. As such, we use Lemmas [1](#A4.EGx29) and [2](#Thmlemma2) to obtain the following expression for $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})$:

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}y].$ (15) $$

We substitute the expressions for $\epsilon^{*}(z_{t},t)$ from Eq. [14](#A1.E14) and $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})$ from Eq. [15](#A1.E15), respectively, in addition to replacing $y$ with $Ax^{*}$, to obtain $z_{t-1}$ in terms of $z_{t}$:

$$ $\displaystyle z_{t-1}=$ $\displaystyle\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\left[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}y\right]$ $\displaystyle+\sqrt{1-\bar{\alpha}_{t-1}}(-\sqrt{1-\bar{\alpha}_{t}}(\sqrt{\bar{\alpha}_{t}}\mu-z_{t})),$ $\displaystyle z_{t-1}=$ $\displaystyle\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\left[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}y\right]$ $\displaystyle-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mu+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}z_{t},$ $\displaystyle z_{t-1}=$ $\displaystyle\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\left[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}Ax^{*}\right]$ $\displaystyle-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mu+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}z_{t},$ $\displaystyle z_{t-1}=$ $\displaystyle\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right]z_{t}$ $\displaystyle+\left[\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}-\bar{\alpha}_{t}\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right]\mu$ $\displaystyle-\left[\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mathbf{I}_{n}\right]\mu+\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}Ax^{*},$ $\displaystyle z_{t-1}=$ $\displaystyle\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right]z_{t}$ $\displaystyle+\left[(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right]\mu-\left[\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mathbf{I}_{n}\right]\mu$ $\displaystyle+\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}Ax^{*}.$ $$

On further simplification, we have

$$ $\displaystyle z_{t-1}=$ $\displaystyle K_{t}z_{t}+E_{t}\mu-F_{t}\mu+\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}Ax^{*},$ $$

where we have the following matrix definitions,

$$ $\displaystyle K_{t}=\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right],$ (16) $$

$$ $\displaystyle E_{t}=\left[(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right],$ (17) $$

$$ $\displaystyle F_{t}=\left[\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mathbf{I}_{n}\right].$ (18) $$

The goal is to obtain the upper bound for $\|x^{\mathrm{gen}}-x^{*}\|_{2}$. Note that $\|x^{\mathrm{gen}}-x^{*}\|_{2}=\|z_{0}-x^{*}\|_{2}$. So, first, we subtract $x^{*}$ from both sides to obtain

$$ $\displaystyle z_{t-1}-x^{*}=$ $\displaystyle K_{t}z_{t}+E_{t}\mu-F_{t}\mu+\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}Ax^{*}-x^{*}.$ $$

Further, we add and subtract $K_{t}x^{*}$ to the right side to obtain

$$ $\displaystyle z_{t-1}-x^{*}=$ $\displaystyle K_{t}z_{t}-K_{t}x^{*}+E_{t}\mu-F_{t}\mu+\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}Ax^{*}-x^{*}+K_{t}x^{*}.$ $$

We further simplify the above expression to obtain

$$ $\displaystyle z_{t-1}-x^{*}=$ $\displaystyle K_{t}\left(z_{t}-x^{*}\right)+E_{t}\mu-F_{t}\mu+K_{t}x^{*}+D_{t}x^{*},$ $$

where the matrix definition of $D_{t}$ is

$$ $\displaystyle D_{t}=\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}A-\mathbf{I}_{n}.$ (19) $$

Now, we obtain the expression for $\|z_{t-1}-x^{*}\|_{2}$ in terms of $\|z_{t}-x^{*}\|_{2}$.

$$ $\displaystyle\|z_{t-1}-x^{*}\|_{2}=\|K_{t}(z_{t}-x^{*})+E_{t}\mu-F_{t}\mu+K_{t}x^{*}+D_{t}x^{*}\|_{2}.$ $$

Applying the triangle inequality repeatedly, we get

$$ $\displaystyle\|z_{t-1}-x^{*}\|_{2}\leq\|K_{t}(z_{t}-x^{*})\|_{2}+\|K_{t}x^{*}\|_{2}+\|D_{t}x^{*}\|_{2}+\|E_{t}\mu\|_{2}+\|F_{t}\mu\|_{2}.$ (20) $$

Before obtaining the upperbound for $\|z_{0}-x^{*}\|$, for $\gamma(t)>0$, we will first show that $\|K_{t}\|_{2},\|D_{t}\|_{2},\|E_{t}\|_{2},\|F_{t}\|_{2}<1\ \forall\ t\ \in[1,T].$ Here $\|K_{t}\|_{2}$ refers to the spectral norm of the matrix $K_{t}$. To show this, we establish a few relationships that will be the recurring theme used in proving that $\|K_{t}\|_{2},\|D_{t}\|_{2},\|E_{t}\|_{2},\|F_{t}\|_{2}<1\ \forall\ t\ \in[1,T].$

The spectral norm of the matrix M is defined as $\|M\|_{2}=\max_{x\neq\mathrm{0}}\frac{\|Mx\|_{2}}{\|x\|_{2}}$. From this definition, we immediately note the following two inequalities.

- •
$\|Mx\|_{2}\leq\|M\|_{2}\|x\|_{2}$ as $\|M\|_{2}=\max_{x\neq\mathrm{0}}\frac{\|Mx\|_{2}}{\|x\|_{2}}$.
- •
$\|MN\|_{2}=\max_{x\neq\mathrm{0}}\frac{\|MNx\|_{2}}{\|x\|_{2}}\leq\max_{x\neq\mathrm{0}}\frac{\|M\|_{2}\|Nx\|_{2}}{\|x\|_{2}}\leq\max_{x\neq\mathrm{0}}\frac{\|M\|_{2}\|N\|_{2}\|x\|_{2}}{\|x\|_{2}}=\|M\|_{2}\|N\|_{2}.$

Further, we note that the following are well-established properties for spectral norms and positive definite matrices. Consider a positive definite matrix $M$, *i.e.,* $M\succ 0$.

- •
$\|M\|_{2}$ is equal to the largest eigen value of $M$, *i.e.,* $\lambda_{\mathrm{max}}(M)$.
- •
$\|M^{-1}\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(M)}$ as the eigenvalues of $M^{-1}$ are the reciprocal of the eigenvalues of $M$.
- •
$\|-M\|_{2}=\|M\|_{2}$.

We refer the readers to Lemmas [3](#Thmlemma3), [8](#Thmlemma8), and [10](#Thmlemma10), where we show that $\|K_{t}\|_{2},\|E_{t}\|_{2},\|F_{t}\|_{2}<1\ \forall\ t\ \in[1,T]$, if $\gamma(t)>0$.

Similarly, Lemma [6](#Thmlemma6) shows that $\|D_{t}\|_{2}<1\ \forall\ t\ \in[1,T]$, if $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}{A^{T}A}}$.

We first apply the inequality $\|Mx\|_{2}\leq\|M\|_{2}\|x\|_{2}$ to simplify Eq. [20](#A1.E20) as follows.

$$ $\displaystyle\|z_{t-1}-x^{*}\|_{2}\leq\|K_{t}\|_{2}\|z_{t}-x^{*}\|_{2}+\|K_{t}x^{*}\|_{2}+\|D_{t}x^{*}\|_{2}+\|E_{t}\mu\|_{2}+\|F_{t}\mu\|_{2}.$ (21) $$

Therefore, we can recursively obtain the upper bound for $\|z_{t}-x^{*}\|_{2}$ in terms of $\|z_{T}-x^{*}\|_{2}$. This process, repeated $T$ times, provides the upper bound for $\|z_{0}-x^{*}\|_{2}$.

$$ $\displaystyle\|z_{0}-x^{*}\|_{2}$ $\displaystyle\leq\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T}\|_{2}\|(z_{T}-x^{*})\|_{2}$ $\displaystyle\qquad+(\|K_{1}\|_{2}+\|K_{1}\|_{2}\|K_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|K_{T}\|_{2})\|x^{*}\|_{2}$ $\displaystyle\qquad+(\|D_{1}\|_{2}+\|K_{1}\|_{2}\|D_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|D_{T}\|_{2})\|x^{*}\|_{2}$ $\displaystyle\qquad+(\|E_{1}\|_{2}+\|K_{1}\|_{2}\|E_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|E_{T}\|_{2})\|\mu\|_{2}$ $\displaystyle\qquad+(\|F_{1}\|_{2}+\|K_{1}\|_{2}\|F_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|F_{T}\|_{2})\|\mu\|_{2}.$ (22) $$

Let $\lambda_{k}=\max_{t}\left(\|K_{1}\|_{2},\|K_{2}\|_{2},\dots,\|K_{T}\|_{2}\right)$. Since for $\gamma(t)>0$, $\|K_{1}\|_{2},\dots,\|K_{T}\|_{2}<1$, we note that $\lambda_{k}<1$.

Therefore, $\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T}\|_{2}$ can be upper bounded by $\lambda_{k}^{T}$.

Additionally, note that $\|K_{1}\|_{2}\|K_{2}\|_{2}\leq\|K_{1}\|_{2}$ as $\|K_{2}\|_{2}<1$. Therefore, $(\|K_{1}\|_{2}+\|K_{1}\|_{2}\|K_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|K_{T}\|_{2})$ can be upper bounded by $T\|K_{1}\|_{2}$.

Similarly, $(\|K_{1}\|_{2}\|D_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|D_{T}\|_{2})$ can be upperbounded by $(T-1)\|K_{1}\|_{2}$.

The same applies to $(\|K_{1}\|_{2}\|E_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|E_{T}\|_{2})$ and $(\|K_{1}\|_{2}\|F_{2}\|_{2}+\dots+\|K_{1}\|_{2}\|K_{2}\|_{2}\dots\|K_{T-1}\|_{2}\|F_{T}\|_{2})$.

Therefore, the upper bound in Eq. [22](#A1.E22) can be simplified as

$$ $\displaystyle\|z_{0}-x^{*}\|$ $\displaystyle\leq\lambda_{k}^{T}\|(z_{T}-x^{*})\|_{2}+T\|K_{1}\|_{2}\|x^{*}\|_{2}+(\|D_{1}\|_{2}+(T-1)\|K_{1}\|_{2})\|x^{*}\|_{2}$ $\displaystyle\qquad+(\|E_{1}\|_{2}+(T-1)\|K_{1}\|_{2})\|\mu\|_{2}+(\|F_{1}\|_{2}+(T-1)\|K_{1}\|_{2})\|\mu\|_{2}.$ (23) $$

Consequently, in Lemmas [4](#A4.EGx38), [7](#A4.EGx63), [9](#A4.EGx72), [10](#Thmlemma10), we show

$$ $\displaystyle\|K_{1}\|_{2}$ $\displaystyle\leq\frac{\sqrt{\bar{\alpha}_{1}}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}<1\ \text{if}\ \gamma(1)>0,$ $\displaystyle\|D_{1}\|_{2}$ $\displaystyle\leq\begin{matrix}\frac{1}{\gamma(1)\lambda_{\mathrm{min}}(A^{T}A))-1}<1\ \text{if}\ \gamma(1)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)},\end{matrix}$ $\displaystyle\|E_{1}\|_{2}$ $\displaystyle\leq\frac{1-\bar{\alpha}_{1}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}<1\ \text{if}\ \gamma(1)>0,$ $\displaystyle\|F_{1}\|_{2}$ $\displaystyle=0.$ (24) $$

For our choice of $\gamma(1)=\frac{2kT}{\lambda_{\mathrm{min}}(A^{T}A)}$, we first note that $\gamma(1)>0$ and $\gamma(1)>\frac{2}{\lambda_{\mathrm{max}}(A^{T}A)}$ for $k>1$. Therefore, we can rewrite the above inequalities as

$$ $\displaystyle\|K_{1}\|_{2}$ $\displaystyle\leq\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT},$ $\displaystyle\|D_{1}\|_{2}$ $\displaystyle\leq\begin{matrix}\frac{1}{2kT-1},\end{matrix}$ $\displaystyle\|E_{1}\|_{2}$ $\displaystyle\leq\frac{1-\bar{\alpha}_{1}}{1+2kT},$ $\displaystyle\|F_{1}\|_{2}$ $\displaystyle=0.$ (25) $$

Therefore, Eq. [23](#A1.E23) can be upper bounded using Eq. [25](#A1.E25) as shown below:

$$ $\displaystyle\|z_{0}-x^{*}\|_{2}$ $\displaystyle\leq\lambda_{k}^{T}\|(z_{T}-x^{*})\|_{2}+T\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|x^{*}\|_{2}+\left(\frac{1}{2kT-1}\right)\|x^{*}\|_{2}+$ (26) $\displaystyle\qquad\left(\frac{1-\bar{\alpha}_{1}}{1+2kT}\right)\|\mu\|_{2}+(T-1)\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|x^{*}\|_{2}+2(T-1)\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|\mu\|_{2}.$ $$

As $T\rightarrow\infty$, we observe the following:

$$ $\displaystyle\lim_{T\rightarrow\infty}\lambda_{k}^{T}\|(z_{T}-x^{*})\|_{2}$ $\displaystyle=0\ \ \ \left(\lambda_{k}<1\right),$ $\displaystyle\lim_{T\rightarrow\infty}T\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|x^{*}\|_{2}$ $\displaystyle=\left(\frac{\sqrt{\bar{\alpha}_{1}}}{2k}\right)\|x^{*}\|_{2}\ \ \ (\text{if}\ k>0),$ $\displaystyle\lim_{T\rightarrow\infty}\left(\frac{1}{2kT-1}\right)\|x^{*}\|_{2}$ $\displaystyle=0,$ $\displaystyle\lim_{T\rightarrow\infty}\left(\frac{1-\bar{\alpha}_{1}}{1+2kT}\right)\|\mu\|_{2}$ $\displaystyle=0,$ $\displaystyle\lim_{T\rightarrow\infty}(T-1)\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|x^{*}\|_{2}$ $\displaystyle=\left(\frac{\sqrt{\bar{\alpha}_{1}}}{2k}\right)\|x^{*}\|_{2}\ \ \ (\text{if}\ k>0),$ $\displaystyle\lim_{T\rightarrow\infty}2(T-1)\left(\frac{\sqrt{\bar{\alpha}_{1}}}{1+2kT}\right)\|\mu\|_{2}$ $\displaystyle=\left(\frac{\sqrt{\bar{\alpha}_{1}}}{k}\right)\|\mu\|_{2}\ \ \ (\text{if}\ k>0).$ $$

Therefore, in the limit $T\rightarrow\infty$, we have

$$ $\displaystyle\|z_{0}-x^{*}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{k}\left(\|x^{*}\|_{2}+\|\mu\|_{2}\right)\ \text{or},$ $\displaystyle\|x^{\mathrm{gen}}-x^{*}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{k}\left(\|x^{*}\|_{2}+\|\mu\|_{2}\right).$ $$

###### Lemma 1 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]$. The optimal denoiser $\epsilon^{*}(z_{t},t)$ is given by

$$ $\displaystyle\epsilon^{*}(z_{t},t)=-\sqrt{1-\bar{\alpha}_{t}}(\sqrt{\bar{\alpha}_{t}}\mu-z_{t}).$ $$

###### Proof.

We first observe the distribution of $z_{t}$.

For the diffusion forward process, we know that $z_{t}=\sqrt{\bar{\alpha}_{t}}z_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon$, where $\epsilon\sim\mathcal{N}(\mathbf{0}_{n},\mathbf{I}_{n})$.

Note that $z_{0}$ is a sample from the Gaussian distribution $\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$.

Consequently, we note that $z_{t}$ is a sample from the Gaussian distribution $\mathcal{N}(\sqrt{\bar{\alpha}_{t}}\mu+\mathrm{0}_{n},\bar{\alpha}_{t}\mathbf{I}_{n}+(1-\bar{\alpha}_{t})\mathbf{I}_{n})$. On simplification, we note that $z_{t}$ is a sample from $\mathcal{N}(\sqrt{\bar{\alpha}_{t}}\mu,\mathbf{I}_{n})$.

We denote the PDF of $z_{t}$’s marginal distribution as $q_{t}(z_{t})$.

Since we are using the optimal denoiser, the reverse process PDF at $t$, induced by the optimal denoiser, $p_{*,t}(z_{t})$ is the same as the forward process PDF at $t$, which is $q_{t}(z_{t})$.

Here, note that in Sec. [3](#S2.F3), we denote the reverse process PDF as $p_{\theta,t}$, where the reverse process is governed by the denoiser $\epsilon_{\theta}$. We replace this notation with $p_{*,t}(z_{t})$ as we are using the optimal denoiser.

Therefore, the score function at $t$ is given by $\nabla_{z_{t}}\mathrm{log}\ p_{*,t}(z_{t})=\nabla_{z_{t}}\mathrm{log}\ q_{t}(z_{t})$.

The score function for the Gaussian distribution $q_{t}(z_{t})$ with mean $\sqrt{\bar{\alpha}_{t}}\mu$ and covariance matrix $\mathbf{I}_{n}$, *i.e.,* $\nabla_{z_{t}}(\mathrm{log}\ q_{t}(z_{t}))$ is given by $\sqrt{\bar{\alpha}_{t}}\mu-z_{t}$.

Finally, Luo () shows that for the diffusion step $t$, the optimal denoiser can be obtained from the score function using the following expression:

$$ $\epsilon^{*}(z_{t},t)=-\sqrt{1-\bar{\alpha}_{t}}\nabla_{z_{t}}\mathrm{log}\ q_{t}(z_{t})\ \ \Rightarrow\ \ \epsilon^{*}(z_{t},t)=-\sqrt{1-\bar{\alpha}_{t}}(\sqrt{\bar{\alpha}_{t}}\mu-z_{t}).$ (27) $$

∎

###### Lemma 2 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]$. The projected posterior mean estimate, $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$, from the projection step in line 5 of Algorithm [1](#alg1) is given by

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=[I+\gamma(t)A^{T}A]^{-1}[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}y],$ $$

where the penalty coefficient from Algorithm [1](#alg1), $\gamma(t)>0\ \forall\ t\in[1,\dots,T].$

###### Proof.

We start with the unconstrained minimization in line 5 of Algorithm [1](#alg1), given by

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=\operatorname*{arg\,min}_{z}\frac{1}{2}\left(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}+\gamma(t)\|y-Az\|_{2}^{2}\right).$ $$

Note that we replaced the penalty function $\Pi(z)$ with $\|y-Az\|_{2}^{2}$, as we are required to generate a sample that satisfies the constraint $y=Az$.

Since the objective function is convex with respect to $z$, we obtain the global minimum by setting the gradient with respect to $z$ to 0, *i.e.,*

$$ $\displaystyle\nabla_{z}\left(\frac{1}{2}\left(\|z-\hat{z}_{0}(z_{t};\epsilon_{\theta})\|_{2}^{2}+\gamma(t)\|y-Az\|_{2}^{2}\right)\right)$ $\displaystyle=0,$ $\displaystyle\nabla_{z}\left(\frac{1}{2}\left(z^{T}z-2z^{T}\hat{z}_{0}(z_{t};\epsilon_{\theta})+\hat{z}_{0}(z_{t};\epsilon_{\theta})^{T}\hat{z}_{0}(z_{t};\epsilon_{\theta})\right)\right)+\gamma(t)\nabla_{z}\left(\frac{1}{2}\|y-Az\|_{2}^{2}\right)$ $\displaystyle=0,$ $\displaystyle z-\hat{z}_{0}(z_{t};\epsilon_{\theta})+\gamma(t)\nabla_{z}\left(\frac{1}{2}\|y-Az\|_{2}^{2}\right)$ $\displaystyle=0,$ $\displaystyle z-\hat{z}_{0}(z_{t};\epsilon_{\theta})+\gamma(t)\nabla_{z}\left(\frac{1}{2}\left(y^{T}y+z^{T}A^{T}Az-2y^{T}Az\right)\right)$ $\displaystyle=0,$ $\displaystyle z-\hat{z}_{0}(z_{t};\epsilon_{\theta})+\gamma(t)\left(A^{T}Az-A^{T}y\right)$ $\displaystyle=0,$ $\displaystyle\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]z-\left(\hat{z}_{0}(z_{t};\epsilon_{\theta})+\gamma(t)A^{T}y\right)$ $\displaystyle=0.$ $$

Solving this, we obtain the following expression for $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})$:

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}(\hat{z}_{0}(z_{t};\epsilon_{\theta})+\gamma(t)A^{T}y).$ $$

Note that the inverse of $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]$ exists as $A^{T}A\succ 0$ (from Assumption [2](#Thmassumption2)) and $\gamma(t)>0$, which ensures $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]\succ 0$.
Further, substituting the expression for $\hat{z}_{0}(z_{t};\epsilon_{\theta})$, we obtain

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon_{\theta})=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}\left[\frac{z_{t}-\sqrt{1-\bar{\alpha}_{t}}\epsilon_{\theta}(z_{t},t)}{\sqrt{\bar{\alpha}_{t}}}+\gamma(t)A^{T}y\right].$ $$

Given that $P_{\mathrm{data}}=\mathcal{N}(\mathbf{\mu},\mathbf{I}_{n})$, for the $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, we use the expression for the optimal denoiser $\epsilon^{*}(z_{t},t)$ (check Eq. [27](#A1.E27)) in place of $\epsilon_{\theta}(z_{t},t)$ to obtain

$$ $\displaystyle\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})$ $\displaystyle=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}\left[\frac{z_{t}+(1-\bar{\alpha}_{t})(\sqrt{\bar{\alpha}_{t}}\mu-z_{t})}{\sqrt{\bar{\alpha}_{t}}}+\gamma(t)A^{T}y\right],$ $\displaystyle\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})$ $\displaystyle=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}\left[\frac{z_{t}+\sqrt{\bar{\alpha}_{t}}\mu-z_{t}-\bar{\alpha}_{t}\sqrt{\bar{\alpha}_{t}}\mu+\bar{\alpha}_{t}z_{t}}{\sqrt{\bar{\alpha}_{t}}}+\gamma(t)A^{T}y\right].$ $$

This can be finally simplified to obtain the expression

$$ $\hat{z}_{0,\mathrm{pr}}(z_{t};\epsilon^{*})=[\mathbf{I}_{n}+\gamma(t)A^{T}A]^{-1}[\mu-\bar{\alpha}_{t}\mu+\sqrt{\bar{\alpha}_{t}}z_{t}+\gamma(t)A^{T}y].$ $$

∎

###### Lemma 3 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1].$ If $\ \bar{\alpha}_{t}<\bar{\alpha}_{t-1}$ and the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(t)>0\ \forall\ t\in[1,T]$, the spectral norm of the matrix $K_{t}$, $\|K_{t}\|_{2}$, with $K_{t}$ as defined in Eq. [16](#A1.E16), is less than 1.

###### Proof.

We want to show that

$$ $\displaystyle\|K_{t}\|_{2}=\left\|\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right]\right\|_{2}<1.$ $$

The spectral norm follows the triangle inequality. Therefore, after simplifying the expression with triangle inequality, we need to show

$$ $\displaystyle\left\|\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right\|_{2}+\left\|\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right\|_{2}<1.$ $$

We note that for $\gamma(t)>0$, $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]\succ 0$, and therefore $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\succ 0$. Similarly, $\mathbf{I}_{n}\succ 0$.

Further, we use the identities that if $M\succ 0$, then $\|M\|_{2}=\lambda_{\mathrm{max}}(M)$, $\|M^{-1}\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(M)}$, and $\|cM\|_{2}=|c|\|M\|_{2}$.

Therefore, $\|\mathbf{I}_{n}\|_{2}=1$, $\|\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}$. Further, note that $\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\geq 0$ and $\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\geq 0$. Substituting these, the inequality simplifies to

$$ $\displaystyle\frac{\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}<1.$ $$

Therefore, it is sufficient to show that

$$ $\displaystyle\frac{\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}<1-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}.$ $$

For any diffusion process with noise coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, where $\bar{\alpha}_{t}>\bar{\alpha}_{t-1}\ \forall\ t\in[1,T]$, Lemma [5](#Thmlemma5) shows that $\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\leq 1-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}$. Therefore, it is sufficient to show that $\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])>1$.

To proceed further, we use the Weyl’s inequality Horn & Johnson (), which states that for any two real symmetric matrices $P\in\mathbb{R}^{n\times n}$ and $Q\in\mathbb{R}^{n\times n}$, if the eigenvalues are represented as $\lambda_{\mathrm{max}}(P)=\lambda_{1}(P)>=\lambda_{2}(P)\dots>=\lambda_{n}(P)=\lambda_{\mathrm{min}}(P)$, and $\lambda_{\mathrm{max}}(Q)=\lambda_{1}(Q)>=\lambda_{2}(Q)\dots>=\lambda_{n}(Q)=\lambda_{\mathrm{min}}(Q)$, then we have the following inequality:

$$ $\lambda_{i}(P)+\lambda_{j}(Q)\leq\lambda_{i+j-n}(P+Q).$ (28) $$

For $i=j=n$, we have $\lambda_{\mathrm{min}}(P)+\lambda_{\mathrm{min}}(Q)\leq\lambda_{\mathrm{min}}(P+Q)$.

For $P=\mathbf{I}_{n}$ and $Q=\gamma(t)A^{T}A$ with $\gamma(t)>0$, this inequality can be exploited as both these matrices are real and symmetric.
Therefore, we have

$$ $\displaystyle\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])$ $\displaystyle\geq\lambda_{\mathrm{min}}(\mathbf{I}_{n})+\lambda_{\mathrm{min}}(\gamma(t)A^{T}A),$ (29) $\displaystyle\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])$ $\displaystyle\geq 1+\gamma(t)\lambda_{\mathrm{min}}(A^{T}A).$ (30) $$

Note that now it is sufficient to show $1+\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)>1$. For $\gamma(t)>0$, this inequality holds true as $\lambda_{\mathrm{min}}(A^{T}A)>0$ ($A^{T}A\succ 0$). Therefore,

$$ $\displaystyle\left\|\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right]\right\|_{2}<1.\$ $$

∎

###### Lemma 4 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]$. If $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}$ and the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(t)>0\ \forall\ t\in[1,T]$, $\|K_{1}\|_{2}$ with $K_{t}$ as defined in Eq. [16](#A1.E16) is given by

$$ $\displaystyle\|K_{1}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}.$ (31) $$

###### Proof.

We want to find an upper bound for

$$ $\displaystyle\|K_{t}\|_{2}=\left\|\left[\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right]\right\|_{2}.$ $$

Applying the triangle inequality for spectral norm, we get

$$ $\displaystyle\left\|K_{t}\right\|_{2}\leq\left\|\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right\|_{2}+\left\|\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\mathbf{I}_{n}\right\|_{2}.$ $$

We use the same simplifications shown in Lemma [3](#Thmlemma3) to obtain

$$ $\displaystyle\|K_{t}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}+\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}.$ $$

For $t=1$, we know that $\bar{\alpha}_{t-1}=\bar{\alpha}_{0}=1$. Therefore, we obtain

$$ $\displaystyle\|K_{1}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{\lambda_{\mathrm{min}}([\mathbf{I}_{n}+\gamma(1)A^{T}A])}.$ $$

Further, the denominator can be lower bounded using Weyl’s inequality, as shown in Eq. [30](#A1.E30). Therefore, we obtain

$$ $\displaystyle\|K_{1}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{\lambda_{\mathrm{min}}([\mathbf{I}_{n}+\gamma(1)A^{T}A])}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}.$ $$

Hence, we have shown that

$$ $\displaystyle\|K_{1}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{1}}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}.$ $$

∎

###### Lemma 5 .

For any $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]\ \forall t\in[1,T]$, if $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}$, then

$$ $\displaystyle\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}<1-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}.$ $$

###### Proof.

Squaring on both sides, we get

$$ $\displaystyle\bar{\alpha}_{t-1}\bar{\alpha}_{t}<1+(1-\bar{\alpha}_{t-1})(1-\bar{\alpha}_{t})-2\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}.$ $$

After further simplification, we have to show

$$ $\displaystyle\bar{\alpha}_{t-1}\bar{\alpha}_{t}<(1-\bar{\alpha}_{t})+(1-\bar{\alpha}_{t-1})+\bar{\alpha}_{t-1}\bar{\alpha}_{t}-2\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}},$ $\displaystyle 0<(1-\bar{\alpha}_{t})+(1-\bar{\alpha}_{t-1})-2\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}},$ $\displaystyle 0<(\sqrt{1-\bar{\alpha}_{t-1}}-\sqrt{1-\bar{\alpha}_{t}})^{2}.$ $$

Since $\bar{\alpha}_{t}\neq\bar{\alpha}_{t-1}$, we know that $\sqrt{1-\bar{\alpha}_{t-1}}\neq\sqrt{1-\bar{\alpha}_{t}}$. Therefore $(\sqrt{1-\bar{\alpha}_{t-1}}-\sqrt{1-\bar{\alpha}_{t}})^{2}>0$. Therefore, we conclude that

$$ $\displaystyle\sqrt{\bar{\alpha}_{t-1}}\sqrt{\bar{\alpha}_{t}}<1-\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}.$ $$

Note that this clearly holds for the edge case $t=1$, where we have $\sqrt{\bar{\alpha}_{1}}<1$, and for $t=T$, where we have $0<1-\sqrt{1-\bar{\alpha}_{T-1}}$. For the choices of $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, these clearly hold true.
∎

###### Lemma 6 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]\ \forall\ t\in[1,T]$. For the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$, $\|D_{t}\|_{2}$, with $D_{t}$ as defined in Eq. [19](#A1.E19), is less than 1.

###### Proof.

Note that the matrix $D_{t}$ is given by,

$$ $\displaystyle D_{t}=\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}A-\mathbf{I}_{n}.$ $$

Using the matrix inversion identity, $(AB)^{-1}=B^{-1}A^{-1}$, we rewrite $D_{t}$ as follows.

$$ $\displaystyle D_{t}=\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\left(A^{T}A\right)^{-1}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]\right]^{-1}-\mathbf{I}_{n}.$ $\displaystyle D_{t}=\sqrt{\bar{\alpha}_{t-1}}\left[\frac{\left(A^{T}A\right)^{-1}}{\gamma(t)}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]\right]^{-1}-\mathbf{I}_{n}.$ $\displaystyle D_{t}=\sqrt{\bar{\alpha}_{t-1}}\left[\frac{\left(A^{T}A\right)^{-1}}{\gamma(t)}+\mathbf{I}_{n}\right]^{-1}-\mathbf{I}_{n}.$ $$

We observe that the choice of $\gamma(t)$ is greater than 0. More precisely, $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$.
Now, if $\|-\frac{(A^{T}A)^{-1}}{\gamma(t)}\|_{2}<1$, then we can apply the Neumann’s series for matrix inversion, which states that if $\|M\|_{2}<1$, then

$$ $\displaystyle[\mathbf{I}_{n}-M]^{-1}=\sum_{i=0}^{\infty}M^{i}.$ (32) $$

First, note that $\|-A\|_{2}=\|A\|_{2}$. Therefore, $\|-\frac{(A^{T}A)^{-1}}{\gamma(t)}\|_{2}=\|\frac{(A^{T}A)^{-1}}{\gamma(t)}\|_{2}$ for $\gamma(t)>0$. From the theorem statement, $\gamma(t)>0.$

Additionally, we know that $\|\frac{(A^{T}A)^{-1}}{\gamma(t)}\|_{2}=\lambda_{\mathrm{max}}\left(\frac{(A^{T}A)^{-1}}{\gamma(t)}\right)=\frac{1}{\gamma(t)\lambda_{\mathrm{min}}\left(A^{T}A\right)}$.

Therefore, it is enough to show that $\frac{1}{\gamma(t)\lambda_{\mathrm{min}}\left(A^{T}A\right)}<1$ to apply the Neumann’s series.

However, we know that $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}((A^{T}A)^{-1})}$. Therefore, we observe that $\frac{1}{\gamma(t)\lambda_{\mathrm{min}}\left(A^{T}A\right)}<\frac{1}{2}<1$.

Thus, we have shown that $\|\frac{(A^{T}A)^{-1}}{\gamma(t)}\|_{2}<1$. Therefore, using Eq. [32](#A1.E32), we get

$$ $\displaystyle\left[\mathbf{I}_{n}-\left(-\frac{(A^{T}A)^{-1}}{\gamma(t)}\right)\right]^{-1}=\sum_{i=0}^{\infty}\left(\frac{(-A^{T}A)^{-1}}{\gamma(t)}\right)^{i}=\mathbf{I}_{n}+\sum_{i=1}^{\infty}\left(\frac{(-1)^{i}(A^{T}A)^{-i}}{\gamma(t)^{i}}\right).$ (33) $$

The last equality stems from the fact that for any matrix $M\in\mathbb{R}^{n\times n}$, $M^{0}=\mathbf{I}_{n}$. Substituting this expression for the second term in $D_{t}$, we get

$$ $\displaystyle D_{t}=\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left(\frac{(-1)^{i}(A^{T}A)^{-i}}{\gamma(t)^{i}}\right)\right)+\sqrt{\bar{\alpha}_{t-1}}\mathbf{I}_{n}-\mathbf{I}_{n}.$ $$

On further simplification, we have

$$ $\displaystyle D_{t}=\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left(\frac{(-1)^{i}(A^{T}A)^{-i}}{\gamma(t)^{i}}\right)\right)-\left(1-\sqrt{\bar{\alpha}_{t-1}}\right)\mathbf{I}_{n}.$ $$

Computing the spectral norm and using the triangle inequality, we get

$$ $\displaystyle\|D_{t}\|_{2}$ $\displaystyle=\left\|\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left(\frac{(-1)^{i}(A^{T}A)^{-i}}{\gamma(t)^{i}}\right)\right)-\left(1-\sqrt{\bar{\alpha}_{t-1}}\right)\mathbf{I}_{n}\right\|_{2},$ $\displaystyle\leq\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left\|\frac{(-1)^{i}(A^{T}A)^{-i}}{\gamma(t)^{i}}\right\|_{2}\right)+\left\|(1-\sqrt{\bar{\alpha}_{t-1}})\mathbf{I}_{n}\right\|_{2}.$ $$

The inequality arises from the triangle inequality for spectral norms. Note that each of the matrices within the summation is either positive definite or negative definite, and the spectral norms of all these matrices can be represented as $\left\|\frac{(A^{T}A)^{-i}}{\gamma(t)^{i}}\right\|_{2}$.
Therefore, we get

$$ $\displaystyle\|D_{t}\|_{2}\leq\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left\|\frac{(A^{T}A)^{-i}}{\gamma(t)^{i}}\right\|_{2}\right)+\left(1-\sqrt{\bar{\alpha}_{t-1}}\right).$ $$

Using the inequality $\|MN\|_{2}\leq\|M\|_{2}\|N\|_{2}$ multiple times, we get the following:

$$ $\displaystyle\left\|\frac{(A^{T}A)^{-i}}{\gamma(t)^{i}}\right\|_{2}\leq\frac{1}{\gamma(t)^{i}}\left(\left\|(A^{T}A)^{-1}\right\|_{2}\right)^{i}.$ $$

Additionally, for the above equation, we used $\|cM\|_{2}=|c|\|M\|_{2}$. Here, $c$ is $\gamma(t)$, which is greater than 0.
Since $A^{T}A\succ 0$, we have $\left\|(A^{T}A)^{-1}\right\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(A^{T}A)}$. Therefore, we have the following inequality:

$$ $\displaystyle\left\|\frac{(A^{T}A)^{-i}}{\gamma(t)^{i}}\right\|_{2}\leq\frac{1}{\gamma(t)^{i}}\left(\frac{1}{\lambda_{\mathrm{min}}(A^{T}A)}\right)^{i}=\frac{1}{\left(\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)\right)^{i}}.$ $$

Using this to upper bound $\|D_{t}\|_{2}$, we get

$$ $\displaystyle\|D_{t}\|_{2}\leq\sqrt{\bar{\alpha}_{t-1}}\left(\sum_{i=1}^{\infty}\left(\frac{1}{\left(\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)\right)^{i}}\right)\right)+\left(1-\sqrt{\bar{\alpha}_{t-1}}\right).$ $$

Finally, the summation of an infinite geometric series of the form $a+a^{2}+\dots$, where $a<1$ is $\frac{a}{1-a}$. Here, note that we have $\gamma(t)>\frac{1}{\lambda_{\mathrm{min}}(A^{T}A)}$. Therefore, $\frac{1}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)}<1$. Finally, we have

$$ $\displaystyle\sum_{i=1}^{\infty}\left(\frac{1}{\gamma(t)^{i}(\lambda_{\mathrm{min}}(A^{T}A))^{i}}\right)=\frac{\frac{1}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)}}{1-\frac{1}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)}}=\frac{1}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1}.$ $$

So, we obtain

$$ $\displaystyle\|D_{t}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{t-1}}}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1}+\left(1-\sqrt{\bar{\alpha}_{t-1}}\right).$ (34) $$

Now, for $\|D_{t}\|_{2}<1$, we need to show

$$ $\displaystyle\frac{\sqrt{\bar{\alpha}_{t-1}}}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1}+\left(1-\sqrt{\bar{\alpha}_{t-1}}\right)<1,\text{or}$ $\displaystyle\frac{\sqrt{\bar{\alpha}_{t-1}}}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1}<\sqrt{\bar{\alpha}_{t-1}}.$ $$

This simplifies to showing $\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1>1$, which is true if $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$. And, from the statement of the lemma, we know that $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$.

Therefore, we have shown that $\|D_{t}\|_{2}<1$ for $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$.
∎

###### Lemma 7 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]\ \forall\ t\in[0,T]$. For the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(1)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$, $\|D_{1}\|_{2}$, with $D_{t}$ as defined in Eq. [19](#A1.E19), is upper bounded by

$$ $\displaystyle\|D_{1}\|_{2}$ $\displaystyle\leq\begin{matrix}\frac{1}{\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)-1}.\end{matrix}$ $$

###### Proof.

Note that the matrix $D_{t}$ is given by,

$$ $\displaystyle D_{t}=\gamma(t)\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}A^{T}A-\mathbf{I}_{n}.$ $$

From Eq. [34](#A1.E34) in Lemma [6](#Thmlemma6), we know that if $\gamma(t)>\frac{1}{\lambda_{\mathrm{min}}(A^{T}A)}$,

$$ $\displaystyle\|D_{t}\|_{2}\leq\frac{\sqrt{\bar{\alpha}_{t-1}}}{\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)-1}+\left(1-\sqrt{\bar{\alpha}_{t-1}}\right).$ $$

From the lemma, we know that $\gamma(t)>\frac{2}{\lambda_{\mathrm{min}}(A^{T}A)}$. Therefore, we use Eq. [34](#A1.E34) and substitute for $t=1$ and $\bar{\alpha}_{0}=1$, we get

$$ $\displaystyle\|D_{1}\|_{2}\leq\frac{1}{\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)-1}.$ $$

∎

###### Lemma 8 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]$. If $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}\forall\ t\in[1,T]$ with the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(t)>0$, $\|E_{t}\|_{2}<1$ where $E_{t}$ is defined as in Eq. [17](#A1.E17).

###### Proof.

We know that the matrix $E_{t}$ is defined as

$$ $\displaystyle E_{t}=\left[(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right].$ $$

First, we use the identity $\|cM\|_{2}=|c|\|M\|_{2}$, where $c$ is any real number, we need to show

$$ $\displaystyle(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\left\|\left[\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right]\right\|_{2}<1.$ $$

Note that $(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\geq 0$. Further, for $\gamma(t)>0$, $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]\succ 0$, and therefore $\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\succ 0$.

We use the identity that for $M\succ 0,\|M^{-1}\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(M)}$.

Therefore, $\|\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\|_{2}=\frac{1}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}$. We use this expression to simplify the inequality as

$$ $\displaystyle\frac{(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}}{\lambda_{\mathrm{min}}(\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right])}<1.$ $$

We use to Weyl’s inequality (check Eq. [30](#A1.E30)) to lower bound the denominator and thereby upper bound the left side. Therefore, it is sufficient to show

$$ $\displaystyle\frac{(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}}{1+\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)}<1.$ $$

We observe that the numerator $(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}$ is always less than 1. However, we know that the denominator $1+\gamma(t)\lambda_{\mathrm{min}}(A^{T}A)$ is strictly greater than 1 for $\gamma(t)>0$ as $(A^{T}A)^{-1}$ exists and $\lambda_{\mathrm{min}}(A^{T}A)>0$. Therefore, the left side is always less than 1. This leads to

$$ $\displaystyle\left\|(1-\bar{\alpha}_{t})\sqrt{\bar{\alpha}_{t-1}}\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}\right\|_{2}<1.$ $$

∎

###### Lemma 9 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1],$ If $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}\forall\ t\in[1,T]$ with the penalty coefficients from Algorithm [1](#alg1) given by $\gamma(t)>0$, $\|E_{1}\|_{2}$, with $E_{t}$ defined as in Eq. [17](#A1.E17), is upper bounded by

$$ $\displaystyle\sigma_{\mathrm{max}}(E_{1})\leq\frac{1-\bar{\alpha}_{1}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}.$ $$

###### Proof.

We know that $E_{t}$ is given by

$$ $\displaystyle E_{t}=\sqrt{\bar{\alpha}_{t-1}}(1-\bar{\alpha}_{t})\left[\mathbf{I}_{n}+\gamma(t)A^{T}A\right]^{-1}.$ $$

We first substitute for $t=1$ and $\sqrt{\bar{\alpha}}_{0}=1$

$$ $\displaystyle E_{1}=(1-\bar{\alpha}_{1})[\mathbf{I}_{n}+\gamma(1)A^{T}A]^{-1}.$ $$

We use the identity $\|cM\|_{2}=|c|\|M\|_{2}$, where $c$ is any real number, to get

$$ $\displaystyle\left\|E_{1}\right\|_{2}=(1-\bar{\alpha}_{1})\left\|[\mathbf{I}_{n}+\gamma(1)A^{T}A]^{-1}\right\|_{2}.$ $$

Here, note that $(1-\bar{\alpha}_{1})\geq 0$. Similar to Lemma [8](#Thmlemma8), we can rewrite the spectral norm as

$$ $\displaystyle\left\|E_{1}\right\|_{2}=\frac{1-\bar{\alpha}_{1}}{\lambda_{\mathrm{min}}([\mathbf{I}_{n}+\gamma(1)A^{T}A])}.$ $$

Again, using Weyl’s inequality and performing similar modifications as in Lemma [8](#Thmlemma8), we obtain the following upper bound for the spectral norm

$$ $\displaystyle\left\|E_{1}\right\|_{2}\leq\frac{1-\bar{\alpha}_{1}}{1+\gamma(1)\lambda_{\mathrm{min}}(A^{T}A)}.$ $$

∎

###### Lemma 10 .

Suppose Assumption [2](#Thmassumption2) holds. Consider a $T$-step diffusion process with coefficients $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$ such that $\bar{\alpha}_{0}=1$, $\bar{\alpha}_{T}=0$, $\bar{\alpha}_{t}\in[0,1]$. If $\bar{\alpha}_{t}<\bar{\alpha}_{t-1}\ \forall t\in[0,T]$, $\|F_{t}\|_{2}$, with $F_{t}$ as defined in Eq. [18](#A1.E18), is less than 1. Additionally, $\|F_{1}\|_{2}$ is 0.

###### Proof.

Note that $F_{t}$ is given by the expression,

$$ $\displaystyle F_{t}=\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\mathbf{I}_{n}.$ $$

First, we use the identity $\|cM\|_{2}=|c|\|M\|_{2}$, where $c$ is any real number. Therefore, we need to show

$$ $\displaystyle\|F_{t}\|_{2}=\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}\left\|\mathbf{I}_{n}\right\|_{2}<1.$ $$

For the given conditions on $\bar{\alpha}_{0},\dots,\bar{\alpha}_{T}$, we observe that at least one of the terms in $\sqrt{1-\bar{\alpha}_{t-1}}\sqrt{1-\bar{\alpha}_{t}}\sqrt{\bar{\alpha}_{t}}$ is always less than 1. Therefore $\|F_{t}\|_{2}<1$. And, since $\bar{\alpha}_{0}=1$, for $F_{1}$, we have $\sqrt{1-\bar{\alpha}_{0}}=0$. Therefore, $F_{1}$ is a null matrix and $\|F_{1}\|_{2}=0$.
∎

###### Lemma 1 .

###### Proof.

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

###### Lemma 7 .

###### Proof.

###### Lemma 8 .

###### Proof.

###### Lemma 9 .

###### Proof.

###### Lemma 10 .

###### Proof.

## Appendix B Metrics

For the FTSD and J-FTSD metrics, we train the time series and condition encoders using the procedure given in Narasimhan et al. (). For FTSD, we only train the time series encoder using supervised contrastive loss to maximize the similarity of time series chunks that belong to the same sample. For J-FTSD, we perform contrastive learning training in a CLIP-like manner to maximize the similarity between time series and corresponding paired metadata, as explained in Narasimhan et al. (). We use Informer models as the encoders. Additionally, just as in the case of (Paul et al., ; Narasimhan et al., ), we observe that the approaches corresponding to the lowest values of FD metrics have the lowest TSTR and DTW scores, and the highest SSIM scores. This further validates the correctness of the FTSD and J-FTSD metrics used for evaluation.

We sourced the implementations of DTW and SSIM from the public domain. For SSIM, we used 1D uniform filters from SciPY Virtanen et al. (). We set the values of $C_{1}$ and $C_{2}$ to $1^{-4}$ and $9^{-4}$.

For the constraint violation magnitude, we computed the violation for each constraint, excluding the allowable constraint violation budget.

## Appendix C Datasets

We compared CPS against the existing baselines for six settings - air quality (Chen, ), air quality conditional (Chen, ), traffic (Hogue, ), traffic conditional (Hogue, ), stocks Yoon et al. (), and waveforms. The training and testing splits for the air quality and traffic datasets are taken from Narasimhan et al. (). We additionally evaluated the constrained generation approaches on the stocks and the waveforms datasets. We used the preprocessing scripts provided by Yoon et al. () for the stocks dataset. The waveforms dataset was synthetically generated. We generated $~{}64,000$ sinusoidal waveforms of varying amplitudes, phases, and frequencies. The amplitude varies from $0.1$ to $1.0$. The phase varies from 0 0 to $2\pi$. The frequency limits were chosen based on the Nyquist criterion. The generators and the GAN models were trained on this dataset. However, for the TSTR metrics, we created a subset of this dataset with $~{}16,000$ samples. All the datasets except the waveforms dataset were standard normalized.

The Air Quality dataset is a multivariate dataset with six channels. The total number of train, val, and test samples are 12166, 1537, and 1525, respectively.

The Traffic dataset is univariate. The total train, val, and test samples are 1604, 200, and 201, respectively.

The Stocks dataset is a multivariate dataset with six channels. The total train, val, and test samples are 2871, 358, and 360, respectively.

The truncated form of the waveforms dataset used for evaluation consists of 13320, 1665, and 1665 train, val, and test samples, respectively.

The horizon length used for all datasets is 96.

## Appendix D Implementation

In this section, we will describe the implementation details of our approach, each baseline, trained models, metrics, etc.

### D.1 Diffusion Model Architecture

We utilize the Time Weaver-CSDI denoiser for all the diffusion models used in this work. The training hyperparameters and the model parameters are precisely the same as indicated in (Narasimhan et al., ). The total number of residual layers is 10 for all the experiments. Further, we used 200 denoising steps with a linear noise schedule for the diffusion process. All the baselines and CPS use the same base diffusion model with the Time Weaver-CSDI denoiser backbone.

We use 256 channels in each residual layer, with 16-dimensional vectors representing each time series channel. The diffusion time step input embedding is a 256-dimensional vector. Further, the metadata encoder has an embedding size of 256 for the conditional case. The metadata encoder has two attention layers with eight attention heads. All our experiments use a learning rate of $10^{-4}$. Our training procedure and the hyperparameters are precisely the same as the values in Narasimhan et al. ().

### D.2 Constrained Posterior Sampling Implementation

For the CPS implementation, we use CVXPY Diamond & Boyd (). We first implement the constraint violation function with the violation threshold set to $0.005$ for all the constraints except the bounds like argmax, argmin, OHLC, and the trend constraint. For example, consider the mean constraint. The constraint violation function for this constraint is implemented as $\max\left(\left|\frac{1}{L}\left(\sum_{u=1}^{L}c(u)\right)-\mu_{c}\right|-0.005,0\right)$, where $L$ is the time series horizon. We do not provide the constraint violation threshold for the bounds. Though the allowable constraint violation threshold is $0.01$, we performed the projection step with a constraint violation threshold of $0.005$ to ensure that the sample strictly lies within the constraint set. We use the same choice of $\gamma(t)\ \forall t\in[1,T]$ as described in Sec. [3](#S3). However, we clip the value of $\gamma(t)$ to $100,000$ after certain denoising steps, as the CVXPY solvers cannot handle extremely high values of $\gamma(t)$. We note that this clipping usually occurs after 150 denoising steps.

### D.3 Baseline Implementation

This section will explain all the details about the baseline implementations. Specifically, we use two baselines - Constrained Optimization Problem (COP) and Guided DiffTime. We note that both approaches were proposed in (Coletta et al., ). However, the implementation of these approaches is not publicly available. Based on the details provided in (Coletta et al., ), we have implemented the baselines for comparison against CPS.

#### D.3.1 Constrained Optimization Problem Implementation

The Constrained Optimization Problem, COP, has two variants. These are referred to as COP and COP-FineTuning, respectively. In COP, we perturb a randomly selected sample from the training and validation datasets. In COP-FineTuning, we perturb the sample generated from the Time Weaver-CSDI diffusion model.

Note that (Coletta et al., ) suggests to extract statistical features to be imposed as distributional constraints. For example, Coletta et al. () suggests extracting autocorrelation features for the stocks dataset. However, since it is practically impossible to list all the statistical features for each dataset to obtain the distributional constraints, Coletta et al. () suggests the use of the critic function from a Wasserstein GAN (Arjovsky et al., ). The details of the GAN training are summarized below.

COP has two objectives - maximize the $l_{2}$ distance from a randomly selected sample from the training and maximize the critic value from a Wasserstein GAN.

Similarly, COP FineTuning has two objectives - minimize the $l_{2}$ distance from a generated sample and maximize the critic value from a Wasserstein GAN.

We optimize for these objectives while ensuring constraint satisfaction.

As suggested in (Coletta et al., ), we use the SLSQP solver from SciPy Virtanen et al. (). Unlike (Coletta et al., ), which performs piecewise optimization, we note that all the constraints used in our work are global. Therefore, piecewise optimization is very suboptimal. For example, it is suboptimal to break a time series into chunks and perform optimization for each piece when the goal is to generate a sample with a specific mean value. This is also pointed out in (Coletta et al., ). Therefore, we perform COP for the whole time series at once. We consider two budgets - 0.005 and 0.01. This is similar to Coletta et al. (). However, unlike their approach, we stop with 0.01 as the allowable constraint violation in our case is 0.01 for all methods.

We used a weight of 0.1 for the critic’s objective. We noticed that for different values (1.0,0.1,0.01) of this weight, there was very little change in the DTW and the SSIM metrics.

#### D.3.2 Critic Function Implementation

Coletta et al. () suggest using the critic function in a Wasserstein GAN (Arjovsky et al., ) to enforce realism in the COP approach. Therefore, we used the WaveGAN (Donahue et al., ) implementation from (Alcaraz & Strodthoff, ). The implementation from (Alcaraz & Strodthoff, ) has the gradient penalty loss, an improved training procedure to enforce the required Lipschitz continuity for the critic function. Additionally, the WaveGAN training with gradient penalty has been implemented (Alcaraz & Strodthoff, ) for generating time series samples for the ECG domain. Therefore, we use their implementation to obtain the critic function for the COP baseline. The number of parameters is adjusted such that the diffusion model and the GAN model have a comparable number of parameters.

Similar to the diffusion model, we used the same architecture and training hyperparameters for all the datasets and experimental settings. Specifically, we trained the WaveGAN model with a learning rate of $10^{-4}$ for all the datasets. The input to the generator is a 48-dimensional random vector. Additionally, we ensured that the total number of parameters was equally distributed between the generator and the discriminator to prevent either of the models from overpowering the other.

#### D.3.3 Guided DiffTime Implementation

We use the same Time Weaver-CSDI denoiser as in the case of CPS. For the guidance weight, we experimented with the following weights - $\left(0.00001,0.0001,0.001,0.01,0.1,1.0\right)$. We chose the best guidance weight based on the constraint violation rate. Note that we used the same guidance weight for all individual constraints. Using PyTorch, we implemented all the constraints mentioned in Sec. [4](#S4). Additionally, we augmented the Guided DiffTime approach with the DiffTime algorithm for fixed values. In other words, after each step of denoising followed by guidance update, we enforced the fixed value constraints, as specified in (Coletta et al., ). This applies to the values at argmax, argmin, 1, 24, 48, 72, and 96 timestamps.

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 3
presentation: 4
questions: see limitations above
soundness: 3
strengths: - The paper provides a solid theoretical basis for CPS, including proofs of convergence and constraint satisfaction.
- Experiments on real-world datasets demonstrate CPS's ability to generate high-quality and constraint-compliant time series.
- It seems like that the framework can handle a large number of constraints, and all at sampling time.
summary: The paper introduces Constrained Posterior Sampling, a diffusion-based method for generating time series that adhere to hard constraints. The approach is theoretically grounded and the empirical evaluations demonstrate CPS's effectiveness in producing high-quality, constraint-compliant time series. 
However, the novelty of CPS is somewhat limited due to its conceptual similarities with existing methods (see below for extended discussion) .
weaknesses: - The key weakness of this approach is novelty. CPS shares significant conceptual similarities with existing methods. 
In particular :
[A] also uses a projection step within diffusion processes to enforce constraints, and consider arbitrary constraint sets, while providing theoretical guarantees for convex constraint sets. The algorithm proposed in [A] is also very similar to the proposed one. 
[B] implements hard constraints on the outputs of autoencoders for linear constraints, as those considered in this paper. 
[C] uses “mirror mappings” which apply to convex constraint sets. 
[D] further extend the classes of constraints can be handled handling constraints represented by convex polytopes. 

- Given the above, the paper lacks a thorough comparison with existing methods and this omission makes it difficult to assess CPS's relative advantages or limitations.

- The projection step in CPS can be computationally intensive, especially for high-dimensional constraints, how can this be handled? 

[A] Jacob K Christopher, Stephen Baek, Ferdinando Fioretto. Constrained Synthesis with Projected Diffusion Models. NeurIPS 2024 - ArXiv: https://arxiv.org/abs/2402.03559v1 (first appeared on ArXiv on Feb 2024)

[B] Thomas Frerix, Matthias Nießner, and Daniel Cremers. Homogeneous linear inequality constraints for neural network activations. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 748–749, 2020.

[C] Guan-Horng Liu, Tianrong Chen, Evangelos Theodorou, and Molei Tao. Mirror diffusion models for constrained and watermarked generation. NeurIPS 2023. ((first appeared on ArXiv on Oct 2023)

[D] Nic Fishman, Leo Klarner, Emile Mathieu, Michael Hutchinson, and Valentin De Bortoli. Metropolis sampling for constrained diffusion models. NeurIPS 2023. ((first appeared on ArXiv on July 2023)

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 4
questions: •	Even though they propose new method, but the reviewer think the method is ot novel. It poses regularize similar to (Coletta et al., 2024), but with convex optimization strategy.

•	The author mentioned that the guided gradient is not enough for restricting the constraints. Could you clarify the main differences of two methods? In terms of convex optimization with global optimum, two of them should be same with proper set of hyperparameters.

•	Furthermore, the authors mentioned about the constraints sets as convex. However, the objective $||z-\hat{z}_{0,pr}||^2$ that the author proposed might not be optimum for non-convex optimization.

•  The comparison includes only (Coletta et al., 2024). Could you consider additional comparison methods, such as the instance-aware guidance strategy in [1]?

•	In Theorem 2, why does $k$ approach infinity as $T$ goes to infinity? More explanation on the convergence analysis would be helpful.

•	Please correct minor errors. Such as generation(Yoon et al.) in line 140 and Recently (Coletta ~) -> use citet.

[1] DIFFUSION-TS: INTERPRETABLE DIFFUSION FOR GENERAL TIME SERIES GENERATION, ICLR 2024.
soundness: 2
strengths: • The paper investigate the very important and practical issues in time series generation, which encompasses generating particular examples with specific condition. This is the unique properties in time series generation compared to vision tasks.

• The authors propose a convex optimization based methods which is much more effective and simple for time series than existing methods.

• The paper conducts experiments on various datasets, including well-known ones, demonstrating the effectiveness of the constrained synthesis.

• The authors establish a theoretical background of their methods.
summary: This paper investigates generating time series data samples with specific constraints. The authors propose a training-free constrained posterior sampling method based on $z-\hat{z}_{0,pr}$ optimization problem within objective function. Furthermore, they theoretically justify the use of imposing condition as a conditional sampling. They empirically prove their idea across various settings, including different setups and datasets.
weaknesses: Please refer to the Questions section.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: - In Section 3, should $\hat{z}_0$, and various related symbols, have the subscript $t$? If not, why is 0 used there?
- What does the form of the distribution that CPS samples from, given in Theorem 1, tell us?
- How many constraints are imposed in Table 1 and Figure 6?
- How badly are the generated time series in Figure 6 violating the constraints?
- How is it possible that CPS is effectively copying samples (see traffic conditional and stocks datasets in Figure 6) from the test set that the generator has presumably not seen during training? Are these time series really characterised so precisely by a relatively small number of numbers?
soundness: 4
strengths: The paper is very well-written, easy to read, and makes clear arguments. The proposed method CPS is simple, and seems easy to apply. It may even be applicable to other domains, like tabular data, where meaningful constraints can be specified as simple inequalities. The results are very good, with CPS outperforming baseline methods in almost all experiments.
summary: The paper proposes a method called CPS to sample time series from a diffusion model with given constraints. The trick is projecting the intermediate values in the diffusion model sampling process towards the constraint set. In later iterations of the sampling process, the diffusion model corrects the projected time series so it stays close to the data generating distribution. The paper derives some theoretical properties of their sampler, including a convergence theorem in the case of Gaussian data. The paper also contains experiments on several real and toy datasets, where CPS outperforms baselines.
weaknesses: I'm not sure how meaningful the similarity metrics comparing the generated sample and the test sample that satisfies the all constraints is, since an optimal generator may not get the best possible score. This is because any finite set of constraints is likely to leave some freedom to generate different outputs, which would not be identical to the reference sample from the test set. This is shown in Table 1, where none of COP, COP-FT and CPS violate any constraints on several datasets, but have different values for the DTW and SSIM metrics.

A comparison with Loss DiffTime would still be interesting, since it would show how beneficial training specifically for a given set of constraints is.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: **Q1** Page 5, Lines 240-244 "Our key intuition is to design $\gamma(t)$ as a strictly **decreasing** function of $t$ that takes small values for the initial denoising steps (t close to T) and tends to $\infty$ for the final denoising steps". Does the function form of $gamma(t)$ affect the result? Is it possible to select others?

**Q2** The main contribution is the adaptive projection step that solves the penalized optimization task. From my understanding, this is more like a soft constraint and it becomes a hard one when $t=1$. However, in practice, how to choose $T$?

**Q3** In every denoising step, the protection step solves a regularized optimization task. However, the optimization problem might be difficult if some non-linear constraints are chosen. Is it guaranteed that the optimal or local optimal solution can be obtained?

**Q4** In Theorem 2, what is design parameter k? It seems missing in the algorithm.
soundness: 4
strengths: **1** This paper has extensive experiments to show the effectiveness of their method.

**2** A theoretical result is provided for validating the effectiveness under a linear constraint case. 

**3** Proper baselines are chosen for comparison and comprehensive literature review is provided.

**4** The proposed method basically provides some interesting ideas of handling constrained generation in diffusion models. I think the method can be extended to other scenarios.
summary: In this paper, the authors proposed a diffusion-based sampling algorithm - Constrained Posterior Sampling (CPS), to generate time series data meeting certain constraints. The basic idea is solving a regularized optimization task in every denoising step and the penalty term  will be strengthened such that the resulting data fall within the constrained set. Overall, I think this paper is well-written, has clear presentation and motivation. I will consider raising my score if the authors handle my questions properly.
weaknesses: **1** In each step, a constrained optimization task needs to be solved, which can be computationally intensive. Additionally, the current framework is more like as a general framework without imposing specific restrictions on the constraints. It would be beneficial to develop a universal algorithm that aims to minimize the constrained optimization for a general class of constraints, even if the optimal solution is not achieved.

**2**  The fidelity of generated data is not properly evaluated in this paper. I suggest using classification method to for fidelity evaluation. That is using a classifier to discriminate real and synthetic data and use the testing error as fidelity metric.

Meta-review: This paper considers the constrained sampling problem from a diffusion model. It has two orthogonal components, the development of the constrained posterior sampling (CPS) algorithm and its application to time series. The proposed CPS algorithm leverages proximal optimization algorithms with a carefully designed varying objective to gradually impose the hard constraints. One major weakness pointed out by the reviewers is the lack of thorough comparisons to existing methods. The authors have added comparison experiments to several methods including PDM, Diffusion-TS and Loss-DiffTime. Still, some baselines are missing. Constrained sampling can be viewed as a special case of conditional sampling and inverse problem [1]. For instance, Stochastic Control Guidance (SCG) [2] can also be used to achieve constrained sampling. In addition to experimental comparison, the paper can also benefit from a detailed comparison in methodology level to highlight its contribution. Finally, the position of the paper can be improved. It seems CPS is not restricted to time series generation. If it has superior performance in general, experiments should be conducted in other domains as well to make the paper stronger. If it only works for time series, detailed discussion should be included to explain the reasons. Overall, this is a borderline paper. The authors are encouraged to prepare a comprehensive revision based on these comments. 
[1] A Survey on Diffusion Models for Inverse Problems
[2] Symbolic Music Generation with Non-Differentiable Rule Guided Diffusion

Decision: Reject