# 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=
-->
Sign in with Wallet
Connect another wallet