  In this paper we prove that the space of two parameter, matrix-valued BMO
functions can be characterized by considering iterated commutators with the
Hilbert transform. Specifically, we prove that $$\| B \|_{BMO} \lesssim \|
[[M_B, H_1],H_2] \|_{L^2(\mathbb{R}^2;\mathbb{C}^d) \rightarrow
L^2(\mathbb{R}^2;\mathbb{C}^d)} \lesssim \| B \|_{BMO}.$$ The upper estimate
relies on Petermichl's representation of the Hilbert transform as an average of
dyadic shifts, and the boundedness of certain paraproduct operators, while the
lower bound follows Ferguson and Lacey's proof for the scalar case.
