Joint estimation of low-rank components and connectivity graph in high-dimensional graph signals: Application to brain imaging

Highlights

• We formulate the problem of joint estimation of low-rank components and connectivity graph for distorted high-dimensional graph signal.

• We propose two approaches to this joint estimation problem. One approach applies alternating optimization to solve the related subproblems. Another approach solves the joint estimation directly.

• We provide analysis to understand the impact of inexact graph on the low- rank components estimation accuracy, and we apply the analysis result to inform the choice of regularization parameter when the graph is inexact.

• We perform extensive experiments on synthetic and real data, especially brain imaging (MEG) data, to illustrate and validate our proposed algo- rithms.

Abstract

Given high-dimensional, graph-smooth and grossly-corrupted signals, this paper presents the algorithm to simultaneously estimate the intrinsic low-rank components and the underlying graph from high-dimensional noisy graph signal. We assume that the perturbation on low-rank components is sparse and signals are smooth on unknown underlying graph. The proposed algorithm learns the low-rank components by exploiting estimate of the graph, and refines the graph estimation with learned low-rank components. We propose two solutions to this problem: One applies alternating optimization to solve the subproblems. The other solves the problem directly. Furthermore, we analyze the impact of inexact graph on low-rank components estimation to justify our approach. We conduct extensive experiments on the proposed algorithm with synthetic data and real brain imaging data, Magnetoencephalography (MEG) and compare it with state-of-the-art methods. In particular, our algorithm is applied to estimate the low-rank components for classifying MEG signals evoked when a subject views a face or non-face image. We observe that our proposed algorithm is competitive in estimating low-rank components, adequately capturing intrinsic task-related information in lower-dimensional representation. This leads to improved classification performance. Furthermore, we notice that our estimated graph indicates brain active regions for the visual activity that are consistent with neuroscientific findings.

Introduction

We consider the problem of uncovering intrinsic low-rank components of high-dimensional data. This data has a large number of features per sample. We assume the intrinsic representations of the data lie on some low-dimensional subspace, referred to as low-rank components. Furthermore, we assume the data is defined on the vertices of some underlying graph and changes smoothly between the connected vertices, referred to as graph-smoothness [1], [2]. However, in many problems, the underlying graph is unknown or inexact [3], [4], [5] and needs to be estimated from the data. Therefore, in this work, we propose to apply graph-smoothness property as a bridge to combine graph learning and low-rank components estimation together. As will be discussed, our proposed algorithm estimates low-rank components and learns the underlying graph simultaneously.

High-dimensional data is common in many engineering areas such as image and video processing, biomedical imaging, computer networks, and transportation networks. As a specific application of our proposed algorithm, we are interested in automated analysis of brain imaging data. In many cases, the goal is to find the spatio-temporal neural signature of a task, by performing classification on cortical activation evoked by different stimuli [6], [7], [8], [9], [10], [11]. Common brain imaging techniques are Electroencephalography (EEG) and Magnetoencephalography (MEG). These measurements are high-dimensional spatio-temporal data. Furthermore, the measurements are degraded by various types of noise (e.g., sensor noise, ambient magnetic field noise, etc.) and some noise models are complicated (non-Gaussian). The high-dimensionality and noise limit both speed and accuracy of signal analysis, which may result in unreliable signature modeling for classification. The high-dimensionality of these signals also increases complexity of classifier. Combination of a complex classifier and availability of few data samples (due to time, cost, or study limitations) can easily lead to model overfitting. Thus, for a reliable study of brain imaging data, there is a need for a dimensionality reduction method that ensures inclusion of task-related information.

Note that it has been recognized that there are patterns of anatomical links, statistical dependencies or causal interactions between distinct units within a nervous system [12]. Some techniques have also been developed to estimate this brain connectivity graph [13], [14]. However, this task is complicated. In many cases, the estimated graph may not be accurate. Thus, it is necessary to have a method to refine the graph.

Our work addresses the low-rank components estimation and graph learning simultaneously. Previous works have studied these two problems independently [15], [16]. For low-rank components estimation, several linear or nonlinear methods have been proposed making use of the graph Laplacian of the sample affinity graphs. Specifically, given the high-dimensional datasets that consist of $n$ $p$-dimensional data points, these methods make use of the sample affinity graphs of $n$ vertices to capture the geometry of the data points and aim to preserve the structure. In Laplacian Eigenmap (LE) [17], dimensionality is reduced by making the nodes that are close in the high-dimensional space to be close in the low-dimensional space. Locality Preserving Projections (LPP) [18] is based on LE, but assumes that there is a vector of transformation. The vector is calculated as the results of LE. Local Linear Embedding (LLE) [19] aims to recover the low-dimensional structure using neighboring points relationships, which are calculated from the high-dimensional data. Multidimensional Scaling (MDS) [20] aims to recover the low-dimensional structure by minimizing the difference between the distance of two nodes in high-dimensional space and the distance of these two nodes in low-dimensional space. Isomap [21] is a nonlinear generalization of MDS. These methods can be viewed as kernel PCA on specially constructed graph [22]. In addition, various approaches have been proposed to incorporate spectral graph regularization [23], [24], [25] to uncover the low-rank components. In these works, graphs are fixed in their algorithms, precomputed from the noisy input data. These works do not consider inaccuracy in the graphs. In contrast, our proposed algorithms refine the graph to improve low-rank components estimation at the same time. Some low-rank components estimation algorithms have been applied for brain imaging data. Principal Components Analysis (PCA) and Independent Components Analysis (ICA) [26] are commonly-used algorithms to estimate the latent representations from noisy brain imaging data [27]. Maximum-likelihood estimation based algorithm [28] has been applied to estimate low-rank signals for MEG/EEG. Robust Principal Components Analysis (RPCA) has been applied in fMRI data to improve the classification accuracy [29].

Graph estimation or graph learning has been an active research topic recently [2], [30]. Different algorithms have been proposed based on various assumptions on the graph. The work of Spyrou and Escudero [31] assumes that original data is a four-mode tensor and it applies graph regularized tensor factorization of a four-mode tensor to denoise data. The work of Maretic et al. [5] assumes that graph signals are sparse and signals are a sparse linear combination of some atoms in a structured graph dictionary. We refer their method as Sparse Graph Dictionary method (SGDict) and we compare it with our proposed method in graph estimation accuracy. The work of Chepuri et al. [32] assumes that topology is sparse and is represented with a sparse edge selection vector. In particular, [32] assumes the number of edges is known and learns the connectivities from the data. It is related to [33], which will be introduced in the next paragraph. Note that [32] only estimates edge existence. However, [33] estimates both edge existence and edge weights. The work of Berger et al. [34] assumes that signals are smooth on the graph, which is also similar to [33]. But it controls the node degree in the graph. Compared with our proposed algorithms, this work only considers data with Gaussian noise. The methods proposed by Pavez and Ortega [4], Rabbat [35] and Hu et al. [36] also assume that signal is smooth on the graph. In these works, data, considered as random vectors, is modeled by a Gaussian Markov Random Field (GMRF) distribution, and its corresponding precision matrix is estimated as graph matrix. The works of Kao et al. [37] and Kalofolias [3] combine sparsity and smoothness assumption together and assume that signals are smooth on the sparse graph. The work of Bars et al. [38] assumes signals are smooth on the graph and bandlimited, which imply sparsity of signal representation in the spectral domain. Study on dynamic graph (topology changes with time) estimation is another active research topic. The work of Villafañe-Delgado et al. [39] proposes the dynamic Graph Fourier Transform (DGFT) to deal with dynamic graph signals with different Laplacian matrices. The work of Kalofolias et al. [40] proposes an algorithm to learn the dynamic graphs based on the assumption that changes in graph are smooth along the time. The work of Yamada et al. [41] infers the time-varying graphs based on the assumption that changes in graph are sparse along the time. The work of Shen et al. [42] proposes a data model based on the autoregressive model to model the time-varying graph signals and learns the graph from the data. Another data model is used by Shafipour et al. [43] and Wai et al. [44]. In this model, graph signals are generated by local diffusion dynamics on the graph. Shafipour et al. [43] infer the graph from non-stationary graph signals. Wai et al. [44] apply low-pass graph filter to detect communities from data with Gaussian noise.

The work of Dong et al. [33] is most related to our work. While its focus is to learn the connectivity graph topology, it also estimates some noise-free data of the input data as output. Gaussian noise model and Frobenius norm optimization are employed in [33] to make it suitable for the problem when noise is Gaussian with small amplitudes. As will be discussed, our work assumes another signal model with a low-rank component plus a sparse perturbation. Experiment results suggest that our method outperforms other methods for high-dimensional graph data corrupted by complicated noise, such as brain imaging signals. Furthermore, our work analyzes the impact of inexact graph on low-rank components estimation.

Graph signal processing tools have been applied to some brain imaging tasks [45]. In [46], graph Fourier transform is applied to decompose brain signals into low, medium, and high-frequency components for analysis of functional brain networks properties. The work of Liu et al. [47] uses eigenvectors of graph Laplacian matrix to approximate intrinsic subspace of high-dimensional brain imaging data. They evaluate different brain connectivity estimations to compute the graph. Estimation of the brain connectivity graph using a Gaussian noise model has been proposed in [36]. The work of Hu et al. [48] introduces a graph regression model for learning structural brain connectivity of Alzheimer’s disease from PiB-PET imaging data. The work of Griffa et al. [49] proposes a method to capture the network of spatio-temporal connectivity integrating neural oscillations and anatomical connectivity from diffusion fMRI. The work of Medaglia et al. [50] applies the graph Fourier transform to distill functional brain signals into aligned and liberal components. The work of Guo et al. [51] applies deep convolutional networks on graph signals for brain signal analysis. Besides, some works [52], [53], [54], [55] make use of graph smoothness regularization to learn the latent manifold in magnetic resonance imaging (dMRI) reconstruction application. The goal of reconstruction is to recover high-dimensional or high resolution data from undersampled inputs. Graph signal processing has also been shown to be useful in image compression [56], temperature data [57], wireless sensor data [58] and Internet traffic analysis [59]. A few signal features motivated by graph signal processing have also been proposed [60], [61].

In this paper, we propose algorithms to estimate low-rank components and learn the underlying graph simultaneously. Based on a signal model with sparse perturbation, the proposed algorithm estimates low-rank components of data using graph smoothness assumption and refines the graph with estimated low-rank components. Our specific contributions are:

• We formulate the problem of joint estimation of low-rank components and connectivity graph for distorted high-dimensional graph signal.

• We propose two approaches to this joint estimation problem. One approach applies alternating optimization to solve the related subproblems. Another approach solves the joint estimation directly.

• We provide analysis to understand the impact of inexact graph on the low-rank components estimation accuracy, and we apply the analysis result to inform the choice of regularization parameter when the graph is inexact.

• We perform extensive experiments on synthetic and real data, especially brain imaging (MEG) data, to illustrate and validate our proposed algorithms.

The joint estimation was first formulated and studied in our previous conference paper [62]. This work is an extension of [62]. Compared with [62], we have included substantial new contributions: a new method to address joint estimation, a study of the impact of graph distortion on low-rank component estimation, extensive and new experiments to understand the subtleties of joint estimation, and MEG experiments with consideration of field spreading and relevant frequency bands.

The rest of the paper is organized as follows. The proposed algorithm and its two solutions are presented in Section 2. Then, a detailed analysis on the impact of graph distortion is described in Section 3. Detailed experiment results on synthetic data are provided in Section 4 and results on real brain imaging data as well as MNIST image data are presented in Section 5. Section 6 concludes the work.