--- title: "Multidimensional Scaling with Distance Bounds" author: "Jan de Leeuw" date: "Version 12, January 16, 2017" output: pdf_document: keep_tex: yes number_sections: yes toc: yes toc_depth: 3 html_document: keep_md: yes number_sections: yes toc: yes toc_depth: 3 fontsize: 12pt graphics: yes bibliography: updown.bib abstract: We give an algorithm, with R code, to minimize the multidimensional scaling stress loss function under the condition that some or all of the fitted distances are between given positive upper and lower bounds. This paper combines theory, algorithms, code, and results of @deleeuw_E_17b and @deleeuw_E_17c. --- {r function_code, echo = FALSE} source("updown.R")  {r packages, echo = FALSE} options (digits = 10) library (captioner)  {r equation_numbering, echo = FALSE} figure_nums <- captioner (prefix = "Figure") lemma_nums <- captioner (prefix = "Lemma") theorem_nums <- captioner (prefix = "Theorem") lemma_nums (name= "feasibility", caption ="[Feasibility]", display = FALSE) lemma_nums (name= "monotonicity", caption ="[Monotonicity]", display = FALSE) theorem_nums (name= "convergence", caption ="[Convergence]", display = FALSE) theorem_nums (name= "necessary", caption ="[Necessary]", display = FALSE) figure_nums (name = "equi_regular_conf", caption = "Equidistance Data, Unconstrained Left, Lower Bounds Right", display = FALSE) figure_nums (name = "equi_regular_both", caption = "Equidistance Data, Unconstrained Left, Lower Bounds Right", display = FALSE) figure_nums (name = "equi_regular_equal_1", caption = "Equidistance Data, Distances Constrained to One", display = FALSE) figure_nums (name = "equi_regular_equal_2", caption = "Equidistance Data, Distances Constrained to One", display = FALSE) figure_nums (name = "gruijter_unconstrained", caption = "De Gruijter Data, Unconstrained Solution", display = FALSE) figure_nums (name = "gruijter_from_below", caption = "De Gruijter Data, MDS from Below", display = FALSE) figure_nums (name = "gruijter_from_above", caption = "De Gruijter Data, MDS from Above", display = FALSE) figure_nums (name = "gruijter_interval", caption = "De Gruijter Data, Distance Intervals", display = FALSE)  Note: This is a working paper which will be expanded/updated frequently. All suggestions for improvement are welcome. The directory [gifi.stat.ucla.edu/updown](http://gifi.stat.ucla.edu/updown) has a pdf version, the complete Rmd file with all code chunks, the bib file, and the R source code. #Introduction In this paper we study the minimization of *stress*, defined as $$\mathbf{stress}(X):=\frac14\sum_{i=1}^n\sum_{j=1}^n w_{ij}(\delta_{ij}-d_{ij}(X))^2$$ over $n\times p$ *configuration matrices* $X$. Here $W=\{w_{ij}\}$ and $\Delta=\{\delta_{ij}\}$ are given symmetric, hollow, and non-negative matrices of, respectively, *weights* and *dissimilarities*. The $d_{ij}(X)$ are Euclidean distances between the points with coordinates in the rows of $X$. Thus $$d_{ij}^2(X):=(x_i-x_j)'(x_i-x_j)=(e_i-e_j)'XX'(e_i-e_j)=\mathbf{tr}\ X'A_{ij}X.$$ Here $e_i$ and $e_j$ are unit vectors, with all elements equal to zero except for the single element $i$ or $j$, which is equal to one. Also $A_{ij}:=(e_i-e_j)(e_i-e_j)'$. The problem of minimizing **stress** is well-studied (see, for example, @borg_groenen_05). In this paper we analyze the problem of minimizing **stress** over all configurations $X$ for which $\alpha_{ij}\leq d_{ij}(X)\leq\beta_{ij}$ for some or all pairs $i,j$, where the $\alpha_{ij}$ and $\beta_{ij}$ are given bounds satisfying $0\leq\alpha_{ij}<\beta_{ij}\leq+\infty$. We suppose the constraints are strongly consistent, i.e. there is at least one $n\times p$ configuration matrix $X$ such that \alpha_{ij} {r equi_confplot, fig.align = "center", echo = FALSE} par(mfrow=c(1,2), pty = "s") plot(h2ax, type = "n", col = "RED", xlab = "dim 1", ylab = "dim 2") text (h2a$x, as.character(1:nrow (h2a$x)), col = "RED") plot(h2b$x, type = "n", col = "RED", xlab = "dim 1", ylab = "dim 2") text (h2b$x, as.character(1:nrow (h2b$x)), col = "RED") for (i in c(2,3,6,8,9)) lines (matrix(c(h2b$x[1, ], h2b$x[i, ]), 2, 2, byrow = TRUE))  r figure_nums("equi_regular_conf") In the second analysis, now using eight objects with equal dissimilarties, we require$d_{ij}(X)\geq 1$for all$i,j$with$|i-j|=1$, and$d_{ij}(X)\leq 2$for all$i,j$with$|i-j|=2$. {r equi_equal_data, echo = FALSE, cache = TRUE} x <- matrix(c(0,0,0,1,0,2,1,2,2,2,2,1,2,0,1,0),8,2, byrow =TRUE) delta <- rep (1, 28) bndup <- matrix(Inf, 8, 8) bndlw <- matrix(0, 8, 8) bndlw[abs(outer(1:8,1:8,"-"))==1] <- 1.0 bndup[abs(outer(1:8,1:8,"-"))==2] <- 2.0 bndlw <- triangleFromSymmetric(bndlw, diagonal = FALSE) bndup <- triangleFromSymmetric(bndup, diagonal = FALSE) h2e <- smacofUpDown(delta, xini = x, verbose = FALSE) h2c <- smacofUpDown(delta, xini = 2.1 * h2e$x, bndup = bndup, bndlw = bndlw, verbose = FALSE)  In the unconstrained case we find **stress** equal to r h2e$stress after r h2e$itel iterations, in the constrained case we find r h2c$stress after r h2c$itel iterations. The lower bound constraints are all active, which means the distances between successive points are all one (indicated with lines in figure r figure_nums("equi_regular_both", display = "n"). The upper bound constraints are all inactive.
{r equi_equal_plot, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2), pty="s") plot(h2e$x, type = "n", col = "RED", xlab = "dim 1", ylab = "dim 2") text (h2e$x, as.character(1:nrow (h2e$x)), col = "RED") plot(h2c$x, type = "n", col = "RED", xlab = "dim 1", ylab = "dim 2") text (h2c$x, as.character(1:nrow (h2c$x)), col = "RED") for (i in 1:7) lines (matrix(c(h2c$x[i, ], h2c$x[i + 1, ]), 2, 2, byrow = TRUE)) 
r figure_nums("equi_regular_both")

So far we have assumed $\alpha_{ij}<\beta_{ij}$ for all $i,j$. It is possible, however, to let $\alpha_{ij}=\beta_{ij}$, which means we require $d_{ij}(X)=\alpha_{ij}=\beta_{ij}$. Note the complete set of constraints still needs to be consistent, and the initial configuration must still be feasible. In our example we require all seven distances between successive points to be equal to one. {r equi_equal_d, echo = FALSE, cache = TRUE} bndup <- matrix(Inf, 8, 8) bndlw <- matrix(0, 8, 8) bndlw[abs(outer(1:8,1:8,"-"))==1] <- 1.0 bndup[abs(outer(1:8,1:8,"-"))==1] <- 1.0 bndlw <- triangleFromSymmetric(bndlw, diagonal = FALSE) bndup <- triangleFromSymmetric(bndup, diagonal = FALSE) h2d <- smacofUpDown(delta, xini = x, bndup = bndup, bndlw = bndlw, itmax = 2000, verbose = FALSE)  Our algorithm in this case becomes extremely slow. It stops at the upper bound of r h2d$itel iterations, with **stress** r h2d$stress. Of course all constraints are now active. The Lagrange multipliers are large, serving as penalty parameters to force the equality constraints. {r equi_equal_d_out, echo = FALSE} mprint(h2dmultipliers)  {r equi_equal_d_plot, fig.align = "center", echo = FALSE} par(pty ="s") plot(h2dx, xlab = "dim 1", ylab = "dim 2", type = "n") text(h2d$x, as.character(1:8), col = "RED") for (i in 1:7) lines (matrix(c(h2d$x[i, ], h2d$x[i + 1, ]), 2, 2, byrow = TRUE))  r figure_nums("equi_regular_equal_1") Note, however, that the previous analysis where we merly required$d_{ij}(X)\geq 1$for$|i-j|=1$produced a solution feasible for the current problem with **stress** equal to r h2c$stress. This makes it interesting see what goes on if we increase the upper bound of the number of iterations even more, say, to 10000. {r equi_equal_e, echo = FALSE, cache = TRUE} h2f <- smacofUpDown(delta, xini = x, bndup = bndup, bndlw = bndlw, itmax = 10000, verbose = FALSE)  We now find **stress** r h2f$stress after r h2f$itel iterations. The Lagrange multipliers for this solution are much smaller (first seven correspond with upper bounds, second seven with lower bounds). What basically happens is that for each $(i,j)$ with $|i-j|=1$ either the upper or the lower constraint is seen to be violated and gets a zero Lagrange multiplier. The corresponding Lagrange multiplier for the corresponding other constraint is positive, because it is interpreted as being satisfied. This is a consequence of using floating point with limited precision, and the difference between $\alpha_{ij}=0.99\cdots$ and $\beta_{ij}=1.00\cdots$. {r equi_equal_e_out_1, echo = FALSE} mprint(h2f$multipliers[1:7], d = 10) mprint(h2f$multipliers[8:14], d = 10)  The successive distances are {r equi_equal_e_out_2, echo = FALSE} mprint(symmetricFromTriangle(h2fdist, diagonal = FALSE)[outer(1:8,1:8,"-")==1], d = 10)  {r equi_equal_e_plot, fig.align = "center", echo = FALSE} par(pty ="s") plot(h2fx, xlab = "dim 1", ylab = "dim 2", type = "n") text(h2f$x, as.character(1:8), col = "RED") for (i in 1:7) lines (matrix(c(h2f$x[i, ], h2f$x[i + 1, ]), 2, 2, byrow = TRUE))  r figure_nums("equi_regular_equal_2") ##Dutch Political Parties 1967 As the next illustration we use data from @degruijter_67, with average dissimilarity judgments between Dutch political parties in 1967. The data are {r poldist-data, echo = FALSE} poldist <- structure( c( 5.63, 5.27, 4.6, 4.8, 7.54, 6.73, 7.18, 6.17, 6.72, 5.64, 6.22, 5.12, 4.59, 7.22, 5.47, 5.46, 4.97, 8.13, 7.55, 6.9, 4.67, 3.2, 7.84, 6.73, 7.28, 6.13, 7.8, 7.08, 6.96, 6.04, 4.08, 6.34, 7.42, 6.88, 6.36, 7.36 ), Labels = c("KVP", "PvdA", "VVD", "ARP", "CHU", "CPN", "PSP", "BP", "D66"), Size = 9L, call = quote(as.dist.default(m = polpar)), class = "dist", Diag = FALSE, Upper = FALSE ) print (poldist) poldata <- as.vector (poldist) polnames <- attr (poldist, "Labels")  {r politics_comp, echo = FALSE, cache = TRUE} xini <- torgerson (symmetricFromTriangle(poldata, diagonal = FALSE), p = 2) h1a <- smacofUpDown (poldata, xini = xini) h1b <- smacofUpDown (poldata, xini = .1 * xini, bndup = poldata) h1c <- smacofUpDown (poldata, xini = 10 * xini, bndlw = poldata)  First, three different but comparable anlyses are done. The first does not impose restriction, the second is *MDS from below*, which requires$d_{ij}(X)\leq\delta_{ij}$for all$i,j$. And the third is *MDS from above*, for which$d_{ij}(X)\geq\delta_{ij}$for all$i,j$. The configurations are quite similar, except for the position of D'66, which at the time was a novelty in Dutch politics. The value of **stress** at the solutions is, respectively, r h1a$stress, r h1b$stress, and r h1c$stress.
{r politics_plot_1, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2), pty="s") plot(h1a$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h1a$x, polnames, col = "RED") plot(h1a$delta, h1a$dist, col = "RED", xlab = "dissimilarities", ylab = "distances") abline(0, 1) 
r figure_nums("gruijter_unconstrained")

{r politics_plot_2, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2), pty="s") plot(h1b$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h1b$x, polnames, col = "RED") plot(h1b$delta, h1b$dist, col = "RED", xlab = "dissimilarities", ylab = "distances") abline(0, 1) 
r figure_nums("gruijter_from_below")

{r politics_plot_3, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2), pty="s") plot(h1c$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h1c$x, polnames, col = "RED") plot(h1c$delta, h1c$dist, col = "RED", xlab = "dissimilarities", ylab = "distances") abline(0, 1) 
r figure_nums("gruijter_from_above")

{r gruijter_interval, echo = FALSE, cache = TRUE} xpol <- structure(c(-4.84538498357043e-16, -2.54322587848901, -3.89644810362578, -3.42647895687901, -1.35322222513677, 1.35322222513677, 3.426478956879, 3.89644810362578, 2.54322587848901, -3.95655709625336, -3.03089857746785, -0.687048929599558, 1.97827854812668, 3.71794750706741, 3.71794750706741, 1.97827854812668, -0.687048929599557, -3.03089857746785), .Dim = c(9L, 2L)) bndlw <- rep (2, 36) bndup <- rep (8, 36) h1d <- smacofUpDown(poldata, xini = xpol, bndlw = bndlw, bndup = bndup, verbose = FALSE)  The same data are used in the next analysis, which requires $2\leq d_{ij}(X)\leq 8$ for all $i,j$. We converge to **stress** r h1d$stress in r h1d$itel iterations. There are eight active constrains, with four distances equal to the lower bound and four equal to the upper bound. The configuration in figure r figure_nums("gruijter_interval", display = "n") shows the distances equal to the lower bound in blue and those equal to the upper bound in red. The active constraints are also clearly visible in the Shepard plot in figure r figure_nums("gruijter_interval", display = "n").
{r gruijter_interval_plot, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2), pty="s") plot(h1d$x, type = "n", xlab = "dim 1", ylab = "dim 2") text (h1d$x, polnames) lines (matrix(c(h1d$x[1, ], h1d$x[3, ]), 2, 2, byrow = TRUE), col = "BLUE") lines (matrix(c(h1d$x[3, ], h1d$x[4, ]), 2, 2, byrow = TRUE), col = "BLUE") lines (matrix(c(h1d$x[4, ], h1d$x[5, ]), 2, 2, byrow = TRUE), col = "BLUE") lines (matrix(c(h1d$x[6, ], h1d$x[7, ]), 2, 2, byrow = TRUE), col = "BLUE") lines (matrix(c(h1d$x[1, ], h1d$x[6, ]), 2, 2, byrow = TRUE), col = "RED") lines (matrix(c(h1d$x[3, ], h1d$x[6, ]), 2, 2, byrow = TRUE), col = "RED") lines (matrix(c(h1d$x[4, ], h1d$x[8, ]), 2, 2, byrow = TRUE), col = "RED") lines (matrix(c(h1d$x[5, ], h1d$x[8, ]), 2, 2, byrow = TRUE), col = "RED") plot(h1d$delta, h1d$dist, col = "RED", xlab = "dissimilarities", ylab = "distances") abline (h = 2) abline (h = 8) abline (0, 1) 
r figure_nums("gruijter_interval")

#Appendix: Code ##updown.R {r file_auxilary, code = readLines("updown.R")}  ##auxilary.R {r file_auxilary, code = readLines("auxilary.R")}  ##mdsUtils.R {r file_auxilary, code = readLines("mdsUtils.R")}  ##smacofUpDown.R {r file_auxilary, code = readLines("smacofUpDown.R")}  #References