Keywords: Large-scale graph visualization, matrix-based representation, node-link diagrams, multi-scale navigation, matrix ordering, graph aggregation Visualization de graphe de grande taille, représentation matricielle, diagramme noeud et lien, navigation multi-échelle, réordonnancement de matrice, agrégation de graphe

This article presents the Zoomable Adjacency Matrix Explorer (ZAME), a visualization tool for exploring networks at a scale of millions of nodes and tens of millions of edges. ZAME presents an adjacency matrix graph representation aggregated at multiple scales. It allows analysts to explore a graph at many levels, zooming and panning with interactive performance from the most summary to the most detailed views.

Several components work together to make this performance possible and the results meaningful. A "pyramid" of aggregated views paged on demand to OpenGL GPU shader programs supports smooth multiscale browsing in huge datasets. Efficient matrix ordering algorithms group related elements to make the views meaningful. Using ZAME, we can explore the entire French Wikipedia, over 500,000 articles and 6,000,000 links, with interactive performance on standard consumer-level computer hardware.

Introduction

Few Internet users today are unaware of the Wikipedia online encyclopedia project. Written exclusively by volunteers on the Internet, Wikipedia has become the ninth most popular site on the Web since its inception in 2001. As of this writing, its versions in 250 languages contain more than six million articles, over 1.6 million in the English version alone, making it the largest encyclopedia ever assembled. Most users of Wikipedia today read individual articles or series of articles as in a traditional encyclopedia. However, some users would like to view and analyze the collection as a whole, such as Wikipedia administrators or sociology researchers who want to study its structure and connectivity.

If this can be done easily, more-typical Wikipedia readers might be able to profit from a graphical overview showing where the articles they have read fit in a grouping of related articles. These overviews might include visual clues about the character of the article to help them decide where to go next, such as article length, number of accesses or citations, update activity, authorship, or perhaps rating by some future online editorial agency that the reader trusts. Such tools could enable whole new patterns of exploration that many people will use to study Wikipedia.

This pattern is not limited to Wikipedia, but is part of two general trends. Increasingly, analysts in many fields need to explore and understand large, complex networks with millions of vertices and millions or billion of edges. Also, where appropriate tools have been developed, ordinary users with personal computers today explore datasets on grand scalesthe geography of our entire planet, or all the stocks of the New York exchange-that until recently were the domain of specialists.

For example, we are collaborating with sociologists seeking to understand the dynamics of large online social networks. The social network of the Wikipedia project links tens of thousands of writers to hundreds of thousands or millions of articles; open-source software projects such as the Debian Linux distribution link hundreds of developers to thousands of software modules. Currently they have only statistical tools to test their models of these networks, because interactively exploring such large networks remains a challenge for existing information visualization tools.

In this article, we present the Zoomable Adjacency Matrix Explorer (ZAME), the first system that permits meaningful visual navigation of the whole of Wikipedia, from an overview of the entire corpus down to the individual articles, seamlessly and with interactive performance. ZAME is based on a zoomable multiscale adjacency matrix representation of Wikipedia articles and their internal links. It implements three original features: a fast automatic reordering mechanism to find a good layout; special indices supporting a rich variety of data aggregations; and GPU-accelerated rendering with programmble shaders to deliver interactive and smooth framerates.

ZAME can be used with any network including social networks and bioinformatics networks. The system provides two different scales of zoom-geometric zoom as well as detail zoom-allowing the user to independently change both the amount of data visible on the screen as well as the size allocation of each data point. In addition, we use special visual representations of aggregated edges implemented in shader programs to make the aggregation (and its contents, if desired) evident to the user.

The remainder of this article is organized as follows: we first present a survey of related work. This is followed by a description of ZAME's features and design. We conclude with a discussion of its current implementation and performance results.

INRIA

Related Work

The most crucial attribute of a graph visualization is its readability, its effectiveness at conveying the information required by the user tasks [START_REF] Ghoniem | On the readability of graphs using nodelink and matrix-based representations: a controlled experiment and statistical analysis[END_REF]. For example, Ghoniem et al. proposed a taxonomy of graph tasks and evaluated the effectiveness of node-link and adjacency matrix graph representations to support them in [START_REF] Plaisant | Task taxonomy for graph visualization[END_REF]. They found that graph visualization readability depends strongly on the graph's size (number of nodes and edges) and density (average edges per node). Different visualizations had better readability for different graph sizes and densities, as well as for different tasks.

Several recent efforts visualize large graphs with aggregated representations that present node-link and/or matrix graphs at multiple levels of aggregation. The following sections describe several important ones.

Node-Link Visualization of Large Networks

Node-link diagrams can effectively visualize networks of about one million vertices if they are relatively sparse. Hachul and Jünger [START_REF] Hachul | An experimental comparison of fast algorithms for drawing general large graphs[END_REF] compared six large-scale graph drawing algorithms for 29 example graphs, some synthetic and some real, of which the largest had 143,437 vertices and 409,593 edges. Of the six algorithms, only three scaled well: HDE [START_REF] Harel | Graph drawing by high-dimensional embedding[END_REF], FM 3 [START_REF] Hachul | Drawing large graphs with a potential-field-based multilevel algorithm (extended abstract)[END_REF] and to some extent GRIP [START_REF] Gajer | Grip: Graph drawing with intelligent placement[END_REF]. However, the densities of the sample graphs were small, typically less than 4.0. When the density or the size grows, some dimensionality reduction mechanism is needed to maintain readablity.

Hierarchical aggregation allows larger graphs to be visualized and navigated, assuming that there is an algorithm for finding suitable aggregations at each level in a reasonable time. An aggregation is suitable if each aggregated level can be visualized effectively as a node-link diagram, and if navigation between levels has sufficient visual continuity for users to maintain their mental map of the whole network and avoid getting lost. No automated strategy published to date can select an appropriate algorithm for an arbitrary network, but there are many successful aggregation algorithms for specific categories of graphs. For example, Auber et al. [START_REF] Auber | Multiscale visualization of small world networks[END_REF] present effective algorithms for aggregating and visualizing the important class of networks (including the global Internet and many social networks) known as small-world networks [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF], whose characteristics include power-law degree distribution, high clustering coefficient and small diameter. Systems such as Tulip offer multiple clustering algorithms and are designed to permit smooth navigation on large aggregated networks [START_REF] Auber | Tulip : A huge graph visualisation framework[END_REF].

Gansner et al. propose another method involving a topological fisheye that is usable when a correct 2D layout can be computed on a large graph [START_REF] Gansner | Topological fisheye views for visualizing large graphs[END_REF]. After the network is laid out, it is topologically simplified to reduce the level of detail at coarser resolutions. The fisheye selects a focus node and displays the full detailed network around it, showing the remainder of the network at increasingly coarse resolutions for nodes farther away from the focus. This technique preserves users' mental maps of the whole graph while magnifying (distorting) the network around the focus point. However, effectively laying out an arbitrary large, dense graph remains an open problem for node-link diagrams. All the methods require a good global initial layout, which can be very expensive to compute.

Wattenberg [START_REF] Wattenberg | Visual exploration of multivariate graphs[END_REF] describes a method for aggregating networks according to attributes on their vertices. The aggregation is only computed according to the attribute values, much like pivot tables in spreadsheet calculators or data cubes in OLAP databases. This approach works best when the values are categorical or numerical with a low cardinality. The article only refers to categorical attributes on the vertices. When no categorical attribute is suitable for computing the pivots, as with the Wikipedia hypertext link network visualized in this paper, this approach is not effective-the whole network would be displayed as a single point.

Overall, the node-link representation has two major weaknesses [START_REF] Ghoniem | On the readability of graphs using nodelink and matrix-based representations: a controlled experiment and statistical analysis[END_REF]: (i) it copes poorly with dense networks, and (ii) without a good layout, it requires aggregation methods to reduce the density enough to be readable. Because these methods are very dataset dependent, current node-link visualization systems leave the choice of aggregation and layout to the users, who therefore need considerable knowledge and experience to get good results.

Matrix Visualization of Large Networks

Several recent articles have used the alternative adjacency matrix representation to visualize large networks. Abello and van Ham demonstrated [START_REF] Abello | MatrixZoom: A visual interface to semi-external graphs[END_REF] the effectiveness of the matrix representation coupled with a hierarchical aggregation mechanism to visualize and navigate in networks too large to fit in main memory. Their approach is based on the computation of a hierarchy on the network displayed as a tree in the rows and columns of the aggregated matrix representation. The aggregation is a hierarchical clustering in which items are grouped but not ordered. It is computed according to memory constraints as well as semantic ones. The users operate on the tree to navigate and understand the portion of the network they are currently viewing.

Navigation in this approach is constrained by the hierarchical aggregation: users navigate in a tree that has been computed in advance. The main challenge is to find an aggregation algorithm that is both fast and that produces a hierarchy meaningful to the user. Unfortunately, this choice is typically dataset-dependent. Without a meaningful hierarchy, users navigate in clusters containing unrelated entries and cannot make sense of what they see.

Henry and Fekete in [START_REF] Henry | MatrixExplorer: a dual-representation system to explore social networks[END_REF] proposed several methods based on reordering the rows and columns of the matrix as opposed to just clustering similar nodes. They describe two methods based on approximate Traveling Salesman Problem (TSP) solutions, which are computed on the similarity of connection patterns and not on the network itself. One method uses a TSP solver directly; the other initially computes a hierarchical clustering and then reorders the leaves using a constrained TSP. Both algorithms yield orderings that reveal clusters, outliers and interesting visual patterns based on the roles of vertices and edges. Because the matrix is ordered, not just clustered, both navigation (panning) between and within clusters of similar items may reveal useful structure. Unfortunately, both reordering methods are at best quadratic; since they need to compute the full distance matrix between all the vertices. Therefore, it is difficult to scale them to hundreds of thousands or millions of vertices. 3 The Zoomable Adjacency Matrix Explorer ZAME is a network visualization system based on a multiscale adjacency matrix representation. Attributes associated with the vertices and edges can be mapped to visual attributes of the visualization, such as color, transparency, labels, border width, size, etc., using selectable schemes to aggregate each attribute appropriately at higher levels. ZAME integrates all the views of a large network, from the most general overview down to the details, and provides multiscale navigation techniques that operate at all levels. Therefore, it can be used for tasks at all levels, from understanding the network's overall structure, to exploring the distribution of results to a content-based search query at multiple levels of aggregation, to performing analysis tasks such as finding cliques or the most central vertices. Additional relevant graph tasks can be found in [START_REF] Plaisant | Task taxonomy for graph visualization[END_REF].

INRIA

Interactive System

Representation

The main technical challenge of the system is managing the huge scale of such graphson the order of a million vertices for the Wikipedia database-while delivering the real-time framerates necessary for smooth interaction. Furthermore, the graph must be laid out (the matrix nodes reordered) to group similar nodes so that the visualization becomes readable, i.e. in such a way that patterns emerge and conclusions can be drawn from the data.

To achieve these goals, our tool consists of four different components (see Figure 2): a hierarchical data structure for storing the graph and its attributes in an aggregated format;

an automatic reordering component for computing a meaningful order for the matrix to support visual analysis;

an accelerated rendering mechanism for efficiently displaying and caching a massive graph dataset; and a set of navigation techniques for exploring the graph.

The following sections describe these components in detail.

Multiscale Data Aggregation

To support multiscale graph exploration with interactive response, we designed an index structure tailored to its visualization abstraction, the 3-dimensional binary pyramid of "detail levels" shown in Figure 3. Users pan across the surface of one detail level; they zoom up and down between detail levels (along perspective lines meeting at the tip of the pyramid, so features match between zoom levels). Detail level zero of this abstraction, the bottom level of the pyramid, is the adjacency matrix of the raw data with the nodes arranged (according to the reordering permutation) on the rows and columns at the edge and the edges between them indicated at the row/column intersections on the level's surface. Every detail level above the base has half the length and width (number of nodes) and a quarter the intersection squares (possible edges) as the level below it. Therefore, each node at a higher, more summary detail level represents two nodes at the level below it (except the last node on a level, in the case where the lower level had an odd number of nodes), and each intersection square indicating a possible edge at this level represents four possible edges at the level below it. Combining nodes and edges for higher detail levels require that we define how to aggregate the underlying data and its corresponding visual representation. See Sections 3.3 and 3.5 for these definitions. The original network at detail level 0 is stored as a standard graph in the InfoVis Toolkit (IVTK) [START_REF] Fekete | The InfoVis Toolkit[END_REF]. The specialized index structure is maintained in "zoomable" equivalents specialized for the pyramid index structure.

A standard IVTK graph is implemented by a pair of tables, one for its vertices and another for its edges. Each row represents one vertex or edge. It contains attributes in "internal columns" that maintain the graph topology. For each vertex, there is a first and last edge for each of two linked edge lists for its outgoing and incoming edges. Each of these lists is maintained in a pair of columns in the edge table: the "next edge" and "previous edge" for the outgoing and incoming edge lists. To complete the topology, two more columns of the edge table store its first and second vertex (source and sink edge vertices for directed graphs). The doubly-linked edge lists are an optimization for fast edge removal, but the back-link can be omitted to reduce memory consumption if necessary. Thus, the IVTK needs to store four numbers for each vertex, and either six or four numbers for each edge. Its current implementation uses CERN COLT [START_REF]The Colt project[END_REF] primitive integer arrays that have very little overhead; so in a 32-bit Java implementation today, the memory consumed is 16 bytes for each vertex and 16-24 bytes for each edge.

A zoomable graph has the same basic structure. However, instead of single-layer vertex and edge tables, it uses specialized "zoomable" tables with multilevel indices. Also, it maintains some invariants that accelerate important operations. Below is a description of the basic operations: iterator(int level) returns an iterator over all the elements of a specified aggregation level.

The zoomable vertex table refers to the original graph's vertex table (its related table.) Each aggregated vertex at level 0 refers to one vertex in the related table, except that their order is changed by a "permutation" structure that implements the reordering described in 3.2. The numbering of aggregated vertices at level 1 begins immediately after those at level 0, 2 after 1, and so forth (except that an odd number at any level is rounded up by one). Each pair of vertices at at level n is aggregated by a single vertex at level n + 1. So, given the number of a vertex at any level, this simple numbering scheme makes it straightforward to compute corresponding vertices at levels above and below it. It is also straightforward to calculate the size of the index: it is bounded by the series 1/2+1/4+. . . < 1; so all the index levels at most double the size of the table. Like vertices at level 0, an aggregated vertex has pointers (edge numbers) to its lists of out-and in-edges, but these refer to aggregated edges.

Unfortunately, the zoomable edge table index is much more complicated to build and maintain, because the set of possible edges is not dense. The number of possible edges between N nodes is on the order of N 2 ; so in 32-bit signed arithmetic, calculations based on a simple enumeration of possible edges analogous to those used for vertex indices would overflow for more than 2 15 (32K) vertices. As with the vertex table, edges at level n follow edges at level n -1. However, because the corresponding edge numbers at each level cannot be calculated, they must be stored. Achieving fast access requires several optimizations. To compute the level of an edge, we perform a binary search in a vector containing the starting index of each level. At each level, the outgoing edges are stored in order by their first (source) vertex: the outgoing edges of vertex n follow the outgoing edges of vertex n-1. Similarly, all the outgoing edges for one source vertex are sorted in order of their second (sink) vertex. This arrangment makes it very fast to search for an edge given its vertices using two levels of binary search. This optimization allowed us to omit the previous and next edge columns for aggregated edges.

Despite memory optimizations such as using binary search in edge lists ordered at two levels to replace next-edge and previous-edge lists, aggregated edge indices are still very costly in terms of memory, several times larger than the original data. We still need to maintain a separate incoming edge list. We also need an extra column to store the "super edge" (that is, the corresponding edge in the level above). This appears wasteful; since the "super vertex" corresponding with this endpoint could be calculated and the super edge found by a binary search in the edge list of the super vertex. Wikipedia has around six million edges; so 24 Mb are required just to add this one column to the base level, and the total including all its index levels is around five times more. However, we still chose to store this information, because the basic operation of aggregating edge attributes requires going through each edge from level 0 up and accumulating the aggregated results on the super edge. Without direct access, the complexity of this operation would be n × log(n) instead of n, and it would require tens of minutes instead of minutes.

Because the amount of information used for aggregated indices on a huge file such as Wikipedia exceeds the virtual memory capacity of a 32-bit Java Virtual Machine (JVM), we implemented a paging mechanism that allocates columns of memory in fairly large fixed-size chunks that are retrieved from disk when needed. Fortunately, the memory layout of the zoomable aggregated graph is very well suited to paging. Most operations are performed in vertex-and edge-order on a specified level; so they tend to use consecutive indexes likely to be allocated nearby on disk.

The total size of the aggregated edge table depends dramatically on the quality of the ordering. A good ordering groups both edges and non-edges; so multiple nearby edges and non-edges aggregate at each step and the size of successive index levels rapidly diminishes. The worst case is a "pepper and salt" pattern of widely-spread edges that do not aggregate INRIA significantly for many levels, resulting in an aggregated edge table that can be 4 to 8 times larger than the original edge table. Because Wikipedia averages only 12 links per page, without a good ordering it has this kind of profile. Our current reordering methods, though imperfect, improve this to about a factor of 5.

Reordering

Of all the algorithms described in the literature on matrix reordering, few are sub-quadratic. We experimented with those based on linear dimension reduction such as Principal Component Analysis (PCA) or Correspondance Analysis (CA) and greedy TSP algorithms such as the Nearest-Neighbor TSP (NNTSP) heuristic.

High-Dimensional Embedding

PCA and Correspondence Analysis were used effectively by Chauchat and Risson [START_REF] Chauchat | a new method and a software integrating Jacques BERTIN's Graphics and Multidimensional Data Analysis Methods[END_REF] for reordering matrices. Their matrices were small enough to permit computing the eigenvectors directly, but this is obviously infeasible for hundreds of thousands or millions of points. Harel and Koren describe a modified PCA method called "High Dimensional Embedding" [START_REF] Harel | Graph drawing by high-dimensional embedding[END_REF] that can efficiently lay out very large graphs by computing only k-dimensional vectors where k is typically 50. This method is designed for laying out node-link diagrams by using the two or three first components of the PCA for positioning of the vertices. We used it and improved it for reordering matrices.

The solution proposed by Harel and Koren consists of choosing a set of pivot vertices that the algorithm tries to place near the outer edges of the graph, and to use the graph distances to these pivots as the coordinates of each vertex. When 50 pivots are chosen, these coordinates are in 50 dimensions. PCA is computed on these dimensions and the eigenvectors are computed using a power-iteration that converges very quickly in real cases. Their main algorithm is listed in Figure 9. The BFS algorithm is a simple Breadth First Search computed when no edge weights exist; otherwise, it is replaced by a Dijkstra shortest path computation.

Figure 4) compares the results of HDE versus TSP on a relatively small social network dataset. Although HDE's results (right side) are noticeably worse, the calculation is much faster, requiring just a fraction of second instead of 30 seconds required by the standard TSP solver that produced the figure on the left. Although the TSP solver clearly could not be applied to the far larger Wikipedia dataset, HDE reordering is quite fast, but again, the local quality of the ordering is poor.

More seriously, the original algorithm can choose pivots in a very ineffective way. Consider a large connected network with two very distant star-shaped vertices V 1 and V 2 , much more distant than any other vertices except the children of these two vertices. HDE will pick a random pivot first, then take the farthest vertex from the pivot, say around V 1 . The next vertex will then be a pivot around V 2 . In turn, the next will be another vertex around V 1 etc. until all the children of V 1 and V 2 are enumerated, producing a very biased distribution of pivots (see Figure 5). This is not merely a theoretical problem; something similar actually RR n 6163 To avoid this bias, we modified HDE to penalize edges progressively according to the number of times they already participate in a path between existing pivots (Figure 9). This penalty encourages the algorithm to place new pivots in regions of the graph not traversed INRIA by paths between previous pivots, which hopefully represent important different features of its overall structure.

This algorithm avoids the pathological distribution of pivots, but it does not solve the fundamental question of how many pivots are required to adequately represent a large graph, which is still an open question.

Nearest-Neighbor TSP Approximation

Although the Traveling Salesman Problem is NP-complete, it has good approximations in many cases. When the edges are weighted and the weight values are not too similar, the Nearest-Neighbor TSP (NNTSP) approximation algorithms can be can be effective [START_REF] Gutin | The Traveling Salesman Problem and Its Variations[END_REF]. An initial ordering can be computed in linear time by limiting the search distance. In our NNTSP algorithm we limit the distance to 6 since the number of nodes examined grows by a factor of 12 (the average number of links per page in Wikipedia) on each iteration, and we felt that pages more than 6 hops away were unlikely to be good matches.

Because hypertext links are unweighted, we computed various dissimilarity functions between the source and destination pages of the link to provide them.

To date, we have tried three link-weighting or distance-calculation algorithms. In order of increasing quality and running time, they are dissimilarity of adjacency patterns, distance in the HDE pivot space and distance in a randomly-selected pivot space. (Dis)similarity of adjacency patterns is determined by merging the set of destination vertices between the edge sets of a pair of source vertices. Visually, an ordering that optimizes this measure will tend to maximize the number of groups of aligned squares on an adjacency matrix representation, which is good for perceptual analysis. This measure can be calculated efficiently-in time linear in the average number of vertices per node-so it was the first one we tried.

NNTSP using this distance measure does a reasonable job of aligning almost all the vertices of Wikipedia-about 97%, approximately the first 460,000 of the 474,594 pages in its largest connected component that we actually analyze (since graph distance cannot be computed for unconnected components). However, at the end of the NNTSP algorithm, when relatively few unvisited vertices remain to be placed, the chance of finding pages with common links is small. Therefore, since the adjacency pattern measure is likely to be zero for all of the remaining pages, even when it is computed for neighbors of neighbors of neighbors . . . many cannot discriminate between them and the search for a match becomes very long.

The second and third methods we tried differ from the first by the adjacency pattern dissimilarity metric. First, they find some number (currently 10) of HDE pivots, respectively either by Harel and Koren's method [START_REF] Harel | Graph drawing by high-dimensional embedding[END_REF] or randomly. They compute the distances between every graph point and each pivot, which takes about 20 seconds per pivot for Wikipedia. Then, each graph distance is approximated by the Triangle Inequality: the true distance between vertices is less than the sum of the distances between the two vertces and any pivot. The shortest distance is used (corresponding to the closest pivot). Dissimilarity is computed on the distance-to-pivots coordinate space and nearest-neighbor is also computed in that space. These methods find better orderings for the graph than the first (which gives up and accepts a low-quality ordering for the last 3% of the points) at the expense of several minutes longer calculation time.

Results

The HDE method provides the best overview of the graph at a high, heavily-aggregated level, with the majority of connected nodes grouped in the upper left corner. In contrast, the NNTSP algorithm orders the Wikipedia data reasonably well globally, but dramatically better locally. The global order is actually not as good as HDE's, as reflected by a poorer compression factor of the indicesd-NNTSP leaves more "salt and pepper" vertices. However, its ordering is dramatically better at a small scale. Basically, it finds many large groups of closely-related pages, groups them in the order and arranges them reasonably. It is the default reordering method we apply.

We provide some examples of the reordering performance of the three algorithms proposed in this paper to show the difference in results on both global as well as local level. See Figure 10 for overview and detail views of the results of reordering the same matrix using the HDE, NNTSP, and NNTSP2 algorithms.

Aggregating Data Attributes

To make views of data at a higher level meaningful, attributes of the original vertices (such as article names, creation date, or number of edits) and edges (such as link weights) must be combined in meaningful ways, with rapid access to details including the original data values. This requires considerable flexibility in order to accommodate a range of data types and semantics as well as user intentions. There is a new trend in the classification community to use so-called "symbolic data analysis" for richer aggregation [START_REF] Billard | Symbolic Data Analysis: Conceptual Statistics and Data Mining[END_REF]. We summarize the principle and explain how we map symbolic data to visualizations.

It is important to make it evident to the user when a particular cell is representing an aggregated, as opposed to an original, attribute. Section 3.5 describes some of the visual representations we employ for this purpose. Typically, if enough display space is available, a histogram can faithfully visualize the aggregated values for each item. If less space is available, a min/max range or Tukey diagram can be used. In the worst case, where only one or a few pixels are available, the aggregated value can be used to modulate the color shade of the whole cell.

Categorical Attributes

Categorical attributes-including Boolean values-have a cardinality: the number of categories. At the non-aggregated level, a categorical variable can hold one categorical value, such as a US state. When aggregating categorical values, we compute a distribution, i.e. the count of each item aggregated per category.

INRIA

Numerical Attributes

Numerical attributes may be aggregated by various methods such as mean, maximum or minimum. They can also be transformed into categorical values by binning the values in intervals. Numerical data permits a wide range of analysis such as calculating standard deviations and other statistics. ZAME internally computes a discrete distribution for the aggregated values using a bin width computed according to [START_REF] Scott | On optimal and data-based histograms[END_REF]. We also keep track of the mean, extreme, and median values.

Nominal Attributes

Unlike even ordinal data, there is no inherent relationship between nominal attributes such as article names, authors, or subject titles. Unfortunately, they are often vital to understanding what elements of the visual representation refer to. Like numeric attributes, nominal attributes of specific datasets can be aggregated using special methods such as concatenation, finding common words, or sampling representative labels. Methods such as excentric labels [START_REF] Fekete | Excentric labeling: dynamic neighborhood labeling for data visualization[END_REF] may be used to display many more labels on demand.

Visualization

For the visualization component of our matrix navigation tool, we are given an elusive challenge: to efficiently render a matrix representation of a large-scale aggregated graph structure consisting of thousands if not millions of nodes and edges. The rendering needs to be efficient enough to afford interactive exploration of the graph with a minimum of latency.

We can immediately make an important observation: for matrices of this magnitude, it is the screen resolution that imposes a limitation on the amount of visible entities. In other words, there is no point in ever drawing entities that are smaller than a pixel at the current level of geometric zoom. In fact, the user will often want the entities to be a great deal larger than that at any given point in time. This works in our favor and significantly limits the depth we need to traverse into the aggregated graph structure in order to render a single view of the matrix.

At the same time, we must recognize that accessing the aggregated graph structure may be a costly operation and one which is not guaranteed to finish in a timely manner. Clearly, in order to achieve real-time framerates, we must decouple the rendering and interaction from the data storage.

Our system solves this problem by utilizing a tile management component that is responsible for caching individual tiles of the matrix at different detail levels. Missing tiles are temporarily exchanged for coarser tiles of a lower detail level until the correct tile is loaded and available. The scheme even allows for predictive tile loading policies, which try to predict user navigation depending on history and preload tiles that may be needed in the future.

In the following text, we describe these aspects of tile management and predictive loading in more depth. We also describe our use of programmable shaders and textures for efficiently rendering matrix visualizations.

Basic Rendering

Instead of attempting to actively fill and stroke each cell representing an edge in our adjacency matrix, we use 2D textures for storing tiled parts of the matrix in video memory. Textures are well-suited for this purpose since they are regular array structures, just like matrix visualizations. Also, they are accessible to programmable vertex and fragment shaders running on the GPU (graphical processing unit) of modern 3D graphics cards. This allows us to avoid sending excessive geometry data to the GPU and instead render a few large triangles with procedural textures representing the matrix.

Our basic matrix rendering shader accesses the texture information for a given position and discards the fragment if no edge is present. If there is an edge there, the data stored in the texture is used to render a color or visual representation (see Section 3.4.1). Stroking is performed automatically by detecting whenever a pixel belongs to the outline of the edge-in this case, black is drawn instead of the color from the visual representation of the edge.

Tile Management

Adjacency matrices may represent millions of nodes on a side, and thus storing the full matrix in texture memory is impossible. Rather, we conceptually split the full matrix into tiles of a fixed size and focus on providing an efficient texture loading and caching mechanism for individual tiles depending on user navigation.

In our implementation, we preallocate a fixed pool of tiles of a given size. We use an LRU cache to keep track of which tiles are in use and which can be recycled. As the user pans and zooms through the matrix, previously cached tiles can be retrieved from memory and drawn efficiently without further cost. Tiles which are not in the cache must be fetched from the aggregate graph structure-this is done in a background thread that keeps recycling and building new tiles using the cache and the tile pool.

While an unknown tile is being loaded in the background thread, the tile manager uses an imposter tile. Typically, imposter tiles are found by stepping upwards in detail zoom levels until a coarser tile covering the area of the requested tile is eventually found in the cache. The imposter is thus an aggregation of higher zoom levels and therefore not a perfectly correct view of the tile, but it is sufficient until the real tile has finished loading.

Predictive Tile Loading

Beyond responding to direct requests from the rendering, our tile caching mechanism can also attempt to predictively load tiles based on an interaction history over time. For example, if the user is increasing the detail zoom level of the visualization, we may try to preload a number of lower-level tiles in anticipation of the user continuing this operation. Alternatively, if the user is panning in one direction, we may try to preload tiles in this direction to make the interaction smoother and more correct.

INRIA
Our tile manager implementation supports an optional predictive tile loading policy that plugs into the background thread of the tile manager. Depending on the near history of requested tiles, the policy can choose to add additional tile requests to the command queue. Furthermore, the visualization itself can give hints to the policy on the user's interaction.

Aggregated Visual Representations

By employing programmable fragment shaders to render procedural textures representing matrix tiles, we get access to a whole new set of functionality at nearly no extra rendering cost. In our system, we use this capability to render special visual representations for aggregated edges. As indicated in Section 3.3, we can use these to give the user an indication of the data that has been aggregated to form a particular edge.

Currently, we support the following such visual representations (see Figure 7 for examples of these):

Standard color shade. Single color to show occupancy, or a two-color ramp scale to indicate the value.

Average. Computed average value of aggregated edges shown as a "watermark"value in the cell. Each representation has been implemented as a separate fragment shader and can easily be exchanged. Furthermore, new representations can also be added. Depending on the availability and interpretation of the data contained in the tile textures, the user can therefore switch between any of these representations at will and with no performance cost.

Figure 8 shows a general overview of the fragment shaders used in our system. The texture representing the matrix tile is first accessed to see whether there is an edge to draw at all; if not, the fragment is discarded and nothing is drawn. The next step is to check whether the current fragment resides on the outer border of a cell, in which case the fragment is part of the stroke and black is output. Finally, the last step depends on the actual visual representation chosen, and determines the color of the fragment depending on its position in INRIA the cell. The output color can either be the currently active OpenGL color for flat shading, or a ramp color scale indexed using the edge data.

Navigation

In the previous sections, we have used the terms geometric zoom and detail zoom without defining them properly. Navigation techniques for the ZAME system operate on these two properties as well as the visual substrate on which the adjacency matrix is laid out:

Geometric zoom encodes the position and dimensions of the currently visible viewport on the visual substrate.

Detail zoom describes the current level of detail of the adjacency matrix.

In other words, the viewport defined by the geometric zoom governs which part of the matrix is mapped to the physical window on the user's screen. This is a continuous measure. The detail zoom, on the other hand, governs how much detail is shown in the window, i.e. at which discrete level in the hierarchical pyramid structure we are drawing the matrix. Since the hierarchy has discrete aggregation levels, detail zoom is also a discrete measure. ZAME provides all of the basic navigation and interaction techniques of a graph visualization tool (see Table 1 for an overview of navigation techniques in the system). Users can pan around in the visualization by grabbing and dragging the visual canvas itself, or by manipulating the scrollbars.

Operation

Geometric and detail zoom are often coupled using a mapping function from user context to zoom level. A typical navigation technique might for instance be designed so that the amount of visible detail depends on the available display space in the window. In the ZAME system, this mapping is implemented by a zoom policy that can be configured by the user. We currently only support a policy that tries to maintain a desired cell size given the available display space, but we will explore more complex policies tailored to the user's current task in the future.

Beyond this coupled zooming operation, ZAME also supports direct geometric or detail zooming. In this way, the tool provides drill-down functionality to give details-on-demand on any edge or vertex of the matrix visualization.

Dataset

Nodes Edges Load Reorder (secs) (secs) InfoVis04 1,000 1,000 10 30 Protein-protein graph 100,000 1,000,000 10 30 Wikipedia (Fr) 500,000 6,000,000 50 50

Table 2: Performance measurements for standard graph datasets processed and visualized using ZAME.

Results and Discussion

Here we give some brief results on our implementation and performance of the matrix navigation tool.

Implementation

Our implementation is built in Java using only standard libraries and toolkits. Rendering is performed using the JOGL 1.0.0 with OpenGL 2.0 and the OpenGL Shading Language (GLSL). The implementation is built on the InfoVis Toolkit [START_REF] Fekete | The InfoVis Toolkit[END_REF] and will be made publicly available as an extension module to this software.

Performance Measurements

Performance measurements of the different phases of the ZAME system for several graph datasets are presented in Table 2. The measurements were conducted on a Intel Core T M 2, 2.13 GHz computer with 2 GB of RAM and an NVIDIA GeForce FX 7800 graphics card with 128 MB of video memory. For the navigation, the visualization window was maximized at 1680×1200 resolution.

Conclusions and Future Work

This article has described our tool for navigating in massive networks on the scale of millions of nodes and edges. It allows smooth interactive visualization of the link structure of the French Wikipedia, half a million articles connected by six million links. This article describes the technical innovations we used: a fast automatic reordering mechanism to find a good layout; special indices supporting a rich variety of data aggregations; and GPU-accelerated rendering with shader programs to deliver interactive and smooth framerates.

INRIA

In the future, we intend to continue exploring the problem of multiscale navigation and how to provide powerful yet easy-to-use interaction techniques for this task. More specifically, we are interested in exploring the human aspects of detail zoom versus geometric zoom and suitable policies for coupling these together. We are also interested in more closely studying the user utility of the aggregate visual representations introduced in this work.

Figure 1 :

 1 Figure 1: Screenshot of ZAME showing a min-max aggregation of numerical values.

Figure 2 :

 2 Figure 2: Overview of the components of ZAME.

Figure 3 :

 3 Figure 3: Conceptual structure of the aggregated graph data structure.

Figure 4 :

 4 Figure 4: Results of the algorithms TSP (left) and HDE (right) applied to a social network

Figure 5 :

 5 Figure 5: Pathological HDE pivot selection for the Wikipedia dataset.

Figure 6 :

 6 Figure 6: Matrix visualization rendering pipeline.

Figure 7 :Figure 8 :

 78 Figure 7: Seven different visual representations for aggregated edges.

HDE

Figure 10 :INRIAFigure 11 :

 1011 Figure 10: Overview (left) and detail (right) views of reordering using the NNTSP2 algorithm.INRIA

Table 1 :

 1 Overview of ZAME interaction techniques.

		Input	Effect
	zooming	MMB + mouse y	zooming
	panning	RMB + mouse x/y	viewport position
	det-zooming mouse wheel	detail zoom
	geo-zooming shift + MMB + mouse y geo zoom
	selection	LMB + mouse x/y	cell select

RR n 6163

Acknowledgments

This work was supported in part by the French Autograph and SEVEN projects funded by ANR.

Function HighDimDraw(G(V = 1, . . . , n, E), m) % This function finds an m-dimensional layout of G: | % Initialize the penalized edge length E to 1 for each edge

| % Penalize all edges in shortest path from pi to pi+1

2 end for end for return X 1 , X 2 , . . . , X m end for


Vincent Drouard
email: vincent.drouard@inria.fr

Silèye Ba

Radu Horaud

Switching Linear Inverse-Regression Model for Tracking Head Pose *

We propose to estimate the head-pose angles (pitch, yaw, and roll) by simultaneously predicting the pose parameters from observed high-dimensional feature vectors, and tracking these parameters over time. This is achieved by embedding a Gaussian mixture of linear inverse-regression model into a dynamic Bayesian model. The use of a switching Kalman filter (SKF) enables a principled way of carrying out this embedding. The SKF governs the temporal predictive distribution of the pose parameters (modeled as continuous latent variables) conditioned by the discrete variables associated with the mixture of linear inverse-regression formulation. We formally derive the equations of the proposed switching linear regression model, we propose an approximation that is both identifiable and computationally tractable, we design an EM procedure to estimate the SKF parameters in closed-form, and we carry out experiments and comparisons with other methods using recently released datasets.

Introduction

Recent advances in computer vision have demonstrated the relevance of representing images and image regions with feature vectors lying in high-dimensional feature spaces, e.g. SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], SURF [START_REF] Bay | Speeded-up robust features (surf)[END_REF], and any of their variants, or CNN-based features which may be used in conjunction with regression [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF][START_REF] Sharif Razavian | Cnn features off-the-shelf: an astounding baseline for recognition[END_REF] and tracking [START_REF] Ma | Hierarchical convolutional features for visual tracking[END_REF]. The rationale of representing image regions with high-dimensional feature vectors is that the latter supposedly embed hidden information, such as identity or pose. For example, in the case of face analysis, one can infer both face recognition and face orientation from such features. In the case of face orientation, or head pose, the task consists of extracting a low-dimensional parameterization, i.e. pitch, yaw and roll, from the high-dimensional feature space -a parameterization of the head-pose manifold. Not surprisingly, some of the best performing headpose estimation methods rely either on dimensionality reduction followed by regression, [START_REF] Srinivasan | Head pose estimation using view based eigenspaces[END_REF][START_REF] Raytchev | Head pose estimation by nonlinear manifold learning[END_REF][START_REF] Hu | Head pose estimation by non-linear embedding and mapping[END_REF][START_REF] Li | Query driven localized linear discriminant models for head pose estimation[END_REF][START_REF] Benabdelkader | Robust head pose estimation using supervised manifold learning[END_REF][START_REF] Foytik | A two-layer framework for piecewise linear manifold-based head pose estimation[END_REF][START_REF] Sundararajan | Head pose estimation in the wild using approximate view manifolds[END_REF], or on high-dimensional-to-low-dimensional regression, e.g. [START_REF] Murphy-Chutorian | Head pose estimation for driver assistance systems: A robust algorithm and experimental evaluation[END_REF][START_REF] Marin-Jimenez | Detecting people looking at each other in videos[END_REF][START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF][START_REF] Drouard | Head pose estimation via probabilistic highdimensional regression[END_REF].

Nevertheless, these feature-based approaches estimate the head-pose parameters from one image and are not designed to track the parameters over an image sequence. The observed feature vectors contain more than just head pose information, e.g. variabilities in illumination, appearance, shape, identity, background, clutter, etc. Moreover, errors in face localization are inherent, i.e. the bounding-box needed to extract the feature vector is not always correctly aligned with the face itself. All these time-varying phenomena may induce large oscillations and inconsistencies in the estimation of the parameter values, yielding non-smooth temporal trajectories.

In this paper we propose to simultaneously predict headpose parameters from observed feature vectors and to track these parameters over time, based on embedding regression into a dynamic Bayesian model. Without loss of generality, we adopt a HOG-based description of faces, hence we need to predict a low-dimensional output (head-pose parameters) from a high-dimensional input (HOG vectors). To solve the latter we train the regression of [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF] which is a generative Gaussian mixture of linear regression model. The proposed dynamic model is based on the switching Kalman filter (SKF) formulation. The proposed SKF governs the temporal predictive distribution of the pose parameters (which are continuous latent variables) conditioned by the state variables (the discrete latent variables associated with the mixture of linear regression formulation). The rationale of plugging regression into a dynamic Bayesian model is that the latter filters the prediction of the former while taking full advantage of the rich generative regression formulation.

We formally derive a switching dynamic Bayesian model, we devise an approximation of this model that is both identifiable and computationally tractable, we design an EM procedure that estimate the parameters, and we carry out experiments and comparisons with other methods. The principle of the tracker is summarized in Fig. 1 and an example is shown in Fig. 2.

The remainder of the paper is organized as follows. Section 2 discusses the related work. Section 3 summarizes the mixture of linear regression method used to predict pose parameters form the observed data. Section 4 describes in detail the proposed dynamic Bayesian model and Section 5 formally derives closed-form formulas for estimating the model parameters, i.e. model training. Experiments and comparisons with other methods are described in Section 6. Section 7 draws some conclusions. 1

Related Work

Head-pose tracking has been an actively investigated topic; head-pose estimation and tracking methods were surveyed [START_REF] Murphy-Chutorian | Head pose estimation in computer vision: a survey[END_REF]. Many approaches rely on extracting facial landmarks, on tracking these landmarks over the image sequence and on estimating a rigid transformation between consecutive images, e.g. [START_REF] Gee | Fast visual tracking by temporal consensus[END_REF][START_REF] Uřičář | Detector of facial landmarks learned by the structured output SVM[END_REF], or between consecutive image pairs, e.g. [START_REF] Yang | Model-based head pose tracking with stereovision[END_REF]. Similarly, [START_REF] Maurer | Tracking and learning graphs and pose on image sequences of faces[END_REF] builds a face graph based on the landmarks and tracks this graph over the image sequence. Another landmark-based approach [START_REF] Yao | Using affine correspondence to estimate 3-d facial pose[END_REF] consists of using a 3D model of a generic face that embeds model-centered coordinates of facial landmarks, e.g. nose tip, eyes, lip corners, etc. The model is first fitted to the face detected in the first image and then fitted to the subsequent faces by tracking the landmarks. These methods heavily rely on landmark detection and tracking as well as on the robust estimation of the 2D-landmark-to-3D-landmark rigid transformation, i.e. the pose parameters. Therefore these methods are limited to frontal views of faces, because the landmarks are partially or totally occluded in side views of faces. Moreover, they track the facial landmarks instead of the pose parameters, hence they do not yield smooth pose trajectories. The advantage of the proposed method is that it relies neither on facial landmark detection nor on landmark tracking. The proposed method, once trained based on pairs of HOG descriptors and pose parameters, can deal with side views of faces, unlike landmark-based methods.

Head-pose tracking was also addressed using sampling methods based on particle filters, which allow to sample the temporal predictive distribution e.g. [START_REF] Ba | A probabilistic framework for joint head tracking and pose estimation[END_REF]. A principled way of combining a latent-variable temporal filter with the observed data is an important issue. In [START_REF] Tu | Accurate head pose tracking in low resolution video[END_REF] it is proposed to extract a high-dimensional feature vector from a face and then to apply PCA to reduce its dimensionality. This assumes that the high-dimensional to low-dimensional mapping is linear (which may not be the case) and it does not guarantee that the PCA output contains pose information. Particle filtering can also be combined with a 3D deformable model and with facial landmarks, e.g. [START_REF] Dornaika | Head and facial animation tracking using appearance-adaptive models and particle filters[END_REF][START_REF] Taheri | Joint albedo estimation and pose tracking from video[END_REF]. As already outlined, landmark extraction is not always possible. The advantage of the proposed method over these particle-filter trackers is both theoretical and methodological: the feature-space to parameter-space mapping is combined with a dynamic model, and the estimation of the model parameters yields closed-form EM procedures. Moreover, the proposed SKF approximation, which makes the model computationally tractable, amounts to a variational approximation, thus yielding an extremely efficient runtime method.

Switching state space models have also been used to solve tracking problems. For example, [START_REF] Ghahramani | Switching state-space models[END_REF], [START_REF] Oh | Learning and inference in parametric switching linear dynamic systems[END_REF] and [START_REF] Kooij | A nonparametric hierarchical model to discover behavior dynamics from tracks[END_REF] show that the use of switching linear models helps tracking. In [START_REF] Pavlovic | Learning switching linear models of human motion[END_REF] switching models are applied for tracking people in videos in order to obtain motion-capture data, and three different approaches for inferring the parameters are compared, namely the Viterbi algorithm, variational inference, and the generalized pseudo Bayesian algorithm of order 2 (GPB2). The reported results obtained with these three approaches are quite similar. Viterbi has the lowest complexity, GPB2 yields the smoothest parameter trajectories, while the variational inference achieves a good compromise between low complexity and smooth trajectories.

The proposed method combines high-dimensional to low-dimensional mixture of linear regressions with a switching state-space model. In practice we adopt the GPB2 algorithm which, in combination with the generative regression model, yields closed-form expressions for the estimation of the tracked parameters. Hence, it is more efficient than sampling techniques which are often used in conjunction with generative tracking methods.

Mixture of Linear Inverse Regressions

In this section we summarize the mixture of linear inverse regressions of [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF], which is named Gaussian locally linear mapping (GLLiM). GLLiM interchanges the roles of the input (high dimensional) and of the output (low dimensional), such that a low-to-high regression is being learned. The immediate consequence of this inverse regression strategy is a dramatic reduction in the number of model parameters, thus facilitating the task of training.

Let X and Y be two random variables, and let x and y denote their realizations, where X ∈ R L is the lowdimensional output (pose parameters) and Y ∈ R D (D L) is the high-dimensional input (feature vectors). Once trained, the goal is to predict x given both an input y and

× p(x|Z = k; θ)p(Z = k; θ), (1)
where θ denotes the model parameters. Assuming that e k is a zero-mean Gaussian variable with a diagonal covariance matrix Σ k ∈ R D×D , we obtain that

p(y|x, Z = k; θ) = N (y; A k x + b k , Σ k ). (2)
If we further assume that X follows a mixture of Gaussians via the same assignment Z = k, we can write that

p(x|Z = k; θ) = N (x; c k , Γ k ), p(Z = k; θ) = π k , ( 3
)
where

c k ∈ R L , Γ k ∈ R L×L and K k=1 π k = 1. Note that this representation induces a partition of R L into K regions R k ,
= {c k , Γ k , π k , A k , b k , Σ k } K k=1 .
The model parameters can be estimated via an EM algorithm. The expectation step computes the responsibilities, namely p(Z n = k|x n , y n ; θ (old) ), given the old parameter values θ (old) , while the maximization step computes new parameter values via maximization of the expected complete-data log-likelihood function, namely θ (new) = argmax E[log p(x, y, Z|θ (old) )], which yields a closed-form solution [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF]. Initial responsibilities are obtained by fitting a K-component GMM to the lowdimensional data {x n } N n=1 . Once the inverse parameter vector θ is estimated one obtains a low-dimensional to high-dimensional inverse predictive distribution [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF]. The high-dimensional to lowdimensional forward predictive distribution has a closedform expression:

p(x|y; θ * ) = K k=1 π * k N (y; c * k , Γ * k ) K j=1 π * j N (y; c * j .Γ * j ) N (x; A * k y+b * k , Σ * k ), (4)
which is a Gaussian mixture fully defined by the forward parameters

θ * = {c * k , Γ * k , π * k , A * k , b * k , Σ * k } K k=1
, that can be obtained in analytically from the inverse parameters:

c * k = A k c k + b k , Γ * k = Σ k + A k Γ k A k , π * k = π k , (5) A * k = Σ * k A k Σ -1 k , b * k = Σ * k (Γ -1 k c k -A k Σ -1 k b k ), (6) Σ * k = (Γ -1 k + A k Σ -1 k A k ) -1 . (7)

The Dynamic Bayesian Model

The main difference between the probabilistic regression model outlined in Section 3 and the proposed temporal model is that the conditional distribution p(x|y) is replaced with p (x t |y 1:t ), where t is the time index. The proposed graphical model is shown on Fig. 3, where Z t is the discrete latent variable associated with the Gaussian mixture of linear regression, X t and Y t are the latent head pose and the observed high-dimensional feature vector at t, respectively. Using Bayes rule and marginalization we obtain:

(X t = x t |Y 1:t = y 1:t ) = K j=1 K i=1 Xt-1 1 p(y t ) ×p(x t , x t-1 , Z t = j, Z t-1 = i, y t |y 1:t-1 )dX t-1 . (8)
Under the Markovian assumption and using the conditional independencies associated with the proposed graphical model of Fig. 3, we can write the term inside the integral of (8) as follows:

p(x t , x t-1 , Z t = j, Z t-1 = i, y t |y 1:t-1 ) = p (y t |x t , Z t = j) p (x t |x t-1 , Z t = j) × p (Z t = j|Z t-1 = i) p x t-1 , Z t-1 = i|y 1:t-1 . (9)
The first probability on the right hand side of this equation, p (y t |x t , Z t = j) is the Gaussian distribution introduced in (2). The main difference between the static model and the dynamic model is that (3) is replaced with:

p(x t |x t-1 , Z t = j) = N (x t |C j x t-1 , Q j ), (10)
p(Z t = j|Z t-1 = i) = τ ij . ( 11
)
The parameters of the temporal model will be jointly denoted by φ:

φ = {C j , Q j , τ ij , i, j = 1 . . . K}. (12)
By substituting (2), [START_REF] Drouard | Head pose estimation via probabilistic highdimensional regression[END_REF], and ( 11) into (9), by using basic properties of Gaussian distributions (Gaussian product and Gaussian integral), e.g. [START_REF] Bishop | Pattern recognition and machine learning[END_REF], and after some derivations, one can show that (8) can be written as a K 2 -component GMM:

p(x t |y 1:t ; ψ t|t-1 ) = K i=1 K j=1 π ij t|t-1 N (x t |µ ij t|t-1 , W ij t|t-1 ), (13)
which in turn can be approximated with another Kcomponent GMM, namely:

p(x t |y 1:t ; λ t ) ≈ K j=1 ρ j t N (x t |η j t , V j t ). ( 14
)
The parameters of the these two Gaussian mixtures are denoted with

ψ t|t-1 = {π ij t|t-1 , µ ij t|t-1 , W ij t|t-1 , i, j = 1 . . . K} (15)
and with

λ t = {ρ j t , η j t , V j t , j = 1 . . . K}. (16)
The K-component GMM approximation (13) of the K 2component GMM ( 14) is necessary in order to avoid an exponential grow of the number of components, hence it guarantees the computational tractability of the temporal model. As discussed in [START_REF] Murphy | Switching kalman filters[END_REF], three approaches have been proposed to avoid the number of components to explode. Our approximation is based on the generalized pseudo Bayesian algorithm (GPBa) of order 2 (GPB2) which, according to [START_REF] Pavlovic | Learning switching linear models of human motion[END_REF], yields a smooth output.

One interesting consequence of replacing (3) with [START_REF] Drouard | Head pose estimation via probabilistic highdimensional regression[END_REF] is that the parameter set θ in Section 3 is replaced with a reduced parameter set θ r = {A j , b j , Σ j } K j=1 . Consequently, the formulae in ( 5) and ( 6) are simplified:

A * j = Σ * j A j Σ -1 j , b * j = -A * j b j , Σ * j = (A j Σ -1 j A j ) -1 . ( 17
)
It can be shown that the parameter set (15) can be written as a function of θ r , λ t-1 and φ:

W ij t|t-1 = Σ * j -1 + P ij t-1 -1 , (18)
µ ij t|t-1 = W ij t|t-1 Σ * j -1 A * j y t + b * j + P ij t-1 C j η j t-1 , (19)
π ij t|t-1 ∝ r ij t|t-1 = ρ i t-1 τ ij N (d ij t|t-1 |0, S ij t|t-1 ), (20)
where P ij t-1 , d ij t|t-1 , and S ij t|t-1 are defined by:

P ij t-1 = Q j + C j V i t-1 C j -1 , (21)
d ij t|t-1 = y t -A j (C j η i t-1 ) -b j , (22)
S ij t|t-1 = Σ j + A j (Q j + C j V i t-1 C j )A j . (23)
Eq. ( 21) defines the covariance of the prediction variable dynamics, Eq. ( 22) is the difference between the observation and the predicted observation, given η t-1 , and Eq. ( 23) defines the associated covariance matrix.

The mean ( 19) can be seen as a "weighted" linear combination of the dynamical prediction C j η j t-1 and of the prediction based on observation A * j y t + b * j , where the "weights" are covariance matrices. Thus the confidence related to the covariance matrices defines the weights of the dynamical prediction and the observation prediction in the final estimation. Eq. ( 18) is the associated covariance matrix, which is the inverse of the sum of the precision matrix of the temporal prediction P ij t-1 and precision matrix Σ * j -1

of the observation y t . The GMM proportions in Eq. ( 20) are defined as a product between three terms: the proportions of the i th components at t -1, ρ i t-1 , the switching filter transition probabilities τ ij , and

N (d ij t|t-1 |0, S ij t|t-1 ).
Using the mixture reduction scheme explained in [START_REF] Salmond | Mixture reduction algorithms for point and extended object tracking in clutter[END_REF], the parameters of the K-component GMM λ t can now be evaluated from the parameters of the K 2 -component GMM ψ t|t-1 , with the following formulas:

η j t = K i=1 πij t|t-1 µ ij t|t-1 , (24)
V j t = K i=1 πij t|t-1 W ij t|t-1 + (µ ij t|t-1 -η j t )(µ ij t|t-1 -η j t ) , (25)
ρ j t = K i=1 π ij t|t-1 , with πij t|t-1 = π ij t|t-1 / k k=1 π kj t|t-1 . ( 26
)

Estimating the Model Parameters

This section describes the estimation of the parameters φ in [START_REF] Foytik | A two-layer framework for piecewise linear manifold-based head pose estimation[END_REF] via learning (please consult [START_REF] Murphy | Switching kalman filters[END_REF] for a detailed description). We remind that the estimation of the regression parameters θ is described in detail in [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF] and summarized in Section 3. We use an EM procedure to estimate the parameters C j , Q j K j=1 . During the E-step we compute the complete-data expected log-likelihood. First, we express the complete-data log-likelihood L as follows:

L = log p (x 1:T , Z 1:T , y 1:T ; θ) = T t=1 K j=1 δ(t, j) log N (y t ; A j x t + b j , Σ j ) + T t=2 K j=1 δ(t, j) log N (x t ; C j x t-1 , Q j ) + T t=2 K j=1 K i=1 δ(t, j)δ(t -1, i) log τ ij + log p (x 1 , Z 1 ) , (27)
where δ(t, q) is equal to 1 if Z t = q and 0 otherwise. To complete the E-step we evaluate the above expected loglikelihood. The M-step maximizes it with respect to the parameters that must be estimated. Thus we obtain the following formulas for the estimation of the parameters φ:

C j = T t=2 p (Z t = j|y t ) E x t x T t-1 × T t=2 p (Z t = j|y t ) E x t-1 x T t-1 -1 , (28)
Q j = 1 T t=2 p (Z t = j|y t ) × T t=2 p (Z t |y t ) E x t x T t -C j E x t-1 x T t , (29)
where:

E x t x T t-1 = V t,t-1 + η t η t-1 , E x t-1 x T t-1 = V t-1 + η t-1 η t-1 .
Parameters η t , η t-1 , V t-1 and V t,t-1 are obtained using the smoothing statistics p (x t |y 1:T ) which is expressed as follows:

p(x t |y 1:T ) = K j=1 K i=1 p (x t , Z t = j, Z t+1 = i|y 1:T ) = K j=1 K i=1 p (x t |Z t = j, Z t+1 = i, y 1:T ) ×p (Z t = j, Z t+1 = i|y 1:T ) , (30)
where

p(x t |Z t = j, Z t+1 = i, y 1:T ) =p (x t |Z t = j, y 1:t ) p y t+1:T |Z t+1 = i, x t . ( 31
)
The first term of (31) p (x t |Z t = j, y 1:t ) is the forward recursion and is defined in Section 4 as N x t |η j t , V j t . The second term p y t+1:T |Z t+1 = i, x t defines the backward recursion:

p(y t+1:T |Z t+1 = i, x t ) = Xt+1 p (x t+1 |x t , Z t+1 = i) ×p y t+1 |x t+1 , Z t+1 = i p y t+2:T |x t+1 dX t+1 , (32)
where p y t+2:

T |x t+1 = N x t+1 |η b t+1 , V b t+1 .
The second term of Eq. ( 30) can be decomposed as follows, using [START_REF] Kim | Dynamic linear models with markov-switching[END_REF]:

p (Z t = j, Z t+1 = i|y 1:T ) p (Z t+1 = i|y 1:T ) p (Z t+1 = i|Z t = j) p (Z t = j|y 1:t ) p (Z t+1 = i|y 1:t ) , (33)
where:

p (Z t = j|y 1:t ) = ρ j t (34)
p (Z t+1 = i|y 1:t ) = K j=1 p (Z t+1 = i|Z t = j) p (Z t = j|y 1:t ) (35)
p (Z t+1 = i|y 1:T ) = K j=1 p (Z t+1 = i, Z t+2 = j|y 1:T ) (36)
As outlined in Section 4, the number of component increases [START_REF] Oh | Learning and inference in parametric switching linear dynamic systems[END_REF], hence we merge the Gaussians twice: first over Z t+1 to obtain a mixture of K Gaussian components with mean η jb t , covariance V jb t and proportions p (Z t = j|y 1:T ) = ρ jb t , second over Z t and thus we obtain a single Gaussian component, with mean η t and covariance V t .

To estimate the transition matrix {τ ij } K i,j=1 we employ the Lagrange multiplier method to maximize the loglikelihood with respect to τ ij , hence we obtain the following expression for the transition probabilities:

τ ij = T t=2 p (Z t = j|y 1:T ) p (Z t-1 = i|y 1:T ) T t=2 p (Z t-1 = i|y 1:T ) (37)

Experiments

To gauge the performance of the proposed method we used two datasets: the Biwi-Kinect head-pose dataset [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF] and the EYEDIAP dataset [START_REF] Mora | Eyediap: a database for the development and evaluation of gaze estimation algorithms from rgb and rgb-d cameras[END_REF]. Biwi-Kinect comprises 24 videos of 20 different people (16 men and 4 women) recorded with a Kinect camera. During the recordings people were asked to move their heads freely in front of the camera. 3D head pose (pitch, yaw, and roll angles) annotations are automatically and accurately provided for each video frame using the face-shift software. The angle values range from -60 • to 60 • for pitch, -75 • to 75 • for yaw and -20 • to 20 • for roll. The dataset provides RGB and depth images as well as the calibration matrices. The 3D nose positions are provided as well. EYEDIAP is a dataset for gaze and head-pose estimation. It provides 94 videos of 16 people recorded using different configurations, such as static and turning heads. The dataset provides RGB videos (in both HD and GVA quality) and depth videos with the associated calibrations matrices. For each video, annotations of both head-pose and gaze are provided for each frame. The angle values range from -40 • to 40 • for pitch and -50 • to 50 • for yaw. In our experiments we only used the RGB images where people are looking at a moving object (named A FT M in the dataset).

The proposed model, referred to as HPE SKF (headpose estimation based on SKF) is compared to the following methods: (i) a landmark-based approach that uses the facial landmark localization method of [START_REF] Uřičář | Detector of facial landmarks learned by the structured output SVM[END_REF] (Flandmarks) combined with 2D-to-3D landmark-based pose estimation method, namely the PnP (perspective n-point) algorithm available with OpenCV, (ii) a depth head model based on method [START_REF] Mora | Gaze estimation from multimodal kinect data[END_REF] that learns a 3D head model using 16 manually annotated facial landmarks on multiple frames to learn the model, and ICP (iterative closest point algorithm) to estimate the transformation (rotation and translation) of the head pose between two consecutive frames in order to track the pose over time, (iii) the regression-based method of [START_REF] Drouard | Head pose estimation via probabilistic highdimensional regression[END_REF] which is referred to as HPE-GLLiM, and (iv) the regression method [START_REF] Drouard | Head pose estimation via probabilistic highdimensional regression[END_REF] combined with a standard Kalman filter [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Bishop | Pattern recognition and machine learning[END_REF]. Both (i) and (ii) perform tracking. For evaluation, on the Biwi-Kinect dataset we used the leave-one-out protocol: all the data related to one person are put aside and the remaining data of the other persons are used to train the model.

Using the EYEDIAP dataset, we compared our method against the baseline method of [START_REF] Mora | Gaze estimation from multimodal kinect data[END_REF] which requires to learn a person based on his/her 3D head model. In this case we did not apply the leave-one-out protocol but instead we learned a person based head-pose regression model (θ) on a subsample of the frames associated with each person. Then, for each person we estimate the tracking parameters (φ) over the whole video. Thus, we obtain a person-based tracking model for head pose. As a measure of performance, we use the mean absolute deviation between the ground-truth and the estimated pose.

Face regions are extracted from images with a face detector [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF]. This detector is efficient and robust with both frontal-and side-views of faces. Using the detection we run a face tracker using particle filtering to extract a face at each frames of the videos. Nevertheless, the obtained face regions are noisy, i.e. the bounding boxes are not always nicely aligned onto the faces. This yields extremely realistic input data for the tested methods, unlike other head-pose benchmarks that use manually extracted bounding boxes, e.g. using the nose tip as the bounding box center. From each face region thus detected, we extract HOG features with different cell resolutions (in a pyramid-like fashion), namely 32 × 32, 16 × 16 and 8 × 8 pixels, with block size of 2 × 2 cells and 8 bins to quantize the gradient orientation. This results in feature vectors of size D = 1888.

The regression parameters θ (Section 3) and the SKF parameters φ Eq. ( 12) are learned separately. First, the regression parameters θ are estimated using the EM algorithm described in [START_REF] Deleforge | High-dimensional regression with Gaussian mixtures and partially-latent response variables[END_REF]. Second, the filtering parameters φ are estimated using the method described in Section 5. The covariances are initialized with identity matrices and the transition matrix {τ ij } K i,j=1 is initialized with the Battacharrya distance [START_REF] Bhattachayya | On a measure of divergence between two statistical population defined by their population distributions[END_REF] between two subspaces obtained using the parameters θ of the low dimensional space defined by Z. The results obtained with the Biwi-Kinect and EYE-DIAP datasets are summarized in Table 1 and Table 2, respectively, namely the average and standard deviations of the absolute error between the head-pose estimated values and the ground-truth values. The proposed tracking method improves head-pose parameter tracking, compared to all the other methods. For example, the average error for the yaw angles in Table 1 is of 8.6 • while all other methods yield an error larger than 12 • . We observe the same behavior for the pitch angle. Moreover, our method also reduces the standard deviation. Compared to the landmark-based method of [START_REF] Uřičář | Detector of facial landmarks learned by the structured output SVM[END_REF], our method is able to provide an estimation for each test input, whereas the method based on landmarks is unable to provide an output when some of the landmarks are not visible due to extreme head orientations. In this case [START_REF] Uřičář | Detector of facial landmarks learned by the structured output SVM[END_REF] yields very large errors, e.g. first row of Table 1. We also note that the proposed HPE SKF method performs much better than a standard Kalman filter. On Table 2 the global average error and standard deviations are lower than Table 1 due to the person-based model learning. The behavior shown in Table 2 follows the same as the one of Table 1, i.e. our tracking method decreases the head pose estimation average error and standard deviation with respect to the HPE GLLiM and the classic Kalman filtering. [START_REF] Mora | Gaze estimation from multimodal kinect data[END_REF] (that released the EYEDIAP dataset) achieves better results for pitch, while our method obtains better results for yaw.

Conclusions

Head-pose trackers have the advantage of combining information from both past and present, and hence they avoid oscillations between consecutive estimations solely based on independent observations, as it is the case with many head-pose estimation methods, e.g. Fig. 2 and Fig. 4. Overall, the output of our tracking method is both more accurate and smoother than the output of several head-pose methods (with and without tracking). Moreover, noisy observations, e.g. due to badly aligned bounding boxes or to partial occlusions, do not impact too much the proposed tracker because the temporal model, once properly trained, does not allow oscillations between consecutive estimations. Furthermore, our approach does not only reduce both the estimation error and the standard deviation, it smoothes the estimation of head poses over the whole video. This is very useful for other temporal tasks, such as the estimation of eye gaze or of the visual focus of attention [START_REF] Massé | Simultaneous estimation of gaze direction and visual focus of attention for multi-personto-robot interaction[END_REF].

The method presented in this paper is very general and it is limited neither to head pose nor to HOG features. It can be applied to estimate and track the pose parameters of objects of all kinds, provided that their image appearance varies as a function of their 3D orientation in a consistent way. Therefore, one can combine our tracker with other types of features, such as features obtained by training a neural network and by substituting the last layer of the latter with the generative regression model used above.

Figure 1 .

 1 Figure 1. The method starts by learning a mixture of linear regression that allows the prediction of a head-pose from a HOG vector estimated from the bounding box of a face. Hence, Eq. (4) (Section 3) is applied at t -1 (top) and at t (bottom) and head poses are thus predicted, they are denoted A on the figure. Notice that, because of various perturbations in the data and of inherent flaws in face detection, the two predictions use two different affine transformations and hence they are associated with two different Gaussian components in the mixture, e.g. magenta and green on the figure. The proposed dynamic model combines the temporal prediction of the filter from t -1 to t, denoted B on the figure, with the pose predicted at t, to yield a filtered pose estimate, denoted C on the figure. The mixture of linear regression is plugged in the SKF model in a principled way.

Figure 2 .

 2 Figure 2. Yaw angles predicted with the mixture of linear regression method [10] (top sequence and red plot) and with the proposed method (bottom sequence and green plot). The ground-truth yaw trajectory is shown in blue.

  where R k is the region where the transformation (A k , b k ) is most likely invoked, e.g. Fig. 1. This model is described by the inverse parameter set θ

Figure 3 .

 3 Figure 3. Proposed graphical model

Figure 4 .

 4 Figure 4. Comparison between the estimated pitch angle (top) and yaw angle (bottom) with three different methods: HPE GLLiM (red), Kalman Filter (blue), and the proposed HPE SKF (green) for the Biwi-Kinect dataset.

Table 1 .

 1 Average (Avg.) and standard deviation (Std.) of the absolute error (in degrees) for the pitch, yaw and roll angles (when applicable) on the Biwi-Kinect dataset. Head bounding boxes are extracted using a face detection algorithm.

		Pitch	Yaw	Roll
	Methods Avg.	Std. Avg.	Std. Avg. Std.
	[39]	13.12 10.79 21.1 14.16 -	-
	[10]	14.35 13.73 12.52 13.52 10.89 9.82
	[1, 6]	27.43 17.61 15.13 12.09 14.62 9.94
	HPE SKF 10.03 8.73 8.6	7.21 8.48 8.01

Table 2 .

 2 Average (Avg.) and standard deviation (Std.) of the absolute error (in degrees) for the pitch and yaw angles on the EYE-DIAP dataset. Head bounding boxes are extracted using a face detection algorithm combined with a face tracker.

		Pitch	Yaw
	Methods	Avg.	Std.	Avg.	Std.
	Funes Mora et al. [26] 4.17	5.59	6.89 14.42
	HPE GLLiM [10]	7.94	9.23 10.62 11.95
	Kalman filter [1, 6]	23.17 18.67 25.55 21.22
	HPE SKF	5.34	8.30	6.68	9.76

Supplementary material for this paper can be found at https://team.inria.fr/perception/research/ head-pose-tracking/

* Funding from the European Union FP7 ERC Advanced Grant VHIA (#340113) is greatly acknowledged.