# An Improved Exact Algorithm for Least-Squares Unidimensional Scaling

Authors: Gintaras Palubeckis
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Givenn objects and an n×n symmetric dissimilarity matrix D with zero main diagonal
and nonnegative off-diagonal entries, the least-squares unidimensional scaling problem
asks to find an arrangement of objects along a straight line such that the pairwise distances
between them reflect dissimilarities represented by the matrix D. In this paper, we propose an improved branch-and-bound algorithm for solving this problem. The main
ingredients of the algorithm include an innovative upper bounding technique relying on
the linear assignment model and a dominance test which allows considerably reducing the
redundancy in the enumeration process. An initial lower bound for the algorithm is provided
by an iterated tabu search heuristic. To enhance the performance of this heuristic
we develop an efficient method for exploring the pairwise interchange neighborhood of a
solution in the search space. The basic principle and formulas of the method are also used
in the implementation of the dominance test. We report computational results for both
randomly generated and real-life based problem instances. In particular, we were able to
solve to guaranteed optimality the problem defined by a 36×36 Morse code dissimilarity matrix.


---

## Body

## 1. Introduction

Least-squares (orL2) unidimensional scaling is an important optimization problem in the field of combinatorial data analysis. It can be stated as follows. Suppose that there are n objects and, in addition, there are pairwise dissimilarity data available for them. These data form an n×n symmetric dissimilarity matrix, D=(dij), with zero main diagonal and nonnegative off-diagonal entries. The problem is to arrange the objects along a straight line such that the pairwise distances between them reflect dissimilarities given by the matrix D. Mathematically, the problem can be expressed as
(1)minΦ(x)=∑i=1n-1‍∑j=i+1n‍(dij-|xi-xj|)2,
where x denotes the vector (x1,…,xn) and xi, i∈{1,…,n}, is the coordinate for the object i on the line. We note that two (or even more) objects may have the same coordinate. In the literature, the double sum given by (1) is called a stress function [1].Unidimensional scaling and its generalization, multidimensional scaling [1, 2], are dimensionality reduction techniques for exploratory data analysis and data mining [3]. The applicability of these techniques has been reported in a diversity of fields. In recent years, unidimensional scaling has found new applications in several domains including genetics, psychology, and medicine. Caraux and Pinloche [4] developed an environment to graphically explore gene expression data. This environment offers unidimensional scaling methods, along with some other popular techniques. Ng et al. [5] considered the problem of constructing linkage disequilibrium maps for human genome. They proposed a method for this problem which is centered on the use of a specific type of the unidimensional scaling model. Armstrong et al. [6] presented a study in which they evaluated interest profile differences across gender, racial-ethnic group, and employment status. The authors used two measures of profile similarity, the so-called Q correlation and Euclidean distances. In the former case, unidimensional scaling appeared to be a good model for interpreting the obtained results. Hubert and Steinley [7] applied unidimensional scaling to analyze proximity data on the United States Supreme Court Justices. They have found that agreement among the justices was better represented by this approach than using a hierarchical classification method. Dahabiah et al. [8, 9] proposed a data mining technique based on a unidimensional scaling algorithm. They applied this technique on a digestive endoscope database consisting of a large number of pathologies. Ge et al. [10] developed a method for unsupervised classification of remote sensing images. Unidimensional scaling plays a pivotal role in this classifier. Bihn et al. [11] examined the diversity of ant assemblages in the Atlantic Forest of Southern Brazil. To visualize nestedness of ant genera distribution, a unidimensional scaling approach was adopted.From (1), we see that in order to solve the unidimensional scaling problem (UDSP for short) we need to find not only a permutation of objects, but also their coordinates preserving the order of the objects as given by the permutation. Thus, at first glance, it might seem that the UDSP is a much more difficult problem than other linear arrangement problems. Fortunately, this is not the case. As shown by Defays [12], the UDSP can be reduced to the following combinatorial optimization problem:
(2)maxp∈ΠF(p)=∑k=1n‍(∑l=1k-1‍dp(k)p(l)-∑l=k+1n‍dp(k)p(l))2,
where Π is the set of all permutations of W={1,…,n}, that is, the set of all vectors (p(1),p(2),…,p(n)) such that p(k)∈W, k=1,…,n, and p(k)≠p(l), k=1,…,n-1, l=k+1,…,n. Letting Si, i∈W, denote the sum ∑j=1,j≠indij, we can write
(3)F(p)=∑k=1n‍(2∑l=1k-1‍dp(k)p(l)-Sp(k))2=∑k=1n‍(2∑l=k+1n‍dp(k)p(l)-Sp(k))2.
Thus, the UDSP is essentially a problem over all permutations of n objects. With an optimal solution p* to (2) at hand, we can compute the optimal coordinates for the objects by setting xp*(1)=0 and using the following relationship [13] for k+1=2,…,n:
(4)xp*(k+1)=xp*(k)+(∑l=k+1ndp*(k)p*(l)-∑l=1k-1dp*(k)p*(l)=xp*(k)+-∑l=k+2ndp*(k+1)p*(l)+∑l=1kdp*(k+1)p*(l))n-1.There is an important body of the literature on algorithms for solving the UDSP. Since the feasible solutions to the problem are linear arrangements of objects, a natural approach to the UDSP is to apply a general dynamic programming paradigm. An algorithm based on this paradigm was proposed by Hubert and Arabie [14]. However, as pointed out in the papers by Hubert et al. [15] and Brusco and Stahl [13], the dynamic programming algorithm requires very large memory for storing the partial solutions and, in fact, runs out of memory long before running out of time. As reported in these papers, the largest UDSP instances solved using dynamic programming were of sizes 25 and 26. An alternative approach for solving the UDSP is to use the branch-and-bound method. The first algorithm following this approach is that of Defays [12]. In his algorithm, during branching, the objects are assigned in an alternating fashion to either the leftmost free position or the rightmost free position in the n-permutation under construction. For example, at level 4 of the search tree, the four selected objects are p(1), p(n), p(2), and p(n-1), and there are no object assigned to positions 3,…,n-2 in the permutation p. To fathom partial solutions, the algorithm uses a symmetry test, an adjacency test, and a bound test. Some more comments on the Defays approach were provided in the paper by Brusco and Stahl [13]. The authors of that paper also proposed two improved branch-and-bound algorithms for unidimensional scaling. The first of them adopts the above-mentioned branching strategy borrowed from the Defays method. Meanwhile, in the second algorithm of Brusco and Stahl, an object chosen for branching is always assigned to the leftmost free position. Brusco and Stahl [13, 16] presented improved bounding procedures and an interchange test expanding the adjacency test used by Defays [12]. These components of their branch-and-bound algorithms are briefly sketched in Section 2 and Section 4, respectively. Brusco and Stahl [13] have reported computational experience for both randomly generated and empirical dissimilarity matrices. The results suggest that their branch-and-bound implementations are superior to the Defays algorithm. Moreover, the results demonstrate that these implementations are able to solve UDSP instances that are beyond the capabilities of the dynamic programming method.There are also a number of heuristic algorithms for unidimensional scaling. These include the smoothing technique [17], the iterative quadratic assignment heuristic [15], and simulated annealing adaptations [18, 19]. A heuristic implementation of dynamic programming for escaping local optima has been presented by Brusco et al. [20]. Computational results for problem instances with up to 200 objects were reported in [18–20].The focus of this paper is on developing a branch-and-bound algorithm for the UDSP. The primary intention is to solve exactly some larger problems than those solved by previous methods. The significant features of our algorithm are the use of an innovative upper bounding procedure relying on the linear assignment model and the enhancement of the interchange test from [13]. In the current paper, this enhancement is called a dominance test. A good initial solution and thus lower bound are obtained by running an iterated tabu search algorithm at the root node of the search tree. To make this algorithm more efficient, we provide a novel mechanism and underlying formulas for searching the pairwise interchange neighborhood of a solution in O(n2) time, which is an improvement over a standard O(n3)-time approach from the literature. The same formulas are also used in the dominance test. We present the results of empirical evaluation of our branch-and-bound algorithm and compare it against the better of the two branch-and-bound algorithms of Brusco and Stahl [13]. The largest UDSP instance that we were able to solve with our algorithm was the problem defined by a 36×36 Morse code dissimilarity matrix.The remainder of this paper is organized as follows. In Section2, we present an upper bound on solutions of the UDSP. In Section 3, we outline an effective tabu search procedure for unidimensional scaling and propose a method for exploring the pairwise interchange neighborhood. Section 4 gives a description of the dominance test for pruning nodes in the search. This test is used in our branch-and-bound algorithm presented in Section 5. The results of computational experiments are reported in Section 6. Finally, Section 7 concludes the paper.


## 2. Upper Bounds

In this section, we first briefly recall the bounds used by Brusco and Stahl [13] in their second branch-and-bound algorithm. Then we present a new upper bounding mechanism for the UDSP. We discuss the technical details of the procedure for computing our bound. We choose a variation of this procedure which yields slightly weakened bounds but is quite fast.Like in the second algorithm of Brusco and Stahl, in our implementation of the branch-and-bound method partial solutions are constructed by placing each new object in the leftmost free position. For a partial solutionp=(p(1),…,p(r)), r∈{1,…,n-1}, we denote by W(p) the set of the objects in p. Thus W(p)={p(1),…,p(r)}, and W-(p)=W∖W(p) is the set of unselected objects. At the root node of the search tree, p is empty, and W-(p) is equal to W. The first upper bound used by Brusco and Stahl, when adopted to p, takes the following form:
(5)UBS≔Fr(p)+∑i∈W-(p)‍Si2,
where Fr(p) denotes the sum of the terms in (3) corresponding to assigned objects. Thus, Fr(p)=∑k=1r‍(Sp(k)-2∑l=1k-1‍dp(k)p(l))2.To describe the second bound, we defineD^=(d^ij) to be an n×(n-1) matrix obtained from D by sorting the nondiagonal entries of each row of D in nondecreasing order. Let ηik=(Si-2∑j=1k-1‍d^ij)2 and μk=max1⩽i⩽nηik for k=1,…,⌈n/2⌉. For other values of k, μk is obtained from the equality μk=μn-k+1, k∈{1,…,n}. The second upper bound proposed in [13] is given by the following expression:
(6)Fr(p)+∑k=r+1n‍μk.The previous outlined bounds are computationally attractive but not strong enough when used to prune the search space. In order to get a tighter upper bound, we resort to the linear assignment model(7)max∑i=1n-r‍∑j=1n-r‍eijyijs.t.∑j=1n-r‍yij=1,i=1,…,n-r,∑i=1n-r‍yij=1,j=1,…,n-r,yij∈{0,1},i,j=1,…,n-r,
where eij, i,j=1,…,n-r are entries of an (n-r)×(n-r) profit matrix E with rows and columns indexed by the objects in the set W-(p). Next, we will construct the matrix E. For a permutation p∈Π, we denote by fp(k) (resp., fp(k)′) the sum of dissimilarities between object p(k) and all objects preceding (resp., following) p(k) in p. Thus, fp(k)=∑l=1k-1‍dp(k)p(l), fp(k)′=∑l=k+1n‍dp(k)p(l). From (3), we have
(8)F(p)=∑k=1n‍(2max(fp(k),fp(k)′)-Sp(k))2,
and 2max(fp(k),fp(k)′)⩾Sp(k)=fp(k)+fp(k)′ for each k=1,…,n. Suppose now that p=(p(1),…,p(r)) is a partial solution and Π(p) is the set of all n-permutations that can be obtained by extending p. Assume, for the sake of simplicity, that W-(p)={1,…,n-r}. In order to construct a suitable matrix E, we have to bound fp(k) and fp(k)′ in (8) from above. To this end, for any i∈W-(p), let (di1′,…,di,n-r-1′) denote the list obtained by sorting dij, j∈W-(p)∖{i}, in nonincreasing order. To get the required bounds, we make use of the sums γij=∑l=1j‍di1′, i∈W-(p), j=0,1,…,n-r-1, and ρik=∑l=1r‍dip(l)+γi,k-r-1, i∈W-(p), k=r+1,…,n. Suppose that p^ is any permutation in Π(p) and p^(k)=i for a position k∈{r+1,…,n}. Then it is easy to see that fp^(k)=fi⩽ρik and fp^(k)′=fi′⩽γi,n-k. Indeed, in the case of fp^(k), for example, we have fp^(k)=∑l=1r‍dip(l)+∑l=r+1k-1‍dip^(l), and the last term in this expression is less than or equal to γi,k-r-1. We define E to be the matrix with entries eij=(2max(ρi,r+j,γi,n-r-j)-Si)2, i,j=1,…,n-r. Let U* denote the optimal value of the linear assignment problem (LAP) (7) with the profit matrix E. From the definition of E, we immediately obtain an upper bound on the optimal value of the objective function F.Proposition 1.
For a partial solutionp=(p(1),…,p(r)), U:=Fr(p)+U*⩾maxp^∈Π(p)F(p^).From the definition ofρ and γ values, it can be observed that, with the increase of k, ρik increases, while γi,n-k decreases. Therefore, when constructing the ith row of the matrix E, it is expedient to first find the largest k∈{r+1,…,n} for which ρik⩽γi,n-k. We denote such k by k~ and assume that k~=r if ρi,r+1>γi,n-r-1. Then ei,k-r=(2γi,n-k-Si)2 for k=r+1,…,k~, and ei,k-r=(2ρik-Si)2 for k=k~+1,…,n.The most time-consuming step in the computation of the upper boundU is solving the linear assignment problem. In order to provide a computationally less expensive bounding method we replace U* with the value of a feasible solution to the dual of the LAP (7). To get the bound, we take αi=max1⩽k⩽n-reik=ei,n-r=Si2, i=1,…,n-r, βj=max1⩽i⩽n-r(eij-αi), j=1,…,n-r. Clearly, αi+βj⩾eij for each pair i,j=1,…,n-r. Thus the vector (α1,…,αn-r,β1,…,βn-r) satisfies the constraints of the dual of the LAP instance of interest. Let U~=∑i=1n-r‍αi+∑j=1n-r‍βj. Then, from the duality theory of linear programming, we can establish the following upper bound for F(p^).Proposition 2.
For a partial solutionp=(p(1),…,p(r)), U′:=Fr(p)+U~⩾maxp^∈Π(p)F(p^).Proof.
The inequality follows from Proposition1 and the fact that U~⩾U*.One point that needs to be emphasized is thatU′ strengthens the first bound proposed by Brusco and Stahl.Proposition 3.
U
′
⩽
U
B
S.Proof.
Indeed,Fr(p)+∑i∈W-(p)‍Si2=Fr(p)+∑i=1n-r‍αi⩾Fr(p)+U~. The equality in the previous chain of expressions comes from the choice of αi and our assumption on W-(p), whereas the inequality is due to the fact that βj⩽0 for each j=1,…,n-r.As an example, we consider the following4×4 dissimilarity matrix borrowed from Defays [12]:
(9)[-4234-7227-4324-].
The same matrix was also used by Brusco and Stahl [13] to illustrate their bounding procedures. From Table 2 in their paper, we find that the solution p=(2,4,1,3) is optimal for this instance of the UDSP. The optimal value is equal to 388. We compute the proposed upper bound for the whole Defays problem. In this case, W-(p)={1,…,4} and Fr(p)=0. The matrix E is constructed by noting that k~=2 for each i=1,…,4. It takes the following form:
(10)[812525811698181169169818116981252581].
Using the expressions for αi and βj, we find that α1=α4=81, α2=α3=169, β1=β4=0, and β2=β3=-56. Thus, U~=∑i=14‍αi+∑j=14‍βj=500-112=388=maxp∈ΠF(p). Hence, our bound for the Defays instance is tight. For comparison, the bound (5) is equal to ∑i=14‍αi=500. To obtain the second bound of Brusco and Stahl, we first compute the μ values. Since μ1=μ4=169 and μ2=μ3=81, the bound (6) (taking the form ∑i=14‍μi) is also equal to 500.


## 3. Tabu Search

As it is well known that the effectiveness of the bound test in branch-and-bound algorithms depends on the value of the best solution found throughout the search process. A good strategy is to create a high-quality solution before the search is started. For this task, various heuristics can be applied. Our choice was to use the tabu search (TS) technique. It should be noted that it is not the purpose of this paper to investigate heuristic algorithms for solving the UDSP. Therefore, we restrict ourselves merely to presenting an outline of our implementation of tabu search.Forp∈Π, let pkm denote the permutation obtained from p by interchanging the objects in positions k and m. The set N2(p)={pkm∣k,m=1,…,n,k<m} is the 2-interchange neighborhood of p. At each iteration, the tabu search algorithm performs an exploration of the neighborhood N2(p). A subset of solutions in N2(p) are evaluated by computing Δ(p,k,m)=F(pkm)-F(p). In order to get better quality solutions we apply TS iteratively. This technique, called iterated tabu search (ITS), consists of the following steps: (1)
initialize the search with some permutation of the objects;(2)
apply a tabu search procedure. Check if the termination condition is satisfied. If so, then stop. Otherwise proceed to(3);(3)
apply a solution perturbation procedure. Return to(2).In our implementation, the search is started from a randomly generated permutation. Also, some data are prepared to be used in the computation of the differencesΔ(p,k,m). We will defer this issue until later in this section. In step (2) of ITS, at each iteration of the TS procedure, the neighborhood N2(p) of the current solution p is explored. The procedure selects a permutation pkm∈N2(p), which either improves the best available objective function value so far or, if this does not happen, has the largest value of Δ(p,i,j) among all permutations pij∈N2(p) for which the objects in positions i and j are not in the tabu list. The permutation pkm is then accepted as a new current solution. The number of iterations in the tabu search phase is bounded above by a parameter of the algorithm. The goal of the solution perturbation phase is to generate a permutation which is employed as the starting point for the next invocation of the tabu search procedure. This is achieved by taking the current permutation and performing a sequence of pairwise interchanges of objects. Throughout this process, each object is allowed to change its position in the permutation at most once. At each step, first a set of feasible pairs of permutation positions with the largest Δ values is formed. Then a pair, say (k,m), is chosen randomly from this set, and objects in positions k and m are switched. We omit some of the details of the algorithm. Similar implementations of the ITS technique for other optimization problems are thoroughly described, for example, in [21, 22]. We should note that also other enhancements of the generic tabu search strategy could be used for our purposes. These include, for example, multipath adaptive tabu search [23] and bacterial foraging tabu search [24] algorithms.In the rest of this section, we concentrate on how to efficiently compute the differencesΔ(p,k,m). Certainly, this issue is central to the performance of the tabu search procedure. But more importantly, evaluation of pairwise interchanges of objects stands at the heart of the dominance test, which is one of the crucial components of our algorithm. As we will see in the next section, the computation of the Δ values is an operation the dominance test is based on.Brusco [19] derived the following formula to calculate the difference between F(pkm) and F(p):
(11)Δ(p,k,m)=(4∑l=k+1m‍dp(k)p(l))((∑l=k+1m‍dp(k)p(l))+2fp(k)-Sp(k))+(4∑l=km-1‍dp(l)p(m))((∑l=km-1‍dp(l)p(m))-2fp(m)+Sp(m))+∑l=k+1m-1‍4(dp(l)p(m)-dp(k)p(l))((dp(l)p(m)-dp(k)p(l))+∑l=k+1m-1‍4(dp(l)p(m)-dp(k)p(l))+2fp(l)-Sp(l)(dp(l)p(m)-dp(k)p(l))).
Using (11), Δ(p,k,m) for all permutations in N2(p) can be computed in O(n3) time. However, this approach is actually too slow to be efficient. In order to describe a more elaborate approach we introduce three auxiliary n×n matrices A=(aiq), B=(biq), and C=(cij) defined for a permutation p∈Π. Let us consider an object i=p(s). The entries of the ith row of the matrices A and B are defined as follows: if q<s, then aiq=∑l=qs-1‍dip(l), biq=∑l=qs-1‍dip(l)(dip(l)+2fp(l)-Sp(l)); if q>s, then aiq=∑l=s+1q‍dip(l), biq=∑l=s+1q‍dip(l)(-dip(l)+2fp(l)-Sp(l)); if q=s, then aiq=biq=0. Turning now to the matrix C, we suppose that j=p(t). Then the entries of C are given by the following formula:
(12)cij={∑l=min(s,t)+1max(s,t)-1‍dip(l)djp(l)if|s-t|>10if|s-t|=1.
Unlike A and B, the matrix C is symmetric. We use the matrices A, B, and C to cast (11) in a form more suitable for both the TS heuristic and the dominance test.Proposition 4.
Forp∈Π, k∈{1,…,n-1}, and m∈{k+1,…,n},
(13)Δ(p,k,m)=4ap(k)m(ap(k)m+2ap(k),1-Sp(k))+4ap(m)k(ap(m)k-2ap(m),1+Sp(m))+4(bp(m),k+1-bp(k),m-1-2cp(k)p(m)).Proof.
First, notice that the sums∑l=k+1m‍dp(k)p(l) and ∑l=km-1‍dp(l)p(m) in (11) are equal, respectively, to ap(k)m and ap(m)k. Next, it can be seen that fp(k)=ap(k),1 and fp(m)=ap(m),1. Finally, the third term in (11), ignoring the factor of 4, can be rewritten as
(14)∑l=k+1m-1‍dp(l)p(m)(dp(l)p(m)+2fp(l)-Sp(l))-∑l=k+1m-1‍dp(l)p(m)dp(k)p(l)-∑l=k+1m-1‍dp(k)p(l)dp(l)p(m)-∑l=k+1m-1‍dp(k)p(l)(-dp(k)p(l)+2fp(l)-Sp(l)).
It is easy to see that the first sum in (14) is equal to bp(m),k+1, each of the next two sums is equal to cp(k)p(m), and the last sum is equal to bp(k),m-1. After simple manipulations, we arrive at (13).Having the matricesA,B,andC and using (13), the 2-interchange neighborhood of p can be explored in O(n2) time. The question that remains to be addressed is how efficiently the entries of these matrices can be updated after the replacement of p by a permutation p′ chosen from the neighborhood N2(p). Suppose that p′ is obtained from p by interchanging the objects in positions k and m. Then, we can write the following formulas that are involved in updating the matrices A and B:
(15)aiq:=ai,q-1+dip′(q),(16)aiq:=ai,q+1+dip′(q),(17)biq:=bi,q-1+dip′(q)(2ap′(q),1-Sp′(q)-dip′(q)),(18)biq:=bi,q+1+dip′(q)(2ap′(q),1-Sp′(q)+dip′(q)).
From now on, we will assume that the positions s and t of the objects i and j are defined with respect to the permutation p′. More formally, we have i=p′(s), j=p′(t). The procedure for updating A,B, and C is summarized in the following: (1)
for each objecti∈{1,…,n} do the following:
(1.1)
ifs<k, then compute aiq, q=k,…,m-1, by (15) and biq, q=k,…,n, by (17);(1.2)
ifs>m, then compute aiq, q=m,m-1,…,k+1, by (16) and biq, q=m,m-1,…,1, by (18);(1.3)
ifk<s<m, then compute aiq by (15) for q=m,…,n and by (16) for q=k,k-1,…,1;(1.4)
ifs=k or s=m, then first set ais:=0 and afterwards compute aiq by (15) for q=s+1,…,n and by (16) for q=s-1,s-2,…,1;(1.5)
ifk⩽s⩽m, then first set bis:=0 and afterwards compute biq by (17) for q=s+1,…,n and by (18) for q=s-1,s-2,…,1;(2)
for each pair of objectsi,j∈{1,…,n} such that s<t do the following:
(2.1)
ifs<k<t<m, then set cij:=cij+dip′(k)djp′(k)-dip′(m)djp′(m);(2.2)
ifk<s<m<t, then set cij:=cij+dip′(m)djp′(m)-dip′(k)djp′(k);(2.3)
if|{s,t}∩{k,m}|=1, then compute cij anew using (12);(2.4)
ifcij has been updated in (2.1)–(2.3), then set cji:=cij.It can be noted that for some pairs of objectsi and j the entry cij remains unchanged. There are five such cases: (1)s,t<k; (2)k<s,t<m; (3)s,t>m; (4)s<k<m<t; (5){s,t}={k,m}.We now evaluate the time complexity of the described procedure. Each of Steps(1.1) to (1.5) requires O(n) time. Thus, the complexity of Step (1) is O(n2). The number of operations in Steps (2.1), (2.2), and (2.4) is O(1), and that in Step (2.3) is O(n). However, the latter is executed only O(n) times. Therefore, the overall complexity of Step (2), and hence of the entire procedure, is O(n2).Formula (13) and the previous described procedure for updating matrices A, B, and C lie at the core of our implementation of the TS method. The matrices A, B, and C can be initialized by applying the formulas (15)–(18) and (12) with respect to the starting permutation. The time complexity of initialization is O(n2) for A and B, and O(n3) for C. The auxiliary matrices A, B, and C and the expression (13) for Δ(p,k,m) are also used in the dominance test presented in the next section. This test plays a major role in pruning the search tree during the branching and bounding process.


## 4. Dominance Test

As is well known, the use of dominance tests in enumerative methods allows reducing the search space of the optimization problem of interest. A comprehensive review of dominance tests (or rules) in combinatorial optimization can be found in the paper by Jouglet and Carlier [25]. For the problem we are dealing with in this study, Brusco and Stahl [13] proposed a dominance test for discarding partial solutions, which relies on performing pairwise interchanges of objects. To formulate it, let us assume that p=(p(1),…,p(r)), r<n, is a partial solution to the UDSP. The test checks whether Δ(p,k,r)>0 for at least one of the positions k∈{1,…,r-1}. If this condition is satisfied, then we say that the test fails. In this case, p is thrown away. Otherwise, p needs to be extended by placing an unassigned object in position r+1. The first condition in our dominance test is the same as in the test suggested by Brusco and Stahl. To compute Δ(p,k,r), we use subsets of entries of the auxiliary matrices A, B, and C defined in the previous section. Since only the sign of Δ(p,k,r) matters, the factor of 4 can be ignored in (13), and thus Δ′(p,k,r)=Δ(p,k,r)/4 instead of Δ(p,k,r) can be considered. From (13) we obtain
(19)Δ′(p,k,r)=ap(k)r(ap(k)r+2ap(k),1-Sp(k))Δ′(p,k,r)=+ap(r)k(ap(r)k-2ap(r),1+Sp(r))Δ′(p,k,r)=+bp(r),k+1-bp(k),r-1-2cp(k)p(r).We also take advantage of the case whereΔ′(p,k,r)=0. We use the lexicographic rule for breaking ties. The test outputs a “Fail” verdict if the equality Δ′(p,k,r)=0 is found to hold true for some k∈{1,…,r-1} such that p(k)>p(r). Indeed, p=(p(1),…,p(k-1),p(k),p(k+1),…,p(r)) can be discarded because both p and partial solution p′=(p(1),…,p(k-1),p(r),p(k+1),…,p(k)) are equally good, but p is lexicographically larger than p′.We note that the entries of the matricesA, B, and C used in (19) can be easily obtained from a subset of entries of A, B, and C computed earlier in the search process. Consider, for example, cp(k)p(r) in (19). It can be assumed that k∈{1,…,r-2} because if k=r-1, then cp(k)p(r)=0. For simplicity, let p(r)=i and p(r-1)=j. To avoid ambiguity, let us denote by C′ the matrix C that is obtained during processing of the partial solution (p(1),…,p(r-2)). The entry cp(k)i′ of C′ is calculated by (temporarily) placing the object i in position r-1. Now, referring to the definition of C, it is easy to see that cp(k)i=cp(k)i′+dp(k)jdij. Similarly, using formulas of types (15) and (17), we can obtain the required entries of the matrices A and B. Thus (19) can be computed in constant time. The complexity of computing Δ′(p,k,r) for all k is O(n). The dominance test used in [13] relies on (11) and has complexity O(n2). So our implementation of the dominance test is significantly more time-efficient than that proposed by Brusco and Stahl.In order to make the test even more efficient we, in addition, decided to evaluate those partial solutions that can be obtained fromp by relocating the object p(r) from position r to position k∈{1,…,r-2} (here position k=r-1 is excluded since in this case relocating would become a pairwise interchange of objects). Thus, the test compares p against partial solutions pk=(p(1),…,p(k-1),p(r),p(k),…,p(r-1)), k=1,…,r-2. Now suppose that p and pk are extended to n-permutations p~∈Π and, respectively, p~k∈Π such that, for each l∈{r+1,…,n}, the object in position l of p~ coincides with that in position l of p~k. Let us define δrk to be the difference between F(p~k) and F(p~). This difference can be computed as follows:
(20)δrk=(4∑l=kr-1‍dp(r)p(l))((∑l=kr-1‍dp(r)p(l))-2fp(r)+Sp(r))δrk=+4∑l=kr-1‍dp(r)p(l)(dp(r)p(l)+2fp(l)-Sp(l)).
The previous expression is similar to the formula (6) in [19]. For completeness, we give its derivation in the Appendix. We can rewrite (20) using the matrices A and B.Proposition 5.
Forp∈Π, k∈{1,…,n-1}, and r∈{k+1,…,n},
(21)δrk=4ap(r)k(ap(r)k-2ap(r),1+Sp(r))+4bp(r)k.The partial solutionp is dominated by pk if δrk>0. Upon establishing this fact for some k∈{1,…,r-2}, the test immediately stops with the rejection of p. Of course, like in the case of (13), a factor of 4 in (21) can be dropped. The computational overhead of including the condition δrk>0 in the test is minimal because all the quantities appearing on the right-hand side of (21), except bp(r),1, are also used in the calculation of Δ′(p,k,r) (see (19)). Since the loop over k needs to be executed the time complexity of this part of the test is O(n). Certainly, evaluation of relocating the last added object in p gives more power to the test. Thus our dominance test is an improvement over the test described in the Brusco and Stahl [13] paper.


## 5. A Branch-and-Bound Algorithm

Our approach to the UDSP is to start with a good initial solution provided by the iterated tabu search heuristic invoked in the initialization phase of the algorithm and then proceed with the main phase where, using the branch-and-bound technique, a search tree is built. For the purpose of pruning branches of the search tree, the algorithm involves three tests:(1) the upper bound test (based on U′); (2) the dominance test; (3) the symmetry test. The first two of them have already been described in the prior sections. The last one is the most simple and cheap to apply. It hinges on the fact that the values of the objective function F at a permutation p∈Π and its reverse p~ defined by p~(i)=p(n-i+1), i=1,…,n, are the same. Therefore, it is very natural to demand that some fixed object in W, say w, being among those objects that are assigned to the first half of the positions in p. More formally, the condition is i⩽⌈n/2⌉ where i is such that p(i)=w. If this condition is violated for p=(p(1),…,p(⌈n/2⌉)), then the symmetry test fails, and p is discarded. Because of the requirement to fulfill the previous condition, some care should be taken in applying the dominance test. In the case of Δ′(p,k,r)=0 and p(k)>p(r) (see the paragraph following (19)), this test in our implementation of the algorithm gives a “Fail” verdict only if either r⩽⌈n/2⌉ or p(k)≠w.Through the preliminary computational experiments, we were convinced that the upper bound test was quite weak when applied to smaller partial solutions. In other words, the time overhead for computing linear assignment-based bounds did not pay off for partial permutationsp=(p(1),…,p(r)) of length r that is small or modest compared to n. Therefore, it is reasonable to compute the upper bound U′ (described in Section 2) only when r⩾λ(D), where λ(D) is a threshold value depending on the UDSP instance matrix D and thus also on n. The problem we face is how to determine an appropriate value for λ(D). We combat this problem with a simple procedure consisting of two steps. In the first step, it randomly generates a relatively small number of partial permutations p1,…,pτ, each of length ⌊n/2⌋. For each permutation in this sample, the procedure computes the upper bound U′. Let U1′,…,Uτ′ be the values obtained. We are interested in the normalized differences between these values and F(pheur), where pheur is the solution delivered by the ITS algorithm. So the procedure calculates the sum Z=∑i=1τ‍100(Ui′-F(pheur))/F(pheur) and the average R=Z/τ of such differences. In the second step, using R, it determines the threshold value. This is done in the following way: if R<R1, then λ(D)=⌊n/3⌋; if R1⩽R<R2, then λ(D)=⌊n/2⌋; if R>R2, then λ(D)=⌊3n/4⌋. The described procedure is executed at the initialization step of the algorithm. Its parameters are R1, R2, and the sample size τ.Armed with tools for pruning the search space, we can construct a branch-and-bound algorithm for solving the UDSP. In order to make the description of the algorithm more readable we divide it into two parts, namely, the main procedure B&B (branch-and-bound) and an auxiliary procedure SelectObject. The former builds a search tree with nodes representing partial solutions. The task of SelectObject is either to choose an object, which can be used to extend the current permutation, or, if the current node becomes fathomed, to perform a backtracking step. In the description of the algorithm given in the following, the set of objects that can be used to create new branches is denoted byV. The main procedure goes as follows.B&B(1)
Apply the ITS algorithm to get an initial solutionp* to a given instance of the UDSP.SetF*≔F(p*).(2)
Calculateλ(D).(3)
Initialize the search tree by creating the root node with the empty setV attached to it.SetL≔0 (L denotes the tree level the currently considered node belongs to).(4)
Setr≔L.(5)
Call SelectObject(r). It returns L and, if L=r, additionally some unassigned object (let it be denoted by v). If L<0, then go to (7). Otherwise, check whether L<r. If so, then return to (4) (perform a backtracking step). If not (in which case L=r), then proceed to (6).(6)
Create a new node with emptyV. Set p(r+1)≔v. Increment L by 1, and go to (4).(7)
Calculate the coordinates for the objects by applying formula (4) with respect to the permutation p* and stop.Steps(1) and (2) of B&B comprise the initialization phase of the algorithm. The initial solution to the problem is produced using ITS. This solution is saved as the best solution found so far, denoted by p*. If later a better solution is found, p* is updated. This is done within the SelectObject procedure. In the final step, B&B uses p* to calculate the optimal coordinates for the objects. The branch-and-bound phase starts at Step (3) of B&B. Each node of the search tree is assigned a set V. The partial solution corresponding to the current node is denoted by p. The objects of V are used to extend p. Upon creation of a node (Step (3) for root and Step (6) for other nodes), the set V is empty. In SelectObject, V is filled with objects that are promising candidates to perform the branching operation on the corresponding node. SelectObject also decides which action has to be taken next. If L=r, then the action is “branch,” whereas if L<r, the action is “backtrack.” Step (4) helps to discriminate between these two cases. If L=r, then a new node in the next level of the search tree is generated. At that point, the partial permutation p is extended by assigning the selected object v to the position r+1 in p (see Step (6)). If, however, L<r and r is positive, then backtracking is performed by repeating Steps (4) and (5). We thus see that the search process is controlled by the procedure SelectObject. This procedure can be stated as follows.SelectObject(r) (1)
Check whether the setV attached to the current node is empty. If so, then go to (2). If not, then consider its subset V′={i∈V∣iisunmarkedandui(p)>F*}. If V′ is nonempty, then mark an arbitrary object v∈V′, and return with it as well as with L=r. Otherwise go to (6).(2)
Ifr>0, then go to (3). Otherwise, set V≔{1,…,n}, ui(p)≔∞, i=1,…,n, and go to (4).(3)
Ifr=n-1, then go to (5). Otherwise, for each i∈W-(p), perform the following operations. Temporarily set p(r+1)≔i. Apply the symmetry test, the dominance test, and, if r⩾λ(D), also the bound test to the resulting partial solution p′=(p(1),…,p(r+1)). If it passes all the tests, then do the following. Include i into V. If r⩾λ(D), assign to ui(p) the value of the upper bound U′ computed for the subproblem defined by the partial solution p′. If, however, r<λ(D), then set ui(p)≔∞. If after processing all objects in W-(p), the set V remains empty, then go to (6). Otherwise proceed to (4).(4)
Mark an arbitrary objectv∈V, and return with it as well as with L=r.(5)
LetW-(p)={j}. Calculate the value of the solution p=(p(1),…,p(n-1),p(n)=j). If F(p)>F*, then set p*:=p and F*≔F(p).(6)
Ifr>0, then drop the rth element of p. Return with L:=r-1.In the previous description,W-(p) denotes the set of unassigned objects, as before. The objects included in V⊆W-(p) are considered as candidates to be assigned to the leftmost available position in p. At the root node, V=W-(p). In all other cases, except when |W-(p)|=1, an object becomes a member of the set V only if it passes two or, for longer p, three tests. The third of them, the upper bound test, amounts to check if U′>F*. This condition is checked again on the second and subsequent calls to SelectObject for the same node of the search tree. This is done in Step (1) for each unmarked object i∈V. Such a strategy is reasonable because during the branch-and-bound process the value of F* may increase, and the previously satisfied constraint U′>F* may become violated for some i∈V. To be able to identify such objects, the algorithm stores the computed upper bound U′ as ui(p) at Step (3) of SelectObject.A closer look at the procedure SelectObject shows that two scenarios are possible:(1) when this procedure is invoked following the creation of a new node of the tree (which is done either in Step (3) or in Step (6) of B&B); (2) when a call to the procedure is made after a backtracking step. The set V initially is empty, and therefore Steps (2) to (5) of SelectObject become active only in the first of these two scenarios. Depending on the current level of the tree, r, one of the three cases is considered. If the current node is the root, that is, r=0, then Step (2) is executed. Since in this case the dominance test is not applicable, the algorithm uses the set V comprising all the objects under consideration. An arbitrary object selected in Step (4) is then used to create the first branch emanating from the root of the tree. If r∈[1,…,n-2], then the algorithm proceeds to Step (3). The task of this step is to evaluate all partial solutions that can be obtained from the current partial permutation p by assigning an object from W-(p) to the (r+1) th position in p. In order not to lose at least one optimal solution it is sufficient to keep only those extensions of p for which all the tests were successful. Such partial solutions are memorized by the set V. One of them is chosen in Step (4) to be used (in Step (6) of B&B) to generate the first son of the current node. If, however, no promising partial solution has been found and thus the set V is empty, then the current node can be discarded, and the search continued from its parent node. Finally, if r=n-1, then Step (5) is executed in which the permutation p is completed by placing the remaining unassigned object in the last position of p. Its value is compared with that of the best solution found so far. The node corresponding to p is a leaf and, therefore, is fathomed. In the second of the above-mentioned two scenarios (Step (1) of the procedure), the already existing nonempty set V is scanned. The previously selected objects in V have been marked. They are ignored. For each remaining object i, it is checked whether ui(p)>F*. If so, then i is included in the set V′. If the resulting set V′ is empty, then the current node of the tree becomes fathomed, and the procedure goes to Step (6). Otherwise, the procedure chooses a branching object and passes it to its caller B&B. Our implementation of the algorithm adopts a very simple enumeration strategy: branch upon an arbitrary object in V′ (or V when Step (4) is performed).It is useful to note that the tests in Step(3) are applied in increasing order of time complexity. The symmetry test is the cheapest of them, whereas the upper bound test is the most expensive to perform but still effective.


## 6. Computational Results

The main goal of our numerical experiments was to provide some measure of the performance of the developed B&B algorithm as well as to demonstrate the benefits of using the new upper bound and efficient dominance test throughout the solution process. For comparison purposes, we also implemented the better of the two branch-and-bound algorithms of Brusco and Stahl [13]. In their study, this algorithm is denoted as BB3. When referring to it in this paper, we will use the same name. In order to make the comparison more fair, we decided to start BB3 by invoking the ITS heuristic. Thus, the initial solution from which the branch-and-bound phase of each of the algorithms starts is expected to be of the same quality.


### 6.1. Experimental Setup

Both our algorithm and BB3 have been coded in the C programming language, and all the tests have been carried out on a PC with an Intel Core 2 Duo CPU running at 3.0 GHz. As a testbed for the algorithms, three empirical dissimilarity matrices from the literature in addition to a number of randomly generated ones were considered. The random matrices vary in size from 20 to 30 objects. All off-diagonal entries are drawn uniformly at random from the interval[1,100]. The matrices, of course, are symmetric.We also used three dissimilarity matrices constructed from empirical data found in the literature. The first matrix is based on Morse code confusion data collected by Rothkopf [26]. In his experiment, the subjects were presented pairs of Morse code signals and asked to respond whether they thought the signals were the same or different. These answers were translated to the entries of the similarity matrix with rows and columns corresponding to the 36 alphanumeric characters (26 letters in the Latin alphabet and 10 digits in the decimal number system). We should mention that we found small discrepancies between Morse code similarity matrices given in various publications and electronic media. We have decided to take the matrix from [27], where it is presented with a reference to Shepard [28]. Let this matrix be denoted by M=(mij). Its entries are integers from the interval [0,100]. Like in the study of Brusco and Stahl [13], the corresponding symmetric dissimilarity matrix D=(dij) was obtained by using the formula dij=200-mij-mji for all off-diagonal entries and setting the diagonal of D to zero. We considered two UDSP instances constructed using Morse code confusion data—one defined by the 26×26 submatrix of D with rows and columns labeled by the letters and another defined by the full matrix D.The second empirical matrix was taken from Groenen and Heiser [29]. We denote it by G=(gij). Its entry gij gives the number of citations in journal i to journal j. Thus the rows of G correspond to citing journals, and the columns correspond to cited journals. The dissimilarity matrix D was derived from G in two steps, like in [13]. First, a symmetric similarity matrix, H=(hij), was calculated by using the formulas hij=hji=(gij/∑l=1,l≠in‍gil)+(gji/∑l=1,l≠jn‍gjl), i,j=1,…,n, i<j, and hii=0, i=1,…,n. Then D was constructed with entries dij=⌊100(maxl,khlk-hij)⌋, i,j=1,…,n, i≠j, and dii=0, i=1,…,n.The third empirical matrix used in our experiments was taken from Hubert et al. [30]. Its entries represent the dissimilarities between pairs of food items. As mentioned by Hubert et al. [30], this 45×45 food-item matrix, Dfood, was constructed based on data provided by Ross and Murphy [31]. For our computational tests, we considered submatrices of Dfood defined by the first n∈{20,21,…,35} food items. All the dissimilarity matrices we just described as well as the source code of our algorithm are publicly available at http://www.proin.ktu.lt/~gintaras/uds.html.Based on the results of preliminary experiments with the B&B algorithm, we have fixed the values of its parameters:τ=10, R1=0, and R2=5. In the symmetry test, w was fixed at 1. The cutoff time for ITS was 1 second.


### 6.2. Results

The first experiment was conducted on a series of randomly generated full density dissimilarity matrices. The results of B&B and BB3 for them are summarized in Table1. The dimension of the matrix D is encoded in the instance name displayed in the first column. The second column of the table shows, for each instance, the optimal value of the objective function F* and, in parentheses, two related metrics. The first of them is the raw value of the stress function Φ(x). The second number in parentheses is the normalized value of the stress function that is calculated as Φ(x)/∑i<j‍dij. The next two columns list the results of our algorithm: the number of nodes in the search tree and the CPU time reported in the form hours : minutes : seconds or minutes : seconds or seconds. The last two columns give these measures for BB3.Table 1
Performance on random problem instances.
Inst.
Objective function value
B&B
BB3
Number of nodes
Time
Number of nodes
Time
p20
9053452 (179755.4; 0.284)
1264871
4.4
1526725
5.6
p21
9887394 (187927.7; 0.285)
3304678
10.6
4169099
14.5
p22
12548324 (224540.6; 0.282)
6210311
20.5
8017112
29.0
p23
14632006 (239619.8; 0.274)
12979353
45.0
16624663
1:06.7
p24
16211210 (281585.9; 0.294)
31820463
1:53.3
42435561
2:55.5
p25
15877026 (312670.0; 0.330)
103208158
6:44.6
137525505
9:52.4
p26
20837730 (347147.8; 0.302)
200528002
13:35.4
267603807
20:43.6
p27
21464870 (358595.8; 0.311)
287794935
20:37.7
392405143
32:36.5
p28
23938264 (396136.0; 0.317)
581602877
44:09.3
812649051
1:13:05.4
p29
28613820 (453778.6; 0.315)
2047228373
2:50:57.2
3044704554
4:38:57.0
p30
30359134 (464556.9; 0.315)
3555327537
5:09:09.9
5216804247
8:55:04.6Table 2
Performance on empirical dissimilarity matrices.
Inst.
Objective function value
B&B
BB3
Number of  nodes
Time
Number of nodes
Time
M26
180577772 (1944889.1; 0.219)
108400003
7:43.2
201264785
16:14.3
J28
60653172 (716703.3; 0.249)
2019598593
3:48:00.6
6780517853
10:39:06.6
F20
19344116 (105143.2; 0.098)
236990
2.7
1971840
7.0
F21
22694178 (122997.9; 0.102)
568998
5.2
4716777
16.2
F22
26496098 (149594.9; 0.110)
1178415
10.6
9993241
35.7
F23
30853716 (175768.1; 0.116)
2524703
25.2
24357098
1:37.8
F24
35075372 (217471.2; 0.130)
6727692
1:00.1
51593000
3:30.4
F25
40091544 (249100.2; 0.134)
14560522
2:28.3
139580891
10:01.4
F26
45507780 (288604.8; 0.142)
40465339
6:55.3
298597411
23:18.8
F27
51859210 (333682.5; 0.148)
116161004
19:29.8
670492972
55:57.8
F28
58712130 (378864.2; 0.153)
253278124
44:35.8
1366152797
2:03:26.5
F29
66198704 (420398.1; 0.156)
530219761
1:42:26.5
3616257778
5:32:08.4
F30
74310052 (462935.3; 0.157)
972413816
3:20:56.6
9360343618
15:56:41.2As it can be seen from the table, B&B is superior to BB3. We observe, for example, that the decrease in CPU time over BB3 is around40% for the largest three instances. Also, our algorithm builds smaller search trees than BB3.In our second experiment, we opt for three empirical dissimilarity matrices mentioned in Section6.1. Table 2 summarizes the results of B&B and BB3 on the UDSP instances defined by these matrices. Its structure is the same as that of Table 1. In the instance names, the letter “M” stands for the Morse code matrix, “J” stands for the journal citations matrix, and “F” stands for the food item data. As before, the number following the letter indicates the dimension of the problem matrix. Comparing results in Tables 1 and 2, we find that for empirical data the superiority of our algorithm over BB3 is more pronounced. Basically, B&B performs significantly better than BB3 for each of the three data sets. For example, for the food-item dissimilarity matrix of size 30×30, B&B builds almost 10 times smaller search tree and is about 5 times faster than BB3.Table3 presents the results of the performance of B&B on the full Morse code dissimilarity matrix as well as on larger submatrices of the food-item dissimilarity matrix. Running BB3 on these UDSP instances is too expensive, so this algorithm is not included in the table. By analyzing the results in Tables 2 and 3, we find that, for i∈{22,…,35}, the CPU time taken by B&B to solve the problem with the i×i food-item submatrix is between 196 (for i=30) and 282 (for i=27) percent of that needed to solve the problem with the (i-1)×(i-1) food-item submatrix. For the submatrix of size 35×35, the CPU time taken by the algorithm was more than 14 days. In fact, a solution of value 120526370 was found by the ITS heuristic in less than 1 s. The rest of the time was spent on proving its optimality. We also see from Table 3 that the UDSP with the full Morse code dissimilarity matrix required nearly 21 days to solve to proven optimality using our approach.Table 3
Performance of B&B on larger empirical dissimilarity matrices.
Instance
Objective function value
Number of  nodes
Time
F31
82221440 (519445.5; 0.164)
2469803824
8:57:03.8
F32
91389180 (570324.1; 0.166)
5040451368
19:44:12.7
F33
100367800 (641239.5; 0.174)
14318454604
55:24:20.0
F34
110065196 (710405.4; 0.180)
33035826208
129:36:50.0
F35
120526370 (784156.4; 0.185)
82986984742
342:56:45.6
M36
489111494 (4068295.6; 0.230)
243227368848
494:51:17.5In Table4, we display optimal solutions for the two largest UDSP instances we have tried. The second and fourth columns of the table contain the optimal permutations for these instances. The coordinates for objects, calculated using (4), are given in the third and fifth columns. The number in parenthesis in the fourth column shows the order in which the food items appear within Table 5.1 in the book by Hubert et al. [30]. From Table 4, we observe that only a half of the signals representing digits in the Morse code are placed at the end of the scale. Other such signals are positioned between signals encoding letters. We also see that shorter signals tend to occupy positions at the beginning of the permutation p.Table 4
Optimal solutions for the full Morse code dissimilarity matrix and the 35× 35 food-item dissimilarity matrix.
i
Morse code full
Food item data
p
(
i
)
x
p
(
i
)
p
(
i
)
x
p
(
i
)
1
(E)•
0.00000
orange (3)
0.00000
2
(T) −
5.83333
watermelon (2)
0.28571
3
(I)••
24.77778
apple (1)
0.28571
4
(A)•-
34.08333
banana (4)
1.28571
5
(N)-•
44.11111
pineapple (5)
1.62857
6
(M) − −
52.33333
lettuce (6)
21.25714
7
(S)•••
70.30556
broccoli (7)
22.31429
8
(U)••-
83.13889
carrots (8)
22.54286
9
(R)•-•
93.47222
corn (9)
22.88571
10
(W)•--
102.97222
onions (10)
25.11429
11
(H)••••
114.44444
potato (11)
29.97143
12
(D)-••
124.30556
rice (12)
51.62857
13
(K)-•-
131.52778
spaghetti (19)
58.97143
14
(V)•••-
145.00000
bread (13)
64.48571
15
(5)•••••
152.44444
bagel (14)
69.85714
16
(4)••••-
159.91667
cereal (16)
72.11429
17
(F)••-•
170.75000
oatmeal (15)
72.51429
18
(L)•-••
182.25000
pancake (18)
76.37143
19
(B)-•••
189.52778
muffin (17)
77.80000
20
(X)-••-
199.55556
crackers (20)
88.97143
21
(6)-••••
205.66667
granola bar (21)
93.00000
22
(3)•••--
220.13889
pretzels (22)
100.74286
23
(C)-•-•
229.47222
nuts (24)
104.37143
24
(Y)-•--
238.41667
popcorn (23)
107.62857
25
(7)--•••
249.30556
potato chips (25)
112.48571
26
(Z)--••
254.88889
doughnuts (26)
120.08571
27
(Q)--•-
264.38889
pizza (31)
126.51429
28
(P)•--•
270.83333
cookies (27)
136.88571
29
(J)•---
282.94444
chocolate bar (29)
139.40000
30
(G)--•
292.47222
cake (28)
141.05714
31
(O) − − −
300.22222
pie (30)
143.54286
32
(2)••---
310.50000
ice cream (32)
152.00000
33
(8)---••
320.36111
yogurt (33)
157.11429
34
(1)•----
333.50000
butter (34)
161.80000
35
(9)----•
341.86111
cheese (35)
165.08571
36
(0)-----
350.27778In Table5, we evaluate the effectiveness of tests used to fathom partial solutions. Computational results are reported for a selection of UDSP instances. The second, third, and fourth columns of the table present 100Nsym/N, 100Ndom/N, and 100NUB/N, where Nsym (respectively, Ndom and NUB) is the number of partial permutations that are fathomed due to symmetry test (respectively, due to dominance test and due to upper bound test) in Step (3) of SelectObject, and N=Nsym+Ndom+NUB. The last three columns give the same set of statistics for the BB3 algorithm. Inspection of Table 5 reveals that Ndom is much larger than both Nsym and NUB. It can also be seen that Ndom values are consistently higher with BB3 than when using our approach. Comparing the other two tests, we observe that Nsym is larger than NUB in all cases except when applying B&B to the food-item dissimilarity submatrices of size 20×20 and 25×25. Generally, we notice that the upper bound test is most successful for the UDSP instances defined by the food-item matrix.Table 5
Relative performance of the tests.
Instance
B & B
BB3
Sym.
Dom.
UB
Sym.
Dom.
UB
p20
11.58a
88.24b
0.18c
11.11a
88.73b
0.16c
p25
8.86a
90.87b
0.27c
8.48a
91.43b
0.09c
p30
7.01a
92.84b
0.15c
6.29a
93.64b
0.07c
M26
11.40a
88.47b
0.13c
10.44a
89.52b
0.04c
J28
8.55a
89.23b
2.22c
5.68a
94.26b
0.06c
F20
8.01a
75.63b
16.36c
13.25a
85.42b
1.33c
F25
7.04a
81.27b
11.69c
10.31a
89.39b
0.30c
F30
8.95a
83.91b
7.14c
7.74a
92.16b
0.10c
F35
8.53a
87.23b
4.24c
—a
—b
—c
M36
8.82a
91.14b
0.04c
—a
—b
—c
aPercentage of partial solutions fathomed due to symmetry test.
bPercentage of partial solutions fathomed due to dominance test.
cPercentage of partial solutions fathomed due to bound test.We end this section with a discussion of two issues regarding the use of heuristic algorithms for the UDSP. From a practical point of view, heuristics, of course, are preferable to exact methods as they are fast and can be applied to large problem instances. If carefully designed, they are able to produce solutions of very high quality. The development of exact methods, on the other hand, is a theoretical avenue for coping with different problems. In the field of optimization, exact algorithms, for a problem of interest, provide not only a solution but also a certificate of its optimality. Therefore, one of possible applications of exact methods is to use them in the process of assessment of heuristic techniques. Specifically, optimal values can serve as a reference point with which the results of heuristics can be compared. However, obtaining an optimality certificate in many cases, including the UDSP, is very time consuming. In order to see the time differences between finding an optimal solution for UDSP with the help of heuristics and proving its optimality (as it is done by the branch-and-bound method described in this paper), we run the two heuristic algorithms, iterated tabu search and simulated annealing (SA) of Brusco [19], on all the problem instances used in the main experiment. For two random instances, namely, p21 and p25, SA was trapped in a local maximum. Therefore, we have enhanced our implementation of SA by allowing it to run in a multi start mode. This variation of SA located an optimal solution for the above-mentioned instances in the third and, respectively, second restarts. Both tested heuristics have proven to be very fast. In particular, ITS took in total only 2.95 seconds to find the global optima for all 30 instances we experimented with. Performing the same task by SA required 5.44 seconds. Such times are in huge contrast to what is seen in Tables 1–3. Basically, the long computation times recorded in these tables are the price we have to pay for obtaining not only a solution but also a guarantee of its optimality. We did not perform additional computational experiments with ITS and SA because investigation of heuristics for the UDSP is outside the scope of this paper.Another issue concerns the use of heuristics for providing initial solutions for the branch-and-bound technique. We have employed the iterated tabu search procedure for this purpose. Our choice was influenced by the fact that we could adopt in ITS the same formulas as those used in the dominance test, which is a key component of our approach. Obviously, there are other heuristics one can apply in place of ITS. As it follows from the discussion in the preceding paragraph, a good candidate for this is the simulated annealing algorithm of Brusco [19]. Suppose that B&B is run twice on the same instance from different initial solutions, both of which are optimal. Then it should be clear that in both cases the same number of tree nodes is explored. This essentially means that any fast heuristic capable of providing optimal solutions is perfect to be used in the initialization step of B&B. An alternative strategy is to use in this step randomly generated permutations. The effects of this simple strategy are mixed. For example, for problem instances F25 and F27, the computation time increased from 148.3 and, respectively, 1169.8 seconds (see Table 2) to 413.3 and, respectively, 2062.5 seconds. However, for example, for p25 and M26, a somewhat different picture was observed. In the case of random initial solution, B&B took 407.5 and 468.7 seconds for these instances, respectively. The computation times in Table 1 for p25 and Table 2 for M26 are only marginally shorter. One possible explanation of this lies in the fact that, as Table 5 indicates, the percentage of partial solutions fathomed due to bound test for both p25 and M26 is quite low. However, in general, it is advantageous to use a heuristic for obtaining an initial solution for the branch-and-bound algorithm, especially because this solution can be produced at a negligible cost with respect to the rest of the computation.


## 6.1. Experimental Setup

Both our algorithm and BB3 have been coded in the C programming language, and all the tests have been carried out on a PC with an Intel Core 2 Duo CPU running at 3.0 GHz. As a testbed for the algorithms, three empirical dissimilarity matrices from the literature in addition to a number of randomly generated ones were considered. The random matrices vary in size from 20 to 30 objects. All off-diagonal entries are drawn uniformly at random from the interval[1,100]. The matrices, of course, are symmetric.We also used three dissimilarity matrices constructed from empirical data found in the literature. The first matrix is based on Morse code confusion data collected by Rothkopf [26]. In his experiment, the subjects were presented pairs of Morse code signals and asked to respond whether they thought the signals were the same or different. These answers were translated to the entries of the similarity matrix with rows and columns corresponding to the 36 alphanumeric characters (26 letters in the Latin alphabet and 10 digits in the decimal number system). We should mention that we found small discrepancies between Morse code similarity matrices given in various publications and electronic media. We have decided to take the matrix from [27], where it is presented with a reference to Shepard [28]. Let this matrix be denoted by M=(mij). Its entries are integers from the interval [0,100]. Like in the study of Brusco and Stahl [13], the corresponding symmetric dissimilarity matrix D=(dij) was obtained by using the formula dij=200-mij-mji for all off-diagonal entries and setting the diagonal of D to zero. We considered two UDSP instances constructed using Morse code confusion data—one defined by the 26×26 submatrix of D with rows and columns labeled by the letters and another defined by the full matrix D.The second empirical matrix was taken from Groenen and Heiser [29]. We denote it by G=(gij). Its entry gij gives the number of citations in journal i to journal j. Thus the rows of G correspond to citing journals, and the columns correspond to cited journals. The dissimilarity matrix D was derived from G in two steps, like in [13]. First, a symmetric similarity matrix, H=(hij), was calculated by using the formulas hij=hji=(gij/∑l=1,l≠in‍gil)+(gji/∑l=1,l≠jn‍gjl), i,j=1,…,n, i<j, and hii=0, i=1,…,n. Then D was constructed with entries dij=⌊100(maxl,khlk-hij)⌋, i,j=1,…,n, i≠j, and dii=0, i=1,…,n.The third empirical matrix used in our experiments was taken from Hubert et al. [30]. Its entries represent the dissimilarities between pairs of food items. As mentioned by Hubert et al. [30], this 45×45 food-item matrix, Dfood, was constructed based on data provided by Ross and Murphy [31]. For our computational tests, we considered submatrices of Dfood defined by the first n∈{20,21,…,35} food items. All the dissimilarity matrices we just described as well as the source code of our algorithm are publicly available at http://www.proin.ktu.lt/~gintaras/uds.html.Based on the results of preliminary experiments with the B&B algorithm, we have fixed the values of its parameters:τ=10, R1=0, and R2=5. In the symmetry test, w was fixed at 1. The cutoff time for ITS was 1 second.


## 6.2. Results

The first experiment was conducted on a series of randomly generated full density dissimilarity matrices. The results of B&B and BB3 for them are summarized in Table1. The dimension of the matrix D is encoded in the instance name displayed in the first column. The second column of the table shows, for each instance, the optimal value of the objective function F* and, in parentheses, two related metrics. The first of them is the raw value of the stress function Φ(x). The second number in parentheses is the normalized value of the stress function that is calculated as Φ(x)/∑i<j‍dij. The next two columns list the results of our algorithm: the number of nodes in the search tree and the CPU time reported in the form hours : minutes : seconds or minutes : seconds or seconds. The last two columns give these measures for BB3.Table 1
Performance on random problem instances.
Inst.
Objective function value
B&B
BB3
Number of nodes
Time
Number of nodes
Time
p20
9053452 (179755.4; 0.284)
1264871
4.4
1526725
5.6
p21
9887394 (187927.7; 0.285)
3304678
10.6
4169099
14.5
p22
12548324 (224540.6; 0.282)
6210311
20.5
8017112
29.0
p23
14632006 (239619.8; 0.274)
12979353
45.0
16624663
1:06.7
p24
16211210 (281585.9; 0.294)
31820463
1:53.3
42435561
2:55.5
p25
15877026 (312670.0; 0.330)
103208158
6:44.6
137525505
9:52.4
p26
20837730 (347147.8; 0.302)
200528002
13:35.4
267603807
20:43.6
p27
21464870 (358595.8; 0.311)
287794935
20:37.7
392405143
32:36.5
p28
23938264 (396136.0; 0.317)
581602877
44:09.3
812649051
1:13:05.4
p29
28613820 (453778.6; 0.315)
2047228373
2:50:57.2
3044704554
4:38:57.0
p30
30359134 (464556.9; 0.315)
3555327537
5:09:09.9
5216804247
8:55:04.6Table 2
Performance on empirical dissimilarity matrices.
Inst.
Objective function value
B&B
BB3
Number of  nodes
Time
Number of nodes
Time
M26
180577772 (1944889.1; 0.219)
108400003
7:43.2
201264785
16:14.3
J28
60653172 (716703.3; 0.249)
2019598593
3:48:00.6
6780517853
10:39:06.6
F20
19344116 (105143.2; 0.098)
236990
2.7
1971840
7.0
F21
22694178 (122997.9; 0.102)
568998
5.2
4716777
16.2
F22
26496098 (149594.9; 0.110)
1178415
10.6
9993241
35.7
F23
30853716 (175768.1; 0.116)
2524703
25.2
24357098
1:37.8
F24
35075372 (217471.2; 0.130)
6727692
1:00.1
51593000
3:30.4
F25
40091544 (249100.2; 0.134)
14560522
2:28.3
139580891
10:01.4
F26
45507780 (288604.8; 0.142)
40465339
6:55.3
298597411
23:18.8
F27
51859210 (333682.5; 0.148)
116161004
19:29.8
670492972
55:57.8
F28
58712130 (378864.2; 0.153)
253278124
44:35.8
1366152797
2:03:26.5
F29
66198704 (420398.1; 0.156)
530219761
1:42:26.5
3616257778
5:32:08.4
F30
74310052 (462935.3; 0.157)
972413816
3:20:56.6
9360343618
15:56:41.2As it can be seen from the table, B&B is superior to BB3. We observe, for example, that the decrease in CPU time over BB3 is around40% for the largest three instances. Also, our algorithm builds smaller search trees than BB3.In our second experiment, we opt for three empirical dissimilarity matrices mentioned in Section6.1. Table 2 summarizes the results of B&B and BB3 on the UDSP instances defined by these matrices. Its structure is the same as that of Table 1. In the instance names, the letter “M” stands for the Morse code matrix, “J” stands for the journal citations matrix, and “F” stands for the food item data. As before, the number following the letter indicates the dimension of the problem matrix. Comparing results in Tables 1 and 2, we find that for empirical data the superiority of our algorithm over BB3 is more pronounced. Basically, B&B performs significantly better than BB3 for each of the three data sets. For example, for the food-item dissimilarity matrix of size 30×30, B&B builds almost 10 times smaller search tree and is about 5 times faster than BB3.Table3 presents the results of the performance of B&B on the full Morse code dissimilarity matrix as well as on larger submatrices of the food-item dissimilarity matrix. Running BB3 on these UDSP instances is too expensive, so this algorithm is not included in the table. By analyzing the results in Tables 2 and 3, we find that, for i∈{22,…,35}, the CPU time taken by B&B to solve the problem with the i×i food-item submatrix is between 196 (for i=30) and 282 (for i=27) percent of that needed to solve the problem with the (i-1)×(i-1) food-item submatrix. For the submatrix of size 35×35, the CPU time taken by the algorithm was more than 14 days. In fact, a solution of value 120526370 was found by the ITS heuristic in less than 1 s. The rest of the time was spent on proving its optimality. We also see from Table 3 that the UDSP with the full Morse code dissimilarity matrix required nearly 21 days to solve to proven optimality using our approach.Table 3
Performance of B&B on larger empirical dissimilarity matrices.
Instance
Objective function value
Number of  nodes
Time
F31
82221440 (519445.5; 0.164)
2469803824
8:57:03.8
F32
91389180 (570324.1; 0.166)
5040451368
19:44:12.7
F33
100367800 (641239.5; 0.174)
14318454604
55:24:20.0
F34
110065196 (710405.4; 0.180)
33035826208
129:36:50.0
F35
120526370 (784156.4; 0.185)
82986984742
342:56:45.6
M36
489111494 (4068295.6; 0.230)
243227368848
494:51:17.5In Table4, we display optimal solutions for the two largest UDSP instances we have tried. The second and fourth columns of the table contain the optimal permutations for these instances. The coordinates for objects, calculated using (4), are given in the third and fifth columns. The number in parenthesis in the fourth column shows the order in which the food items appear within Table 5.1 in the book by Hubert et al. [30]. From Table 4, we observe that only a half of the signals representing digits in the Morse code are placed at the end of the scale. Other such signals are positioned between signals encoding letters. We also see that shorter signals tend to occupy positions at the beginning of the permutation p.Table 4
Optimal solutions for the full Morse code dissimilarity matrix and the 35× 35 food-item dissimilarity matrix.
i
Morse code full
Food item data
p
(
i
)
x
p
(
i
)
p
(
i
)
x
p
(
i
)
1
(E)•
0.00000
orange (3)
0.00000
2
(T) −
5.83333
watermelon (2)
0.28571
3
(I)••
24.77778
apple (1)
0.28571
4
(A)•-
34.08333
banana (4)
1.28571
5
(N)-•
44.11111
pineapple (5)
1.62857
6
(M) − −
52.33333
lettuce (6)
21.25714
7
(S)•••
70.30556
broccoli (7)
22.31429
8
(U)••-
83.13889
carrots (8)
22.54286
9
(R)•-•
93.47222
corn (9)
22.88571
10
(W)•--
102.97222
onions (10)
25.11429
11
(H)••••
114.44444
potato (11)
29.97143
12
(D)-••
124.30556
rice (12)
51.62857
13
(K)-•-
131.52778
spaghetti (19)
58.97143
14
(V)•••-
145.00000
bread (13)
64.48571
15
(5)•••••
152.44444
bagel (14)
69.85714
16
(4)••••-
159.91667
cereal (16)
72.11429
17
(F)••-•
170.75000
oatmeal (15)
72.51429
18
(L)•-••
182.25000
pancake (18)
76.37143
19
(B)-•••
189.52778
muffin (17)
77.80000
20
(X)-••-
199.55556
crackers (20)
88.97143
21
(6)-••••
205.66667
granola bar (21)
93.00000
22
(3)•••--
220.13889
pretzels (22)
100.74286
23
(C)-•-•
229.47222
nuts (24)
104.37143
24
(Y)-•--
238.41667
popcorn (23)
107.62857
25
(7)--•••
249.30556
potato chips (25)
112.48571
26
(Z)--••
254.88889
doughnuts (26)
120.08571
27
(Q)--•-
264.38889
pizza (31)
126.51429
28
(P)•--•
270.83333
cookies (27)
136.88571
29
(J)•---
282.94444
chocolate bar (29)
139.40000
30
(G)--•
292.47222
cake (28)
141.05714
31
(O) − − −
300.22222
pie (30)
143.54286
32
(2)••---
310.50000
ice cream (32)
152.00000
33
(8)---••
320.36111
yogurt (33)
157.11429
34
(1)•----
333.50000
butter (34)
161.80000
35
(9)----•
341.86111
cheese (35)
165.08571
36
(0)-----
350.27778In Table5, we evaluate the effectiveness of tests used to fathom partial solutions. Computational results are reported for a selection of UDSP instances. The second, third, and fourth columns of the table present 100Nsym/N, 100Ndom/N, and 100NUB/N, where Nsym (respectively, Ndom and NUB) is the number of partial permutations that are fathomed due to symmetry test (respectively, due to dominance test and due to upper bound test) in Step (3) of SelectObject, and N=Nsym+Ndom+NUB. The last three columns give the same set of statistics for the BB3 algorithm. Inspection of Table 5 reveals that Ndom is much larger than both Nsym and NUB. It can also be seen that Ndom values are consistently higher with BB3 than when using our approach. Comparing the other two tests, we observe that Nsym is larger than NUB in all cases except when applying B&B to the food-item dissimilarity submatrices of size 20×20 and 25×25. Generally, we notice that the upper bound test is most successful for the UDSP instances defined by the food-item matrix.Table 5
Relative performance of the tests.
Instance
B & B
BB3
Sym.
Dom.
UB
Sym.
Dom.
UB
p20
11.58a
88.24b
0.18c
11.11a
88.73b
0.16c
p25
8.86a
90.87b
0.27c
8.48a
91.43b
0.09c
p30
7.01a
92.84b
0.15c
6.29a
93.64b
0.07c
M26
11.40a
88.47b
0.13c
10.44a
89.52b
0.04c
J28
8.55a
89.23b
2.22c
5.68a
94.26b
0.06c
F20
8.01a
75.63b
16.36c
13.25a
85.42b
1.33c
F25
7.04a
81.27b
11.69c
10.31a
89.39b
0.30c
F30
8.95a
83.91b
7.14c
7.74a
92.16b
0.10c
F35
8.53a
87.23b
4.24c
—a
—b
—c
M36
8.82a
91.14b
0.04c
—a
—b
—c
aPercentage of partial solutions fathomed due to symmetry test.
bPercentage of partial solutions fathomed due to dominance test.
cPercentage of partial solutions fathomed due to bound test.We end this section with a discussion of two issues regarding the use of heuristic algorithms for the UDSP. From a practical point of view, heuristics, of course, are preferable to exact methods as they are fast and can be applied to large problem instances. If carefully designed, they are able to produce solutions of very high quality. The development of exact methods, on the other hand, is a theoretical avenue for coping with different problems. In the field of optimization, exact algorithms, for a problem of interest, provide not only a solution but also a certificate of its optimality. Therefore, one of possible applications of exact methods is to use them in the process of assessment of heuristic techniques. Specifically, optimal values can serve as a reference point with which the results of heuristics can be compared. However, obtaining an optimality certificate in many cases, including the UDSP, is very time consuming. In order to see the time differences between finding an optimal solution for UDSP with the help of heuristics and proving its optimality (as it is done by the branch-and-bound method described in this paper), we run the two heuristic algorithms, iterated tabu search and simulated annealing (SA) of Brusco [19], on all the problem instances used in the main experiment. For two random instances, namely, p21 and p25, SA was trapped in a local maximum. Therefore, we have enhanced our implementation of SA by allowing it to run in a multi start mode. This variation of SA located an optimal solution for the above-mentioned instances in the third and, respectively, second restarts. Both tested heuristics have proven to be very fast. In particular, ITS took in total only 2.95 seconds to find the global optima for all 30 instances we experimented with. Performing the same task by SA required 5.44 seconds. Such times are in huge contrast to what is seen in Tables 1–3. Basically, the long computation times recorded in these tables are the price we have to pay for obtaining not only a solution but also a guarantee of its optimality. We did not perform additional computational experiments with ITS and SA because investigation of heuristics for the UDSP is outside the scope of this paper.Another issue concerns the use of heuristics for providing initial solutions for the branch-and-bound technique. We have employed the iterated tabu search procedure for this purpose. Our choice was influenced by the fact that we could adopt in ITS the same formulas as those used in the dominance test, which is a key component of our approach. Obviously, there are other heuristics one can apply in place of ITS. As it follows from the discussion in the preceding paragraph, a good candidate for this is the simulated annealing algorithm of Brusco [19]. Suppose that B&B is run twice on the same instance from different initial solutions, both of which are optimal. Then it should be clear that in both cases the same number of tree nodes is explored. This essentially means that any fast heuristic capable of providing optimal solutions is perfect to be used in the initialization step of B&B. An alternative strategy is to use in this step randomly generated permutations. The effects of this simple strategy are mixed. For example, for problem instances F25 and F27, the computation time increased from 148.3 and, respectively, 1169.8 seconds (see Table 2) to 413.3 and, respectively, 2062.5 seconds. However, for example, for p25 and M26, a somewhat different picture was observed. In the case of random initial solution, B&B took 407.5 and 468.7 seconds for these instances, respectively. The computation times in Table 1 for p25 and Table 2 for M26 are only marginally shorter. One possible explanation of this lies in the fact that, as Table 5 indicates, the percentage of partial solutions fathomed due to bound test for both p25 and M26 is quite low. However, in general, it is advantageous to use a heuristic for obtaining an initial solution for the branch-and-bound algorithm, especially because this solution can be produced at a negligible cost with respect to the rest of the computation.


## 7. Conclusions

In this paper we have presented a branch-and-bound algorithm for the least-squares unidimensional scaling problem. The algorithm incorporates a LAP-based upper bound test as well as a dominance test which allows reducing the redundancy in the search process drastically. The results of computational experiments indicate that the algorithm can be used to obtain provably optimal solutions for randomly generated dissimilarity matrices of size up to30×30 and for empirical dissimilarity matrices of size up to 35×35. In particular, for the first time, the UDSP instance with the 36×36 Morse code dissimilarity matrix has been solved to guaranteed optimality. Another important instance is defined by the 45×45 food-item dissimilarity matrix Dfood. It is a great challenge to the optimization community to develop an exact algorithm that will solve this instance. Our branch-and-bound algorithm is limited to submatrices of Dfood of size around 35×35. We also remark that optimal solutions for all UDSP instances used in our experiments can be found by applying heuristic algorithms such as iterated tabu search and simulated annealing procedures. The total time taken by each of these procedures to reach the global optima is only a few seconds. An advantage of the branch-and-bound algorithm is its ability to certify the optimality of the solution obtained in a rigorous manner.Before closing, there are two more points worth of mention. First, we proposed a simple, yet effective, mechanism for making a decision about when it is useful to apply a bound test to the current partial solution and when such a test most probably will fail to discard this solution. For that purpose, a sample of partial solutions is evaluated in the preparation step of B&B. We believe that similar mechanisms may be operative in branch-and-bound algorithms for other difficult combinatorial optimization problems. Second, we developed an efficient procedure for exploring the pairwise interchange neighborhood of a solution in the search space. This procedure can be useful on its own in the design of heuristic algorithms for unidimensional scaling, especially those based on the local search principle.



---

*Source: 890589-2013-05-07.xml*