EMBL Centre for Network Analysis workshop 9-11 December 2015



Data integration and mining with graphs:


Tensor factorizations for dynamic graphs



Jean-Karim Hériché

Ellenberg lab

Cell Biology and Biophysics unit










 As gene networks become available for multiple time points and experimental conditions, there's a need for principled ways of detecting patterns and their evolution in such dynamic graphs. We will introduce a tensor representation of dynamic graphs and tensor factorization methods as ways of investigating patterns in these graphs.




 Many biological processes are driven by systems that evolve over time. New technologies make the dynamics of such systems more accessible with increasingly higher time resolution. In particular, data on gene networks can be acquired for multiple time points and experimental conditions. Detecting patterns and how they evolve over time and/or across experimental conditions in such multidimensional data is becoming a common task. For example, one may have performed AP-MS experiments of several proteins at different cell cycle or developmental stages and be interested in identifying protein complexes and how they change between time points.


Dynamic graphs

 - What is a dynamic graph ?

        A dynamic graph can be seen as multiple instances of a graph whose edges change between two instances. This definition includes time-varying graphs whose edges change as a function of time but also multigraphs (graphs in which nodes are connected with edges of different types) if you segregate each edge type into a separate instance of the graph. Examples of dynamic graphs are protein interaction networks and gene co-expression networks obtained at different cell cycle or developmental stages or under different experimental conditions (e.g. control vs treatment vs disease). With this definition, a dynamic graph can therefore be represented as a series of adjacency matrices.


 - Different approaches to mining dynamic graphs

        It is common to treat dynamic graphs in one of three ways:


All the above treat time or experimental conditions separately from the nodes. A global approach treating time/conditions and nodes simultaneously would better capture evolving patterns. One of the first such approach did so by extending the modularity coefficient [1] to multiple instances of a graph [2,3] by defining links between graph instances. Here, we'll introduce another global approach which treats genes, time and conditions as multiple dimensions of the data.



 - What is a tensor ?

  There are different definitions of tensors used in different fields but they all come down to the same concept that tensors are a generalization of scalars, vectors and matrices: a scalar is a rank zero tensor, a vector is a rank one tensor, a matrix a rank two tensor... Therefore for most practical purposes, a tensor is a multidimensional array. It's easy to see that a dynamic graph represented as a series of adjacency matrices is a three-dimensional tensor.


 - Terminology and notation:

 Although the notation and terminology are not fully standardized, many papers now follow the conventions in [4] which we also adopt here. A tensor will be denoted using underlined bold face e.g. X to distinguish it from the matrix X.

 The order of a tensor is its number of dimensions also called modes (or sometimes ways). The equivalent of matrix rows and columns are fibers and when used as vectors, they are assumed to be column vectors. Slices are 2D sections of a tensor. See Figure 1 for a graphical representation of these terms for third-order tensors.

Figure 1: Tensor terminology

 - Working with tensors:

 Most of the tasks we're interested in rely on transforming tensors into matrices using a process called matricization or unfolding. Matricization consists in rearranging the elements of the tensor into a matrix. The most commonly used way of doing this is called mode-n matricization which takes the mode-n fibers of the tensors and arranges them as columns of a matrix (Figure 2). The mode-n unfolding of Z is denoted Z(n). The order of columns is not important as long as it is consistent across calculations. Sometimes it can be useful to represent a tensor as a vector, in a process called vectorization. Ordering of the elements in the vector is also not important as long as it is consistent.

 Although tensors can be multiplied together, we will only need the mode-n multiplication of a tensor with a matrix in which the mode-n fibers of the tensor are multiplied by the matrix. This is written Z = Y xn V and can be expressed in matricized form: Z(n) = VY(n). See Figure 3 for an intuitive graphical example.

        Working with matricized tensors also makes use of various matrix products: the Kronecker product, the Khatri-Rao product and the Hadamard products. See paragraphs 1 and 2 of [5] and chapter 1.4 of [6] for a more detailed introduction to tensor algebra.

Figure 2: Matricization of a tensor
Figure 3: Mode-n multiplication with a matrix

Graph nodes clustering by matrix factorization

        Nodes of static graphs can be clustered by applying clustering methods to the adjacency matrix. Among these methods, matrix factorization approaches have already been successfully applied to some computational biology questions (see [7] for a review on non-negative matrix factorization (NMF) in computational biology, [8,9] for examples of spectral clustering). Spectral clustering and NMF solve a problem of the form: A = W * H + E. A more general form of factorization is the singular value decomposition (SVD): A = U * S *  VT where columns of U and VT are constrained to be orthogonal. SVD corresponds to the decomposition of the matrix into the sum of outer products of vectors (Figure 4)

Figure 4: Singular value decomposition

Tensor factorization

 Tensor factorization can be seen as a natural extension of SVD to an arbitrary number of dimensions (Figure 5).

Figure 5: A higher order SVD
The general form is the Tucker decomposition which factorizes a tensor into a small core tensor and a set of matrices (See Figure 6). See [5] for a general introduction and [6] for non-negative factorizations.

        The Tucker decomposition can be written in terms of matrix unfolding of the tensor X (of size I x J x K):




where ⊗ is the Kronecker product and A, B and C are factor matrices of size I x R, J x S and K x T respectively where R,S and T are the numbers of factors desired for each dimension. The factor matrices are usually constrained to have orthogonal columns.

Like PCA and other matrix decompositions, the Tucker decomposition is not unique as the core can be modified without affecting the quality of the solution as long as the inverse modification is applied to the matrices. Restricting the core tensor to a diagonal gives the CP (or PARAFAC) decomposition (Figure 7, which is equivalent to Figure 5) in which the same R factors are shared across the matrices. The CP decomposition can also be expressed in terms of matrices:

X(1)A(C B)T  

X(2)B(C A)T

X(3)C(B A)T

where ⊙ is the Khatri-Rao product. A form of non-negative tensor factorisation is obtained by constraining A, B and C to have non-negative entries.

Figure 6: Tucker decomposition of a third order tensor
Figure 7: The CP decomposition

Non-negative tensor factorisation of dynamic graphs

 This can be thought of as an extension of spectral clustering to more than one adjacency matrix. Our dynamic graph is represented by a third order tensor G of genes x genes x time points (or experimental conditions) in which each frontal slice is the genes x genes adjacency matrix of a weighted undirected graph. The edge weights are also assumed to be non-negative (≥ 0). Using a non-negative CP decomposition, we can factorize G into non-negative factor matrices A, B and C. Because the graph is undirected, the adjacency matrices are symmetric which leads to A = B. A is a genes x factors matrix and C is a factors x time points matrix. Therefore, if we think of the factors as representing clusters, elements of A can be interpreted as gene cluster memberships and elements of C as measuring the temporal activity of the clusters.



- MATLAB toolboxes:  tensor, nway (compatible with Octave)

- R packages*: PTak, ThreeWay, rTensor

- perl:  Algorithms::Cube (https://github.com/jkh1/Perlstuff)

- python*: scikit-tensor (https://github.com/mnick/scikit-tensor)


* non-negativity constraint is not available in R and python packages (December 2015).


Example in R

## Load required package



## Where the data is

## This is a simulated dynamic graph with 30 nodes and 6 time points

dataDir <- "~/Documents/Bio-IT/courses/Data_integration_and_mining_with_graphs/dynamic_graph"


## Get list of adjacency matrices

## Each matrix is in a csv file

fileList <- dir(path=dataDir,pattern = ".csv")


## Read all adjacency matrices into an array

A <- array(as.numeric(NA),dim=c(30,30,6)) # There are 6 matrices of size 30 x 30

for (i in 1:length(fileList)){

    A[,,i] <- as.matrix(read.delim(file.path(dataDir,fileList[i]), sep = ';', header=TRUE, row.names=1))



## Convert list of adjacency matrices to a tensor

G <- as.tensor(A)


## Perform CP factorization

cpG <- cp(G,num_components = 4)

## View the factor matrices

## Note that the first two are identical (up to numerical accuracy)

## due to the fact that the graph adjacency matrices are symmetric





Example in perl


use strict;

use warnings;

use Algorithms::Matrix;

use Algorithms::Cube;


# Where the data is

my $dataDir = "$ENV{'HOME'}/Documents/Bio-IT/courses/Data_integration_and_mining_with_graphs/dynamic_graph";


# Number of time points

my $T = 6;


# Read dynamic graph from files as series of adjacency matrices

my @G;

foreach my $t(0..$T-1) {

  my ($rowLabels,$colLabels);

  ($G[$t], $rowLabels, $colLabels) =



# Create Cube (i.e. third order tensor) object from matrices

my $G = Algorithms::Cube->new(@G);


# CP factorization with non-negativity constraint

my $r = 4; # Number of components

my ($A,$B,$C,@norms) = $G->nncp($r,symmetry=>1,norms=>1);




1. Newman MEJ, Girvan M (2004) Finding and evaluating community structure in networks. Phys. Rev. E 69, 026113.

2. Mucha PJ, Richardson T, Macon K, Porter MA, Onnela JP (2010) Community structure in timedependent, multiscale, and multiplex networks. Science 328: 876–878.

3. Bassett DS, Porter MA, Wymbs NF, Grafton ST, Carlson JM, et al. (2013) Robust detection of dynamic community structure in networks. Chaos 23: 013142.

4. Kiers HAL (2000) http://www.wiley.com/legacy/wileychi/chemometrics/582_ftp.pdf. J. Chemometrics 14: 105-122.

5. Kolda TG, Bader BW (2009) Tensor decompositions and applications. SIAM Rev., 51(3), 455–500.

6. Cichocki A, Zdunek R, Phan AH, Amari SI (2009) Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. John Wiley & Sons.

7. Devarajan K (2008) Nonnegative matrix factorization: an analytical and interpretive tool in computational biology. PLoS Comput Biol 4(7):e1000029.

8. Kluger Y, Basri R, Chang JT, et al. (2003) Spectral biclustering of microarray data: coclustering genes and conditions. Genome Res. 13(4):703–16.

9. Imangaliyev S, Keijser B, Crielaard W, Tsivtsivadze E (2015) Personalized microbial network inference via co-regularized spectral clustering. Methods 83:28-35.