Estimates of large genetic covariance matrices are commonly obtained by pooling results from a series of analyses of small subsets of traits. Procedures available to pool the part-estimates differ in their efficacy in accounting for unequal accuracies of estimates and sampling correlations, and ensuring that pooled matrices are within the parameter space. We propose a maximum likelihood (ML) approach to combine estimates, treating sets from individual part-analyses as matrices of mean squares and cross-products from independent families. This facilitates simultaneous pooling of estimates for all sources of variation considered, readily allows for weighted estimation or a given structure of the pooled matrices, and provides a framework for regularized estimation by penalizing the likelihood. A simulation study is presented, comparing the quality of combined estimates for several procedures, including truncation or shrinkage of either canonical or individual matrix eigen-values, iterative summation of expanded part matrices, and the ML approach, considering a range of penalties. Shrinking eigen-values of individual matrices towards their mean reduced losses in the pooled estimates, but substantially increased proportional losses in their phenotypic counterparts and thus yielded estimates differing most from corresponding full multivariate analyses of all traits. Assuming a simple pseudo-pedigree structure when combining estimates for all random effects simultaneously using ML allowed sampling correlations between estimates of different components from the same part-analysis to be approximated sufficiently to yield pooled matrices closest to full multivariate results, with little change in phenotypic components. Imposing a mild penalty to shrink matrices for random effects towards their sum proved highly advantageous, markedly reducing losses in estimates and more than compensating for the reduction in efficiency of using the data inherent in analyses by parts. Penalized ML provides a flexible alternative to current methods for pooling estimates from part-analyses with good sampling properties, and should be adopted more widely.