Padé Based-Magnetic Resonance Spectroscopy (MRS)

Article type: Research Article

Authors: Belkić, Dževad

Affiliations: Karolinska Institute, Medical Radiation Physics, P.O. Box 260, S-171 76 Stockholm, Sweden

Abstract: This review is on quantum-mechanical signal processing based upon the Padé approximant (PA). We presently link the PA and the Lanczos algorithm to design the Padé-Lanczos approximant (PLA). The PLA is operationalised with the recursive algorithm called the fast Padé transform (FPT) for both parametric and nonparametric estimations of spectra. The FPT for any given power series is defined by the unique quotient of two polynomials. This processor provides a meaningful result even when the original expansion diverges. It can significantly accelerate slowly converging sequences/series. As opposed to a single polynomial, e.g. the fast Fourier transform (FFT), the FPT can analytically continue general functions outside their definition domains. Moreover, we show that the FPT is an efficient solver of generalised eigenproblems e.g. the quantum-mechanical evolution/relaxation matrix U comprised of auto-correlation functions. These generic functions can be either computed theoretically or measured experimentally. Such a concept, put forward as a computational tool, surpasses its initial purpose. Indeed autocorrelation functions represent a veritable alternative formulation of quantum mechanics. This is not just because all the major observables e.g. complete energy spectra, local density of states, quantal rate constants, etc, are expressible through the auto-correlation functions. It is also because these and other observables could be given completely in terms of some appropriate, relatively small informational parts that can be singled out and analysed separately from the unwanted/redundant remainder of the full data set of auto-correlation functions. The needed dimensionality reduction of original large problems treated by the FPT can be achieved by e.g. windowing using the band-limited decimation. Alternatively, as done in this work, the Lanczos tridiagonalisation can be employed yielding sparse Jacobi matrices in terms of the Lanczos coupling parameters {αn,βn} that have their very important physical interpretations. The FPT is naturally ingrained in the Schrodinger picture of quantum mechanics and in the total time-independent Green function for the studied system. This yields a versatile framework for a unified treatment of spectroscopy and collision within signal processing and quantum mechanics. In the quantum-mechanical method of nearest neighbours or tight bindings, we use time signal points {cn} as the only input data to derive the exact analytical expressions for FPT, the continued fractions, the Lanczos polynomials {Pn(ω),Qn(ω)}, the couplings {αn,βn}, the Lanczos state vectors ψn, the total wave function Υ(ω) at any frequency ω and the ‘Hamiltonian’ operator Ωˆ. Within signal processing we also analyse ordinary difference equations, base-quotient concept, iterative relaxations, nonlinear accelerations, rational approximations, regularisations of spurious roots and the methods of moments. The analysis is performed within the context of the so-called harmonic inversion problem, or equivalently, spectral decomposition or quantification problem. This is an inverse/reconstruction problem which is aimed at retrieving uniquely the whole information from the input time signal {cn}, such that the output data contain the number K, position R(ωk), height |dk|, width I(ωk) and phase Arg(dk) of every complex harmonic including all possible degeneracies/multiplicities of normal mode frequencies {ωk} due to overlapping resonances. The FPT has the unique virtues of unequivocally providing all these parameters for the nondegenerate and degenerate spectra with the least computational effort by attaining the needed accuracy via robust and stable signal processing. This could be fully accomplished in only two simple steps: (1) solving a system of linear equations to extract the expansion/prediction coefficients of the characteristic/secular polynomial QK(ω) directly from the cn’s and (2) rooting QK(ω) to obtain all the sought complex frequencies {ωk} (1⩽k⩽K). Solely step (1) is needed in the nonparametric FPT, which is given by the polynomial quotient PK(ω)/Qk(ω). This yields the shape spectrum in the FPT at any ω. In the quantification problem, the parametric FPT is required by carrying out step (2) which gives all the K complex frequencies {ωk}. The corresponding complex amplitudes {dk} are obtained from the two separate explicit expressions for the nondegenerate and degenerate spectrum as the residues of PK(ω)/Qk(ω) taken at one single frequency ω=ωk. This produces the K isolated Lorentzian resonances in a nondegenerate spectrum if the roots {ωk} of Qk(ω) are all distinct (no multiplicities). When some of the roots {ωk} of Qk(ω) coincide i.e. have multiplicities, a degenerate non-Lorentzian spectrum PK(ω)/Qk(ω) is obtained in the FPT which describes all isolated and completely or partially overlapping resonances. By contrast, other parametric methods, e.g. linear predictor (LP) or the Hankel-Lanczos singular value decomposition (HLSVD), are much more computationally demanding as they obtain the set {dk} by solving another system of linear equations using all the elements {ωk} (1⩽k⩽K), where errors in some of the ωk’s might cause severe inaccuracies in the dk’s. Moreover, to arrive at the dk’s the method of filter diagonalisation (FD) performs additional computations for the full state vectors {Υk} from the generalised eigenproblem of the evolution/relaxation matrix i.e. the Hankel data matrix. Relative to the linear FFT, the advantages of the nonlinear FPT are in the significantly increased resolution for the same signal length or in the same resolution for shorter signal lengths. Furthermore, as opposed to other nonlinear methods (FD, etc) that typically undergo wild oscillations before they eventually stabilise, the FPT exhibits a strikingly stable convergence with the increased signal length. This is the present finding via an illustrative spectral analysis of a time signal, which has been experimentally measured via Magnetic Resonance Spectroscopy (MRS) at the magnetic field strength of 7T encoded from the brain of a healthy volunteer. We also analyse unphysical (spurious or extraneous) resonances in parametric estimation of spectra. Spurious peaks stem from noise corruption of the time signal and from the so-called overdetermination problem. If the signal of the length N happens to have less than N/2 peaks, the problem becomes algebraically overdetermined, since there are more equations than unknowns. This leads to singular values associated with false peaks that represent spurious resonances. In the PA, extraneous roots could appear in both the numerator and denominator polynomial. Spurious roots in the denominator polynomial of the PA are undesirable, since they lead to unphysical spikes in the Padé spectrum. Spurious roots in the numerator polynomial of the PA are also unwelcome, since they fill in the valleys with unphysical antiresonances and this destroys the phase minimum as well as the uniqueness feature of the PA. We examine this important problem by using the so-called constrained root reflection, which is an analytical procedure for regularising spurious roots. First, we separate unequivocally the genuine from spurious resonances in the Padé power spectrum, which is itself the Padé-Chebyshev approximant (PCA). Second, the unstable PCA containing diverging and converging exponentials is properly stabilised. This is done by a special root reflection, which reverses the sign of diverging exponentials so that they are relocated on the side of genuine resonances. Such a procedure is accomplished under the constraint that the parameter and the shape spectra of the PA are the same. The ensuing constrained root reflection has its physical justification in preservation of the total energy of the signal via the Parseval identity. The resulting method is called the Padé-Schur approximant (PSA) which, as a stable estimator, possesses only converging exponentials. The PSA describes the Padé power spectrum as the unique ratio of two Schur polynomials whose computed roots are all adequately regularised by being placed on the side of physical resonances. In this way, rather than attempting to eliminate or reduce the noise content from the measured data, as has often been attempted previously with a potential risk of losing weak genuine spectral features, the PSA processes noise as measured together with the physical signal.

Keywords: Padé approximant, continued fractions, Lanczos algorithm, signal processing

Keywords: 30B70, 30E05, 33D45, 41A21, 62Fxx, 65Bxx, 94A12

DOI: 10.3233/JCM-2003-3405

Journal: Journal of Computational Methods in Sciences and Engineering, vol. 3, no. 4, pp. 563-733, 2003