Given a multi-input, multi-output linear time-invariant system, equations have been derived to characterize a system of lower dimension which optimally approximates the given system with respect to the sum of squared errors of the impulse response. Here expressions are developed for the derivative of the penalty function with respect to all parameters of the approximating system, thus enabling the use of gradient methods to minimize approximation errors.