Decompose Commutators
05AB1E, 36 29 bytes
ÄàD(Ÿsgãsgãã.ΔD2FsN._`øδ*O}-Q

-7 bytes thanks to @ovs in another answer of mine.
Extremely slow since I'm using a brute-force approach with three cartesian products. But performance is irrelevant for code-golf challenges I guess (would love to see this same challenge as fastest-code or fastest-algorithm).
Try it online (only works for matrices of size 2x2, and even then usually times out if the absolute maximum of the integers in the matrix \$M\$ is too large..)
Verify whether the output \$[A,B]\$ results in the input \$M\$.
Explanation:
Assumes all the integers used in the resulting matrices \$[A,B]\$ will be within the range \$[-\max(\lvert M\rvert), \max(\lvert M\rvert)]\$. If you know any way to prove or disprove this assumption, let me know.
Note that there are of course matrices \$[A, B]\$ with integers larger than the integers of the resulting matrix \$M\$. But this \$M\$ will have other possible \$[A, B]\$ which are within this min-max range. I.e. \$[A,B]=\$ [[[1,2],[3,4]], [[-1,-2],[-3,-3]]] will have \$M=\$ [[0,2],[-3,0]], but using this \$M\$ to find a possible \$[A,B]\$ will result in [[[3,2],[3,3]], [[2,2],[3,3]]], which are all within the range \$[-3,3]\$.

Let me start with an explanation of 2FsN._`øδ*O}-, which will calculate \$M = AB-BA\$ of a given matrix-pair \$[A,B]\$ (which is the same code as in the verifier TIO-link above).
2F       # Loop 2 times:
  s      #  Swap the top two values on the stack
         #  (so the matrix-pair is at the top again after the first iteration)
   N._   #  Rotate the pair of matrices the (0-based) loop-index amount of times
         #  (to reverse the order from [A,B] to [B,A] in the second iteration)
      `  #  Push both matrices separated to the stack
  ø      #  Zip/transpose the top matrix; swapping rows/columns
   δ     #  Apply double-vectorized:
    *    #   Multiply
     O   #  And then sum each inner row together
 }       # After the loop:
  -      # Subtract the values at the same positions in the two matrices from each other

As for the rest of the program to create all possible matrices and brute-force over them:
Ä        # Convert each integer in the (implicit) input-matrix `M` to its absolute value
 à       # Pop and push the flattened maximum of these absolute values of `M`
  D(     # Duplicate it, and negate the maximum
    Ÿ    # Pop both and push a list in the range [max, -max]
s        # Swap to get the input matrix again
 g       # Pop and push its length to get the matrix dimension
  ã      # Create the cartesian product of the [max, -max]-list that many times
         # (i.e. if the input-matrix has a dimension of 3, we'll now have all possible 
         #  triplets of integers within the [max, -max]-range)
sgã      # Do the same, to create all possible combinations of these triplets
         # (so we now have all possible matrices of the same dimensions as the input-matrix,
         #  using integers within the [max, -max]-range)
ã        # Cartesian product yet again (defaults to 2) to create all possible pairs of matrices
.Δ       # Then find the first pair of matrices which is truthy for:
  D      #  Duplicate the current pair
   2FsN._`øδ*O}-
         #  Calculate M as mentioned above
    Q    #  And check if its equal to the (implicit) input-matrix
         # (after which this matrix-pair is output implicitly as result)

Just to give an idea of how inefficient this brute-force approach is:
The given example input-matrix \$M=\$ [[11,12,12],[1,1,1],[1,2,3]] of the challenge description, will create all possible matrix-pairs of the same 3x3 dimensions with integers in the range \$[-12,12]\$. This means the list of matrix-pairs will contain \$(((12^2+1)^3)^3)^2=802830827198685151406498569488525390625\$ matrix-pairs.
And the given example input-matrix \$M=\$ [[-11811,-9700,-2937],[3,14,15],[2,71,82]] would result in \$(((11811^2+1)^3)^3)^2=400239694511528196808960294025932709663435907226315389994060779277222573307615163703296535923733538440194289670070277830289553752923910908747710464\$ matrix-pairs..