Consider the linear mixed model with normal random effects, $$Y = X\beta + Zv + \epsilon$$ where \(\beta\) and \(v\) are vectors of fixed and random effects. One of most popular methods to solve the Henderson's Mixed Model Equation related to the problem is EM-type algorithm. Its drawback, however, comes from repetitive matrix inversion during recursive estimation steps. Kim (2017) proposed a novel method of avoiding such difficulty, letting the estimation more fast, stable, and scalable.