Conjugate gradient variants for ℓp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\ell}_{p}$$\end{document}-regularized image reconstruction in low-field MRI

We consider the MRI physics in a low-field MRI scanner, in which permanent magnets are used to generate a magnetic field in the millitesla range. A model describing the relationship between measured signal and image is derived, resulting in an ill-posed inverse problem. In order to solve it, a regularization penalty is added to the least-squares minimization problem. We generalize the conjugate gradient minimal error (CGME) algorithm to the weighted and regularized least-squares problem. Analysis of the convergence of generalized CGME (GCGME) and the classical generalized conjugate gradient least squares (GCGLS) shows that GCGME can be expected to converge faster for ill-conditioned regularization matrices. The ℓp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\ell}_{p}$$\end{document}-regularized problem is solved using iterative reweighted least squares for p=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=1$$\end{document} and p=12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=\frac{1}{2}$$\end{document}, with both cases leading to an increasingly ill-conditioned regularization matrix. Numerical results show that GCGME needs a significantly lower number of iterations to converge than GCGLS.


Introduction
In low-field magnetic resonance imaging (MRI), magnetic field strengths in the millitesla (mT) range are used to visualize the internal structure of the human body. In traditional MRI scanners, magnetic field strengths of several tesla are the norm. While these high-field MRI scanners yield images of excellent quality, their cost, size and infrastructure demands make them unattainable for developing countries. Therefore, the design of low-field MRI scanners is of great clinical relevance. This research is part of a project that aims toward creating an inexpensive low-field MRI scanner using a Halbach cylinder that can be used for medical purposes. A Halbach cylinder is a configuration of permanent magnets that generates a magnetic field inside the cylinder and a very weak, or in the ideal case, no magnetic field outside of it. Imaging can be done by making use of the variations in the magnetic field. However, the resulting reconstruction problem is very ill-posed. This is due to the nonlinearity of the magnetic field inside the Halbach cylinder that we consider. This field leads to non-bijective mappings and potentially gives rise to aliasing artifacts in the solution. Additionally, in the center of the cylinder, there is very little variation in the field, limiting the spatial resolution in that area. Another complication we face is low signal-to-noise ratios. Nevertheless, in a similar project, Cooley et al. [6] have shown that it is possible to reconstruct magnetic resonance images given signals obtained with a device based on a Halbach cylinder, using a simplified signal model in which similar assumptions are made as in high-field MRI. In this paper, we revisit the underlying physics and formulate the general signal model for MRI without making these assumptions.
Regularization is required to limit the influence of noise on the solution of the image reconstruction problem as much as possible. In this paper, we reformulate the weighted and regularized least-squares problem such that the conjugate gradient minimal error (CGME) method (see for example [2]) can be used to solve it for nontrivial covariance and regularization matrices, filling a gap in existing literature as far as we know. We do this by deriving the Schur complement equation for the residual. A similar approach is taken by Orban and Arioli [25] to derive generalizations of the Golub-Kahan algorithm. Using these algorithms, they formulate generalizations of LSQR, Craig's method and LSMR (see [2]) for the general regularization problem. We explain in which cases generalized CGME (GCGME) may have an advantage over generalized conjugate gradient least squares (GCGLS). Additionally, we apply GCGME to MRI data with different types of regularization.
The present paper results from our efforts to address the challenges of low-field MRI using advanced image processing. It is interdisciplinary in nature, with an emphasis on image reconstruction techniques. The contributions of this paper include a signal model for low-field MRI that does not rely on any field assumptions as encountered in highfield MRI. Also, a new generalization of the conjugate gradient method is presented for the weighted and regularized least-squares problem, including an analysis of when this generalization is expected to perform best. Although we focus on a low-field MRI setting, this algorithm is generally applicable to p -regularized least-squares problems.

Low-field MRI
In magnetic resonance imaging (MRI), the internal structure of the body is made visible by measuring a voltage signal that is induced by time variations of the transverse magnetization within a body part of interest. Based on this measured signal, an image of the spin density of different tissue types may be obtained.
To be specific, first the body part of interest is placed in a static magnetic field ⃗ B = B 0 (⃗ r) ⃗ i x that is oriented in the x-direction in our Halbach measurement setup (see Fig. 1a) with a position-dependent x-component B 0 = B 0 (⃗ r) . A net magnetization will be induced that is oriented in the same direction as the static magnetic field. In the above expression, = 267 × 10 6 rad s −1 T −1 is the proton gyromagnetic ratio, ℏ = 1.055 × 10 −34 m 2 kg s −1 is Planck's constant divided by 2 , k B = 1.381 × 10 −23 m 2 kg s −2 K −1 is Boltzmann's constant, and T is the temperature in kelvin.
Subsequently, a radiofrequency pulse is emitted to tip the magnetization toward the transverse yz-plane. After this pulse has been switched off (in our model at t = 0 ), the magnetization rotates about the static magnetic field with a precessional frequency (also known as the Larmor frequency) given by and will relax back to its equilibrium given by Eq. (1). During this process, an electromagnetic field is generated that can be locally measured outside the body using a receiver coil. This measured signal is amplified, demodulated, and low-pass filtered, and for the resulting signal, we have the signal model [23]: where is the domain occupied by the body part of interest, T 2 (⃗ r) is the transverse relaxation time, c(⃗ r) is the socalled coil sensitivity with amplification included, M ⟂ (⃗ r, 0) is the transverse magnetization at t = 0 , and is the difference between the Larmor frequency and the demodulation frequency that is used. For this demodulation frequency, we take the frequency that corresponds to the static magnetic field at the center of our imaging domain.
Furthermore, using Eq. (2) in the expression for M 0 , we have and since the initial transverse magnetization M ⟂ (⃗ r, 0) is proportional to M 0 (⃗ r) , we can also write our signal model as: where it is understood that all remaining proportionality constants have been incorporated in the coil sensitivity c(⃗ r) . Conventionally, the spatial dependence of is ignored. Therefore, the 2 term usually does not appear in MRI literature. However, we incorporate it into our model because of the relatively large inhomogeneities in the magnetic field we are considering. We remark that Eq. (5) is a general MRI signal model, but it is more suitable for low-field MRI because the assumptions made for high-field MRI (namely, a very strong and homogeneous magnetic field) do not hold for low field. Ignoring T 2 relaxation, the final signal model becomes The measurements taken in an MRI scanner consist of noisy samples of the signal given by Eq. (6): where b i denotes the ith sample of the signal, measured at time t i . L is the number of time samples, and e i is the measurement error.

Model-based image reconstruction
In high-field MRI, the magnetic field is manipulated in such a way that Eq. (6) constitutes a Fourier transform. The resulting linear problem is well posed, and the image can be efficiently obtained using an inverse FFT. However, in low-field MRI, the magnetic field is usually strongly inhomogeneous, which prevents us from using standard FFT routines. Model-based image reconstruction can be applied instead [10].
In order to estimate ( ) , we write it as a finite series expansion of the form: where (⋅) denotes the object basis function, j is the center of the jth basis function and x j are the coefficients. Usually, rectangular basis functions are used, in which case N is the number of pixels. Combining Eqs. (6) and (8) yields where When the basis functions are highly localized, a "center of pixel" approximation can be used: Here, x y is the pixel size and z is the thickness of the slice that is being imaged. Combining Eqs. (7) and (8) yields one system of equations: where the elements of are described by Eq. (11). This problem is ill-posed due to the nature of the magnetic field that is present within the Halbach cylinder. As shown in Fig. 2, the field has a high degree of symmetry. The precessional frequency depends linearly on the magnitude of the field, which means that several pixels will correspond to the same frequency. Therefore, using only one measured signal, it is impossible to determine the contribution of each pixel to the signal. By rotating the object to be imaged and hence obtaining a multitude of different signals corresponding to different rotations of the same object, we plan to mitigate this problem. The same approach was taken by Cooley et al. [6] ( Table 1).

Methodology
The model that is used to reconstruct is given by the linear system of Eq. (12). We can attempt to solve for by finding a solution to the least-squares problem This can be done by applying the conjugate gradient method introduced by Hestenes and Stiefel in 1952 [19] to the normal equations with H denoting the Hermitian transpose of . The conjugate gradient method tailored to Eq. (14) was proposed in [19] and is usually referred to as conjugate gradient for least squares (CGLS). The difference with the standard conjugate gradient method lies in the increased stability of the CGLS method. A review of the literature reveals that this method is known by other names as well.
H = H , Table 1 Overview of the main matrices and vectors used in this work N is the total number of pixels in the image and M is the total number of acquired data points. F can have an arbitrary number of rows, which determines the dimensionality of D

Symbol Description Dimensions
On the other hand, the second normal equations can be solved using the conjugate gradient method as well. In the literature, this is usually called conjugate gradient minimal error (CGME). However, in [1] it is called conjugate gradient normal error (CGNE), while [30] uses the term Craig's method. It was introduced by Craig in 1955 [7]. CGLS and CGME are discussed by Björck in [2], Hanke in [14] and Saad in [29]. While CGLS minimizes the residual = − in the 2 norm over the Krylov sub- , CGME minimizes the error (over the same subspace). The main drawback of this latter method is that, in theory, it only works for consistent problems for which ∈ ( ) . This means that the method is of limited use for most problems in practice, because the presence of noise renders the system inconsistent. In [21], this problem is circumvented by defining an operator that projects onto the column space of . Subsequently, = can be solved using CGME. The obvious disadvantage of this method is that has to be calculated and stored.

Regularization of the problem
Regularization of an ill-posed problem aims to make the problem less sensitive to noise by taking into account additional information, i.e., it aims at turning an ill-posed problem into a well-posed one. Like many iterative methods, both CGLS and CGME have a regularizing effect if the iterating procedure is stopped early: keeping the number of iterations low keeps the noise from corrupting the result too much. If a large number of iterations is used, noise can have a very strong effect on the solution. The regularizing properties of CGLS were established by Nemirovskii in [24] and are discussed in [2,9,14], among others. CGME's regularizing effect was shown by Hanke in [15]. However, we are interested in what Hansen [17] calls general-form Tikhonov regularization, i.e., adding a regularization term to minimization problem (13), leading to where is a weighting matrix, and is a Hermitian positive definite matrix. Using a CG algorithm to solve Eq. (16) is a natural choice [10]. The CG method is often used to solve image reconstruction problems in MRI when a conventional Fourier model is insufficient (see for example [11,27,34]). Additionally, it is used as a building block for other algorithms used in MRI by Pruessman [26], Ramani and Fessler [28] and Ye et al. [38], among others. It is straightforward to generalize CGLS to regularized and weighted (15) least-squares problems of the form of Eq. (16). In this case, because of the well-posedness of the resulting minimization problem, the noise does not influence the solution as much as when Eq. (13) is considered and increasing the number of iterations does not lead to a noisier solution. In this paper, we will use = −1 , where is the covariance matrix of the noise: For our application, the noise can be considered to be white, which means that = . However, for completeness, we consider the general case. In case = , Eq. (17) reduces to a minimization problem with standard Tikhonov regularization [36]. The optimal value of the regularization parameter is usually unknown. An approach that is often used to find a suitable value is the L-curve method [16]. By taking the gradient and setting it equal to , the normal equations are obtained: Again, the conjugate gradient method can be used to solve Eq. (18). We will use the term GCGLS (generalized CGLS) to refer to the conjugate gradient method applied to the normal Eq. (18). Saunders [30] extended Craig's method, which is mathematically equivalent to CGME, to the regularized leastsquares problem with = and = . He introduces an additional variable and considers the constrained minimization problem By defining ̃ = √ = − , he shows that this constrained minimization problem is equivalent to = is consistent, and hence, Eq. (19) can be solved using CGME. Unfortunately, no advantages to using CGME were found. Note that such a reformulization is necessary because the standard way of including the regularization matrix = , by simply solving the so-called damped least-squares problem (20) � √ � = � � using CGME, is not possible, due to the inconsistency of the system. Reformulation of CGME for general-form regularization can be achieved using a Schur complement approach as will be shown as follows. Again, we consider Eq. (17). We introduce the variable = −1 ( − ), and we note that || − || 2 −1 = || || 2 . Then, minimization problem (17) can be formulated as a constrained minimization problem: and using the technique of Lagrange multipliers, we find that If we eliminate from Eq. (23), the original normal Eq. (18) is obtained, whereas if we assume is invertible and we subsequently eliminate , we end up with a different set of equations. As mentioned before, the first option leads to the GCGLS method. The latter approach leads to the GCGME method.

GCGLS
By applying the conjugate gradient method to Eq. (18) and making some adjustments to increase stability (see [2] for details), the GCGLS algorithm is obtained: Here, M is the total number of data points measured and N is the number of pixels in the image. The residual of the normal Eq. (18) is denoted by k . We remark that the vectors on the left side can be overwritten by the vectors on the right. Only eight vectors have to be stored, namely , , , , , , and −1 . Note that the recursion for k+1 is included to avoid an extra multiplication with . It can be ignored in case = . In this algorithm, only three matrix-vector multiplications are carried out per iteration: k+1 , H k and k . Additionally, one system with has to be solved (if ≠ ). A slightly different formulation of the GCGLS algorithm can be found in [34].

If
is invertible, can be eliminated from Eq. (23), yielding Subsequently, can be obtained from as: In [25], Arioli and Orban derive a generalization of Craig's method [7] based on Schur complement (24). Below, we formulate a similar generalization of the CGME method applied to this system. We are not aware this generalization of CGME has been formulated elsewhere.
in which max ( ) and min ( ) are the largest and smallest eigenvalue of , respectively. In this section, we bound the condition numbers of the two Schur complement matrices in Eqs. (18) and (24) to gain insight into when GCGME can be expected to perform better than GCGLS, and vice versa. Given two HPD matrices and , the following bound on the condition number holds: This inequality follows from Weyl's theorem [37], which states that for eigenvalues of Hermitian matrices and , the following holds: Here, i ( ) denotes any eigenvalue of the matrix . For GCGLS, we have that and, using the following inequalities with max ( ) the largest singular value of , we get that Analogously, for CGME, we have Here, k is the residual of the normal Eq. (24). Note that the original CGME algorithm can be recovered from the generalized CGME algorithm given above by taking 1 = and = , the zero matrix. Only seven vectors have to be stored, namely , , , , , −1 and . Like GCGLS, GCGME needs four matrix operations per iteration: k , H k , −1 k and −1 k . We remark that there is an essential difference between GCGLS and GCGME. GCGLS iterates for the solution vector , and the equality k = −1 ( − k ) is explicitly imposed. The equality k = 1 −1 H k is not enforced and is only (approximately) satisfied after convergence. GCGME, on the other hand, iterates for k . The equality k = 1 −1 H k is enforced, while k = −1 ( − k ) is only satisfied approximately after convergence.

Convergence of GCGLS and GCGME
The convergence of the conjugate gradient method depends on the condition number of the system matrix. Again, suppose that CG is used to solve the system = for the unknown vector , where is a Hermitian positive definite (HPD) matrix and is a known vector. Then, the following classical convergence bound holds [2]: where 2 ( ) is the 2 -norm condition number of , which, for HPD matrices, is equal to and using similar manipulations as above we obtain These inequalities indicate that if GCGLS can be expected to perform best, and that if GCGME should be preferred. This latter situation may occur when the regularization term is minimized in the p -norm with p ∈ (0, 1] , as we will discuss in the next section.

Types of regularization
Instead of an 2 -penalty, we will consider the more general case of an p penalty with p ∈ (0, 2] . Then, the minimization problem becomes A vast literature regarding this 2 p minimization problem is available. In for example [3,4,20,22], this problem is solved using a majorization-minimization approach. In this work, we will focus on the classical approach using iterative reweighted least squares (IRLS), also known as iterative reweighted norm (IRN), see for example [2], for solving minimization problem (37), in which GCGLS and GCGME can be used as building blocks. Their performances will be compared. We choose the IRLS algorithm for three reasons: its simplicity, the fact that it is a wellknown technique and that in this algorithm, the regularization matrix changes in each iteration, which makes it especially interesting for us, because we can test whether GCGME indeed performs better in case Eq. (36) holds. This work is not meant to evaluate the performance of IRLS as a solver for Eq. (37), and we do not compare it with other methods. For completeness, however, we do mention that we could also have chosen to evaluate both approaches as a building block of the split Bregman method [12] for the 1 -regularized problem, for example. In [4], Chan and Liang use CG as a building block for their half-quadratic algorithm that solves Eq. (37) as well. A comparison between GCGLS and GCGME could be carried out in this context too.
IRLS is an iterative method that can solve an p -regularized minimization problem by reducing it to a sequence of 2 -regularized minimization problems. Note that for a vector of length N, so Furthermore, is some regularizing matrix. Note that Eq. (37) can be rewritten as: where and | | is the element-wise modulus of . This is simply another instance of minimization problem (17), with = H . However, now depends on . So, when the kth iterate k is known, k+1 is found as follows: where This is repeated until convergence. Furthermore, in Eq. (43), is a small number that is added to the denumerator to prevent division by zero. We will use = 10 −6 . We observe that in each IRLS iteration, we simply encounter an instance of minimization problem (17) again with k = H k , which can be solved using either GCGLS or GCGME. When carrying out calculations with −1 k , we will use Due to the sparsity-inducing property of the p penalty when p ≤ 1 (see for example [8]), −1 k = diag | k | will contain an increasing number of entries nearly equal to zero. In cases where is an invertible matrix, When GCGME is used, we can take advantage of this structure, instead of calculating k and working with its inverse. Moreover, when is an orthogonal matrix, no additional computations are necessary to compute inverses. The regularization matrix = H k will become illconditioned when elements of k become small. Therefore, we expect that, when combined with IRLS, GCGME will perform better than GCGLS for p ≤ 1 . Numerical experiments are carried out to investigate this further.

Different choices for p
We will minimize the following 1 -regularized least-squares problem and the 1∕2 -regularized least-squares problem to obtain approximations to the optimal solution . For a general , this results in the following two minimization problems: and We note that in the latter case, the objective function is not convex which means that the obtained solution does not necessarily correspond to a global minimum, see for example [5]. For each of these two minimization problems, we will consider two different regularization operators.

Regularizing using the identity matrix
First, we set = . In case the 1 penalty is used, the minimization problem reduces to This is known as least absolute shrinkage and selection operator (LASSO) regularization which was first introduced by Tibshirani in [35]. If the regularization parameter is set to a sufficiently high value, the resulting solution will be sparse. The same holds for the 1∕2 -regularized minimization problem: The rationale behind choosing this type of regularization is the fact that the intensity of many pixels in MRI images is equal to 0. In both cases ( p = 1 and p = 1∕2 ), the regu- This is especially useful for GCGME, because calculating the product of −1 and a vector is trivial in this case.

Regularizing using first-order differences
Additionally, we consider the case where is a first-order difference matrix that calculates the values of the jumps between each pair of neighboring pixels. Suppose our  image consists of n × n pixels. If we define the 1D firstorder difference operator 1D ∈ ℝ n×n the 2D first-order difference matrix is given by: where ⊗ denotes the Kronecker product. This type of regularization is known as anisotropic total variation regularization. A reason for choosing = is that neighboring pixels are very likely to have the same values in MR images. This is due to the fact that neighboring pixels tend to represent the same tissue. However, is not a square matrix, which means that, in the 1 case, k has to be calculated explicitly and then inverted when GCGME is used. Although this makes regularization with first-order differences in combination with GCGME less attractive than with GCGLS, we do include this technique to investigate the relative reconstruction quality of this widely used regularization method. The resulting minimization problems are equal to Eqs. (45) and (46) with = :

Four different minimization problems
We will investigate all four minimization problems (47), (48), (51) and (52). Since the least-squares term is the same in all four minimization problems, the difference between them lies in the penalty term used, as summarized in Table 2. In each of the four cases, we will use both GCGLS and GCGME to compare their rate of convergence.

Numerical simulations
For our simulations, we use a simulated magnetic field as shown in Fig. 1a. (We also have access to a measured field map, but it is measured on a very coarse grid, making it unsuitable for our purposes.) The magnetic field within the FoV of 14 cm by 14 cm is clearly inhomogeneous, as shown in Fig. 2. The magnetic field has an approximately quadrupolar profile. This is because the Halbach cylinder is designed to generate a field that is as uniform as possible. However, due to practical limitations, such as the finite length of the cylinder, this uniformity cannot be attained, leading to a quadratic residual field profile. See for example [6,18]. We do not use a switched linear gradient coil, as is done in conventional MRI. Instead, the inhomogeneous background field is used for readout encoding. For a thorough exploration of the use of non-bijective encoding maps in MRI, we refer to [13,18,[31][32][33]. Performing slice selection in the presence of a nonhomogeneous background field is nontrivial, but this complication is ignored here. We assume that the entire measured signal originates from one slice. We simulate the signal generation inside the Halbach cylinder using Eqs. (11) and (12). The dwell time is set to t = 5 × 10 −6 , and the readout window is 0.5 ms, leading to 101 data

B [mT]
x (cm)  points per measurement. Additionally, the field is rotated by 5° after each individual measurement, so in order to cover a full circle, 72 different angles are considered. We note that this is similar to a radial frequency-domain trajectory dataset in conventional MRI. In [18], quadrupolar fields are used to generate such a dataset. However, the field we are using is only approximately quadrupolar, so it is not a true radial frequency-domain trajectory experiment. The system consists of 72 × 101 = 7272 equations. The numerical phantom of 64 × 64 pixels is shown in Fig. 1b, resulting in a matrix of size 7272 × 4096 . We assume that the repetition time T R is long enough for the magnetization vector to relax back to its equilibrium. Also, the echo time is assumed to be so short as to make T 2 -weighting negligible.
Since the background field is almost homogeneous in the center, as shown in Fig. 2, we decided to place the object of interest in the numerical phantom off-center. Within a homogeneous region in the field, distinguishing between the different pixels is impossible. Another obstacle in the reconstruction process is the fact that the background field is almost symmetrical in both the x-and the y-axis, potentially leading to aliasing artifacts in the lower half of the image (because the object of interest is placed in the upper half of the image). We could reconstruct by leaving out all the columns in matrix corresponding to the pixels in the lower half of the image. Another way of circumventing this problem is by using several receiver coils with different sensitivity maps to break the symmetry of the problem [18,33]. However, we choose not to take these approaches, so we can see how severe these artifacts are for the different objective functions.
The coil sensitivity c is assumed to be constant, so it is left out of the calculations. White Gaussian noise is added, so the covariance matrix is simply the identity matrix. We assume an SNR of 20. The numerical experiments are carried out using MATLAB version 2015a. Often, CG is stopped once the residual is small enough. However, GCGLS and GCGME are solving different normal equations, so the residuals are different for both methods. Therefore, a comparison using such a stopping criterion would not be fair. Instead, a fixed number of CG iterations is used per IRLS iteration. The value of the regularization parameter is chosen heuristically. The number of IRLS iterations is set to 10. We consider both 10 and 1000 CG iterations per IRLS iteration. The initial guess 0 in GCGLS (and 0 in GCGME) is the zero vector. During the first IRLS iteration, we set = , which means that = * . After the first IRLS iteration, we calculate the weight matrix according to Eq. (43). We use warm starts, i.e., we use the final value of our iterate k (or k for GCGME) of the previous IRLS iteration as an initial guess for the next IRLS iteration. Table 3 shows the parameters that were chosen for all four different minimization problems. The regularization parameter was chosen heuristically in each case.

Results and discussion
All resulting images are shown in Fig. 3. We note that in all cases (except perhaps the ‖ ‖ 1∕2 1∕2 one), GCGME yields a result that resembles the original more than GCGLS does. GCGLS tends to yield aliasing artifacts in the lower half of the image. This effect is less pronounced for the GCGME results, especially when ‖ ‖ 1∕2 1∕2 is used as the penalty term. The objective function value is plotted as a function of the iteration number in Fig. 4. We see that GCGME attains a lower objective function value in all cases. However, both methods should in theory converge to the same value for the ‖ ‖ 1 -and ‖ ‖ 1 -penalty terms. Evidently, GCGLS has not converged yet. If we increase the number of CG iterations to 1000, GCGLS and GCGME converge to the same result, as can be seen in Appendix B. The GCGME result is the same, whether 10 or 1000 CG iterations are carried out, which means that GCGME has already converged in the first case. However, GCGLS needs a  significantly larger number of iterations to converge. In case = , GCGLS and GCGME both need 0.069 s per iteration. When = , GCGME needs slightly more time per iteration than GCGLS: 0.072 versus 0.069 s.

Discussion of the results
GCGLS needs a large number of CG iterations to converge, while for GCGME, this number is low (typically, 10 is sufficient). This can be explained by the observation that as we get closer to the solution, many elements of the vector | k | 2−p will converge to zero, due to the sparsity-enforcing properties of the p penalty when p ≤ 1 . Therefore, −1 k = diag | k | 2−p will contain an increasing number of very small entries, which means that the matrix k = H k will become more and more ill-conditioned as the number of IRLS iterations grows. That means that, after a few IRLS iterations, 2 ( k ) ≫ 2 ( ) will hold, in which case GCGME performs better than GCGLS, which is consistent with our results.
It is interesting to note that when the number of CG iterations for GCGLS is set to 10, GCGLS appears to have reached convergence after 4-5 IRLS iterations, yielding an image with aliasing artifacts in the form of an additional shape in the lower half of the image, as well as regions of intensity in the corners of the image. However, convergence is not actually attained yet. The number of CG iterations needs to be increased to a 1000 before convergence is reached.
We observe that the ‖ ‖ 1∕2 1∕2 penalty is best at repressing the aliasing artifacts in the lower half of the image.

Conclusion
We formulated a general MRI signal model describing the relationship between measured signal and image which is more suitable for low-field MRI because the assumptions that are usually made in high-field MRI do not hold here. The discretized version yields a linear system of equations that is very ill-posed. Regularization is needed to obtain a reasonable solution. We considered the weighted and regularized least-squares problem. A second set of normal equations was derived, which allowed us to generalize the conjugate gradient minimal error (CGME) method to include nontrivial weighting and regularization matrices. We compared our GCGME method to the classical GCGLS method by applying both to data simulated using our signal model. Different regularization operators were considered: the identity matrix and the anisotropic total variation operator that determines the size of the jumps between neighboring pixels. The regularization term was measured in the 1 -norm and the 1 2 -norm, and iterative reweighted least squares (IRLS) was used to solve the resulting minimization problems. In each IRLS iteration, an 2 -regularized minimization problem was solved using GCGLS or GCGME. GCGME converges much faster than GCGLS, due to the regularization matrix becoming increasingly ill-conditioned as the number of IRLS iterations increases. This makes GCGME the preferred algorithm for our application.