ACCELERATED MAGNETIC RESONANCESPECTROSCOPY WITH VANDERMONDE FACTORIZATION XiaoboQu⋆1 Jiaxi Ying1 Jian-FengCai2 Zhong Chen1 7 1 1 DepartmentofElectronicScience, XiamenUniversity,Xiamen,China 0 2 DepartmentofMathematics,Hong KongUniversityofScience andTechnology,HongKongSAR, China 2 n a J ABSTRACT peaks inevitably results in loss of sparsity since more non- 4 zerosarepresentedinthespectrum[5]. ThislimitationofCS Multi-dimensional magnetic resonance spectroscopy is an 2 leadtothedegenerationofreconstructionperformance. importanttoolforstudyingmolecularstructures,interactions ] anddynamicsinbio-engineering. Thedataacquisitiontime, Anotherlineofworkisconcernedwiththelowrankness h however,isrelativelylongandnon-uniformsamplingcanbe of MRS signals in the time domain. The time domain p applied to reduce this time. To obtain the full spectrum,a signal of MRS is called the free induction decay (FID), - d reconstruction method with Vandermonde factorization is whichisgenerallyapproximatedbyasumofafewdecaying e proposed.This method explores the general signal property exponential functions. By transforming the FID into a so m in magnetic resonancespectroscopy: Its time domain signal called Hankel matrix, the number of exponentionals will s. is approximated by a sum of a few exponentials. Results equal to the rank of this matrix [6]. Assuming the number c on synthetic and realistic data show that the new approach of spectral peaks is much smaller than the FID data points, si canachievefaithfulspectrumreconstructionandoutperforms the low rank Hankel matrix completion (LRHMC) [5] was y state-of-the-artlowrankHankelmatrixmethod. proposedtorecovertheFIDsignalsintheNUSMRS.Unlike h CS,whichseeksthesparsityofthespectruminthefrequency p IndexTerms— MagneticResonanceSpectroscopy,Van- domain and encounters problems to represent the signals [ dermondefactorization,Hankelmatrix,Lowrank with fast decay, LRHMC tried to minimize the number of 1 peaks. ExperimentalresultsonsimulatedandrealMRSdata v 1. INTRODUCTION show that broad peaks can be recovered by LRHMC much 7 better than the l norm minimization on the spectrum [5]. 1 1 0 Magnetic resonance spectroscopy (MRS) has been regarded The recovery condition of low rank Hankel matrix under 7 as an indispensable tool in studying a molecular structure uniformlyrandom samplingand Gaussian randomencoding 0 and dynamics or interactions of biopolymers in chemistry canbefoundin[7]and[8],respectively. . 1 and biology. However, the duration of a multi-dimensional In LRHMC, the connection between exponentials and 0 MRS experimentis proportionalto the numberof measured matrixstructuresisstillnotfullyexplored. Forexample,itis 7 datapointsandincreasesrapidlywithspectralresolutionand 1 stillunclearhowthesignalsubspaceisrelatedtoeachspecific dimensionality. The non-uniformsampling (NUS) approach : exponential function. Taking this into account, we propose v offers a general solution for a dramatic reduction in mea- to use the Vandermonde structure of the Hankel matrix i X surement time. Reconstructing the full signal from a non- to reconstruct the FID. These exponentials are explicitly uniformlysampledsignalisessentialforthenextstepofdata r presented in the Vandermondedecompositionof the Hankel a analysis. Thereconstructionmaybesuccessfulbyexploiting matrix. This nice property allows enforcing more matrix the inherent structure of the signal in the time or frequency structures in the signal model thus has great potential to domains. improvethe spectrumreconstruction. Experimentresultson One line of work is concernedwith the sparsity of MRS syntheticdataandrealMRSdatashowthatthenewapproach in the frequency domain. Compressed sensing (CS) [1] requires significantly fewer measurements than LRHMC to suggests that if the signal enjoys a sparse representation in achieveafaithfulreconstructionoftheFID. some transform domain, it is possible to recover a signal The rest of this paper is organized as follows. Section even when the number of samples is far below its ambient 2 introduces related work. Section 3 presents the proposed dimension.CShasbeendemonstratedasaneffectivetoolfor Hankel matrix completion with Vandermonde factorization reconstructingNUSspectrumbyassumingthatthespectrum (HVaF)method. Section4presentsthenumericalresultson is sparse in the frequencydomain[2, 3, 4]. However,broad bothsimulatedandreal-worlddata. Section5concludesthis *Correspondingauthor:XiaoboQu([email protected]) workanddiscussessomefuturework. 2. RELATEDWORK Note that we choose the dimension of y to be odd simply to yield a squared matrix Hy. Actually, our algorithm and LetRbeaHankeloperatorwhichmapsavectorx ∈ Cn to resultsdonotrelyonthedimensionbeingodd. a Hankel matrix Rx ∈ Cn1×n2 with n1 +n2 = n+1 as In this paper, we formulate the signal reconstruction as follows Hankel matrix completion with Vandermonde factorization. Specifically, we aim to find x, the Hankel matrix of which [Rx] =x , ∀i∈{0,...,n −1}, j ∈{0,...,n −1}, ij i+j 1 2 canbefactorizedintotwoVandermondematrices,i.e.,Hx= Inparticular,wedenotetheHankeloperatorbyHinsteadof UVH with U and V Vandermonde matrices. To force Rinthecasen =n . U and V to be Vandermonde matrices is equal to restrict 1 2 The LRHMC [5] is based on the observations that the each column of U and V to contain single component of HankelmatrixRxconstructedbytheFIDxislowrankifthe exponentialfunction. AccordingtoKroneckerstheorem[9], numberofspectralpeaksismuchsmallerthanthedatapoints we can equally constrain the Hankel matrix of each column in the whole spectrum. Hence, the reconstruction problem ofUandVtoberank1. Thereforetheoptimization(3)can canbeformulatedasthelowrankmatrixcompletionproblem beequivalentlyrewrittenas min kHxk + λky−Dxk2, (1) find x,U,V (4) x ∗ 2 2 s.t. rank(R(U ))=1, rank(R(V ))=1, where x is the acquired NUS FID data, D is an operator (:,r) (:,r) of the NUS schedule, k·k∗ is the nuclear norm defined as a Hx=UVH, y=Dx, forallr ∈{1,...,K}. sumofmatrixsingularvalues.TheefficiencyofLRHMChas Giventhatitisdifficulttodevelopareliableandcompu- beenverifiedonnumericalsimulationsandrealMRSdata[5]. tationalalgorithmtosolvetherank-constrainedproblemand However,wewillpresentinSection3thattheHankelmatrix themeasurementsinMRS experimentsareusuallycontami- of the FID has Vandermonde factorization and experiment natedbyboundednoise,werelax(4)andsolvethefollowing results in Section 4 show that the new approach exploiting probleminstead Vandermonde factorization can achieve much better recon- structionthanLRHMCfromthesameNUSdata. K λ min (kRQ Uk +kRQ Vk )+ ky−Dxk2. (5) U,V,xX r ∗ r ∗ 2 2 3. HANKELMATRIXCOMPLETIONWITH r=1 VANDERMONDEFACTORIZATION s.t.Hx=UVH. ThegeneralsignalpropertyinMRSthatFIDcanbeapprox- To solve (5), we develop an algorithm based on half imated by a sum of a few decaying exponentials has been quadratic methods with continuation for its advantage in widely acknowledged[5, 6]. Let vector y ∈ C2N−1 be the handling multivariable optimization [10]. We introduce the completeFID term Rx−UVH 2 to keep Rx close enough to UVH (cid:13) (cid:13)F instea(cid:13)d of addressin(cid:13)g the constraint Rx = UVH directly, R y = d e−k∆t/τr+i2πfrk∆t, (2) andproposethefollowingoptimization, k X r r=1 K where d , τ and f are the complex amplitude, decay time min (kRQ Uk +kRQ Vk )+ β Hx−UVH 2 andfreqruenrcy,resperctively,ofthek-thexponential. U,V,xXr=1 r ∗ r ∗ 2 (cid:13)(cid:13) (cid:13)(cid:13)F Letdefinez = e−∆t/τr+i2πfr∆t. Itisobservedthatthe λ HankelmatrixorftheFIDy∈C2N−1admitsaVandermonde + 2 ky−Dxk22. factorization (6) 1 ··· 1 Whenβ → ∞, the solutionto (6)is approachingto (5). To Hy= z..1 ··..· z..R c1 ... s(6o)l,vaen(d6)t,hseonm(6e)aiusxrielfioarrymvualaritaebdlienstoarteheinftorolldouwciendgienqtouimvaoldenelt . . . zN−1 ··· zN−1 cR (3) form: 1 R 11... zz...1 ········· zz1NN...−−11 . Um,Vin,xXrK=1(kBrλk∗+kCrk∗)+ β2 (cid:13)(cid:13)Hx−UVH(cid:13)(cid:13)2F (7) R R + ky−Dxk2, 2 2 Obviously,(3)canbeeasilyrewrittenasaformoftheproduct oftwofactormatricesandthispaperfocusonthelatterform. s.t. B =RQ U, C =RQ V. r r r r In (7), the first two terms are non-smoothbut separable, andtheothertermsaresmooth,whichmakesitrelativeeasier Ground than (6) in developingnumericalalgorithms. For a givenβ, wesolve(7)byAlternatingDirectionMethodofMultipliers (ADMM)[5]. The details of the new algorithm are omitted hereduetolimitedspace. 4. NUMERICALEXPERIMENTS In this section, we will evaluate the performance of the proposed HVaF on synthetic data and real MRS. The state- of-the-art LRHMC [5] is compared with HVaF. The NUS is Poisson-gap sampling [11]. We empirically observe that HVaFis veryrobustto K and thuswe set K to be 64 in all (a) theexperimentsofthispaper. 4.1. Syntheticdata Fig.1showscomparisonsbetweenasimulatedfullysampled referencespectrumanditsNUSreconstructionsobtainedus- ingtheLRHMCandHVaFalgorithms.Thesimulatedsignals with 127pointsare generatedaccordingto (2). The relative leastnormalizederror(RLNE)isdefinedbykx−yk /kyk , 2 2 wherexandyarethereconstructedsignalandtruesignal. It is observed that when sampling rate is relatively high, both LRHMC and HVaF can achieve a reconstruction with low RLNE;however,assamplingratedecreases,HVaFobtainsa significantlylowerRLNEthanLRHMC.ThusHVaFrequires (b) fewer NUS samples to achieve faithful reconstructions than LRHMC. In particular, the line shape of the low-intensity Fig.1: Reconstructionsofthesyntheticspectrumcontaining peakcanbepreservedbetterbyHVaFthanLRHMC. five peaks. Fig. 1(a) The reconstructions by LRHMC and HVaF from 25% NUS. Fig. 1(b) The reconstruction RLNE 4.2. RealisticMRSdata withrespecttodifferentsamplingrates. Theerrorbarsarethe standard deviations of the RLNE over 100 NUS resampling HereweapplytheproposedHVaFtorecovera2-Dspectrum trials. withthesize512×255fromtheNUSFID.Fig2showsaNUS 2D 1H-15N HSQC spectrum of the intrinsically disordered cytosolic domain of human CD79b protein from the B-cell extra assumptions on the signal compared with LRHMC. receptor (More details about the MRS data and experiment Themethodallowsthemorereductioninmeasurementtime setting can be foundin [5]). Obviously, the HVaF produces andachievesmorefaithfulreconstructionthanLRHMC.The the more faithful recovery (Fig. 2(c)) of the ground truth Vandermondefactorization method will also be extended to (Fig. 2(a)) than LRHMC (Fig. 2(b)). To clearly compare thehigh-dimensionalMRSreconstruction[12]. the reconstructedresults, we presentone of the slices in the reconstruction. Fig. 2(d)illustratesthatseverallow-intensity peaks are notably compromised in the LRHMC spectrum, 6. ACKNOWLEDGMENTS whileHVaFsucceedstorecovermostofthepeaks. This work was supported by National Natural Science 5. CONCLUSIONANDDISCUSSION FoundationofChina (61571380,61302174andU1632274), Fundamental Research Funds for the Central Universities We develop the HVaF reconstruction as a new technique (20720150109), Natural Science Foundation of Fujian to obtain high-quality spectra from non-uniform sampling Province of China (2015J01346, 2016J05205), Important measurements. HVaFexplorestheVandermondestructureof Joint Research Project on Major Diseases of Xiamen City theHankelmatrixconstructedbytheFIDinsteadofthelow (3502Z20149032). The authorsthank Vladislav and Maxim rankstructureinLRHMC.NotethatHVaFdoesnotintroduce forsharingtheMRSdatausedinpaper[5]. 15N 15N (ppm) Ground Truth (ppm) HVaF LRHMC 110 110 115 115 120 120 (cid:52)(cid:82)(cid:85)(cid:84)(cid:72) (cid:40)(cid:54)(cid:65)(cid:38) (cid:44)(cid:50)(cid:40)(cid:45) 125 125 8.6 8.4 8.2 8 8.6 8.4 8.2 8 8.6 8.4 8.2 81 H (ppm) (a) (b) (c) (d) Fig. 2: The LRHMC and HVaF reconstructions in the 2D 1H-15N HSQC spectrum experiment. (a) the uniformly-sampled spectrum;(b)and(c)andaretheLRHMCandHVaFreconstructionsusing17%sampleddata,respectively.(d)thesliceofthe reconstructionslocatedatthe8.25ppminthedimensionof1H. 7. REFERENCES [7] Y. Chen and Y. Chi, “Robust spectral compressed sensingviastructuredmatrixcompletion,”IEEETrans- [1] E. J. Cande`s, J. Romberg, and T. Tao, “Robust actions on Information Theory, vol. 60, no. 10, pp. uncertainty principles: exact signal reconstruction 6576–6601,2014. from highly incomplete frequency information,” IEEE [8] J.-F.Cai,X.Qu,W.Xu,andG.-B.Ye,“Robustrecovery TransactionsonInformationTheory,vol.52,no.2,pp. of complex exponential signals from random gaussian 489–509,2006. projectionsvialowrankHankelmatrixreconstruction,” [2] X. Qu, X. Cao, D. Guo, and Z. Chen, “Compressed AppliedandComputationalHarmonicAnalysis,vol.41, sensing for sparse magnetic resonance spectroscopy,” no.2,pp.470–490,2016. in International Society for Magnetic Resonance in [9] F. Andersson, M. Carlsson, J.-Y. Tourneret, and Medicine 18th Scientific Meeting, 2010, Conference H. Wendt, “A new frequency estimation method for Proceedings,p.3371. equallyandunequallyspaceddata,”IEEETransactions on Signal Processing, vol. 62, no. 21, pp. 5761–5774, [3] K. Kazimierczuk and V. Y. Orekhov, “Accelerated 2014. NMR spectroscopy by using compressed sensing,” Angewandte Chemie-International Edition, vol. 50, [10] X. Qu, Y. Hou, F. Lam, D. Guo, J. Zhong, and no.24,pp.5556–5559,2011. Z. Chen, “Magnetic resonance image reconstruction from undersampledmeasurementsusing a patch-based [4] X. Qu, D. Guo, X. Cao, S. Cai, and Z. Chen, nonlocal operator,” Medical Image Analysis, vol. 18, “Reconstruction of self-sparse 2D NMR spectra from no.6,pp.843–856,2014. undersampleddata inthe indirectdimension,”Sensors, vol.11,no.9,pp.8888–8909,2011. [11] S. G. Hyberts, K. Takeuchi, and G. Wagner, “Poisson- gap sampling and forward maximum entropy recon- [5] X.Qu, M.Mayzel,J.-F.Cai, Z.Chen,andV.Orekhov, structionforenhancingtheresolutionandsensitivityof “Accelerated NMR spectroscopy with low-rank recon- protein nmr data,” Journal of the American Chemical struction,” Angewandte Chemie-International Edition, Society,vol.132,no.7,pp.2145–2147,2010. vol.54,no.3,pp.852–854,2015. [12] J. Ying, H. Lu, Q. Wei, J.-F. Cai, D. Guo, J. Wu, [6] H. M. Nguyen, X. Peng, M. N. Do, and Z.-P. Liang, Z. Chenc, and X. Qu, “Hankel matrix nuclear “Denoising MR spectroscopic imaging datawith low- norm regularized tensor completion for N-dimensional rank approximations,” IEEE Transactions on Biomedi- exponential signals,” https://arxiv.org/abs/1604.02100, calEngineering,vol.60,no.1,pp.78–89,2013. 2016.