K <- 25L
D <- diag(1 + runif(K))
a <- runif(1)
W <- round(matrix(runif(K^2), K, K))
W[upper.tri(W, TRUE)] <- 0
W <- W + t(W)
Lambda <- Harman23.cor$cov[1:4,1:4]

mult <-
  function(D, a, W, Lambda, y) {
    temp <- Lambda %*% matrix(y, nrow = ncol(Lambda))
    # optimize this particularly if D - a * W is sparse!!!
    scalar <- t(y) %*% c(apply(D - a * W, 1, "%*%", t(temp)))
    return(scalar[1])
  }

all.equal(c( t(y) %*% kronecker(D - a * W, Lambda) %*% y ), 
          mult(D, a, W, Lambda, y))