# Efficient method (maybe parallel method) to calculate a 256$\times$256 matrix?

I want to calculate a matrix of shape 256$$\times$$256 as follows

ParallelTable[Tr[emat[[i]].IntegrateUnitaryHaar[UV.emat[[j]].ConjugateTranspose[UV],{u,15},{v,1}]]/16,{i,1,256},{j,1,256}]


where emat[[i]] is a matrix. I have already used ParallelTable to execute the code parallelly. But it seems the program is still very inefficient. I wonder if there are some more efficient methods to execute the code?

Edit Here is some more information about the code that I want to speed up. I create an emat to help me construct the target matrix. emat is defined as

emat=ConstantArray[0,{256,16,16}];
(*Calculate Bases*)
Do[emat[[(i-1)*64+(j-1)*16+(k-1)*4+l]]=KroneckerProduct[KroneckerProduct[KroneckerProduct[Subscript[\[Sigma], i],Subscript[\[Sigma], j]],Subscript[\[Sigma], k]],Subscript[\[Sigma], l]],{i,1,4},{j,1,4},{k,1,4},{l,1,4}];


where Sigma are pauli matrices defined as

Subscript[\[Sigma], 1] = ( {
{0, 1},
{1, 0}
} ); Subscript[\[Sigma], 2] = ( {
{0, -I},
{I, 0}
} ); Subscript[\[Sigma], 3] = ( {
{1, 0},
{0, -1}
} ); Subscript[\[Sigma], 4] = ( {
{1, 0},
{0, 1}
} );


Before creating the target matrix, I need some more materials I dubbed U and V matrices. Code as follows

U = Array[Subscript[u, #1, #2] &, {15, 15}];
V = Array[Subscript[v, #1, #2] &, {1, 1}];
UV = ArrayFlatten[{{V, 0}, {0, U}}];
changebasisu = IdentityMatrix[hilbertdim];
changebasisu[[1, 1]] = 1/Sqrt[2]; changebasisu[[1, 16]] = 1/Sqrt[2];
changebasisu[[16, 1]] = 1/Sqrt[2]; changebasisu[[16, 16]] = -1/Sqrt[2];
UV = changebasisu . UV . ConjugateTranspose[changebasisu];


Finally I can write the formula for the elements of the target matrix to be

ParallelTable[Tr[emat[[i]].IntegrateUnitaryHaar[UV.emat[[j]].ConjugateTranspose[UV],{u,15},{v,1}]]/16,{i,1,256},{j,1,256}]


where IntegrateUnitaryHaar is a function to calculate integration with respect to the Haar measure on the unitary groups. To execute this function, you need to install a related package inside Mathematica and using code << IntU inside Mathematica.

The execution result of $Version on my Mathematica is 13.1.0 for Microsoft Windows (64-bit) (June 16, 2022). I've tried to execute the following code to calculate only the first line of the target matrix, with AbsoluteTiming function as ParallelDo[pmat[[1,j]]=Tr[emat[[1]].IntegrateUnitaryHaar[UV.emat[[j]].ConjugateTranspose[UV],{u,15},{v,1}]]/16,{j,1,256}];//AbsoluteTiming  The result is {884.955,Null} which is very slow! Since I think every element in essence can be calculated in an absolutely uncorrelated way. So I really expect parallelly execute the code to dramatically slow down the running time. But I didn't solve it yet, is it I just didn't do it right? The total code showing as follows(the code needs to download the package and then install it inside Mathematica. I temporarily don't know how to add it directly to the code below.) Clear["Global*"]; Subscript[\[Sigma], 1] = ( { {0, 1}, {1, 0} } ); Subscript[\[Sigma], 2] = ( { {0, -I}, {I, 0} } ); Subscript[\[Sigma], 3] = ( { {1, 0}, {0, -1} } ); Subscript[\[Sigma], 4] = ( { {1, 0}, {0, 1} } ); emat = ConstantArray[0, {256, 16, 16}]; (*Calculate Bases*) Do[emat[[(i - 1)*64 + (j - 1)*16 + (k - 1)*4 + l]] = KroneckerProduct[ KroneckerProduct[ KroneckerProduct[Subscript[\[Sigma], i], Subscript[\[Sigma], j]], Subscript[\[Sigma], k]], Subscript[\[Sigma], l]], {i, 1, 4}, {j, 1, 4}, {k, 1, 4}, {l, 1, 4}]; (*target matrix*) pmat = ConstantArray[0, {256, 256}]; (*Prepare package for doing haar integral*) << IntU; (*Routine to do haar integral*) U = Array[Subscript[u, #1, #2] &, {15, 15}]; V = Array[Subscript[v, #1, #2] &, {1, 1}]; UV = ArrayFlatten[{{V, 0}, {0, U}}]; changebasisu = IdentityMatrix[16]; changebasisu[[1, 1]] = 1/Sqrt[2]; changebasisu[[1, 16]] = 1/Sqrt[2]; changebasisu[[16, 1]] = 1/Sqrt[2]; changebasisu[[16, 16]] = -1/Sqrt[2]; UV = changebasisu . UV . ConjugateTranspose[changebasisu]; pmat = ParallelTable[ pmat[[i, j]] = Tr[emat[[i]] . IntegrateUnitaryHaar[ UV . emat[[j]] . ConjugateTranspose[UV], {u, 15}, {v, 1}]]/ 16, {i, 1, 256}, {j, 1, 256}]; // AbsoluteTiming pmat  • To make the before/after comparison more effective you can include timing information from AbsoluteTiming. Also, include the output of $Version to your post. What type of entries does emat have? Please include emat and UV. In short, please present a minimal example and include definitions.
– Syed
Mar 31, 2023 at 6:11
• @Syed Thanks for pointing out those. I edited my post. Mar 31, 2023 at 6:38
• You are recomputing IntegrateUnitaryHaar[UV.emat[[j]].ConjugateTranspose[UV],{u,15},{v,1}] for every value of i instead of memoizing it. Mar 31, 2023 at 8:07
• Use sparse Pauli matrices for speed: Clear[σ, emat]; σ[i_] := σ[i] = SparseArray[PauliMatrix[i]]; emat[i_Integer, j_Integer, k_Integer, l_Integer] := emat[i, j, k, l] = KroneckerProduct[σ[i], σ[j], σ[k], σ[l]] will give you memoized, sparse, quick access to the emat values without forcing a pre-compute. Mar 31, 2023 at 8:08
• Can you please edit your post so that the code can be copy-pasted into a notebook and will actually run? Please don't expect people to piece it together by themselves. Mar 31, 2023 at 8:19

I tried your first code and found out that your emat tensor is nothing else but Flatten[TensorProduct[#,#,#]&[PauliMatrix[{1,2,3}], 3]

The execution observed via task manager in Windows is showing an excessive use of the virtual memory on the disk.

I have a SSD as system disk and a backup HDD. The HDD is used and the prediction of Timing shown by Parallel amounts to hours.

The simple reason seems to be the unbeleavable amount of symbolic expression trees written and read from and to disk for any tiny step.

If RAM memory holds the complete execution tree there is a time factor of 0.1 with respect to a SSD depending on the interface.

I finally finish the calculation in finite time. The method is based on the comments of @Roman. Codes showing as follows

Clear["Global*"];
\[Sigma][i_] := \[Sigma][i] = SparseArray[PauliMatrix[i]];
emat = SparseArray[ConstantArray[0, {256, 16, 16}]];
Do[emat[[(i - 1)*64 + (j - 1)*16 + (k - 1)*4 + l]] =
KroneckerProduct[\[Sigma][i], \[Sigma][j], \[Sigma][k], \[Sigma][
l]], {i, 1, 4}, {j, 1, 4}, {k, 1, 4}, {l, 1, 4}];
(*Prepare package for doing haar integral*)
<< IntU;
(*Routine to do haar integral*)
U = Array[Subscript[u, #1, #2] &, {15, 15}];
V = Array[Subscript[v, #1, #2] &, {1, 1}];
UV = ArrayFlatten[{{V, 0}, {0, U}}];
changebasisu = IdentityMatrix[16];
changebasisu[[1, 1]] = 1/Sqrt[2]; changebasisu[[1, 16]] = 1/Sqrt[2];
changebasisu[[16, 1]] = 1/Sqrt[2]; changebasisu[[16, 16]] = -1/Sqrt[2];
UV = changebasisu . UV . ConjugateTranspose[changebasisu];
integralemat[j_Integer] :=
integralemat[j] =
IntegrateUnitaryHaar[
UV . emat[[j]] . ConjugateTranspose[UV], {u, 15}, {v, 1}];
pmat = SparseArray[ConstantArray[0, {256, 256}]];
Do[pmat[[i, j]] = Tr[emat[[i]] . integralemat[j]]/16, {i, 1, 256}, {j,
1, 256}]; // AbsoluteTiming


The result of AbsoluteTiming function is {890.54, Null}`. But this time it's the time for the total target matrix instead of one row.

• I tried your first code and found out that your emat is nothing else but Flatten[TensorProduct[#,#,#]&[PauliMatrix[{1,2,3}], 3] Mar 31, 2023 at 16:48