Square Root Of Matrix Research Articles

Square-root second-order extended Kalman filter and its application in target motion analysis

In this study, a new non-linear filter based on second-order extended Kalman technology and square-root Kalman algorithm is proposed and applied to the problem of target motion analysis (TMA). By orthogonalising the state vectors with the square root of covariance matrix, the authors reduce the computational complexity in the second-order extended Kalman algorithm, and then we deduce a second-order extended Kalman algorithm based on square-root Kalman technology, which is more numerically stable and can guarantee positive and semi-definiteness of the state covariance. The authors apply the new algorithm to a TMA problem, and prove with both numerical simulation and actual experimentation that their algorithm performs well in the significant non-linear systems.

Fast Generalized Eigenvector Tracking Based on the Power Method

This letter develops a power method-based algorithm for tracking generalized eigenvectors when stochastic signals having unknown correlation matrices are observed. The proposed approach is based on the fact that the generalized eigenvalue problem is reduced to a standard eigenvalue problem, for which the power method can be easily applied by changing the metric of the vector space. The difficulty of applying the power method to this problem lies in the computational load of the inverse of the square root of a matrix at every update. The proposed algorithm can avoid the computation of eigenvalue decomposition to obtain the inverse of the matrix square root. Experimental results show that the proposed tracking algorithm gives a similar performance in tracking to the direct computation of singular value decomposition (SVD), while the computation order of the proposed algorithm is lower than the SVD.

Generalized indirect covariance NMR formalism for establishment of multidimensional spin correlations.

Multidimensional nuclear magnetic resonance (NMR) experiments measure spin-spin correlations, which provide important information about bond connectivities and molecular structure. However, direct observation of certain kinds of correlations can be very time-consuming due to limitations in sensitivity and resolution. Covariance NMR derives correlations between spins via the calculation of a (symmetric) covariance matrix, from which a matrix-square root produces a spectrum with enhanced resolution. Recently, the covariance concept has been adopted to the reconstruction of nonsymmetric spectra from pairs of 2D spectra that have a frequency dimension in common. Since the unsymmetric covariance NMR procedure lacks the matrix-square root step, it does not suppress relay effects and thereby may generate false positive signals due to chemical shift degeneracy. A generalized covariance formalism is presented here that embeds unsymmetric covariance processing within the context of the regular covariance transform. It permits the construction of unsymmetric covariance NMR spectra subjected to arbitrary matrix functions, such as the square root, with improved spectral properties. This formalism extends the domain of covariance NMR to include the reconstruction of nonsymmetric NMR spectra at resolutions or sensitivities that are superior to the ones achievable by direct measurements.

Ensemble covariances adaptively localized with ECO-RAP. Part 2: a strategy for the atmosphere

Part 1’s localization method, Ensemble COrrelations Raised to A Power (ECO-RAP), is incorporated into a Local Ensemble Transform Kalman Filter (LETKF). Because brute force incorporation would be too expensive, we demonstrate a factorization property for Part 1’s Covariances Adaptively Localized with ECO-rap (CALECO) forecast error covariance matrix that, together with other simplifications, reduces the cost. The property inexpensively provides a large CALECO ensemble whose covariance is the CALECO matrix. Each member of the CALECO ensemble is an element-wise product between one raw ensemble member and one column of the square root of the ECO-RAP matrix. The LETKF is applied to the CALECO ensemble rather than the raw ensemble. The approach enables the update of large numbers of variables within each observation volume at little additional computational cost. Under plausible assumptions, this makes the CALECO and standard LETKF costs similar. The CALECO LETKF does not require artificial observation error inflation or vertically confined observation volumes both of which confound the assimilation of non-local observations such as satellite observations. Using a 27 member ensemble from a global NumericalWeather Prediction (NWP) system, we depict four-dimensional (4-D) flow-adaptive error covariance localization and test the ability of the CALECO LETKF to reduce analysis error.

Ensemble covariances adaptively localized with ECO-RAP. Part 2: a strategy for the atmosphere

Part 1’s localization method, Ensemble COrrelations Raised to A Power (ECO-RAP), is incorporated into a Local Ensemble Transform Kalman Filter (LETKF). Because brute force incorporation would be too expensive, we demonstrate a factorization property for Part 1’s Covariances Adaptively Localized with ECO-rap (CALECO) forecast error covariance matrix that, together with other simplifications, reduces the cost. The property inexpensively provides a large CALECO ensemble whose covariance is the CALECO matrix. Each member of the CALECO ensemble is an element-wise product between one raw ensemble member and one column of the square root of the ECO-RAP matrix. The LETKF is applied to the CALECO ensemble rather than the raw ensemble. The approach enables the update of large numbers of variables within each observation volume at little additional computational cost. Under plausible assumptions, this makes the CALECO and standard LETKF costs similar. The CALECO LETKF does not require artificial observation error inflation or vertically confined observation volumes both of which confound the assimilation of non-local observations such as satellite observations. Using a 27 member ensemble from a global NumericalWeather Prediction (NWP) system, we depict four-dimensional (4-D) flow-adaptive error covariance localization and test the ability of the CALECO LETKF to reduce analysis error.

Ensemble covariances adaptively localized with ECO-RAP. Part 2: a strategy for the atmosphere

Part 1’s localization method, Ensemble COrrelations Raised to A Power (ECO-RAP), is incorporated into a Local Ensemble Transform Kalman Filter (LETKF). Because brute force incorporation would be too expensive, we demonstrate a factorization property for Part 1’s Covariances Adaptively Localized with ECO-rap (CALECO) forecast error covariance matrix that, together with other simplifications, reduces the cost. The property inexpensively provides a large CALECO ensemble whose covariance is the CALECO matrix. Each member of the CALECO ensemble is an element-wise product between one raw ensemble member and one column of the square root of the ECO-RAP matrix. The LETKF is applied to the CALECO ensemble rather than the raw ensemble. The approach enables the update of large numbers of variables within each observation volume at little additional computational cost. Under plausible assumptions, this makes the CALECO and standard LETKF costs similar. The CALECO LETKF does not require artificial observation error inflation or vertically confined observation volumes both of which confound the assimilation of non-local observations such as satellite observations. Using a 27 member ensemble from a global NumericalWeather Prediction (NWP) system, we depict four-dimensional (4-D) flow-adaptive error covariance localization and test the ability of the CALECO LETKF to reduce analysis error.

No eigenvalues outside the support of the limiting empirical spectral distribution of a separable covariance matrix

We consider a class of matrices of the form C n = ( 1 / N ) A n 1 / 2 X n B n X n ∗ × A n 1 / 2 , where X n is an n × N matrix consisting of i.i.d. standardized complex entries, A n 1 / 2 is a nonnegative definite square root of the nonnegative definite Hermitian matrix A n , and B n is diagonal with nonnegative diagonal entries. Under the assumption that the distributions of the eigenvalues of A n and B n converge to proper probability distributions as n N → c ∈ ( 0 , ∞ ) , the empirical spectral distribution of C n converges a.s. to a non-random limit. We show that, under appropriate conditions on the eigenvalues of A n and B n , with probability 1, there will be no eigenvalues in any closed interval outside the support of the limiting distribution, for sufficiently large n . The problem is motivated by applications in spatio-temporal statistics and wireless communications.

Newton's method with exact line search for the square root of a matrix

We show how to incorporate the line search into Newton's method and simplify solution procedure for obtaining the square root of a matrix. The line search improves the global convergence properties of Newton's method and ameliorates the stability of iteration procedure. Some numerical examples are given to show the effectiveness of the new method.

Unbiased ensemble square root filters

There is an established framework that describes a class of ensemble Kalman filter algorithms as square root filters (SRFs). These algorithms produce analyses by updating a state estimate and a square root of the ensemble covariance matrix. The matrix square root of the forecast covariance is post-multiplied by another matrix to give a matrix square root of the analysis covariance. The choice of post-multiplier is not unique and can be multiplied by any orthogonal matrix to give another scheme that is also a SRF. Not all filters of this type bear the desired relationship to the forecast ensemble: the analysis ensemble mean may not be equal to the analysis state estimate and consequently there may be an accompanying shortfall in the spread of the analysis ensemble as expressed by the ensemble covariance matrix. This points to the need for a restricted version of the notion of an ensemble SRF, which we call an unbiased ensemble SRF. This paper provides a generic set of necessary and sufficient conditions for the scheme to be centred on the analysis state estimate (unbiased). A few of these results have already been published elsewhere in the literature; this paper brings these together with new results and provides simple proofs and examples, as well as a mathematical description of the set of unbiased ensemble SRFs. Many (but not all) published ensemble SRF algorithms satisfy the criteria that we establish. While these conditions are not a cure-all and cannot deal with independent sources of bias such as model and observation errors, they should be useful to designers of ensemble SRFs in the future.

3-D Full-Vectorial Analysis of Strong Optical Waveguide Discontinuities Using Pade Approximants

A full-vectorial 3-D numerical method applicable to high-index contrast optical waveguide discontinuities is presented. Rigorous treatment of the longitudinal boundary condition is incorporated in the formulation. The square root of the characteristic matrix is approximated using Padeacute approximants which results in an efficient implementation. The biconjugate gradient stabilized method is utilized to iteratively calculate the reflected and transmitted fields. A preconditioner is proposed which results in reduced number of iterations. The proposed method is applied to various optical waveguide facets exhibiting strong transverse and longitudinal refractive index discontinuities. In all cases, the modal reflectivities of the fundamental TE-Like and TM-Like modes are calculated for both the full-vectorial and the semi-vectorial formulations. Significant difference in the calculated modal reflectivity is seen between the full and semi-vectorial models. The error in the power balance remains low in the full-vectorial case irrespective of the waveguide dimensions. However, in the semi-vectorial case, the error in the power balance is found to increase when the waveguide width is reduced

Fast adaptive LDA using quasi-Newton algorithm

A new adaptive algorithm for linear discriminant analysis (LDA) based on the quasi-Newton optimization technique is presented. The proposed algorithm uses the secant method for adaptive computation of the inverse Hessian matrix and the Newton–Raphson method for optimal estimation of the step size at each iteration. Current adaptive method, based on the Newton–Raphson optimization technique, uses a direct calculation of the inverse Hessian matrix, which can be both laborious to calculate and invert for systems with large number of dimensions. The new algorithm has the advantage of automatic optimal selection of the step size using the current data samples and also adaptive computation of the inverse Hessian matrix that overcomes its sensitivity to data condition. Based on the new adaptive algorithm, we present a self-organizing neural network for adaptive computation of the square root of the inverse covariance matrix and use it in cascaded form with a principal component analysis (PCA) network for LDA. Experimental results demonstrated fast convergence and lower computational cost of the new algorithm compared to the adaptive gradient descent and Newton–Raphson LDA algorithms, respectively and justified its advantages for on-line pattern recognition applications with stationary and non-stationary multidimensional input data.

Computing the square roots of matrices with central symmetry

For computing square roots of a nonsingular matrix A, which are functions of A, two well known fast and stable algorithms, which are based on the Schur decomposition of A, were proposed by Björk and Hammarling [A. Björk, S. Hammarling, A Schur method for the square root of a matrix, Linear Alg. Appl. 52/53 (1983) 127–140], for square roots of general complex matrices, and by Higham [N.J. Higham, Computing real square roots of a real matrix, Linear Alg. Appl. 88/89 (1987) 405–430], for real square roots of real matrices. In this paper we further consider (the computation of) the square roots of matrices with central symmetry. We first investigate the structure of the square roots of these matrices and then develop several algorithms for computing the square roots. We show that our algorithms ensure significant savings in computational costs as compared to the use of standard algorithms for arbitrary matrices.

An efficient algorithm for double‐difference tomography and location in heterogeneous media, with an application to the Kilauea volcano

Improving our understanding of crustal processes requires a better knowledge of the geometry and the position of geological bodies. In this study we have designed a method based upon double‐difference relocation and tomography to image, as accurately as possible, a heterogeneous medium containing seismogenic objects. Our approach consisted not only of incorporating double difference in tomography but also partly in revisiting tomographic schemes for choosing accurate and stable numerical strategies, adapted to the use of cross‐spectral time delays. We used a finite difference solution to the eikonal equation for travel time computation and a Tarantola‐Valette approach for both the classical and double‐difference three‐dimensional tomographic inversion to find accurate earthquake locations and seismic velocity estimates. We estimated efficiently the square root of the inverse model's covariance matrix in the case of a Gaussian correlation function. It allows the use of correlation length and a priori model variance criteria to determine the optimal solution. Double‐difference relocation of similar earthquakes is performed in the optimal velocity model, making absolute and relative locations less biased by the velocity model. Double‐difference tomography is achieved by using high‐accuracy time delay measurements. These algorithms have been applied to earthquake data recorded in the vicinity of Kilauea and Mauna Loa volcanoes for imaging the volcanic structures. Stable and detailed velocity models are obtained: the regional tomography unambiguously highlights the structure of the island of Hawaii and the double‐difference tomography shows a detailed image of the southern Kilauea caldera–upper east rift zone magmatic complex.

Iterative refinement method for the approximate factorization of a matrix inverse

An iterative refinement technique for the sparse approximate inverse factor $Z$ of a symmetric matrix $S$, where $Z{Z}^{T}\ensuremath{\approx}{S}^{\ensuremath{-}1}$ is proposed. Its efficiency is demonstrated for the approximate inverse square root and Cholesky factorization of the overlap matrix $S$ in the generalized eigenvalue problem occurring in linear scaling electronic structure theory using a localized, nonorthogonal basis set. A version of the refinement method useful for local perturbations of the overlap matrix is also proposed. In addition to its applications in electronic structure theory, the scheme may be useful for iterative refinement of preconditioners.

Novel Finite Element Analysis of Optical Waveguide Discontinuity Problems

In this paper, a novel finite-element method (FEM) to rigorously and efficiently solve the optical waveguide discontinuity problems is presented. Instead of performing the time-consuming modal solutions on both discontinuity sides, the square root of the characteristic matrix is efficiently approximated using Taylor's series expansion, and then the interface boundary conditions are enforced at the discontinuity plane to solve for the reflected and transmitted fields. The high numerical precision and effectiveness of the proposed method is demonstrated through the analysis of various discontinuity problems, and the excellent agreement of the results obtained using the present finite element method and those obtained using other rigorous approaches in the literature.

CLT for linear spectral statistics of large-dimensional sample covariance matrices

Let $B_n=(1/N)T_n^{1/2}X_nX_n^*T_n^{1/2}$ where $X_n=(X_{ij})$ is $n\times N$ with i.i.d. complex standardized entries having finite fourth moment, and $T_n^{1/2}$ is a Hermitian square root of the nonnegative definite Hermitian matrix $T_n$. The limiting behavior, as $n\to\infty$ with $n/N$ approaching a positive constant, of functionals of the eigenvalues of $B_n$, where each is given equal weight, is studied. Due to the limiting behavior of the empirical spectral distribution of $B_n$, it is known that these linear spectral statistics converges a.s. to a nonrandom quantity. This paper shows their rate of convergence to be $1/n$ by proving, after proper scaling, that they form a tight sequence. Moreover, if $\expp X^2_{11}=0$ and $\expp|X_{11}|^4=2$, or if $X_{11}$ and $T_n$ are real and $\expp X_{11}^4=3$, they are shown to have Gaussian limits.

Computing the square root and logarithm of a real P-orthogonal matrix

A real square matrix A is P-orthogonal if A T PA= P; P is a fixed real nonsingular matrix, but most of the results in this work require that it is symmetric positive definite or P T =P −1, P 2=±I . The class of P-orthogonal matrices includes, for instance, orthogonal and symplectic matrices as particular cases. We present an efficient iterative method for computing the P-orthogonal factor in the generalized polar decomposition, which generalizes the well-known Newton's method for the standard polar decomposition. A connection between Newton's method for the matrix square root and polar iterates brings out a new iterative method for computing the principal square root of a P-orthogonal matrix. One important feature of this method is that, when P is symmetric positive definite, it allows us to restore the P-orthogonal property of the exact square root by computing the nearest P-orthogonal matrix. We also analyse the problem of finding the nearest P-symmetric and P-skew-symmetric matrices. New bounds and new estimates for the Padé error of the matrix logarithm are given in order to improve the existing Briggs–Padé algorithms and adapt them to P-orthogonal matrices. Special attention will be paid to the orthogonal case.

Accelerated Stokesian dynamics: Brownian motion

A new Stokesian dynamics (SD) algorithm for Brownian suspensions is presented. The implementation is based on the recently developed accelerated Stokesian dynamics (ASD) simulation method [Sierou and Brady, J. Fluid Mech. 448, 115 (2001)] for non-Brownian particles. As in ASD, the many-body long-range hydrodynamic interactions are computed using fast Fourier transforms, and the resistance matrix is inverted iteratively, in order to keep the computational cost O(N log N). A fast method for computing the Brownian forces acting on the particles is applied by splitting them into near- and far-field contributions to avoid the O(N3) computation of the square root of the full resistance matrix. For the near-field part, representing the forces as a sum of pairwise contributions reduces the cost to O(N); and for the far-field part, a Chebyshev polynomial approximation for the inverse of the square root of the mobility matrix results in an O(N1.25 log N) computational cost. The overall scaling of the method is thus roughly of O(N1.25 log N) and makes possible the simulation of large systems, which are necessary for studying long-time dynamical properties and/or polydispersity effects in colloidal dispersions. In this work the method is applied to study the rheology of concentrated colloidal suspensions, and results are compared with conventional SD. Also, a faster approximate method is presented and its accuracy discussed.

Fast adaptive algorithms and networks for class-separability features

In this article, we introduce accelerated algorithms for linear discriminant analysis (LDA) and feature extraction from unimodal multiclass Gaussian data. Current adaptive methods based on the gradient descent optimization technique use a fixed or a monotonically decreasing step size in each iteration, which results in a slow convergence rate. Here, we use a variable step size, optimally computed in each iteration using the steepest descent method, in order to accelerate the convergence of the algorithm. Based on the new adaptive algorithm, we present a self-organizing neural network for adaptive computation of the square root of the inverse covariance matrix ( Σ −1/2) and use it (i) in a network for optimal feature extraction from Gaussian data and (ii) in cascaded form with a principal component analysis network for LDA. Experimental results demonstrate fast convergence and high stability of the algorithm and justify its advantages for on-line pattern recognition applications with stationary and non-stationary input data.

Using SEM Programs to Perform Matrix Manipulations and Data Simulation

A didactic discussion of a simple method for Cholesky decomposition and obtaining a symmetric square root of a positive definite covariance matrix using popular structural equation modeling programs is discussed. The applicability of such software for purposes of accomplishing some matrix manipulations is demonstrated. The approach can also be used to generate multivariate normal data. The outlined procedure is illustrated with a numerical example.