Code covered by the BSD License  

Highlights from
MMX - Multithreaded matrix operations on N-D matrices

5.0

5.0 | 3 ratings Rate this file 25 Downloads (last 30 days) File Size: 460 KB File ID: #37515
image thumbnail

MMX - Multithreaded matrix operations on N-D matrices

by Yuval

 

16 Jul 2012

N-D matrix operations. Like Boettcher's ndfun or Tursa's mtimesx, only faster.

| Watch this File

File Information
Description

MMX treats an N-D matrix of double precision values as a set of pages
of 2D matrices, and performs various matrix operations on those pages.
MMX uses multithreading over the higher dimensions to achieve good
performance. Full singleton expansion is available for most operations.

C = MMX('mult', A, B) is equivalent to the matlab loop
for i=1:N,
    C(:,:,i) = A(:,:,i) * B(:,:,i);
end
Singleton expansion is enabled on all dimensions so for example if
A = randn(5,4,3,10,1);
B = randn(4,6,3,1 ,6);
C = zeros(5,6,3,10,6);
then C = mmx('mult',A,B) equivalent to
for i = 1:3
   for j = 1:10
      for k = 1:6
         C(:,:,i,j,k) = A(:,:,i,j,1) * B(:,:,i,1,k);
      end
   end
end

C = MMX('mult', A, B, mod) and where mod is a modifier string, will
transpose one or both of A and B. Possible values for mod are
'tn', 'nt' and 'tt' where 't' stands for 'transposed' and 'n' for
'not-transposed'. For example
>> size(mmx('mult',randn(4,2),randn(4,2),'tn'))
ans = 2 2

C = MMX('square', A, []) will perform C = A*A'
C = MMX('square', A, [],'t') will perform C = A'*A

C = MMX('square', A, B) will perform C = 0.5*(A*B'+B*A')
C = MMX('square', A, B, 't') will perform C = 0.5*(A'*B+B'*A)

C = MMX('chol', A, []) will perform C = chol(A)

C = MMX('backslash', A, B) will perform C = A\B
Unlike other commands, 'backslash' does not support singleton
expansion. If A is square, mmx will use LU factorization, otherwise it
will use QR factorization. In the underdetermined case, (i.e. when
size(A,1) < size(A,2)), mmx will give the least-norm solution which
is equivalent to C = pinv(A)*B, unlike matlab's mldivide.

C = MMX('backslash', A, B, 'U') or MMX('backslash', A, B, 'L') will
perform C = A\B assuming that A is upper or lower triangular,
respectively.

C = MMX('backslash', A, B, 'P') will perform C = A\B assuming that A
is symmetric-positive-definite.

MMX(n) does thread control: mmx will automatically start a number of
threads equal to the number of available processors, however the
number can be set manually to n using the command mmx(n). mmx(0) will
clear the threads from memory.

IMPORTANT NOTE: The functions which assume special types of square
matrices as input ('chol' and 'backslash' for 'U','L' or 'P'
modifiers) do not check that the inputs are indeed what you say they
are, and produce no error if they are not. Caveat computator.

COMPILATION: To compile run 'build_mmx'. Type 'help build_mmx' to read
about compilation issues and options

Acknowledgements

Mtimesx Fast Matrix Multiply With Multi Dimensional Support inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.14 (R2012a)
Other requirements Works better if you can link it to a single-threaded BLAS library
Tags for This File  
Everyone's Tags
backslash, multiple multiplication, multiply, multiprod, multithreaded
Tags I've Applied
Add New Tags Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (10)
04 Jun 2013 Alexandros Iliopoulos

Dear Yuval,

Indeed, I'm happy to report that the numerical errors were not due to your package, after all.

Interestingly, the errors I was having were not due to MMX giving me the pseudo-inverse solution (which is what I wanted), but rather they were due to MATLAB's pinv() function, whose output I passed to MMX. In particular, I found out that MATLAB does not use symmetric-matrix-specific algorithms for functions such as pinv, svd, or sqrt, which may result in output matrices having non-zero imaginary terms. Of course, this did mess up the MMX output (as noted earlier), but the input was wrong to begin with.

After going around that problem, MMX worked like a charm!

04 Jun 2013 Alexandros Iliopoulos  
04 Jun 2013 Yuval

Alexandros,

Thank you so much for your feedback. This is not a bug. As described in the documentation, for rank deficient matrices mmx gives the least-norm (pseudoinverse) solution rather than the maximally-sparse solution. You didn't specify how your matrices were generated but clearly they are not 'in general position'. Try using random A's and you will see no discrepancy (to machine precision). Alternatively, compare mmx's output to pinv(A)*B.

03 Jun 2013 Alexandros Iliopoulos

Dear Yuval,

It appears to me that your package suffers from numerical artifacts, at least with respect to 'BACKSLASH' computations.

To check, I ran something like this:
% A:3x3x1850, B:3x1x1850
C = mmx( 'backslash', A, B );
D = zeros( size( C ) );
for k = 1 : size(D,3)
D(:,:,k) = A(:,:,k) \ B(:,:,k);
end
E = sum( (C - D).^2, 1 );
stem( permute( E, [1 3 2] ) )

It is readily visible that E is zero everywhere except for 4-5 indices, where the error is pretty large.

Unfortunately, I haven't been able to consruct a MWE with rand() data or something similar. Please PM me if you want me to send you the data that gave me the errors.

I hope that you have the time to work on this issue (assuming that I haven't done something silly), as I think that this is a very handy package to have! I'd be willing to help if I can.

03 Jun 2013 Alexandros Iliopoulos  
12 Apr 2013 Yuval

Thanks Adam!
Please spread the word.
Adding complex number support should be pretty easy.
If I see that mmx gets more traction in the community and that more people want it i'd be happy to do it myself.

As it stands it seems that it's buried in a dark corner of file exchange...

12 Apr 2013 Adam Cooman

The package works great. It all compiled nicely and is a lot faster than everything else out there. Also, the backslash option is really usefull.

Watch out when passing complex numbers though, the function only works on the real part of the data. Would it be possible to extend the function to complex matrices too?

19 Sep 2012 Yuval

Jan, look in the comments inside build_mmx.m

11 Sep 2012 Jan Simon

Do you have explicite suggestions for "linking to a single-threaded BLAS library"?

16 Jul 2012 Yuval

This package's homepage is here
http://www.cs.washington.edu/people/postdocs/tassa/code/#mmx

Contact us