R/fista_LpS.R

In kevinbai92/LSvarEstimate: Low Rank Plus Sparse VAR Model Estimation

#' low-rank plus sparse structure matrix estimation function solved by FISTA algorithm
#' @param A design matrix
#' @param b target vector
#' @param alpha a positive constant the measuring the degree of separating sparse from low-rank, default value is 0.25
#' @param lambda tuning parameter for sparse component
#' @param mu tuning parameter for low-rank component
#' @param niter the number of iterations for FISTA algorithm
#' @param backtracking boolean argument, indicating whether use backtracking method or not
#' @param x.true the true coefficient
#' @return A list consisting of:
#' \item{sparse.comp}{the estimated sparse component}
#' \item{lr.comp}{the estimated low rank component}
#' \item{obj.val}{the values of objective function}
#' \item{rel.err}{the relative error, if we know the true transition matrix}
#' @details This function is the main estimation function. It provides the estimators for sparse and low rank components simultaneously and use the constraint space to separate sparse from low rank. The main theoretical results can be found in the papers: Basu Sumanta, Xianqi Li, George Michalidis (2019).
#' @examples
#' A <- matrix(stats::rnorm(10000), nrow = 400, ncol = 25)
#' b <- matrix(stats::rnorm(10000), nrow = 400, ncol = 25)
#' lambda <- 0.1
#' mu <- 0.5
#' niter <- 1000
#' x.true <- diag(25)
#' fista_LpS(A, b, alpha = 0.25, lambda = lambda, mu = mu, niter = niter, TRUE, x.true = x.true)
#' @references
#' Basu, Sumanta, Xianqi Li, and George Michailidis. "Low rank and structured modeling of high-dimensional vector autoregressions." IEEE Transactions on Signal Processing 67.5 (2019): 1207-1222.
#' @export


fista_LpS <- function(A, b, alpha = 0.25, lambda, mu, niter, backtracking = TRUE, x.true){
    tnew = t <- 1
    p <- dim(A)[2]
    x1 <- matrix(0, nrow = p, ncol = p)
    xnew1 = xnew2 <- x1
    y1 = y2 <- x1
    AtA <- t(A) %*% A
    Atb <- t(A) %*% b

    if(backtracking == TRUE){
        L <- norm(A, "F")^2 / 5
        gamma <- 2
    }else{
        L <- norm(A, "F")^2
    }

    obj.val = rel.err <- c()
    for(i in 1:niter){
        if(backtracking == TRUE){
            L.bar <- L
            found <- FALSE
            while(found == FALSE){
                y <- y1 + y2
                prox1 <- prox_sparse_func(y1, y, A, b, 2*L.bar, lambda, AtA, Atb)
                prox2 <- prox_nuclear_func(y2, y, A, b, 2*L.bar, mu, AtA, Atb)

                ### Restricted solution space
                for(j in 1:p){
                    for(k in 1:p){
                        if(abs(prox2[j,k]) > alpha){
                            prox2[j,k] <- alpha * sign(prox2[j,k])
                        }
                    }
                }

                prox <- prox1 + prox2
                if(f_func(prox, A, b) <= Q_func(prox, y, A, b, L.bar, AtA, Atb)){
                    found <- TRUE
                }else{
                    L.bar <- L.bar * gamma
                }
            }
            L <- L.bar
        }
        x1 <- xnew1
        x2 <- xnew2
        xnew1 <- prox1
        xnew2 <- prox2
        t = tnew
        tnew <- (1 + sqrt(1 + 4*t^2))/2
        y1 <- xnew1 + (t - 1) / tnew * (xnew1 - x1)
        y2 <- xnew2 + (t - 1) / tnew * (xnew2 - x2)
        xnew <- xnew1 + xnew2

        obj.val <- c(obj.val, f_func(xnew, A, b) + sparse_pen(xnew1, lambda) + nuclear_pen(xnew2, mu))
        rel.err <- c(rel.err, norm(xnew - x.true, "F") / norm(x.true, "F"))
    }
    return(list(sparse.comp = xnew1, lr.comp = xnew2, obj.val = obj.val, rel.err = rel.err))
}