```{r timeloop10, echo = FALSE, cache = TRUE} m <- 100 n <- 10 wv <- dv <- rep (1, n * (n - 1) / 2) wm <- dm <- 1 - diag (n) set.seed(12345) user1 <- rep (0, m) user2 <- rep (0, m) user3 <- rep (0, m) user4 <- rep (0, m) for (j in 1:m) { xold <- rnorm(2*n) t1 <- system.time(h1 <- smacofR(wm, dm, 2, xold = matrix(xold, n, 2), eps = 1e-15, itmax = 100), gcFirst = TRUE) t2 <- system.time(h2 <- smacofRC(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t3 <- system.time(h2 <- smacofRCU(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t4 <- system.time(h2 <- smacofRCU64(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) user1[j] <- t1[1] user2[j] <- t2[1] user3[j] <- t3[1] user4[j] <- t4[1] } ``` The four medians for $n=10$ are `r median(user1)` for `smacofR()`, `r median(user2)` for `smacofRC()`, `r median(user3)` for `smacofRCU()`, and `r median(user4)` for `smacofRCU64()`. This still gives roughly the same conclusions on relative times as for $n=4$.

```{r timeloop10_out, fig.align="center", echo = FALSE} boxplot(list(smacofR = user1, smacofRC = user2, smacofRCU = user3, smacofRCU64 = user4)) ```

```{r timeloop50, echo = FALSE, cache = TRUE} m <- 100 n <- 50 wv <- dv <- rep (1, n * (n - 1) / 2) wm <- dm <- 1 - diag (n) set.seed(12345) user1 <- rep (0, m) user2 <- rep (0, m) user3 <- rep (0, m) user4 <- rep (0, m) for (j in 1:m) { xold <- rnorm(2*n) t1 <- system.time(h1 <- smacofR(wm, dm, 2, xold = matrix(xold, n, 2), eps = 1e-15, itmax = 100), gcFirst = TRUE) t2 <- system.time(h2 <- smacofRC(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t3 <- system.time(h2 <- smacofRCU(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t4 <- system.time(h2 <- smacofRCU64(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) user1[j] <- t1[1] user2[j] <- t2[1] user3[j] <- t3[1] user4[j] <- t4[1] } ``` The medians for $n=50$ are `r median(user1)` for `smacofR()`, `r median(user2)` for `smacofRC()`, `r median(user3)` for `smacofRCU()`, and `r median(user4)` 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()`.

```{r timeloop50_out, fig.align="center", echo = FALSE} boxplot(list(smacofR = user1, smacofRC = user2, smacofRCU = user3, smacofRCU64 = user4)) ```

```{r timeloop100, echo = FALSE, cache = TRUE} m <- 100 n <- 100 wv <- dv <- rep (1, n * (n - 1) / 2) wm <- dm <- 1 - diag (n) set.seed(12345) user1 <- rep (0, m) user2 <- rep (0, m) user3 <- rep (0, m) user4 <- rep (0, m) for (j in 1:m) { xold <- rnorm(2*n) t1 <- system.time(h1 <- smacofR(wm, dm, 2, xold = matrix(xold, n, 2), eps = 1e-15, itmax = 100), gcFirst = TRUE) t2 <- system.time(h2 <- smacofRC(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t3 <- system.time(h2 <- smacofRCU(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t4 <- system.time(h2 <- smacofRCU64(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) user1[j] <- t1[1] user2[j] <- t2[1] user3[j] <- t3[1] user4[j] <- t4[1] } ``` For $n=100$ things begin to change. The medians are `r median(user1)` for `smacofR()`, `r median(user2)` for `smacofRC()`, `r median(user3)` for `smacofRCU()`, and `r median(user4)` 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()`.

```{r timeloop100_out, fig.align="center", echo = FALSE} boxplot(list(smacofR = user1, smacofRC = user2, smacofRCU = user3, smacofRCU64 = user4)) ```

```{r timeloop250, echo = FALSE, cache = TRUE} m <- 100 n <- 250 wv <- dv <- rep (1, n * (n - 1) / 2) wm <- dm <- 1 - diag (n) set.seed(12345) user1 <- rep (0, m) user2 <- rep (0, m) user3 <- rep (0, m) user4 <- rep (0, m) for (j in 1:m) { xold <- rnorm(2*n) t1 <- system.time(h1 <- smacofR(wm, dm, 2, xold = matrix(xold, n, 2), eps = 1e-15, itmax = 100), gcFirst = TRUE) t2 <- system.time(h2 <- smacofRC(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t3 <- system.time(h2 <- smacofRCU(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) t4 <- system.time(h2 <- smacofRCU64(wv, dv, 2, xold = xold, eps = 1e-15, itmax = 100), gcFirst = TRUE) user1[j] <- t1[1] user2[j] <- t2[1] user3[j] <- t3[1] user4[j] <- t4[1] } ``` The medians for $n=250$ are `r median(user1)` for `smacofR()`, `r median(user2)` for `smacofRC()`, `r median(user3)` for `smacofRCU()`, and `r median(user4)` 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()`.

```{r timeloop250_out, fig.align="center", echo = FALSE} boxplot(list(smacofR = user1, smacofRC = user2, smacofRCU = user3, smacofRCU64 = user4)) ```

#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 file_auxilary, code = readLines("smacofR.R"), size = 'footnotesize'} ``` ##R Glue ###smacofRC.R ```{r file_auxilary, code = readLines("smacofRC.R"), size = 'footnotesize'} ``` ###utilsRC.R ```{r file_auxilary, code = readLines("utilsRC.R")} ``` ###jacobiRC.R ```{r file_auxilary, code = readLines("jacobiRC.R")} ``` ###lapackeRC.R ```{r file_auxilary, code = readLines("lapackeRC.R")} ``` ##C code ###smacof.h ```{c file_auxilary, code = readLines("smacof.h"), eval = FALSE} ``` ###smacof.c ```{c file_auxilary, code = readLines("smacof.c"), eval = FALSE} ``` ###utils.c ```{c file_auxilary, code = readLines("utils.c"), eval = FALSE} ``` ###jacobi.c ```{c file_auxilary, code = readLines("jacobi.c"), eval = FALSE} ``` ###lapacke.c ```{c file_auxilary, code = readLines("lapacke.c"), eval = FALSE} ``` #References