Embracing tensors

[Laurent S]

Hi all,

I think Julia has great potential, and can, unlike Matlab, rid itself of its tunnel vision on matrices and embrace tensors from the ground up. Tensors, and tensor decompositions (such as those which generalize the matrix SVD/EVD/...), have been gaining traction in recent years, and supporting them properly would certainly increase Julia's adoption rate looking forward. Below is a list of some of my tensor-related suggestions/comments:

• repmat currently allows replication along the first two dimensions, the second argument could instead define by which multiplier to replicate along a vector of modes.
  Possibly even consider changing the name to "rep" or "reparr", since the method also applies to vectors and tensors, not just matrices. On a related note, I'm very glad Julia has broadcasting operators, this makes my life so much easier.
• Although I can't find a reference in the documentation, it seems Julia uses the work "rank" internally to describe the number of dimensions an array (or tensor) has. This could be confusing, rank already has a different definition for matrices, and disjoins into two different concepts for tensors (i.e., rank and multilinear rank). For tensors, the number of dimensions is referred to as the tensor's order in the literature.
• The implementation of rank should allow the user to set a tolerance, preferably a relative tolerance (unlike Matlab's absolute tolerance).
• For tensors, rank could return the tensor's multilinear rank, which is a vector of ranks (one rank for each mode). Or, it could be a separate function (e.g., mlrank). The matrix case is just a special case of mlrank where the column rank equals the row rank.
• If A and B are matrices of the proper dimensions, A*B and B*A are well-defined. Now let B be a tensor (say, three-dimensional), then A*B and B*A are still defined in the multilinear algebra sense, but not yet in Julia. For the former, the effect is that each column vector B[:,j,k] of B is mapped to C[:,j,k] = A*B[:,j,k]. For the latter, each row vector B[i,:,k] of B is mapped to C[i,:,k] = B[i,:,k]*A. It is really a straightforward generalization of the matrix-matrix multiplication to matrix-tensor multiplication.
• Related to the previous point, the * operator would only allow multiplying a tensor from the left and from the right. In general, B could also be multiplied by a matrix A along its third (or nth) dimension. I'm not sure what the best way to include such functionality is. A tmmul(Tensor,Matrix,dim) function perhaps? Or a new operator?
• Julia already has the Kronecker product. There is a related product called the Khatri--Rao product, which would also be nice to have since it is an essential operation in one of the most prominent tensor decompositions.
• normfro(T) and norm(T,"normfro") where T is a third- or higher-order tensor should return norm(T[:]) (currently returns an error).
• Thank you for the many nice improvements to basic functions like sum and broadcasting operators, I love that I can do things like sum(T,2:ndims(T)) .+ randn(1,5).
  However, here is an inconsistency I noticed (or is it desired behaviour?):
  A = randn(3,3)
  size(A[:,3]) # is a (3,) -- This is fine, because
  size(A[:,3:3]) # is (3,1) -- I can still have it retain singleton dimensions if I want.
  size(sum(A,2:2)) # is (3,1) -- As expected.
  size(sum(A,2)) # is (3,1) -- Why is this not (3,), as with the indexing examples?
  

Thanks for all the great work. I'm very much looking forward to interactive plotting.

Best,

Laurent

[Tim Holy]

Looks like a good list, I'm sure I'd get some value from some of these as
well. I suspect there's little reason not to have most or all of them, except
for the minor point that someone actually has to write the code!

Best,
--Tim

On Wednesday, May 29, 2013 03:10:36 PM Laurent S wrote:
> Hi all,
>
> I think Julia has great potential, and can, unlike Matlab, rid itself of
> its tunnel vision on matrices and embrace tensors from the ground up.
> Tensors, and tensor decompositions (such as those which generalize the
> matrix SVD/EVD/...), have been gaining traction in recent years, and
> supporting them properly would certainly increase Julia's adoption rate
> looking forward. Below is a list of some of my tensor-related
> suggestions/comments:
>

> - repmat currently allows replication along the first two dimensions,

> the second argument could instead define by which multiplier to replicate
> along a vector of modes.
> Possibly even consider changing the name to "rep" or "reparr", since the
> method also applies to vectors and tensors, not just matrices. On a
> related note, I'm very glad Julia has broadcasting operators, this makes my
> life so much easier.

> - Although I can't find a reference in the documentation, it seems Julia

> uses the work "rank" internally to describe the number of dimensions an
> array (or tensor) has. This could be confusing, rank already has a
> different definition for matrices, and disjoins into two different
> concepts for tensors (i.e., rank and multilinear rank). For tensors, the
> number of dimensions is referred to as the tensor's order in the

> literature<http://www.sandia.gov/~tgkolda/pubs/bibtgkfiles/TensorReview.pdf
> > .
> - The implementation of rank should allow the user to set a

> tolerance, preferably a relative tolerance (unlike Matlab's absolute
> tolerance).

> - For tensors, rank could return the tensor's multilinear rank, which is

> a vector of ranks (one rank for each mode). Or, it could be a separate
> function (e.g., mlrank). The matrix case is just a special case of mlrank
> where the column rank equals the row rank.

> - If A and B are matrices of the proper dimensions, A*B and B*A are

> well-defined. Now let B be a tensor (say, three-dimensional), then A*B
> and B*A are still defined in the multilinear algebra sense, but not yet in
> Julia. For the former, the effect is that each column vector B[:,j,k] of B
> is mapped to C[:,j,k] = A*B[:,j,k]. For the latter, each row vector
> B[i,:,k] of B is mapped to C[i,:,k] = B[i,:,k]*A. It is really a
> straightforward generalization of the matrix-matrix multiplication to
> matrix-tensor multiplication.

> - Related to the previous point, the * operator would only allow

> multiplying a tensor from the left and from the right. In general, B
> could also be multiplied by a matrix A along its third (or nth) dimension.
> I'm not sure what the best way to include such functionality is. A
> tmmul(Tensor,Matrix,dim) function perhaps? Or a new operator?

> - Julia already has the Kronecker product. There is a related product

> called the Khatri--Rao product, which would also be nice to have since it
> is an essential operation in one of the most prominent tensor
> decompositions.

> - normfro(T) and norm(T,"normfro") where T is a third- or higher-order

> tensor should return norm(T[:]) (currently returns an error).

> - Thank you for the many nice improvements to basic functions like sum

[Stefan Karpinski]

Agreed, these all seem like good suggestions. Would you be willing to create issues for these to get the ball rolling?

[Laurent S]

Sure, I'll start putting them on GitHub tomorrow.

Laurent

[Simon Byrne]

On Wednesday, 29 May 2013 23:10:36 UTC+1, Laurent S wrote:

• If A and B are matrices of the proper dimensions, A*B and B*A are well-defined. Now let B be a tensor (say, three-dimensional), then A*B and B*A are still defined in the multilinear algebra sense, but not yet in Julia. For the former, the effect is that each column vector B[:,j,k] of B is mapped to C[:,j,k] = A*B[:,j,k]. For the latter, each row vector B[i,:,k] of B is mapped to C[i,:,k] = B[i,:,k]*A. It is really a straightforward generalization of the matrix-matrix multiplication to matrix-tensor multiplication.
• Related to the previous point, the * operator would only allow multiplying a tensor from the left and from the right. In general, B could also be multiplied by a matrix A along its third (or nth) dimension. I'm not sure what the best way to include such functionality is. A tmmul(Tensor,Matrix,dim) function perhaps? Or a new operator?

I think what would be really useful for tensors would be a way of using Einstein sum notation. Numpy has the `einsum` function, however I dug into the code once, and it seemed that for tensors of order >= 3, it always ends up using simple loops, whereas in some cases it would be more efficient to use BLAS routines, as certain configurations are equivalent to reshaped matrix matrix multiplications (this was a while ago, things might have changed since then). 

I've been playing around with writing a macro to do this, so you could use something like

%einsum X[i,j] = A[i,k,l] * B[k,l,j] + Y[j] * Z[i]

which would do the relevant reshapings, but it's a bit tricky...

Simon

[Stefan Karpinski]

Long ago, I played around with Einstein summation notation for multidimensional comprehensions, but we ended up going with something much simpler instead. Doesn't mean we shouldn't have something like that though, accessible via a macro.

[Gustavo Lacerda]

If you look back on the archive, I started a small thread proposing a Tensor package.  It may be worth reading.

G

--
--
Gustavo Lacerda
http://www.optimizelife.com