Towards an Efficient Use of the BLAS Library for Multilinear Tensor Contractions

The main motivations behind the Basic Linear Algebra Subroutines were modularity, efficiency, and portability through standardization [10]. With the advent of architectures with a hierarchical memory, BLAS was extended with its levels 2 and 3, raising the granularity level to allow an efficient use of the memory hierarchy—by amortizing the costly data movement with the calculation of a larger number of floating point operations. With respect to the full potential of a processor, the efficiency of BLAS 1, 2, and 3 routines is roughly 5%, 20% and 90+%, respectively. It is widely accepted that for large enough operands, these values are nearly perfect. Moreover, the scalability of BLAS 1 and 2 kernels is rather limited, while that of BLAS 3 is typically close to perfect. As demonstrated by its universal usage, BLAS succeeded in providing a portability layer between the computing architecture and both numerical libraries and simulation codes.

In addition to the reference library,⁶⁶6Available at http://www.netlib.org/blas/. nowadays many implementations of BLAS exist, including hand-tuned [12, 11], automatically-tuned [23, 24], and versions developed by processors manufacturers [13, 15, 14]. Despite their excellent efficiency, such libraries might not be well suited for tensor computations. For instance, in the following it will become apparent that the design of the BLAS interface falls short when matrices are originated by slicing higher dimensional objects; specifically, costly transpositions are caused by the constraint that matrices have to be stored so that one of the dimensions is contiguous in memory. Recent efforts aim at lifting some of these shortcomings [25].

A list of BLAS related definitions follows.

• •

  Stride: distance in memory between two adjacent elements of a vector or a matrix. BLAS requires matrices to be stored by columns, that is, with vertical $stride=1$.

• •

  Leading dimension: distance between two adjacent row-elements in a matrix stored by columns. In BLAS routines, leading dimensions are specified through arguments.

• •

  GEMM: $C\leftarrow\alpha A^{\bullet}B^{\bullet}+\beta C$; kernel for general matrix-matrix multiplication (BLAS 3). $A,B$ and $C$ are matrices; $\alpha$ and $\beta$ are scalars; $M^{\bullet}$ indicates that matrix $M$ can be used with/without transposition.

• •

  GEMV: $y\leftarrow\alpha A^{\bullet}x+\beta y$; kernel for general matrix-vector multiplication (BLAS 2). $A$ is a matrix, $x$ and $y$ are vectors, and $\alpha$ and $\beta$ are scalars; $A^{\bullet}$ indicates that matrix $A$ can be used with/without transposition.

• •

  GER: $C\leftarrow\alpha xy^{H}+C$; outer product (BLAS 2). $C$ is a matrix, $x$ and $y$ are vectors, and $\alpha$ is a scalar.

• •

  DOT: $\alpha\leftarrow x^{H}y$; inner product (BLAS 1). $x$ and $y$ are vectors, and $\alpha$ is a scalar.

In order to illustrate our methodology with a first well-known example, and thus facilitate the following discussion, here we consider one of the most basic contractions, the multiplication of matrices.

To the reader familiar with the jargon based on Einstein convention, given two 2-dimensional tensors $A^{ih}$ and $B_{hj}$, we want to contract the $h$ index. In linear algebra terms, the matrices $A$ and $B$ have to be multiplied together.

For the sake of this example, we assume that only level 1 and 2 BLAS are available, while level 3 is not; in other words, we assume that there exists no routine that takes (pointers to) $A$ and $B$ as input, and directly returns the matrix $C$ as output. Therefore, we decompose the operation in terms of level 1 and 2 kernels by slicing the operands, exactly the same way we will later decompose a contraction between tensors in terms of matrix-matrix multiplications.

In the following, we provide a number of alternatives for slicing the matrix product. Matrices $A,B$ and $C$ are of size $m\times k$, $k\times n$, and $m\times n$, respectively. Also, $F^{i}$ and $F_{j}$ respectively denote the $i$-th row and $j$-th column of the matrix $F$.

1. 1.

   $A$ is sliced horizontally. The slicing of $i$, the first index of $A$, leads to a decomposition of $C$ as a sequence of independent GEMVs: $C:=\forall_{i}\ a^{i}B$

   $$C:=AB=\left[\begin{array}[]{c}a^{1}\\ \hline\cr\vdots\\ \hline\cr a^{m}\end{array}\right]B=\left[\begin{array}[]{c}a^{1}B\\ \hline\cr\vdots\\ \hline\cr a^{m}B\end{array}\right].$$

2. 2.

   $B$ is sliced vertically. The slicing of $j$, the second index of $B$, also leads to a decomposition of $C$ as a sequence of independent GEMVs: $C:=\forall_{j}\ Ab_{j}$

   $$C:=AB=A[b_{1}|b_{2}|\dots|b_{n}]=[Ab_{1}|Ab_{2}|\dots|Ab_{n}].$$

3. 3.

   $A$ is sliced vertically and $B$ horizontally. The slicing of $h$, the contracted index, leads to expressing $C$ as a sum of GERs: $C:=\sum_{h}a_{h}b^{h}$

   $$C:=AB=[a_{1}|\dots|a_{k}]\left[\begin{array}[]{c}b^{1}\\ \hline\cr\vdots\\ \hline\cr b^{k}\end{array}\right]=a_{1}b^{1}+a_{2}b^{2}+\dots+a_{k}b^{k}.$$

4. 4.

   $A$ is sliced horizontally and $B$ vertically. If both indices $i$ and $j$ are sliced, $C$ is decomposed as a two-dimensional sequence of DOTs: $C:=\forall_{i}\forall_{j}\ a^{i}b_{j}$

   $$C:=AB=\left[\begin{array}[]{c}a^{1}\\ \hline\cr\vdots\\ \hline\cr a^{m}\end{array}\right][b_{1}|\dots|b_{n}]=\left[\begin{array}[]{c | c | c}a^{1}b_{1}&\dots&a^{1}b_{n}\\ \hline\cr\vdots&\ddots&\\ \hline\cr a^{m}b_{1}&&a^{m}b_{n}\end{array}\right].$$

While mathematically these four expressions compute exactly the same matrix product, in terms of performance they might differ substantially. From this perspective, slicings 1–3 lead to level 2 BLAS kernels (GEMVs for 1 and 2, GERs for 3) and are thus to be preferred to slicing 4, which instead leads to level 1 kernels (DOT). In the following sections we will follow the same procedure for higher-dimensional tensors and multiple contractions.

Due to their high-dimensionality, tensors are usually presented together with their indices. As commonly done in physics related literature, in the rest of the paper we use greek letters to denote indices. We now illustrate some concepts which are not strictly necessary to understand the methodology introduced in the following sections, but that are important to establish a connection with the language commonly used by physicists and chemists.

As explained in Appendix A, there is a substantial difference between upper (contravariant) and lower (covariant) indices. In general, the connection between covariant and contravariant indices is given by a special tensor $g_{\mu\nu}$ called “metric”. For example, an upper index of a vector $x$ can be lowered by multiplying the vector with the metric tensor (in itself having properties similar to those of a symmetric matrix):

Notice that the right hand side is written without explicit reference to the sum over the $\mu$ index. This convention is usually known as the “Einstein notation”: Repeated indices come in couples (one upper and one lower index) and are summed over their entire range. We refer to this operation as a “contraction”. It is important to keep in mind that contractions happen only between equal indices lying on opposite positions. Following this simple rule, matrix-vector multiplication is written as

Since a tensor can have more than one upper or lower index, we interpret lower indices as generalized “column modes” and the upper indices as generalized “row modes”. In general, it is common to refer to a $k$-dimensional tensor having $i$ row modes and $j$ column modes (with $i+j=k$) as a ($i,j$)-tensor or an “order-$k$ tensor”. Finally, since relative order among indices is important, the position of an index should be always accounted for; this is of particular importance in light of the results we present.

Having established the basic notations and conventions (see Appendix A for more details), we are ready to study numerical multilinear operations between tensors. We begin by looking at an example where only one index ($\gamma$) is contracted. In the language of linear algebra, the operation can be described as an inner product between a (3,0)-tensor and a (2,0)-tensor

In accordance with the convention, the summation over $\gamma$ can only take place once the index is lowered in one of the two tensors. For example, $\gamma$ can be lowered in $T$ as

The inner product (1) can then be written as a left multiplication of the last column-mode of $T$ with the first row-mode of $S$. Since $\gamma$ could have been lowered in $S$ (instead of $T$), mathematically this operation is exactly equal to a right multiplication between the last row-mode of $T$ and the first column-mode of $S$⁷⁷7 Once the tensor indices are given explicitly, the order of the operands does not affect the outcome: $T\indices{{}^{\alpha\beta}_{\gamma}}S^{\gamma\rho}\equiv S^{\gamma\rho}T\indices{{}^{\alpha\beta}_{\gamma}}$.

The inner product (3) constitutes only one way of contracting one index of $T$ with one of $S$. In general, depending on which index is contracted, there are six mathematically distinct cases:

Since the index sizes can differ substantially, the distinction among all these cases is not just formal. Additionally, a tensor can be symmetric with respect to the interchange of the first and third index but not with respect to the interchange between the first and the second. Consequently the mutual ordering among indices is quite important prompting to distinguish among cases with different ordering. Thus not only matters which index is contracted, but also in which position each index is. Ultimately, a tensor is a numerical object stored in the physical memory of a computer. We will see that index structure and storage strategies play a fundamental role in the optimal use of the BLAS library.

Here we discuss how the procedure used in Sec. (2.2) to slice matrices can be generalized to tensors, to produce objects manageable by BLAS. In general, an N-dimensional tensor is sliced by (N-1)-dimensional objects, and every slice is itself a tensor. Notice that the slicing language shifts the focus from the parts of the partitioned tensor to the partitioning object itself. As we will see in the next section, this shift in focus is at the core of our analysis.

To allow an easy visualization, let us again consider the contraction case (i) of Eq. (4). In order to reduce this operation to a sequence of matrix-matrix multiplications, the (2,1)-tensor $T\indices{{}^{\alpha\beta}_{\sigma}}$ needs to be sliced by a 2-dimensional object; such a slicing can be realized in several ways. For the sake of simplicity, let us represent $T$ as a 3-dimensional cube, and consider its slicing by one or more planes. Figure 1 shows three (out of seven) possible slicings of $T$ by different choices of planes.

• (a)

  The cube is sliced by a plane defined by the modes $\alpha$ and $\sigma$.

• (b)

  The cube is sliced by a plane defined by the modes $\beta$ and $\sigma$.

• (c)

  This double slicing is defined by combining the two previous slicing, (a) and (b).

A simple yet effective way of describing the slicing of a tensor is by means of a vector $\vec{s}$ that indicates the indices that are sliced:

In this language, the previous three cases correspond to

• (a)

  $\vec{s}_{a}(T)=(s_{\alpha},s_{\beta},s_{\sigma})=$ (0, 1, 0);

• (b)

  $\vec{s}_{b}(T)=$ (1, 0, 0);

• (c)

  $\vec{s}_{c}(T)=$ (1, 1, 0); this is the sum of the previous two slicings, $\vec{s}_{a}(T)+\vec{s}_{b}(T)=\vec{s}_{c}(T)$.

In the first case, $s_{\beta}=1$ means that the index $\beta$ is decomposed in a sequence of as many slices as the range of $\beta$. Accordingly, $s_{\alpha}$ and $s_{\sigma}$ take value 0, i.e., the slicing is made by planes defined by the modes $\alpha$ and $\sigma$ along the $\beta$ direction. Note that any combination of slicings results in an object of dimensionality $N-\|\vec{s}(T)\|_{1}$. For instance, slicing (c) produces a set of $3-2=1$-dimensional objects.

We now relate the slicing to mappings onto BLAS routines. To this end, we have to specify how the tensors are stored in memory. In the rest of the paper we always assume the first index of each tensor to have $stride=1$. The same way in Sec. 2.2 we decomposed the contraction between $A$ and $B$, in the following we decompose the contraction $R\indices{{}^{\alpha\beta\rho}}:=T\indices{{}^{\alpha\beta}_{\sigma}}\;S^{\sigma\rho}$ in three different ways. Here the tensors $T$, $S$ and $R$ are of size $m\times n\times k$, $k\times\ell$, and $m\times n\times\ell$, respectively. (Then the slices $R^{\alpha i\rho}$ are of size $m\times\ell$, and $T\indices{{}^{\alpha i}_{\sigma}}$ are of size $m\times k$.)

• (a)

  $\vec{s}_{a}(T)=$ (0, 1, 0) leads to a sequence of matrix products directly mappable onto GEMMs since the dimensions with stride=1 remain unsliced: $R=\forall_{i}\ T_{i}S$.
• (b)

  $\vec{s}_{b}(T)=$ (1, 0, 0) also leads to a sequence of matrix products; in this case, however, to enable the use of GEMM, each slice of $T$ must be copied on a separate memory region: $R=\forall_{i}\ (T_{i})_{cp}\ S$.
• (c)

  $\vec{s}_{c}(T)=$ (1, 1, 0) leads to a sequence of GEMVs: $R=\forall_{i}\forall_{j}\ T_{ij}S$.
  

Once again, these three slicings lead to mathematically equivalent algorithms for computing the tensor $R$. In terms of performance, (a) is to be preferred, as it directly maps onto matrix products. Although slicing (b) also reduces the problem to a sequence of matrix-matrix products, since the slices have both dimensions with stride $>1$, these cannot be directly mapped onto GEMM; this means that a costly preliminary transposition is necessary. Finally, slicing (c) should be avoided, as it forms $R$ by means of matrix-vector products, thus relying on the slow BLAS 2 kernel.

As illustrated in Sec. 2.4, a slicing transforms a contraction between two generic tensors $T$ and $S$ into a sequence (or a sum) of contractions between the slices of $T$ and those of $S$. To allow a mapping onto BLAS 3, both the slices of $S$ and $T$ need to conform to the GEMM interface (see Sec. 2.1). In the first place, they have to be matrices; hence $S$ and $T$ have to have at least $N=2$ indices, and if $N>2$, the tensor has to be sliced along $N-2$ of them. Moreover, the slices of $S$ and $T$ must have exactly one contracted and one free (non-contracted) index; in other words, we are looking for two-dimensional slices that share one dimension, the same way in the contraction $C:=AB$, the matrices $A\in\mathbf{R}^{m\times k}$ and $B\in\mathbf{R}^{k\times n}$ share the $k$-dimension, while $m$ and $n$ are free. Finally, the BLAS interface does not allow the use of matrices with stride $>1$ in both dimensions; therefore it is not possible to slice the indices of $S$ and $T$ that correspond to stride 1 modes without compromising the use of GEMM.

The previous considerations can be formalized as three requirements. Define $N(X)$ as the number of modes of tensor $X$, and let $R:=TS$ be a contraction.

• R1.

  The modes of $T$ and $S$ with stride 1 must not be sliced.

• R2.

  $T$ and $S$ have to be sliced along $N(T)-2$ and $N(S)-2$ modes, respectively. (In a compact language, $\|\vec{s}(T)\|_{1}=N(T)-2$ and $\|\vec{s}(S)\|_{1}=N(S)-2$). Notice that whenever a contracted index is sliced, both $T$ and $S$ are affected.

• R3.

  The slices must have exactly one free and one contracted index. Therefore, it must be avoided to slice only all the contracted indices or only all the free ones.

Despite their simplicity, these three requirements apply universally to contractions $R=TS$ of any order between tensors of generic dimensionality. If they are all satisfied, then it is possible to compute $R$ via GEMMs; this is for instance the case of Example (a) in Sec. 2.4.

A precise classification of contractions is given in the next section. Let us now briefly discuss the scenarios in which exactly one of the requirements does not hold.

• F1.

  If requirement R1 is not satisfied (the resulting slice does not have any mode with stride 1), the use of GEMM is still possible by copying each slice in a separate portion of memory. We refer to this operation as a COPY + GEMM, an example of whom is given by case (b) in Sec. 2.4.

• F2.

  If $\|\vec{s}(T)\|_{1}>N(T)-2$ or $\|\vec{s}(S)\|_{1}>N(S)-2$, then it is not possible to use BLAS 3 subroutines; in this case, one has to rely on routines from level 2 or level 1, such as GEMV and GER, or even on no BLAS at all. An example is provided by case (c) in Sec. 2.4.

• F3.

  As for F2, if all the contracted indices or all the free indices are sliced, one can use only BLAS 2, BLAS 1, or even no BLAS at all. An example is given by the contraction $R\indices{{}^{\alpha\beta\rho}}:=T\indices{{}^{\alpha\beta}_{\sigma}}\;S^{\sigma\rho}$ with $\vec{s}(T)=(0,0,1)$: the computation reduces to a sequence of element-wise operations.

We now use the requirements for cataloging all contractions in a small set of classes; most importantly, for each class we provide a simple recipe, indicating how to achieve an optimal use of BLAS. The aim is to empower computational scientists with a precise and easy-to-implement tool: Ultimately, these recipes shed light on the importance of choosing the correct strategy for generating and storing tensors in memory.

Given two tensors $t_{1}$ and $t_{2}$, and the number of contractions $p$, let us define $\Delta(t_{i})$ as $N(t_{i})-p$, with $i=1,2$ ($\Delta$ counts the number of indices of $t_{i}$ minus the number of contracted indices, that is, the number of free indices). The following classification is entirely based on $\Delta(t_{1})$ and $\Delta(t_{2})$.

Independently of the order of the tensors $t_{1}$ and $t_{2}$ and of the number of contractions, our classification by means of $\Delta(t_{i})$ is a general tool to easily determine whether GEMM is usable or not. For a given operation, the user can immediately determine the class to which it belongs. If both $\Delta(t_{1})$ and $\Delta(t_{2})$ equal 0, it is only possible to use DOT. If exactly one between $\Delta(t_{1})$ and $\Delta(t_{1})$ equals 0, then BLAS 2 operations are usable. Finally, the most interesting scenario is when both $\Delta(t_{1})$ and $\Delta(t_{2})$ are greater or equal to 1: In this case there is always the opportunity of using GEMM, provided the tensors can be generated or stored in memory suitably. In this respect, the recipes also provide guidelines on how tensors should be arranged for maximum performance.

If the user has no freedom in terms of generation and storage of the tensors, GEMM is still amenable, but at the cost of an overhead caused by copying and transposing. Whenever possible, this option should be avoided. Nonetheless, in general COPY + GEMM still performs better than BLAS 2 kernels; moreover, the transpositions can be optimized to reduce the number of memory accesses and therefore the total overhead. In the next section, we present the performance that one can expect for each of the classes.

We consider the contractions $R\indices{{}^{\alpha}_{\eta}}:=A^{\alpha\beta\gamma}B_{\eta\beta\gamma}$ and $R\indices{{}^{\beta}_{\eta}}:=A^{\alpha\beta\gamma}B_{\gamma\alpha\eta}$, and the operations resulting from the slicings depicted in Fig. 2 and 3.¹⁰¹⁰10 Since $A$ and $B$ are cubic tensors, these two contractions perform the same number of arithmetic operations. Moreover, if the tensors in the second contraction are stored in such a way that a) $A$ is transposed so that the first and the second indices are exchanged, and b) $B$ is transposed so that the first and the third indices are exchanged, then they are also computing the same quantity, i.e., $R\indices{{}^{\alpha}_{\eta}}=R\indices{{}^{\beta}_{\eta}}$. Specifically, in Figg. 4 and 5 we show the execution times for the computation of $R$ when using GEMM (in two different ways), COPY + GEMM, GEMV, and GER.

As expected, the two slicings leading to GEMM outperform considerably those that instead make use of BLAS 2 routines. In fact, for tensors of size 150 and larger, each routine attains the efficiency level mentioned in the introduction: 90+% of the theoretical peak for GEMM, and between 10% to 20% of peak for BLAS 2 routines. This is an example of the clear benefits due to the use of BLAS-3 routines instead of those from BLAS 2 and 1. Additionally, the plots show that when COPY + GEMM is used, for tensors of size 200 or larger the overhead due to copying becomes apparent.¹¹¹¹11 While it is possible to optimize the copying, amortizing the cost over multiple slices, here we present the timings for a simple slice-by-slice copy, as a typical user would implement. In short, the results are in agreement with the discussion from Sec. 3.3: Whenever the contraction falls in the third class, the use of GEMM is preferable. If the dimensions are small, there is little difference between COPY + GEMM and GEMM; for mid and large sizes instead, it is always convenient to store the tensors so as to enable the direct use of GEMM, thus avoiding any overhead (class 3.2).

Figure 5 illustrates the results for a parallel execution on all the eight cores of tesla, using multi-threaded BLAS. It is not surprising to observe the increase of the gap between GEMM and BLAS 2 based routines. Moreover, comparing Fig. 4 and Fig. 5, one can also appreciate the larger performance gap between GEMM and COPY + GEMM, confirming the fact that on multi-threaded architectures it is even more crucial to carefully store the data to avoid unnecessary overheads. In fact, as explained in Sec. 3.3, the tensors with $\Delta(t_{i})\geq 1$ should be stored in memory so that the dimension with stride 1 is not one of the contracted indices; this choice prevents the use of COPY + GEMM and the consequent loss in performance. Finally, we point out that even though the two GEMM-based algorithms traverse the tensors in different ways, they attain the exact same efficiency.

In conclusion, with this first example we confirmed that when performing double contractions on square tensors, any slicing that leads to GEMM  constitutes the most efficient choice, The situation does not change as the order of the tensors and the number of contracted indices increase.

We shift the focus to real world scenarios in which the tensors are rectangular and possibly highly skewed. The objective is to test the validity of our recipes in situations, inspired by scientific applications, more general than contractions between square tensors. For this purpose, we concentrate on double contractions between 4-dimensional tensors. We present here two numerical examples: The first arises in the Coupled-Cluster method, while the second is motivated by covariant operations appearing in General Relativity. All these tests refer to the tensors $X$ and $Y$ presented in Sec. 3.2.3, and the slicings described for Eq. (8). In particular, $R^{\alpha\beta\gamma\delta}$ was computed via four algorithms (resulting from four different slicings), using GEMM, COPY + GEMM, DOT, and GER, respectively.¹²¹²12In this case the indices contracted remain the same, consequently $R^{\alpha\beta\gamma\delta}:=X^{i\alpha j\beta}Y\indices{{}_{i}^{\gamma}{}_{j}^{\delta}}$ and $X^{i\alpha j\beta}Y\indices{{}_{j}^{\gamma}{}_{i}^{\delta}}$ involve exactly the same number of operations.