In this paper, we consider the identification of matrix state-space models (MSSM) of the following form: X(k + 1) = A(2)X(k) A(1)(T) + B2U(k) B-1(T) Y(k) = C2X(k) C-1(T) + E(k) for all time dependent quantities and matrices of appropriate dimensions. Due to the large size of these matrices, vectorization does not allow the use of standard multivariable subspace methods such as N4SID or MOESP. In this paper, the resulting Kronecker structure that appears in the system matrices due to vectorization is exploited for developing a scalable subspace-like identification approach. This approach consists of first estimating the Markov parameters associated to the MSSM via the solution of a regularized bilinear least-squares problem that is solved in a globally convergent manner. Second, a bilinear low-rank minimization problem is tackled which allows to write a three-dimensional low-rank tensor and consequently to estimate the state sequence and the lower-dimensional matrices A(1), A(2), B-1, B-2, C-1, C-2. A numerical example on a large-scale adaptive optics system demonstrates the ability of the algorithm to handle the identification of state-space models within the class of Kronecker structured matrices in a scalable manner which results in more compact models.