In this paper, a fast recursive algorithm with three grade QR decomposition for multi-variable transfer function matrix identification is presented. The algorithm avoids the repeated processing of the measured data of subsystems identification and a parameter estimation that minimizes the overall loss function can be obtained without expanding of parameter vector dimensions. The simulation indicates the validity of the algorithm.