Untrained neural network embedded Fourier phase retrieval from few
  measurements
Fourier phase retrieval (FPR) is a challenging task widely used in various
applications. It involves recovering an unknown signal from its Fourier
phaseless measurements. FPR with few measurements is important for reducing
time and hardware costs, but it suffers from serious ill-posedness. Recently,
untrained neural networks have offered new approaches by introducing learned
priors to alleviate the ill-posedness without requiring any external data.
However, they may not be ideal for reconstructing fine details in images and
can be computationally expensive. This paper proposes an untrained neural
network (NN) embedded algorithm based on the alternating direction method of
multipliers (ADMM) framework to solve FPR with few measurements. Specifically,
we use a generative network to represent the image to be recovered, which
confines the image to the space defined by the network structure. To improve
the ability to represent high-frequency information, total variation (TV)
regularization is imposed to facilitate the recovery of local structures in the
image. Furthermore, to reduce the computational cost mainly caused by the
parameter updates of the untrained NN, we develop an accelerated algorithm that
adaptively trades off between explicit and implicit regularization.
Experimental results indicate that the proposed algorithm outperforms existing
untrained NN-based algorithms with fewer computational resources and even
performs competitively against trained NN-based algorithms.
Liyuan Ma, Hongxia Wang, Ningyi Leng, Ziyang Yuan
# Untrained neural network embedded Fourier phase retrieval from few measurements

###### Abstract

Fourier phase retrieval (FPR) is a challenging task widely used in various applications. It involves recovering an unknown signal from its Fourier phaseless measurements. FPR with few measurements is important for reducing time and hardware costs, but it suffers from serious ill-posedness. Recently, untrained neural networks have offered new approaches by introducing learned priors to alleviate the ill-posedness without requiring any external data. However, they may not be ideal for reconstructing fine details in images and can be computationally expensive. This paper proposes an untrained neural network (NN) embedded algorithm based on the alternating direction method of multipliers (ADMM) framework to solve FPR with few measurements. Specifically, we use a generative network to represent the image to be recovered, which confines the image to the space defined by the network structure. To improve the ability to represent high-frequency information, total variation (TV) regularization is imposed to facilitate the recovery of local structures in the image. Furthermore, to reduce the computational cost mainly caused by the parameter updates of the untrained NN, we develop an accelerated algorithm that adaptively trades off between explicit and implicit regularization. Experimental results indicate that the proposed algorithm outperforms existing untrained NN-based algorithms with fewer computational resources and even performs competitively against trained NN-based algorithms.

## 1 Introduction

In optics, most detectors can only record the magnitude or intensity of the signal, while losing the phase information [1]. Fourier phase retrieval (FPR) is
an important inverse problem that seeks to recover an unknown signal \(\mathbf{x}^{*}\in\mathbb{R}^{n}\) from its Fourier magnitude \(\mathbf{b}\in\mathbb{R}^{m}\). The problem can be described as

\[\text{Find }\mathbf{x}\in\mathbb{R}^{n}\text{ s.t. }\left|\mathcal{F}\mathbf{x} \right|+\boldsymbol{\delta}=\mathbf{b}, \tag{1}\]

where \(\mathcal{F}\) represents the Fourier transform operator, \(\left|\cdot\right|\) denotes the element-wise absolute value operator, and \(\boldsymbol{\delta}\) is the additive noise. Note that a 2D image can also be expressed as \(\mathbf{x}^{*}\in\mathbb{R}^{n}\) with \(n_{1}\times n_{2}\) pixels, and the corresponding Fourier magnitude as \(\mathbf{b}\in\mathbb{R}^{m}\) with \(m_{1}\times m_{2}\) elements by a lexicographical order, where \(n_{1}\cdot n_{2}=n\) and \(m_{1}\cdot m_{2}=m\). FPR arises in many applications, such as X-ray crystallography, ptychography, and diffraction imaging [2, 3, 4, 5].

FPR is an ill-posed inverse problem due to the non-uniqueness of its solution. Even though \(m\geq 2n-1\), there can be up to \(2^{n-2}\) solutions for Eq. (1), along with global phase shift, conjugate inversion, and spatial shift [6]. Most existing literature [7, 8, 9, 10] conducts simulation experiments under the oversampled condition of \(m\geq 2n-1\). FPR with few measurements is seriously ill-posed, as the uniqueness of the solution lacks theoretical analysis under the condition of \(m<2n-1\). Classical algorithms such as Gerchberg-Saxton (GS) [11] and hybrid input and output (HIO) [12] often perform poorly under this condition. However, there are still many scenarios [13, 14] where FPR with few measurements is common, especially when the detector's resolution is limited.

To overcome the ill-posedness of FPR, there are various regularization methods that introduce priors of the signal to FPR models. One of the common models is built as

\[\hat{\mathbf{x}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathbb{R}^{n}}\,\frac {1}{2m}\|\mathbf{b}-\left|\mathcal{F}\mathbf{x}\right|\|_{2}^{2}+\alpha R( \mathbf{x}), \tag{2}\]

where the first and second terms are data fidelity and regularization terms respectively, and the regularization parameter \(\alpha>0\) pursues the tradeoff between these two terms. Explicit regularizers, such as \(\delta_{\Omega}(\mathbf{x})\), \(\|\mathbf{x}\|_{1}\), and \(\|\mathbf{x}\|_{TV}\)[15, 16, 17], often capture certain properties of signals, such as support domain constraints, sparsity, and edge preservation. However, they are usually insufficient for representing complex information in images.

Recently, implicit regularizers based on learning have been proposed to represent more complicated priors of images. According to whether labeled data is required for training, this approach can be divided into supervised and unsupervised learning. [18] uses supervised learning with pairs of data \(\{\mathbf{b}_{i},\mathbf{x}_{i}\}\) to train an end-to-end neural network (NN) that represents an implicit prior. The Plug-and-Play (PnP) framework [19] leverages a pre-trained denoiser to implicitly represent a denoising prior. prDeep [8] adapts the Regularization by Denoising (RED) framework \(R(\mathbf{x})=\frac{1}{2}\mathbf{x}^{T}\left(\mathbf{x}-D(\mathbf{x})\right)\) to solve phase retrieval (PR) problems. RED penalizes the residual difference between the image with its denoised self and the correlations between the image with the residual, where \(D(\mathbf{x})\) denotes the pre-trained DnCNN [20] denoiser. In addition, trained generative networks [21, 22, 23] learn the distribution of images from labeled datasets to construct a generative prior. However, all of these supervised learning-based
algorithms require acquiring a large amount of data with ground truth images. This can be expensive or even impossible in many fields. Moreover, limited generalization can easily occur when the distribution of the training and test data is inconsistent.

Unsupervised learning offers new approaches to impose implicit regularization without requiring labeled data. Among them, untrained NNs [24] propose a new paradigm for training randomly initialized NNs using only a single observation, without any external data. Untrained NNs, such as Deep Image Prior (DIP) [25] and its variant Deep Decoder (DD) [26], demonstrate that their network architecture can capture statistical priors of an image from the observation itself. Net-PGD [27] imposes an untrained generative prior by constraining \(\mathbf{x}\) to the generation space represented by a DD-based network, and solves Gaussian PR using the projected gradient descent (PGD) framework. It is beneficial for Gaussian PR with few measurements compared to hand-crafted priors. However, it is not ideal for solving more challenging FPR. DeepMMSE [10] is an unsupervised learning-based algorithm which utilizes an untrained generative network with dropout to approximate the minimum mean squared error (MMSE) estimator of the image. Although it performs well in terms of FPR, its computational cost is quite high due to numerous iterations and a large number of network parameters. When \(\mathbf{x}\) is a \(128\times 128\) grayscale image, \(\mathbf{b}\) has a size of \(256\times 256\), and the network structures are set as specified in section 3, the number of network parameters in DeepMMSE is 708402, which is almost seven times that of Net-PGD, which has 108160 parameters. This makes it impractical to use when computational resources are limited. Therefore, it is necessary to explore unsupervised learning-based FPR algorithms that can achieve both fast and high-quality reconstruction.

To achieve fast and high-quality reconstruction of FPR using unsupervised learning, we adopt a DD-based network to impose a generative prior for FPR due to its latent ability to reconstruct the image from few measurements. However, unlike Net-PGD, we use the superior alternating direction method of multipliers

Figure 1: Comparison of different algorithms based on untrained NNs conducting FPR on _Cameraman_. The networks do not require any external data for training. We present heat maps of reconstructed errors, which show the difference between the reconstructed image and the ground truth. The color scheme on the right side of the heat map represents the range of error values. The computational times for DD, Net-PGD, DeepMMSE and our algorithm are 58.48s, 68.23s, 739.99s, and 36.02s, respectively.
(ADMM) [7, 28] framework. As shown in Fig. 1, the DD algorithm, which employs a DD-based network to impose a generative prior for FPR within the ADMM framework, outperforms Net-PGD in the PGD framework. In addition, multiple studies [29, 30] indicate that untrained NNs tend to fit smooth images and struggle to recover high-frequency structures. This is further supported by Fig. 1, which shows that three algorithms based on untrained NNs have large reconstruction errors on the high-frequency structures in the image. To address this problem, we combine the DD-based network with total variation (TV) regularization, which is beneficial for preserving the local structures of images [31, 32]. Based on the ADMM framework, we propose a vanilla TV-regularized FPR algorithm with DD, which is called Vanilla-TRAD.

In addition, we develop an accelerated algorithm for achieving fast reconstruction of FPR. In Vanilla-TRAD, the generative network is trained to find a solution within its range for each iteration of FPR. This results in a high computational cost. We have observed that the later iterations of the generative network are unnecessary since they yield solutions with little improvement in quality. Therefore, we introduce an acceleration strategy inspired by the hybrid steepest descent (HSD) method [33, 34]. This strategy dynamically adjusts whether to use DD in the iterations, which improves the reconstruction speed. The proposed algorithm is called Accelerated-TRAD. Fig. 1 illustrates the benefits of our proposed Accelerated-TRAD in terms of reconstruction quality and computational time. The acceleration strategy offers a new approach to adaptively trading off between explicit and implicit regularization.

The performance of Vanilla-TRAD and Accelerated-TRAD is extensively evaluated under various settings. The experimental results of FPR with different measurement lengths and noise levels show that Accelerated-TRAD outperforms existing untrained NN-based algorithms with less computational cost and performs competitively against trained NN-based algorithms. Analysis of the parameters shows that Accelerated-TRAD is not sensitive to parameter selection in the acceleration strategy.

The rest of this paper is organized as follows. In section 2, Vanilla-TRAD and Accelerated-TRAD are introduced. In Section 3, we compare the proposed algorithms to state-of-the-art algorithms under various settings. We also analyze the selection of regularizers and parameters for the proposed algorithms. Section 4 is the conclusion. The code is available at [https://github.com/Liyuan-2000/TRAD.git](https://github.com/Liyuan-2000/TRAD.git).
The proposed methods

### Vanilla TV-regularized FPR algorithm under untrained generative prior

We consider a regularized FPR problem under untrained generative prior. According to Eq. (2), the problem can be described as

\[\begin{split}&\min_{\mathbf{x},\boldsymbol{\theta}}\;\frac{1}{2m} \|\mathbf{b}-\left|\mathcal{F}\mathbf{x}\right|\|_{2}^{2}+\alpha R(\mathbf{x} ),\\ &\mathrm{s.t.}\quad\mathbf{G}(\boldsymbol{\theta})-\mathbf{x}=0, \end{split} \tag{3}\]

where \(R(\mathbf{x})\) is an explicit regularizer that enforces desired properties on \(\mathbf{x}\), while \(\mathbf{G}(\boldsymbol{\theta})\) is an untrained generative network with parameters \(\boldsymbol{\theta}\) that restricts \(\mathbf{x}\) in \(Range(\mathbf{G})\).

An untrained generative network can implicitly capture priors of an image \(\mathbf{x}\) through its structure without external data [25, 26, 35]. In this paper, a fixed low-dimensional latent code \(\mathbf{Z}_{0}\in\mathbb{R}^{d_{0}\times c_{0}}\) is chosen as the input of network. The output is \(\mathbf{G}(\boldsymbol{\theta})\in\mathbb{R}^{d_{J}\times c_{out}}\), where \(c_{out}=1\) for a grayscale image, \(c_{out}=3\) for an RGB image, and \(d_{J}\cdot c_{out}=n\). The structure of the network is illustrated by Fig. 2. It is composed of \(J\) layers, each of which consists of \(1\times 1\) convolutions, a ReLU activation function \(\mathrm{relu}(\cdot)\), a channel normalization operator \(\mathrm{cn}(\cdot)\) and an upsampling operator \(\mathbf{U}_{j}\in\mathbb{R}^{d_{j+1}\times d_{j}}\). That is to say,

\[\mathbf{Z}_{j+1}=\mathbf{U}_{j}\mathrm{cn}\left(\mathrm{relu}\left(\mathbf{Z} _{j}\mathbf{W}_{j}\right)\right),j=0,1,\cdots,J-1, \tag{4}\]

where \(\mathbf{W}_{j}\in\mathbb{R}^{c_{j}\times c_{j+1}}\) denotes the weight matrix corresponding to the \(1\times 1\) convolutions. The output layer of the network includes the \(1\times 1\) convolutions and a sigmoid activation function \(\mathrm{sigmoid}(\cdot)\). It is expressed as:

\[\mathbf{G}(\boldsymbol{\theta})=\mathrm{sigmoid}\left(\mathbf{Z}_{J}\mathbf{ W}_{J}\right), \tag{5}\]

where \(\mathbf{W}_{J}\in\mathbb{R}^{c_{J}\times c_{out}}\). All parameters in the network can be vectored as \(\boldsymbol{\theta}=\mathrm{vec}\left(\mathbf{W}_{0},\mathbf{W}_{1},\cdots, \mathbf{W}_{J}\right)\) and \(\mathrm{vec}(\cdot)\) is the vectorization of a matrix. For brevity, we utilize \(\{c_{0},c_{1},\cdots,c_{J}\}\) to represent the structure of a \(J\)-layer network. To reduce the cost of computation and storage, we adopt a network [27] similar to DD [26], where the sizes of the channels satisfy \(d_{0}<d_{1}<\cdots<d_{J}\).

As mentioned before, while untrained generative networks are useful for reducing the ill-posedness of PR with few measurements, they struggle to recover high-frequency structures since they tend to fit smooth images [29]. We overcome this limitation by choosing a proper regularizer \(R(\mathbf{x})\). For instance, \(R(\mathbf{x})=\|\mathbf{x}\|_{TV}\) encourages the preservation of local structures [31]. The TV norm \(\|\mathbf{x}\|_{TV}\) can be either anisotropic or isotropic, i.e., \(\|\nabla\mathbf{x}\|_{1}\) or \(\|\nabla\mathbf{x}\|_{2}\). Other regularization terms that protect high-frequency characteristics can also be considered for selection.

Eq. (3) is a non-convex optimization problem with multiple local minima. This makes it challenging for general gradient-based algorithms to find the global
optimal solution. To solve Eq. (3), we use the ADMM [7, 36, 37] framework, which mitigates the difficulty by dividing the problem into several subproblems that can be solved more easily. Note that the data fidelity term in Eq. (3) is not Lipschitz differentiable, making it difficult to design efficient algorithms that use gradients. To address this, we replace it with its smoothed version [38]:

\[f(\mathbf{u})=\frac{1}{2m}\|\sqrt{\mathbf{b}^{2}+\varepsilon\mathbf{1}}-\sqrt{ |\mathcal{F}\mathbf{u}|^{2}+\varepsilon\mathbf{1}}\|_{2}^{2}, \tag{6}\]

where \(\varepsilon>0\) represents the penalization parameter, and \(\mathbf{1}\in\mathbb{R}^{m}\) denotes a vector whose elements are all ones. The gradient of \(f(\mathbf{u})\) is

\[\nabla f(\mathbf{u})=\mathbf{u}-\mathcal{F}^{-1}\left(\frac{\sqrt{\mathbf{b}^ {2}+\varepsilon\mathbf{1}}}{\sqrt{|\mathcal{F}\mathbf{u}|^{2}+\varepsilon \mathbf{1}}}\odot\mathcal{F}\mathbf{u}\right), \tag{7}\]

where \(\odot\) denotes the Hadamard product and \(\mathcal{F}^{-1}\) denotes the inverse Fourier transform operator.

Now we reformulate Eq. (3) into

\[\begin{split}&\min_{\mathbf{x},\boldsymbol{\theta}}\ f\left(\mathbf{G}(\boldsymbol{\theta}) \right)+\alpha\|\mathbf{x}\|_{TV},\\ &\text{s.t.}\quad\mathbf{G}(\boldsymbol{\theta})-\mathcal{P} \mathbf{x}=0,\end{split} \tag{8}\]

where \(\mathcal{P}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) is a zero-padding operator, which appends \(m-n\) zeros after the last element of \(\mathbf{x}\). Introducing \(\mathcal{P}\) aims to execute the Fourier transform using the Fast Fourier transform (FFT) algorithm.

By means of the augmented Lagrangian multiplier method, we transform Eq. (8) into an unconstrained optimization problem:

\[\underset{\boldsymbol{\eta}}{\text{maxim}}\min_{\boldsymbol{\theta},\mathbf{x }}\ L_{\rho}(\boldsymbol{\theta},\mathbf{x},\boldsymbol{\eta})=f(\mathbf{G}( \boldsymbol{\theta}))+\alpha\|\mathbf{x}\|_{TV}+\langle\boldsymbol{\eta}, \mathbf{G}(\boldsymbol{\theta})-\mathcal{P}\mathbf{x}\rangle+\frac{\rho}{2}\| \mathbf{G}(\boldsymbol{\theta})-\mathcal{P}\mathbf{x}\|_{2}^{2}, \tag{9}\]

where \(\boldsymbol{\eta}\in\mathbb{R}^{m}\) is the Lagrangian multiplier associated with the equality constraint, and \(\rho>0\) is the coefficient of the quadratic penalty term. Following

Figure 2: Illustration of the untrained generative network used in this paper.
the ADMM framework, we can minimize \(L_{\rho}(\mathbf{\theta},\mathbf{x},\mathbf{\eta})\) with respect to \(\mathbf{\theta}\), \(\mathbf{x}\), and \(\mathbf{\eta}\) in an alternating way. Given initial variables \(\mathbf{x}_{0}\), \(\mathbf{\eta}_{0}\) and \(k=0\),

\[\mathbf{\theta}_{k+1} =\operatorname*{arg\,min}_{\mathbf{\theta}}\;\left\{f(\mathbf{G}(\bm {\theta}))+\frac{\rho}{2}\|\mathbf{G}(\mathbf{\theta})-\mathcal{P}\mathbf{x}_{k}+ \frac{\mathbf{\eta}_{k}}{\rho}\|_{2}^{2}\right\}, \tag{10a}\] \[\mathbf{x}_{k+1} =\operatorname*{arg\,min}_{\mathbf{x}}\;\left\{\alpha\|\mathbf{x }\|_{TV}+\frac{\rho}{2}\|\mathbf{G}(\mathbf{\theta}_{k+1})-\mathcal{P}\mathbf{x}+ \frac{\mathbf{\eta}_{k}}{\rho}\|_{2}^{2}\right\},\] (10b) \[\mathbf{\eta}_{k+1} =\mathbf{\eta}_{k}+\rho\left(\mathbf{G}(\mathbf{\theta}_{k+1})-\mathcal{ P}\mathbf{x}_{k+1}\right). \tag{10c}\]

The first subproblem (10a) is used to update the network parameters. We solve it in two steps. First, we calculate a vector \(\mathbf{u}\in\mathbb{R}^{m}\), which acts a substitute for \(\mathbf{G}(\mathbf{\theta})\) in Eq. (10a). Since \(f\) is differentiable, we make a first-order approximation of \(f\) at \(\mathcal{P}\mathbf{x}_{k}-\frac{\mathbf{\eta}_{k}}{\rho}\) and then calculate \(\mathbf{u}_{k+1}\) using the Karush-Kuhn-Tucker (KKT) condition:

\[\begin{split}\mathbf{u}_{k+1}&=\operatorname*{arg\, min}_{\mathbf{u}}\left\{f(\mathbf{u})+\frac{\rho}{2}\|\mathbf{u}-\mathcal{P} \mathbf{x}^{k}+\frac{\mathbf{\eta}_{k}}{\rho}\|^{2}\right\}\\ &=\operatorname*{arg\,min}_{\mathbf{u}}\;f\left(\mathcal{P} \mathbf{x}_{k}-\frac{\mathbf{\eta}_{k}}{\rho}\right)+\left\langle\nabla f\left( \mathcal{P}\mathbf{x}_{k}-\frac{\mathbf{\eta}_{k}}{\rho}\right),\mathbf{u}- \mathcal{P}\mathbf{x}_{k}+\frac{\mathbf{\eta}_{k}}{\rho}\right\rangle\\ &+\frac{\rho}{2}\|\mathbf{u}-\mathcal{P}\mathbf{x}_{k}+\frac{\bm {\eta}_{k}}{\rho}\|^{2}\\ &=\mathcal{P}\mathbf{x}_{k}-\frac{\mathbf{\eta}_{k}}{\rho}-\frac{1}{ \rho}\nabla f\left(\mathcal{P}\mathbf{x}_{k}-\frac{\mathbf{\eta}_{k}}{\rho}\right).\end{split} \tag{11}\]

The second step is to project \(\mathbf{u}_{k+1}\) into \(Range(\mathbf{G})\) by searching for a suitable \(\mathbf{\theta}_{k+1}\), which serves as an internal loop at the \(k+1\)-th iteration.

\[\mathbf{\theta}_{k+1}=\operatorname*{arg\,min}_{\mathbf{\theta}}\;\|\mathbf{G}(\mathbf{ \theta})-\mathbf{u}_{k+1}\|_{2}^{2}. \tag{12}\]

The detailed process for solving Eq. (12) will be explained in subsection 2.2. For convenience, we denote \(\mathbf{v}_{k+1}=\mathbf{G}\left(\mathbf{\theta}_{k+1}\right)\) and substitute \(\mathbf{v}_{k+1}\) for \(\mathbf{G}\left(\mathbf{\theta}_{k+1}\right)\) in Eq. (10b) and Eq. (10c).

Now consider the second subproblem (10b). We make a first-order approximation of \(\|\cdot\|_{TV}\) at \(\mathcal{P}^{-1}\left(\mathbf{v}_{k+1}+\frac{\mathbf{\eta}_{k}}{\rho}\right)\) again by the gradient \(\nabla\|\cdot\|_{TV}\) and then calculate \(\mathbf{x}_{k+1}\) using the KKT condition. The result is

\[\mathbf{x}_{k+1}=\mathcal{P}^{-1}\left(\mathbf{v}_{k+1}+\frac{\mathbf{\eta}_{k}}{ \rho}\right)-\frac{\alpha}{\rho}\Delta\left[\mathcal{P}^{-1}\left(\mathbf{v}_ {k+1}+\frac{\mathbf{\eta}_{k}}{\rho}\right)\right], \tag{13}\]

where \(\Delta\) denotes the Laplace operator due to \(\nabla\|\mathbf{x}\|_{TV}=\operatorname{div}\left(\frac{\nabla\mathbf{x}}{\| \mathbf{x}\|}\right)=\frac{1}{\|\mathbf{x}\|}\Delta\mathbf{x}\) and \(\operatorname{div}(\cdot)\) denotes the divergence.

As to the initialization, we just set \(\mathbf{x}_{0}=\mathcal{F}^{-1}(\mathbf{b})\) and \(\mathbf{\eta}_{0}=\mathbf{0}\). Thus, we have developed an algorithm to solve FPR, called the Vanilla TV-regularized FPR algorithm with DD (Vanilla-TRAD). The detailed steps are summarized in Algorithm 1.
### Accelerated TV-regularized FPR algorithm under untrained generative prior

Vanilla-TRAD effectively combines TV regularization and the untrained generative prior through the ADMM framework. However, its computational cost may be expensive because solving Eq. (12) is required to update the network parameters for each \(k\), even though parameter reduction has already been considered in the selection of network structure. In this subsection, we propose an acceleration algorithm to improve the efficiency of Vanilla-TRAD.

Before introducing the acceleration algorithm, we first discuss the specific implementation of Eq. (12) in Vanilla-TRAD. The loss function is \(g(\mathbf{\theta})=\|\mathbf{G}(\mathbf{\theta})-\mathbf{u}_{k+1}\|_{2}^{2}\). Given the initial network parameters \(\mathbf{\theta}_{0}\), learning rate \(\gamma_{0}\), and the number of internal loops \(l_{0}\), the parameters can be updated by minimizing \(g(\mathbf{\theta})\) using gradient descent,

\[\mathbf{\theta}_{k+1}^{(j+1)}=\mathbf{\theta}_{k+1}^{(j)}-\gamma_{k}\nabla g\left(\bm {\theta}_{k+1}^{(j)}\right),j=0,1,\cdots,l_{k}-1, \tag{14}\]

where \(\mathbf{\theta}_{k+1}^{(0)}=\mathbf{\theta}_{k}\) and \(\mathbf{\theta}_{k+1}=\mathbf{\theta}_{k+1}^{(l_{k})}\). Variants of gradient descent, such as SGD [39] and Adam [40], can also be used. The learning rate \(\gamma_{k}\) decreases as \(k\) increases, following the formula \(\gamma_{k}=\gamma_{0}\beta^{\lfloor k/\kappa_{1}\rfloor}\), where \(0<\beta<1\) denotes the decay factor and \(\kappa_{1}\) represents the epoch of decay. The number of internal loops, denoted by \(l_{k}\), is calculated as \(l_{k}=\text{round}\left(l_{0}\zeta^{\lfloor k/\kappa_{2}\rfloor}\right)\), where \(\zeta>1\) denotes the growth factor, \(\kappa_{2}\) represents the epoch of growth, and \(\text{round}(\cdot)\) is the integer-valued function. Note that the number of internal loops \(l_{k}\) gradually increases with \(k\), so the computational time for each iteration also increases with \(k\).

Referring to the HSD method [34], we modify the 4th line of Vanilla-TRAD to a weighted combination of \(\mathbf{u}_{k+1}\) and \(\mathbf{G}(\mathbf{\theta}_{k+1})\).

\[\mathbf{v}_{k+1}=\mu_{k}\mathbf{G}\left(\mathbf{\theta}_{k+1}\right)+(1-\mu_{k}) \mathbf{u}_{k+1}. \tag{15}\]

During the early iterations, we expect to play the role of the untrained generative prior to obtain a relatively stable intermediate solution. Once stability
is reached, the weight \(\mu_{k}\) can be reduced to 0, since it is unnecessary to project \(\mathbf{u}_{k+1}\) into \(Range(\mathbf{G})\) when it is close to the optimal solution. Thus, the weight \(\mu_{k}\) can be designed to be 1 in the early iterations and decrease to 0 later. We choose a specific formula given by

\[\mu_{k}=\exp\left\{-\left(\frac{\max\left\{0,k-\kappa_{3}\right\}}{\lambda} \right)^{2}\right\}, \tag{16}\]

where \(\kappa_{3}\) denotes the iteration step at which the weight \(\mu_{k}\) begins to decrease, and \(\lambda\) represents the rate of decrease. As shown in Fig. 3, a small \(\kappa_{3}\) indicates early weight decay, while a small \(\lambda\) indicates fast weight decay.

```
0:\(\mathbf{b},\mathbf{x}_{0}=\mathcal{F}^{-1}(\mathbf{b}),\boldsymbol{\eta}_{0}= \mathbf{0},\rho>0,\varepsilon>0\), \(K\), \(\boldsymbol{\theta}_{0}\), \(\gamma_{0}\), \(l_{0}\), \(\beta\), \(\zeta\), \(\lambda\), \(\kappa_{1}\), \(\kappa_{2}\), \(\kappa_{3}\)
0:\(\mathbf{x}_{K}\)
```
1:for\(k=0\) : \(K-1\)
2:\(\mathbf{u}_{k+1}=\mathcal{P}\mathbf{x}_{k}-\frac{\eta_{k}}{\rho}-\frac{1}{ \rho}\nabla f\left(\mathcal{P}\mathbf{x}_{k}-\frac{\eta_{k}}{\rho}\right)\)
3:\(\boldsymbol{\theta}_{k+1}^{(0)}=\boldsymbol{\theta}_{k}\)
4:\(\gamma_{k}=\gamma_{0}\beta^{\left\lfloor k/\kappa_{1}\right\rfloor}\)
5:\(l_{k}=\text{round}\left(l_{0}\zeta^{\left\lfloor k/\kappa_{2}\right\rfloor}\right)\)
6:for\(j=0\) : \(l_{k}-1\)
7:\(\boldsymbol{\theta}_{k+1}^{(j+1)}=\boldsymbol{\theta}_{k+1}^{(j)}-\gamma_{k} \nabla g\left(\boldsymbol{\theta}_{k+1}^{(j)}\right)\)
8:end
9:\(\boldsymbol{\theta}_{k+1}=\boldsymbol{\theta}_{k+1}^{(l_{k})}\)
10:\(\mu_{k}=\exp\left\{-\left(\frac{\max\{0,k-\kappa_{3}\}}{\lambda}\right)^{2}\right\}\)
11:\(\mathbf{v}_{k+1}=\mu_{k}\mathbf{G}\left(\boldsymbol{\theta}_{k+1}\right)+(1- \mu_{k})\mathbf{u}_{k+1}\)
12:\(\mathbf{x}_{k+1}=\mathcal{P}^{-1}\left(\mathbf{v}_{k+1}+\frac{\boldsymbol{\eta} _{k}}{\rho}\right)-\frac{\alpha}{\rho}\Delta\left[\mathcal{P}^{-1}\left( \mathbf{v}_{k+1}+\frac{\boldsymbol{\eta}_{k}}{\rho}\right)\right]\)
13:\(\boldsymbol{\eta}_{k+1}=\boldsymbol{\eta}_{k}+\rho\left(\mathbf{v}_{k+1}- \mathcal{P}\mathbf{x}_{k+1}\right)\)
14:end ```

**Algorithm 2** Accelerated-TRAD

Note that Accelerated-TRAD removes the network in later iterations. This benefits both the preservation of high-frequency structures in images and the
reduction of computational cost. The proposed accelerated algorithm, referred to as Accelerated-TRAD, is shown in Algorithm 2.

**Discussions.** 1. Accelerated-TRAD is an accelerated version of Vanilla-TRAD. For the first \(\kappa_{3}\) iterations, the format of Accelerated-TRAD is consistent with that of Vanilla-TRAD. In the case of \(\mu_{k}\equiv 1\), Accelerated-TRAD is equivalent to Vanilla-TRAD.

2. The main complexity of the proposed algorithms lies in the internal loops implemented by \(\mathbf{G}(\boldsymbol{\theta})\). In Accelerated-TRAD, the removal of the network in some iterations can effectively reduce computational complexity. Additionally, removing the network in later iterations is more effective than doing so earlier. This is because there are more internal loops in the later iterations, resulting in a higher time cost.

3. It seems that Accelerated-TRAD has added two parameters compared to Vanilla-TRAD. In fact, they can be simply chosen as \(\kappa_{3}=1000\) and \(\lambda=10\). It will be demonstrated through experiments in subsection 3.4 that the algorithm is not sensitive to parameter selection.

## 3 Experiment

Performance evaluation is conducted on FPR at different measurement lengths and noise levels. We compare the proposed algorithms with state-of-the-art algorithms on PR, especially learning-based methods, including prDeep [8], NetPGD [27] and DeepMMSE [10].

* prDeep [8]: A supervised learning algorithm based on the pre-trained DnCNN [20] is used to solve PR problems with different measurement models and noise levels. The algorithm integrates the DnCNN denoiser into the RED framework with the HIO initialization.
* Net-PGD [27]: An unsupervised learning-based algorithm combines the projected gradient descent (PGD) algorithm with a DD-based network to solve Gaussian PR. The algorithm implements gradient descent using external loops to produce the estimated solution, which is then projected into the network space using internal loops.
* DeepMMSE [10]: An unsupervised learning-based algorithm utilizes an untrained generative network with dropout to approximate the minimum mean squared error (MMSE) estimator of the image in PR. The algorithm also requires the HIO initialization, as described in [8].

The experimental data, shown in Fig. 4, consists of \(128\times 128\) grayscale images cropped or resized from standard \(256\times 256\) images, including 2 natural images, 2 remote sensing images, and 2 microscopic images. We evaluate the performance of the compared algorithms using the Peak Signal to Noise Ratio (PSNR), Structural Similarity (SSIM) of the reconstructed results, as well as the time cost of calculation. Due to the randomization factors involved in initializing
and updating network parameters, the algorithms may produce different results in different runs. To ensure a fair comparison, each algorithm is run five times. The average PSNR, SSIM and time cost of the five reconstructed results are then calculated.

The parameters in the proposed algorithms are as follows. We use isotropic TV norm in our experiments. The initial learning rate is \(\gamma_{0}=0.005\), and the decay factor and epoch are \(\beta=0.5\) and \(\kappa_{1}=500\). The initial number of internal loops is \(l_{0}=5\), and the growth factor and epoch are \(\zeta=1.2\) and \(\kappa_{2}=500\). The parameters \(\mathbf{\theta}_{0}\) in the network are initialized according to He initialization [41]. The elements of the latent code \(\mathbf{Z}_{0}\) are generated from the Gaussian distribution \(\mathcal{N}(0,0.01)\). Additionally, we set \(\rho=1\), \(\epsilon=0.001\) and \(\alpha=\frac{1}{884}\) by default. The structure of the network is designed to \(\{128,128,128,128\}\), which is a 3-layer network. The parameters in Accelerated-TRAD are set to \(\kappa_{3}=1000\) and \(\lambda=10\), which will be analyzed in subsection 3.4.

Following [10], the network input of DeepMMSE are generated from the Gaussian distribution \(\mathcal{N}(0,0.1)\). For Net-PGD, its network structure and internal loops are consistent with our algorithms. The learning rate of its external loops is initialized to 0.5, and the decay factor and epoch are set to 0.7 and 500, respectively. \(\sigma_{w}\) in prDeep is set to 0.1 when no noise is added to the Fourier magnitude. Both Net-PGD and the proposed algorithms are conducted 2000 iterations while DeepMMSE is conducted 50000 iterations. prDeep is executed for 200 iterations four times, once for each of the denoiser networks trained at standard deviations of 50, 40, 20, and 10. All NNs adopt the Adam optimizer. All algorithms are conducted on NVIDIA A-100 GPUs. Among them, prDeep uses MATLAB R2021a, while the other algorithms use the Pytorch framework with Python 3.9.

Figure 4: The data used in the experiment. First column: 2 natural images. Second column: 2 remote sensing images. Third column: 2 microscopic images.
### FPR at different sampling ratios

This subsection compares algorithms for FPR at six sampling ratios ranging from 1.5: 0.1: 2.0. The sampling ratio, denoted as \(r\), represents the ratio of the measurement length to image length in each dimension. For 2D images with equal width and height, \(r=\sqrt{m}/\sqrt{n}\). Table 1 shows the quantitative results for the compared algorithms at sampling ratios of 1.7 and 1.9. Additionally, a visual comparison between the algorithms at six sampling ratios can be seen in Fig. 5. For each algorithm, the images with the highest PSNR in five runs are shown.

Table 1 and Fig. 5 indicate that the reconstructed results of Accelerated-TRAD have almost the highest PSNR, SSIM and visual quality. Among all unsupervised learning-based algorithms, Accelerated-TRAD performs the best in terms of both the quality of the reconstructed results and computational time. It also performs competitively against supervised learning-based prDeep with only a minimal increase in computational time, no more than 10 seconds.

### FPR from noisy measurements

This subsection compares the robustness of algorithms for FPR. The measurement model with noise is:

\[\mathbf{b}=\left|\mathcal{F}\mathbf{x}\right|+\boldsymbol{\delta}, \tag{17}\]

where \(\boldsymbol{\delta}\sim\mathcal{N}(0,\sigma^{2})\) stands for Gaussian noise with a standard deviation of \(\sigma\).

When the sampling ratios are fixed, we compare algorithms at different noise levels. The noise level is represented by the Signal Noise Ratio (SNR), which is defined as

\[\text{SNR}=20\text{log}_{10}\frac{\text{Var}\left(\left|\mathcal{F}\mathbf{x} \right|\right)}{\sigma^{2}}, \tag{18}\]

where \(\text{Var}\left(\left|\mathcal{F}\mathbf{x}\right|\right)\) is the variance of \(\left|\mathcal{F}\mathbf{x}\right|\). A low SNR indicates a high level of noise.

\begin{table}
\begin{tabular}{l l l l l l l l} \hline \hline \multirow{2}{*}{\(r\)=1.7} & \multicolumn{6}{c}{PSNR(dB)/SSIM} & Time (s) \\  & \multicolumn{1}{c}{Cameraman} & Stream and Bridge & \multicolumn{1}{c}{Stadem} & \multicolumn{1}{c}{Tuson} & \multicolumn{1}{c}{Butterfly Nebula} & \multicolumn{1}{c}{Pillars of Creation} & Total \\ \hline prDeep [8] & 0.30/0.95 & 35.64/0.88 & 24.26/0.64 & 21.35/0.34 & **22.62/0.52** & 25.62/0.66 & **29.78** \\ Net-PGD [27] & 14.87/0.28 & 13.24/0.25 & 13.42/0.22 & 19.00/0.28 & 19.86/0.38 & 18.71/0.31 & 57.21 \\ DeepMISE [10] & 38.53/0.97 & 20.24/0.53 & 25.22/0.75 & 24.57/**0.39** & 19.34/0.42 & 20.86/0.48 & 78.486 \\ Vanilla-TRAD & 25.92/0.89 & 22.66/0.70 & 21.76/0.63 & 21.53/0.37 & 21.59/0.48 & 24.58/0.61 & 65.24 \\ Accerrelated-TRAD & **51.12/1.00** & **45.50/1.00** & **31.03/0.79** & **25.09/0.57** & 19.53/0.45 & **37.67/0.92** & **31.03** \\ \hline \multirow{2}{*}{\(r\)=1.9} & \multirow{2}{*}{Cameraman} & Stream and Bridge & \multicolumn{1}{c}{Stadem} & \multicolumn{1}{c}{Tuson} & \multicolumn{1}{c}{Butterfly Nebula} & \multicolumn{1}{c}{Pillars of Creation} & Total \\ \cline{2-7}
\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline  & \multicolumn{3}{c}{PSNR(dB)/SSIM (\(r\)=1.6)} & \multicolumn{3}{c}{PSNR(dB)/SSIM (\(r\)=1.8)} \\  & SNR=30 & SNR=20 & SNR=10 & SNR=30 & SNR=20 & SNR=10 \\ \hline prDeep [8] & 19.34/0.30 & 18.80/0.29 & 18.00/0.22 & 18.88/0.30 & 19.57/0.31 & 18.75/0.25 \\ Net-PGD [27] & 18.94/0.31 & 18.41/0.23 & 17.24/0.15 & 19.17/0.35 & 18.70/0.28 & 18.14/0.18 \\ DeepMMSE [10] & 19.29/0.29 & 19.16/0.24 & 17.65/0.16 & 20.46/0.33 & 19.28/0.25 & 18.63/0.19 \\ Vanilla-TRAD & **21.09/0.41** & 19.71/0.28 & **19.54/0.31** & **22.56/0.47** & **20.79/0.35** & 20.09/0.34 \\ Accelerated-TRAD & 20.33/0.33 & **20.68/0.36** & 19.05/0.28 & 21.54/0.36 & 20.45/0.32 & **20.49/0.35** \\ \hline \hline \end{tabular}
\end{table}
Table 2: The PSNR and SSIM of the reconstructed results by different algorithms on _Pillars of Creation_ at different noise levels. The best and second best results at each column are **boldfaced** and underlined respectively.

Figure 5: The reconstructed results of different algorithms at various sampling ratios, along with the corresponding PSNR(dB)/SSIM for each image annotated below.
Table 2 displays the PSNR and SSIM of the reconstructed results using different algorithms at SNRs of 30, 20, and 10 for sampling ratios of 1.6 and 1.8. It is evident that the proposed algorithms are more robust than others. The PSNR of reconstructed results in Vanilla-TRAD or Accelerated-TRAD is at least 1dB higher than that in the other three algorithms.

### Why combine TV with DD?

We construct two numerical experiments to explore the reasons for combining TV with DD. (1) To analyze the effectiveness of combining explicit and implicit regularization, we compare the proposed algorithms with related algorithms that have no regularization, TV regularization only, and DD regularization only. (2) To determine the necessity of combining TV and DD, we replace TV with explicit sparse or low-rank regularizer, and DD with SIREN [42] or DIP [25]. We then compare these alternatives to Accelerated-TRAD. All compared algorithms are executed under the ADMM framework. The average PSNR of six reconstructed images in the experimental data is calculated.

In the first experiment, the quantitative results of five algorithms at six sampling ratios are shown in Fig. 6. A visual comparison of the reconstructed _Cameraman_ at the sampling ratio of 1.7 can be seen in Fig. 7. Accelerated-TRAD performs the best and accurately recovers both high- and low-frequency information. While explicit TV regularization is effective for FPR at high sampling ratios, it may cause artifacts. In contrast, implicit DD regularization is good for FPR at low sampling ratios but leads to oversmoothed patterns. Vanilla-TRAD, which combines DD and TV regularization, can slightly alleviate this issue and improve the PSNR of reconstructed results. The accelerated technique used in Accelerated-TRAD not only reduces the time cost but also effectively plays the role of the two priors.

In the second experiment, the TV regularization is replaced with sparse regularization by \(l_{1}\) norm \(\|\mathbf{x}\|_{1}\) and low-rank regularization by nuclear norm \(\|\mathbf{x}\|_{*}\), respectively. Additionally, the DD is replaced with SIREN and DIP,

Figure 6: Quantitative results of our algorithms as well as related algorithms with no regularization, TV regularization only, and DD regularization only.
respectively. The corresponding network structures are described as follows. The learning rate and number of iterations are consistent with Accelerated-TRAD.

* SIREN [42]: An implicitly defined network is utilized to represent continuous signals and their derivatives. The network is composed of multilayer perceptrons (MLPs) with periodic sinusoidal activation. It consists of four hidden layers with 128 channels, takes image coordinates as input, and outputs pixel values.
* DIP [25]: An untrained generative network uses a U-Net-like "hourglass" architecture with skip connections. It includes upsampling and downsampling layers with channel numbers and kernel sizes of {16, 32, 64, 128, 128} and {3, 3, 3, 3, 3}, respectively. Upsampling and downsampling are implemented using nearest neighbor interpolation and convolution strides. Skip connections have channel numbers and kernel sizes of {4, 4, 4, 4, 4} and {1, 1, 1, 1}, respectively. The network input follows a uniform distribution of \(\mathcal{U}(0,0.1)\).

Fig. 8 shows the quantitative results of Accelerated-TRAD and its alternatives using different regularization combinations. It is apparent that Accelerated-TRAD, which combines TV and DD, outperforms the other regularization combinations. Additionally, the number of parameters in DD is only one-tenth that in DIP, and slightly higher than that in SIREN.

While experimental results have shown the reasons for combining TV and

Figure 7: Visual comparison of our algorithms as well as related algorithms with no regularization, TV regularization only, and DD regularization only at the sampling ratio of 1.7.
DD, the underlying mechanism behind this combination has yet to be theoretically explained. This is worth further research.

### The impact of parameters

This subsection examines the impact of the parameters \(\lambda\) and \(\kappa_{3}\) in Accelerated-TRAD. Firstly, we explore the influence of \(\kappa_{3}\), which represents the iteration step at which the weight \(\mu_{k}\) begins to decrease. We use Vanilla-TRAD to perform FPR on the experimental data, due to the consistency of Vanilla-TRAD and Accelerated-TRAD in the first \(\kappa_{3}\) iterations. Fig. 9 shows the PSNR of the reconstructed images during the iterations. We observe that the PSNR of almost all images stabilizes when \(k>500\). This indicates that Accelerated-TRAD is not sensitive to \(\kappa_{3}\), as long as \(\kappa_{3}>500\). Furthermore, \(\kappa_{3}\) generally does not need to change with different \(\mathbf{x}\).

Next, we consider the impact of \(\lambda\), which determines the decrease rate of \(\mu_{k}\). When \(\kappa_{3}=1000\), we apply Accelerated-TRAD to _Pillars of Creation_ at \(\lambda\) values

Figure 8: Quantitative results of Accelerated-TRAD and its alternatives using other explicit and implicit regularization combinations. The number of parameters in DD, SIREN, and DIP are 108160, 50049, and 1037193, respectively.

Figure 9: PSNR of reconstructed results during the iterations of Vanilla-TRAD at the sampling ratio of 2.0.
of 10, 100, and 200, respectively, and compare with Vanilla-TRAD. The PSNR of the reconstructed images during the iterations is shown in Fig. 10. It can be observed that Accelerated-TRAD is superior to Vanilla-TRAD in terms of both reconstructed quality and computational time. In Accelerated-TRAD, the PSNR of reconstructed results increases at different rates when \(\lambda\) takes different values, but eventually reaches similar results. A small \(\lambda\) tends to result in a low computational time, due to the rapid decrease of the weight \(\mu_{k}\).

## 4 Conclusion

This paper proposes an untrained NN-based algorithm called Vanilla-TRAD that effectively combines explicit TV regularization and implicit untrained generative prior with the ADMM framework, along with its accelerated version, Accelerated-TRAD. The untrained generative prior is beneficial for FPR with few measurements, while TV regularization is beneficial for recovering high-frequency information in the image. Numerical experiments confirm the effectiveness and necessity of combining the two priors. The acceleration technique in Accelerated-TRAD reduces the computational cost and enhances the role of the two priors in improving the reconstructed quality. Extensive experiments demonstrate that Accelerated-TRAD achieves fast and high-quality reconstruction of FPR under various settings. This paper presents a new approach to FPR under limited measurements and computational resources. In the future, we plan to explore a general paradigm that combines the prior based on image gradient with the untrained generative prior.

Figure 10: Comparison of Vanilla-TRAD and Accelerated-TRAD with different \(\lambda\). The algorithms perform FPR on _Pillars of Creation_ at the sampling ratio of 2.0. The computational times for Vanilla-TRAD and Accelerated-TRAD with \(\lambda=10\), \(\lambda=100\) and \(\lambda=200\) are 70.77s, 35.37s, 41.22s and 47.67s, respectively.