The concept and methods are well accepted in subspace identification of linear multivariable systems. However with regards to bilinear systems, a major drawback of most of the subspace identification methods is to induce enormous dimension of the data matrices, which grows exponentially with the increase of the model order. Accordingly huge storage and computation burden have prohibited the use of subspace identification methods for bilinear systems in the modeling of many industrial bilinear processes. In this paper, a computationally efficient subspace identification procedure for bilinear systems is proposed to provide a solution to tackle the computational difficulties. The new square data matrices are formatted with much smaller dimension than traditional approaches and the QR factorization is replaced with a fast Cholesky factorization based on displacement structure theory. A fast subspace decomposition (FSD) is developed to replace traditional singular value decomposition (SVD) algorithm, then Kalman state estimates can be extracted from a large space with less computation time. Finally, two case studies, a typical bilinear dynamic plant and a real nonlinear process Continuous Stirred Tank Reactor (CSTR), are presented to show the efficiency of this identification method