# Efficiently computing the product of a multi-dimensional matrix (or tensor) and vectors

Update: Thank you very much for all of you who answered below. I'm studying each answer now. In the long term, I'm more interested in solutions that work for sparse tensors (sorry I should have mentioned this earlier), but I can use your ideas as a temporary solution (I will need to add more RAM for larger size problems.)

Let $T=(T_{ijk})$ be an $n\times n\times n$ tensor and $x=(x_1,x_2,\ldots,x_n)$, $y=(y_1,y_2,\ldots,y_n)$, $z=(z_1,z_2,\ldots,z_n)$ be $n$-dimensional vectors.

(Edit: as suggested by @WolfgangBangerth I should not use the $\otimes$ symbol that's usually used for outer product. I replaced it by $\otimes_i$, which means now the product of a tensor and a vector at mode $i$.)

The products \begin{align} a &= T\otimes_2 y\otimes_3 z\\ b &=T\otimes_3 z\otimes_1 x \\ c &=T\otimes_1 x\otimes_2 y \end{align} are $n$-dimensional vectors defined by: \begin{align} a_i &= \sum_{1\le j\le n}\sum_{1\le k\le n} T_{ijk}y_jz_k,\quad i=1,\ldots,n \\ b_j &= \sum_{1\le i\le n}\sum_{1\le k\le n} T_{ijk}x_iz_k,\quad i=1,\ldots,n \\ c_k &= \sum_{1\le i\le n}\sum_{1\le j\le n} T_{ijk}x_iy_j,\quad i=1,\ldots,n. \end{align}

I have an algorithm doing:

Repeat:

1. Compute $a$, update $x=f(a)$.
2. Compute $b$, update $y=g(b)$.
3. Compute $c$, update $z=h(c)$.

I would like to ask for a way (or toolboxes/libraries that help me) to efficiently compute $a,b,c$ at each iterations, without doing loops.

Language: C++ is preferred but Matlab is also acceptable if easier.

Thank you very much in advance for your help!

• I think you should just drop the $\otimes$ symbol for the product -- it is usually used for outer products, not inner products. – Wolfgang Bangerth Oct 29 '16 at 2:01
• @WolfgangBangerth: You are right, thanks. I edited it. – Khue Oct 29 '16 at 13:46

I did a little test with the Tensor package in Eigen C++ (release candidate 1 of version 3.3) http://eigen.tuxfamily.org/index.php?title=Main_Page

The path to the Tensor include file shows this package as "unsupported" and, presumably, this will still be the case in version 3.3. I do know that this package is rapidly reaching the "supported" stage. For more information on this, you could contact the Eigen developers directly if that is a concern.

My test code is pasted below. On my system, the Eigen versions of the tensor products you show below are around 3X faster than the brute force versions. The syntax of the operation seems fairly convenient to me.

#include <boost/timer.hpp>
#include <unsupported/Eigen/CXX11/Tensor>

const int n = 300;
Eigen::Tensor<double, 3> T(n, n, n);
Eigen::Tensor<double, 1> y(n), z(n);
int ii = 1;
for (int i = 0; i < n; i++) {
y(i) = i+1;
z(i) = i + 3;
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
T(i, j, k) = ii++;
}
}
}
//cout << T << endl;
typedef Eigen::IndexPair<int> IP;
Eigen::array<IP, 1> d10 = { IP(1, 0) }; // second index of T
boost::timer timer;
Eigen::Tensor<double, 1>  Tyz = T.contract(y, d10).contract(z, d10);
//cout << Tyz << endl;
printf("elapsed time=%8.3f\n", timer.elapsed());

// brute force
Eigen::VectorXd v(n);
v.setZero();
timer.restart();
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
v(i) += T(i,j,k)*y(j)*z(k);
}
}
}
//cout << v.transpose() << endl;
printf("elapsed time=%8.3f\n", timer.elapsed());

• Thanks a lot. I have compiled your code (Eigen 3.3 RC1) and got an error: 'd10' was not declared in this scope. I got a few errors involving the file CXX11Meta.h as well. Could you please check? And it seems to me that Eigen does not support sparse tensors yet. Sorry I should have mentioned that I'm more interested in sparse tensors :( – Khue Oct 29 '16 at 13:19
• Sorry about that! Should have checked the code again after editing! I changed this line in the code above to fix the error: Eigen::array<IP, 1> d10 = { IP(1, 0) }; // second index of T. Yes, I believe you will need a C++ compiler that supports C++11 to compile parts of the Tensor package (not true in general for the rest of Eigen). I used MS VS 2013 on windows so the C++11 support doesn't have to be complete. Eigen has quite good support for sparse matrices so hopefully sometime in the future they will support sparse tensors. – Bill Greene Oct 29 '16 at 13:59
• I have g++4.8 that supports C++11. I will contact the developpers. Thanks, Bill. – Khue Oct 30 '16 at 19:48

In Matlab you can do these operations in a vectorized way using the commands reshape, shiftdim, and permute. The essential idea is that contraction of a tensor with a vector is equivalent to matrix multiplication of that vector with an unfolded version of the tensor. For the first example in the question, the command is:

a = z*reshape(y*reshape(shiftdim(T,1), n, n*n), n, n);


where $y$, $z$, and $c$ are row vectors. This computes:

$$a_i = \sum_{j=1}^n \sum_{k=1}^n T_{ijk} y_j z_k.$$

The reason for using row vectors and multiplying them from the right (rather than using column vectors and multiplying them from the left) is that MATLAB stores tensors in column-first order, so when reshaping and multiplying, you access elements of the tensor in the same order as it is stored in memory. This way it is more efficient in terms of cache usage and so on.

Explanation:

Let's break it down as follows:

(5) a = z*                                              ;
(4)       reshape(                                 n, n)
(3)               y*
(2)                 reshape(               n, n*n)
(1)                         shiftdim(T,1)


(1) The shiftdim(T,1) shifts all the modes to the left, wrapping around. That is, $$T_{ijk} \rightarrow T_{jki}.$$

(2) The innermost reshape unfolds the 3-tensor into a $n$-by-$n^2$ matrix where the $j$ mode is exposed on the left: $$\underbrace{T_{jki}}_{3-\text{tensor}} \rightarrow \underbrace{T_{(j)(ki)}}_{\text{matrix}}.$$

(3) Multiplying on the left by $y$ contracts $y$ with the $j$ mode, creating a temporary row vector of length $n^2$, lets call it $p$, which is a vectorized version of the tensor after contracting it with the $y$.

$$p_{ki} = \sum_j y_j T_{jki}$$

(4) The second reshape refolds $p$ into a matrix to expose the $k$ mode:

$$\underbrace{p_{(ki)}}_{\text{row vector}} \rightarrow \underbrace{p_{(k)(i)}}_{\text{matrix}}.$$

(5) The final multiplication by $z$ contracts $v$ with the $k$'th mode:

$$a_k = \sum_k z_k p_{ki}$$

In general, you can perform any sequence of tensor-vector or tensor-tensor contractions by unfolding the tensors into matrices that expose the modes being contracted over, then performing matrix multiplication, then folding the result back up into a tensor again.

The vectorized operations in Matlab that allow you to do this are shiftdim and permute to shift or permute the modes of the tensor (to bring the desired modes to the front or the back, or do the reverse later), and reshape to unfold and fold the tensor.

In Python with numpy, you can use the built in commands tensordot or einsum, which do this under the hood.

Timing example:

Here's a timing example you can run to get you an idea of how much faster this is:

n = 250;
T = randn(n,n,n);
y = randn(1,n);
z = randn(1,n);
disp('Vectorized:')
tic
a1 = z*reshape(y*reshape(shiftdim(T,1), n, n*n), n, n);
toc

a2 = zeros(1,n);
disp('Unvectorized:')
tic
for ii=1:n
for jj=1:n
for kk=1:n
a2(ii) = a2(ii) + T(ii,jj,kk)*y(jj)*z(kk);
end
end
end
toc

disp('Error:')
disp(norm(a2-a1))


On my laptop this displays the following results:

Vectorized:
Elapsed time is 0.056645 seconds.
Unvectorized:
Elapsed time is 1.781048 seconds.
Error:
3.1182e-11

• Thanks a lot, Nick! Detailed, well explained answer! Does it work for sparse tensors? – Khue Oct 29 '16 at 13:45
• I don't think Matlab supports sparse tensors. The Sandia tensor toolbox, as mentioned on Wolfgang Bangerth's answer, does support sparse tensors, though unfortunately you have to fill out some forms and agree to a license to download it. It is written by Tammy Kolda, Brett Bader, and collaborators, who you may recognize as the authors of the famous SIAM review paper "Tensor Decompositions and Applications", sandia.gov/~tgkolda/pubs/pubfiles/TensorReview.pdf – Nick Alger Oct 29 '16 at 20:31
• Have just tested and it turns out that using that Matlab Toolbox is still slower than brute force C++ (in the sparse case) :( For dense tensors I think your solution is good enough. I accepted @BillGreene's answer because he was the first to answer (and to give a solution in C++). I wish I could accept more than one answer :P Thanks again. And +1. – Khue Oct 30 '16 at 19:54
• Incidentally, the idea of folding and unfolding tensors along successive modes is a key principle behind "Tensor Trains", wherein one unfolds the tensor along the first mode, does a low rank approximation, refolds/unfolds the right low rank factor, does a low rank approximation to it, and so on until there are no modes left. – Nick Alger Nov 1 '16 at 0:42
• Thank you for the great explanation! For future readers: A couple of other einsum implementations in MATLAB are shown here. Regarding multidimensional sparse arrays, there's a FEX submission that can be mentioned in this context - ndSparse. An example of using ndSparse for Kronecker products can be found here. – Dev-iL Aug 28 '19 at 18:05

I think you should look at Tammy Kolda's "Tensor Toolbox" for matlab. It has many of the kind of operations you are looking for implemented in efficient ways.

• Thanks. I have tried it but got an error when doing the multiplication (using the ttv function): Error using * MTIMES is not fully supported for integer classes. I guess I should contact the author for that. – Khue Oct 29 '16 at 13:34
• I intend to try Tensorlab as well: tensorlab.net. – Khue Oct 29 '16 at 13:42
• @khue -- yes, that's a question for the authors (after reading the manual, of course). – Wolfgang Bangerth Oct 31 '16 at 19:23
• Thanks. I could manage to use those toolboxes (by adding several double(.) to their code to bypass the integer errors). However, as I have already replied to @NickAlger, brute force C++ is still faster in my experiments. – Khue Nov 1 '16 at 11:20