PCA by ALS

Description

This program approximates the n x nvar input matrix Y
by the product of n x ndim F and nvar x ndim A matrices as
Y = F %*% t(A) + E

When the optimal scaling of Y is employed, the result should be close to the one from PRINCIPALS of Young, Takane and de Leeuw (1978).

Usage

pca_by_als(
  Y,
  ndim = 2,
  mlevel = 4,
  maxiter = 1000,
  epsd = 1e-06,
  SQUAREM = 0,
  nSQUAREM = 1,
  maxalpha = -1,
  always = 1,
  reset1 = 1,
  reset2 = 2,
  print = 1,
  plot = 0
)

Arguments

Details

The puropose of this function is to demonstrate how the ALS (Alternating Least Squares) and inline SQUAREM can be utilized in PCA context.

The parameter set (F,A) is divided into F and A and the conditional minimization of rss w.r.t. F and w.r.t. A will be alternated until convergence.

When requested, the optimal scaling of each column of Y will be performed as a part of ALS iteration.

The program may not be very efficient but it reduces rss at each phase of iteration.

This function finds F and A matrices of rank ndim and the optimal transformation of each column of the data matrix Y which minimize sum over j of( (Yhat[,j]-(F%*%t(A))[,j])^2 / var(Yhat[,j]) )
where Yhat=t(Y) is the column-wise optimally scaled Y.

The actual minimization is done by normalizing Yhat so that var(Yhat[,j])=1 and minimizing
rss=sum( (Yhat-F%*%A)^2 )
by alternating least squares..

The core part of this function is as follows:

# Repeat the following 4 phases until convergence.

# 1 Given Yhat and F, estimate A
A <- t( solve(t(F)%*%F)%*%t(F)%*%Yhat )
# 2 Given Yhat and F, estimate F
F <- t( solve(t(A)%*%A)%*%t(A)%*%t(Yhat) )
# 3 Given F and A, estimate each column of Yhat
Yhat[,j] <- optimally scaled Y[,j] which minimizes rss.
# 4 Column standardize Yhat.
# Check convergence: rss should decrease during the iteration.
#

Note that the column standardized Y will be used as the initial Yhat.

Note that when mlevel=0 rss will be minimized without normalization of Yhat. Therefore, the result should be equievalent to PCA after suitable normalization of F and A.

Value

A list of

  F   n x ndim     componet score matrix whose dispersion matrix is I
  A   nvar x ndim  componet loadings matrix satisfying t(A)%*%A = diagonal
  Yhat n x nvar    optimally scaled data matrix
  rmse             = sqrt( rss/(n*nvar) )
  mlevel           transformation level

References

Young, W. F., Takane, Y. and de Leeuw, J. (1978)
The principal components of mixed measurement level multivariate data:
An alternating least squares method with optimal scaling features.
Psychometrika volume 43, pages 279-281

Examples

#
# Comparison of pca by als and native prcomp
#
#
# generate column centered Y
seed <- 1701
set.seed(seed)
n <- 50; nvar <- 6
Y <- matrix(rnorm(n*nvar),n,nvar)
Y <- scale( Y, center=TRUE )
attributes(Y)$`scale.center`=NULL

# pca by als (no optimal transformation)
ndim <- 2
res <- pca_by_als( Y, ndim, SQUAREM=0, always=1, print=1, mlevel=0 )

# normalize the result so that Disp(F)=I and t(A)%*%A=diagonal
FA <- normalize_fa( res$F, res$A )
F1 <- FA$F; A1 <- FA$A
rmse1 <- sqrt(ssq(Y-F1%*%t(A1))/n/nvar)

# pca by native prcomp: Y == F2 %*% t(A2), t(F2)%*%F2 = diag, t(A2)%*%A2 = I
resprc <- prcomp( Y, rank=ndim )
sqrteva <- resprc$sdev[1:ndim]*sqrt((n-1))   # vardef adjusted, eva of Y'Y
A2 <- resprc$rotation
# get F2 by solving Y=F2%*%t(A2) or t(Y)=A2%*%t(F2)  for F2
F2 <- Y%*%A2%*%solve(t(A2)%*%A2)

# from prcomp to pca:
# Note that sqrt eva from native prcomp is the eva of Y'Y, not Y'Y/n
A21 <- A2%*%diag(sqrteva)/sqrt(n)
F21 <- F2%*%diag(sqrt(n)/sqrteva)
rmse21 <- sqrt(ssq(Y-F21%*%t(A21))/n/nvar)

# FA <- normalize_fa( F21, A21 )
# F21n <- FA$F; A21n <- FA$A

Print(A1,A21)
Print(rmse1,rmse21)



# Simulating PRINCIPALS

# generate F according to the FA model
seed <- 1701
set.seed(seed)
n <- 200; nvar <- 9
ndim0 <- 3
F <- matrix(rnorm(n*ndim0),n,ndim0)
A <- matrix(rnorm(nvar*ndim0), nvar,ndim0)
A[,3] <- 0.5*A[,3]

# transform Y
Y <- F%*%t(A) + 0.2*matrix(rnorm(n*nvar),n,nvar)
 # transform Y
 Y[,2]=(Y[,2]-min(Y[,2]))^3
 Y[,3]=(Y[,3]-min(Y[,3]))^3
 Y[,4]=(Y[,4]-min(Y[,4]))^(1/3)
 Y[,5]=(Y[,5]-min(Y[,5]))^(1/3)
Y <- scale(Y, center=TRUE, scale=TRUE)

# round Y
Y <- round(10*Y,0)

# PRINCIPALS
ndim <- 2
res5 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=5)
res4 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=4)
res3 <- pca_by_als(Y, ndim, SQUAREM=3, always=1, print=1, plot=2, mlevel=3)
Print(res5$mlevel, res5$rmse)
Print(res4$mlevel, res4$rmse)
Print(res3$mlevel, res3$rmse)