HD 131399 | 2015-06-12 | Ab |  0.845 $\pm$ 0.004 |  195.6 $\pm$ 0.2 |  ${}^{\text{(d)}}$(7.94 $\pm$ 2.10)$\times 10^{-6}$ |  (9.36 $\pm$ 8.01)$\times 10^{-6}$ |  (1.05 $\pm$ 0.17)$\times 10^{-5}$ |  8.23$\times 10^{-6}$
HD 131399 | 2016-05-07 | Ab |  0.825 $\pm$ 0.003 |  194.2 $\pm$ 0.2 |  ${}^{\text{(d)}}$(7.66 $\pm$ 2.21)$\times 10^{-6}$ |  (7.58 $\pm$ 2.10)$\times 10^{-6}$ |  (7.92 $\pm$ 1.37)$\times 10^{-6}$ |  7.86$\times 10^{-6}$
HIP 65426 | 2017-02-09 | b |  0.834 $\pm$ 0.001 |  149.9 $\pm$ 0.1 |  ${}^{\text{(e)}}$(9.91 $\pm$ 7.38)$\times 10^{-5}$ |  (6.20 $\pm$ 0.46)$\times 10^{-5}$ |  (5.87 $\pm$ 0.17)$\times 10^{-5}$ |  6.17$\times 10^{-5}$
HIP 65426 | 2018-05-13 | b |  0.825 $\pm$ 0.001 |  149.9 $\pm$ 0.1 |  ${}^{\text{(a)}}$(8.47 $\pm$ 2.24)$\times 10^{-5}$ |  (9.08 $\pm$ 2.44)$\times 10^{-5}$ |  (8.98 $\pm$ 0.09)$\times 10^{-5}$ |  9.06$\times 10^{-5}$
|  |  |  |  |  ${}^{\text{(f)}}$(8.39 $\pm$ 2.22)$\times 10^{-5}$ |  |  |
HIP 72192 | 2015-06-11 | CC1 |  4.632 $\pm$ 0.012 |  166.7 $\pm$ 0.1 | – |  (2.94 $\pm$ 2.31)$\times 10^{-6}$ |  (3.24 $\pm$ 1.07)$\times 10^{-6}$ |  2.88$\times 10^{-6}$
HIP 72192 | 2015-06-11 | CC2 |  4.891 $\pm$ 0.006 |  342.8 $\pm$ 0.1 | – |  (5.44 $\pm$ 1.85)$\times 10^{-6}$ |  (6.07 $\pm$ 1.10)$\times 10^{-6}$ |  5.95$\times 10^{-6}$

In this section, we compare the photometry estimations obtained with the
tested algorithms on some real known sources present in the considered
datasets. When available, we also compare them with measurements published in
the literature that have been obtained either with the SPHERE Data Center
implementations of TLOCI and PCA or with ANDROMEDA (for HD 95086 b only). For
PCA (VIP), the reported uncertainties are computed from the residual combined
image (i.e., after subtraction of the on-axis PSF, derotation, and stacking)
in an annulus at the same angular separation than the source of interest. For
PACO, the photometry is estimated by accounting for the spatial correlations
of the nuisance component. As concerns the proposed algorithm, and as for its
detection stage, the control of the uncertainties is an intricate task given
that the core of the deep model works as a black-box. For that reason, we did
not produce estimations of the standard-deviation, which is left for future
work. For the experiments we conduct with PCA (VIP), PACO and the proposed
method, we use the same pre-reduction of the raw observations and the same
off-axis PSF template that has been measured prior the science observations.
It is not ensured that these specificity also hold for photometry estimates
extracted from the literature, which could induce (unknown) systematic
variations that are not taken into account. Similarly, given the variability
of the observing conditions between different observations, absolute
comparisons between multiple epochs of the same target should be done with
caution.

Table 8 reports the photometry estimations on known real sources. Even in the
absence of absolute ground-truth, we can make some relative comparisons for a
given dataset and for a given source. In that view, photometry estimates
produced by the proposed method are compatible with published and/or obtained
results with PCA (VIP) and PACO for almost all sources. The largest
discrepancy is for CC5 of HIP 88399 (2018-04-11), where a factor six lies
between the estimates from PCA (VIP) and the proposed algorithm. Very likely,
the PCA (VIP) estimate is not reliable as the source is located at large
angular separation, i.e. in an area of the field of view where PCA (VIP) is
prone to large errors (see Fig. 17). Besides, PACO and the proposed method
lead to quite close estimations, which also supports the previous claim.

Photometry estimates of CC-12 that we identified from HD 95086 datasets in
Sect. 4.2.3 are consistent among the different algorithms. They are also
consistent between the two epochs where CC-12 has been detected. For CC-11
that we identified from the same datasets, we note a discrepancy by about 40%
between the estimates obtained (i) by the proposed method, and (ii) by the
other tested algorithms. However, this two groups of estimations are
consistent between the 2015 and 2018 epochs. The discrepancy could be
attributed to the presence of a bright background source (denoted by CC-8 in
Fig. 16) that could lead to an overestimation of the photometry with
algorithms that do not rely on a training step with synthetic sources. In any
case, these two candidate point-like sources should be considered with
caution.

##  5 Conclusion

We have described the key principles of a new algorithm for detecting and
characterizing point-like sources at high contrast from ADI observations. The
detection stage combines the statistics-based model of PACO with deep learning
in a three step procedure: (i) the data are centered and whitened using the
PACO framework, (ii) a CNN is trained to detect synthetic sources from the
pre-processed images, and (iii) a detection map is inferred. While the CNN
itself works as a black-box approach, the proposed method encompasses prior
domain knowledge such as the apparent motion of sources and the expected shape
of the exoplanetary signal inside the ADI datasets. More importantly, the
proposed detection approach capitalizes on the statistical model of the
nuisance component embedded in PACO to improve the stationarity and the
contrast during a pre-processing step. Once a candidate source has been
identified, its photometry can be estimated using a dedicated deep learning
module, also trained in a supervised fashion. Tested on eleven SPHERE-IRDIS
datasets, the proposed detection method performs on-par with or better than
PACO and other standard algorithms of the field like PCA in terms of
precision-recall trade-off.

We are currently working on the extension of the proposed algorithm for the
joint processing of multi-spectral datasets such as the ones provided by the
SPHERE-IFS instrument using the angular plus spectral differential imaging
technique. Besides, we are currently investigating the three main limitations
of the proposed approach: (i) the lack of control of the uncertainties, (ii)
the data-dependence of the learning procedure which does not take benefits
from multiple observations to build a more general and robust model of the
nuisance component, (iii) its task-dependence which is not adapted to
reconstruct spatially resolved objects like circumstellar disks. Besides, we
would like to incorporate within our deep models some meta-data (e.g.,
monitoring of the observing conditions, telemetry of the adaptive optics),
with the aim to further improve their sensitivity and robustness.

## Acknowledgements

This project is supported in part by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (COBREX;
grant agreement n° 885593). The work of TB and JP was supported in part by the
Inria/NYU collaboration, the Louis Vuitton/ENS chair on artificial
intelligence and the French government under management of Agence Nationale de
la Recherche as part of the Investissements d’avenir program, reference
ANR19-P3IA0001 (PRAIRIE 3IA Institute). The work of JM was supported in part
by the ERC grant number 714381 (SOLARIS project) and by ANR 3IA MIAI@Grenoble-
Alpes (ANR-19-P3IA-0003).

OF, TB, JM, JP, ML, and AML conceived and designed the method as well as the
analysis presented in this paper. OF and TB developed, tested, and implemented
the algorithm. OF and ML selected the raw data. ML pre-reduced them through
the SPHERE Data Centre. OF and TB performed the analysis of the data. OF, TB,
JM, JP, ML, and AML wrote the manuscript.

## Data Availability

The raw data used in this article are freely available on the ESO archive
facility at http://archive.eso.org/eso/eso_archive_main.html. They were pre-
reduced with the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble),
PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice), Observatoire de
Paris/LESIA (Paris), and Observatoire de Lyon/CRAL (Lyon, France). The
resulting pre-processed datasets will be shared based on reasonable request to
the corresponding author.

## References

  * Allard et al. (2003) Allard F., Guillot T., Ludwig H.-G., Hauschildt P. H., Schweitzer A., Alexander D. R., Ferguson J. W., 2003, in Symposium-International Astronomical Union. pp 325–332
  * Allard et al. (2007) Allard F., Allard N. F., Homeier D., Kielkopf J., McCaughrean M. J., Spiegelman F., 2007, Astronomy & Astrophysics, 474, L21
  * Amara & Quanz (2012) Amara A., Quanz S. P., 2012, Monthly Notices of the Royal Astronomical Society, 427, 948
  * Badrinarayanan et al. (2017) Badrinarayanan V., Kendall A., Cipolla R., 2017, IEEE transactions on pattern analysis and machine intelligence, 39, 2481
  * Beuzit et al. (2019) Beuzit J.-L., et al., 2019, Astronomy & Astrophysics, 631, A155
  * Boucaud et al. (2020) Boucaud A., et al., 2020, Monthly Notices of the Royal Astronomical Society, 491, 2481
  * Bowler (2016) Bowler B. P., 2016, Publications of the Astronomical Society of the Pacific, 128, 102001
  * Burrows et al. (1997) Burrows A., et al., 1997, The Astrophysical Journal, 491, 856
  * Cabayol et al. (2021) Cabayol L., et al., 2021, Monthly Notices of the Royal Astronomical Society, 506, 4048
  * Cantalloube et al. (2015) Cantalloube F., et al., 2015, Astronomy & Astrophysics, 582, A89
  * Cantalloube et al. (2020) Cantalloube F., et al., 2020, in Adaptive Optics Systems VII. pp 1027–1062
  * Carbillet et al. (2011) Carbillet M., et al., 2011, Experimental Astronomy, 30, 39
  * Castellá et al. (2016) Castellá B. F., et al., 2016, in Adaptive Optics Systems V. pp 697–710
  * Chabrier et al. (2000) Chabrier G., Baraffe I., Allard F., Hauschildt P., 2000, The Astrophysical Journal, 542, 464
  * Chauvin et al. (2004) Chauvin G., Lagrange A.-M., Dumas C., Zuckerman B., Mouillet D., Song I., Beuzit J.-L., Lowrance P., 2004, Astronomy & Astrophysics, 425, L29
  * Chauvin et al. (2005) Chauvin G., et al., 2005, Astronomy & Astrophysics, 438, L29
  * Chauvin et al. (2017) Chauvin G., et al., 2017, Astronomy & Astrophysics, 605, L9
  * Chauvin et al. (2018) Chauvin G., et al., 2018, Astronomy & Astrophysics, 617, A76
  * Cheetham et al. (2019) Cheetham A., et al., 2019, Astronomy & Astrophysics, 622, A80
  * Chen et al. (2010) Chen Y., Wiesel A., Eldar Y. C., Hero A. O., 2010, IEEE Transactions on Signal Processing, 58, 5016
  * Chomez et al. (2022) Chomez A., et al., 12/2022, submitted to Astronomy & Astrophysics
  * Conte et al. (1995) Conte E., Lops M., Ricci G., 1995, IEEE Transactions on Aerospace and Electronic Systems, 31, 617
  * Currie et al. (2012a) Currie T., Fukagawa M., Thalmann C., Matsumura S., Plavchan P., 2012a, The Astrophysical Journal Letters, 755, L34
  * Currie et al. (2012b) Currie T., et al., 2012b, The Astrophysical Journal Letters, 760, L32
  * Daglayan et al. (2022) Daglayan H., Vary S., Cantalloube F., Absil P.-A., Absil O., 2022, arXiv preprint arXiv:2210.10609
  * Dahlqvist et al. (2020) Dahlqvist C.-H., Cantalloube F., Absil O., 2020, Astronomy & Astrophysics, 633, A95
  * Dahlqvist et al. (2021a) Dahlqvist C.-H., Louppe G., Absil O., 2021a, Astronomy & Astrophysics, 646, A49
  * Dahlqvist et al. (2021b) Dahlqvist C.-H., Cantalloube F., Absil O., 2021b, Astronomy & Astrophysics, 656, A54
  * Delorme et al. (2017) Delorme P., et al., 2017, in Annual meeting of the French Society of Astronomy and Astrophysics.
  * Desgrange et al. (2022) Desgrange C., et al., 2022, arXiv preprint arXiv:2206.00425, 639, A113
  * Desidera et al. (2021) Desidera S., et al., 2021, Astronomy & Astrophysics, 651, A70
  * Dohlen et al. (2008) Dohlen K., Saisse M., Origne A., Moreaux G., Fabron C., Zamkotsian F., Lanzoni P., Lemarquis F., 2008, in SPIE Astronomical Telescopes + Instrumentation. p. 701859
  * Fergus et al. (2014) Fergus R., Hogg D. W., Oppenheimer R., Brenner D., Pueyo L., 2014, The Astrophysical Journal, 794, 161
  * Flasseur et al. (2018a) Flasseur O., Denis L., Thiébaut É., Langlois M., 2018a, in IEEE International Conference on Image Processing. pp 2735–2739
  * Flasseur et al. (2018b) Flasseur O., Denis L., Thiébaut É., Langlois M., 2018b, Astronomy & Astrophysics, 618, A138
  * Flasseur et al. (2018c) Flasseur O., Denis L., Thiébaut É. M., Langlois M., 2018c, in SPIE Astronomical Telescopes + Instrumentation. p. 107032R
  * Flasseur et al. (2020a) Flasseur O., Denis L., Thiébaut É., Langlois M., 2020a, Astronomy & Astrophysics, 634, A2
  * Flasseur et al. (2020b) Flasseur O., Denis L., Thiébaut É., Langlois M., 2020b, Astronomy & Astrophysics, 637, A9
  * Flasseur et al. (2021) Flasseur O., Thé S., Denis L., Thiébaut É., Langlois M., 2021, A&A, 651, A62
  * Flasseur et al. (2022) Flasseur O., Bodrito T., Mairal J., Ponce J., Langlois M., Lagrange A.-M., 2022, in Adaptive Optics Systems VIII. pp 1154–1167
  * Galicher et al. (2018) Galicher R., et al., 2018, Astronomy & Astrophysics, 615, A92
  * Gawlikowski et al. (2021) Gawlikowski J., et al., 2021, arXiv preprint arXiv:2107.03342
  * Gebhard et al. (2022) Gebhard T. D., Bonse M. J., Quanz S. P., Schölkopf B., 2022, arXiv preprint arXiv:2204.03439
  * Gonzalez et al. (2016) Gonzalez C. G., Absil O., Absil P.-A., Van Droogenbroeck M., Mawet D., Surdej J., 2016, Astronomy & Astrophysics, 589, A54
  * Gonzalez et al. (2017) Gonzalez C. A. G., et al., 2017, The Astronomical Journal, 154, 12pp
  * Gonzalez et al. (2018) Gonzalez C., Absil O., Van Droogenbroeck M., 2018, Astronomy & Astrophysics, 613, A71
  * Goodfellow et al. (2014) Goodfellow I., Pouget-Abadie J., Mirza M., Xu B., Warde-Farley D., Ozair S., Courville A., Bengio Y., 2014, Advances in neural information processing systems, 27
  * He et al. (2015) He K., Zhang X., Ren S., Sun J., 2015, in Proceedings of the IEEE international conference on computer vision. pp 1026–1034
  * He et al. (2016) He K., Zhang X., Ren S., Sun J., 2016, in Proceedings of the IEEE conference on computer vision and pattern recognition. pp 770–778
  * Huertas-Company & Lanusse (2022) Huertas-Company M., Lanusse F., 2022, arXiv preprint arXiv:2210.01813
  * Hüllermeier & Waegeman (2021) Hüllermeier E., Waegeman W., 2021, Machine Learning, 110, 457
  * Jensen-Clem et al. (2017) Jensen-Clem R., et al., 2017, The Astronomical Journal, 155, 19
  * Jovanovic et al. (2015) Jovanovic N., et al., 2015, Publications of the Astronomical Society of the Pacific, 127, 890
  * Kay (1993) Kay S. M., 1993, Fundamentals of statistical signal processing: estimation theory. Prentice-Hall, Inc.
  * Kendall et al. (1948) Kendall M. G., Stuart A., Ord J. K., 1948, The advanced theory of statistics. Vol. 1, JSTOR
  * Keppler et al. (2018) Keppler M., et al., 2018, Astronomy & Astrophysics, 617, A44
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014, arXiv preprint arXiv:1412.6980
  * Lafrenière et al. (2007) Lafrenière D., Marois C., Doyon R., Nadeau D., Artigau E., 2007, The Astrophysical Journal, 660, 770
  * Lagrange et al. (2009) Lagrange A.-M., et al., 2009, Astronomy & Astrophysics, 493, L21
  * Langlois et al. (2021) Langlois M., et al., 2021, Astronomy & Astrophysics, 651, A71
  * Ledoit & Wolf (2004) Ledoit O., Wolf M., 2004, Journal of Multivariate Analysis, 88, 365
  * Macintosh et al. (2014) Macintosh B., et al., 2014, Proceedings of the National Academy of Sciences, 111, 12661
  * Macintosh et al. (2015) Macintosh B., et al., 2015, Science, 350, 64
  * Marois et al. (2006) Marois C., Lafrenière D., Doyon R., Macintosh B., Nadeau D., 2006, The Astrophysical Journal, 641, 556
  * Marois et al. (2008) Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J., Lafrenière D., Doyon R., 2008, Science, 322, 1348
  * Marois et al. (2013) Marois C., Correia C., Véran J.-P., Currie T., 2013, Proceedings of the International Astronomical Union, 8, 48
  * Marois et al. (2014) Marois C., Correia C., Galicher R., Ingraham P., Macintosh B., Currie T., De Rosa R., 2014, in SPIE Astronomical Intrumentation + Telescopes. p. 91480U
  * Mawet et al. (2014) Mawet D., et al., 2014, The Astrophysical Journal, 792, 97
  * Mesa et al. (2019) Mesa D., et al., 2019, Monthly Notices of the Royal Astronomical Society, 488, 37
  * Milletari et al. (2016) Milletari F., Navab N., Ahmadi S.-A., 2016, in 2016 fourth international conference on 3D vision (3DV). pp 565–571
  * Morzinski et al. (2014) Morzinski K. M., et al., 2014, in Adaptive Optics Systems IV. p. 914804
  * Mugnier et al. (2009) Mugnier L. M., Cornia A., Sauvage J.-F., Rousset G., Fusco T., Védrenne N., 2009, Journal of the Optical Society of America A, 26, 1326
  * Nielsen et al. (2012) Nielsen E. L., et al., 2012, The Astrophysical Journal, 750, 53
  * Nielsen et al. (2017) Nielsen E. L., et al., 2017, The Astronomical Journal, 154, 218
  * Nielsen et al. (2019) Nielsen E. L., et al., 2019, The Astronomical Journal, 158, 13
  * Nocedal & Wright (1999) Nocedal J., Wright S. J., 1999, Numerical optimization. Springer
  * Pairet et al. (2019) Pairet B., Cantalloube F., Gomez Gonzalez C. A., Absil O., Jacques L., 2019, Monthly Notices of the Royal Astronomical Society, 487, 2262
  * Paszke et al. (2019) Paszke A., et al., 2019, Curran Associates, Inc., pp 8024–8035
  * Pavlov et al. (2008) Pavlov A., Möller-Nilsson O., Feldt M., Henning T., Beuzit J.-L., Mouillet D., 2008, in SPIE Astronomical Telescopes + Instrumentation. p. 701939
  * Pueyo (2018) Pueyo L., 2018, Handbook of Exoplanets, pp 705–765
  * Rameau et al. (2013a) Rameau J., et al., 2013a, The Astrophysical journal letters, 772, L15
  * Rameau et al. (2013b) Rameau J., et al., 2013b, The Astrophysical Journal Letters, 779, L26
  * Reddi et al. (2019) Reddi S. J., Kale S., Kumar S., 2019, arXiv preprint arXiv:1904.09237
  * Ronneberger et al. (2015) Ronneberger O., Fischer P., Brox T., 2015, in International Conference on Medical image computing and computer-assisted intervention. pp 234–241
  * Ruffio et al. (2017) Ruffio J.-B., et al., 2017, The Astrophysical Journal, 842, 14
  * Samland et al. (2021) Samland M., Bouwman J., Hogg D., Brandner W., Henning T., Janson M., 2021, Astronomy & Astrophysics, 646, A24
  * Santos (2008) Santos N. C., 2008, New Astronomy Reviews, 52, 154
  * Schneider et al. (2011) Schneider J., Dedieu C., Le Sidaner P., Savalle R., Zolotukhin I., 2011, Astronomy & Astrophysics, 532, A79
  * Simonyan & Zisserman (2015) Simonyan K., Zisserman A., 2015
  * Soummer et al. (2012) Soummer R., Pueyo L., Larkin J., 2012, The Astrophysical Journal Letters, 755, L28
  * Sudre et al. (2017) Sudre C. H., Li W., Vercauteren T., Ourselin S., Jorge Cardoso M., 2017, in , Deep learning in medical image analysis and multimodal learning for clinical decision support. Springer, pp 240–248
  * Traub & Oppenheimer (2010) Traub W. A., Oppenheimer B. R., 2010, Exoplanets, pp 111–156
  * Vigan et al. (2010) Vigan A., Moutou C., Langlois M., Allard F., Boccaletti A., Carbillet M., Mouillet D., Smith I., 2010, Monthly Notices of the Royal Astronomical Society, 407, 71
  * Vigan et al. (2021) Vigan A., et al., 2021, Astronomy & Astrophysics, 651, A72
  * Wagner et al. (2016) Wagner K., Apai D., Kasper M., Kratter K., McClure M., Robberto M., Beuzit J.-L., 2016, Science, 353, 673
  * Wahhaj et al. (2015) Wahhaj Z., et al., 2015, Astronomy & Astrophysics, 581, A24
  * Wahhaj et al. (2021) Wahhaj Z., et al., 2021, arXiv preprint arXiv:2101.08268
  * Wainwright & Simoncelli (2000) Wainwright M. J., Simoncelli E. P., 2000, in Advances in Neural Information Processing Systems. pp 855–861
  * Wang et al. (2020) Wang L., Wang C., Sun Z., Chen S., 2020, IEEE Access, 8, 167939
  * Wertz et al. (2017) Wertz O., Absil O., González C. G., Milli J., Girard J. H., Mawet D., Pueyo L., 2017, Astronomy & Astrophysics, 598, A83
  * Yalniz et al. (2019) Yalniz I. Z., Jégou H., Chen K., Paluri M., Mahajan D., 2019, arXiv preprint arXiv:1905.00546
  * Yip et al. (2019) Yip K. H., et al., 2019, in Joint European Conference on Machine Learning and Knowledge Discovery in Databases. pp 322–338

##  Appendix A Discussion about the statistical model and the whitening of the
observations

###  A.1 Refinement of the statistical model

Table 9: Comparison between a multi-variate Gaussian model (Eqs. (2)) and a GSM model (Eqs. (16)) for the statistical modeling of the nuisance component in the pre-processing step of the proposed detection approach. Reported scores are AUC of F1R score (best when close to 1) as an overall measurement of the precision-recall trade-off of the underlying detector. Mean results and standard-deviation are obtained on the eleven SPHERE-IRDIS datasets considered in this work, see Sect. 4 for the recording logs. Ang. Sep. (”) | Gaussian model (default) | GSM model (variant)
---|---|---
$\bm{\left[0;2\right]}$” | $0.88\pm 0.03$ | $0.90\pm 0.03$
$\bm{\left[2;4\right]}$” | $0.90\pm 0.05$ | $0.91\pm 0.04$
$\bm{\left[4;6\right]}$” | $0.89\pm 0.04$ | $0.91\pm 0.01$
$\bm{\left[6;7\right]}$” | $0.88\pm 0.06$ | $0.90\pm 0.03$

Concerning the statistical model of the nuisance component, the multi-variate
Gaussian assumption described in Sect. 2.1.1 is a convenient approximation
leading to closed-form expressions for the underlying estimators. However,
this formulation neglects all temporal fluctuations of the data. In Flasseur
et al. (2020a, b), we consider a refinement of this model using a multi-
variate Gaussian scaled mixture (GSM; Conte et al. (1995); Wainwright &
Simoncelli (2000)). It amounts to model the distribution of patch
$\bm{f}_{n,t}$ centered around pixel $n$ at time $t$ by a temporally weighted
multi-variate Gaussian
$\mathcal{N}(\bm{m}_{n},\sigma_{n,t}^{2}\,\mathbf{C}_{n})$. Under this model,
the sample estimates
$\\{\widehat{\bm{m}}_{n};\widehat{\sigma}_{n,t}^{2};\widehat{\mathbf{S}}_{n}\\}$
of the local mean, of the temporal scaling factors, and of the covariance
coming from the maximum likelihood are obtained with a fixed-point iterative
method as follows (Flasseur et al., 2020a):

$\begin{cases}\widehat{\bm{m}}_{n}=\frac{1}{\sum\limits_{t=1}^{T}1/\hat{\sigma}_{n,t}^{2}}\cdot\sum\limits_{t=1}^{T}\frac{1}{\hat{\sigma}_{n,t}^{2}}\mathbf{E}_{n,t}\,\bm{r}\,\in\mathbb{R}^{K}\,,\vspace{1mm}\\\
\widehat{\sigma}_{n,t}^{2}=\frac{1}{K}(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})^{\top}\,\widehat{\mathbf{S}}_{n}^{-1}\,(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})\,\in\mathbb{R}_{+}\,,\vspace{1mm}\\\
\widehat{\mathbf{S}}_{n}=\frac{1}{T}\sum\limits_{t=1}^{T}\frac{1}{\hat{\sigma}_{n,t}^{2}}(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})(\mathbf{E}_{n,t}\,\bm{r}-\widehat{\bm{m}}_{n})^{\top}\,\in\mathbb{R}^{K\times
K}\,.\end{cases}$ (16)

The regularized covariance $\widehat{\mathbf{C}}_{n}$ is estimated on the fly
by shrinkage of $\widehat{\mathbf{S}}_{n}$ as defined in Eq. (5), with
$Q=\left(\sum_{t=1}^{T}1/\widehat{\sigma}_{n,t}^{2}\right)^{2}\big{/}\left(\sum_{t=1}^{T}1/\widehat{\sigma}_{n,t}^{4}\right)$
the equivalent number of patches involved in the computation of
$\widehat{\mathbf{S}}_{n}$ in the presence of the weighting factors
$\\{\widehat{\sigma}_{n,t}^{2}\\}_{t=1:T}$, see Flasseur et al. (2020a). The
scaling factors $\\{\widehat{\sigma}_{n,t}^{2}\\}_{t=1:T}$ can be interpreted
as the local variance of the the residual data (i.e., after centering and
whitening). This approach allows to identify and to neutralize outliers,
taking the form both of spatially interpolated defective pixels and of local
areas displaying fluctuations on some temporal frames larger than on other
ones. These properties transfer to the estimators of the statistics of the
nuisance component with an improved robustness against bad data, thus leading
to a better detection sensitivity and characterization accuracy.

We tested to integrate a GSM model in the pre-processing step of the proposed
algorithm. Table 9 compares, with the same procedure as described in Sect.
2.2.3, the detection performance in terms of precision and recall of the
proposed algorithm using either a multi-variate Gaussian model (Eqs. (2)) or a
GSM model (Eqs. (16)) in the pre-processing step. It illustrates that, as for
PACO, the GSM model leads to an improved detection sensitivity, with an
improved stability over datasets (i.e., smaller error bars). However, the gain
is significantly smaller than the typical gain, reaching 10 to 15%, obtained
when substituting the multi-variate Gaussian assumption with a GSM model in
the PACO algorithm (Flasseur et al., 2020a). This observation can be
attributed to the presence of the additional learning stage of our proposed
method that, hopefully, partially correcting for the approximate fidelity
(with respect to the observations) of the statistical model embedded in the
pre-processing step.

Given the increased computational burden of the proposed detection approach
with a GSM model by a factor between 10 and 30 (which represents the typical
number of iterations needed to reach convergence of the estimators (16)), we
recommend, for practical reasons, to use the standard version of deep PACO
embedding a multi-variate Gaussian model, as defined in Sect. 2.1.1. This is
also the choice we made for the presentation of the results in the main core
of this paper. The alternative version of the algorithm embedding a GSM model
can be reserved to refine, in a second step, the reduction of some datasets
with ambiguous detections. When comparing with state-of-the-art detection
algorithms, we use for PACO the version embedding a GSM model to compare
fairly the proposed method against the best setting of existing methods.

###  A.2 Refinement of the whitening of the observations

Table 10: Comparison between non-normalized (Eqs. (7)) and normalized (Eqs. (17)) outputs produced by the pre-processing step of the proposed detection approach. Reported scores are AUC of F1R score (best when close to 1) as an overall measurement of the precision-recall trade-off of the underlying detector. Mean results and standard-deviation are obtained on the eleven SPHERE-IRDIS datasets considered in this work, see Sect. 4 for the recording logs. Ang. Sep. (”) | Not normalized (default) | Normalized (variant)
---|---|---
$\bm{\left[0;2\right]}$ | $0.88\pm 0.03$ | $0.88\pm 0.04$
$\bm{\left[2;4\right]}$ | $0.90\pm 0.05$ | $0.92\pm 0.05$
$\bm{\left[4;6\right]}$ | $0.89\pm 0.04$ | $0.92\pm 0.03$
$\bm{\left[6;7\right]}$ | $0.88\pm 0.06$ | $0.91\pm 0.05$
Figure 20: (a): Illustration of two data frames impacted by a wind-driven halo
near the star. The direction of main elongation of the halo is symbolized by
the white arrow. (b): Detection map obtained on the corresponding dataset with
the proposed detection algorithm embedding the whitening procedure either
defined in Sect. 2.1.2 (i.e., without normalization by the whitened off-axis
PSF) or defined in Appendix A.2 (i.e., with normalization by the whitened off-
axis PSF). The four known exoplanets (HR 8799 b, c, d, e) are circled in
green. Dataset: HR 8799 (2015-07-04), see Sect. 4 for the observation logs.

The pre-processing step (centering and local whitening) described in Sect.
2.1.2 is computationally quite efficient since it involves, for non-
overlapping square patches, only $\lfloor N/K\rceil$ ($\simeq 10^{4}$ for the
SPHERE-IRDIS instrument) matrix multiplications of size $K\times K$. However,
it has two drawbacks: (i) it leads to some spatial discontinuities between
adjacent patches, (ii) it does not account for the spatially non-stationary
transformation induced by the whitening process on the off-axis PSF in terms
of shape and of intensity. In practice, we observe that this approach
sometimes lacks of robustness for observations recorded under medium to bad
conditions, i.e., where the unmixing between the nuisance component and the
sought objects is even more difficult. Limitation (i) can be addressed by
considering overlapping patches (e.g., with a patch stride of one pixel).
Limitation (ii) can be addressed by adding an output term explicitly
accounting for the shape and intensity transformation induced by the whitening
process. In this context, the pre-processed images $\widetilde{\bm{r}}$ in
$\mathbb{R}^{N\times 2}$ are now formed by the concatenation of two images
$\widetilde{\bm{a}}$ in $\mathbb{R}^{N}$ and $\widetilde{\bm{b}}$ in
$\mathbb{R}^{N}$ defined by:

$\begin{cases}\widetilde{b}_{n^{\prime}}=\frac{\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\bm{h}^{\top}\,\widehat{\mathbf{C}}^{-1}_{n}\,(\bm{r}_{n}-\widehat{\bm{m}}_{n})\right]_{n^{\prime}}}{q_{n^{\prime}}}\,,\forall
n^{\prime}\in\llbracket 1;N\rrbracket\,,\\\
\widetilde{a}_{n^{\prime}}=\frac{\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\bm{h}^{\top}\,\widehat{\mathbf{C}}^{-1}_{n}\,\bm{h}\right]_{n^{\prime}}}{q_{n^{\prime}}}\,,\forall
n^{\prime}\in\llbracket 1;N\rrbracket\,,\end{cases}$ (17)

with $q_{n^{\prime}}$ the number of patches averaged at each location
$n^{\prime}$ of the field of view, as defined in Eq. (6). For overlapping
squared patches with a unit patch stride, $q_{n^{\prime}}$ is equal to $K$
almost everywhere, excepted on the borders of the field of view where it
progressively tends to zero. The term $\bm{b}$ can be interpreted as the
correlation between the whitened off-axis PSF and the centered plus whitened
observations, while $\bm{a}$ is a normalization term representing the auto-
correlation of the whitened off-axis PSF. It can be noted that for each pixel
$n^{\prime}$ of the field of view, the ratio
$\widetilde{b}_{n^{\prime}}/\sqrt{\widetilde{a}_{n^{\prime}}}$ (respectively,
the ratio $\widetilde{b}_{n^{\prime}}/\widetilde{a}_{n^{\prime}}$) corresponds
to the S/N of detection (respectively, to the source flux) that would be
estimated by PACO at pixel $n^{\prime}$ (Flasseur et al., 2018b). To prevent
the deep model (built from the outputs of the pre-processing step) to learn
only the mapping
$\widetilde{\bm{r}}\rightarrow\widetilde{\bm{b}}/\sqrt{\widetilde{\bm{a}}}$,
we resort to a residual learning procedure. It consists in evaluating the loss
function (see Sect. 2.2.3) jointly on the detection map produced by the
proposed algorithm and also on the PACO S/N map
$\widetilde{\bm{b}}/\sqrt{\widetilde{\bm{a}}}$ computed on the fly. This
strategy explicit rewards the deep model to perform better than PACO.

Table 10 compares, with the same procedure as described in Sect. 2.2.3, the
detection performance in terms of precision and recall of the proposed
algorithm using either the whitening procedure described in Sect. 2.1.2, and
the whitening procedure described in this Appendix. It shows that for the
eleven SPHERE-IRDIS datasets we study in details in this work, the variant
approach accounting for the transformation induced by the whitening process on
the off-axis PSF leads only to a slight improvement of the overall detection
performance. This is due to the fact that none of the eleven considered
datasets was recorded under bad observing conditions. Figure 20 gives a
qualitative comparison between the two approaches for a dataset of HR 8799
(see Sect. 4 for the recording logs) impacted by the wind-driven halo effect.
The variant procedure described in this appendix allows to avoid numerous
evident false alarms occurring when the normalization of the pre-processed
frames is omitted.

Compared to the whitening procedure described in Sect. 2.1.2, the variant
described in this appendix has a computational burden increased by a factor
$4\times K$, i.e. by typically a factor between 200 and 500 for the VLT/SPHERE
instrument. For practical reasons, we recommend to use the standard whitening
procedure by default, as defined in Sect. 2.1.2. This is also the choice we
made for the presentation of the results in the main core of this paper. The
alternative version of the algorithm accounting for the whitening of the off-
axis PSF can be reserved to refine, in a second step, the reduction of some
datasets for which obvious and numerous false alarms are experienced at
inference time.

##  Appendix B Main settings of the detection and characterization stages of
the proposed algorithm

Table 11 gives a quick-look summary of the main settings for the proposed
algorithm, both for the detection and for the characterization modules. Fields
are classified in three categories: pre-processing, generation of the training
set, and deep learning model.

Table 11: Summary of the main settings used for the detection and characterization modules of the proposed algorithm. ${}^{\text{(a)}}$We recall that the training sets are generated on the fly for the detection stage, i.e., the notion of epochs is used only to schedule the learning rate. For the characterization stage, the term epoch is used in a classical meaning so that the generated patches are seen multiple times (corresponding to the number of epochs) by the network. parameters | detection | characterization
---|---|---
$\blacktriangleright$ pre-processing
input | observations $\bm{r}\in\mathbb{R}^{N\times T}$ | observations $\bm{r}\in\mathbb{R}^{N\times T}$
strategy for known sources (training only) | temporal shuffling | source masking + temporal shuffling
temporal centering | yes | yes
spatial whitening | yes | no
whitening patch shape | square | –
whitening patch area $K$ |  $\in\llbracket 7^{2};12^{2}\rrbracket$ pixels (automatic) | –
whitening patch stride |  $\lfloor\sqrt{K}\,\rceil$ pixels (default) or 1 pixel (variants of App. A) | –
whitening statistical model | multi-variate Gaussian (default) or GSM (variant of App. A.1) | –
whitening output quantity |  $\widehat{\mathbf{L}}_{n}^{\top}\,(\bm{r}_{n}-\widehat{\bm{m}}_{n})$ s.t. $\widehat{\mathbf{L}}_{n}\,\widehat{\mathbf{L}}_{n}^{\top}=\widehat{\mathbf{C}}_{n}^{-1}\,,$ $\forall n\in\mathbb{P}$ (default) or | –
|  $\text{concat.}\left(\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\bm{h}^{\top}\,\widehat{\mathbf{C}}^{-1}_{n}\,(\bm{r}_{n}-\widehat{\bm{m}}_{n})\right]_{n^{\prime}};\left[\sum\limits_{n\in\mathbb{P}}\mathbf{E}_{n}^{\top}\,\bm{h}^{\top}\,\widehat{\mathbf{C}}^{-1}_{n}\,\bm{h}\right]_{n^{\prime}}\right)\,,\forall n^{\prime}\in\llbracket 1;N\rrbracket$ (variant of Appendix A.2) |
parallactic derotation | yes | yes
temporal collapsing | no | yes (after injections if any)
implantation | GPUs or CPUs (parallelized) | GPUs or CPUs (parallelized)
output | pre-processed observations $\widetilde{\bm{r}}\in\mathbb{R}^{N\times T}$ | pre-processed observations $\widetilde{\bm{r}}\in\mathbb{R}^{N\times T}$
$\blacktriangleright$ generation of the training set
input | pre-processed observations $\widetilde{\bm{r}}\in\mathbb{R}^{N\times T}$ | pre-processed observations $\widetilde{\bm{r}}\in\mathbb{R}^{N\times T}$
pre-processing update on injection arcs | yes | yes
data augmentation | yes | yes
total number $P$ of sources | $\in\llbracket 25,000;50,000\rrbracket$ | $\simeq 40,000$
number $S$ of training sets | $\in\llbracket 500;1,000\rrbracket$ | $\simeq 8,000$
number $P^{\left[s\right]}$ of sources per set $s$ | $\in\llbracket 1,10\rrbracket$ | $\in\llbracket 1;10\rrbracket$
location $\phi$ of sources | uniform in polar system | uniform in Cartesian system
contrast $\alpha$ of sources | uniform in $\left[3\widehat{\sigma}_{\phi_{p}}^{\texttt{PACO}};12\widehat{\sigma}_{\phi_{p}}^{\texttt{PACO}}\right]$ | uniform in $\left[1\times 10^{-6};3\times 10^{-5}\right]$ (default)
cropping patch area $J$ | – |  $31^{2}$ pixels
implantation | CPUs (parallelized) | CPUs (parallelized)
output | sets of pre-processed observations with injections $\\{\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{\left[s\right]}\in\mathbb{R}^{N\times T}\\}_{s=1:S}$ | sets of pre-processed patches with injections $\\{\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}^{[p]}\in\mathbb{R}^{J}\\}_{p=1:P}$
$\blacktriangleright$ deep learning model
input | a set of pre-processed observations with injections $\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{3.8315pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.68205pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{r}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{1.91574pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{r}\hss$\crcr}}}\limits}^{\left[s\right]}\in\mathbb{R}^{N\times T}$ | a pre-processed patch with injection $\mathchoice{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{4.02501pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.8175pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\bm{p}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern 0.80002pt$\hss\leavevmode\resizebox{2.0125pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\bm{p}\hss$\crcr}}}\limits}^{[p]}\in\mathbb{R}^{J}$
architecture | U-Net (Res-Net 18 backbone) | custom VGG-like
task | pixel-wise classification | regression
number of weights |  $\simeq 11$ millions |  $\simeq 1.2$ millions
pre-trained weights | no | no
optimization loss | Dice2 (overlap measure) | absolute relative error
optimizer | AMSGrad | Adam
validation metric | F1R score (precision / recall trade-off) | absolute relative error
batch size | 1 (set of pre-processed observations with injections) | 1024 (pre-processed patches with injections)
weight decay | $10^{-5}$ | 0
initial learning rate | $10^{-3}$ | $10^{-3}$
learning rate scheduling | yes (-10% every 10 epochs) | yes (-70% every 50 epochs)
number of epochs${}^{\text{(a)}}$ |  $\in\llbracket 1;100\rrbracket$ (on the fly) | 300 (fixed)
implantation | GPUs | GPUs
output | detection map $\widehat{\bm{y}}\in\left[0;1\right]^{M}$ | photometry estimates $\widehat{\alpha}\in\mathbb{R}_{+}$

See pages 1- of deep_paco_mnras_2022_supplementary_material_compressed.pdf


# Regularization Theory of the Analytic Deep Prior Approach

Clemens Arndt Center for Industrial Mathematics, University of Bremen, Germany
<EMAIL_ADDRESS>
(March 2022)

###### Abstract

The analytic deep prior (ADP) approach was recently introduced for the
theoretical analysis of deep image prior (DIP) methods with special network
architectures. In this paper, we prove that ADP is in fact equivalent to
classical variational Ivanov methods for solving ill-posed inverse problems.
Besides, we propose a new variant which incorporates the strategy of early
stopping into the ADP model. For both variants, we show how classical
regularization properties (existence, stability, convergence) can be obtained
under common assumptions.

##  1 Introduction

In particular the field of image processing (e.g. denoising, deblurring) is a
constant source for challenging inverse problems. The restoration of a
corrupted image is typically ill-posed, so regularization techniques are
needed to obtain a natural looking result. In other words, the restoration
method should incorporate some prior knowledge about the appearance of natural
images. However, dependent on the application it can be very difficult to give
a mathematically exact definition of what natural looking images are. This
makes it hard to encode such prior knowledge in a penalty term for classical
variational regularization approaches (e.g. TV regularization [30, 8]).

However, deep learning methods with convolutional neural networks have proven
to be quite successful in generating and restoring images [22, 16, 35]. One
reason for that is the use of appropriate training data, but [25] shows that
just the architecture of an untrained network can already serve as an image
prior. The so-called deep image prior (DIP) approach consists in optimizing
the weights of a neural network $\varphi_{\theta}$ to minimize the loss
function

$\frac{1}{2}\|A\varphi_{\theta}(z)-y^{\delta}\|^{2}$ (1.1)

for some forward operator $A$ and noisy data $y^{\delta}$ (the network’s input
$z$ is randomly chosen and kept fixed). Although no training data and no
penalty functional is used, DIP produces remarkable results in different image
processing tasks, as can be seen in [25]. Even challenging problems like
sparse angle computed tomography [3] or compressive sensing [19] can be solved
this way.

Developing regularization theory for deep learning methods is of high interest
[2]. A very prosporous approach is to combine classical theory with deep
learning (e.g. [26, 29]). The number of papers which analyze DIP from a
theoretical point of view is also growing. In [18] a functional is
constructed, which measures the ability of the neural network
$\varphi_{\theta}$ to approximate an arbitrary image. This functional can then
be used as a penalty term in a classical variational method. The authors of
[32] analyze how fast a DIP network approximates the low-frequency and high-
frequency components of the target image. By controlling this so-called
spectral bias, overfitting is avoided. In [20] the ability to denoise images
is attributed to convolutional layers, which are faster in fitting smooth
images than noisy ones. The role and the choice of hyperparameters for DIP
approaches is described in [34]. A Bayesian perspective is presented in [9],
where DIP is interpreted as a Gaussian process.

The choice of architecture is crucial for applications of DIP. Generative
neural networks are a natural choice due to their ability to reproduce natural
looking images. But the authors of [13] took a LISTA-like network [17] instead
to develop the so-called analytic deep prior (ADP) approach. This may not lead
to a better practical performance of DIP, but it’s the foundation for an
interesting theory. The main aspect consists in interpreting the training of a
neural network as the optimization of a Tikhonov functional. There is an
analogy to [1], where the penalty term for a Tikhonov functional is optimized.
But in contrast to that, the focus of [13] is on the forward operator inside
the functional (see section 2).

This work summarizes deeper investigations of the ADP model. The main result
(Theorem 3.3) is an equivalence between the ADP approach and classical Ivanov
methods [36]. Out of this follows a complete analysis of the regularization
properties of ADP including the existence of solutions, stability of
reconstructions and convergence towards the ground truth for vanishing noise.

In practical applications of DIP, gradient descent and early stopping is used
to minimize the loss function (1.1). Thus, a global (or at least a local)
minimum is not reached in general. While this fact was not considered in the
theoretical derivation of ADP, we propose a new variant (called ADP-$\beta$)
which incorporates the effect of early stopping into the model (section 3.2).
We also analyze the regularization properties of this new approach.

In section 4 we compare different numerical ways to compute ADP and DIP (with
a LISTA-like architecture) solutions of simple inverse problems.111Code
available at https://gitlab.informatik.uni-
bremen.de/carndt/analytic_deep_prior We find that numerical solutions of both
methods are mostly similar to each other, which is important for using the ADP
theory for interpretations of DIP. But there can also be observed some
interesting disparities between the different numerical ways. This illustrates
a crucial difference between the analytical definition of DIP as a
minimization problem and the numerical implementation as a gradient descent
iteration.

##  2 Preliminaries and methods

We consider an inverse problem based on the operator equation

$Ax^{\dagger}=y^{\dagger}$ (2.1)

where we want to recover the unknown ground truth $x^{\dagger}$ as good as
possible. The data $y^{\dagger}$ is typically not known exactly, but we have
only access to noisy data $y^{\delta}$.

######  Assumption 2.1.

We make the following assumptions for the inverse problem (2.1).

  * •

Let $X,Y$ be Hilbert spaces and $A\in L(X,Y)$.

  * •

There exists $x^{\dagger}\in X$ and for a given $\delta>0$, it holds
$\|y^{\delta}-y^{\dagger}\|\leqslant\delta$ for $y^{\delta}\in Y$.

  * •

Let $R\colon X\to[0,\infty]$ be a convex, coercive and weakly lower
semicontinuous functional with $R\not\equiv\infty$.

We recall the definition of Bregman distances, which we will use in Theorem
3.12 for a convergence result, similar to the ones in [6, 21].

######  Definition 2.2 (Bregman distance).

For a convex functional $R\colon X\to[0,\infty]$ with subdifferential
$\partial R$ and $\tilde{x},x\in X$, the Bregman distance is defined as the
set

$D_{R}(\tilde{x},x)=\\{R(\tilde{x})-R(x)-\langle
p,\tilde{x}-x\rangle\,|\,p\in\partial R(x)\\}.$ (2.2)

The DIP approach (introduced by [25]) for the inverse problem (2.1) consists
in solving

$\min_{\theta}\frac{1}{2}\|A\varphi_{\theta}(z)-y^{\delta}\|^{2}$ (2.3)

via a gradient descent w.r.t. the parameters $\theta$ of a neural network
$\varphi_{\theta}$, as already described in the introduction. Despite the use
of a neural network, DIP is a model-based approach and not data-based. To
derive the ADP approach, we have to make two assumptions (see [13] for
details).

The first one is choosing $\varphi_{\theta}$ to be a LISTA-like network [17],
which consists of several layers of the form

$x^{l+1}=S_{\alpha\lambda}(x^{l}-\lambda B^{*}(Bx^{l}-y^{\delta})),$ (2.4)

where $B=\theta$ is the trainable parameter. Originally, this architecture is
inspired by ISTA [12], an algorithm for finding sparse solutions of the
inverse problem (2.1), and $S_{\alpha\lambda}$ is the shrinkage function. More
general, we can choose $S_{\alpha\lambda}$ to be the proximal mapping of a
penalty functional $R$. Then, (2.4) equals a proximal forward-backward
splitting algorithm [11, Theorem 3.4] which converges to the solution of the
minimization problem

$\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x).$ (2.5)

The second assumption is letting the number of layers tend to infinity. This
might be difficult in practice (see section 4), but it causes the output
$\varphi_{\theta}(z)$ of the network to be a solution of (2.5). Therefore the
ADP model (introduced by [13]) is defined as

$\displaystyle\begin{split}&\min_{B\in
L(X,Y)}\frac{1}{2}\|Ax(B)-y^{\delta}\|^{2}\\\ \text{s.t.}\quad
x(B)=\,&\mathrm{arg}\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha
R(x).\end{split}$ (2.6)

While DIP is about optimizing the weights of a neural network, ADP is about
optimizing the forward operator in a Tikhonov functional. If we add an
additional regularization term for the operator $B$, we get the (new)
ADP-$\beta$ model

$\begin{split}\min_{B\in
L(X,Y)}\frac{1}{2}\|Ax(B)-y^{\delta}\|^{2}+\beta\|B-A\|^{2}\\\
\text{s.t.}\quad x(B)=\mathrm{arg}\min_{x\in
X}\,\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x).\end{split}$ (2.7)

The reason for this modification will be explained in section 3.2.

To guarantee uniqueness of $x(B)$, the functional $R$ should be strictly
convex, but this is not always required. If we assume $R$ even to be strongly
convex, $x(B)$ depends continuously on $B$ as the following theorem states. It
will be useful for proving existence and stability results for ADP-$\beta$.

######  Theorem 2.3.

Let $R\colon X\to[0,\infty]$ be a strongly convex, coercive and weakly lower
semicontinuous functional. Then

$x(B)=\mathrm{arg}\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x)$
(2.8)

depends continuously on $B\in L(X,Y)$.

The proof can be found in the appendix A.1.

##  3 Theoretical results

###  3.1 Equivalence to classical methods

DIP solutions of inverse problems are naturally restricted to be the output of
a neural network. Analogously, only elements of the set

$U_{\alpha R}=\left\\{\hat{x}\in X\,\bigg{|}\,\exists B\in
L(X,Y):\hat{x}=\mathrm{arg}\min_{x\in
X}\,\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x)\right\\}$ (3.1)

can be solutions of the ADP approach. By definition

$\displaystyle\min_{x\in U_{\alpha R}}\frac{1}{2}\|Ax-y^{\delta}\|^{2}$ (3.2)

is equivalent to the original ADP problem (2.6). To get a better understanding
of this minimization problem we investigate $U_{\alpha R}$. It will turn out
that the set $U_{\alpha R}$ can be characterized in a much easier way, even
without using an operator $B\in L(X,Y)$. For this purpose, we formulate the
following lemmas.

######  Lemma 3.1.

Let $R\colon X\to[0,\infty]$ be a convex, coercive and weakly lower
semicontinuous functional and $\hat{x}\in X$, $y^{\delta}\in Y$,
$y^{\delta}\neq 0$ and $\alpha>0$ be arbitrary. If there exists $v\in\partial
R(\hat{x})$ such that

$\alpha\langle v,\hat{x}\rangle\leqslant\frac{\|y^{\delta}\|^{2}}{4}$ (3.3)

holds, then there exists a linear operator $B\in L(X,Y)$ which fulfills

$\hat{x}=\mathrm{arg}\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha
R(x).$ (3.4)

The proof can be found in the appendix A.2.

######  Lemma 3.2.

Let $R\colon X\to[0,\infty]$ be a convex, coercive and weakly lower
semicontinuous functional and $\hat{x}\in X$, $y^{\delta}\in Y$, $\alpha>0$ be
arbitrary. If for every $v\in\partial R(\hat{x})$

$\alpha\langle v,\hat{x}\rangle>\frac{\|y^{\delta}\|^{2}}{4}$ (3.5)

holds, then there exists no linear operator $B\in L(X,Y)$ which fulfills

$\hat{x}=\mathrm{arg}\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha
R(x).$ (3.6)

The proof can be found in the appendix A.3. For given $y^{\delta}\in Y$,
$\hat{x}\in X$ and a penalty term $R$, these lemmas state whether there exists
a linear forward operator $B\colon X\to Y$ such that $\hat{x}$ is the Tikhonov
solution w.r.t. $R$ of the inverse problem w.r.t. $y^{\delta}$. As a
consequence, we can write the ADP minimization problem with a much simpler
side constraint.

######  Theorem 3.3.

Let Assumption 2.1 hold. Then, for all $y^{\delta}\in Y$, $y^{\delta}\neq 0$,
$\alpha>0$ the formulation

$\displaystyle\begin{split}&\min_{x\in X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}\\\
\text{s.t.}\quad&\exists v\in\partial R(x):\quad\alpha\langle
v,x\rangle\leqslant\frac{\|y^{\delta}\|^{2}}{4}\end{split}$ (3.7)

is equivalent to the ADP-Problem (2.6).

###### Proof.

According to Lemma 3.1 and Lemma 3.2 there exists a linear operator $B\in
L(X,Y)$ such that

$\hat{x}=\mathrm{arg}\min_{x\in X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha
R(x)$ (3.8)

if and only if $\hat{x}$ fulfills the side constraint of (3.7). ∎

######  Remark 3.4.

For the standard Tikhonov penalty term $R(x)=\frac{1}{2}\|x\|^{2}$, it holds
$\partial R(x)=x$. In this case we get

$\displaystyle\begin{split}\min_{x\in X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}\\\
\text{s.t.}\quad\|x\|^{2}\leqslant\frac{\|y^{\delta}\|^{2}}{4\alpha}\end{split}$
(3.9)

as an equivalent formulation of the ADP problem (2.6). For
$r=\|y^{\delta}\|^{2}/(4\alpha)$, this equals the Ivanov regularization method

$\begin{split}\min_{x\in X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}\\\
\text{s.t.}\quad\|x\|^{2}\leqslant r.\end{split}$ (3.10)

As [36] shows, this method is in fact equivalent to the Tikhonov method

$\min_{x\in
X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}+\frac{\tilde{\alpha}}{2}\|x\|^{2}$ (3.11)

for some $\tilde{\alpha}$ dependent on $y^{\delta}$ and $r$. We note that the
Tikhonov parameter $\tilde{\alpha}$ may be equal to zero and in particular it
differs from the parameter $\alpha$ of the ADP problem (see section 3.2).

######  Remark 3.5.

There are also cases in which the side constraint of (3.7) defines a non-
convex feasible set. Then, the ADP problem is more difficult to solve. We give
a simple two-dimensional example with the penalty term
$R\colon\mathbb{R}^{2}\to[0,\infty)$,

$R(x_{1},x_{2})=\begin{cases}3\cdot|x_{1}-5|&\text{for
}3\cdot|x_{1}-5|\geqslant|x_{2}|,\\\ |x_{2}|&\text{for
}|x_{2}|>3\cdot|x_{1}-5|.\end{cases}$ (3.12)

This functional has a non-centered minimum at $(5,0)^{\mathrm{T}}$ and the
absolute value of its gradient $|\partial R(x)|$ is strongly dependent on the
direction. Because of these properties, it’s easy to show that the term
$\langle v,x\rangle$, $v\in\partial R(x)$ in the side constraint of (3.7) is
non-convex w.r.t. $x\in\mathbb{R}^{2}$.

###  3.2 Parameter choice and early stopping

By construction of the ADP model, we expect it in application to act like DIP.
But in the previous section it turned out that ADP behaves in fact equivalent
to classical methods like Tikhonov’s. When we apply ADP to an inverse problem,
the question arises whether ADP can also deliver something that is “new” and
not equivalent to a Tikhonov solution. This section presents, how the model
has to be changed to produce ADP solutions that are more similar to DIP
solutions. In the same time, we derive a strategy for choosing the parameter
$\alpha$ of the ADP model.

When we compare the ADP method

$\begin{split}&\min_{B\in L(X,Y)}\frac{1}{2}\|Ax(B)-y^{\delta}\|^{2}\\\
\text{s.t.}\quad&x(B)=\mathrm{arg}\min_{x\in
X}\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\frac{\alpha_{\text{ADP}}}{2}\|x\|^{2}\end{split}$
(3.13)

to the equivalent (see Remark 3.4) Tikhonov method

$\min_{x\in
X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}+\frac{\tilde{\alpha}}{2}\|x\|^{2},$ (3.14)

we have to make sure not to confuse the parameters $\alpha_{\text{ADP}}$ and
$\tilde{\alpha}$ of both models with each other. At first we state the
following relation between these parameters.

######  Lemma 3.6.

If the solutions of (3.13) and (3.14) coincide,
$\tilde{\alpha}\leqslant\alpha_{\text{ADP}}$ holds. Equality of the parameters
could only occur if $y^{\delta}$ was in the kernel of $A^{*}$ or a singular
vector of $A$.

The proof can be found in the appendix A.4. In general, we can assume that
$\tilde{\alpha}<\alpha_{\text{ADP}}$ holds. So in any application it makes
sense to choose the ADP parameter greater than one would choose the parameter
of a Tikhonov model. But independent of the parameter choice, the ADP solution
will always be equivalent to a Tikhonov solution (Remark 3.4). To make ADP
more similar to DIP, we apply early stopping [15, section 7.8]. This strategy
is often used in the application of DIP but wasn’t considered for the ADP
model yet.

For a given inverse problem, we could solve the ADP problem (3.13) with a
gradient descent algorithm w.r.t. the operator $B$ (see section 4 for details)
and terminate this iteration early. Taking $B^{0}=A$ as initial value leads by
definition to $x(B^{0})$ being equal to the Tikhonov solution w.r.t. the
parameter $\alpha_{\text{ADP}}$. We assume the iteration to converge
successfully towards the minimizer $\hat{x}$ of (3.13). Since $\hat{x}$ is
also the minimizer of (3.14), the limit of the iteration is also a Tikhonov
solution but w.r.t. the parameter $\tilde{\alpha}$. Because of
$\tilde{\alpha}<\alpha_{\text{ADP}}$, the starting solution $x(B^{0})$ is a
stronger regularized Tikhonov solution than the limit $\hat{x}$ of the
iteration (see figure 1).

Figure 1: Comparison of ADP solutions during gradient descent to the Tikhonov
method (orange: at the start of the gradient descent, red: by early stopping,
gray: limit of the gradient descent). The forward operator is an integration
like in (4.1) and the data $y^{\delta}$ is contaminated with gaussion noise
(PSNR=40). The Tikhonov parameter and the stopping criterion for ADP were
chosen a posteriori to achieve the most accurate reconstructions.

If we apply early stopping, we take some $x(B^{k})$ in between, which in
general does not equal a Tikhonov solution w.r.t. $A$ (see figure 1). This
strategy makes sense if we expect $x(B^{k})$ to be a better solution than (the
Tikhonov solutions) $x(B^{0})$ and $\hat{x}$. That could be the case if
$x(B^{0})$ is slightly over-regularized and $\hat{x}$ is slightly under-
regularized. Because then, the optimal regularization would lay in between.

Now, we come back to the parameter choice. If we have a criterion for
estimating a suitable Tikhonov parameter $\alpha_{\text{Tik}}$ for a given
inverse problem, we should try to choose $\alpha_{\text{ADP}}$ in a way that

$\tilde{\alpha}<\alpha_{\text{Tik}}<\alpha_{\text{ADP}}$ (3.15)

holds. Because then, $x(B^{0})$ will be slightly over-regularized and
$\hat{x}$ slightly under-regularized, as proposed.

In the example of figure 1, we see that the ADP solution $x(B^{k})$, obtained
with early stopping, is a better approximation for the ground truth
$x^{\dagger}$ than the most accurate Tikhonov solution, which corresponds to
$\alpha_{\text{Tik}}$. But this result is strongly dependent on the particular
inverse problem. The Tikhonov method is optimal for data that is normally
distributed. If the given distribution differs from that, it is theoretically
possible that ADP with early stopping produces a better solution than the
Tikhonov method.

Finally, we want to include the early stopping strategy directly into the ADP
model to be able to investigate its effect on the regularization of inverse
problems. Early stopping enforces the iterated variable to stay close to the
initial value. Because of $B^{0}=A$, we can expect $\|B^{k}-A\|$ to be small
for small $k$. This leads to using $\|B^{k}-A\|$ as an additional penalty term
in the ADP problem, which has a similar effect as early stopping [4, section
2.3], [33, section 4]. What comes out is the ADP-$\beta$ model

$\begin{split}\min_{B\in
L(X,Y)}\frac{1}{2}\|Ax(B)-y^{\delta}\|^{2}+\beta\|B-A\|^{2}\\\
\text{s.t.}\quad x(B)=\mathrm{arg}\min_{x\in
X}\,\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x).\end{split}$ (3.16)

###  3.3 Properties of ADP

The equivalence between ADP and the Ivanov method (with general convex penalty
term $R$), shown in section 3.1, allows to obtain some regularization
properties (existence, stability, convergence) for ADP. We suppose that
Assumption 2.1 holds. Besides the functional

$\tilde{R}(x)=\min_{v\in\partial R(x)}\langle v,x\rangle$ (3.17)

is assumed to be well-defined. Because then, the side constraint of (3.7) can
be formulated as

$\tilde{R}(x)\leqslant\frac{\|y^{\delta}\|^{2}}{4\alpha}.$ (3.18)

Due to $\langle v,x\rangle=R(x)+R^{*}(v)$ for $v\in\partial R(x)$, where
$R^{*}$ denotes the convex conjugated functional, coercivity of $R$ implies
coercivity of $\tilde{R}$.

######  Remark 3.7 (Existence).

There exists a solution of the ADP problem (2.6) if the functional
$\tilde{R}$, defined in (3.17), is weakly lower semicontinuous. This follows
from the equivalence theorem 3.3 and [37, Theorem 2.1] about the existence of
Ivanov solutions.

Uniqueness of solutions and stability w.r.t. the data $y^{\delta}$ is less
trivial. First, the right hand side of the side constraint (3.18) is dependent
on $y^{\delta}$, which isn’t the case for ordinary Ivanov problems. Secondly,
we know from Remark 3.5, that the constraint (3.18) does not always define a
convex feasible set. Nevertheless, for the special case
$R(x)=\frac{1}{2}\|x\|^{2}$ we can obtain a convenient stability result. In
this case, $\tilde{R}(x)=\|x\|^{2}$ is a strictly convex functional.
Additionally, if the given inverse problem is ill-posed, we can assume the ADP
solutions to fulfill the constraint (3.18) with equality. Under these
conditions, the following theorem provides stability of ADP.

######  Theorem 3.8 (Stability).

For $R(x)=\frac{1}{2}\|x\|^{2}$, let $(y_{k})\subset Y$ be a sequence with
${y_{k}\to\hat{y}\in Y}$ and assume that the corresponding ADP solutions
$x_{k}$, $\hat{x}$ are unique and fulfill the side constraint in (3.9) with
equality. Then ADP is stable, which means $x_{k}\to\hat{x}$.

The proof can be found in the appendix A.5.

To obtain a convergence result for ADP, it makes sense to use standard
convergence theorems, either of the Tikhonov method [21, Theorem 4.4] or of
the Ivanov method [23, Theorem 2.5], [31, Theorem 3]. They differ especially
in the source conditions they require for the ground truth $x^{\dagger}$ and
in the parameter choice rules. If we assume $\tilde{R}$ to be convex, by [36,
Theorem 2] and the equivalence theorem 3.3, the Tikhonov problem

$\min_{x\in X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}+\tilde{\alpha}\tilde{R}(x)$
(3.19)

is equivalent to the ADP formulation (3.7) for suitable chosen
$\tilde{\alpha}\geqslant 0$.

######  Remark 3.9 (Convergence).

Because of the equivalence between (3.7) and (3.19), the convergence of ADP
solutions $x_{\alpha}^{\delta}$ to $x^{\dagger}$ for vanishing $\delta$ w.r.t.
the Bregman distance can be directly derived from Tikhonov convergence
theorems. But the ADP parameter $\alpha$ does not coincide with the Tikhonov
parameter $\tilde{\alpha}$. That’s why, for ADP we do not get an explicit
parameter choice rule like $\alpha\sim\delta$. Besides, a source condition for
$x^{\dagger}$ has to be fulfilled by the functional $\tilde{R}$ (defined in
(3.17)) and not by the penalty term $R$.

###  3.4 Properties of ADP-

For proving the existence of solutions of variational regularization schemes,
[21, Theorem 3.1] provides a useful framework. If we want to apply this for
ADP-$\beta$, it has to be ensured that $B\mapsto x(B)$ is weak-weak continuous
[21, Assumptions 2.1]. But unfortunately, in general this is not the case.

To obtain convenient regularization properties anyway, we restrict to
${X=Y=L^{2}(\Omega)}$ with $\Omega\subset\mathbb{R}^{n}$. In this setting, we
consider a forward operator $A\colon X\to Y$ that can be parametrized by a
function $f\in L^{p}(\Omega)$, $p\in[1,\infty)$. More precisely, we take a
continuous, bilinear operator ${T\colon L^{p}(\Omega)\times X\to Y}$ and
define

$Ax=T(f,x).$ (3.20)

The same parametrization of operators by functions is used in [5]. One typical
example would be a convolutional operator $T(f,x)=f*x$.

The crucial idea is the additional restriction $f\in W^{1,p}(\Omega)$ to take
advantage of the compact embedding of Sobolev spaces $W^{1,p}(\Omega)\subset
L^{p}(\Omega)$. A similar strategy is used in [24] for achieving weak-weak
continuity of the forward operator.

We define the parametrized ADP-$\beta$ approach as

$\displaystyle\begin{split}\min_{g\in
W^{1,p}}\frac{1}{2}\|T(f,x_{g})-y^{\delta}\|_{L^{2}}^{2}&+\beta\|f-g\|_{W^{1,p}}^{2}\\\
\text{s.t.}\quad x_{g}=\mathrm{arg}\min_{x\in
L^{2}}\,\frac{1}{2}\|T(g,x)-y^{\delta}\|_{L^{2}}^{2}&+\alpha R(x).\end{split}$
(3.21)

In particular, this can be interpreted as a Tikhonov method for solving the
nonlinear inverse problem $F(g^{\dagger})=y^{\dagger}$ with the forward
operator $F\colon W^{1,p}(\Omega)\to Y$, $F(g)=T(f,x_{g})$.

######  Remark 3.10 (Existence).

The forward operator $F$ is weak-strong continuous if the penalty term $R$ is
strongly convex. This holds, because weak convergence $g_{k}\rightharpoonup g$
w.r.t. $W^{1,p}(\Omega)$ implies convergence by norm in $L^{p}(\Omega)$, by
Theorem 2.3 the convergence of $x_{g_{k}}\to x_{g}$ follows, and the bilinear
operator $T$ is continuous. This is more than enough to fulfill the
assumptions of [21, Theorem 3.1], which provides the existence of a solution
of (3.21).

A weak stability result for the parametrized ADP-$\beta$ method could be
directly obtained from [21, Theorem 3.2]. But this particular framework even
allows to prove strong stability.

######  Theorem 3.11 (Stability).

For $p=2$, let $R$ be a strongly convex penalty term, $(y_{k})\subset Y$ a
convergent sequence with $y_{k}\to\hat{y}$ and $(g_{k})\subset
W^{1,p}(\Omega)$ the corresponding solutions of the ADP-$\beta$ problem
(3.21). Then, $(g_{k})$ has a convergent subsequence and the limit of each
subsequence is an ADP-$\beta$ solution corresponding to $\hat{y}$.

The proof can be found in the appendix A.6.

While proving existence and stability of ADP-$\beta$-solutions required a
smart parametrization and the use of compact embeddings, a convergence theorem
(w.r.t. the Bregman distance) can be proven for the general formulation (2.7).
Similar to classical results like [21, Theorem 4.4] or [6, Theorem 2], we need
to assume a source condition

$\exists w\in Y:\quad A^{*}w\in\partial R(x^{\dagger}).$ (3.22)

The parameter $\beta$ turns out to be really helpful for obtaining a
convergence result.

######  Theorem 3.12 (Convergence).

Let Assumption 2.1 hold, $x^{\dagger}$ be an $R$-minimizing solution of (2.1)
which fulfills the source condition (3.22) and assume there exist ADP-$\beta$
solutions $\hat{x}_{\alpha}^{\delta}$ of (2.7). If $\alpha$ is chosen
proportional to $\delta$, there exists $d\in
D_{R}(\hat{x}_{\alpha}^{\delta},x^{\dagger})$ which fulfills $d=O(\delta)$.

The proof can be found in the appendix A.7.

##  4 Numerical Computations

Aim. We want to see whether there is a similarity between ADP and DIP also on
the numerical side. The ADP approach is based on the idea of using a LISTA
network in a DIP method. Usually LISTA architectures contain round about ten
layers, but ADP is motivated with a network of infinite depth (see section 2).
To derive the ADP model, the output of this infinite network is then replaced
by the solution of a minimization problem. So the question arises, whether
numerically computed ADP solutions of an inverse problem are yet similar to
solutions obtained via a DIP with LISTA architecture.

In this section, we present algorithms for the computation of ADP solutions
and we compare them with DIP solutions. In doing so, the focus is not on the
performance of the methods (in comparison to other state-of-the-art
reconstruction algorithms) but on the similarity of the different solutions.

Methods. From Theorem 3.3, we know that the ADP problem is equivalent to an
Ivanov problem. This creates a possibility to compute ADP solutions easily,
fast and almost exactly (we call this method ADP Ivanov). In contrast to that,
it is more difficult to realize a LISTA architecture with infinite depth. But
there are at least two possibilities to simulate such a network.

The first idea (Algorithm 1: DIP LISTA $L=\infty$) is to begin with a network
$\varphi_{B}$ of ten layers and to increase the network depth during the
training process of the DIP. This is done implicitly with a simple trick. In
each training step, the network’s input is set to be the network’s output of
the previous step [13, Appendix 3]. So the original input will pass through
more and more layers and in each step the last ten layers are optimized (via
backpropagation).

initialize $B_{0}$, $z_{0}$ (e.g. $B_{0}=A$, $z_{0}=\text{random noise}$);

for  _$k=0,1,...$_ do

$z_{k+1}=\varphi_{B_{k}}(z_{k})$;

$\mathrm{loss}_{k}=\frac{1}{2}\|A\varphi_{B_{k}}(z_{k})-y^{\delta}\|^{2}$;

$B_{k+1}=\mathrm{update}(\nabla_{B_{k}}\mathrm{loss}_{k})$;

end for

return $z_{k}$

Algorithm 1 DIP LISTA L=$\infty$

initialize $B_{0}$ (e.g. $B_{0}=A$);

for  _$k=0,1,...$_ do

1\. Calculate $x(B_{k})$ with fixed point iteration

2\. IFT provides: $\nabla_{B_{k}}x(B_{k})$;

3\. Update:

$\mathrm{loss}_{k}=\frac{1}{2}\|Ax(B_{k})-y^{\delta}\|^{2}$;

$B_{k+1}=\mathrm{update}(\nabla_{B_{k}}\mathrm{loss}_{k})$;

end for

return $x(B_{k})$

Algorithm 2 ADP IFT

The second idea (Algorithm 2: ADP IFT) is to compute $x(B)$ from (2.6) with a
classical algorithm like ISTA. After that, one can compute the gradient of
$x(B)$ w.r.t. $B$ (see the proof of [13, Lemma 4.1]) via the implicit function
theorem (IFT). Thus, backpropagation through a big amount of layers is
avoided.

For the standard DIP approach, we use a LISTA-like architecture of ten Layers
(DIP LISTA $L=10$) and optimize the weights via backpropagation. So, in total
we compare four different methods (ADP Ivanov, ADP IFT, DIP LISTA $L=\infty$,
DIP LISTA, $L=10$). Since solving the Ivanov problem results in the exact ADP
solution, we use this as a reference for the other three methods (for which we
don’t have convergence guarantees).

In all methods we use the elastic net functional [38]
$R(x)=\alpha_{1}\|x\|_{1}+\frac{\alpha_{2}}{2}\|x\|^{2}$ as a penalty term. So
there is one parameter for $\ell^{1}$-regularization (leads to sparsity) and
one parameter for $\ell^{2}$-regularization (leads to stability and
smoothness). In the LISTA-architecture, this is realized by substracting the
gradient of the $\ell^{2}$-term before applying the activation function.

Setting. We consider two different artificial inverse problems (inversion of
the integration operator and a deconvolution) on $L^{2}(I)$ for an interval
$I\subset\mathbb{R}$. The forward operators are

$(A_{1}x)(t)=\int_{0}^{t}x(s)\,\mathrm{d}s\quad\text{ and }\quad A_{2}x=g*x,$
(4.1)

$g$ being a Gaussian function. Both of them lead to ill-posed inverse
problems. We chose three different ground truth functions and created data by
applying the forward operators and adding normally distributed random noise.
This leads to six examples in total, which is enough for some basic
observations. Figure 2 shows the reconstructions corresponding to the
integration operator $A_{1}$. The three rows contain the three different
ground truth functions and each column contains a different method. For
comparison, the actual ADP solution (ADP Ivanov) and the ground truth is
displayed in every plot. Since we are only interested in finding similarities
and disparities between the solutions of the different methods, the choice of
the regularization parameters plays a minor role. So, we took the same values
$\alpha_{1}$, $\alpha_{2}$ for each method and simply chose them a posteriori
for each example to minimize the $L^{2}$-error between reconstructions and
ground truth. Figure 3 shows the analogous results for the deconvolution
problem (forward operator $A_{2}$).

Figure 2: Computation of ADP and DIP reconstructions via the IFT, via a DIP
with LISTA network ($L=\infty$ and $L=10$) and via the equivalent Ivanov
problem. The forward operator is $A_{1}$ (integration). The given data has
PSNR=40 due to additiv Gaussian noise. The regularization parameters
$\alpha_{1}$, $\alpha_{2}$ are chosen for each example (row) separately but
are the same for each method (column).

Observations. From these experiments, we can make the following observations.
There is a significant difference between using $L=10$ or $L=\infty$ layers in
a LISTA network. With an infinite number of layers, the reconstructions are
looking more realistic. The results of the DIP LISTA L=$\infty$ (Algorithm 1)
method and of the IFT method (Algorithm 2) are always looking quite similar.
This was expected because both of them simulate an infinitely deep LISTA
network. Differences are probably due to the different ways the gradients are
computed or due to to slow convergence of the methods.

Figure 3: Computation of ADP and DIP reconstructions via the IFT, via a DIP
with LISTA network ($L=\infty$ and $L=10$) and via the equivalent Ivanov
problem. The forward operator is $A_{2}$ (convolution). The given data has
PSNR=45 due to additiv Gaussian noise. The regularization parameters
$\alpha_{1}$, $\alpha_{2}$ are chosen for each example (row) separately but
are the same for each method (column).

In most of the cases, the reconstructions of these both methods are looking
quite similar to the actual ADP solution. But sometimes they contain artifacts
(e.g. the peaks in figure 3, third row). It seems that there are some spots
which are hard to reconstruct for the DIP methods and others are rather
simple. Besides, the ADP problem (2.6) is not a convex minimization problem
w.r.t. $B$. So there is no guarantee for the methods which do gradient descent
(DIP LISTA $L=\infty$ and the IFT method) to converge towards the global
minimizer. Figure 4 shows that the reconstructions of these methods are indeed
dependent on the initial value $B_{0}$ of the algorithms. In contrast to that,
the Ivanov problem from Theorem 3.3 is convex (with the elastic net penalty
term $R$). That’s probably why the actual ADP solutions are the only ones
which never contain strange artifacts and the only ones that are always quite
good reconstructions of the ground truth.

Figure 4: The same setting and methods as in figure 3 but with different
initial values $B_{0}$.

The easiest possibility to slightly improve the reconstrution quality is to
apply early stopping. In doing so, the most severe artifacts in the
reconstructions can be diminished. This case corresponds to the ADP-$\beta$
approach (see section 3.2), whose additional convex term $\beta\|B-A\|^{2}$ is
a numerical advantage because it stabilizes the gradient descent for finding
the minimizer. Indeed, adding the gradient of the $\beta$-term to the update
in Algorithm 2 (ADP IFT) can also diminish the severe artifacts. But we do not
include experimental results about this, since the most interesting part is
the comparison with the equivalent Ivanov problem, which doesn’t exist for
ADP-$\beta$.

The main conclusion is the numerical verification of the derivation of the ADP
problem from the DIP approach. It is possible to use the theoretical analysis
of the ADP problem for interpretations of the DIP approach because of the
similarity between the reconstructions from the different numerical methods.
However, the examples from figure 4 illustrate that DIP can be formulated as a
minimization problem (2.3) but a numerical computed solution is not
automatically a global minimizer of this problem. If early stopping is used,
it is probably not even a local minimizer. Hence, there is a significant
difference between the theoretical definition and the practical implementation
of DIP.

##  5 Conclusion

ADP and ADP-$\beta$ were introduced as methods for solving ill-posed inverse
problems in a typical Hilbert space setting (Assumption 2.1). Both of them are
motivated by considering DIP with a LISTA-like architecture. The main result
is an equivalence of ADP to the classical method of Ivanov regularization.

We have proven existence, stability and convergence results for both ADP and
ADP-$\beta$. The obtained regularization properties are comparable to the ones
of classical methods like Tikhonov’s. In principal, these results can be
transferred to DIP with LISTA-like networks. But due to non-convexity of the
DIP minimization problem, numerically computed DIP solutions can differ
significantly from exact ADP solutions, although they are similar in many
cases. We conclude that theoretical analyses of the DIP approach should
consider the whole optimization process and not only the properties of the
minimizer.

One very important part is the early stopping of the DIP optimization process.
In the ADP setting, we incorporated this strategy with an additional penalty
term, which resulted in the ADP-$\beta$ model. The effect of this
regularization can be seen by comparing the convergence theorems of ADP and
ADP-$\beta$. Theorem 3.12 provides a parameter choice rule
($\alpha\sim\delta$) for ADP-$\beta$, which is a big advantage over ADP.

A generalization of the ADP regularization results to DIP with general
convolutional neural networks (CNNs) would be very desirable. The LISTA
architecture was suitable because of its similarity to proximal splitting
algorithms and the possibility to interpret the output as a solution of a
variational problem. Finding similar connections for general CNNs is harder.
However in [7], CNNs are used to model proximal mappings and in [28], CNNs are
interpreted as algorithms for sparse coding. Besides, [10] asserts that most
common activation functions are in fact proximal mappings and they establish a
theory for characterizing the fixed point sets of neural networks as solutions
of variational inequalities. These directions could provide ideas for possible
future extensions.

## Acknowledgements

I want to thank Dr. Daniel Otero Baguer, Prof. Peter Maaß, Dr. Tobias Kluth
and many more colleagues from the University of Bremen and Dr. Yury Korolev
from the University of Cambridge for helpful advice and feedback.

##  Appendix A Proofs of theoretical results

###  A.1 Theorem 2.3

Continuity of $B\mapsto x(B)$.

###### Proof.

Let $(B_{k})\subset L(X,Y)$ be a sequence of operators with $B_{k}\to B\in
L(X,Y)$. At first, we mention that the sequence $(x(B_{k}))$ is bounded
because

$\begin{split}\alpha
R(x(B_{k}))\leqslant\frac{1}{2}\|B_{k}x(B_{k})-y^{\delta}\|^{2}+\alpha
R(x(B_{k}))\leqslant\frac{1}{2}\|B_{k}x(B)-y^{\delta}\|^{2}+\alpha
R(x(B))\end{split}$ (A.1)

holds. Further, we can estimate

$\displaystyle\begin{split}&\quad\,\,\frac{1}{2}\|Bx(B_{k})-y^{\delta}\|^{2}+\alpha
R(x(B_{k}))\\\
&=\frac{1}{2}\|B_{k}x(B_{k})-y^{\delta}+(B-B_{k})\,x(B_{k})\|^{2}+\alpha
R(x(B_{k}))\\\
&\leqslant\frac{1}{2}(\|B_{k}x(B_{k})-y^{\delta}\|+\|(B-B_{k})\,x(B_{k})\|)^{2}+\alpha
R(x(B_{k}))\\\ &=\frac{1}{2}\|B_{k}x(B_{k})-y^{\delta}\|^{2}+\alpha
R(x(B_{k}))\\\
&\,+\|(B-B_{k})\,x(B_{k})\|\cdot\left(\|B_{k}x(B_{k})-y^{\delta}\|+\frac{1}{2}\|(B-B_{k})\,x(B_{k})\|\right).\end{split}$
(A.2)

Because of the boundedness of $(x(B_{k}))$ and the convergence $B_{k}\to B$,
the term

$\|(B-B_{k})\,x(B_{k})\|\cdot\left(\|B_{k}x(B_{k})-y^{\delta}\|+\frac{1}{2}\|(B-B_{k})\,x(B_{k})\|\right)$
(A.3)

converges to zero. For the remaining terms, we can estimate

$\frac{1}{2}\|B_{k}x(B_{k})-y^{\delta}\|^{2}+\alpha
R(x(B_{k}))\leqslant\frac{1}{2}\|B_{k}x(B)-y^{\delta}\|^{2}+\alpha R(x(B)),$
(A.4)

and $\|B_{k}x(B)-y^{\delta}\|^{2}$ converges to $\|Bx(B)-y^{\delta}\|^{2}$. So
$(x(B_{k}))$ is a minimizing sequence of the strongly convex functional
$\frac{1}{2}\|Bx-y^{\delta}\|^{2}+\alpha R(x)$ because $x(B)$ is the
minimizer. By [27, Theorem 1], the minimizing sequence converges to the
minimizer $x(B)$. ∎

###  A.2 Lemma 3.1

First part of the equivalence theorem for ADP to Ivanov problems.

###### Proof.

Let $\hat{x},v,y^{\delta},\alpha$ and $R$ be given according to the
assumptions. We have to find a linear operator $B$ such that

$-B^{*}(B\hat{x}-y^{\delta})\in\alpha\partial R(\hat{x}).$ (A.5)

holds. Because of $v\in\partial R(\hat{x})$, we just try to solve the equation

$-B^{*}(B\hat{x}-y^{\delta})=\alpha v$ (A.6)

for $B$. If $\hat{x}=0$, solving would be trivial. Otherwise, we can decompose
$v$ into

$v=\mu\hat{x}+v_{\bot}\qquad\text{s.t. }\langle v_{\bot},\hat{x}\rangle=0.$
(A.7)

Accordingly it is $\mu={\langle v,\hat{x}\rangle}/{\|\hat{x}\|^{2}}$. With
that, we can write the equation from above as

$B^{*}B\hat{x}+\alpha\mu\hat{x}+\alpha v_{\bot}=B^{*}y^{\delta}.$ (A.8)

We consider a linear operator $B\colon X\to Y$ of the form

$Bx=(\sigma_{1}\langle x,\hat{x}\rangle+\sigma_{2}\langle
x,v_{\bot}\rangle)\cdot y^{\delta}$ (A.9)

with two coefficients $\sigma_{1}$ and $\sigma_{2}$ to be determined later.
Then, the adjoint operator is given by

$B^{*}y=\langle y,y^{\delta}\rangle(\sigma_{1}\hat{x}+\sigma_{2}v_{\bot})$
(A.10)

and it holds

$\displaystyle B^{*}B\hat{x}$
$\displaystyle=B^{*}((\sigma_{1}\|\hat{x}\|^{2})\cdot
y^{\delta})=\sigma_{1}^{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}\hat{x}+\sigma_{1}\sigma_{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}v_{\bot},$
(A.11) $\displaystyle B^{*}y^{\delta}$
$\displaystyle=\sigma_{1}\|y^{\delta}\|^{2}\hat{x}+\sigma_{2}\|y^{\delta}\|^{2}v_{\bot}.$
(A.12)

To fulfill (A.8), we have to solve

$\sigma_{1}^{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}\hat{x}+\sigma_{1}\sigma_{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}v_{\bot}+\alpha\mu\hat{x}+\alpha
v_{\bot}=\sigma_{1}\|y^{\delta}\|^{2}\hat{x}+\sigma_{2}\|y^{\delta}\|^{2}v_{\bot}.$
(A.13)

Because $\hat{x}$ and $v_{\bot}$ are orthogonal to each other, we get the two
equations

$\displaystyle\sigma_{1}^{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}+\alpha\mu=\sigma_{1}\|y^{\delta}\|^{2},$
(A.14)
$\displaystyle\sigma_{1}\sigma_{2}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}+\alpha=\sigma_{2}\|y^{\delta}\|^{2}.$
(A.15)

Notice that (A.15) and the coefficient $\sigma_{2}$ could be ignored if
$v_{\bot}=0$ held.

Equation (A.14) can be solved for $\sigma_{1}$ with a quadratic formula, which
leads to

$\sigma_{1}=\frac{1}{2\|\hat{x}\|^{2}}\pm\sqrt{\frac{1}{4\|\hat{x}\|^{4}}-\frac{\alpha\mu}{\|\hat{x}\|^{2}\|y^{\delta}\|^{2}}}.$
(A.16)

Accordingly,

$\frac{\alpha\mu}{\|y^{\delta}\|^{2}}\leqslant\frac{1}{4\|\hat{x}\|^{2}}$
(A.17)

must hold to get real solutions. We know from above that $\mu={\langle
v,\hat{x}\rangle}/{\|\hat{x}\|^{2}}$. If we insert this, we will see that this
matches exactly the assumptions of the lemma.

Now, equation (A.15) has to be solved for $\sigma_{2}$. Excluding $\sigma_{2}$
leads to

$\sigma_{2}(\sigma_{1}\|\hat{x}\|^{2}\|y^{\delta}\|^{2}-\|y^{\delta}\|^{2})+\alpha=0.$
(A.18)

If the term inside of the parenthesis doesn’t equal zero, there will exist a
solution $\sigma_{2}$. If the term equaled zero, equation (A.14) would lead to
$\mu=0$. But in this case, we could choose $\sigma_{1}=0$ (the quadratic
formula allows two solutions), and then it is no problem to find a solution
for $\sigma_{2}$, too.

By finding solutions for $\sigma_{1}$ and $\sigma_{2}$, we showed that the
operator $B$ defined in (A.9) solves equation (A.6). So the lemma is proved. ∎

###  A.3 Lemma 3.2

Second part of equivalence theorem for ADP to Ivanov problems.

###### Proof.

Let $\hat{x},y^{\delta},\alpha$ and $R$ be given according to the assumptions.
Assume there exists a linear operator $B$ such that

$0\in B^{*}(B\hat{x}-y^{\delta})+\alpha\partial R(\hat{x})$ (A.19)

holds. It follows

$v:=-\frac{1}{\alpha}B^{*}(B\hat{x}-y^{\delta})\in\partial R(\hat{x}).$ (A.20)

We can calculate $\alpha\langle v,\hat{x}\rangle=-\|B\hat{x}\|^{2}+\langle
y^{\delta},B\hat{x}\rangle$. So according to the assumptions,

$-\|B\hat{x}\|^{2}+\langle
y^{\delta},B\hat{x}\rangle>\frac{\|y^{\delta}\|^{2}}{4}$ (A.21)

must hold. But with Young’s inequality, we get

$-\|B\hat{x}\|^{2}+\langle
y^{\delta},B\hat{x}\rangle\leqslant-\|B\hat{x}\|^{2}+\frac{1}{4}\|y^{\delta}\|^{2}+\|B\hat{x}\|^{2}=\frac{\|y^{\delta}\|^{2}}{4}.$
(A.22)

Obviously, this is a contradiction. That’s why such an operator $B$ can’t
exist. ∎

###  A.4 Lemma 3.6

Relation between the ADP parameter and the Tikhonov parameter of the
equivalent problem.

###### Proof.

Let $\hat{x}$ be the solution of the ADP problem (3.13). Because of the
equivalence to the Tikhonov method, $\hat{x}$ is the solution of (3.14) in the
same time. Besides, $x(A)$ is the Tikhonov solution w.r.t. the parameter
$\alpha_{\text{ADP}}$ of the inverse problem. Because of the minimizing
properties of $\hat{x}$ and $x(A)$,

$\displaystyle\frac{1}{2}\|A\hat{x}-y^{\delta}\|^{2}$
$\displaystyle\leqslant\frac{1}{2}\|Ax(A)-y^{\delta}\|^{2},$ (A.23)
$\displaystyle\frac{1}{2}\|Ax(A)-y^{\delta}\|^{2}+\frac{\alpha_{\text{ADP}}}{2}\|x(A)\|^{2}$
$\displaystyle\leqslant\frac{1}{2}\|A\hat{x}-y^{\delta}\|^{2}+\frac{\alpha_{\text{ADP}}}{2}\|\hat{x}\|^{2}$
(A.24)

holds. It follows $\|x(A)\|^{2}\leqslant\|\hat{x}\|^{2}$. Both $x(A)$ and
$\hat{x}$ are Tikhonov solutions of the same problem (only with different
parameters). So $\tilde{\alpha}\leqslant\alpha_{\text{ADP}}$ must hold because
the norm of $\hat{x}$ is greater (or equal) than the norm of $x(A)$.

Now, we assume $\tilde{\alpha}=\alpha_{\text{ADP}}>0$. By Remark 3.4, the
problems

$\displaystyle\min_{x\in
X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}+\frac{\alpha_{\text{ADP}}}{2}\|x\|^{2},$
(A.25) $\displaystyle\min_{x\in
X}\frac{1}{2}\|Ax-y^{\delta}\|^{2}\quad\text{s.t.
}\|x\|^{2}\leqslant\frac{\|y^{\delta}\|^{2}}{4\alpha_{\text{ADP}}}$ (A.26)

are equivalent. The solution $\hat{x}$ fulfills

$(A^{*}A+\alpha_{\text{ADP}}\cdot\mathrm{Id})\hat{x}=A^{*}y^{\delta}$ (A.27)

and we assume $A^{*}y^{\delta}\neq 0$. Then, $\hat{x}$ must fulfill the side
constraint of (A.26) with equality, otherwise $A^{*}A\hat{x}=A^{*}y^{\delta}$
would hold, which is a contradiction. Accordingly we get
$\alpha_{\text{ADP}}\|\hat{x}\|^{2}=\|y^{\delta}\|^{2}/4$ and by computing the
inner product of (A.27) with $\hat{x}$, it follows

$\|A\hat{x}\|^{2}+\frac{\|y^{\delta}\|^{2}}{4}=\|A\hat{x}\|^{2}+\alpha_{\text{ADP}}\|\hat{x}\|^{2}=\langle
A\hat{x},y^{\delta}\rangle.$ (A.28)

If we then apply the Cauchy-Schwarz and Young’s inequality, we get

$\langle
A\hat{x},y^{\delta}\rangle\leqslant\|A\hat{x}\|\cdot\|y^{\delta}\|\leqslant\|A\hat{x}\|^{2}+\frac{\|y^{\delta}\|^{2}}{4},$
(A.29)

which means these inequalities must in fact hold as equalities. Therefore,
$A\hat{x}$ and $y^{\delta}$ must be linear dependent (Cauchy-Schwarz) and
$2\|A\hat{x}\|=\|y^{\delta}\|$ must hold (Young). It follows

$A\hat{x}=\frac{1}{2}y^{\delta}.$ (A.30)

We can plug this into (A.27) and get

$\alpha_{\text{ADP}}\hat{x}=\frac{1}{2}A^{*}y^{\delta}.$ (A.31)

Accordingly $AA^{*}y^{\delta}=\alpha_{\text{ADP}}y^{\delta}$ holds, so
$y^{\delta}$ is a singular vector of $A$. ∎

###  A.5 Theorem 3.8

Stability of the ADP approach.

###### Proof.

We follow some of the ideas of the proofs of [14, Theorem 2.1] and [31,
Theorem 2].

Let $x_{k}$ and $\hat{x}$ be unique solutions of (3.9) for
$y^{\delta}=y_{k},\hat{y}$ with $y_{k}\to\hat{y}$. The sequence $(x_{k})$ is
bounded, so there exists a weakly convergent subsequence $(x_{k_{l}})$,
$x_{k_{l}}\rightharpoonup x_{\infty}$. For arbitrary $\varepsilon>0$ and $x\in
X$ with $\|x\|^{2}\leqslant\|\hat{y}\|^{2}\cdot(4\alpha)^{-1}-\varepsilon$, it
holds

$\|Ax_{\infty}-\hat{y}\|\leqslant\liminf_{l\to\infty}\|Ax_{k_{l}}-y_{k_{l}}\|\leqslant\lim_{l\to\infty}\|Ax-
y_{k_{l}}\|=\|Ax-\hat{y}\|$ (A.32)

because $x_{k_{l}}$ minimizes the ADP problem w.r.t. $y_{k_{l}}$ and $x$
fulfills the side constraint for $l$ big enough. With $\varepsilon\to 0$ and
because of the uniqueness of the solutions, we obtain $x_{\infty}=\hat{x}$.
Arguing with a subsequence of a subsequence leads to the weak convergence
$x_{k}\rightharpoonup\hat{x}$ of the whole sequence.

According to the assumptions, it holds
$\|x_{k}\|^{2}=\|\hat{y_{k}}\|^{2}\cdot(4\alpha)^{-1}$. So $y_{k}\to\hat{y}$
implies ${\|x_{k}\|\to\|\hat{x}\|}$ and together with the weak convergence, we
finally obtain $x_{k}\to\hat{x}$. ∎

###  A.6 Theorem 3.11

Stability of the ADP-$\beta$ approach.

###### Proof.

First, we note that the sequence $(g_{k})$ is bounded in $W^{1,2}(\Omega)$.
Hence, there exists at least one weakly convergent subsequence. For any
subsequence with $g_{k}\rightharpoonup\hat{g}$, it holds

$\displaystyle T(f,x_{g_{k}})-y_{k}\,$
$\displaystyle\to\,T(f,x_{\hat{g}})-\hat{y}$ (A.33)

because of the arguments from Remark 3.10. For arbitrary $g\in
W^{1,2}(\Omega)$,

$\begin{split}&\quad\,\frac{1}{2}\|T(f,x_{\hat{g}})-\hat{y}\|_{L^{2}}^{2}+\beta\|\hat{g}-f\|_{W^{1,2}}^{2}\\\
&\leqslant\liminf_{k\to\infty}\frac{1}{2}\|T(f,x_{g_{k}})-y_{k}\|_{L^{2}}^{2}+\beta\|g_{k}-f\|_{W^{1,2}}^{2}\\\
&\leqslant\lim_{k\to\infty}\frac{1}{2}\|T(f,x_{g})-y_{k}\|_{L^{2}}^{2}+\beta\|g-f\|_{W^{1,2}}^{2}\\\
&=\frac{1}{2}\|T(f,x_{g})-\hat{y}\|_{L^{2}}^{2}+\beta\|g-f\|_{W^{1,2}}^{2}\end{split}$
(A.34)

holds because of the minimizing property of $g_{k}$ w.r.t. $y_{k}$. Hence,
$\hat{g}$ is a minimizer of (3.21) w.r.t $\hat{y}$. If we choose $g=\hat{g}$,
the first and the last line in (LABEL:eq_stab_minimizer) coincide, and we get

$\lim_{k\to\infty}\frac{1}{2}\|T(f,x_{g_{k}})-y_{k}\|_{L^{2}}^{2}+\beta\|g_{k}-f\|_{W^{1,2}}^{2}=\frac{1}{2}\|T(f,x_{\hat{g}})-\hat{y}\|_{L^{2}}^{2}+\beta\|\hat{g}-f\|_{W^{1,2}}^{2}.$
(A.35)

If follows
$\lim_{k\to\infty}\|g_{k}-f\|_{W^{1,2}}^{2}=\|\hat{g}-f\|_{W^{1,2}}^{2}$.
Hence, $(g_{k})$ converges by norm to $\hat{g}$. ∎

###  A.7 Theorem 3.12

Convergence of the ADP-$\beta$ approach.

###### Proof.

According to (3.22), we can choose
$d=R(\hat{x}_{\alpha}^{\delta})-R(x^{\dagger})-\langle
A^{*}w,\hat{x}_{\alpha}^{\delta}-x^{\dagger}\rangle$ and there exists an
operator ${\hat{B}\in L(X,Y)}$ that fulfills
$\hat{x}_{\alpha}^{\delta}=x(\hat{B})$.

Because of the minimizing property of $\hat{x}_{\alpha}^{\delta}$,

$\alpha
R(\hat{x}_{\alpha}^{\delta})\leqslant\frac{1}{2}\|\hat{B}\hat{x}_{\alpha}^{\delta}-y^{\delta}\|^{2}+\alpha
R(\hat{x}_{\alpha}^{\delta})\leqslant\frac{1}{2}\|\hat{B}x^{\dagger}-y^{\delta}\|^{2}+\alpha
R(x^{\dagger}).$ (A.36)

holds. If follows

$\displaystyle\begin{split}d&=R(\hat{x}_{\alpha}^{\delta})-R(x^{\dagger})-\langle
A^{*}w,\hat{x}_{\alpha}^{\delta}-x^{\dagger}\rangle\leqslant\frac{1}{2\alpha}\|\hat{B}x^{\dagger}-y^{\delta}\|^{2}-\langle
w,A\hat{x}_{\alpha}^{\delta}-y^{\dagger}\rangle\\\
&\leqslant\frac{1}{2\alpha}\left(\|\hat{B}x^{\dagger}-Ax^{\dagger}\|+\|y^{\dagger}-y^{\delta}\|\right)^{2}+\|w\|\|A\hat{x}_{\alpha}^{\delta}-y^{\dagger}\|\\\
&\leqslant\frac{1}{2\alpha}\left(\|x^{\dagger}\|\|\hat{B}-A\|+\delta\right)^{2}+\|w\|\|A\hat{x}_{\alpha}^{\delta}-y^{\dagger}\|.\end{split}$
(A.37)

We will show $\|\hat{B}-A\|=O(\delta)$ and
$\|A\hat{x}_{\alpha}^{\delta}-y^{\dagger}\|=O(\delta)$ to deduce $d=O(\delta)$
for $\alpha$ chosen proportional to $\delta$. Because of the minimizing
property of $\hat{B}$, we get

$\beta\cdot\|\hat{B}-A\|^{2}\leqslant\frac{1}{2}\|Ax(\hat{B})-y^{\delta}\|^{2}+\beta\cdot\|\hat{B}-A\|^{2}\leqslant\frac{1}{2}\|Ax(A)-y^{\delta}\|^{2}.$
(A.38)

From standard convergence results of the Tikhonov method [21, Theorem 4.4] or
[6, Theorem 2], we get $\|Ax(A)-y^{\delta}\|=O(\delta)$. So
$\|\hat{B}-A\|=O(\delta)$ holds.

Besides,

$\displaystyle\|A\hat{x}_{\alpha}^{\delta}-y^{\dagger}\|\leqslant\|Ax(\hat{B})-y^{\delta}\|+\|y^{\delta}-y^{\dagger}\|\leqslant\|Ax(A)-y^{\delta}\|+\delta$
(A.39)

holds and we can use $\|Ax(A)-y^{\delta}\|=O(\delta)$ again. So $d=O(\delta)$
follows. ∎

## References

  * [1] G. S. Alberti, E. De Vito, M. Lassas, L. Ratti, and M. Santacesaria.  Learning the optimal Tikhonov regularizer for inverse problems.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [2] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [3] D. O. Baguer, J. Leuschner, and M. Schmidt.  Computed tomography reconstruction using deep image prior and learned reconstruction methods.  Inverse Problems, 36(9):094004, 2020.
  * [4] C. Bishop.  Regularization and complexity control in feed-forward networks, pages 141–148.  EC2 et Cie, 1995.  International Conference on Artificial Neural Networks ICANN’95.
  * [5] I. Bleyer and R. Ramlau.  A double regularization approach for inverse problems with noisy data and inexact operator.  Inverse Problems, 29:025004, 2013.
  * [6] M. Burger and S. Osher.  Convergence rates of convex variational regularization.  Inverse Problems, 20(5):1411–1421, 2004.
  * [7] E. Celledoni, M. J. Ehrhardt, C. Etmann, B. Owren, C.-B. Schönlieb, and F. Sherry.  Equivariant neural networks for inverse problems.  Inverse Problems, 37(8):085006, 2021.
  * [8] A. Chambolle and T. Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40:120–145, 2011.
  * [9] Z. Cheng, M. Gadelha, S. Maji, and D. Sheldon.  A Bayesian Perspective on the Deep Image Prior.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [10] P. L. Combettes and J.-C. Pesquet.  Deep neural network structures solving variational inequalities.  Set-Valued and Variational Analysis, 28:491–518, 2020.
  * [11] P. L. Combettes and V. R. Wajs.  Signal recovery by proximal forward-backward splitting.  Multiscale Modeling & Simulation, 4(4):1168–1200, 2005.
  * [12] I. Daubechies, M. Defrise, and C. De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics, 57(11):1413–1457, 2004.
  * [13] S. Dittmer, T. Kluth, P. Maass, and D. Otero Baguer.  Regularization by Architecture: A Deep Prior Approach for Inverse Problems.  Journal of Mathematical Imaging and Vision, 62:456–470, 2020.
  * [14] H. W. Engl, K. Kunisch, and A. Neubauer.  Convergence rates for Tikhonov regularisation of non-linear ill-posed problems.  Inverse Problems, 5(4):523–540, 1989.
  * [15] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [16] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
  * [17] K. Gregor and Y. LeCun.  Learning fast approximations of sparse coding.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, page 399–406, Madison, WI, USA, 2010. Omnipress.
  * [18] A. Habring and M. Holler.  A generative variational model for inverse problems in imaging.  SIAM Journal on Mathematics of Data Science, 4(1):306–335, 2022\.
  * [19] R. Heckel and M. Soltanolkotabi.  Compressive sensing with un-trained neural networks: Gradient descent finds a smooth approximation.  In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 4149–4158. PMLR, 2020.
  * [20] R. Heckel and M. Soltanolkotabi.  Denoising and regularization via exploiting the structural bias of convolutional generators.  In International Conference on Learning Representations, 2020.
  * [21] B. Hofmann, B. Kaltenbacher, C. Pöschl, and O. Scherzer.  A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators.  Inverse Problems, 23(3):987–1010, 2007.
  * [22] V. Jain and S. Seung.  Natural image denoising with convolutional networks.  In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2009.
  * [23] B. Kaltenbacher and A. Klassen.  On convergence and convergence rates for Ivanov and Morozov regularization and application to some parameter identification problems in elliptic PDEs.  Inverse Problems, 34(5):055008, 2018.
  * [24] T. Kluth, C. Bathke, M. Jiang, and P. Maass.  Joint super-resolution image reconstruction and parameter identification in imaging operator: analysis of bilinear operator equations, numerical solution, and application to magnetic particle imaging.  Inverse Problems, 36(12):124006, 2020.
  * [25] V. Lempitsky, A. Vedaldi, and D. Ulyanov.  Deep Image Prior.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 9446–9454, 2018.
  * [26] H. Li, J. Schwab, S. Antholzer, and M. Haltmeier.  NETT: solving inverse problems with deep neural networks.  Inverse Problems, 36(6):065005, jun 2020.
  * [27] C. G. Looney.  Convergence of minimizing sequences.  Journal of Mathematical Analysis and Applications, 61(3):835–840, 1977.
  * [28] V. Papyan, Y. Romano, J. Sulam, and M. Elad.  Theoretical foundations of deep learning via sparse representations: A multilayer sparse model and its connection to convolutional neural networks.  IEEE Signal Processing Magazine, 35(4):72–89, 2018.
  * [29] Y. Romano, M. Elad, and P. Milanfar.  The little engine that could: Regularization by denoising (red).  SIAM Journal on Imaging Sciences, 10(4):1804–1844, 2017.
  * [30] L. I. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, 60(1):259–268, 1992.
  * [31] T. I. Seidman and C. R. Vogel.  Well posedness and convergence of some regularisation methods for non-linear ill posed problems.  Inverse Problems, 5(2):227–238, 1989.
  * [32] Z. Shi, P. Mettes, S. Maji, and C. G. M. Snoek.  On Measuring and Controlling the Spectral Bias of the Deep Image Prior.  International Journal of Computer Vision, 2022.
  * [33] J. Sjoberg and L. Ljung.  Overtraining, regularization, and searching for minimum with application to neural networks.  International Journal of Control, 62, 1994.
  * [34] Y. Sun, H. Zhao, and J. Scarlett.  On architecture selection for linear inverse problems with untrained neural networks.  Entropy, 23:1481, 2021.
  * [35] Y. Tai, J. Yang, and X. Liu.  Image super-resolution via deep recursive residual network.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2790–2798, 2017.
  * [36] V. Vasin.  Relationship of several variational methods for the approximate solution of ill-posed problems.  Mathematical notes of the Academy of Sciences of the USSR, 7:161–165, 1970.
  * [37] C. R. Vogel.  A constrained least squares regularization method for nonlinear ill-posed problems.  SIAM Journal on Control and Optimization, 28(1):34–49, 1990.
  * [38] H. Zou and T. Hastie.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 67(2):301–320, 2005.