Examples of GN and NR

Description

Examples of GN and NR in the context of Logistic Regression. Multidimensional Scaling, Factor Analysis and PCA.

Examples

## Not run: 

#
# For comparison reason, the following package is required:
library(nlsr)
#



#####################################################################
#
#     ML and Bayes MAP Estimation of Logistic Regression Parameters
#            by Newton-Raphson Method
#
#####################################################################
#
# probability (logistic function with matrix X, no 1.7)
# X is the regressor matrix, beta is the regression coefficient vector
# Returns a vector of length nrow(X).
Prob <- function( beta, X ){
 if( is.matrix(X) | is.data.frame(X) ) z=X%*%beta
 else z=X*beta
 P=1/(1+exp(-z))
 return( P )
} # end of P
#
#
#
# This is the objective function to be minimized:
# Minus log likelihood
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a scalar.
#
mllh <- function( beta, data ){
 # minus log likelihood: multiple regressor matrix X
 locx=which(regexpr("x",colnames(data)) > 0)
 X=as.matrix(data[,locx]); f=data$f; r=data$r
 p=Prob(beta,X)
 llh=sum( r*log(p)+(f-r)*log(1-p) )
 return( -llh )
} # end of mllh
#
#
# First derivative of mllh wrt beta.
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a vector of size length(beta).
#
dmllh_a <- function( beta, data ){
 # analytic first derivative of mllh: multiple regressor matrix X
 locx=which(regexpr("x",colnames(data)) > 0)
 X=as.matrix(data[,locx]); f=data$f; r=data$r
 p=Prob(beta,X)
 dab=colSums(c(r-f*p)*X)
 # cat(".da.")
 return( -dab )
} # end of dmllh_a
#
#
dmllh_n <- function( beta, data, mllh=NULL ){
 # numeric first derivative of mllh: multiple regressor matrix X
 # Not Used.
 res=JacobianMat( beta, mllh, data=data)
 return( res )
} # end of dmllh_n
#
#
#
# Second derivative of mllh wrt beta.
# parameter beta is the first argument and the second argument data
# will be passed to this via ... parameter.
# Returns a square matrix of size length(beta).
#
d2mllh_nn <- function( param, data, mllh=NULL ){
 # numeric second derivative matrix using numeric first derivative of mllh
 # dmllh_n comes from the parent environment
 # Not Used.
 res=JacobianMat( param, dmllh_n, data=data, mllh=mllh )
 # cat(".d2an.")
 return( res )
} # end of d2mllh_nn
#
#
d2mllh_an <- function( param, data ){
 # numeric second derivative matrix using analytic first derivative of mllh
 # dmllh_a comes from the parent environment
 # Not Used.
 res=JacobianMat( param, dmllh_a, data=data )
 # cat(".d2an.")
 return( res )
} # end of d2mllh_an
#
#
d2mllh_aa <- function( param, data ){
 # analytic second derivative matrix of mllh
 # not yet available
 locx=which(regexpr("x",colnames(data)) > 0)
 X=as.matrix(data[,locx]); f=data$f; r=data$r
 p=Prob(param,X)
 dPdbetaa=c(p*(1-p))*X
 H=t(f*X)%*%dPdbetaa
 return(H)
} # end of d2mllh_aa
#
#
#
#
#
#
# generate data
#
#
# seed for random numbers
seed=1701
set.seed(seed)
#
# true value of beta
beta=c(-1,2,3)
# # of parameters = # of regressors = # of columns of X including intercept.
nq=length(beta)
#
# # of trials per obs
ff=c(1,2,3)               # will be recycled.
# ff=1
#
# generate regressor matrix X
n=50
#
# regressor matrix as 1 and normal random variable
X=cbind(1, matrix(rnorm(n*(nq-1)),n,nq-1) )
colnames(X)=paste("x",(1:nq)-1,sep="")
#
# true response probability
p=Prob(beta,X)
#
# generate f=# of trials, r=# of successes, y=sample proportion
f=sample( ff, n, rep=1 )
r=mapply( function(size,prob) rbinom(1,size,prob) , f, p )
y=r/f
#
# data frame
data=data.frame( id=1:n, X, f, r, y, p )
print(data)
#
# initial value
# beta0=c(0,1,1)
beta0=rev(beta)
#
#
#
# Estimation Starts Here.
#
#
#
# Arguments to nlminb or NR
#
#  First "beta0" argument is the starting value of the parameter.
#  Second argument "mllh" is the objective function to be minimized
#  which must have "beta0" as its first argument.
#  The additional argument to this function can be passed through ... .
#  or, the objects in the environment in which the objective function
#  "mllh" is defined can be accessed from within the objective function.
#  In the example below, data, which contains f and r,
#  is passed to "mllh"  via ... argument of nlminb or NR.
#
#  Note that, unlike nlm or nlfb which require the gradient to be returned
#  as an attribute "gradient" of the function value,
#  nlminb or NR uses the gradient or hessian functions which return
#   the gradient or hessian as their values.
#
#
# mle by native nlminb
# no analytic gradient and hessian
# Ones like dmllh_n or d2mllh_nn must be used in nlminb.
resnlminb=nlminb( beta0, mllh, data=data, control=list(trace=0) )
# with analytic gradient
resnlminb2=nlminb( beta0, mllh, data=data, control=list(trace=0)
                   , gradient=dmllh_a )
# with analytic gradient and Hessian
resnlminb3=nlminb( beta0, mllh, data=data, control=list(trace=0)
                   , gradient=dmllh_a, hessian=d2mllh_aa )
#
#
# solution by NR
# no analytic gradient and hessian
# Functions similar to dmllh_n or d2mllh_nn will be created by in NR.
resNR=NR( x=beta0, mllh,  data=data, print=2 )
# with analytic gradient
resNR2=NR( x=beta0, mllh,  data=data, gradient=dmllh_a, print=2 )
# with analytic gradient and Hessian
resNR3=NR( x=beta0, mllh,  data=data, gradient=dmllh_a, hessian=d2mllh_aa
           , print=2 )
#
#
# comparison of results
p1=resnlminb$par; p2=resnlminb2$par; p3=resnlminb3$par
pNR1=resNR$par; pNR2=resNR2$par; pNR3=resNR3$par
cat("\n\nComparizon of the Results:  Parameters\n")
Print(p1,p2,p3)
Print(pNR1,pNR2,pNR3)
#
# comparison of gradient
g1=dmllh_a( p1, data ); g2=dmllh_a( p2, data ); g3=dmllh_a( p3, data )
gNR1=dmllh_a( pNR1, data ); gNR2=dmllh_a( pNR2, data )
gNR3=dmllh_a( pNR3, data )
cat("\n\nComparison of the Results:  Gradients\n")
Print(g1,g2,g3)
Print(gNR1,gNR2,gNR3)
#
#
#
#
#
# Bayes MAP
#
#
# prior
#
# The prior distribution of beta is N(mu, inv(invSigma)), where
# mu is the prior mean vector and invSigma is the prior dispersion matrix.
#
# In addition to data which contains (r,f), mu and invSigma will be
# passed to the objective function "mlogpostpdf" via ... argument
# of NR or nlminb.
#
#
logpriorpdf <- function( beta, mu, invSigma ){
 # logarithm of multivariate normal pdf (constant terms omitted)
 lpdf=-0.5*c(beta-mu)%*%invSigma%*%c(beta-mu) + 0.5*det(invSigma)
 return(lpdf)
} # logpriorpdf
#
# popsterior pdf is proportional to the likelihood x prior pdf
mlogpostpdf <- function( beta, data, mu, invSigma ){
 # minus log posterior density
 # log post pdf of beta is  log likelihood  +  log prior pdf
 lppdf= -mllh( beta, data ) +  logpriorpdf( beta, mu, invSigma )
 return( -lppdf )
} # end of mlogpostpdf
#
#
#
# prior constant (hyper parameters) vague prior N(0,100 I)
mu=rep(0,length(beta0)); invSigma=0.01*diag(length(beta0))
#
# MAP estimates
resNR4=NR( x=beta0, mlogpostpdf,  data=data, mu=mu, invSigma=invSigma  )
pNR4=resNR4$par;
resnlminb4=nlminb( beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma
                   , control=list(trace=0) )
p4=resnlminb4$par
#
# evaluate the final gradient numerically
gNR4=JacobianMat( pNR4, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
g4=JacobianMat( p4, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
#
#
# prior constant (hyper parameters) sharp prior N(0,I)
mu=rep(0,length(beta0)); invSigma=diag(length(beta0))
#
# MAP estimates
resNR5=NR( x=beta0, mlogpostpdf,  data=data, mu=mu, invSigma=invSigma  )
pNR5=resNR5$par
resnlminb5=nlminb( beta0, mlogpostpdf, data=data, mu=mu, invSigma=invSigma
                   , control=list(trace=0) )
p5=resnlminb5$par
#
# evaluate the final gradient numerically
gNR5=JacobianMat( pNR5, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
g5=JacobianMat( p5, mlogpostpdf, data=data, mu=mu, invSigma=invSigma )
#
cat("\n\nComparison\n")
Print(p1,p4,p5)
Print(pNR1,pNR4,pNR5)
Print(gNR4,g4,gNR5,g5)
#
#
#####################################################################
#
#    End of
#     ML and Bayes MAP Estimation of Logistic Regression Parameters
#            by Newton-Raphson Method
#
#####################################################################
#
#
#
#
#
#
#
#
#
#####################################################################
#
#                 Multidimensional Scaling
#
#####################################################################
#
#
residmds <- function( vecX ){
 # residual of metric mds:  vech(O-D)
 # Shin-ichi Mayeakawa
 # 20210713,20220713,14,22
 #
 # Note that ns, O, and locLH come from the environment where
 # this function was defined.
 #
 # convert vecX to matrix X
 X=matrix(vecX,ns)
 #
 # Euclidean distance
 XX=rowSums(X*X)
 XX11=XX%*%matrix(1,1,ns)
 D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
 #
 # lower half of O-D
 resid=c( (O-D)[locLH] )
 #
 return( resid )
 #
} # end of residmds
#
#
residmds2 <- function( vecX ){
 # residual of metric mds:  vech(D-O) to be used by nlfb
 # Shin-ichi Mayeakawa
 # 20210713,20220713,14,22
 #
 return( -residmds( vecX ) )
} # end of residmds2
#
#
rssmds <- function( vecX ){
 # residual of metric mds:  ssq(vech(O-D))
 # Shin-ichi Mayeakawa
 # 20210713,20220713,14,22
 #
 return( sum( residmds(vecX)^2 ) )
} # end of rssmds
#
#
Jacmds <- function( vecX ){
 # Numeric Jacobian of vech(O-D), i.e. -vech(D)
 # Shin-ichi Mayeakawa
 # 20210713,20220713,14,22
 F=JacobianMat( vecX, residmds )
 return( F )
} # end of gradmds
#
#
Jacmds_a <- function( vecX ){
 # Analytic Jacobian of vech(D)
 # Shin-ichi Mayeakawa
 # 20220723
 # convert vecX to matrix X
 X=matrix(vecX,ns)
 # Euclidean distance
 XX=rowSums(X*X)
 XX11=XX%*%matrix(1,1,ns)
 D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
 vecD=D[locLH]
 #
 # This is the Jacobian of vech(D) wrt vec(X)
 F=matrix(0,npair,ns*ndim)
 #
 if(0){
  k=0
  for( i in 2:ns ){
   for( j in 1:(i-1) ){
    k=k+1
    for( a in 1:ndim ){
     F[k,(a-1)*ns+i]=(X[i,a]-X[j,a])/D[i,j]
     F[k,(a-1)*ns+j]=-(X[i,a]-X[j,a])/D[i,j]
    }
   }
  }
 } # en of skip
 #
 GX=G%*%X
 for( k in 1:npair ){
  for( a in 1:ndim ){
   ximxja=GX[k,a]
   F[k,(a-1)*ns+locvLH[k,1]]=ximxja/vecD[k]
   F[k,(a-1)*ns+locvLH[k,2]]=-ximxja/vecD[k]
  }
 }
 #
 return( F )
 #
} # end of Jacmds_a
#
#
HessQmds_a <- function( vecX ){
 # approximate Hessian calculated from analytic Jacobian of vech(D)
 # Shin-ichi Mayeakawa
 # 20220722,23
 F=Jacmds_a( vecX )
 H=t(F)%*%F
 H=H + 0.01*min(abs(diag(H)))*diag(nrow(H))
 return( 2*H )
} # end of HessQmds_a
#
#
gradmds <- function( vecX ){
 # numeric gradient of rssmds
 # Shin-ichi Mayeakawa
 # 20210713,20220713,14,22,23
 res=JacobianMat( vecX, rssmds )
 return( res )
} # end of gradmds
#
#
gradmds_a <- function( vecX ){
 # analytic gradient of rssmds
 # Shin-ichi Mayeakawa
 # 20220723
 res=c( -2*residmds(vecX)%*%Jacmds_a(vecX) )
 return( res )
} # end of gradmds_a
#
#
normalizemds <- function( vecX ){
 # normalize X in each of GN/NR iteration
 # Shin-ichi Mayekawa
 # 20220723
 #
 # convert vecX to matrix X
 X=matrix(vecX,ns)
 # column center X
 return( c( normalize_config(X,center=1) ))
} # end of normalizemds
#
#
#
#
# generate test data
#
# X=as.matrix( expand.grid(-3:3, -3:3) )
X=as.matrix( expand.grid(-2:2, -2:2) )
# X=as.matrix( expand.grid(-1:1,-1:1) )
ns=nrow(X); ndim=ncol(X)
sname=paste("s",1:ns,sep="")
#
set.seed(1701)
sigmaE=0.5
XX=rowSums(X*X)
XX11=XX%*%matrix(1,1,ns)
D=sqrt( abs( XX11-2*X%*%t(X)+t(XX11) ) )
O=D+sigmaE*rnorm(ns*ns)
O=round( 0.5*(O+t(O))-diag(diag(O)), 2 )
locLH=lower.tri(O, diag=FALSE)
npair=sum(locLH)
# locvLH=vechindex(ns,1,0,type=2)
locvLH=cbind(row(O)[locLH],col(O)[locLH])
# G matrix for  paired comparison
G=matrix(0,npair,ns)
for( k in 1:npair ){
 G[k,locvLH[k,1]]=1; G[k,locvLH[k,2]]=-1
}
#
#
#
# # of dimensions
ndim=2
#
#
# initial by classic mds
J=diag(ns)-matrix(1,ns,ns)/ns
B=-0.5*J%*%O^2%*%J
ev=eigen(B)
P=ev$vectors[,1:ndim]
L=diag(ev$values[1:ndim])
Xcmds=P%*%sqrt(L)
#
X0=Xcmds
rss0=rssmds( c(X0) )
g0=gradmds( c(X0) )
title=paste("CMDS: rss =", round(rss0,3))
plot(X0, main=title)
text( X0, sname, pos=4 )
#
#
#
# Gauss-Newton
#
# by nlfb
resnlfb=nlfb(c(X0),residmds2, control=list(femax=1000,japprox="jacentral"))
X1=matrix(resnlfb$coefficients,ns)
rss1=rssmds( c(X1) )
g1=gradmds( c(X1) )
title=paste("GN by nlfb: rss =", round(rss1,3))
plot(X1, main=title)
text( X1, sname, pos=4 )
#
#
# by GN
resGN=GN( c(X0), residmds, maxiter=100, normalize=normalizemds )
X2=matrix(resGN$par,ns)
rss2=rssmds( c(X2) )
g2=gradmds( c(X2) )
title=paste("GN by GN: rss =", round(rss2,3))
plot(X2, main=title)
text( X2, sname, pos=4 )
#
#
# by GN with analytic gradient
resGN2=GN( c(X0), residmds, maxiter=100, normalize=normalizemds
           , jacobian=Jacmds_a )
X21=matrix(resGN2$par,ns)
rss21=rssmds( c(X21) )
g21=gradmds( c(X21) )
title=paste("GN by GN with analytic gradient: rss =", round(rss2,3))
plot(X21, main=title)
text( X21, sname, pos=4 )
#
Print(rss0,rss1,rss2, rss21)
#
#
#
#
# Newton-Raphson
#
# by nlminb
resnlminb=nlminb( c(X0), rssmds, control=list(trace=0) )
X3=matrix(resnlminb$par,ns)
rss3=rssmds( c(X3) )
g3=gradmds( c(X3) )
title=paste("NR by nlminb: rss =", round(rss3,3))
plot(X3, main=title)
text( X3, sname, pos=4 )
#
#
# by NR
resNR=NR( c(X0), rssmds, maxiter=1000, flipd=1, normalize=normalizemds )
X4=matrix(resNR$par,ns)
rss4=rssmds( c(X4) )
g4=gradmds( c(X4) )
title=paste("NR by NR: rss =", round(rss4,3))
plot(X4, main=title)
text( X4, sname, pos=4 )
#
Print(rss0, rss1, rss2, rss3, rss4)
#
#
# by NR by supplying analytic gradient and quasi Hessian.
resNR2=NR( c(X0), rssmds, maxiter=1000, flipd=1
           , normalize=normalizemds, grad=gradmds_a, hessian=HessQmds_a )
#
#
X41=matrix(resNR2$par,ns)
rss41=rssmds( c(X41) )
g41=gradmds( c(X41) )
title=paste("NR by NR with quasi H: rss =", round(rss41,3))
plot(X41, main=title)
text( X41, sname, pos=4 )
#
Print(rss1, rss2, rss3, rss4, rss41)
#
#
#
#
if(0){
 #
 library(lazy.mds)
 #
 # wmdsGN
 resmds=wmdsGN( O, ndim=2, mlevel=9, maxiter=1000, print=1 )
 Xw=resmds$X
 rssw=rssmds( c(Xw) )
 gw=gradmds( c(Xw) )
 title=paste("GV by wmdsGN: rss =", round(rssw,3))
 plot(Xw, main=title)
 text( Xw, sname, pos=4 )
 #
 Print(rss0, rss1, rss2, rss3, rss4, rssw)
 #
 Print(max(abs(g0)), max(abs(g1)),max(abs(g2))
       ,max(abs(g3)),max(abs(g4)),max(abs(gw)))
 #
 #
}
#
#
#
#####################################################################
#
#                End of Multidimensional Scaling
#
#####################################################################
#
#
#
#
#
#
#
#
#
#####################################################################
#
#                Factor Analysis
#
#####################################################################
#
#
#
critfa <- function( param, logpsi=0, method="ML", LH=0 ){
 # evaluation of ML or LS criterion function: param=c( c(Lambda), psi )
 # Shin-ichi Mayekawa
 # 20210705,20220723
 # uselogpsiasparam -> logpsi: 20220728
 #
 #  S, logdets, ndim and nvar are from the parent env.
 #
 #
 Lambda=matrix(param[1:(nvar*ndim)],nvar)
 psi=param[-(1:(nvar*ndim))]
 if( logpsi ) psi=exp(psi)
 Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
 #
 if( method == "ML" ){
  invSigma=solve(Sigma)
  crit=tr(S%*%invSigma)+log(det(Sigma)) - logdetS - nvar
 }
 else{
  if( LH == 1 ) crit=sum((S-Sigma)[locLH]^2)
  else crit=sum((S-Sigma)^2)
 }
 #
 return( crit )
 #
} # end of critfa
#
#
#
#
critfapsi <- function( psi, logpsi=0, method="ML", LH=0 ){
 # evaluation of ML or LS criterion function: param=psi
 # Shin-ichi Mayekawa
 # 20210707
 # uselogpsiasparam -> logpsi: 20220728
 #
 #  S, logdets, ndim and nvar are from the parent env.
 #
 #
 # psi or log(psi)
 # If the parameter is log(psi), recover original as exp(log(psi)).
 if( logpsi ) psi=exp(psi)
 #
 psi=c(psi)
 #
 if( method == "ML" ){
  # conditional MLE of Lambda given psi
  Psimh=diag(1/sqrt(psi))
  ev=eigen( Psimh%*%S%*%Psimh )
  P=ev$vectors[,1:ndim]
  D=ev$values[1:ndim]
  Lambda=sqrt(psi)*P%*%diag(sqrt(D-1), nrow=ndim)%*%invPhimh
  # log likelihood ratio stat
  Sigma=Lambda%*%t(Lambda)+diag(psi)
  invSigma=solve(Sigma)
  crit=tr(S%*%invSigma)+log(det(Sigma)) - logdetS - nvar
 }
 else{
  # conditional LSE of Lambda given psi
  ev=eigen( S-diag(psi) )
  P=ev$vectors[,1:ndim]
  Gamma=ev$values[1:ndim]
  Lambda=P%*%diag(sqrt(Gamma))
  # rss
  Sigma=Lambda%*%t(Lambda)+diag(psi)
  if( LH == 1 ) crit=sum((S-Sigma)[locLH]^2)
  else crit=sum((S-Sigma)^2)
 }
 #
 return( crit )
 #
} # end of critfapsi
#
#
#
#
gradfa_a <- function(  param, logpsi=0, method="ML", LH=0 ){
 # analytic gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
 # Shin-ichi Mayekawa
 # 20220724,25
 # avoid Jacfa_a: 20220726dnc
 # correction for rss_LH: 20220726dnc,0728
 # uselogpsiasparam -> logpsi: 20220728
 #
 #
 Lambda=matrix(param[1:(nvar*ndim)],nvar)
 psi=param[-(1:(nvar*ndim))]
 if( logpsi ) psi=exp(psi)
 Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
 #
 if( method == "LS" ){
  SS=S-Sigma
  # dLambda=0.5*(-4*diag(diag(SS))%*%Lambda -4*SS%*%Lambda )
  if( LH == 0 ) dLambda=-4*SS%*%Lambda%*%Phi
  else dLambda=-2*((diag(diag(SS))+SS)%*%Lambda%*%Phi )
  dpsi=-2*diag(SS)
  if( logpsi == 1 ) dpsi=dpsi*psi
  g=c( c(dLambda), dpsi )
 }
 else{
  invSigma=matSwp(Sigma)
  iSSmSiS=invSigma%*%(Sigma-S)%*%invSigma
  dLambda=2*iSSmSiS%*%Lambda%*%Phi
  dpsi=diag(iSSmSiS)
  if( logpsi == 1 ) dpsi=dpsi*psi
  g=c(c(dLambda),dpsi)
 }
 #
 return( g )
 #
} # gradfa_a
#
#
gradfa <- function( param, logpsi=0, method="ML" ){
 # numerical gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
 # Shin-ichi Mayekawa
 # 20220724,25
 # uselogpsiasparam -> logpsi: 20220728
 #
 return( JacobianMat( param, critfa, method=method
                      , logpsi=logpsi ) )
} # end of gradfa
#
#
#
HessQfa_a <- function( param, logpsi=0, method="LS", LH=0 ){
 # analytic quasi Hessian of LS or ML critfa wrt param=c(c(Lambda,psi)
 # Shin-ichi Mayekawa
 # 20220725
 # LH: 20220728
 # uselogpsiasparam -> logpsi: 20220728
 #
 if( method == "ML" ){
  cat("\nerror(HessQfa) This is for method=LS only.\n")
  stop()
 }
 F=Jacfa_a( param, logpsi=logpsi, LH=LH )
 H=t(F)%*%F
 H=H + 0.01*min(abs(diag(H)))*diag(nrow(H))
 return( 2*H )
} # end of HessQfa_a
#
#
#
gradfapsi <- function( psiorlogpsi, logpsi=0, method="ML" ){
 # numerical gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
 # Shin-ichi Mayekawa
 # 20220726
 # uselogpsiasparam -> logpsi: 20220728
 #
 # This works with  LH=0 and LH=1
 #
 return( JacobianMat( psiorlogpsi, critfapsi, method=method
                      , logpsi=logpsi ) )
} # end of gradfapsi
#
#
gradfapsi_a <- function( psiorlogpsi, logpsi=0, method="ML" ){
 # analytic gradient of LS or ML critfapsi wrt param=psi or log(psi)
 # Shin-ichi Mayekawa
 # 20220726
 # uselogpsiasparam -> logpsi: 20220728
 #
 # This is for LH=0
 #
 #
 # psi or log(psi)
 # If the parameter is log(psi), recover original as exp(log(psi)).
 if( logpsi ) psi=exp(psiorlogpsi)
 else psi=psiorlogpsi
 psi=c(psi)
 #
 if( method == "ML" ){
  # conditional MLE of Lambda given psi
  Psimh=diag(1/sqrt(psi))
  ev=eigen( Psimh%*%S%*%Psimh )
  Omega2=ev$vectors[,(ndim+1):nvar]^2
  theta=ev$values[(ndim+1):nvar]
  o2t=Omega2%*%diag(1-theta)
  g=c( rowSums(o2t)/psi )
  if( logpsi == 1 ) g=g*psi
 }
 else{
  # conditional LSE of Lambda given psi
  ev=eigen( S-diag(psi) )
  P2=ev$vectors[,(ndim+1):nvar]^2
  gamma=ev$values[(ndim+1):nvar]
  g=c(-2*rowSums(P2%*%diag(gamma)) )
  # above is correct if rss=ssq(S-Sigma), not if rss=ssq((S-Sigma)[locLH])
  if( logpsi == 1 ) g=g*psi
 }
 #
 return(g)
 #
} # end of gradfapsi_a
#
#
#
#
#
#
residfa <- function( param, logpsi=0 ){
 # residual of LS fa vec(S-Sigma)
 # Shin-ichi Mayekawa
 # 20220728
 # uselogpsiasparam -> logpsi: 20220728
 #
 #
 Lambda=matrix(param[1:(nvar*ndim)],nvar)
 psi=param[-(1:(nvar*ndim))]
 if( logpsi ) psi=exp(psi)
 Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
 #
 return( c(S-Sigma) )
 #
} # end of residfa
#
#
residfa2 <- function( param, logpsi=0 ){
 # residual of LS fa vec(Sigma-S)
 # Shin-ichi Mayekawa
 # 20220728
 # uselogpsiasparam -> logpsi: 20220728
 #
 #
 return( -residfa( param, logpsi=logpsi ) )
 #
} # end of residfa2
#
#
#
Jacfa_a <- function( param, logpsi=0, LH=0 ){
 # analytic Jacobian of vech(Sigma) or vec(Sigma) wrt param=c(c(Lambda,psi)
 # Shin-ichi Mayekawa
 # 20220724,25
 # LH=0: 20220728
 # uselogpsiasparam -> logpsi: 20220728
 #
 #
 Lambda=matrix(param[1:(nvar*ndim)],nvar)
 psi=param[-(1:(nvar*ndim))]
 if( logpsi ) psi=exp(psi)
 Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
 #
 #
 if( LH == 1 ){
  npair=nvar*(nvar-1)/2
  # First, calculate  d vec(Sigma) / d vec(Lambda')
  F=matrix(0,npair,nvar*ndim+nvar)
  #
  for( k in 1:npair ){
   i=locvLH[k,1]; j=locvLH[k,2]
   if( i == j ){
    F[k,ndim*(i-1)+(1:ndim)]=2*Lambda[i,]
    if( logpsi == 0 ) F[k,ndim*nvar+i]=1
    else F[k,ndim*nvar+i]=psi[i]
   }
   else{
    F[k,ndim*(i-1)+(1:ndim)]=Lambda[j,]
    F[k,ndim*(j-1)+(1:ndim)]=Lambda[i,]
   }
  }
 } # end of LH
 else{
  npair=nvar*nvar
  # First, calculate  d vec(Sigma) / d vec(Lambda')
  F=matrix(0,npair,nvar*ndim+nvar)
  #
  for( k in 1:npair ){
   # col-major rollout
   j=trunc((k-1)/nvar)+1; i=((k-1)%%nvar)+1
   if( i == j ){
    F[k,ndim*(i-1)+(1:ndim)]=2*Lambda[i,]
    if( logpsi == 0 ) F[k,ndim*nvar+i]=1
    else F[k,ndim*nvar+i]=psi[i]
   }
   else{
    F[k,ndim*(i-1)+(1:ndim)]=Lambda[j,]
    F[k,ndim*(j-1)+(1:ndim)]=Lambda[i,]
   }
  }
 }  # end of LH!=0
 #
 # F0=matrix(as.numeric(abs(F)>0),npair)
 #
 # swap the columns
 #  d vec(S) / d vec(Lambda') -> d vec(S) / d vec(Lambda)
 # od=c(t(matrix(1:(nvar*ndim),ndim,nvar)))]
 od=matindex(ndim,nvar, type=1, byrow=1)
 F[,1:(nvar*ndim)]=F[,od]
 #
 # F1=matrix(as.numeric(abs(F)>0),npair)
 #
 return( F )
 #
} # end of Jacfa_a
#
#
#
Jacfa_aa <- function( param, logpsi=0, LH=0 ){
 # returns analytic Jacobian of vec(S-Sigma) as gradient attribute
 # Shin-ichi Mayekawa
 # 20220728
 # uselogpsiasparam -> logpsi: 20220728
 #
 res=Jacfa_a( param, logpsi=logpsi, LH=LH )
 res0=0
 attributes(res0)$gradient=res
 return( res0 )
} # end of Jacfa_aa
#
#
gradfa_Ja <- function(  param, logpsi=0, method="ML", LH=0 ){
 # analytic gradient of LS or ML critfa wrt param=c(c(Lambda,psi)
 # If method="LS" Jacfa_a is used.
 # Shin-ichi Mayekawa
 # 20220724,25
 # LH: 20220728
 # uselogpsiasparam -> logpsi: 20220728
 # Should be the same as gradfa_a
 #
 #
 Lambda=matrix(param[1:(nvar*ndim)],nvar)
 psi=param[-(1:(nvar*ndim))]
 if( logpsi ) psi=exp(psi)
 Sigma=Lambda%*%Phi%*%t(Lambda)+diag(psi)
 #
 if( method == "LS" ){
  # by Jacobian:  applicable for vec(S) and vech(S)
  F=Jacfa_a( param, logpsi=logpsi, LH=LH )
  # gradient
  g=c( -2*t(F)%*%(S-Sigma)[locLH] )
 }
 else{
  invSigma=matSwp(Sigma)
  iSSmSiS=invSigma%*%(Sigma-S)%*%invSigma
  dLambda=2*iSSmSiS%*%Lambda%*%Phi
  dpsi=diag(iSSmSiS)
  if( logpsi == 1 ) dpsi=dpsi*psi
  g=c(c(dLambda),dpsi)
 }
 #
 return( g )
 #
} # gradfa_Ja
#
#
#
#
#
#
#
#
# generate data
seed <- 1701
set.seed(seed)
#
nvar <- 20;  ndim0 <- 3
ps <- 0.2
df <- 500
#
A <- matrix(runif(nvar*ndim0),nvar)
Sigma <- A%*%t(A)
Sigma <- Sigma+ps*diag(nvar)
dS <- sqrt(diag(Sigma))
Sigma <- diag(1/dS)%*%Sigma%*%diag(1/dS)
#
# correlation matrix from Wishart random matrix
S <- rWishart( 1, df, Sigma )
S <- S[,,1]/df
dS <- sqrt(diag(S))
S <- diag(1/dS)%*%S%*%diag(1/dS)
#
# const to be used later
#
# log determinant of S
logdetS=log(det(S))
#
# lower half of S:   NOT used.
locLH=lower.tri(S,diag=TRUE)
# vecS=S[locLH]
npair=sum(locLH)
locvLH=cbind(row(S)[locLH],col(S)[locLH])
#
#
#
# # of dimensions to be used
ndim=2
#
#
#
# initial values by pca
ev=eigen( S )
Lambda0=ev$vectors[,1:ndim]%*%diag(sqrt(ev$values[1:ndim]), nrow=ndim)
psi0=pmax(0.001, diag(S-Lambda0%*%t(Lambda0)) )
#
# factor correlation matrix (constant)
Phi=diag(ndim)
invPhi=diag(ndim)
invPhimh=diag(ndim)
#
#
# initial value
param0=c(c(Lambda0),psi0)
param00=c(c(Lambda0),log(psi0))
#
#
#
#
#
#         Maximum Likelihood Solution
#
#
#
#
#
#
# ML fa by nlminb: param=c( c(Lambda),psi )
resnlminb=nlminb( param0, critfa, control=list(trace=0)
                  , lower=1e-3, upper=1-1e-3 )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
# ML fa by nlminb with analytic gradient: param=c( c(Lambda),psi )
resnlminb=nlminb( param0, critfa, control=list(trace=0)
                  , lower=1e-3, upper=1-1e-3, gradient=gradfa_a )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
#
# ML fa by NR:  param=c( c(Lambda),psi )
resNR=NR( param0, critfa, logpsi=0, maxiter=50, flipd=1 )
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb,critNR)
#
# ML fa by NR with analytic gradient:  param=c( c(Lambda),psi )
resNR=NR( param0, critfa, logpsi=0, maxiter=50, flipd=1
          , gradient=gradfa_a)
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb,critNR)
#
#
#
#
# ML fa by nlminb: param=c( c(Lambda),log(psi) )
resnlminb=nlminb( param00, critfa, control=list(trace=0)
                  , logpsi=1 )
critnlminb=resnlminb$objective
param=resnlminb$par
LambdaML=matrix(param[1:(nvar*ndim)],nvar)
psiML=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb)
#
# ML fa by nlminb with analytic gradient: param=c(c(Lambda),log(psi))
resnlminb=nlminb( param00, critfa, control=list(trace=0)
                  , logpsi=1, gradient=gradfa_a )
critnlminb=resnlminb$objective
param=resnlminb$par
LambdaML=matrix(param[1:(nvar*ndim)],nvar)
psiML=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb)
#
#
# ML fa by NR: param=c( c(Lambda),log(psi) )
resNR=NR( param00, critfa, logpsi=1, maxiter=50, flipd=1 )
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb,critNR)
#
#
# ML fa by NR with analytic gradient: param=c( c(Lambda),log(psi) )
resNR=NR( param00, critfa, logpsi=1, maxiter=50, flipd=1
          , gradient=gradfa_a )
critNR=resNR$objective
param=resNR$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlminb,critNR)
#
#
#
#
#
#
#
# ML fa by nlminb:  param=psi
# with constraints on psi. (no effect)
resnlminb=nlminb( psi0, critfapsi, control=list(trace=0)
                  , lower=0.02, upper=0.98 )
critnlminb=resnlminb$objective
psinlminb=resnlminb$par
Print(critnlminb,psinlminb)
#
g1=gradfapsi( resnlminb$par )
g2=gradfapsi_a( resnlminb$par )
Print(g1,g2)
#
#
# ML fa by nlminb:  param=log(psi)
resnlminb=nlminb( log(psi0), critfapsi, control=list(trace=0)
                  , logpsi=1 )
critnlminb=resnlminb$objective
psinlminb=exp( resnlminb$par )
Print(critnlminb,psinlminb)
#
g1=gradfapsi( resnlminb$par, logpsi=1 )
g2=gradfapsi_a( resnlminb$par, logpsi=1 )
Print(g1,g2)
#
#
# ML fa by NR:  param=psi
resNR=NR( psi0, critfapsi, maxiter=50, flipd=1 )
critNR=resNR$objective
psiNR=resNR$par
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par )
g2=gradfapsi_a( resNR$par )
Print(g1,g2)
#
#
# ML fa by NR with analytic gradient:  param=psi
resNR=NR( psi0, critfapsi, maxiter=50, flipd=1, gradient=gradfapsi_a )
critNR=resNR$objective
psiNR=resNR$par
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par )
g2=gradfapsi_a( resNR$par )
Print(g1,g2)
#
#
#
# ML fa by NR:  param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50, flipd=1 )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1 )
g2=gradfapsi_a( resNR$par, logpsi=1 )
Print(g1,g2)
#
#
# ML fa by NR with analytic gradient:  param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50, flipd=1
          , gradient=gradfapsi_a )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1 )
g2=gradfapsi_a( resNR$par, logpsi=1 )
Print(g1,g2)
#
#
Print(critnlminb,critNR)
#
#
#
#
#
#
#    Least Squares Solution by Newton-Raphson Method
#
#
#
#
#
#
# LS fa by nlminb:  param=c(c(Lambda),psi)
resnlminb=nlminb( param0, critfa, control=list(trace=0)
                  , logpsi=0, method="LS" )
critnlminb=resnlminb$objective
param=resnlminb$par
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlminb)
#
g1=gradfa( param, method="LS" )
g2=gradfa_a( param, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlminb:  param=psi
resnlminb=nlminb( psi0, critfapsi, control=list(trace=0)
                  , method="LS" )
critnlminb=resnlminb$objective
Print(critnlminb,psinlminb)
g1=gradfapsi( resnlminb$par, method="LS" )
g2=gradfapsi_a( resnlminb$par, method="LS" )
Print(g1,g2)
#
#
#
# LS fa by nlminb with reparametrization:  param=log(psi)
resnlminb=nlminb( log(psi0), critfapsi, control=list(trace=0)
                  , logpsi=1, method="LS" )
critnlminb=resnlminb$objective
psinlminb=exp( resnlminb$par )
Print(critnlminb,psinlminb)
g1=gradfapsi( resnlminb$par, logpsi=1, method="LS" )
g2=gradfapsi_a( resnlminb$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by NR with analytic g and H:  param=c(c(Lambda),psi)
resNR=NR( param0, critfa, logpsi=0, maxiter=500
          , flipd=1, method="LS", gradient=gradfa_a, hessian=HessQfa_a )
critNR=resNR$objective
param=resNR$par
Print(critNR,critnlminb)
#
g1=gradfa( resNR$par, logpsi=0, method="LS" )
g2=gradfa_a( resNR$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
#
# LS by NR with analytic g and H:
#        param=c(c(Lambda),log(psi))
resNR=NR( param00, critfa, logpsi=1, maxiter=500
          , flipd=1, method="LS", gradient=gradfa_a, hessian=HessQfa_a )
critNR=resNR$objective
param=resNR$par
Print(critNR,critnlminb)
#
g1=gradfa( resNR$par, logpsi=1, method="LS" )
g2=gradfa_a( resNR$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by NR:  param=log(psi)
resNR=NR( log(psi0), critfapsi, logpsi=1, maxiter=50
          , flipd=1, method="LS" )
critNR=resNR$objective
psiNR=exp( resNR$par )
Print(critNR,critnlminb,psiNR)
#
g1=gradfapsi( resNR$par, logpsi=1, method="LS" )
g2=gradfapsi_a( resNR$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
#
#
#    Least Squares Solution by Gauss-Newton Method
#
#
#
#
# LS fa by nlfb  param=c(c(Lambda),psi)
resnlfb=nlfb( param0,residfa2,control=list(femax=1000,japprox="jacentral"))
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlfb)
#
g1=gradfa( param, logpsi=0, method="LS" )
g2=gradfa_a( param, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlfb with analytic jacobian:  param=c(c(Lambda),psi)
resnlfb=nlfb( param0, residfa2, jacfn=Jacfa_aa, control=list(femax=1000) )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critnlfb)
#
g1=gradfa( param, logpsi=0, method="LS" )
g2=gradfa_a( param, logpsi=0, method="LS" )
Print(g1,g2)
#
#
#
# LS fa by nlfb  param=c(c(Lambda),log(psi))
resnlfb=nlfb( param00, residfa2
              , control=list(femax=1000,japprox="jacentral"), logpsi=1 )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlfb)
#
g1=gradfa( param, logpsi=1, method="LS" )
g2=gradfa_a( param, logpsi=1, method="LS" )
Print(g1,g2)
#
#
# LS fa by nlfb with analytic jacobian:  param=c(c(Lambda),log(psi))
resnlfb=nlfb( param00, residfa2, control=list(femax=1000)
              , logpsi=1, jacfn=Jacfa_aa )
param=resnlfb$coefficients
critnlfb=resnlfb$ssquares
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critnlfb)
#
g1=gradfa( param, logpsi=1, method="LS" )
g2=gradfa_a( param, logpsi=1, method="LS" )
Print(g1,g2)
#



# LS fa by GN: param=c(c(Lambda),psi)
resGN=GN( param0, residfa, logpsi=0, maxiter=50 )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=0, method="LS" )
g2=gradfa_a( resGN$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by GN with analytic jacobian: param=c(c(Lambda),psi)
resGN=GN( param0, residfa, logpsi=0, maxiter=50
          , jacobian=Jacfa_a )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=param[-(1:(nvar*ndim))]
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=0, method="LS" )
g2=gradfa_a( resGN$par, logpsi=0, method="LS" )
Print(g1,g2)
#
#
# LS fa by GN: param=c(c(Lambda),log(psi))
resGN=GN( param00, residfa, logpsi=1, maxiter=50 )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=1, method="LS" )
g2=gradfa_a( resGN$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#
#
# LS fa by GN with analytic jacobian: param=c(c(Lambda),log(psi))
resGN=GN( param00, residfa, logpsi=1, maxiter=50, jacobian=Jacfa_a )
param=resGN$par
critGN=resGN$objective
Lambda=matrix(param[1:(nvar*ndim)],nvar)
psi=exp( param[-(1:(nvar*ndim))] )
Print(critGN, critnlfb)
#
g1=gradfa( resGN$par, logpsi=1, method="LS" )
g2=gradfa_a( resGN$par, logpsi=1, method="LS" )
Print(g1,g2)
#
#
#

#####################################################################
#
#                End of Factor Analysis
#
#####################################################################
#
#
#
#
#
#
#
#
#
#
#####################################################################
#
#       PCA by Newton-Raphson and Gauss-Newton method
#
#####################################################################
#
#
#
# Given n x nvar Y matrix, find n x ndim F and nvar x ndim A such that
#  rss = tr( t(Y - F \
# is minimized where ndim << min(n,nvar).
#
# The solution using svd or eigen decomposition is well known
# but here we try to minimize rss directly w.r.t c( c(F), c(A) )
# by Newton-Raphson or Gauss-Newton method.
#
# It's not efficient at all but works.
#
#
#
# We use nlfb function of nlsr package for comparison.
#
library(nlsr)
#
#
rssPCA <- function( vecPCA, ndim=2, Y=0 ){
 # rss of PCA: ssq(Y-Yhat)=ssq(Y-F%*%t(A))
 # recover F and A from vecPCA and evaluate rss
 n=nrow(Y); nvar=ncol(Y)
 F=matrix(vecPCA[1:(n*ndim)],n,ndim)
 A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
 rss=ssq(Y-F%*%t(A))
 return( rss )
} # end of rssPCA
#
#
drssPCA <- function( vecPCA, ndim=2, Y=0 ){
 # gradient of rss of PCA
 # recover F and A from vecPCA and evaluate rss
 n=nrow(Y); nvar=ncol(Y)
 F=matrix(vecPCA[1:(n*ndim)],n,ndim)
 A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
 dF=2*(F%*%t(A)%*%A-Y%*%A)
 dA=2*(A%*%t(F)%*%F-t(Y)%*%F)
 res=c(c(dF),c(dA))
 return( res )
} # end of drssPCA
#
#
ymyhatPCA <- function( vecPCA, ndim=2, Y=0 ){
 # residual as c(Y-Yhat)
 # recover F and A from vecPCA and evaluate vec(Y-F%*%t(A))
 n=nrow(Y); nvar=ncol(Y)
 F=matrix(vecPCA[1:(n*ndim)],n,ndim)
 A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
 resid=c(Y-F%*%t(A))
 return( resid )
} # end of ymyhatPCA
#
#
yhatmyPCA <- function( vecPCA, ndim=2, Y=0 ){
 # residual as c(Yhat-Y) to be used by nlfb
 # recover F and A from vecPCA and evaluate vec(F%*%t(A)-Y)
 return( -ymyhatPCA( vecPCA, ndim=ndim, Y=Y ) )
} # end of yhatmyPCA
#
#
JacPCA <- function( vecPCA, ndim=2, Y=0 ){
 # Jacobian: c( vec( d c(Yhat) / d F ), vec( d c(Yhat) / d A ) )
 # recover F and A from vecPCA and evaluate vec(Y-F%*%t(A))
 n=nrow(Y); nvar=ncol(Y)
 F=matrix(vecPCA[1:(n*ndim)],n,ndim)
 A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
 #
 # Jac
 temp=NULL
 for( i in 1: n ) temp=append(temp,list(A))
 temp=Reduce( b_diag, temp )
 # dF=temp[,c(t(matrix(1:(ndim*n),ndim,n)))]
 # dF=dF[c(t(matrix(1:(nvar*n),nvar,n))),]
 dF=temp[,matindex(ndim,n, type=1, byrow=1)]
 dF=dF[matindex(nvar,n, type=1, byrow=1),]
 temp=NULL
 for( j in 1: nvar ) temp=append(temp,list(F))
 temp=Reduce( b_diag, temp )
 # dA=temp[,c(t(matrix(1:(ndim*nvar),ndim,nvar)))]
 dA=temp[,matindex(ndim,nvar, type=1, byrow=1)]
 #
 res=cbind( dF, dA )
 #
 return( res )
} # end of JacPCA
#
#
#
normalizePCA <- function( vecPCA, ndim=2, Y=0 ){
 # normalize X in each of GN/NR iteration
 # so that F'F=n I  and A'A=diag and A[1,] > 0
 # 20220723,30
 #
 # convert vecX to matrix X
 n=nrow(Y); nvar=ncol(Y)
 F=matrix(vecPCA[1:(n*ndim)],n,ndim)
 A=matrix(vecPCA[-(1:(n*ndim))],nvar,ndim)
 #
 # normalize
 res=normalize_config(F,center=0, rotate=1, A=A)
 #
 F=res$F; A=res$A
 for( a in 1:ndim ){
  if( A[1,a] < 0 ){
   A[,a]=-A[,a]
   F[,a]=-F[,a]
  }
 }
 res=list(F=F,A=A)
 #
 return( unlist(res) )
} # end of normalizePCA
#
#
#
#
# generate data
seed <- 1701
set.seed(seed)
n <- 50; nvar <- 9
Y <- matrix(rnorm(n*nvar),n)%*%(diag(1/(1:nvar)))
# centering
Y=scale( Y, center=TRUE, scale=FALSE )
#
#
#
# initial value
ndim=2
#
set.seed(1701)
F0=matrix(rnorm(n*ndim),n,ndim)
A0=t( solve(t(F0)%*%F0)%*%t(F0)%*%Y )
vecPCA0=c(c(F0),c(A0))
rss0=ssq(Y-F0%*%t(A0))
#
#
#
#
# Estimation starts here.
#
# Using all of F and A as parameters
# normalizing the solution after each NR iteration.
# with or without analytic gradient or Jacobian.
#
#
#
# Newton-Raphson
#
#
#
# by nlminb
cat("\n\n**** Starting PCA by NR by nlminb with analytic gradient ****\n\n")
resnlm=nlminb( vecPCA0, rssPCA, control=list(trace=0), Y=Y, ndim=ndim
               , gradient=drssPCA )
vecPCAnlm=resnlm$par
F=matrix(vecPCAnlm[1:(n*ndim)],n,ndim)
A=matrix(vecPCAnlm[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
Fnlm=resnorm$F; Anlm=resnorm$A
rssnlm=ssq(Y-Fnlm%*%t(Anlm))
gradnlm=drssPCA( vecPCAnlm, ndim=ndim, Y=Y )
Print(rss0,rssnlm, min(abs(gradnlm)))
#Print(A1,Anlm)
#
#
# by NR
cat("\n\n**** Starting PCA by NR by NR with analytic gradient **** \n\n")
resNR=NR( vecPCA0, rssPCA, Y=Y, ndim=ndim, normalize=normalizePCA
          , gradient=drssPCA,  flipd=1, maxiter2=3, eps=-9, print=2 )
vecPCANR=resNR$par
F=matrix(vecPCANR[1:(n*ndim)],n,ndim)
A=matrix(vecPCANR[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FNR=resnorm$F; ANR=resnorm$A
rssNR=ssq(Y-FNR%*%t(ANR))
gradNR=drssPCA( vecPCANR, ndim=ndim, Y=Y )
Print(rss0,rssNR, min(abs(gradNR)))
# Print(A1,ANR)
#
#
#
# Gauss-Newton
#
#
#
# by nlsr
cat("\n\n**** Startin PCA by GN by nlfb ****\n\n")
resnlsr=nlfb( vecPCA0, yhatmyPCA, Y=Y, ndim=ndim
              , control=list(femax=500,japprox="jacentral") )
vecPCAnlsr=resnlsr$coefficients
F=matrix(vecPCAnlsr[1:(n*ndim)],n,ndim)
A=matrix(vecPCAnlsr[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
Fnlsr=resnorm$F; Anlsr=resnorm$A
rssnlsr=ssq(Y-Fnlsr%*%t(Anlsr))
gradnlsr=drssPCA( vecPCAnlsr, ndim=ndim, Y=Y )
Print(rss0,rssnlsr, min(abs(gradnlsr)))
#
#
#
# by GN
cat("\n\n**** Starting PCA by GN  w/o analytic Jacobian ****\n\n")
resGN=GN( vecPCA0, ymyhatPCA, Y=Y, ndim=ndim, normalize=normalizePCA
          , print=2 )
vecPCAGN=resGN$par
F=matrix(vecPCAGN[1:(n*ndim)],n,ndim)
A=matrix(vecPCAGN[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FGN=resnorm$F; AGN=resnorm$A
rssGN=ssq(Y-FGN%*%t(AGN))
gradGN=drssPCA( vecPCAGN, ndim=ndim, Y=Y )
Print(rss0,rssGN, min(abs(gradGN)))
#
#
# by GN
cat("\n\n**** Starting PCA by GN  with analytic Jacobian ****\n\n")
resGN=GN( vecPCA0, ymyhatPCA, Y=Y, ndim=ndim, normalize=normalizePCA
          , jacobian=JacPCA, print=2 )
vecPCAGN=resGN$par
F=matrix(vecPCAGN[1:(n*ndim)],n,ndim)
A=matrix(vecPCAGN[-(1:(n*ndim))],nvar,ndim)
resnorm=normalize_config(F,center=0, rotate=1, A=A)
FGN=resnorm$F; AGN=resnorm$A
rssGN2=ssq(Y-FGN%*%t(AGN))
gradGN2=drssPCA( vecPCAGN, ndim=ndim, Y=Y )
Print(rss0,rssGN2, min(abs(gradGN2)))
#
#
#
# final results
# rss
rss=cbind(rssNR,rssnlm, rssGN,rssGN2, rssnlsr)
cat("\n\n Comparison of final rss\n")
Print(rss)
grad=cbind( min(abs(gradnlm)),min(abs(gradNR)),min(abs(gradGN))
            , min(abs(gradGN2)), min(abs(gradnlsr)) )
colnames(grad)=c("NR",   "nlm",  "GN",   "GN2",  "nlsr")
cat("\n\n Comparison of final gradient\n")
Print(grad)
#
#
#####################################################################
#
#       End of PCA by Newton-Raphson and Gauss-Newton method
#
#####################################################################
#



## End(Not run)