21
$\begingroup$

Consider this simple matrix number multiplication:

lth = 200;
mtx = RandomReal[{0, 1}, {lth, lth}];
ls = RandomReal[{0, 1}, {lth}];

Et = Function[{t}, Sin[(π t)/20] Sin[2 t]];
Etc = Compile[{{t, _Real}}, Et[t], 
   CompilationOptions -> {"InlineCompiledFunctions" -> True, 
     "InlineExternalDefinitions" -> True}];


Table[Etc[t]*mtx;, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.244659, Null} *)

It seems that this is very slow, compared to matrix vector multiplication:

Table[mtx.ls;, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.038151, Null} *)

Moreover, the matrix addition also seems to be slow

Table[mtx + mtx;, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.153648, Null} *)

Question: So why the matrix number multiplication and addition so much slower than the matrix vector multiplication? Are there ways to speed them up?

I'm using 10.3 on OS X 10.11.4.


Edit

The slowness of the matrix number multiplication doesn't seem to come from Etc, for example:

Table[Etc[t];, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.001614, Null} *)

Table[1.*mtx;, {t, 0., 20, 0.01}]; // AbsoluteTiming
(* {0.235871, Null} *)

Edit 2

Here is a comparison to Matlab:

lth=200;
mtx=rand(lth);
ls=rand(lth,1);


tic;
for t=0:0.01:20
mtx2=1.*mtx;
end
toc


tic;
for t=0:0.01:20
mtx2=mtx*ls;
end
toc
Elapsed time is 0.034530 seconds.
Elapsed time is 0.015745 seconds.

Mathematica is as fast as Matlab in matrix vector multiplication, but about 7X slower in the matrix number multiplication.


Edit 3

Here are more detailed comparisons between Mathematica and Matlab for matrix number multiplication and matrix vector multiplication, for different problem size.

Compare α*A, where α is a scalar and A is a matrix

enter image description here

Here is a comparison of for the spare array method suggested by Anton Antonov:

enter image description here

It seems that even using the sparse array method, Mathematica is still way too slow than Matlab in this operation, especially looking at cases with large matrix size.

Compare A.v, where v is a vector and A is a matrix

enter image description here

It seems that Mathematica has a much closer performance to Matlab in the matrix vector mulplitcation than that of the matrix number multiplication.

I'm happy about the <2X slowness in the matrix vector multiplication, but the slowness of the matrix number multiplication seems to me like something is not working properly, given that this operation is so fundamental that it should have been optimized at an early stage.

code used for comparison

(*===Mathematica functions ====*)
matrixScaling[lth_] := Module[{mtx},
  mtx = RandomReal[{0, 1}, {lth, lth}];
  First@AbsoluteTiming[Table[1.*mtx;, {2000}];]/2000/lth
  ]
matrixScaling2[lth_] := Module[{mtx, sp},
  mtx = RandomReal[{0, 1}, {lth, lth}];
  First@AbsoluteTiming[
      Table[sp = SparseArray[{_, _} -> 1., Dimensions[mtx]]; 
        sp*mtx;, {2000}];]/2000/lth
  ]
matrixDot[lth_] := Module[{mtx, ls},
  mtx = RandomReal[{0, 1}, {lth, lth}];
  ls = RandomReal[{0, 1}, {lth}];
  First@AbsoluteTiming[Table[mtx.ls;, {2000}];]/2000/lth
  ]


(*===Matlab functions ====*)
Needs["MATLink`"]
OpenMATLAB[]

FilePrint["~/Documents/MATLAB/matrix_scale.m"]
(*
function time = matrix_scale( lth )

mtx=rand(lth);
tic;
for i=0:2000
mtx2=1.*mtx;
end
time=toc;
time=time/lth/2000;
end
*)

FilePrint["~/Documents/MATLAB/matrix_dot.m"]
(*
function time = matrix_dot( lth )

mtx=rand(lth);
ls=rand(lth,1);
tic;
for i=0:2000
mtx2=mtx*ls;
end
time=toc;
time=time/lth/2000;
end
*)

(*====comparison ====*)
nls = Range[100, 1000, 100];
timeScaleMMA = matrixScaling /@ nls;
timeScaleMMA2 = matrixScaling2 /@ nls;
timeDotMMA = matrixDot /@ nls;

matrixScaleMatlab = MFunction["matrix_scale"];
matrixDotMatlab = MFunction["matrix_dot"];
timeScaleMatlab = matrixScaleMatlab /@ nls;
timeDotMatlab = matrixDotMatlab /@ nls;

(*===plot results ===*)

BarChart[1.*^6 Transpose@{timeScaleMMA, timeScaleMatlab}, 
 ChartLabels -> {Range[100, 1000, 100], None}, 
 PlotTheme -> "Detailed", 
 FrameLabel -> {"matrix size", "time (μs)"}, 
 ChartLegends -> {"Mathematica", "Matlab"}, ImageSize -> 400, 
 AspectRatio -> 1/GoldenRatio]

BarChart[1.*^6 Transpose@{timeScaleMMA, timeScaleMMA2, 
    timeScaleMatlab}, ChartLabels -> {Range[100, 1000, 100], None}, 
 PlotTheme -> "Detailed", 
 FrameLabel -> {"matrix size", "time (μs)"}, 
 ChartLegends -> {"Mathematica", "Mathematica sparse array", 
   "Matlab"}, ImageSize -> Medium]

BarChart[1.*^6 Transpose@{timeDotMMA, timeDotMatlab}, 
 ChartLabels -> {Range[100, 1000, 100], None}, 
 PlotTheme -> "Detailed", 
 FrameLabel -> {"matrix size", "time (μs)"}, 
 ChartLegends -> {"Mathematica", "Matlab"}, ImageSize -> 400, 
 AspectRatio -> 1/GoldenRatio]
$\endgroup$
7
  • $\begingroup$ The slowness of the first Table is probably due to the calling of the function Etc, which calculates sinusoids, and not to the multiplication by the matrix. $\endgroup$
    – bill s
    May 14, 2016 at 13:33
  • $\begingroup$ @bills it doesn't seem to come from calling Etc, see update. $\endgroup$ May 14, 2016 at 15:03
  • $\begingroup$ The matrix number multiplication is 2.2 times slower using version 9.0.1 than it is using version 10.4.1. $\endgroup$
    – Karsten7
    May 15, 2016 at 1:07
  • $\begingroup$ On my MacBook Pro, 3GHz Intel Core i7 with MMa 10.0.2 the matrix multiplication is three times slower where matrix vector mlutiplication is slightly faster than in your case. I found that the first m.v takes .6x10^-6 where the subsequent calcualtions take .12x10^-6 so memoization seems to take place where this appears not to be the case with the matrix multiplication which has a consistent timing of approx .7x10^-5. $\endgroup$
    – Sander
    May 15, 2016 at 6:31
  • $\begingroup$ @Sander I'm guessing the extra time for the first run is due to loading some fast library function(for example BLAS) in the matrix vector multiplication. That's also consistent with the fact that the matrix vector multiplication is comparable with Matlab speed, since it seems that both of then use BLAS. $\endgroup$ May 15, 2016 at 6:49

4 Answers 4

11
$\begingroup$

Mr.Wizard proposed that the slowness is due to the copying of the data. It seems that this intuition is correct. We can use the two different arguments passing mechanism in LibraryLink to test this conjecture.

In LibraryLink, there are the "copied" passing and the "shared" passing. The copied passing copies the data from Mathematica to the library function while the shared passing shares the data between them. If the slowness really comes from the copying of the data, then we would expect a huge performance difference between these two way of passing arguments.

This plot compares the original Mathematica calculation, LibraryLink copied passing and LibraryLink shared passing for calculating 1.*A, where A is a matrix. We can see that indeed the copied passing scales very close to Mathematica performance, while the shared passing scales more like Matlab performance.

enter image description here


Code for the tests

Here are the code I used for the tests.

Fortran code that multiply a matrix and a number

!matscale.f90
subroutine matscale(num,a,n) bind(c)
implicit none
integer      :: n
real (kind=8) :: num
real (kind=8) :: a(n,n)
a=num*a(:,:)
end subroutine matscale

LibraryLink wrapper for copied passing

//array_copy.cpp
#include "/Applications/Mathematica.app/Contents/SystemFiles/IncludeFiles/C/WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
  return 0;
}


extern "C"{ void matscale(double* num, double a[], int* n);}

EXTERN_C DLLEXPORT int mymatscale(WolframLibraryData libData, mint \
Argc, MArgument* Args, MArgument Res)
{


double  num = MArgument_getReal(Args[0]);
MTensor _ta = MArgument_getMTensor(Args[1]);
double* a = libData->MTensor_getRealData(_ta);
int  n = MArgument_getInteger(Args[2]);


matscale(&num, a, &n);


MArgument_setInteger(Res, 0);
return LIBRARY_NO_ERROR;
}

LibraryLink wrapper for shared passing

The only difference is that we have an an extra libData->MTensor_disownAll(_ta); at the end:

//array_share.cpp
#include "/Applications/Mathematica.app/Contents/SystemFiles/IncludeFiles/C/WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
  return 0;
}


extern "C"{ void matscale(double* num, double a[], int* n);}

EXTERN_C DLLEXPORT int mymatscale(WolframLibraryData libData, mint \
Argc, MArgument* Args, MArgument Res)
{


double  num = MArgument_getReal(Args[0]);
MTensor _ta = MArgument_getMTensor(Args[1]);
double* a = libData->MTensor_getRealData(_ta);
int  n = MArgument_getInteger(Args[2]);


matscale(&num, a, &n);

libData->MTensor_disownAll(_ta);

MArgument_setInteger(Res, 0);
return LIBRARY_NO_ERROR;
}

Compilation of the library

ifort -fPIC -c -r8 matscale.f90
icc -fPIC -c array_copy.cpp
ifort -fPIC -lstdc++ -dynamiclib matscale.o array_copy.o -o array_copy.dylib
icc -fPIC -c array_share.cpp
ifort -fPIC -lstdc++ -dynamiclib matscale.o array_share.o -o array_share.dylib

Load library function

mtxCopy = 
  LibraryFunctionLoad[
   "array_copy.dylib", 
   "mymatscale", {Real, {Real, 2, Automatic}, Integer}, Integer];

mtxShare = 
  LibraryFunctionLoad[
   "array_share.dylib", 
   "mymatscale", {Real, {Real, 2, "Shared"}, Integer}, Integer];

Perform the test

matrixScaling[lth_] := 
 Module[{mtx}, mtx = RandomReal[{0, 1}, {lth, lth}];
  First@AbsoluteTiming[Table[1.*mtx;, {2000}];]/2000/lth]

matrixScaling2[lth_] := 
 Module[{mtx, \[Alpha] = 1.}, mtx = RandomReal[{0, 1}, {lth, lth}];
  First@AbsoluteTiming[Table[mtxCopy[\[Alpha], mtx, lth];, {2000}];]/
    2000/lth]

matrixScaling3[lth_] := 
 Module[{mtx, \[Alpha] = 1.}, mtx = RandomReal[{0, 1}, {lth, lth}];
  First@AbsoluteTiming[
      Table[mtxShare[\[Alpha], mtx, lth];, {2000}];]/2000/lth]

nls = Range[100, 1000, 100];

timeScaleMMA = matrixScaling /@ nls;
timeScaleMMA2 = matrixScaling2 /@ nls;
timeScaleMMA3 = matrixScaling3 /@ nls;

BarChart[1.*^6 Transpose@{timeScaleMMA, timeScaleMMA2, timeScaleMMA3},
  ChartLabels -> {Range[100, 1000, 100], None}, 
 PlotTheme -> "Detailed", 
 FrameLabel -> {"matrix size", "time (\[Mu]s)"}, 
 ChartLegends -> {"Mathematica", "LibraryLink-Copy", 
   "LibraryLink-Share"}, ImageSize -> Medium]
$\endgroup$
3
  • $\begingroup$ Sorry I left the work for you to do, but thanks for doing it. I'm glad my idea proved useful. $\endgroup$
    – Mr.Wizard
    May 16, 2016 at 14:35
  • $\begingroup$ Nice response and discussion. Please just provide an edited version of your response or write how we can do scalar-to-matrix multiplications quickly from now on. I mean, how we can call LibraryLink->"Shared" for such multiplications as fast as possible in general environments. $\endgroup$
    – Faz
    May 16, 2016 at 16:35
  • $\begingroup$ @FazlollahSoleymani This discussion is aimed at finding the reason of the slowness. I think it's not easy for me to write an easy-to-use solution in general environments, which seems to require non-trivial design efforts. I will be very happy to see someone with better knowledge than me to come up with a canonical solution. $\endgroup$ May 16, 2016 at 21:18
4
$\begingroup$

I believe that on packed arrays both Dot and Times are performed by external libraries, e.g. Intel MKL, and that following Mathematica's paradigm of immutability the library does not act directly upon the original array but rather a copy. I conjecture that this copying or transport is the cause of the slow-down that you observe and that within Mathematica Dot is faster primarily because there is less to copy.

At the moment I cannot think of a way to test this conjecture; can anyone else?

$\endgroup$
2
  • 1
    $\begingroup$ Rectangular matrices? For me, with mtx = RandomReal[{0, 1}, {100000, 1}]; ls = RandomReal[{0, 1}, {1}];, Times is faster than Dot. Varying the number of rows and columns results in a various results, although I have not systematically explored the behavior. $\endgroup$
    – Michael E2
    May 15, 2016 at 12:38
  • 1
    $\begingroup$ It seems that your conjecture is correct, we can test it using the two ways of passing arguments in LibraryLink. See my answer. $\endgroup$ May 16, 2016 at 14:29
3
$\begingroup$

So why the matrix number multiplication and addition so much slower than the matrix vector multiplication? Are there ways to speed them up?

It seems that if we use SparseArray for this computation we can get ~3 times speed-up:

n = 2000;
tres =
 Table[
  (Print[lth];
   mtx = RandomReal[{0, 1}, {lth, lth}];
   {lth, (AbsoluteTiming[
         Table[Etc[t]*mtx;, {t, 0., 20, 20/(n - 1)}]][[1]] - 
       AbsoluteTiming[Table[{};, {t, 0., 20, 20/(n - 1)}]][[1]])/
     n, (AbsoluteTiming[
         Table[(sp = SparseArray[{{_, _} -> Etc[t]}, Dimensions[mtx]];
            sp*mtx;), {t, 0., 20, 20/(n - 1)}]][[1]] - 
       AbsoluteTiming[Table[{};, {t, 0., 20, 20/(n - 1)}]][[1]])/
     n}), {lth, Range[100, 1000, 100]}]

Quartiles[tres[[All, 3]]/tres[[All, 2]]]

(* {0.347113, 0.443113, 0.905894} *)

ListPlot[Transpose[{tres[[All, 1]], tres[[All, 3]]/tres[[All, 2]]}], 
 PlotRange -> All, PlotTheme -> "Detailed"]

enter image description here

Here is the code I used to do timing exploration:

Needs["Developer`"]
{lth, n} = {400, 2000};
mtx = RandomReal[{0, 1}, {lth, lth}];
spmtx = SparseArray[mtx];
pmtx = Developer`ToPackedArray[mtx];
(AbsoluteTiming[
     sp = SparseArray[{{_, _} -> 4.}, Dimensions[mtx]]; 
     Do[sp*mtx;, {n}]][[1]]
   - AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[
     sp = SparseArray[{{_, _} -> 4.}, Dimensions[mtx]]; 
     Do[sp*pmtx;, {n}]][[1]]
   - AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[Do[4.*mtx;, {n}]][[1]] - 
   AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[Do[mtx*4.;, {n}]][[1]] - 
   AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[Do[4.*spmtx;, {n}]][[1]] - 
   AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[Do[4.*pmtx;, {n}]][[1]] - 
   AbsoluteTiming[Do[{}, {n}]][[1]])/n
(AbsoluteTiming[Table[Etc[t]*mtx;, {t, 0., 20, 20/(n - 1)}]][[1]] - 
   AbsoluteTiming[Table[{};, {t, 0., 20, 20/(n - 1)}]][[1]])/n
(AbsoluteTiming[
     Table[sp = SparseArray[{{_, _} -> Etc[t]}, Dimensions[mtx]]; 
      sp*mtx;, {t, 0., 20, 20/(n - 1)}]][[1]] - 
   AbsoluteTiming[Table[{};, {t, 0., 20, 20/(n - 1)}]][[1]])/n
$\endgroup$
14
  • 1
    $\begingroup$ Yes, I'm interested in matrix of size range from 200 to 2000. We are working on high harmonic generation from atoms, which requires large-scale simulations sometimes. $\endgroup$ May 15, 2016 at 1:35
  • 1
    $\begingroup$ Are you aware of subroutines in BLAS or MKL that can perform matrix number multiplication? I only found xSCAL that does the vector number multiplication, but there are no counterparts for matrix. $\endgroup$ May 15, 2016 at 1:56
  • 1
    $\begingroup$ @xslittlegrass BLAS is based on FORTRAN. Matrices in FORTRAN can be seen as vectors. (Unfolding is column-wise). So axpy is the one to apply. $\endgroup$ May 15, 2016 at 8:40
  • 1
    $\begingroup$ @TheDoctor Here is one of the classical reference. $\endgroup$ May 30, 2016 at 15:52
  • 1
    $\begingroup$ @xslittlegrass: this reference mentions that the time-dependent Schrodinger equation is integrated directly on a numerical grid. I would have thought that one could use Chebyshev-polynomial expansion and the fast-Fourier-transform algorithm, modified for atoms exposed to an intense laser field. See, e.g., Electronic wave propagation with Mathematica (DOI: dx.doi.org/10.1016/S0010-4655(00)00196-X) $\endgroup$
    – TheDoctor
    Jun 10, 2016 at 7:56
1
$\begingroup$

Not a full answer, but there are other tools available, such as Outer and TensorProduct. For example, compare

AbsoluteTiming[(mat1 = Table[Etc[t] mtx, {t, 0., 20, 0.01}]);]

with

AbsoluteTiming[(mat2 = TensorProduct[Et[Range[0., 20, 0.01]], mtx]);]

TensorProduct is over twice as fast. The answers are identical:

Norm[Flatten[mat1 - mat2]] // Chop
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.