This paper addresses the problem of estimating the parameters in a multivariable bilinear model on the basis of observed input-output data. The main contribution is to develop, analyse, and empirically study new techniques for computing a Maximum-Likelihood based solution. In particular, the emphasis here is on developing practical methods that are illustrated to be numerically reliable, robust to choice of initialisation point, and numerically efficient in terms of how computation and memory requirements scale relative to problem size. This results in new methods that can be reliably deployed on systems of non-trivial state, input and output dimension. Underlying these developments is a new approach (in this context) of employing the Expectation-Maximisation method as a means for robust and gradient free computation of the Maximum-Likelihood solution.