sparsePcaArd <- function(
				x,	# Data matrix (rows: variables; column; samples)
				k,	# Factor dimension	
				a,	# Shape parameter for the inverse Gamma prior for the scale variable [gamma] (Archambreau et al. NIPS, 2008)
				b,	# Scale parameter for the inverse Gamma prior for the scale variable [gamma] (Archambreau et al. NIPS, 2008)
				itrmax = 1000	# Max number of iteration
					)
 {

	p <- dim( x )[1]	# Dimension of data
	n <- dim( x )[2]	# Dimension of data
 	unit <- diag( rep( 1, length = k ) )	# Identity matrix

 # Parameters involved in posterior on presion inducing sparsity of the loading matrix
 	chi <- 2 * b				# Inverse Gaussian parameter 1 (Eq(18) of Archambreau et al. NIPS, 2008) 
 	omega <- - a / ( p * k ) + 1 / 2	# Inverse Gaussian parameter 2 (Eq(18) of Archambreau et al. NIPS, 2008)

 # Set initial paramaeters
 	cmat <- x %*% t(x)					# Squared-sum matrix of data 
 	eig <- eigen( cmat, symmetric = TRUE )		# Eigenvalue decomposition of data
 	load <- eig$vectors[ , 1:k] 				# Factor loading matrix (k dominat eigenvectors)
 	psi <- 0.01* sum( diag( cmat ) ) / ( n * p ) 	# Noise variances (scalar)
 
 # Loop for variational mean field approximation 	
	for( itr in 1:itrmax ){

		print("--- Iteration ---")
		print( itr )

		# Parameters of the posterior distribution on k latent factor variables 
		# Eq(17) of Archambreau et al. NIPS, 2008
 		barS <- solve( t( load ) %*% load / psi + unit )	# Posterior covariance matrix 
 		barZ <- barS %*% t( load ) %*% x / psi 			# Posterior mean

		# Parameters involved in the posterior distribution on the scale parameters
		# Eq(17) of Archambreau et al. NIPS, 2008
 		phi <- load^2

		# Eq(23) of Archambreau et al, NIPS, 2008
 		denom <- apply( sqrt( phi*chi ), MARGIN=c(1,2), besselK, nu=omega )
 		numer <- apply( sqrt( phi*chi ), MARGIN=c(1,2), besselK, nu=(omega+1) )
 		gammaExp <- sqrt(chi/phi)*numer/denom	# Expectation of gamma
		gammaExp[ is.nan(gammaExp)==TRUE ] <- 0
		gammaExp[ gammaExp >= 10^10 ] <- 10^10

		# Sparsity configuration matrix 
		sparseInd <- matrix( 1, ncol=k, nrow=p )	# Initialization
		sparseInd[ gammaExp >= 10^10 ] <- 0		# If the expected gamma is large enough, the correponding element of loading matrix become zero.  

		# Upper-thresholding for the expected gamma
		gammaExp[gammaExp >= 1e15] <- 1e15

 		# Some sufficient statistics 
 		Txz <- x %*% t( barZ ) 
 		Tzz <- barZ %*% t( barZ )
 		VarZ <- n*barS + Tzz
 
		# Eq(20) of Archambreau et al. NIPS, 2008
		#print("--- Expected Gamma ---")
		#print( gammaExp )
 		for(i in 1:p){
			#print(i)	
			prj <- solve ( diag( gammaExp[i,] ) + VarZ / psi  )
			load[i,] <- prj %*% Txz[i,] / psi
 		}

		# Eq(21) of Archambreau et al. NIPS, 2008
 		tmp <- t( load )%*%load
 		psi <- ( sum( diag( x%*%t(x) - 2 * Txz %*% t(load) ) ) + sum( diag(VarZ%*%tmp) ) ) / ( n * p )

		#print( "--- Loadings ---" )
 		#print( round( load, digits=5 ) )
		#print( "--- Noise Variance ---")
 		#print( round( psi, digits=5 ) )
		degs.ARD <- which( abs( load ) >= 10^{-7}, arr.ind = FALSE )
		degs.ARD <- length( degs.ARD ) / ( p * k )
		print( "--- Degree of Sparseness ---" )
 		print( round( degs.ARD, digits=5 ) )


 	}

 	return( list( 
			load = load, 
			psi = psi, 
			gammaExpect = gammaExp,
			sparseInd = sparseInd 
			) 
		)

 }