Nothing
#' Smooth NPD to Nearest PSD or PD Matrix
#'
#' Smoothing an indefinite matrix to a PSD matrix via theory described by Lurie
#' and Goldberg
#'
#'
#' @param R Indefinite Matrix.
#' @param start.val Optional vector of start values for Cholesky factor of S.
#' @param Wghts An optional matrix of weights such that the objective function
#' minimizes wij(rij - sij)^2, where wij is Wghts[i,j].
#' @param PD Logical (default = FALSE). If PD = TRUE then the objective
#' function will smooth the least squares solution to insure Positive
#' Definitness.
#' @param Penalty A scalar weight to scale the Lagrangian multiplier. Default =
#' 50000.
#' @param eps A small value to add to zero eigenvalues if smoothed matrix must
#' be PD. Default = 1e-07.
#' @return \item{RLG}{Lurie Goldberg smoothed matrix.} \item{RKB}{Knol and
#' Berger smoothed matrix.} \item{convergence}{0 = converged solution, 1 =
#' convergence failure.} \item{start.val}{Vector of start.values.}
#' \item{gr}{Analytic gradient at solution.} \item{Penalty}{Scalar used to
#' scale the Lagrange multiplier.} \item{PD}{User-supplied value of PD.}
#' \item{Wghts}{Weights used to scale the squared euclidean distances.}
#' \item{eps}{Value added to zero eigenvalue to produce PD matrix.}
#' @author Niels Waller
#' @keywords statistics
#' @export
#' @examples
#'
#' data(BadRLG)
#'
#' out<-smoothLG(R = BadRLG, Penalty = 50000)
#' cat("\nGradient at solution:", out$gr,"\n")
#' cat("\nNearest Correlation Matrix\n")
#' print( round(out$RLG,8) )
#'
#' ################################
#' ## Rousseeuw Molenbergh example
#' data(BadRRM)
#'
#' out <- smoothLG(R = BadRRM, PD=TRUE)
#' cat("\nGradient at solution:", out$gr,"\n")
#' cat("\nNearest Correlation Matrix\n")
#' print( round(out$RLG,8) )
#'
#' ## Weights for the weighted solution
#' W <- matrix(c(1, 1, .5,
#' 1, 1, 1,
#' .5, 1, 1), nrow = 3, ncol = 3)
#' tmp <- smoothLG(R = BadRRM, PD = TRUE, eps=.001)
#' cat("\nGradient at solution:", out$gr,"\n")
#' cat("\nNearest Correlation Matrix\n")
#' print( round(out$RLG,8) )
#' print( eigen(out$RLG)$val )
#'
#' ## Rousseeuw Molenbergh
#' ## non symmetric matrix
#' T <- matrix(c(.8, -.9, -.9,
#' -1.2, 1.1, .3,
#' -.8, .4, .9), nrow = 3, ncol = 3,byrow=TRUE)
#' out <- smoothLG(R = T, PD = FALSE, eps=.001)
#'
#' cat("\nGradient at solution:", out$gr,"\n")
#' cat("\nNearest Correlation Matrix\n")
#' print( round(out$RLG,8) )
#'
smoothLG <- function(R, start.val = NULL, Wghts = NULL, PD = FALSE,
Penalty = 50000, eps=1e-07){
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
## smoothLG: a function for smoothing an indefinite matrix
## to a PSD matrix via theory described by Lurie and Goldberg
## Niels Waller
## February 14, 2016
## February 29, 2016 added var names
## June 30, 2016 added error checks
##
## Arguments:
## R Indefinite Matrix
## start.val optional start values for Cholesky factor of S
## Wghts an optional matrix of weights such that
## fnc minimizes wij(rij - sij)^2
## where wij is Wght[i,j]
## PD logical. if PD = TRUE the fnc will smooth
## the least squares solution to insure PD
## Penality Scalar weight for the Lagrangian multiplier.
## Defaut = 50000
## eps value to add to zero eigenvalues if smoothed matrix must be PD
##
## Value:
## RLG Lurie Goldberg smooth matrix
## RKB Knol and Berger smoothed matrix
## convergence 0 = converged, 1 = failure
## start.val start.values
## gr analytic gradient at solution
## Penalty
## PD user-supplied value of PD
## Wghts
## eps
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
## add some error checkig code to check that input
## matrices are symmetric etc.
# if(!isSymmetric(R)) stop("\nInput matrix not symmetric\n")
if( (PD!=TRUE) & (PD!=FALSE) ) stop("\nPD should be a logical (TRUE or FALSE)\n")
## Declarations and definitions
LG = 1 #start val for Lagrange multiplier
Nvar <- ncol(R)
lowtriN<-Nvar*(Nvar+1)/2
I <- diag(Nvar)
## if weights not supplied all weights are 1
if(is.null(Wghts)) W <-matrix(1,Nvar,Nvar)
else W <- Wghts
## Z is a matrix of Zeros
Z <-matrix(0,Nvar,Nvar)
## define trace function
TR <- function(x) sum(diag(x))
## define custom diag function
DIAG <- function(x) diag(diag(x))
## corSmooth: a function to smooth an ID matrix to a PD matrix
## by a function described by Knol and Berger
## This function is used to generate start values
corSmooth <-function (R, eps = .01) #1e+08 * .Machine$double.eps)
{
KDK <- eigen(R)
Dplus <- D <- KDK$values
Dplus[D <= 0] <- eps
Dplus <- diag(Dplus)
K <- KDK$vectors
Rtmp <- K %*% Dplus %*% t(K)
invSqrtDiagRtmp <- diag(1/sqrt(diag(Rtmp)))
invSqrtDiagRtmp %*% Rtmp %*% invSqrtDiagRtmp
}
## smooth R by Knol and Berger method
## Rkb = R Knol and Berger
RKB<-corSmooth(R)
## Compute Cholesky of RS for start values
if( is.null(start.val) ){
start.val <-c((t(chol(RKB)))[lower.tri(RKB, diag=T)][-1], LG)
}
## position of the Lagrange multiplier in the
## list of parameter estimates
LagrangeNum <- length(start.val)
## makeL create a lower-triangular matrix of the same
## dimension as R. This will be updated during optimization.
makeL <- function(x){
L <- I
lower <- lower.tri(L, diag = TRUE)
L[lower] <- c(1,x[-LagrangeNum])
L
} #End makeL
## fnc to minimize (1/2) squared euclidean distance
## between S and R.
fnc <- function(x){
L <- makeL(x)
S <- L%*% t(L)
.5*TR( (sqrt(W) * (R - S)) %*% ( sqrt(W) *(R - S)) ) + x[LagrangeNum]* Penalty* TR( DIAG(S-I) %*% DIAG(S-I) )
} #End fnc
###############################################
## Compute gradient of fnc
grvec<-rep(0,lowtriN-1)
grFNC <- function(x){
L <- makeL(x)
Lt <- t(L)
S <- L %*% t(L)
## function to compute gradient element Lij
gradLij <- function(i,j){
Jij<-Jji<-Z
Jij[i,j]<-1
Jji[j,i]<-1
## partial deriv of main function w.r.t. Lij
dfdLij <- TR( (W*(R - S)) %*% -(L%*%Jji + Jij%*%Lt))
## partial deriv of constraint function w.r.t. Lij
dgdLij <- 2 * x[lowtriN]* Penalty * TR( (DIAG(S - I) %*% (L%*%Jji + Jij%*%Lt)) )
dfdLij + dgdLij
} #End gradLij
## fill in gradient vector
k<-0
for(j in 1:Nvar){
for(i in 2:Nvar){
if(i>=j) {
k = k+1
grvec[k]<-gradLij(i,j)
}
}
}
grvec[k+1] <- Penalty * TR( DIAG(S-I)%*%DIAG(S-I) )
grvec
}
##minimize fnc.
out <- optim(par = start.val, fn = fnc,
gr=grFNC,
method="BFGS",
control=list(fnscale=1,
maxit=10000,
reltol=1e-20,
abstol=1e-20))
## check if function converged
if(normF(grFNC(out$par)) > 1e-05) out$convergence <-1
# assembled smoothed matrix
L<-makeL(out$par[-LagrangeNum])
Rsmooth <- L %*% t(L)
## Is smoothed matrix a proper correlation matrix
if( abs(sum(abs(diag(Rsmooth))) - Nvar) > 1e-05 ) out$convergence <-1
if(det(Rsmooth) > 1) out$convergence <-1
if(out$convergence==0){
## Smooth to insure PD matrix?
if(PD) Rsmooth<-corSmooth(Rsmooth, eps)
colnames(Rsmooth) <- rownames(Rsmooth) <- colnames(R)
colnames(RKB) <- rownames(RKB) <- colnames(R)
Rsmooth <- .5 * (Rsmooth + t(Rsmooth))
RKB <- .5 * (RKB + t(RKB))
list(RLG = Rsmooth,
RKB = RKB,
convergence = out$convergence,
par = out$par,
start.val = start.val,
gr=grFNC(out$par),
Penalty = Penalty,
PD = PD,
Wghts = W,
eps = eps
)
}
## if function failed to converge
else list(RLG = 999,
RKB = 999,
convergence = 1,
# par = rep(99, length(start.val)),
start.val = start.val,
gr=NULL,
Penalty = Penalty,
PD = PD,
Wghts=W,
eps=eps
)
} ## End smoothLG
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.