Phase Error Correction in Magnetic Resonance: A Review of Models, Optimization Functions, and Optimizers in Traditional Statistics and Neural Networks

1. Introduction

Magnetic Resonance (MR) techniques, including Nuclear Magnetic Resonance (NMR) spectroscopy and Magnetic Resonance Imaging (MRI), play a central role in fields ranging from chemical analysis to medical diagnostics. These techniques rely on the precise interpretation of phase-sensitive data to ensure accurate results. However, the presence of phase errors—caused by factors such as eddy currents, motion artifacts, and timing discrepancies—poses a significant challenge, often distorting critical information.

To address phase errors, researchers have developed methods that use multiple acquisitions to mitigate distortions. For instance, phase errors induced by eddy currents can be minimized by combining data collected under positive and negative magnetic fields (Bieri et al., 2005; Jiru, 2008). Similarly, in phase-contrast MRI, gradient pulses applied in opposing directions effectively reduce motion-induced phase errors, improving the precision of velocity measurements (Wymer et al., 2020). These techniques illustrate the effectiveness of multi-acquisition strategies in achieving more accurate results.

However, practical constraints often limit the feasibility of collecting multiple data sets or acquisitions, particularly in high-throughput or time-sensitive applications. In these cases, phase error correction must rely on computational methods applied to single data sets, placing greater demands on the robustness and accuracy of automatic correction algorithms. While early efforts, such as Ernst's proposal for automatic phase correction over fifty years ago (Ernst, 1969), showed promise, subsequent developments have revealed significant limitations. Many automatic methods struggle with complex datasets, leading to inconsistencies or even introducing additional errors (Minderhoud et al., 2020). As a result, manual phase correction remains widely used, either in conjunction with automatic corrections (Emwas et al., 2018) or as a standalone method (Canlet et al., 2023). Manual correction is still recommended for achieving more accurate results (Ben-Tal et al., 2022).

The persistence of manual correction underscores the need for more reliable automatic methods. Recent research has turned its focus toward innovative solutions, exploring approaches such as advanced statistical models and machine learning techniques. These efforts aim to enhance the accuracy of automatic corrections and minimize reliance on manual interventions. For example, the integration of neural network-based models has shown promise in tackling nonlinear phase errors, offering potential pathways to fully automated solutions.

This review seeks to address these evolving challenges by providing a comprehensive analysis of current phase correction methodologies across MR modalities, including NMR, MRI, and others. Unlike previous studies that focused narrowly on specific techniques (Binczyk et al., 2015; Haskell et al., 2023), this work takes a broader view, examining the foundational models, optimization functions, and optimizers used in both traditional statistical methods and emerging neural network architectures. By critically evaluating the strengths and limitations of these approaches, this review aims to highlight key opportunities for innovation.

As MR technologies continue to advance, developing effective phase correction methods is essential to fully realize their potential. This review will summarize the state-of-the-art techniques, identify persistent challenges, and propose future directions for research to drive progress in this crucial area.

2.1. Global phase error correction using a uniform phase adjustment value

The simplest model for phase correction involves applying a global phase adjustment across the entire spectrum. This can be expressed by the formula:

$$\theta =c$$

(1)

In this formula, θ represents the phase correction value applied uniformly across the spectrum. The value c is the global adjustment parameter, determined to cancel out the phase error by being equal in magnitude but opposite in sign to the global phase deviation of the spectrum.

The main challenge lies in determining the appropriate value for c. Ernst (Ernst, 1969) introduced a method based on integrals to estimate global phase errors; this method remains widely used today. The global phase correction model can be formulated as follows:

$$\theta =-\mathit{arctan}\left({\displaystyle \frac{\int D\left(x\right)dx}{\int A\left(x\right)dx}}\right)$$

(2)

In this formula, arctan, the inverse tangent function, produces an angle representing the phase error. The variable x represents the frequency- or location-related variable, D(x) denotes the observed dispersion spectrum, and A(x) stands for the observed absorption spectrum. A(x) and D(x) represent the real and imaginary parts of the complex data in the MR frequency domain, which is commonly referred to as the frequency domain for one-dimensional (1D) MR data and as k-space for multi-dimensional MR data. The integrals in the formula are calculated over the entire spectrum, as global phase correction aims to address errors across all spectral data points.

While Ernst’s original paper (Ernst, 1969) did not explicitly include a formula involving the dispersion spectrum, it utilized the Hilbert transform of the observed frequency domain data. In theory, when there is no signal decay over time, a pure dispersion spectrum is equivalent to the Hilbert transform of the corresponding pure absorption spectrum, which shifts the phase by 90 degrees. This 90-degree phase difference is a key concept used for phase adjustment in Džakula’s research (Džakula, 2000).

However, real MR signals are subject to decay and various types of errors over time. Therefore, it is more accurate to use the observed dispersion spectrum directly rather than relying on the Hilbert-transformed spectrum. This approach has been adopted in subsequent research. Nonetheless, when the observed dispersion spectrum is unavailable, the Hilbert-transformed spectrum remains a viable alternative.

2.2. Linear phase correction model with zero and first order parameters

The global phase correction model, while useful, does not account for the variability of phase errors across the spectrum. This limitation arises because phase errors can vary with frequency and are not uniformly distributed. To address these variations, a linear phase correction model with zero and first-order parameters was introduced. This model provides a more sophisticated approach by incorporating frequency-dependent terms.

For one-dimensional (1D) frequency domain data, the phase correction model can be expressed as (Balacco, 1994; Bao et al., 2013; Brown et al., 1989; Chen et al., 2002; E. C. Craig & Marshall, 1988; de Brouwer, 2009; Džakula, 2000; Heuer, 1991; Jiang, 2024a; Jiru, 2008; Montigny et al., 1990; Morris, 2017; Prostko et al., 2022; Rout et al., 2023; Sotak et al., 1984; Van et al., 2011; Worley et al., 2015; Worley & Powers, 2014; Wright et al., 2012):

$$\theta ={\beta }_0+{\beta }_{1}x$$

(3)

In this formula, β₀ and β₁ represent the zero and first-order parameters, respectively. The variable x can be any function related to position or frequency in the observed absorption spectrum, such as a scaled index (Chen et al., 2002).

For two-dimensional (2D) frequency domain data, the phase correction model can be extended as (Feinauer et al., 1997; Larry Bretthorst, 2008; Moussavi et al., 2014):

$$\theta ={\beta }_0+{\beta }_{1}x+{\beta }_{2}y$$

(4)

Here, x and y represent the two spatial frequency directions in k-space, which refers to the frequency domain of MR data with more than one dimension, such as 2D, 3D, or even 4D and beyond.

2.3. Linear phase correction model with higher-order terms

In cases where phase errors exhibit non-linear characteristics, the linear model can be further extended to include higher-order terms. This extension allows for a more accurate representation of complex phase error patterns.

For 1D MR data, a linear phase correction model with higher-order terms can be expressed as (Gan & Hung, 2022; Hardy et al., 2009; Hoffman et al., 1992; Jaroszewicz et al., 2023; Vanvaals & Vangerwen, 1990; Worley & Powers, 2014):

$$\theta ={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\dots$$

(5)

For 2D MR data, the model can be written as (Zheng Chang & Qing-San Xiang, 2005):

$$\theta ={\beta }_0+{\beta }_{1}x+{\beta }_2\mathrm{y}+{\beta }_3{x}^2+{\beta }_4{y}^2+\dots$$

(6)

This approach can also be extended to three-dimensional (3D) and four-dimensional (4D) MR data, providing enhanced flexibility for correcting non-linear phase errors.

2.4. Neural networks with multiple layers

As discussed in previous sections, traditional phase correction models have evolved from simple global phase shifts to more sophisticated linear models that address frequency-dependent errors and incorporate higher-order terms. While these advancements have significantly improved phase correction, they still face limitations when dealing with complex, non-linear phase error patterns. To overcome these challenges, modern techniques have introduced advanced machine learning approaches, such as neural networks.

Neural networks, particularly those with multiple layers, offer a promising alternative by learning complex, non-linear patterns directly from data. Unlike traditional models that rely on predefined mathematical formulation, neural networks can adaptively learn and correct phase errors based on extensive training data. This adaptability allows them to capture intricate relationships between phase errors and spectral features, making them effective for handling non-linear phase error patterns.

To understand neural networks, consider them as a series of interconnected nodes (neurons), where each node performs a mathematical operation—typically a weighted sum followed by an activation function. In this context, a simple linear regression model can be seen as a basic form of a neural network, consisting of just an input layer and an output layer.

However, a typical neural network includes one or more hidden layers, each applying non-linear activation functions to model complex patterns in the data. The network's final output is generated through a sequence of transformations across all layers, with each layer enhancing the model's ability to correct phase errors effectively. Figure 1 illustrates a neural network architecture consisting of five layers: an input layer, three intermediate hidden layers, and an output layer. This structure demonstrates how the network processes data through multiple stages to improve phase correction.

2.4.1. Unsupervised neural network learning

Neural network learning can mainly be classified into supervised and unsupervised approaches, with supervised training being the dominant method. The primary difference lies in the use of ground truth data: supervised learning requires ground truth phase values to be included in the MR dataset during training, enabling the trained model to predict and adjust phase values for new datasets. In contrast, unsupervised learning identifies patterns and structures in the data without relying on ground truth values.

In the context of phase correction, unsupervised learning becomes particularly valuable since phase errors are not directly observable, and these statistical models we have discussed earlier are typically applied within an unsupervised learning framework.

Now, let’s delve deeper into how unsupervised neural network learning can be applied to phase correction. Based on the ideas from Shamaei (Shamaei et al., 2023), Figure 2 illustrates a neural network architecture designed for unsupervised phase correction using deep autoencoders. This architecture comprises several nonlinear hidden layers and includes two main components: an encoder and a decoder.

The encoder processes MR signals in complex format, converting them into a simpler form that the network can more easily interpret. It begins with a dropout layer that randomly removes some of data points during each training step. This is followed by six convolutional blocks and a fully connected layer, which combines the features learned by the previous layers to make final predictions.

Each convolutional block includes a one-dimensional convolutional layer (Conv1D) followed by a Rectified Linear Unit (ReLU) activation function.

Convolution is just to merge two functions, illustrating how one function is altered by the other. One-dimensional convolution can be expressed as (Talvila, 2009):

$$\left(f\ast g\right)\left(t\right)=f\left(t\right)\ast g\left(t\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f\left(\tau \right)g\left(t-\tau \right)d\tau$$

(7)

In digital convolution, a small matrix, known as a filter or kernel, is applied to a larger input matrix by sliding it across the input and performing element-wise multiplications. The results of these multiplications are summed to generate a single number for each position of the filter, creating an output matrix or feature map. Unlike simple matrix multiplication, which involves dot products of rows and columns, convolution captures local patterns by repeatedly applying the filter across the entire input, making it essential for tasks like image processing in neural networks.

In 1D convolution (Conv1D), a one-dimensional input data matrix is convolved with a one-dimensional filter of shorter length, producing a one-dimensional output that highlights local patterns in the input sequence. In 2D convolution (Conv2D), a two-dimensional input data matrix is convolved with a two-dimensional filter of fewer elements, generating a two-dimensional output that captures spatial features and patterns within the image.

Each convolutional block in Figure 2 also includes an activation function. The ReLU activation function filters the data by allowing only positive values to pass through and can be expressed as

$$f\left(x\right)=\mathrm{m}\mathrm{a}\mathrm{x}(0, x)$$

(8)

These Conv1D and ReLU blocks work together to extract and process relevant features from the input data. The incorporation of these components enables the network to handle the complexities of phase correction more efficiently.

Following the encoder, the decoder uses the information provided by the encoder to reconstruct a Lorentzian lineshape of a specific peak in the input signal using the latent space parameters.

2.4.2. Supervised neural network learning

In unsupervised learning, there is no ground truth needed for phase values or absorption spectra without phase errors (Shamaei et al., 2023). The Lorentzian lineshape is used as a reference model for phase adjustments, as shown in Figure 2. When ground truth phase values or absorption spectra without phase errors are available, they can be used for comparison, enabling supervised neural network learning. This approach follows the structure in Figure 1 and the architectural pattern depicted in Figure 2.

Synthetic data are commonly used to establish ground truth for supervised learning. For instance, Jurek (Jurek et al., 2023) utilized the BrainWeb MRI Simulator (https://www.mcgill.ca/bic/neuroinformatics/software-neuroinformatics-tools-and-analysis-brainweb-mri-simulator) with random phase shifts to achieve this goal. The phase correction architecture, based on SRCNN (Super-Resolution Using Deep Convolutional Networks) (C. Dong et al., 2016), includes an input layer, three hidden convolutional layers with non-linear mapping, and an output layer that performs convolutions twice to produce the real and imaginary parts of the phase-corrected image as separate channels. This structure resembles the five layers shown in Figure 1, although the number of nodes differs significantly.

Another example is the FID-A software (Simpson et al., 2017), used in supervised learning by Bugler (Bugler et al., 2024) to obtain ground truth. This approach also involves an input layer, an output layer, and multiple hidden layers. Except for the output layer's convolution, all other convolutional and dense layers use the ReLU activation function, which is independently applied to both the real and imaginary components of the data. A dense layer, or fully connected layer, integrates all information to make decisions or predictions based on identified patterns.

2.4.3. Semi-supervised neural network learning

Semi-supervised neural network learning combines elements of supervised and unsupervised learning to enhance model performance, particularly when ground truth data is limited. This approach is especially useful in the context of phase correction, where it utilizes both labelled and unlabelled data.

An example of this method can be found in Ma et al. (Ma et al., 2024), where the semi-supervised learning framework consists of two components: supervised learning and unsupervised learning. The supervised learning component, similar to the approach used by Bugler et al. (Bugler et al., 2024), relies on simulated data generated with the FID-A software (Simpson et al., 2017). In contrast, the unsupervised learning component uses an in vivo dataset extracted from the Big GABA repository (Mikkelsen et al., 2017), which serves as template spectra. These templates are modified by introducing additional offsets to simulate real-world variations. The results from these two components are integrated to enhance overall model performance.

However, it is important to note that the unsupervised component still involves some supervision, as it utilizes a form of pseudo ground truth. This suggests that the approach is, in fact, supervised learning rather than semi-supervised, as the authors claimed (Ma et al., 2024).

Regarding the architecture in Ma’s study (Ma et al., 2024), the CNN-SR model (convolutional neural network-based spectral registration) is employed for both training processes. Like other neural network models discussed earlier, CNN-SR includes an input layer, an output layer, and several internal hidden layers. Specifically, CNN-SR comprises convolutional layers, batch normalization layers, max-pooling layers, and dense (fully connected) layers. Each hidden layer is followed by a ReLU activation function, while the final output layer uses a linear activation function.

2.5. Supervised neural network combined with high-order polynomial models

In the pursuit of more accurate and robust phase correction, integrating traditional linear models with advanced neural network architectures presents a compelling approach. Neural networks excel at capturing complex, non-linear relationships, while linear models provide simplicity and interpretability. This combination enhances phase correction performance by leveraging the strengths of both techniques.

For instance, You et al. (You et al., 2022) applied a 3D multichannel U-Net CNN in combination with a third-order polynomial regression model. The 3D U-Net architecture, originally introduced by Ronneberger et al. (Ronneberger et al., 2015), features a contracting path on the left for encoding and an expansive path on the right for decoding, forming a U-shaped structure with more than 30 layers. In their study, the U-Net used the hyperbolic tangent (tanh) activation function instead of the more common ReLU. This activation function ensures smooth, bounded transformations by mapping inputs to the range [-1, 1], facilitating the capture of intricate patterns in phase correction. The tanh function is mathematically expressed as

$$\tanh \left(x\right)={\displaystyle \frac{e^x-e^{-x}}{e^x+e^{-x}}}$$

(9)

Ground truth data is essential for training the 3D U-Net, as with all supervised neural networks. In You et al. (You et al., 2022), manually corrected phase values served as the ground truth. After phase correction by the 3D CNN, the refined phase values were further adjusted using a third-order polynomial least-squares regression. This two-step process combined the neural network’s ability to learn complex relationships with the precision of polynomial regression, effectively mitigating spurious local errors and improving accuracy.

Similarly, Srinivas et al. (Srinivas et al., 2023) employed the same 3D U-Net architecture and tanh activation function as You et al. (You et al., 2022). While the specifics of their ground truth data generation were not fully detailed, it appears to follow a similar strategy. Furthermore, the post-processing step involving third-order polynomial regression of the 3D U-Net inferences was incorporated into the architecture diagram.

This section 2 has outlined the main approaches to phase error correction, starting with simple global adjustments and advancing to more complex linear models and neural network-based methods. These models form the foundation of current phase correction strategies and demonstrate the increasing sophistication in addressing phase errors. Their strengths and limitations will be critically evaluated in later sections, particularly in Section 5.1: Challenges in Phase Error Correction Models.

3. Optimization Functions for Phase Correction

Phase errors and the corresponding adjustments are not directly observable. Consequently, all phase correction models—whether traditional statistical or neural network-based—require an optimization process to determine appropriate parameters. This is necessary unless phase errors have already been estimated by other methods, such as manual phase estimation, DISPA (Dispersion vs Absorption) approaches (E. C. Craig & Marshall, 1988; Herring & Phillips, 1984; Hoffman et al., 1992; A. G.; R. Marshall D. C., 1978; Sotak et al., 1984; WACHTER E. A et al., 1989), or the PAMPAS (Phase Angle Measurement from Peak AreaS) method (Džakula, 2000). When phase errors are already estimated, a simple regression model can be applied without further optimization. However, when phase errors are not pre-estimated, this section discusses several optimization functions that can be used for phase error correction.

3.1. Minimization of integral of dispersion spectrum

Consider a simple linear model given by formula (3):

$$\theta ={\beta }_0+{\beta }_{1}x$$

Now, consider formula (2):

$$\theta =-\mathit{arctan}\left({\displaystyle \frac{\int D\left(x\right)dx}{\int A\left(x\right)dx}}\right)$$

Using this, we can formulate an optimization function for phase correction. The following optimization function is conceptually similar to, but presented slightly differently from, the so-called Ernst method described in Binczyk’s paper (Binczyk et al., 2015):

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\int S_D\left(x,{\beta }_0,{\beta }_1\right)dx$$

(10)

Binczyk referred to $S_D\left(x,{\beta }_0, {\beta }_1\right)$ as the magnitude of the dispersion spectrum after phase correction using a linear model with two parameters. Thus, $S_D\left(x,{\beta }_0, {\beta }_1\right)$ can be considered the absolute value of the dispersion spectrum, as it represents the magnitude component only. However, the integral of the absolute dispersion does not necessarily approach zero; in fact, it might actually reach a maximum when there is no phase error. The integral of the absolute dispersion spectrum does reach a minimum when the dispersion peaks are in a Lorentzian shape, which could occur if the dispersion is actually the absorption spectrum or the negative of the absorption spectrum. Therefore, it might be better to consider $D\left(x,{\beta }_0, {\beta }_1\right)$, representing the intensity rather than the absolute intensity, in the dispersion spectrum. The question remains: even with this consideration, is formula (10), used as one of the objective functions in Jaroszewicz’s study (Jaroszewicz et al., 2023), an effective optimization function?

3.2. Minimization of absolute integral of dispersion spectrum

We initially believed that formula (10) represented the optimization function based on Ernst’s paper (Ernst, 1969), as claimed in Binczyk’s paper (Binczyk et al., 2015). However, upon closer examination, this might not be the case. Ernst only stated, “The dispersion mode signal under nonsaturating conditions has a vanishing integral, for slow as well as for arbitrarily fast passage, whereas the integral of the absorption mode signal has a non-zero value” (Ernst, 1969). Clearly, Ernst did not explicitly mention an optimization function; he merely described the characteristics of MR data. Even if we interpret Ernst’s statement as implying optimization, he did not specify whether two individual functions or a combined optimization function should be used.

More importantly, even if we develop an optimization function based on Ernst’s statement—'The dispersion mode signal under nonsaturating conditions has a vanishing integral'—the function should converge to zero rather than merely seeking a minimum. For instance, although -20 is smaller than -2, -2 is closer to zero. On the other hand, if we consider $S_D\left(x,{\beta }_0, {\beta }_1\right)$ in formula (10) as the absolute dispersion spectrum, it diverges from Ernst’s original concept.

Returning to formula (2): $\theta =-\mathit{arctan}\left({\displaystyle \frac{\int D\left(x\right)dx}{\int A\left(x\right)dx}}\right)$, we observe that θ reaches its minimum when the numerator converges to zero, not simply when it reaches a minimum value.

Therefore, the following optimization function, which minimizes the absolute integral of the dispersion spectrum, might make more sense:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}|\int D\left(x,{\beta }_0,{\beta }_1\right)dx|$$

(11)

3.3. Maximization of R ratio

In his original paper, Ernst (Ernst, 1969) did suggest an optimization function: “It is a rather complicated function of the phase angle α, such that an iterative search procedure must be used to determine the proper phase setting for maximum R.” Here, R is defined as the ratio of the “maximum positive signal excursion to maximum negative signal excursion” (Ernst, 1969), and is mathematically expressed as:

$$\widehat{\theta } \underset{\left(\theta \right)}{=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}} \left(R\right)\underset{\left(\theta \right)}{=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}} \left({\displaystyle \frac{\mathrm{m}\mathrm{a}\mathrm{x}\left(A\left(x,\theta \right)\right)}{|\mathrm{m}\mathrm{i}\mathrm{n}(A\left(x,\theta \right)\left)\right|}}\right)$$

(12)

Here, $A\left(x, \theta \right)$ represents the absorption spectrum after phase correction with the value θ, where x is the variable related to frequency or location. This approach aims to find the optimal global phase adjustment for correcting phase errors.

3.4. Minimization of the integral of absolute absorption spectrum

The integral of the absolute absorption spectrum reaches its minimum when there is no phase error. Therefore, an optimization function—the minimization of the integral of the absolute absorption spectrum—can be formulated as:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\int \left|A\left(x,{\beta }_0,{\beta }_1\right)\right|dx$$

(13)

Here, $A\left(x,{\beta }_0, {\beta }_1\right)$ represents the absorption spectrum after phase correction using a linear model with two parameters, where x is the variable related to frequency or location. This optimization function is also referred to as area minimization (de Brouwer, 2009). Although de Brower did not explicitly provide the formula, the concept can be inferred from the MATLAB code in the supplemental file (de Brouwer, 2009).

3.5. Maximization of the integral of the absorption spectrum

Conversely, we can consider another property of the absorption spectrum: the integral of the absorption spectrum reaches its maximum when there is no phase error. This can be mathematically expressed as:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}\int A\left(x,{\beta }_0,{\beta }_1\right)dx$$

(14)

This method was originally suggested by Nelson and Brown (Nelson & Brown, 1989), who also proposed a variation that applies it exclusively to peak regions. It is also the main objective function in Jaroszewicz’s study (Jaroszewicz et al., 2023).

Hardy et al. (Hardy et al., 2009) extended this approach to three-dimensional gradient-echo images by searching for the parameters of a linear model with third-order polynomial terms for phase error correction. In this case, the model involves three position variables and up to 20 parameters: one zero-order term, three first-order terms, six second-order terms, and ten third-order terms. Despite the complexity, the optimization still involves maximizing the integral of the real part of the images, though it is performed in the image domain rather than the frequency domain. Similar to NMR data, MRI raw data are initially in the time domain, and the MRI signals can be analyzed in the frequency domain, commonly referred to as k-space. However, MRI frequency domain data (k-space) can be further transformed into the image domain, where the visual format directly represents the reconstructed spatial information we see in the final images.

3.6. Minimization of the sum of squares of differences between heights at peak-region maxima in the absorption spectrum and the magnitude spectrum

This method was also proposed by Nelson and Brown (Nelson & Brown, 1989). The magnitude spectrum can be calculated using the following equation:

$$M\left(x\right)=\sqrt{{A\left(x\right)}^2+{D\left(x\right)}^2}$$

(15)

Here, x represents the frequency or location-related variable, and M(x), A(x), and D(x) are the magnitude, absorption, and dispersion spectra respectively.

The rationale is that peak heights of absorption and magnitude should be the same when there is no phase error. The optimization function can be written as:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{argmin}{\sum }_{peakmaxima}{(\mathrm{A}\left(x,{\beta }_0,{\beta }_1\right)-\mathrm{M}(x\left)\right)}^2$$

(16)

Note that M(x) is independent of the phase.

3.7. Minimization of the sum of squares of negative values within peak regions in the absorption spectrum

This function aims to minimize the sum of squares of negative values within the peak regions of the absorption spectrum. The rationale is that after phase correction, the negative values within the peak regions should be minimized. This method was also proposed by Nelson and Brown (Nelson & Brown, 1989).

Formula:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\sum }_{peakregions}{\left(\min \left(0,\mathrm{A}\left(x,{\beta }_0,{\beta }_1\right)\right)\right)}^2$$

(17)

3.8. Minimization of entropy with negative peak penalty

Chen (Chen et al., 2002) proposed an optimization function for phase error correction based on entropy minimization. Entropy minimization is a widely used approach in various machine learning applications.

Entropy measures the average uncertainty associated with a set of states and their corresponding probabilities. The general formula for entropy can be written as (Ellerman, 2017):

$$H\left(X\right)=-{\sum }_{x\in X }p\left(x\right)logp\left(x\right)$$

(18)

where p(x) represents the probability of a specific realization x of the random variable X.

For continuous random variables, the concept of entropy is extended to differential entropy, where the summation is replaced by an integral over the probability density function.

While this formula is conceptually straightforward, there is no direct association between MR data and probabilities. To address this, Chen et al. (Chen et al., 2002) proposed replacing p(x) in formula (18) with normalized derivatives of MR data. Additionally, a penalty term is incorporated. The resulting minimization function, based on Chen’s suggestion, is expressed as:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\{-{\sum }_i[{\displaystyle \frac{{|A}_i^m|}{{\sum }_i{|A}_i^m|}}\mathit{ln}\left({\displaystyle \frac{{|A}_i^m|}{{\sum }_i{|A}_i^m|}}\right)]+\mathsf{\gamma }[{\sum }_{i}I\left(A_i<0\right)A_i^2\left]\right\}$$

(19)

Here, $A_i$ a shorthand notation for $A_{(i, {\beta }_0, {\beta }_1)}$, which represents the i-th data point in the absorption spectrum after applying linear phase correction with two parameters β₀​ and β₁. The parameter m denotes the order of the derivative and can range from 1 to 4. The parameter γ is a penalty factor, and I() is an indicator function.

Following Chen’s idea as presented in formula (19), Binczyk proposed a simplified version of this optimization function (Binczyk et al., 2015):

$$\left(\widehat{{\beta }_0}, \widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\{-{\sum }_i\left[{|A}_i\right|\mathit{ln}\left({|A}_i|\right)]+$$

(20)

In this version, the absolute value of the phased absorption is used instead of the normalized absolute value of the derivative of the phased absorption. Additionally, no detailed information is provided about the penalty term.

3.9. Negative peak penalty after normalization

Instead of treating the negative peak as a penalty term in (19) as proposed by Chen (Chen et al., 2002) , de Brower (de Brouwer, 2009) opted to use the penalty component alone as an optimization function, with a slightly modified formula. This can be expressed by the following formula based on the MATLAB code provided in the supplementary file (de Brouwer, 2009):

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\sum }_i{(\left|A_{(i,{\beta }_0,{\beta }_1)}\right|-A_{(i,{\beta }_0,{\beta }_1)})}^2$$

(21)

While de Brower (de Brouwer, 2009) suggested performing normalization before minimization, a constant scaling factor does not impact the minimization process, although it might affect the calculation speed.

3.10. Mean squared error

Mean squared error (MSE) is one of the most commonly used loss functions, with a basic formula that is straightforward and can be expressed as:

$$MSE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2$$

(22)

Here, n represents the number of data points, $y_i$ is the i-th observed value, and ​${\widehat{y}}_i$ is the estimated value for ​$y_i$. Without dividing by n, this formula is known as the sum of squared errors (SSE) or is often referred to as the L2 error. This measure can also be used for phase error correction. For example, it was used as the main component of the loss function to compare the fit and spectra in an extension of the cross-entropy stochastic optimization method (Ravanbakhsh et al., 2010). In this case, the other component of the loss function was the absolute error (Ravanbakhsh et al., 2010).

While the concept of MSE is straightforward, applying it to data with phase errors can be challenging. Phase errors distort the data, making it difficult to compare with data that are free from such errors. To effectively use MSE for phase correction optimization, it is essential to have ground truth data without phase errors. This can be achieved using synthetic data (Jurek et al., 2023; Vanvaals & Vangerwen, 1990), previously acquired reference spectra of pure solid-state components (Prostko et al., 2022), manual phasing of the data (Srinivas et al., 2023; You et al., 2022), or specific mathematical models such as Lorentzian, Gaussian, and Voigt functions (Shamaei et al., 2023).

3.11. Stein's unbiased risk estimate minimization

Since it is difficult to apply MSE without extra ground truth data, Stein's Unbiased Risk Estimate (SURE) was introduced for phase correction optimization process. SURE is an estimate of the mean squared error proposed by Stein (Stein, 1981). It has been applied to MRI data for adaptive phase correction (Pizzolato et al., 2020) and it operates in the image domain. The corresponding loss function is defined as (Pizzolato et al., 2020):

$$\phi \left(I_0,\lambda \right)={\displaystyle \frac{1}{2XY}}{\left|\left|I_0-\rho \left(I_0,\lambda \right)\right|\right|}^2-{\widehat{\sigma }}^2+{\displaystyle \frac{{\widehat{\sigma }}^2}{XY}}({div}_{\mathfrak{R}}\left[\rho \left(I_0,\lambda \right)\right]+{div}_{\mathfrak{I}}\left[\rho \left(I_0,\lambda \right)\right])$$

(23)

Here, the 2D image consists of X data points in the x direction, and Y data points in the y direction, resulting in X×Y pixels. $I_{0 }$is the origin image. $\rho \left(I_0, \lambda \right)$ represents a regularized or phase-corrected version of the image. ${\widehat{\sigma }}^2$ denotes the noise variance. div is a divergence operator. ℜ[⋅] and ℑ[⋅] represent the real and imaginary parts of the function, respectively. The divergence terms can be estimated using a Monte Carlo approach.

3.12. Maximization of refocusing ratio

For 2D image domain data, where phase errors may be related to both the x- and y- directions, the analysis could focus solely on the x direction. As a result, the second-order linear model for phase correction remains the same as in formula (5), excluding any additional terms denoted by “+...”:

$$\theta ={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2$$

Chang and Xiang (Zheng Chang & Qing-San Xiang, 2005) proposed the following optimization function for phase correction in the image domain, which is slightly reformatted here for clarity:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1},\widehat{{\beta }_2}\right)=\underset{({\beta }_0,{\beta }_1,{\beta }_2)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}{\displaystyle \frac{\left|{\sum }_{\left(x,y\right)}S(x,y,{\beta }_0,{\beta }_1,{\beta }_2)\right|}{{\sum }_{\left(x,y\right)}\left|S\left(x,y,{\beta }_0,{\beta }_1,{\beta }_2\right)\right|}}$$

(24)

While this function may seem complex, it is actually quite straightforward. $S(x,y,{\beta }_0,{\beta }_1, {\beta }_2$) refers to the real part of the image data for each pixel after phase correction using a second-order linear model. The ratio of the absolute value of the summation to the summation of absolute values is known as the refocusing ratio (R) (Chang and Xiang, 2005). The R value ranges from 0 to 1, with R=1 indicating no phase error.

3.13. Posterior probability maximization

Similar to section 3.12, this approach also works in the 2D image domain but utilizes a different phase correction model, which follows the same form as in formula (4):

$$\theta ={\beta }_0+{\beta }_{1}x+{\beta }_{2}y$$

Bretthorst (Larry Bretthorst, 2008) proposed simplified posterior probability functions for ${\beta }_1$and ${\beta }_2$​, which take the form of Student's t-distribution. Thus, this t-distribution is used to estimate the first-order parameters. Additionally, the zero-order phase parameter that maximizes the posterior probability is estimated using the following formula (Larry Bretthorst, 2008):

$${\left({\beta }_0\right)}_{est}=-{\displaystyle \frac{1}{2}}{\tan }^{-1}({\displaystyle \frac{X}{W}})$$

(25)

where

$$\begin{gathered} X=2{\sum }_j{\sum }_l{F}_{Rjl}{F}_{Ijl}, \\ W={\sum }_j{\sum }_l[F_{Rjl}^2-F_{Ijl}^2] \end{gathered}$$

Here, $F_{Rjl }$and $F_{Ijl}$ refer to the real and imaginary components of a pixel at the j-th and l-th intersection in the two directions of the image domain, while, ${\sigma }^2$ represents the noise variance.

3.14. Maximum likelihood estimation (MLE) based cost function

To address phase errors induced by motion in 3D imaging, and assuming a linear phase correction model, Van et al. (Van et al., 2011) proposed an MLE-based cost function:

$$\begin{gathered} (\widehat{\mathit{a}},\widehat{a_0})=\underset{(\mathit{a},a_0)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left(R\left(\mathit{a},a_0\right)\right) \\ =\underset{(\mathit{a},a_0)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\sum }_{\mathit{r}}{|e^{j\mathrm{\angle }\widetilde{I_b}\left(\mathit{r}\right)}-e^{j\left[\mathrm{\angle }{I}_0\left(\mathit{r}\right)+\left(\mathit{a}·\mathit{r}\right)+a_0\right)]}|}^2 \end{gathered}$$

(26)

Here, $\mathit{a}$ refers to a vector of the first-order parameters (slopes) of the linear model for the three different dimensions in the image, and $a_0$ is the intercept (zero order or global parameter) of the linear model. The variable $\mathit{r}$ represents image space across multiple dimensions, j = $\sqrt{-1}$ (commonly represented as i in mathematics), ∠ is the phase extraction operator, $\widetilde{I_b}\left(\mathit{r}\right)$ is the image with phase error, and $I_0\left(\mathit{r}\right)$ is the ideal image without distortion.

Although formula (26) appears complicated, it essentially minimizes the sum of squared errors, but in a nonlinear fashion due to the nature of complex image data.

3.15. Maximization of Pearson correlation between magnitude and real spectra

For isolated peaks, an ideal real spectrum should have peak heights that match those of the corresponding magnitude spectrum. The two spectra are positively correlated: as the magnitude spectrum increases, the real spectrum also increases, and as the magnitude spectrum decreases, the real spectrum decreases. Building on this relationship, Wright et al. (Wright et al., 2012) proposed maximizing the Pearson correlation coefficient between the real spectrum and the magnitude spectrum. This is done to optimize zero- and first-order parameters in a linear phase correction model for each of the two specific regions separately: the choline/creatine region (3.36–2.86 ppm) and the citrate region (2.86–2.37 ppm). These regions are essential for phase correction, as they are used as internal references in the application. The method relies on the presence of these metabolites in all samples to ensure effective phase and frequency correction, which is critical for accurate analysis in prostate MR Spectroscopic Imaging (MRSI). The formula can be written as:

$$\left(\widehat{{\beta }_0},\widehat{{\beta }_1}\right)=\underset{({\beta }_0,{\beta }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}{\displaystyle \frac{{\sum }_i(A_{(i,{\beta }_0,{\beta }_1)}-\overline{A_{({\beta }_0,{\beta }_1)}})(M_{(i,{\beta }_0,{\beta }_1)}-\overline{M_{({\beta }_0,{\beta }_1)}})}{\sqrt{{\sum }_i{(A_{(i,{\beta }_0,{\beta }_1)}-\overline{A_{({\beta }_0,{\beta }_1)}})}^2{(M_{\left(i,{\beta }_0,{\beta }_1\right)}-\overline{M_{\left({\beta }_0,{\beta }_1\right)}})}^2}}}$$

(27)

Here, A and M represent the real and magnitude spectra, respectively.

3.16. Tail height minimization and penalty for negative peaks

Bao et al. (Bao et al., 2013) proposed a two-step process for phase error correction: coarse tuning and fine tuning. In the first step, coarse tuning, the following tail height minimization objective function (OF) is used (Bao et al., 2013):

$$OF={\sum }_{i=1}^{PN}abs\{{\sum }_{k=0}^{M}R[peak\left(i\right).Start-k]-{\sum }_{k=0}^{M}R[peak\left(i\right).End+k]\}$$

(28)

Here, R refers to the real part of a spectrum, PN is the number of identified peaks, and M denotes the number of chosen points, where 1 ≤ M ≤ 4.

In the second step, fine tuning is based on minimizing the sum of squared errors, with the penalty function (PF) defined as (Bao et al., 2013):

$$PF={\sum }_{i=1}^N{\{TempR\left[i\right]-abs(TempR\left[i\right]\left)\right\}}^2$$

(29)

Where:

$$TempR\left[i\right]=\left\{\begin{array}{c}R\left(i\right), default situation\\ 0, if the ith point is in distorted peaks\\ -R\left(i\right), if the ith point is in negative peaks\end{array}\right.$$

Here, R(i) indicates the i-th data point in the real part of the spectrum.

3.17. Sum of squared error minimization with scaling

Worley and Powers (Worley & Powers, 2014) proposed an optimization function to simultaneously search for parameters in a linear phase correction model and a scaling factor. The minimization function is derived from Worley and Powers (2014):

$$\left(\widehat{b},\widehat{{\varphi }_0},\widehat{{\varphi }_1}\right)=\underset{(b,{\varphi }_0,{\varphi }_1)}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\sum }_{j=1}^N{|b·s\left({\mathrm{ɷ}}_j\right)e^{i\left({\varphi }_0+{\varphi }_1{\mathrm{ɷ}}_j\right)}-r({\mathrm{ɷ}}_j\left)\right|}^2$$

(30)

In this formula, b is the scaling factor, ϕ₀ and ϕ₁ are the zero-order and first-order parameters for a linear phase correction model, respectively, and i = $\sqrt{-1}$. The index j refers to a point in the spectrum with length N, s(ɷ_j) is the j-th point of the real part of mean-centered row in data matrix X, r(ɷ_j) is the corresponding point in a reference spectrum, and ɷ_j represents the frequency at the j-th point.

If the scaling factor b is ignored, formula (30) reduces to the sum of squared errors (SSE).

3.18. Phase difference minimization

This optimization function is specifically designed for MRI data with multiple scans. Due to the long scan time, these images often contain phase errors caused by patient movement relative to the first scan. To correct these phase errors, a simple optimization function can be applied to the region of interest (ROI) (Xia & Chen, 2021):

$${\varphi }^{\mathrm{\prime }}(k_x,k_y)=\underset{\left(\varphi (k_x,k_y)\right)}{\mathrm{argmin}}(\arg \left[G_j^{\mathrm{\prime }}(k_x,k_y)\right]-\arg \left[G_j(k_x,k_y)\right])$$

(31)

Here, ϕ′(k_x, k_y) is the estimated phase adjustment value at location (k_x,k_y). The term arg[G_j′(k_x, k_y)] and arg[G_j(k_x, k_y) represent the phase values (angles) between two iterations. The initial value of G_j(k_x, k_y) is set as the first scan’s k-space data, including the phase shift due to movement. Through iterative adjustments, G_j(k_x,k_y) is updated until the phase difference is minimized, producing a corrected image.

3.19. Mean absolute error

The mean absolute error (MAE) is a straightforward metric, and its general formula can be expressed as:

$$MAE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n|y_i-{\widehat{y}}_i|$$

(32)

where n represents the number of data points, $y_i$ is the i-th observed value, and ​${\widehat{y}}_i$ is the estimated value for ​$y_i$. Without dividing by n, this formula is known as the sum of absolute errors (SAE) or is often referred to as the L1 error.

Directly measuring the phase shift MAE presents a significant challenge because the phase shift itself is neither known nor directly observable. This limitation is common to all phase correction models, but it is particularly relevant here. Traditional phase error correction methods do not depend on knowing the true value of the phase shift. Instead, they use indirect formulas to estimate phase. To overcome this limitation, Bugler et al. (Bugler et al., 2024) trained a convolutional neural network (CNN) model using simulations where the ground truth of the phase shift was known. This approach relies on the availability of ground truth; without it, the method may not be applicable.

3.20. Linear combination of mean absolute error

The MAE loss function can be extended to a linear combination of multiple MAE loss functions for different variables. Ma et al. (Ma et al., 2024) proposed such a loss function, defined as follows:

$$\begin{gathered} Loss=100\times \left(RealSpectraMAE+ImaginarySpectraMAE\right) \\ +{\displaystyle \frac{FrequencyMAE}{10}}+{\displaystyle \frac{PhaeMAE}{20}} \end{gathered}$$

(33)

The paper reports that the values 100, 10, and 20 in the loss function were empirically chosen to balance error types by weighting frequency and phase MAEs relative to the real and imaginary parts of the spectra.

This approach requires ground truth data to calculate multiple types of MAEs. The real and imaginary MAEs are computed by comparing the template spectra from the BIG GABA database

(Mikkelsen et al., 2017), which serve as approximate ground truth values. The frequency and phase MAEs are calculated using comparisons to simulated ground truth, providing a reference for phase and frequency errors.

This section 3 has reviewed a wide range of optimization functions used for phase error correction, ranging from traditional methods like integral minimization to more advanced techniques, such as entropy-based penalties and maximum likelihood estimation. Their strengths and limitations will be critically evaluated in Section 5.2: Challenges in Phase Error Correction Optimization Functions.

4. Optimizers in Traditional Statistics and Neural Networks

As discussed earlier, phase errors and the corresponding phase values that require adjustment are not directly observable. Therefore, except in rare cases where phase errors have already been estimated, each phase correction model—whether based on traditional statistical methods or neural networks—requires an optimization process to determine the appropriate parameters.

To perform this optimization, both optimization functions and optimizers are needed. While considerable attention has been paid to developing optimization functions tailored to specific challenges, the optimizers that execute these functions often receive less focus, particularly within traditional statistical methods.

An optimization function represents the goal, while an optimizer is the computational procedure used to achieve that goal. Optimizers are algorithms designed to iteratively adjust model parameters to find the optimal solution. The selection of an optimizer plays a crucial role in the optimization process, affecting aspects like convergence speed and the ability to avoid getting stuck in local minima or maxima, and the overall accuracy of the model. In the following sections, we will examine the optimizers used for phase error correction in both traditional statistical models and neural networks.

4.1. False position

The False Position method, also known as Regula Falsi, is an ancient technique that combines elements of the bisection method and linear interpolation to iteratively refine a solution. This optimization method starts with an initial “false” position—a guess that serves as an approximation rather than the exact solution. The process involves evaluating the function at the endpoints of the interval and using these values to linearly interpolate a new point within the interval. This new point becomes one of the endpoints for the next iteration, gradually narrowing the interval and approaching the root. Despite its age, this method has continued to find applications, including the zero-order equality (ZOE) method for phase correction (Balacco, 1994).

4.2. Golden section bisection

The golden section bisection method is an optimization technique designed to find the minimum or maximum of a unimodal function. It builds on the principles of bisection, similar to the false position method, but with a key difference: it utilizes the golden ratio, approximately 1.618, to divide the search interval into segments that maintain a consistent ratio. By leveraging the golden ratio, this method efficiently narrows down the search space and has been used in various applications. For example, in adaptive phase correction (APC), this method was employed to reduce the number of evaluations in the phase correction optimization process (Pizzolato et al., 2020).

4.3. Simplex

The Simplex algorithm is a well-known method for solving linear programming problems, first proposed by George Dantzig in 1947 (Gass, 2011). The algorithm iteratively explores the vertices of the feasible region to identify the optimal solution. More specifically, the Simplex method moves from one vertex of the feasible region to an adjacent one in order to maximize or minimize the objective function. The algorithm has been widely applied in various optimization processes, including phase error correction (Brown et al., 1989; Chen et al., 2002; de Brouwer, 2009).

4.4. Nelder-Mead

Unlike the Simplex algorithm, which is used for linear programming, the Nelder-Mead method (Nelder & Mead, 1965) is a direct search optimization technique that utilizes an n-simplex, which is an n-dimensional polytope with n+1 vertices. This method is particularly useful for optimizing non-differentiable functions. It iteratively adjusts the simplex based on function evaluations at its vertices, moving towards the optimum without requiring gradient information. Although the method itself does not require derivatives, the original paper includes a procedure for estimating the Hessian matrix near the minimum, which can be useful for statistical estimation problems. The Hessian matrix is a square matrix of second-order partial derivatives of a scalar-valued function, which describes the local curvature of the function. The Nelder-Mead method has been cited for its effectiveness in various applications, including phase error correction (Binczyk et al., 2015), and in fact, this is the default method for the widely used R function optim().

4.5. Powell’s method

In Brown’s phase correction paper (Brown et al., 1989), it was mentioned that one of the optimizers was “Powell’s steepest descent.” However, it appears that “Powell’s steepest descent” is not a recognized algorithm; rather, Powell’s method and steepest descent are distinct methods. In this section, we will focus on Powell’s method and discuss the steepest descent method in the following section.

Powell's method was introduced by Powell in 1964 (Powell, 1964). Unlike many modern optimization techniques, Powell's method is a gradient-free algorithm that minimizes functions without the need for derivative calculations (Vassiliadis & Conejeros, 2001). It achieves minimization by iteratively performing line searches along a set of directions, updating these directions based on the results, all without relying on any gradient or derivative information.

4.6. Steepest descent

Remarkably, the steepest descent method, developed by Augustin-Louis Cauchy in 1847, remains widely used today (Andrei, 2022; Petrova & Solov’ev, 1997). Cauchy, renowned for his work in mathematics, introduced this gradient-based technique as a practical approach for optimizing functions. It involves selecting a descent direction based on the negative gradient and finding an optimal step size via line search. Despite its simplicity, it remains fundamental in numerical optimization and continues to influence modern methods.

4.7. Quasi-Newton

The Newton-Raphson algorithm, widely used in traditional statistics, refines model parameters iteratively using both first and second derivatives, offering high accuracy with few iterations due to its quadratic convergence rate. However, computing the Hessian matrix can be costly for large-scale problems.

Quasi-Newton methods address this by approximating the inverse Hessian matrix using only first-order derivatives. Notable methods include the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm and its limited-memory variant, L-BFGS. These methods balance computational efficiency and convergence speed, making them suitable for complex problems where exact Hessian calculations are impractical. Quasi-Newton methods have been applied in phase error correction (Brown et al., 1989; Hardy et al., 2009).

4.8. Hypersphere

Hypersphere optimization techniques constrain the search space within a hypersphere, which is a high-dimensional analogue of a circle or sphere. These methods perform angular searches orthogonal to the direction of steepest descent, reducing the risk of inefficient oscillatory iterations across ridges by stabilizing the search path (Brown et al., 1989). Recent advancements (Kumar et al., 2022) have further refined these techniques, introducing improved spherical search strategies with local distribution-induced self-adaptation, making them effective for non-convex optimization problems, both with and without constraints.

4.9. Levenberg-Marquardt nonlinear least squares

The gradient descent method iteratively adjusts parameters to minimize the objective function by moving in the direction of steepest descent. The Gauss-Newton method approximates the Hessian matrix by using the Jacobian (the partial derivatives of the residuals (or errors)) and ignoring the second-order partial derivatives. The Levenberg-Marquardt algorithm combines the directionality of gradient descent with the efficiency of the Gauss-Newton method, making it particularly effective for nonlinear objective functions. It updates parameters by adjusting the step size based on a damping factor, balancing the trade-off between speed and stability during convergence. This approach is especially useful for problems with complex error surfaces, where other methods might struggle with slow convergence or instability. The Levenberg-Marquardt method was employed in the optimization process for phase error correction (Worley et al., 2015; Worley & Powers, 2014).

4.10. Trust-region-reflective

The trust-region-reflective method iteratively approximates the objective function within a “trust region”, an area where the model is expected to closely match the true objective. The region's size adjusts based on the model's performance: it expands if predictions are accurate and contracts if not. This approach prevents overly aggressive updates while allowing steady progress toward the optimal solution. It's particularly effective for complex, nonlinear objective functions with constraints, making it a reliable choice for challenging optimization problems. The Matlab lsqnonlin() function, which offers this algorithm, was applied in phase error correction (Van et al., 2011).

4.11. Subplex/Sbplx

The Subplex method (Rowan, 1990) was derived from the Nelder-Mead n-simplex method, with Sbplx being a re-implementation of Subplex used in phase error correction (Prostko et al., 2022). Subplex improves upon the standard Nelder-Mead n-simplex approach by breaking the optimization problem into smaller, more manageable subproblems, which can be solved more efficiently. It combines local exploration with global refinement, concentrating on promising regions of the parameter space while remaining adaptable to explore new areas. This blend of local and global search enables Subplex to effectively navigate complex, multimodal objective functions where traditional Nelder-Mead n-simplex methods might become trapped in local optima.

4.12. Stochastic gradient descent (SGD)

PyTorch, a Python library derived from the Lua-based Torch, includes the optimization package torch.optim. The package has been applied in phase error correction (Shamaei et al., 2023), although the specific optimizer used was not mentioned. Among the various optimizers PyTorch offers, SGD and Adam are commonly highlighted in its documentation, alongside options like LBFGS, which we discussed earlier.

Stochastic Gradient Descent (SGD) is a fundamental optimizer for neural networks, particularly effective for large-scale and high-dimensional data. It updates parameters iteratively based on a randomly selected subset of data, making SGD faster and more scalable, but it also introduces noise into the optimization process. Consequently, SGD is commonly combined with the standard backpropagation algorithm (LeCun 1998, Dong 2016) and has been applied in phase error correction (Jurek et al., 2023).

Adam is another method for addressing the noise introduced by SGD and will be discussed in the next section.

4.13. Adaptive moment estimation (Adam)

To address the noise introduced by SGD and improve convergence, several enhancements have been developed, including momentum, learning rate schedules, and adaptive learning rates used in methods like AdaGrad, RMSprop, and Adam.

AdaGrad (Adaptive Gradient Algorithm) and RMSProp (Root Mean Square Propagation) are adaptive optimization algorithms that adjust learning rates based on past gradient information. RMSProp reduces the learning rate for frequently occurring features, while AdaGrad adapts the learning rate individually for each parameter. Adam (Adaptive Moment Estimation) combines the benefits of RMSProp and AdaGrad by using both the exponentially decaying average of past gradients and squared gradients, which enhances convergence speed and stability.

Adam has become one of the most widely used optimizers in deep learning, including for phase error correction (Shamaei et al., 2023; Srinivas et al., 2023; You et al., 2022), due to its robustness and efficiency. It performs well in a variety of settings and often requires minimal hyperparameter tuning to achieve good results.

4.14. Bayesian optimization

Bayesian optimization employs probabilistic models and Bayes’ theorem to strategically determine where to sample next, with the goal of achieving optimal results with as few evaluations as possible. The approach involves building a probabilistic model of the objective function and using this model to guide the search process, effectively balancing the exploration of uncertain areas with the exploitation of known promising regions. Bayesian optimization has been applied in phase error correction (Larry Bretthorst, 2008).

This section 4 has examined the optimizers employed in phase error correction, encompassing both traditional statistical methods, such as the Nelder-Mead algorithm and quasi-Newton, and modern approaches like stochastic gradient descent and adaptive moment estimation. These optimization techniques are essential for fine-tuning models and achieving accurate phase corrections. Their strengths and limitations will be critically evaluated in later sections, particularly in Section 5.3: Challenges in Phase Error Correction Optimizers.

5. Challenges

Despite the wide variety of models, optimization functions, and optimizers developed for MR data phase error correction, not all phase error correction in practice rely on fully automated approaches. A recent review by Ben-Tal et al. (Ben-Tal et al., 2022) even suggested that fully manual phase correction should be used, as it provides more accurate results, highlighting the ongoing limitations of automated approaches. So, why does automatic phase error correction still underperform in certain cases? In the following sections, we explore the key challenges from three perspectives: models, optimization functions, and optimizers.

5.1. Challenges in phase error correction models

To correct phase errors, a model is needed when phase error estimates from other approaches are unavailable. This model—whether based on statistical methods or machine learning—estimates the phase adjustment values required for correction. Once determined, these values are applied to the observed data to adjust the phase and eliminate the errors.

Among the five modeling approaches discussed in Section "2. Phase Correction Models," global phase error correction with a uniform phase adjustment value is the simplest. However, this model is too simplistic to handle phase errors, as phase errors are not unform in real-world MR data. A linear model, which includes both zero-order and first-order parameters, is the most commonly used for phase error correction. It introduces a first-order parameter to adjust for frequency-related phase errors, but it assumes a linear relationship between phase errors and frequency, which cannot adequately address non-linear phase errors.

To manage non-linear phase errors, higher-order terms are introduced into the linear model. While these terms can capture some non-linear phase errors, the non-linear pattern of phase errors cannot be fully addressed by simple quadratic or cubic terms.

Neural networks with multiple layers can handle complex, non-linear patterns, but training often requires ground truth data, typically from simulations, template spectra, or manually phased spectra. Simulation-trained models may struggle with real-world data due to the gap between simulated conditions and actual variability. Template spectra, while standardized, often lack the diversity of real-world signals, leading to poor generalization when models encounter noise, biological variability, or machine-specific differences. Manually phased spectra rely on the expert's skill, and even slight human errors can introduce phase inaccuracies, affecting model quality and consistency. These limitations can weaken the performance and robustness of neural networks trained on such data.

While unsupervised learning does not require ground truth data, it operates under the strong assumption that all peaks follow Lorentzian line shapes (Shamaei et al., 2023), which may not always hold true for real-world data.

Semi-supervised learning, which aims to balance the strengths of both supervised and unsupervised methods, has the potential to overcome these limitations. However, as discussed earlier, Ma’s study (Ma et al., 2024) is not a true example of semi-supervised learning. Instead, it is a supervised approach.

Recent research has added a step involving a linear model with higher-order terms following neural network-based training, making the final model more interpretable and smoothing out potential irregularities from neural network learning. This combination of traditional statistics and modern machine learning has been promising. Although the phased data are not statistically different from manually corrected data, they still fall slightly short, leaving room for improvement, even with sophisticated 3D U-nets (Shamaei et al., 2023; You et al., 2022).

5.2. Challenges in phase error correction optimization functions

Since phase errors and phase adjustment values are not directly observable, developing a model—whether it is a statistical model or a machine learning architecture—requires an optimization function. Section 3, “Optimization Functions for Phase Correction,” discusses 20 optimization functions that have been applied to phase error correction. Here, we would like to explore the challenges associated with these optimization functions.

5.2.1. Integral of the imaginary component

Minimizing the integral of the imaginary component (dispersion spectrum), as shown in formula (10), follows the approach by Binczyk et al. (Binczyk et al., 2015). However, this approach does not fully align with the characteristics of an ideal dispersion spectrum. To address this, we propose minimizing the absolute integral of the dispersion spectrum. The ideal dispersion spectrum, without phase errors, has an integral of zero, with one skewed positive peak and one skewed negative peak. When phase errors are present, the integral deviates from zero, becoming either positive or negative. Therefore, while the integral itself is neither a minimum nor a maximum, minimizing the absolute integral could provide an effective solution (formula 11).

However, minimizing an absolute function poses a significant challenge due to its non-differentiability at zero. The derivative of an absolute value function is undefined where the function crosses zero, complicating the use of gradient-based optimization methods. This issue is particularly problematic when trying to minimize the absolute integral to zero, as the function becomes non-smooth, making it difficult to compute gradients and achieve convergence.

5.2.2. Integral of the real component

Minimizing the integral of the absolute value of the real component (as described in formula (13)) considers only absolute values. While using absolute values may complicate differentiation, it is essential to note that the integral will never converge to zero; instead, it approaches a minimum positive value when the absorption spectrum is ideal. However, it is important to recognize that a reversed negative peak, which is 180 degrees out of phase with the ideal peak, can also reach this same minimal value.

Conversely, maximizing the integral of the absorption spectrum (as indicated in formula (14)) does not involve absolute values, making it a differentiable function. More significantly, since an ideal absorption spectrum devoid of phase errors contains no negative values, the integral naturally attains its maximum when the absorption is optimal. Therefore, formula (14) serves as an effective optimization function.

Additionally, the formula for maximizing the refocusing ratio in the image domain is presented in formula (24). The refocusing ratio is defined as the ratio of the absolute value of the sum of the real component to the sum of the absolute real component. In this context, the numerator corresponds to the absolute value from formula (14), while the denominator refers to formula (13). The objective is to balance maximizing the absolute sum while managing the sum of the absolute values. While this approach appears logical, it can still yield negative intensity values, as an image with a 180-degree phase error can achieve the maximum refocusing ratio. In such cases, the phase shift inverts the real component without altering its absolute intensity.

5.2.3. Peak height

The maximization of the R ratio (formula 12) presents some challenges. In essence, the numerator represents the maximum positive peak height, while the denominator represents the smallest absolute height among the negative peaks in the absorption spectrum. Since an ideal absorption spectrum should have no negative peaks, the denominator would ideally be zero. However, as previously discussed, the derivative at zero introduces complications. Furthermore, because the formula is a ratio, a logarithmic transformation is commonly applied to facilitate its maximization, but log(0) is undefined. Thus, the R ratio, as defined in formula (12), is not a suitable maximization function. Additionally, it only considers two peaks in the entire spectrum.

5.2.4. Entropy

Minimizing entropy with a negative peak penalty, as outlined in formula (19), offers a balanced approach by accounting for all data points while simultaneously controlling negative values. However, since entropy is typically applied to probability functions, there remains the question of whether scaled absolute derivatives of the absorption spectrum can legitimately be treated as probability functions, especially when high-order derivatives are involved.

Formula (20) is a variation of formula (19), also aiming to minimize entropy with a negative peak penalty, but instead of using derivatives, it applies absolute absorption values directly. While formula (19) raises the issue of whether the scaled absolute derivatives can be considered probability functions, they at least fall within the [0, 1] range required for probabilities. In contrast, absolute absorption values do not adhere to this range.

It is noteworthy that the minimization process remains unaffected by constant factors or additive constants. Therefore, even if the values do not strictly satisfy the probability range, minimizing entropy still holds if we assume the scaled absolute absorption spectrum can be considered a probability function. This assumption may be more reasonable than using derivatives, as an ideal absorption peak resembles a Lorentzian function, which can be interpreted as a probability function.

Following this reasoning, if we ignore the scaling factor in the first term of formula (20) while leaving the second term unchanged, the penalty term in formula (20) —which penalizes the sum of squared negative absorption values—loses some of its impact, thereby reducing the overall effectiveness of the penalty.

Nevertheless, both formulas (19) and (20) still raise the question of whether it is reasonable to apply the concept of entropy to absolute absorption values or to the scaled absolute values of their derivatives.

5.2.5. Squared errors

Formula (16) is a variation of the sum of squared errors. The objective is to minimize the difference between absorption peak heights and their corresponding magnitude peak heights, based on the assumption that these heights should be equal. While this method is effective for isolated peaks, it may not be as accurate for situations where the ideal absorption peak height does not correspond with the magnetic peak height due to overlapping. Moreover, this approach tends to overlook most data points that are not local maxima.

Formula (17) involves minimizing the sum of squared negative values in the absorption spectrum. This is based on the premise that an ideal absorption spectrum should not contain negative values, making this formula a variation of the sum of squared errors with negative values considered as errors. Although this method is effective for controlling negative values, it does not account for positive values in the spectrum.

Another related idea is to use the negative peak penalty term from formula (19) on its own to develop an optimization function, resulting in formula (21). Instead of using an indicator, formula (21) directly subtracts the absorption value from its absolute value. If the value is positive, the result is 0; if the value is negative, it results in a positive value, which is considered an error, as an ideal absorption spectrum has no negative values. Looking at formula (21), it essentially sums the squared errors if negative values are treated as errors, making it effectively the same as formula (17). As discussed earlier, while this method is effective at controlling negative values, it does not account for the positive values in the spectrum.

The sum of squared errors can be converted to the mean squared error by dividing it by the number of terms used to calculate the sum, as shown in formula (22), which is a standard statistical measure. The error is measured between the observed value and its fit, estimate, or reference value. Ideally, the reference would be the ground truth; however, in real-world MR data, such a ground truth does not exist. As a result, researchers often rely on synthetic data, template spectra, manually phased data, or specific mathematical models to serve as references. However, as discussed in Section 5.1, these approaches limit the effectiveness of the mean squared error approach.

To address this complicated situation in mean squared errors in MR data without ground truth, Stein's unbiased risk estimate based adaptive phase error correction method was introduced (Pizzolato et al., 2020), and the loss function is shown in formula (23). The main component of the loss function is the sum of squared errors, which measures the difference between the original image and the phase-corrected image. The effectiveness of this approach depends on how closely the original image approximates a phase-error-free image. If the original image deviates significantly from an ideal phase-error-free image, the formula may not perform effectively.

One might wonder why the mean spectrum of real component spectra is not used as a reference for phase error correction, as it often is in regular statistical analyses. The reason is that, unless the spectra are already free of phase errors, their mean will not produce a phase-error-free spectrum—even with a large sample size. This is because most phase errors are systematic biases and do not follow a normal distribution. These errors are influenced by factors such as machine settings, acquisition conditions, or experimental setups, and they persist across multiple measurements. Consequently, averaging spectra does not eliminate these systematic errors in the same way that averaging can mitigate random noise.

The maximum likelihood estimation-based cost function, also referred to as a minimization function, is shown in formula (26). Despite its seemingly complex form, it essentially represents the sum of squared errors. This method aims to correct phase errors by fitting a model to the distorted image data. However, a significant challenge is obtaining an ideal, distortion-free reference image, which is often unavailable in real-world situations.

The second step of Bao’s method (Bao et al., 2013) also employs the sum of squared errors, as shown in formula (29). This approach is quite similar to the one in formula (21), but with a key modification: it minimizes the sum of squared negative values for non-distorted positive peaks, as well as non-distorted negative peaks after reversing them. However, it disregards the distorted peaks altogether. The drawback of this approach is clear—by ignoring the distorted peaks, which should receive more attention, this method is likely to leave residual phase errors even after the phase correction process.

To address the lack of ground truth, Worley and Powers (Worley & Powers, 2014) used the mean spectrum of the dataset as a reference spectrum, combining phase error correction with a scaling step, as described in formula (30). The goal is to minimize the sum of squared errors between the scaled spectrum and the reference spectrum. However, as previously discussed, most phase errors arise from systematic biases rather than random errors. Minderhoud et al. (Minderhoud et al., 2020) demonstrated that phase errors can vary significantly between machines, indicating they are primarily systematic and machine-specific. As a result, when spectra are acquired from the same machine, averaging them cannot eliminate these systematic errors. Therefore, while formula (30) may align each spectrum more closely with the mean spectrum, it is unlikely to correct non-random phase errors.

5.2.6. Bayesian optimization

Through a complex Bayesian calculation, Bretthorst (Larry Bretthorst, 2008) derived a posterior probability function in the form of a t-distribution. This function was used to estimate the two first-order phase correction parameters in the image domain, while the formula for estimating the zero-order phase correction parameter is given in formula (25). The challenge lies in whether all the assumptions made during the derivation process can be satisfied with real-world data.

5.2.7. Pearson correlation

Wright et al. (Wright et al., 2012) proposed maximizing the correlation between the real component spectrum and the magnitude spectrum for two specific peak ranges, as shown in formula (27). While this approach is feasible, it requires that the two specific metabolites be included in all samples. A major limitation of this method is that focusing on just these two ranges may not provide sufficient phase correction for other regions of the spectrum, as much of the data is disregarded. Additionally, an important question arises: does maximizing this correlation truly yield a phase-corrected real component spectrum? It’s important to note that the relationship between the magnitude spectrum ($\sqrt{{Real}^2+{Imaginary}^2}$) and the real component is nonlinear, which raises doubts about whether maximizing the Pearson correlation between them can effectively achieve a phase-free real component spectrum.

5.2.8. Phase difference

A phase difference-based method is generally challenging because a phase error-free real component is not directly observable, making it difficult to establish a true ground truth for phase values. However, in specific cases, such as correcting phase errors caused by patient movement in MRI data, it is possible to minimize the phase differences between the current scan and the previous scan immediately before it, as outlined in formula (31). In MRI, even slight patient movement during image acquisition can introduce phase errors due to shifts in the spatial encoding of signals, impacting the accuracy of the reconstructed images. This approach uses the initial scan as a reference for iterative corrections. However, if the reference scan itself contains significant phase errors due to movement or other factors, these errors could be propagated rather than corrected, potentially limiting the effectiveness of the method.

5.2.9. Absolute errors

Similar to squared errors, absolute errors are generally challenging to apply in phase error correction. However, in simulation-based studies where ground truth phase values are available, minimizing the mean absolute error in phase can be used to train a convolutional neural network. (Bugler et al., 2024). As with squared errors, applying simulation-based training to real-world data may lead to a gap between simulated and actual conditions, potentially limiting its effectiveness for phase correction.

Ma et al. (Ma et al., 2024) extended this approach by combining multiple mean absolute errors across the real spectrum, imaginary spectrum, frequency, and phase, as shown in formula (33). However, this method also relies on simulations and template spectra due to the need for ground truth, and thus faces similar limitations as the Bugler et al. approach (Bugler et al., 2024). Additionally, the weighting coefficients in Ma et al. (Ma et al., 2024) appear somewhat arbitrary, with 2000 times more weight given to the real and imaginary spectra based on template spectra from BIG GABA (Mikkelsen et al., 2017) than to phase based on simulations, thereby reducing the contribution of mean absolute errors in phase correction.

5.3. Challenges in phase error correction optimizers

Optimizers are algorithms designed to adjust model parameters iteratively in order to minimize or maximize a specific objective function. While various optimizers have been applied to phase error correction, each comes with its own set of limitations.

Both the false position and golden section bisection methods can converge under the right conditions with well-behaved functions. However, their assumptions and fixed search behaviors may make them less effective for functions with multiple extrema, as these characteristics can slow convergence or prevent the methods from finding the global optimum.

The Simplex algorithm is effective for linear programming problems but struggles with non-linear or non-smooth problems, as it is specifically designed for linear objectives and feasible regions. It may fail to find optimal solutions when non-linear constraints or objective functions are involved.

The Nelder-Mead method is a robust optimizer for non-differentiable functions and does not require gradient information. However, it can suffer from slow convergence and is highly sensitive to initial conditions, especially in complex or high-dimensional optimization problems.

Powell’s method, another derivative-free optimizer, avoids gradient calculations by performing line searches along iteratively updated directions. It is generally considered robust for smooth, unconstrained optimization but can struggle with non-smooth functions. Conversely, the steepest descent method, which relies on gradients, can be sensitive to the learning rate or step size and often struggles with ill-conditioned problems where the function’s curvature varies significantly.

Quasi-Newton methods face the challenge of memory requirements, as they need to store approximations of the Hessian matrix, which can be demanding for large-scale problems. The computational expense of updating these approximations can also impact efficiency in high-dimensional spaces. While L-BFGS reduces memory usage with a limited-memory approach, it is typically very efficient for large-scale problems, though it may compromise convergence speed compared to full-memory methods.

Hypersphere optimization, involving angular searches, can help avoid inefficient oscillations but can still face issues with local minima and determining how constraints interact with spherical boundaries. Additionally, it may involve computational complexity and scaling challenges, especially in high-dimensional spaces.

The Levenberg-Marquardt algorithm is heavily dependent on the initial guess, which can affect its convergence to local minima. It is also computationally expensive due to iterative updates and Jacobian computations. While it balances speed and stability through a damping factor, careful tuning is required to optimize performance and avoid convergence issues.

The trust-region-reflective method can be computationally demanding due to the need to solve subproblems at each iteration to define and adjust the trust region. The choice of trust region size is critical, as an inappropriate size can lead to slow or inefficient convergence. The method’s computational cost and parameter tuning can impact its performance, especially in complex, nonlinear problems.

Subplex/Sbplx decomposes the problem into smaller, simpler optimization tasks, which can result in increased computational overhead due to multiple local searches. Additionally, managing the trade-off between local refinement and global exploration can impact both the efficiency and convergence of the method, particularly in complex, high-dimensional problems.

Due to the stochastic nature of SGD, updates can be noisy, potentially causing the optimizer to oscillate around minima, preventing it from settling at the optimum. Choosing the right learning rate is challenging: if the learning rate is too high, it can lead to instability, while if it is too low, convergence may become slow. Recent developments often incorporate adaptive learning rate methods or momentum-based approaches to address these sensitivity issues.

While Adam often converges quickly and efficiently, it may sometimes find solutions that generalize less well than other optimizers like SGD. Its generalization depends on the problem and how it handles sharp vs. flat minima. Adam might not always converge to the global minimum, particularly in complex, non-convex optimization problems. It is also important to note that Adam's performance can be sensitive to hyperparameter settings, and its computational and memory requirements can be substantial.

Bayesian optimization is effective for costly evaluations but struggles with model selection and tuning complexity. It can be computationally demanding, especially in high-dimensional spaces or with very noisy functions, and requires careful balancing of exploration and exploitation to optimize performance.

6. Conclusion

In this review, we explored three statistical models, several machine learning architectures, 20 optimization functions, and 14 optimizers used for phase error correction. We also examined the challenges each of these methods presents.

Statistical models remain the dominant approach for phase error correction, with the linear model—utilizing zero and first-order parameters—being the most commonly employed. In some cases, higher-order terms are added to enhance accuracy. Recently, machine learning architectures, particularly convolutional neural networks (CNNs), have gained traction in phase error correction due to their capacity for handling complex patterns in data.

Although a wide range of optimization functions are used in phase error correction, they generally fall into two categories:

MR data-specific characteristics, such as phase difference, peak height, the integral of the real component, and the integral of the imaginary component.

Common statistical approaches, including Pearson correlation, absolute errors, squared errors, entropy, and Bayesian methods.

Among these, squared error is the most widely used component in statistical optimization functions. For MR data-specific optimizations, real component integrals are often favored due to their clearer patterns in phase-error-free signals. However, considering multiple MR-specific features simultaneously could potentially improve accuracy in phase correction.

The main challenge in phase error correction, particularly in optimization functions based on statistical measures like squared error, is the absence of ground truth in real-world data. This is a significant hurdle for both statistical models and machine learning methods, including neural networks. Without a clear ground truth, models may struggle to generalize effectively, especially in diverse, noisy environments.

In terms of optimizers, Nelder-Mead is the most commonly applied for statistical model optimization in phase error correction, largely because it is the default method in the optim() function in R. However, many studies do not explicitly mention the optimizer used. On the machine learning side, Adam is the most frequently utilized optimizer. Both methods, however, come with limitations—Nelder-Mead can be sensitive to local minima, and Adam can struggle with generalization, particularly in non-convex objective functions.

From our own practice, Quasi-Newton optimizers like BFGS and L-BFGS have shown strong performance across both statistical and machine learning models. While memory requirements can be a challenge—particularly for BFGS—this issue has lessened in importance with modern computing advancements, making these methods more viable for large-scale applications.

Given the extensive array of models, optimization functions, and optimizers currently available, is there still a need for further development? The answer remains affirmative. Despite the progress made, significant room for improvement remains in fully automatic phase error correction. As Ben-Tal (2022) pointed out, manual correction is still often recommended for achieving the highest accuracy.

Our ongoing research seeks to depart from the conventional optimization-based approaches to phase error correction. Instead, we are exploring a novel method that does not rely on any optimization process, and thus eliminates the need for optimization functions or optimizers altogether. This approach has the potential to serve as a ground truth to guide machine learning models, especially in the context of phase error correction.

For the future of automatic, AI-based phase error correction, which is seamlessly integrated with other MR data preprocessing steps, the collection of MR data from diverse machines, facilities, and environments is crucial. Optimization functions must account for these variances, and ideally, the process would no longer rely on ground truth data. This can be achieved by leveraging MR-specific characteristics more effectively.

The underlying motivation for this review is our belief in the power of ensemble learning—bridging traditional statistical models and machine learning architectures. By dynamically selecting appropriate optimization functions and optimizers tailored to specific datasets, we believe it is possible to achieve fully automated and highly accurate phase error correction. This is the direction we see for the future.