---
# System prepended metadata

title: Normal Equations

---


# Normal Equations

### 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.

#### Data

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


```julia
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])
```




<table class="data-frame"><thead><tr><th></th><th>x</th><th>y</th></tr><tr><th></th><th>Int64</th><th>Float64</th></tr></thead><tbody><p>10 rows × 2 columns</p><tr><th>1</th><td>1</td><td>1.1</td></tr><tr><th>2</th><td>1</td><td>1.2</td></tr><tr><th>3</th><td>2</td><td>1.9</td></tr><tr><th>4</th><td>2</td><td>1.2</td></tr><tr><th>5</th><td>2</td><td>2.0</td></tr><tr><th>6</th><td>2</td><td>1.7</td></tr><tr><th>7</th><td>3</td><td>1.0</td></tr><tr><th>8</th><td>3</td><td>1.7</td></tr><tr><th>9</th><td>4</td><td>1.1</td></tr><tr><th>10</th><td>1</td><td>1.7</td></tr></tbody></table>



### One way model

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

$$
y_{ij} = \mu + \alpha_i + e_{ij}
$$

is

$$
\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}
$$

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

### Computing X'X as the product of full-stored X' and X


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


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

      0.000011 seconds (11 allocations: 352 bytes)



```julia
X = [ones(n) X]
```




    10×5 Array{Float64,2}:
     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




```julia
y = data[:y];
```


```julia
lhs = X'X
rhs = X'y;
```


```julia
sol = lhs\rhs
```




    5-element Array{Float64,1}:
     -0.04583333333333304
      1.3791666666666667 
      1.745833333333333  
      1.3958333333333333 
      1.1458333333333333 



#### Verify solution 


```julia
[lhs*sol rhs]
```




    5×2 Array{Float64,2}:
     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.

### Computing X'X as the product of sparse-stored X' and X


```julia
using SparseArrays
X = sparse(1:n,data[:x],1.0);
```


```julia
X =[ones(n) X];
```


```julia
y = sparse(y);
```


```julia
lhs = X'X
rhs = X'y;
```

### [Supplemental Note: Computing full-stored or sparse-stored X'X without matrix multiplication](SupplementalNote/1.SupplementalNote.ipynb)


# Exercise

This exercise will provide an opportunity to practice your matrix calculation skills using Julia in the context of comparison of efficiencies of different ways to calculate the left-hand-side of the normal equations.

Ues the @time macro to measure performance of functions below

1. Compute the left-hand-side of the normal equations as the product of full-stored $\mathbf{X}^{T}$ and $\mathbf{X}$
2. Compute the left-hand-side of the normal equations as the product of sparse-stored $\mathbf{X}^{T}$ and $\mathbf{X}$

#### Data


```julia
using Distributions,DataFrames
```


```julia
nobs,nlevels = 1_000_000,1000
x = sample(1:nlevels,nobs)
α = randn(nlevels)
y = [α[i] .+ randn() for i in x];
```


```julia
data = DataFrame(x=x,y=y)
first(data,5)
```




<table class="data-frame"><thead><tr><th></th><th>x</th><th>y</th></tr><tr><th></th><th>Int64</th><th>Float64</th></tr></thead><tbody><p>5 rows × 2 columns</p><tr><th>1</th><td>740</td><td>0.140606</td></tr><tr><th>2</th><td>419</td><td>-0.197064</td></tr><tr><th>3</th><td>368</td><td>-0.109038</td></tr><tr><th>4</th><td>426</td><td>-0.478222</td></tr><tr><th>5</th><td>531</td><td>0.00339639</td></tr></tbody></table>



#### Answer 1. Computing X'X as product of full-stored X


```julia
X = zeros(nobs,nlevels);
```


```julia
@time for i = 1:nobs
    j = data[:x][i]
    X[i,j] = 1.0
end
```

      2.637775 seconds (3.51 M allocations: 69.172 MiB, 4.34% gc time)



```julia
@time lhs = X'X;
```

     10.474817 seconds (943.22 k allocations: 52.955 MiB, 0.55% gc time)


#### Answer 2. Computing X'X as product of sparse-stored X


```julia
using SparseArrays
ii = 1:nobs
@time X = sparse(ii,data[:x],1.0);
```

      0.137251 seconds (351.10 k allocations: 62.857 MiB, 8.09% gc time)



```julia
@time lhs = X'X;
```

      0.322128 seconds (290.65 k allocations: 48.350 MiB, 1.94% gc time)


<!--stackedit_data:
eyJoaXN0b3J5IjpbNzkyNDcxMDA4XX0=
-->