The Vectorizer program for matlab tensor calculations
download (an executable; compiled on a utoronto cslab machine)
Problem description
It is often necessary to compute sums of products of high-dimensional matrices (tensors) in matlab.
Examples
- C(i)=sum_j A(i,j) B(j)
- C(i,j)=sum_k A(i,k)B(k,j)
The above expressions are easily implemented in matlab, as simple matrix multiplications. However, there are many simple cases that
matlab cannot handle easily:
- C(i,j,w)=sum_k A(i,k,w)B(k,j,w)
- C(i,j,:)=A(i,:)-B(j,:)
- C(a,b,c,d,e,f,g)=sum_k A(a,k,b) * sum_j (B(c,d,j,k)-D(j,k,f,g))*prod_{l,m} W(k,j,a,b,c,l,m)*2
The last expression is a useful example that describes the expressions in question. In general, we are
talking about nested sums products of matrices along certain dimensions, although the matrices could
also be added to one another. Matlab, unfortunately, does not support easy computation of such expressions.
The general recipe for computing such expressions is to recursively repmat every matrix, permute its dimensions,
perform element-wise multiplication (or addition or subtraction). This is precisely what the Vectorizer does.
It accepts expressions written in a rather easy to understand way, and outputs matlab code to compute
the answer. It is written in Haskell and the source code
is provided, compiled using ghc.
The above examples are translated to the Vectorizer language. They are similar, except
that the ranges of each variables must be specified. This is because being a syntactic program,
it does not know the dimensions of the matrices which are needed in order to reshape and repmat the matrices.
If a similar program is ever implemented inside matlab then specifying the ranges of variables will not be necessary.
The examples below should be analyzed,
because they reveal the operation of the program. It is not hard, since the program
does the intuitively obvious thing (to me, at least).
Warning: never use * in the vectorizer; always use .*. The reason is that the program always assumes that
all operations are elementwise, and treats them as "syntactic things"; it does not care what symbol you give it for the operation,
be it +, -, *, .*, .+, or anything else. Because of that, if * is used instead of .*, the code will almost surely produce an error,
but it may sometimes run incorrectly.
- C[size_i](i)=sum(j=size_j|A(i,j).*B(j))
produces
C=sum(A.*repmat(reshape(B,[1,size_j]),[size_i,1]),2);
- C[n,m](i,j)=sum(k=l|A(i,k).*B(k,j))
produces
C=sum(repmat(reshape(A,[n,1,l]),[1,m,1]).*repmat(reshape(permute(B,[2,1]),[1,m,l]),[n,1,1]),3);
- C[size_i,size_j,size_w](i,j,w)=sum(k=size_k|A(i,k,w).*B(k,j,w))
produces
C=permute(sum(repmat(reshape(permute(A,[1,3,2]),[size_i,size_w,1,size_k]),[1,1,size_j,1]).*
repmat(reshape(permute(B,[3,2,1]),[1,size_w,size_j,size_k]),[size_i,1,1,1]),4),[1,3,2]);
- C[size_i,size_j,yo](i,j,:)=A(i,:)-B(j,:)
produces
C=repmat(reshape(A,[size_i,1,yo]),[1,size_j,1])-repmat(reshape(B,[1,size_j,yo]),[size_i,1,1]);
- C[s_a,s_b,s_c,s_d,s_e,s_f,s_g](a,b,c,d,e,f,g)=sum(k=s_k|A(a,k,b) .*
sum(j=s_j|(B(c,d,j,k)-D(j,k,f,g)).* prod(l=s_l,m=s_m| W(k,j,a,b,c,l,m).*2())))
produces
(Note the brackets after the 2! Their absence could be the source of mysterious syntax errors)
C=repmat(reshape(sum(repmat(reshape(permute(A,[1,3,2]),[s_a,s_b,1,1,1,1,s_k]),[1,1,s_c,s_d,s_f,s_g,1]).*
permute(sum((repmat(reshape(permute(B,[1,2,4,3]),[s_c,s_d,s_k,1,1,1,1,s_j]),[1,1,1,s_f,s_g,s_a,s_b,1])-
repmat(reshape(permute(D,[2,3,4,1]),[1,1,s_k,s_f,s_g,1,1,s_j]),[s_c,s_d,1,1,1,s_a,s_b,1])).*
repmat(reshape(permute(prod(prod(W.*2,7),6),[5,1,3,4,2]),[s_c,1,s_k,1,1,s_a,s_b,s_j]),
[1,s_d,1,s_f,s_g,1,1,1]),8),[6,7,1,2,4,5,3]),7),[s_a,s_b,s_c,s_d,1,s_f,s_g]),[1,1,1,1,s_e,1,1]);
Explanation of the operation of the program
The program operates on a syntactic level, so the user must provide the program with the dimensions of
the matrix. This is done by the declaration C[m,n,...](i,j,...); In this case, the variable i
ranges from 1 to m,
j ranges from 1 to n, etc.
When we write sum(i=m, j=n,...|...), then i also ranges from 1 to m,
and the variable j ranges from 1 to n.
When a matrix indexing operation is written as A(i,j,...), it is assumed that
its first dimension is of size m and second
is of size n, depending on the range of the indecies. Every matrix is merely a symbol and the same letter
A can be used in different manners, although this is misleading.
In addition, the program does not understand references such as A(i,j,i,...). Writing this
will probably result in an error, or worse, incorrect code.
Some simple optimizations are implemented as well, namely, that prevent from making too large repmats
of the matrices in the inner-most sums. Unfortunately, I forgot the details.
It goes without saying that there is no guarantee that in some obscure case this program will fail
to produce an answer that does what expected, although I never encountered such a case.
Important Notes
- The program always outputs misleading syntax errors when it finds some; they are never to be believed.
- The program is only sometimes tolerant for white spaces, so they could also be the source of
otherwise mysterious syntax errors.
In other words, do not use whitespaces.
- The characters _ and : are legal variable names, so things like C(i,j,:)
are allowed, making the code look closer to that of matlab. This is identical to C(i,j,x)
if : is replaced by x throughout the code.
- If you are using this program and getting very peculiar syntax errors, please send me an email
and I may try to fix them.
home