Sławomir Bąk 
email: slawomir.bak@inria.fr
  
Krzysztof Kurowski 
email: krzysztof.kurowski@man.poznan.pl
  
Krystyna Napierała 
  
Human Re-identification System On Highly Parallel GPU and CPU Architectures

Keywords: Re-identification, Covariance Matrix, Generalized Eigenvalues, High Performance Computing, GPU

The paper presents a new approach to the human reindetification problem using covariance features. In many cases, a distance operator between signatures, based on generalized eigenvalues, has to be computed efficiently, especially once the real-time response time is expected from the system. This is a challenging problem as many procedures are in fact computationally intensive tasks and must be repeated constantly. To deal with this problem we have successfully designed and then tested a new video surveillance system. To obtain the required high efficiency we took the advantage of highly parallel computing architectures such as FPGA, GPU and CPU units to perform calculations. However, we had to propose a new GPU-based implementation of the distance operator for querying the example database. Thus, in this paper we present experimental evaluation of the proposed solution in the light of the database response time depending on its size.

Introduction

Human re-identification is one of the most challenging and important problems in computer vision and pattern recognition. The re-identification problem can be defined as a determination whether a given person of interest has already been observed over a network of cameras. This issue (also called the person re-identification problem) can be considered on different levels depending on information cues currently available in the system. Biometrics such as face, iris or gait can be used to recognize identities. Nevertheless, in most video surveillance scenarios such detailed information is not available due to a low video resolution or a difficult segmentation (crowded environments, e.g. airports, metro stations). Therefore a robust modeling of a global appearance of an individual is necessary to re-identify a given person of interest. In these identification techniques (named appearance-based approaches) clothing is the most reliable information about an identity of an individual (there is an assumption that individuals wear the same clothes between different sightings). A model of appearance has to handle differences in illumination, pose and camera parameters to allow matching appearances of the same individual observed in different cameras. High accuracy of re-idenfication approaches can only be achieved using an appearance representation based on descriptors which are invariant across different camera views. Recently, a covariance descriptor [START_REF] Tuzel | Region covariance: A fast descriptor for detection and classification[END_REF] has proved its effectiveness in recognition [START_REF] Bak | Person re-identification using spatial covariance regions of human body parts[END_REF] and classification approaches [START_REF] Tosato | Multi-class classification on riemannian manifolds for video surveillance[END_REF]. It has been shown that the performance of the covariance features is superior to other methods as rotation and illumination changes are absorbed by the covariance matrix [START_REF] Bak | Person re-identification using spatial covariance regions of human body parts[END_REF]. Moreover, integral images [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] used for fast covariance computation make this descriptor very efficient concerning extraction of the covariance.

However, as covariance matrices does not lay on Euclidean space, there is necessary to apply complex differential geometry to compute a distance between two covariances. The distance operator is computationally heavy as it requires solving the generalized eigenvalues problem. As a consequence, matching of the covariance descriptors is slower than matching of other computer vision descriptors, which are usually represented by vectors laying on Euclidean space. This often makes covariance-based approaches difficult to apply in real-time systems in spite of their effectiveness. Moreover, in the person re-identication problem there is usually a large number of candidate matches which makes the issue much more challenging (concerning the matching accuracy as well as the response-time of a database of human appearances). Hence, we propose a new hybrid architecture based on Graphics Processing Units (GPUs) and Field Programmable Gate Arrays (FPGAs) accelerators to take advantage of high performance computing and make covariance-based approaches applicable to large-scale databases. This paper makes the following contributions:

-We describe a new GPU-and FPGA-based architecture for the person reidentification problem. This architecture can be easily adjusted to more general video surveillance problems (such as object recognition or object classification)(Section 3). -We propose an implementation for finding generalized eigenvalues and eigenvectors for distance operator, using NVIDIA GPU architecture (Section 5).

We evaluate our approach in Section 6 before concluding the paper.

Motivations and Related Work

There is an increasing demand for effective surveillance systems, i.e. systems that can perform low-cost, low-power and high-speed operations. Recently, recognition problems became one of the most important tasks in video surveillance.

As recognition is an extremely difficult task, the existing approaches are computationally heavy. Thus, a new high performance architectures are necessary to apply these approaches to real-time systems. In [START_REF] Sriram | Design-space exploration of biologically-inspired visual object recognition algorithms using cpus, gpus, and fpgas[END_REF] biologically-inspired algorithms are adjusted to GPU to perform large-scale object recognition. Similarly, in [START_REF] Poli | Processing neocognitron of face recognition on high performance environment based on gpu with cuda architecture[END_REF] GPU-based neural network is presented to recognize human faces. We offer a surveillance system based on FPGA and GPU architectures giving much better computing facilities in comparison with traditional CPU-based systems. The most demanding part of the system -the generalized eigenvalues calculations -is performed using the GPU architecture. There are already a few approaches of finding eigenvalues of symmetric matrices using GPUs [START_REF] Lessig | Eigenvalue computation with cuda[END_REF], but these implementations are focused on computation of the eigenvalues of large matrices (the implementations are optimized for matrices larger than 1024). In contrary to these approaches, we concern a domain-specific problem where it is necessary to solve the generalized eigenvalues problem for a large number of small matrices (the covariance descriptor is mostly represented by square matrix of size between 8 and 16). To the best of our knowledge, we are the first to propose an implementation for finding the generalized eigenvalues and eigenvectors of a large amount of small matrices using NVIDIA GPU cards.

System Architecture

Our GPU-and FPGA-based surveillance system for the person re-identification problem (Fig. 1) is designed to assign the tasks to the most suitable architectures.

The system consists of a network of cameras with embedded FPGA units, which preprocess a video stream (image denoising, object detection and classification). FPGA units are well suited for video processing tasks usually exposing a high level of parallelism. Therefore, only the necessary information (blobs of interest) is sent in the network, without the need of transfering the whole video stream to the central unit. The partially transformed data is collected on a central unit with a CPU and GPU processor. In our approach a human appearance is represented by a set of covariance matrices extracted on different resolutions from detected body parts [START_REF] Bak | Person re-identification using spatial covariance regions of human body parts[END_REF]. In total, the human signature is represented by 26 covariance matrices of the size 11. Covariance signature can be computed on a CPU as there exist an efficient way to extract covariances using integral images.

Signature is then stored on a GPU unit which serves as a database of signatures.

The advantage of such a solution is that, first, CPU unit is offloaded from storing this information, and second, when the distance is calculated, the data is already stored on a proper unit. We use a Tesla S1070 with four GPU units, on which there is 4GB of available global memory for each unit. Taking into account the free space needed for calculations for a query to a database, about 200,000 signatures can be stored in the database on one unit, which is sufficient for the purposes of our system. The user (administrator of a surveillance system) can select an object of interest and query the database with a new signature. The new signature is compared using the distance operator to all signatures in the database. Distance operator is calculated in parallel directly on the GPU unit.

The result is a list of the most similar signatures.

The most important part in the system is the database stored on GPU and the calculation of distance on this architecture. The preprocessing on FPGA is not that crucial for making the system real-time, therefore we show only how to calculate the distance on GPU (Sec. 5), and we evaluate experimentally the time of the database response depending on its size (Sec. 6).

Covariance Descriptor

In [START_REF] Tuzel | Region covariance: A fast descriptor for detection and classification[END_REF] the covariance of d-features has been proposed to characterize a region of interest. Now, we introduce the geodesic distance definition proposed by [START_REF] Förstner | A metric for covariance matrices[END_REF] as its computation is the main topic of our work. The distance between two covariance is defined as

ρ(C i , C j ) = d k=1 ln 2 λ k (C i , C j ) (1) 
where

λ k (C i , C j ) k=1...d are the generalized eigenvalues of C i and C j , determined by λ k C i x k -C j x k = 0, k = 1 . . . d and x k = 0 are the generalized eigenvectors.
The complexity of finding generalized eigenvalues increases rapidly with the size of d. The efficient methods work for dimensions below 5. Often, a dimension of covariance matrix has to be decreased to conform requirements of a real-time systems at the expense of the accuracy. Hence, as the distance computation is the main bottleneck in covariance-based approaches we decided to take advantage of GPU to speed up the computation of the generalized eigenvalues.

GPU Implementation

GPU Architecture

GPU is a high performance computing unit in which the emphasis is put on the computing units instead of data caching and control flow units in contrast with CPU. The memory is not cached, so to obtain the maximum performance it is crucial that the programmer ensures the coalesced memory accesses. In our setup we use the NVIDIA Tesla 1070S with 4 graphic cards, each consiting of 30 SMs (Streaming Multiprocessors) with 8 scalar processors. The total number of available processors is therefore equal to 240x4 GPU. The threads work in a SIMD model (Single Instruction Multiple Data) and are grouped into blocks. All threads of a block reside on the same processor core. Threads within one block are split into basic scheduling units called warps, consisting of 32 threads. A warp executes one common instruction at a time, so full efficiency is realized when all 32 threads of a warp agree on their execution path. On GPU there are different types of memory. The two most important are shared memory and global memory. Shared memory is allocated per block and shared among all threads of a block. It is organized in banks, and if the data is accessed so that each thread of a block accesses a different bank, shared memory is as fast as registers. Global memory is a slow off-chip memory, which can be accessed by both CPU and GPU. It is also used to synchronize data between threads in different blocks. In order to have a low latency, it should be accessed in a coalesced fashion.

Generalized Eigenvalues Problem

The most computationally heavy component of calculating the distance between the signatures is the calculation of generalized eigenvalues of positive definite symmetric matrices A and B (two corresponding covariance matrices from the compared signatures). It is defined by the equation Ax = λBx, where λ denotes a vector of eigenvalues and x is the eigenvector. This equation can be decomposed to the equation

(L -1 AL -T )L T x = λ(L T x) (2) 
where L is the upper triangular matrix calculated as B = LL T . We can notice that the decomposed Eq. ( 2) already corresponds to original eigenvalues problem. We solve the generalized eigenvalues problem in 4 steps: Let us note that there are many methods to calculate the eigenvalues of a symmetric matrix. On the CPU we use a Jacobi algorithm, which does not need to perform the tridiagonalization first, however this algorithm is not easy to parallelize. For the GPU implementation we first tridiagonalize the matrix and then we use the bisection algorithm, which can be parallelized much more efficiently.

Porting the Algorithm to the GPU Architecture

Two sources of parallelism can be used in the algorithm. The first one is obvious -a comparison of two signatures involves comparing n pairs of covariance matrices, which can be naturally processed in parallel. The second one involves the parallelism extracted from each step of the equation performed on a given pair of covariance matrices. We will now briefly describe how we use the parallelism of each of the procedures to efficiently port it to the GPU architecture.

In Cholesky Decomposition of n matrices, the formula to calculate each element of the L upper triangular matrix is given as

L j,j = B j,j - j-1 k=1 L 2 j,k , L i,j = 1 L j,j (A i,j - j-1 k=1 L i,k L j,k ) f or i > j (3)
The matrix B is loaded from global memory to shared memory in a coalesced way. After the calculations, the matrix L is loaded back to global memory, also assuring coalescence. For the calculations we use the Cholesky-Crout algorithm, which starts from the upper left corner of the matrix L and proceeds to calculate the matrix column by column, as processing data in column order ensures the memory coalescence. To calculate an element of a column, only the elements from the columns to the left are needed. So, all the elements in one column can be processed in parallel. For our matrices of size 11, we can therefore use 11 threads to calculate each column in parallel, using 11 iterations to calculate the whole matrix. The natural decomposition of a problem would be to assign one matrix to a block. However, 11 threads is much less than the size of a warp, which is the smallest scheduling unit. Assigning two matrices to a block enables to process two matrices without the increase of processing time, as 22 threads still form only one warp scheduled in one operation. More matrices however will lead to the creation of more warps, which in our specific situation (calculation of each matrix is independent) makes it similar to the creation of more blocks -it does not matter if in the next step another warp or another block is scheduled. So in our implementation there are n/2 blocks, each block consists of 22 threads and calculates two matrices. Our preliminary experiments showed that indeed, inreasing size of a block from 11 to 22 decreased the time of calculating n matrices twice, while futher increase of the size of a block (up to 110 with 10 matrices assigned) did not bring further improvement. Forward Substitution solves the equation Lx = A. Two forward substitutions solve the equation C = (L -1 AL -T )L T X. We process them together to avoid copying the intermediate data to the global memory. An element x m,k of a matrix

x is calculated as x m,k = Am- P m-1 i=1 Lm,ix i,k Lm,m
. All the columns of x are calculated independently. To calculate an element in the column, one could calculate the sum in the equation using map-reduce model, however for small matrices of size 11 it does not save much calculations so we calculate them sequentially. As a result, for one matrix we need 11 threads, and following our observations from Cholesky Distance procedure, we assign 22 threads to one block. Tridiagonalization is the part of the generalized eigenvalues algorithm which has the lowest level of paralellism. To tridiagonalize the symmetric matrix we use the Householder transformation [START_REF] Householder | Unitary triangularization of a nonsymmetric matrix[END_REF]. In this algorithm, n -2 iterations need to be performed sequentially; in each iteration the appropriate elements in the k-th row and column are zeroed. In one iteration some of the computations, such as matrix multiplication and vector-vector multiplications can be parallelised. To do this, we use n × n threads for each matrix. We also tested the version in which there are only 11 threads, each calculating one column of the multiplied matrices, but it was less efficient. As 121 threads is more than a size of a warp, we can assign one matrix per block without loosing the efficiency. Bisection Algorithm is used to calculate the eigenvalues of a symmetric tridiagonal matrix C. A detailed description of the bisection algorithm can be found e.g. in [START_REF] Demmel | Applied numerical linear algebra[END_REF]. This algorithm finds all the eigenvalues of a matrix with a given approximation. The core function of this algorithm is the Count() procedure returning the number of eigenvalues present in a given interval. The algorithm starts with the initial interval constructed using Gerschgorin's theorem. Then, it is divided in two and Count() procedure returns a number of eigenvalues in each subset. If it uquals zero, the node is abandoned, otherwise it is further subdivided into two subsets, until the size of a subset is not bigger than the assumed approximation. For our purposes it is enough to use the approximation equal to 10e-6. The main source of parallelisation comes from the fact that the Count() function in each node can be calculated independently. As only 11 eigenvalues can be found, on each level of the binary tree there will be only up to 11 nodes containing at least one eigenvalue. We therefore use 11 threads to calculate one matrix. There are also some other sources of parallelism in the calculation of Count() function itself, but we will not present it here. Again, we assign 22 threads to a block, calculating 2 matrices. Gerschgorin's procedure cannot be paralellized, so we process it sequentially, but a parallelism is obtained by calculating all the procedures for different matrices in parallel.

Experimental results

In our experimental setup, we calculated the performance for matrices number ranging from 10 to 3000. As we described in Section 3, we assume that the database of signatures is stored directly on the GPU. In our time estimation, we do not take into account the time of the data transfer, because the reference signatures already reside in the device memory, and the time of transferring the query signature is negligible. Below we present the speedup obtained with comparison to the optimal version on the CPU -that is, a version with Jacobi algorithm (implementation from LTI library), as on CPU this version is faster than calculation of tridiagonalization and then bisection algorithm. The results are presented in Fig. 2(a). One can easily note that the speedup grows with the number of matrices, and reaches its maximum (66) from about 1500 matrices (corresponding to about 50 signatures). Distributing the database of signatures equally between the GPU cards, and performing the calculations for a query signature on all the GPU cards in parallel would result in further speedup improvements. In these tests we used oneTesla node consisting of 4 GPUs. Table 1 presents the time of the component procedures for 5 numbers of matrices. The Cholesky and Forward Substitution times are the lowest, while the biggest impact on the total time comes from the tridiagonalization procedure, as it has the lowest degree of parallelism. Nevertheless, for example for the size of 3000, the total time on the GPU is equal to 3.76ms, while the time on the CPU is equal to 240ms. In Fig. 2(a) one can notice the "stairs" (decrease of speedup) which occur in regular intervals every 480 matrices. This is a result of a similar effect observed in component functions. Let us analyse it on the time performance of the Bisection Algorithm (Fig. 2(b)). One can clearly see that the time rises suddenly every 480 matrices. This is correlated to the number of multiprocessors on the GPU cards. The architecture is the most efficiently used when all the processors (on Tesla S1070 -240 processors) have some blocks assigned. As in our model each block calculates two matrices, the architecture is used most efficiently when k*480 matrices are processed. Otherwise some processors remain idle. The worst case is when k*480+1 matrices are processed -then, the time is almost equal to processing (k+1)*480 matrices, which results in sudden decrease of speedup in these points.

Conclusions

In this paper we demonstrate that it is possible to improve significantly the performance of example video surveillance procedures, in particular the distance operator for querying the database. We observe that in our approach the speedup grows with the number of matrices. Moreover, we can easily distribute demanding calculations of signatures on many GPU units and obtain very good scalability of the system. Although GPU-based approach of four routines of generalized eigenvalues problem have been already proposed, they are optimized for solving only one large matrix. In our approach we use many very small matrices, which can be efficiently stored in a shared memory onto many GPUs. This requires different memory alignments, memory access optimization and simplification of some sub-procedures, but as tested experimentally could yield much better performance.

Fig. 1 .

 1 Fig. 1. The GPU-based architecture for the person re-identification.

1 .

 1 Calculate the Cholesky Decomposition to solve the equation B = LL T 2. Calculate twice the Forward Substitution to solve C = (L -1 AL -T )L T x 3. Tridiagonalize the symmetric matrix C to prepare for solving the eigenvalues problem 4. Use the Bisection Algorithm to find the eigenvalues of a symmetric tridiagonal matrix C

Fig. 2 .

 2 Fig. 2. Speedup for finding generalized eigenvalues and time of Bisection algorithm.

Table 1 .

 1 Time[ms] of component procedures

	N	Cholesky Forward.S Tridiag. Bisect. Total.GPU Total.CPU
	200	0.037	0.035	0.140 0.147	0.359	16
	400	0.049	0.071	0.272 0.187	0.579	32
	600	0.082	0.078	0.389 0.325	0.874	48
	800	0.093	0.112	0.522 0.352	1.08	64
	1000 0.126	0.120	0.657 0.505	1.41	80