Efficient Decomposition-Based Algorithms for $\ell_{1}$-Regularized Inverse Problems with Column-Orthogonal and Kronecker Product Matrices

We are interested in solving discrete ill-posed inverse problems of the form ${\bf A}{\bf x}\approx{\bf b}$, where ${\bf A}\in\mathbb{R}^{m\times n}$, $m\geq n$, ${\bf b}\in\mathbb{R}^{m}$, and ${\bf x}\in\mathbb{R}^{n}$. It is assumed that the matrix ${\bf A}$ has full column rank but is ill-conditioned, and ${\bf b}$ is contaminated by additive Gaussian noise: ${\bf b}={\bf A}{\bf x}+{\bf e}$, where ${\bf e}$ is a Gaussian noise vector. Throughout this paper, we normalize the problem by dividing ${\bf A}$ and ${\bf b}$ by the standard deviation of ${\bf e}$. This corresponds to whitening the noise in the data. Due to the ill-posedness of this problem, the solution of the normal equations, ${\bf x}=({\bf A}^{\top}{\bf A})^{-1}{\bf A}^{\top}{\bf b}$, is poor, and therefore, we will apply regularization. With generalized Tikhonov regularization [48], ${\bf x}$ is selected to solve the minimization problem

where $\lambda$ is a regularization parameter and ${\bf L}\in\mathbb{R}^{p\times n}$ is a regularization matrix. ${\bf L}$ is often selected as the discretization of some derivative operator [28]. We also assume that

where $\mbox{Null}({\bf C})$ denotes the null space of the matrix ${\bf C}$.

Instead of considering the Tikhonov problem Eq. 1, we will consider a sparsity-preserving problem

where $\mu$ is a regularization parameter. We will solve Eq. 3 using two iterative methods, split Bregman (SB) [24] and Majorization-Minimization (MM) [31, 32]. Both methods share a minimization sub-problem of the same form [47].

Our focus is on Kronecker Product (KP) matrices ${\bf A}$ and ${\bf L}$, that is, matrices that can be written as ${\bf A}={\bf A}_{1}\otimes{\bf A}_{2}$ and ${\bf L}={\bf L}_{1}\otimes{\bf L}_{2},$ where $\otimes$ denotes the Kronecker product, which is defined for matrices ${\bf C}\in\mathbb{R}^{m\times n}$ and ${\bf D}\in\mathbb{R}^{p\times q}$ by

Thus, we are interested in solving problem Eq. 3 of the specific form

KP matrices ${\bf A}$ and ${\bf L}$ arise in many applications. In image deblurring problems, ${\bf A}$ is a KP matrix, provided that the blur can be separated in the horizontal and vertical directions [29]. Another case where ${\bf A}$ is a KP matrix is in 2D nuclear magnetic resonance (NMR) problems [6, 23], where the goal is to reconstruct a map of relaxation times from a probed material. ${\bf L}$ is a KP matrix with certain forms of regularization, such as when orthonormal framelets [44] or wavelets [12] are used in two dimensions. Both are tight frames, namely, they form redundant orthogonal bases for $\mathbb{R}^{n}$. This redundancy allows for robust representations where the loss of information is more tolerable [11]. When used as regularization in Eq. 3, framelet [8, 10, 11, 54] and wavelet-based [18] methods enforce sparsity on the solution within the tight frame.

In [8], tight frames are used for regularization in a general $\ell_{p}$-$\ell_{q}$ problem without the assumptions of having KP matrices. Other minimization problems with KP matrices ${\bf A}$ or ${\bf L}$ have been studied before. Least-squares problems of the form

have been solved efficiently using the singular value and QR decompositions [19, 20]. The special case where the forward operator has the form ${\bf A}={\bf I}\otimes{\bf A}_{2}$, which is equivalent to a matrix inverse problem, has also been analyzed [26]. For least-squares problems of the form

the generalized singular value decompositions (GSVDs) of $\{{\bf A}_{1},{\bf B}_{1}\}$ and $\{{\bf A}_{2},{\bf B}_{2}\}$ can be used to obtain a joint decomposition and solve the problem [3, 35, 52]. In particular, Givens rotations can be used to diagonalize the matrix containing the generalized singular values, which then diagonalizes the problem. A regularized problem of the form

occurs in adaptive optics [3, 35]. The structure of this problem can be utilized to rewrite it as a least-squares problem with a preconditioner. Another problem that utilizes KP matrices is the constrained least-squares problem

which is considered in [4]. In this case, both the forward operator and the constraint are KP matrices, and a null space method can be applied to the matrix equations. The solution is unique under certain circumstances, which can be analyzed using the GSVD [4].

In cases where ${\bf A}$ is not a KP matrix, there are methods to approximate it by a Kronecker product ${\bf A}\approx{\bf C}\otimes{\bf D}$ so that the benefits of the KP structure can still be utilized. The nearest Kronecker product can be found by minimizing

Van Loan and Pitsianis [53] developed a method for solving this problem, and in [42], Pitsianis extended it to solve the problem

These KP approximations have been applied to the forward matrix ${\bf A}$ in imaging problems [7, 22, 33, 34].

Main Contributions: We utilize joint decompositions that exploit the KP structure, computed at the first iteration only, to solve Eq. 5 with SB and MM where the regularization parameter is selected at each iteration as in [47]. We introduce a family of regularization operators that are column orthogonal. This family includes framelet and wavelet matrices. We show that with these matrices, we can obtain the GSVDs needed for the joint decompositions directly from the singular value decompositions (SVDs) of ${\bf A}_{1}$ and ${\bf A}_{2}$, allowing us to use framelets and wavelets without computing a GSVD. Other regularization matrices, such as finite-difference operators, would require computing GSVDs to implement these same decomposition-based algorithms. This makes using framelet and wavelet operators efficient for solving $\ell_{1}$-regularized problems with SB and MM. We also show that the solution of the generalized Tikhonov problem with a fixed regularization parameter will generate the same solution using framelets and wavelets.

This paper is organized as follows. In Section 2, we review two iterative methods for solving Eq. 3, split Bregman and Majorization-Minimization. In Section 3, the GSVD is presented and we discuss how it can be used to obtain a joint decomposition that can then be utilized within the iterative methods. Theoretical results are also presented for when ${\bf L}_{1}$ and ${\bf L}_{2}$ are column orthogonal, which allow this joint decomposition to be obtained from the SVD rather than the GSVD. In Section 4, we introduce a family of column-orthogonal regularization operators that include the framelet and wavelet regularization matrices. The iterative methods are then tested on image deblurring examples in Section 5. Although the results presented in Section 3 apply to a general ${\bf A}$, the examples are presented for cases with $m\geq n$. Conclusions are presented in Section 6.

In this section, we will consider solving Eq. 3 using two iterative methods: the split Bregman (SB) method and the Majorization-Minimization (MM) method. Both methods only require computing two GSVDs, which can be reused for every iteration.

Under the assumption that both ${\bf A}$ and ${\bf L}$ are KP matrices, we can solve the minimization problem Eq. 18 occurring for both the SB and MM methods by using a joint decomposition resulting from the individual GSVDs of $\{{\bf A}_{1},{\bf L}_{1}\}$ and $\{{\bf A}_{2},{\bf L}_{2}\}$. The GSVD is a joint decomposition of a matrix pair $\{{\bf A},{\bf L}\}$ introduced by Van Loan [51]. It has been formulated in different ways [30, 39], but for this paper, we use the definition of the GSVD given in [30].

The ratios $\gamma_{i}=\upsilon_{i}/\mu_{i}$ in the GSVD are known as the generalized singular values (GSVs) [27] and are sorted in decreasing order in Lemma 3.1. If $\mu_{i}=0$ and $\upsilon_{i}=1$, then the corresponding GSV $\gamma_{i}$ is infinite. The columns of ${\bf U}$ that correspond to $\upsilon_{i}=0$ span $\mbox{Null}({\bf A})$ while the columns of ${\bf V}$ that correspond to $\mu_{i}=0$ span $\mbox{Null}({\bf L})$.

From the GSVDs of the factors $\{{\bf A}_{1},{\bf L}_{1}\}$ and $\{{\bf A}_{2},{\bf L}_{2}\}$,

we can construct a joint decomposition of ${\bf A}$ and ${\bf L}$ using the properties of the Kronecker product [52]. This joint decomposition is given by

Let us assume, without loss of generality, that ${\bf x}$ is square, meaning that $\text{array}({\bf x})={\bf X}\in\mathbb{R}^{n\times n}$, where $\text{array}({\bf x})$ is an operator that reshapes the vector ${\bf x}\in\mathbb{R}^{n^{2}\times 1}$ into an $n\times n$ matrix. We will also assume that ${\bf A}_{1},{\bf A}_{2},{\bf L}_{1},{\bf L}_{2}$ have at least as many rows as columns. With assumption Eq. 2, $t=n$ in Lemma 3.1 and the solution to Eq. 18 is then given as

Since $t=n$, ${\bf Z}_{1}$, ${\bf Z}_{2}$, and ${\bf Z}={\bf Z}_{1}\otimes{\bf Z}_{2}$ are invertible, and we may set ${\bf Z}_{1}^{\top}={\bf Y}_{1}^{-1}$ and ${\bf Z}_{2}^{\top}={\bf Y}_{2}^{-1}$. Let $\tilde{{\bf U}}_{1}$ and $\tilde{{\bf U}}_{2}$ be the first $n$ columns of ${\bf U}_{1}$ and ${\bf U}_{2}$, respectively, and let $\tilde{{\bf V}}_{1}$ and $\tilde{{\bf V}}_{2}$ be the last $n$ columns of ${\bf V}_{1}$ and ${\bf V}_{2}$, respectively. Define also $\tilde{{\bf B}}=\tilde{{\bf U}}_{2}^{\top}{\bf B}\tilde{{\bf U}}_{1}$ and $\tilde{{\bf H}}^{(k)}=\tilde{{\bf V}}_{2}^{\top}{\bf H}^{(k)}\tilde{{\bf V}}_{1}$, where ${\bf B}=\text{array}({\bf b})$ and ${\bf H}^{(k)}=\text{array}({\bf h}^{(k)})$. Then, using Eq. 22 and

where $\text{vec}({\bf F})$ is the vectorization of the matrix ${\bf F}$, the solution in Eq. 23 can be written as

where $\odot$ denotes element-wise multiplication and $\tilde{\bm{\Phi}}$ and $\tilde{\bm{\Psi}}$ are given by

Here, $\upsilon_{i}$ and $\mu_{i}$ refer to the (potentially) non-zero entry in column $i$ of $\tilde{\bm{\Upsilon}}_{1}\otimes\tilde{\bm{\Upsilon}}_{2}$ or $\tilde{\bf M}_{1}\otimes\tilde{\bf M}_{2}$.

In both SB and MM, ${\bf L}{\bf x}$ must also be computed at each iteration. This computation can be written using Eq. 24 as

where $\hat{{\bf M}}=\text{array}(\left[\mu_{1},\dots,\mu_{n^{2}}\right]^{\top})$.

In the $\chi^{2}$ dof test, we also compute ${\bf L}^{\dagger}_{\bf A}{\bf h}^{(k)}$. Defining $\overline{{\bf m}}\in\mathbb{R}^{n^{2}}$ by

and $\overline{{\bf M}}=\text{array}(\overline{{\bf m}})$, we have that

Notice that Eq. 25, Eq. 26, and Eq. 27 allow each step in the iterative methods to be computed using only the individual GSVDs in Eq. 21.

This section presents a family of regularization operators ${\bf L}$, which we adopted for the solution of Eq. 3. Specifically, we will consider a two-level framelet analysis operator [13, 44] and the Daubechies D4 discrete wavelet transform [12].

Consider a matrix that corresponds to a one-dimensional (1D) transform of the form

and its two-dimensional (2D) transform consists of the form

where ${\bf T}_{ij}={\bf T}_{i}\otimes{\bf T}_{j}$ for $i,j=1,\dots,t.$ This then yields the 2D-transform-regularized problem

To solve Eq. 34 we observe that ${\bf L}_{T}^{1D}\otimes{\bf L}_{T}^{1D}={\bf P}_{T}{\bf L}_{T}^{2D}$, where ${\bf P}_{T}$ is a permutation matrix. Since

for any permutation matrix ${\bf P}$, we have that

As a result, solving