9. Exercise: Computing full-stored or sparse-stored X’X without matrix multiplication#

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

9.1. Data#

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

	x	y
	Int64	Float64
1	1	1.1
2	1	1.2
3	2	1.9
4	2	1.2
5	2	2.0
6	2	1.7
7	3	1.0
8	3	1.7
9	4	1.1
10	1	1.7

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

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

\[\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 that any row of \(\mathbf{X}\) contains only two non-zero elements. These correspond to \(\mu\) and \(\alpha_i\) in the model. Recall that \(\mathbf{x}_i\) denotes row \(i\) of \(\mathbf{X}\). The first element of \(\mathbf{x}_i\) corresponds to \(\mu\). Thus, all \(\mathbf{x}_i\) will contain a “1” in this position. The second element of \(\mathbf{x}_i\) corresponds to \(\alpha_1\). Thus, \(\mathbf{x}_1\), \(\mathbf{x}_2\) and \(\mathbf{x}_{10}\) contain a “1” in this position because observations 1, 2 and 10 are from treatment 1. So, the contribution from the first observation, for example, to the \(\mathbf{X'X}\) matrix is

\[\begin{split} \mathbf{x}_1\mathbf{x}'_1 = \begin{bmatrix} 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} \begin{bmatrix} 1 & 1 & 0 & 0 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} \end{split}\]

The contributions from observations 2 and 10 would be identical to this because \(\mathbf{x}_2\) and \(\mathbf{x}_{10}\) are identical to \(\mathbf{x}_1\). In general, set \(pos\_\mu = 1\), which is the column position of the “1” in any \(\mathbf{x}_i'\) corresponding to \(\mu\), and set \(pos\_\alpha\) equal to the column position for the “1” in any \(\mathbf{x}_i'\) corresponding to \(\alpha_i\). Then, the positions of the contributions to \(\mathbf{X'X}\) from any observations are: (\(pos\_\mu\),\(pos\_\mu\)), (\(pos\_\mu\),\(pos\_\alpha\)), (\(pos\_\alpha\),\(pos\_\mu\)) and (\(pos\_\alpha\),\(pos\_\alpha\)). Further, in the one-way model, the contribution to each of these positions is a “1”. So, \(\mathbf{X'X}\) can be constructed efficiently by setting \(pos\_\alpha=1+ilevel\), where ilevel is the level of the factor A for observation \(i\), and adding “1” to positions (\(pos\_\mu\),\(pos\_\mu\)), (\(pos\_\mu\),\(pos\_\alpha\)), (\(pos\_\alpha\),\(pos\_\mu\)) and (\(pos\_\alpha\),\(pos\_\alpha\)) in \(\mathbf{X'X}\) for each observation in the data file. Similarly, \(\mathbf{X'y}\) can be constructed efficiently by adding \(\mathbf{y}_i\) to positions \(pos\_\mu\) and \(pos\_\alpha\) in \(\mathbf{X'y}\). This strategy is used in the program given below.

9.1.1. Number of levels for \(\alpha\)#

nlevels = length(unique(data[!,:x]));

9.1.2. Make X’X and X’y#

p = nlevels + 1
lhs = fill(0.0,p,p)
rhs = fill(0.0,p);
for i in 1:size(data,1)
    pos_μ              = 1
    pos_α              = 1 + data[i,:x]
    lhs[pos_μ,pos_μ]  += 1.0
    lhs[pos_μ,pos_α]  += 1.0
    lhs[pos_α,pos_μ]  += 1.0
    lhs[pos_α,pos_α]  += 1.0
    
    y            = data[i,:y]
    rhs[pos_μ]  += y
    rhs[pos_α]  += y   
end

lhs

5×5 Matrix{Float64}:
0  3.0  4.0  2.0  1.0
0  3.0  0.0  0.0  0.0
0  0.0  4.0  0.0  0.0
0  0.0  0.0  2.0  0.0
0  0.0  0.0  0.0  1.0

rhs

5-element Vector{Float64}:
599999999999998
0
8
7
1

9.1.3. Solution#

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

sol = lhs\rhs

5-element Vector{Float64}:
  1.7458333333333331
 -0.4125000000000002
 -0.04583333333333354
 -0.3958333333333334
 -0.6458333333333334

9.1.4. Verify solution#

[lhs*sol rhs]

5×2 Matrix{Float64}:
6  14.6
0   4.0
8   6.8
7   2.7
1   1.1

9.2. comparison of efficiencies of different ways to calculate X’X#

9.3. Big Example#

using Distributions

9.3.1. Generate data#

n       = 1_000_000
p       = 1000
levels  = sample(1:p,n)
α       = randn(p)
y       = [α[i] .+ randn() for i in levels];

9.4. Computing X’X as product of full-stored (X’ and X)#

X = fill(0.0,(n,p+1));

@time for i = 1:n
    j      = levels[i] + 1
    X[i,1] = 1.0
    X[i,j] = 1.0
end

  0.169930 seconds (3.98 M allocations: 75.968 MiB, 2.12% gc time)

@time lhs = X'X;

  7.706762 seconds (2.46 M allocations: 129.639 MiB, 6.73% compilation time)

9.5. Computing full-stored X’X without matrix multiplication#

lhs = fill(0.0,p,p);

@time for i in 1:length(levels)
    pos_α               = levels[i]
    lhs[pos_α,pos_α]   += 1.0
end

  0.155171 seconds (5.49 M allocations: 99.004 MiB)

9.6. Computing sparse-stored X’X as product of sparse-stored (X’ and X)#

using SparseArrays, LinearAlgebra
ii = 1:n
@time X = sparse(ii,levels,1.0);

  0.175950 seconds (359.91 k allocations: 65.146 MiB, 85.41% compilation time)

X = [ones(n) X]

1000000×1001 SparseMatrixCSC{Float64, Int64} with 2000000 stored entries:
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿
⣿

@time lhs = X'X;

  0.408823 seconds (419.04 k allocations: 76.342 MiB, 38.05% compilation time)

lhs

1001×1001 SparseMatrixCSC{Float64, Int64} with 3001 stored entries:
⡟⢍⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉
⡇⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀
⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄

rhs = X'y;

QRLhs = qr(lhs)
sol = QRLhs\rhs
[lhs*sol rhs]

1001×2 Matrix{Float64}:
 21777.6      21777.6
    47.8103      47.8103
   633.311      633.311
 -1072.6      -1072.6
     3.79227      3.79227
  -934.483     -934.483
  2224.34      2224.34
   458.346      458.346
  1531.21      1531.21
 -1023.51     -1023.51
    76.7795      76.7795
 -1028.66     -1028.66
  -428.257     -428.257
     ⋮        
  -668.414     -668.414
  -149.372     -149.372
   432.522      432.522
  -617.182     -617.182
  1041.06      1041.06
   465.445      465.445
  -377.921     -377.921
   569.776      569.776
  -541.006     -541.006
  -317.002     -317.002
  -250.474     -250.474
  1137.16      1137.16

9.7. Computing sparse-stored X’X without matrix multiplication#

lhs = spzeros(p,p)

1000×1000 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀

@time for i in 1:length(levels)
    pos_α  = levels[i]
    lhs[pos_α,pos_α]   += 1.0
end

  0.235890 seconds (6.55 M allocations: 147.876 MiB, 20.27% compilation time)

Modern Programming in Julia for Genomic Prediction

Exercise: Computing full-stored or sparse-stored X’X without matrix multiplication

Contents