R/sparse.R

In evgam: Generalised Additive Extreme Value Models

.joinSmooth <- function(lst, sparse = FALSE) {
  nms <- c("S", "first.para", "last.para", "rank", "null.space.dim", "id")
  nbi <- sapply(lst, function(x) x$nb)
  starts <- cumsum(c(0, nbi))
  if (length(lst) == 1) {
    lst <- lst[[1]]$smooth
  } else {
    lst <- lapply(lst, function(x) x$smooth)
    smooths <- sapply(lst, length) > 0
    for (i in which(smooths)) {
      for (j in seq_along(lst[[i]])) {
        lst[[i]][[j]]$first.para <- lst[[i]][[j]]$first.para + starts[i]
        lst[[i]][[j]]$last.para <- lst[[i]][[j]]$last.para + starts[i]
      }
    }
    lst <- unlist(lst, recursive=FALSE)
  }
  out <- lapply(lst, function(x) subset(x, names(x) %in% nms))
  lst2 <- list()
  nb <- tail(starts, 1)
  for (i in seq_along(lst)) {
    temp <- list()
    ind <- lst[[i]]$first.para:lst[[i]]$last.para
    for (j in seq_along(lst[[i]]$S)) {
      if (!sparse) {
        temp2 <- matrix(0, nb, nb)
      } else {
        temp2 <- Matrix::Matrix(0, nb, nb, sparse = TRUE)
      }
      temp2[ind, ind] <- lst[[i]]$S[[j]]
      temp[[j]] <- temp2
    }
    lst2[[i]] <- temp
  }
  lst2 <- unlist(lst2, recursive=FALSE)
  for (i in seq_along(out)) {
    if (length(out[[i]]$null.space.dim) == 0)
      out[[i]]$null.space.dim <- 0
  }
  if (sparse)
    lst2 <- lapply(lst2, as, 'CsparseMatrix')
  Sreps <- sapply(lst, function(x) length(x$S))
  possible <- seq_len(sum(Sreps))
  ids <- lapply(lapply(out, '[[', 'id'), as.integer)
  got <- as.integer(unlist(ids))
  possible <- possible[!possible %in% got]
  for (i in seq_along(ids)) {
    if (length(ids[[i]]) == 0) {
      ids[[i]] <- possible[1:Sreps[i]]
      possible <- possible[!possible %in% unlist(ids)]
    } else {
      if (length(ids[[i]]) < Sreps[i])
        ids[[i]] <- rep(ids[[i]], Sreps[i])
    }
  }
  ids <- unlist(ids)
  uids <- unique(ids)
  nsp <- length(uids)
  A <- matrix(0, nsp, length(ids))
  A[cbind(ids, 1:ncol(A))] <- 1
  attr(out, 'rho_rep') <- order(ids)
  labs <- unlist(lapply(lst, function(x) paste0(x$label, sapply(x$margin, '[[', 'label'))))
  attr(out, 'label') <- labs[attr(out, 'rho_rep')]
  attr(out, "Sl") <- lst2#[attr(out, 'rho_rep')]
  attr(out, "nb") <- nb
  attr(out, 'nsp') <- length(uids)
  attr(out, 'A') <- A
  out
}

.gH_dense <- function(x, likdata, sandwich=FALSE, deriv=2) {
  nX <- length(likdata$X)
  if (nX == 1) {
    out <- .gH1(x, likdata$X[[1]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
  } else {
    if (nX == 2) {
      out <- .gH2(x, likdata$X[[1]], likdata$X[[2]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
    } else {
      if (nX == 3) {
        out <- .gH3(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
      } else {
        if (nX == 4) { # added with evgam_0.1.2 (05/04/2020)
          out <- .gH4(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
        } else {
          if (nX == 5) { # added with evgam_0.1.5 (29/06/2021)
            out <- .gH5(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$X[[5]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
          } else {
            if (nX == 6) { # added with evgam_1.0.1 (31/08/2022)
              out <- .gH6(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$X[[5]], likdata$X[[6]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
            } else {
              stop("Number of model parameters not in {1, 2, 3, 4, 5, 6}")
            }
          }
        }
      }
    }
  }
  out
}

.gH_sparse <- function(x, likdata, sandwich=FALSE, deriv=2) {
  nX <- length(likdata$X)
  if (nX == 1) {
    out <- .gHsp1(x, likdata$X[[1]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
  } else {
    if (nX == 2) {
      out <- .gHsp2(x, likdata$X[[1]], likdata$X[[2]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
    } else {
      if (nX == 3) {
        out <- .gHsp3(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
      } else {
        if (nX == 4) { # added with evgam_0.1.2 (05/04/2020)
          out <- .gHsp4(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
        } else {
          if (nX == 5) { # added with evgam_0.1.5 (29/06/2021)
            out <- .gHsp5(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$X[[5]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
          } else {
            if (nX == 6) { # added with evgam_1.0.1 (31/08/2022)
              out <- .gHsp6(x, likdata$X[[1]], likdata$X[[2]], likdata$X[[3]], likdata$X[[4]], likdata$X[[5]], likdata$X[[6]], likdata$dupid, likdata$duplicate, as.integer(sandwich), deriv)
            } else {
              stop("Number of model parameters not in {1, 2, 3, 4, 5, 6}")
            }
          }
        }
      }
    }
  }
  out
}

.gH <- function(x, likdata, sandwich=FALSE, deriv=2) {
  if (!likdata$sparse) {
    out <- .gH_dense(x, likdata, sandwich, deriv)
  } else {
    out <- .gH_sparse(x, likdata, sandwich, deriv)
  }
  out
}

.gH.nopen.dense <- function(pars, likdata, likfns, sandwich=FALSE, deriv=2) {
  pars <- as.vector(likdata$compmode + likdata$CH %*% (as.vector(pars) - likdata$compmode))
  if ('sandwich' %in% names(formals(likfns$d120))) {
    temp <- likfns$d120(pars, likdata, sandwich)
  } else {
    temp <- likfns$d120(pars, likdata)
  }
  if (is.null(likdata$agg))
    temp <- .gH(temp, likdata, sandwich, deriv)
  temp[[1]] <- likdata$k * temp[[1]]
  temp[[1]] <- t(temp[[1]] %*% likdata$CH)
  if (deriv > 1) {
    temp[[2]] <- likdata$k * temp[[2]]
    temp[[2]] <- crossprod(likdata$CH, temp[[2]]) %*% likdata$CH
    attr(temp, "PP") <- temp[[2]] / norm(temp[[2]], "F")
  }
  temp
}

.gH.nopen.sparse <- function(pars, likdata, likfns, sandwich=FALSE, deriv=2) {
  pars <- as.vector(likdata$compmode + likdata$CH %*% (as.vector(pars) - likdata$compmode))
  if ('sandwich' %in% names(formals(likfns$d120))) {
    temp <- likfns$d120(pars, likdata, sandwich)
  } else {
    temp <- likfns$d120(pars, likdata)
  }
  if (is.null(likdata$agg))
    temp <- .gH(temp, likdata, sandwich, deriv)
  temp[[1]] <- likdata$k * temp[[1]]
  temp[[1]] <- as.matrix(temp[[1]] %*% likdata$CH)
  temp[[1]] <- t(temp[[1]])
  if (deriv > 1) {
    temp[[2]] <- likdata$k * temp[[2]]
    temp[[2]] <- Matrix::crossprod(likdata$CH, temp[[2]]) %*% likdata$CH
    # attr(temp, "PP") <- temp[[2]] / Matrix::norm(temp[[2]], "F")
  }
  temp
}

.gH.nopen <- function(pars, likdata, likfns, sandwich=FALSE, deriv=2) {
  if (!likdata$sparse) {
    out <- .gH.nopen.dense(pars, likdata, likfns, sandwich, deriv)
  } else {
    out <- .gH.nopen.sparse(pars, likdata, likfns, sandwich, deriv)
  }
out
}

.gH.pen.dense <- function(pars, likdata, likfns, deriv=2) {
  temp <- .gH.nopen(pars, likdata, likfns, deriv=deriv)
  temp[[1]] <- temp[[1]] + crossprod(pars, likdata$S)[1, ]
  if (deriv > 1) {
    if (norm(likdata$S, "F") > 0) {
      attr(temp, "PP") <- attr(temp, "PP") + likdata$S / norm(likdata$S, "F")
    } else {
      attr(temp, "PP") <- attr(temp, "PP") + likdata$S
    }
    temp[[2]] <- temp[[2]] + likdata$S
  }
  temp
}

.gH.pen.sparse <- function(pars, likdata, likfns, deriv=2) {
  temp <- .gH.nopen(pars, likdata, likfns, deriv=deriv)
  # if (!likdata$sparse) {
  #   temp[[1]] <- temp[[1]] + base::crossprod(pars, likdata$S)[1, ]
  # } else {
    temp[[1]] <- temp[[1]] + Matrix::crossprod(pars, likdata$S)[1, ]
  # }
  if (deriv > 1) {
    # attr(temp, "PP") <- attr(temp, "PP") + likdata$S / Matrix::norm(likdata$S, "F")
    temp[[2]] <- temp[[2]] + likdata$S
  }
  temp
}

.gH.pen <- function(pars, likdata, likfns, deriv=2) {
  if (likdata$sparse) {
    out <- .gH.pen.sparse(pars, likdata, likfns, deriv)
  } else {
    out <- .gH.pen.dense(pars, likdata, likfns, deriv)
  }
  out
}

.search.dir.dense <- function(g, H, kept=NULL) {
  if (is.null(kept))
    kept <- !logical(length(g))
  H0 <- H
  g0 <- g
  g <- g[kept]
  H <- H[kept, kept, drop=FALSE]
  if (any(!is.finite(g)))
    stop("Some gradient non-finite")
  if (any(!is.finite(H)))
    stop("Some Hessian non-finite")
  H2 <- .precondition(H)
  H2 <- .perturb(H2)
  R <- attr(H2, "chol")
  d <- attr(H2, "d")
  out <- numeric(length(kept))
  piv <- ipiv <- attr(R, "pivot")
  ipiv[piv] <- seq_len(length(piv))
  out[kept] <- .precond_solve_dense(attr(H2, "chol"), g)
  g0[!kept] <- 0
  attr(out, "gradient") <- g0
  attr(out, "Hessian") <- H0
  attr(out, "cholHessian") <- R
  attr(out, "rankHessian") <- attr(H2, 'rank0')
  out
}

.search.dir.sparse <- function(g, H, kept=NULL) {
  if (is.null(kept))
    kept <- !logical(length(g))
  H0 <- H
  g0 <- g
  g <- g[kept]
  H <- H[kept, kept, drop=FALSE]
  if (any(!is.finite(g)))
    stop("Some gradient non-finite")
  # if (any(!is.finite(H)))
  #   stop("Some Hessian non-finite")
  out <- numeric(length(kept))
  H2 <- .precondition(H)
  H2 <- .perturb(H2)
  attr(attr(H2, 'chol'), 'H2') <- H2
  out[kept] <- .precond_solve_sparse(attr(H2, 'chol'), g)
  g0[!kept] <- 0
  attr(out, "gradient") <- g0
  attr(out, "Hessian") <- H2
  attr(out, "cholHessian") <- attr(H2, 'chol')
  attr(out, "rankHessian") <- attr(H2, 'rank0')
  out
}

.search.dir <- function(g, H, kept=NULL) {
  if (inherits(H, 'matrix')) {
    out <- .search.dir.dense(g, H, kept)
  } else {
    out <- .search.dir.sparse(g, H, kept)
  }
  out
}

.Hdata_dense <- function(H) {
  out <- list(H0=H)
  H2 <- .precondition(H)
  H2 <- .perturb(H2)
  out$H <- H2
  out$dH <- attr(H2, "d")
  out$cH <- attr(H2, "chol")
  # out$iH <- MASS::ginv(out$H0)
  out$iH <- .precond_solve(out$cH)
  out$kept <- !logical(nrow(H))
  out
}

.Hdata_sparse <- function(H) {
  out <- list(H0=H)
  H2 <- .precondition(H)
  H2 <- .perturb(H2)
  out$H <- H2
  out$dH <- attr(H2, "d")
  out$cH <- attr(H2, "chol")
  # out$iH <- Matrix::solve(out$H0)
  # t1 <- Matrix::tcrossprod(Matrix::solve(out$cH, Matrix::Diagonal(nrow(H)))) * Matrix::tcrossprod(out$dH)
  # range(t1 - out$iH)
  # t2 <- Matrix::tcrossprod(Matrix::solve(out$cH, Matrix::Diagonal(nrow(H))))
  # range(t2 - out$iH)
  # out$dH * Matrix::t(Matrix::chol2inv(out$cH) * out$dH)# out$iH <- Matrix::solve(out$cH, Matrix::Diagonal(out$dH), transpose=TRUE))
  # out$iH <- out$dH * Matrix::t(Matrix::chol2inv(out$cH) * out$dH)# out$iH <- Matrix::solve(out$cH, Matrix::Diagonal(out$dH), transpose=TRUE))
  out$iH <- NULL#out$dH * Matrix::t(Matrix::chol2inv(out$cH) * out$dH)# out$iH <- Matrix::solve(out$cH, Matrix::Diagonal(out$dH), transpose=TRUE))
  out$kept <- !logical(nrow(H))
  out
}

.Hdata <- function(H) {
  if (inherits(H, 'matrix')) {
    out <- .Hdata_dense(H)
  } else {
    out <- .Hdata_sparse(H)
  }
  out
}

.perturb_dense <- function(A) {
  d <- attr(A, "d")
  eps <- 1e-12
  cholA <- suppressWarnings(chol(A, pivot=TRUE))
  rk <- rk0 <- attr(cholA, "rank")
  while(rk < nrow(A)) {
    diag(A) <- diag(A) + eps
    cholA <- suppressWarnings(chol(A, pivot=TRUE))
    rk <- attr(cholA, "rank")
    eps <- 1e2 * eps
  }
  attr(cholA, 'd') <- d
  attr(A, "chol") <- cholA
  attr(A, "rank") <- rk
  attr(A, "rank0") <- rk0
  return(A)
}

.perturb_sparse <- function(A) {
  d <- attr(A, 'd')
  eps <- 1e-12
  cholA <- suppressWarnings(Matrix::Cholesky(A, LDL = TRUE))
  rk0 <- sum(Matrix::diag(cholA) > 1e-8)#length(d)#Matrix::rankMatrix(A, method = 'qr')
  # cholA <- suppressWarnings(try(Matrix::Cholesky(A, LDL = FALSE)))
  # cholA <- suppressWarnings(try(Matrix::chol(A), silent = TRUE))
  while(any(Matrix::diag(cholA) <= 1e-8)) {
  # while(inherits(cholA, 'try-error')) {
    Matrix::diag(A) <- Matrix::diag(A) + eps
    cholA <- suppressWarnings(Matrix::Cholesky(A, LDL = TRUE))
    # cholA <- suppressWarnings(try(Matrix::Cholesky(A, LDL = FALSE)))
    # cholA <- suppressWarnings(try(Matrix::chol(A), silent = TRUE))
    #rk <- Matrix::rankMatrix(A, method = 'qr')
    eps <- 1e2 * eps
  }
  # cholA <- Matrix::Cholesky(A)
  attr(cholA, 'd') <- d
  attr(A, "chol") <- cholA
  # attr(A, "chol2") <- Matrix::Cholesky(A)
  attr(A, "rank") <- length(d)#Matrix::rankMatrix(A, method = 'qr')
  attr(A, "rank0") <- rk0
  return(A)
}

.perturb <- function(A) {
  if (inherits(A, 'matrix')) {
    out <- .perturb_dense(A)
  } else {
    out <- .perturb_sparse(A)
  }
  out
}

.precondition_dense <- function(H) {
  d <- 1/sqrt(abs(pmax(diag(H), 1e-8)))
  out <- d * t(t(H) * d)
  attr(out, "d") <- d
  out
}

.precondition_sparse <- function(H) {
  d <- 1/sqrt(abs(pmax(Matrix::diag(H), 1e-8)))
  out <- d * Matrix::t(Matrix::t(H) * d)
  attr(out, "d") <- d
  out
}

.precondition <- function(H) {
  if (inherits(H, 'matrix')) {
    out <- .precondition_dense(H)
  } else {
    out <- .precondition_sparse(H)
  }
  out
}

# .precond_solve_dense <- function(L, x) {
#   d <- attr(L, "d")
#   x <- d * as.matrix(x)
#   piv <- ipiv <- attr(L, "pivot")
#   ipiv[piv] <- seq_len(nrow(L))
#   d * backsolve(L, backsolve(L, x[piv, , drop=FALSE], upper.tri=TRUE, transpose=TRUE))[ipiv, , drop=FALSE]
# }

.precond_solve_dense <- function(L, x) {
  d <- attr(L, "d")
  piv <- ipiv <- attr(L, "pivot")
  ipiv[piv] <- seq_len(nrow(L))
  if (is.null(x)) {
    x <- diag(length(d))
    out <- tcrossprod(backsolve(L, x))[ipiv, ipiv] * tcrossprod(d)
  } else {
    x <- d * as.matrix(x)
    out <- d * backsolve(L, backsolve(L, x[piv, , drop=FALSE], upper.tri=TRUE, transpose=TRUE))[ipiv, , drop=FALSE]
  }
  out
}

.precond_solve_sparse <- function(L, x) {
  d <- attr(L, 'd')
  if (is.null(x)) {
    out <- Matrix::tcrossprod(Matrix::solve(L, Matrix::Diagonal(nrow(L)))) * Matrix::tcrossprod(d)
  } else {
    # out <- d * Matrix::solve(L, Matrix::solve(Matrix::t(L), x * d))
    out <- d * Matrix::solve(L, x * d)
  }
  out
}

.precond_solve <- function(L, x = NULL) {
  if (is.null(attr(L, 'd')))
    attr(L, 'd') <- rep(1, nrow(L))
  if (inherits(L, 'matrix')) {
    out <- .precond_solve_dense(L, x)
  } else {
    out <- .precond_solve_sparse(L, x)
  }
  out
}

.d0logdetH_dense <- function(x) {
  if (is.null(x$cholHessian)) {
    out <- as.vector(determinant(x$Hessian)[[1]])
  } else {
    cH <- x$cholHessian
    out <- 2 * sum(log(diag(cH)))
    out <- out - 2 * sum(log(attr(cH, "d")))
  }
  list(d0 = out)
}

.d0logdetH_sparse <- function(x) {
  if (is.null(x$cholHessian)) {
    out <- as.vector(Matrix::determinant(x$Hessian)[[1]])
  } else {
    cH <- x$cholHessian
    out <- sum(log(Matrix::diag(cH)))
    out <- out - 2 * sum(log(attr(cH, "d")))
  }
  list(d0 = out)
}

.d0logdetH <- function(x, sparse) {
  if (!sparse) {
    out <- .d0logdetH_dense(x)
  } else {
    out <- .d0logdetH_sparse(x)
  }
  out
}

.d1beta <- function(lsp, beta, spSl, H) {
  out <- list(d0 = beta)
  spSlb <- sapply(spSl, function(x) x %*% beta)
  if (inherits(spSlb, 'list'))
    spSlb <- sapply(spSlb, as.vector)
  out$d1 <- -.precond_solve(H$cH, spSlb)
  out$spSlb <- spSlb
  out
}

.choltr <- function(x, b) {
  if (inherits(x, 'matrix')) {
    sum(diag(.precond_solve(x, b)))
  } else {
    sum(Matrix::diag(.precond_solve(x, b)))
  }
}

.d1H0_diag <- function(dbeta, likdata, likfns, H) {

  X <- likdata$X
  idpars <- likdata$idpars
  CH <- likdata$CH

  nb <- nrow(dbeta$d1)
  nsp <- ncol(dbeta$d1)
  nX <- length(X)
  n <- nrow(X[[1]])

  ind <- .indices(nX)

  beta <- likdata$compmode + likdata$CH %*% (dbeta$d0 - likdata$compmode)

  GH <- likdata$k * likfns$d340(beta, likdata)

  dbeta$d1 <- CH %*% dbeta$d1

  d1eta <- lapply(seq_len(nX), function(i) X[[i]] %*% dbeta$d1[idpars == i, , drop=FALSE])

  trd1H <- numeric(nsp)
  
  for (l in 1:nsp) {

    if (!likdata$sparse) {
      d1Hl <- matrix(0, nb, nb)
    } else {
      d1Hl <- Matrix::Matrix(0, nrow = nb, ncol = nb, sparse = TRUE)
    }

    for (i in 1:nX) {
      for (j in 1:nX) {
        v <- numeric(n)
        for (k in 1:nX)
          v <- v + GH[,ind$i3[i, j, k]] * d1eta[[k]][, l]
        if (!likdata$sparse) {
          d1Hl[idpars == i, idpars == j] <- crossprod(X[[i]], X[[j]] * v)
        } else {
          d1Hl[idpars == i, idpars == j] <- Matrix::crossprod(X[[i]], X[[j]] * v)
        }
      }  
    }

    trd1H[l] <- .choltr(H$cH, d1Hl)
  
  }

  list(d1 = trd1H)

}

.d1H0 <- function(dbeta, likdata, likfns, H) {
  
  X <- likdata$X
  idpars <- likdata$idpars
  CH <- likdata$CH
  
  nb <- nrow(dbeta$d1)
  nsp <- ncol(dbeta$d1)
  nX <- length(X)
  n <- nrow(X[[1]])
  
  ind <- .indices(nX)
  
  beta <- likdata$compmode + likdata$CH %*% (dbeta$d0 - likdata$compmode)
  
  GH <- likdata$k * likfns$d340(beta, likdata)
  
  dbeta$d1 <- CH %*% dbeta$d1
  
  d1eta <- lapply(seq_len(nX), function(i) X[[i]] %*% dbeta$d1[idpars == i, , drop=FALSE])
  
  trd1H <- numeric(nsp)
  
  d1H <- list()
  
  for (l in 1:nsp) {
    
    if (!likdata$sparse) {
      d1Hl <- matrix(0, nb, nb)
    } else {
      d1Hl <- Matrix::Matrix(0, nrow = nb, ncol = nb, sparse = TRUE)
    }
    
    for (i in 1:nX) {
      for (j in 1:nX) {
        v <- numeric(n)
        for (k in 1:nX)
          v <- v + GH[,ind$i3[i, j, k]] * d1eta[[k]][, l]
        if (!likdata$sparse) {
          d1Hl[idpars == i, idpars == j] <- crossprod(X[[i]], X[[j]] * v)
        } else {
          d1Hl[idpars == i, idpars == j] <- Matrix::crossprod(X[[i]], X[[j]] * v)
        }
      }  
    }
    
    d1H[[l]] <- d1Hl
    
  }
  
  list(d1 = d1H)
  
}

.d1logdetH <- function(dbeta, likdata, likfns, spSl, H) {
  d1 <- .d1H0_diag(dbeta, likdata, likfns, H)$d1
  d1 <- d1 + sapply(spSl, function(x) .choltr(H$cH, x))
  list(d1 = d1, dbeta = dbeta)
}

.d2logdetH <- function(dbeta, likdata, likfns, spSl, H) {
  d1 <- .d1H0(dbeta, likdata, likfns)$d1
  out <- .d2H0_diag(dbeta, likdata, likfns, H)$d2
  nsp <- length(d1)
  for (l in 1:nsp) {
    d1[[l]] <- d1[[l]] + spSl[[l]]
    out[l, l] <- out[l, l] + sum(diag(.precond_solve(H$cH, spSl[[l]])))
    d1[[l]] <- .precond_solve(H$cH, d1[[l]])
  }
  for (l in 1:nsp) {
    for (m in 1:nsp) {
      out[l, m] <- sum(diag(d1[[l]] %*% d1[[m]])) - out[l, m]
    }
  }
  list(d2 = out)
}

.d12logdetH <- function(dbeta, likdata, likfns, spSl, H) {
  d1 <- .d1H0(dbeta, likdata, likfns)$d1
  out <- .d2H0_eigen(dbeta, likdata, likfns, H)$d2
  nsp <- length(d1)
  for (l in 1:nsp) {
    d1[[l]] <- d1[[l]] + spSl[[l]]
    out[l, l] <- out[l, l] + sum(diag(.precond_solve(H$cH, spSl[[l]])))
    d1[[l]] <- .precond_solve(H$cH, d1[[l]])
  }
  for (l in 1:nsp) {
    for (m in 1:nsp) {
      out[l, m] <- out[l, m] - sum(t(d1[[l]]) * d1[[m]])#sum(diag(d1[[l]] %*% d1[[m]]))
    }
  }
  list(d1 = sapply(d1, function(x) sum(diag(x))), d2 = out)
}

.d12logdetH_diag <- function(dbeta, likdata, likfns, spSl, H) {
  d1 <- .d1H0(dbeta, likdata, likfns)$d1
  out <- .d2H0_diag(dbeta, likdata, likfns, H)$d2
  nsp <- length(d1)
  for (l in 1:nsp) {
    d1[[l]] <- d1[[l]] + spSl[[l]]
    out[l] <- out[l] + sum(diag(.precond_solve(H$cH, spSl[[l]])))
    d1[[l]] <- .precond_solve(H$cH, d1[[l]])
    out[l] <- out[l] - sum(t(d1[[l]]) * d1[[l]])
  }
  list(d1 = sapply(d1, function(x) sum(diag(x))), d2 = out)
}

.pivchol_rmvn_sparse <- function(n, mu, Sig) {
  Ch <- Matrix::Cholesky(Sig, permutation = TRUE, LDL = FALSE)
  P <- as(Ch, "pMatrix")
  L <- as(Ch, "Matrix")
  Z <- matrix(rnorm(n * ncol(Sig)), ncol(Sig))
  as.matrix(Matrix::t(P) %*% L %*% Z) + as.vector(mu)
}