Benchmarking Solutions to Mrode (4th Ed.): BLUPF90, R, and Julia

The publication of the 4th edition of Linear Models for the Prediction of the Genetic Merit of Animals by R. A. Mrode and I. Pocrnic (2023, CABI) provides a comprehensive pedagogical reference for linear mixed models in quantitative genetics and animal breeding. Accompanying the text are three open-source companion solution repositories implementing the worked examples across the textbook:

  1. BLUPF90 Suite (Examples_BLUPF90): Parameter and data setups for the legacy Fortran 90 binaries developed by Ignacy Misztal and colleagues at the University of Georgia.
  2. R Scripts (Examples_R): Educational scripts using standard R packages (Matrix, pedigreemm, lme4, and the tidyverse).
  3. Julia Package (Mrode4edExamples.jl / Examples-Julia): A native, type-stable implementation built on the domain-specific libraries RelationshipMatrices.jl and Breeding.jl.

Below is an empirical benchmark and architectural evaluation comparing the computational paradigms, execution performance, numerical reliability, and algorithmic scaling of these three solution suites.

Benchmark Comparison

Execution times and memory footprints were measured across 14 representative models on an Intel Core i9 system running Linux x86_64:

  • Julia (Warm): Median in-memory execution time measured via BenchmarkTools.jl (excluding JIT compilation latency) and heap memory allocations.
  • BLUPF90 (blupf90+): Median wall-clock execution of the compiled Fortran binary (including process startup, Intel oneAPI MKL initialization, sparse assembly, and text disk I/O).
  • R (Rscript): Median wall-clock execution of standalone R script execution (including R interpreter startup, package dependencies, and matrix operations).

Benchmark runtime comparison across Mrode 4th edition models

Detailed Micro-Benchmark Table

Chapter / Model Paradigm & Effects Julia (Warm) Julia Memory BLUPF90 (blupf90+) R (Rscript)
Ex 4.1: Animal Model Univariate additive animal model, PEV, and reliabilities 0.422 ms 230.0 KB 186.33 ms 2,250.4 ms
Ex 4.2: Sire Model Random sire evaluation 0.180 ms 72.7 KB 184.64 ms 2,168.9 ms
Ex 4.3: Reduced Animal Model RAM with non-parent absorption and back-solving 0.169 ms 88.0 KB 193.77 ms 2,178.6 ms
Ex 4.4: Genetic Groups Phantom parent / unknown parent groups (QP) 0.456 ms 206.4 KB 181.67 ms Failed (Issue #1)
Ex 5.1: Repeatability Model Permanent environmental effects ($pe$) and DYD 0.236 ms 121.4 KB 245.83 ms 1,156.0 ms
Ex 5.2: Common Environment Maternal common environmental effect ($c^2$) 0.193 ms 114.4 KB 237.46 ms 1,163.8 ms
Ex 6.1: Bivariate Model Multi-trait with missing records and unequal designs 0.486 ms 333.5 KB 200.28 ms 1,167.1 ms
Ex 8.1: Maternal Animal Model Direct genetic + maternal genetic ($r_{am}$) + maternal $pe$ 0.236 ms 186.7 KB 251.64 ms 1,179.4 ms
Ex 9.1: Social Interaction Model Associative/social interaction effects and group housing 0.411 ms 383.2 KB 296.11 ms 1,179.1 ms
Ex 10.2: Random Regression Legendre polynomials on additive genetic and $pe$ 0.505 ms 650.9 KB 340.33 ms Failed (Issue #1)
Ex 11.3: GBLUP Genomic relationship matrix ($G$, VanRaden 2008) 0.257 ms 160.5 KB 732.17 ms 1,176.2 ms
Ex 12.1: Single-step GBLUP $H^{-1}$ joint matrix, blending and tuning $G$ to $A_{22}$ 0.528 ms 266.9 KB 707.94 ms 2,059.7 ms
Ex 12.2: ssGBLUP with APY Algorithm for Proven and Young ($G^{-1}_{\text{APY}}$) 0.328 ms 192.0 KB 836.80 ms Failed (Directory bug)
Ex 13.1: Pedigree Dominance Inversion of dominance relationship matrix ($D^{-1}$) 0.318 ms 135.3 KB 242.12 ms 1,715.2 ms

Architectural Characteristics

flowchart TD
    subgraph BLUPF90 ["BLUPF90 Suite (Fortran 90 + MKL)"]
        F1["*.par Parameter File"] --> F2["Disk I/O (*.dat, *.ped, *.gen)"]
        F2 --> F3["Sparse Hash / IJA Assembly"]
        F3 --> F4["FSPAK Direct Cholesky / PCG"]
        F4 --> F5["Disk Outputs (solutions, acc_bf90)"]
    end

    subgraph RSuite ["R Scripts (Pedagogical)"]
        R1["Data Frames & Vectors"] --> R2["Dense model.matrix()"]
        R2 --> R3["crossprod() & explicit block rbind/cbind"]
        R3 --> R4["solve(LHS, RHS) & solve(LHS) for PEV"]
        R4 --> R5["In-memory Matrices (O(N³) scaling)"]
    end

    subgraph JuliaSuite ["Julia Package (Breeding.jl Ecosystem)"]
        J1["DataFrames & Genotype Matrices"] --> J2["Sparse incidence() & ainv(ped)"]
        J3["Symbolic MME Engine: mme()"]
        J2 --> J3
        J3 --> J4["Sparse Solver: solve_mme()"]
        J3 --> J5["Takahashi Sparse PEV / reliability()"]
    end

1. BLUPF90 Suite: Production Out-of-Core Engine

  • Computational Framework: Compiled Fortran 90 binaries linked against Intel oneAPI MKL OpenMP routines and the FSPAK sparse Cholesky package.
  • Workflow: Strictly file-driven. Requires parameter files (*.par) declaring column positions, effect types (cross, cov), and covariance structures. Inputs are space-delimited text files (*.dat, *.ped, *.gen, *.xrf), and outputs are written to disk (solutions, acc_bf90, fspak90.ord).
  • Benchmark Profile: Individual model execution takes 180 ms to 840 ms. On textbook-scale problems ($N = 5$ to $30$ animals), runtime is almost entirely overhead: process spawning, shared library loading (libmkl_intel_lp64.so), text file parsing, and disk writes.
  • Strengths & Limitations: While inefficient for rapid interactive prototyping or tiny pedagogical models, BLUPF90’s sparse hash table assembly and iteration-on-data algorithms are industrial standards capable of scaling to millions of animals and hundreds of traits without memory exhaustion. However, algorithmic internals (e.g., REML iterations, Gibbs sampling loops, and iterative solver steps) are enclosed in opaque executables rather than exposed as modular routines.

2. R Suite: Pedagogical Matrix Algebra

  • Computational Framework: Standalone scripts utilizing Matrix, pedigreemm, lme4, and tidyverse.
  • Workflow: Constructs mixed model equations via explicit block matrix operations:

$$\mathbf{LHS} = \begin{bmatrix} \mathbf{X}’\mathbf{X} & \mathbf{X}’\mathbf{Z} \\ \mathbf{Z}’\mathbf{X} & \mathbf{Z}’\mathbf{Z} + \alpha \mathbf{A}^{-1} \end{bmatrix}$$

  • Benchmark Profile: Script execution takes 1,150 ms to 2,250 ms per invocation. Runtime is dominated by R package startup (loading tidyverse + pedigreemm + lme4 requires 1.1–1.8 s).
  • Strengths & Limitations: Clear pedagogical exposition of direct textbook matrix formulas. However, its memory and computational scaling are prohibitive:
    • Reliance on dense matrix expansion (as.matrix(), model.matrix()) and dense inverses (solve(LHS)) imposes $O(N^3)$ computational time and $O(N^2)$ memory scaling.
    • Calculation of prediction error variances (PEV) and individual reliabilities requires inverting the entire coefficient matrix $\mathbf{C} = \mathbf{LHS}^{-1}$.

3. Julia Package (Mrode4edExamples.jl): Native, Type-Stable Modern Design

  • Computational Framework: Pure, type-stable Julia built on top of RelationshipMatrices.jl and Breeding.jl.
  • Workflow: High-level, allocation-aware mixed model engine using sparse incidence matrices and native sparse relationship matrices (ainv(), hinv(), nrm(), grm()).
  • Benchmark Profile:
    • Warm in-memory execution: 0.17 ms to 0.53 ms per example with minimal heap allocations (70–650 KB).
    • Cold execution: Incurs a Time-to-First-Execution (TTFX) compilation cost (~56 s compilation overhead out of 78.4 s for all 42 examples across the textbook; subsequent replays take ~22 s for the entire book).
  • Algorithmic Sophistication:
    • Sparse PEV / Reliability: Computes reliabilities and PEVs via sparse Takahashi equations (reliability(m, 1)) without forming or inverting the full coefficient matrix $\mathbf{C}$.
    • Colleau’s Algorithm: Extracts submatrices $\mathbf{A}_{22}$ for single-step evaluations without dense formation of $\mathbf{A}$.
    • Complete Algorithmic Coverage: Implements exact underlying mechanisms in native code: Henderson’s rules, APY core/non-core recursions, breed-of-origin partial relationship matrices, Fisher scoring for threshold models, Method 3 / canonical sums of squares, AI- and EM-REML, univariate/multivariate Gibbs samplers, and iterative solvers (Jacobi, Gauss-Seidel, iteration on data, PCG).

Codebase Health & Reproducibility Audit

1. Known Runtime Errors in the R Suite

Testing the published R scripts revealed three distinct categories of runtime failures, as documented in Examples_R Issue #1:

  1. Matrix Dimension Dropping in Chapter 4 (Chapter_04_3.R):
    • Subsetting the right-hand side vector RHS = RHS[-c(11), ] without drop = FALSE simplifies RHS into a 1D vector. Consequently, solve(LHS, RHS) returns a dimensionless vector rather than a matrix, causing subsequent rownames(solutions_QP) <- ... calls to fail with: 
      Error in rownames<-(*tmp*, value = ...): attempt to set 'rownames' on an object with no dimensions
    • Fix: Preserve matrix dimensions: RHS = RHS[-c(11), , drop = FALSE].
  2. Declaration Order Bug in Chapter 10 (Chapter_10.R):
    • Line 61 attempts to factorize animal IDs using levels = pedX@label before pedX and the pedigree data frame are declared (lines 78–88), halting execution with Error: object 'pedX' not found.
    • Fix: Reorder the script to define the pedigree prior to incidence matrix setup.
  3. Deprecated Dependency in Chapter 13 (Chapter_13_2.R):
    • Line 14 calls read_table2(), which was removed in modern versions of the readr package (readr >= 2.0.0), producing could not find function "read_table2".
    • Fix: Replace with read_table().
  4. Directory Coupling in Chapter 12 (Chapter_12_2.R):
    • The script hardcodes source("APY_fja.R"), causing execution to fail if the process working directory is not strictly set to Content/Chapter_12/.

2. Discrepancy Tracking in the Julia Suite

The Julia suite (Mrode4edExamples.jl) incorporates continuous verification against the textbook’s printed values (42 test suites, 303 automated numerical checks). Where printed textbook values cannot be reproduced due to errata, the codebase checks the corrected mathematical solution:

  • Example 6.3: DYD calculation for sire 1 in the book erroneously uses $\text{YD}_{41} = 25.7$ instead of $201 - 175.7 = 25.3$.
  • Example 6.4: Animal 6 solutions in the book were printed with inverted signs (a copy of animal 5).
  • Example 6.5: MACE bull 3 has a typographical error in $\mathbf{W}_1$ (spurious $(1,2)$ element), and bull 2’s PC omits the grandson’s MGD group.
  • Example 10.2: The printed quadratic form $\mathbf{k}’\mathbf{P}\mathbf{k} = 323,463$ deviates from the book’s own stated $\mathbf{t}$ and $\mathbf{P}$ matrices ($320,126$).
  • Example 11.8: Typo in pseudo-SNPs for animal 5 in block 2 propagates through $\mathbf{G}_h$ and $2\sum pq = 4.42$; the correlation printed as $0.999$ is actually $0.992$.
  • Example 17.7: AI-REML $\hat{a}_4$ solution is $+0.254$, correcting the sign printed in the book.

Summary Takeaway

  • For Research, Prototyping, and Education: Mrode4edExamples.jl is currently the most robust and computationally capable companion. It combines the clean mathematical syntax of modern scientific computing with production-grade algorithms (sparse MME assembly, Takahashi equations, Colleau recursions), executing worked models in under 1 millisecond per run.
  • For Industrial Operations: BLUPF90 remains the established workhorse for massive, out-of-core animal breeding evaluations where datasets scale into millions of pedigrees and genotypes.
  • For R Users: The R companion scripts offer straightforward linear algebra walkthroughs, but require the fixes from Issue #1 and should not be used as templates for realistic breeding dataset sizes due to dense matrix memory scaling.