map_b <- function (delta, x, w, beta) { s <- tcrossprod (x) d <- squareDist (s) r <- -2 * w * (delta - d) diag (r) <- -rowSums (r) return (s + r / beta) } map_gamma <- function (b, p) { e <- slanczos (b, p) ea <- e\$values ev <- e\$vectors return (ev[, 1:p] %*% diag (sqrt (pmax (0, ea[1:p])))) } map_phi <- function (delta, alpha, w = 1 - diag (nrow (delta)), p = 2, beta = 4 * nrow (delta), basis = matrixBasis (nrow(delta), p)) { n <- nrow (delta) x <- matrix (basis %*% alpha, n, p) b <- map_b (delta, x, w, beta) x <- map_gamma (b, p) return (drop (as.vector(x) %*% basis)) } map_db <- function (x, y, w, beta) { n <- nrow (x) da <- tcrossprod (x, y) + tcrossprod (y, x) dd <- matrix (0, n, n) for (i in 1:n) for (j in 1:n) dd <- dd + w[i, j] * sum ((x[i,] - x[j,]) * (y[i,] - y[j,])) * aij (i, j, n) return (da - 2 * dd / beta) } map_dgamma <- function (b, db, p) { n <- nrow (b) e <- eigen (b) l <- e\$values k <- e\$vectors xi <- crossprod (k, db %*% k) dg <- matrix (0, n, p) for (s in 1:p) { dgs <- (1 / (2 * sqrt (l[s]))) * xi[s, s] * k[, s] for (t in 1:n) { if (t == s) next dgs <- dgs - (sqrt (l[s]) / (l[t] - l[s])) * xi[s, t] * k[, t] } dg[, s] <- dgs } return (dg) } map_dphi <- function (delta, alpha, epsilon, w = 1 - diag (nrow (delta)), p = 2, beta = 4 * nrow (delta), basis = matrixBasis (nrow(delta), p)) { n <- nrow (delta) x <- matrix (basis %*% alpha, n, p) y <- matrix (basis %*% epsilon, n, p) b <- map_b (delta, x, w, beta) db <- map_db (x, y, w, beta) dg <- map_dgamma (b, db, p) return (drop (as.vector(dg) %*% basis)) } janJacobian <- function (delta, alpha, w = 1 - diag (nrow (delta)), p = 2, beta = 4 * nrow (delta), basis = matrixBasis (nrow(delta), p)) { n <- nrow (delta) m <- p * (n - 1) j <- matrix (0, m, m) for (i in 1:m) { j[i, ] <- map_dphi (delta, alpha, ei (i, m), w, p, beta, basis) } return (j) } power <- function (delta, x, w = 1 - diag (nrow (x)), beta = 4 * nrow (x), itmax = 1000, eps = 1e-10, verbose = FALSE, basis = matrixBasis (nrow(x), ncol(x))) { n <- nrow (x) p <- ncol (x) set.seed (12345) y <- drop (as.vector (x) %*% basis) eold <- ei (1, ncol (basis)) lold <- -Inf itel <- 1 repeat { enew <- map_dphi (delta, y, eold, w, p, beta, basis) lnew <- sum (eold * enew) enew <- enew / mnorm (enew) edif <- 1 - abs (sum (eold * enew)) if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( lold, digits = 10, width = 15, format = "f" ), formatC ( lnew, digits = 10, width = 15, format = "f" ), formatC ( edif, digits = 10, width = 15, format = "f" ), "\n" ) } if ((itel == itmax) || (edif < eps)) break eold <- enew lold <- lnew itel <- itel + 1 } return (list (e = enew, l = lnew, itel = itel)) }