Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille1.fr

Luís Paquete
email: paquete@dei.uc.pt

Marco Simões
email: msimoes@student.dei.uc.pt

José R Figueira
email: jose.figueira@mines.inpl-nancy.fr

Connectedness and Local Search for Bicriteria Knapsack Problems

This article reports an experimental study on a given structural property of connectedness of optimal solutions for two variants of the bicriteria knapsack problem. A local search algorithm that explores this property is then proposed and its performance is compared against exact algorithms in terms of running time and number of optimal solutions found. The experimental results indicate that this simple local search algorithm is able to find a representative set of optimal solutions in most of the cases, and in much less time than exact approaches.

Introduction

Stochastic local search algorithms have been applied successfully to many multicriteria combinatorial optimization (MCO) problems. It is widely accepted that their sound performance is related to some structural properties of the solution space that allow local search procedures to find reasonably good quality solutions in an effective manner. However, little is known about which properties these are and how they can affect the performance of this class of algorithms.

In this article, the notion of connectedness of the set of optimal solutions for MCO problems (a.k.a. the efficient set ) [START_REF] Ehrgott | Connectedness of efficient solutions in multiple criteria combinatorial optimization[END_REF] is analyzed from an experimental point of view, and related to the performance of a particular class of multicriteria local search algorithms [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. For a given efficient set of a MCO problem instance, a graph can be constructed such that each node represents an optimal solution and an edge connects two nodes if the corresponding optimal solutions are neighbors for a given neighborhood structure. The efficient set is connected with respect to that neighborhood structure if the underlying graph is also connected, that is, there is a path between any pair of nodes. If the efficient set is connected and the neighborhood structure is tractable from a computational point of view, local search algorithms would be able to find the efficient set in a very effective manner [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF], by starting from at least one efficient solution. However, worst-case results have shown that the efficient set for many MCO problems is not connected in general with respect to different neighborhood structures [START_REF] Ehrgott | Connectedness of efficient solutions in multiple criteria combinatorial optimization[END_REF][START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF], except for very few particular cases [START_REF] Gorski | On a particular case of the multi-criteria unconstrained optimization problem[END_REF][START_REF] Da Silva | Geometrical configuration of the Pareto fronteir of the bi-criteria 0-1-knapsack problem[END_REF]. Moreover, some recent results indicate that approximate solutions that are obtained from independent metaheuristic runs on the bicriteria traveling salesman problem are strongly clustered with respect to small-sized neighborhood structures [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF].

The degree of connectedness is here investigated experimentally for two variants of the bicriteria knapsack problem: The bicriteria unconstrained optimization problem (BUOP) and the bicriteria knapsack problem with bounded cardinality (BKP-BC). Both problems are NP-hard and intractable in the general case [START_REF] Ehrgott | Multicriteria optimization[END_REF]. The experimental results suggest that the efficient set for the two problems above is very often connected with respect to elementary neighborhood structures, despite of the negative results reported in the literature [START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF]. Based on these positive findings, a local search algorithm is proposed and its performance is compared with that of multicriteria dynamic programming algorithms in terms of running-time and number of optimal solutions found. A special technique is introduced that allows the early termination of the complete neighborhood exploration without harming algorithmic performance in terms of solution quality.

The article is organized as follows. Section 2 and Section 3 give the connectedness results for BUOP and BKP-BC, respectively. Moreover, for each of the sections above, the local search algorithms and the exact approaches for the corresponding problem are introduced. Numerical results are shown on a large set of random instances of different structure and size. Finally, Section 4 presents conclusions and further work.

The Bicriteria Unconstrained Optimization Problem

This section introduces a variant of the classical knapsack problem where the capacity constraint is transformed into an additional criterion to be optimized. Then, an experimental analysis on the connectedness property of the corresponding efficient set is reported, as well as the performance of a local search algorithm on several instances of the problem.

Problem Definition

The original (single-criterion) 0/1 knapsack problem is formulated as follows:

max n i=1 p i x i s.t. n i=1 w i x i ≤ W (1)
where p = (p 1 , p 2 , . . . , p j , . . . , p n ) is the profit vector, p j representing the amount of profit on item j, j = 1, . . . , n, and x = (x 1 , x 2 , . . . , x j , . . . , x n ) with x j = 1 if the item j is included in the subset of selected items (knapsack) and x j = 0 otherwise; w = (w 1 , w 2 , . . . , w j , . . . , w n ) is the weight vector, w j representing the amount of investment on item j, j = 1, . . . , n; and W is the overall amount available or budget. The sum of profits and the sum of weights of a given solution x are denoted by p(x) and w(x), respectively.

By transforming the capacity constraint of Problem (1) into a criterion, the following bicriteria unconstrained optimization problem [START_REF] Ehrgott | Multicriteria optimization[END_REF] is obtained: max (p(x), -w(x)) .

(

A proper meaning to the operator "max" above is given as follows. Let X denote the set of feasible solutions of Problem [START_REF] Ehrgott | Multicriteria optimization[END_REF]. The image of the feasible solutions when using the vector maximizing function of Problem (2) defines the feasible region in the criteria space, denoted here by Z ⊆ N 2 . A feasible solution x ∈ X is efficient if there does not exist another feasible solution x ′ ∈ X such that p(x ′ ) ≥ p(x) and w(x ′ ) ≤ w(x), with at least one strict inequality in one of above (or (p(x ′ ), -w(x ′ )) ≥ (p(x), -w(x))). A vector z ∈ Z is nondominated if there is some efficient solution x such that z = (p(x), -w(x)). A vector z ∈ Z dominates a vector z ′ ∈ Z (or z ′ is dominated by z) if z ≥ z ′ holds; if neither z ≥ z ′ nor z ′ ≥ z holds, then both are (mutually) nondominated. The set of all efficient solutions and the set of nondominated vectors are called the efficient set and the nondominated set, respectively. The usual goal of MCO is to find a minimal complete set, that is, the smallest subset of the efficient set whose image coincides with the nondominated set. This subset may not be unique.

Connectedness Analysis

This section describes an experimental analysis for investigating the influence of problem size and degree of conflict between the two criteria on the connectedness property of the efficient set for BUOP. A multicriteria dynamic programming (MDP-BUOP) algorithm is implemented to compute the efficient set. This algorithm consists of the first phase of the Nemhauser-Ullman algorithm for the single-criterion 0/1 knapsack problem [START_REF] Nemhauser | Discrete dynamic programming and capital allocation[END_REF]. It has shown to be theoretically efficient for several input data distributions [START_REF] Beier | Probabilistic analysis of knapsack core algorithms[END_REF]. The MDP-BUOP sequential process consists of n stages. At any stage i, the algorithm generates a set S i of states, which represents a set of promising feasible solutions made up of the first i items, i = 1, . . . , n. A state s = (s p , s w ) ∈ S i represents a feasible solution of profits s p and weight s w . The MDP-BUOP algorithm follows the recursion:

S i := vmax {(s p + p i , s w -w i ), s ∈ S i-1 }
for i = 1, . . . , n, with the basis case S 0 := (0, 0). Operator "vmax" returns the states that are nondominated in S i . At the last stage n, the set S n corresponds to the nondominated set.

In order to obtain the efficient set with the MDP-BUOP algorithm, a binary string is generated with each new state and updated accordingly during the sequential process. For this reason, the implementation keeps states with the same component values. The removal of dominated states at each stage is performed by the algorithm of Kung et al. [START_REF] Kung | On finding the maxima of a set of vectors[END_REF]. Only two sets of states are maintained during the overall sequential process since, at any stage i > 0, only set S i-1 is required.

Table 1. Percentage of BUOP instances connected with respect to the 1-flip neighborhood (all instances are connected with respect to the 1-flip-exchange neighborhood) and average size of the efficient set for each instance size and data correlation (%con and avgs, respectively). These two sets are implemented as height-balanced binary search trees in order to allow logarithmic-time operations. If the implementation does not terminate before one hour of CPU-time, or if RAM resources available are exceeded, the run is cancelled and the output is omitted. The code is written in C++. BUOP instances are defined with two parameters: problem size (n) and correlation between profit and weight vectors (ρ). Both parameters affect the size of the efficient set. The positive (resp. negative) data correlation will increase (resp. decrease) the degree of conflict between the two criteria. The size of the instances ranges from 300 to 3000 and the correlations are ρ ∈ {-0.8, -0.4, 0.0, 0.4, 0.8}. Profit and weight integer values are generated randomly according to a uniform distribution in [1, M/n], where M denotes the maximum possible integer value. The generation of correlated data follows the procedure given by Verel et al. [START_REF] Verel | Analyzing the effect of objective correlation on the efficient set of MNK-landscapes[END_REF]. For each problem size and each correlation degree, 30 different and independent instances are randomly generated.

ρ = -0.8 ρ = -0.4 ρ = 0.0 ρ = 0.4 ρ = 0.
Two neighborhood structures are considered for the experimental analysis: the 1-flip and 1-flip-exchange neighborhoods. Two feasible solutions are 1flip neighbors if they differ exactly on one assignment. In other words, a given neighbor can be reached by adding or removing one item from a given solution. Hence, this neighborhood structure is directly related to the Hamming distance between binary strings. The 1-flip-exchange neighborhood is an extension of the neighborhood above. Two feasible solutions are 1-flip-exchange neighbors if one can be obtained from the other by exchanging two items, adding one item, or removing one item. The size of the neighborhood is linear with n for the 1flip neighborhood structure, while it is quadratic in the case of 1-flip-exchange. Both neighborhoods coincide with the neighborhood structures used by Gorski et al. [START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF] for a similar class of problems.

For the connectedness analysis, MDP-BUOP outputs the efficient set for every instance. For each neighborhood structure, an adjacency matrix is built, indicating whether each two efficient solutions are neighbors or not. Based on this matrix, the connectedness of the corresponding graph is tested. Since this

Algorithm 1 Pareto Local Search Input: n ∈ N, p, w ∈ N n , x0 := 0. Output: Set VT 1: VF := {x0} 2: VT := ∅ 3: while VF = ∅ do 4: select x ⋆ from VF 5: VF := VF \ {x ⋆ } 6: VT := VT ∪ {x ⋆ } 7: for all x ∈ N (x ⋆ ) do 8: if {x ′ | x ′ ∈ VF ∪ VT , (p(x ′ ), -w(x ′ )) ≥ (p(x), -w(x))} = ∅ then 9: VF := {x ′ | x ′ ∈ VF , (p(x), -w(x)) ≥ (p(x ′ ), -w(x ′ ))} 10: VT := {x ′ | x ′ ∈ VT , (p(x), -w(x)) ≥ (p(x ′ ), -w(x ′ ))} 11:
VF := VF ∪ {x} 12: return VT analysis involves a large usage of memory resources (more than 2Gb for largesize instances), only results for a limited number of instance sizes and correlation values are presented. Table 1 gives the percentage of instances with the efficient set that is connected with respect to the 1-flip neighborhood, as well as the average size of the efficient set, rounded to the nearest integer. Although the correlation in the input data influences the size of the efficient set, it does not seem to affect the connectedness results. However, the proportion of instances with a connected efficient set slightly decreases with the increase of the instance size. Finally, for those set of instances, the efficient set is always connected with respect to the 1-flip-exchange neighborhood.

Local Search for BUOP

The local search algorithm for BUOP (PLS-BUOP) is based on the Pareto Local Search [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. The pseudo-code is given in Algorithm 1. For simplification purpose, it is assumed that all feasible solutions have a distinct image in the criterion space. Two archives of nondominated solutions are maintained, V T and V F , respectively. Archive V T contains the set of solutions whose neighborhood has already been explored, while V F contains the remaining ones. PLS-BUOP starts with an efficient solution that initializes the archive. Then, at each iteration, a solution is chosen from V F , and its neighborhood is explored (N (x) denotes the neighborhood of a given solution x ∈ X). All the nondominated neighboring solutions are used to update V F and dominated solutions are discarded from V F and V T . The algorithm terminates once V F is empty. This algorithm stops naturally when a Pareto local optimum set is found, it does not cycle, and if connectedness of the efficient set holds, it is able to identify it by starting from at least one efficient solution (with a proper change of the conditions in Algorithm 1) [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. In the following, some particular details of the implementation are described. -Initialization. The algorithm starts with the efficient solution x 0 = 0. This solution maps, in the criteria space, to an extreme point from the nondominated set, with null profit and weight values. -Selection. At each iteration of the algorithm, the solution from V F with the smallest profit value is selected (Algorithm 1, line 4). -1-flip neighborhood exploration. In order to perform the 1-flip neighborhood exploration in an efficient manner, a preliminary step ranks the items into different layers, with respect to the dominance-depth ranking [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF]. Two different ranks are required to cover the cases of adding and deleting an item. This step is illustrated in Fig. 1. Let L i denote the set of items in the i-th layer and let u ∈ L i and v ∈ L j , j > i. Then, for the case of the addition, it holds that (p v , -w v ) ≥ (p u , -w u ). Therefore, the neighborhood exploration starts by examining the items in the first layer, and proceeds with the items in the subsequent layers. Within the same layer, the exploration follows the nondecreasing order of the weights. The exploration stops after verifying that no item of a given layer belongs to the current solution. Indeed, any neighboring solution constructed from the subsequent layers is dominated by at least one neighboring solution built from this layer. For the case of deletion, a similar reasoning applies with the corresponding changes. -1-flip-exchange neighborhood exploration. For this neighborhood, the same reasoning as above applies when adding or removing an item. In order to exchange pairs of items more efficiently, the following pre-processing procedure is applied: First, for each item i ∈ {1, . . . , n}, a tuple records the profit and weight difference with respect to each different item j. Then, the n -1 tuples are sorted in terms of dominance-depth ranking. When considering the exchange of item i with another item, the exploration follows the order given by the dominance-depth ranking, with the nondecreasing order of the weighs for tuples within the same layer; the exploration stops once the items corresponding to all tuples of a given layer can be exchanged with i. The exploration is iterated for every item in the knapsack. -Data structures. Archives V F and V T are implemented as height-balanced binary search trees. The removal of dominated solutions follows the algorithm of Kung et al. [START_REF] Kung | On finding the maxima of a set of vectors[END_REF].

MDP-BUOP and two PLS-BUOP versions, each one for each neighborhood structure, were run on the instances described in Section 2.2. The implementation of MDP-BUOP was modified in order to output a minimal complete set. All the algorithms share the same programming language, data structures, compilation options and were run in the same machine with a time bound of one hour. The PC used for the experiments was an iMac with Mac OS X v10.5.5 (Leopard), 2.4 GHz Intel Core with 4MB L2 Cache and 2GB SDRAM. All codes were compiled with g++ version 4.0.1 using the -O3 flag. The CPU-time taken by MDP-BUOP and PLS-BUOP with 1-flip neighborhood for instances with correlation ρ = -0.4 and ρ = 0.4 is reported in Fig. 2. In order to distinguish the differences of performance, the CPU-time is presented in log scale. Similar relationship between performance of both implementations was obtained for the remaining correlation values. Table 2 presents the percentage of efficient solutions and the total number of solutions returned by PLS-BUOP using 1flip neighborhood, averaged over results obtained in 30 instances for each size and correlation. Since PLS-BUOP using the 1-flip-exchange neighborhood took always more time than MDP-BUOP, the results of the former are omitted.

Although the experimental analysis indicated the existence of instances with unconnected efficient set with respect to the 1-flip neighborhood (see Table 1), the results show that the local search approach using the same notion of neighborhood is able to identify a minimal complete set in many cases. For other instances, it leads to the identification of more than 99.9% of the efficient set. Furthermore, the local search performs very efficiently in comparison to the dynamic programming approach in terms of computational time. Indeed, the larger the instance size, the larger the gap between PLS-BUOP and MDP-BUOP in Table 2. Average percentage of efficient solutions and average number of solutions found by PLS-BUOP using the 1-flip neighborhood (%ef and avgs, respectively). ρ = -0.8 ρ = -0.4 ρ = 0.0 ρ = 0.4 ρ = 0.8 size %ef avgs %ef avgs %ef avgs %ef avgs %ef avgs 300 100.0 8 965 99.9 10 838 100.0 14 501.2 100.0 24 702 99.9 55 159 600 99.9 32 242 100.0 39 056 100.0 52 675.4 100.0 92 717 100.0 228 427 900 99.9 69 208 99.9 83 357 99.9 112 597.3 99.9 201 642 99.9 495 153 1200 100.0 118 483 99.9 144 259 99.9 195 395.6 100.0 338 903 99.9 823 611 1500 99.9 179 932 100.0 217 230 100.0 301 844.6 99.9 526 029 99.9 1 352 816 1800 99.9 252 972 100.0 308 373 99.9 423 450.0 99.9 716 598 99.9 1 818 607 2100 99.9 337 443 99.9 409 770 99.9 563 839.5 99.9 969 068 99.9 2 431 713 2400 100.0 434 333 99.9 530 843 100.0 719 388.3 99. [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF] terms of CPU time. However, PLS-BUOP appears to be much more efficient for negatively correlated data, while MDP-BUOP was not able to solve all the instances for positively correlated profit and weight values.

The Bicriteria Knapsack Problem with Bounded Cardinality

This section reports a similar analysis on a variant of the previous problem, where an additional cardinality constraint is considered.

Problem Definition

The bicriteria knapsack problem with bounded cardinality is a variant of the BUOP that is obtained by limiting the number of chosen items by a cardinality bound (k). The formulation of BKP-BC is as follows: max (p(x), -w(x))

s.t. n i=1 x i ≤ k (3)
In the problem above, the operator "max" follows the same meaning as in Problem [START_REF] Ehrgott | Multicriteria optimization[END_REF]. The same terminology and notation will be used for this problem.

Connectedness Analysis

The efficient set is computed by means of a multicriteria dynamic programming (MDP-BKP-BC) algorithm that extends the MDP-BUOP algorithm given in Section 2.2. This algorithm has k × n stages. At each stage, the algorithm generates the set T (i,j) of states, which represents a set of potential efficient solutions made up of the first i items, i = 1, . . . , n, with cardinality j, j = 1, . . . , k.

A state t = (t p , t w ) ∈ T (i,j) represents a feasible solution of profits t p and weight t w . This approach follows the recurrence relation:

T (i,j) := vmax T (i-1,j) ∪ (t s + p i , t w -w i ), t ∈ T (i-1,j -1)
for i = 1, . . . , n and j = 1, . . . , k, with the basis cases T (i,0) := (0, 0), for i = 0, . . . , n and T (0,j) := 0) for j = 0, . . . , k. The nondominated set of states is given by the set vmax

T (n,0) ∪ • • • ∪ T (n,k) .
The implementation follows the same principles presented in Section 2.2 for the MDP-BUOP algorithm. However, the implementation of MDP-BKP-BC has to keep 2(k +1) sets during the overall process, since at each state i, the sets T (i,j) and T (i-1,j) , for j = 0, . . . , k are required. For this reason, MDP-BKP-BC should take more time than MDP-BUOP for the same instance size.

A similar set of instances to those used for the previous problem is generated. In addition to the instance size (n) and to the data correlation (ρ), several different values for the cardinality bound (k) are considered: n/10, n/5, and n/2. For each problem size, correlation and cardinality bound, 30 different and independent instances are generated randomly, as explained in Section 2.2. The two neighborhood structures described in Section 2.2 were also considered for this problem. The use of the 1-flip-exchange neighborhood for this problem is motivated by the conjecture that many efficient solutions may have the same cardinality, due to the additional cardinality constraint.

Due to limited memory resources, results were only obtained for instances of limited size (these instances can be inferred from Table 3). Differently from the connectedness results obtained in the first problem (see Table 1), no instance with an efficient set that is connected with respect to the 1-flip neighborhood was found. However, the efficient set for all the instances were connected with respect to the 1-flip-exchange neighborhood. These results corroborate those of Gorski et al. [START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF] for much smaller instances (up to 100 items).

Local Search for BKP-BC

Given the positive results reported in the previous section, a local search (PLS-BKP-BC) was developed under the same reasoning of PLS-BUOP (see Section 2.3). The only difference is in the neighborhood exploration since a maximum number of items has to be ensured in the solution when considering the possibility of adding an item.

MDP-BKP-BC and two versions of PLS-BKP-BC were run in the same instances defined in the previous section. The experiments were performed in a computer cluster with 6 nodes, each with an AMD Phenom II X6 processor with 3.2GHz, 3 and 6 MB L2 and L3 Cache, respectively, and 12 MB DDR3 SDRAM. The operating system was Ubuntu 8.04 LTS. Both codes were compiled with gcc version 4.2.4 using the -O3 flag.

The CPU-time taken by the three approaches is plotted in Fig. 3 for four different instance parameter settings. In order to distinguish the differences of performance, the CPU-time is presented in log scale. Clearly, both PLS-BKP-BC versions take much less time than MDP-BKP-BC to terminate. Similar results hold for the remaining instances. Table 3 reports the percentage of efficient solutions found by the two versions of PLS-BKP-BC. As expected from the connectedness analysis reported in the previous section, PLS-BKP-BC using the 1-flip-exchange neighborhood was always able to find a minimal complete set for all instances up to 1800 items in less than one hour of CPU-time. PLS-BKP-BC using 1-flip neighborhood is, in many cases, able to find more than half of a minimal complete set in less than one second.

Concluding Remarks

This article describes an experimental analysis on the structure of the efficient set, in terms of connectedness, for two MCO problems. Despite of the negative results reported in the literature for similar problems [START_REF] Ehrgott | Connectedness of efficient solutions in multiple criteria combinatorial optimization[END_REF][START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF], the experimental analysis for the problems investigated in this paper are quite promising. For both bicriteria versions of the unconstrained optimization problem and knapsack problem with a bounded cardinality constraint, the experiments suggest that small-sized neighborhood structures give rise to connected efficient sets quite Table 3. Average percentage of efficient solutions found by PLS-BKP-BC using the 1-flip and the 1-flip-exchange neighborhood structures (%ef1 and %ef2, respectively). Although the large number of connected instances motivates the use of local search algorithms, it is still an open question whether those approaches are efficient enough as compared to exact algorithms. The experimental analysis reported in this article gives a clear positive answer for the second problem. For the first problem, without cardinality constraint, some preliminary results indicated that the local search proposed in this article under the same neighborhood that (empirically) provides connectedness would not be worthwhile in terms of running time. Still, using a smaller neighborhood structure allows the same algorithm to find more than 99.9% of the efficient set in a significantly less amount of time than the exact approach.

The simplicity of the local search approach proposed in this article is very appealing for implementation purpose, needs no definition of parameters and requires a minimum number of modifications in order to be applied to other type of knapsack problems. For instance, the same principles can be applied to the multicriteria knapsack problem (with several maximizing profit criteria and one capacity constraint) by ignoring (or penalizing) infeasible neighboring solutions. However, finding appropriate definitions of neighborhoods that give rise to a large number of connected efficient sets for knapsack problems with capacity constraints is still under investigation.

A natural question is whether it is possible to derive analytical results for MCO problems that would prove connectedness by assuming some structure or distribution on the input data. Connections with neighborhood structures arising in the context of the linear programming formulation of the MCO problem may provide further insights [START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF]. Moreover, the derivation of bounds on the run-time of multiobjective evolutionary algorithms that start from efficient solutions are also of interest. For instance, the size of the efficient set for the first problem is polynomially bounded for many input data distributions [START_REF] Beier | Probabilistic analysis of knapsack core algorithms[END_REF]. This suggests that a polynomial run-time bound may be achieved by such type of algorithms.

Fig. 1 .

 1 Fig. 1. Initial dominance-depth ranking of items for the exploration of the neighborhood (left: ranking for adding, right: ranking for deleting).

Fig. 2 .

 2 Fig. 2. CPU time in seconds (average value and standard deviation, given in log scale) of MDP-BUOP and PLS-BUOP using the 1-flip neighborhood for instances with data correlation ρ = -0.4 (left) and ρ = 0.4 (right).

Fig. 3 .

 3 Fig. 3. CPU time in seconds (average value and standard deviation, given in log scale) of the three approaches for instances with parameters k = n/10 (top) and k = n/5 (bottom), ρ = -0.4 (left) and ρ = 0.4 (right).

  card. %ef1 %ef2 %ef1 %ef2 %ef1 %ef2 %ef1 %ef2 %ef1 %ef2 300 n/10 72.6 100.0 54.1 100.0 42.0 100.0 44.2 100.0 47.7 100.0 n/5 83.9 100.0 63.0 100.0 53.5 100.0 54.4 100.0 56.3 100.0 n/2 95.6 100.0 88.6 100.0 85.3 100.0 83.1 100.0 84.7 100.0 600 n/10 75.3 100.0 50.9 100.0 41.6 100.0 41.6 100.0 47.6 100.0 n/5 84.2 100.0 63.3 100.0 51.3 100.0 50.0 100.0 55.3 100.0 n/2 95.4 100.0 88.3 100.0 83.1 100.0 83.4 100.0 --900 n/10 75.6 100.0 50.1 100.0 39.6 100.0 40.5 100.0 47.0 100.0 n/5 83.1 100.0 61.0 100.0 50.4 100.0 51.5 100.0 --1200 n/10 75.6 100.0 50.5 100.0 40.2 100.0 39.9 independently of the size and of the structure of input data. In fact, it is not yet clear what structure of the input data may generate, in general, an unconnected efficient set under the 1-flip-exchange neighborhood that was used in this article.

  8 size %con avgs %con avgs %con avgs %con avgs %con

										avgs
	300 100.0	8 965 96.7 10 838 100.0 14 501 100.0 24 702 93.3	55 159
	600 96.7 32 242 100.0 39 056 100.0 52 675 96.7 92 717 90.0 22 8427
	900 90.0 69 208 93.3 83 357 86.7 112 598 76.7 201 642 60.0 495 154
	1200 96.7 118 483 90.0 144 259 86.7 195 396 83.3 338 903	-	823 612
	1500 90.0 179 932 93.3 217 230 86.7 301 845	-	526 029	-	1 352 817
	1800 83.3 252 972 93.3 308 373	-	423 450	-	716 598	-	1 818 608
	2100 76.7 337 443	-	409 771	-	563 840	-	969 069	-	2 431 715

Acknowledgements. The authors acknowledge Jochen Gorski for the discussion on the main topic of this article. This work was partially supported by the Portuguese Foundation for Science and Technology (PTDC/EIA-CCO/098674/2008) and the project "Connectedness and Local Search for Multiobjective Combinatorial Optimization" funded by the Deutscher Akademischer Austausch Dienst and Conselho de Reitores das Universidades Portuguesas. The third author acknowledges a CEG-IST grant from Instituto Superior Técnico (PTDC/GES/73853/2006).


Roland Badeau
email: roland.badeau@enst.fr.

Rémy Boyer
email: remy.boyer@lss.supelec.fr

MULTILINEAR SINGULAR VALUE DECOMPOSITION FOR STRUCTURED TENSORS

Keywords: Multilinear SVD, fast algorithms, structured and unstructured tensors AMS subject classifications. 15A15, 15A09, 15A23

The Higher-Order SVD (HOSVD) is a generalization of the Singular Value Decomposition (SVD) to higher-order tensors (i.e. arrays with more than two indices) and plays an important role in various domains. Unfortunately, this decomposition is computationally demanding. Indeed, the HOSVD of a third-order tensor involves the computation of the SVD of three matrices, which are referred to as "modes", or "matrix unfoldings". In this paper, we present fast algorithms for computing the full and the rank-truncated HOSVD of third-order structured (symmetric, Toeplitz and Hankel) tensors. These algorithms are derived by considering two specific ways to unfold a structured tensor, leading to structured matrix unfoldings whose SVD can be efficiently computed 1 .

1. Introduction. The subject of multilinear decomposition is now mature [START_REF] Comon | Tensor Decompositions, state of the art and applications[END_REF][START_REF] Sidiropoulos | Low-Rank Decomposition of Multi-Way Arrays: A Signal Processing Perspective[END_REF]. There are essentially two families. The first one is known under the name of CANDECOMP/PARAFAC (CANonical DECOMPosition or PARAllel FACtors model) and was independently proposed in [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition[END_REF][START_REF] Harshman | PARAFAC: Parallel factor analysis[END_REF]. This decomposition is very useful in several applications and is linked to the tensor rank [START_REF] Howell | Global properties of tensor rank[END_REF]. The second one is related to the multidimensional rank [START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF] and is known under the name of Tucker decomposition [START_REF] Tucker | Some Mathematical Notes on Three-mode Factor Analysis[END_REF]. This decomposition is a more general form which is often used. Orthogonality constraints are not required in the general Tucker decomposition but if needed one can refer to the Higher-Order Singular Value Decomposition (HOSVD) [START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF] or multilinear SVD.

The HOSVD is a generalization of the SVD to higher-order tensors (ie. arrays with more than two indices). This decomposition plays an important role in various domains, such as harmonic retrieval [START_REF] Papy | Exponential data fitting using multilinear algebra. The single-channel and the multichannel case[END_REF], image processing [START_REF] Huang | Multi-way methods in image analysisrelationships and applications[END_REF], telecommunications, biomedical applications (magnetic resonance imaging and electrocardiography), web search [START_REF] Sun | CubeSVD: A Novel Approach to Personalized Web Search[END_REF], computer facial recognition [START_REF] Wang | Facial Expression Decomposition[END_REF], handwriting analysis [START_REF] Savas | Analyses and Tests of Handwritten Digit Recognition Algorithms[END_REF], statistical methods involving Independent Component Analysis (ICA) [START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF].

In [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF], it was shown that the HOSVD of a third-order tensor involves the computation of the SVD of three matrices called modes, leading to a high computational cost. A first approach for reducing the complexity of tensor-based methods consists in a dimensionality reduction: only the principal components of the HOSVD are calculated, leading to the rank-truncated HOSVD. In this paper, we present a standard and fast algorithm for calculating the full and the rank-truncated HOSVD, which only computes the left factors of the three SVD's. Next, we focus on structured tensors, such as symmetric and Toeplitz tensors, which naturally arise in signal processing methods involving higher-order statistics [START_REF]Fast Reliable Algorithms for Matrices with Structure[END_REF]Chapter 9], and Hankel tensors [START_REF] Papy | Exponential data fitting using multilinear algebra. The single-channel and the multichannel case[END_REF], introduced in the context of the Harmonic Retrieval problem [START_REF] Liu | Almost Sure Identifiability of Constant Modulus Multidimensional Harmonic Retrieval[END_REF], which is at the heart of many signal processing applications. To the best of our knowledge, there are no specific HOSVD algorithms proposed in the literature for exploiting tensors structures. In this paper however, we show that such tensors can be efficiently decomposed. We first observe that standard unfoldings [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF] do not present a particularly noticeable structure even in the case of structured tensors. Consequently, we introduce two different ways to unfold a structured tensor which clarify the link between structured modes and structured tensors. By doing this, we can exploit fast product techniques [START_REF] Golub | Matrix computations[END_REF]. A second point of this work concerns Hankel and symmetric tensors. The modes of these structured tensors are column-redundant so it is possible to reduce the computational cost of the HOSVD algorithm by taking the redundant structure of each mode into account. Finally, our fastest implementation of the rank-truncated HOSVD (dedicated to Hankel tensors) has a quasilinear complexity with respects to the tensor dimension.

Note that for applications involving very large tensor dimensions, an even lower complexity may be required. In this case, one may be interested in rank-revealing tensor decompositions which can be computed faster than the rank-truncated HOSVD. Such an approach is developed in [START_REF] Oseledets | Tucker dimensionality reduction in linear time[END_REF], based on cross approximation techniques which are derived from LU factorizations [START_REF] Tyrtyshnikov | Incomplete cross approximation in the mosaic-skeleton method[END_REF][START_REF] Bebendorf | Approximation of boundary element matrices[END_REF]. An algorithm is proposed which provides a Tucker-like low rank approximation of unstructured cube tensors, the complexity of which is linear with respects to the tensor dimension in many cases [START_REF] Oseledets | Tucker dimensionality reduction in linear time[END_REF]. This linear complexity is nevertheless obtained via an approximated rank reduction, in comparison with an exact Tucker decomposition such as the HOSVD.

2. Preliminaries in multilinear algebra. We present some basic definitions in the context of third-order tensor algebra. These definitions can be extended to order greater than three and we refer the interested reader to [START_REF] Comon | Tensor Decompositions, state of the art and applications[END_REF][START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF] for instance.

Tucker's product. The Tucker's product, also called s-mode product, of a thirdorder complex-valued tensor A ∈ C I1×I2×I3 by a matrix B ∈ C Js×Is for s ∈ [1 : 3] is defined according to:

[A × 1 B] j1i2i3 = I1-1 i1=0 [A] i1i2i3 [B] j1i1 , (2.1)
[A × 2 B] i1j2i3 = I2-1 i2=0 [A] i1i2i3 [B] j2i2 , (2.2)
[A × 3 B] i1i2j3 = I3-1 i3=0 [A] i1i2i3 [B] j3i3 , (2.3)
where we denoted the entries of A by [A] i1i2i3 with i s ∈ {0 . . . I s -1}. We have the following properties:

A × s B × s ′ C = (A × s B) × s ′ C = (A × s ′ C) × s B, (2.4)
(A × s B) × s C = A × s (BC). (2.5)
Mode of a tensor. There are several ways to represent an I 1 × I 2 × I 3 third-order complex-valued tensor A as a collection of matrices.

Definition 2.1. The modes (also called "matrix unfoldings") A 1 , A 2 , A 3 are usually defined as follows:

[A 1 ] i1,i2I3+i3 = [A] i1i2i3 , (2.6) [A 2 ] i2,i3I1+i1 = [A] i1i2i3 , (2.7
)

[A 3 ] i3,i1I2+i2 = [A] i1i2i3 . (2.8)
These matrices are of dimension (I 1 × I 2 I 3 ), (I 2 × I 3 I 1 ), (I 3 × I 1 I 2 ), respectively. The dimensions of the vector spaces generated by the columns of the modes of A are called column rank (or 1-mode rank) R 1 , row rank (or 2-mode rank) R 2 and 3-mode rank R 3 , respectively.

Multilinear SVD (HOSVD).

Theorem 2.2 (Third-Order SVD [START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF][START_REF] Tucker | Some Mathematical Notes on Three-mode Factor Analysis[END_REF]). Every I 1 × I 2 × I 3 tensor A can be written as the product:

A = S × 1 U (1) × 2 U (2) × 3 U (3)
(2.9)

where × s represents the This decomposition is a generalization of the SVD because the diagonality of the matrix containing the singular values, in the matrix case, is a special case of allorthogonality. Also, the HOSVD of a second-order tensor (matrix) yields the matrix SVD, up to trivial indeterminacies. The matrix of s-mode singular factors, U (s) , can be found as the matrix of left singular vectors of the mode A s , defined in (2.6)-(2.8). The s-mode singular values correspond to the singular values of this matrix unfolding. Note that the s-mode singular factors of a tensor, corresponding to the nonzero smode singular values, form an orthonormal basis for its s-mode vector subspace, like in the matrix case.

The core tensor S can then be computed (if needed) by bringing the matrices of s-mode singular factors to the left side of equation (2.9):

S = A × 1 U (1) H × 2 U (2) H × 3 U (3) H (2.10)
where (.) H denotes the conjugate transpose. Mode decompositions. Expression (2.9) can be written in terms of modes as follows:

A 1 = U (1) S 1 U (3) ⊗ U (2) H , A 2 = U (2) S 2 U (3) ⊗ U (1) H , A 3 = U (3) S 3 U (1) ⊗ U (2) H ,
where ⊗ denotes the Kronecker product and S 1 , S 2 and S 3 denote respectively the first, second and third modes of the core tensor S.

HOSVD algorithm for unstructured tensors.

In this section we present an efficient implementation of the HOSVD in the general framework of unstructured tensors, from which our fast algorithms for structured tensors will be derived in section 4. Let I = 1 3 (I 1 + I 2 + I 3 ). The computational costs of the various algorithms presented below are related to the flop (floating point operation) count. For example, a dot product of I-dimensional vectors involves 2I flops (I multiplications plus I additions).

The calculation of the HOSVD of tensor A requires the computation, for all s ∈ [1 : 3], of the left factor U (s) in the full SVD of matrix A s , as defined above. Note that in many applications, we are interested in computing the HOSVD truncated at ranks (M 1 , M 2 , M 3 ), which means that we only compute the M s first columns of the matrix U (s) (M s is often supposed to be much lower than I s ). We will suppose throughout this paper that this possibly truncated SVD is computed by means of the orthogonal iteration method, although other algorithms such as the Golub-Reinsch SVD and R-SVD [7, pp. 253-254] could also be applied. When computing only the n × r left factor U in the rank r-truncated SVD of an n × m matrix A with n < m, the orthogonal iteration method consists in recursively computing the n × r matrix

B i = A(A H U i-1
), involving 2r matrix / vector products, and the QR factorization [7, pp. 410-411]. Thus the computational cost of one iteration is 2r c(n, m) + 2r 2 n flops, where c(n, m) = 2nm is the cost of 1 matrix / vector product, and 2r 2 n is the cost of 1 QR factorization [7, pp. 231-232]. Besides, the s-mode A s has n = I s rows and m = s ′ =s I s ′ columns. Assuming that s ′ =s I s ′ is much greater than I s , the dominant cost of one iteration for computing U (s) is 4M s I 1 I 2 I 3 flops. Finally, the Tucker product (2.10) can be computed by folding for instance its first mode given by

B i = U i R i of this n × r matrix
S 1 = U (1) H A 1 U (3) ⊗ U (2) .
(

A fast implementation of equation (2.11) was proposed in [START_REF] Anderson | Improving the speed of multi-way algorithms: Part I. Tucker3[END_REF], whose complexity is

6M s I 1 I 2 I 3 flops, where M = 1 3 (M 1 + M 2 + M 3 ).
Note that the computation of the Tucker product is generally not needed in applications, this is why it will be omitted in the following developments.

The computational cost of the full and rank-truncated HOSVD is summarized in table 2.1 (the full HOSVD is the same as the rank-truncated HOSVD with M s = I s for all s = 1, 2, 3). In this table and below, the global cost is provided as a maximum w.r.t. I 1 , I 2 , I 3 , under the constraint I 1 + I 2 + I 3 = 3I. In particular, the maximal complexity per iteration is obtained for cube tensors (I 1 = I 2 = I 3 = I) and equals 12M I 3 .

Operation Cost per iteration SVD of A 1 4M 1 I 1 I 2 I 3 SVD of A 2 4M 2 I 1 I 2 I 3 SVD of A 3 4M 3 I 1 I 2 I 3 Global cost 12M I 3
3. Structured tensors and reordered tensor modes.

In this section, we present three tensor structures which are usual in many applications. Next, we introduce new reordered tensor modes which clarify the link between structured tensors and structured modes.

3.1. Structured tensors. Definition 3.1 (Toeplitz tensors). A Toeplitz tensor is a structured tensor which satisfies the following property: for all i 1 ∈ {0 . . .

I 1 -1}, i 2 ∈ {0 . . . I 2 -1}, i 3 ∈ {0 . . . I 3 -1}, ∀k ∈ {0 . . . min(I 1 -i 1 , I 2 -i 2 , I 3 -i 3 ) -1}, [A] i1+k,i2+k,i3+k = [A] i1i2i3 .
Below, any permutation of 3 elements will be denoted π = (π 1 , π 2 , π 3 ) where π 1 , π 2 , π 3 ∈ {1, 2, 3}, according to the following definition

π : (i 1 , i 2 , i 3 ) → (i π1 , i π2 , i π3 ).
Definition 3.2 (Symmetric tensors). A cube (I × I × I) tensor A which is unchanged by any permutation π is called a symmetric tensor:

∀i 1 , i 2 , i 3 ∈ {0, . . . , I -1}, [A] π(i1,i2,i3) = [A] i1i2i3 .
Example 1 (Fast higher PCA for real moment and cumulant). The HOSVD can be viewed (cf. reference [START_REF] Kroonenberg | Three-Mode Principal Component Analysis: Theory and Applications[END_REF] ) as a higher Principal Component Analysis (PCA). This technique is often used as a data dimensional reduction for moment and cumulant tensors [START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF]. Third-order moment and cumulant tensors are defined according to

[M] t1t2t3 = E{x(t 1 )x(t 2 )x(t 3 )}, (3.1)
[C] t1t2t3 = E{x(t 1 )x(t 2 )x(t 3 )} + 2E{x(t 1 )}E{x(t 2 )}E{x(t 3 )} -E{x(t 1 )}E{x(t 2 )x(t 3 )} -E{x(t 2 )}E{x(t 1 )x(t 3 )} -E{x(t 3 )}E{x(t 1 )x(t 2 )}. (3.
2) where t 1 , t 2 , t 3 ∈ {0 . . . I -1}, and x(t) is a real random process.

Moment and cumulant tensors, defined in (3.1) and (3.2), are symmetric tensors according to definition 3.2. The proof is straightforward and can be generalized to larger orders [START_REF] Comon | Tensor Decompositions, state of the art and applications[END_REF]. Moreover, if x(t) is a third-order stationary process, the moment and cumulant tensors defined in (3.1) and (3.2) are third-order Toeplitz tensors according to definition 3.2. Indeed, if x(t) is a stationary process, its probability distribution is invariant to temporal translations. This property implies

[C] t+i1,t+i2,t+i3 = [C] i1i2i3 .

Definition 3.3 (Hankel tensors). A Hankel tensor is a structured tensor whose coefficients

[A] i1i2i3 only depend on i 1 + i 2 + i 3 .
Note that a cube Hankel tensor is symmetric. Hankel tensors were introduced in [START_REF] Papy | Exponential data fitting using multilinear algebra. The single-channel and the multichannel case[END_REF] in the context of the Harmonic Retrieval problem [START_REF] Liu | Almost Sure Identifiability of Constant Modulus Multidimensional Harmonic Retrieval[END_REF]. This problem is at the heart of many signal processing applications.

Example 2 (Definition and properties of the harmonic model).

We consider the complex harmonic model defined according to:

x n = M m=1 α m z n m , for n ∈ [0 : N -1] (3.3)
where N is the analysis duration and M is the known number of components, z m = e δm+iφm is called the m th pole of x n where i = √ -1, φ m is called the m th angularfrequency belonging to (-π, π], and δ m is the m th damping factor. In the sequel, we assume that all the poles are distinct. In addition, α m = a m e iφm is the non-zero m th complex amplitude, ie., a m = 0, ∀m. Besides, we define the Vandermonde matrices Z (I1) , Z (I2) and Z (I3) associated to model (3.3), according to [Z (Is) ] n,m = z n m , and we assume that M ≤ min(I 1 , I 2 , I 3 ). Then the Hankel tensor

[A] i1i2i3 = x (i1+i2+i3) associated to model (3.3) is diagonalizable according to A = D × 1 Z (I1) × 2 Z (I2) × 3 Z (I3)
where × i denotes the i-th Tucker's product and

[D] jkℓ = α j if j = k = ℓ 0 otherwise
is a hyper-cubic M ×M ×M super-diagonal core tensor. As a consequence, the Hankel tensor A is a rank-(M, M, M ) tensor. Following standard subspace-based parametric estimation methods, the harmonic model can then be estimated by computing the rank M -truncated HOSVD of tensor A [START_REF] Papy | Exponential data fitting using multilinear algebra. The single-channel and the multichannel case[END_REF].

Modes of structured tensors.

As mentioned in the introduction, standard unfoldings of structured tensors [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF] do not present a particularly noticeable structure. Consequently, we introduce in this section two different ways to unfold a structured tensor which clarify the link between structured modes and structured tensors.

Example 3. Consider the 4 × 4 × 4 symmetric and Toeplitz tensor [A] ijk = 3 * max(i, j, k) -sum(i, j, k). The classical 1-mode is formed of 4 symmetric submatrices:

A 1 =         0 2 4 6 2 1 3 5 4 3 2 4 6 5 4 3         2 1 3 5 1 0 2 4 3 2 1 3 5 4 3 2         4 3 2 4 3 2 1 3 2 1 0 2 4 3 2 1         6 5 4 3 5 4 3 2 4 3 2 1 3 2 1 0         .
However, by permuting its columns, we define an other mode A ′ 1 (referred to below as the type-1 reordered tensor mode), which is formed of 4 Toeplitz matrices T 0 , T 1 , T 2 and T 3 (referred to below as the type-1 oblique sub-matrices):

A ′ 1 =           6 5 4 3     T3     4 5 3 4 2 3 4 2     T2     2 3 4 1 2 3 3 1 2 5 3 1     T1     0 1 2 3 2 0 1 2 4 2 0 1 6 4 2 0     T0     2 3 4 1 2 3 3 1 2 5 3 1     T1     4 5 3 4 2 3 4 2     T2     6 5 4 3     T3       .
It can be noted that this reordered mode satisfies an axial blockwise symmetry with respect to its central oblique sub-matrix T 0 . Obviously, the left singular factor U (1) in the SVD of A 1 is the same as the left singular factor in the SVD of A ′ 1 , since both matrices have the same columns. However, we will show below that the SVD of A ′ 1 can be computed efficiently, by exploiting the Toeplitz structure of the oblique sub-matrices T k .

Example 4. Consider the 4×4×4 Hankel tensor [A] ijk = i+j +k. The standard 1-mode is formed of 4 Hankel sub-matrices: However, by permuting its columns, we define an other mode A ′′ 1 (referred to below as the type-2 reordered tensor mode), which is formed of 7 rank-1 matrices R 0 . . . R 6 (referred to below as the type-2 oblique sub-matrices):

A 1 =         0 1 2 3 1 2 3 4 2 3 4 5 3 4 5 6        
A ′′ 1 =           0 1 2 3     R0     1 1 2 2 3 3 4 4     R1     2
    R5     6 7 8 9     R6       .
Again, the left singular factor in the SVD of A 1 is the same as the left singular factor in the SVD of A ′′ 1 , since both matrices have the same columns. However, we will show below that the SVD of A ′′ 1 can be computed efficiently, by exploiting the rank-1 structure of the oblique sub-matrices R k , and the Hankel structure of the matrix obtained by removing the repeated columns in A ′′ 1 . In the following section, the type-1 and type-2 oblique sub-matrices will be defined in the general case by "slicing" a third-order tensor according to 2 different oblique directions, as shown in Fig. 3.1.

3.2.1.

Oblique sub-matrices of a tensor. Definition 3.4 (Type-1 and type-2 oblique sub-matrices of a tensor). For any permutation π, the oblique sub-matrices of a tensor A are defined as follows:

• For all k ∈ {0 . . . I π3 -1}, let J

(π2,π3) (k) = min(I π2 , I π3 -k). The coefficients of the k th type-1 oblique

I π1 × J (1) (π2,π3) (k) sub-matrix of A are [T (π) k ] ij = [A] π -1 (i,j,k+j) (3.4)
where

0 ≤ i ≤ I π1 -1 and 0 ≤ j ≤ J (1) (π2,π3) (k) -1. • For all k ∈ {0 . . . I π2 + I π3 -2}, let J (2) (π2,π3) (k) = min(I π2 , I π3 , 1 + k, I π2 + I π3 -1 -k).
(3.5)

The coefficients of the I π1 × J

(2)

(π2,π3) (k) type-2 oblique sub-matrix of A are [R (π) k ] ij = [A] π -1 (i, max(k-Iπ 3 +1,0)+j, min(k,Iπ 3 -1)-j) (3.6)
where 0 ≤ i ≤ I π1 -1 and 0 ≤ j ≤ J

(2)

(π2,π3) (k) -1. Y I π3 - ? I π2 I π1 - T (π) 0 T (π) 1 - T (π) Iπ 3 -1 - . . . A ? T (π1,π3,π2) 1 ? T (π1,π3,π2) I2-1 • • • - ? Y A I π2 I π3 I π1 - R (π) 0 R (π) 1 - R (π) Iπ 3 -1 - . . . 6 R (π) Iπ 3 6 R (π) Iπ 2 +Iπ 3 -2 • • • Type-1 oblique sub-matrices
Type-2 oblique sub-matrices k ] ij = [A] i,j,k+j , which does not depend on π. 2. Applying equation (3.4) to i + 1 and j + 1 (for all 0 ≤ i < I π1 -1 and 0 ≤ j < J

(1) (3.6) shows that for all permutation π and index k ∈ {0 . . .

(π2,π3) (k) -1) yields [T (π) k ] i+1,j+1 = [A] π -1 (i+1,j+1,k+j+1) . However, since the tensor A is Toeplitz, [A] π -1 (i+1,j+1,k+j+1) = [A] π -1 (i,j,k+j) = [T (π) k ] ij . Therefore [T (π) k ] i+1,j+1 = [T (π) k ] ij , which means that the matrix T k is Toeplitz. 3. If [A] i1i2i3 is of the form [A] i1i2i3 = x (i1+i2+i3) , equation
I π2 + I π3 -2}, [R (π) k ] ij = x (i+k)
does not depend on j.

Reordered tensor modes.

Below we introduce the type-1 and type-2 reordered tensor modes, formed by concatenating the type-1 and type-2 oblique sub-matrices.

Definition 3.6.

The type-1 reordered tensor modes are defined by concatenating the type-1 oblique sub-matrices:

• A ′ 1 is the I 1 × (I 2 I 3 ) matrix [T (1,3,2) I2-1 , . . . , T (1,3,2) 0 = T (1,2,3) 0 , . . . , T (1,2,3)
I3-1 ],
• the weighting factors d

(s) k = 1 if 0 ≤ k < I 2 if I ≤ k < J .
In this way, the cost of the (rank-truncated) HOSVD is reduced to that of the (rank-truncated) SVD of H s D s , which is 2M I 3 flops per iteration. In particular, it can be noted that the compression and weighting of the modes lead to a complexity 6 times as low as that of the algorithm in table 2.1.

Algorithms for Hankel tensors. In the case of (I 1 ×I 2 ×I 3 ) Hankel tensors, [A] i1i2i3 is of the form [A] i1i2i3 = x (i1+i2+i3) , and we proved in point 3 of proposition 3.5 that for all permutation π and index k ∈ {0 . . .

I π2 + I π3 -2}, [R (π) k ] ij = x (i+k) .
In particular, all columns of the type-2 oblique sub-matrix R (π) k are equal. Therefore for each s-mode we can define

• the non-redundant Hankel matrix

H s (i, k) = x (i+k) , of dimension I s ×J s with J s = ( s ′ =s I s ′ ) -1;
• the weighting factors d

(s) k = J (2) {π s ′ } s ′ =s (k) = min({I s ′ } s ′ =s , 1 + k, J s -k) (here d (s)
k is the number of columns of the k th oblique sub-matrix of the s-mode, defined in equation (3.5)). It can be noted that the weighting function 1+k → d The fast SVD-based algorithm for computing the full or rank-truncated HOSVD of the Hankel tensor A is summarized in table 4.1. The compression and weighting of

d (s) k =        1 + k if 1 ≤ 1 + k < min({I s ′ } s ′ =s ) min({I s ′ } s ′ =s ) if min({I s ′ } s ′ =s ) ≤ 1 + k ≤ max({I s ′ } s ′ =s ) J s -k if max({I s ′ } s ′ =s ) < 1 + k ≤ J s 0 elsewhere. (4.1) 6
- max({I s ′ } s ′ =s )) min({I s ′ } s ′ =s )) Js min({I s ′ } s ′ =s ) 1 d (s) k 1

Operation

Cost per iteration SVD of

H 1 D 1 4M 1 I 1 (I 2 + I 3 ) SVD of H 2 D 2 4M 2 I 2 (I 1 + I 3 ) SVD of H 3 D 3 4M 3 I 3 (I 1 + I 2 ) Global cost 24M I 2
the modes allow a reduction of the complexity of one order of magnitude w.r.t. the algorithm in table 2.1. If additionally the Hankel tensor is cube (I 1 = I 2 = I 3 = I), then it is symmetric, and the three modes are equal. In this case, the global complexity is reduced 8M I 2 flops.

Operation

Cost per iteration SVD of

A ′ 1 2M1(90I 2 log 2 (I) + M1I1) SVD of A ′ 2 2M2(90I 2 log 2 (I) + M2I2) SVD of A ′ 3 2M3(90I 2 log 2 (I) + M3I3) Global cost 6M 90I 2 log 2 (I)
4.2. Algorithms exploiting the Toeplitz or Hankel structure. In the above developments, we assumed that the rank-r rank-truncated SVD of an n × m matrix with n < m was computed by means of the orthogonal iteration method [7, pp. 410-411], which consists in recursively performing 2r matrix / vector products and 1 QR factorization of an n × r matrix (a full SVD corresponds to the case r = n). We mentioned that the computational cost of one iteration is 2r c(n, m) + 2r 2 n flops, where c(n, m) is the cost of 1 matrix / vector product, and 2r 2 n is the cost of 1 QR factorization [7, pp. 231-232].

In the following, we will focus on the HOSVD of Toeplitz or Hankel tensors, which can be computed efficiently, using fast matrix / vector products. Indeed, the computational cost of a product between a p × q Toeplitz or Hankel matrix and a vector can be reduced from 2pq flops to 15(p + q) log 2 (p + q) flops, by means of Fast Fourier Transforms (FFT) [7, pp. 188-191,201-202].

Algorithms for Toeplitz tensors. In the case of Toeplitz tensors, we mentioned in point 2 of proposition 3.5 that for all permutation π and index k ∈ {0 . . . I π3 -1}, the type-1 oblique sub-matrix T (π) k is Toeplitz. Therefore the oblique modes A ′ s are formed of Toeplitz blocks. As a consequence, the computational cost of the multiplication of A ′ s by a vector of appropriate dimension can be reduced from 2I 3 flops to 90I 2 log 2 (I) flops 2 . By introducing those fast products into the orthogonal iteration method, the cost of the (rank-truncated) SVD of A ′ s is reduced to 2M s (90I 2 log 2 (I) + M s I s ) per iteration.

The fast algorithm for computing the full or rank-truncated HOSVD of a Toeplitz tensor is summarized in table 4.2. If additionally the tensor A is symmetric, then the three modes are equal. Moreover, as shown in section 4.1, the SVD of A ′ s can be replaced by that of H s D s , where the I × I(I+1) 2 matrix H s is also block-Toeplitz. Therefore the cost of the SVD of H s D s is half that of the SVD of A ′ s . As a consequence, the compression and weighting of the modes lead to a complexity 6 times as low as that of the fast HOSVD algorithm in table 4.2.

Algorithms for Hankel tensors. In the case of Hankel tensors, we noted in section 4.1 that the HOSVD could be obtained by computing the SVD of the matrices H s D s , where each compressed mode H s is a Hankel matrix (H s (i, k) = x (i+k) ). Therefore we can again use fast matrix-vector products to further reduce the complexity. More precisely, the computational cost of the multiplication of the I s ×(( s ′ =s I s ′ )-1) Hankel matrix H s by a vector of appropriate dimension can be reduced from 4I 2 flops to 45I log 2 (I) flops, by means of FFT's. By introducing those fast products into the orthogonal iteration method, the cost of the SVD of H s D s is reduced to 2M s (45I log 2 (I) + M s I s ) per iteration.

The ultra-fast algorithm for computing the full or rank-truncated HOSVD of a Hankel tensor is summarized in table 4.3. Its global cost is provided as a maximum over M 1 , M 2 , M 3 , under the constraint M 1 + M 2 + M 3 = 3M . It can be noted that the cost due to the fast matrix / vector products, and the cost due to the QR factorizations can be of the same order of magnitude if M = O(log 2 (I)). If additionally the Hankel tensor is cube (I 1 = I 2 = I 3 = I), then it is symmetric, and the three modes are equal. In this case, the global complexity is three times as low as that of the ultra-fast HOSVD algorithm in table 4.3.

Comparison of the complexities.

The overall costs of the various HOSVD algorithms presented above are summarized in table 4.4 (sorted in decreasing order of complexity). Note that only the complexity upper bounds are given in this table, and that the calculation of the tensor S is not included. Besides, it can be noted that the FFT-based HOSVD algorithms are not always the fastest, because of the high constants in table 4.4. The best choice for computing the HOSVD actually depends on I, and possibly on M (the dominant cost of all algorithms is linear w.r.t M , except that of the ultra-fast algorithms for Hankel tensors). Fig. 4.2 represents the different complexities for M = 10. From this figure we can draw general remarks, which actually stand for any value of parameter M :

• the best algorithm for computing the HOSVD of Toeplitz (resp. symmetric Toeplitz) tensors is that dedicated to such tensors if I 400, or that dedicated to unstructured (resp. symmetric) tensors otherwise; • the best algorithms for computing the HOSVD of symmetric, Hankel and cube Hankel tensors are always those dedicated to such tensors. In other respects, the comparison between the fast and ultra-fast computations of the HOSVD for Hankel and cube Hankel tensors are sensitive to parameter M , as can be noted in table 4.4. Our simulations showed that:

• for small values of M (M ≪ I), the ultra-fast algorithm is faster if I 70;

• for moderate values of M (M ≃ I/2), the ultra-fast algorithm is faster if I 80; • for large values of M (M ≃ I), the ultra-fast algorithm is faster if I 100.

Conclusions.

In this paper, we proposed to decrease the computational cost of the full or rank-truncated HOSVD, which is basically O(M I 3 ), by exploiting the structure of symmetric, Toeplitz, and Hankel tensors. For symmetric and Hankel tensors, our solution is based on the fact that the HOSVD can be reduced to the SVD of three non-redundant (no column are repeated) matrices whose columns are

i

  , are the s-mode singular values of A and the columns of U (s) are the s-mode singular factors.

Fig. 3 . 1 . 1 .

 311 Fig. 3.1. Type-1 and type-2 oblique sub-matrices of a tensor

  (s) k (plotted in Fig. 4.1) is piecewise linear:

+ kFig. 4 . 1 .

 41 Fig. 4.1. Weighting function d (s) k for Hankel tensors

2

  Under the constraint I 1 +I 2 +I 3 = 3I, the maximum cost is obtained for cube tensors (I 1 = I 2 = I 3 = I). Besides, left or right multiplying an I ×k oblique sub-matrix T (π) I-k by a vector of appropriate dimension normally involves 2Ik flops. This complexity is reduced to 15(I + k) log 2 (I + k) flops by means of FFT's. Therefore left or right multiplying the block-Toeplitz matrix A ′ s by a vector of appropriate dimension normally involves 2 I-1 k=0 2Ik ∼ 2I 3 flops, or 2 I-1 k=0 15(I + k) log 2 (I + k) ∼ 90I 2 log 2 (I) by means of FFT's.

(

  the cost corresponds to a single iteration of the orthogonal iteration method) Operation Cost per iteration SVD of H1D1 2M1(45I log 2 (I) + M1I1) SVD of H2D2 2M2(45I log 2 (I) + M2I2) SVD of H3D3 2M3(45I log 2 (I) + M3I3) Global cost 6M (45I log 2 (I) + M I)

(Fig. 4 . 2 .

 42 Fig. 4.2. Flops count vs. size I for M = 10.

  Tucker s-mode product[START_REF] Lathauwer | Signal Processing based on Multilinear Algebra[END_REF], U (s) is an unitary I s × I s matrix and S is an all-orthogonal and ordered I 1 × I 2 × I 3 tensor. All-orthogonality means that the matrices S is=α , obtained by fixing the s th index to α, are mutually orthogonal with respect to (w.r.t.) the standard inner product. Ordering means that S is=0

	S is=1	. . .	S is=Is-1 for all possible values of s. The Frobenius-norms S is=i ,
	symbolized by σ

Table 2 .

 2 1 HOSVD Algorithm for unstructured tensors(the cost corresponds to a single iteration of the orthogonal iteration method)

Table 4 . 1

 41 Fast HOSVD algorithms for Hankel tensors (the cost corresponds to a single iteration of the orthogonal iteration method)

Table 4 . 2

 42 Fast HOSVD algorithm for Toeplitz tensors (the cost corresponds to a single iteration of the orthogonal iteration method)

Table 4 . 3

 43 Ultra-fast HOSVD algorithm for Hankel tensors

Table 4 . 4

 44 Complexities of the HOSVD algorithms

• A ′ 2 is the I 2 × (I 3 I 1 ) matrix [T (2,1,3)

I3-1 , . . . , T

(2,1,3) 0 = T

(2,3,1) 0 , . . . , T

(2,3,1)

I1-1 , . . . , T

(3,2,1) 0 = T (3,1,2) 0 , . . . , T

I2-1 ]. In the same way, the type-2 reordered tensor modes are defined by concatenating the type-2 oblique sub-matrices:

I1+I2-2 ]. Proposition 3.7.

1. For all s = 1, 2, 3, the mode A s and the reordered modes A ′ s , A ′′ s admit the same singular values and left singular vectors. 2. If A is a symmetric tensor, then

, and this unique mode admits an axial blockwise symmetry w.r.t. its central oblique sub-matrix. Proof.

1. For all s = 1, 2, 3, the columns of the reordered modes A ′ s and A ′′ s form a permutation of the columns of the mode A s defined in section 2. 2. This is a corollary of point 2 in proposition 3.5.

Fast algorithms for computing the HOSVD of structured tensors.

In this section, the reordered tensor modes introduced above are used to efficiently compute the HOSVD of structured tensors. The first improvement consists in exploiting the column-redundancy of symmetric and Hankel tensors. To further reduce the computational cost, we then exploit the fast matrix-vector product techniques specific to Toeplitz and Hankel matrices.

Algorithms exploiting column-redundancy.

Here we suppose that the s-mode of tensor A is redundant, e.g. some columns of the s-mode are equal (this is the case of symmetric and Hankel tensors for instance). We aim at exploiting this redundancy in order to efficiently implement the HOSVD of A. Toward this end, we define the I s ×J s matrix H s as the matrix obtained by removing the repeated columns in the s-mode (J s ≤ s ′ =s I s ′ ), and we denote d (s) k the number of occurrences of the k th column of H s in the s-mode. Then we consider the I s × I s correlation matrix of the s-mode: C (s) = A s A s H . It is clear that this matrix can be factorized as

where

(if the s-mode is not redundant, then we define H s as the s-mode itself and D s is defined as the J s × J s identity matrix). As a consequence, the M s highest singular values and left singular vectors of the s-mode of dimensions I s × s ′ =s I s ′ are the same as those of the smaller I s × J s matrix H s D s . Algorithms for symmetric tensors. In the case of (I × I × I) symmetric tensors, we proved in point 2 of proposition 3.7 that A ′ 1 = A ′ 2 = A ′ 3 , and that this unique mode admits an axial blockwise symmetry. Therefore we can define

• the non-redundant matrix H s = [T ], ∀s ∈ {1, 2, 3}, of dimension I × J with J = I(I + 1)/2;