Speed up large matrix multiplications with or without high precision?

Questions
,
1

I have several matrices (2D matrices) sizes can go up to 30000*30000 and I need to do matrix multiplications. However, even if I use 16 Digits the calculation take extremely long compare to MATLAB. I am not expecting Mathematica to be faster but not 100 times slower.

So for MATLAB I create

A=rand(300,300);
tic
A*A;
toc

and it takes only 0.004 seconds with format long ( double precision)

For Mathematica

a = RandomReal[{-10, 10}, {300, 300}, WorkingPrecision -> 15];
AbsoluteTiming[f = a.a;]

taking 13 seconds and there is no time difference if I increase the WorkingPrecision to 30. But more interestingly if I put the working precision to 6 digits (single precision) time consumed is same. An example.

a = RandomReal[{-10, 10}, {300, 300}, WorkingPrecision -> 6];
AbsoluteTiming[f = N[a].N[a];]
AbsoluteTiming[g = a.a;]
Max[Abs[f - g]]

The difference between g and f far being zero. I don’t understand why I get different results in here.

I can’t imagine doing such a calculation with really large matrices. I am certain that something is not right in my Mathematica approach. What is the way to speed up matrix multiplication in Mathematica.

2

In[3]:= a = RandomReal[{-10, 10}, {300, 300}];
AbsoluteTiming[f = a.a;]

Out[4]= {0.135985, Null}

In[5]:= a = RandomReal[{-10, 10}, {300, 300}];
AbsoluteTiming[f = a.a;]

Out[6]= {0.00157647, Null}

In[7]:= a = RandomReal[{-10, 10}, {300, 300}];
AbsoluteTiming[f = a.a;]

Out[8]= {0.00154512, Null}

I suspect that some functions run slower the first time they are used.

3

that is an interesting observation, however, what you did corresponds to

a = RandomReal[{-10, 10}, {300, 300}];
AbsoluteTiming[f = N[a].N[a];]

just working with 6 digits not even double precision.

4

In[15]:= a = RandomReal[{-10., 10.}, {300, 300}];
a.a; // AbsoluteTiming
$Version

Out[16]= {0.0046302, Null}

Out[17]= "10.4.1 for Microsoft Windows (64-bit) (April 11, 2016)"
5

If you want to use more digits than 6 as in

a = RandomReal[{-10., 10.}, {300, 300},WorkingPrecision → 16]

again time is 13 sec. That is the part is confusing me. Why even double precision (which is physical) does take that long.

6

I think that WolkingPrecision → 16 means that 16 digits are exact, i.e. Mma is working with more digits than 16 (more than two words) to guarantee specified exactitude. When working with double precision, only double words are used. It’s only a supposition, not a certainty…

Apropos,

In[4]:= a = RandomReal[{-10., 10.}, {10000, 10000}];
a.a; // AbsoluteTiming
$Version

Out[5]= {91.8452, Null}

Out[6]= "10.4.1 for Microsoft Windows (64-bit) (April 11, 2016)"
7

I think being exact is not an easy concept. My problem, I want to work with double precision (at least) and when I am working with double precision I expect a similar performance from Mathematica that of MATLAB for simple (highly parallelizable) calculation like matrix multiplication.

a = RandomReal[{-10, 10}, {10000, 10000}];
AbsoluteTiming[f = a.a;]

this takes around 17 seconds. It is 6 digits, same as your example. With MATLAB double precision calculation takes 17 sec.
10.3.1 for Mac OS X x86 (64-bit) (December 9, 2015)

8

What’s wrong with Machine Precision? It’s the default and close to 16 digits.

In[1]:= $MachinePrecision

Out[1]= 15.9546
9

You are right. In Mathematica we deal with:

Precision[x]
which gives the effective number of digits of precision in the number x,

and

Accuracy[x]
which gives the effective number of digits to the right of the decimal point in the number x.

There is also:

MachinePrecision - a symbol used to indicate machine-number precision.

The last is somewhat appropriate to traditional interpretation.

10

When you use:

a = RandomReal[{-10, 10}, {300, 300}, WorkingPrecision -> n];

for whatever n, the numbers that are generated are so-called arbitrary precision numbers. When you don’t use a WorkingPrecision option:

a = RandomReal[{-10, 10}, {300, 300}];

it will create machineprecision numbers:

MachineNumberQ[a[[1,1]]]

Evaluates to True. While in the former case (independent of n) it is False.

Machine precision numbers do not track the precisions of numbers, while arbitrary precision does. This makes a HUGE difference! On most modern machine a MachinePrecision number will be a ‘double precision floating point’:

$MachinePrecision
15.9546  (53 Log10[2])

By default this will be used. Note that Mathematica only SHOWS 6 digits by default on screen! i.e. RandomReal[] gives 0.531239 (but in reality it has all 53 bits set randomly)…

11

FYI:

a = RandomReal[{-10, 10}, {5000, 5000}];
AbsoluteTiming[a.a;]

OR

A=rand(5000,5000);tic;A*A;toc

gives 3 seconds on both of my machines, just milliseconds difference… Presumably because both use BLAS (afaik).

12

Well, there is nothing wrong with machine precision but also I know several problems more is required. Also one migh ask, why Wolfram developed extended precision ? I am trying find the best use of the tool that I have.

13

Thank you for the reply. I was getting similar times for the example that you also showed but I did not know that Mathematica was only showing 6 digits instead of machine precision. How can I make displaying full digits?

From your response I can see that it does not matter 6 for 60 digits with WorkingPrecision there is a cost of keeping track and that is quite huge.

14

Arbitrary precision is sometimes useful for certain applications. In most cases however double precision is more than enough for most applications. The main problem is that it can’t use the very efficient multiplication algorithms in the CPU, and have to resort to different algorithm with overhead. The WorkingPrecision does matter in performance:

a = RandomReal[{-10, 10}, {50, 50}, WorkingPrecision -> 6];
AbsoluteTiming[a.a;]
a = RandomReal[{-10, 10}, {50, 50}, WorkingPrecision -> 60];
AbsoluteTiming[a.a;]
a = RandomReal[{-10, 10}, {50, 50}, WorkingPrecision -> 600];
AbsoluteTiming[a.a;]
a = RandomReal[{-10, 10}, {50, 50}, WorkingPrecision -> 6000];
AbsoluteTiming[a.a;]

But only for really large amounts of digits…

You can use e.g.

NumberForm[%, 20]

to see all the digits. Or change the default in the preferences:

ScreenShot2016-07-12at20.04.43.png
ScreenShot2016-07-12at20.04.43.png1636×1612 138 KB

In the slightly highlighted field…

15

Basically by default it uses MachinePrecision, good for 16 digits of precision (on most modern computers). If you want more, we can, using what Mathematica calls ‘arbitrary precision’. This however prohibits very efficient CPU function calls that work very nicely with doubles…