Cross-product Penalized Component Analysis (XCAN)

TL;DR

XCAN is a flexible sparse matrix factorization method that balances variance maximization and structural preservation, enabling versatile data exploration across disciplines.

Contribution

Introduces XCAN, a novel cross-product penalized component analysis method that extends previous PCA frameworks for improved data exploration.

Findings

Effective in diverse data exploration tasks

Balances variance and structure in factorization

Applicable across multiple disciplines

Abstract

Matrix factorization methods are extensively employed to understand complex data. In this paper, we introduce the cross-product penalized component analysis (XCAN), a sparse matrix factorization based on the optimization of a loss function that allows a trade-off between variance maximization and structural preservation. The approach is based on previous developments, notably (i) the Sparse Principal Component Analysis (SPCA) framework based on the LASSO, (ii) extensions of SPCA to constrain both modes of the factorization, like co-clustering or the Penalized Matrix Decomposition (PMD), and (iii) the Group-wise Principal Component Analysis (GPCA) method. The result is a flexible modeling approach that can be used for data exploration in a large variety of problems. We demonstrate its use with applications from different disciplines.

Equations

$$\mathbf{X}=\hat{\mathbf{T}}{\hat{\mathbf{P}}}^T+\mathbf{E},$$

$$\{\hat{\mathbf{P}},\hat{\mathbf{T}}\}=\underset{\mathbf{P},\mathbf{T}}{\mathrm{argmin}}\Vert \mathbf{X}-\mathbf{T}{\mathbf{P}}^T{\Vert }_F^2$$

$${\hat{\mathbf{p}}}^{SL}=\underset{\mathbf{p}}{\mathrm{argmax}}\Vert \mathbf{Xp}{\Vert }_2^{2}s.t.\Vert \mathbf{p}{\Vert }_1\le c,\Vert \mathbf{p}{\Vert }_2^2=1$$

$$\{{\hat{\mathbf{P}}}^{SP},{\hat{\mathbf{Q}}}^{SP}\}=\underset{\mathbf{P},\mathbf{Q}}{\mathrm{argmin}}\Vert \text{X}-\mathbf{XP}{\mathbf{Q}}^T{\Vert }_F^2+{\lambda }_2\sum _{h=1}^H\Vert {\mathbf{p}}_h{\Vert }_2^2+{\lambda }_1\sum _{h=1}^H\Vert {\mathbf{p}}_h{\Vert }_{1}s.t.{\mathbf{Q}}^T\mathbf{Q}=\mathbf{I}$$

$$\{{\hat{\mathbf{p}}}^P,{\hat{\mathbf{u}}}^P\}=\underset{\mathbf{p},\mathbf{u}}{\mathrm{argmax}}{\mathbf{u}}^T\mathbf{Xp}s.t.\Vert \mathbf{u}{\Vert }_1\le c_1\Vert \mathbf{p}{\Vert }_1\le c_2,\Vert \mathbf{u}{\Vert }_2^2\le 1,\Vert \mathbf{p}{\Vert }_2^2\le 1$$

$${\hat{d}}^P=({\hat{\mathbf{u}}}^P{)}^T\mathbf{X}{\hat{\mathbf{p}}}^P$$

$$\mathbf{X}=\mathbf{X}(\mathbf{I}-{\hat{d}}^P{\hat{\mathbf{u}}}^P({\hat{\mathbf{p}}}^P{)}^T)$$

$$\{{\hat{\mathbf{P}}}^C,{\hat{\mathbf{T}}}^C\}=\underset{\mathbf{P},\mathbf{T}}{\mathrm{argmin}}\Vert \text{X}-\text{T}{\mathbf{P}}^T{\Vert }_F^2+\lambda \sum _{h=1}^H\Vert {\mathbf{p}}_h{\Vert }_1+\lambda \sum _{h=1}^H\Vert {\mathbf{t}}_h{\Vert }_1$$

$$\mathbf{X}={\hat{\mathbf{U}}}^X{\hat{\mathbf{S}}}^X({\hat{\mathbf{P}}}^X{)}^T+\mathbf{E}$$

$$\{{\hat{\mathbf{P}}}^X,{\hat{\mathbf{S}}}^X,{\hat{\mathbf{U}}}^X\}=\underset{\mathbf{P},\mathbf{S},\mathbf{U}}{\mathrm{argmin}}\Vert \text{X}-\text{U}\mathbf{S}{\mathbf{P}}^T{\Vert }_F^2+{\lambda }_r{F}_r+{\lambda }_c{F}_{c}s.t.\Vert {\mathbf{u}}_h{\Vert }_2^2\le 1,\Vert {\mathbf{p}}_h{\Vert }_2^2\le 1$$

$$F_r=\sum _{h=1}^H\Vert ({\mathbf{u}}_h{\mathbf{u}}_h^T)\oslash \mathbf{XXt}{\Vert }_F^2,$$

$$F_c=\sum _{h=1}^H\Vert ({\mathbf{p}}_h{\mathbf{p}}_h^T)\oslash \mathbf{XtX}{\Vert }_F^2,$$

$$\mathbf{XtX}={\mathbf{X}}^T\mathbf{X}\oslash \Vert diag({\mathbf{X}}^T\mathbf{X}){\Vert }_2$$

$$\mathbf{XXt}=\mathbf{X}{\mathbf{X}}^T\oslash \Vert diag(\mathbf{X}{\mathbf{X}}^T){\Vert }_2$$

$$\{{\hat{\mathbf{P}}}^X,{\hat{\mathbf{S}}}^X,{\hat{\mathbf{U}}}^X\}=\underset{\mathbf{P},\mathbf{S},\mathbf{U}}{\mathrm{argmin}}\Vert \text{X}-\text{U}\mathbf{S}{\mathbf{P}}^T{\Vert }_F^2+{\lambda }_0{F}_0+{\lambda }_c{F}_c+{\lambda }_r{F}_r,$$

$$F_0=\sum _{h=1}^H(\Vert {\mathbf{p}}_h{\Vert }_2^2-1{)}^2+(\Vert {\mathbf{u}}_h{\Vert }_2^2-1{)}^2$$

$$\mathbf{X}=\left[\begin{array}[]{cc}\mathbf{X}_{1}&\mathbf{0}\\ \mathbf{0}&\mathbf{X}_{2}\\ \mathbf{0}&2\mathbf{X}_{3}\\ \end{array}\right]+0.15\cdot\mathbb{N}(\mathbf{0},\mathbf{I}),$$

$${\mathbf{X}}_{nc}=1{\hat{\mathbf{p}}}_0^T+\hat{\mathbf{T}}{\hat{\mathbf{P}}}^T+\mathbf{E}s.t.\hat{\mathbf{T}}\ge 0,\hat{\mathbf{P}}\ge 0$$

$$\frac{\partial F}{\partial U}=-2(\text{X}-\text{U}\mathbf{S}{\mathbf{P}}^T)\text{P}\mathbf{S}+{\lambda }_0\frac{\partial {\mathbf{F}}_0}{\partial U}+{\lambda }_r\frac{\partial {\mathbf{F}}_{\mathbf{r}}}{\partial U}$$

$$\frac{\partial F}{\partial S}=-2{\text{U}}^T(\text{X}-\text{U}\mathbf{S}{\mathbf{P}}^T)\text{P}$$

$$\frac{\partial F}{\partial P}=-2(\text{X}-\text{U}\mathbf{S}{\mathbf{P}}^T{)}^T\text{U}\mathbf{S}+{\lambda }_0\frac{\partial {\mathbf{F}}_0}{\partial P}+{\lambda }_c\frac{\partial {\mathbf{F}}_{\mathbf{c}}}{\partial P}$$

$$\frac{\partial F_0}{\partial u_{hi}}=4u_{hi}({\mathbf{u}}_h^T{\mathbf{u}}_h-1)$$

$$\frac{\partial F_r}{\partial u_{hi}}=4u_{hi}{\mathbf{n}}_{hi}^T{\mathbf{n}}_{hi},$$

$$\frac{\partial F_0}{\partial p_{hj}}=4p_{hj}({\mathbf{p}}_h^T{\mathbf{p}}_h-1)$$

$$\frac{\partial F_c}{\partial p_{hj}}=4p_{hj}{\mathbf{m}}_{hj}^T{\mathbf{m}}_{hj},$$

Full text

1 Introduction

Matrix factorization methods, which decompose a matrix into a product of factors, are extensively employed to understand complex data. Factors are often useful for highlighting and interpreting special observations (outliers), clusters of similar observations, groups of related variables, and crossed relationships between observations and variables.

Principal component analysis (PCA) [1, 2] is a key tool for that purpose. PCA factorizes a matrix using the criterion of maximizing variance. The PCA model follows the expression:

[TABLE]

where $\mathbf{X}$ is a $N\times M$ data matrix, $\mathbf{\hat{T}}(N\times K)$ is the score matrix containing the projection of the objects onto the principal components (PCs) sub-space, $\mathbf{\hat{P}}(M\times K)$ is the loading matrix containing the linear combination of the variables represented in each of the PCs, and $\mathbf{E}(N\times M)$ is the matrix of residuals. PCA satisfies:

[TABLE]

where $\|\cdot\|_{F}$ stands for the Frobenius norm. Depending on the area of application, loading vectors are constrained to unit length, in order to leave the data variance in the scores to ease interpretation. One interesting property of PCA is that loading vectors can be computed simultaneously or sequentially with exactly the same parameter estimates. This property is a consequence of the components being orthogonal, and while that leads to nice mathematical properties, it seldom reflects the underlying biological or chemical reality.

When used for interpretation by exploring the PCs, PCA has a major shortcoming: The PCs are linear combinations of all the variables, and often combine unrelated sources of variance [3]. On the other hand, for easier interpretation, it is desirable to find factorizations that correspond to a limited number of original variables. This can be achieved by means of rotation [1] or sparse methods like sparse principal component analysis (SPCA) [4, 5]. Another approach is to constrain loadings to agree with the structure of the correlation matrix, like in the group-wise principal component analysis (GPCA) [6].

While previous methods focus on the mode of the variables (columns) of the data, PCA shows exactly the same limitation for interpreting the mode of the observations: i) score vectors are typically non-sparse, and ii) unrelated observations can provide very similar scores. Extensions that apply the sparsity idea to both modes already exist, like some variants of co-clustering [7] or the penalized matrix decomposition (PMD) [8].

In this paper, we introduce the cross-product penalized component analysis (XCAN). XCAN is a matrix factorization based on a loss function that allows a trade-off between variance maximization and structural preservation, aimed at solving the aforementioned problems. The result is a flexible modeling approach that can be used for data exploration in a large variety of problems. We will demonstrate its use with applications from different disciplines.

The rest of the paper is organized as follows. Section 2 introduces the methods on which XCAN is based. Section 3 presents the XCAN algorithm. Section 4 illustrates XCAN through four case studies, including simulated as well as real data from different application fields. Conclusions and future work are discussed in Section 5.

2 Related methods

The proposed XCAN method is inspired by previous developments, notably (i) the SPCA framework based on the lasso (least absolute shrinkage and selection operator), (ii) extensions of SPCA to constrain both modes of the factorization, like co-clustering or the PMD, and (iii) GPCA. From SPCA, we inherit the approach of defining a set of meta-parameters to define the loss function as a trade-off. This trade-off is between captured variance and structural penalties, following the GPCA strategy. The new loss function is intended to reflect both the structure among observations and among variables, in a similar way as in co-clustering or PMD.

2.1 Sparse PCA

The SPCA idea is grounded on the work of several authors more than two decades ago, as described in [9]. There are various versions of SPCA, most based on the modification of the PCA loss in eq. (2) by including sparsity-inducing constraints or penalties with the $L_{0}$ or $L_{1}$ norms [10]. The $L_{0}$-norm of a vector refers to the number of non-zero elements in the vector, and the $L_{1}$-norm of a vector computes the sum of the absolute values of the vector entries. The application of the $L_{1}$-norm in a regression setting was originally called the least absolute shrinkage and selection operator (lasso) [11]. In this section, we focus on the lasso versions of SPCA for its widespread use in model interpretability [12].

The SCoTLASS algorithm [4] incorporates the lasso in the PCA calibration as follows:

[TABLE]

where $\mathbf{\hat{p}}^{SL}$ is the resulting sparse loading, and $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$ refer to the $L_{1}$ and $L_{2}$ norms, respectively. To obtain successive components, the SCoTLASS optimization constrains the second and further sparse loadings to be orthogonal to the rest.

The SCoTLASS criterion is very computational demanding [13] and has numerical limitations. For this reason, Zou et al. [5] introduced an alternative formulation to generate sparse components based on regularized regression, using a criterion close to the (naive555The regular elastic net makes use of a scaling factor between the lasso and the ridge penalty.) elastic net [14], which is a combination of both the lasso and the ridge ($L_{2}$-norm) penalties:

[TABLE]

where we distinguish between sparse loadings $\mathbf{\hat{P}}^{SP}$ (sometimes also referred to as weights in similar modeling frameworks) and orthonormal loadings $\mathbf{\hat{Q}}^{SP}$, $\mathbf{p}_{h}$ represents the $h$th column vector in $\mathbf{P}$, with $H$ the number of components. The solution proposed for eq. (4) is a biconvex optimization where sparse weights and orthogonal loadings are obtained using an alternating approach. The algorithm is simultaneous, in the sense that all components are computed in the same alternating iteration. An alternative sequential variant is defined in [15]. A particular solution of SPCA in eq. (4) is when $\lambda_{2}=\infty$, which is the most popular choice when the number of columns in the data is much higher than the number of rows. Then, the sparse loadings can be computed by soft-thresholding, simplifying and improving the efficiency of the computation: $g_{\lambda}(p)=\operatorname{sign}(p)\left(|p|-\lambda\right)_{+}$.

2.2 Extensions of sparsity to both modes

Witten et al. [8] propose a new sparse algorithm referred to as the Penalized Matrix Decomposition (PMD). It can be used to constrain both the number of observations and variables contributing to each factor using soft-thresholding. The PMD follows:

[TABLE]

The corresponding pseudo-singular value is then obtained as:

[TABLE]

After each component is obtained, projection deflation is performed as:

[TABLE]

The authors show that this solution, when only applying sparsity to the loadings, is connected with SCoTLASS and SPCA [8].

A similar approach was introduced in [7] with the goal of performing coclustering:

[TABLE]

where $\mathbf{t}_{h}$ represents the $h$th column vector in $\mathbf{T}$, and the same penalty variable is used in both modes. This loss is used within an alternating optimization, which, unlike PDM, produces all components in a single run.

2.3 Group-wise PCA

GPCA limits each component to represent the variability of a single group of variables, which in turn allows to interpret each component independently (provided the variables are not overlapping). GPCA starts with the identification of a set of $K$ (possibly overlapping) groups of correlated variables. This is achieved by applying thresholding in a pseudo-correlation matrix $\mathbf{M}(M\times M)$ of the data. In the original formulation of GPCA, the MEDA approach (Missing-data for Exploratory Data analysis) [16] was implemented to obtain $\mathbf{M}$.

Besides, GPCA may be simpler to use in practice than sparse methods based on the lasso, because by inspecting $\mathbf{M}$ we can often identify suitable values for the threshold or even when the GPCA model is not appropriate at all. In comparison, the main challenge when using sparse methods is to find suitable values for meta-parameters like $c$ (3), $\lambda_{1}$ and $\lambda_{2}$ (4), $c_{1}$ and $c_{2}$ (5) or $\lambda$ (8). This advantage, however, comes at a price. While the capability to reflect the structure in the map $\mathbf{M}$ is an appealing property, GPCA is an “all or nothing” approach, meaning that a variable is either in a group or not, and different GPCA models can be obtained for very similar values of the threshold.

3 Cross-product Penalized Component Analysis (XCAN)

With XCAN, we would like to make the most of the advantages of sparse methods in one or the two modes and combine that with the idea behind GPCA. Like PMD, XCAN factorizes a matrix into three matrices:

[TABLE]

The loss function for the XCAN factorization is as follows:

[TABLE]

where the first part is the actual model, with $\mathbf{S}$ constrained to be diagonal, and $F_{c}$ and $F_{r}$ are defined to constrain the structure of the model. The meta-parameters $\lambda_{c}$ and $\lambda_{r}$ control the level of these penalties. We define the penalties as follows:

[TABLE]

[TABLE]

where $\mathbf{u}_{h}$ and $\mathbf{p}_{h}$ are the $h$th column vectors in $\mathbf{U}$ and $\mathbf{P}$, respectively, $H$ is the total number of XCAN components (XCs), $\oslash$ is the Hadamard (element-wise) division. $\mathbf{XtX}(M\times M)$ and $\mathbf{XXt}(N\times N)$ denote maps of relationship between variables and observations, respectively, and are given as inputs. To avoid numerical problems in the divisions, values below a threshold in $\mathbf{XtX}$ and $\mathbf{XXt}$ are set to that threshold. In this paper, we fixed the value of this threshold to 0.01 in all experiments.

3.1 Cross-product matrices and rationale

We will generally refer to $\mathbf{XtX}$ and $\mathbf{XXt}$ as cross-product matrices, because several of their possible definitions can be computed as cross-products, and we expect them to be symmetric. This is also the reason for the name XCAN, where “X” stands for “cross”. For instance, $\mathbf{XtX}$ may be set to the correlation matrix of $\mathbf{X}$ and $\mathbf{XXt}$ to the correlation matrix of $\mathbf{X}^{T}$, respectively. $\mathbf{XtX}$ can also be set to any map $\mathbf{M}$ used in GPCA.

The rationale behind the definition of the loss in (10)-(11) is as follows. $\mathbf{XtX}$ and $\mathbf{XXt}$ contain the correlation structure in the variables and observations, respectively. Values close to 0 in those matrices identify unrelated variables or observations. Using an element-wise division, we prevent unrelated elements to be part of the same component. That way, we obtain sparse components that agree with the structure enforced by the input cross-product matrices.

The choice of matrices $\mathbf{XtX}$ and $\mathbf{XXt}$ is principal in XCAN, and should be carefully done taking into account the goal of the analysis, in a similar way to when selecting the data pre-processing. For instance, it is customary to mean-center data for some types of analyses, but not for others (e.g., for spectra). Although not compulsory, using matrices $\mathbf{XtX}$ and $\mathbf{XXt}$ that are consistent with the data pre-processing is expected to provide more coherent results. By taking that into account, we use the following definition of cross-product in our experiments in this paper, inspired by Pearson’s correlation:

[TABLE]

[TABLE]

The advantage of these definitions is that they do not require data to be column-wise or row-wise mean-centered, like correlation matrices do.

Cross-product matrices can be conveniently post-processed to modify the behavior of XCAN, which adds flexibility to the modeling approach. For instance, a possible post-processing operation is thresholding. As discussed before, GPCA was defined with the suitable property that its meta-parameter (i.e., the threshold) can be set upon visual inspection of the map $\mathbf{M}$. Similarly, we can inherit this idea in XCAN by thresholding $\mathbf{XtX}$ and $\mathbf{XXt}$, in order to discard minor correlations from the analysis. This results in useful means to further impose sparsity, as in GPCA. We can also use thresholding in such a way that only positive correlations are kept in $\mathbf{XtX}$ and $\mathbf{XXt}$. This can be useful to derive loadings and scores vectors in XCAN where all non-zero elements have the same sign, something in line with non-negativity constraints. In the examples, we will use this thresholding idea or non-negativity constraints, when suitable.

Since cross-product matrices can be considered as a way to include structural penalties in both modes of the data in XCAN, we can also use them in different ways, e.g., as in chemometrics literature, to impose smoothness, to connect different data sets of same individuals or variables (data fusion) or to include apriori information into a model. Studying those applications in detail is out of the scope of this paper. We will, however, show an application of XCAN by incorporating the class labels of the samples into the model.

3.2 Algorithmic Approach

The XCAN model is fitted to the data by solving for all components simultaneously using gradient-based all-at-once optimization. To constrain the vectors in $\mathbf{P}^{X}$ and $\mathbf{U}^{X}$ to unit length, like in the SVD, a suitable way is to redefine the loss as:

[TABLE]

where

[TABLE]

If $\lambda_{0}$ is set to a sufficiently large value, this additional term in the loss will serve the purpose of normalizing factors $\mathbf{a}_{h}$ and $\mathbf{b}_{h}$. In our experiments in Section 4, we set $\lambda_{0}=1$.

We solve eq. (15) by computing the partial derivatives of the loss function with respect to $\mathbf{U},\mathbf{S}$ and $\mathbf{P}$ (as given in the Appendix), constructing the gradient, and then using a gradient-based optimization algorithm. In our experiments, we use the Poblano Toolbox [17], that has several unconstrained gradient-based optimization algorithms such as the nonlinear conjugate gradient (NCG), and limited-memory BFGS (L-BFGS) [18]. When non-negativity constraints are desired, we can also use the limited-memory BFGS with bound constraints (LBFGS-B)666We use the implementation at https://github.com/stephenbeckr/L-BFGS-B-C..

It should be noted that the XCAN optimization can have many local minima. In the experiments, we initialized the algorithm using the PCA solution, but ideally, multiple random starts would help with the local minima problem.

4 Experiments

4.1 Simulation

We start with a simulated experiment to demonstrate the properties and flexibility of XCAN. We simulate three data sets $\mathbf{X}_{1}$, $\mathbf{X}_{2}$ and $\mathbf{X}_{3}$, each with 5 observations and 5 variables. Each data set is generated with high correlation between the variables, using the simuleMV tool [19]. We, then, construct $\mathbf{X}$ as follows:

[TABLE]

where $\mathbb{N}(\mathbf{0},\mathbf{I})$ denotes noise randomly drawn from the normal distribution with 0 mean and unit variance. Finally, $\mathbf{XtX}$ follows eq. (13) and $\mathbf{XXt}$ eq. (14).

Matrices $\mathbf{XtX}$ and $\mathbf{XXt}$ for one of the simulations are shown in Figure 1. We can see that variables (Figure 1(a)) are approximately grouped in two groups of 5 and observations (Figure 1(b)) describe three major groups.

Figure 2 shows the result of applying XCAN with three components and for different values of the meta-parameters $\lambda_{c}$ and $\lambda_{r}$. The figure shows an upper bar plot with the scores, computed as $\mathbf{\hat{T}}^{X}=\mathbf{\hat{U}}^{X}\mathbf{\hat{S}}^{X}$, and a lower bar plot with the loadings $\mathbf{\hat{P}}^{X}$. Figure 2(a) shows regular PCA, since the structural penalties are deactivated, which is used as a baseline. In such setting, each loading/score vector contains information about all variables/observations, respectively. Figures 2(b), 2(c) and 2(d) show the application of XCAN with structural penalties in the loadings, scores, and both loadings and scores, respectively. In all cases, the XCAN model works as expected, and the variance remains reasonably close to the PCA model. This shows that the penalties meet the true data structure at a minor price in terms of explained variance.

In order to see the performance of XCAN with different number of components, we have also compared the models with one, two, and three XCs in the example, with the structural penalties activated in loadings and scores. Results can be inspected by comparing Figures 3 and 2(d). We observe that the XCAN model of the simulated data is very stable.

As discussed before, we can apply thresholding in XCAN. If we set all values in $\mathbf{XtX}$ and $\mathbf{XXt}$ with magnitude less than 0.5, i.e., within the interval $(-0.5,0.5)$, to zero, we obtain the matrices in Figure 4, where most of the spurious correlations are discarded. Figure 5 compares 3-component XCAN models with and without thresholding in the cross-product matrices. We can see that thresholding can be effective in terms of imposing sparsity. In the current example, we can achieve a sparser model using lower values of the meta-parameters and capture higher variance.

4.2 Animal Data

This case study makes use of a toy data set created to illustrate the co-clustering algorithm in [7]. The data contains 34 observations, most of them representing animals, each with 17 attributes, including binary and continuous attributes. The data values are non-negative. Binary attributes are used to describe if an animal has (1) or does not have (0) a specific feature. For instance, in the feature ‘Carnivore’, ‘Lion’ contains a 1 and ‘Cow’ a 0.

Scores and loadings plots for the first two PCs of PCA are shown in Figure 6, for auto-scaled data, that is, mean-centered across the animals mode and scaled to unit variance within the attributes. The first PC is dominated by birds, e.g., ‘Eagle’, ‘Blackbird’ or ‘Chicken’, which share features like ‘Wings’, ‘Feathers’ or ‘Has a beak’. Birds are small, and for this reason they are located opposite to big animals, in particular, to the ‘House’, which was included in the data set as an outlier. The second PC contrasts ‘Dangerous’ and ‘Extinct’ animals, most notably the ‘T. Rex’, with those ‘Domesticized’ and ‘Eaten by Caucasians’, which are also correlated with ‘Breathe under water’. Both components are interpretable, but still complex in the sense that they mix different concepts: birds + small vs big, dangerous vs domesticated + fish. As such, this data set is a perfect example of the limitations of PCA reported in the introduction.

To perform the XCAN analysis, we auto-scaled the data. Since the goal is to understand differences among individuals and group them as in (co-)clustering, it makes sense to center the data so that we study the variability among the set of individuals, instead of with respect to some center of coordinates. Scaling seems to equalize the influence of the variables in the model. Since we want each component to be a cluster of similar (and not antagonist) individuals, we constrain $\mathbf{XXt}$ to be non-negative. In particular, every entry in $\mathbf{XXt}$ with a value below 0.5 is hard thresholded to 0, i.e., values within the interval $(-\infty,0.5)$ are thresholded. That way, components will be extra sparse and yield only positive (or only negative) scores. In case a component contains negative scores we simply change its sign. For this analysis, we are not concerned with the signs of the loadings, but we still want loadings to be extra sparse, so we threshold the entries of $\mathbf{XtX}$ using the threshold value of 0.5, i.e., setting every entry with a value within the interval $(-0.5,0.5)$ to 0. This allows negative correlations of -0.5 or lower. Resulting cross-product matrices are shown in Figure 7.

The 9-component XCAN model is shown in Figure 8. We can see that all XCs are sparse in both scores and loadings, and all scores are non-negative, as expected. The model extracts meaningful components, that can be interpreted as co-clusters. The first component represents wild birds, with their corresponding features. Interestingly, ‘Chicken’ is not included in the component. The reason for that can be seen in Figure 9, which represents a zoom of matrix $\mathbf{XXt}$: ‘Chicken’ and ‘Eagle’, even though they share the features in the first component, they are completely different animals according to matrix $\mathbf{XXt}$, showing a correlation close to 0. Thus, XCAN does not place them in the same component. Since ‘Eagle’ is strongly correlated with the other birds, the model selects the former among the group of scores of the first component. The second component represents domesticated animals consumed for food. The third component focuses on the feature of extinct, but only two (‘T.Rex’ and ‘Neanderthal’) out of the five extinct animals (including ‘Mamouth’, ‘Sabre Tiger’ and ‘Triceratops’) show relevant scores. The reason for that is similar to before and also apparent in Figure 9. ’Neanderthal’ and ’Triceratops’ are uncorrelated, and for this reason they cannot share the same component. The rest of the components are more or less self-explanatory.

The combination of sparsity and grouping in components 1 and 3 is coherent with GPCA, but not with PCA/SPCA models, where one would expect that all individuals with feathers and wings, or those extinct, should be placed in the same component. If, for some reason, we are interested only in sparsity but not on the grouping characteristic of XCAN, then traditional sparse methods should be used for the analysis. With XCAN, we explicitly avoid, within a component, objects or features that are not associated in the cross-product matrices.

In conclusion, we observe that the XCAN model performs exactly as expected: components are sparse and reflect accurately the structure in the cross-product matrices.

4.3 Real Data

4.3.1 Vast Challenge for Cybersecurity

The increase in the number of cybersecurity incidents, coupled with the shortage of specialized professionals, has created a need for efficient data analysis tools to support the detection, triaging and analysis of incidents [20]. In particular, anomaly-based Intrusion Detection Systems (IDS) [21] are fundamental resources to unveil new attack strategies. A large number of intrusion detection approaches based on PCA have been proposed in the last two decades [22, 23, 24, 25, 26]. In a recent paper [27], the multivariate detection approach was extended to GPCA.

With GPCA, we can identify anomalies in the data following a straightforward approach. components simplify the interpretation by reducing the number of variables to examine. Components can be interpreted one at a time, following a very similar workflow to the one security analysts use in traditional software tools. Since XCAN inherits the features of GPCA, we will explore its performance in the cybersecurity domain. In comparison to GPCA, XCAN can also impose sparsity in the rows, which in the cybersecurity domain typically corresponds to time-resolved data. This is useful to speed up the analysis of an incident, so that the analyst can focus on a few points in time to troubleshoot the problem.

The data of the present case study comes from the VAST 2012 2nd mini-challenge [28] and was captured in a corporate network during a time frame of two days. During that time period, a botnet compromised the network, causing performance problems and the emergence of spyware. The raw data is parsed in a total of 265 features at 1 minute intervals, yielding a 2345 x 265 data matrix. More details can be found in [29]. Cross-product matrices following eq. (13) and eq. (14) are shown in Figure 10. The plots illustrate that applying a group-wise constrained model is reasonable, both in observations and features, since the cross-product matrices contain squares of high correlation.

We used the following strategy for the application of XCAN in this problem domain. Cross-product matrices were thresholded to positive values above 0.7 ($\mathbf{XtX}$) and 0.3 ($\mathbf{XXt}$). Data was auto-scaled for similar reasons as in the previous case study. The mean centering is useful to detect anomalies as deviations from the mean. The scaling normalizes the relevance of the variables in the model.

Figure 11 compares the first four components by GPCA (left column) and the 4-component model by XCAN with sparse loadings (middle column) and both sparse scores and loadings (right column). We re-ordered the components of XCAN to match those in GPCA. To improve the visualization for security analysts, we included statistical control limits in the scores of each component, so that anomalies can be easily identified. Both GPCA and XCAN provide useful results for anomaly detection, and components can be easily interpreted one-at-a-time. For instance, the first component of GPCA (upper-left figure) describes a number of anomalies around sampling time 400, which are related to the variables in the loadings. So a brief description of the location in time of the anomalies, and their diagnosis (the variables related to the abnormal behavior) is obtained in a single plot. In comparison to PCA (e.g., see [26]), this approach for anomaly detection is much simpler and easy to understand for security analysts.

We can also see that GPCA and XCAN with sparse loadings are very similar. If we apply also sparsity in the rows in XCAN, a reduced set of observations are identified for each component, allowing the analyst to focus on a subset of time points to proceed with a more detailed forensic analysis.

4.3.2 Olive Data

Oil samples from 24 brands and several types (olive, corn, sesame, etc.) were obtained. Infrared spectra were measured using a Nicolet 5-DX FT-IR system. Each spectrum consisted of 1556 measurements from 3600 to 600 cm-1 of which two regions were used in the analysis here in accordance with the original publication [30]. The resulting spectra, with dimension $80\times 188$, are shown in Figure 12, with different colors representing different brands of oil.

One meaningful way to model spectra with matrix factorization is:

[TABLE]

where $\mathbf{\hat{T}}$ and $\mathbf{\hat{P}}$ are non-negative, and $\mathbf{\hat{p}}_{0}$ takes the roll of a baseline. For this example, we will follow the same approach with XCAN, and also incorporate the baseline and non-negativity constraints on $\mathbf{\hat{U}}^{X}$, $\mathbf{\hat{S}}^{X}$ and $\mathbf{\hat{T}}^{X}$. The result with $\lambda_{c}=\lambda_{r}=0$ is shown in Figure 13. The top shows the baseline, contained in $\mathbf{\hat{p}}_{0}$. We observe that components are generally not sparse.

Since in this case study the data is not centered and it is far from the origin of coordinates, the cross-product expressions in eqs. (13) and (14) are impractical: they are almost completely filled with 1s. For this reason, we post-processed cross-product matrices by subtracting a baseline computed with the minimum values of each column (row) for $\mathbf{XtX}$ ($\mathbf{XXt}$). This approach gives us the cross-product matrices in Figure 14(a) and (b). We can see that the cross-product matrix $\mathbf{XXt}$ still does not indicate any groups of observations. For this reason, we will not use this cross-product matrix with XCAN. Instead, here, we will perform a different analysis to illustrate another potential use of XCAN with data from samples, which has class information, and define $\mathbf{XXt}$ based on the classes of the samples as in Figure 14(c). Each matrix entry contains 1 if the corresponding observations belong to the same class, or 0 otherwise. This illustrates how we can use external information to influence the XCAN outcome, instead of using the cross-product matrix definitions given earlier.

The 4-component XCAN model constrained only with $\mathbf{XtX}$, along with the baseline, is shown in Figure 15. The model is similar to the one maximizing variance but with sparse loadings. Figure 16 shows the same analysis using the observation class map within XCAN. Now each component shows the particularities of a different oil brand. Notice this works as a description of the oil brands, not a discrimination. For instance, the first component shows the pattern of the yellow class, without an attempt to distinguish it from the other brands. It might not be a complete description, though. For instance, look at the fourth component. It is focused on the blue class of oil, which is different from the rest in the last interval of wavelengths. However, XCAN only shows one wavelength in this component, because this wavelength is uncorrelated to the others where the blue brand shows its peculiarities. This is interesting information that can improve interpretation but cannot be found with similar sparse methods.

5 Conclusion

In this paper, we introduce the Cross-product penalized component analysis (XCAN) method and illustrate its use with examples in different application domains. XCAN combines variance maximization and structural penalties, that are specified in the form of cross-product matrices and result in sparse matrix factorizations. This provides a flexible modeling framework to explore complex data and enhance the structural information. We plan to extend the application of XCAN to a variety of problems, in particular to data fusion and the derivation of gray models with a-priori information.

Appendix: Partial Derivatives

Let us call $F$ the loss function in eq. (15). The partial derivatives are as follows:

[TABLE]

[TABLE]

[TABLE]

where all partial derivatives with respect to $\mathbf{U}$, $\mathbf{D}$ and $\mathbf{P}$ are of the same dimension of the corresponding matrices. The partial derivatives of the terms in the loss functions $F_{0}$, $F_{r}$ and $F_{c}$ can be calculated element-wise, using the symmetric structure of $\mathbf{XtX}$ and $\mathbf{XXt}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $\mathbf{n}_{hi}=\mathbf{u}_{h}\oslash\mathbf{x}\mathbf{x}t_{i}$, being $\mathbf{x}\mathbf{x}t_{i}$ the $i$-th row of $\mathbf{XXt}$, and $\mathbf{m}_{hj}=\mathbf{p}_{h}\oslash\mathbf{x}t\mathbf{x}_{j}$, being $\mathbf{x}t\mathbf{x}_{j}$ the $j$-th row of $\mathbf{XtX}$.

Acknowledgement

This work is partly supported by the Spanish Ministry of Economy and Competitiveness and FEDER funds through project TIN2017-83494-R and the “Plan Propio de la Universidad de Granada”.