1
#' Check of KKT Conditions for Linear Mixed Model
2
#'
3
#' @description This function checks the KKT conditions
4
#' @param x rotated x. Should be U^T X, where U is the matrix of eigenvectors
5
#'   and X contains the first column of ones for the intercept. x should be a
6
#'   mtrix of dimension n x (p+1). These are outputted by the constructor
7
#'   functions. See \code{\link{ggmix_data_object}} for details
8
#' @param y rotated y. Should be U^T Y, where U is the matrix of eigenvectors
9
#'   and Y is the response.
10
#' @param eigenvalues non-zero eigenvalues of the kinship matrix, or the square
11
#'   of the singular values of the matrix used to construct the kinship matrix
12
#' @inheritParams ggmix
13
#' @inheritParams logliklasso
14
#' @param tol.kkt Tolerance for determining if an entry of the subgradient is
15
#'   zero
16
#' @rdname kkt_check
17
#' @note \code{grr_sigma2} and \code{grr_beta0} are functions for the gradient
18
#'   of sigma2 and beta0, respectively
19
kkt_check <- function(eta, sigma2, beta, eigenvalues, x, y, nt,
20
                                    lambda, tol.kkt = 1e-3) {
21

22
  # eta = eta_next; sigma2 = sigma2_next;
23
  # beta = beta_next;
24
  # # beta = as.matrix(c(beta0_next,rep(0,p)))
25
  # eigenvalues = Lambda;
26
  # x = utx;
27
  # y = uty; nt = n;
28
  # lambda = lambda; tol.kkt = 0.000001
29
  # =======================
30

31 4
  di <- 1 + eta * (eigenvalues - 1)
32 4
  wi <- (1 / sigma2) * (1 / di)
33

34

35
  # scale the weights to sum to nvars
36 4
  wi_scaled <- as.vector(wi) / sum(as.vector(wi)) * nt
37
  # wi_mat <- diag(wi_scaled)
38

39
  # KKT for beta0
40
  # grr_beta0(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt)
41
  # kkt_beta0 <- abs(grr_beta0(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt)) < tol.kkt
42

43 4
  kkt_beta0 <- grr_beta0(
44 4
    eta = eta, sigma2 = sigma2, beta = beta,
45 4
    eigenvalues = eigenvalues, x = x, y = y, nt = nt
46
  )
47

48
  # KKT for eta
49
  # gr_eta_lasso_fullrank(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt)
50
  # kkt_eta <- gr_eta_lasso_fullrank(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) < tol.kkt
51 4
  kkt_eta <- gr_eta_lasso_fullrank(
52 4
    eta = eta, sigma2 = sigma2, beta = beta,
53 4
    eigenvalues = eigenvalues, x = x, y = y, nt = nt
54
  )
55

56
  # KKT for sigma2
57
  # grr_sigma2(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt)
58
  # kkt_sigma2 <- grr_sigma2(eta = eta, sigma2 = sigma2, beta = beta, eigenvalues = eigenvalues, x = x, y = y, nt = nt) < tol.kkt
59 4
  kkt_sigma2 <- grr_sigma2(
60 4
    eta = eta, sigma2 = sigma2, beta = beta,
61 4
    eigenvalues = eigenvalues, x = x, y = y, nt = nt
62
  )
63

64
  # KKT for beta
65
  # g0 <- (1 / nt) * crossprod(x[,-1, drop = F], wi) %*% (y - x %*% beta) / (colSums(sweep(x[,-1, drop = F]^2, MARGIN = 1, wi_vec, '*')))
66

67
  # KKT for beta
68
  # g0 <- (1 / sum(wi_scaled)) * crossprod(x[,-1, drop = F] * wi_scaled, (y - x %*% beta ))
69 4
  g0 <- (1 / sum(wi)) * crossprod(x[, -1, drop = F] * wi, (y - x %*% beta))
70

71
  # this gives same result as g0
72
  # g1 <- colSums((1 / nt) * sweep(sweep(x[,-1], MARGIN = 1, wi_vec, '*'), MARGIN = 1, drop((y - x %*% beta)),'*'))
73
  # plot(drop(g0),g1)
74
  # abline(a=0,b=1)
75

76 4
  g <- g0 - lambda * sign(beta[-1])
77

78
  # this is for when beta=0 and should be between -1 and 1
79 4
  gg <- g0 / lambda
80

81
  # which of the betas are non-zero
82 4
  oo <- abs(beta[-1]) > 0
83

84
  # if all betas are 0 then set to TRUE, else abs(g[oo]) will give error since 'oo' is all FALSE
85
  # kkt_beta_nonzero <- if (all(!oo)) TRUE else max(abs(g[oo])) < tol.kkt
86
  # kkt_beta_subgr1 <- min(gg[!oo]) > -1
87
  # kkt_beta_subgr2 <- max(gg[!oo]) < 1
88

89 4
  kkt_beta_nonzero <- if (all(!oo)) 0 else sum(abs(g[oo]) > tol.kkt)
90 4
  kkt_beta_subgr <- sum(abs(gg[!oo]) > 1)
91
  # if (sum(abs(g[oo]) > tol.kkt) > 0) plot(abs(g[oo]))
92

93 4
  return(c(
94 4
    kkt_beta0 = kkt_beta0,
95 4
    kkt_eta = kkt_eta,
96 4
    kkt_sigma2 = kkt_sigma2,
97 4
    kkt_beta_nonzero = kkt_beta_nonzero,
98 4
    kkt_beta_subgr = kkt_beta_subgr
99
  ))
100
}
101

102
#' @rdname kkt_check
103
grr_sigma2 <- function(eta, sigma2, beta, eigenvalues, x, y, nt) {
104 4
  di <- 1 + eta * (eigenvalues - 1)
105

106 4
  sigma2 - (1 / nt) * sum((((y - x %*% beta)^2) / di))
107
}
108

109
#' @rdname kkt_check
110
grr_beta0 <- function(eta, sigma2, beta, eigenvalues, x, y, nt) {
111 4
  di <- 1 + eta * (eigenvalues - 1)
112 4
  wi <- (1 / sigma2) * diag(1 / di)
113

114 4
  as.numeric(crossprod(x[, 1, drop = FALSE], wi) %*% (y - x %*% beta))
115
}

Read our documentation on viewing source code .

Loading