Sparse Hierachical Extrapolated Parametric Methods for Cortical Data Analysis

TL;DR

This paper introduces a novel approach that leverages the hierarchical structure of icosahedral meshes in cortical data analysis, significantly accelerating data factorization methods like sparse dictionary learning.

Contribution

It demonstrates how to exploit the hierarchical mesh structure for faster cortical data factorization and compares various methods on multiple neuroimaging datasets.

Findings

Exploiting mesh hierarchy yields several orders of magnitude speedup.

Different sparsity norms and extrapolation methods impact performance.

The approach is validated on structural and functional MRI datasets.

Abstract

Many neuroimaging studies focus on the cortex, in order to benefit from better signal to noise ratios and reduced computational burden. Cortical data are usually projected onto a reference mesh, where subsequent analyses are carried out. Several multiscale approaches have been proposed for analyzing these surface data, such as spherical harmonics and graph wavelets. As far as we know, however, the hierarchical structure of the template icosahedral meshes used by most neuroimaging software has never been exploited for cortical data factorization. In this paper, we demonstrate how the structure of the ubiquitous icosahedral meshes can be exploited by data factorization methods such as sparse dictionary learning, and we assess the optimization speed-up offered by extrapolation methods in this context. By testing different sparsity-inducing norms, extrapolation methods, and factorization schemes, we compare the performances of eleven methods for analyzing four datasets: two structural and two functional MRI datasets obtained by processing the data publicly available for the hundred unrelated subjects of the Human Connectome Project. Our results demonstrate that, depending on the level of details requested, a speedup of several orders of magnitudes can be obtained.

Equations

$$X\approx DBC$$

$$L=X^{T}D,K=D^{T}D,M=L^{T}L$$

$$(DL)$$

$$S(z,\alpha {)}_{ij}=sign(z_{ij})\max \left(0,|z_{ij}|-\alpha \right)$$

$$(PNNMF)$$

$$(PPNMF)\min ||X-DB{B}^T{D}^{T}X|{|}_2^{2}B\ge 0$$

$$\text{by repeating}B\leftarrow B\odot \left(\frac{1(n_k,n_d)}{2}+\left[MB\right]\oslash \left[\left(KB{B}^{T}M+MB{B}^{T}K\right)B\right]\right)$$

$$e=\left[\left[L^T-KBC\right]\odot \left[L^T-KBC\right]\right]1(n_s,1)$$

$${\Delta }^{+}$$

$${\Delta }^{-}$$

$$B^{t+1}$$

$$y^t$$

$$B^{t+1}$$

$$\tau$$

$$y^t$$

$$B^{t+1}$$

$$\tau$$

$$\displaystyle f(x)=\left\{\begin{array}[]{ll}exp\left(-\frac{{cos}^{-1}(\langle x_{o},x\rangle)}{\pi\sigma}\right)&~{}~{}~{}~{}~{}\frac{{cos}^{-1}(\langle x_{o},x\rangle)}{\pi\sigma}\leq\tau\\ 0&~{}~{}~{}~{}~{}\texttt{otherwise}\end{array}\right.$$

$$G=G^{+}-G^{-}{G}^{+}\ge 0,G^{-}\ge 0$$

$$Y_{i,j}\leftarrow Y_{i,j}-\alpha \left(G_{i,j}^{+}-G_{i,j}^{-}\right)$$

$${\alpha }_{i,j}=\frac{Y_{i,j}}{G_{i,j}^{+}}$$

$$Y_{i,j}\leftarrow \frac{Y_{i,j}}{G_{i,j}^{+}}\left(G_{i,j}^{+}-G_{i,j}^{+}+G_{i,j}^{-}\right)=\frac{Y_{i,j}}{G_{i,j}^{+}}{G}_{i,j}^{-}$$

$$Y\leftarrow Y\odot G^{-}\oslash G^{+}$$

$$||Y|{|}_2^2=⟨Y,Y⟩$$

$$L$$

$$K$$

$$M$$

$$(DL)$$

$$e=\frac{1}{2}\left[⟨X,X⟩+⟨DBC,DBC⟩-2⟨X,DBC⟩\right]$$

$$\frac{\partial e(B+\tau F)}{\partial \tau }{|}_{\tau =0}$$

$$B\leftarrow B-\eta \left(KDBC{C}^T-L^T{C}^T\right)$$

$$\frac{\partial e(C+\tau F)}{\partial \tau }{|}_{\tau =0}$$

$$(DL)$$

$$(PPNMF)\min ||X-DB{B}^T{D}^{T}X|{|}_2^{2}B\ge 0$$

$$\frac{\partial e(B+\tau F)}{\partial \tau }{|}_{\tau =0}$$

Full text

1 Introduction

Several modalities, such as EEG and MEG, are not able to image deep brain structures. MRI modalities benefit from better signal-to-noise ratios at the surface of the brain. For these practical reasons, many neuroimaging studies focus on the cortex. Cortical data are usually processed indepedently for the two hemispheres. The data of each hemisphere are projected onto a reference mesh, where subsequent analyses are carried out. Several multiscale approaches have been proposed for analyzing these surface data, such as spherical harmonics and graph wavelets [9, 14, 6]. As far as we know, however, these tools have not been exploited for accelerating data factorization schemes, such as nonnegative matrix factorization [8, 13, 10] and sparse dictionary learning [1, 2].

In this paper, we accelerate cortical data factorization by exploiting the hierarchical structure of icosahedral meshes and by investigating the effects of a novel extrapolation scheme adapted to nonnegative factorizations. Our results demonstrate that, depending on the level of details requested, a speedup of several orders of magnitudes can be obtained. The remainder of the paper is organized as follows. In section 2 we present the four factorization scheme considered for this work, we explain how mesh structure was exploited, how the factorizations were initialized and how they were accelerated by extrapolation. Section 3 presents our experimental results, obtained with two structural and functional datasets generated from the data available for the hundred unrelated HCP subjects [11]. Discussions conclude the paper.

2 Methods

2.1 Factorizations

In this work, we compare four factorizing schemes: a variant of sparse dictionary learning [2], two Non-negative Matrix factorizations (NNMF) [8, 3] and a projected Non-negative Matrix Factorization (PNMF) [13]. Our goal is to decompose a data matrix $X$ of size $n_{f}\times n_{s}$, containing $n_{f}$ positive measurements acquired for $n_{s}$ subjects, as a product $BC$ between $n_{d}$ basis vectors, stored in a matrix $B$ of size $n_{f}\times n_{d}$ and loadings $C$ of size $n_{d}\times n_{s}$. We assume that the locations of the $n_{f}$ measures are in bijection with the faces of an icosahedral mesh ${\cal M}$ generated from the icosahedron ${\cal M}o$ as explained in the next section. In the context of neuroimaging the number of subjects $n_{s}$ is usually of the order of a few hundreds. The dimension $n_{f}$ is generally much larger, for instance 327680 faces for the largest HCP and freesurfer cortical meshes [11, 4]. This large dimension significantly slows down the computations and multiplies the local minima which could trap the alternating minimization scheme commonly used for solving factorization problems.

In this work, we propose to reduce the spatial dimension by introducing a positive design matrix $D$ of size $n_{f}\times n_{k}$ and decomposing $X$ as follows:

[TABLE]

$D$ can be interpreted as an additional set of positive cortical basis. We start with a restricted number of coarse cortical maps, which we gradually refine to let the optimization focus on the cortical regions where decomposition errors are large.

For sparse dictionary learning, the L1 norm of the matrices $B$ and $C$ is penalized at the same time as the L2 decomposition error. Starting from a random initialization, $B$ and $C$ are iteratively updated by an alternating proximal gradient descent [2]. Under the following notations

[TABLE]

the parametric dictionary learning problem solved and the alternating minimization scheme, known as PALM [2], write:

[TABLE]

where $\lambda$ is a constant sparsity parameter, we set $\eta=10^{-1}/||L||_{2}$, and the proximal operator $S$ applies a soft thresholding to matrix components:

[TABLE]

NNMF frameworks impose a positivity constraint on $B$ and $C$ and are usually optimized via multiplicative updates [8]. This constraint is often sufficient for generating sparse and non overlapping basis [8, 13, 3]. We compared two “parametric” NNMF frameworks: a framework where the square of the Frobenius norm of $B$ and $C$ is penalized to alleviate the ambiguity problem [3] and a framework generating sparse $B$ and $C$ [3]:

[TABLE]

where $\odot$ denotes the entrywise product, $\oslash$ the entrywise division, $\left[x\right]_{+}=max(x,0)$ and $\mathbf{1}(n,m)$ is the matrix of ones of size as $n\times m$.

In addition, we considered a projective NMF (PNMF) scheme [13]. PNMF generates the loadings by projecting the data on the basis. As a result, only $B$ needs to be determined and no penalization is required for handling ambiguity issues. We obtained the following parametric scheme after reducing the amplitude of the updates by two for stabilizing the optimization [3]:

[TABLE]

2.2 Hierarchical Optimization on Icosahedral Meshes

Icosahedral meshes are generated by iteratively subdividing the faces of an icosahedron in four smaller triangles, as illustrated in figure 1. This procedure generates almost perfectly regular meshes. Spherical harmonics and graph wavelets allow to exploit the hierarchical structure of these meshes [9, 14, 6]. In order to preserve positivity when factorizing cortical data, we adopted an approach more primitive but closely related to spherical wavelets [9]. More precisely, we initialize $D$ by concatenating twenty rotated versions of the same positive cortical map, centered on the faces of the original icosahedron. A thousand time, we initialize $B$ and $C$ randomly and run a thousand iterations of our factorization method. The best pair $(B,C)$ is then gradually refined by a procedure preserving the association between the column of D and the faces of a icosahedral mesh gradually refined. Each refinement consists of two steps. An average local error is first computed for each column of D by projecting the reconstruction error $X-DBC$ using $D$:

[TABLE]

Then, the faces associated with the worst errors are subdivided, D is updated by rotating the removed columns and scaling their support by half accordingly, the matrices $K,L,M$ required for the factorization are updated, and $B$ is refined by dividing by four and replicating four time the rows corresponding to the columns removed from D. Local errors $e$ and $B$ updates are not exact because $D$ is not an orthogonal basis. However, these updates are very efficient and constitutes good approaximations when the overlap between the cortical maps in $D$ is limited.

2.3 Initialization and Design Matrices

For dictionary learning, $B$ and $C$ were initialized by sampling from ${\cal N}(0,1)$. For PPNMF schemes, $B$ was then scaled by the inverse of the Frobenius norm of $L$ and replaced by its absolute value. For NMF schemes, $C$ was initialized uniformly with $\frac{1}{n_{d}}\mathbf{1}(n_{d},n_{s})$ and $B$ was generated by averaging random selections of five colums of $L$[7]. The following positive functions

[TABLE]

where $\langle.,.\rangle$ is the standard inner product, were evaluated for $(\sigma=0.015,\tau=3)$, for each of the $n_{f}$ face centers $x$ of ${\cal M}$ and for $x_{o}$ center of the twenty original ${\cal M}o$ icosahedron faces. We built an initialization $D_{O}$ for $D$ by concatenating these twenty cortical maps, which means that optimization started at $n_{k}=20$.

2.4 Extrapolation

Extrapolation was introduced for accelerating the convergence of iterative soft-thresholding algorithms [1]. Extrapolation exploits previous gradient steps for extending the current one, as shown in figure 2. For NMF schemes, an additional projection is necessary to preserve positivity. However, this projection generates zeros which can trap the optimization. For this reason, we investigated the log extrapolation illustrated in figure 2. This novel extrapolation, which was bounded in amplitude and run after ten standard updates to prevent instability, corresponds to an extrapolation of the logarithm of the matrix components.

3 Results

3.1 Data and Parameters

We validated our methods using the data of the hundred unrelated HCP subjects [11]. The outliers of each subject myelin and cortical thickness maps were removed by limiting the absolute difference with the median to 4.4478 median absolute deviations. We generated a regional homogeneity map (reHo) [15] for each subject by first bandpass-filtering between 0.05 and 0.1 Hz the rs-fMRI processed with the ICA+FIX pipeline with MSMAll registration [5]. The timeseries obtained were then normalized to zero mean and unit L2 norm, concatenated, and reHo was measured for neighborhoods of three edges of radius [15]. The amplitude of low frequency fluctuations (ALFF) [16] was on the contrary measured for each subject scan separately, for the frequency band 0.05-0.1 Hz. All these positive maps were projected onto the fsaverage5 mesh using the transformation provided on Caret website [12]. For all the experiments sparsity level $\lambda$ was set to $5$ for dictionary learning, $1/2$ for sparse NMF schemes. $\lambda$ was set to $1/||L||_{2}$ for alleviating ambuiguity issues with non-sparse NMF schemes.

3.2 Computational Time and Extrapolation

We measured the maximal speed up by comparing the computational time required for running the first million of optimization steps with the small design matrix $D_{O}$ with the computational time required for running the same algorithm at the original fsaverage5 resolution. The results presented in figure 3 correspond to speedups between $2.4\times 10^{3}$ and $6.25\times 10^{5}$ for the most time consuming projected NMF schemes. These speedups of three to six orders of magnitude grant us the possibility to generate a good initial pair $(B,C)$ by running a large number of random initializations.

We compared the extrapolation strategies by measuring the reconstruction error $||X-DBC||_{2}^{2}$ when running our eleven algorithms ten times for factorizing the left hemisphere myelin data. $D$ was set to $D_{O}$ and the algorithms were run for a thousand iterations. The median errors reported in figure 4 demonstrate that both extrapolation approaches significantly speed up the convergence. However, the log-extrapolation is less likely to get trapped in local minima, and always reaches a slightly lower energy than state of the art extrapolation.

3.3 Methods Comparison

Figure 5 presents the L1 norm of B and C obtained for all our algorithm and datasets. Ten steps of refinement were conducted. At each step, five faces were subdivided and the algorithms run for ten optimization steps. Our results suggest that the tradeoff between the sparsity of B and C is different for the parameters selected. Because projected NMF do not control the sparsity of the loadings, PNMF basis tend to be very sparse but the projected loadings are not. The other factorization schemes balanced the L1 norms of B and C. For the parameters selected, the basis generated by dictionary learning were slightly sparser. We illustrate in figure 6 the results obtained when decomposing the myelin data using LE-PNNMF for a larger number of iterations. The basis obtained nicely decompose the map of large data variability into weakly overlapping components. The refinement had focused accordingly.

4 Discussions

In this paper we exploit the structure of the icosahedral meshes commonly used in neuroimaging for accelerating optimization tasks such as data factorization. We compare four factorizations schemes and investigate the use of extrapolation for further reducing computational time. Our experiments with structural and functional data acquired by the Human Connectome Project demonstrate that our approach is particularly interesting for processing massive datasets.

Appendix A Proofs

All the optimization schemes implement an alternative gradient descent. They alternatively minimize the basis $B$ and the loadings $C$ by performing a gradient descent, which depends on the scheme. For the Dictionary Learning, a step of proximal gradient is performed [2, 1] as explained in the next section.

For the non-negative schemes, a multiplicative update is performed [8]. This update adapts the step size of the gradient for each matrix component independently, so that the matrix components remain positive. More precisely, let $Y$ denote $B$ or $C$ and $G$ the gradient of the differentiable part of the objective function. $G$ is first decomposed into its positive and negative parts:

[TABLE]

A standard gradient update with positive step size $\alpha$ would be, for any component $Y_{i,j}$:

[TABLE]

For nonnegative factorization, the positive step size $\alpha$ is set independently for each component as follows [8]:

[TABLE]

which yields

[TABLE]

This update maintains the positivity of $Y$ and can be expressed in a more elegant form using componentwise products and divisions:

[TABLE]

All the gradients were derived using Gâteaux derivatives, as explained in the next sections. We introduce the notation $\langle.,.\rangle$ for the standard inner product, hence:

[TABLE]

We recall the following notations introduced in the paper:

[TABLE]

Appendix B Dictionary Learning

The differenciable part of the objective of the parametric dictionary learning:

[TABLE]

can be expressed as follows:

[TABLE]

As a result, for $B$, the Gâteaux derivative of $e$ in the direction $F$ writes:

[TABLE]

Hence, a gradient descent of step $\eta$ on $B$ writes:

[TABLE]

Applying a componentwise soft-thresholding $S$ by $\eta\lambda$ [2, 1] leads to the optimization step describe in Section 2.1.

The derivation for $C$ is very similar:

[TABLE]

As a result, our parametric dictionary learning writes:

[TABLE]

Appendix C PPNMF

For our parametric PNMF scheme [13]

[TABLE]

the Gâteaux derivative in the direction $F$ writes:

[TABLE]

which leads to the following multiplicative update:

[TABLE]

Shrinking the amplitude of this update by two as follows:

[TABLE]

leads to the update presented in Section 2.1:

[TABLE]

Appendix D PNNMF

For the parametric NMF scheme

[TABLE]

the Gâteaux derivative with respect to $B$ in the direction $F$ writes:

[TABLE]

Following [8, 3], we join the gradient term originating from the penalty with the negative part of the gradient, and we project the components of the matrix back to positive real numbers. These operations yield the following multiplicative update:

[TABLE]

The same derivation yields for $C$:

[TABLE]

hence the update:

[TABLE]

and we obtain the PNNMF scheme presented in Section 2.1.

Appendix E SPNNMF

The derivation of the SPNNMF scheme follows exactly the derivation of the PNNMF. Since $B$ and $C$, the derivative of the L1 norms is a matrix of ones of the same size as the updated matrix.