PGA <- function(x,y, max_iter){    # X should be centered. i.e., mean 0
	n <- length(y)
	x <- as.matrix(x)	
	x <- cbind(rep(1,n),x)
	xdim <- dim(x)
	if (xdim[1] != n) 
		stop("check dim of x")
	p=xdim[2]
	beta_current=rep(0,p)
	beta_all<-matrix(rep(0,p*max_iter),p,max_iter)	
	j<-rep(0,max_iter)

	e=y
	e_new=e
	e_new2=rep(0,max_iter)
	e2_new=rep(0,max_iter)
	e2_prev=0
	
	for (iter in (1:max_iter)){
		corrs=(t(e_new)%*%x)^2/colSums(x^2)
		j[iter]=which.max(corrs)
		beta_sub=t(e_new)%*%x[,j[iter]]/sum(x[,j[iter]]^2)
		beta_current[j[iter]]=beta_current[j[iter]] + beta_sub    ## add lambda?
		beta_all[,iter] = beta_current
		e <- e-beta_sub*x[,j[iter]]                               ##  also needs modify?
	
		e_new=e
		e2=e^2
	
		e_new2[iter]=mean(e_new^2)
		e2_new[iter]=mean(e2)
		if (abs(e2_new[iter]-e2_prev)<1e-2)   ## stopping criteria
			break
		e2_prev=e2_new[iter]
	}

	return(list(coeff_final=beta_current,coeff_all=beta_all))
}