lrank <- function (x, p = 2) { e <- eigen (x) v <- e\$vectors[, 1:p] return (tcrossprod (x %*% v, v)) } perturb <- function (b, e, p = 2) { n <- nrow (b) eb <- eigen (b) vb <- eb\$vectors ab <- eb\$values if (is.vector(e) && (length (e) == 2)) e <- outer (vb[, e[1]], vb[, e[2]]) + outer (vb[, e[2]], vb[, e[1]]) zeta <- crossprod(vb, e %*% vb) m1 <- m2 <- matrix (0, n, n) for (s in 1:p) { m2 <- m2 + zeta [s, s] * outer (vb[, s], vb[, s]) for (t in 1:n) { if (t == s) next m1 <- m1 + (ab[s] / (ab[t] - ab[s])) * zeta[s, t] * (outer (vb[, t], vb[, s]) + outer (vb[, s], vb[, t])) } } return (list ( ee = e, mm = m2 - m1 )) } power <- function (diss, c, w = 1 - diag (nrow (delta)), beta = 4 * nrow (delta), itmax=1000, eps = 1e-10, verbose = FALSE) { n <- nrow (diss) set.seed (12345) d <- squareDist (c) r <- -2 * w * (diss - d) diag(r) <- -rowSums (r) b <- c + r / beta eold <- doubleCenter (crossprod (matrix (rnorm (n * 100), 100, n)) / 100) eold <- eold / mnorm (eold) itel <- 1 sold <- -Inf repeat { dold <- squareDist (eold) f <- - 2 * w * dold diag(f) <- -rowSums (f) g <- eold - f / beta enew <- doubleCenter (perturb (b, g)\$mm) snew <- sum (eold * enew) enew <- enew / mnorm (enew) if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( sold, digits = 10, width = 15, format = "f" ), formatC ( snew, digits = 10, width = 15, format = "f" ), "\n" ) } if ((itel == itmax) || (abs (snew - sold) < eps)) break itel <- itel + 1 eold <- enew sold <- snew } return (list (e = enew, s = snew, beta = beta)) } basisMatrix <- function (diss, c, beta = 4 * nrow (diss), w = 1 - diag (nrow (diss))) { n <- nrow (diss) d <- squareDist (c) r <- -2 * w * (diss - d) diag(r) <- -rowSums (r) b <- c + r / beta m <- n * (n - 1) / 2 h <- matrix (0, m, m) dd <- matrix (0, n, n) dcbasis <- matrix (0, n^2, m) ij <- 1 for (i in 1:(n-1)) for (j in (i+1):n) { dcbasis [, ij] <- as.vector (aij (i, j, n)) ij <- ij + 1 } dcbasis <- svd (dcbasis)\$u ij <- 1 for (i in 2:n) for (j in 1:(i-1)) { e <- matrix (dcbasis [, ij], n, n) dd <- squareDist (e) f <- - 2 * w * dd diag(f) <- -rowSums (f) g <- e - f / beta p <- perturb (b, g)\$mm kl <- 1 for (k in 2:n) for (l in 1:(k-1)) { f <- matrix (dcbasis [, kl], n, n) h[ij, kl] <- sum (f * p) kl <- kl + 1 } ij <- ij + 1 } return (h) }