4 FSSS

In this section we implement the main algorithm, Feature Subspace Stability Selection (FSSS), which is described in Algorithm 1 of the paper.

Input: Features \(X\in\mathbb{R}^{n\times p}\), response variable \(y\in\mathbb{R}^n\), stability threshold \(\alpha\in(1/2, 1)\), even integer \(B\), and a base variable selection procedure
Subsampling: Obtain estimates \(\{\widehat{S}^{(2\ell-1)}, \widehat{S}^{(2\ell)}\}_{l=1}^{B/2}\) from applying base procedure to complementary subsamples, and compute \(\mathcal{P}_{\rm avg} = \frac{1}{B} \sum_{l=1}^B \mathcal{P}_{\widehat S^{(\ell)}}\).
Sequentially add variables: Set \(\widehat{S} \gets \varnothing\) and perform:
- 1. Let \(\mathcal{F} \gets \{j\in\{1, \dots, p\}\setminus \widehat{S} : \text{tr}( \mathcal{P}_{v} \mathcal{P}_{\rm avg}) \geq \alpha \text{ for }v = \mathcal{P}_{\widehat{S}^\perp}X_j\}\).
- 1. If \(|\mathcal{F}| > 0\), randomly select one \(j\) from \(\mathcal{F}\), and:
  - 1. If \(\pi(\hat{S}\cup \{j\}) \geq \alpha\), go to step (c);
  - 1. Otherwise, let \(\mathcal{F} \gets \mathcal{F} \setminus \{j\}\), and repeat (b) until either (b)(i) is satisfied, or \(\mathcal{F}\) becomes empty.
- 1. If \(|\mathcal{F}| = 0\), return \(\widehat{S}\). Otherwise, let \(\widehat{S} \gets \widehat{S} \cup \{j\}\), and repeat (a)–(c).
Output: The selection set \(\widehat{S}\).

The subsampling part of this algorithm was implemented in the previous section. What remains is the part where we sequentially add variables. There are two helper functions needed, which we define next:

#' FSSS algorithm
#' 
#' Algorithm 1. A function that implements the FSSS algorithm
#' 
#' @param X the n by p design matrix
#' @param bags a list returned from l0_subsampling
#' @param alpha a value in (0,1), the threshold for stability with default value 0.9
#' @param type a character, either "greedy" (greedy algorithm) or "any-path" 
#' @param verbose 0 or 1, logical value for outputting the selection path
#' 
#' @return a vector representing the selection set
#' 
#' @export
FSSS <- function(X, bags, alpha = 0.9, type = "greedy", verbose = 1) {
  Pavg = bags$Pavg
  base_lst = bags$base_lst
  d = ncol(X)
  
  if(nrow(Pavg) != nrow(X) ) stop("The Pavg has different dimension with X!")
  if(!(type %in% c("greedy", "any-path"))) stop("The type of FSSS must be either 'greedy' or 'any-path'!")
  
  node_par = c()
  orth_par = NULL
  stability_output = 0
  
  while(TRUE) {
    if(length(node_par) == d) {
      break
    }
    node_chil = lapply(setdiff(1:d, node_par), function(i) {
      c(node_par, i)
    })
    if (is.null(orth_par)) {
      dir_chil = lapply(node_chil, function(x) {
        X[,x]
      })
    } else {
      dir_chil = lapply(node_chil, function(x) { # find the gs-orthogonal vector
        Xj = X[, setdiff(unlist(x), node_par)]
        terms = lapply(seq_len(ncol(orth_par)), function(i) {
          sum(Xj * orth_par[,i]) / sum(orth_par[,i]^2) *  matrix(orth_par[,i])
        })
        rs = Xj - apply( matrix( unlist(terms), ncol = ncol(orth_par) ), 1, sum)
        return(rs)
      })
    }
    psi_chil = sapply(dir_chil, function(x) {
      if(norm(unlist(x), "2") < 1e-7) {
        return(1)
      } else {
        return( 1 - t(unlist(x)) %*% Pavg %*% unlist(x) / norm(unlist(x), "2")^2 )
      }
    })
    
    if(type == "greedy") {
      # for the greedy path
      idx = which.min(psi_chil)
      stability = find_stab(base_lst, X[,node_chil[[idx]], drop = FALSE])
    } else {
      # for any path
      candidate = which(psi_chil < 1 - alpha)
      while(length(candidate) > 0) {
        idx = ifelse(length(candidate) != 1, sample(candidate, 1), candidate)
        stability = find_stab(base_lst, X[,node_chil[[idx]], drop = FALSE])
        if( stability >= alpha & psi_chil[idx] <= 1 - alpha ) {
          break
        } else {
          candidate = setdiff(candidate, idx)
        }
      }
      
      if(length(candidate) == 0) break
    }
    
    # add feature to selection set
    if( psi_chil[idx] <= 1 - alpha & stability >= alpha) {
      orth_par = cbind( orth_par, dir_chil[[idx]])
      node_par = node_chil[[idx]]
      if(verbose == 1) cat("current selection set:", node_par, "\n")
    } else {
      break
    }
  }
  
  return(node_par)
}

We can try running this 10 times on our example:

Selection_set = list()
while(length(Selection_set) < 10) {
  Selection_set_obj = FSSS(ex$X, ex$bags, alpha = 0.85, type = "any-path", verbose = 0)
  S = sort(Selection_set_obj)
  if(!(list(S) %in% Selection_set) ) {
    Selection_set = c(Selection_set, list(S))
  }
  cat("finished", length(Selection_set), "\n")
}

## finished 1 
## finished 2 
## finished 3 
## finished 4 
## finished 5 
## finished 6 
## finished 7 
## finished 7 
## finished 8 
## finished 8 
## finished 8 
## finished 8 
## finished 9 
## finished 9 
## finished 9 
## finished 10