σ0 = PauliMatrix[0]; (*{{1, 0}, {0, 1}}*)
σ1 = PauliMatrix[1]; (*{{0, 1}, {1, 0}}*)
σ2 = PauliMatrix[2]; (*{{0, -I}, {I, 0}}*)
σ3 = PauliMatrix[3]; (*{{1, 0}, {0, -1}}*)

O2 = {{0, 0}, {0, 0}};

γ0 = ArrayFlatten[{{σ0, O2}, {O2, σ0}}];
γ1 = ArrayFlatten[{{O2, σ1}, {-σ1, O2}}];
γ2 = ArrayFlatten[{{O2, σ2}, {-σ2, O2}}];
γ3 = ArrayFlatten[{{O2, σ3}, {-σ3, O2}}];

MatrixForm[{{O2, σ1}, {-σ1, O2}}]

sigma[M__] := Module[{submat, basis, fb, nrm, cf},
               submat = Partition[M, {2, 2}, 2]; (*2x2 submatrices*)
               basis = {σ0, σ1, σ2, σ3};
               fb = Flatten[basis, {{1}, {2, 3}}];
               nrm = Diagonal[fb.ConjugateTranspose[fb]]; (*norm of basis*)
               cf[m_] := Flatten[m].ConjugateTranspose[fb]/nrm; (*coeff for σ*)
               Map[{"I", "σ1", "σ2", "σ3"}.cf[#] &, submat, {2}]]

sigma[γ2] // MatrixForm

sigma[γ1.γ2] // MatrixForm

M = ArrayFlatten@Outer[Times, σ2, σ3];
MatrixForm[M]
sigma[M] // MatrixForm