$\ell_{1}$-Regularized ICA: A Novel Method for Analysis of Task-related fMRI Data

In the introduction, we mainly focus on the application of ICA to fMRI. In practice, ICA is applied not only to biological data analysis but also to other various problems such as image processing or finance. Therefore, we formulate our proposed ICA in general setting for the application to various problems.

Originally, ICA is proposed as one of the blind source separation methods. The objective of all ICA methods is to estimate latent independent source $\bm{S}=[\bm{s}_{1},\ldots,\bm{s}_{K}]^{T}$ and mixing matrix $\bm{A}$ simultaneously without any prior knowledge on them. In many algorithms of ICA, non-gaussianity is the quantity to be maximized.

Among various ICA methods, FastICA [17] is a widely used algorithm. In FastICA, factorized matrices $\bm{A},\bm{S}$ are estimated as follows.

1. 1.

   The observed matrix $\bm{X}$ is prewhitened by $\bm{Q}\in\mathbb{R}^{K\times p}$ as $\bm{Z}=\bm{QX}$, where $\bm{Q}$ denotes the prewhitening matrix and $\bm{Z}$ denotes the prewhitened observed matrix.

2. 2.

   Next, the projection matrix $\bm{W}\in\mathbb{R}^{K\times K}$ is evaluated, which maximizes non-gaussianity among rows of the signal source matrix $\bm{S}=\bm{W}^{T}\bm{Z}$. The projection matrix $\bm{W}$ is frequently called unmixing matrix. The cost function for non-gaussianity is defined by

   $$J_{f}(\bm{w}_{i})=\frac{1}{N^{2}}\left\{\sum_{j=1}^{N}G(\bm{w}_{i}^{T}\bm{z}_{j})-\sum_{j=1}^{N}G(\nu_{j})\right\}^{2}.$$

   (1)

   In the above equation, $G(.)$ is a nonlinear function and defined by $G(y)=-e^{-y^{2}/2}$. The variable $\nu_{j}$ $(j\in\{1,\ldots,N\}$) is a normal random variable with zero mean and unit variance. In maximizing $J_{f}(\bm{w}_{i})$, a fixed point algorithm in equation (2) is used.

   $$\bm{w}_{i}\leftarrow\frac{1}{N}\sum_{j=1}^{N}\bm{z}_{j}G^{\prime}(\bm{w}_{i}^{T}\bm{z}_{j})-G^{\prime\prime}(\bm{w}_{i}^{T}\bm{z}_{j})\bm{w}_{i}.$$

   (2)

   In addition, the orthogonal constraint $\bm{W}^{T}\bm{W}=\bm{I}_{K}$ is imposed on the projection matrix $\bm{W}$, where $\bm{I}_{K}$ is $K$-dimensional identity matrix. In practice, the column vectors of $\bm{W}$ are mutually orthogonalized by Gram-Schmidt method, where the column vectors are updated as follows,

   $$\bm{w}_{i}\leftarrow\bm{w}_{i}-\sum_{j=1}^{i-1}(\bm{w}_{i}^{T}\bm{w}_{j})\bm{w}_{j}.$$

   (3)

3. 3.

   The mixing matrix $\bm{A}$ is calculated as $\bm{A}=\bm{Q}^{\dagger}\bm{W}$, where $\bm{Q}^{\dagger}$ is Moore-Penrose pseudo-inverse matrix of $\bm{Q}$. Accordingly, the observed matrix $\bm{X}$, the mixing matrix $\bm{A}$, and signal source matrix $\bm{S}$ are related by $\bm{X}\approx\bm{A}\bm{S}$, which is shown by the orthogonality of $\bm{W}$.

Before going into our proposed method, we should mention several related works regarding ICA with sparse constraint. First, there is a method called SparseICA [18], [19], which was introduced for blind source separation in image processing in real-world application [20]. In this method, independent source matrix $\bm{S}$ is computed as the product of the sparse coefficient matrix $\bm{C}$ and the dictionary matrix $\bm{\Phi}$ by $\bm{S}=\bm{C}\bm{\Phi}$.

Next, ICA with regularization terms was proposed for the linear non-gaussian acyclic model [21]. In this method, the regularization term for the separation matrix defined by $\bm{W}^{T}\bm{Q}$ is introduced. It should be noted that the construction of sparse matrix in ICA is similar to our proposed method explained in the following section. However, the sparse constraint on the matrix is different. In this method, the sparse constraint is imposed on the separation matrix $\bm{W}^{T}\bm{Q}$, while the estimated mixture matrix $\bm{A}$ is sparse in our method. Similar studies with sparse separation matrix $\bm{W}^{T}\bm{Q}$ can also be found in other references [22], [23].

There is also other related work [24], where the sparse constraint is imposed on the mixing matrix $\bm{A}$ like our method. Additionally, the objective of the work is to identify sparse neuronal networks like ours. However, there is a difference in the methodology: this method is based on Bayesian framework and the algorithm is a stochastic one. More precisely, the mixing matrix $\bm{A}$ and parameters in the model are assumed to follow some probabilistic distributions and optimized by Markov chain Monte Carlo method. In contrast, our method is a deterministic one, where the cost function is minimized by DC algorithm.

In our method, the $i$-th vector in the projection matrix $\bm{w}_{i}$ minimizing cost function in equation (4) must be evaluated.

Here $\|\cdot\|_{p}$ is $\ell_{p}$-norm, $\widetilde{J_{f}}(\bm{w}_{i})=-{J_{f}}(\bm{w}_{i})$, and $\bm{Q}^{\#}\in\mathbb{R}^{p\times K}$ is given by the minimizer in equation (5).

with $\|\cdot\|_{\rm F}$ being Frobenuis norm and $\|\cdot\|_{0}$ being the matrix $\ell_{0}$-norm or the number of non-zero elements in the matrix, respectively. By definition, the minimizer matrix $\bm{Q}^{\#}$ is close to the original $\bm{Q}^{\dagger}$ and sparse with only $k_{0}$ nonzero elements. We need the sparse estimated mixture matrix $\bm{A}=\bm{Q}^{\#}\bm{W}$, whose sparsity depends on the sparsity of $\bm{Q}^{\#}$ [25]. Hence, approximate sparse inverse whitening matrix $\bm{Q}^{\#}$ is favorable rather than the original inverse whitening matrix $\bm{Q}^{\dagger}$ in our proposed method. In practical computation, the algorithm of orthogonal matching pursuit [26], [27] is applied to obtain each row of $\bm{Q}^{\#}$.

We apply widely-used DC algorithm for minimization of the cost function. DC algorithm is applicable to the minimization problem, where the cost functions $J(\bm{w})$ is represented by the difference of two convex functions $g(\bm{w}),h(\bm{w})$ as

Note that Hessians of the two convex functions must be bounded above. In DC algorithm, the following two processes are repeated until convergence or maximum iteration number $M$ is reached.

1. (i)

   Update ${}_{\nabla}\bm{h}^{(t)}$ by choosing one element from the subgradient set of $h(\bm{w})$, namely ${}_{\nabla}\bm{h}^{(t)}\in\partial h(\bm{w}^{(t)})$ with $\bm{w}^{(t)}$ being the $t$-th update of $\bm{w}$.

2. (ii)

   Update $\bm{w}^{(t+1)}$ by $\bm{w}^{(t+1)}=\mathop{\rm arg~{}min}\limits_{\bm{w}}\{g(\bm{w})-({}_{\nabla}\bm{h}^{(t)})^{T}\bm{w}\}$.

It is considered that the convergent solution by DC algorithm is often a global optimal solution. Therefore, DC algorithm is one of the effective approaches for nonconvex optimization.

For application of DC algorithm to minimization in our method, the cost function $J_{o}(\bm{w}_{i})$ must be expressed by the difference of two convex functions. First, Lipschitz constant $L$ of $\widetilde{J_{f}}(\bm{w}_{i})$ is introduced, which can be defined because $\widetilde{J_{f}}(\bm{w}_{i})$ is secondary differentiable. Then, the function $(L/2)\|\bm{w}_{i}\|_{2}^{2}-\widetilde{J_{f}}(\bm{w}_{i})$ is easily proved to be convex by the definition of Lipschitz constant $L$. Moreover, regularization term $\|\bm{Q}^{\#}\bm{w}_{i}\|_{1}$ is clearly convex, and the sum of two convex functions $(L/2)\|\bm{w}_{i}\|_{2}^{2}+\alpha\|\bm{Q}^{\#}\bm{w}_{i}\|_{1}$ is also convex. Therefore, $J_{o}(\bm{w}_{i})$ can be represented by the difference of two convex functions $g(\bm{w}_{i}),h(\bm{w}_{i})$ like $J_{o}(\bm{w}_{i})=g(\bm{w}_{i})-h(\bm{w}_{i})$, where the two convex functions are defined in equations (7) and (8).

Note that $g(\bm{w}_{i})$ and $h(\bm{w}_{i})$ are continuous functions in our $\ell_{1}$-regularized ICA.

In DC algorithm, computation of Lipschitz constant $L$ needs maximum eigenvalue of Hessian $\nabla^{2}\widetilde{J_{f}}(\bm{w}_{i})$, whose computational complexity is $O(K^{3})$. In practice, $L$ is evaluated by backtracking line search algorithm [28]. In this algorithm, we prepare the tentative Lipschitz constant $l^{(t)}$ in the $t$-th update, then $l^{(t)}$ is accepted as the final Lipschitz constant if the criterion of inequality (9) is satisfied for $\zeta\in(0,1)$. Otherwise, $l^{(t)}$ is rescaled by $l^{(t)}\leftarrow\eta l^{(t)}$ with the constant $\eta>1$, then inequality (9) is checked again.

After evaluation of Lipschitz constant, we move on to the formulation of DC algorithm in our proposed method. In DC algorithm, subderivative $\partial h(\bm{w}_{i})$ is necessary for the calculation of ${}_{\nabla}\bm{h}^{(t)}$ in the process (i), which is given by equation (10),

Then, the process (ii) for the update of $\bm{w}_{i}^{(t)}$ is represented as follows.

The minimization problem in equation (11) can be regarded as a generalized lasso problem [29], whose solution cannot be obtained analytically. Hence, we attempt to find the solution numerically by the method of alternating directions method of multipliers (ADMM) [30]. In the application of ADMM, $\bm{w}_{i}^{(t+1)}$ is computed by minimizing augmented Lagrangian $\mathcal{L}_{\rho}(\bm{w},\bm{\gamma},\bm{\tau})$,

The variables $\bm{w},\bm{\gamma},\bm{\tau}$ minimizing $\mathcal{L}_{\rho}(\bm{w},\bm{\gamma},\bm{\tau})$ are alternately updated by equations (13), (14) and (15) until convergence.

where superscript $[\theta]$ means the variable of the $\theta$-th update in ADMM iteration. Update of $\bm{w}^{[\theta+1]}$ can be rewritten by putting the solution of $\partial\mathcal{L}_{\rho}(\bm{w},\bm{\gamma}^{[\theta]},\bm{\tau}^{[\theta]})/\partial\bm{w}=\bm{0}$, which leads to

Similarly, by solving equation $\partial\mathcal{L}_{\rho}(\bm{w}^{[\theta+1]},\bm{\gamma},\bm{\tau}^{[\theta]})/\partial\bm{\gamma}=\bm{0}$ analytically, update of the $j$-th element in $\bm{\gamma}^{[\theta+1]}$ is represented as follows,

where $(\bm{Q}^{\#})_{j,}$ means the $j$-th row vector of $\bm{Q}^{\#}$ and $\mathrm{S}_{a}(x)$ is soft thresholding function defined in equation (18),

Finally, by combining these results, our proposed algorithm is obtained, which is summarized in Algorithm 1.

We can derive a convergence condition for our proposed algorithm for $\ell_{1}$-regularized ICA. The theorem and the proof of the convergence condition are given in the following.

For general DC algorithm, the condition of non-increasing function $J_{o}(\bm{w}_{i}^{(t+1)})\leq J_{o}(\bm{w}_{i}^{(t)})$ must be satisfied for the updated variable $\bm{w}_{i}^{(t+1)}$ [31]. When this condition is satisfied, it is guaranteed that the sequence $\{\bm{w}^{(t)}\}_{t\in\mathbb{Z}_{\geq 0}}$ converges to a point $\bm{w}^{*}$ satisfying the condition for subgradient as

The significant difference of DC algorithm in our $\ell_{1}$-regularized ICA is orthogonalization and normalization after variable update, which requires a slight change in the proof for convergence condition of DC algorithm. For this purpose, we should distinguish the variable after orthogonalization and normalization, which is denoted by $\widehat{\bm{w}}_{i}^{(t+1)}$. With this variable, we can show that our proposed algorithm converges to a point in equation (19) if the inequality $J_{o}(\widehat{\bm{w}}_{i}^{(t+1)})\leq J_{o}(\bm{w}_{i}^{(t)})$ is satisfied.

The convergence theorem for our proposed algorithm is stated as follows. It should be commented that the convergence condition in the theorem is a sufficient condition, which can be relaxed.

For validity of our proposed method, we apply our method to synthetic data to evaluate its performance. In addition, we compare the performance of our method with widely-used FastICA algorithm. In our numerical experiments, we use scikit-learn library (version 1.2.2) for Python (version 3.10.11).