An Improved Exact Algorithm for Least-Squares Unidimensional Scaling

Abstract

Given n objects and an 𝑛×𝑛 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.

1. Introduction

Least-squares (or 𝐿₂) unidimensional scaling is an important optimization problem in the field of combinatorial data analysis. It can be stated as follows. Suppose that there are 𝑛 objects and, in addition, there are pairwise dissimilarity data available for them. These data form an 𝑛×𝑛 symmetric dissimilarity matrix, 𝐷=(𝑑_𝑖𝑗), 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 𝐷. Mathematically, the problem can be expressed as minΦ(𝑥)=^𝑛−1𝑛_{𝑖=1𝑗=𝑖+1}𝑑_𝑖𝑗−𝑥_𝑖−𝑥_𝑗², (1) where 𝑥 denotes the vector (𝑥₁,…,𝑥_𝑛) and 𝑥_𝑖, 𝑖∈{1,…,𝑛}, is the coordinate for the object 𝑖 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 𝑄 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: max_𝑝∈Π𝐹𝑝=^𝑛_𝑘=1^𝑘−1_𝑙=1𝑑_{𝑝(𝑘)𝑝(𝑙)}−^𝑛_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}², (2) where Π is the set of all permutations of 𝑊={1,…,𝑛}, that is, the set of all vectors (𝑝(1),𝑝(2),…,𝑝(𝑛)) such that 𝑝(𝑘)∈𝑊, 𝑘=1,…,𝑛, and 𝑝(𝑘)≠𝑝(𝑙), 𝑘=1,…,𝑛−1, 𝑙=𝑘+1,…,𝑛. Letting 𝑆_𝑖, 𝑖∈𝑊, denote the sum ∑^𝑛_{𝑗=1,𝑗≠𝑖}𝑑_𝑖𝑗, we can write 𝐹𝑝=^𝑛_𝑘=12^𝑘−1_𝑙=1𝑑_{𝑝(𝑘)𝑝(𝑙)}−𝑆_𝑝(𝑘)²=^𝑛_𝑘=12^𝑛_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}−𝑆_𝑝(𝑘)². (3) Thus, the UDSP is essentially a problem over all permutations of 𝑛 objects. With an optimal solution 𝑝^∗ to (2) at hand, we can compute the optimal coordinates for the objects by setting 𝑥_𝑝∗(1)=0 and using the following relationship [13] for 𝑘+1=2,…,𝑛: 𝑥_{𝑝∗(𝑘+1)}=𝑥_{𝑝∗(𝑘)}+^𝑛_𝑙=𝑘+1𝑑_{𝑝∗(𝑘)𝑝∗(𝑙)}−^𝑘−1_𝑙=1𝑑_{𝑝∗(𝑘)𝑝∗(𝑙)}−^𝑛_𝑙=𝑘+2𝑑_{𝑝∗(𝑘+1)𝑝∗(𝑙)}+^𝑘_𝑙=1𝑑_{𝑝∗(𝑘+1)𝑝∗(𝑙)}𝑛⁻¹. (4)

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 𝑛-permutation under construction. For example, at level 4 of the search tree, the four selected objects are 𝑝(1), 𝑝(𝑛), 𝑝(2), and 𝑝(𝑛−1), and there are no object assigned to positions 3,…,𝑛−2 in the permutation 𝑝. 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 𝑂(𝑛²) time, which is an improvement over a standard 𝑂(𝑛³)-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 Section 2, 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 solution 𝑝=(𝑝(1),…,𝑝(𝑟)), 𝑟∈{1,…,𝑛−1}, we denote by 𝑊(𝑝) the set of the objects in 𝑝. Thus 𝑊(𝑝)={𝑝(1),…,𝑝(𝑟)}, and 𝑊(𝑝)=𝑊⧵𝑊(𝑝) is the set of unselected objects. At the root node of the search tree, 𝑝 is empty, and 𝑊(𝑝) is equal to 𝑊. The first upper bound used by Brusco and Stahl, when adopted to 𝑝, takes the following form: 𝑈_BS∶=𝐹_𝑟𝑝+_{𝑖∈𝑊(𝑝)}𝑆²_𝑖, (5) where 𝐹_𝑟(𝑝) denotes the sum of the terms in (3) corresponding to assigned objects. Thus, 𝐹_𝑟(𝑝)=∑^𝑟_𝑘=1(𝑆_𝑝(𝑘)−2∑^𝑘−1_𝑙=1𝑑_{𝑝(𝑘)𝑝(𝑙)})².

To describe the second bound, we define 𝐷=(𝑑_𝑖𝑗) to be an 𝑛×(𝑛−1) matrix obtained from 𝐷 by sorting the nondiagonal entries of each row of 𝐷 in nondecreasing order. Let 𝜂_𝑖𝑘=(𝑆_𝑖−2∑^𝑘−1_𝑗=1𝑑_𝑖𝑗)² and 𝜇_𝑘=max_{1⩽𝑖⩽𝑛}𝜂_𝑖𝑘 for 𝑘=1,…,⌈𝑛/2⌉. For other values of 𝑘, 𝜇_𝑘 is obtained from the equality 𝜇_𝑘=𝜇_{𝑛−𝑘+1}, 𝑘∈{1,…,𝑛}. The second upper bound proposed in [13] is given by the following expression: 𝐹_𝑟𝑝+^𝑛_𝑘=𝑟+1𝜇_𝑘. (6)

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 max^𝑛−𝑟_𝑖=1^𝑛−𝑟_𝑗=1𝑒_𝑖𝑗𝑦_𝑖𝑗s.t.^𝑛−𝑟_𝑗=1𝑦_𝑖𝑗=1,𝑖=1,…,𝑛−𝑟,_{𝑛−𝑟𝑖=1}𝑦_𝑖𝑗𝑦=1,𝑗=1,…,𝑛−𝑟,_𝑖𝑗∈{0,1},𝑖,𝑗=1,…,𝑛−𝑟, (7) where 𝑒_𝑖𝑗, 𝑖,𝑗=1,…,𝑛−𝑟 are entries of an (𝑛−𝑟)×(𝑛−𝑟) profit matrix 𝐸 with rows and columns indexed by the objects in the set 𝑊(𝑝). Next, we will construct the matrix 𝐸. For a permutation 𝑝∈Π, we denote by 𝑓_𝑝(𝑘) (resp., 𝑓^′_𝑝(𝑘)) the sum of dissimilarities between object 𝑝(𝑘) and all objects preceding (resp., following) 𝑝(𝑘) in 𝑝. Thus, 𝑓_𝑝(𝑘)=∑^𝑘−1_𝑙=1𝑑_{𝑝(𝑘)𝑝(𝑙)}, 𝑓^′_𝑝(𝑘)=∑^𝑛_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}. From (3), we have 𝐹𝑝=^𝑛_𝑘=12max𝑓_𝑝(𝑘),𝑓^′_𝑝(𝑘)−𝑆_𝑝(𝑘)², (8) and 2max(𝑓_𝑝(𝑘),𝑓^′_𝑝(𝑘))⩾𝑆_𝑝(𝑘)=𝑓_𝑝(𝑘)+𝑓^′_𝑝(𝑘) for each 𝑘=1,…,𝑛. Suppose now that 𝑝=(𝑝(1),…,𝑝(𝑟)) is a partial solution and Π(𝑝) is the set of all 𝑛-permutations that can be obtained by extending 𝑝. Assume, for the sake of simplicity, that 𝑊(𝑝)={1,…,𝑛−𝑟}. In order to construct a suitable matrix 𝐸, we have to bound 𝑓_𝑝(𝑘) and 𝑓^′_𝑝(𝑘) in (8) from above. To this end, for any 𝑖∈𝑊(𝑝), let (𝑑^′_𝑖1,…,𝑑^′_{𝑖,𝑛−𝑟−1}) denote the list obtained by sorting 𝑑_𝑖𝑗, 𝑗∈𝑊(𝑝)⧵{𝑖}, in nonincreasing order. To get the required bounds, we make use of the sums 𝛾_𝑖𝑗=∑^𝑗_𝑙=1𝑑^′_𝑖1, 𝑖∈𝑊(𝑝), 𝑗=0,1,…,𝑛−𝑟−1, and 𝜌_𝑖𝑘=∑^𝑟_𝑙=1𝑑_{𝑖𝑝(𝑙)}+𝛾_{𝑖,𝑘−𝑟−1}, 𝑖∈𝑊(𝑝), 𝑘=𝑟+1,…,𝑛. Suppose that 𝑝 is any permutation in Π(𝑝) and 𝑝(𝑘)=𝑖 for a position 𝑘∈{𝑟+1,…,𝑛}. Then it is easy to see that 𝑓𝑝_(𝑘)=𝑓_𝑖⩽𝜌_𝑖𝑘 and 𝑓𝑝^′_(𝑘)=𝑓^′_𝑖⩽𝛾_{𝑖,𝑛−𝑘}. Indeed, in the case of 𝑓𝑝_(𝑘), for example, we have 𝑓𝑝_(𝑘)=∑^𝑟_𝑙=1𝑑_{𝑖𝑝(𝑙)}+∑^𝑘−1_𝑙=𝑟+1𝑑_𝑖𝑝_(𝑙), and the last term in this expression is less than or equal to 𝛾_{𝑖,𝑘−𝑟−1}. We define 𝐸 to be the matrix with entries 𝑒_𝑖𝑗=(2max(𝜌_{𝑖,𝑟+𝑗},𝛾_{𝑖,𝑛−𝑟−𝑗})−𝑆_𝑖)², 𝑖,𝑗=1,…,𝑛−𝑟. Let 𝑈^∗ denote the optimal value of the linear assignment problem (LAP) (7) with the profit matrix 𝐸. From the definition of 𝐸, we immediately obtain an upper bound on the optimal value of the objective function 𝐹.

Proposition 1. For a partial solution 𝑝=(𝑝(1),…,𝑝(𝑟)), 𝑈∶=𝐹_𝑟(𝑝)+𝑈^∗⩾max𝑝_∈Π(𝑝)𝐹(𝑝).

From the definition of 𝜌 and 𝛾 values, it can be observed that, with the increase of 𝑘, 𝜌_𝑖𝑘 increases, while 𝛾_{𝑖,𝑛−𝑘} decreases. Therefore, when constructing the 𝑖th row of the matrix 𝐸, it is expedient to first find the largest 𝑘∈{𝑟+1,…,𝑛} for which 𝜌_𝑖𝑘⩽𝛾_{𝑖,𝑛−𝑘}. We denote such 𝑘 by 𝑘 and assume that 𝑘=𝑟 if 𝜌_𝑖,𝑟+1>𝛾_{𝑖,𝑛−𝑟−1}. Then 𝑒_{𝑖,𝑘−𝑟}=(2𝛾_{𝑖,𝑛−𝑘}−𝑆_𝑖)² for 𝑘=𝑟+1,…,𝑘, and 𝑒_{𝑖,𝑘−𝑟}=(2𝜌_𝑖𝑘−𝑆_𝑖)² for 𝑘=𝑘+1,…,𝑛.

The most time-consuming step in the computation of the upper bound 𝑈 is solving the linear assignment problem. In order to provide a computationally less expensive bounding method we replace 𝑈^∗ with the value of a feasible solution to the dual of the LAP (7). To get the bound, we take 𝛼_𝑖=max_{1⩽𝑘⩽𝑛−𝑟}𝑒_𝑖𝑘=𝑒_{𝑖,𝑛−𝑟}=𝑆²_𝑖, 𝑖=1,…,𝑛−𝑟, 𝛽_𝑗=max_{1⩽𝑖⩽𝑛−𝑟}(𝑒_𝑖𝑗−𝛼_𝑖), 𝑗=1,…,𝑛−𝑟. Clearly, 𝛼_𝑖+𝛽_𝑗⩾𝑒_𝑖𝑗 for each pair 𝑖,𝑗=1,…,𝑛−𝑟. Thus the vector (𝛼₁,…,𝛼_𝑛−𝑟,𝛽₁,…,𝛽_𝑛−𝑟) satisfies the constraints of the dual of the LAP instance of interest. Let 𝑈=∑^𝑛−𝑟_𝑖=1𝛼_𝑖+∑^𝑛−𝑟_𝑗=1𝛽_𝑗. Then, from the duality theory of linear programming, we can establish the following upper bound for 𝐹(𝑝).

Proposition 2. For a partial solution 𝑝=(𝑝(1),…,𝑝(𝑟)), 𝑈^′∶=𝐹_𝑟(𝑝)+𝑈⩾max𝑝_∈Π(𝑝)𝐹(𝑝).

Proof. The inequality follows from Proposition 1 and the fact that 𝑈⩾𝑈^∗.

One point that needs to be emphasized is that 𝑈^′ strengthens the first bound proposed by Brusco and Stahl.

Proposition 3. 𝑈^′⩽𝑈_𝐵𝑆.

Proof. Indeed, 𝐹_𝑟(𝑝)+∑_{𝑖∈𝑊(𝑝)}𝑆²_𝑖=𝐹_𝑟(𝑝)+∑^𝑛−𝑟_𝑖=1𝛼_𝑖⩾𝐹_𝑟(𝑝)+𝑈. The equality in the previous chain of expressions comes from the choice of 𝛼_𝑖 and our assumption on 𝑊(𝑝), whereas the inequality is due to the fact that 𝛽_𝑗⩽0 for each 𝑗=1,…,𝑛−𝑟.

As an example, we consider the following 4×4 dissimilarity matrix borrowed from Defays [12]: −4234−7227−4324−. (9) 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 𝑝=(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, 𝑊(𝑝)={1,…,4} and 𝐹_𝑟(𝑝)=0. The matrix 𝐸 is constructed by noting that 𝑘=2 for each 𝑖=1,…,4. It takes the following form: 812525811698181169169818116981252581. (10) Using the expressions for 𝛼_𝑖 and 𝛽_𝑗, we find that 𝛼₁=𝛼₄=81, 𝛼₂=𝛼₃=169, 𝛽₁=𝛽₄=0, and 𝛽₂=𝛽₃=−56. Thus, 𝑈=∑⁴_𝑖=1𝛼_𝑖+∑⁴_𝑗=1𝛽_𝑗=500−112=388=max_𝑝∈Π𝐹(𝑝). Hence, our bound for the Defays instance is tight. For comparison, the bound (5) is equal to ∑⁴_𝑖=1𝛼_𝑖=500. To obtain the second bound of Brusco and Stahl, we first compute the 𝜇 values. Since 𝜇₁=𝜇₄=169 and 𝜇₂=𝜇₃=81, the bound (6) (taking the form ∑⁴_𝑖=1𝜇_𝑖) 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.

For 𝑝∈Π, let 𝑝_𝑘𝑚 denote the permutation obtained from 𝑝 by interchanging the objects in positions 𝑘 and 𝑚. The set 𝑁₂(𝑝)={𝑝_𝑘𝑚∣𝑘,𝑚=1,…,𝑛,𝑘<𝑚} is the 2-interchange neighborhood of 𝑝. At each iteration, the tabu search algorithm performs an exploration of the neighborhood 𝑁₂(𝑝). A subset of solutions in 𝑁₂(𝑝) are evaluated by computing Δ(𝑝,𝑘,𝑚)=𝐹(𝑝_𝑘𝑚)−𝐹(𝑝). 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 Δ(𝑝,𝑘,𝑚). We will defer this issue until later in this section. In step (2) of ITS, at each iteration of the TS procedure, the neighborhood 𝑁₂(𝑝) of the current solution 𝑝 is explored. The procedure selects a permutation 𝑝_𝑘𝑚∈𝑁₂(𝑝), which either improves the best available objective function value so far or, if this does not happen, has the largest value of Δ(𝑝,𝑖,𝑗) among all permutations 𝑝_𝑖𝑗∈𝑁₂(𝑝) for which the objects in positions 𝑖 and 𝑗 are not in the tabu list. The permutation 𝑝_𝑘𝑚 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 (𝑘,𝑚), is chosen randomly from this set, and objects in positions 𝑘 and 𝑚 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 Δ(𝑝,𝑘,𝑚). 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 𝐹(𝑝_𝑘𝑚) and 𝐹(𝑝): Δ=4𝑝,𝑘,𝑚_{𝑚𝑙=𝑘+1}𝑑_{𝑝(𝑘)𝑝(𝑙)}^𝑚_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}+2𝑓_𝑝(𝑘)−𝑆_𝑝(𝑘)+4^𝑚−1_𝑙=𝑘𝑑_{𝑝(𝑙)𝑝(𝑚)}^𝑚−1_𝑙=𝑘𝑑_{𝑝(𝑙)𝑝(𝑚)}−2𝑓_𝑝(𝑚)+𝑆_𝑝(𝑚)+^𝑚−1_𝑙=𝑘+14𝑑_{𝑝(𝑙)𝑝(𝑚)}−𝑑_{𝑝(𝑘)𝑝(𝑙)}𝑑_{𝑝(𝑙)𝑝(𝑚)}−𝑑_{𝑝(𝑘)𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙). (11) Using (11), Δ(𝑝,𝑘,𝑚) for all permutations in 𝑁₂(𝑝) can be computed in 𝑂(𝑛³) time. However, this approach is actually too slow to be efficient. In order to describe a more elaborate approach we introduce three auxiliary 𝑛×𝑛 matrices 𝐴=(𝑎_𝑖𝑞), 𝐵=(𝑏_𝑖𝑞), and 𝐶=(𝑐_𝑖𝑗) defined for a permutation 𝑝∈Π. Let us consider an object 𝑖=𝑝(𝑠). The entries of the 𝑖th row of the matrices 𝐴 and 𝐵 are defined as follows: if 𝑞<𝑠, then 𝑎_𝑖𝑞=∑^𝑠−1_𝑙=𝑞𝑑_{𝑖𝑝(𝑙)}, 𝑏_𝑖𝑞=∑^𝑠−1_𝑙=𝑞𝑑_{𝑖𝑝(𝑙)}(𝑑_{𝑖𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙)); if 𝑞>𝑠, then 𝑎_𝑖𝑞=∑^𝑞_𝑙=𝑠+1𝑑_{𝑖𝑝(𝑙)}, 𝑏_𝑖𝑞=∑^𝑞_𝑙=𝑠+1𝑑_{𝑖𝑝(𝑙)}(−𝑑_{𝑖𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙)); if 𝑞=𝑠, then 𝑎_𝑖𝑞=𝑏_𝑖𝑞=0. Turning now to the matrix 𝐶, we suppose that 𝑗=𝑝(𝑡). Then the entries of 𝐶 are given by the following formula: 𝑐_𝑖𝑗=^{max(𝑠,𝑡)−1}_()𝑙=m∑_{in𝑠,𝑡+1}𝑑_{𝑖𝑝(𝑙)}𝑑_{𝑗𝑝(𝑙)}||||if𝑠−𝑡>10if𝑠−𝑡=1. (12) Unlike 𝐴 and 𝐵, the matrix 𝐶 is symmetric. We use the matrices 𝐴, 𝐵, and 𝐶 to cast (11) in a form more suitable for both the TS heuristic and the dominance test.

Proposition 4. For 𝑝∈Π, 𝑘∈{1,…,𝑛−1}, and 𝑚∈{𝑘+1,…,𝑛}, Δ𝑝,𝑘,𝑚=4𝑎_{𝑝(𝑘)𝑚}𝑎_{𝑝(𝑘)𝑚}+2𝑎_𝑝(𝑘),1−𝑆_𝑝(𝑘)+4𝑎_{𝑝(𝑚)𝑘}𝑎_{𝑝(𝑚)𝑘}−2𝑎_𝑝(𝑚),1+𝑆_𝑝(𝑚)+4𝑏_{𝑝(𝑚),𝑘+1}−𝑏_{𝑝(𝑘),𝑚−1}−2𝑐_{𝑝(𝑘)𝑝(𝑚)}. (13)

Proof. First, notice that the sums ∑^𝑚_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)} and ∑^𝑚−1_𝑙=𝑘𝑑_{𝑝(𝑙)𝑝(𝑚)} in (11) are equal, respectively, to 𝑎_{𝑝(𝑘)𝑚} and 𝑎_{𝑝(𝑚)𝑘}. Next, it can be seen that 𝑓_𝑝(𝑘)=𝑎_𝑝(𝑘),1 and 𝑓_𝑝(𝑚)=𝑎_𝑝(𝑚),1. Finally, the third term in (11), ignoring the factor of 4, can be rewritten as ^𝑚−1_𝑙=𝑘+1𝑑_{𝑝(𝑙)𝑝(𝑚)}𝑑_{𝑝(𝑙)𝑝(𝑚)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙)−^𝑚−1_𝑙=𝑘+1𝑑_{𝑝(𝑙)𝑝(𝑚)}𝑑_{𝑝(𝑘)𝑝(𝑙)}−^𝑚−1_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}𝑑_{𝑝(𝑙)𝑝(𝑚)}−^𝑚−1_𝑙=𝑘+1𝑑_{𝑝(𝑘)𝑝(𝑙)}−𝑑_{𝑝(𝑘)𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙). (14) It is easy to see that the first sum in (14) is equal to 𝑏_{𝑝(𝑚),𝑘+1}, each of the next two sums is equal to 𝑐_{𝑝(𝑘)𝑝(𝑚)}, and the last sum is equal to 𝑏_{𝑝(𝑘),𝑚−1}. After simple manipulations, we arrive at (13).

Having the matrices 𝐴,𝐵,and𝐶 and using (13), the 2-interchange neighborhood of 𝑝 can be explored in 𝑂(𝑛²) time. The question that remains to be addressed is how efficiently the entries of these matrices can be updated after the replacement of 𝑝 by a permutation 𝑝^′ chosen from the neighborhood 𝑁₂(𝑝). Suppose that 𝑝^′ is obtained from 𝑝 by interchanging the objects in positions 𝑘 and 𝑚. Then, we can write the following formulas that are involved in updating the matrices 𝐴 and 𝐵: 𝑎_𝑖𝑞∶=𝑎_{𝑖,𝑞−1}+𝑑_{𝑖𝑝′(𝑞)},𝑎_𝑖𝑞∶=𝑎_𝑖,𝑞+1+𝑑_{𝑖𝑝′(𝑞)},𝑏_𝑖𝑞∶=𝑏_{𝑖,𝑞−1}+𝑑_{𝑖𝑝′(𝑞)}2𝑎_{𝑝′(𝑞),1}−𝑆_{𝑝′(𝑞)}−𝑑_{𝑖𝑝′(𝑞)},𝑏_𝑖𝑞∶=𝑏_𝑖,𝑞+1+𝑑_{𝑖𝑝′(𝑞)}2𝑎_{𝑝′(𝑞),1}−𝑆_{𝑝′(𝑞)}+𝑑_{𝑖𝑝′(𝑞)}. (15)(16)(17)(18) From now on, we will assume that the positions 𝑠 and 𝑡 of the objects 𝑖 and 𝑗 are defined with respect to the permutation 𝑝^′. More formally, we have 𝑖=𝑝^′(𝑠), 𝑗=𝑝^′(𝑡). The procedure for updating 𝐴,𝐵, and 𝐶 is summarized in the following: (1)for each object 𝑖∈{1,…,𝑛} do the following:  (1.1) if 𝑠<𝑘, then compute 𝑎_𝑖𝑞, 𝑞=𝑘,…,𝑚−1, by (15) and 𝑏_𝑖𝑞, 𝑞=𝑘,…,𝑛, by (17);  (1.2) if 𝑠>𝑚, then compute 𝑎_𝑖𝑞, 𝑞=𝑚,𝑚−1,…,𝑘+1, by (16) and 𝑏_𝑖𝑞, 𝑞=𝑚,𝑚−1,…,1, by (18);  (1.3) if 𝑘<𝑠<𝑚, then compute 𝑎_𝑖𝑞 by (15) for 𝑞=𝑚,…,𝑛 and by (16) for 𝑞=𝑘,𝑘−1,…,1;  (1.4) if 𝑠=𝑘 or 𝑠=𝑚, then first set 𝑎_𝑖𝑠∶=0 and afterwards compute 𝑎_𝑖𝑞 by (15) for 𝑞=𝑠+1,…,𝑛 and by (16) for 𝑞=𝑠−1,𝑠−2,…,1;  (1.5) if 𝑘⩽𝑠⩽𝑚, then first set 𝑏_𝑖𝑠∶=0 and afterwards compute 𝑏_𝑖𝑞 by (17) for 𝑞=𝑠+1,…,𝑛 and by (18) for 𝑞=𝑠−1,𝑠−2,…,1; (2)for each pair of objects 𝑖,𝑗∈{1,…,𝑛} such that 𝑠<𝑡 do the following:  (2.1)if 𝑠<𝑘<𝑡<𝑚, then set 𝑐_𝑖𝑗∶=𝑐_𝑖𝑗+𝑑_{𝑖𝑝′(𝑘)}𝑑_{𝑗𝑝′(𝑘)}−𝑑_{𝑖𝑝′(𝑚)}𝑑_{𝑗𝑝′(𝑚)};  (2.2)if 𝑘<𝑠<𝑚<𝑡, then set 𝑐_𝑖𝑗∶=𝑐_𝑖𝑗+𝑑_{𝑖𝑝′(𝑚)}𝑑_{𝑗𝑝′(𝑚)}−𝑑_{𝑖𝑝′(𝑘)}𝑑_{𝑗𝑝′(𝑘)};  (2.3)if |{𝑠,𝑡}∩{𝑘,𝑚}|=1, then compute 𝑐_𝑖𝑗 anew using (12);  (2.4)if 𝑐_𝑖𝑗 has been updated in (2.1)--(2.3), then set 𝑐_𝑗𝑖∶=𝑐_𝑖𝑗.

It can be noted that for some pairs of objects 𝑖 and 𝑗 the entry 𝑐_𝑖𝑗 remains unchanged. There are five such cases: (1)𝑠,𝑡<𝑘; (2)𝑘<𝑠,𝑡<𝑚; (3)𝑠,𝑡>𝑚; (4)𝑠<𝑘<𝑚<𝑡; (5){𝑠,𝑡}={𝑘,𝑚}.

We now evaluate the time complexity of the described procedure. Each of Steps (1.1) to (1.5) requires 𝑂(𝑛) time. Thus, the complexity of Step (1) is 𝑂(𝑛²). The number of operations in Steps (2.1), (2.2), and (2.4) is 𝑂(1), and that in Step (2.3) is 𝑂(𝑛). However, the latter is executed only 𝑂(𝑛) times. Therefore, the overall complexity of Step (2), and hence of the entire procedure, is 𝑂(𝑛²).

Formula (13) and the previous described procedure for updating matrices 𝐴, 𝐵, and 𝐶 lie at the core of our implementation of the TS method. The matrices 𝐴, 𝐵, and 𝐶 can be initialized by applying the formulas (15)–(18) and (12) with respect to the starting permutation. The time complexity of initialization is 𝑂(𝑛²) for 𝐴 and 𝐵, and 𝑂(𝑛³) for 𝐶. The auxiliary matrices 𝐴, 𝐵, and 𝐶 and the expression (13) for Δ(𝑝,𝑘,𝑚) 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 𝑝=(𝑝(1),…,𝑝(𝑟)), 𝑟<𝑛, is a partial solution to the UDSP. The test checks whether Δ(𝑝,𝑘,𝑟)>0 for at least one of the positions 𝑘∈{1,…,𝑟−1}. If this condition is satisfied, then we say that the test fails. In this case, 𝑝 is thrown away. Otherwise, 𝑝 needs to be extended by placing an unassigned object in position 𝑟+1. The first condition in our dominance test is the same as in the test suggested by Brusco and Stahl. To compute Δ(𝑝,𝑘,𝑟), we use subsets of entries of the auxiliary matrices 𝐴, 𝐵, and 𝐶 defined in the previous section. Since only the sign of Δ(𝑝,𝑘,𝑟) matters, the factor of 4 can be ignored in (13), and thus Δ^′(𝑝,𝑘,𝑟)=Δ(𝑝,𝑘,𝑟)/4 instead of Δ(𝑝,𝑘,𝑟) can be considered. From (13) we obtain Δ^′𝑝,𝑘,𝑟=𝑎_{𝑝(𝑘)𝑟}𝑎_{𝑝(𝑘)𝑟}+2𝑎_𝑝(𝑘),1−𝑆_𝑝(𝑘)+𝑎_{𝑝(𝑟)𝑘}𝑎_{𝑝(𝑟)𝑘}−2𝑎_𝑝(𝑟),1+𝑆_𝑝(𝑟)+𝑏_{𝑝(𝑟),𝑘+1}−𝑏_{𝑝(𝑘),𝑟−1}−2𝑐_{𝑝(𝑘)𝑝(𝑟)}. (19)

We also take advantage of the case where Δ^′(𝑝,𝑘,𝑟)=0. We use the lexicographic rule for breaking ties. The test outputs a “Fail” verdict if the equality Δ^′(𝑝,𝑘,𝑟)=0 is found to hold true for some 𝑘∈{1,…,𝑟−1} such that 𝑝(𝑘)>𝑝(𝑟). Indeed, 𝑝=(𝑝(1),…,𝑝(𝑘−1),𝑝(𝑘),𝑝(𝑘+1),…,𝑝(𝑟)) can be discarded because both 𝑝 and partial solution 𝑝^′=(𝑝(1),…,𝑝(𝑘−1),𝑝(𝑟),𝑝(𝑘+1),…,𝑝(𝑘)) are equally good, but 𝑝 is lexicographically larger than 𝑝^′.

We note that the entries of the matrices 𝐴, 𝐵, and 𝐶 used in (19) can be easily obtained from a subset of entries of 𝐴, 𝐵, and 𝐶 computed earlier in the search process. Consider, for example, 𝑐_{𝑝(𝑘)𝑝(𝑟)} in (19). It can be assumed that 𝑘∈{1,…,𝑟−2} because if 𝑘=𝑟−1, then 𝑐_{𝑝(𝑘)𝑝(𝑟)}=0. For simplicity, let 𝑝(𝑟)=𝑖 and 𝑝(𝑟−1)=𝑗. To avoid ambiguity, let us denote by 𝐶^′ the matrix 𝐶 that is obtained during processing of the partial solution (𝑝(1),…,𝑝(𝑟−2)). The entry 𝑐^′_{𝑝(𝑘)𝑖} of 𝐶^′ is calculated by (temporarily) placing the object 𝑖 in position 𝑟−1. Now, referring to the definition of 𝐶, it is easy to see that 𝑐_{𝑝(𝑘)𝑖}=𝑐^′_{𝑝(𝑘)𝑖}+𝑑_{𝑝(𝑘)𝑗}𝑑_𝑖𝑗. Similarly, using formulas of types (15) and (17), we can obtain the required entries of the matrices 𝐴 and 𝐵. Thus (19) can be computed in constant time. The complexity of computing Δ^′(𝑝,𝑘,𝑟) for all 𝑘 is 𝑂(𝑛). The dominance test used in [13] relies on (11) and has complexity 𝑂(𝑛²). 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 from 𝑝 by relocating the object 𝑝(𝑟) from position 𝑟 to position 𝑘∈{1,…,𝑟−2} (here position 𝑘=𝑟−1 is excluded since in this case relocating would become a pairwise interchange of objects). Thus, the test compares 𝑝 against partial solutions 𝑝_𝑘=(𝑝(1),…,𝑝(𝑘−1),𝑝(𝑟),𝑝(𝑘),…,𝑝(𝑟−1)), 𝑘=1,…,𝑟−2. Now suppose that 𝑝 and 𝑝_𝑘 are extended to 𝑛-permutations 𝑝∈Π and, respectively, 𝑝_𝑘∈Π such that, for each 𝑙∈{𝑟+1,…,𝑛}, the object in position 𝑙 of 𝑝 coincides with that in position 𝑙 of 𝑝_𝑘. Let us define 𝛿_𝑟𝑘 to be the difference between 𝐹(𝑝_𝑘) and 𝐹(𝑝). This difference can be computed as follows: 𝛿_𝑟𝑘=4^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}−2𝑓_𝑝(𝑟)+𝑆_𝑝(𝑟)+4^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}𝑑_{𝑝(𝑟)𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙). (20) 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 𝐴 and 𝐵.

Proposition 5. For 𝑝∈Π, 𝑘∈{1,…,𝑛−1}, and 𝑟∈{𝑘+1,…,𝑛}, 𝛿_𝑟𝑘=4𝑎_{𝑝(𝑟)𝑘}𝑎_{𝑝(𝑟)𝑘}−2𝑎_𝑝(𝑟),1+𝑆_𝑝(𝑟)+4𝑏_{𝑝(𝑟)𝑘}. (21)

The partial solution 𝑝 is dominated by 𝑝_𝑘 if 𝛿_𝑟𝑘>0. Upon establishing this fact for some 𝑘∈{1,…,𝑟−2}, the test immediately stops with the rejection of 𝑝. Of course, like in the case of (13), a factor of 4 in (21) can be dropped. The computational overhead of including the condition 𝛿_𝑟𝑘>0 in the test is minimal because all the quantities appearing on the right-hand side of (21), except 𝑏_𝑝(𝑟),1, are also used in the calculation of Δ^′(𝑝,𝑘,𝑟) (see (19)). Since the loop over 𝑘 needs to be executed the time complexity of this part of the test is 𝑂(𝑛). Certainly, evaluation of relocating the last added object in 𝑝 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 𝑈^′); (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 𝐹 at a permutation 𝑝∈Π and its reverse 𝑝 defined by 𝑝(𝑖)=𝑝(𝑛−𝑖+1), 𝑖=1,…,𝑛, are the same. Therefore, it is very natural to demand that some fixed object in 𝑊, say 𝑤, being among those objects that are assigned to the first half of the positions in 𝑝. More formally, the condition is 𝑖⩽⌈𝑛/2⌉ where 𝑖 is such that 𝑝(𝑖)=𝑤. If this condition is violated for 𝑝=(𝑝(1),…,𝑝(⌈𝑛/2⌉)), then the symmetry test fails, and 𝑝 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 Δ^′(𝑝,𝑘,𝑟)=0 and 𝑝(𝑘)>𝑝(𝑟) (see the paragraph following (19)), this test in our implementation of the algorithm gives a “Fail” verdict only if either 𝑟⩽⌈𝑛/2⌉ or 𝑝(𝑘)≠𝑤.

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 permutations 𝑝=(𝑝(1),…,𝑝(𝑟)) of length 𝑟 that is small or modest compared to 𝑛. Therefore, it is reasonable to compute the upper bound 𝑈^′ (described in Section 2) only when 𝑟⩾𝜆(𝐷), where 𝜆(𝐷) is a threshold value depending on the UDSP instance matrix 𝐷 and thus also on 𝑛. The problem we face is how to determine an appropriate value for 𝜆(𝐷). 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 𝑝₁,…,𝑝_𝜏, each of length ⌊𝑛/2⌋. For each permutation in this sample, the procedure computes the upper bound 𝑈^′. Let 𝑈^′₁,…,𝑈^′_𝜏 be the values obtained. We are interested in the normalized differences between these values and 𝐹(𝑝_heur), where 𝑝_heur is the solution delivered by the ITS algorithm. So the procedure calculates the sum 𝑍=∑^𝜏_𝑖=1100(𝑈^′_𝑖−𝐹(𝑝_heur))/𝐹(𝑝_heur) and the average 𝑅=𝑍/𝜏 of such differences. In the second step, using 𝑅, it determines the threshold value. This is done in the following way: if 𝑅<𝑅₁, then 𝜆(𝐷)=⌊𝑛/3⌋; if 𝑅₁⩽𝑅<𝑅₂, then 𝜆(𝐷)=⌊𝑛/2⌋; if 𝑅>𝑅₂, then 𝜆(𝐷)=⌊3𝑛/4⌋. The described procedure is executed at the initialization step of the algorithm. Its parameters are 𝑅₁, 𝑅₂, 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 by 𝑉. The main procedure goes as follows.

B&B  (1) Apply the ITS algorithm to get an initial solution 𝑝^∗ to a given instance of the UDSP.   Set 𝐹^∗∶=𝐹(𝑝^∗).  (2) Calculate 𝜆(𝐷).  (3)Initialize the search tree by creating the root node with the empty set 𝑉 attached to it.  Set 𝐿∶=0 (𝐿 denotes the tree level the currently considered node belongs to).  (4) Set 𝑟∶=𝐿.  (5) Call SelectObject(𝑟). It returns 𝐿 and, if 𝐿=𝑟, additionally some unassigned object (let it be denoted by 𝑣). If 𝐿<0, then go to (7). Otherwise, check whether 𝐿<𝑟. If so, then return to (4) (perform a backtracking step). If not (in which case 𝐿=𝑟), then proceed to (6).  (6) Create a new node with empty 𝑉. Set 𝑝(𝑟+1)∶=𝑣. Increment 𝐿 by 1, and go to (4).  (7)Calculate the coordinates for the objects by applying formula (4) with respect to the permutation 𝑝^∗ 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 𝑝^∗. If later a better solution is found, 𝑝^∗ is updated. This is done within the SelectObject procedure. In the final step, B&B uses 𝑝^∗ 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 𝑉. The partial solution corresponding to the current node is denoted by 𝑝. The objects of 𝑉 are used to extend 𝑝. Upon creation of a node (Step (3) for root and Step (6) for other nodes), the set 𝑉 is empty. In SelectObject, 𝑉 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 𝐿=𝑟, then the action is “branch,” whereas if 𝐿<𝑟, the action is “backtrack.” Step (4) helps to discriminate between these two cases. If 𝐿=𝑟, then a new node in the next level of the search tree is generated. At that point, the partial permutation 𝑝 is extended by assigning the selected object 𝑣 to the position 𝑟+1 in 𝑝 (see Step (6)). If, however, 𝐿<𝑟 and 𝑟 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(𝑟) (1)Check whether the set 𝑉 attached to the current node is empty. If so, then go to (2). If not, then consider its subset 𝑉^′={𝑖∈𝑉∣𝑖isunmarkedand𝑢_𝑖(𝑝)>𝐹^∗}. If 𝑉^′ is nonempty, then mark an arbitrary object 𝑣∈𝑉^′, and return with it as well as with 𝐿=𝑟. Otherwise go to (6). (2)If 𝑟>0, then go to (3). Otherwise, set 𝑉∶={1,…,𝑛}, 𝑢_𝑖(𝑝)∶=∞, 𝑖=1,…,𝑛, and go to (4). (3)If 𝑟=𝑛−1, then go to (5). Otherwise, for each 𝑖∈𝑊(𝑝), perform the following operations. Temporarily set 𝑝(𝑟+1)∶=𝑖. Apply the symmetry test, the dominance test, and, if 𝑟⩾𝜆(𝐷), also the bound test to the resulting partial solution 𝑝^′=(𝑝(1),…,𝑝(𝑟+1)). If it passes all the tests, then do the following. Include 𝑖 into 𝑉. If 𝑟⩾𝜆(𝐷), assign to 𝑢_𝑖(𝑝) the value of the upper bound 𝑈^′ computed for the subproblem defined by the partial solution 𝑝^′. If, however, 𝑟<𝜆(𝐷), then set 𝑢_𝑖(𝑝)∶=∞. If after processing all objects in 𝑊(𝑝), the set 𝑉 remains empty, then go to (6). Otherwise proceed to (4). (4)Mark an arbitrary object 𝑣∈𝑉, and return with it as well as with 𝐿=𝑟. (5)Let 𝑊(𝑝)={𝑗}. Calculate the value of the solution 𝑝=(𝑝(1),…,𝑝(𝑛−1),𝑝(𝑛)=𝑗). If 𝐹(𝑝)>𝐹^∗, then set 𝑝^∗∶=𝑝 and 𝐹^∗∶=𝐹(𝑝). (6)If 𝑟>0, then drop the 𝑟th element of 𝑝. Return with 𝐿∶=𝑟−1.

In the previous description, 𝑊(𝑝) denotes the set of unassigned objects, as before. The objects included in 𝑉⊆𝑊(𝑝) are considered as candidates to be assigned to the leftmost available position in 𝑝. At the root node, 𝑉=𝑊(𝑝). In all other cases, except when |𝑊(𝑝)|=1, an object becomes a member of the set 𝑉 only if it passes two or, for longer 𝑝, three tests. The third of them, the upper bound test, amounts to check if 𝑈^′>𝐹^∗. 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 𝑖∈𝑉. Such a strategy is reasonable because during the branch-and-bound process the value of 𝐹^∗ may increase, and the previously satisfied constraint 𝑈^′>𝐹^∗ may become violated for some 𝑖∈𝑉. To be able to identify such objects, the algorithm stores the computed upper bound 𝑈^′ as 𝑢_𝑖(𝑝) 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 𝑉 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, 𝑟, one of the three cases is considered. If the current node is the root, that is, 𝑟=0, then Step (2) is executed. Since in this case the dominance test is not applicable, the algorithm uses the set 𝑉 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 𝑟∈[1,…,𝑛−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 𝑝 by assigning an object from 𝑊(𝑝) to the (𝑟+1) th position in 𝑝. In order not to lose at least one optimal solution it is sufficient to keep only those extensions of 𝑝 for which all the tests were successful. Such partial solutions are memorized by the set 𝑉. 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 𝑉 is empty, then the current node can be discarded, and the search continued from its parent node. Finally, if 𝑟=𝑛−1, then Step (5) is executed in which the permutation 𝑝 is completed by placing the remaining unassigned object in the last position of 𝑝. Its value is compared with that of the best solution found so far. The node corresponding to 𝑝 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 𝑉 is scanned. The previously selected objects in 𝑉 have been marked. They are ignored. For each remaining object 𝑖, it is checked whether 𝑢_𝑖(𝑝)>𝐹^∗. If so, then 𝑖 is included in the set 𝑉^′. If the resulting set 𝑉^′ 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 𝑉^′ (or 𝑉 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 𝑀=(𝑚_𝑖𝑗). Its entries are integers from the interval [0,100]. Like in the study of Brusco and Stahl [13], the corresponding symmetric dissimilarity matrix 𝐷=(𝑑_𝑖𝑗) was obtained by using the formula 𝑑_𝑖𝑗=200−𝑚_𝑖𝑗−𝑚_𝑗𝑖 for all off-diagonal entries and setting the diagonal of 𝐷 to zero. We considered two UDSP instances constructed using Morse code confusion data—one defined by the 26×26 submatrix of 𝐷 with rows and columns labeled by the letters and another defined by the full matrix 𝐷.

The second empirical matrix was taken from Groenen and Heiser [29]. We denote it by 𝐺=(𝑔_𝑖𝑗). Its entry 𝑔_𝑖𝑗 gives the number of citations in journal 𝑖 to journal 𝑗. Thus the rows of 𝐺 correspond to citing journals, and the columns correspond to cited journals. The dissimilarity matrix 𝐷 was derived from 𝐺 in two steps, like in [13]. First, a symmetric similarity matrix, 𝐻=(ℎ_𝑖𝑗), was calculated by using the formulas ℎ_𝑖𝑗=ℎ_𝑗𝑖=(𝑔_𝑖𝑗/∑^𝑛_{𝑙=1,l≠𝑖}𝑔_𝑖𝑙)+(𝑔_𝑗𝑖/∑^𝑛_{𝑙=1,𝑙≠𝑗}𝑔_𝑗𝑙), 𝑖,𝑗=1,…,𝑛, 𝑖<𝑗, and ℎ_𝑖𝑖=0, 𝑖=1,…,𝑛. Then 𝐷 was constructed with entries 𝑑_𝑖𝑗=⌊100(max_𝑙,𝑘ℎ_𝑙𝑘−ℎ_𝑖𝑗)⌋, 𝑖,𝑗=1,…,𝑛, 𝑖≠𝑗, and 𝑑_𝑖𝑖=0, 𝑖=1,…,𝑛.

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, 𝐷_food, was constructed based on data provided by Ross and Murphy [31]. For our computational tests, we considered submatrices of 𝐷_food defined by the first 𝑛∈{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, 𝑅₁=0, and 𝑅₂=5. In the symmetry test, 𝑤 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 Table 1. The dimension of the matrix 𝐷 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 𝐹^∗ and, in parentheses, two related metrics. The first of them is the raw value of the stress function Φ(𝑥). The second number in parentheses is the normalized value of the stress function that is calculated as Φ(𝑥)/∑_𝑖<𝑗𝑑_𝑖𝑗. 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.

Table 2 

Performance on empirical dissimilarity matrices.

As 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 around 40% 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 Section 6.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.

Table 3 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 𝑖∈{22,…,35}, the CPU time taken by B&B to solve the problem with the 𝑖×𝑖 food-item submatrix is between 196 (for 𝑖=30) and 282 (for 𝑖=27) percent of that needed to solve the problem with the (𝑖−1)×(𝑖−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.

In Table 4, 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 𝑝.

Table 4 

Optimal solutions for the full Morse code dissimilarity matrix and the 35 × 35 food-item dissimilarity matrix.

In Table 5, 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 100𝑁_sym/𝑁, 100𝑁_dom/𝑁, and 100𝑁_UB/𝑁, where 𝑁_sym (respectively, 𝑁_dom and 𝑁_UB) 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 𝑁=𝑁_sym+𝑁_dom+𝑁_UB. The last three columns give the same set of statistics for the BB3 algorithm. Inspection of Table 5 reveals that 𝑁_dom is much larger than both 𝑁_sym and 𝑁_UB. It can also be seen that 𝑁_dom values are consistently higher with BB3 than when using our approach. Comparing the other two tests, we observe that 𝑁_sym is larger than 𝑁_UB 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.

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 to 30×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 𝐷_food. 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 𝐷_food 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.

Appendix

Derivation of (20)

We will deduce an expression for 𝛿_𝑟𝑘, 1⩽𝑘<𝑟⩽𝑛, assuming that 𝑝 is a solution to the whole problem, that is, 𝑝∈Π. Relocating the object 𝑝(𝑟) from position 𝑟 to position 𝑘 gives the permutation 𝑝^′=(𝑝(1),…,𝑝(𝑘−1),𝑝(𝑟),𝑝(𝑘),…,𝑝(𝑛)). By using the definition of 𝐹 we can write 𝛿_𝑟𝑘=𝐹𝑝^′𝑝=−𝐹_{𝑘−1𝑙=1}𝑑_{𝑝(𝑟)𝑝(𝑙)}−^𝑛_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}²−^𝑟−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}−^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}²+^𝑟−1_𝑙=𝑘^𝑙−1_𝑗=1𝑑_{𝑝(𝑙)𝑝(𝑗)}−^𝑛_𝑗=𝑙+1𝑑_{𝑝(𝑙)𝑝(𝑗)}+2𝑑_{𝑝(𝑟)𝑝(𝑙)}²−^𝑟−1_𝑙=𝑘^𝑙−1_𝑗=1𝑑_{𝑝(𝑙)𝑝(𝑗)}−^𝑛_𝑗=𝑙+1𝑑_{𝑝(𝑙)𝑝(𝑗)}². (A.1) Let us denote the first two terms in (A.1) by 𝛿^′ and the remaining two terms by 𝛿^′′=∑^𝑟−1_𝑙=𝑘𝛿_𝑙^′′. It is easy to see that 𝛿^′ stems from moving the object 𝑝(𝑟) to position 𝑘 and 𝛿^′′ stems from shifting objects 𝑝(𝑙), 𝑙=𝑘,…,𝑟−1, by one position to the right. To get a simpler expression for 𝛿^′, we perform the following manipulations: 𝛿^′=^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}²−2^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}²−^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}²+2^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}−^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}². (A.2)

Continuing, we get 𝛿^′=^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}²−2^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}−2^𝑘−1_𝑙=1𝑑^𝑛_{𝑝(𝑟)𝑝(𝑙)𝑙=𝑟+1}𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}²+2^𝑟−1_𝑙=𝑘𝑑^𝑛_{𝑝(𝑟)𝑝(𝑙)𝑙=𝑟+1}𝑑_{𝑝(𝑟)𝑝(𝑙)}+^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}²−^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}²−2^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}−^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}²+2^𝑘−1_𝑙=1𝑑^𝑛_{𝑝(𝑟)𝑝(𝑙)𝑙=𝑟+1}𝑑_{𝑝(𝑟)𝑝(𝑙)}+2^𝑟−1_𝑙=𝑘𝑑^𝑛_{𝑝(𝑟)𝑝(𝑙)𝑙=𝑟+1}𝑑_{𝑝(𝑟)𝑝(𝑙)}−^𝑛_𝑙=𝑟+1𝑑_{𝑝(𝑟)𝑝(𝑙)}². (A.3)

Finally, 𝛿^′=−4^𝑘−1_𝑙=1𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}+4^𝑟−1_𝑙=𝑘𝑑^𝑛_{𝑝(𝑟)𝑝(𝑙)𝑙=𝑟+1}𝑑_{𝑝(𝑟)𝑝(𝑙)}=−4𝑓_𝑝(𝑟)−^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}+4^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}𝑆_𝑝(𝑟)−𝑓_𝑝(𝑟)=4^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑟−1_𝑙=𝑘𝑑_{𝑝(𝑟)𝑝(𝑙)}−2𝑓_𝑝(𝑟)+𝑆_𝑝(𝑟). (A.4)

Similarly, for 𝛿_𝑙^′′ we have 𝛿_𝑙^′′=^𝑙−1_𝑗=1𝑑_{𝑝(𝑙)𝑝(𝑗)}−^𝑛_𝑗=𝑙+1𝑑_{𝑝(𝑙)𝑝(𝑗)}²+4𝑑_{𝑝(𝑟)𝑝(𝑙)}^𝑙−1_𝑗=1𝑑_{𝑝(𝑙)𝑝(𝑗)}−^𝑛_𝑗=𝑙+1𝑑_{𝑝(𝑙)𝑝(𝑗)}+4𝑑²_{𝑝(𝑟)𝑝(𝑙)}−^𝑙−1_𝑗=1𝑑_{𝑝(𝑙)𝑝(𝑗)}−^𝑛_𝑗=𝑙+1𝑑_{𝑝(𝑙)𝑝(𝑗)}²=4𝑑_{𝑝(𝑟)𝑝(𝑙)}𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙)+𝑓_𝑝(𝑙)+4𝑑²_{𝑝(𝑟)𝑝(𝑙)}=4𝑑_{𝑝(𝑟)𝑝(𝑙)}𝑑_{𝑝(𝑟)𝑝(𝑙)}+2𝑓_𝑝(𝑙)−𝑆_𝑝(𝑙). (A.5) By summing up (A.4) and (A.5) for 𝑙=𝑘,…,𝑟−1 we obtain (20).