Accelerating Magnetic Resonance $\text{T}_{1\rho}$ Mapping Using Simultaneously Spatial Patch-based and Parametric Group-based Low-rank Tensors (SMART)

Scalars are denoted by italic letters (e.g., $b$). Vectors are denoted by bold lowercase letters (e.g., b ). Matrices are denoted by bold capital letters (e.g., B). Tensors of order three or higher are denoted by Euler script letters (e.g., $\mathcal{X}$). Indexed scalar $x_{ijk}$ denotes the $(i,j,k)$th element of third-order tensor $\mathcal{X}$. Operators are denoted by capital italic letters (e.g., $C$). $\ C^{T}$ is used to denote the adjoint operator of $C$.

The mode-$i$ matricization (unfolding) of a tensor $\mathcal{X}\in\mathbb{C}^{N_{1}\times N_{2}\times\ldots\times N_{d}}$ is defined as $\textbf{X}_{(i)}\in\mathbb{C}^{N_{i}\times\left(N_{1}\ldots N_{i-1}N_{i+1}\ldots N_{d}\right)}$, arranging the data along the $N_{i}$th dimension to be the columns of $\textbf{X}_{(i)}$. The opposite operation of tensor matricization is the folding operation, which arranges the elements of the matrix $\textbf{X}_{(i)}$ into the $d$th order tensor $\mathcal{X}$. Tensor multiplication is more complex than matrix multiplication. In this study, only the tensor $i$-mode product[35] is considered, which is defined as the multiplication (e.g., $\mathcal{X}\times_{i}\mathbf{U}$) of a tensor $\mathcal{X}$ with a matrix $\textbf{U}\in\mathbb{C}^{J\times N_{i}}$ in mode $i$. Elementwise, we have $w_{n_{1}\cdots n_{i-1}jn_{i+1}\cdots n_{d}}=\sum_{n_{i}=1}^{N_{i}}x_{n_{1}\cdots n_{i-1}n_{i}n_{i+1}\cdots n_{d}}u_{jn_{i}}$, where $\mathcal{W}=\left(\mathcal{X}\times{}_{i}\mathbf{U}\right)\in\mathbb{C}^{N_{1}\times\ldots\times N_{i-1}\times J\times N_{i+1}\times\ldots\times N_{d}}$. The $i$-mode product can also be calculated by the matrix multiplication $\mathbf{W}_{(i)}=\mathbf{U}\mathbf{X}_{(i)}$.

Let $\textbf{X}\in\mathbb{C}^{N_{\textit{voxel}}\times N_{\textit{TSL}}}$ be the $\text{T}_{1\rho}$-weighted image series to be reconstructed, where $\textit{N}_{\textit{voxel}}$ denotes the voxel number of the $\text{T}_{1\rho}$-weighted image, and $\textit{N}_{\textit{TSL}}$ is the number of TSLs. $\mathbf{Y}\in\mathbb{C}^{N_{\textit{voxel}}\times N_{c}\times N_{\textit{TSL}}}$ denotes the corresponding k-space measurement, where $N_{c}$ denotes the coil number. The forward model for $\text{T}_{1\rho}$ mapping is given by

$$\mathbf{Y}=E\mathbf{X}+\zeta,$$

(1)

where $\zeta$ denotes the measurement noise and $E$ denotes the encoding operator[36, 37] given by $E=AFS$, where $A$ is the undersampling operator that undersamples $k$-space data for each $\text{T}_{1\rho}$-weighted image. $F$ denotes the Fourier transform operator. $S$ denotes an operator which multiplies the coil sensitivity map by each $\text{T}_{1\rho}$-weighted image coil-by-coil. The general formulation for recovering $\mathbf{X}$ from its undersampled measurements can be formulated as

$$\begin{gathered}\underset{\textbf{X}}{\operatorname{\arg\min}}\frac{1}{2}\left\|E\mathbf{X}-\mathbf{Y}\right\|_{F}^{2}+\lambda R(\mathbf{X})\end{gathered},$$

(2)

where $\left\|\cdot\right\|_{F}$ denotes the Frobenius norm, $R(\mathbf{X})$ denotes a combination of regularizations, and $\lambda$ is the regularization parameter.

In this study, the SMART reconstruction method assumes that the $T_{1\rho}$-weighted image series $\mathbf{X}$ can be expressed as a high‐order low‐rank representation on a spatial patch-based tensor and a parametric group-based tensor (shown in Fig. 1). The problem in (2) can be formulated as a constrained optimization on the two high‐order low‐rank tensors.

The SMART method performance is affected by the parameters of the spatial and parametric tensor construction.These include patch size b, maximum similar patch number $N_{p,max}$, normalized $l_{2}$-norm distance threshold $\lambda_{m}$ in the spatial tensor construction, and the group number $\textit{N}_{g}$ in the parametric tensor construction. These parameters should be carefully tuned to obtain the best reconstruction performance.

$\rm T_{1\rho}$ maps were obtained by fitting the reconstructed $\rm T_{1\rho}$-weighted images with different TSLs pixel-by-pixel:

$$M_{k}=M_{0}\exp\left(-t_{k}/T_{1\rho}\right)_{k=1,2,\ldots,N_{\textit{TSL }}},$$

(15)

where $M_{0}$ denotes the baseline image intensity without applying a spin-lock pulse, and $M_{k}$ is the signal intensity for the kth TSL image. $\text{T}_{1\rho}$ map was estimated using the nonlinear least-squares fitting method with the Levenberg–Marquardt algorithm[48] from the reconstructed $\text{T}_{1\rho}$-weighted images.

Numerical phantom, real phantom, and in vivo experiments were performed to evaluate the tissue clustering method performance. All MR scans involved were performed on 3T scanners. The details of data acquisition are shown in the supplementary information and Table S1.

The proposed method was evaluated on the $\text{T}_{1\rho}$ mapping datasets of the brain collected from six volunteers on 3T scanners[20, 53]. The local institutional review board approved the experiments, and informed consent was obtained from each volunteer. The details of data acquisition are shown in supplementary information and Table S1.

For the 2D datasets, the fully sampled k-space data were retrospectively (Retro) undersampled along the ky dimension using a pseudo-random undersampling pattern [54] with acceleration factors (R) = 4 and 6. For the 3D datasets, the fully sampled k-space datasets were retrospectively undersampled using Poisson disk random[55] patterns with R = 6.76 and 9.04. The sampling masks for each TSL were different. For the prospective (Pro) study , two 2D datasets were acquired from one volunteer (27 years old) with R = 4.48 and 5.76, respectively. Fully sampled data were also acquired as a reference in the prospective experiment.

SMART and four state-of-the-art methods were used to reconstruct the undersampled k-space data. These methods include the high-dimensionality undersampled patch-based reconstruction method (PROST) [45] using the spatial patch-based low-rank tensor, low-rank plus sparse (L + S) method[36] enforcing the spatial global low-rank and sparsity of the image matrix, locally low-rank method (LLR)[13] using the spatial patch-based low-rank matrix, and model-driven low rank and sparsity priors method (MORASA)[29] which enforces the global low-rankness of the matrix in both spatial and temporal directions and sparsity of the image matrix. In addition, a modified PROST reconstruction method (PROST + HM) was also compared with SMART to verify the effectiveness of the parametric group-based low-rank tensor in image reconstruction. This jointly uses the spatial patch-based low-rank tensor and a parametric low-rank matrix by replacing the patch-based parameter tensor in SMART with a Hankel matrix. Experiments with high acceleration factors of 10.2 and 11.7 in the 2D scenario and 11.26 and 13.21 in the 3D scenario were also performed to verify the effectiveness of the proposed method.

$\text{T}_{1\rho}$ quantification was assessed using the normalized root mean square error (nRMSE)[49]. The quality of the reconstructed $\text{T}_{1\rho}$-weighted images was assessed based on the nRMSE, high-frequency error norm (HFEN)[56, 57], structural similarity index measure (SSIM)[58], and peak signal-to-noise ratio (PSNR).

In this subsection, several different values of the parameters related to the spatial and parametric tensor construction mentioned in the previous section were set and used for 2D image reconstruction. The nRMSE values for the reconstructions of SMART with R = 4 and R = 6 using different parameters were compared.

The representative zero-filling reconstruction of the Retro 2D dataset with R = 6, the intermediate reconstruction at the 12th iteration, the parameter histogram, the estimated $\text{T}_{1\rho}$ map, the corresponding tissue groups, and the ranks of tissue groups were also compared to demonstrate the necessity of updating tissue maps in the iteration.

All reconstructions, $\text{T}_{1\rho}$ estimation, and analyses were performed in MATLAB 2017b (MathWorks, Natick, MA, USA) on an HP workstation with 500 GB DDR4 RAM and 32 cores (two 16-core Intel Xeon E5-2660 2.6 GHz processors). The reconstruction parameters used in all methods are presented in the supplementary information Table S2. The tensor transforms, including the folding and unfolding operators, and the tensor multiplication were implemented using the MATLAB tensor toolbox by Brett et al. from Sandia National Laboratories¹¹1http://www.tensortoolbox.org/. The HOSVD using algorithm 1 of the supplementary information was implemented based on the transforms from the above toolbox.

We conducted an ablation study to analyze the spatial and parametric tensors’ effects on the reconstruction. Two models are constructed as follows:
model 1:

$$\underset{\mathbf{X}}{\arg\min}\frac{1}{2}\|E\mathbf{X}-\mathbf{Y}\|_{F}^{2}+\lambda_{1}\sum_{i}\left\|\mathcal{T}_{i}\right\|_{*}\text{ s.t. }\mathcal{T}_{i}=P_{i}(\textbf{X})$$

(17)

model 2:

$$\underset{\mathbf{X}}{\arg\min}\frac{1}{2}\|E\textbf{X}-\mathbf{Y}\|_{F}^{2}+\lambda_{2}\sum_{j}\left\|\mathcal{Z}_{j}\right\|_{*}\text{ s.t. }\mathcal{Z}_{j}={H}_{j}(\mathbf{X})$$

(18)

Model 1 only includes the spatial tensor, and model 2 only includes the parametric tensor. Model 1 is the same as the PROST method. Similarly, we applied ADMM to solve the optimization problem in model 2. These two models were used to reconstruct images from the retrospectively undersampled in vivo dataset with R = 4 and 6.

$\text{T}_{1\rho}$-weighted images at TSL = 80 ms of the numerical phantom with different SNRs are shown in Fig.2(a) to represent low, medium, and high SNRs. Fig.2(b) shows the rank curves of the numerical and real phantom, and in vivo datasets. The first row of Fig.2(b) is for the fully sampled dataset, while the second is for reconstructing the undersampled data.

Fig. 6 shows the prospective reconstructed $\text{T}_{1\rho}$-weighted images (at TSL = 40 ms) from another volunteer and the magnified images using the SMART, PROST + HM, PROST, L + S, LLR, and MORASA methods. Visual artifacts can be observed in the magnified images of reconstructions at all accelerating factors using the L + S and LLR methods. The images reconstructed using the PROST + HM and PROST methods were blurred compared to those reconstructed using the SMART method. Some image details (namely, the blood vessel marked in a yellow arrow) can narrowly be seen in the reconstructions using the PROST + HM and MORASA methods and disappeared in the reconstruction using the PROST method at R = 5.76. In contrast, the image details were well preserved using the SMART method. The $\text{T}_{1\rho}$ maps reconstructed using the methods above are shown in supplementary information Fig. S4. Similar conclusions can be drawn from the $\text{T}_{1\rho}$ maps.

Fig.9 shows the effects of four parameters of the spatial and parametric tensor construction in 2D imaging. As shown in Fig. 9(a), the smallest nRMSE is achieved at $b=9$. When the patch size decreases from 9 to 3, the nRMSE increases rapidly due to the small image patches’ limited spatial modeling ability. When the patch size is larger than 9, the reconstruction performance slowly degrades since the discrepancy among the patches increases, causing the reduced low rankness of the constructed spatial tensor. The nRMSE of reconstruction decreases with the increase of $N_{p,max}$, as shown in Fig. 9(c), but the trend slows down when $N_{p,max}\geq 30$. Considering the computational complexity, we chose $N_{p,max}=30$ in this study. Fig. 9(d) shows that $\lambda_{m}$ needs to be larger than 0.15; otherwise, the nRMSE will significantly increase. As shown in Fig. 9(b), changing $N_{g}$ has little effect on the nRMSE of the reconstruction for R = 4. In contrast, large $N_{g}$ shows improved performance for the higher acceleration factor of 6. Therefore, $N_{g}=60$ was used in this study.

Supplementary information Fig. S7 shows the representative zero-filling reconstruction with R = 6, the intermediate reconstruction at the 12th iteration, the parameter histogram, the estimated $\text{T}_{1\rho}$ map, the corresponding tissue groups, and the ranks of tissue groups. The initial $\text{T}_{1\rho}$-weighted image of the zero-filling reconstruction shows strong aliasing artifact, and the $\text{T}_{1\rho}$ map estimated from zero-filling images is quite blurred. The quality of $\text{T}_{1\rho}$-weighted image and $\text{T}_{1\rho}$ map was significantly improved at the 12th iteration. After iteration, the distribution of the parameter histogram also changed (shown in Fig. S7(b)). In the figure of tissue groups, the pixel intensities along the pixel number direction are not as uniform as those at the 12th iteration. The ranks in several groups are reduced from 2 to 1.

Fig. 10 shows the reconstructed images for models 1 and 2 with R = 4 and 6, respectively. For the results of model 1, the aliasing artifacts due to undersampling can be well removed from the reconstructed images by applying the spatial tensor, but the reconstructed images are blurry. In contrast, the sharpness of images reconstructed by model 2 is improved, but aliasing artifacts still exist. These results are in line with our hypothesis that the parametric group-based tensor can help improve image sharpness and preserve more image details in the reconstructions.

In this study, all the reconstruction parameters included in the ADMM and CG algorithms were selected empirically. The relative change in the solution between two consecutive iterations was used to show the convergence numerically and is defined as follows:

$$\left\|\mathbf{X}_{\textit{iter}}-\mathbf{X}_{\textit{iter-1 }}\right\|_{F}/\left\|\mathbf{X}_{\textit{iter-1 }}\right\|_{F}$$

(19)

where $\mathbf{X}_{\textit{iter}}$ denotes the reconstructed image at the iterth ADMM iteration. Fig. 11 shows the relative changes of 2D imaging with R = 4, 6, 10.2, and 11.7. The relative change rapidly stabilizes within a few iteration numbers and has no significant alterations when the iteration number reaches 15. Additionally, we changed the CG iteration number from 7 to 20 to test its effects on the reconstruction. The nRMSE variations of the reconstructed images with different CG iteration numbers were less than 1%, indicating the CG iteration number may have little effect on reconstruction after several iterations.

We used several tricks to speed up the reconstruction in our experiments. For the patch extraction step, a patch offset with three was used to reduce the number of searched patches. Also, the patch extraction step was implemented in a graphics processing unit (GPU). The time consumption of this step was reduced from 9.5s to 0.14s for the 2D patch extraction and from 600 mins to 171s for the 3D patch extraction. The memory footprints were 411 MB and 821 MB for the 2D and 3D patch extractions. The $\text{T}_{1\rho}$ estimation for the histogram analysis of the tissue clustering was applied every three iterations to speed up the reconstruction. Additionally, parallel computing was applied in the HOSVD denoising step of tensors $\mathcal{T}$ and $\mathcal{Z}$ for each spatial and parametric tensor. The reconstruction times of the six methods are listed in supplementary information Table S3.

For 3D imaging, the images can be reconstructed slice-by-slice using the SMART method, where the 2D patch extraction is implemented. Compared to the 2D patch extraction, the correlation information between the slices can be utilized to improve the reconstruction. Therefore, in this study, the 3D patch extraction was selected in the SMART method. Fig. 12 shows a slice of the reconstructed images at R = 9.04 using the 2D and 3D patch extraction, respectively. The reference images from the fully sampled data and magnified images are also shown. As shown in the magnified images, the 3D patch extraction displays better image resolution and detail, particularly in the regions indicated by the red frames.

Previous studies [59, 60] have utilized convolutional sparse coding (CSC) to reconstruct images from undersampled data. Thanh et al.[59] proposed a filter-based dictionary learning method using a 3D CSC to recover the high-frequency information of the MRI images. The CSC employed a sparse representation of the entire image computed by the sum of a set of convolutions with dictionary filters instead of dividing an image into overlapped patches. It can solve high computational complexity and long-running time problems due to patch-based approaches and avoid artifacts caused by patch aggregation. However, filters must be carefully chosen, and artifacts may appear in some reconstructed images due to inaccurate filters[60].

Low rankness has been widely used in fast MR parametric mapping. The annihilating filter-based low-rank Hankel matrix method (ALOHA)[30] exploits the transform domain sparsity and the coil sensitivity smoothness, representing a low-rank Hankel structured matrix in the weighted $k$-space domain. The prior information used in ALOHA differs from the image domain redundancy utilized in our study. Specifically, the spatial and temporal redundancies are mainly represented as a summation of exponentials and the Hankel matrix at each voxel with a low-rank property. Multi-scale low-rank with different patch sizes was also used to capture global and local information. It has a more compact signal representation and may further improve MR reconstruction. Frank Ong et al.[61, 62] have applied multi-scale low-rank in face recognition preprocessing and dynamic contrast-enhanced MRI. However, when constructing multiple low-rank tensors, the memory footprint is large, and time consumption is increased several times compared with a single patch size. Another problem is that the multi-scale low-rank decomposition is not translation invariant; shifting the input changes the resulting decomposition. This translation variant nature often creates blocking artifacts near the block boundaries, which can be visually jarring. Introducing overlapping partitions of the patches so that the overall algorithm is translation invariant can remove these artifacts. The proposed SMART method uses a single patch size with overlapped patches considering the limits of memory and reconstruction time.

In this study, the setup of the proposed method is non-convex, and the patch selection operator is time-varying. The convergence of the whole algorithm cannot be completely guaranteed. However, based on our experimental results, the proposed method works well despite the lack of a convergence guarantee. We found that after several iterations, the relative change of reconstructions between two consecutive iterations and the patch extract selection varies slightly. Generally, this estimation based on the patch structure remains consistent along with the increasing iteration. In addition, the block matching algorithm for patch extraction is heuristic, and the current selection strategy may not present the best clusters. In future work, patches obtained through a sliding operation throughout the data may be used to construct the high order tensor. The local and nonlocal redundancies and similarities between the contrast images can be exploited using the tensor dictionary learning method by jointly applying the low-rank constraints on the dictionaries, and the sparse constraints on the core coefficient tensors [63]. This may further improves the reconstruction performance.

The iterative method’s long reconstruction time significantly hinders the proposed method’s clinical application. However, there are two potential solutions to this problem. One solution is using parallel computing via GPU. We have already implemented the patch extraction in GPU. Although the reconstruction time is still too long for clinical use, we plan to improve the speed of the SMART method by leveraging parallel computing technology, such as multi-GPU and cluster systems. The other solution is the deep learning-based reconstruction method. In recent years, the unrolling-based strategy has established a bridge from traditional iterative CS reconstruction to deep neural network design [64]. Previous studies have shown that an unrolling-based deep network significantly reduced the reconstruction time and improved the reconstruction quality through the learnable regularizations and parameters. However, deep learning-based reconstruction requires a considerable amount of high-quality training data, and the network performance may degrade due to the motion artifact during the data acquisition. Training data collection is critical and will be the bottleneck of the deep learning-based reconstruction method. In our future work, we will try both strategies to promote the clinical use of the proposed method.

The 2D dataset was collected on a 3T scanner (TIM TRIO,Siemens, Erlangen, Germany) using a 12-channel head coil from three volunteers (age: 24 ± 2 years). The 3D dataset was collected on a 3T scanner (uMR 790, United Imaging Healthcare, Shanghai, China) using a 32-channel head coil from two healthy volunteers (age: 26 ± 1 years). $\text{T}_{1\rho}$-weighted images were acquired using the fast spin echo (FSE) sequence for 2D imaging and the modulated flip angle technique in refocused imaging with an extended echo train (MATRIX) sequence for 3D imaging, both using a spin-lock pulse at the beginning of the sequence to generate the $\text{T}_{1\rho}$-weighting [20, 53]. The main imaging parameters that used in the 2D and 3D imaging applications are listed in Supplementary information Table S1.

The 2D phantom dataset was also collected a 3T scanner (TIM TRIO,Siemens, Erlangen, Germany) using a 12-channel head coil. The imaging sequence was the same as used in the in vivo data acquisition. To obtain an accurate $\text{T}_{1\rho}$ estimation of the phantom, a longer TR of 4000 ms was used. Therefore, the scan time of the phantom data acquisition was longer than that of the in vivo data acquisition.

The reconstruction parameters that used the SMART method are listed in Supplementary information Table S2. In the L + S method, the singular-value thresholding and soft-thresholding algorithms were used to solve for the L and S components, respectively. The ratio for singular-value thresholding was set as 0.02, whereas that for soft-thresholding was set as [0.02, 0.025, 0.025, 0.035, and 0.035] to achieve optimal performance. In the LLR method, the iteration numbers were set as 40, 50, and 60 for R = 4 and 6, respectively. The block size was set as $8\times 8$, and the ratio of the largest singular value used for the threshold of the singular value decomposition was initially set as 0.03, reduced to 0.01 after 10 iterations, and finally reduced to 0.005 in the final 10 iterations. The parameters for PROST and PROST + HM were the same as those for SMART. In the MORASA, the rank of the Casorati matrix was 2. The threshold for the wavelet coefficients and that of the Hankel matrix were 0.01 and 0.03, respectively.

Supplementary Information Fig.S1(a-c) shows the $\text{T}_{1\rho}$-weighted image of the real phantom at TSL = 40 ms, the obtained $\text{T}_{1\rho}$ map, and the estimated $\text{T}_{1\rho}$ values of the nine vials, respectively.

Supplementary information Fig. S4 shows the corresponding $\text{T}_{1\rho}$ maps reconstructed using the SMART, PROST + HM, PROST, L + S, LLR, and MORASA methods and the error maps between the reconstructed $\text{T}_{1\rho}$ map and the reference map. The $\text{T}_{1\rho}$ maps reconstructed using the SMART method exhibited superior reconstruction performance compared to those reconstructed using the other four methods, with the smallest nRMSEs and error maps for both acceleration factors.

To clearly demonstrate the necessity of updating tissue maps in the iteration, the representative zero-filled reconstruction with R = 6, the intermediate results at the 12th iteration , the parameter histogram, the estimated $\text{T}_{1\rho}$ map, and the corresponding tissue groups and the ranks of tissue groups are included in the Supplemental information Fig. S7. To show the tissue group, we randomly selected 10 of them and chose 10 pixels from each of the group randomly, since the total tissue groups are too much to be shown. Then, an image block with the size of $10\times 5$ can be constructed for each tissue group, where 10 is the pixel number chosen and 5 is the TSL number. The image blocks are shown as an example of the tissue groups. To calculate the rank of a tissue group, we construct the low-rank matrix with all pixels in the tissue group. The size of the matrix is $N_{p}\times 5$ where $N_{p}$ is the pixel number in the group and 5 is the TSL number. The rank value of the matrix was calculated using the SVD method. The rank curves of all tissue groups are shown.

The reconstruction time for each method used in the 2D (R = 4 ) and 3D (R = 6.76) imaging applications is listed in Supplementary information Table S3.

The algorithms for HOSVD and the SMART method are shown in Algorithms 1 and 2, respectively.

The convergence of the proposed method is still an open problem[18, 39]. Although empirical evidence suggests that the reconstruction algorithm has very strong convergence behavior, the satisfying proof is still under study. In this study, we analyze the time complexity of the SMART method using the $\mathcal{O}$ notation, and we give a weak convergence result with the following Theorem :

Theorem 1: Suppose Assumption (A1) in the Proof section holds, and the multipliers are bounded. Then, we have

$$\min_{i\in[N]}\left\|\left(\textbf{X}^{*},\mathcal{T}^{*},\mathcal{Z}^{*},\alpha_{1}^{*},\alpha_{2}^{*}\right)\right\|_{F}^{2}\leq\mathcal{O}\left(\frac{1}{N}\right)$$

(20)

where $\left(\textbf{X}^{*},\mathcal{T}^{*},\mathcal{Z}^{*},\alpha_{1}^{*},\alpha_{2}^{*}\right)$ denotes the subgradient of the augmented Lagrangian function $L^{i}\left(\textbf{X},\mathcal{T},\mathcal{Z},\alpha_{1},\alpha_{2}\right)$ and $N$ is the total number of iterations.

The proof of Theorem 1 is described in the next section.

In this section, we provide proofs on the following assumptions:

(A1) The sequence of paired proximal operators $(\Gamma^{n},\Gamma^{n+1})$ is asymptotically nonexpansive with a sequence { $\epsilon^{n+1}$ }, e.g.,

$$\|\Gamma^{n}(u)-\Gamma^{n+1}(u)\|_{F}^{2}\leq(1+\epsilon^{n+1})\|u-v\|_{F}^{2},\forall{u,v,n}$$

(21)

The assumption (A1) is a standard assumption for analyzing the convergence of the learned iterative algorithms, which also can be found in [65, 66].

Considering the operator $P_{i}$ is performed at every iteration in, let $\Gamma=[V(\mathcal{T}_{1}),\ldots,V(\mathcal{T}_{i}),\ldots,V(\mathcal{Z}_{1}),\ldots,V(\mathcal{Z}_{j}),\ldots]^{T}$, $\alpha=[V(\alpha_{11}),\ldots,V(\alpha_{1i}),\ldots,V(\alpha_{21}),\ldots,V(\alpha_{2j}),\ldots]^{T}$,and $Q=[P^{n}_{1},\ldots,P^{n}_{i},\ldots,H_{1},\ldots,H_{j},\ldots]$, where $V$ denotes an operator which vectorizes the tensor. Then (7) can be expressed as

$$\begin{gathered}L^{n}\left(\textbf{X},\Gamma,\alpha\right):=\frac{1}{2}\left\|E\textbf{X}-\textbf{Y}\right\|_{F}^{2}+\lambda_{1}\sum_{i}\left\|\mathcal{T}_{i}\right\|_{*}+\lambda_{2}\sum_{i}\left\|\mathcal{Z}_{j}\right\|_{*}\\ -\left\langle\alpha,\Gamma-Q^{n}(\textbf{X})\right\rangle+\frac{\mu}{2}\left\|\Gamma-Q^{n}(\mathbf{X})\right\|^{2}_{2}\end{gathered}$$

(22)

According to the $\frac{\mu}{2}$-strongly convex of $L^{n}(\textbf{X}^{n},\Gamma,\alpha^{n})$ with respect to $\Gamma$, we have

$$\begin{gathered}L^{n}(\textbf{X}^{n},\Gamma^{n},\alpha^{n})-L^{n}(\textbf{X}^{n},\Gamma^{n+1},\alpha^{n})\geq\frac{\mu}{2}\|\Gamma^{n}-\Gamma^{n+1}\|_{2}^{2}\end{gathered}$$

(23)

Since $L^{n}(\textbf{X},\Gamma^{n+1},\alpha^{n})$ is also $\frac{\mu}{2}$-strongly convex with respect to $\mathbf{X}$, we have

$$\begin{gathered}L^{n}(\mathbf{X}^{n},\Gamma^{n+1},\alpha^{n})-L^{n}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n})\\ \geq\frac{\mu}{2}\|\mathbf{X}^{n+1}-\mathbf{X}^{n}\|^{2}_{F}\end{gathered}$$

(24)

According to the update rule of the multiplier $\alpha$, we have

$$\begin{gathered}L^{n}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n})-L^{n}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n+1})\\ =-\frac{1}{\mu}\|\alpha^{n+1}-\alpha^{n}\|_{2}^{2}\end{gathered}$$

(25)

On the other hand, we have

	$\displaystyle L^{n}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n+1})-L^{n+1}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n+1})$		(26)
	$\displaystyle\geq-\left\langle\alpha^{n+1},Q^{n+1}(\mathbf{X}^{n+1})-Q^{n}(\mathbf{X}^{n+1})\right\rangle$
	$\displaystyle+\frac{\mu}{2}\|Q^{n+1}(\mathbf{X}^{n+1})-Q^{n}(\mathbf{X}^{n+1})\|_{2}^{2}=0$

where the sencond equality comes from the non-expansion assumption (A1) of $Q^{n}$. With respect to the first-order optimization condition for the update on $\mathbf{X}$, that is, $\nabla_{\mathbf{X}}L^{n}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n})=0$, we have

$$\begin{gathered}0=E^{*}(E\mathbf{X}^{n+1}-\mathbf{Y})+(Q^{n})^{*}\alpha^{n}-\mu(Q^{n})^{*}(Q^{n}(\mathbf{X}^{n+1}-\Gamma^{n+1}))\\ =E^{*}(E\mathbf{X}^{n+1}-\mathbf{Y})+(Q^{n})^{*}\alpha^{n+1}\end{gathered}$$

(27)

where $E^{*}$ and $Q^{*}$ represent the Hermitian adjoint of operators $E$ and $Q$, respectively. Then we have

	$\displaystyle\underline{\sigma}_{Q^{n}}\left\|\alpha^{n+1}-\alpha^{n}\right\|_{2}\leq\|E^{*}E(\mathbf{X}^{n+1}-\mathbf{X}^{n})\|_{F}$		(28)
	$\displaystyle\leq\bar{\sigma}_{E^{*}E}\left\|\mathbf{X}^{n+1}-\mathbf{X}^{n}\right\|_{F}$

where $\underline{\sigma}_{Q^{n}}$ denotes the minimum singular value of $Q^{n}$, and $\bar{\sigma}_{E^{*}E}$ denotes the maximum singular value of $E^{*}E$.

Combining (23), (24), (25), (26), (28), we have

	$\displaystyle L^{n}(\mathbf{X}^{n},\Gamma^{n},\alpha^{n})-L^{n+1}(\mathbf{X}^{n+1},\Gamma^{n+1},\alpha^{n+1})$		(29)
	$\displaystyle\geq(\frac{\mu}{2}-\frac{\bar{\sigma}^{2}_{E^{*}E}}{\underline{\sigma}^{2}_{Q^{n}}})\|\mathbf{X}^{n+1}-\mathbf{X}^{n}\|^{2}_{F}+\frac{\mu}{2}\|\Gamma^{n+1}-\Gamma^{n}\|^{2}_{2}$

Summing both sides of the above inequality from 0 to $N$, we have

		$\displaystyle\frac{N}{2}\left[\left(\frac{\mu}{2}-\frac{\bar{\sigma}_{E^{*}E}^{2}}{\underline{\sigma}_{Q^{n}}^{2}}\right)\left\|\mathbf{X}^{n+1}-\mathbf{X}^{n}\right\|_{F}^{2}+\frac{\mu}{2}\left\|\Gamma^{n+1}-\Gamma^{n}\right\|_{2}^{2}\right]$		(30)
		$\displaystyle\leq L^{0}\left(\mathbf{X}^{0},\Gamma^{0},\alpha^{0}\right)-L^{N}\left(\mathbf{X}^{N},\Gamma^{N},\alpha^{N}\right)<+\infty$

On the other hand, let $\left(\mathbf{X}_{n+1}^{*},{\Gamma}_{n+1}^{*},{\alpha}_{n+1}^{*}\right)\in\partial L\left(\mathbf{X}^{n+1},{\Gamma}^{n+1},{\alpha}^{n+1}\right)$, we have

$$\left\{\begin{array}[]{l}{\Gamma}_{n+1}^{*}\in{\alpha}^{n+1}-\mu\left(Q^{n+1}\mathbf{X}^{n+1}-{\Gamma}^{n+1}\right)+\lambda\sum_{p}\partial\left\|\Gamma_{p}^{n+1}\right\|_{*}\\ \mathbf{X}_{n+1}^{*}=E^{*}\left(E\mathbf{X}^{n+1}-\mathbf{Y}\right)-Q^{n+1}{\alpha}^{n+1}\\ +\mu{(Q^{n+1})}^{*}\left(Q^{n+1}\mathbf{X}^{n+1}-{\Gamma}^{n+1}\right)\\ {\alpha}_{n+1}^{*}={\Gamma}^{n+1}-Q^{n+1}\mathbf{X}^{n+1}\end{array}\right.$$

(31)

where $p$ is equal to the sum of $i$ and $j$. According to the iterative rule of ADMM algorithm, we have

$\displaystyle\left\{\begin{array}[]{l}0\in{\alpha}^{n}-\mu\left(Q^{n}\mathbf{X}^{n}-\Gamma^{n+1}\right)+\lambda\sum_{p}\partial\left\|\Gamma_{p}^{n+1}\right\|_{*}\\ 0=E^{*}\left(E\mathbf{X}^{n+1}-\mathbf{Y}\right)-Q^{n}\alpha^{n}\\ +\mu(Q^{n})^{*}\left(Q^{n+1}\mathbf{X}^{n+1}-{\Gamma}^{n+1}\right)\\ 0=\alpha_{n+1}-\alpha_{n}+\mu({\Gamma}^{n+1}-Q^{n}\mathbf{X}^{n+1})\end{array}\right.$

(32)

Combing (31) and (32), we have

$$\left\{\begin{array}[]{l}{\Gamma}_{n+1}^{*}={\alpha}^{n}-{\alpha}^{n+1}-\mu\left(Q^{n+1}\boldsymbol{X}^{n+1}-Q^{n}\boldsymbol{X}^{n}\right)\\ \boldsymbol{X}_{n+1}^{*}=Q^{n}{\alpha}^{n}-Q^{n+1}{\alpha}^{n+1}\\ {\alpha}_{n+1}^{*}=\frac{1}{\mu}\left({\alpha}^{n}-{\alpha}^{n+1}\right)\end{array}\right.$$

(33)

Then, we have

		$\displaystyle\left\|\left({\Gamma}_{n+1}^{*},\boldsymbol{X}_{n+1}^{*},{\alpha}_{n+1}^{*}\right)\right\|_{F}^{2}$		(34)
	$\displaystyle\leq$	$\displaystyle\left(2+(1+\epsilon^{n})+\frac{1}{\mu^{2}}\right)\left\|{\alpha}^{n}-{\alpha}^{n+1}\right\|_{2}^{2}$
		$\displaystyle+2\mu^{2}(1+\epsilon^{n})\left\|\boldsymbol{X}^{n+1}-\boldsymbol{X}^{n}\right\|_{F}^{2}$
	$\displaystyle\leq$	$\displaystyle\left[\frac{\bar{\sigma}_{E^{*}E}^{2}}{\underline{\sigma}_{{Q}^{n}}^{2}}\left(3+\epsilon^{n}+\frac{1}{\mu^{2}}\right)+2\mu^{2}(1+\epsilon^{n})\right]\left\|\boldsymbol{X}^{n+1}-\boldsymbol{X}^{n}\right\|_{F}^{2}$

where ${\bar{\sigma}}_{Q^{n}}$ denotes the largest singular values of $Q^{n}$, the first inequality is due to the nonexpensive assumption of $Q^{n}$, and the second inequality is due to the relation (28). By the convexity of norm, we have

		$\displaystyle\min_{i\in[N]}\left\|\left({\Gamma}_{i}^{*},\boldsymbol{X}_{i}^{*},{\alpha}_{i}^{*}\right)\right\|_{F}^{2}$		(35)
	$\displaystyle\leq$	$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\left\|\left({\Gamma}_{n}^{*},X_{n}^{*},{\alpha}_{n}^{*}\right)\right\|_{F}^{2}$
	$\displaystyle\leq$	$\displaystyle\left[\frac{\bar{\sigma}_{E^{*}E}^{2}}{\underline{\sigma}_{{Q}^{n}}^{2}}\left(3+\epsilon^{n}+\frac{1}{\mu^{2}}\right)+2\mu^{2}(1+\epsilon^{n})\right]$
		$\displaystyle\cdot\frac{1}{N}\sum_{n=1}^{N}\left\|\boldsymbol{X}^{n+1}-\boldsymbol{X}^{n}\right\|_{F}^{2}$
	$\displaystyle\leq$	$\displaystyle\left[\frac{\bar{\sigma}_{E^{*}E}^{2}}{\underline{\sigma}_{{Q}^{n}}^{2}}\left(3+\epsilon^{n}+\frac{1}{\mu^{2}}\right)+2\mu^{2}(1+\epsilon^{n})\right]\left(\frac{\mu}{2}-\frac{\bar{\sigma}_{E^{*}E}^{2}}{\mu\underline{\sigma}_{{Q^{n}}}^{2}}\right)^{-1}$
		$\displaystyle\cdot\frac{L\left(\boldsymbol{X}^{0},{\Gamma}^{0},{\alpha}^{0}\right)-L\left(\boldsymbol{X}^{N},{\Gamma}^{N},{\alpha}^{N}\right)}{N}$

From aforementioned analysis, we know that we can choose $\mu=\mathcal{O}(\frac{\bar{\sigma}_{E^{*}E}}{\underline{\sigma}_{Q}})$ which implies that

$$\min_{i\in[N]}\left\|\left({\Gamma}_{i}^{*},\boldsymbol{X}_{i}^{*},{\alpha}_{i}^{*}\right)\right\|_{F}^{2}\leq\mathcal{O}\left(\frac{1}{N}\right)$$

(36)