{'fraction_non_alphanumeric': 0.042300466638385555, 'fraction_numerical': 0.010423610690261197, 'mean_word_length': 4.340453074433657, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bilevel learning of l1-regularizers with closed-form gradients(BLORC)

We present a method for supervised learning of sparsity-promoting regularizers, a key ingredient in many modern signal reconstruction problems. The parameters of the regularizer are learned to minimize the mean squared error of reconstruction on a training set of ground truth signal and measurement pairs. Training involves solving a challenging bilevel optimization problem with a nonsmooth lower-level objective. We derive an expression for the gradient of the training loss using the implicit closed-form solution of the lower-level variational problem given by its dual problem, and provide an accompanying gradient descent algorithm (dubbed BLORC) to minimize the loss. Our experiments on simple natural images and for denoising 1D signals show that the proposed method can learn meaningful operators and the analytical gradients calculated are faster than standard automatic differentiation methods. While the approach we present is applied to denoising, we believe that it can be adapted to a wide-variety of inverse problems with linear measurement models, thus giving it applicability in a wide range of scenarios.


INTRODUCTION
Regularized image reconstruction has been used in many applications in imaging and signal processing. Typically, we are interested in solving the optimization problem of the form: where A ∈ R m×n is a linear forward model or measurement matrix, R : R n → R is a regularization functional, x ∈ R n denotes the underlying signal that we aim to recover, and y ∈ R m denotes its measurements. A prominent theme in designing the regularization functional, R, has been that of sparsity [1]: the idea that the reconstructed signal x should admit a sparse (having a small number of nonzero elements) Emails representation in some domain. While the 0 "norm" is a canonical measure of sparsity, but due to its non-convexity, its relaxed version, the 1 norm is often used in practice. A common and successful approach to promoting sparsity in signal reconstructions is to use R(x) = W x 1 , where W ∈ R k×n is an analysis dictionary or sparsifying transform. While there are several choices for W that work well in practice (e.g., wavelets, finite differences, discrete cosine transform (DCT)), several authors have also sought to learn W from data, an approach called dictionary (or transform) learning [2][3][4].
Learning the regularizer in a supervised fashion has been often cast as a bilevel problem in earlier works [5][6][7][8][9]. One of the methods to solve the bilevel problem was to replace the lower-level problem by its first order optimality conditions obtained through KKT conditions and then solving the bilevel problem using an implicit differentiation method. However, sparsity-promoting regularizers typically involve nonsmooth penalties such as the 1 norm. Some works [10][11][12][13] used relaxed versions of the 1 norm to perform implicit differentiation. However, these smooth approximations introduce errors in calculating the gradients close to the non-differentiable points. We propose a direct approach to solving the bilevel formulation by replacing the lower-level variational formulation by an implicit closed-form expression. While we demonstrate supervised learning of transforms for denoising (A = I), our approach can be extended to general inverse problems.
Automatic differentiation packages like PyTorch [14] and CVXPY's Diffcp [15][16][17] allow for the exact evaluation of the Jacobian of an arbitrarily complicated differentiable function, by partitioning the function into a sequence of simple operations, which are by themselves trivially differentiable. We compare our approach of deriving gradients by implicit differentiation of a closed-form expression with the auto-differentiation approaches in terms of accuracy and time of computation. To our knowledge, no work has to date solved the bilevel sparse signal reconstruction problem by exploiting implicit closedform expressions or without making relaxations. Extensive experimental results on 1D signals and images demonstrate that our method learns meaningful transforms taking into account both the dataset and the task at hand, and with less runtime than the autodifferentiation methods.

BILEVEL FORMULATION OF SUPERVISED LEARNING AND ITS ANALYSIS
The general bilevel sparsifying transform learning formulation for denoising is: The bilevel learning formulation (2) consists of two problems: the upper-level and the lower-level problem. The upperlevel problem minimizes the cost Q(W ) with respect to the parameters W , where Q(W ) compares the ground truth x t with the parameterized reconstruction output x * t (W ) obtained from the lower-level optimization problem. The lower-level problem is a variational reconstruction problem to obtain the parameterized reconstruction solution x * t (W ) from y t . In this supervised setting, the training samples (x t , y t ) are given and denote a clean image and its noisy version. Our objective is to learn the transform matrix W so that the lower level reconstruction using this sparsity operator is as close as possible to the ground truth. It is important to note two things in the above bilevel problem: first, no constraint is imposed on W during learning; and second, the algorithm learns the scaling of the regularization penalty, hence a separate scalar regularization strength is not required (i.e., β||W x|| 1 = ||βW x|| 1 for β ≥ 0). As a first step towards solving (2), we derive an implicit closed-form expression for the lower-level problem.

Closed-form solution obtained by duality
Consider the lower-level functional with x, y ∈ R n and W ∈ R k×n . It is strictly convex in x (the first term is strictly convex, while the second is convex), and therefore, has a unique global minimizer. Thus, we can write x * (W ) = arg min x J (x, W ), without the possibility of the minimizer not existing or being a nonsingleton. Note that although x * depends on both y and W , the y-dependence is not relevant for this derivation and we will not include it explicitly.
A key component in deriving the closed-form expression for the minimizer of (3) is that we need to know the sign pattern of W x * (W ). So, the closed-form expression is an implicit equation (4), where the reconstruction x * (W ) is dependent on sign(W x * ) (where sign(z) i is defined to be -1 when z i < 0; 0 when z i = 0; and 1 when z i > 0), which itself is a function of x * . So, to obtain the closed-form expression, we first solve for x * using any well-known iterative minimization algorithm like ADMM or PGD [18] and use it to compute the sign pattern. It is natural to question the approach of using an implicit closed-form equation if there exist well-known iterative minimization algorithms like ADMM and PGD to solve for x * . But we note that our objective in this step is not to find x * , but rather to obtain a closed-form expression that allows us to take gradients with respect to W .
Let the nonzero pattern s denote R sign(W x * ) and let W 0 and W ± contain the rows of W , whose indices are given by the set k =0 and k =0 , respectively. Then, the closed-form expression for x * (W ) = arg min x 1 2 x − y 2 2 + W x 1 is obtained from Lagrangian dual analysis [19] as where P N (W0) is the projector matrix onto the nullspace of W 0 and is given by , and R is a row selection matrix that helps retain only the nonzero elements of the sign vector sign(W x * ).

Gradient calculations
Once, we have a closed-form expression of the lower level problem, hour next step is to compute the derivative of x * (W ) with respect to W using matrix calculus [20]. Note that in order for the gradients to exist, the sign pattern vector sign(W x * ) has to remain constant in an open set containing W . Only then the closed-form expression for x * (W ) is valid in each region where sign(W x * ) is constant. As a brief example of this, consider the scalar denoising problem x * (w) = arg min x 1 2 (x − y) 2 + |wx| and s(w) = sign(wx * (w)). Assuming that y ≥ 0, one can show that x * (w) = y − |w| when y − |w| ≥ 0 and 0 otherwise. As a result, s((0, y)) = 1, s((−y, 0)) = −1, and s((−∞, −y] ∪ 0 ∪ [y, ∞)) = 0 is piecewise constant; a similar result holds when y ≤ 0. Thus Q(w) = (x * (w) − x t ) 2 is smooth except at w = 0, −y, y. So, there exist intervals in the domain of w, where Q(w) is differentiable with respect to w. This can be shown for higher dimensional cases as well.
Theorem 2 (Differential of the closed-form). Let c(W ) = sign(W x * (W )) denote a sign vector. Then if c(W ) is a constant vector in an open neighbourhood containing W , then the gradients of the closed-form expression in (4) wrt W exist and the form of the differential [20] is given by: where C = W † 0 ∂W 0 P N (W0) , and ∂x * , ∂W 0 , and ∂W ± are the differentials of x * ,W 0 , and W ± , respectively. So based on Theorem 2, we present the final expression of the gradient of the upper level cost U (W ) with respect to W .
be the upper-level cost function that is smooth with respect to the intermediate reconstruction x * (W ), then by the chain rule in matrix calculus, we derive an expression for the gradient of the cost U (W ) with respect to W 0 and W ± (denoted by ∇ W0 U and ∇ W± U ) as

BLORC ALGORITHM
We propose an iterative Bilevel Learning Of 1 Regularizers with Closed-form gradients (BLORC) algorithm for (2). We perform minibatch gradient descent to learn the transform matrix W from supervised training pairs (x t , y t ). In all our experiments, we start from W 0 = I n×n , the identity matrix. At the start of each epoch, the training pairs are randomly shuffled to remove any kind of dataset bias. Each sample in a batch is processed as follows. (a) given the measurement y t and current W , the lower-level reconstruction problem is solved iteratively using ADMM to obtain an estimate of x * t (W ). (b) the sign vector sign(W x * t (W )) is obtained after hard-thresholding W x * t (W ) with a small parameter γ. The sign vector is of paramount importance as it decides the rowsplit of W into W 0 and W ± . Hence, we run ADMM for a few thousands of iterations to ensure convergence of the sign pattern. (c) Then, using (7) and (8), we obtain the gradient of the upper level cost on a single training pair, i.e., ∇ W U which is the row-concatenation of ∇ W0 U and ∇ W± U . Averaging the gradients over the samples in the batch yields the minibatch gradient.
At the end of each batch, we update the matrix W based on the learning rate α and the mini-batch gradient. The updated W is used in the next batch in step (a) above. Our method of gradient calculation can also be extended when the training pairs are image patches instead of 1D signals. For image denoising experiments, image patches of size √ n × √ n with an overlap stride of r were extracted. Then the 2D patches were converted to 1D arrays as a pre-processing (first) step. In the end, the rows of the learned W were reshaped to look like convolutional filter patches. Except this first and last step, all the intermediate steps were the same for both 1D signals and image patches.

NUMERICAL EXPERIMENTS
To demonstrate how well BLORC learns the transform matrix W , we perform denoising experiments on synthetic 1D signals.
These synthetic 1D signals were chosen to be sparse with respect to a specific transform that also provides a baseline to compare our learned transforms with.

Denoising experiments
In our experiment, we generate M = 4000 training pairs (x t , y t ), where the x t 's are piece-wise constant signals of length n = 64 (with peak value normalized to 1) and y t 's are noisy versions with additive i.i.d. Gaussian noise with standard deviation σ = 0.1. Figure 1(a) shows a single pair of such (x t , y t ) . We perform minibatch gradient descent with batch-size B = 100 and run the algorithm for E = 750 epochs. The learning rate was chosen to be α = 10 −4 and the sign threshold was γ = 10 −3 . The learned transform is shown in Figure 1(b). We repeated the experiment with the same parameters but with x t 's chosen as smoothly varying signals of different harmonics that are sparse in the discrete cosine transform (DCT) domain as in Figure 1(c). The learned transform is row-rearranged such that it has maximum correlation with the 1D-DCT matrix and is shown in Figure 1 (d).
What is interesting in the learned transforms of Figures 1(b) and 1(d) is that they exactly capture the same intuition as the standard finite difference transform and the 1D-DCT transform, respectively, but in addition, they also have some novel features learned for the denoising task. Experimental results suggest that the learned transforms perform better than the standard transforms for denoising on a test-set. The average PSNR for 20 piecewise-constant test signals denoised using the learned W was 26.2 dB whereas that using the standard transform was 25.4 dB, with the PSNR of the noisy signals being 19.5 dB.
Extending the BLORC algorithm to images or more specifically to image patches, we learn reasonable transforms as well. We chose images (normalized) of size 256 × 256 with directional patterns (vertical strips and diagonal strips) and generated their noisy versions with i.i.d. Gaussian noise with σ = 0.1. We extracted image patches of size 8 × 8 with an overlap stride of 7. The reason for choosing a large stride was to ensure that distinct training pairs are obtained. As a pre-processing step, the image patches of size 8 × 8 were vectorized to size n = 64. All the parameters and the intermediate steps of the minibatch gradient descent are unchanged from the 1D signal experiment above. The learned convolutional filters in Figs. 2(b) and 2(d) look orthogonal to their image counterparts in Figs. 2(a) and 2(c), respectively, which is the expected output. In particular, these transforms are learned not to sparsify the training pairs but to minimize the gap between reconstructions and the ground truth.

Accuracy of gradient computations
While the BLORC algorithm uses an explicit form of the gradient, in automatic differentiation approaches, the task is divided into a sequence of differentiable operations as computational graphs on which backpropagation is performed through the chain rule. This division into a sequence of operations and calculating gradients for each node of the graph can take significant time, which can be bypassed if an explicit form of the gradient already exists that connects both the upper-level and the lower-level problem. This advantage in time can be crucial for larger datasets and batch-sizes. As a demonstration of this, we calculate and plot ∇ W U for a single training pair (x t , y t ) using three methods with the W = I initialization in Figure 3. In the first method, we use the direct expressions in (7) and (8) to get ∇ W U , which we denote as "By BLORC" in Figure 3. For the second method, we used the CVXPY package to run an iteration of optimization over (2) for the same training pair and obtained the gradient. Finally, as the third method, we used PyTorch to calculate the gradient of our closed-form expression. As the baseline for our comparisons, we calculated the numerical gradient of the upper level problem in (2) by noting the incremental change in the cost for incremental changes in each element of the matrix W . The errors for the three methods in Table 1 have been calculated with the ground truth value being set to the one from the numerical gradient method. The analytical form of the gradient in BLORC makes itt faster and more accurate compared to automatic differentiation approaches as is evident from Table 1. In Figure 3, the gradient ∇ W U calculated using all three methods encapsulates the feature of the single training pair (x t , y t ) but the one calculated using BLORC was more accurate overall as shown in Table 1.

CONCLUSION
In this paper, we solve the bilevel problem of learning sparsity regularizers by analytically calculating the gradients of an implicit closed-form expression. The learned regularizers are optimal for the task at hand (here, denoising) and outperform known regularizers on test data for the task. Our mathematical    analysis of bilevel reconstruction may lay the cornerstone for learning sparsifying (including deep) regularizers for inverse problems with general forward operators A, which is of major interest to the computational imaging community.