Monte Carlo Sampling Offers Faster Route to Genomic Prediction Accuracy in Hybrid Breeding

Every breeding program runs on a quiet promise: that the numbers attached to each plant or animal, the estimated genetic values that decide which lines move forward and which are discarded, actually mean something. Measuring how much those numbers can be trusted has long been one of the most computationally punishing tasks in quantitative genetics. A new study published in Theoretical and Applied Genetics shows that a family of simulation-based techniques, originally developed for livestock evaluations, can deliver trustworthy measures of prediction accuracy in the kind of multi-trait genomic models that dominate modern hybrid crop breeding, without the brutal matrix algebra that exact methods demand.

The research, led by Antero Heikkilä of the University of Jyväskylä together with colleagues at the Natural Resources Institute Finland, the University of Helsinki, and collaborators, tackled a problem that grows exponentially worse as breeding datasets expand. In genomic best linear unbiased prediction, or GBLUP, the gold standard for measuring the precision of an estimated genetic value is the prediction error variance, or PEV. This quantity describes how far a prediction is expected to stray from the true genetic value under repeated sampling, and it flows directly from the inverse of the coefficient matrix of the mixed model equations. Reliability, the squared correlation between true and estimated genetic value, follows from the PEV and the genetic variance. The catch is that inverting that coefficient matrix scales cubically with its size, and the matrix grows with both the number of individuals and the number of traits being modeled simultaneously.

For a modest dataset, a modern computer can grind through the inversion and produce exact PEVs for every individual. But breeding programs increasingly juggle tens of thousands of genotyped lines, multiple heterotic pools, and several traits measured across environments, all of which must be treated as correlated character states in a multiple-trait framework. The coefficient matrix for such a model can reach dimensions in the hundreds of thousands, at which point exact inversion becomes not merely slow but practically impossible, both in time and in memory. Storing a 160,000-dimensional coefficient matrix in double precision alone requires roughly 205 gigabytes, and the space complexity grows with the square of the number of random effects multiplied by the square of the number of traits.

The Finnish team’s answer was to sidestep the inversion entirely by borrowing a strategy from the animal breeding literature. The core idea, first laid out by García-Cortés and colleagues in the mid-1990s and refined by Hickey and colleagues in 2009, is deceptively simple: instead of computing PEVs analytically, simulate genetic values from their assumed distributions many times, solve the mixed model equations for each simulated dataset, and compare the simulated true values with their predictions. The discrepancies across these Monte Carlo replicates carry the information needed to approximate the PEV. Because the method never requires the inverse of the coefficient matrix, it can lean on fast iterative solvers such as the preconditioned conjugate gradient algorithm, which scales far more gently with problem size.

What makes the new study distinctive is its application to a multiple-trait GBLUP model tailored to hybrid breeding, a setting quite different from the single-trait animal models where these formulas were born. The researchers built a model with three random effect groups: the general combining abilities of parent lines from two heterotic pools, and the specific combining abilities of their crosses. General combining ability captures the additive genetic effects a parent transmits to progeny, while specific combining ability captures non-additive effects such as dominance and epistasis, which are especially important in hybrid crops like maize. The team simulated a realistic maize scenario with two parent pools of 1,000 doubled-haploid lines each, 10,000 F1 hybrids, two traits resembling grain yield and disease resistance, and two environments, yielding four correlated character states and 48,000 genetic values whose accuracy needed assessing.

The simulation machinery came in two flavors. Marker-based Monte Carlo generates genetic values directly from centered SNP marker matrices, requiring computational effort that scales linearly with the number of genotyped individuals and markers. Cholesky Monte Carlo instead factors the genomic relationship matrices and multiplies the resulting triangular factors by random standard normal draws, an approach that becomes cheaper when the number of genotyped individuals is smaller than twice the number of markers. For the specific combining ability group, the researchers faced a combinatorial obstacle: the full interaction marker matrix for all possible crosses would have dimensions of millions by tens of thousands. They circumvented this with a Hadamard-product construction that approximates the Kronecker structure of the cross relationship matrix, and verified empirically that the approximation left the results untouched.

Four previously published approximation formulas were put through their paces. PEV1 and PEV2, both from García-Cortés and colleagues, take complementary views of the same quantity, with PEV1 comparing total genetic variance against the mean squared prediction and PEV2 averaging the squared differences between simulated and predicted values directly. PEV3 combines the two as a variance-weighted average, and NF2, proposed by Hickey’s group, blends PEV2 with the empirical variance of the predictions. The results were strikingly consistent: with 500 Monte Carlo samples, correlations between exact and approximated reliabilities exceeded 0.99 for three of the four formulas, root mean squared errors peaked at 0.047, and regression slopes sat close to one with intercepts near zero. Crucially, the methods proved unbiased across all groups and traits.

The formulas were not interchangeable, however, and the pattern of their differences told a coherent statistical story. PEV1 approximated low reliabilities well but grew noisy as reliability climbed, while PEV2 showed the opposite behavior, excelling at high reliabilities and struggling at the low end. PEV3 and NF2, by pooling information, delivered stable performance across the whole reliability range. The team traced these patterns to the asymptotic sampling variances of each estimator, showing that the realized sampling variances across their simulations matched theoretical predictions. Their practical advice to practitioners is to report root mean squared error and maximum absolute difference rather than relying on correlations alone, which can look flattering on a wide reliability scale yet mask individual errors exceeding 0.2.

The computational benchmarks revealed where the Monte Carlo approach earns its keep. For the largest coefficient matrix tested, at dimension 160,000, exact Cholesky inversion took nearly three hours, while a single solve with the preconditioned conjugate gradient algorithm took under two minutes. Because Monte Carlo replicates are independent, the workload distributes perfectly across computing nodes: with 50 nodes and 500 samples, the researchers estimated roughly 35 minutes of wall-clock time, more than four times faster than the exact method, and the advantage widens as datasets grow since inversion scales cubically while iterative solving scales closer to quadratically. The iterative framework also supports out-of-core computing, reading data from disk rather than holding the full coefficient matrix in memory.

The authors are candid about the method’s limits. The approach assumes that genetic and residual variance components are known without error, whereas in practice they must be estimated first, and any misspecification propagates directly into the accuracy approximations. Solving the equations hundreds of times remains a genuine burden for the very largest datasets, even with parallelization. Yet the study demonstrates something breeders have needed: that decades-old sampling formulas, designed for animal models, transfer cleanly to the multiple-trait, non-additive world of hybrid plant breeding. As genomic datasets balloon, the ability to quantify prediction confidence cheaply and reliably may prove as consequential as the predictions themselves, giving breeders a clearer view not just of which crosses look promising, but of how much they should believe the numbers saying so.

Subject of Research: Monte Carlo approximation of prediction error variances and reliabilities in multiple-trait genomic prediction for hybrid breeding

Article Title: Approximating prediction error variances and reliabilities in a multiple-trait genomic prediction model using Monte Carlo sampling

Article References: Heikkilä, A., Strandén, I., Lidauer, M. H., Nordhausen, K., & Taskinen, S. (2026). Approximating prediction error variances and reliabilities in a multiple-trait genomic prediction model using Monte Carlo sampling. Theoretical and Applied Genetics, 139(9), Article 250. https://doi.org/10.1007/s00122-026-05336-0

Image Credits: AI Generated

DOI: 10.1007/s00122-026-05336-0

Keywords: genomic prediction, GBLUP, Monte Carlo sampling, prediction error variance, reliability, hybrid breeding, maize, mixed model equations, plant breeding, quantitative genetics, preconditioned conjugate gradient, combining ability