Closed-form EM algorithm matrix manipulations

Indexing the kernel matrrix list Hl

All of the above would be performed in kernL() so ind1 and ind2 would be available in environment. We now enter the indxFn() function. Set k = 1. First find the index for which variable k has interactions (in relation to the positioning in intr). Variable 1 appears in columns 1 and 2 of intr.

k <- 1
ind.int1 <- intr[1, ] == k; ind.int2 <- intr[2, ] == k
(ind.int <- which(ind.int1 | ind.int2))
## [1] 1 2

Which of the Hadamard products (i.e. intr) involve variable k, and where are the relevant Hadamard products in relation to the full index? The reason for the formula below is that the Hadamard products are calculated in the same order that appears in intr, and we add p because the first p elements are the p kernel matrices.

(k.int <- ind.int + p)
## [1] 4 5

Which variables have interaction with variable k? When I wrote this, I was thinking which of the \(\lambda_j\) need to be multipled by \(\lambda_k\), \(j \neq k\)? In other words, what are the other half of the pair of variable k in the matrix intr?

(k.int.lam <- c(intr[1, ][ind.int2], intr[2, ][ind.int1]))
## [1] 2 3

Next I simply call nok the indices of all variables excluding k. sidenote: I am beginning to think that k.int.lam == nok.

(nok <- (1:p)[-k])
## [1] 2 3

Finally, these are the indices of the Hadamard products which do not involve variable k.

(k.noint <- which(!(ind.int1 | ind.int2)) + p)
## [1] 6

Indices for \(\mathbf P_k \mathbf R_k + (\mathbf P_k \mathbf R_k)^\top\)

We have a list called H2l which contains all possible two-way terms Hi %*% Hj + Hj %*% Hi, i,j = 1,...,h which arises as a result of squaring H = H1 + ... + Hh. It is efficient to calculate these two-way terms once at the beginning and then recall them as needed. For our example, the entries of this list consist of

ind <- paste(ind1, ind2, sep = "x")
names(ind) <- as.character(1:length(ind))
ind
##     1     2     3 
## "1x2" "1x3" "2x3"

For clarity, in this example we rename the entries of ind to reflect the three unique scale parameters, as follows:

ind.tmp
##     1     2     3 
## "1x2" "1x3" "2x3"

To remind ourselves, the matrices \(\mathbf P_k\), \(\mathbf R_k\) and the product \(\mathbf P_k \mathbf R_k\) are defined as

\[ \begin{align} \mathbf P_k &= \mathbf H_k + \sum_{j\in M}\lambda_j\mathbf H_{kj} \\ \mathbf R_k &= \sum_{j\neq k} \lambda_j \mathbf H_j \\ \mathbf P_k \mathbf R_k &= \sum_{j\neq k} \lambda_j \mathbf H_k\mathbf H_j + \sum_{j\in M} \sum_{j'\neq k} \lambda_j\lambda_{j'} \mathbf H_{kj} \mathbf H_{j'} \end{align} \]

For now, ignore the scale parameters in the formulae above, and just concentrate on the kernel matrices. For the first part of \(\mathbf P_k \mathbf R_k\), we require the matrix product indices where variables \(k\) is multiplied with all other variables except itself. We call this za.

(za <- which((ind1 %in% k & ind2 %in% nok) | (ind2 %in% k & ind1 %in% nok)))
## [1] 1 2

# Check
ind[za]
##     1     2 
## "1x2" "1x3"

For the second part, it is a double sum of the products of the Hadamard matrices involving variable k, and all the kernel matrices except k. We have already coded the indices as k.int and nok respectively. Note that these indices are in relation to the full index 1, 2, ..., 6. In R, we can use the expand.grid() function to create a data frame from all possible combinations of k.int and nok, which would give us the index of the double sum. We then find the positions of these coordinates in ind1 and ind2.

(grid.PR <- expand.grid(k.int, nok))
##   Var1 Var2
## 1    4    2
## 2    5    2
## 3    4    3
## 4    5    3
(zb <- which((ind1 %in% grid.PR[,1] & ind2 %in% grid.PR[,2]) |
            (ind2 %in% grid.PR[,1] & ind1 %in% grid.PR[,2])
))
## integer(0)

Finally, what’s left is to take care of the scale parameters. Our scale parameters are contained in the vector of length q=6 called lambda. For the first part, that is simply the index nok. For the second part, we need to find the indices using expand.grid() again, but this time using the indices k.int.lam and nok. k.int.lam would give us the indices of the scale parameters which have interactions with k.

(nok)
## [1] 2 3
(grid.PR.lam <- expand.grid(k.int.lam, nok))
##   Var1 Var2
## 1    2    2
## 2    3    2
## 3    2    3
## 4    3    3

The required product \(\mathbf P_k \mathbf R_k\) is then, in a manner of speaking,

sum(lambda[PR.lam index] * H2l[PR index])

Indices for \(\mathbf P_k \mathbf U_k +(\mathbf P_k \mathbf U_k)^\top\)

The formulae of interest are

\[ \begin{align} \mathbf P_k &= \mathbf H_k + \sum_{j\in M}\lambda_j\mathbf H_{kj} \\ \mathbf U_k &= \mathop{\sum\sum}_{k',j \in M \ \& \ k'\neq k} \lambda_{k'} \lambda_j \mathbf H_{k'j} \\ \mathbf P_k \mathbf U_k &= \mathop{\sum\sum}_{k',j \in M \ \& \ k'\neq k} \lambda_j \lambda_{k'} \mathbf H_k \mathbf H_{k'j} + \mathop{\sum\sum}_{j,k',j' \in M \ \& \ k'\neq k} \lambda_j \lambda_{k'} \lambda_{j'} \mathbf H_{kj} \mathbf H_{k'j'} \end{align} \]

The idea is similar to the above, albeit the indices can be a bit confusing. For the first part of the sum, we need the indices of the double sum involving k.noint paired with k. Recall that k.noint are the indices of the Hadamard products which do not involve variable k (mathematically, it is the set \(\{(k',j) \in M : k' \neq k\}\)). We then find the corresponding index from ind1 and ind2.

(grid.PU1 <- expand.grid(k, k.noint))
##   Var1 Var2
## 1    1    6
(zc <- which((ind1 %in% grid.PU1[,1] & ind2 %in% grid.PU1[,2]) |
             (ind2 %in% grid.PU1[,1] & ind1 %in% grid.PU1[,2])))
## integer(0)

# Check
ind[zc]
## named character(0)
ind.tmp[zc]
## named character(0)

For the second part, we use expand.grid() to find the indices for the double sum involving k.int (the Hadamard products involving k) and k.noint.

(grid.PU2 <- expand.grid(k.int, k.noint))
##   Var1 Var2
## 1    4    6
## 2    5    6
(zd <- which((ind1 %in% grid.PU2[,1] & ind2 %in% grid.PU2[,2]) |
             (ind2 %in% grid.PU2[,1] & ind1 %in% grid.PU2[,2])))
## integer(0)

# Check
ind[zd]
## named character(0)
ind.tmp[zd]
## named character(0)

Finally, we need to take care of the scale parameters. For the first part, we only require the index from k.noint, while for the second part we need the combinations of k.int.lam (indices of the scale parameters which have interactions with k) and k.noint.

(k.noint)
## [1] 6
(grid.PU.lam <- expand.grid(k.int.lam, k.noint))
##   Var1 Var2
## 1    2    6
## 2    3    6

Calculation of \(\mathbf S_k\)

The matrix \(\mathbf S_k\) is given by the formula

\[ \mathbf S_k = \mathbf P_k\mathbf R_k + (\mathbf P_k\mathbf R_k)^\top + \mathbf P_k\mathbf U_k + (\mathbf P_k\mathbf U_k)^\top. \]

As \(\mathbf P_k\mathbf R_k\) and \(\mathbf P_k\mathbf U_k\) are made up linearly of two-way matrix products of the kernel matrices which are stored in H2l, all we need is to add the right entries of H2l together (and not forgetting the respective scale parameters). The indices are given by the function indxFn().

indB <- indxFn(1)
indB$PRU  # = c(za, zc, zb, zd) i.e. index of Hl to sum together
## [1] 1 2
rbind(indB$PRU.lam1, indB$PRU.lam2)  # index of the lambdas to cross-product with Hl
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    0    0    0    2    3    2    3    2    3
## [2,]    2    3    6    2    2    3    3    6    6

For k=1, we calculate \(\mathbf S_1\) as

lambda.PRU <- c(rep(1, sum(indB$PRU.lam1 == 0)), lambda[indB$PRU.lam1])
lambda.PRU <- lambda.PRU * lambda[indB$PRU.lam2]
S <- Reduce("+", mapply("*", H2l[indB$PRU], lambda.PRU, SIMPLIFY = FALSE))

This is an efficient way to calculate \(\mathbf S_k\) by simply recalling the already multiplied matrices. In the EM algorithm, we would need to repeat this calculation for each k=1,...,l and also for each EM step t=1,2,....

Note that \(\mathbf S_k\) is calculated this way only if there are parsimonious interactions. When only a single scale paramters is used (e.g. using one.lam = TRUE), then \(\mathbf S_k = 0\). When no interactions are present, then \(\mathbf U_k = 0\), and we only need \(\mathbf P_k \mathbf R_k + (\mathbf P_k \mathbf R_k)^\top\). But this becomes easier because \(\mathbf P_k = \mathbf H_k\) as there are no Hadamard interactions.

Also realise that we never calculate \(\mathbf R_k\) and \(\mathbf U_k\) explicitly, because only \(\mathbf S_k\) is required in the closed form expression of \(\lambda_k^{(t+1)}\).

Efficient calculation of \(\mathbf P_k^2\)

When (parsimonious) interactions are present, \(\mathbf P_k^2\) is given by

This sum is made of two parts. The first is by adding up relevant squared kernel matrices and Hadamard products. We can collate these matrix products into a list of length h = l + m called Hsql. The second part comes from H2l as we have discussed above. Now it is a matter of summing up the right parts.

For k=1, the first part is getting the indices of the squared terms correctly. This is easy as we have already obtained this previously.

(Psq <- c(k, k.int))
## [1] 1 4 5
(Psq.lam <- k.int.lam)
## [1] 2 3

For the second part, we use combn() to generate all possible two-way combinations of c(k, k.int). This would give us the indices for the sums in the second part above. The corresponding scale paramters lambda are obtained the same way, but from all possible combinations of c(0, k.int.lam). The two columns of grid.Psq.lam give the index for which lambda needs to be multiplied. An entry of 0 means that only the other non-zero column entry is used, e.g. for grid.Psq.lam[1,], we multiply 1 * lambda[2]; for grid.Psq.lam[3,], we multiply lambda[2] * lambda[3].

(grid.Psq <- t(combn(c(k, k.int), 2)))
##      [,1] [,2]
## [1,]    1    4
## [2,]    1    5
## [3,]    4    5
(ze <- which((ind1 %in% grid.Psq[,1] & ind2 %in% grid.Psq[,2]) |
             (ind2 %in% grid.Psq[,1] & ind1 %in% grid.Psq[,2])))
## integer(0)

# Check
ind[ze]
## named character(0)
ind.tmp[ze]
## named character(0)

grid.Psq.lam <- NULL
if (length(k.int.lam) > 0) grid.Psq.lam <- t(combn(c(0, k.int.lam), 2))
grid.Psq.lam
##      [,1] [,2]
## [1,]    0    2
## [2,]    0    3
## [3,]    2    3

The code, which is found in the function for kernL(), is given by

# First part of sum
Psql[[k]] <<- Reduce("+", mapply("*", Hsql[indB$Psq],
                                 c(1, lambda[indB$Psq.lam] ^ 2),
                                 SIMPLIFY = FALSE))

# Second part of sum
lambda.P2 <- c(rep(1, sum(indB$P2.lam1 == 0)), lambda[indB$P2.lam1])
lambda.P2 <- lambda.P2 * lambda[indB$P2.lam2]
Psql[[k]] <<- Psql[[k]] + Reduce("+", mapply("*", H2l[indB$P2],
                                             lambda.P2,
                                             SIMPLIFY = FALSE))

Note that in cases where there are no interactions (or with non-parsimonious interactions), then \(\mathbf P_k^2 = \mathbf H_k^2\) and does not depend on \(\lambda_k\), thus can be calculated once and stored.