RANK-ONE AND SPARSE MATRIX DECOMPOSITION FOR DYNAMIC MRI

. We introduce a rank-one and sparse matrix decomposition model for dynamic magnetic resonance imaging (MRI). Since (cid:96) p -norm (0 < p < 1) is generally nonconvex, nonsmooth, non-Lipschitz, we propose reweighted (cid:96) 1 norm to surrogate (cid:96) p -norm. Based on this, we put forward a modiﬁed alterna- tive direction method. Numerical experiments are also given to illustrate the eﬃciency of our algorithm.


1.
Introduction. Dynamic magnetic resonance imaging (MRI) reconstructs a temporal series of images to resolve the motion or the variation of the imaged object. It is often used in cardiac, perfusion, functional, and gastro-intestinal imaging. However, achieving high spatio-temporal resolutions is challenging and inspiring in dynamic MRI due to the hardware limitations and the risk of peripheral nerve stimulation. This is because the data in MRI, samples in k-space of the spatial Fourier transform of the object, are acquired sequentially in time. A variety of MR techniques therefore aim to reduce the number of data required for accurate reconstruction.
Recently, R. Otazo et al. [12] applied low-rank and sparse model into MRI, and considered the following optimization problem min rank(L) + β T S 0 where · 0 is 0 -norm, which counts the number of nonzero entries. β ∈ R is a real positive weighting parameter to balance the weights of rank and sparsity, and d is the undersampled k − t data. L ∈ C m×n is a low-rank matrix, which is better at capturing the global structure. S ∈ C m×n is a sparse matrix, which can capture the local structure. Generally, we assume that the dynamic component S has a sparse representation in some known sparsifying transformation T , hence the idea of minimizing T S and not S itself. For acquisition with multiple receiver coils, E is given by the frame-by-frame multicoil encoding operator, which performs a multiplication by coil sensitivities followed by a Fourier transform according to the sampling pattern.

XIANCHAO XIU AND LINGCHEN KONG
It will be nice if one can solve (1), however, it is NP-hard due to the combinatorial nature of the 0 -norm. Instead one typically resorts to its relaxations that will be solvable in polynomial time. Fortunately, Chandrasekaran et al. [5] suggested to replace rank with nuclear norm and 0 norm with 1 norm, hence, (1) can be reformulated as follows where · * is the nuclear norm (defined as the sum of all singular values) of a matrix, and · 1 is the 1 norm in componentwise (defined as the sum of absolute values of all entries). In fact, 1 norm is a loose approximation of 0 norm and often leads to an overpenalized problem, consequently, some further improvements are required. Among such efforts, a very natural improvement is the suggestion of the use of p (0 < p < 1) norm, which is defined as x p = ( n i=1 |x i | p ) 1/p for any x ∈ C n . Even though the p norm is generally nonconvex, nonsmooth, non-Lipschitz, and NP-hard (see, e.g. [6,9]), it is shown in [8,7] that under some assumptions, the sequence generated by iterative reweighted 1 or 2 minimization algorithms converges to the sparest solution, which is also the global minimizer. Compared with the convex 1 norm regularization, it has significant advantages in many fields such as finance, econometrics and signal processing. Sequently, X. Li et al. [10] considered the following classical low-rank and sparse matrix decomposition model for surveillance video min rank(L) + S 0 and suggested to use rank-one matrix to surrogate low-rank matrix because the background component has relative small changes over a short period, that is, where u ∈ R m , 1 ∈ R n , and L = u1 T . Inspired by their work, in this paper, we focus our attention on the rank-one and sparse matrix decomposition model For classical low rank and sparse matrix decomposition model (2), Wright et al. [13] adopted the iterative thresholding technique to solve. However, the iterative thresholding scheme converges extremely slowly. Lin et al. [11] presented the novel augmented Lagrange multipliers, and established its convergence analysis. Building on the above work, Yuan and Yang [16] showed alternating direction method of multipliers (ADMM) to handle it. After that, ADMM began to be revisited frequently due to its success in the emerging application of structured convex optimization problems, like statistical learning [3], graphical model selection [2], matrix completion [11]. Xiu et al. [15] proposed modified iterative reweighted algorithm (MIRL1), which worked well in surveillance video. An interesting question occurs: if we can combine ADMM with MIRL1 to deal with dynamic MRI?
We will give an affirmative answer in this paper, and the rest of this paper is organized as follows. We will analyze (3) and derive modified alternative direction method in Section 2. Section 3 presents experimental results, including comparisons with other algorithm and computational performance. Finally, Section 4 concludes our paper.
2. Model and Algorithm. We begin this section with model analysis. It is easy to write the Lagrangian function associated with (3), with λ ∈ R being a tuning parameter. Now, we show how to apply modified alternating direction method (MADM) to solve (4). We first divide this optimization problem into three subproblems: where M 0 = E * d and M k can be derived from the third subproblem. Here, we call these subproblems in (5) as u-subproblem, S-subproblem, M-subproblem, respectively.
• S-subproblem seems much more troublesome. As we know, T S p p (0 < p < 1) is not Lipschitz continuous at some points, and this poses the difficulty to design algorithms for solving (4). Much of the recent research has concentrated on the reweighted 1 minimization problems [1,4,17], and shown that it performs well in recovering sparse signal. Thus, in this paper, for a fixed u k , the reweighted 1 version for S-subproblem is revised to where • denotes the Hadamard product, W is the reweighted matrix and 0 < W ij ≤ 1 for entries in W .
Following the similar idea of [15], we can derive the following modified iterative reweighted 1 algorithm (MIRL1) for the S-subproblem (6).
By exploiting the method from Literature [18] to update the weight, which is different from classic iterative reweighted 1 algorithm [4], we take W τ as
• M-subproblem ensures data consistency, which is Based on the above arguments, the framework of modified alternative direction method (MADM) for (4) is obtained as below.  [15] for more details.
3. Numerical Experiments. In this section, we conduct numerical experiments to compare the performance of the MADM and IST in [12]. All experiments are performed in MATLAB version 8.2 with Intel Core I5 2.60GHz CPU and 8GB of RAM. We test cardiac perfusion data, which are available at at the web site cai2r.net/resources/software/ls-reconstruction-matlab-code. Some example frames of the test data can be found in Figure 1.
In this paper, we choose parameter λ = 0.01 as suggested in [16], stopping criterion ε 0 = 5e − 4. The following Figures 1-4 show the computational results. Figure 1 shows the computational results of cardiac perfusion by IST, Figures 2-4 show the computational results by MADM with p = 0.1, 0.5, 0.9, respectively. In these figures, the first row is original frames, the second row is reconstructed background parts, and the third is reconstructed moving parts. From the numerical experiments, we see that: 1) MADM is easier to find the pathological parts than IST. See the third row of Figures 1-4; 2) The smaller p we choose, the more obvious interested parts we get. See the third row of . Conclusions. In this paper, we introduced a new rank-one and sparse matrix decomposition for dynamic MRI. Rather than handling the nonconvex p -norm directly, we applied the convex reweighted 1 -norm. Combining MIRLI and ADM, we developed MADM for (4), and illustrated the effectiveness of the proposed algorithm.