R/als.mg.r

In gesca: Generalized Structured Component Analysis (GSCA)

als.mg <- function (Z, W0, A0, W, A, V, I, PHT, nvar, nlv, ng, itmax, ceps)
{
	#---------------------------------------------------
	# ALS algorithm for GSCA (multiple group analysis)
	# Heungsun Hwang, Sunmee Kim
	# No calcuation of Kronecker products  
	# Last revised Jan 29, 2016
	# Elements of A and W can be fixed to constants (0 > & < 1)
	# Z = block-diagonal matrix of "normalized" data for all groups
	#---------------------------------------------------

	sizea <- ncol(A)
	aindex0 <- which(A0 >= 1)
	Psi <- Z%*%V%*%I
	Gamma <- Z%*%W
	it <- 0			# iteration counter
	imp <- 100000	# initial improvement
	f0 <- 10^10 	# initial function value

	while (it <= itmax && imp > ceps) {
		it <- it + 1
		# step 1: update A
		for (t in 1:sizea) {
			if ( !all(A0[,t] == 0) ) {
				H1 <- diag(1,sizea)
				H1[t,t] <- 0
				aindex <- which(A0[,t] >= 1) 	# free parameters
				if (length(aindex) != 0) {
					a <- A[,t,drop=FALSE]
					a[aindex] <- 0					# if fixed values, not all zeros
					e <- matrix(0,1,sizea)
					e[t] <- 1
					Y <- Psi - Gamma%*%(A%*%H1 + a%*%e)	# Revised 2017-09-28
					X <- Gamma[,aindex,drop=FALSE]
					A[aindex,t] <- solve(t(X)%*%X, t(X)%*%Y%*%t(e))
				}
			}
		}		
		vecA <- A[aindex0]
		A[aindex0] <- PHT%*%vecA
		
		# step 2: update W
		kk = 0
		for (g in 1:ng) {
			k <- kk + 1
			kk <- kk + nlv
			p <- g*nvar + (g-1)*nlv
			s <- 0
			for (j in k:kk) {
				s <- s + 1
				windex <- which(W0[,j] == 99)
				w <- W[,j,drop=FALSE]
				w[windex] <- 0 # if fixed values, not all zeros
				beta <- I[p+s,,drop=FALSE] - A[j,,drop=FALSE]
				H2 <- diag(1,ng*nlv)
				H2[j,j] <- 0
				H3 <- diag(1,sizea*ng)
				H3[p+s, p+s] <- 0
				Delta <- W%*%H2%*%A - V%*%H3%*%I - w%*%beta
				Zp <- Z[,windex]
				temptheta <- solve(t(Zp)%*%Zp,t(Zp))%*%(Z%*%Delta)%*%t(beta)
				theta <- t(solve(t(beta%*%t(beta)),t(temptheta)))
				zw <- Zp%*%theta
				theta <- theta%*%(1/sqrt(t(zw)%*%zw))
				W[windex,j] <- theta
				V[windex,p+s] <- theta
			}
		}
		
		Gamma <- Z%*%W
		Psi <- Z%*%V%*%I
		dif <- Psi-Gamma%*%A
		f <- sum(diag(t(dif)%*%dif))
		imp <- f0-f
		info <- c(it,f,f0,imp)
		f0 <- f	
	}
	
	output.als.mg <- list(W = W, A = A, Psi = Psi, Gamma = Gamma, f = f, it = it, imp = imp)
	output.als.mg
}