In this paper we consider functions in the Hardy space H defined in the unit disc of matrix-valued functions. We show that it is possible, as in the scalar case, to decompose those functions as linear combinations of suitably modified matrix-valued Blaschke products, 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.