  Let $B$ be a locally integrable matrix function, $W$ a matrix A${}_p$ weight
with $1 < p < \infty$, and $T$ be any of the Riesz transforms. We will
characterize the boundedness of the commutator $[T, B]$ on $L^p(W)$ in terms of
the membership of $B$ in a natural matrix weighted BMO space. To do this, we
will characterize the boundedness of dyadic paraproducts on $L^p(W)$ via a new
matrix weighted Carleson embedding theorem. Finally, we will use some of the
ideas from these proofs to (among other things) obtain quantitative weighted
norm inequalities for these operators and also use them to prove sharp $L^2$
bounds for the Christ/Goldberg matrix weighted maximal function associated with
matrix A${}_2$ weights.
