  In this paper we consider functions in the Hardy space $\mathbf{H}_2^{p\times
q}$ defined in the unit disc of matrix-valued. We show that it is possible, as
in the scalar case, to decompose those functions as linear combinations of
suitably modified matrix-valued Blaschke product, in an adaptive way. The
procedure is based on a generalization to the matrix-valued case of the maximum
selection principle which involves not only selections of suitable points in
the unit disc but also suitable orthogonal projections. We show that the
maximum selection principle gives rise to a convergent algorithm. Finally, we
discuss the case of real-valued signals.
