Visualization of complex dynamic datasets by means of mathematical optimization

Omega 86 (2019) 125–136 Contents lists available at ScienceDirect Omega journal homepage: www.elsevier.com/locate/omega Visualization of complex dynamic datasets by means of mathematical optimization Emilio Carrizosa a , Vanesa Guerrero b , ∗, Dolores Romero Morales c a Instituto de Matemáticas de la Universidad de Sevilla (IMUS), Seville, Spain b Department of Statistics, Universidad Carlos III de Madrid, Getafe, Spain c Copenhagen Business School, Frederiksberg, Denmark a r t i c l e i n f o Article history: Received 27 April 2017 Accepted 23 July 2018 Available online 26 July 2018 Keywords: Visualization Dynamic magnitude Multidimensional scaling Difference of convex optimization a b s t r a c t In this paper we propose an optimization model and a solution approach to visualize datasets which are made up of individuals observed along different time periods. These individuals have attached a timedependent magnitude and a dissimilarity measure, which may vary over time. Difference of convex optimization techniques, namely, the so-called Difference of Convex Algorithm, and nonconvex quadratic binary optimization techniques are used to heuristically solve the optimization model and develop this visualization framework. This way, the so-called Dynamic Visualization Map is obtained, in which the individuals are represented by geometric objects chosen from a catalogue. A Dynamic Visualization Map faithfully represents the dynamic magnitude by means of the areas of the objects, while it trades off three different goodness of fit criteria, namely the correct match of the dissimilarities between the individuals and the distances between the objects representing them, the spreading of such objects in the visual region, and the preservation of the mental map by ensuring smooth transitions along snapshots. Our procedure is successfully tested on dynamic geographic and linguistic datasets. © 2018 Elsevier Ltd. All rights reserved. 1. Introduction Revealing and interpreting the underlying structures in complex datasets is a challenge that analysts have to face. Such demanding tasks arise in diverse contexts as health care, [8] , risk management, [3,59] or text mining, [4,60] . Information Visualization arises as a discipline to give answers to such demanding tasks by designing suitable visualization frameworks, which bring hidden patterns to light and enhance interpretability [7,21,45] . Mathematical Optimization has broadly contributed to its development in terms of modeling and algorithmic approaches [8,16–19,29,30,51] , but many problems still deserve further attention. This paper contributes to the literature concerning applications of Mathematical Optimization to Information Visualization. In particular, we focus on one of the most relevant data types according to [57] , namely, dynamic multidimensional data, which consist of observations whose multiple attributes change over time [22,26,50] . This work deals with datasets which consist of N individuals, V = { v 1 , . . . , v N } , for which a magnitude or weight has been historically observed during T time periods, and a measure of proximity between individ- ∗Corresponding author. E-mail addresses: [email protected] (E. Carrizosa), v[email protected] (V. Guerrero), [email protected] (D. Romero Morales). uals, given through a time-dependent dissimilarity matrix, is also given [36,49] . Visualizing individuals with attached dissimilarities has been historically done by Multidimensional Scaling (MDS) [12] , whose aim is to represent the individuals in the dataset as points in a low dimensional space (usually R 2 ) in such a way that the distances between the points approximate the given dissimilarities. A straightforward approach to visualize dynamic multidimensional datasets, involving dissimilarities observed along T time periods, would consist of executing T independent MDS, yielding T snapshots, one per period. Nevertheless, this approach might yield difficult-to-interpret visualizations [34,61] , especially when the dissimilarities change abruptly in consecutive periods or, since MDS results are invariant under rotations and reflections, and thus the snapshots may turn upside-down. In order to illustrate this statement, let us consider the dataset consisting of N = 13 stock market indices, studied along T = 200 time periods. Dissimilarities are measured through the correlation between stock market indices as studied in [34] . The plots for time periods t = 30 and t = 31 , obtained with the isoMDS() function in R, [54] , are given in Fig. 1 . As the reader can observe, a visual effort is required to read from snapshot t = 30 to t = 31 , since stock market indices appear rotated and further from each other. This example calls for the construction of visualization frameworks which trade off a faithfull https://doi.org/10.1016/j.omega.2018.07.008 0305-0483/© 2018 Elsevier Ltd. All rights reserved. 126 E. Carrizosa et al. / Omega 86 (2019) 125–136 Fig. 1. Illustration of the lack of preservation of mental map in two consecutive snapshots when MDS is run independently in the stock market indices dataset in [34] . representation of dissimilarities and a preservation of the mental map [48] , i.e., the transitions in the layouts in two consecutive time periods should be smooth, in the sense that the objects representing the individuals do not suffer big displacements in the shift from one period to the next one. On top of the challenge posed by the fact that data are time-varying, we also want to visualize the magnitude attached at each time period to each individual. To do this, a first approach might consist of executing an MDS, and then replacing points by geometric objects objects, say, discs or rectangles, centered at the MDS points, having areas proportional to the magnitude values [38] . Now, the sizes of the objects play a key role in the computation of distances between them. Therefore, the changes in the perception of distances, induced by the sizes of these (a posteriori) depicted objects, might yield misleading conclusions about the dissimilarities between the individuals. In [18] , the authors developed a generalization of MDS, based on a Mathematical Optimization model, to visualize static data. That framework simultaneously incorporates the information about dissimilarities and a magnitude, which were optimally rescaled within user’s given bounds, by (a priori) deciding which geometric object (disc, rectangle, etc.) represents each individual. In that work, dissimilarities were reproduced as distances between the geometric objects representing each individual. See e.g. Fig. 2 for an example in which N = 27 Danish words, depicted as rectangles, whose areas represent the importance of the word in the Danish news in 1995, and distances between rectangles reflect the semantic relatedness between the words they represent (the more related, the smaller their dissimilarity is). However, the visualization framework in Carrizosa et al. [18] suffers from the same lack of smoothness than MDS when dealing with time-varying data, namely the location of the objects can abruptly change in consecutive periods when each snapshot is built independently. This work presents a novel application of mixed integer nonconvex optimization to Information Visualization. We propose a new Mathematical Optimization model and a heuristic algorithm as solution approach to visualize dynamic datasets involving a magnitude and dissimilarities. We design a one-stage procedure, which involves a non-trivial generalization of the approach presented in [18] . This new visualization framework simultaneously Fig. 2. Visualization of Danish words by means of different rectangles, scaled according to their frequencies, and distances representing words semantic relatedness, [18] . builds a collection of T snapshots, each representing a time period, in which the individuals under consideration are depicted as geometric objects located in a visualization region , whose areas represent the magnitude or weights and the distance between the objects depict the dissimilarities. The novelty of our model with respect to the one in [18] is twofold. First, the preservation of the mental map is incorporated into the optimization model by pursuing smooth transitions between two consecutive snapshots. Second, we assume that each individual has attached a catalogue of E. Carrizosa et al. / Omega 86 (2019) 125–136 127 Fig. 3. Importance of the choice of the geometric objects depicting each individual. candidate objects to be depicted, and the choice of the object becomes a decision of the model, making it more flexible. For instance, consider a set of N = 3 individuals V = { v 1 , v 2 , v 3 } which could be depicted by means of two different rectangles: one whose basis has two units and one unit of height, and its 90 °rotation. Two possible representations are depicted in Fig. 3 . Observe that, whereas the rectangles appear collapsed in Fig. 3 (left), a different choice from the catalogue makes the visualization clearer, Fig. 3 (right), in which different rectangles have been chosen for individuals v 2 and v 3 . Although there have been some studies, see e.g. [27] , suggesting that people place more importance on objects that are horizontally oriented than when appearing as vertically oriented, we consider that tuning the angle used to represent each individual is helpful to visualize them in a more meaningful and clear manner. Both types of information, namely dissimilarities and a magnitude, are reproduced regardless of the orientation of the objects. The optimization model proposed to visualize the complex dynamic dataset described above can be formulated as a nonconvex Mixed Integer Nonlinear Problem (MINLP), which is handled through a heuristic algorithm, wich is based on an alternating strategy. On one hand, Difference of Convex (DC) optimization techniques are used to optimize the continuous variables involved in the problem [40,41] . Our proposal is a follow-up work to [18] , where the scaling variables to proportionally depict the information have become constants. This choice simplifies the expressions of the DC functions involved in the optimization problem, and thus, reduces the number of parameters to be chosen by the user. On the other hand, nonconvex quadratic binary optimization tools are exploited to handle the combinatorial stage of the alternating procedure [14] . Our approach is clearly different from existing techniques in the literature, mostly ad-hoc multi-stage procedures, which often depend on user’s manual tuning or exploit the nature of the data, and thus cannot be used for arbitrary datasets. Some examples are found in Graph Drawing [5,6,44,46,59,61] , geographical applications [2,47] and text visualizations, which have word clouds as their main instrument [20,25,37] . Word clouds are usually generated using ad-hoc multi-stage procedures, aiming to avoid empty spaces and overlap between words, which are written horizontally and/or vertically and scaled according to a magnitude. Our visualization framework is based on an optimization model and a solution approach. Therefore, whereas traditional approaches to build word clouds are not versatile enough to include further information, our model is able to build visuals for a linguistic dataset which incorporate additional features, either as new criteria or constraints, and different encodings for the words, which can enrich the visuals according to the user’s needs. In particular, we consider the representation of dissimilarities as well as their temporal evolution, jointly with the visualization of a dynamic magnitude. See also [1,33] for further applications. The remainder of the paper is organized as follows. Section 2 is devoted to describe the model to visualize the dynamic complex dataset under study. In Section 3 , we present a solution approach based on DC optimization tools and nonconvex quadratic binary optimization. Some computational tests, involving data of different nature, are included in Section 4 . Section 5 contains some conclusions and future lines of research. Finally, the Appendix closes the paper with some technical details. 2. Dynamic Visualization Map: the model In what follows, the problem of visualizing a dynamic dataset by means of geometric objects is formally stated and written as a mathematical optimization program. Let ⊆R 2 be a visualization region, which acts as the computer screen. Let us consider a set of N individuals, V = { v 1 , . . . , v N } , which have been observed over a time horizon of T time periods. For each t = 1 , . . . , T , let V ( t ) ⊆V be the subset of individuals to be represented in time period t , and let | V ( t )| be its cardinality. The elements in V ( t ) have attached a magnitude ω(t) = ω i,t i ∈ V (t) ∈ R | V (t) | + and a dissimilarity measure δ(t) = δij,t i,j∈ V (t) ∈ R | V (t) | ×| V (t) | + . Let B i = B 1 i , . . . , B s i i be a catalogue of geometric objects, called in what follows reference objects , which are assumed to be closed convex sets, centered in the origin of the coordinate system in R 2 . The elements in B i are the candidates to represent each v i ∈ V , although just one of them is chosen for the whole time horizon. Let τ, r p i,t ∈ R + be real positive numbers, which scale the area of the reference objects in B i . The scaling of the reference objects is made in such a way that the area of r p i,t B p i is equal to ω i, t . Besides this scaling, τis a positive parameter to be chosen by the user, which rescales all the objects in all periods in order to make sure they fit into . The dynamic dataset described above is visualized by means of a collection of T snapshots, each containing the individuals in V ( t ) depicted as geometric objects, such as discs or rectangles. In order to properly visualize the magnitudes and dissimilarities attached to the data, five conditions are considered: (C1) Each individual v i ∈ V is represented by means of the same reference object, chosen from the corresponding catalogue B i , throughout the whole time horizon. (C2) In each time period t , the area of the geometric object used to represent each individual in V ( t ) is proportional to its magnitude ω( t ). (C3) In each time period t , the distance between the geometric objects representing the individuals in V ( t ) resemble the dissimilarities δ( t ). (C4) In each time period, the geometric objects are spread over the visualization region . (C5) The transitions between two consecutive snapshots are smooth. In what follows, we introduce a mathematical optimization model which considers (C1) and (C2) as hard conditions, whereas the violation of conditions (C3)–(C5) is minimized. A visualization framework satisfying these conditions is called in what follows a Dynamic Visualization Map . Let x = (x p i ) i =1 , ... ,N,p=1 , ... ,s i be decision variables defined as x p i = 1 if individual v i is represented by B p i ∈ B i 0 otherwise , and let c i,t ∈ R 2 be continuous variables, which translate the reference objects in B i and determine their positions. In other words, if the reference object B p i ∈ B i is chosen to represent individual v i , then v i will be represented at time period t by the convex 128 E. Carrizosa et al. / Omega 86 (2019) 125–136 body c i,t + τr p i,t B p i . Therefore, constructing a Dynamic Visualization Map is stated as a Mixed Integer Nonlinear Optimization Problem (MINLP), whose aim is to find the choice x of reference objects and the values of the translation vectors c 1 , 1 , ... , c N,T to obtain a good fit in criteria (C3)–(C5) modeled through an objective function F to be detailed later. Hence, one has to solve a problem of the form: min c 1 , 1 , ... , c N,T , x F (c 1 , 1 , ... , c N,T , x ) s . t .  p=1 , ... ,s i x p i = 1 , i = 1 , . . . , N, (DyV iMap) c i,t + τr p i,t x p i B p i ⊆, i = 1 , . . . , N, p = 1 , . . . , s i , t = 1 , . . . , T , c i,t ∈ R 2 , i = 1 , . . . , N;t = 1 , . . . , T , x p i ∈ { 0 , 1 } , i = 1 , . . . , N, p = 1 , . . . , s i . The first constraint in ( DyViMap ) ensures condition (C1) is satisfied, namely, for each individual, only one reference object is chosen, among the candidates in its catalogue, to represent the individual in all time periods. The second constraint when x p i = 1 ensures that the whole geometric object representing v i , constructed by translating and scaling the reference object B p i , must fit into . In particular, this means that c i, t ∈ . The second constraint when x p i = 0 does not add any new information, given the observation we just made on c i, t . Observe that setting r p i,t so that the area of r p i,t B p i is equal to ω i, t , condition (C2) is satisfied, independently of the choice of reference object. Finally, the type of the variables is modeled through the third and fourth constraints. Observe that there are 2 ×N ×T continuous variables and N  i =1 , ... ,N s i binary variables. In order to illustrate the role of the second constraint in ( DyViMap ), we present an example. Let be the unit square, namely = [0 , 1] 2 . The catalogues for the individuals v i , i = 1 , ... , N, in V are defined as B i = B 1 i , B 2 i ,where B 1 i is a square with side of length  and B 2 i is a ball with radius of length ρ. Therefore, for each individual there are two binary variables: x 1 i and x 2 i . On one hand, the second constraint in ( DyViMap ) when p = 1 becomes 0 ≤c i,t ±τr 1 i,t x 1 i  2 ,  2 ≤1 , (1) which ensures that the scaled squared does lie inside when x 1 i = 1 , and c i, t does lie too when x 1 i = 0 . Observe that in case B 1 i were a rectangle, we would distinguish among its two different length sides yielding an analogous constraint. On the other hand, the second constraint in ( DyViMap ) when p = 2 becomes 0 ≤c i,t ±τr 1 i,t x 2 i ( ρ, ρ) ≤1 , (2) which ensures that the scaled disc does lie inside when x 2 i = 1 , and c i, t does lie too when x 2 i = 0 . The objective function in ( DyViMap ), F , models the violation of conditions (C3)–(C5) considering a weighted sum of three functions, F MDS , F spread and F smooth , through a vector λ=(λ1 , λ2 , λ3 ) , such that λk ≥0 and  3 k =1 λk = 1 , yielding F = λ1 F MDS + λ2 F spread + λ3 F smooth . The first term, F MDS , measures the discrepancy between the given dissimilarities and the distances between the objects (condition (C3)), by considering the STRESS expression in MDS [12,23] . The second one, F spread , quantifies the spread of the objects in the visualization region (condition (C4)). Finally, F smooth models the smoothness in the transition between two consecutive periods (condition (C5)). In order to measure the distance between two geometric objects in the t -th time period, namely the distance between c i,t + τr p i,t B p i and c j,t + τr q j,t B q j , we consider the infimum distance between two closed convex sets, which is a convex function, [35] . Let  · denote the Euclidean norm, then the infimum distance between c i,t + τr p i,t B p i and c j,t + τr q j,t B q j is defined as d τr p i,t B p i ;τr q j,t B q j : R 2 ×R 2 −→ R + (c i,t , c j,t )  −→ inf b i,t ∈B p i b j,t ∈B q j  c i,t + τr p i,t b i,t −c j,t + τr q j,t b j,t  , Then, the expressions of F MDS , F spread and F smooth are: F MDS (c 1 , 1 , . . . , c N,T , x ) = T  t=1  i,j∈ V (t)  p=1 , ... ,s i q =1 , ... ,s j  d τr p i,t B p i ;τr q j,t B q j (c i,t , c j,t ) −κδij,t  2 x p i x q j , F spread (c 1 , 1 , . . . , c N,T , x ) = − T  t=1  i,j∈ V (t)  p=1 , ... ,s i q =1 , ... ,s j d 2 τr p i,t B p i ;τr q j,t B q j (c i,t , c j,t ) x p i x q j F smooth (c 1 , 1 , . . . , c N,T ) = T −1  t=1  i =1 , ... ,N  c i,t −c i,t+1  2 . Note that F MDS and F spread are generalizations of those presented in [18] for the particular case in which the dataset is not dynamic, i.e., just one time period is considered. The expression of F MDS includes a positive parameter κ(to be chosen by the user), which scales the dissimilarities to make them comparable with the distance between objects measured by means of the infimum distance. The spread criterion, modeled through F spread , aims to separate the geometric objects representing the individuals as much as possible by means of the squared infimum distance between them. The preservation of the mental map is modeled through F smooth , which imposes that the locations of the objects, given by their translation vectors, do not suffer big changes from one time period to the next one, [61] , and by the fact that an individual is depicted by means of the same reference object throughout the whole time horizon, since x does not depend on the time period. The aim of the presented model is to obtain a trade-off between the criteria involved to enhance the interpretability of the dynamic complex data structure under consideration. 3. Dynamic Visualization Map: the algorithmic approach Section 2 states the problem of building a Dynamic Visualization Map as a Mixed Integer Nonlinear Problem (MINLP). This section is devoted to present a solution approach to solve ( DyViMap ), in which continuous and binary variables are optimized in an alternating procedure: the choice x of reference objects depicting the individuals, belonging to their corresponding catalogue, is optimized for translations c 1 , 1 , . . . , c N,T fixed, then the translation vectors c 1 , 1 , . . . , c N,T are optimized for x fixed, and the process is repeated until a stopping criterion is satisfied. This solution approach is a heuristic motivated from the theoretical work and the numerical results obtained in [18] . On one hand, observe that if the continuous variables c 1 , 1 , . . . , c N,T in ( DyViMap ) are fixed, the resulting problem is a nonconvex binary quadratic optimization problem with assignment constraints of the form min x  i,j∈ V  p=1 , ... ,s i q =1 , ... ,s j a pq ij (c) x p i x q j s . t .  s i p=1 x p i = 1 , i = 1 , . . . , N, (DyV iMap) c c i,t + τr p i,t x p i B p i ⊆, i = 1 , . . . , N, p = 1 , . . . , s i , t = 1 , . . . , T , x p i ∈ { 0 , 1 } , i = 1 , . . . , N, p = 1 , . . . , s i , E. Carrizosa et al. / Omega 86 (2019) 125–136 129 where the coefficients a pq ij (c) take the form a pq ij (c) =  t :i,j∈ V (t ) λ1  d τr p i,t B p i ;τr q j,t B q j (c i,t , c j,t ) −κδij,t  2 −λ2 d 2 τr p i,t B p i ;τr q j,t B q j (c i,t , c j,t )  Problem ( DyViMap ) c can thus be solved by standard MINLP Global Optimization solvers. However, when, at most, two reference objects are in the catalogue of each individual ( s i = 2 , for all i = 1 , . . . , N), the problem ( DyViMap ) c can be rewritten as an unconstrained convex quadratic 0–1 problem. In this particular case x 2 i = 1 −x 1 i for all i = 1 , . . . , N, and then, rising up the diagonal of the quadratic form in the objective until it is positive semidefinite, we can obtain an equivalent convex objective function, [9] . With this rising-up classical trick, the unconstrained convex quadratic 0– 1 reformulation holds and Mixed Integer Quadratic Programming solvers can be used instead. Observe than in this particular case, the number of binary variables is reduced to N . On the other hand, for fixed values of x , ( DyViMap ) becomes a nonlinear continuous optimization problem, which involves a difference of convex (DC) function to be minimized. Indeed, when x is fixed, F is DC since F MDS is DC (squared difference of a positive convex function and a positive real number), F smooth is convex, F spread is concave, and λk ≥0, k = 1 , 2 , 3 . It is worth noting that, since x is fixed, each individual is allocated to one reference object, and the model is analogous to the single-object case studied in [18] with the additional convex term F smooth in the objective. Thus, the Difference of Convex Algorithm (DCA) is a suitable tool to find good quality solutions which requires a DC decomposition of the objective function. The performance of the DCA strongly depends on the choice of the DC decomposition, [10,11,28] . In this work, as in [18,39,53] , we seek a DC decomposition of F , with fixed x , whose expression is formed by a quadratic separable convex function minus a convex function, as stated in Proposition 1 . Proposition 1. For a given x , function F can be expressed as a DC function, F = u −(u −F ) , where the quadratic separable convex function u is given by u = T  t=1  i,j∈ V (t) 2 max { λ1 −λ2 , 0 }  c i,t  2 +  c j,t  2 + λ1 κ2 δ2 ij,t  +2 λ3 N  i =1   c i, 1  2 +  c i,T  2 + 2  t=2 , ... ,T −1  c i,t  2  Proof. See Appendix.  Roughly speaking, DCA consists of an iterative process in which a sequence of convex programs are solved. At each iteration, the concave part is replaced by its affine majorization at a certain feasible point, and the resulting convex problem is then solved. For further details about the DCA and its application see [40–43,52] . Thanks to the DC decomposition of F given in Proposition 1 , for x fixed, one needs to solve N ×T convex quadratic problems with simple constraints: min c i,t M i,t  c i,t  2 −c  i,t γ¯ c i,t  s . t . c i,t + τr p i,t x p i B p i ⊆, (DyV iMap) DCA i,t c i,t ∈ R 2 , for scalars M i,t ∈ R + , which follow from the coefficients that multiply each term in the u part (after grouping terms) in Proposition 1 . Vectors γ¯ c i,t ∈ R 2 are subgradients of the function u −F evaluated in the locations obtained in the previous iterations of DCA, ¯ c = ¯ c 1 , 1 , . . . , ¯ c N,T . When has an amenable form, for instance a box or a disc, the optimal solution of problems (DyV iMap) DCA i,t , i = 1 , . . . , N, t = 1 , . . . , T , can be readily obtained by differentiating and equating the gradient to zero (considering correctly the constraints). In this case, obviously, running times will be strongly reduced since the convex optimization problems to be solved at each stage of DCA have a closed expression for their optimal value. The DCA scheme for solving ( DyViMap ) with fixed x is outlined in Algorithm 1 . Algorithm 1 DCA scheme for ( DyViMap ) with fixed x . Input: c ini = c ini 1 , 1 , . . . , c ini N,T , such that c ini i,t + τr p i,t x p i B p i ⊆, i = 1 , . . . , N, t = 1 , . . . , T . 1: ¯ c ← c ini 2: repeat 3: Compute γ¯ c i,t ∈ ∂ ( u −F ) ( ¯ c ) ; 4: Compute c = c 1 , 1 , . . . , c N,T as the solution of the problem (DyV iMap) DCA i,t , for all i = 1 , . . . , N, t = 1 , . . . , T ; 5: ¯ c ← c; 6: until stop condition is met. Output: ¯ c = ¯ c 1 , 1 , . . . , ¯ c N,T . Summarizing, solving the problem ( DyViMap ) by means of an alternating algorithm which optimizes its continuous and binary variables, respectively, requires the call to a DCA subroutine in the first case ( Algorithm 1 to optimize c 1 , 1 , . . . , c N,T ) and an integer nonconvex quadratic solver to optimize x . The alternating routine to solve ( DyViMap ) is given in Algorithm 2 . Algorithm 2 Alternating scheme for ( DyViMap ). Input: x ini = (x ini ) p i ∈ { 0 , 1 } S , where S =  N k =1 s k , and c ini = c ini 1 , 1 , . . . , c ini N,T , such that c ini i,t + τr p i,t (x ini ) p i B p i ⊆, i = 1 , . . . , N, p = 1 , . . . , s i , t = 1 , . . . , T . 1: ¯ c ← c ini ; 2: ¯ x ← x ini ; 3: repeat 4: ¯ c ← Algorithm 1 ( ¯ c ) ; 5: ¯ x ← solve (DyV iMap) c ; 6: until stop condition is met. Output: ¯ c = ¯ c 1 , 1 , . . . , ¯ c N,T , ¯ x = ¯ x p i , i = 1 , . . . , N, p = 1 , . . . , s i . 4. Computational experience The methodology proposed in Section 3 is illustrated in two datasets, for which one or two reference objects per individual are considered, as usually done in the literature. Algorithm 2 has been coded in AMPL, [31] , and the quadratic binary problems have been solved with CPLEX 12.6 [24] . The computational experiments have been carried out on a PC Intel® Core TM i7-2600K, 16GB of RAM. Since problem ( DyViMap ) is expected to be multimodal, since it extends standard MDS (known to be multimodal [58,62] ), DCA may get stuck in local optima. For this reason, Algorithm 2 has been embedded in a multistart routine. We set = [0 , 1] 2 , and λ1 = 0 . 7 , λ2 = 0 . 2 and λ3 = 0 . 1 . The choices of τand κare made dependent on the dataset following the expressions in (3) and (4) , respectively. τ= 1 max t N  i =1 ω i,t ·0 . 05 (3) 130 E. Carrizosa et al. / Omega 86 (2019) 125–136 Fig. 4. U.S. dataset visualization in T = 1 . κ= N(N −1) T T  t=1  i,j=1 , ... ,N i  = j δij,t ·0 . 30 (4) We test the methodology proposed in Section 3 for the problem ( DyViMap ) by visualizing the evolution of the population in the states forming the U.S. across three time periods. The U.S. dataset consists of N = 50 individuals, the states forming the country, for which the population in T = 3 time periods (three years: 1890, 1950 and 2010), ω( t ), has been recorded, [13] . The dissimilarities, δ( t ), depict the geodesic distance between the centroids of the states. These dissimilarities have been computed in R, [54] , running the gdist function to the latitude/longitud coordinates given in state.center data. The dataset measurement is incomplete, in the sense that there are no data available for Alaska (AK) and Hawaii (HI) in 1890. In other words, V (1) consists of all states excepting AK and HI, whereas V (2) and V (3) contain the 50 states. The unit square in R 2 is considered as the reference object for all the states, this is B i = { [0 , 1] 2 }, i = 1 , . . . , N. Since there exists only one reference object per individual, variables x are known (fixed) a priori and the alternating strategy in Algorithm 2 reduces to Algorithm 1 , which finds the 2 ×N ×T = 300 coordinates of the translation vectors c i, t . The maximum number of iterations in Algorithm 1 is set to 100 and the number of iterations of the multistart routine to 50. Initial values of the translation vectors, c ini are uniformly generated in , taking into account the feasible region. Figs. 4–6 show the Dynamic Visualization Map representing the U.S. dataset. Besides the good fitting of dissimilarities, we observe how the squares representing the states are spread over region and they do not suffer big shifts from one period to another. The alternating algorithm for the problem ( DyViMap ), Algorithm 2 , is tested in a linguistic dataset. The dataset consists of the most popular words arising in Danish news around Fig. 5. U.S. dataset visualization in T = 2 . Fig. 6. U.S. dataset visualization in T = 3 . the topic of immigration by year between 1995 and 2015, T = 21 . The relevance of each word per year, ω( t ), is measured by means of the term frequency inverse document frequency (tf-idf) weighting factor, whereas δ( t ) depicts the semantic relatedness of pairs of words using the cosine vector similarity formula, [55] . Then, this similarity is converted into a dissimilarity, values ranging between 0 and 1, with 0 the most similar and 1, the most different. We have a corpus V of N = 203 words, and we set E. Carrizosa et al. / Omega 86 (2019) 125–136 131 Fig. 7. Linguistic dataset visualization in T = 1 , ... , 6 . 132 E. Carrizosa et al. / Omega 86 (2019) 125–136 Fig. 8. Linguistic dataset visualization in T = 7 , ... , 12 . E. Carrizosa et al. / Omega 86 (2019) 125–136 133 Fig. 9. Linguistic dataset visualization in T = 13 , ... , 18 .