## Data Acquisition

The 5D EP-COSI data was acquired on a Siemens 3T Skyra scanner (Siemens Healthineer, Erlangen, Germany) with a dedicated “receive” 24-channel phase-array breast coil and a body “transmit” coil (FOV: 160×160×120mm3, 1.5 mL voxel volume, TR/TE were 1500/35 ms). 64 t1 points sampled were used along F1 (15) with a spectral bandwidth (SW) of 1250 Hz and 512 complex t2 points, and a SW of 1190 Hz along F2. A three-pulse sequence (27) was employed before the global water suppression. A non-water suppressed scan with one t1 point was acquired for eddy current phase correction and coil combination (28). Two spatial and one spectral dimensions (ky,kz,t1) were non-uniformly sampled with an exponentially-weighted sampling density along t1 and gaussian sampling density along the ky-kz plane for an acceleration factor of 8. The total scan time was 28 minutes and 48 seconds.

## Covariance and Inner-Product COSY processing

After the GS-CS reconstruction, a hybrid spectral-spatial data matrix \(\mathbf{D}\)(x, y, z, t2, t1\()\) was outputted where x, y and z were the Fourier transform of kx, ky and kz dimensions in k-space. After Fourier transforming the direct spectral dimension (t2), we get the mixed time-frequency matrix \(\mathbf{A}\)(F2, t1) for every special location (x = 1,2,3…,16; y = 1,2,3…,16; z = 1,2,3…,8), with a stack of 1D spectra. Fourier transforming the indirect dimension t1 of \(\mathbf{A}\) yielded the conventional Fourier transformed COSY spectrum \(\mathbf{S}\)(F2, F1\()\) of size (512×64).

A covariance transform was instead obtained from **A** in the following manner (25).

Step 1: Make matrix \(\mathbf{A}\) offset free by subtracting the average 1D spectrum (\({\mathbf{A}}_{\text{a}\text{v}\text{g}}\)) from it.

\(\stackrel{\sim}{\mathbf{A}} = \mathbf{A} – {\mathbf{A}}_{\text{a}\text{v}\text{g}}\) . This 1D spectrum is formed by averaging over the t1 dimension in \(\mathbf{A}\).

Step 2: Apply a singular value decomposition (SVD) (30) to the transposed mixed time-frequency matrix \({\stackrel{\sim}{\mathbf{A}}}^{\text{T}}\) = \(\stackrel{\sim}{\mathbf{U}}\).\(\stackrel{\sim}{\mathbf{W}}\).\({\stackrel{\sim}{\mathbf{V}}}^{\text{T}}\), where \(\stackrel{\sim}{\mathbf{U}}\) and \(\stackrel{\sim}{\mathbf{V}}\) are the singular vectors and \(\stackrel{\sim}{\mathbf{W}}\) is the diagonal matrix with singular values as its diagonal elements.

Step 3: The final covariance transformed spectrum is then calculated as,

, which gives a high-resolution spectrum of size (512×512), with the spectral width along indirect spectral dimension same as that of the direct spectral dimension.

The similarity between covariance and Fourier-transformed spectra is based on the Parseval’s theorem as described in (20). It follows that the covariance spectrum \(\mathbf{C}= {\stackrel{\sim}{\mathbf{A}}}^{\text{T}}.\mathbf{A}\) correspond to the 2D spectrum squared, i.e., \({\mathbf{S}}^{2}\), provided that the \({\mathbf{A}}_{\text{a}\text{v}\text{g}}\) in the calculation of covariance vanishes. However, it is shown later that the condition of vanishing average might not hold near the spectral center due to the relatively slow spin precision, as well as when reducing the number of t1 points (26). It is shown that this problem can be overcome by discarding the average. i.e, by choosing \(\stackrel{\sim}{\mathbf{A}} = \mathbf{A}\) in step 1. This is called the inner-product covariance transform where the Parseval’s theorem assures the correspondence between a matrix \(\mathbf{I}= {\mathbf{A}}^{\text{T}}.\mathbf{A}\) and \({\mathbf{S}}^{2}\), even if the average doesn’t vanish (26). High resolution IP covariance matrix is therefore calculated as

\({\mathbf{S}}_{\mathbf{I}\mathbf{P}} = {\left({\mathbf{A}}^{\text{T}}.\mathbf{A}\right)}^{1/2} = \mathbf{U}.\mathbf{W}.{\mathbf{U}}^{\text{T}}\) , where \(\mathbf{U}\), \(\mathbf{V}\) and \(\mathbf{W}\) are the singular vectors singular values of \({\mathbf{A}}^{\text{T}}\).

Even though high resolution can be achieved by CT while also minimizing t1-ridging, reduced number of t1 points can introduce spurious correlations in the CT/IP spectrum. However, the deterministic nature of the sampling scheme and resonant positions enable us to mask out these spurious correlations while retaining the advantage in resolution (31). This may be approached in multiple ways. One option is to determine these spurious correlations based on their intensity, frequencies, and t1 sampling as shown in (31). Another option is to use a map of true resonance frequencies obtained via high resolution FFT as a reference to identify spurious correlations. In this work, we used the resonance frequency locations identified from a prior knowledge based synthetic spectra to remove the spurious correlations.