16
$\begingroup$

I have some data generated from some program, and it appears that matrix multiplication on these data are about 10 times slower than on some random data:

Get["tb.dat"];

xls = Range[-500, 500, 1000/(1000 - 1)] // N;

Re[Conjugate[#].(xls*#)] & /@ tb; // AbsoluteTiming
(* {0.067147, Null} *)

tb2 = RandomComplex[{0., 1. + I}, Dimensions[tb]];

Re[Conjugate[#].(xls*#)] & /@ tb2; // AbsoluteTiming
(* {0.004564, Null} *)

So what are the reasons for the performance problem with my data?

The data files are here (6MB).

Update

Here is the same problem with the packed array. Since Save unpacks data, so I have to use DumpSave to demonstrate the problem:

DumpGet["tb3PK.dump"];

Developer`PackedArrayQ@tb3PK
(* True *)

xls = Range[-500, 500, 1000/(1000 - 1)] // N;

Re[Conjugate[#].(xls*#)] & /@ tb3PK; // AbsoluteTiming
(* {0.065747, Null} *)

tb2 = RandomComplex[{0., 1. + I}, Dimensions[tb3PK]];

Re[Conjugate[#].(xls*#)] & /@ tb2; // AbsoluteTiming
(* {0.003482, Null} *)

The file tb3PK.dump is here . I'm using OS X 10 and Mathematica version 10.0.1

Improve this question
$\endgroup$
12
  • 7
    $\begingroup$ It is because your data is not packed, but data generated by RandomComplex is. Try this: {{Developer`PackedArrayQ@tb, ByteCount@tb},{Developer`PackedArrayQ@tb2, ByteCount@tb2}}. However according to what I have learned about packed arrays, your data should be packable because Equal @@ Map[Head, tb, {2}] is true, i.e. all of your elements are of the same type. However, Developer`PackedArrayQ@Developer`ToPackedArray[tb] returns False. I am not writing this as an answer as I hope that an answer will address how to pack your data. $\endgroup$
    – C. E. ♦
    Commented Jan 9, 2015 at 4:35
  • $\begingroup$ @Pickett No wonder tb3=Developer``ToPackedArray[tb] is as slow as the unpacked data. I didn't know that Developer`ToPackedArray could fail to pack an array without giving any error message. $\endgroup$
    – xslittlegrass
    Commented Jan 9, 2015 at 5:04
  • 1
    $\begingroup$ @MichaelE2 but why the problem still exists even when I set SetSystemOptions["CatchMachineUnderflow" -> False], it seems to me if the problem is caused by the auto convention to high precision, turn off it would resolve the problem ? $\endgroup$
    – xslittlegrass
    Commented Jan 9, 2015 at 19:55
  • 1
    $\begingroup$ Well, I'm not sure. Setting it to False speeds it up, but only by a factor of 4. Perhaps there's a cost to underflow even when it is not "caught." Someone else is going to have to answer this. Sorry. $\endgroup$
    – Michael E2
    Commented Jan 9, 2015 at 20:42
  • 1
    $\begingroup$ Get on its own should work. DumpGet has no doc page, though there exists an error page. $\endgroup$
    – Sjoerd C. de Vries
    Commented Jan 9, 2015 at 22:54

2 Answers 2

Reset to default
13
+50
$\begingroup$

For readers who didn't read all the comments, the slowdown is due to a lack of packing of tb, whereas RandomReal returns packed arrays when more than 250 elements are generated. The reason why packing tb fails is because some elements have different precision than others, and (I think?) ToPackedArray requires arrays to be of homogeneous type.

To fix this, chop and multiply by 1. + 0. I, then pack:

<< Developer`
tbP = ToPackedArray[(1. + 0. I) Chop@tb];
PackedArrayQ[tb2]
Re[Conjugate[#].(xls*#)] & /@ tb; // AbsoluteTiming
Re[Conjugate[#].(xls*#)] & /@ tbP; // AbsoluteTiming

producing

True
{0.031200, Null}
{0., Null}

Does anyone know a better way of taking an array composed of heterogeneous-precision real and complex numbers, and converting it into an array of homogeneous-precision complex numbers? The above trick works, but I feel like there ought to be a better way.

Update

To address the Update in your question, I'm stumped. Comparing the DumpSave output with the solution I used above:

DumpGet["tb3PK.dump"];
Re[Conjugate[#].(xls*#)] & /@ tbP; // AbsoluteTiming
Re[Conjugate[#].(xls*#)] & /@ tb3PK; // AbsoluteTiming
Tally[Precision /@ Flatten[tbP]]
Tally[Precision /@ Flatten[tb3PK]]
ByteCount[tbP]
ByteCount[tb3PK]

giving

{0., Null}
{0.040000, Null}
{{MachinePrecision, 100000}}
{{MachinePrecision, 100000}}
1600152
1600152

In other words, the DumpSave version is indeed slower as you observed, despite having the same bytecount and precision, and I'm not sure why. You may want to ask that as a separate question.

Improve this answer
$\endgroup$
9
  • $\begingroup$ Could use N[array]. $\endgroup$
    – Daniel Lichtblau
    Commented Jan 9, 2015 at 15:36
  • $\begingroup$ I think the numbers in the tb are all MachinePrecision, for example And @@ ((Precision[#] == MachinePrecision &) /@ tb) returns True. Also I'm interested in the log spectrum of the data, where small numbers are important to me, so I would prefer to not Chop if possible. $\endgroup$
    – xslittlegrass
    Commented Jan 9, 2015 at 15:37
  • $\begingroup$ @DanielLichtblau: Using N[tb] and then packing fails to pack the array, because the result is still heterogeneous-precision. $\endgroup$
    – DumpsterDoofus
    Commented Jan 9, 2015 at 15:47
  • $\begingroup$ @xslittlegrass: Tally[Precision /@ Flatten[tb]] shows that most, but not all elements are MachinePrecision. $\endgroup$
    – DumpsterDoofus
    Commented Jan 9, 2015 at 15:50
  • 1
    $\begingroup$ @xslittlegrass 15.954589770191005 == MachinePrecision does return True, but when using as precision, they're different, the latter represents machine precision, while the former represents arbitrary precision, they lead to 2 different algorithm… well, it's a little hard for me to explain in my own words, but there're many related posts in this site, for example this one. $\endgroup$
    – xzczd ♦
    Commented Jan 10, 2015 at 5:08
2
$\begingroup$

I propose a silly workaround, instead of a workable explanation:

Get["tb.dat"];

xls = Range[-500, 500, 1000/(1000 - 1)] // N;
test[tb_] := (Re[Conjugate[#].(xls*#)] & /@ tb;) // AbsoluteTiming;
test[RandomComplex[{0., 1. + I}, Dimensions[tb]]]
(* {0.002002, Null} *)
test[tb]
(* {0.040038, Null} *)
test[tb + ConstantArray[0. + 0. I, Dimensions@tb]]
(* {0.003003, Null} *)
test[tb (1.0 + 0. I)]
(* {0.002002, Null} *)

This tricks Mathematica into doing the right thing. I find it often (too often) helps a lot to multiply with 1., or just add 0.. From my tests it doesn't seem to matter if you multiply by (1.0 + 0. I) or add a constant array, you get the same performance boost. It also doesn't seem to matter if the array is packed or not.

Improve this answer
$\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.