We discuss identification of discrete-time bilinear state space systems with multiple inputs and multiple outputs. Subspace identification methods for bilinear systems suffer from the curse of dimensionality. Already for relatively low order systems, the matrices involved become so large that the method cannot be used in practice. We have modified the subspace algorithm such that it reduces the dimension of the matrices involved. Only the rows that have the largest influence on the model are selected; the remaining rows are discarded. This obviously leads to an approximation error. The initial model that we get from the subspace method is optimized using the principle of separable least squares. According to this principle, we can first solve for the matrices that enter non-linearly in the output error criterion and then obtain the others by solving a linear least squares problem.