This paper considers the linear approximation and identification of multi-input multi-output (MIMO) Wiener-Hammerstein systems, or LNL systems. Evaluating the input-output cross-covariance matrix of the MIMO LNL system for Gaussian inputs, we show that the best linear approximation of the MIMO LNL system in the mean square sense can be obtained by the orthogonal projection (ORT) subspace identification method. For each allocation of the poles of the best linear approximation between the two linear subsystems, the unknown parameters in the numerators of the linear subsystems and the coefficients of a basis function expansion of the nonlinearity are estimated by applying the separable least-squares. The best LNL system is the one that gives the minimum mean square output error. A numerical example is included to show the feasibility of the present approach. (C) 2016 Elsevier Ltd. All rights reserved.