Summary: It has been demonstrated by \textit{N. T. Young} [NATO ASI Series F34 (1987)] that given a stable matrix-valued function \(G_ 0(s)\) and a nonnegative integer \(k\), there exists a unique superoptimal approximation \(\Phi(s)\) with no more than \(k\) poles in the left half plane that minimizes the sequence \((s^ \infty_ 1(G_ 0+\Phi),\;s^ \infty_ 2(G_ 0+\Phi),\dots)\), with respect to lexicographic ordering, where \(s^ \infty_ i(G_ 0+\Phi):=\sup_ \omega[s_ i(G_ 0+\Phi)(j\omega)]\) and \(s_ i(\cdot)\) are the singular values in descending order of magnitude. This paper presents a constructive state- space algorithm that evaluates the superoptimal approximating matrix function. The procedure recursively minimizes each frequency-dependent singular value with the aid of all-pass transformations constructed from the \(k\)th Schmidt pairs of a sequence of Hankel operators. The algorithm may be stopped after an arbitrary number of, say, \(\ell\leq \min(m,p)\) steps. The representation formula at the \(\ell\)th stage will characterize all matrix functions that have \(\leq k\) poles in the left half plane and that minimize \(s^ \infty_ 1(G_ 0+\Phi),\dots,s^ \infty_ \ell(G_ 0+\Phi)\).