# Speedup matrix number multiplication

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

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

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

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]

• 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. – bill s May 14 '16 at 13:33
• @bills it doesn't seem to come from calling Etc, see update. – xslittlegrass May 14 '16 at 15:03
• The matrix number multiplication is 2.2 times slower using version 9.0.1 than it is using version 10.4.1. – Karsten 7. May 15 '16 at 1:07
• 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. – Sander May 15 '16 at 6:31
• @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. – xslittlegrass May 15 '16 at 6:49

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.

## 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


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

mtxShare =
"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)"},

• Sorry I left the work for you to do, but thanks for doing it. I'm glad my idea proved useful. – Mr.Wizard May 16 '16 at 14:35
• 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. – Fazlollah May 16 '16 at 16:35
• @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. – xslittlegrass May 16 '16 at 21:18

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?

• 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. – Michael E2 May 15 '16 at 12:38
• It seems that your conjecture is correct, we can test it using the two ways of passing arguments in LibraryLink. See my answer. – xslittlegrass May 16 '16 at 14:29

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"]


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 = DeveloperToPackedArray[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

• 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. – xslittlegrass May 15 '16 at 1:35
• 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. – xslittlegrass May 15 '16 at 1:56
• @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. – Anton Antonov May 15 '16 at 8:40
• @TheDoctor Here is one of the classical reference. – xslittlegrass May 30 '16 at 15:52
• @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) – TheDoctor Jun 10 '16 at 7:56

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
`