weirdPauliMatricesDecomposition.m

generatePair[pauliIndices_List, tuple_List] := {
   BitXor[
    ReplaceAll[
     pauliIndices, {3 -> 0, 2 -> 1}
     ],
    tuple
    ],
   (-1)^Total@tuple[[Flatten@Position[pauliIndices, 2 | 3]]] 
   };

matrixElementsFromPauliTuple[matrix_, pauliIndices_] := Table[
     {
        #[[2]],
        {FromDigits[tuple, 2] + 1, FromDigits[#[[1]], 2] + 1}
        } &@generatePair[pauliIndices, tuple],
     {tuple, Tuples[{0, 1}, Length@pauliIndices]}
     ] // Transpose // Total@Times[
      matrix[[Sequence @@ #]] & /@ #[[2]],
      #[[1]]
      ] &;

coefficientsFromMatrix[matrix_] := Table[
    matrixElementsFromPauliTuple[matrix,
     IntegerDigits[idx - 1, 4, Log[2, Length@matrix]]
     ],
    {idx, 4^Log[2, Length@matrix]}
    ]/Length@matrix;

matrixElementsFromPauliTuple[
 IdentityMatrix@4,
 {0, 0}
 ]

coefficientsFromMatrix[IdentityMatrix@4]
coefficientsFromMatrix[IdentityMatrix@8] // Length

alternative "standard" method:

numQubits = 7;
mat = RandomReal[{-1, 1}, 2^numQubits {1, 1}];
(Tr@Dot[
      KroneckerProduct @@ (SparseArray@*PauliMatrix /@ IntegerDigits[#, 4, numQubits]),
      mat
      ]) & /@ Range[0, 4^numQubits - 1] // Dimensions

Roughly 4s to run with 7 qubits. Roughly 35s for 8 qubits.

direct method with LinearSolve:

{sparsePaulis[0], sparsePaulis[1], sparsePaulis[2], sparsePaulis[3]} = SparseArray /@ PauliMatrix /@ {0, 1, 2, 3};
integerToPauliTuple[integer_, numQubits_] := IntegerDigits[integer - 1, 4, numQubits];
matrixOfPaulis[numQubits_] := matrixOfPaulis[numQubits] = Transpose@SparseArray@Table[
      KroneckerProduct @@ (sparsePaulis /@ integerToPauliTuple[idx, numQubits]) // Flatten,
      {idx, 4^numQubits}
];
randomMat[numQubits_] := RandomReal[{-1, 1}, 2^numQubits {1, 1}];

LinearSolve[matrixOfPaulis@8, Flatten@randomMat@8]

This manages to solve up to 8 qubits, in roughly 8 seconds. Hangs for more than a minute for 9 qubits.