The use of derivative-free methods to give maximum likelihood estimates of bivariate (co)variance parameters is illustrated. An algorithm is given to estimate the four variance and two covariance parameters of two random effects associated with each trait when both traits are measured on all animals and the same fixed and random model hold for both traits. By reparameterising in terms of canonical heritabilities and a transformation matrix, a six-dimensional problem is reduced to a two-dimensional problem. It is shown how to derive the estimate of the transformation matrix given the values for the canonical heritabilities. Maximization is then only over the two dimensions of canonical heritabilities and the transformation matrix. The use and the properties of the method are illustrated with examples from simulated selection experiment data.