R/penutils.R In fdaSP: Sparse Functional Data Analysis Methods

Defines functions getGroupMat_OVGLASSOpenfun

```# Authors: Mauro Bernardi
#          Department Statistical Sciences
#      	   Via Cesare Battisti, 241
#          E-mail: mauro.bernardi@unipd.it

# Last change: July 20, 2020

# :::::::::::::::::::::::::::::::::::::::::::::::::
# Functions
# 0. lambdamax (wrap)
# 1. getlambda_LASSO
# 2. getlambda_ELNET
# 3. getlambda_OVGLASSO
# 4. penfun
# 5. getGroupMat_OVGLASSO

# 4. penfun
penfun <- function(method, p, q, splOrd, regpars) {

# get the number of groups
if (q == 1) {
B <- p - splOrd + 1
} else {
B <- (p - splOrd + 1) * (q - splOrd + 1)
}

# lasso
if (method == "lasso") {
Cmat         <- diag(1L, p * q)
Pmat         <- Cmat
Cvec         <- rep(1L, p * q)
grMat        <- NULL
grMatWeights <- diag(1L, p * q)
norme        <- NULL
}

# ridge
if (method == "ridge") {
Cmat         <- diag(1L, p * q)
Pmat         <- Cmat
Cvec         <- rep(1L, p * q)
grMat        <- NULL
grMatWeights <- diag(1L, p * q)
norme        <- NULL
}

# ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# overlap group lasso (function-to-function)
# no overall group
if ((method == "ogl") & (q > 1)) {

# get the group matrices
grMat        <- getGroupMat_OVGLASSO(p = p, q = q, splOrd = splOrd)
grMatWeights <- 1.0 / apply(grMat, 2, sum)

if (is.null(regpars)) {
Cpen <- NULL
Cmat <- NULL
Cvec <- NULL
norme <- NULL
} else {

# get the weights
norme_     <- grMat %*% (regpars^2)

# Set the norm to a given threshold
norme <- fSet2GivenTolerance(x = norme_, n = B, toler = .Machine\$double.eps)

Cvec       <- as.vector(t(grMat) %*% (1.0 / (2.0 * sqrt(norme))))
Cpen       <- diag(1, p * q)
diag(Cpen) <- Cvec
Cmat       <- invvec(Cvec, p, q)
}
}

# ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# overlap group lasso (function-to-function)
# overall group
if ((method == "ogl&1") & (q > 1)) {
# get the group matrices
grMat        <- getGroupMat_OVGLASSO(p = p, q = q, splOrd = splOrd)
grMat        <- rbind(grMat, rep(1L, p * q))
grMatWeights <- 1.0 / apply(grMat, 2, sum)

# get the weights
if (is.null(regpars)) {
Cpen  <- NULL
Cmat  <- NULL
Cvec  <- NULL
norme <- NULL
} else {
# get the weights
# Cvec   <- 1.0 / (2.0 * sqrt(crossprod(regpars)))
norme_ <- grMat %*% (regpars^2)

# Set the norm to a given threshold
norme <- fSet2GivenTolerance(x = norme_, n = B, toler = .Machine\$double.eps)

#Cvec       <- rep(Cvec, p*q) + as.vector(t(grMat) %*% (1.0 / (2.0 * sqrt(norme))))
Cvec       <- as.vector(t(grMat) %*% (1.0 / (2.0 * sqrt(norme))))
Cpen       <- diag(1, p * q)
diag(Cpen) <- Cvec
Cmat       <- invvec(Cvec, p, q)
}
}

# ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# overlap group lasso (function-to-scalar)
# no overall group
if ((method == "ogl") & (q == 1)) {

# get the group matrices
grMat        <- getGroupMat_OVGLASSO(p = p, q = 1, splOrd = splOrd)
grMatWeights <- 1.0 / apply(grMat, 2, sum)

if (is.null(regpars)) {
Cpen <- NULL
Cmat <- NULL
Cvec <- NULL
norme <- NULL
} else {

# get the weights
norme_    <- grMat %*% (regpars^2)

# Set the norm to a given threshold
norme <- fSet2GivenTolerance(x = norme_, n = B, toler = .Machine\$double.eps)

Cvec       <- as.vector(t(grMat) %*% (1.0 / (2.0 * sqrt(norme))))
Cpen       <- diag(1, p)
diag(Cpen) <- Cvec
Cmat       <- Cpen
}
}

# ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# overlap group lasso (function-to-scalar)
# no overall group
if ((method == "ogl&1") & (q == 1)) {

# get the group matrices
grMat        <- getGroupMat_OVGLASSO(p = p, q = 1, splOrd = splOrd)
grMat        <- rbind(grMat, rep(1L, p))
grMatWeights <- 1.0 / apply(grMat, 2, sum)

if (is.null(regpars)) {
Cpen  <- NULL
Cmat  <- NULL
Cvec  <- NULL
norme <- NULL
} else {

# get the weights
norme_    <- grMat %*% (regpars^2)

# Set the norm to a given threshold
norme <- fSet2GivenTolerance(x = norme_, n = B, toler = .Machine\$double.eps)

Cvec       <- as.vector(t(grMat) %*% (1.0 / (2.0 * sqrt(norme))))
Cpen       <- diag(1, p)
diag(Cpen) <- Cvec
Cmat       <- Cpen
}
}

# get the output
res              <- NULL
res\$Cpen         <- Cpen
res\$Cmat         <- Cmat
res\$Cvec         <- Cvec
res\$grMat        <- grMat
res\$grMatWeights <- grMatWeights
res\$regp.grnorm  <- norme

# Return output
return(res)
}

# 5. getGroupMat_OVGLASSO
# get groups of overlap group LASSO
getGroupMat_OVGLASSO <- function(p, q, splOrd = c(3, 4)) {

if (q > 1) {
# get the number of groups
B     <- (p - splOrd + 1) * (q - splOrd + 1)
grMat <- matrix(0, B, p * q)

# get the groups matrix
if ((splOrd == 3) || (splOrd == 4)) {
# spline of order 3
if (splOrd == 3) {
cnt <- 1
for (i in 0:(p-splOrd)) {
for (j in 1:(q-splOrd+1)){
ind             <- c(q*i+j, q*i+j+1, q*i+j+2, p+q*i+j, p+q*i+j+1, p+q*i+j+2, 2*p+q*i+j, 2*p+q*i+j+1, 2*p+q*i+j+2)
grMat[cnt, ind] <- 1
cnt             <- cnt + 1
}
}
}

# spline of order 4
if (splOrd == 4) {
# cnt <- 1
# for (i in 0:(p-splOrd)) {
#   for (j in 1:(q-splOrd+1)) {
#     ind             <- c(q*i+j, q*i+j+1, q*i+j+2, q*i+j+3, p+q*i+j, p+q*i+j+1, p+q*i+j+2, p+q*i+j+3, 2*p+q*i+j, 2*p+q*i+j+1, 2*p+q*i+j+2, 2*p+q*i+j+3, 3*p+q*i+j, 3*p+q*i+j+1, 3*p+q*i+j+2, 3*p+q*i+j+3)
#     grMat[cnt, ind] <- 1
#     cnt             <- cnt + 1
#   }
# }
# cnt <- 1
# for (i in 0:(p - splOrd)) {
#   init_row <- i * q
#   for (j in 1:(q - splOrd + 1)) {
#     init_col <- j
#     v1       <- init_row + init_col + 0:(splOrd-1)
#     v2       <- q + init_row + init_col + 0:(splOrd-1)
#     v3       <- 2*q + init_row + init_col + 0:(splOrd-1)
#     v4       <- 3*q + init_row + init_col + 0:(splOrd-1)
#     ind      <- c(v1, v2, v3, v4)
#
#     grMat[cnt, ind] <- 1
#     cnt             <- cnt + 1
#   }
# }
B <- (p - splOrd + 1) * (q - splOrd + 1)
grMat <- matrix(0, B, p * q)
cnt <- 1
for (j in 1:(q - splOrd + 1)) {
init_col <- j
for (i in 0:(p - splOrd)) {
init_row <- i*q

v1  <- init_row + init_col + 0:3
v2  <- q + init_row + init_col + 0:3
v3  <- 2*q + init_row + init_col + 0:3
v4  <- 3*q + init_row + init_col + 0:3
ind <- c(v1, v2, v3, v4)

grMat[cnt, ind] <- 1
cnt <- cnt + 1
}
}
}
} else {
cat("Error in getGroupMat_OVGLASSO, the spline order in not equat to 3 or 4!\n")
}
} else if (q == 1) {
# get the number of groups
B     <- (p - splOrd + 1)
grMat <- matrix(0, B, p)

# get the groups matrix
if ((splOrd == 3) || (splOrd == 4)) {
# spline of order 3
if (splOrd == 3) {
cnt <- 1
for (i in 0:(p-splOrd)) {
ind             <- c(i+1, i+2, i+3)
grMat[cnt, ind] <- 1
cnt             <- cnt + 1
}
}

# spline of order 4
if (splOrd == 4) {
cnt <- 1
for (i in 0:(p-splOrd)) {
ind             <- c(i+1, i+2, i+3, i+4)
grMat[cnt, ind] <- 1
cnt             <- cnt + 1
}
}
} else {
cat("Error in getGroupMat_OVGLASSO, the spline order in not equat to 3 or 4!\n")
}
} else {
cat("q is negative or equal to zero in getGroupMat_OVGLASSO!\n")
}

# return output
return(grMat)
}
```

Try the fdaSP package in your browser

Any scripts or data that you put into this service are public.

fdaSP documentation built on Oct. 6, 2023, 1:07 a.m.