Algebraic methods and computational strategies for pseudoinverse-based MR image reconstruction (Pinv-Recon)

Abstract

Image reconstruction in Magnetic Resonance Imaging (MRI) is fundamentally a linear inverse problem, such that the image can be recovered via explicit pseudoinversion of the encoding matrix by solving —a method referred to here as Pinv-Recon. While the benefits of this approach were acknowledged in early studies, the field has historically favored fast Fourier transforms (FFT) and iterative techniques due to perceived computational limitations of the pseudoinversion approach. This work revisits Pinv-Recon in the context of modern hardware, software, and optimized linear algebra routines. We compare various matrix inversion strategies, assess regularization effects, and demonstrate incorporation of advanced encoding physics into a unified reconstruction framework. While hardware advances have already significantly reduced computation time compared to earlier studies, our work further demonstrates that leveraging Cholesky decomposition leads to a two-order-of-magnitude improvement in computational efficiency over previous Singular Value Decomposition-based implementations. Moreover, we demonstrate the versatility of Pinv-Recon on diverse in vivo datasets encompassing a range of encoding schemes, starting with low- to medium-resolution functional and metabolic imaging and extending to high-resolution cases. Our findings establish Pinv-Recon as a versatile and robust reconstruction framework that aligns with the increasing emphasis on open-source and reproducible MRI research.

Introduction

Magnetic resonance imaging (MRI) is a cornerstone of modern medical diagnostics, providing unparalleled soft tissue contrast and functional information. Its widespread adoption, with over 100 million scans performed annually¹, demonstrates its clinical importance. Central to the effectiveness of MRI is the process of image reconstruction^2,3, which converts the raw k-space data into diagnostic images. MRI reconstruction has traditionally relied on the fast Fourier transform (FFT) to rapidly transform uniformly sampled Cartesian k-space data into images. However, FFT-based methods struggle with non-Cartesian sampling patterns, parallel imaging (multi-coil data), and other complex encoding schemes, necessitating more sophisticated reconstruction algorithms^{4, 5–6}. As a result, iterative reconstruction techniques and other advanced methods have been developed to address these challenges, improving image quality while managing computational complexity.

From a mathematical perspective, MR image reconstruction can be understood as a linear inverse problem, where the unknown image is encoded into the measured k-space data via an encoding matrix, ⁷. Direct inversion of this encoding matrix, while theoretically possible, has historically been deemed computationally prohibitive in terms of both memory requirements and computational time, because the complexity scales with , where N is the number of elements on one side of Encode^{5,8, 9–10}. This is in contrast to the computational complexity of of the FFT. Early attempts, for example Van de Walle et al. in 2000¹¹, explored pseudoinverse-based reconstruction but were constrained by the computing resources of that era. Pruessmann’s seminal 2001 paper on sensitivity encoding with arbitrary k-space trajectories also illustrated the challenges of direct matrix inversion on early-2000s hardware and advocated for more efficient iterative solutions⁷. Since then, few studies have revisited the computational burden of direct matrix inversion in MRI.

Recent advances in computational hardware and software have significantly altered the landscape of numerical linear algebra computations. The transition to 64-bit operating systems in the early 2000s enabled the handling of very large arrays (exceeding the 4 GB memory limit of 32-bit systems)¹², and compatible software packages have emerged in the subsequent years. Progress in multi-threaded central processing units (CPUs), parallel computing platforms, powerful graphical processing units (GPUs), and tensor processing units (TPUs) have driven the rise in artificial intelligence and deep learning, which in turn have prompted further developments in these areas. Combined with highly optimized numerical linear algebra libraries, these have made large-scale matrix operations, such as direct pseudoinversion of the encoding matrix, increasingly practical for MR image reconstruction applications.

In this work, we investigate algebraic methods and computational strategies for direct pseudoinversion of the MR encoding matrix, an approach we term ‘Pinv-Recon’. By leveraging modern computational resources, we demonstrate that Pinv-Recon can be a viable alternative to conventional reconstruction methods, especially for low- to medium-resolution imaging applications. This approach offers several benefits, including simplicity of implementation, inherently interpretable regularization, the ability to account for multiple intertwined encoding mechanisms and distortions in a single framework, and support for computing relevant image reconstruction metrics, such as the spatial response function (SRF) and the noise matrix.

The remainder of this paper is organized as follows. Algebraic Background covers potential matrix decomposition methods for directly computing the pseudoinverse of the encoding matrix for MR image reconstruction, expanding beyond the previously suggested method of singular value decomposition (SVD). Computational Evaluation and Numerical Simulations evaluates Pinv-Recon in terms of computational feasibility, the impact of different regularization approaches, and its versatility in incorporating various encoding mechanisms and distortions. Exemplary Applications include Pinv-Recon of a range of 2D and 3D imaging experiments, including hyperpolarized Carbon-13 brain, hyperpolarized Xenon-129 lung, proton knee, proton brain (the ISMRM CG-SENSE challenge dataset¹³), and proton abdominal imaging. The final section provides a discussion of the advantages, limitations, and potential future directions of Pinv-Recon.

Algebraic background

Inspired by the notation introduced by Pruessmann et al. in Encoding and reconstruction in parallel MRI¹⁴, we first formalize the MR image reconstruction problem and then outline possible solutions, including various methods for pseudoinversion of the encoding matrix. We use descriptive variable names in equations for clarity and readability of the algebraic formulation. MRI can be modeled as a linear forward mapping, or “encoding”, of the unknown image onto the raw k-space data, mediated by an encoding matrix, :

1

where both and have been flattened into 1D vectors. is a dense, complex matrix that relates each measured data point to all unknown image voxels. In its simplest form, Encode considers only basic gradient encoding (Fourier encoding) and reduces to a Fourier matrix, with each element written as:

2

where is the imaginary unit, is the th k-space sampling point, and is the th spatial location (voxel position). The k-space points, and spatial locations, can both be discretized to arbitrary grids, including Cartesian, non-Cartesian, uniform, non-uniform, overdiscretized, and masked ones (See Supplementary Material 1).

MRI reconstruction can be considered the inverse problem of this forward mapping, i.e., solving for the unknown image given the measured data and the known matrix that encapsulates all the encoding physics and geometry. Image reconstruction can be performed by first explicitly inverting Encode to find the reconstruction matrix ( and then applying it to the data by a single matrix-vector multiplication:

3

This assumes that Encode is full rank and thus invertible, or has been regularized so that a pseudoinverse (i.e. Moore-Penrose pseudoinverse, denoted by ) exists. Alternatively, the image can be obtained without explicitly forming , for example via iterative algorithms.

The following sections highlight several commonly used matrix pseudoinversion methods and linear system solvers. Among these, only truncated SVD-based pseudoinversion and iterative solvers have been commonly considered in the context of MR reconstruction prior to this work^8,9,11,15.

Direct inversion methods

The true inverse (denoted by ) of a matrix exists only if the matrix is square and of full rank. However, is typically rectangular (total number of measurements by total number of unknowns ) in most MR applications. Direct inversion methods decompose the full rectangular encoding matrix to compute the Moore–Penrose pseudoinverse, which is a generalized solution for inverting such rectangular matrices¹⁶. We highlight two common approaches below.

Singular value decomposition: Singular value decomposition (SVD)¹⁷ factorizes and inverts the encoding matrix as:

4

where and are unitary matrices (of dimensions and , respectively) and is a diagonal matrix containing the singular values of Encode. The set of singular values (the SVD spectrum) provides a useful visual impression of the overall conditioning of Encode. are the truncated versions of each matrix respectively. Truncated SVD (tSVD) regularizes and ensures stability in the inversion by removing small singular values. The ratio of the largest and the remaining smallest singular value in is then referred to as the condition number ( ). Typically, the truncation threshold is set to capture a predetermined percentage of the energy of the original data, such as 95%. Note that reconstruction matrices with different condition numbers can be realized from the same original untruncated SVD decomposition.

QR decomposition: QR decomposition¹⁸ factors a matrix into orthogonal and triangular factors and provides the pseudoinverse as:

5

where is an matrix with orthonormal columns and is an upper triangular matrix. Q and R may also be truncated for regularization. QR decomposition is typically faster than SVD and maintains good numerical stability, but does not provide the information of the SVD spectrum.

Normal-equation-based methods

Instead of decomposing the full rectangular Encode, the normal equations provide a method to calculate the Moore-Penrose pseudoinverse for both over- and under- determined problems:

• Overdetermined problems ( ): where the left-pseudoinverse solution gives providing the least-squares solution for the image (i.e., minimizing ), and

• Underdetermined problems ( ): where the right-pseudoinverse solution gives providing the minimum-norm solution (i.e., ).

The two cases can be derived from one another by simple Hermitian transpose operations; for simplicity, we focus on the more typical overdetermined case going forward. In this case, the term is referred to as the Gram matrix, which is Hermitian and positive semi-definite.

Normal-equation-based methods avoid directly inverting the raw rectangular matrix. Instead, they decompose the smaller, square Gram matrix to derive the pseudoinverse, which is computationally more efficient than decomposing the full rectangular matrix. However, it can be problematic for ill-conditioned problems because the condition number of is squared compared to that of . This issue can be mitigated by applying appropriate Tikhonov regularization^19,20, which improves the conditioning of the matrix, similar to singular value truncation. Additionally, Tikhonov regularization helps ensure positive-definiteness, a required property for certain matrix decomposition methods. It is implemented by adding a small to the diagonal of the Gram matrix:

6

For very large , the solution approaches a simple conjugate-phase reconstruction²¹, whereas smaller values more fully account for density compensation and deapodization. The choice of can be guided by principles like Morozov’s discrepancy principle, which suggests selecting such that the residual is on the order of the noise level in the data²².

Eigenvalue decomposition: Eigenvalue decomposition¹⁷ applies to the Gram matrix and gives Recon using the normal equation:

7

where the columns of are eigenvectors and is a diagonal matrix of eigenvalues. Note that for a Hermitian matrix. Eigenvalue decomposition is a special case of the SVD applied to a Hermitian matrix (or a symmetric matrix in the real case), thus resulting in only one base P. The eigenvalues of the Gram matrix are equivalent to the singular values of Encode squared. Regularization can be applied by truncating the eigenvalues, by adding to the eigenvalues for Tikhonov regularization, or by combining both.

Cholesky decomposition: Cholesky decomposition¹⁷ is known for being one of the most computationally efficient matrix decomposition algorithms, serving as an attractive choice for Pinv-Recon. It is a specialized method for square, Hermitian positive-definite matrices, factorizing the Gram matrix and giving Recon as:

8

where L is a lower triangular matrix. The inversion of L is computationally efficient because of its triangular structure.

Block-wise matrix inversion

When working with large encoding matrices that exceed available GPU memory, or when parallel computing resources can be exploited, block-wise matrix decomposition can be used¹⁷. This involves partitioning Encode or the Gram matrix into smaller blocks, and then applying the algorithm (e.g., factorization, inversion) at the block level, sequentially or in parallel. Block-wise methods are already commonly used internally in linear algebra libraries for performance. In this work, we use the block-wise Cholesky decomposition to overcome GPU memory limitations²³. Note this involves some copying overhead.

Direct solution without explicit inversion

The matrix inversion methods described above explicitly calculate the pseudoinverse of Encode to obtain the reconstruction operator Recon, and we will refer to them collectively as Pinv-Recon. Other approaches exist which directly obtain image without calculating Recon. For example, directly solving equation (1) with simple matrix division of data by Encode, adding Tikhonov regularization to ensure numerical stability. In practice, this can be implemented using a linear solver (e.g., MATLAB’s backslash operator ), which implicitly performs a matrix factorization (such as QR or Cholesky-based factorization) to find the least-squares solution without explicitly returning the pseudoinverse.

Iterative algebraic reconstruction techniques, such as conjugate-gradient (CG) methods, may also utilize the explicit encoding matrix and solve for image without finding Recon, and have been widely explored due to their computational efficiency¹⁷. For example, it has been explored by Kadah and Hu in 1998⁹, by Wilm et al. in 2011²⁴, and by Li et al. in 2015²⁵, incorporating various encoding mechanisms such as distortion and coil sensitivity which will be further discussed in Computational Evaluation and Numerical Simulations - Generalized Encoding. We will hereafter refer to this approach as . Alternative iterative algorithms also exist, for example LSQR²⁶.

Image metrics derived from Recon

Recon, which is generally not calculated in modern MR image reconstruction algorithms or in direct solution methods such as , contains useful information for the calculation of relevant image metrics^14,27. The spatial response function (SRF) measures the proportion of a pixel represented in the reconstructed image. It effectively calculates the forward encoding step followed by the backward projection, and accounts for the fact that each image value does not exclusively reflect the magnetization at the corresponding position, but also exhibits contamination from surrounding signals. Using Recon and Encode, the SRF matrix (SRF) of all image pixels can be calculated by simple matrix-matrix multiplication:

9

In the optimal scenario, SRF approaches identity. An SRF map can be obtained by reshaping the diagonal elements of SRF to the image size, where a value of 1 indicates the signal from that pixel is fully represented in the final image. Recon also allows simple calculation of the noise matrix X, which reflects the noise propagated during reconstruction:

10

where is the noise covariance matrix.

Computational evaluation and numerical simulations

Computational time

Initial work was undertaken to explore the memory requirements and computation times for decomposing reconstruction matrices using different methods. Random complex single-precision matrices were generated to simulate encoding matrices for 2D image matrix sizes (MTX) 32 32 up to 320 320 (with equivalent MTX for 1D and 3D images noted in the x-axis of Fig. 1b). Figure 1a illustrates Pinv-Recon of the Shepp-Logan phantom for generic 1D, 2D and 3D random sampling. The matrix decomposition times using SVD, QR, eigenvalue, and Cholesky decomposition in MATLAB (The Mathworks, MA, USA) were measured on a high-performance workstation (Intel® Xeon® Gold 6448 CPU, 2 32 cores, 1 TB RAM; NVIDIA® A100 80 GB GPU), similar to those available on modern MR scanners. The results are summarized in Fig. 1b,c. All computations in MATLAB by default utilized multi-threaded CPU processing across the available CPU cores with the LAPACK (Linear Algebra Package) library²⁸, providing a strong baseline for performance comparisons with GPU acceleration. GPU acceleration was achieved by assigning Encode to the GPU using gpuArray in MATLAB, the timing of which was not included in the results but can be considered negligible relative to the computation time. Computational times required on a less powerful mobile workstation are also shown in Supplementary Material 2.

Fig. 1 [Images not available. See PDF.]

(a) Illustrations of Pinv-Recon using generic 1D, 2D and 3D trajectories. (b) Plot of computational times for SVD, computing Encode Encode, and Cholesky decomposition for a range of matrix sizes using CPU only (solid lines) and with GPU acceleration (dotted lines) (c) Corresponding table of computational times also including QR and Eigenvalue decomposition. Asterisks (*) indicate decomposition using block-wise Cholesky decomposition.

The computational time required for a 2 GB encoding matrix (approximately equivalent to a MTX 128 128 2D image with full Fourier encoding) ranged from under a second to several minutes, depending on the inversion method used. Among the four methods, Cholesky decomposition was fastest, followed by eigenvalue decomposition, QR, and then SVD. This ranking is expected from their algorithmic cost. All methods exhibited scaling approximately proportional to the cube of the matrix dimension ( ). Use of the GPU provided additional speedups for the larger matrices (acceleration ranging from 1.7- to 7.8-fold for matrices above 2 GB), though the benefit varied with the decomposition method and matrix size. For very small matrices, the GPU-associated overheads and latencies dominate, and a CPU implementation is more efficient. A block-wise Cholesky algorithm was employed to invert matrices larger than the GPU memory (indicated by asterisks in Fig. 1c). While it accelerates computation relative to the CPU, the associated overhead reduces the speedup compared to a non-blockwise approach.

A few additional points are worth noting. For normal-equation-based approaches, the formation of the Gram matrix (computing ) determines the lower bound of the computational time, and therefore is included in Fig. 1b,c for reference. The computational time of iterative solvers is largely determined by factors that influence convergence, such as the matrix’s condition number and the chosen stopping criteria. While these methods can avoid expensive matrix–matrix operations by relying primarily on matrix–vector multiplications, their runtime still grows with the number of iterations required. For conjugate gradient methods that incorporate the FFT, the complexity scales as , offering significantly greater speed-ups for large matrices compared to smaller ones. The decomposition times for non-Cartesian trajectories (circular coverage of k-space) or masked images (spatial restriction to object coverage) will be shorter than images of equivalent MTX listed in Fig. 1c, which assume equal and , i.e. full Fourier encoding, for simplicity. For example, a MTX 128 128 image sampled with 15776-point spiral, masked to a circular field-of-view (covering 78% of image space) results in an encoding matrix of only = 15776 12780, and a Gram matrix of 12780 12780, requiring much shorter decomposition time than for = 16384 16384 (128 128 ). The present computational evaluation assumes single precision. If double precision (64-bit) were used instead, memory usage and computation time will increase, approximately doubling in many cases; conversely, emerging half-precision (16-bit) or mixed-precision arithmetic could further accelerate these computations at the cost of reduced numerical stability.

Regularization

Regularization is an essential step in pseudoinverse reconstruction, determining the trade-off between noise amplification and resolution. We compared tSVD versus Tikhonov regularization using a Shepp–Logan phantom simulation. A 128 128 Shepp–Logan phantom was forward-encoded with an 8-arm spiral trajectory (Fig. 2a) and reconstructed using Pinv-Recon under varying noise and regularization settings. The Shepp-Logan images had pixel intensities ranging from 0 to 1, and random complex Gaussian noise was added in k-space at noise levels of 0, 0.1*noise, 1*noise, and 10*noise (where noise=randn(size(Data))+1i*randn(size(Data))/2). These noise levels were chosen so as to be visible in the image. tSVD reconstructions were performed while varying the energy cutoff (from minimal to near-total truncation), and Tikhonov reconstructions were performed while varying (from to times the largest eigenvalue).

Fig. 2 [Images not available. See PDF.]

Shepp-Logan Phantom simulations showing the effects of different regularization approaches. (a) The MTX 128 128 spiral trajectory used for simulating the forward encoding and inverse reconstruction. (b) The SVD spectrum of the encoding matrix, with vertical lines showing the truncation of the SVD at different values. The square root of the singular values of Encode Encode after Tikhonov Regularization at corresponding values are plotted in black dotted lines. (c, e) Shepp-Logan images with different noise levels reconstructed with decreasing regularization from left to right with tSVD and Tikhonov regularization respectively. (d, f) Graphs of MSEs against regularization for varying noise levels, with solid lines showing the MSEs between the original Shepp-Logan image and the reconstructed Shepp-Logan image, and the dotted lines showing the MSEs between purely reconstructing the noise and the ideal Shepp-Logan image. Extent of regularization is expressed for tSVD as the cumulative energy truncated ( ), and for Tikhonov regularization as the regularization parameter divided by the largest eigenvalue ( ).

The log of the SVD spectrum (Fig. 2b) is shown, normalized to the maximum singular value. Truncation of the SVD at , and are highlighted with the vertical gray grid lines. The square root of the singular values of Encode Encode after Tikhonov regularization at corresponding values are plotted in black dotted lines. Figure 2c,e show reconstructions using tSVD and Tikhonov regularization, respectively, with increasing noise levels (top to bottom) and decreasing regularization (left to right). Figure 2d,f display mean squared errors (MSE) for different noise levels. The y-axis markers highlight convergence points: (error between the Shepp-Logan image and a zero image) and (error between the Shepp-Logan image and an ‘ideal’ reference Shepp-Logan image accounting for the circular k-space field of view).

Insufficient regularization amplifies noise. In low-noise scenarios, minimal truncation suffices for a low MSE. As noise increases, stronger regularization is needed to prevent divergence from the reference image. However, excessive truncation degrades resolution, causing higher MSEs. In SVD truncation, only low-resolution components remain, resulting in a blurry image (Fig. 2c, left). Over-regularizing with Tikhonov regularization (large ) results in conjugate phase reconstruction, where the reconstructed image converges to a zero image (Fig. 2e, left).

In intermediate regimes, both regularization methods converged to similar reconstruction errors, and all matrix decomposition methods yielded virtually identical results when using the same Tikhonov regularization . The normalized MSEs between them are on the order of , indicating that the choice of regularization outweighs the choice of decomposition in terms of image reconstruction variability.

Generalized encoding

We demonstrated the versatility of Pinv-Recon using a numerical Shepp–Logan phantom by incorporating multiple encoding components into . Ideally, Encode can account for all relevant desired encoding mechanisms as well as undesired encoding distortions.

Desired encoding mechanisms include spatial localization via Fourier encoding (Eq. (2), used in the majority of MR imaging), which can be complemented by advanced techniques such as parallel imaging to exploit spatially varying coil sensitivity profiles, and spectral encoding methods like Dixon’s chemical shift encoding using multiple echo times^14,29. These can be mathematically represented as individual factors in the encoding matrix. For example,

• Coil sensitivity encoding: where is the spatial sensitivity profile of the th receive coil, and

• Multi-echo chemical shift encoding: where is the chemical shift frequency of the th species (with the th echo time).

Undesired encoding distortions inevitably occur in addition to these intended encoding mechanisms, and should be accounted for. Common distortions include:

• off-resonance: where is the local field inhomogeneity at voxel , and is the time of the th readout³⁰,

• Spatiotemporal gradient imperfections: where the true spatiotemporal k-space location replaces the nominal k-space coordinates ^{31, 32, 33–34}, and

• Patient motion: where is a rotation matrix and is a translation vector describing the motion (e.g. due to patient movement) at time . This simplifies to a constant rigid-body shift if the imaging volume is not centered at the gradient isocenter but does not move during imaging.

Each of these factors can be incorporated as a multiplicative term modifying the encoding matrix or as a replacement for nominal k-space or image space coordinates. These encoding components can multiply together to form the net encoding relationship without breaking the linear mapping between image and data, allowing Pinv-Recon to solve for all these effects in one system.

Figure 3 demonstrates incrementally including all of the above mechanisms and distortions, with the encoding matrix in the final column being: where the k-space vector and the spatial coordinates were modified as and respectively to include gradient imperfections and patient motion. The data was also undersampled (acceleration factor 2 2) with the inlcusion of the coil sensitivity encoding in the column 4. Row 3 shows reconstruction using Encode accounting for Gradient encoding only, and row 4 shows reconstruction with all relevant encoding accounted for.

Fig. 3 [Images not available. See PDF.]

Illustration of additional encoding included in Pinv-Recon as a single unified framework. From top to bottom: encoding mechanism, SVD spectrum of encoding matrix, reconstructed image accounting for gradient encoding only, reconstructed image accounting for gradient and other encoding. From left to right: spiral gradient encoding + spatial off-center positioning, + gradient nonlinearities, + off-resonance, + coil sensitivity encoding (with 2x2 undersampling of the data), + chemical shift separation (of water and fat).

In addition to the encoding examples highlighted here, other forms of desired and undesired encoding may also be included^{35, 36, 37–38}, for example simultaneous multi-slice imaging and spatially varying relaxation decay.

One memory handling strategy for constructing an encoding matrix with coil sensitivity is to perform it in a procedural manner for each coil. Consider a dataset acquired with coils, each coil with k-space points and image space points. The resulting full matrix is of dimensions ( ) , and the Gram matrix has dimensions . Rather than constructing and multiplying in its entirety, can be computed for each coil and added together (equation (11)), significantly reducing computational time.

11

Then, instead of preforming a large matrix-matrix multiplication between the inverted Gram matrix and (equation (6)), can be computed for each coil individually then summed together, before multiplying with the inverted Gram matrix to obtain the reconstructed image.

Exemplary applications

Several emerging MRI applications feature small to medium matrix sizes and could benefit from non-FFT-based reconstruction, making them good candidates for Pinv-Recon. Examples include hyperpolarized MRI³⁹, diffusion MRI⁴⁰, certain functional MRI acquisitions⁴¹, multi-nuclear spectroscopic imaging⁴², and low-field MRI⁴³. These often have low SNR and use non-Cartesian trajectories or require distortion corrections, which are challenging for traditional FFT-based pipelines but are naturally handled in the explicit encoding matrix formulation. In the following subsections, we demonstrate Pinv-Recon in various experimental scenarios, beginning with smaller image sizes, and then pushing to larger matrix sizes to show the possibilities of Pinv-Recon. Additional examples are shown in Supplementary Materials 3 and 4.

Low resolution 2D spectral-spatial spiral (hyperpolarized carbon-13 MRI)

As shown in Fig. 4, we first applied Pinv-Recon to an 8-channel hyperpolarized C metabolic imaging dataset acquired with a 2D spiral trajectory (MTX 20 20)⁴⁴.

Fig. 4 [Images not available. See PDF.]

Application of Pinv-Recon on multichannel hyperpolarized carbon-13 data. From left to right: Pyruvate, Lactate, Bicarbonate. From top to bottom: root-sum-of-squares coil combination, Pinv-Recon with coil sensitivity map included, Pinv-Recon with coil sensitivity map included on an overdiscretized spatial grid (MTX 64 64), mean SNR for the first two methods (in blue and red respectively), SRFs and SVD spectra from the reconstruction with and without coil sensitivity included respectively.

Data was reconstructed using (1) Pinv-Recon with gradient encoding only and root-sum-of-squares coil combination, (2) Pinv-Recon with coil sensitivity encoding included, as well as (3) Pinv-Recon with coil sensitivity encoding and overdiscretization to a higher-resolution spatial grid of MTX 64 64. Since the coil sensitivity is spatially-varying but time-invariant, the reconstruction matrix has to be calculated for each of the eight individual slices, but can be used to reconstruct the images across all time points and metabolites in eight matrix-matrix multiplications. Including the coil sensitivity profiles allowed for a more SNR-optimal combination than simple root-sum-of-squares coil combination. Overdiscretized reconstruction can be straightforwardly implemented by defining a finer range of r (Eq. (2)), which is much simpler compared to over-discretized reconstruction using gridding, often applied by interpolating the data to a finer Cartesian grid^4,45. Pinv-Recon also allowed the calculation of the SRF and noise maps demonstrating that inclusion of the coil sensitivities improves the SRF. Moreover, Pinv-Recon provides inherent and straightforward denoising by regularization, either by truncation of the SVD or by Tikhonov regularization, which is especially useful for SNR-limited techniques such as hyperpolarized carbon-13 imaging⁴⁶.

Medium resolution 3D stack-of-spirals (hyperpolarized Xenon-129 MRI)

As shown in Fig. 5, a hyperpolarized Xenon-129 dataset (acquired using a golden-angle stack-of-spirals with MTX 80 80 41) demonstrates Pinv-Recon at a larger matrix size and with a separable Z-dimension. This is an example of non-Cartesian in-plane sampling combined with Cartesian phase-encoding along Z with a single slab-selection pulse. In reconstruction, the uniformly sampled Z-dimension was processed separately: Pinv-Recon was applied only along Z without regularization (equivalent to a standard 1D FFT), while a separate Pinv-Recon was applied to the non-uniformly sampled in-plane 2D spiral. Reconstruction at the native resolution and over-discretized to twice the resolution can be performed at reasonable compute times of just a couple of seconds using Cholesky decomposition.

Fig. 5 [Images not available. See PDF.]

(a) Reconstruction of stack-of-spirals data using, from left to right, Pinv-Recon, , NUFFT, NUFFT with regularization, and gridding. Gridding and NUFFT results in streaking artifacts that can be reduced with regularization, whilst Pinv-Recon eliminates these artifacts and boosts SNR. First row shows the images in the original scale, and second row shows them in log scale. (b) The golden angle spiral in one Z phase-encode of the stack-of-spirals trajectory (c) Top row: the SRFs for tSVD regularization and Tikhonov regularization respectively. Bottom row: the corresponding noise matrices (d) The singular values of Encode, with the singular values of the Gram matrix after Tikhonov regularization shown by the dotted line.

Pinv-Recon was also compared against several methods: (1) , (2) non-uniform FFT (NUFFT) reconstruction⁴⁷ (3) NUFFT reconstruction with additional regularization, and (4) conventional gridding. The benefit of inherent regularization in Pinv-Recon is highlighted, allowing better reconstruction of ill-conditioned undersampled data⁴⁸. The golden-angle spiral trajectory had unevenly-spaced interleaves, creating voids in k-space that translate to density compensation factors larger than unity. This creates an ill-posed problem which gridding algorithms may not be able to solve properly without regularization. Streaking artifacts that were present in the gridding and the NUFFT reconstructed images can be minimized with additional regularization. However, this is inherently accounted for in the Pinv-Recon and implementations with Tikhonov regularization included.

The SRF (Fig. 5c, top row) and the noise matrix (Fig. 5c, bottom row) for the reconstruction of the hyperpolarized Xenon image can be calculated for tSVD and Tikhonov regularization respectively, showing that the image is less well-represented and less noisy towards the center of the image with this level of Tikhonov regularization.

High resolution 2D Cartesian VDPD ( -weighted Knee)

We further demonstrate the efficiency of Pinv-Recon for separable imaging dimensions, using a high resolution 18-channel 2D Variable Density Poisson Disk (VDPD) knee dataset (MTX 384 336), which is evenly sampled in the readout dimension, but variably sampled in the phase encoding dimension (Fig. 6). The dataset undersamples k-space in the phase encoding dimension with an acceleration factor (R) of around 3. Along the frequency encoding line, each phase encode is reconstructed separately, accounting for coil sensitivity.

Fig. 6 [Images not available. See PDF.]

High Resolution -weighted knee image reconstructed with Pinv-Recon. (a) Plot of the k-space trajectory, showing variable sampling along the phase encoding dimension from left to right (b) Reconstructed image using gradient encoding only (c) Reconstructed image using gradient and coil sensitivity encoding (d) Eigenvalue spectrum, for gradient encoding only (blue) and for gradient encoding with coil sensitivity encoding (black) (e) SRF for reconstruction with gradient encoding only (f) SRF for reconstruction with gradient encoding and coil sensitivity encoding.

We show that by incorporating the coil sensitivity profiles, a higher resolution and more accurate image can be reconstructed by ameliorating pseudo-random sub-sampling artifacts (Fig. 6c). From the corresponding SRF map with values close to 1 (Fig. 6f), it can be seen that almost full nominal image resolution can be recovered, in contrast to the gradient-encoding-only SRF map with values of around 0.3 (Fig. 6e).

High resolution 2D radial with undersampling (CG-SENSE challenge dataset)

We explored the application of Pinv-Recon to even larger matrix sizes, demonstrating its computational feasibility using the widely recognized open-source radial brain dataset from the ISMRM CG-SENSE reproducibility challenge¹³. Using this dataset, we show that a MTX 320 320 image can be reconstructed in just a few minutes, with the possibility of including non-Cartesian SENSE reconstruction for undersampled data. The size of the encoding matrix can be reduced by masking out image voxels beyond the field-of-view of interest. A circular mask with diameter equal to the side of the image was used initially, providing an image to which a signal threshold can be applied to generate an even smaller mask for subsequent reconstructions. This approach reduces the number of unknowns without discarding image content.

Instead of using forward and inverse gridding or NUFFT to iteratively reconstruct the image, we demonstrate the computational feasibility of Pinv-Recon as an alternative to non-Cartesian SENSE by simply incorporating the coil sensitivity into the encoding matrix and pseudoinverting it. Reconstruction with acceleration factors of R = 1 to 4 are shown in Fig. 7. gives a similar image to Pinv-Recon but does not straightforwardly provide the SRFs and noise maps (columns 4 and 5 in Fig. 7). Especially in the context of non-Cartesian undersampled reconstruction, these metrics provide valuable information about the optimality of encoding and reconstruction. Using the SVD, whilst computationally more expensive than other matrix decomposition methods, can still be useful for shedding light onto the conditioning of the problem, as shown in the rightmost panel of Fig. 7.

Fig. 7 [Images not available. See PDF.]

CG-SENSE reproducibility challenge dataset. From top to bottom: undersampling factors from 1 to 4. From left to right: reconstruction using CG-SENSE reference implementation in MATLAB, reconstruction using Pinv-Recon with coil sensitivity encoding, reconstruction using , SRF maps, noise map, SVD spectra.

High resolution 3D stack-of-stars (abdominal imaging)

Separable 3D Pinv-Recon was demonstrated using a high-resolution abdominal 3D stack-of-stars spoiled gradient echo acquisition (Fig. 8) which has MTX 320 320 144. Similar to the ISMRM CG-SENSE reproducibility challenge, only image pixels in the masked region are included in reconstruction. The Z phase-encoding dimension was factored out via a 1D FFT, and the XY dimensions were reconstructed using Pinv-Recon, forming a 3D image which is shown in the leftmost column in axial, coronal and sagittal views.

Fig. 8 [Images not available. See PDF.]

Left: Axial, coronal and saggital views of the 3D stack-of-stars abdominal imaging dataset. In the middle four columns, top: varying the number of spokes used in reconstruction, middle: varying the condition number by changing the tSVD threshold, bottom: varying the number of iterations in . Right: SVD spectrum, SRF map, and noise map.

The potential for real-time image reconstruction, i.e., updating the reconstructed image during scanning, is demonstrated in this dataset. Unlike gridding and iterative reconstruction methods, Pinv-Recon is data agnostic in that the reconstruction matrix can be calculated simply using the known trajectory, even before data is acquired. Once data acquisition finishes, only a simple matrix-vector reconstruction remains. The top row of Fig. 8 shows images obtained using increasing fractions of the available data spokes from left to right, indicating the potential for real-time updates during imaging. The acquisition time of certain MRI scans can last up to several minutes (6 minutes for this dataset), and real-time updates can be useful for radiographers or scanner operators to monitor the acquisition.

Comparison of Pinv-Recon and is shown in the bottom two rows of Fig. 8. The middle row shows images obtained using varying thresholds for Pinv-Recon with tSVD, demonstrating the trade-off between resolution and SNR offered by varying the condition number , also illustrated in the SVD spectrum on the right. The bottom row shows images reconstructed using . As is common with CG-based iterative methods, image quality first improves, peaks and then degrades with an increasing number of iterations.

Discussion

With the drastic increase in computational power over the past two decades⁴⁹, inverting large encoding matrices has become much more feasible. For example, Van de Walle et al. (2000) reported approximately 75 minutes to reconstruct a MTX 64 64 image ( = 64 4033)¹¹; today, the same size matrix can be decomposed in seconds (approximately 2 s with SVD or less than 0.1 s with Cholesky on our workstation). In our experiments, a 2 GB encoding matrix (roughly equivalent to a MTX 128 128 image with full Fourier encoding) was inverted via Cholesky Decomposition in about 0.4 s. Memory limitations, a major barrier in the past, are far less restrictive on modern 64-bit systems that can address terabytes of RAM^12,50. Looking forward, ongoing developments in high-performance computing architectures, for example larger GPU memory⁵¹, extreme parallelism with AI chips⁵² and half-precision or mixed-precision arithmetic⁵³, promise to further extend the capabilities of Pinv-Recon. Mixed-precision computations, in particular, offers an interesting approach to speed-up computations by using lower precision for initial calculations but storing results at higher precision, increasing computational efficiency whilst limiting the compromise on accuracy⁵⁴. Moreover, modern randomized algorithms for low-rank matrix approximation, such as the svdsketch function in MATLAB⁵⁵, are emerging as potential ways to approximate the pseudoinverse more quickly by sketching the dominant subspace. While such approaches were not used in this work, they represent promising avenues for handling extremely large encoding matrices by trading off some accuracy for speed.

Table 1 summarizes the various algebraic methods and computational strategies mentioned, their applicable scenarios, as well as their advantages and limitations. Among the matrix inversion methods examined, SVD is the most robust and provides valuable insight (via the SVD spectrum) into the conditioning of the problem^11,56. However, SVD is computationally expensive. Our exploration of alternative factorizations revealed that many of these methods, though not previously considered in MRI reconstruction, can offer substantial speed-ups. In particular, Cholesky decomposition was observed to be over 100 faster than SVD in some cases. These methods operate on the normal equations (decomposing the Gram matrix), effectively reducing the problem size to the smaller dimension of Encode. Depending on the problem size relative to the available compute resources, different computational strategies may be considered. For very small encoding matrices, CPU computations are preferred. GPU acceleration is useful for larger encoding matrices, and block-wise inversion on the GPU is preferred for even larger encoding matrices that do not fit on the GPU memory. Other strategies to better manage the inversion include masking voxels outside region-of-interest, computing in single precision, reducing the problem size by separating uniformly sampled dimensions, and procedural construction of the Gram matrix with coil sensitivity encoding.Table 1

Table summarizing algebraic methods and computational strategies, their applicable scenarios, advantages and limitations respectively. Notes: (i) “Small/medium/large” refer to the size of the encoding problem relative to available CPU/GPU memory and bandwidth. (ii) Methods that explicitly compute the reconstruction matrix ( ) enable direct calculation of SRF and noise matrices, the reusing of Recon, and real-time feedback during acquisition. (iii) Normal-equation / Gram-matrix methods square the singular values (eigenvalues of ), worsening conditioning; Tikhonov regularization ( ) is recommended to ensure positive-definiteness and improve stability.

Method/strategy	Applicable scenarios (by problem size)	Advantages	Limitations
Algebraic methods
SVD of Encode (tSVD)	Small to medium encoding matrices when interpretability is desired and sufficient memory available.	Provides singular-value spectrum (conditioning insight); robust to rank deficiency; Recon explicitly calculated; intuitive regularization via truncation.	Highest compute and memory cost among direct methods; less suitable for very large problems; GPU gains modest for small matrices.
Cholesky of	Medium to large problems where speed is critical and Gram matrix fits in memory.	Very fast among direct solvers; operates on smaller N*N Gram matrix; Recon explicitly calculated; well-optimized BLAS/LAPACK/GPU.	Works on Gram matrix (squared conditioning); requires regularization for positive semi-definite matrix; no singular-value spectrum.
Eigenvalue decomposition of	Small to medium problems when eigenspectrum of Gram is useful for analysis/tuning.	Eigenspectrum ( ) aids interpretability; operates on smaller Gram matrix; Recon explicitly calculated; supports truncation and/or Tikhonov in eigenbasis.	Slower and memory-heavier than Cholesky; still suffers squared conditioning.
QR of Encode	Small to medium matrices when faster factorization than SVD is preferred.	Numerically stable; Recon explicitly calculated; Low memory footprint; scalable; Flexible regularization; avoids explicit inversion.	No singular-value information; can be slower than normal-equation methods.
Iterative solvers on normal equations (e.g., CG, LSQR)	Very large problems where explicit Recon is not needed or memory is constrained.	Numerically stable; scalable; Flexible regularization; avoids explicit inversion; early stopping possible.	No explicit Recon (SRF/noise not immediate); parameter tuning required; convergence may slow for ill-conditioned systems; runtime depends on iterations.
Computational strategies
Block-wise inversion (Encode or )	Problems that do not fit into device memory; multi-GPU/cluster or out-of-core settings.	Enables out-of-core and parallel decomposition; leverages hardware to scale up.	Overhead and implementation complexity; limited benefit for small problems.
Masking voxels outside region-of-interest	All sizes, especially large images.	Significantly reduces the size of Encode.	Limited benefit for smaller problems.
Precision choice (single/mixed precision)	All sizes, especially memory-limited cases.	Cuts memory and bandwidth; often shortens runtime on modern accelerators.	Potentially reduced numerical robustness for ill-conditioned problems.
Multi-threaded CPU	All sizes where CPU cores are available.	Exploits available cores; stable, predictable; more available memory.	May lag optimized GPU for very large problems; performance depends on library setup.
GPU acceleration	Medium to large problems that fit in GPU memory; benefit grows with size.	Significant speedups for large dense factorizations; batched operations possible.	For very small problems, kernel overhead dominates; very large problems can exceed memory; CPU-GPU transfers can bottleneck.
Problem separability via FFT (decouple fully sampled dimensions)	Whenever >= 1 dimension is uniformly and fully sampled.	Reduces effective matrix dimension before inversion; shrinks memory and runtime.	Often when separability holds; requires careful scaling/density-compensation handling.
Procedural/coil-wise Gram matrix construction	Large matrices with coil sensitivity encoding.	Reduce memory consumption by avoiding large matrix-matrix multiplications.	Added implementation complexity; speed-up negligible for small problems or when structure does not factor cleanly.

Perhaps the greatest strength of Pinv-Recon is its simplicity and versatility. It seamlessly accommodates non-Cartesian trajectories and a variety of corrections that would otherwise require complex dedicated algorithms^{5,57, 58–59}. Conventional reconstruction pipelines often involve multiple stages (e.g., density compensation^60,61, trajectory correction, off-resonance modulation, motion gating, coil combination, deapodization, gradient nonlinearity correction), each adding complexity and variability. Many of these algorithms cannot be easily combined with one another. In Pinv-Recon, all these effects are incorporated into the single encoding matrix, thereby simplifying the reconstruction workflow and improving reproducibility. Iterative methods can, in principle, also incorporate such factors, as exemplified by previous work using the explicit encoding matrix in conjugate-gradient methods, for example Kadah and Hu’s Algebraic Reconstruction for Magnetic Resonance Imaging Under Inhomogeneity⁹ and Wilm et al.’s Higher order reconstruction for MRI in the presence of spatiotemporal field perturbations²⁴. Leveraging modern compute resources, Pinv-Recon implements this more simply, and also calculates Recon to provide, by simple matrix-matrix multiplication, additional metrics such as SRF maps and noise amplification maps that can serve as quality assurance measures and guide future sequence design^{27,62, 63, 64–65}. Unlike iterative reconstructions, which require full data before processing, Pinv-Recon can pre-compute Recon and reuse it across acquisitions with the same encoding, offering the option of real-time updates during acquisition and leaving only a matrix-vector multiplication after acquisition. Regularization—implemented transparently via tSVD, Tikhonov regularization, or both—further reduces noise and stabilizes ill-conditioned, undersampled problems, as demonstrated in our Xenon-129 dataset. Pinv-Recon can be further extended to non-linear reconstructions⁶⁶, for example compressed-sensing problems, where k-space is randomly subsampled⁶⁷, as well as reconstruction using implicit neural representation⁶⁸.

Recent work related to pseudoinversion in MRI reconstruction has also investigated various matrix decomposition techniques to facilitate efficient inversion, including QR decomposition and Jacobi SVD^69,70. These methods typically focus on unaliasing undersampled Cartesian images by inverting much smaller coil sensitivity encoding matrices after an initial reconstruction using gradient encoding. In contrast, the present work directly inverts the full encoding matrix (including gradient and coil sensitivity terms) and extends this framework to non-Cartesian sampling schemes.

With the move towards open-source and reproducible MR acquisition and reconstruction, as evidenced by the growing popularity of the vendor-neutral pulse sequence programming environment Pulseq, Pinv-Recon shows great promise as a simple, versatile, and reliable reconstruction method^71,72. An example implementation of Pinv-Recon for Pulseq is now available and can be found online⁷³. Researchers developing novel reconstruction techniques may use Pinv-Recon as a reference reconstruction method since it provides the generic linear-least squares solution (or the minimum 2-norm solution for underdetermined cases). Whilst the simulations and characterization in this manuscript were run in MATLAB for its efficient matrix operations, we have also included Python-equivalent code in the GitHub repository linked in the Data Availability Statement, to provide an option to run the code in an open-source programming language. Moreover, Pinv-Recon can also be used as an instructive educational tool for students familiarizing themselves with various aspects of MR physics, exploiting its modularity to visualize each encoding mechanism.

We acknowledge that Pinv-Recon is not intended to replace all FFT-based or iterative methods; rather, it excels in small to medium matrix size applications (up to approximately 256 256 in-plane) where its computational burden is manageable and its advantages in flexibility and reproducibility are most pronounced. For very high-resolution imaging, iterative methods may remain more practical.

In conclusion, we have shown that direct pseudoinversion of the encoding matrix—once considered impractical—is now within reach for many MR reconstruction problems. Pinv-Recon leverages modern computational resources and optimized linear algebra routines to provide a simple, linear, and extendable reconstruction framework that is particularly well-suited for MRI applications with small to medium matrix sizes and specialized encoding requirements. The presented framework can unify multiple linear encoding strategies and correction mechanisms into a single, comprehensive algebraic model. This one-step approach not only computes the image directly but also yields valuable image metrics such as spatial response functions and noise maps, and has the potential for real-time feedback during imaging, making it a valuable tool for both research and clinical applications.

Methods

Computational evaluation and numerical simulations

The computational evaluation and numerical simulations were performed as described in Section 3, using Matlab 2023b (The MathWorks, MA, USA). The MATLAB implementation can be also found in the Git repository linked in the Data Availability Statement.

Applications

All reconstructions were performed using the equations outlined in Section 2. Details of the acquisition of each dataset are as follows:

Hyperpolarized carbon-13 dataset: A healthy volunteer was scanned after IV injection of 35 mL ofhyperpolarized [1- C]pyruvate (polarized for approximately 4 hours in the SPINlab (GE Healthcare, WI, USA)). The acquisition was a spectral-spatial single-arm spiral sequence covering a 240 240 mm field-of-view with eight 20 mm slices; flip angles were for pyruvate, for lactate, and for bicarbonate, acquired at 4 s intervals. The scan was performed on a GE 3T Premier (GE Healthcare, WI, USA) system using an 8-channel transmit-receive C/ H head coil (Rapid Biomedical, Rimpar, Germany). Coil sensitivity maps were estimated via polynomial smoothing of pyruvate images from each individual coil and normalized by the sum-of-squares of all coils. This project was approved by the University of Nottingham Medical School Ethics Committee (Ref: 416-1911 BRAIN DNP) and was conducted in accordance with the standards set by the latest version of the Declaration of Helsinki, except for registration in a database. Informed, written consent was obtained from the volunteer before participation.

Hyperoplarized xenon-129 dataset: A healthy 34-year-old female participant inhaled 1 L of enriched Xe gas (polarized for approximately 30 minutes). Imaging was performed on a 3T GE Premier scanner with a flexible Xe chest coil (PulseTeq, Cobham, UK) using a stack-of-spirals 3D sequence covering a 400 400 200 mm field-of-view. The sequence employed 16 spiral interleaves with golden-angle rotation, a bandwidth of 250 kHz, TR = 15.6 ms, TE = 2.3 ms, and a hard excitation pulse yielding an effective flip angle of approximately per volume. Participants gave informed written consent, and ethics approval was provided by the University of Oxford Central University Research Ethics Approval (Reference: R77150/RE001).

Pinv-Recon of the dataset was compared with conjugate gradient iterative reconstruction with 11 iterations, NUFFT using the Michigan Image Reconstruction Toolbox (MIRT)⁴⁷(density compensation calculated using the Pipe and Menon method⁷⁴, NUFFT object calculated with minmax interpolator with Kaiser-Bessel scaling, kernel width of 6, table-based interpolation with a oversampling factor) without and with additional regularization with quadratic penalized weighted least squares (10 iterations, with a beta of ), and conventional gridding (Kaiser-Bessel filtering⁴ with kernel width of 3 and oversampling to 1.125 of the original grid size, density compensation factor determined using Voronoi triangulation).

T₂-weighted knee: The knee dataset was acquired with a dedicated knee coil with 18 receive channels using a 2D fast recovery fast spin echo (frFSE) acquisition with fat saturation. The acquisition parameters were TR = 6 s, TE = 108 ms. The trajectory was a 1D variable density poison disc (VDPD) sampling with image matrix size of 384 in the readout direction and 336 in the phase encoding direction. The phase encoding dimension was acquired with 3.23 acceleration (104/336), and the center of k-space (-12:12) was fully sampled for coil sensitivity calibration.

ISMRM CG-SENSE challenge dataset: The ISMRM CG-SENSE Challenge dataset was obtained from CG-SENSE revisited: Results from the first ISMRM reproducibility challenge¹³.

Stack-of-stars abdominal dataset: For the 3D stack-of-stars dataset (using the Liver Acquisition with Volume Acceleration LAVA-Star sequence), the acquisition parameters were FOV = 40.0 40.0 31.68 cm, resolution = 1.25 1.25 2.2 mm, TE/TR = 1.488/3.276 ms, FA = 12 , BW = ±62.5 kHz, with interleaved fat suppression, resulting in a scan time of 5 min 58 sec.

Author contributions

K.Y.: writing - original draft (lead); conceptualization (supporting); data collection (equal); methodology (equal); software (equal); formal analysis (equal); writing - review and editing (equal). C.T.: conceptualization (supporting); writing - review and editing (equal). R.F.S.: conceptualization (supporting); writing - review and editing (equal). B.W.: visualization; writing - review and editing (equal). A.M.: data collection (equal). S.S.: data collection (equal); writing - review and editing (equal). P.M.: data collection (equal). D.A.: data collection (equal). F.V.G.: funding acquisition (equal); supervision. D.J.T.: funding acquisition (equal); supervision (equal); writing - review and editing (equal). J.T.G.: funding acquisition (equal); supervision (equal); writing - review and editing (equal). F.W.: conceptualization (lead); writing - original draft (supporting); methodology (equal); software (equal); formal analysis (equal); supervision (equal); writing - review and editing (equal). All authors read and approved the final manuscript.

Funding

K.Y. acknowledges an Oxford-Medical Research Council Doctoral Training Partnership iCASE award, the Oxford-Radcliffe Scholarship, and GE HealthCare for graduate funding. D.J.T. was funded by a British Heart Foundation Senior Basic Science Research Fellowship (FS/19/18/34252). J.T.G. is funded by the Oxford Biomedical Research Centre, and acknowledges the BHF Centre of Research Excellence, University of Oxford, and the UK Medical Research Council for funding. The hyperpolarized Xenon-129 datasets were funded by the British Heart Foundation (RE/18/3/34214) and the hyperpolarized Carbon-13 dataset was funded by Nottingham Life cycle 5.

Data availability

The code supporting the numerical simulations in this manuscript, as well as some example datasets can be found at https://github.com/univ39/Pinv-Recon.

Declarations

Competing interests

C.T. is an employee of The Mathworks. R.F.S. and F.W. are employees of GE HealthCare. S.S. is an Editorial Board Member of Scientific Reports. The remaining authors declare no relevant conflicts of interest.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1038/s41598-025-21929-z.

Publisher’s note

References

1. Webb, A; Obungoloch, J. Five steps to make MRI scanners more affordable to the world. Nature; 2023; 615, pp. 391-393.2023Natur.615.391W1:CAS:528:DC%2BB3sXltlGrs7w%3D [DOI: https://dx.doi.org/10.1038/d41586-023-00759-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36918678]

2. Lauterbur, PC. Image formation by induced local interactions: Examples employing nuclear magnetic resonance. Nature; 1973; 242, pp. 190-191.1973Natur.242.190L1:CAS:528:DyaE3sXhs1Kqu74%3D [DOI: https://dx.doi.org/10.1038/242190a0]

3. Edelman, R.R. The history of MR imaging as seen through the pages of radiology. Radiology. https://doi.org/10.1148/radiol.14140706.

4. Beatty, PJ; Nishimura, DG; Pauly, JM. Rapid gridding reconstruction with a minimal oversampling ratio. IEEE Trans. Med. Imaging; 2005; 24, pp. 799-808.2005ITMI..24.799B [DOI: https://dx.doi.org/10.1109/TMI.2005.848376] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15959939]

5. Wright, KL; Hamilton, JI; Griswold, MA; Gulani, V; Seiberlich, N. Non-Cartesian parallel imaging reconstruction. J. Magnetic Resonance Imaging; 2014; 40, pp. 1022-1040. [DOI: https://dx.doi.org/10.1002/jmri.24521]

6. Man, L-C; Pauly, JM; Macovski, A. Multifrequency interpolation for fast off-resonance correction. Magnetic Resonance Med.; 1997; 37, pp. 785-792.1:STN:280:DyaK2s3nt1Wjuw%3D%3D [DOI: https://dx.doi.org/10.1002/mrm.1910370523]

7. Pruessmann, KP; Weiger, M; Börnert, P; Boesiger, P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magnetic Resonance Med.; 2001; 46, pp. 638-651.1:STN:280:DC%2BD3MrksVSgsQ%3D%3D [DOI: https://dx.doi.org/10.1002/mrm.1241]

8. Chen, D-Q; Marr, RB; Lauterbur, PC. Reconstruction from NMR data acquired with imaging gradients having arbitrary time dependence. IEEE Trans. Med. Imaging; 1986; 5, pp. 162-164.1986ITMI..5.162C1:STN:280:DC%2BD1c%2Fos1CmsQ%3D%3D [DOI: https://dx.doi.org/10.1109/TMI.1986.4307765] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18244002]

9. Kadah, Y; Hu, X. Algebraic reconstruction for magnetic resonance imaging under B/sub 0/ inhomogeneity. IEEE Trans. Med. Imaging; 1998; 17, pp. 362-370.1998ITMI..17.362K1:STN:280:DyaK1cvgsFejtQ%3D%3D [DOI: https://dx.doi.org/10.1109/42.712126] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9735900]

10. Cruz, G., Yaman, B., Akçakaya, M., Doneva, M. & Prieto, C. Chapter 2 - MRI Reconstruction as an Inverse Problem. In Akçakaya, M., Doneva, M. & Prieto, C. (eds.) Advances in Magnetic Resonance Technology and Applications, 7 of Magnetic Resonance Image Reconstruction, 37–57, https://doi.org/10.1016/B978-0-12-822726-8.00011-7 (Academic Press, 2022).

11. Walle, R et al. Reconstruction of MR images from data acquired on a general nonregular grid by pseudoinverse calculation. IEEE Trans. Med. Imaging; 2000; 19, pp. 1160-1167.2000ITMI..19.1160V [DOI: https://dx.doi.org/10.1109/42.897806] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11212364]

12. Wang, S.P. Computer Memory and Storage. In Wang, S.P. (ed.) Computer Architecture and Organization: Fundamentals and Architecture Security, 45–69, https://doi.org/10.1007/978-981-16-5662-0_3 (Springer, Singapore, 2021).

13. Maier, O et al. CG-SENSE revisited: Results from the first ISMRM reproducibility challenge. Magnetic Resonance Med.; 2021; 85, pp. 1821-1839. [DOI: https://dx.doi.org/10.1002/mrm.28569]

14. Pruessmann, KP. Encoding and reconstruction in parallel MRI. NMR Biomed.; 2006; 19, pp. 288-299. [DOI: https://dx.doi.org/10.1002/nbm.1042] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16705635]

15. Kirchner, T; Fillmer, A; Tsao, J; Pruessmann, KP; Henning, A. Reduction of voxel bleeding in highly accelerated parallel 1H MRSI by direct control of the spatial response function. Magnetic Resonance Med.; 2015; 73, pp. 469-480.1:CAS:528:DC%2BC2MXht1Kltrw%3D [DOI: https://dx.doi.org/10.1002/mrm.25185]

16. Barata, JCA; Hussein, MS. The Moore-Penrose Pseudoinverse: A Tutorial Review of the Theory. Braz. J. Phys.; 2012; 42, pp. 146-165.2012BrJPh.42.146B [DOI: https://dx.doi.org/10.1007/s13538-011-0052-z]

17. Golub, G.H. & Loan, C. F.V. Matrix Computations (JHU Press, 2013).

18. Gander, W. Algorithms for the QR decomposition. Res. Rep.; 1980; 80, pp. 1251-1268.

19. Tikhonov, A.N. & Arsenin, V. I. V.I. Solutions of Ill-Posed Problems. Scripta Series in Mathematics (Winston and distributed solely by Halsted Press, 1977).

20. Hansen, PC; O’Leary, DP. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput.; 1993; 14, pp. 1487-1503.1241596 [DOI: https://dx.doi.org/10.1137/0914086]

21. Pauly, J., Chapter 5 - Reconstruction of Non-Cartesian Data (2005). Available at: https://users.fmrib.ox.ac.uk/~karla/reading_group/lecture_notes/AdvRecon_Pauly_read.pdf.

22. Chapter 5: Tikhonov regularization. In Linear and Nonlinear Inverse Problems with Practical Applications, Computational Science & Engineering, 63–82, https://doi.org/10.1137/1.9781611972344.ch5 (Society for Industrial and Applied Mathematics, 2012).

23. Chen, J., Jin, Z., Shi, Q., Qiu, J. & Liu, W. Block algorithm and its implementation for Cholesky factorization. ICCGI21 (2013).

24. Wilm, BJ; Barmet, C; Pavan, M; Pruessmann, KP. Higher order reconstruction for MRI in the presence of spatiotemporal field perturbations. Magnet. Resonance Med.; 2011; 65, pp. 1690-1701. [DOI: https://dx.doi.org/10.1002/mrm.22767]

25. Li, S et al. Algebraic reconstruction technique for parallel imaging reconstruction of undersampled radial data: Application to cardiac cine. Magnetic Resonance Med.; 2015; 73, pp. 1643-1653. [DOI: https://dx.doi.org/10.1002/mrm.25265]

26. Paige, CC; Saunders, MA. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Softw.; 1982; 8, pp. 43-71.661121 [DOI: https://dx.doi.org/10.1145/355984.355989]

27. Sánchez-González, J et al. Minimum-norm reconstruction for sensitivity-encoded magnetic resonance spectroscopic imaging. Magnetic Resonance Med.; 2006; 55, pp. 287-295. [DOI: https://dx.doi.org/10.1002/mrm.20758]

28. Mathworks, Lapack in matlab (2025). Available at: https://uk.mathworks.com/help/matlab/math/lapack-in-matlab.html.

29. Wiesinger, F et al. IDEAL spiral CSI for dynamic metabolic MR imaging of hyperpolarized [1-13C]pyruvate. Magnetic Resonance Med.; 2012; 68, pp. 8-16.1:CAS:528:DC%2BC38XotlCnsrs%3D [DOI: https://dx.doi.org/10.1002/mrm.23212]

30. Haskell, MW; Nielsen, J-F; Noll, DC. Off-resonance artifact correction for MRI: A review. NMR Biomed.; 2023; 36, [DOI: https://dx.doi.org/10.1002/nbm.4867] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36326709]e4867.

31. Duyn, J.H., Yang, Y., Frank, J.A. & van der Veen, J.W. Simple correction method for k-space trajectory deviations in MRI. Journal of Magnetic Resonance (San Diego, Calif.: 1997)132, 150–153, https://doi.org/10.1006/jmre.1998.1396 (1998).

32. Rahmer, J; Mazurkewitz, P; Börnert, P; Nielsen, T. Rapid acquisition of the 3D MRI gradient impulse response function using a simple phantom measurement. Magnetic Resonance Med.; 2019; 82, pp. 2146-2159. [DOI: https://dx.doi.org/10.1002/mrm.27902]

33. Barmet, C; Zanche, ND; Pruessmann, KP. Spatiotemporal magnetic field monitoring for MR. Magnetic Resonance Med.; 2008; 60, pp. 187-197. [DOI: https://dx.doi.org/10.1002/mrm.21603]

34. Sipilä, P; Greding, S; Wachutka, G; Wiesinger, F. 2H transmit-receive NMR probes for magnetic field monitoring in MRI. Magnetic Resonance Med.; 2011; 65, pp. 1498-1506.1:CAS:528:DC%2BC3MXkvVartr0%3D [DOI: https://dx.doi.org/10.1002/mrm.22741]

35. King, K.F. & Angelos, L. SENSE image quality improvement using matrix regularization. In Proceedings of the 9th Annual Meeting of ISMRM, Glasgow, Scotland, 1771 (2001).

36. Maudsley, AA. Fourier imaging using rf phase encoding. Magnetic Resonance Med.; 1986; 3, pp. 768-777.1:STN:280:DyaL2s%2FmtlGqtg%3D%3D [DOI: https://dx.doi.org/10.1002/mrm.1910030512]

37. Kartäusch, R et al. Spatial phase encoding exploiting the Bloch-Siegert shift effect. Magnetic Resonance Mater. Phys. Biol. Med.; 2014; 27, pp. 363-371. [DOI: https://dx.doi.org/10.1007/s10334-013-0417-0]

38. Sacolick, LI; Wiesinger, F; Hancu, I; Vogel, MW. B1 mapping by Bloch-Siegert shift. Magnetic Resonance Med.; 2010; 63, pp. 1315-1322. [DOI: https://dx.doi.org/10.1002/mrm.22357]

39. Larson, PEZ; Gordon, JW. Hyperpolarized metabolic MRI–acquisition, reconstruction, and analysis methods. Metabolites; 2021; 11, 386.1:CAS:528:DC%2BB3MXhslWgu77P [DOI: https://dx.doi.org/10.3390/metabo11060386] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34198574][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8231874]

40. Kiselev, VG. Fundamentals of diffusion MRI physics. NMR Biomed.; 2017; 30, [DOI: https://dx.doi.org/10.1002/nbm.3602] e3602.

41. Hernandez-Garcia, L; Lahiri, A; Schollenberger, J. Recent progress in ASL. NeuroImage; 2019; 187, pp. 3-16. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2017.12.095] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29305164]

42. Wei, Y et al. Multi-nuclear magnetic resonance spectroscopy: State of the art and future directions. Insights Into Imaging; 2022; 13, 135. [DOI: https://dx.doi.org/10.1186/s13244-022-01262-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35976510][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9382599]

43. Arnold, TC; Freeman, CW; Litt, B; Stein, JM. Low-field MRI: Clinical promise and challenges. J. Magnetic Resonance Imaging; 2023; 57, pp. 25-44. [DOI: https://dx.doi.org/10.1002/jmri.28408]

44. Schulte, RF et al. Saturation-recovery metabolic-exchange rate imaging with hyperpolarized [1-13C] pyruvate using spectral-spatial excitation. Magnetic Resonance Med.; 2013; 69, pp. 1209-1216.2013mri.book...S1:CAS:528:DC%2BC3sXmtlWnsLg%3D [DOI: https://dx.doi.org/10.1002/mrm.24353]

45. Schomberg, H; Timmer, J. The gridding method for image reconstruction by Fourier transformation. IEEE Trans. Med. Imaging; 1995; 14, pp. 596-607.1995ITMI..14.596S1:STN:280:DC%2BD1c%2FmtlOjtQ%3D%3D [DOI: https://dx.doi.org/10.1109/42.414625] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18215864]

46. Kim, Y et al. Denoising of hyperpolarized 13C MR images of the human brain using patch-based higher-order singular value decomposition. Magnetic Resonance Med.; 2021; 86, pp. 2497-2511.1:CAS:528:DC%2BB3MXhslGjt7jP [DOI: https://dx.doi.org/10.1002/mrm.28887]

47. Fessler, J.A., Michigan Image Reconstruction Toolbox (MIRT) (0000). Available at: https://web.eecs.umich.edu/~fessler/code/.

48. Uecker, M et al. Real-time MRI at a resolution of 20 ms. NMR Biomed.; 2010; 23, pp. 986-994. [DOI: https://dx.doi.org/10.1002/nbm.1585] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20799371]

49. Moore, GE. Cramming more components onto integrated circuits. Electronics; 1965; 38, pp. 1-14.

50. Mathworks, Table of Matlab release features (2024). Available at: https://uk.mathworks.com/matlabcentral/discussions/general/847991-table-of-matlab-release-features.

51. NVIDIA, NVIDIA Blackwell Architecture (2025). Available at: https://www.nvidia.com/en-us/data-center/technologies/blackwell-architecture/.

52. Cerebras, Product - Chip. (2024). Available at: https://cerebras.ai/product-chip/

53. Mathworks, What is half precision? (2025). Available at: https://uk.mathworks.com/help/coder/ug/what-is-half-precision.html.

54. Higham, NJ; Mary, T. Mixed precision algorithms in numerical linear algebra. Acta Numerica; 2022; 31, pp. 347-414.4436588 [DOI: https://dx.doi.org/10.1017/S0962492922000022]

55. Mathworks, Svdsketch - compute svd of low-rank matrix sketch (2025). Available at: https://uk.mathworks.com/help/matlab/ref/svdsketch.html.

56. Desplanques, B; Cornelis, D; Achten, E; Walle, R; Lemahieu, I. Iterative reconstruction of magnetic resonance images from arbitrary samples in k-space. IEEE Trans. Nuclear Sci.; 2002; 49, pp. 2268-2273.2002ITNS..49.2268D [DOI: https://dx.doi.org/10.1109/TNS.2002.803673]

57. Zhu, K et al. Hybrid-space SENSE reconstruction for simultaneous multi-slice MRI. IEEE Trans. Med. Imaging; 2016; 35, pp. 1824-1836.2016ITMI..35.1824Z [DOI: https://dx.doi.org/10.1109/TMI.2016.2531635] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26915118][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4988924]

58. Harkins, KD; Does, MD; Grissom, WA. Iterative method for Predistortion of MRI gradient waveforms. IEEE Trans. Med. Imaging; 2014; 33, pp. 1641-1647.2014ITMI..33.1641H [DOI: https://dx.doi.org/10.1109/TMI.2014.2320987] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24801945][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4128553]

59. Uecker, M et al. Espirit—An eigenvalue approach to autocalibrating parallel mri: Where sense meets grappa. Magnetic Resonance Med.; 2014; 71, pp. 990-1001. [DOI: https://dx.doi.org/10.1002/mrm.24751]

60. Rasche, V; Proksa, R; Sinkus, R; Börnert, P; Eggers, H. Resampling of data between arbitrary grids using convolution interpolation. IEEE Trans. Med. Imaging; 1999; 18, pp. 385-392.1999ITMI..18.385R1:STN:280:DyaK1MzktlWruw%3D%3D [DOI: https://dx.doi.org/10.1109/42.774166] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10416800]

61. Tan, H. & Zheng, Y. Point spread function optimization for MRI reconstruction. In Proceedings. (ICASSP ’05). IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005., 2, ii/477–ii/480 Vol. 2, https://doi.org/10.1109/ICASSP.2005.1415445 (2005).

62. Pruessmann, KP; Weiger, M; Scheidegger, MB; Boesiger, P. SENSE: Sensitivity encoding for fast MRI. Magnetic Resonance Med.; 1999; 42, pp. 952-962.1:STN:280:DC%2BD3c%2FhtFCruw%3D%3D [DOI: https://dx.doi.org/10.1002/(SICI)1522-2594(199911)42:5<952::AID-MRM16>3.0.CO;2-S]

63. Roemer, PB; Edelstein, WA; Hayes, CE; Souza, SP; Mueller, OM. The NMR phased array. Magnetic Resonance Med.; 1990; 16, pp. 192-225.1:STN:280:DyaK3M%2FpsVKhsg%3D%3D [DOI: https://dx.doi.org/10.1002/mrm.1910160203]

64. Kellman, P; McVeigh, ER. Image reconstruction in SNR units: A general method for SNR measurement. Magnetic Resonance Med.; 2005; 54, pp. 1439-1447. [DOI: https://dx.doi.org/10.1002/mrm.20713]

65. Tunnicliffe, E.M., Graves, M.J. & Robson, M.D. Use of the Noise Covariance Matrix in Array Coil Quality Assurance. In Proceedings of the International Society for Magnetic Resonance in Medicine (Montreal, Quebec, Canada, 2011).

66. Uecker, M; Hohage, T; Block, KT; Frahm, J. Image reconstruction by regularized nonlinear inversion–Joint estimation of coil sensitivities and image content. Magnetic Resonance Med.; 2008; 60, pp. 674-682. [DOI: https://dx.doi.org/10.1002/mrm.21691]

67. Lustig, M; Donoho, D; Pauly, JM. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magnetic Resonance Med.; 2007; 58, pp. 1182-1195. [DOI: https://dx.doi.org/10.1002/mrm.21391]

68. Zhu, Y; Liu, Y; Zhang, Y; Liang, D. Implicit neural representation for medical image reconstruction. Phys. Med. Biol.; 2025; [DOI: https://dx.doi.org/10.1088/1361-6560/addfa5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/41285110][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC12087020]

69. Ullah, I et al. QR-decomposition based SENSE reconstruction using parallel architecture. Comput. Biol. Med.; 2018; 95, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.compbiomed.2018.01.013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29428871]

70. Qazi, SA; Saeed, A; Nasir, S; Omer, H. Singular value decomposition using Jacobi Algorithm in pMRI and CS. Appl. Magnetic Resonance; 2017; 48, pp. 461-471. [DOI: https://dx.doi.org/10.1007/s00723-017-0874-0]

71. Layton, KJ et al. Pulseq: A rapid and hardware-independent pulse sequence prototyping framework. Magnetic Resonance Med.; 2017; 77, pp. 1544-1552.1:CAS:528:DC%2BC2sXkt1SnsLY%3D [DOI: https://dx.doi.org/10.1002/mrm.26235]

72. Veldmann, M et al. Open-source MR imaging and reconstruction workflow. Magnetic Resonance Med.; 2022; 88, pp. 2395-2407. [DOI: https://dx.doi.org/10.1002/mrm.29384]

73. Nielsen, J.-F. SequenceExamples-GE/pge2/spiral at main HarmonizedMRI/SequenceExamples-GE (2025). Available at: https://github.com/HarmonizedMRI/SequenceExamples-GE/tree/main/pge2/spiral.

74. Pipe, J. G. & Menon, P. Sampling density compensation in MRI: Rationale and an iterative numerical solution. Magnetic Resonance Med.41, 179–186. (1999).