Low Rank Magnetic Resonance Fingerprinting

Magnetic Resonance Fingerprinting (MRF) is a relatively new approach that provides quantitative MRI measures using randomized acquisition. Extraction of physical quantitative tissue parameters is performed off-line, based on acquisition with varying parameters and a dictionary generated according to the Bloch equations. MRF uses hundreds of radio frequency (RF) excitation pulses for acquisition, and therefore high under-sampling ratio in the sampling domain (k-space) is required for reasonable scanning time. This under-sampling causes spatial artifacts that hamper the ability to accurately estimate the tissue's quantitative values. In this work, we introduce a new approach for quantitative MRI using MRF, called magnetic resonance Fingerprinting with LOw Rank (FLOR). We exploit the low rank property of the concatenated temporal imaging contrasts, on top of the fact that the MRF signal is sparsely represented in the generated dictionary domain. We present an iterative scheme that consists of a gradient step followed by a low rank projection using the singular value decomposition. Experiments on real MRI data, acquired using a spirally-sampled MRF FISP sequence, demonstrate improved resolution compared to other compressed-sensing based methods for MRF at 5% sampling ratio.


I. INTRODUCTION
Quantitative Magnetic Resonance Imaging (QMRI) is widely used to measure tissue's intrinsic spin parameters such as the spin-lattice magnetic relaxation time (T1) and the spinspin magnetic relaxation time (T2) [1]. Since tissue relaxation times vary in disease, QMRI enables the diagnosis of different pathologies, including multiple sclerosis (MS), Alzheimer, Parkinson, epilepsy and cancer [2], [3], [4], [5], [6], [7]. In addition, the knowledge of tissue relaxation times allows generation of many clinical MR imaging contrasts (such as FLAIR and STIR) off-line, and saves significant amount of scanning time.
Despite the advantages of QMRI, clinical MRI today mainly consists of weighted images. Values in weighted MR imaging are given in arbitrary units, since the signal strength is influenced by both intrinsic parameters (such as relaxation times and concentration of hydrogen atoms) and non-intrinsic ones. Non-intrinsic parameters include transmit and receive coils sensitivities, patient position in the scanner, vendor based scanner specific parameters, and local temperature. This is why weighted MRI images lack quantitative information and as a result, different materials may exhibit similar or identical values. In addition, weighted MRI contrast values vary between different follow-up scans of the same patient. This fact may impair disease monitoring, if based solely on those images. To date, weighted MRI scans are more common than QMRI in the clinic, due to the extremely long time required for QMRI acquisition using conventional techniques [8], [9], [10].
A plethora of methods have been proposed for QMRI. Earlier approaches are based on a series of spin echo (SE) or inversion recovery (IR) images with varying repetition times (TR) and echo times (TE) to evaluate each magnetic parameter (T1 or T2) separately. After acquisition, the curve of intensities for each pixel is matched to the expected magnetic signal, representing the appropriate magnetic tissue parameters [8]. Over time, accelerated methods for QMRI were suggested, such as the popular driven equilibrium single pulse observation of T1 (DESPOT1) [9] or T2 (DESPOT2) [10] and the IR TrueFISP for simultaneous recovery of T1 and T2 quantitative maps [11], [1]. Both techniques do not require long waiting times between excitations to get to an equilibrium state, and therefore they are significantly faster. However, obtaining accurate and high resolution QMRI in a reasonable clinical scanning time is still very challenging.
Recently, an approach for QMRI called magnetic resonance fingerprinting (MRF) has been proposed [12]. MRF uses pseudo-randomized acquisitions to generate many different imaging contrasts, acquired at a high under-sampling ratio. It exploits the different acquisition parameters over time to produce a temporal signature, a "fingerprint" for each material under investigation. By matching this unique signature to a pregenerated set of simulated patterns, the quantitative parameters can be extracted off-line. This approach saves valuable scan time compared to previous methods for accelerated QMRI, and has been proven to be efficient and reliable [12].
MRF utilizes the fact that each tissue responds differently to a given quasi-random pulse sequence. By varying the acquisition parameters (e.g. repetition time (TR), echo time (TE), and radio frequency flip angle (FA)), unique signals are generated from different tissues. After acquisition, a pattern recognition algorithm is used to match the acquired signal to an entry from a dictionary of possible tissue candidates. The dictionary entries are created by simulating the tissue's response to the sequence for a range of T1 and T2 parameter values, using the Bloch equations. The resulting dictionary contains the temporal signatures of various simulated materials, given the pseudo-random pulse sequence. The quantitative parameters, such as the tissue's T1 and T2 relaxation times, can be retrieved from the data by matching the signature acquired to the most correlated entry in the dictionary.
In MRI, data is acquired in Fourier conjugate of the spatial image domain (a.k.a. k-space), where the acquisition time In MRF, we acquire many under-sampled images over time.
When reconstructing by zero filling, the under-sampled data is blurred and contains aliasing artifacts.
of a high resolution, single contrast 3D MRI lasts a few minutes. Since MRF is based on rapid acquisition of hundreds of different contrasts, severe under-sampling is needed in the k-space domain to obtain the temporal resolution required for MRF. In the original MRF paper, reconstruction from under-sampled k-space data is preformed by zero-filling of the missing k-space values [12]. Figure 1 demonstrates the effect of fully sampled versus under-sampled data, acquired with spiral trajectories and recovered by zero-filling. It can be seen that the under-sampled data is blurred and introduces aliasing artifacts. Figure 2 illustrates the noise and under-sampling artifacts of a representative brain voxel intensity as function of time, acquired with an MRF approach based on the fast imaging with steady state precession (FISP) sequence [13]. Clearly, under-sampling also introduces a substantial level of noise in the time domain.
While in the original MRF paper reliable results were presented based on zero-filling reconstruction, exploiting the nature of the signal for improved recovery leads to better performance at higher under-sampling ratios. Therefore, later works added sparse representation and compressed sensing (CS) methods [14], [15]. Davies et al. [16] developed an approach called BLIP (BLoch response recovery via Iterative Projection) that is based on gradient descent and projection onto the dictionary subspace. Wang et al. [17] suggested a method that exploits sparsity in the wavelet domain of each imaging contrast, together with changing the acquisition trajectories between time stamps. Although these techniques have shown better resolution with smaller error compared to the original MRF approach, they do not exploit the temporal similarity across the time domain, which is intrinsic to the MRF acquisition scheme. Such similarity is exploited, via modelling the data as low-rank, in many other MRI applications with high temporal resolution, including cardiac imaging [18] and functional MRI [19].
An initial work describing the implementation of a lowrank model for MRF was developed recently [20] and has also been presented by us in [21]. Here, we apply an improved reconstruction algorithm and include mathematical justification of our approach, together with under-sampling using spiral trajectories. In particular, we exploit the low-rank property of the temporal domain of MRF, via an iterative scheme that consists of a gradient step followed by a projection onto the subspace spanned by the dictionary elements in order to constraint the structure of the tissue behaviour simulated in the dictionary. The estimated images are then decomposed using the singular value decomposition (SVD) and the singular values are soft-thresholded to reduce the nuclear norm value in every step. We show that our algorithm, called magnetic resonance Fingerprinting with LOw Rank constraint (FLOR), increases the accuracy of the resulting parameter maps from highly under-sampled data, compared to previous CS-based methods. Specifically, we demonstrate lower mean square error of the quantitative maps using only 5% of the data, and with no additional computation time.
This paper is organized as follows. Section II describes the MRF problem and provides a review of common reconstruction methods, followed by our low-rank based approach. Section III compares our results to previous MRF algorithms, using realistic MRI data of a human subject. Sections IV and V discuss experimental results using retrospectively undersampled MR data, followed by conclusions.

A. Problem formulation
MRF data consists of multiple frames, acquired in the image's conjugate Fourier domain (a.k.a k-space), where each frame results from different acquisition parameters. We stack the measurements into a Q × L matrix Y, where L is the number of frames and Q is the number of k-space samples in each frame. Every column in Y is an under-sampled Fourier transform of an image frame, X :,i : where F u {·} denotes an under-sampled 2D Fourier transform.
The row X j,: represents the temporal signature of a single pixel (assumed to correspond to a single tissue). The signature depends on the tissue's relaxation times, T1 and T2, and its proton density (PD), grouped as a row vector: Each column, X :,i represents a response image acquired at a single time point with different acquisition parameters, stacked as a column vector: where TR and TE are the repetition time and time to echo and FA represents the flip angle of the RF pulse. Therefore, where f {·} represents the Bloch equations. Note that we omit the off resonance parameter (which appeared in the original MRF paper [12]), since the sequence used in our experiments is derived from the FISP sequence, which is insensitive to off resonance effects [13].
The goal in MRF is to recover, from the measurements Y, the imaging contrasts X and the underlying quantitative parameters of each pixel defined in (2), under the assumptions that every pixel in the image contains a single type of tissue and that Θ 2 is known.
Recovery is performed by defining a dictionary that consists of simulating the signal generated from M tissues using the Bloch equations (represented as M different combinations of T1 and T2 relaxation times), when the length-L acquisition sequence defined in (3) is used. As a result, we obtain a dictionary D of dimensions M × L (M > L as the number of simulated tissues is greater than the sequence length). The PD is not simulated in the dictionary, as it is the gain used to match the Bloch simulation performed on a single spin to the signal obtained from a pixel containing multiple spins. It can be easily determined after the T1 and T2 maps are known. After successful recovery of X, each row in X is matched to a single row in the dictionary, and T1 and T2 are estimated as those used to generate the matched dictionary row. Each dictionary signature has its own unique T1 and T2 values stored in a look up table (LUT), represented as the matrix LUT of dimensions M × 2.

B. Previous Methods
The approach suggested in the original MRF paper [12] is described in Algorithm 1, and uses matched filtering to match dictionary items to the acquired data. In the algorithm, F H {·} is a zero-filling 2D inverse Fourier transform (performed by filling the unknown frequencies in the Fourier domain with zeros and applying the 2D inverse Fourier transform). The parameters, k j are the matching dictionary indices, j is a spatial index and i is the temporal index, representing the ith frame in the acquisition. The parameter maps are extracted from LUT, which holds the values of T1 and T2 for each k j . This approach does not incorporate sparse based reconstruction, which has been proven to be very successful in MRI applications based on under-sampled data [22], [23], [24].

Algorithm 1 Original MRF algorithm Input:
A set of under-sampled k-space images: Y A pre simulated dictionary: D An appropriate look up table: LUT Output: Magnetic parameter maps: T 1 , T 2 , P D Compute for every i and j: Davies et al. [16] suggested a method incorporating sparsity of the data in the dictionary domain (i.e. each pixel is represented by at most one dictionary item), referred to as the BLoch response recovery via Iterative Projection (BLIP) algorithm. This approach is based on the Projected Landweber Algorithm (PLA) which is an extension of the popular iterative hard thresholding method. BLIP (described here as Algorithm 2) consists of iterating between two main steps: A gradient step that enforces consistency with the measurements, and a projection that matches each row of X to a single dictionary atom.

Input:
A set of under-sampled k-space images: Y A pre simulated dictionary: D An appropriate look up table: LUT Output: Magnetic parameter maps: T 1 , T 2 , P D Initialization: µ, X 0 = 0 Iterate until convergence: • Gradient step for every i: Z n+1 onto the dictionary subspace for every j: Restore maps for every j: Wang et al. [17] presented a different approach that enforces sparsity in the wavelet domain of each imaging frame, X :,i .
They also suggested to use a distance metric learning for determining which dictionary atom best matches the signature from acquired data. Specifically, they suggested replacing the Euclidean norm with the Mahalanobis distance. Algorithm 3 describes the basic concept of this method, where Ψ indicates a wavelet transform.

Algorithm 3 Wavelet based MRF recovery Input:
A set of under-sampled k-space images: Y A pre simulated dictionary: D An appropriate look up table: LUT Output: Magnetic parameter maps: T 1 , T 2 , P D Initialization: µ, X 0 :,i = F H {Y :,i } for every i Solve the following minimization problem for every i: Project onto the dictionary subspace for every j: Restore maps for every j: BLIP and a few other works that are based on it (e.g. [17]) do not incorporate the temporal similarity across time points, which is a fundamental nature of the MRF sequence. In addition, there is a high degree of similarity across signatures in D. As a result, the imaging contrasts matrix X is typically a low-rank matrix.
To demonstrate the low-rank property of X, we simulated the acquisition of a real brain with known T1, T2 and PD maps of size 128 × 128 with the FISP sequence [13] and L = 500 TRs. These acquisition parameters were chosen based on previous publications in the field of MRF [12], [13]. Note that the general assumption of X being a low-rank matrix holds as long as temporal similarity exist between timeframes in X, and multiple voxels in the image belong to a single tissue, regardless of the specific acquisition parameters. Figure 3 shows the singular values of X. It can be seen that X is indeed low-rank, as most of the data is represented in the highest 15 singular values.
Using low rank as a prior on MRI data has been proven to be extremely efficient in many other MRI-related tasks [25], [19]. Here we show its contribution in recovering highly under-sampled MRF data. This low-rank property of X can be exploited for improved reconstruction using the following optimization problem: where R is a matrix that matches each pixel (X j,: ) with the dictionary signatures. A priori, it would make sense to constraint R to have one sparse rows that contain the corresponding PD value for each row of X, since only a single dictionary item should match an acquired signature. However, in practice we found that superior results are obtained by relaxing this constraint, and allowing X to be comprised of multiple dictionary elements throughout the iterations. This permits non-simulated signatures to be described by a linear combination of simulated ones, and also preserves the convexity of the equality constraint. The one-sparsity is enforced only at the final stage after X is recovered by using a matched filter (MF) in order to extract the parameter maps. For brevity we write the constraint X = RD as X ∈ D where D = {X : N (X) ⊇ N (D)}. The parameter r is the rank of the matrix, defined as a fixed pre-chosen parameter. Since typically r is not known in advance, we consider a regularized version of (4): for some fixed regularization parameter λ.
Problem (5) is not convex due to the rank constraint. We therefore relax this constraint by replacing the rank of X with the nuclear norm X * , defined as the sum of the singular values of X [26]. This results in the relaxed problem: In order to solve (6) we use the incremental subgradient proximal method [27] as described in Appendix A. Our minimization algorithm, referred to as magnetic resonance Fingerprint with LOw Rank (FLOR), is detailed in Algorithm 4, where the parameter λ is chosen experimentally. Note that by setting λ = 0 and enforcing R to have one-sparse rows, FLOR reduces to BLIP [16]. While Algorithm 4 describes the basic implementation of a low-rank based approach to MRF, we can add an improvement that significantly reduces convergence time. The improvement uses the acceleration approach suggested by Nesterov [28] for minimizing a smooth convex function, and its extension for non smooth composite functions of Beck and Teboulle [29], [30]. The modified algorithm is described in Appendix B. Figure 4 shows the reconstruction error of FLOR as the number of iterations varies with and without this addition.

Algorithm 4 FLOR -MRF with LOw Rank
Re D k j , Xj,: III. EXPERIMENTAL RESULTS This section describes MRI experiments that were carried out using a brain scan of a healthy subject. Data was acquired with GE Signa 3T HDXT scanner. The experimental procedures involving human subjects described in this paper were approved by the Institutional Review Board of Tel-Aviv Sourasky Medical Center, Israel. To generate gold standard data we acquired Fast Imaging Employing Steady-state Acquisition (FIESTA) and Spoiled Gradient Recalled Acquisition in Steady State (SPGR) images, at 4 different flip angles (3 • ,5 • ,12 • and 20 • ), and implemented well known DESPOT1 and DESPOT2 [31] algorithms, after improvements as described in [32], to generate T1,T2 and PD quantitative maps. In addition, the FISP pulse sequence was simulated with constant TE of 2ms, random TR values in the range of 11.5-14.5 ms, and a sinusoidal variation of FA (RF pulses) in the range of 0-70 degrees [13].
To simulate noisy samples, we added complex Gaussian zero-mean noise with standard deviation of σ = 0.5 to the kspace data to obtain SNR of 56. 5 [33]. The tuning parameters were experimentally set as µ = 1 and λ = 5. Data was fed as an input to BLIP, wavelet based reconstruction algorithm and the improved FLOR algorithm (those are described as Algorithms 2,3 above and Algorithm 5 in Appendix B). In addition, we performed reconstruction using 100% of the data (without the addition of noise) via conventional MRF (Algorithm 1), to evaluate the inherent error caused by the discretized dictionary. All the iterative algorithms were run until the difference between consecutive iterations was below the same threshold.
For quantitative error analysis, we calculated the normalized MSE (NMSE) between each quantitative map estimation and the gold standard, defined as: where θ i ,θ i represent a gold standard map (such as T1,T2 or PD) and its corresponding reconstructed map (respectively), N is the number of pixels in the map and j is a spatial index. Non-uniform forward Fourier transform was performed using NUFFT [34], and nonuniform inverse Fourier transform was performed by SPURS [35]. The MATLAB code for reproducing the experiment provided in this section can be be found at: http://webee.technion.ac.il/Sites/People/YoninaEldar/software det18.php. In this code, spiral sampling was based on [36]. Figure 5 shows the resulting maps for the recovery of T1, T2 and PD obtained with the various algorithms against the gold-standard (left). The corresponding error maps of each algorithm versus the gold standard are shown in Figure 7. To allow detailed view of the reconstruction results for the reader, Fig. 6 shows a zoomed region for each map. Please note that these Figures also include the results of the original MRF paper, obtained using 100% of the data without the addition of noise. These results are provided to give the reader a baseline   [ms] Fig. 7: Error maps of T1, T2 in milliseconds and PD in arbitrary units reconstructions. Left: reconstruction using conventional MRF from 100% of the noise-free data, followed by wavelet based algorithm, BLIP and low-rank MRF reconstructions with acceleration from 5% of the noised data.

A. Reconstruction results
for comparison, and to get the sense of the quantization artifacts arising from the discretization of the dictionary. It can be seen that both FLOR and BLIP outperform wavelet-based reconstruction results, when using 5% of sampled data. In addition, FLOR provides a significantly lower error compared to BLIP and comparable to the error obtained by the original MRF algorithm using 100% of the data. Due to the very low sampling ratio in our experiments, the conventional MRF approach result using 5% of the data did not provide valuable reconstruction and is therefore omitted in this analysis. We expand the discussion on this issue further in Section IV.A.

B. Parameter sensitivity analysis
To expand the performance analysis, we evaluated the algorithm's results with different parameter settings. We simulated data acquisition at various noise levels and sampling ratios. This analysis was performed with 2D random under-sampling generated from a polynomial distribution of order 4 [22], to save expensive reconstruction time from a non-uniform sampling process. Figure 8 shows the results of these experiments.
It can be seen that while FLOR provides the lowest NMSE, the improvement of FLOR versus BLIP is reduced as the sampling rate grows. As to the noise, FLOR is more robust to noise and the improvement versus BLIP increases as the noise variance grows. This phenomenon occurs since the low-rank mechanism eliminates low singular values that are related to noise, while BLIP does not handle noise reduction explicitly. It can be seen that the wavelet based algorithm demonstrates robustness to noise, yet it maintains higher reconstruction error.
Next, we checked the tolerance of FLOR to the chosen threshold parameter at the low rank projection stage. We set µ = 1, and examined the results as the value of λ varies. Figure 9 demonstrates the NMSE of the reconstructed images X as a function of the chosen λ. It can be seen that FLOR has an optimum value for the parameter λ, yet it has some tolerance for a range of optional values.
The parameter sensitivity analysis behaves as expected for the experiments that examine the noise variance and undersampling ratio (Figure 8). Generally, the NMSE drops as the noise variance decreases or as the sampling ratio increases. There is an unexpected phenomenon in the T2 NMSE vs.  Fig. 8). This scenario demonstrates improved T2 map reconstruction by FLOR and by BLIP, despite the increase in the noise level. This phenomenon is explained by the MF projection, that adds the quantization error from the dictionary. The NMSE shown on those graphs is calculated from the quantitative maps after the matched filter projection. This projection actually hampers the expected decrease in the NMSE as the noise level increases. We can report that the NMSE of the recovered contrasts, X, (before computing the quantized maps) behaves as expected. Moreover, note that although the T1 NMSE vs. sampling ratio graph (the left bottom plot in Fig. 8) demonstrates close values between BLIP and FLOR, FLOR demonstrates better resolution in all images when examined visually.

A. Relation to previous works
Although an initial work that exploits the low rank structure of MRF sequences has been published before as a short conference paper [20], here we propose a different algorithm for solving the problem. Our solution is based on soft-thresholding the singular values [26], which is mathematically justified in Appendix A. We also use the FISP sequence to minimize the effect of B0 inhomogeneity. Moreover, we compare the results of our algorithm to popular CS methods for MRF, and demonstrate superior results. Finally, as opposed to other CS-based papers for MRF, we examine our algorithm on quantitative maps of a real subject.
While BLIP treats the original MRF problem as an ℓ 0 optimization problem, FLOR solves first the relaxed problem of (5) and only then uses MF to extract the magnetic parameters. It leads to some beneficial properties such as convergence guarantee, and the ability to use the accelerated version described in Appendix B. The results of the original MRF paper are obtained by applying Algorithm 1 on our data. Based on our experiments, the original MRF algorithm did not converge to a meaningful result. The maps shown in our figures can be reproduced by using our code, which can be found at: http://webee.technion.ac.il/Sites/People/YoninaEldar/software det18.php.
The original MRF paper [12] describes a multi coils (32 head coils) acquisition scenario, while in this work we simulated reconstruction from a single coil. It significantly reduces the amount of acquired information and explains why our original MRF reconstruction failed to provide meaningful results (but yet, our single coil-based simulations provide a fair comparison between the various algorithms). Fig. 10: Comparison of FLOR recovery with acceleration step (FLOR I), and FLOR recovery with acceleration step and interpolate matching (FLOR II) out of 5% of the data. T1 and T2 maps are in milliseconds and PD is in arbitrary units. As can be seen in the magnified area, FLOR II demonstrates a smoother solution (with less quantization errors) which is more similar to the original gold standard maps.

B. Computational complexity
FLOR is divided into two main components: The first one recovers the imaging contrasts, and the second extracts the parameter maps from the recovered contrasts. The computational burden of FLOR lies in the low-rank projection step, or specifically, in the SVD calculation. This step does not exist in the other algorithms tested in this paper. However, there are efficient fast techniques to calculate the SVD [37], required by FLOR. Moreover, unlike BLIP, FLOR does not require a MF calculation at every iteration. While the total computation time of the algorithms depends on the image size, for images of size 128 × 128 examined in this paper FLOR provides the fastest convergence.
A time consuming step that exists in all algorithms is the non uniform Fourier transform. Previous implementations of CS-based reconstruction algorithms mainly use the inverse NUFFT (iNUFFT) algorithm [34]. Here we use SPURS, which is a faster approach for the non uniform transform [35]. Based on our observations, SPURS provides improved image reconstruction with the same computational complexity compared to iNUFFT. A possible extension of FLOR is to add additional signatures to the dictionary by linear interpolation, in regions where a few candidates from the dictionary match a single signature from the data. Then, we select the dictionary signatures that exhibit high correlation value (the ones above a certain threshold) and average their matching T1 and T2 values. Thanks to the fact that this improvement expands the possible solutions to include ones that do not exist in the dictionary, it may exhibit improved accuracy compared to the conventional matching. We have implemented this improvement, and the results are shown in Figure 10, with corresponding error maps in Figure 11. In those figures, we compare the recovery maps of FLOR without (FLOR I) and with (FLOR II) the proposed improvement. It can be seen that FLOR II improves the results of FLOR I and produces a smoother solution which better fits the gold standard maps.
V. CONCLUSIONS We presented FLOR, a method for high quality reconstruction of quantitative MRI data using MRF, by utilizing the low-rank property of MRF data. Due to the fact that we exploit low-rank on top of the well known sparsity of MRF in the dictionary matching domain, we are able to obtain high quality reconstruction from highly under-sampled data. We provide results that are comparable to fully sampled MRF, using only 5% of the data. In addition, comparison against CSbased methods for MRF shows the added value of low-rank based reconstruction. Future work will verify the robustness of the algorithm for prospectively sampled MRF data and will examine more complicated patch wise recoveries.
VI. APPENDIX A The basic implementation of FLOR, as described in Algorithm 4 in the paper, aims to solve the following optimization problem: where F u is the partial Fourier transform operator, X has dimensions N 2 × L and D = {X : N (X) ⊇ N (D)}. FLOR solves (8) using the incremental proximal method [27], which treats problems of the form: . The function f i (X) is convex and non-differentiable, h i (X) is a convex function and D is a non-empty, closed, and convex subspace. The general step in solving (9) is given by [27, (4.12)-(4.13)]: where g i k ∈ ∂h i k (X k ), µ k is a positive step size, and P D is the projection operator onto D defined as The optimization problem, defined in the update step of X k+1 , is referred to as the proximal gradient calculation of the nondifferentiable f i k , under the constraint X ∈ D.
The solution of (10b) for f (X) = λ X * without the constraint X ∈ D, is the singular value soft-thresholding operator (SVT) [26] defined as: Here Σ r is a diagonal matrix with the non-zero singular values of Z k on its diagonal, U r and V r are the r left and right singular vectors of the SVD of Z k , associated with the r nonzero singular values, and [x] + = max(0, x). In our case, since Z k ∈ D (as follows from (10a)) and the SVT calculation keeps the operand in the same subspace, the constraint X ∈ D can be omitted. Therefore, (10b) reduces to Combining (13), (16) and (14), the incremental subgradientproximal method for solving (8) consists of two updates in each iteration: This constitutes the core of Algorithms 4 and 5 (in Appendix B). In our framework, the step sizes are set to constant, µ k = µ, and λ is chosen experimentally.

VII. APPENDIX B
Algorithm 4 in the paper describes the basic implementation of low rank for MRF. The accelerated version, mentioned in the paper as accelerated FLOR, is described in Algorithm 5. The difference is in the update of X, where in Algorithm 4 it is updated directly from the low rank solution, and now it is updated with a linear combination between this projection and the one from the iteration before.

P D j = max
Re D k j , Xj,: