Obtaining accurate estimates of the genetic covariance matrix ∑G for multivariate data is a fundamental task in quantitative genetics and important for both evolutionary biologists and plant or animal breeders. Classical methods for estimating ∑G are well known to suffer from substantial sampling errors; importantly, its leading eigenvalues are systematically overestimated. This article proposes a framework that exploits information in the phenotypic covariance matrix ∑P in a new way to obtain more accurate estimates of ∑G. The approach focuses on the "canonical heritabilities" (the eigenvalues of ∑P⁻¹∑G), which may be estimated with more precision than those of ∑G because ∑P is estimated more accurately. Our method uses penalized maximum likelihood and shrinkage to reduce bias in estimates of the canonical heritabilities. This in turn can be exploited to get substantial reductions in bias for estimates of the eigenvalues of ∑G and a reduction in sampling errors for estimates of ∑G. Simulations show that improvements are greatest when sample sizes are small and the canonical heritabilities are closely spaced. An application to data from beef cattle demonstrates the efficacy this approach and the effect on estimates of heritabilities and correlations. Penalized estimation is recommended for multivariate analyses involving more than a few traits or problems with limited data.