#' Factor Loading Adaptive Shrinkage (K factors version)
#'
#' flash provide rank K matrix decomposition with greedy algorithm
#'
#' @param Y is the data matrix (N by P)
#' @param K is the max number of factor user want. the output will provide the actual number of factor
#' automaticaly.
#'
#' @details greedy algorithm on the residual matrix to get a rank one matrix decomposition
#'
#' @export greedy
#'
#' @importFrom ashr ash
#'
#' @return list of factor, loading and variance of noise matrix
#' \itemize{
#' \item{\code{l}} {is a K by N matrix for loadings}
#' \item{\code{f}} {is a K by P matrix for factors}
#' }
#' @examples
#' N = 100
#' P = 200
#' Y = matrix(rnorm(N*P,0,1),ncol=P)
#' g = greedy(Y,10)
#'
greedy = function(Y,K){
N = dim(Y)[1]
P = dim(Y)[2]
#initial the residual for the greedy algorithm
residual = Y
# do rank one decomposition
r_flash = flash(residual)
l_temp = r_flash$l
f_temp = r_flash$f
# test whether it is zero
if(sum(l_temp^2)==0 | sum(f_temp^2)==0){
L_out = 0
F_out = 0
return(list(l = L_out,f = F_out))
}else{
L_out = l_temp
F_out = f_temp
#get the new residual
residual = residual - l_temp %*% t(f_temp)
#itreation for the rank K-1
for(k in 1:K){
#rank one for residual
r_flash = flash(residual)
l_temp = r_flash$l
f_temp = r_flash$f
# get the new residual
residual = residual - l_temp %*% t(f_temp)
#check if we should stop at this iteration
if(sum(l_temp^2)==0 | sum(f_temp^2)==0){
break
}else{
# if not go to next step and restore the l and f
L_out = cbind(L_out,l_temp)
F_out = cbind(F_out,f_temp)
}
# for loop can control the number of the factors if needed
}
return(list(l = L_out,f = F_out))
}
# no return here, since return before
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.