torgerson <- function(delta, p = 2) { doubleCenter <- function(x) { n <- dim(x)[1] m <- dim(x)[2] s <- sum(x) / (n * m) xr <- rowSums(x) / m xc <- colSums(x) / n return((x - outer(xr, xc, "+")) + s) } z <- eigen(-doubleCenter((as.matrix (delta) ^ 2) / 2), symmetric = TRUE) v <- pmax(z\$values, 0) return(z\$vectors[, 1:p] %*% diag(sqrt(v[1:p]))) } enorm <- function (x, w = 1) { return (sqrt (sum (w * (x ^ 2)))) } sqdist <- function (x) { s <- tcrossprod (x) v <- diag (s) return (outer (v, v, "+") - 2 * s) } mkBmat <- function (x) { d <- rowSums (x) x <- -x diag (x) <- d return (x) } mkPower <- function (x, r) { n <- nrow (x) return (abs ((x + diag (n)) ^ r) - diag(n)) } rStressMin <- function (delta, w = 1 - diag (nrow (delta)), p = 2, r = 0.5, eps = 1e-8, itmax = 100000, verbose = TRUE) { delta <- delta / enorm (delta, w) itel <- 1 xold <- torgerson (delta, p = p) xold <- xold / enorm (xold) n <- nrow (xold) nn <- diag (n) dold <- sqdist (xold) rold <- sum (w * delta * mkPower (dold, r)) nold <- sum (w * mkPower (dold, 2 * r)) aold <- rold / nold sold <- 1 - 2 * aold * rold + (aold ^ 2) * nold repeat { p1 <- mkPower (dold, r - 1) p2 <- mkPower (dold, (2 * r) - 1) by <- mkBmat (w * delta * p1) cy <- mkBmat (w * p2) ga <- 2 * sum (w * p2) be <- (2 * r - 1) * (2 ^ r) * sum (w * delta) de <- (4 * r - 1) * (4 ^ r) * sum (w) if (r >= 0.5) { my <- (aold * by) - (aold ^ 2) * (cy - de * nn) } if (r < 0.5) { my <- aold * (by - be * nn) - (aold ^ 2) * (cy - ga * nn) } xnew <- my %*% xold xnew <- xnew / enorm (xnew) dnew <- sqdist (xnew) rnew <- sum (w * delta * mkPower (dnew, r)) nnew <- sum (w * mkPower (dnew, 2 * r)) anew <- rnew / nnew snew <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( sold, digits = 10, width = 13, format = "f" ), formatC ( snew, digits = 10, width = 13, format = "f" ), "\n" ) } if ((itel == itmax) || ((sold - snew) < eps)) break () itel <- itel + 1 xold <- xnew dold <- dnew sold <- snew aold <- anew } return (list (x = xnew, alpha = anew, sigma = snew, itel = itel, r = r)) }