# nocov start
# rm(list = ls())
# source("R/sim-data.R")
# l2norm <- function(x) sqrt(sum(x^2))
#
# alpha <- 0.5
# w_e <- 1 ; w_j <- 1
#
# R <- Y - b0
#
# # this is not the right formula with weights.. but ignore for now. need to figure out the weights after
# lambda_max <- (1 / (n * (1 - alpha))) * max( (1 / w_e) * (crossprod(e, R) ),
# max((1 / w_j) * sapply(Phi_j_list, function(i) l2norm(crossprod(i, R)))))
#
# R_without_j <- Y - b0
#
#
#
# design_array <- array(1:8, dim = c(2,2,2))
#
# apply(design_array, c(1,2), FUN = function(i) i * 2)
#
# phi <- design_array[,,1]
#
# all.equal(phi %*% matrix(b1),
# phi %*% b1)
#
#
# Wj <- E * phi
#
# all.equal(phi %*% b1 + 2 * bE * Wj %*% b1,
# (phi + 2 * bE * Wj) %*% b1)
#
dlogit <- function(r, delta) {
dl <- -1 / (1 + exp(r))
dl
}
dls <- function(r, delta) {
dl <- -r
dl
}
# m1 <- gglasso(x = bardet$x, y = bardet$y, group = group1, loss = "ls")
# violations <- gglasso:::KKT(b0 = m1$b0, beta = m1$beta, y = bardet$y, x = bardet$x, lambda = m1$lambda,
# pf = rep(sqrt(5),20), group = group1, thr = 1e-3, loss = "ls")
# names(fit)
#
# tt <- margin(b0 = fit$a0, betaE = fit$bE, beta = fit$beta, gamma = fit$gamma,
# alpha = fit$alpha,
# y = DT$y, phij = fit$Phi_j, xe_phij = fit$XE_Phi_j,
# e = DT$e, df = fit$df, loss = "ls")
# this gives -R = -(Y - hat(Y))
margin <- function(b0, betaE, beta, gamma, alpha, y, phij, xe_phij, e, df, loss = c("ls", "logit")) {
loss <- match.arg(loss)
nobs <- length(as.vector(y))
beta <- as.matrix(beta)
alpha <- as.matrix(alpha)
gamma <- as.matrix(gamma)
b0MAT <- matrix(rep(b0, nobs), nrow = nobs, byrow = TRUE)
# dim(beta)
# dim(b0MAT)
# dim(xe_phij)
# dim(phij)
# dim(alpha)
# dim(beta)
# b0[1:5]
# b0MAT[1:5,1:5]
# browser()
link <- b0MAT + phij %*% beta + matrix(e) %*% matrix(betaE, nrow = 1) + xe_phij %*% alpha
if (loss %in% c("logit")) {
r <- y * link
} else {
r <- as.vector(y) - link
}
fun <- paste("d", loss, sep = "")
# this is a matrix of -R = -(Y - hat(Y)) of dimension nobs x nlambda
dMat <- apply(r, c(1, 2), eval(fun))
# e + dim(xe_phij %*% beta)
# dim(gamma)
#
#
# if (loss %in% c("logit", "sqsvm", "hsvm")) {
# yxdMat <- t(x) %*% (dMat * y)/nobs
# } else yxdMat <- t(x) %*% dMat/nobs
# yxdMat
return(dMat)
}
KKT <- function(b0, betaE, beta, gamma, alpha, y, phij, xe_phij, e, df,
lambda, lambda2, group, we, wj, wje, thr, loss = c("ls", "logit")) {
loss <- match.arg(loss)
bn <- as.integer(max(group))
nobs <- length(as.vector(y))
# this gives -R = -(Y - hat(Y)) for loss = "ls"
dl <- margin(
b0 = b0, betaE = betaE, beta = beta, gamma = gamma, alpha = alpha, y = y,
phij = phij, xe_phij = xe_phij, e = e, df = df, loss = loss
)
# KKT for beta0 -----------------------------------------------------------
# this is the gradient for beta0, this should be of length nlambda
B0 <- t(dl) %*% matrix(1, nrow = nobs) / nobs
ctr <- 0
for (l in 1:length(lambda)) {
if (abs(B0[l, ]) > thr) {
warning("violate at b0 ", B0[l, ], " lambda=", lambda[l], "\n")
ctr <- ctr + 1
}
}
warning("% of violations for beta0", ctr / length(lambda), "\n")
# KKT for betaE -----------------------------------------------------------
# results for betaE
# BE <- matrix(NA, ncol = length(lambda))
ctr <- 0
for (l in 1:length(lambda)) {
xdMat_betaE <- e + rowSums(do.call(cbind, lapply(seq_along(unique(group)), function(j) {
index <- (group == j)
as.matrix(gamma[j, l] * (xe_phij[, index, drop = F] %*% beta[index, l, drop = F]))
})))
dl_norm_betaE <- t(xdMat_betaE) %*% dl[, l, drop = FALSE] / nobs
if (betaE[l] == 0) {
BE <- dl_norm_betaE / (-lambda[l] * (1 - lambda2) * we)
if (abs(BE) > 1 + thr) {
warning("violate at bE = 0", abs(BE), " lambda=", lambda[l], "\n")
ctr <- ctr + 1
}
} else {
BE <- as.vector(dl_norm_betaE + lambda[l] * (1 - lambda2) * we * sign(betaE[l]))
if (abs(BE) > thr) {
warning("violate at bE != 0", abs(BE), " lambda=", lambda[l], "\n")
ctr <- ctr + 1
} # else {warning("no violation at bE != 0", BE, " lambda=",lambda[l], "\n")}
}
}
warning("% of violations for betaE", ctr / length(lambda), "\n")
# return(list(v = ctr/length(lambda)))
# browser()
# KKT for gamma -----------------------------------------------------------
ctr <- 0
for (l in 1:length(lambda)) {
for (g in 1:bn) {
# browser()
ind <- (group == g)
xdMat_gammaj <- as.matrix(betaE[l] * (xe_phij[, ind, drop = F] %*% beta[ind, l, drop = F]))
dl_norm_gammaj <- t(xdMat_gammaj) %*% dl[, l, drop = FALSE] / nobs
if (gamma[g, l] == 0) {
BE <- dl_norm_gammaj / (-lambda[l] * (1 - lambda2) * we)
if (abs(BE) > 1 + thr) {
warning("violate at gamma_j = 0", BE, " lambda=", lambda[l], "\n")
ctr <- ctr + 1
}
} else {
BE <- as.vector(dl_norm_gammaj + lambda[l] * (1 - lambda2) * we * sign(gamma[g, l]))
if (abs(BE) > thr) {
warning("violate at gamma_j != 0", BE, " lambda=", lambda[l], "\n")
ctr <- ctr + 1
} # else {warning("no violation at bE != 0", BE, " lambda=",lambda[l], "\n")}
}
}
}
warning("% of violations for gamma", ctr / length(lambda), "\n")
# KKT for theta -----------------------------------------------------------
# results for beta (aka theta)
# B <- matrix(NA, ncol = length(lambda))
ctr <- 0
for (l in 1:length(lambda)) {
for (g in 1:bn) {
# browser()
ind <- (group == g)
# this is the pre-multiplier of R in eq (17)
xdMat <- phij[, ind, drop = F] + gamma[g, l] * betaE[l] * xe_phij[, ind, drop = F]
# this is the first part of eq (17)
dl_prod <- t(xdMat) %*% dl[, l, drop = FALSE] / nobs
dl_norm <- l2norm(dl_prod)
# l2 norm of beta for the expansion of covariate corresponding to group g
b_norm <- l2norm(beta[ind, l])
if (b_norm != 0) {
AA <- dl_prod + beta[ind, l] * lambda[l] * (1 - lambda2) * wj[g] / b_norm
# warning(AA,"\n")
if (sum(abs(AA)) >= thr) {
warning("violate at bX != 0", sum(abs(AA)), " lambda=", lambda[l], "\n")
ctr <- ctr + 1
}
} else {
BB <- dl_norm - lambda[l] * (1 - lambda2) * wj[g]
# warning(BB,"\n")
if (BB > thr) {
warning("violate at bX = 0", BB, " lambda=", lambda[l], "\n")
ctr <- ctr + 1
}
}
}
}
warning("% of violations for bTheta", ctr / length(lambda), "\n")
return(ctr / length(lambda))
}
# nocov end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.