# Efficient computation of a matrix times its transpose

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 whether m 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
<<LinearAlgebraBLAS

(* 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,
{
{ "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"}
}
]
$$$$

• Quick search gave this similar question for BLAS on stackoverflow. Have you tried calling BLAS from Mathematica? Perhaps this goes in the right direction Jul 17, 2022 at 17:06
• Are you going to use this symmetric product matrix for further calculations? If yes, then there might be even faster methods that avoid calculating this intermediate altogether. For example, if you're doing linear least squares or calculating a Moore–Penrose pseudoinverse, then you don't need the symmetric intermediate matrix and can go straight to the end result. Jul 17, 2022 at 17:21
• @Roman Yes, but mainly just to compare values of certain blocks. @user293787 This is very interesting indeed. Strangely it seems that SYRK takes about 4 times as long to calculate the multiplication with m an 37000x8 matrix compared to the regular m.Transpose[m]. (where m is a PackedArray)
– Gert
Jul 17, 2022 at 19:38
• In your edit, you seem to assume that m is not packed, but say RandomInteger[{-5,5},{10,10}]//DeveloperPackedArrayQ yields true. I use Version 12.3. Jul 19, 2022 at 12:41
• @user293787: I changed the code so that the 'b-matrices' are properly packed/unpacked. The results are now in line with your findings :) I did load LinearAlgebraBLAS but forgot to add it to the code. This is now fixed. Thank you so much for your help and patience.
– Gert
Jul 20, 2022 at 10:06

Edit and Warning: I noticed that my timing may not have been quite correct, or at least has to be interpreted carefully, please see the edit.

As @Roman points out, one should check if one really has to perform this calculation, but I will not question that here. OP mentions size $$37000 \times 8$$ in the comments. I will use size $$10000 \times 10$$ for illustration.

BLAS call. Motivated by this stackoverflow post let me use SYRK from the BLAS library, and see here for other BLAS commands in Mathematica:

(* sample data *)
n=10000; k=10;
m=RandomReal[{-1,1},{n,k}];

(* container for the result, must be a symbol *)
out=ConstantArray[N[0],{n,n}];

(* calculate and store upper triangular part of m.Transpose[m] in out *)
LinearAlgebraBLASSYRK["U","N",1,m,0,out];


Check. Beware that the code above only calculates the upper triangular part of m.Transpose[m]. So let us check that the upper triangular part of out was calculated correctly:

(* calculate in usual way *)
out2=m.Transpose[m];

(* this gives zero *)
Chop[Norm[Flatten[UpperTriangularize[out-out2]]]]


Timing.

RepeatedTiming[LinearAlgebraBLASSYRK["U","N",1,m,0,out];]

RepeatedTiming[out2=m.Transpose[m];]


So there does seem to be a gain.

Timing (Edit). If I reset the out container to zero each time, but still a packed array, then the SYRK calls take longer:

Table[
out=ConstantArray[N[0],{n,n}]; (* packed! *)
First[AbsoluteTiming[LinearAlgebraBLASSYRK["U","N",1,m,0,out];]]
,{10}]//Mean

Note that the timing is only for the SYRK call. I do not know why it takes longer. I can only guess that a BLAS call somehow puts out in a more BLAS-optimized state, better than a mere packed array. Perhaps someone can clarify. In any case, this may not be a problem if one simply does not reset out.
Memory. If m is square, which I understand it may not be in OP's case, and if one can live with destroying m, then I believe one can use m as the output container also, assuming m is square, as in LinearAlgebraBLASSYRK["U","N",1,m,0,m]`. Whether there is some way to directly get a structured array of some sort, I do not know.