An MTL1TV non-convex regularization model for MR Image reconstruction using the alternating direction method of multipliers

1.   Introduction

Magnetic resonance imaging (MRI) is an approach broadly used for clinical diagnosis because of its ability to effectively reflect changes in human organs ^[1,2,3]. However, a long-standing challenge is that the acquisition time of MRI is relatively long. To address this problem, researchers have proposed numerous techniques, such as sparse sampling ^[4,5] and parallel imaging ^[6,7]. After the theory of compressed sensing (CS) ^[8] was proposed, many scholars applied it to recover magnetic resonance (MR) images from highly undersampled measurements ^[9,10]. Research on MR image compression sensing has thus become increasingly important.

In the MRI domain, the data acquisition of MRI can be modeled as follows:

where ${\rm x} \in {\mathbb R} {^n}$ denotes the desired MRI data, ${\rm y} \in {\mathbb R} {^m}$ is the observed undersampled k-space MRI data, and the matrix ${\rm A} = {\rm RF} \in {\mathbb R}{^{m \times n}}(m < n)$. ${\rm R}$ represents the undersampling operator, such as the radial sampling operator and Cartesian sampling operator. ${\rm F}$ is the discrete Fourier transform (DFT) operator, and $\varepsilon \in {\mathbb R}{^m}$ is the additive Gaussian noise. The objective of MR image reconstruction is to recover ${\rm x}$ from ${\rm y}$. Consequently, the reconstruction problem of (1) can be formulated as

in which $\left\| {{\rm y} - {\rm Ax}} \right\|_2^2$ is the data fidelity term, $\lambda > 0$ is the regularization parameter, and $\varphi (x)$ is some regularizing functional exploiting image prior knowledge. In CS-MRI, $\varphi ({\rm x})$ is the standard total variation (TV) norm (or ${l_1}$-based regularization) ^[11], and the corresponding reconstruction model can be expressed as

where ${\left\| {\rm x} \right\|_{TV}} = {\left\| {D{\rm x}} \right\|_1}$ ($D$ is the gradient operator). Although standard TV regularization is convex and widely used, the ${l_1}$-norm regularized model induces bias for large coefficients ^[12] and, hence, lacks the oracle property. To overcome the shortcomings of ${l_1}$-norm regularization, researchers have proposed some non-convex regularization methods based on the capped-L1 (CaP) ^[13], smoothly clipped absolute deviation (SCAD) ^[14], minimax concave (MC) penalty ^{[15,16,17,18]}, arctangent penalty (ATAN) ^[19], logarithm function (Log) ^[20], transformed ${l_1}$ (TL1) ^[21,22], fraction function penalty ^[23], and non-convex graph total variation (GTV) regularization ^[23]. For instance, Liu et al. ^[24] introduced the minimax concave total variation (MCTV) penalty as a regularization term for MR brain image reconstruction. Further, Luo et al. ^[25] used the arctangent function as the non-convex TV regularize term (AtanTV) for MR image reconstruction. The same team ^[26] introduced the non-convex MR image reconstruction model via the SCAD penalty function. Lu et al. ^[27] also used the non-convex Cauchy total variation (CauchyTV) as the regularization term for MR image reconstruction. These studies have shown that non-convex penalties in image restoration usually have better performance than the ${l_1}$-norm (or standard TV).

The contributions of this study are as follows. First, we slightly modify the transformed ${l_1}$ (TL1) function and propose a modified TL1 penalty function (MTL1). Second, inspired by the previous work, we use the MTL1 as the regularization term (MTL1TV) to construct a non-convex regularization model for MR image reconstruction. Third, the proposed model can be solved by the alternating direction method of multipliers (ADMM) ^[28,29,30]. Finally, numerical experiments on several MR images showed that compared with the traditional TV and state of the art models, the performance of the MTL1TV model was significantly improved. These results confirm the effectiveness of the proposed model.

The remainder of this paper is organized as follows. In Section 2, we formally introduce the MTL1 function and present a non-convex MTL1TV model, which is based on the MTL1 function. In Section 3, we utilize the ADMM algorithm to address the proposed non-convex model. The convergence result is provided in Section 4. The experimental results are then presented in Section 5. Finally, we conclude the paper in Section 6.

2.   Preliminaries and MTL1TV MRI model

In this section, we provide the definition of the the modified TL1 penalty function (MTL1) and give an explicit expression of the proximity operator of MTL1. Then, we use MTL1 to define MTL1TV regularization and structure a non-convex MTL1TV model for MR image reconstruction.

2.1.   Preliminaries

Definition 1. The modified transformed ${l_1}$ function (MTL1) is defined as

where parameter $a > 0$. It is obtained by slightly modifying the following TL1 penalty function ^[20]:

Figure 1 shows that the MTL1 is a good alternative to the ${l_1}$ norm. The larger the value of parameter $a$ is, the closer the behavior of MTL1 is to the ${l_1}$ norm.

The proximal operator of ${\phi _a}(x)$ at $t \in {\mathbb R}$ is defined as

Similar to Theorem Ⅲ.1 of ^[20], we can derive that for the MTL1 function, there exists a closed-formed representation of the optimal solution to (2.3).

Theorem 1. The optimal solution ${x^ * } = \arg \min {\theta _{\lambda, t}}(x)$ is a threshold function defined as

where

with $\varphi (t) = \arccos (1 - \frac{{27\lambda {a^2}}}{{2{{(a + \left| t \right|)}^3}}}),$ and the threshold value ${\delta _\lambda }$ satisfies

To prove Theorem 1, we first give the following lemmas.

Lemma 1. Define three parameters

for any parameters $\lambda > 0$ and $a > 0$, then the inequalities $t_1^ * \le t_3^ * \le t_2^ *$ hold. Furthermore, they are equal to $\frac{a}{{2}}$ when $\lambda = \frac{a}{{{2}}}$.

Lemma 2. For any given $t$, the roots of the following two cubic polynomials of $x$ satisfy the following properties:

1) If $t > t_1^ *$, then the cubic polynomial

has three distinct real roots, and the largest root ${x_0}$ is given by ${x_0} = {g_\lambda }(t)$, $\left| {{g_\lambda }(t)} \right| \le \left| t \right|$.

2) If $t < - t_1^ *$, then the cubic polynomial

has also three distinct real roots, and the smallest root ${x_0}$ is given by ${x_0} = {g_\lambda }(t)$.

The proof of Lemma 2 is essentially the same as Lemma Ⅲ.1 in ^[20] and Lemma $8$ in ^[22], with only minor changes required due to the different values of $t_{1}^{*},$ $t_{2}^{*},$ and $t_{3}^{*}$ and different cubic polynomials. The proof of Theorem 1 is similar to Theorem Ⅲ.1 in ^[20] and Lemma $9$ in ^[22]. The detailed proofs of Lemma 2 and Theorem 1 can be found in Appendixes A and B, respectively. Based on the Theorem 1, the following corollary holds.

Corollary 1. For any $a \ge 2\lambda$, the function $\theta_{\lambda, t}(x)$ defined in (2.6) is strictly convex.

The proof of Corollary 1 can be found in Appendix C.

The graph presented in Figure 2(a) shows the plots of the function $\theta_{\lambda, t}(x) = \lambda {\phi _a}(x) + \frac{1}{2}{(x - t)^2}$ with $t = 0$, $\lambda = 5$, and $a = 15, 50$. It can be seen in Figure 2(a) that the function $\theta_{\lambda, t}(x)$ is strictly convex for $a \ge 2\lambda$. However, when $\lambda = 5$, and $a = 0.5, 1$, the graph presented in Figure 2(b) shows that the function $\theta_{\lambda, t}(x)$ is non-convex.

From Theorem 1, the convexity of the function $\theta_{\lambda, t}(x)$ can be ensured by constraining the parameter $a$ as $a \ge 2\lambda$, which means that there exists the unique global minimizer to $\theta_{\lambda, t}(x)$ when $a \ge 2\lambda$. The unique global minimizer to $\theta_{\lambda, t}(x)$ can be expressed in the following theorem.

According to (2.1), we consider another function ${\psi _a}(x)$, which is induced by ${\phi _a}(x)$:

Figure 3 shows the curves of the functions $\left| x \right|$, ${\phi _a}(x)$, ${\psi _a}(x)$ and $\frac{x^{2}}{a}(a = 10)$. Clearly, the function ${\psi _a}(x)$ is a convex, continuously differentiable function that satisfies $0 \le {\psi _a}(x) \le \left| x \right|$ and $\psi_{a}(x)\leq \frac{x^{2}}{a}$ for all $x \in {\mathbb R}$.

Then, the MTL1 function can be written as

To establish our MTL1TV MRI model, we define a multivariate penalty ${\Phi _a}:{\mathbb R} {^n}\to {\mathbb R}$ based on the above penalty ${\phi _a}$. For this purpose, we first define a multivariate generalization form of the corresponding function ${\psi _a}$ in (2.8).

Definition 2. Let ${\psi _a}$ be correspondingly defined in (2.8); then ${\Psi _a}:{\mathbb R} {^n}\to {\mathbb R}$ is defined as

Remark 1. The function ${\Psi _a}$ satisfies $0 \le {\Psi _a}({\rm v}) \le {\left\| {\rm v} \right\|_1}$ and ${\Psi _a}({\rm v}) \le \frac{1}{a}\left\| {\rm v} \right\|_2^2$ for all ${\rm v} \in {\mathbb R}{^n}$.

Then, the multivariate penalty of ${\phi _a}$ is given by

Finally, by replacing ${\rm v}$ with gradient ${\rm Dx}$ in (2.11) (${\rm D}$ is the first-order difference operator/matrix), we obtain the definition of MTL1TV below.

Definition 3 (MTL1TV). The MTL1TV regularizer: ${\left\| x \right\|_{MTL1TV}}:{\mathbb R} {^n}\to {\mathbb R}$,

Now, the MR image reconstruction model can be formulated as

where $\lambda > 0$ is called the regularization parameter. Because the MTL1 penalty function is non-convex, the model in (2.13) is non-convex. However, the global convexity of the objective function can be controlled by adjusting parameter $a$, as shown in Theorem 2 below.

Theorem 2. Let $\lambda > 0$, $a > 0$. Define $F:{\mathbb R} {^n}\to {\mathbb R}$ as

If $a \ge 2\lambda$ and ${\rm D}^{\rm T}{\rm D} \preceq {{\rm A}^{\rm T}}{\rm A}$, then ${F_a}({\rm x})$ is a convex function.

Proof. Substituting (2.12) into (2.14), we have

Note that $\frac{1}{2}\left\| {\rm y} \right\|_2^2$ does not depend on ${\rm x}$ and that ${\rm y}^{\rm T}{\rm Ax}$ is linear in ${\rm x}$. Hence, the function ${F_a}({\rm x})$ is convex if ${F_2}({\rm x})$ is convex, where ${F_2}({\rm x})$ is defined as

The first term is convex if ${{\rm A}^{\rm T}}{\rm A} \succeq \frac{{2\lambda }}{a}{{\rm D}^{\rm T}}{\rm D}$. Since $a \ge 2\lambda$ and are ${\rm D}^{\rm T}{\rm D} \preceq {{\rm A}^{\rm T}}{\rm A}$ given, it follows that the first term is convex. From Remark 1, the last term is convex. Hence, ${F_2}({\rm x})$ is convex. Therefore, ${F_a}({\rm x})$ is convex.

3.   ADMM for MTL1TV model

In this section, the optimization algorithm is presented in detail. From Theorem 2, (2.13) is a convex optimization problem under certain conditions, which can be solved via the convex optimization algorithms. Hence, the ADMM is used to effectively solve the model in (2.13). First, according to the definition of ${\left\| {\rm x} \right\|_{MTL1TV}}$, we rewrite (2.13) as follows:

Next, we introduce the auxiliary variable ${\rm z}$ with the constraint ${\rm z} = {\rm Dx}$. Then, the optimization problem (3.1) is rewritten as

Under the ADMM framework, the augmented Lagrangian function of (3.2) is

where ${\rm w}\in {\mathbb R} {^n}$ is the Lagrange multiplier, and $\beta > 0$ is a penalty parameter. Then, we invoke the ADMM by iterating the variable updates in (3.5) to (3.7):

Now, we give the detailed steps for solving the subproblems in (3.5) and (3.6) alternatively.

Step 1. Update ${{\rm{x}}^{k + 1}}$ with ${\rm{z}^k}$, and keep ${\rm{w}^k}$ fixed. Let the gradient of this objective function be zero, and obtain the following equation:

The matrix $\rm {A = RF}$, where ${\rm F}$ is the Fourier operator with the property $\rm{F^T} = {\rm{F}^{ - 1}}$. Since ${\rm {D^T}D}$ is a cyclic matrix, it can be diagonalized by Fourier transform. Therefore, ${{\rm{x}}^{k + 1}}$ can be solved by using two Fourier transforms:

Step 2. Update ${{\rm{z}}^{k + 1}}$ with ${{\rm{x}}^{k + 1}}$, and keep ${\rm{w}^k}$ fixed. From (3.6), we have

According to Theorem 1, ${{\rm{z}}^{k + 1}}$ is given by the following expression:

Based on the above analysis, the specific algorithm for solving (3.1) can be summarized as Algorithm 1.

4.   Convergence analysis

In this section, we discuss the convergence of the proposed alternating minimization algorithm. Through this paper, we assume that $ker(A)\cap ker(D) = {0}$. The assumption is very reasonable for the imaging problem ^[31,32]. In Algorithm 1, the dominant computation is the steps to solve the two minimization subproblems (3.5) and (3.6). Inspired by ^[31,32], we have the following convergence result for Algorithm 1.

Assume that $({\rm x}^{*}, {\rm z}^{*}, \lambda^{*})$ is a stationary point that satisfies the first order optimality conditions of (3.6):

Obviously, one can easily find that $\rm x^{*}$ satisfies

which is the first order necessary condition of (2.13) about stationary point. Otherwise, by Algorithm 1, each iteration step about the subproblems (3.5)–(3.7) follows:

In the following, we will prove the sequence $\{{\rm x}^{k}, {\rm z}^{k}, \lambda^{k}\}$ generated by the Algorithm 1 has a limit point $({\rm x}^{*}, {\rm z}^{*}, \lambda^{*})$ that satisfies (4.1). The following two lemmas state that the metric $L({\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}, {\rm \beta}^{k})$ and the sequences ${{\rm x}^{k}}$ and $\{{\rm z}^{k}\}$ are bounded.

Lemma 3. Let $L({\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}, {\rm \beta}^{k})$ be the sequence generated by Algorithm 1, then $\{{\rm w}^{k}\}$ and $\{L({\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}, {\rm \beta}^{k})\}$ are bounded.

Proof. First, we prove the Lagrange multiplier ${\rm w}^{k}\in \mathbb{R}^{n}$ is bounded and $\|{\rm w}^{k}\|\leq \sqrt{n}\lambda$. Let $d\in\partial\lambda\phi_{a}(x)$, then, if $x > 0$, $d = \frac{\lambda a^{2}}{(a+x)^{2}}\in (0, \lambda)$. If $x < 0$, $d = \frac{-\lambda a^{2}}{(a-x)^{2}}\in (-\lambda, 0)$. If $x = 0$, $d\in [\lambda, \lambda]$. Therefore, $|d|\leq\lambda.$ On the other hand, combine (4.4) and (4.5), and we have ${\rm w}^{k+1}\in \partial\lambda\Phi_{a}({\rm z}^{k+1}).$ Thus, by $|d|\leq\lambda$, we have $\|{\rm w}^{k+1}\|_{\infty}\leq\lambda.$ Because ${\rm w}^{k+1}$ is finite dimensional, it implies that $\{{\rm w}^{k}\}$ is bounded and $\|{\rm w}^{k}\|\leq \sqrt{n}\lambda < \infty.$ Then, we prove $\{L({\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}, {\rm \beta}^{k})\}$ is bounded. By (4.5), we notice that

From the above formula and (3.7), we have

By the definitions of ${\rm z}^{k+1}$ and ${\rm x}^{k+1}$, we know that

Summing (4.7)–(4.10), we have

Summing up the above inequality from $k = 0$, we obtain

Because $\theta > 1$, let $k\rightarrow \infty$, and we have

On the other hand, from $\lim\limits_{k \to \infty}\beta^{k} = \infty,$ and the fact that $\{{\rm w}^{k+1}\}$ is bounded, we know

and then

By (4.13) and (4.15), we know

So, $\{L({\rm x}^{k+1}, {\rm z}^{k+1}, {\rm w}^{k+1}, {\rm \beta}^{k+1})\}$ is bounded.

Lemma 4. Let $\{({\rm x}^{k}, {\rm z}^{k})\}$ be the sequence generated by Algorithm 1, then the sequences $\{{\rm x}^{k}\}$ and $\{{\rm z}^{k}\}$ are bounded.

Proof. From (4.9) and (4.10), we obtain that $L({\rm x}^{k+1}, {\rm z}^{k+1}, {\rm w}^{k}, {\rm \beta}^{k})$ is upper bounded. Next, we know that $L({\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}, {\rm \beta}^{k})$ is strongly convex about the variable $x$, and the following inequality holds:

Summing (4.7)–(4.9) and (4.17), we have

Summing the above inequality from $k = 0$, we have

which together with (4.16) yields

Otherwise, if $\lim\limits_{k \to \infty}\|{\rm x}^{k+1}-{\rm x}^{k}\|^{2}_{2}\neq0$, then $\sum\limits_{0\leq i\leq k}\frac{c^{i}}{2}\|{\rm x}^{i+1}-{\rm x}^{i}\|^{2}_{2} = +\infty$; thus,

and (4.16), a contradiction. By (4.20), we know that $\{x^{k}\}$ is a Cauchy sequence and, thus, is convergent and bounded.

Furthermore, we show that $\{{\rm z}^k\}$ is bounded. By (3.4) and (4.12), we have

Because $\lim\limits_{k\rightarrow \infty}{\beta^{k} = \infty}$ and $\{{\rm w}^{k}\}$ are bounded, we also know that $\{{\rm z}^{k}\}$ and $\{{\rm Ax}^{k}\}$ are bounded.

Now, we are ready for proving the following convergence result.

Theorem 3. Let $\{{\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}\}$ be the sequence generated by Algorithm 1, then, there exists a subsequence $\{{\rm x}^{k_{j}}, {\rm z}^{k_{j}}, {\rm w}^{k_{j}}\}$, which converges to a stationary point $({\rm x}^{*}, {\rm z}^{*}, {\rm w}^{*})$ and satisfies (4.1).

Proof. By Lemmas 3 and 4, the sequence $\{{\rm x}^{k}, {\rm z}^{k}, {\rm w}^{k}\}$ is bounded, so there exists a subsequence $\{{\rm x}^{k_{j}}, {\rm z}^{k_{j}}, {\rm w}^{k_{j}}\}$ that converges to a cluster point $({\rm x}^{*}, {\rm z}^{*}, {\rm w}^{*}).$ From the lower semi-continuity of $L({\rm x, z, w}, \beta)$,

and as ${\rm z}^{k_{j}+1}$ is a minimizer of function $L({\rm x}^{k_{j}+1}, {\rm z}^{k_{j}+1}, {\rm w}^{k_{j}}, \beta^{k_{j}})$ with respect to ${\rm z}$,

By the above two inequalities, we get

One can immediately verify that

Therefore, $({\rm x}^{*}, {\rm z}^{*}, {\rm w}^{*})$ is a Karush-Kuhn-Tucker point of the problem (3.2).

Remark 2. From Corollary 1, if $a\geq2\lambda$, the two subproblems (3.5) and (3.6) are continuous, coercive, and strictly convex. Consequently, every subproblem has a global solution. Additionally, from Theorem 2, we know that the objective function is a convex function, which implies the stationary point is the global minimum point.

5.   Experiments

In this section, we present experiments conducted to verify the performance of the MTL1TV method in MATLAB 2019a on a laptop with 32 GB RAM and an Intel Core i7-12700H CPU at 2.70 GHz. We compared our proposed method with the standard TV ^[11], MCTV ^[24], and transformed total variation (TTV) ^[33] methods. For fair performance comparison, the parameters used in comparing algorithms have been appropriately tuned to yield a better reconstruction effect. We used three metrics to assess the quality and accuracy of image reconstruction: the relative error (RE), peak signal-to-noise ratio (PSNR), and structural similarity index (SSIM). These metrics are defined as follows:

where $\rm{x}$, ${\rm{x}}^{k}$ are the original image and reconstructed image, respectively, and

where $\mu_{\rm{x}}$ and $\mu_{{\rm{x}}^{k}}$ are the mean of $\rm{x}$ and ${\rm{x}}^{k}$, $k$ is the number of iterations in the algorithm, $\sigma^{2}$ and $\sigma_{{\rm{x}}^{k}}^{2}$ denote the variance of $\rm{x}$ and ${\rm{x}}^{k}$, and $\sigma_{{\rm{x}}{\rm{x}}^{k}}$ denotes the covariance of $\rm{x}$ and ${\rm{x}}^{k}$. The positive constants $C_{1}$ and $C_{2}$ are stabilization constants. Generally, the MR image reconstruction is better if it has a lower RE and a higher PSNR and SSIM.

The sampling templates and the experimental data are shown in Figure 4. The test data (a)–(c) was three typical MR images: Shepp-Logan, Brain, and Brain 2; (c) and (d) are two different brain MR images. The sampling templates (d)–(f) were radial, random, and Cartesian sampling, respectively. The size of sampling masks and MRI data was $256 \times 256$. In the experiments, the stopping criterion was usually as follows:

or the maximum number of iteration $200$.

5.1.   Experiments on noiseless data

The sampling masks and MR data shown in Figure 4 are used to evaluate our proposed MTL1TV method. Table 1 displays the PSNR, SSIM, RE, and CPU time required to reconstruct various MR images using different methods with varying sampling templates. As observed in Table 1, the time required to reconstruct an MR image using TV is shorter than other methods, but its performance is inferior. Furthermore, when compared with TV, MCTV, and TTV, it is evident that MTL1TV exhibits higher PSNR and SSIM values with lower RE. This suggests that MTL1TV offers improved MR image reconstruction quality, simultaneously maintaining a reasonable runtime as shown in Table 1.

The visual comparison of reconstruction results and the error images are shown in Figures 5−7. In Figure 5, the Shepp-Logan phantom is used to demonstrate the performance of the proposed method. A Cartesian sampling at a sampling rate of $34\%$ is employed to compare with the four reconstruction models proposed above. The regularization parameters $\lambda$ for TV, MCTV, TTV, and MTL1TV are set to $0.0001$, $0.005$, $0.001$, and $0.005$, respectively. For TTV and MTL1TV, the parameter $a$ is set to $1$ and $0.05$, respectively. From the Figure 5, it is observed that compared with other methods, MTL1TV has the best reconstruction capability, and the reconstructed image is most similar to the original image in the visual effects, which is reflected by the comparison of the difference images.

Figure 6 displays the reconstruction results of brain data under radial sampling with 10 trajectory lines and a $3\%$ sampling ratio. We set the regularization parameters $\lambda$ of TV, MCTV, TTV, and MTL1TV to $0.000001$, $0.000005$, $0.001$, and $0.001$, respectively. The value of parameter $a$ is set as $0.5$ for TTV, while it is set as $0.1$ for MTL1TV. Visually, there is no significant difference. The residual images illustrate that the proposed method can visually gain better reconstructions.

For the Brain 2 image, random sampling with a $30\%$ sampling rate and a $0.1$ sampling radius was employed to test the image. The regularization parameters $\lambda$ were $0.00001$, $0.00001$, $0.0001$, and $0.0005$. The value of $a$ in TTV was $a = 0.5$, while in MTL1TV, $a = 0.1$. The reconstructed images and errors are shown in Figure 7. It shows that the reconstruction results of the four methods are similar, and MTL1TV is slightly better than the other three methods.

Additionally, the convergence of the proposed method is empirically tested. In Figure 8, we illustrate the behavior of the proposed MTL1TV model in the iteration process. It respectively shows the curves of the RE, PSNR, and SSIM values on Shepp-Logan phantom data under a radial mask with a $3\%$ sampling ratio versus iteration numbers. It is shown that as the number of iterations increased, the RE value generated by MTL1TV gradually decreased, while the PSNR and SSIM values increased, which exhibits the stable convergence trend in the subsequent iterations.

5.2.   The Parameters selection

In the proposed model, there are two important parameters: $\lambda$ and $a$. The $\lambda$ is the regularization parameter, which trades the sparsity and data consistency. The parameter $a$ has an important impact on the properties of regularization term and it also affects the quality of reconstructed images. Therefore, it is crucial to choose appropriate values of $a$ and $\lambda$, which can help us obtain better reconstruction results. The authors in ^[22] claimed that $a = 1$ is the best among the tests on signals, and the author in ^[33] empirically selects $a = 5$, which achieves better results in image reconstruction tasks. To demonstrate how to select parameters $a$ and $\lambda$, we present the reconstruction results of the Shepp-Logan phantom ($256\times256$) data using a Cartesian mask with a sampling rate of $34\%$ as an example.

To begin, we conduct the experiments with the fixed parameter $a = 1$ and $a = 5$, while allowing $\lambda$ to vary. Subsequently, the parameter $\lambda$ in TV, MCTV, and TTV can also be selected accordingly. The quantitative comparison results are listed in Table 2. From Table 2, we can see that when $\lambda = 0.005$, the performance of the proposed model is the "best". Not only that, when $\lambda = 0.0001$ in TV, $\lambda = 0.005$ in MCTV, $\lambda_{1} = 0.001$ in TTV (when $a = 1$), and $\lambda_{2} = 0.0005$ in TTV (when $a = 5$), and the reconstructed metric PSNR achieves the highest value. Therefore, we set the parameter $\lambda = 0.0001$ in TV, and $\lambda = 0.005$ in MCTV.

Next, we conduct the experiment with the fixed parameter $\lambda_{1}$ and $\lambda_{2}$, respectively, to asses the impact of different values of parameter $a$ on reconstruction performance. From Table 3, it can be clearly seen that when the $a = 1$, TTV exhibits excellent performance. Similarly, when $a = 0.05$, MTL1TV achieves the best results. Furthermore, according to Table 3, it can be observed that when $a = 0.05 > 2\lambda_{2} = 0.01$, the convexity of the MTL1TV model is clearly indicated. In Figure 9, it can be seen that as the parameter $a$ assumes different values, the differences between the reconstructed MR image and the original image also vary. Notably, when $a = 0.05$, the error images are the smallest.

5.3.   Experiments on noisy data

To demonstrate the performance of the proposed model in noisy condition, we conduct experiments on sampled data contaminated with additive Gaussian noise. In these experiments, Gaussian noise with standard deviation ($\sigma = 0.02$) was added to the real and imaginary parts of sampled data, respectively. Figure 10 presents that the reconstructed results, along with the corresponding errors, for Shepp-Logan data under radial sampling ($10$ trajectory lines and a $3\%$ sampling ratio). When compared to the reconstruction errors, it is evident that our proposed method effectively suppresses noise and artifacts. Table 4 presents the quantitative results of the three images (in Figure 4) with additive noise under radial sampling ($10$ trajectory lines and a $3\%$ sampling ratio). Even though all models were given identical parameter settings in both noisy and noise-free environments, it was observed that the values of PSNR, SSIM, and RE for three images decreased under noisy conditions, as shown in Table 4. Nevertheless, they still produced results that were comparable to their noise-free counterparts. When compared to other methods, our method consistently delivered superior results, highlighting its robustness.

6.   Conclusions

In this paper, we introduced a non-convex TV regularization model for undersampled MRI reconstruction. In the new model, based on the MTL1 function, the traditional TV is replaced with the MTL1TV norm, thus improving the fitting performance of the model. We used the ADMM algorithm to compute the minimization problem. The experimental results demonstrated the effectiveness and efficiency of the MTL1TV method. In our future research, we seek to delve into the utilization of the proposed methodology in dynamic MRI, tomographic imaging, and electrical impedance tomography, by amalgamating deep learning frameworks with variational models.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by Natural Science Foundation of Hubei Province (2023AFD013), Philosophy and Social Sciences of Educational Commission of Hubei Province of China (22Y109), and Foundation of Hubei Normal University (2022055).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   Proof of Lemma 2

Proof. 1) Define the new variable $\eta = a + x$ and substitute it into polynomial (2.6), then it becomes

whose discriminant is:

Due to $t > t_1^ *$ and $\Delta > 0$, the cubic equation has three distinct real roots. We change variables as $\eta = y + \frac{a}{3} + \frac{t}{3} = x + a$. The relation between $x$ and $y$ is: $x = y - \frac{{2a}}{3} + \frac{t}{3}$. In terms of $t$, the cubic polynomial is turned into a depressed cubic as:

The three roots in trigonometric form are:

where $\varphi = \arccos \left({1 - \frac{{27\lambda {a^2}}}{{2{{(a + t)}^3}}}} \right)$. It is obvious that ${y_0} > {y_1} > {y_2}$ and ${y_2} < 0$. By the relation $x = y - \frac{{2a}}{3} + \frac{t}{3}$, the three roots in variable $x$ are ${x_i} = {y_i} - \frac{{2a}}{3} + \frac{t}{3}$. From these formulas, we can prove ${x_0} > {x_1} > {x_2}$ and ${x_0} \le t$. Then, the largest root is ${x_0}$, i.e., ${x_0} = {y_0} - \frac{{2a}}{3} + \frac{t}{3} = {g_\lambda }(t)$.

2) We set: $\eta = a - x$ and $y = \eta - \frac{a}{3} + \frac{t}{3}$. So, $x = -y + \frac{{2a}}{3} + \frac{t}{3}$. By a similar analysis as in part (1), there are three distinct roots for the polynomial equation: ${x_0} < {x_1} < {x_2}$ with the smallest solution

where $\varphi = \arccos \left({1 - \frac{{27\lambda {a^2}}}{{2{{(a - t)}^3}}}} \right)$. Therefore, ${x_0} = {g_\lambda }(t)$, when $t < - t_1^ *$.

Appendix B.   Proof of Theorom 1

Proof. We discuss $t = 0$, $t > 0$, and $t < 0$ in the following:

$1)$ $t = 0$: In this case, optimization objective function is ${\theta _{\lambda, t}}(x) = \lambda {\phi _a}(x) + \frac{1}{2}{x^2}$. It is true that $\lambda {\phi _a}(x)$ and $\frac{1}{2}{x^2}$ are both increasing for $x > 0$, and decreasing with $x < 0$. Thus, $\theta (0)$ is the unique minimizer for function ${\theta _{\lambda, t}}(x)$. So ${x^ * } = 0$, when $t = 0.$

$2)$ $t > 0$: Since $\frac{1}{2}{(x - t)^2}$ and $\lambda {\phi _a}(x)$ are both decreasing with $x < 0$, our optimal solution will only be obtained at nonnegative values. Thus, we just need to consider $x \ge 0.$ When $x \ge 0$, we obtain

and

It is clear that $\theta^{''}_{\lambda, t}(x)$ is increasing, and $\theta^{''}_{\lambda, t}(0) = 1 - \frac{2\lambda}{a}$ determines the convexity for the function $\theta_{\lambda, t}(x)$. In the following proof, we further discuss the value of ${x^ * }$ by two conditions: $\lambda < \frac{a}{2}$ and $\lambda > \frac{a}{2}.$

${\rm{a}})$ $\lambda < \frac{a}{2}:$ So, we have $\mathop {\inf }\limits_{x > 0}\theta^{''}_{\lambda, t}(x) = \theta^{''}_{\lambda, t}(0+) = 1-\frac{2\lambda}{a}\geq0$, which means function $\theta'_{\lambda, t}(x)$ is increasing for $x\geq0$, with minimum value $\theta'_{\lambda, t}(0) = \lambda-t = t_2^ *-t.$

${\rm{i}})$ When $0\leq x\leq t_2^ *$, $\theta'_{\lambda, t}(x)$ is always positive, thus the optimal value ${x^ * } = 0$.

${\rm{ii}})$ When $x > t_2^ *$, $\theta'_{\lambda, t}(x)$ is first negative, then positive, and $x > t_2^ * > t_1^ *$. The unique positive stationary point ${x^ * }$ of $\theta'_{\lambda, t}(x)$ satisfies equation: $\theta'_{\lambda, t}({x^ * }) = 0$, which implies

According to Lemma 2, the optimal value is ${x^ * } = x_{0} = g_{\lambda}(t).$

Above all, the value of ${x^ * }$ is

under the condition $\lambda \leq \frac{a}{2}.$

b) $\lambda > \frac{a}{2}:$ In this case, due to the sign of $\theta^{''}_{\lambda, t}(x)$, we know that function $\theta^{'}_{\lambda, t}(x)$ is decreasing at first then switches to be increasing. Its minimum point is $\overline{x} = \sqrt[3]{2\lambda a^{2}}-a$ and the least value

Then, $\theta^{'}_{\lambda, t}(x)\geq t_{1}^{*}-t$ with $x\geq 0.$

${\rm{i}})$ When $0\leq t\leq t_{1}^{*}$, function $\theta_{\lambda, t}(x)$ is always increasing. Thus, optimal value is $x^{*} = 0.$

${\rm{ii}})$ When $t\geq t_{2}^{*}$, $\theta_{\lambda, t}(0+)\leq 0$. Thus, the function $\theta_{\lambda, t}(x)$ is first decreasing and then increasing. There is only one positive optimal stationary point, which is also the optimal solution. Using the Lemma 2, we know that $x^{*} = \theta_{\lambda, t}(x).$

${\rm{iii}})$ When $t_{2}^{*} < t < t_{1}^{*}$, $\theta_{\lambda, t}(0+) > 0$. Thus, the function $\theta_{\lambda, t}(x)$ is first increasing, then decreasing and finally increasing, which implies that there are two positive stationary points and the larger one is a local minima. Using Lemma 2, the local minimized point will be $x_{0} = \theta_{\lambda, t}(x),$ the largest root of (1.1). Since $x_{0}-t+\frac{\lambda a^{2}}{(a+x_{0})^{2}} = 0$, which implies $\frac{\lambda a}{a+x_{0}} = \frac{(t-x_{0})(a+x_{0})}{a}$, we have

Define a function $h(t) = t-g_{\lambda}(t)-\frac{a}{2}$.

First, we prove that $t = t_{3}^{*}$ is a solution to $h(t) = 0.$ Since $\lambda > \frac{a}{2}$, $t_{3}^{*} = \sqrt{2\lambda a}-\frac{a}{2} > 0$. Thus:

Further, by using the relation $\cos(\varphi) = 4\cos^{3}(\frac{\varphi}{3})-3\cos(\frac{\varphi}{3})$ and $0\leq\frac{\varphi}{3}\leq\frac{\pi}{3}$, we have

Plugging this formula into $g_{\lambda}(t_{3}^{*})$ shows that $g_{\lambda}(t_{3}^{*}) = \sqrt{2\lambda a}-a = t_{3}^{*}-\frac{a}{2}.$ Therefore, $t_{3}^{*}$ is the root for $h(t)$ and $t_{3}^{*}\in (t_{1}^{*}, t_{2}^{*}).$

Second, we prove that the function $h(t)$ changes sign at point $t = t_{3}^{*}.$ We prefer to discuss it in two cases.

Case 1: $x \in (t_{3}^{*}, t_{2}^{*})$. According to Lemma 2, we know that $g_{\lambda}(t)$ is the largest root of the cubic polynomial $H(y) = y(a+y)^{2}-t(a+y)^{2}+\lambda a^{2}$ under the condition of $t > t_{1}^{*}$.

For function $H(y)$, we have $H(t) = \lambda a^{2}$ and

Due to $t > t_{3}^{*} = \sqrt{2\lambda a}-\frac{a}{2}$ and $H(t-\frac{a}{2}) < 0$, there is a root $y = g_{\lambda}(t)$, such that $y = g_{\lambda}(t) \in (t-\frac{a}{2}, t)$ for the equation $H(y) = 0$. That is, $t-g_{\lambda}(t) < \frac{a}{2}$, and, thus, $h(t) < 0$. Case 2: $x \in (t_{1}^{*}, t_{3}^{*})$. We have $H(t-\frac{a}{2}) > 0$ and $H(t) > 0$. Due to the proof in Lemma 2, one possible situation is that there are two roots $x_{0}$ and $x_{1}$ within interval $(t-\frac{a}{2}, t)$. However, we can exclude this case. This is because

Furthermore, $x_{0}-x_{1} > \frac{a}{2}$ holds for $t > t_{1}^{*} > \frac{a}{2}$ and $\lambda > \frac{a}{2}$. This is a contradiction of the assumption that $x_{0}$ and $x_{1}$ are in $(t-\frac{a}{2}, t)$. Thus, $H(y) = 0$ has no root in $(t-\frac{a}{2}, t)$.

Therefore, $x_{0} = g_{\lambda}(t) < t-\frac{a}{2}$. That is to say, $h(t) > 0$.

From the discussion earlier, it is true that the optimal solution $x^{*} = 0$, if $0 < t\leq t_{3}^{*}$, and $x^{*} = x_{0} = g_{\lambda}(t)$, if $t > t_{3}^{*}.$

To sum up, we know that under the condition $\lambda > \frac{a}{2}$,

3) $t < 0:$ Because

$x^{*}(t) = -x^{*}(-t),$ which implies that the formula obtained when $t > 0$ above can extend to the case: $t < 0$ by odd symmetry.

Summarizing results from all cases, the proof is completed.

Appendix C.   Proof of Corollary 1

Proof. It is clear to see that the function

is differentiable on ${\mathbb R}\backslash \{ 0\}$. For $x \ne 0$, the derivative of $\theta_{\lambda, t}(x)$ is given by

Let us find the range of $a$ for which $\theta_{\lambda, t}(x)$ is convex. Notice that

and

When we set $a \ge 2\lambda$, the function $\theta'_{\lambda, t}$ is increasing, then the function $\theta_{\lambda, t}(x)$ is strictly convex.