Learning Cocoercive Conservative Denoisers via Helmholtz Decomposition for Poisson Imaging Inverse Problems

Dear reviewer aRzT,

We sincerely appreciate your thoughtful review and constructive feedback. Your comments are highly valuable, and we have carefully considered each point in revising the paper. Our detailed responses are provided below.

Responses to weaknesses:

1.) For most of the theoretical development in the paper, it is assumed that $D_\sigma$ is a CoCo operator and guarantees for the convergence of CoCo-ADMM and COCO-PGD are derived. However, the training for $D_\sigma$ only encourages CoCo nature, doesn’t enforce it. It is possible to check the $\gamma$-cocoerciveness via the spectral norm $\|2\gamma J-I\|_*$ during training—however, it is not possible to check if the denoiser is “conservative enough”. The guarantees require conservativeness, however—so it is unclear, if even though the norm $\|J-J^\mathrm{T}\|_*$ is close to zero, if the guarantee holds enough for it to be useful in reality. A natural next question to ask would be—how much could we relax the conservativeness to still ensure good enough convergence guarantees / proximal nature?

We will use the word encourage rather than enforce, because that spectral regularization can not guarantee the properties.

In the experiments, we found that, if the denoiser is not conservative enough, the convergence can still be stable with fine-tuned parameters. This may be because that the Hamiltonian part of the denoiser is useless (and therefore sometimes harmless) for the denoising performance.

To the best of our knowledge, there is no way to quantify how a non-conservative denoiser can still be good enough for being a proximal operator.

We will include the discussions in the revised paper.

2.) In the evaluation on line 270, the paper claims that the PSNR values in Table 1 validate that “cocoercive and conservative properties are less restrictive for deep denoiser”. I’m not sure why this performance improvement can necessarily be thought of as an indication of lesser restriction—the performance improvement could also have been because of the different loss function.

We agree with your comment. We will discuss carefully about the performance improvements. The performance improvement could also have been because of the different loss function.

3.) Appendix 13 claims that this is useful for real-time applications, even though 256x256 images take 5-10 seconds to process—I’m not sure if this claim makes sense.

We will remove the `real-time' expression in the revised paper.

Please let us know if you have more suggestions. We are very willing to revise the paper based on your comments to further improve its quality.

Dear reviewer G4Uk,

We sincerely appreciate your thoughtful review and constructive feedback. Your comments are highly valuable, and we have carefully considered each point in revising the paper. Our detailed responses are provided below.

Responses to weaknesses:

Thank you for this valuable comments. Please note that, we fully agree with your suggestion that ``I highly recommend writing the entire paper under this assumption and simplifying all parts that can be simplified.'' We have revised the paper carefully, such that all settings are in a finite-dimensional space.

But now, please allow us to respond to your questions. As answered below, most of the theoretical results in the paper still $**hold**$ in an infinite-dimensional setting.

1.1) To my mind, the paper jumps a little between great mathematical generality (e.g., talking about Frechet differentials, suggesting that $V$ could be infinite dimensional) and implicitly making assumptions without introducing them, e.g. the Helmholtz decompositions needs a smoothness assumption, any time the transpose symbol is used, I assume we are in a finite dimensional setting or is this denoting the adjoint operator?

Through out all theoretical parts in the paper, we were assuming a general real Hilbert space, which might be infinite-dimensional. We will include the smooth assumption for the Helmholtz decomposition. The transpose symbol means the adjoint operator in a Hilbert space $V$. When the space is finite dimensional like $\mathbb{R}^n$, the adjoint operator is the transpose.

When comes to training, we will have to consider a finite dimensional space of course. So we agree that the paper should be written in a finite dimensional setting.

1.2) Is the spectral norm the operator norm in the infinite-dimensional setting? If the convergence analysis holds in the finite-dimensional case only, I highly recommend writing the entire paper under this assumption and simplifying all parts that can be simplified.

In an infinite-dimensional setting, $\|\cdot\|_*$ is the operator norm. When the space is finite-dimensional, it is just the spectral norm, that is, the largest singular value.

1.3) Do all convergence results really hold in the infinite-dimensional setting as well?

Thank you for this valuable comments. We have revised the paper, such that all results are included in a finite dimensional space.

1.4) What assumptions do we need to make about $K$ in this case?

When $V$ is just a general real Hilbert space, $K$ is a bounded linear operator, that is $K\in \mathscr{B}[V]$. When considering a finite dimensional space, this is just a matrix.

1.4) Moreover, the notion of the subdifferential and its Lipschitz continuity (Theorem 2.2.) should be introduced.

We will include the subdifferential and its Lipschitz continuity in the revised version. Thank you so much for such a detailed comment.

2.) The denoiser is conservative if the Jacobian is symmetric, i.e., $J=J^\mathrm{T}$. The authors propose to regularize $\|J-J^\mathrm{T}\|_*$, but it seems like any norm could be used. This is not discussed in the paper.

Special thanks for this valuable feedback! We choose to use the spectral norm, mainly because that it is has low computational cost. Due to the power iterative method, one can calculate the spectral norm very easily, without actually calculate the Jacobian, see Appendix A.5 for details. For other norms, we do not know whether there exist similiar algorithms to calculate them. We will include this discussion in the revised paper.

3.) In practice, the necessary assumptions are implemented via soft penalties only, such that all provable statements are lost. The demonstration that the assumptions are met is demonstrated "over 100 patches". Rather than being met on average, I recommend reporting the highest norm after actively optimizing for an $x$ at which a high norm is attained (similar to an adversarial example).

Now the highest norm for $\|J-J^\mathrm{T}\|_*$ is reported below. Please note that the highest norm for cocoercivity has already been reported in Table 1, and thus omitted here.

This table will be included in the revised paper upon acceptance.

4.) The penalties require the computation and subsequent differentiation of the spectral norm of the Jacobian of the denoiser. Given that the denoiser is high-dimensional (mapping images to images), the Jacobian is extremely high-dimensional. While the appendix provides details about the computation, I think it is important to include a short statement on the computational expense of the training of the denoiser.

Thank you for this comment. Since the rebuttal window does not allow extra PDF or images, we report by words the computational cost. We first pre-train a DRUNet without spectral regularization. After that, we post-train the DRUNet with CoCo regularization. This is for saving training time.

Compared with the training for DRUNet without any spectral regularization, for each epoch, the computational time for training CoCo is about $12$ times larger than training DRUNet without spectral regularizations. For the GPU memory cost, DRUNet is about $2080$ MB, and CoCo is about $9360$ MB.

5.) Given the small differences in performance of the different methods, it is important to confirm that all training procedures are identical and to possibly report standard deviations.

We will state in the revised paper, that all training procedures for all denoisers are identical. In the revised paper, we will also report standard deviations with different random seed for denoising performance comparison.

Responses to questions:

Please address the question of the setting (finite vs. infinite-dimensional) of (1) under weaknesses and possibly adapt the paper accordingly.

Please refer to the responses to the weaknesses 1.1.

Please answer the question why you chose the spectral norm to enforce symmetry of $J$.

Please refer to the responses to the weaknesses 2.

Do the values in table 1 (right) change if one optimizes for a point $x$ at which the norms are high?

Unfortunately, yes. Thus we agree with your suggestions, and report the maximum norms, see the responses to the weaknesses 3.

How large is the computational overhead in training the denoisers with spectral penalties involving its Jacobian?

Unfortunately, pretty large: the computational cost for CoCo is 12 times larger than the training for DRUNet, see the responses to the weaknesses 4. As a result, we pre-trained a DRUNet first, and then post-train it with CoCo regularizations to save the training time.

Have all (competing) methods been fine-tuned equally? Did you consider multiple runs with different random seeds?

We have tried our best to fine tune each method, such that they reach the best PSNR values. For the random seed, for a fair comparison, we set the same random seed for all methods.

Responses to limitations:

Thanks for your suggestions on the limitations. We will address them in Section 6 in the revised paper.

Please let us know if you have more suggestions to polish the paper and raise the rating. We are very willing to revise the paper based on your comments to further improve its quality.

Thanks for the detailed response! I will increase my rating to an accept.

Just to make sure my comment regarding the highest |J-J^T| was clear: Since the Jacobian depends on the point where you evaluate it, I was thinking about approximately solving \max_x |J(x)-J^\mathrm{T}(x)| (not just sampling x from a batch). Is this what you did for the new table? If so, how did you approximate the solution of the maximization problem?

Dear reviewer UHTR,

Thank you so much for your valuable feedback and insightful comments. We sincerely appreciate your encouraging words and constructive critique; both are extremely valuable. We have given careful consideration to your points and are focused on refining the paper. Our detailed responses are provided below as part of our pursuit of excellence.

Responses to weaknesses:

1.) While the paper proposes cocoerciveness as a relaxed alternative to non-expansiveness, it is unclear why this particular relaxation is theoretically sufficient and practically better. What exactly fails when the operator is not cocoercive, and how tight is this condition?

Special thanks for this comment! This question is very valuable, and we will discuss this not only here, but also in the revised paper.

There basically two ways to use a prior in the context of PnP. If we use it in a forward way, for example to use it in RED, to the best of our knowledge, the weakest assumption is pseudo-contractiveness as in [71]. But if we use it in a backward way, for example this paper uses it as a prox, the weakest assumption (in terms of a weakly convex prior) has been proven to be `Lipschitz + be the gradient of some convex potential' in [22]. [22] established the equivalence between 'the denoiser is a proximal operator of some weakly convex prior', and 'the denoiser is Lipschitz and is the gradient of some convex potential'. Please note that, 'the denoiser is Lipschitz and is the gradient of some convex potential' can be equivalently rewritten as 'the denoiser is conservative and cocoercive': cocoerciveness ensures that the denoiser is Lipschitz and monotone; combined with conservativeness, the denoiser is then a gradient of a convex potential. The gradient of a convex potential is a special case of monotone operators.

The proof for this equivalence was not included in the paper. Currently, in Theorem A.5, we have proved that if the denoiser is conservative and cocoercive, it is a prox of some weakly convex prior. We will include the equivalence in the revised version.

What exactly fails when the operator is not cocoercive, and how tight is this condition?

As stated above, if we want the denoiser to be a prox of some weakly convex prior, the cocoercive and conservativeness are required. The condition is tight, and actually sufficient and necessary.

2.) The theoretical guarantee that the CoCo denoiser is a proximal operator relies on both cocoerciveness and exact conservativeness (i.e., symmetry of the Jacobian). In practice, both are only enforced softly during training. Does this soft enforcement suffice for convergence? What happens when these conditions are only approximately satisfied?

We softly enforce cocoerciveness and conservativeness by spectral regularization technique. The established convergence requires the two properties. In experiments, we found that, if the denoiser is not conservative enough, the convergence can still be stable with fine-tuned parameters. This may be because that the Hamiltonian part of the denoiser is useless (and therefore sometimes harmless) for the denoising performance. We will include this in the revised paper.

3.) The use of Helmholtz decomposition to interpret denoiser geometry is elegant, but relies on assumptions that may not hold in practice. Is the decomposition stable across different architectures and noise levels? How reliable is the connection between the anti-symmetric part and “useless” Hamiltonian directions in high-dimensional deep networks?

The decomposition is stable with different architechtures and noise levels, as long as the network output is differentiable with respect to the input, such that $J$ is computable. When the dimension is high, the connection between the anti-symmetric part and “useless” Hamiltonian directions is actually illustrated briefly in the equation blow Line 172. This means that the inner product of $D(y)-x$ and $\nabla_y \mathcal{H}(y)$ is zero, and thus the two vector is orthogonal. This further means that, if we consider the continuous dynamics, $\dot{y}=-\nabla_y \mathcal{H}(y)$, $\mathcal{H}(y(t))\equiv Constant$. That is, the dynamics preserve the Hamiltonian functional, if the denoiser only has Hamiltonian part. Since the Hamiltonian functional is a typical denoising loss, this means that a Hamiltonian field does not contribute to the denoising.

We will include more discussion on high dimensional cases in the revised paper.

4.) The convergence analysis (e.g., in Theorem 4.1) introduces a parameter $t<t_0(\gamma)$, where $t_0$ is the root of a cubic equation. However, the paper does not give much intuition on what range of $t$ works best in practice, or how sensitive the algorithm is to the choice of $\gamma$ and $t$. A more detailed discussion or ablation would help.

Since the rebuttal window does not allow PDF or inserting images, we only report the PSNR here. We have conducted oblation study with different $\gamma$ and $t$. We consider deblurring for the image `0005.png' from CBSD68 and kernel 8 from Levin's dataset. We set the Poisson noise level to be $50$. We run the algorithm for 500 iterations to ensure the convergence.

When $\gamma=0.25$, $t_0(\gamma)=\frac{1}{3}$. Thus $t<\frac{1}{3}$.

When $\gamma=0.5$, $t_0(\gamma)\approx 0.3761$. Thus $t\le 0.3761$.

From the two tables, it is clear that when $t$ is not too small, the PnP performance is not sensitive to the choice of $t$. However, different $\gamma$ has significant influence on the PSNR performance. This is because that a larger $\gamma$ means a more restrictive assumption on the denoiser, and therefore results in a worse performance. In practice, we recommend to set $t\approx t_0(\gamma)$, with $\gamma=0.25$. No significant denoising improvement is observed when $\gamma<0.25$.

5.) Although the paper presents a rich and rigorous theoretical framework, the presentation is overwhelmingly formal and math-heavy, with long sections dominated by formulas and lemmas. It may be very difficult for readers—especially those not from optimization theory—to grasp the main ideas. Some intuitive explanation, diagrams, or simplified takeaways would be highly beneficial.

Thank you for this valuable feedback. After thorough reading our own paper, we have to admit that, the paper is indeed a little math-heavy. We will revise the paper carefully, and add some intuitive explanation, diagrams, simplified takeaways, and proof sketches to make the paper easier to read.

Please let us know if you have more suggestions to polish the paper and raise the rating. We are very willing to revise the paper based on your comments to further improve its quality.

Dear reviewer 8mzP,

Thank you so much for your professional feedback and suggestions on the originality and theory. We sincerely appreciate your positive and negative comments, which are extremely helpful. We have carefully considered your comments and are committed to refining our paper. Below, you will find our detailed response as we strive for excellence.

Responses to weaknesses:

1.) Cocoercivity comes down to relaxing by a parameter $\gamma$ the firmly non-expansiveness soft constraint of Terris and al. [66].

Cocoercivity is a much broader (and weaker) concept than firm non-expansiveness, and residual non-expansive. In particular, if $\gamma=1$, $\gamma$-cocoercivity is equivalent to firm non-expansiveness; if $\gamma=0.5$, $\gamma$-cocoercivity is equivalent to residual non-expansiveness. Terris et al. proposed to softly enforce the firm non-expansiveness of the denoiser. By doing so, the denoiser is proved to be a resolvent of some monotone operator prior. If further assumed conservativity, this means that the learned prior is convex. However, as reported by many past papers [28,71], enforcing such constraints significantly compromises the denoising performance and PnP restoration. This paper relaxes this contraint to cocoercivity. If the denoiser is $\gamma$-cocoercive with $\gamma\in(0,1)$, the denoiser is actually a resolvent of some weakly monotone operator prior. If further assumed conservativity, this means that the learned prior is weakly convex. Overall, this paper adopted the soft constraints technique by Terris et al., and extended it to a weaker concept.

Please note that, compared to the pioneer work by Terris and Pesquet et al., this paper has the following significant improvements:

(1) The learned prior is weakly convex, enabling much better PnP restoration performance. In Table 2, compared to RMMO-DRS, the proposed methods has about 1dB increase in PSNR.

(2) We proved the convergence of PnP-ADMM with an implicit weakly convex prior, and a non-strongly convex, non-smooth data fidelity prior.

2.) Enforcing conservativeness of the denoiser (i.e. having a denoiser which is a gradient) can also be found in the literature (RED, Gradient Step denoisers). Besides, the proof technique for convergence strongly rely on the prox characterization of [22] which has already been used a lot in Hurault's works.

Gradient Step denoiser (GS-DRUNet) is indeed a great way to enforce the conservativeness. However, please note that, the training and the inference procedure of GS-DRUNet both need to calculate Jacobian-vector product, because there is always a Jacobian term in GS-DRUNet. This increases the training cost.

This paper wants to propose an alternative way to learn a conservative denoiser. We introduced the Helmholtz decomposition to study the denoising geometry, and showed that the Hamiltonian part is useless for denoising, and thus an ideal denoiser should be conservative. Besides, this decomposition also inspires us how to train a conservative denoiser without modifying the network architechture as in GS-DRUNet.

For the proof, we do utilize the great work [22] by Gribonval and Nikolova. Although Hurault has also use it to characterize the denoiser, there is an important difference between our paper and Hurault's paper. The Lipschitz assumption on the denoiser in Hurault's paper is residual non-expansiveness. This assumption is still restrictive, especiall when the noise level is large. In our paper, the assumption is relaxed to be cocoercive with even small $\gamma$. This relaxation $**fully**$ utilized the result in [22]. Mathematically speaking, since the `Lipschitz + be some gradient of a convex potential' can be equivalently rewritten as 'the denoiser is conservative and cocoercive', the proposed CoCo denoiser is a $**tight**$ utilization of the characterization in [22], and therefore provides better PnP performance.

Responses to questions:

1. We apologize for this careless mistake: in fact, in line 85-88 in the introduction, we have already talked about Bregman-PnP by Hurault et al.. However, we cited the wrong paper. We actually wanted to cite the paper `Convergent Bregman Plug-and-Play Image Restoration for Poisson Inverse Problems' published in Neurips. Besides, in the experiments, we have already compared the proposed method with Bregman-PnP (termed as B-PnP in the paper), see Table 2 and Figure 2.

Besides, thank you for your clarification! We totally agree with your statement that, we proposed to use a Gaussian denoiser as a prox rather than a Bregman prox. Compared to Bregman-PnP, the proposed CoCo-PnP has the following advantages:

(1) In the context of PnP, people use Gaussian denoisers mainly because that Gaussian denoisers are easier to get than other type denoisers. However, as stated in the introduction, there is a lacking literature on using Gaussian denoisers as prox to solve Poisson inverse problems. Our paper aims to solve this problem with Gaussian denoiser as a proximal operator of some implicit weakly convex prior.

(2) As reported in Pages 10 and 30 in the Appendix of the paper `Convergent Bregman Plug-and-Play Image Restoration for Poisson Inverse Problems', for large noise levels, Bregman denoiser does not perform well. Besides, when applied to deblurring with Poisson noises, Bregman-PnP does not outperform existing methods. We have also observed this in the experiments, see Table 2 and Figure 2 for example. Compared to Bregman-PnP (termed B-PnP in our paper), the proposed B-PnP performs significantly better. This suggests that using Gaussian denoisers in Poisson inverse problems is more suitable.

2. The Helmholtz decomposition is useful to not only interpret the denoising geometry, but also inspire us the design the loss functions, such that the learned denoiser is encouraged to be conservative. Proposition 2 does only require the denoiser to be Lipschitz and be some gradient of a convex potential. But, please note that, these two requirements are equivalent to that `the denoiser is conservative and cocoercive'. For the proof, please see line 596-600 in our paper, in which we have proved that the denoiser is a gradient of some convex potential.

Intuitively speaking, the `convexity' of the potential mainly comes from the cocoercivity of the denoiser. Please note that cocoercivity actually requires the denoiser to be a monotone operator. The gradient of a convex potential is a special case of monotone operators.

3. Unfortunately, we are not able to compute it. This is because that, we can not access the convex potential of the denoiser. Different from the great pioneering works by Hurault et al., the proposed denoiser has an implicit potential. According to Theorem 3 in [22], or Proposition 3.1 by Hurault et al. in [29], in order to compute the prior, we need to access the potential. Therefore, we chose to report the PSNR and relative error curves in Figures 2-3, 6-7. These figures validates the convergence.

4. For a differentiable data fidelity term, in the context of operator splitting algorithms, there are basically two ways to use it. We can use the fidelity in a backward way, that is, to use its proximal operator like in PnP-ADMM. Since the fidelity is convex, the proximal operator of such fidelity has good property: it is firmly non-expansive. We can also use the fidelity in a forward way, that is, to use its gradient like in PGD or FBS. But since the fidelity is not smooth, we have to smooth it by considering its Moreau envelope. We chose to include this algorithm not only for completedness of the two ways to use the fidelity, but also that because we wanted to show that, sometimes, the first motivation can not be solved by simply smoothing the non-smooth fidelity. In fact, as reported in Table 2, CoCo-ADMM performs better than CoCo-PEGD: this is because that the fidelity is changed to its envelope, and there can not fully reflect the Poisson distribution.

Thank you again for such professional comments on our paper. Please let us know if you have more suggestions. We are very willing to revise our paper based on your comments.