8. Normal Equations#

8.1. Constructing Normal Equations#

Consider the fixed linear model

\[ \mathbf{y} = \mathbf{X\beta} + \mathbf{e}, \]

where \(\mathbf{y}\) is a \(n\times1\) vector of observations, \(\mathbf{X}\) is a \(n\times p\) known matrix, \(\mathbf{\beta}\) is a \(p\times 1\) vector of fixed effects, and \(\mathbf{e}\) is an \(n\times 1\) vector of residuals that are often assumed to be identically and independently distributed with null mean and variance \(\sigma^2_e\).

The normal equations for this model are

\[ (\mathbf{X}'\mathbf{X})\hat{\mathbf{\beta}} = \mathbf{X}'\mathbf{y}, \]

where \(\mathbf{X'X}\) is \(p\times p\) and \(\mathbf{X'y}\) is \(p\times 1\). In breeding applications, \(p\) may be large, but \(\mathbf{X'X}\) is often very sparse, i.e., has very few non-zero elements. Efficient algorithms take advantage of this sparse structure to reduce computing time and storage requirements.

In this section, we will consider how to compute a full-stored \(\mathbf{X'X}\) matrix and \(\mathbf{X'y}\) vector.

  • One way to proceed would be to first build \(\mathbf{X}\) and then get \(\mathbf{X'X}\) by matrix multiplication.

  • Next, we will look more closely at a one-way model and write a simple program for constructing the normal equations for this model.

8.2. Data#

Consider the following data from a hypothetical one-way experiment with four levels of one factor.

using DataFrames
data = DataFrame(x=[1,1,2,2,2,2,3,3,4,1],y=[1.1,1.2,1.9,1.2,2.0,1.7,1.0,1.7,1.1,1.7])

10 rows × 2 columns

xy
Int64Float64
111.1
211.2
321.9
421.2
522.0
621.7
731.0
831.7
941.1
1011.7

8.3. One way model#

The \(\mathbf{X}\) matrix for the one-way model

\[ y_{ij} = \mu + \alpha_i + e_{ij} \]

is

\[\begin{split} \mathbf{X} = \begin{bmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 0 \\ \end{bmatrix} \end{split}\]

Note \(\mathbf{X'X}\) is singular, but a solution can be obtained as follows.

8.4. Computing X’X as the product of full-stored X’ and X#

n = size(data,1)
p = length(unique(data[!,:x]))
X = zeros(n,p);

@time for i = 1:n
    j = data[!,:x][i]
    X[i,j] = 1.0
end

  0.001662 seconds (1.34 k allocations: 74.640 KiB, 99.19% compilation time)

X = [ones(n) X]

10×5 Matrix{Float64}:
 1.0  1.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  0.0  1.0
 1.0  1.0  0.0  0.0  0.0

y = data[!,:y];

lhs = X'X
rhs = X'y;

8.4.1. Solve normal equations#

sol = lhs\rhs

5-element Vector{Float64}:
 0.7500000000000002
 0.583333333333333
 0.9499999999999996
 0.5999999999999999
 0.3499999999999999

8.4.2. Check solution#

[lhs*sol rhs]

5×2 Matrix{Float64}:
 14.6  14.6
  4.0   4.0
  6.8   6.8
  2.7   2.7
  1.1   1.1

In this section, we will consider how to compute a sparse-stored \(\mathbf{X'X}\) matrix and \(\mathbf{X'y}\) vector.

  • One way to proceed would be to first build \(\mathbf{X}\) as a sparse matrix and then get \(\mathbf{X'X}\) by matrix multiplication.

  • Now we will use the same one-way model and write a simple program for constructing the normal equations for this model.

8.5. Computing X’X as the product of sparse-stored X’ and X#

using SparseArrays, LinearAlgebra
X = sparse(1:n,data[!,:x],1.0)
X =[ones(n) X]

10×5 SparseMatrixCSC{Float64, Int64} with 20 stored entries:
 1.0  1.0   ⋅    ⋅    ⋅ 
 1.0  1.0   ⋅    ⋅    ⋅ 
 1.0   ⋅   1.0   ⋅    ⋅ 
 1.0   ⋅   1.0   ⋅    ⋅ 
 1.0   ⋅   1.0   ⋅    ⋅ 
 1.0   ⋅   1.0   ⋅    ⋅ 
 1.0   ⋅    ⋅   1.0   ⋅ 
 1.0   ⋅    ⋅   1.0   ⋅ 
 1.0   ⋅    ⋅    ⋅   1.0
 1.0  1.0   ⋅    ⋅    ⋅ 

Matrix(X)

10×5 Matrix{Float64}:
 1.0  1.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  0.0  1.0
 1.0  1.0  0.0  0.0  0.0

lhs = X'X
rhs = X'y

5-element Vector{Float64}:
 14.599999999999998
  4.0
  6.8
  2.7
  1.1

QRLhs = qr(lhs) 
sol1 = QRLhs\rhs

5-element Vector{Float64}:
 1.0999999999999992
 0.23333333333333364
 0.6000000000000002
 0.2500000000000008
 0.0

[lhs*sol1 rhs]

5×2 Matrix{Float64}:
 14.6  14.6
  4.0   4.0
  6.8   6.8
  2.7   2.7
  1.1   1.1

8.5.1. Estimation of \(\mathbf{k}'\boldsymbol{\beta}\)#

If \(\mathbf{k}'\boldsymbol{\beta}\) is estimable, its estimate is: \(\mathbf{k}'\hat{\boldsymbol{\beta}}\), where \(\hat{\boldsymbol{\beta}}\) is a solution to the normal equations.

8.5.2. Checking estimability of \(\mathbf{k}'\boldsymbol{\beta}\)#

It can be shown that if \(\mathbf{k}'\boldsymbol{\beta}\) is estimable,

\[ (\mathbf{X}'\mathbf{X})(\mathbf{X}'\mathbf{X})^{-}\mathbf{k} = \mathbf{k}. \]

We will see below how this can be done without computing a generalized inverse of \((\mathbf{X}'\mathbf{X})\).

  1. First get a solution to: \((\mathbf{X}'\mathbf{X})\mathbf{b} = \mathbf{k}\).

  2. Then, compute \(\mathbf{a} = (\mathbf{X}'\mathbf{X})\mathbf{b}\).

  3. Check if \(\mathbf{a}\) is identical to \(\mathbf{k}\)

8.5.2.1. Numerical Example#

k = [0.0;1.0;-1.0;0;0.0]

5-element Vector{Float64}:
  0.0
  1.0
 -1.0
  0.0
  0.0

b = lhs\k
round.(lhs*b,digits=3)

5-element Vector{Float64}:
  0.0
  1.0
 -1.0
 -0.0
  0.0