Let m
be a huge matrix and say I want to calculate m.Transpose[m]
.
The result will be a symmetric matrix, and so it feels like the standard matrix multiplication might just do too much work (or does Mathematica recognize such a product and use an optimized algorithm?). The result from m.Transpose[m]
would also be a standard matrix, while a SymmetrizedArray
would probably be more suitable (at least memory-wise). So, therefore, the following
Question: What would be the best way to calculate a product of a matrix with its own transpose, both memory-wise and time-wise?
EDIT: I tried the SYRK
solution and 'Benchmarked' it (see code below). The code was run on a linux machine where Mathematica had privileged access to 4 cores (via https://unix.stackexchange.com/questions/326579/how-to-ensure-exclusive-cpu-availability-for-a-running-process)
The results are the following:
|Null
| Integer m
| Integer Packed m
| Float m
| Float Packed m|
|-|-|-|-|-|
|m.Transpose[m]
| 8.62596| 8.45888| 5.18896| 5.38446|
|SYRK Integer b| 0.055395| 0.057715| 0.055435| 0.055741|
|SYRK Integer Packed b| 0.200872| 0.194433| 0.195106| 0.195942|
|SYRK Float b| 0.055659| 0.055549| 0.055976| 0.05589|
|SYRK Float Packed b| 0.196128| 0.193704| 0.192634| 0.195042|
Conclusions:
SYRK
is indeed much faster but one should take care not to pack the array to which the result is written.- It seems
SYRK
doesn't care whetherm
is packed - It seems
SYRK
doesn't care whether any of the arrays has integers or floats
These conclusions are in line with the answer by user293787.
(* For PackedArray and BLAS functionality *)
<<Developer`
<<LinearAlgebra`BLAS`
(* Matrix Dimensions *)
{ d1, d2 } = { 40000, 10 };
(* Number of random matrices for which to test *)
nLoops = 10;
(* Generate tables of timings *)
tables =
Reap[ # ][[2,1]]& @
Do[ (* loop over Seeds and Sow tables of timings *)
SeedRandom[s];
(* Define the matrices for multiplication *)
mPacked = RandomInteger[ { -100, 100 }, { d1, d2 } ];
m = FromPackedArray @ mPacked;
mFloat = FromPackedArray @ N @ m;
mFloatPacked = ToPackedArray @ mFloat;
testMats = { m, mPacked, mFloat, mFloatPacked };
(* Define the matrices to be filled. Use Hold to prevent memory overflow *)
integerPackedB = Hold @ ConstantArray[ 0, { d1, d1 } ];
integerB = Hold @ FromPackedArray @ ConstantArray[ 0, { d1, d1 } ];
floatPackedB = Hold @ ConstantArray[ 0., { d1, d1 } ];
floatB = Hold @ FromPackedArray[ConstantArray[ 0., { d1, d1 } ] ];
testBMats = Hold @ { integerB, integerPackedB, floatB, floatPackedB };
Do[ (* timings for Mathematica matrix multiplication *)
ClearSystemCache[];
t[1,i] = First @ AbsoluteTiming[ testMats[[i]].Transpose[ testMats[[i]] ]; ],
{i,4}
];
Do[ (* timings for SYRK matrix multiplication *)
Do[
ClearSystemCache[];
b = ReleaseHold[ testBMats[[1,i]] ];
t[j+1,i] = First @ AbsoluteTiming[ SYRK[ "U", "N", 1,testMats[[j]], 0, b ] ];
b=.;,
{i,4}
],
{j,4}
];
Sow[ Array[ t, {5,4} ] ];
ClearAll[t];
,{ s, nLoops } ];
TableForm[
Plus @@ tables / nLoops,
TableHeadings ->
{
{ "m.Transpose[m]", "SYRK Integer b", "SYRK Integer Packed b", "SYRK Float b", "SYRK Float Packed b" },
{"Integer m", "Integer Packed m","Float m", "Float Packed m"}
}
]
```
m
an 37000x8 matrix compared to the regularm.Transpose[m]
. (where m is a PackedArray) $\endgroup$m
is not packed, but sayRandomInteger[{-5,5},{10,10}]//Developer`PackedArrayQ
yields true. I use Version 12.3. $\endgroup$LinearAlgebra`BLAS`
but forgot to add it to the code. This is now fixed. Thank you so much for your help and patience. $\endgroup$