The M-step involving the scale parameters \(\boldsymbol\lambda\) in the EM algorithm for I-prior models can be found in closed-form in the following situations:
lambda being used.For any other models such as ones involving squared terms and
three-way interactions, the M-step can still be solved using numerical
methods such as a downhill simplex method (in R, we use
optim(method = "Nelder-Mead"). Examples:
| Model | parsm |
length(lambda) |
Closed form? |
|---|---|---|---|
y ~ x1 + x2 + x3 |
TRUE |
3 | Yes |
y ~ (x1 + x2 + x3) |
TRUE |
1 | Yes |
y ~ x1 + x2 + x1:x2 |
TRUE |
2 | Yes |
y ~ x1 + x2 + x1:x2 |
FALSE |
3 | Yes |
y ~ (x1 * x2 * x3) |
TRUE |
1 | Yes |
y ~ x1 * x2 * x3 |
FALSE |
7 | Yes |
y ~ x1 + I(x1 ^ 2) |
FALSE |
2 | Yes |
y ~ x1 * x2 * x3 |
TRUE |
3 | No |
y ~ x1 + I(x1 ^ 2) |
TRUE |
1 | No |
In short, the most complex model for which closed-form for
lambda exists is the parsimonious two-way interaction
model. We describe this below.
Assume there are \(p\) covariates, and each of the \(p\) kernel matrices \(\mathbf H_1, \dots, \mathbf H_p\) are calculated using the appropriate kernels, depending on whether the data is continuous or nominal and what effect is desired. Let the number of unique scale parameters be \(l \leq p\). \(l\) could be less than \(p\), which implies that some of the covariates share a scale parameter. For a group of such covariates, the kernel matrix is simply the sum of each kernel matrix, and thus the kernel matrices can be indexed from \(1, \dots, l\) as well. Otherwise \(l = p\).
If two-way interactions are present between any \(k,j \in \{1,\dots,l\}\), then these are also calculated as \(\mathbf H_{kj} = \mathbf H_k \circ \mathbf H_j\) (the Hadamard product). In general, the scaled kernel matrix looks like \[ \mathbf H_{\lambda} = \sum_{k=1}^l \lambda_k \mathbf H_k + \sum_{k,j\in M} \lambda_k\lambda_j \mathbf H_{kj} \] where the set \(M\) is the index of all two way interaction terms between the \(p\) covariates, i.e. \(M\) \(=\) \(\{(k,j):\) \(k \text{ interacts with } j,\) \(\text{ and }\) \(k < j,\) \(\\\) \(\forall k,j=1,\dots,l \}\). Let the number of two-way interactions be \(m=|M|\). The total number of scale parameters is equal to \(q=l+m\) when there are non-parsimonious interactions present, otherwise it is \(q=l\). The non-parsimonious method of interactions assigns a new scale parameter for each of the Hadamard products of interacting kernel matrices. In comparison, the parsimonious method simply multiplies the corresponding scale parameters together.
For a particular \(\lambda_k\), \(k \in \{1,\dots,q\}\), we partition the sum of the kernel matrix into parts which involve \(\lambda_k\) and parts which do not: \[ \begin{aligned} \mathbf H_\lambda &= \overbrace{\lambda_k \left( \mathbf H_k + \sum_{j\in M}\lambda_j\mathbf H_{kj} \right) \vphantom{\mathop{\sum_{j=1}^p}_{j\neq k}}}^{\lambda_k\text{ is here}} + \overbrace{{\mathop{\sum_{j=1}^l}_{j\neq k}\lambda_j \mathbf H_j} + {\mathop{\sum_{k',j \in M}}_{k'\neq k}\lambda_{k'}\lambda_j \mathbf H_{k'j}}}^{\text{no $\lambda_k$ here}} \\ &= \lambda_k\mathbf {P_k} + {\mathbf R_k} + {\mathbf U_k}. \end{aligned} \]
\(\mathbf P_k\) is the kernel matrix \(\mathbf H_k\) plus the sum-product of the interaction kernel matrices with the scale parameters relating to covariate \(k\), i.e. \(\sum_j \lambda_j\mathbf H_{kj}\). \(\mathbf R_k\) is the sum-product of the kernel matrices and scale parameters excluding \(\lambda_k\mathbf H_k\). \(\mathbf U_k\) is the sum of the interaction cross-product terms excluding those relating to covariate \(k\). Thus, the squared kernel matrix is \[ \begin{align} \mathbf H_{\lambda}^2 =& \ \lambda_k^2\mathbf P_k^2 + \lambda_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) \nonumber \\ & + \mathbf R_k^2 + \mathbf U_k^2 + \mathbf R_k\mathbf U_k + \mathbf U_k\mathbf R_k. \end{align} \]
The closed-form solution for the scale parameters in the M-step at iteration \(t\) is \[ \lambda_k^{(t+1)} = \frac{(\mathbf y - \hat{\boldsymbol\alpha})^\top \mathbf P_k \tilde{\mathbf w}^{(t)} - \frac{1}{2} \text{tr} \left[ \mathbf S_k \tilde{\mathbf W}^{(t)} \right]}{\text{tr} \left[ \mathbf P_k^2 \tilde{\mathbf W}^{(t)} \right]} \] where we have defined \(\mathbf S_k = \mathbf P_k\mathbf R_k + \mathbf R_k\mathbf P_k + \mathbf P_k\mathbf U_k + \mathbf U_k\mathbf P_k\), for each \(k = 1,\dots,l\).
For most cases, \(\mathbf P_k\) and \(\mathbf S_k\) only depend on the kernel matrices and not on the scale parameters, so can be calculated once and stored for efficiency. Further, \(\mathbf U_k\) equals zero for most cases except in the parsimonious multiple scale parameter case thus simplifying calculations. In fact, we can avoid the expensive matrix multiplications involved in evaluating \(\mathbf P_k\), its square, and \(\mathbf S_k\), by storing all possible square and two-way multiplications of the kernel matrices \(\mathbf H_1, \dots, \mathbf H_l\) as the relevant calculation of the M-step simply involves a sum-product of these kernel matrices.
intr is always a 2 x m matrix indexing all
the m two-way interactions in the model.
h is the length of the kernel matrix. If there are
p variables, and m two-way interactions, then
h contains the p kernel matrices, and
m Hadamard products between the kernel matrices according
to the intr indices. Thus h = p + m,
regardless of parsimonious or non-parsimonious interactions.
ind1 and ind2 together give the index of
all possible two-way interactions. In a h x h matrix, these
are the row (ind1) and column (ind2) indices
of the upper triangular matrix excluding the diagonal entries.
The goal is to efficiently compute a list of length h
which contains the H.mat_i ^ 2 for
i = 1,...,h. Incidentally, we have used q to
denote the expanded lambda length which includes
interactions and higher order terms, so q=l+m.
indxFn <- function(k) {
# Indexer helper function used to create indices for H2l. Note: intr, ind1 and
# ind2 are created in kernL().
ind.int1 <- intr[1, ] == k; ind.int2 <- intr[2, ] == k # locating var/kernel matrix
ind.int <- which(ind.int1 | ind.int2) # of interactions (out of 1:no.int)
k.int <- ind.int + p # which kernel matrix has interactions involving k
k.int.lam <- c(intr[1, ][ind.int2], intr[2, ][ind.int1]) # which has
# interaction with k?
nok <- (1:p)[-k] # all variables excluding k
k.noint <- which(!(ind.int1 | ind.int2)) + p # the opposite of k.int
# P.mat %*% R.mat + R.mat %*% P.mat indices ----------------------------------
za <- which((ind1 %in% k & ind2 %in% nok) | (ind2 %in% k & ind1 %in% nok))
grid.PR <- expand.grid(k.int, nok)
zb <- which((ind1 %in% grid.PR[,1] & ind2 %in% grid.PR[,2]) |
(ind2 %in% grid.PR[,1] & ind1 %in% grid.PR[,2]))
grid.PR.lam <- expand.grid(k.int.lam, nok)
# P.mat %*% U.mat + U.mat %*% P.mat indices ----------------------------------
grid.PU1 <- expand.grid(k, k.noint)
zc <- which((ind1 %in% grid.PU1[,1] & ind2 %in% grid.PU1[,2]) |
(ind2 %in% grid.PU1[,1] & ind1 %in% grid.PU1[,2]))
grid.PU2 <- expand.grid(k.int, k.noint)
zd <- apply(grid.PU2, 1, findH2, ind1 = ind1, ind2 = ind2)
grid.PU.lam <- expand.grid(k.int.lam, k.noint)
# P.mat %*% P.mat indices ----------------------------------------------------
grid.Psq <- t(combn(c(k, k.int), 2))
ze <- apply(grid.Psq, 1, findH2, ind1 = ind1, ind2 = ind2 )
grid.Psq.lam <- NULL
if (length(k.int.lam) > 0) grid.Psq.lam <- t(combn(c(0, k.int.lam), 2))
list(
k.int = k.int,
k.int.lam = k.int.lam,
PRU = c(za,zc,zb,zd),
PRU.lam1 = c(rep(0, length(nok) + length(k.noint)),
grid.PR.lam[,1], grid.PU.lam[,1]),
PRU.lam2 = c(nok, k.noint, grid.PR.lam[,2], grid.PU.lam[,2]),
Psq = c(k, k.int),
Psq.lam = k.int.lam,
P2 = ze,
P2.lam1 = grid.Psq.lam[,1],
P2.lam2 = grid.Psq.lam[,2]
)
}
findH2 <- function(z, ind1, ind2){
# This function finds position of H2 (cross-product terms of H). Used in
# indxFn()
x <- z[1]; y <- z[2]
which((ind1 == x & ind2 == y) | (ind2 == x & ind1 == y))
}Regression with 3 covariates and two-way interactions between all 3
covariates. Here, p = 3, l = 3 and
h = q = 6. In full, the index of the H.mat is
c(1, 2, 3, 1:2, 1:3, 2:3).
(mod <- kernL(stack.loss ~ . ^ 2, data = stackloss))
## Sample size: 21
## No. of covariates: 3
## Object size: 251.1 kB
##
## Kernel matrices:
## 1 linear [1:21, 1:21] 383 383 285.2 30.8 30.8 ...
## 2 linear [1:21, 1:21] 34.87 34.87 23.06 17.15 5.34 ...
## 3 linear [1:21, 1:21] 7.37 4.65 10.08 1.94 1.94 ...
## 4 linear x linear [1:21, 1:21] 13355 13355 6575 528 164 ...
## 5 linear x linear [1:21, 1:21] 2822 1782.3 2875.1 59.6 59.6 ...
## 6 linear x linear [1:21, 1:21] 256.9 162.2 232.4 33.3 10.4 ...
##
## Hyperparameters to estimate:
## lambda[1], lambda[2], lambda[3], psi
##
## Estimation methods available:
## direct, em, canonical, mixed, fixed
p <- 3The index of all two-way interactions are obtained by the kernel
loader function. It is contained in model$intr. The
following shows the indices of the variables which have two-way
interactions. For example, variable 1 interacts with variable 2,
variable 1 with 3 and finally 2 with 3. This matrix will always have 2
rows, and columns equal to m.
Next, we list out the indices of all possible two-way terms. This is used to compute the cross-product when multilpying out the square of a sum of matrices.
h <- length(mod$Hl)
z <- 1:h
(ind1 <- rep(z, times = (length(z) - 1):0))
## [1] 1 1 2
(ind2 <- unlist(lapply(2:length(z), function(x) c(NA, z)[-(0:x)])))
## [1] 2 3 3HlAll 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 2Which 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.
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?
Next I simply call nok the indices of all variables
excluding k. sidenote: I am beginning to think that
k.int.lam == nok.
Finally, these are the indices of the Hadamard products which do not
involve variable k.
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:
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 3The required product \(\mathbf P_k \mathbf R_k\) is then, in a manner of speaking,
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.
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 6For 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)}\).
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.
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 3The 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.