gfunc <- function (z, x, u, v, w) { return ((z + (w / outer (u, v)) * (x - z)) * sqrt(outer(u,v))) } sfunc <- function (z, u, v) { return (z / sqrt (outer (u, v))) } lfunc <- function (x, p) { l <- as.matrix(eigen (crossprod(x))\$vectors[,1:p]) yy <- x %*% tcrossprod (l) return (yy) } ffunc <- function (z, x, p, u, v, w) { z1 <- gfunc (z, x, u, v, w) z2 <- lfunc (z1, p) z3 <- sfunc (z2, u, v) return (as.vector(z3)) } cfunc <- function (z, x, p, u, v, w) { n <- length (u) m <- length (v) z <- matrix (z, n, m) return (ffunc (z, x, p, u, v, w)) } hfunc <- function (delta, u, v, w) { ((1 - w / outer (u, v)) * delta) * sqrt (outer (u, v)) } lrPerturb <- function(z, delta, x, p, u, v, w) { m <- ncol (x) gz <- gfunc (z, x, u, v, w) hd <- hfunc (delta, u, v, w) gh <- crossprod(gz, hd) + crossprod(hd, gz) hp <- matrix(0, m, m) e <- eigen(crossprod(gz)) lt <- e\$vectors vt <- e\$values lp <- as.matrix(lt[, 1:p]) for (s in 1:p) { vi <- ifelse (vt == vt[s], 0, 1 / (vt - vt[s])) cinv <- lt %*% diag (vi) %*% t(lt) hp <- hp + cinv %*% gh %*% outer (lp[, s], lp[, s]) } dz <- hd %*% tcrossprod(lp) - gz %*% (hp + t(hp)) return (sfunc (dz, u, v)) } lrDerive <- function(z, x, p, u, v, w) { n <- nrow(z) m <- ncol(z) k <- 1 dz <- matrix(0, n * m, n * m) for (j in 1:m) for (i in 1:n) { dx <- matrix(0, n, m) dx[i, j] <- 1 dz[, k] <- as.vector(lrPerturb(z, dx, x, p, u, v, w)) k <- k + 1 } return(dz) }