# Tweaking the SMACOF Engine Jan de Leeuw Version 06, October 10, 2017 --- Note: This is a working paper which will be expanded/updated frequently. All suggestions for improvement are welcome. The directory [gifi.stat.ucla.edu/tweaking](http://gifi.stat.ucla.edu/tweaking) has a pdf version, the bib file, the complete Rmd file with the code chunks, and the R and C source code. #Introduction The smacof algorithm is a majorization (or MM) algorithm to minimize the least squares loss function \label{E:loss} \sigma(X):=\frac12\mathop{\sum\sum}_{1<=i
Figure 1: Boxplot for n = 4

The four medians for $n=10$ are 0.008 for smacofR(), 0.104 for smacofRC(), 0.028 for smacofRCU(), and 0.03 for smacofRCU64(). This still gives roughly the same conclusions on relative times as for $n=4$.
Figure 2: Boxplot for n = 10

The medians for $n=50$ are 0.044 for smacofR(), 0.113 for smacofRC(), 0.033 for smacofRCU(), and 0.034 for smacofRCU64(). smacofR() is still 3 times faster than smacofRC(), but now smacofRCU() has caught up with smacofR(). Time for smacofRCU64() is the same as for smacofRCU().
Figure 3: Boxplot for n = 50

For $n=100$ things begin to change. The medians are 0.194 for smacofR(), 0.14 for smacofRC(), 0.055 for smacofRCU(), and 0.048 for smacofRCU64(). smacofR() is now the slowest of the three, and smacofRCU() remains three times faster than smacofRC(). smacofRCU64() is now about 15\% faster than smacofRCU().
Figure 4: Boxplot for n = 100

The medians for $n=250$ are 1.2945 for smacofR(), 0.325 for smacofRC(), 0.187 for smacofRCU(), and 0.145 for smacofRCU64(). smacofR() is losing ground, and smacofRC() and smacofRCU() are getting closer, although smacofRCU() still has the clear advantage. smacofRCU64() is about 20\% faster than smacofRCU().
Figure 5: Boxplot for n = 250

#Conclusion More extensive comparisons would be useful. We have not used weights, and consequently did not look at multidimensional scaling problems such as unfolding. We did not include non-metric options or individual difference options yet. This will come later. Also note that the code in the appendix includes both pure R and compact storage C versions to compute the Hessian of $\eqref{E:loss}$. Again, this is for future use, to implement variations of Newton's method and to draw pseudo-confidence intervals around the points of the configuration at the solution. There are also various utilities included, plus an interface to a subset of LAPACKE (@lapacke_16). Not all of this is used in the current project. For now it seems that for large $n$ the compact storage routines in C are comparable in speed to pure R smacof, if not faster, especially if the number of .C() calls within an iteration are limited. And they are obviously superior in terms of storage, although these days that is hardly a consideration any more. Ultimately, however, wasting space is just inelegant, and the interesting exercise was how much execution time we had to sacrifice to please our esthetic prejudices. #Appendix: Code ##Pure R code ###smacofR.R r library(MASS) smacofLossR <- function (d, w, delta) { return (sum (w * (delta - d) ^ 2) / 4) } smacofBmatR <- function (d, w, delta) { dd <- ifelse (d == 0, 0, 1 / d) b <- -dd * w * delta diag (b) <- -rowSums (b) return(b) } smacofVmatR <- function (w) { v <- -w diag(v) <- -rowSums(v) return (v) } smacofGuttmanR <- function (x, b, vinv) { return (vinv %*% b %*% x) } smacofGradientR <- function (x, b, v) { return ((v - b) %*% x) } smacofHmatR <- function (x, b, v, d, w, delta) { n <- nrow (x) p <- ncol (x) r <- n * p h <- matrix (0, r, r) dd <- ifelse (d == 0, 0, 1 / d) cc <- w * delta * (dd ^ 3) for (s in 1:p) { for (t in 1:s) { cst <- matrix (0, n, n) for (i in 1:n) { for (j in 1:n) { cst[i, j] <- cc[i, j] * (x[i, s] - x[j, s]) * (x[i, t] - x[j, t]) } } cst <- -cst diag(cst) <- -rowSums(cst) if (s == t) { h[(s - 1) * n + 1:n, (s - 1) * n + 1:n] <- b - cst } else { h[(s - 1) * n + 1:n, (t - 1) * n + 1:n] <- -cst h[(t - 1) * n + 1:n, (s - 1) * n + 1:n] <- -cst } } } return (h) } smacofHessianR <- function (x, b, v, d, w, delta) { n <- nrow (x) p <- ncol (x) h <- -smacofHmatR (x, b, v, d, w, delta) for (s in 1:p) { h[(s - 1) * n + 1:n, (s - 1) * n + 1:n] <- h[(s - 1) * n + 1:n, (s - 1) * n + 1:n] + v } return(h) } smacofInitialR <- function (delta, p) { n <- nrow(delta) delta <- delta ^ 2 rw <- rowSums (delta) / n sw <- sum (delta) / (n ^ 2) h <- -(delta - outer (rw, rw, "+") + sw) / 2 e <- eigen (h) ea <- e$values ev <- e$vector ea <- ifelse (ea > 0, sqrt (abs(ea)), 0)[1:p] return (ev[, 1:p] %*% diag (ea)) } smacofR <- function (w, delta, p, xold = smacofInitialR(delta, p), itmax = 100, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- as.matrix (dist (xold)) sold <- smacofLossR (dold, w, delta) bold <- smacofBmatR (dold, w, delta) vinv <- ginv (smacofVmatR (w)) repeat { xnew <- smacofGuttmanR (xold, bold, vinv) eiff <- max (abs (xold - xnew)) dnew <- as.matrix (dist (xnew)) bnew <- smacofBmatR (dnew, w, delta) snew <- smacofLossR (dnew, w, delta) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) }  ##R Glue ###smacofRC.R r smacofLossRC <- function (d, w, delta) { m <- length (d) h <- .C( "smacofLossC", as.double (d), as.double (w), as.double (delta), as.integer (m), stress = as.double (0) ) return (h$stress) } smacofDistRC <- function (x, p) { n <- round (length (x) / p) m <- n * (n - 1) / 2 h <- .C("smacofDistC", as.double (x), as.integer (n), as.integer(p), dist = as.double (rep(0, m))) return (h$dist) } smacofInitialRC <- function (delta, p) { m <- length (delta) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C( "smacofInitialC", as.double(delta), as.integer(n), as.integer(p), as.double (rep(0, n)), as.double (rep(0, r)), as.double (rep(0, n * n)), as.double (rep (0, n)), x = as.double (rep (0, n * p)) ) return (h$x) } smacofBmatRC <- function (d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C( "smacofBmatC", as.double(d), as.double(w), as.double(delta), as.integer (n), bmat = as.double(rep(0, r)) ) return (h$bmat) } smacofVmatRC <- function (w) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C("smacofVmatC", as.double(w), as.integer (n), vmat = as.double(rep(0, r))) return (h$vmat) } smacofGuttmanRC <- function (x, b, vinv) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) h <- .C( "smacofGuttmanC", as.double(x), as.double (b), as.double (vinv), as.integer(n), as.integer(p), as.double (rep(0, r)), y = as.double (rep(0, r)) ) return (h$y) } smacofGradientRC <- function (x, b, v) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) h <- .C( "smacofGradientC", as.double(x), as.double (b), as.double (v), as.integer(n), as.integer(p), y = as.double (rep(0, r)) ) return (h$y) } smacofHmatRC <- function (x, b, v, d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) q <- n * (n + 1) / 2 r <- length (x) p <- round (r / n) h <- .C( "smacofHmatC", as.double(x), as.double (b), as.double (v), as.double (w), as.double (delta), as.double (d), as.integer(n), as.integer(p), as.double (rep (0, q)), hmat = as.double (rep (0, r * (r + 1) / 2)) ) return (h$hmat) } smacofHessianRC <- function (x, b, v, d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) q <- n * (n + 1) / 2 r <- length (x) p <- round (r / n) h <- .C( "smacofHessianC", as.double(x), as.double (b), as.double (v), as.double (w), as.double (delta), as.double (d), as.integer(n), as.integer(p), as.double (rep (0, q)), hmat = as.double (rep (0, r * (r + 1) / 2)) ) return (h$hmat) } smacofUpdateRC <- function (x, w, delta, b, vinv) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) q <- n * (n - 1) / 2 h <- .C( "smacofUpdateC", as.double (x), as.double (w), as.double (delta), as.double (vinv), as.integer (n), as.integer (p), dnew = as.double (rep(0, q)), bnew = as.double (b), work = as.double (rep (0, n * p)), snew = as.double (0), xnew = as.double (rep (0, n * p)) ) return (h) } smacofUpdateRC64 <- function (x, w, delta, b, vinv) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) q <- n * (n - 1) / 2 h <- .C64( "smacofUpdateC", SIGNATURE = c(rep("double", 4), rep("integer", 2), rep("double", 5)), x = x, w = w, delta = delta, vinv = vinv, n = n, p = p, dnew = numeric_dc (q), bnew = as.double (b), work = numeric_dc (n * p), snew = numeric_dc (1), xnew = numeric_dc (n * p), INTENT = c(rep("r", 6), "w", "rw", rep ("w", 3)), NAOK = TRUE ) return (h) } smacofRC <- function (w, delta, p, xold = smacofInitialRC(delta, p), itmax = 100, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- smacofDistRC (xold, p) sold <- smacofLossRC (dold, w, delta) bold <- smacofBmatRC (dold, w, delta) vinv <- mpowerRC (smacofVmatRC (w), -1) repeat { xnew <- smacofGuttmanRC (xold, bold, vinv) eiff <- max (abs (xold - xnew)) dnew <- smacofDistRC (xnew, p) bnew <- smacofBmatRC (dnew, w, delta) snew <- smacofLossRC (dnew, w, delta) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) } smacofRCU <- function (w, delta, p, xold = smacofInitialRC(delta, p), itmax = 100, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- smacofDistRC (xold, p) sold <- smacofLossRC (dold, w, delta) bold <- smacofBmatRC (dold, w, delta) vinv <- mpowerRC (smacofVmatRC (w), -1) repeat { h <- smacofUpdateRC(xold, w, delta, bold, vinv) xnew <- h$xnew dnew <- h$dnew snew <- h$snew bnew <- h$bnew eiff <- max (abs (xold - xnew)) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) } smacofRCU64 <- function (w, delta, p, xold = smacofInitialRC(delta, p), itmax = 100, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- smacofDistRC (xold, p) sold <- smacofLossRC (dold, w, delta) bold <- smacofBmatRC (dold, w, delta) vinv <- mpowerRC (smacofVmatRC (w), -1) repeat { h <- smacofUpdateRC64(xold, w, delta, bold, vinv) xnew <- h$xnew dnew <- h$dnew snew <- h$snew bnew <- h$bnew eiff <- max (abs (xold - xnew)) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) }  ###utilsRC.R r dorpolRC <- function (n) { h <- .C("dorpol", as.integer(n), as.double (rep(0, n * n))) return (matrix(h[[2]], n, n)) } docentRC <- function (x) { m <- length (x) n <- round((sqrt (1 + 8 * m) - 1) / 2) h <- .C("docent", as.integer (n), as.double (x), as.double (rep(0, m))) return(h$y) } primatRC <- function (n, m, x, width = 6, precision = 4) { h <- .C( "primat", as.integer (n), as.integer (m), as.integer (width), as.integer (precision), as.double (x) ) } pritrlRC <- function (x, width = 6, precision = 4) { m <- length (x) n <- round ((sqrt (1 + 8 * m) + 1) / 2) h <- .C("pritrl", as.integer (n), as.integer (width), as.integer (precision), as.double (x)) } pritruRC <- function (x, width = 6, precision = 4) { m <- length (x) n <- round ((sqrt (1 + 8 * m) - 1) / 2) h <- .C("pritru", as.integer (n), as.integer (width), as.integer (precision), as.double (x)) } priarrRC <- function (n, m, r, x, width = 6, precision = 4) { h <- .C( "primat", as.integer (n), as.integer (m), as.integer (r), as.integer(width), as.integer (precision), as.double (x) ) } mpowerRC <- function (x, p) { m <- length (x) n <- round ((sqrt (1 + 8 * m) - 1) / 2) h <- .C("mpower", as.integer(n), as.double (x), as.double(p), xpow = as.double(rep(0, m))) return (h$xpow) } trimatRC <- function (x) { m <- length (x) n <- round ((sqrt (1 + 8 * m) - 1) / 2) h <- .C("trimat", as.integer(n), as.double (x), y = as.double (rep(0, n * n))) return(h$y) } mattriRC <- function (x) { n <- round (sqrt (length (x))) m <- n * (n + 1) / 2 h <- .C("mattri", as.integer(n), as.double (x), y = as.double (rep(0, m))) return(h$y) } mutrmaRC <- function (n, m, a, x) { h <- .C("mutrma", as.integer(n), as.integer(m), as.double (a), as.double (x), y = as.double (rep(0, n * m))) return (h$y) }  ###jacobiRC.R r jacobi <- function (a, itmax = 100, eps = 1e-10) { m <- length (a) n <- (sqrt (1 + 8 * m) - 1) / 2 i <- 1:n j <- (-(i ^ 2) / 2) + (n + (3 / 2)) * i - n h <- .C( "jacobiC", as.integer (n), a = as.double (a), e = as.double (rep (0, n * n)), as.double (rep(0, n)), as.double (rep(0, n)), as.integer (itmax), as.double (eps) ) return (list (values = h$a[j], vectors = matrix (h$e, n, n))) }  ###lapackeRC.R r dposvR <- function (a, b) { n <- nrow (a) m <- max (ncol (b), 1) h <- .C("dposv", as.integer(n), as.integer(m), as.double (a), as.double (b)) if (is.null(ncol(b))) { return (h[[4]]) } else { return (matrix(h[[4]], n, m)) } } dsyevdR <- function (a) { n <- nrow (a) x <- matrix (0, n, n) h <- .C("dsyevd", as.integer(n), as.double (a), as.double (x)) return (list(values = h[[3]][1:n], vectors = matrix(h[[2]], n, n))) } dgeqrR <- function (x) { n <- nrow (x) m <- ncol (x) h <- .C("dgeqrf", as.integer(n), as.integer(m), z = as.double (x), tau = as.double (rep (0, m))) return (list (z = matrix(h$z, n, m), tau = h$tau)) } dorgqrR <- function (z, tau) { n <- nrow (z) m <- ncol (z) h <- .C("dorgqr", as.integer(n), as.integer(m), as.double (z), as.double (tau)) return (matrix(h[[3]], n, m)) } dorthoR <- function(x) { n <- nrow (x) m <- ncol (x) h <- .C("dortho", as.integer(n), as.integer(m), as.double (x), as.double(rep(0, m))) return (matrix(h[[3]], n, m)) }  ##C code ###smacof.h c #ifndef SMACOF_H #define SMACOF_H #include #include #include #include #include static inline int VINDEX(const int); static inline int MINDEX(const int, const int, const int); static inline int SINDEX(const int, const int, const int); static inline int TINDEX(const int, const int, const int); static inline int AINDEX(const int, const int, const int, const int, const int); static inline double SQUARE(const double); static inline double THIRD(const double); static inline double FOURTH(const double); static inline double FIFTH(const double); static inline double MAX(const double, const double); static inline double MIN(const double, const double); static inline int IMIN(const int, const int); static inline int IMAX(const int, const int); static inline int IMOD(const int, const int); void dposv(const int *, const int *, double *, double *); void dsyevd(const int *, double *, double *); void dgeqrf(const int *, const int *, double *, double *); void dorgqr(const int *, const int *, double *, double *); void dortho(const int *, const int *, double *); void dorpol(const int *, double *); void primat(const int *, const int *, const int *, const int *, const double *); void priarr(const int *, const int *, const int *, const int *, const int *, const double *); void pritrl(const int *, const int *, const int *, const double *); void pritru(const int *, const int *, const int *, const double *); void mpinve(const int *, double *, double *); void mpower(const int *, double *, double *, double *); void mpinvt(const int *, double *, double *); void trimat(const int *, const double *, double *); void mattri(const int *, const double *, double *); void mutrma(const int *, const int *, const double *, const double *, double *); void jacobiC(const int *, double *, double *, double *, double *, const int *, const double *); void smacofDistC(const double *, const int *, const int *, double *); void smacofLossC(const double *, const double *, const double *, const int *, double *); void smacofBmatC(const double *, const double *, const double *, const int *, double *); void smacofVmatC(const double *, const int *, double *); void smacofGuttmanC(const double *, const double *, const double *, const int *, const int *, double *, double *); void smacofGradientC(const double *, const double *, const double *, const int *, const int *, double *, double *); void smacofHmatC(const double *, const double *, const double *, const double *, const double *, const double *, const int *, const int *, double *, double *); void smacofHessianC(const double *, const double *, const double *, const double *, const double *, const double *, const int *, const int *, double *, double *); void smacofInitialC(const double *, const int *, const int *, double *, double *, double *, double *, double *); static inline int VINDEX(const int i) { return i - 1; } static inline int MINDEX(const int i, const int j, const int n) { return (i - 1) + (j - 1) * n; } static inline int AINDEX(const int i, const int j, const int k, const int n, const int m) { return (i - 1) + (j - 1) * n + (k - 1) * n * m; } static inline int SINDEX(const int i, const int j, const int n) { return ((j - 1) * n) - (j * (j - 1) / 2) + (i - j) - 1; } static inline int TINDEX(const int i, const int j, const int n) { return ((j - 1) * n) - ((j - 1) * (j - 2) / 2) + (i - (j - 1)) - 1; } static inline double SQUARE(const double x) { return x * x; } static inline double THIRD(const double x) { return x * x * x; } static inline double FOURTH(const double x) { return x * x * x * x; } static inline double FIFTH(const double x) { return x * x * x * x * x; } static inline double MAX(const double x, const double y) { return (x > y) ? x : y; } static inline double MIN(const double x, const double y) { return (x < y) ? x : y; } static inline int IMAX(const int x, const int y) { return (x > y) ? x : y; } static inline int IMIN(const int x, const int y) { return (x < y) ? x : y; } static inline int IMOD(const int x, const int y) { return (((x % y) == 0) ? y : (x % y)); } #endif /* SMACOF_H */  ###smacof.c c #include "smacof.h" void smacofDistC(const double *x, const int *n, const int *p, double *d) { int k = 1, nn = *n, pp = *p; for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { double dij = 0.0; for (int s = 1; s <= pp; s++) { dij += SQUARE(x[MINDEX(i, s, nn)] - x[MINDEX(j, s, nn)]); } d[VINDEX(k)] = sqrt(dij); k++; } } } void smacofLossC(const double *dist, const double *w, const double *delta, const int *m, double *stress) { int mm = *m; *stress = 0.0; for (int k = 1; k <= mm; k++) { *stress += w[VINDEX(k)] * SQUARE(delta[VINDEX(k)] - dist[VINDEX(k)]); } *stress /= 2.0; return; } void smacofBmatC(const double *dist, const double *w, const double *delta, const int *n, double *b) { int nn = *n; for (int i = 1; i <= nn; i++) { b[TINDEX(i, i, nn)] = 0.0; } for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { int k = SINDEX(i, j, nn); double dinv = (dist[k] == 0.0) ? 0.0 : 1.0 / dist[k]; double elem = w[k] * delta[k] * dinv; b[TINDEX(i, j, nn)] = -elem; b[TINDEX(i, i, nn)] += elem; b[TINDEX(j, j, nn)] += elem; } } return; } void smacofVmatC(const double *w, const int *n, double *v) { int nn = *n; for (int i = 1; i <= nn; i++) { v[TINDEX(i, i, nn)] = 0.0; } for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { double elem = w[SINDEX(i, j, nn)]; v[TINDEX(i, j, nn)] = -elem; v[TINDEX(i, i, nn)] += elem; v[TINDEX(j, j, nn)] += elem; } } return; } void smacofHmatC(const double *x, const double *bmat, const double *vmat, const double *w, const double *delta, const double *dist, const int *n, const int *p, double *work, double *h) { int nn = *n, pp = *p, np = nn * pp; for (int s = 1; s <= pp; s++) { for (int t = 1; t <= s; t++) { for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { work[TINDEX(i, j, nn)] = 0.0; } } for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { double f1 = (x[MINDEX(i, s, nn)] - x[MINDEX(j, s, nn)]); double f2 = (x[MINDEX(i, t, nn)] - x[MINDEX(j, t, nn)]); double f3 = THIRD(dist[SINDEX(i, j, nn)]); double f4 = delta[SINDEX(i, j, nn)] * w[SINDEX(i, j, nn)]; f3 = (f3 < 1e-10) ? 0 : 1 / f3; double elem = f1 * f2 * f4 * f3; work[TINDEX(i, j, nn)] = -elem; work[TINDEX(i, i, nn)] += elem; work[TINDEX(j, j, nn)] += elem; } } if (s == t) { for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { h[TINDEX((s - 1) * nn + i, (s - 1) * nn + j, nn * pp)] = bmat[TINDEX(i, j, nn)] - work[TINDEX(i, j, nn)]; } } } if (s != t) { for (int i = 1; i <= nn; i++) { for (int j = 1; j <= nn; j++) { int ij = IMIN(i, j); int ji = IMAX(i, j); int sj = (s - 1) * nn + j; int si = (t - 1) * nn + i; h[TINDEX(sj, si, np)] = -work[TINDEX(ji, ij, nn)]; } } } } } return; } void smacofGuttmanC(const double *x, const double *bmat, const double *vinv, const int *n, const int *p, double *work, double *y) { (void)mutrma(n, p, bmat, x, work); (void)mutrma(n, p, vinv, work, y); return; } void smacofGradientC(const double *x, const double *bmat, const double *vmat, const int *n, const int *p, double *work, double *y) { int nn = *n, pp = *p; (void)mutrma(n, p, vmat, x, y); (void)mutrma(n, p, bmat, x, work); for (int i = 1; i <= nn; i++) { for (int s = 1; s <= pp; s++) { y[MINDEX(i, s, nn)] -= work[MINDEX(i, s, nn)]; } } return; } void smacofHessianC(const double *x, const double *bmat, const double *vmat, const double *w, const double *delta, const double *dist, const int *n, const int *p, double *work, double *h) { int nn = *n, pp = *p, np = nn * pp; (void)smacofHmatC(x, bmat, vmat, w, delta, dist, n, p, work, h); for (int j = 1; j <= np; j++) { for (int i = j; i <= np; i++) { h[TINDEX(i, j, np)] = -h[TINDEX(i, j, np)]; } } for (int s = 1; s <= pp; s++) { for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { h[TINDEX((s - 1) * nn + i, (s - 1) * nn + j, np)] += vmat[TINDEX(i, j, nn)]; } } } return; } void smacofNewtonC() { return; } void smacofInitialC(const double *delta, const int *n, const int *p, double *work1, double *work2, double *work3, double *work4, double *x) { int nn = *n, pp = *p, itmax = 100; double s, ss = 0.0, eps = 1e-6; for (int i = 1; i <= nn; i++) { work1[VINDEX(i)] = 0.0; for (int j = 1; j <= nn; j++) { if (i == j) continue; int ij = IMAX(i, j); int ji = IMIN(i, j); s = SQUARE(delta[SINDEX(ij, ji, nn)]); ss += s; work1[VINDEX(i)] += s; if (j < i) continue; work2[TINDEX(ij, ji, nn)] = s; } work1[VINDEX(i)] /= (double)nn; } ss /= SQUARE((double)nn); for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { work2[TINDEX(i, j, nn)] -= work1[VINDEX(i)]; work2[TINDEX(i, j, nn)] -= work1[VINDEX(j)]; work2[TINDEX(i, j, nn)] += ss; work2[TINDEX(i, j, nn)] *= -0.5; } } (void)jacobiC(n, work2, work3, work1, work4, &itmax, &eps); for (int i = 1; i <= nn; i++) { for (int j = 1; j <= pp; j++) { s = work2[TINDEX(j, j, nn)]; if (s <= 0) continue; x[MINDEX(i, j, nn)] = work3[MINDEX(i, j, nn)] * sqrt(s); } } return; } void smacofUpdateC(const double *xold, const double *w, const double *delta, const double *vinv, const int *n, const int *p, double *dnew, double *bmat, double *work, double *snew, double *xnew) { int nn = *n, pp = *p, mm = nn * (nn - 1) / 2; (void)mutrma(n, p, bmat, xold, work); (void)mutrma(n, p, vinv, work, xnew); for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { double dij = 0.0; for (int s = 1; s <= pp; s++) { dij += SQUARE(xnew[MINDEX(i, s, nn)] - xnew[MINDEX(j, s, nn)]); } dnew[SINDEX(i, j, nn)] = sqrt(dij); } } for (int i = 1; i <= nn; i++) { bmat[TINDEX(i, i, nn)] = 0.0; } for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { int k = SINDEX(i, j, nn); double dinv = (dnew[k] == 0.0) ? 0.0 : 1.0 / dnew[k]; double elem = w[k] * delta[k] * dinv; bmat[TINDEX(i, j, nn)] = -elem; bmat[TINDEX(i, i, nn)] += elem; bmat[TINDEX(j, j, nn)] += elem; } } *snew = 0.0; for (int k = 1; k <= mm; k++) { *snew += w[VINDEX(k)] * SQUARE(delta[VINDEX(k)] - dnew[VINDEX(k)]); } *snew /= 2.0; }  ###utils.c c #include "smacof.h" // complete set of orthogonal polynomials void dorpol(const int *n, double *p) { for (int i = 1; i <= *n; i++) { p[MINDEX(i, 1, *n)] = 1.0; double di = (double)i; for (int j = 2; j <= *n; j++) { p[MINDEX(i, j, *n)] = p[MINDEX(i, j - 1, *n)] * di; } } (void)dortho(n, n, p); return; } // double centering void docent(const int *n, double *x, double *y) { int nn = *n; double *ssum = (double *)calloc((size_t)nn, sizeof(double)), t = 0.0; for (int i = 1; i <= nn; i++) { double s = 0.0; for (int j = 1; j <= nn; j++) { int ij = IMAX(i, j); int ji = IMIN(j, i); s += x[TINDEX(ij, ji, nn)]; t += x[TINDEX(ij, ji, nn)]; } ssum[VINDEX(i)] = s / ((double)nn); } t /= (double)SQUARE(nn); for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { y[TINDEX(i, j, nn)] = -(x[TINDEX(i, j, nn)] - ssum[VINDEX(i)] - ssum[VINDEX(j)] + t) / 2.0; } } free(ssum); return; } // print a general matrix void primat(const int *n, const int *m, const int *w, const int *p, const double *x) { for (int i = 1; i <= *n; i++) { for (int j = 1; j <= *m; j++) { printf(" %+*.*f ", *w, *p, x[MINDEX(i, j, *n)]); } printf("\n"); } printf("\n\n"); return; } // print strict lower triangle void pritrl(const int *n, const int *w, const int *p, const double *x) { for (int i = 1; i <= *n; i++) { for (int j = 1; j <= i; j++) { if (i == j) { for (int k = 1; k <= *w + 2; k++) { printf("%c", '*'); } } if (i > j) { printf(" %+*.*f ", *w, *p, x[SINDEX(i, j, *n)]); } } printf("\n"); } printf("\n\n"); return; } // print inclusive lower triangle void pritru(const int *n, const int *w, const int *p, const double *x) { for (int i = 1; i <= *n; i++) { for (int j = 1; j <= i; j++) { printf(" %+*.*f ", *w, *p, x[TINDEX(i, j, *n)]); } printf("\n"); } printf("\n\n"); return; } // print general array void priarr(const int *n, const int *m, const int *r, const int *w, const int *p, const double *x) { for (int k = 1; k <= *r; k++) { for (int i = 1; i <= *n; i++) { for (int j = 1; j <= *m; j++) { printf(" %+*.*f ", *w, *p, x[AINDEX(i, j, k, *n, *m)]); } printf("\n"); } printf("\n\n"); } printf("\n\n\n"); return; } // arbitrary power of a matrix void mpower(const int *n, double *x, double *power, double *xpow) { int nn = *n, itmax = 100; double eps = 1e-6; double *e = (double *)calloc((size_t)SQUARE(nn), sizeof(double)); double *oldi = (double *)calloc((size_t)nn, sizeof(double)); double *oldj = (double *)calloc((size_t)nn, sizeof(double)); (void)jacobiC(n, x, e, oldi, oldj, &itmax, &eps); int k = 1; for (int i = 1; i <= nn; i++) { double s = x[VINDEX(k)]; oldi[VINDEX(i)] = (s > 1e-10) ? pow(s, *power) : 0.0; k += (nn - (i - 1)); } for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { double s = 0.0; for (int k = 1; k <= nn; k++) { s += e[MINDEX(i, k, nn)] * e[MINDEX(j, k, nn)] * oldi[VINDEX(k)]; } xpow[TINDEX(i, j, nn)] = s; } } free(e); free(oldi); free(oldj); return; } // inclusive lower triangle to symmetric void trimat(const int *n, const double *x, double *y) { int nn = *n; for (int i = 1; i <= nn; i++) { for (int j = 1; j <= nn; j++) { y[MINDEX(i, j, nn)] = (i >= j) ? x[TINDEX(i, j, nn)] : x[TINDEX(j, i, nn)]; } } return; } // symmetric to inclusive lower triangular void mattri(const int *n, const double *x, double *y) { int nn = *n; for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { y[TINDEX(i, j, nn)] = x[MINDEX(i, j, nn)]; } } return; } // premultiply matrix by symmetric matrix in inclusive triangular storage void mutrma(const int *n, const int *m, const double *a, const double *x, double *y) { int nn = *n, mm = *m; for (int i = 1; i <= nn; i++) { for (int j = 1; j <= mm; j++) { double s = 0.0; for (int k = 1; k <= nn; k++) { int ik = IMAX(i, k); int ki = IMIN(i, k); s += a[TINDEX(ik, ki, nn)] * x[MINDEX(k, j, nn)]; } y[MINDEX(i, j, nn)] = s; } } } // direct sum of two matrices -- can be used recursively void dirsum(const int *n, const int *m, const double *a, const double *b, double *c) { int nn = *n, mm = *m, mn = nn + mm; for (int j = 1; j <= nn; j++) { for (int i = j; i <= nn; i++) { c[TINDEX(i, j, mn)] = a[TINDEX(i, j, nn)]; } } for (int j = 1; j <= mm; j++) { for (int i = j; i <= mm; i++) { c[TINDEX(nn + i, nn + j, mn)] = b[TINDEX(i, j, mm)]; } } return; }  ###jacobi.c c #include "smacof.h" void jacobiC(const int *nn, double *a, double *evec, double *oldi, double *oldj, const int *itmax, const double *eps) { int n = *nn, itel = 1; double d = 0.0, s = 0.0, t = 0.0, u = 0.0, v = 0.0, p = 0.0, q = 0.0, r = 0.0; double fold = 0.0, fnew = 0.0; for (int i = 1; i <= n; i++) { for (int j = 1; j <= n; j++) { evec[MINDEX(i, j, n)] = (i == j) ? 1.0 : 0.0; } } for (int i = 1; i <= n; i++) { fold += SQUARE(a[TINDEX(i, i, n)]); } while (true) { for (int j = 1; j <= n - 1; j++) { for (int i = j + 1; i <= n; i++) { p = a[TINDEX(i, j, n)]; q = a[TINDEX(i, i, n)]; r = a[TINDEX(j, j, n)]; if (fabs(p) < 1e-10) continue; d = (q - r) / 2.0; s = (p < 0) ? -1.0 : 1.0; t = -d / sqrt(SQUARE(d) + SQUARE(p)); u = sqrt((1 + t) / 2); v = s * sqrt((1 - t) / 2); for (int k = 1; k <= n; k++) { int ik = IMIN(i, k); int ki = IMAX(i, k); int jk = IMIN(j, k); int kj = IMAX(j, k); oldi[VINDEX(k)] = a[TINDEX(ki, ik, n)]; oldj[VINDEX(k)] = a[TINDEX(kj, jk, n)]; } for (int k = 1; k <= n; k++) { int ik = IMIN(i, k); int ki = IMAX(i, k); int jk = IMIN(j, k); int kj = IMAX(j, k); a[TINDEX(ki, ik, n)] = u * oldi[VINDEX(k)] - v * oldj[VINDEX(k)]; a[TINDEX(kj, jk, n)] = v * oldi[VINDEX(k)] + u * oldj[VINDEX(k)]; } for (int k = 1; k <= n; k++) { oldi[VINDEX(k)] = evec[MINDEX(k, i, n)]; oldj[VINDEX(k)] = evec[MINDEX(k, j, n)]; evec[MINDEX(k, i, n)] = u * oldi[VINDEX(k)] - v * oldj[VINDEX(k)]; evec[MINDEX(k, j, n)] = v * oldi[VINDEX(k)] + u * oldj[VINDEX(k)]; } a[TINDEX(i, i, n)] = SQUARE(u) * q + SQUARE(v) * r - 2 * u * v * p; a[TINDEX(j, j, n)] = SQUARE(v) * q + SQUARE(u) * r + 2 * u * v * p; a[TINDEX(i, j, n)] = u * v * (q - r) + (SQUARE(u) - SQUARE(v)) * p; } } fnew = 0.0; for (int i = 1; i <= n; i++) { fnew += SQUARE(a[TINDEX(i, i, n)]); } if (((fnew - fold) < *eps) || (itel == *itmax)) break; fold = fnew; itel++; } return; }  ###lapacke.c c #include "smacof.h" void dposv(const int *n, const int *m, double *a, double *b) { lapack_int nn = (lapack_int)*n, mm = (lapack_int)*m; (void)LAPACKE_dposv(LAPACK_COL_MAJOR, 'U', nn, mm, a, nn, b, nn); return; } void dsyevd(const int *n, double *a, double *x) { lapack_int nn = (lapack_int)*n; (void)LAPACKE_dsyevd(LAPACK_COL_MAJOR, 'V', 'U', nn, a, nn, x); return; } void dgeqrf(const int *n, const int *m, double *a, double *tau) { lapack_int nn = (lapack_int)*n, mm = (lapack_int)*m; (void)LAPACKE_dgeqrf(LAPACK_COL_MAJOR, nn, mm, a, nn, tau); return; } void dorgqr(const int *n, const int *m, double *a, double *tau) { lapack_int nn = (lapack_int)*n, mm = (lapack_int)*m; (void)LAPACKE_dorgqr(LAPACK_COL_MAJOR, nn, mm, mm, a, nn, tau); return; } void dortho(const int *n, const int *m, double *a) { lapack_int nn = (lapack_int)*n, mm = (lapack_int)*m; double *tau = calloc((size_t)mm, sizeof(double)); (void)LAPACKE_dgeqrf(LAPACK_COL_MAJOR, nn, mm, a, nn, tau); (void)LAPACKE_dorgqr(LAPACK_COL_MAJOR, nn, mm, mm, a, nn, tau); free(tau); return; }  #References