Recently developed subspace-based system identification (4SID) techniques have opened new routes to the identification of multi-input multi-output systems. The 4SID techniques guarantee convergence, and run faster than the statistically efficient prediction error methods without much performance loss. The resulting computational load of the 4SID techniques is O(NM(2)), where N is the data length and M is the sliding window size. However, the computational burden O(NM(2)) can become prohibitively large as N and M grow large. Noting that the major bottleneck comes from the QR factorization of an M X N data matrix and that the existing 4SID techniques do not exploit the structure of the matrices arising in the identification procedure, we propose a new implementation of the existing 4SID, which reduces the computational burden to O(NM) by exploiting the displacement and low-rank structure of the matrices.