A user independent denoising method for x‐nuclei MRI and MRS

X‐nuclei (also called non‐proton MRI) MRI and spectroscopy are limited by the intrinsic low SNR as compared to conventional proton imaging. Clinical translation of x‐nuclei examination warrants the need of a robust and versatile tool improving image quality for diagnostic use. In this work, we compare a novel denoising method with fewer inputs to the current state‐of‐the‐art denoising method.


INTRODUCTION
Conventional MRI consists of signal from hydrogen ( 1 H) atoms bound in water and fat, often referred to as proton imaging.The MRI scanner can, however, also acquire signal from other nuclei.These so-called x-nuclei are sensitivity limited because of low endogenous concentration, resulting in an inherently low SNR.][3][4][5][6][7][8][9][10] The lower SNR compared to 1 H is a major limiting factor resulting in longer scan times or low image resolution.The SNR of x-nuclei MRI experiments can, however, be improved in numerous ways.Pre-experiment measures such as coil optimization and isotope enrichment (e.g., hyperpolarization) are used to increase SNR, the latter often being a necessity for 2 H 7 and 13 C 11 imaging.3][14] Last, for post-experiment measures, there exists a variety of processing tools that can increase the SNR.A whole field dedicated to image denoising exists, which uses techniques such as filtering, learning models, statistical methods, or combinations of these, to increase the SNR of MR images. 15In this work, the denoising approach to increasing SNR in x-nuclei acquisitions have been explored in numerical simulations and with in vivo data of human brain and heart.Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA), two novel denoising approaches, have been evaluated for the potential of postprocessing techniques to improve x-nuclei imaging and push it toward clinical application.7][18][19][20][21] Other denoising approaches includes transform denoising (e.g., with wavelets) 22 or convolutional neural network [23][24][25][26] based denoising, although these are less robust and versatile than the statistical approaches for multidimensional data. 27,28GL-HOSVD was developed by Zhang et al. 17 in 2017 and was used on 1 H diffusion-weighted images, and has since been demonstrated for its SNR improvement of dynamic hyperpolarized (HP) 13 C images. 18,29GL-HOSVD applies both a global pre-filtering step and a local patch-based denoising, which is performed on stacks of multidimensional data (tensors), for example, 3 spatial, 1 metabolite, and 1 temporal dimension for dynamic HP- 13 C images.This tensor-based approach uses the multidimensional structure of the data to identify patterns.Local and global HOSVD are both performed using various user inputs within the HOSVD algorithm, which is a potential bias compromising the denoising quality of clinical images having a large variety of noise appearances.tMPPCA was recently introduced by Olesen et al. 20 as a tensor-based extension of the original MPPCA method 21,30 to multidimensional denoising.The tensor-based version was shown to be superior for multi-TE diffusion data.tMPPCA is built on HOSVD using the Marchenko-Pastur distribution to objectively estimate the signal rank (e.g., the signal carrying components to keep after dimensionality reduction) as opposed to the user defined thresholds in GL-HOSVD.This implies that if tMPPCA can denoise the data with a similar performance to HOSVD it can be regarded as an improvement over HOSVD, because of the objective nature of the method.
This study compares GL-HOSVD and tMPPCA denoising on a variety of x-nuclei human data, both with and without previous optimization for HOSVD denoising.This includes 4D HP- 13 C dynamic images of the human heart and brain with metabolites such as pyruvate, lactate, bicarbonate, and alanine, 4D 2 H spectral images of a human brain after ingestion of labeled 2 H glucose, 3D 23 Na images of a human brain, and simulated 13 C dynamic images of the human brain.Comparisons include a visual assessment of the denoising methods as to the original data features, as well as an analysis of residual variance and residual distributions compared to Gaussian reference lines.Furthermore, additional quantitative values are compared by adding Gaussian noise and evaluating denoised images' mean-square-error (MSE) and Bland-Altman analysis to quantify agreement with original data.Previously optimized 18 and user-defined threshold parameters for GL-HOSVD were used to evaluate its robustness when compared to the more user independent approach of tMPPCA.

METHODS
are described elsewhere. 17,20A brief recap of the relevant parts will be presented here.Generally, data corrupted by Gaussian noise can be modeled as the following: where Y represents the noisy observation, A the true underlying signal free of noise and N the zero mean Gaussian noise corrupting the true signal to give the noisy observation.
GL-HOSVD denoising consists of two major steps, the global (prefiltering) and local (patch based) HOSVD denoising.In the HOSVD denoising algorithm, the input data Y is first decomposed with HOSVD into a tensor S. The tensor S is then thresholded against a parameter , and all components below  are nullified.Small components in S mainly reflect the noise N in Y (as well as signal below the noise level).The thresholded tensor is denoted Ŝ.After the thresholding, the inverse HOSVD is applied to the tensor Ŝ, and outputs a tensor Â, which is an estimation of the true signal A. In the global HOSVD, the HOSVD denoising algorithm is applied on the entire multidimensional dataset as a single tensor.The global thresholding parameter is set to where k global is a user defined parameter,  the standard deviation of the noise level, and ∏ R i=1 M i the product of the dimensions of the tensor with rank R. Following the global HOSVD, similar patches of the prefiltered data are stacked and the HOSVD denoising algorithm is applied again on these stacks (so-called local HOSVD), before being aggregated to obtain the final output, a denoised dataset.The local thresholding parameter is set to where k local is a user defined parameter.The stacks of patches are constructed based on a user defined patch size, step size, and radius of search window.By adjusting the parameters k global and k local , the level of global and local HOSVD denoising can be controlled.Similarly, by adjusting the patch size, step size, and radius of search window, the performance and computational time of the local HOSVD denoising of the stacks can be balanced.tMPPCA denoising also uses the HOSVD algorithm, but modifies it in a couple of ways.First, instead of thresholding with a user defined , the number of components that are removed are determined by a set of extra steps.After applying the HOSVD, the individual singular values in the tensor S will predominantly correspond to either signal components or noise components.Note, that there is no singular value with "pure signal components" as noise still enters these, and equally no singular value with "pure noise components" as signal variance below the noise level enters these.However, we will be referring to the two as signal carrying components and noise carrying components.In the MPPCA algorithm, the largest singular values (signal carrying components) will be removed one by one in an iterative process.During each removal step, a consistency test to check whether all the remaining components fit the universal Marchenko-Pastur distribution 31 is performed.Once this consistency criterion is fulfilled, the edge of the Marchenko-Pastur distribution has been identified, and only noise carrying components remain.By running this iterative process, the number of signal-carrying components P are determined, and the remaining noise-only components nullified.This implementation requires no user defined threshold and is, therefore, a more objective approach than the GL-HOSVD algorithm.Second, the tensor MPPCA algorithm includes a recursive step, which speeds up performance and improves on the HOSVD algorithm.The decomposition of the tensor Y with HOSVD amounts to a sequence of regular singular value decompositions (SVDs), each operating on a 2D flattening (matrix) of Y .In the flattenings, all but one dimension of the tensor Y have been pooled together.For the HOSVD algorithm, each regular SVD is processed independently.In the tMPPCA algorithm, after processing each SVD, only the signal carrying components P, determined from the Marchenko-Pastur distribution, is carried over to the next flattening of Y .This means that the first SVD ( Ỹ0 ) will be of size M 1 × M 2 × … M R , the next SVD ( Ỹ1 ) of a reduced size P 1 × M 2 × … M R and so on, until reaching the last dimension R after which the inverse HOSVD is performed.This modification to the HOSVD algorithm will decrease the computational cost for each consecutive SVD and improves each proceeding SVD as only signal carrying components P are carried over.The index ordering of which to perform the SVDs is input in the tMPPCA algorithm along with a patch size, as the algorithm can be performed on subparts of the dataset sequentially or globally on the entire dataset.Note, that when 2D data is input into the denoising algorithm, tMPPCA reduces to original MPPCA.

Simulated ground-truth 13 C brain scans
A simulated dynamic 13 C brain phantom was obtained from a previous study. 18This phantom consisted of a 32 with pyruvate and lactate representing the two metabolites.Kinetic pyruvate-to-lactate conversion rates (k PL ) were defined for gray and white matter as 0.015 and 0.0075 s −1 , respectively.The signal evolution over time of the two metabolites is simulated using a two-site kinetic model. 32

23 Na brain scans
One healthy subject recruited through public announcements for another study 33 was used.The local ethics committee approved the study (no.1-10-72-210-21), registered on ClinicalTrials.gov(no.NCT05215938).
A commercial Helmholtz coil pair (PulseTeq, Chobham, UK) was used for 23 Na imaging.The diameter of the coil loops was 20 cm. 23Na images were acquired with a density-adapted radial acquisition technique, as described by Nagel et al., 34 and zero-filled by a factor of 2 in all three dimensions.The read-out bandwidth was 33.3 kHz.A two-compartment salt-water phantom was included in the scan of the patient.Further acquisition details can be found in Vaeggemose et al. 33 The data analyzed consisted of a 60 × 60 × 60 dimensional dataset (N X × N Y × N Z ) and was B 1 -corrected before denoising.

HP-13 C brain scans (site A)
One healthy subject recruited through public announcements for another study was used.The local ethics committee and the Danish Medicines Agency approved the study (EudraCT, 2020-000352-36).The scanner was equipped with a volume transceiver coil (PulseTeq, Chobham, UK).A spectral-spatial excitation pulse was used with a FA of 12  35 which used the same acquisition parameters.The data analyzed consisted of a 40 , where the metabolites were pyruvate, lactate, and bicarbonate (only 20 × 20 × 6 × 40 on the two latter).

HP-13 C brain scans (site B)
Data from one healthy subject from a previous study 18 was used.The scanner was equipped with a custom-build 32-channel receiver for 13 C acquisition in combination with birdcage RF transmit coils.A spectral-spatial excitation pulse was used with a FA of 20 • on pyruvate and 30 • on metabolites (lactate, bicarbonate).A 2D multi-slice EPI sequence was used with an in-plane resolution of 16 × 16 (FOV = 24 cm) for 8 slices (1.5 cm thickness) and 20 time-points (TR = 3 s, for a total of 60 s).The data analyzed consisted of a 16 , where the metabolites were pyruvate, lactate, and bicarbonate.

HP-13 C heart scans
One cardiac patient recruited through the cardiology department for another study 36 was used.The local ethics committee approved the study (EudraCT, 2018-003533-15).A spectral-spatial pulse excitation was used with a FA of 8 • on pyruvate and 90 • on metabolites (lactate, bicarbonate, and alanine).Readout was a single-shot spiral, 27 ms.FOV was 40 cm, cardiac gated with a pair of pyruvate and metabolite per cycle amounting to 120 excitations (20 timepoints per metabolite).The data analyzed consisted of a 50

Denoising parameters
The GL-HOSVD user defined parameters have previously been optimized for simulated and human dynamic HP- 13 C brain data by Kim et al. 18 In our comparison between denoising methods, the same user defined parameters were applied.The parameters consist of a global threshold of 0.4, a local threshold of 0.8, a patch-size of 5, a step-size of 2, and a search window of 6.Additionally, GL-HOSVD requires an estimation of the noise-level, which is typically estimated through the first or last time-point/Z-slice of the dataset, as they only contain noise.For the HP- 13 C data, the GL-HOSVD denoising was applied for each Z-slice individually, whereas the entire 3D/4D datasets were denoised for 23 Na and 2 H data. tMPPCA is input less in regard to the threshold and only requires a patch-size and a specification of the tensor-structure/flattening order.For all data, a patch-size of X/4 and Y/4 were used, where X and Y are the in-plane spatial matrix dimensions of the MRI data.Previously, the optimal performance was reached at a 10 × 10 patch-size for a numerical phantom, 20 which corresponds approximately to X/4 and Y/4 in most of the data analyzed here.The tensor structure of the HP- 13 C data (including simulations) were chosen such that the three spatial dimensions (X, Y, and Z) were grouped in the first dimension, the temporal information in the second dimension and the metabolites (pyruvate, lactate, bicarbonate, and alanine) in the third dimension of the tensor.The HP-13 C (site A) data was denoised separately on each metabolite because of a difference in matrix size.For the 2 H data, we applied a tensor structure consisting of the spectral data in the first dimension, the two spatial dimensions (X, Y) grouped in the second and the last spatial dimension (Z) along the third dimension.Last, for the 23 Na data each of the three spatial dimensions were retained in the tensor structure.Denoising was performed on the reconstructed complex data with no zero-filling.All plots are of the magnitude data.The spectral 2 H data was fitted with AMARES 37 using prior knowledge of deuterated water (DHO), glucose, Glx, and lactate peaks.

Artificially adding noise to data
Gaussian noise was added to the datasets including real and imaginary components and denoised with the respective methods.Comparison of denoising results from the noise-added data and the original data serve as basis for MSE and Bland-Altman analysis.Evaluation of artifact generation or visual dissimilarities was made from the same denoised imaging results.The noise level of the added noise was chosen to give a visual deterioration of the image quality without being indistinguishable from the original.Noise was added using the awgn MATLAB function which adds white Gaussian noise based on an overall target SNR in dB (assuming the entire dataset to be signal).For the simulated ground-truth HP- 13 C brain data, a whole range of noise-levels (−30-15 dB) were explored in more detail with the use of Monte Carlo simulations with 100 random noise repetitions.

Statistics
The two denoising methods' performances were compared quantitatively using Bland-Altman analysis on the data with added noise.The analysis compares original images against images after GL-HOSVD or tMPPCA denoising (and noise-added images).The reproducibility coefficient (RPC) was computed between the pairs of datasets and was used as a metric to determine agreement with original data.The mean of the differences and 95% confidence intervals are marked on the Bland-Altman plots to indicate variation within data.For the simulated ground-truth HP- 13 C brain data MSE and RPC between denoised and truth were computed at various noise-levels and the voxel-wise variance for a selected slice were plotted (n = 100).The mean voxel variance over the entire dataset were also computed.Residuals were computed and plotted as (Original-Denoising) for each dataset and method.The amount of variance removed and its resemblance to Gaussian noise could be evaluated by looking at the distribution of the residuals.The residual distribution p(r) (real-part of complex data) for each dataset were split into 50 equally spaced bins (histcount in MATLAB) and plotted as log(p(r)) as a function of the squared residual bin value r 2 .This construct would ensure a linear relationship for true Gaussian noise.Furthermore, the log of the residual distribution was normalized to the middle of the distribution (p(r = 0) = 1) and the squared residuals was normalized to the mean variance of the two residual matrices ( 2 matrix ).In addition, any squared residual value above 20 was cut off to avoid any single high-value residual to falsely skew the linear relationship.Gaussian reference lines were fitted to the data, and the slope (variance  2 ) used as another metric of comparison.Assuming that all the variance removed corresponds to noise, a larger slope would correspond to a greater noise removal.

Simulated ground-truth 13 C brain scans
A simulated HP- 13 C brain dataset with a dynamic pyruvate and lactate time-series was analyzed.An example of a slice with lactate over time is shown in Figure 1A, where the "true" data had Gaussian noise added before denoising with GL-HOSVD and tMPPCA.At this noise-level, both denoising methods accurately reconstruct the true data.The overall MSE to "true" are 0.1576, 0.0048 and 0.0047 for noise-added, GL-HOSVD and tMPPCA, respectively (for the entire lactate simulation).tMPPCA is more accurate in the low SNR region of the first three time-steps as it reconstructs the true signal to a better degree than GL-HOSVD (MSE to "true" at time = 3-9 s of 0.1560, 0.0036, and 0.0008 for noise-added, GL-HOSVD, and tMP-PCA, respectively).Figure 1B shows the residuals of the two respective denoising methods.
Additionally, MSE and RPC were computed between denoised and "true" data at various noise-levels, see low-to-mid noise-levels (15-0 dB), whereas tMPPCA pulls slightly ahead at high-to-very-high noise-levels (−5 dB to −15 dB).It has been shown before that low-rank denoising methods introduces non-uniform signal dependent variance, which invalidates SNR as a metric of comparison. 38Figure 2E shows image variance for various noise levels and reveals that the two denoising methods do indeed have a signal dependent variance.However, taking the mean voxel variance, it can be seen that tMPPCA is more consistent at higher noise-levels (see Figure 2F).

Spectral 2 H brain scans
Figure 3 shows the comparison of the two denoising methods on 2 H spectral scans.Both methods better depict the brain outline for all the metabolites in the center most slices (5-7).Similar indications of improvement are also seen along the Z-axis, as metabolites become visible in slice 4 and 8 for the denoised data.GL-HOSVD is removing the noise more thoroughly around the brain, especially for the maximal signal slices (5-7), whereas tMPPCA better reveals signal in the outer slices, especially slice 3 and 4 for glucose, Glx, and lactate and the whole of DHO.Note, that the small artifact, in form of significantly increased signal in one voxel (slice 1), on the GL-HOSVD denoising of DHO originates from the AMARES fitting and not the denoising method.

23 Na brain 3D spatial data
Figure 4A shows the denoising of original data with artificially added noise of 23 Na scans of a healthy brain.GL-HOSVD and tMPPCA performs similarly on this data.Both denoising methods remove the background noise and quite accurately recover the original signal.Some of the slices do, however, have a different signal distribution within the brain compared to the original, see for example slice 27, where the ventricle part is showing an overall more intense area than in the original.Both denoising methods produce a slight artifact in form of the salt-water phantom appearing in slices where it does not appear in the original image, see slice 11 and 35 (marked with white arrow).Note, the brain outline in the residuals (Figure 4B), which indicates a bias in both denoising methods.

HP-13 C brain scans (site A)
The HP- 13 C time-series scan of a healthy brain (site A) can be seen in Figure 5A with original, GL-HOSVD, and tMPPCA denoised data.The applied denoising methods improve visual brain outline depiction throughout most of the time-series.In Figure 5B, additional noise has been added before denoising.Both GL-HOSVD and tMP-PCA remove the added noise quite well for the pyruvate.GL-HOSVD denoising gives a slightly smoother look, whereas tMPPCA denoising leaves a slightly grainier look.For the lactate data, both denoising methods significantly improve on the noise-added images.The level of detail of the denoised images even seems to surpass the denoised image of the original.Background noise is removed more efficiently using tMPPCA denoising as compared to GL-HOSVD, providing an overall cleaner look with better resolved inner most parts of the brain (e.g., Figure 5B, time-step 18 s).A slight spatial bias is observed again in the residuals of the high SNR slices.

HP-13 C brain scans (site B)
The HP- 13 C time-series of a healthy brain (site B) can be seen in Figure 6A with original, GL-HOSVD, and tMP-PCA denoised data.Pyruvate, lactate, and bicarbonate time-series were acquired and are all shown here.The same comparison, excluding the tMPPCA denoising, can be found in Kim et al. 18 Little visual improvement is seen after denoising of the already well-resolved pyruvate.The lactate data sees a similar improvement for both denoising methods, with longer brain visibility over the time-series.Denoising of bicarbonate using tMPPCA makes the brain outline seem clearer in the later time-steps compared to GL-HOSVD, see 50 to 62 s.In Figure 6B, additional noise has been added before denoising of lactate.Both denoising methods seem to surpass, even the original in terms of level of detail, which is especially apparent on the latter time-steps, see 50 to 62 s. tMPPCA results in slightly grainier and more anatomically detailed images, which is especially apparent around time-step 20 s.It is again observed that GL-HOSVD results in a slightly smoother look.This is apparent both before and after adding additional noise, see Figure 6A bicarbonate time-step 20 s and Figure 6B lactate time-step 20 s.As observed in the residuals of pyruvate, GL-HOSVD is removing more noise than tMPPCA, but with a spatial bias occurring in the later time-steps.

HP-13 C heart scans
The HP- 13 C time-series of a healthy heart can be seen in Figure 7 with additional noise added before denoising.At first glance, both methods seem to perform very  similarly for all metabolites.However, as observed earlier, GL-HOSVD leads to a smoother look compared to tMPPCA's grainier look.These seem to be the case for lactate, bicarbonate, and alanine (see time-step 19.8 s for each metabolite).The residuals reveal a slight spatial bias for both methods for all metabolites.

Residual distribution plots
A plot of the log of the residual distribution as a function of the squared residual can be seen in Figure 8 for the simulation, 23 Na brain data, and 13 C brain (site A) data along with Gaussian reference lines.The rest of the datasets can be seen in Figure S1.The shapes of the residual distributions generally follow the Gaussian reference lines.tMPPCA has a slightly higher variance removal for the datasets with lower SNR (e.g., downstream metabolites for HP-13 C data).On the other hand, GL-HOSVD has a slightly higher variance removal for the datasets with higher SNR (e.g., pyruvate for HP-13 C data).Overall, tMP-PCA removes slightly more noise with the mean ratio between the two variances of GL-HOSVD and tMPPCA equating to:  2 tMPPCA ∕ 2 GL−HOSVD = 1.046 ± 0.075.

F I G U R E 5
Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoising applied before and after adding artificial noise to data of HP-13 C brain scans (site A).The first 20 time-steps of pyruvate and lactate data following hyperpolarized [1-13 C]-pyruvate injection in the time-series are shown here (TR = 2 s).(A) Denoising of the original time-series with GL-HOSVD and tMPPCA, with residuals under each individual metabolite (denoted "RES").(B) Denoising after adding artificial noise to the time-series, with residuals under each individual metabolite (denoted "RES").The intensities have been normalized to the maximum signal of the pyruvate data in both cases, but the colorbar of pyruvate adjusted to half the maximum signal to skew focus toward brain tissue rather than only the sinus.

Bland-Altman analysis of denoising performance
To compare the two denoising method's performances with an additional quantitative metric, Bland-Altman analysis was performed on the noise-added datasets.The original data would act as a "ground-truth" for comparison.The results of 23 Na brain, HP-13 C brain lactate (site A) and HP- 13 C brain lactate (site B) can be seen in Figure 8.Additional Bland-Altman analyses for HP- 13 C heart and simulation can be found in Figure S2.The results of Figure 9 shows that GL-HOSVD and tMPPCA have similar RPCs for both 13 C brain lactate datasets (0.116 vs. 0.110 for site A GL-HOSVD and tMPPCA, respectively, and 0.089 vs. 0.085 for site B GL-HOSVD and tMPPCA, respectively) as well as the 23

F I G U R E 7
Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoising applied after adding artificial noise to data of HP-13 C heart scans.The first 20 time-steps of pyruvate, lactate, bicarbonate, and alanine in the time-series are shown here.The sequence repetition time was adjusted to the patient's heartrate, which on average was around 100 bpm (TR ∼ 1.8 s).For each heartbeat a scan of pyruvate and one metabolite (lactate, bicarbonate, or alanine) was performed.Note, that a pyruvate time-series were captured for each metabolite, and the one shown here is the one acquired just before the lactate time-series Each metabolite had artificial Gaussian noise added before denoising, and the residuals are plotted underneath each metabolite (denoted "RES").The intensities have been normalized to the maximum signal of the pyruvate data.Note that the signal from the first six time-steps of the alanine data is an artifact caused by pyruvate flowing in.
23 Na Brain noise-added Simulation 13 C Brain Lactate (Site A) 13 C Brain Pyruvate (Site A) The log of the residual distribution as a function of the squared residual for selected data examples, namely the simulation, sodium ( 23 Na) brain data and carbon ( 13 C) brain (site A) data.The variance  2 for each Gaussian reference is shown and is an indication of the amount of noise removed for the given denoising method (e.g., higher  2 , more noise removed).The log of the residual distribution has been normalized to the middle of the distribution (p(r = 0) = 1) and the squared residuals has been normalized to the mean variance of the two residual matrices ( 2 matrix ).Any squared residual value above 20 has been cut off.

F I G U R E 9
Bland-Altman analysis of sodium ( 23 Na) Brain, hyperpolarized carbon (HP-13 C) brain lactate (site A), and HP- 13 C brain lactate (site B) scans, comparing original with denoising of artificially added noise.The Bland-Altman analysis of the remaining data can be found in Figure S2.The reproducibility coefficient as well as the mean of the difference and 95% confidence intervals are shown on the individual analyses.
GL-HOSVD and tMPPCA, respectively), which is in accordance with the visual analysis.The Bland-Altman analyses found in Figure S2 also shows a similar RPC between GL-HOVSD and tMPPCA for the rest of the datasets (mean RPC improvement over noise-added:  GL−HOSVD RPC = 3.09 ± 1.03 and  tMPPCA RPC = 2.83 ± 0.79).

DISCUSSION
In this study, the denoising performance of the two statistical denoising algorithms GL-HOSVD and tMPPCA was compared on a set of different x-nuclei and simulated MRI data.Generally, both methods had a similar performance for the various data presented here, with a few examples of superior denoising with either method.The denoising of the noise-added 23 Na and 13 C images generally resulted in similar overall RPC values (Figure 9 and S2), which reflects the similarity of the denoised noise-added images to the originals.It is, however, to be noted that any noise present in the original image ought to be subject to denoising as well, making it difficult to directly compare the original and the respective denoised images.An example of this can be seen in Figure 5B, where the denoising of the noise-added images results in images that have higher SNR and are better resolved than the original.Only the simulated data would result in a fully accurate quantitative assessment.The MSE and RPC values of tMPPCA denoising showed superiority for the simulated high-to-very-high noise-level examples (Figure 2B,C).Additionally, tMP-PCA denoising of the first few images of the simulated dynamic data (Figure 1, time-step 1-3) did also result in clearer detail than for GL-HOSVD denoising, which was clearly reflected in the specific image MSE values.It was also observed in the residual distribution plots (Figure 8 and S1) that tMPPCA generally had a higher variance removal in the low SNR/high noise data, which is often what is most desirable.GL-HOSVD generally resulted in smoother images and tMPPCA in more grainier images.This mostly seems to be in favor of tMPPCA, as the fine structured anatomical detail often was better resolved than for GL-HOSVD.Although it is to be noted that this apparent finer detail also could originate from unremoved noise.This was especially apparent for the 13 C data, which might be a consequence of the way the two denoising algorithms deal with the additional temporal dimension.Patterns in the residual plots also reveal some incomplete and/or faulty denoising for both algorithms, which might be because of a breakage of the two algorithms respective assumptions about the data.However, as it seems to be mostly incomplete noise removal in the region-of-interest, the denoised images should be trustworthy.The only exception is of the GL-HOSVD denoising of pyruvate in Figure 6A, where some artifacts seem to be generated in the last time-steps.Note, that despite the sharper appearance of the AMARES fitted denoised 2 H spectral images, a bias and underestimation of the metabolite levels might be introduced after denoising, as pointed out by Clarke and Chiew. 38This is something to keep in mind when applying tools like these for spectral (e.g., 2 H data) or kinetic (e.g., HP-13 C data) fitting.
6][47] The two denoising approaches compared in this work both surpass previous generations of HOSVD based denoising approaches such as global HOSVD, 48 local HOSVD, 16 and tensor image enhancement 39 for GL-HOSVD, 17,18,29 and MPPCA 30 for tMPPCA. 20Hence, they seem like the adequate candidates for comparison of state-of-the-art tensor denoising of multidimensional x-nuclei data.One of the major distinctions between the two methods is the signal rank estimation and reliance on inputs in the algorithm.tMPPCA would generally be assumed to be more versatile than GL-HOSVD when expanding the palette of various data to denoise, because the latter relies on several input parameters and has a larger parameter space to optimize.In terms of using denoising as an automatic postprocessing application to push x-nuclei MRI closer to a clinical setting, it seems from this standpoint that tMPPCA would be the best choice.In addition to this, the computational time of tMPPCA denoising was generally much lower than GL-HOSVD denoising with the current parameters (10-40 times faster mostly).This is likely because of the recursive patch-based nature of tMPPCA, where the dimensionality of the data decreases for each successive SVD, which has also been discussed elsewhere. 20All in all, tMPPCA has shown similar denoising performance to GL-HOSVD, with a slight advantage at lower SNR, and is, therefore, deemed the preferred denoising method for a user independent denoising of x-nuclei MRI and MRS data of various structures.

CONCLUSION
In conclusion, we compared two state-of-the-art denoising algorithms, namely GL-HOSVD and tMPPCA, on a range of x-nuclei data including simulated data.tMPPCA and GL-HOSVD had similar performances for the variety of data analyzed here.As tMPPCA is a more versatile and objective method, because of fewer inputs, it is deemed the preferred method for a user independent denoising application in x-nuclei imaging.
denoising method (e.g., higher  2 , more noise removed).The log of the residual distribution has been normalized to the middle of the distribution (p(r = 0) = 1) and the squared residuals has been normalized to the mean variance of the two residual matrices ( 2 matrix ).Any squared residual value above 20 has been cut off.Note that the residual distribution of the 2 H brain data is of the entire spectral dataset, and not the AMARES fitted residuals shown in Figure 3 in the main article.Figure S2.Bland-Altman analysis of HP-13 C heart data, and 13 C lactate simulations, comparing original with denoising of artificially added noise.The Reproducibility Coefficient (RPC) as well as the mean of the difference and 95% confidence intervals are shown on the individual analyses.
Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoising applied after adding artificial noise to simulated ground truth data.(A) The simulated data shown here represents an axial slice of lactate in a brain over time in a typical hyperpolarized carbon pyruvate-based experiment.The noise-level of the artificially added noise amounts to an overall SNR of 10.The scale of the intensity has been cutoff at 1.2 times the maximum signal of the "true" data.The overall mean-square-error to "true" are 0.1576, 0.0048, and 0.0047 for noise-added, GL-HOSVD, and tMPPCA, respectively.(B) Residuals of GL-HOSVD and tMPPCA.

Figure 2 . 2
Figure 2.An axial slice (timepoint = 24 s) of simulated lactate is shown in Figure 2A. Figure 2B,C shows MSE and RPC as a function of the SNR (noise-level in dB) for the entire dataset of all temporal brain slices.Figure 2D shows noise-levels ranging from −15 to 15 dB, with reference in Figure 2A.MSE and RPC reveal a very similar performance between the two denoising methods at

F I G U R E 3
Comparison between AMARES fit of original, Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoised spectral deuterium ( 2 H) brain data.The AMARES fit uses prior knowledge of the four components observed in the spectra, namely deuterated water (DHO), glucose, glutamate/glutamine (Glx) and lactate, the latter three originating from the downstream metabolism of ingested [6,6 ′ -2 H 2 ]-glucose.The intensities have been normalized to the original data, and the residuals plotted underneath each respective metabolite (denoted "RES").The Z-slice direction represents the axial direction of the brain.

4
Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoising applied after adding artificial noise to data of 23 Na brain scans.(A) Every second axial slice from 5 to 43 are shown here.A two-compartment salt-water phantom was included in the scan, which can be seen as the two circles appearing in the top of the image throughout slice 13-39.The intensities have been normalized to the maximum signal in all the data.The white arrows indicate artifacts of the denoising algorithms, as no phantom is present in the given z-slices.(B) Residuals of every fourth axial slice from 5 to 43.
Global-local higher order singular value decomposition (GL-HOSVD) and tensor Marchenko-Pastur principal component analysis (tMPPCA) denoising applied before and after adding artificial noise to data of HP-13 C brain scans (site B).The first 20 time-steps of pyruvate, lactate and bicarbonate data following hyperpolarized [1-13 C]-pyruvate injection in the time-series are shown here (TR = 3 s).(A) Denoising of the original time-series with GL-HOSVD and tMPPCA, with residuals under each individual metabolite (denoted "RES").The intensities have been normalized to 1/5th of the maximum signal of the pyruvate data, to skew focus more toward brain tissue rather than only the sinus.(B) Denoising of lactate data after adding artificial noise to the time-series, with residuals underneath (denoted "RES").The intensities have been normalized to the maximum signal of the lactate data.
Na brain data (0.033 vs. 0 • on pyruvate and 70 • on metabolites (lactate, bicarbonate).Readout was a single-shot spiral, 21 ms for pyruvate and 25 ms for metabolites.Matrix size was 40 × 40 on pyruvate and 20 × 20 on metabolites.FOV was 35 cm, TR = 500 ms, with 40 timepoints per metabolite (120 excitations in total).Further acquisition details can be found inBøgh et al.,