Symmetric Positive Definite Matrices Research Articles

Robustness of lognormal confidence regions for means of symmetric positive definite matrices when applied to mixtures of lognormal distributions

Robustness of lognormal confidence regions for means of symmetric positive definite matrices when applied to mixtures of lognormal distributions

Riemannian Geometry of Functional Connectivity Matrices for Multi-Site Attention-Deficit/Hyperactivity Disorder Data Harmonization.

The use of multi-site datasets in neuroimaging provides neuroscientists with more statistical power to perform their analyses. However, it has been shown that the imaging-site introduces variability in the data that cannot be attributed to biological sources. In this work, we show that functional connectivity matrices derived from resting-state multi-site data contain a significant imaging-site bias. To this aim, we exploited the fact that functional connectivity matrices belong to the manifold of symmetric positive-definite (SPD) matrices, making it possible to operate on them with Riemannian geometry. We hereby propose a geometry-aware harmonization approach, Rigid Log-Euclidean Translation, that accounts for this site bias. Moreover, we adapted other Riemannian-geometric methods designed for other domain adaptation tasks and compared them to our proposal. Based on our results, Rigid Log-Euclidean Translation of multi-site functional connectivity matrices seems to be among the studied methods the most suitable in a clinical setting. This represents an advance with respect to previous functional connectivity data harmonization approaches, which do not respect the geometric constraints imposed by the underlying structure of the manifold. In particular, when applying our proposed method to data from the ADHD-200 dataset, a multi-site dataset built for the study of attention-deficit/hyperactivity disorder, we obtained results that display a remarkable correlation with established pathophysiological findings and, therefore, represent a substantial improvement when compared to the non-harmonization analysis. Thus, we present evidence supporting that harmonization should be extended to other functional neuroimaging datasets and provide a simple geometric method to address it.

EEG emotion recognition based on enhanced SPD matrix and manifold dimensionality reduction

Recently, Riemannian geometry-based pattern recognition has been widely employed to brain computer interface (BCI) researches, providing new idea for emotion recognition based on electroencephalogram (EEG) signals. Although the symmetric positive definite (SPD) matrix manifold constructed from the traditional covariance matrix contains large amount of spatial information, these methods do not perform well to classify and recognize emotions, and the high dimensionality problem still unsolved. Therefore, this paper proposes a new strategy for EEG emotion recognition utilizing Riemannian geometry with the aim of achieving better classification performance. The emotional EEG signals of 32 healthy subjects were from an open-source dataset (DEAP). The wavelet packets were first applied to extract the time-frequency features of the EEG signals, and then the features were used to construct the enhanced SPD matrix. A supervised dimensionality reduction algorithm was then designed on the Riemannian manifold to reduce the high dimensionality of the SPD matrices, gather samples of the same labels together, and separate samples of different labels as much as possible. Finally, the samples were mapped to the tangent space, and the K-nearest neighbors (KNN), Random Forest (RF) and Support Vector Machine (SVM) method were employed for classification. The proposed method achieved an average accuracy of 91.86%, 91.84% on the valence and arousal recognition tasks. Furthermore, we also obtained the superior accuracy of 86.71% on the four-class recognition task, demonstrated the superiority over state-of-the-art emotion recognition methods.

Empirical Bayes identification of stationary processes and approximation of Toeplitz spectra

We study Statistical Consistency of an approximate subspace identification procedure for the infinite dimensional a posteriori model of a frequency estimation problem in an Empirical Bayesian framework. By first imposing a natural uniform prior probability density on the unknown frequencies, the estimation of the hyperparameters of the a priori distribution can be accomplished by a sequence of subspace identification techniques. These techniques exploit the special structure of the covariance matrix of the a posteriori process which is discovered by making a connection with classical results on energy concentration in deterministic signal processing. The convergence of the spectrum of the subspace estimates to the (nonrational) spectral density of the a posteriori process has an analytic counterpart in the approximation of symmetric positive definite Toeplitz matrices by submatrices of finite rank. This is proven by a weak-sense convergence theorem for Toeplitz spectra.

A Robust Algebraic Domain Decomposition Preconditioner for Sparse Normal Equations

Solving the normal equations corresponding to large sparse linear least-squares problems is an important and challenging problem. For very large problems, an iterative solver is needed and, in general, a preconditioner is required to achieve good convergence. In recent years, a number of preconditioners have been proposed. These are largely serial, and reported results demonstrate that none of the commonly used preconditioners for the normal equations matrix is capable of solving all sparse least-squares problems. Our interest is thus in designing new preconditioners for the normal equations that are efficient and robust and can be implemented in parallel. Our proposed preconditioners can be constructed efficiently and algebraically without any knowledge of the problem and without any assumption on the least-squares matrix except that it is sparse. We exploit the structure of the symmetric positive definite normal equations matrix and use the concept of algebraic local symmetric positive semidefinite splittings to introduce two-level Schwarz preconditioners for least-squares problems. The condition number of the preconditioned normal equations matrix is shown to be theoretically bounded independently of the number of subdomains in the splitting. This upper bound can be adjusted using a single parameter $\tau$ that the user can specify. We discuss how the new preconditioners can be implemented on top of the PETSc library using only 150 lines of Fortran, C, or Python code. Problems arising from practical applications are used to compare the performance of the proposed new preconditioner with that of other preconditioners.

Optimizing and Factorizing the Wilson Matrix

The Wilson matrix, W, is a 4 × 4 unimodular symmetric positive definite matrix of integers that has been used as a test matrix since the 1940’s, owing to its mild ill-conditioning. We ask how close W is to being the most ill-conditioned matrix in its class, with or without the requirement of positive definiteness. By exploiting the matrix adjugate and applying various matrix norm bounds from the literature we derive bounds on the condition numbers for the two cases and we compare them with the optimal condition numbers found by exhaustive search. We also investigate the existence of factorizations with Z having integer or rational entries. Drawing on recent research that links the existence of these factorizations to number-theoretic considerations of quadratic forms, we show that W has an integer factor Z and two rational factors, up to signed permutations. This little matrix continues to be a useful example on which to apply existing matrix theory as well as being capable of raising challenging questions that lead to new results.

Exponential and related probability distributions on symmetric matrices

Pursuing the study initiated in Hassairi and Roula (2019), we show in the present paper that the reliability function of a probability distribution on the cone Ω of positive definite symmetric matrices characterizes the distribution without any invariance condition. We also show that the characterization of the exponential probability distribution on Ω by a memoryless property holds without assuming an invariance condition. We then study the connection between the exponential distribution on Ω and the uniform distribution on a bounded interval of Ω. A notion of matrix Pareto distribution is introduced, and it is shown that this distribution possesses the long tail property.

Accelerating Geostatistical Modeling and Prediction With Mixed-Precision Computations: A High-Productivity Approach With PaRSEC

Geostatistical modeling, one of the prime motivating applications for exascale computing, is a technique for predicting desired quantities from geographically distributed data, based on statistical models and optimization of parameters. Spatial data are assumed to possess properties of stationarity or non-stationarity via a kernel fitted to a covariance matrix. A primary workhorse of stationary spatial statistics is Gaussian maximum log-likelihood estimation (MLE), whose central data structure is a dense, symmetric positive definite covariance matrix of the dimension of the number of correlated observations. Two essential operations in MLE are the application of the inverse and evaluation of the determinant of the covariance matrix. These can be rendered through the Cholesky decomposition and triangular solution. In this contribution, we reduce the precision of weakly correlated locations to single- or half- precision based on distance. We thus exploit mathematical structure to migrate MLE to a three-precision approximation that takes advantage of contemporary architectures offering BLAS3-like operations in a single instruction that are extremely fast for reduced precision. We illustrate application-expected accuracy worthy of double-precision from a majority half-precision computation, in a context where uniform single-precision is by itself insufficient. In tackling the complexity and imbalance caused by the mixing of three precisions, we deploy the <monospace>PaRSEC</monospace> runtime system. <monospace>PaRSEC</monospace> delivers on-demand casting of precisions while orchestrating tasks and data movement in a multi-GPU distributed-memory environment within a tile-based Cholesky factorization. Application-expected accuracy is maintained while achieving up to <inline-formula><tex-math notation="LaTeX">$1.59X$</tex-math></inline-formula> by mixing FP64/FP32 operations on 1536 nodes of <monospace>HAWK</monospace> or 4096 nodes of <monospace>Shaheen II</monospace> , and up to <inline-formula><tex-math notation="LaTeX">$2.64X$</tex-math></inline-formula> by mixing FP64/FP32/FP16 operations on 128 nodes of <monospace>Summit</monospace> , relative to FP64-only operations. This translates into up to 4.5, 4.7, and 9.1 (mixed) PFlop/s sustained performance, respectively, demonstrating a synergistic combination of exascale architecture, dynamic runtime software, and algorithmic adaptation applied to challenging environmental problems.

Finite-difference methods for solving 1D Poisson problem

The paper discusses the formulation and analysis of methods for solving the one-dimensional Poisson equation based on finite-difference approximations - an important and very useful tool for the numerical study of differential equations. In fact, this is a classical approximation method based on the expansion of the solution in a Taylor series, based on which the recent progress of theoretical and practical studies allowed increasing the accuracy, stability, and convergence of methods for solving differential equations. Some of the features of this analysis include interesting extensions to classical numerical analysis of initial and boundary value problems. In the first part, a numerical method for solving the one-dimensional Poisson equation is presented, which reduces to solving a system of linear algebraic equations (SLAE) with a banded symmetric positive definite matrix. The well-known tridiagonal matrix algorithm, also known as the Thomas algorithm, is used to solve the SLAEs. The second part presents a solution method based on an analytical representation of the exact inverse matrix of a discretized version of the Poisson equation. Expressions for inverse matrices essentially depend on the types of boundary conditions in the original setting. Variants of inverse matrices for the Poisson equation with different boundary conditions at the ends of the interval under study are presented - the Dirichlet conditions at both ends of the interval, the Dirichlet conditions at one of the ends and Neumann conditions at the other. In all three cases, the coefficients of the inverse matrices are easily found and the algorithm for solving the problem is practically reduced to multiplying the matrix by the vector of the right-hand side.

Learning a discriminative SPD manifold neural network for image set classification

Performing pattern analysis over the symmetric positive definite (SPD) manifold requires specific mathematical computations, characterizing the non-Euclidian property of the involved data points and learning tasks, such as the image set classification problem. Accompanied with the advanced neural networking techniques, several architectures for processing the SPD matrices have recently been studied to obtain fine-grained structured representations. However, existing approaches are challenged by the diversely changing appearance of the data points, begging the question of how to learn invariant representations for improved performance with supportive theories. Therefore, this paper designs two Riemannian operation modules for SPD manifold neural network. Specifically, a Riemannian batch regularization (RBR) layer is firstly proposed for the purpose of training a discriminative manifold-to-manifold transforming network with a novelly-designed metric learning regularization term. The second module realizes the Riemannian pooling operation with geometric computations on the Riemannian manifolds, notably the Riemannian barycenter, metric learning, and Riemannian optimization. Extensive experiments on five benchmarking datasets show the efficacy of the proposed approach.

The finite cell method with least squares stabilized Nitsche boundary conditions

We apply the recently developed least squares stabilized symmetric Nitsche method for enforcement of Dirichlet boundary conditions to the finite cell method. The least squares stabilized Nitsche method in combination with finite cell stabilization leads to a symmetric positive definite stiffness matrix and relies only on elementwise stabilization, which does not lead to additional fill in. We prove a priori error estimates and bounds on the condition numbers.

The geometry of mixed-Euclidean metrics on symmetric positive definite matrices

Several Riemannian metrics and families of Riemannian metrics were defined on the manifold of Symmetric Positive Definite (SPD) matrices. Firstly, we formalize a common general process to define families of metrics: the principle of deformed metrics. We relate the recently introduced family of alpha-Procrustes metrics to the general class of mean kernel metrics by providing a sufficient condition under which elements of the former belong to the latter. Secondly, we focus on the principle of balanced bilinear forms that we recently introduced. We give a new sufficient condition under which the balanced bilinear form is a metric. It allows us to introduce the Mixed-Euclidean (ME) metrics which generalize the Mixed-Power-Euclidean (MPE) metrics. We unveal their link with the (u,v)-divergences and the (α,β)-divergences of information geometry and we provide an explicit formula of the Riemann curvature tensor. We show that the sectional curvature of all ME metrics can take negative values and we show experimentally that the sectional curvature of all MPE metrics but the log-Euclidean, power-Euclidean and power-affine metrics can take positive values.

Entries of the inverses of large positive definite Toeplitz matrices

AbstractThis is an expository paper embarking on the asymptotic behavior of the entries of the inverses of positive definite symmetric Toeplitz matrices as the matrix dimension goes to infinity. We consider the behavior of the entries in neighborhoods of the four corners as well as the density of the distribution of the entries over all of the inverse matrix.

Pipelined Preconditioned Conjugate Gradient Methods for real and complex linear systems for distributed memory architectures

Preconditioned Conjugate Gradient (PCG) is a popular method for solving large and sparse linear systems of equations. The performance of PCG at scale is affected due to the costly global synchronization steps that arise in dot-products on distributed memory systems. Pipelined PCG (PIPECG) removes the costly global synchronization steps from PCG by only executing a single non-blocking allreduce per iteration and overlapping it with independent computations.In our previous work, we have developed a novel pipelined PCG algorithm called PIPECG-OATI (One Allreduce per Two Iterations) for real linear systems which executes a single non-blocking allreduce per two iterations and provides a large overlap of global communication with independent computations at higher number of cores. Our method achieves this overlap by using iteration combination and by introducing new recurrence and non-recurrence computations. We implement optimizations in the PIPECG-OATI method to use cache memory efficiently.In this work, we present PIPECG-OATI-c method for linear systems with complex Hermitian positive definite and complex symmetric matrices. We compare our method with various pipelined CG methods on a variety of problems and demonstrate that our method always gives the least run times. We performed experiments with our method using 20M and 30M unknowns on up to 16K cores and obtained up to 2.48X performance improvement over PCG and 2.14X performance improvement over PIPECG methods. We also experimented with up to 1-billion unknowns on 16K cores, the largest problem size explored for the CG problem, to our knowledge, and obtained about 25% improvement over PCG.

A new type of CGO solutions and its applications in corner scattering

We consider corner scattering for the operator ∇ ⋅ γ(x)∇ + k 2 ρ(x) in , with γ a positive definite symmetric matrix and ρ a positive scalar function. A corner is referred to one that is on the boundary of the (compact) support of γ(x) − I or ρ(x) − 1, where I stands for the identity matrix. We assume that γ is a scalar function in a small neighborhood of the corner. We show that any admissible incident field will be scattered by such corners, which are allowed to be concave. Moreover, we provide a brief discussion on the existence of non-scattering waves when γ − I has a jump across the corner. In order to prove the results, we construct a new type of complex geometric optics solutions.

A Non-Iterative Method for the Difference of Means on the Lie Group of Symmetric Positive-Definite Matrices

A non-iterative method for the difference of means is presented to calculate the log-Euclidean distance between a symmetric positive-definite matrix and the mean matrix on the Lie group of symmetric positive-definite matrices. Although affine-invariant Riemannian metrics have a perfect theoretical framework and avoid the drawbacks of the Euclidean inner product, their complex formulas also lead to sophisticated and time-consuming algorithms. To make up for this limitation, log-Euclidean metrics with simpler formulas and faster calculations are employed in this manuscript. Our new approach is to transform a symmetric positive-definite matrix into a symmetric matrix via logarithmic maps, and then to transform the results back to the Lie group through exponential maps. Moreover, the present method does not need to compute the mean matrix and retains the usual Euclidean operations in the domain of matrix logarithms. In addition, for some randomly generated positive-definite matrices, the method is compared using experiments with that induced by the classical affine-invariant Riemannian metric. Finally, our proposed method is applied to denoise the point clouds with high density noise via the K-means clustering algorithm.

Seismic Data Denoising With Correlation Feature Optimization Via S-Mean

Random noise elimination acts as an important role in seismic data processing. Moreover, protecting and recovering useful subsurface structure information is also significant. In this study, the S-mean which can obtain the geometric mean of the seismic traces on the symmetric positive definite matrix (SPD) manifold is adopted as a nonlinear filter for seismic denoising. Furthermore, S-mean has the best correlation with other elements based on the S-divergence due to the optimization of finding S-mean on SPD manifold. Therefore, the broken correlation features in noisy seismic data are compensated and maintained well, which can be conducive to describe the subsurface structures. Synthetic example and field data applications qualitatively and quantitatively demonstrate the validity and effectiveness of the proposed workflow.

Epileptic Seizure Detection in EEG Signals Using Discriminative Stein Kernel-Based Sparse Representation

The automatic seizure detection in electroencephalogram (EEG) signals is crucial for the monitoring, diagnosis, and treatment of epilepsy. In this study, an intelligent detection framework with the discriminative Stein kernel-based sparse representation (DSK-SR) is constructed to distinguish epileptic EEG signals. Specifically, in the scheme of DSK-SR, EEG samples are presented by symmetric positive definite (SPD) matrices in the form of covariance descriptors (CovDs). Taking into account the non-Euclidean geometry of the Riemannian manifold of SPD matrices, the traditional SR in Euclidean space cannot be applied in its original form on the manifold. To this end, the DSK defined on the manifold can permit us to embed the manifold into a high-dimensional reproducing kernel Hilbert space (RKHS) to perform SR. Then, test samples are sparsely coded over the training sets, and the classification decision is performed by assessing which class generates the minimal reconstructed residuals. Eventually, competitive experimental results on three widely recognized EEG datasets against the state-of-the-art methods demonstrate the efficacy of the proposed DSK-SR in identifying epileptic EEG signals, indicating its powerful application potential in the automatic seizure detection.

Maximum margin Riemannian manifold-based hyperdisk for fault diagnosis of roller bearing with multi-channel fusion covariance matrix

For rotating machinery, the sudden failure of roller bearing would lead to the downtime of the whole system and even catastrophic accidents. Therefore, multiple accelerometers are usually arranged to comprehensively evaluate the health of roller bearing, enhancing the stability and reliability of monitoring results. This paper proposes a novel fault diagnosis framework by utilizing a multi-channel fusion covariance matrix (MFCM) and Riemannian manifold-based hyperdisk. First, 22 statistical features are acquired from each channel data. Then, MFCM is calculated as the fault feature representation of roller bearing to achieve multi-channel feature fusion, where the element of MFCM represents the correlation information between different channels. Finally, since MFCM is a symmetric positive definite (SPD) matrix, lying on a Riemannian manifold, we design a maximum margin Riemannian manifold-based hyperdisk (MMRMHD) classifier to conduct fault classification, where Log-Euclidean metric (LEM) is introduced to calibrate the distribution of MFCMs. Moreover, to further improve the classification ability of nonlinear SPD data, we map MFCMs into a high-dimensional Hilbert space with the LEM-based kernel function and construct a novel kernelized MMRMHD model. The experimental results on two bearing datasets with multi-channel vibration signals demonstrate the effectiveness and superiority of the proposed fault diagnosis framework.

Computationally Efficient Approximations for Matrix-Based Rényi's Entropy

The recently developed matrix based Renyi's entropy enables measurement of information in data simply using the eigenspectrum of symmetric positive semi definite (PSD) matrices in reproducing kernel Hilbert space, without estimation of the underlying data distribution. This intriguing property makes the new information measurement widely adopted in multiple statistical inference and learning tasks. However, the computation of such quantity involves the trace operator on a PSD matrix $G$ to power $\alpha$(i.e., $tr(G^\alpha)$), with a normal complexity of nearly $O(n^3)$, which severely hampers its practical usage when the number of samples (i.e., $n$) is large. In this work, we present computationally efficient approximations to this new entropy functional that can reduce its complexity to even significantly less than $O(n^2)$. To this end, we leverage the recent progress on Randomized Numerical Linear Algebra, developing Taylor, Chebyshev and Lanczos approximations to $tr(G^\alpha)$ for arbitrary values of $\alpha$ by converting it into matrix-vector multiplications problem. We also establish the connection between the matrix-based Renyi's entropy and PSD matrix approximation, which enables exploiting both clustering and block low-rank structure of $G$ to further reduce the computational cost. We theoretically provide approximation accuracy guarantees and illustrate the properties of different approximations. Large-scale experimental evaluations on both synthetic and real-world data corroborate our theoretical findings, showing promising speedup with negligible loss in accuracy.