Rui Zeng

Jiasong Wu

Zhuhong Shao

Yang Chen

Beijing Chen

Senhadji Lotfi

Huazhong Shu

Color Image Classification via Quaternion Principal Component Analysis Network

Keywords: Deep learning, convolutional neural network, quaternion, QPCANet, PCANet, color image classification

The principal component analysis network (PCANet), which is one of the recently proposed deep learning architectures, achieves the state-of-the-art classification accuracy in various datasets and reveals a simple baseline for deep learning networks. However, the performance of PCANet may be degraded when dealing with color images due to the fact that the architecture of PCANet cannot properly utilize the spatial relationship between each color channel in three dimensional color image. In this paper, a quaternion principal component analysis network (QPCANet), which extends PCANet by using quaternion theory, is proposed for color image classification. Comparing with PCANet, the proposed QPCANet takes into account the spatial distribution information of RGB channels in color images and ensures larger amount of intra-class invariance by using quaternion domain representation for color images.

Experiments conducted on different color image datasets such as UC Merced Land Use, Georgia Tech face, CURet and Caltech-101 have revealed that the proposed QPCANet generally achieves higher classification accuracy than PCANet in color image classification task. The experimental results also verify that QPCANet has much better rotation invariance than PCANet when color image dataset contains lots of rotation information and demonstrate even a simple one layer 2 QPCANet may obtain satisfied accuracy when compared with two layer PCANet.

Introduction

In recent years, image classification becomes a more and more active research topic in the field of pattern recognition and computer vision. It is a very challenging task since image content recognition is subject to various changes of image in terms of illumination, rotation, scaling or more complex deformation. To effectively counter the intra-class variability of image content, numerous approaches have been proposed in the past decades. These methods can be generally divided into two groups: manually extracting intrinsic feature [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Jain | Object detection using Gabor filters[END_REF][START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] and learning unsupervised/supervised feature from data of interest [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF].

For the first group, the most commonly used methods are scale-invariant feature transforms (SIFT) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], Gabor features [START_REF] Jain | Object detection using Gabor filters[END_REF], and local binary pattern (LBP) [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. Both of them achieve great success in different kinds of image classification tasks, such as object recognition, texture classification, and face recognition. However, the limitation of this sort of methods is obviously that we need to manually select the proper method when dealing with new image classification task, where usually requires new domain knowledge.

To remedy the limitation of above hand-crafted methods, the second group, which is well known as deep learning approaches [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF][START_REF] Lecun | Deep learning[END_REF][START_REF] Bengio | Learning deep architectures for AI, Foundations and trends® in Machine Learning[END_REF][START_REF] Bengio | Representation learning: a review and new perspectives[END_REF][START_REF] Guo | Deep learning for visual understanding: A review[END_REF], has become a very active research field in recent years.

The objective of deep learning [START_REF] Bengio | Learning deep architectures for AI, Foundations and trends® in Machine Learning[END_REF] is to extract intrinsic features from data by a multi-level architecture, with the hope that higher-level features represent more concise semantics of the data.

Many approaches have been used to build deep learning architectures, which can be roughly classified into four classes by Bengio et al. [START_REF] Bengio | Representation learning: a review and new perspectives[END_REF]: [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] deep learning architectures based on probabilistic models, for example, deep belief network (DBN) [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF] and deep boltzmann machines (DBM) [START_REF] Salakhutdinov | Deep boltzmann machines[END_REF], etc. (2) deep learning architectures based on reconstruction algorithms, for example, deep autoencoder [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] and sparse coding [START_REF] Olshausen | Emergence of Simple-Cell Receptive Field Properties by Learning a Sparse Code for Natural Images[END_REF], etc. (3) deep learning architectures based on manifold learning algorithms, for example, local coordinate coding (LCC) algorithm [START_REF] Yu | Nonlinear Learning Using Local Coordinate Coding[END_REF][START_REF] Lin | Deep Coding Network[END_REF], etc.

(4) deep learning architectures based on multilayer neural networks, in which convolutional neural network (CNN) [START_REF] Lecun | Backpropagation applied to handwritten zip code recognition[END_REF][START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Russakovsky | Imagenet large scale visual recognition challenge[END_REF][START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF][START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF][START_REF] Szegedy | Going deeper with convolutions[END_REF][START_REF] He | Deep Residual Learning for Image Recognition[END_REF][START_REF] Xu | A deep convolutional neural network for segmenting and classifying epithelial and stromal regions in histopathological images[END_REF][START_REF] Qin | An empirical convolutional neural network approach for semantic relation classification[END_REF][START_REF] Guo | Convolutional feature learning and Hybrid CNN-HMM for scene number recognition[END_REF], first proposed by LeCun et al. [START_REF] Lecun | Backpropagation applied to handwritten zip code recognition[END_REF][START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], gradually becomes the mainstream deep learning structure for image classification after that AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] won the champion of ImageNet [START_REF] Russakovsky | Imagenet large scale visual recognition challenge[END_REF] large scale visual recognition challenge of 2012. Then, Overfeat [START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF],

VGGNet [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF], GoogleNet [START_REF] Szegedy | Going deeper with convolutions[END_REF], and ResNet [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] were further proposed and obtained much lower error rates than AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] in the ImageNet challenge in the year from 2013 to 2015. Generally speaking, multilayer CNN [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF][START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF][START_REF] Szegedy | Going deeper with convolutions[END_REF][START_REF] He | Deep Residual Learning for Image Recognition[END_REF] with many parameters is very suitable for the classification of large scale image datasets, for example, ImageNet database [START_REF] Russakovsky | Imagenet large scale visual recognition challenge[END_REF]. However, when the size of image datasets is not very large, two or three-layer convolutional networks with few parameters may be enough for a good classification performance. For example, Mallat and Bruna [START_REF] Mallat | Group invariant scattering[END_REF][START_REF] Bruna | Invariant scattering convolution networks[END_REF] proposed a mathematical justified wavelet scattering network (ScatNet) whose convolutional layer, nonlinear layer, pooling layer are constructed by prefixed complex wavelets, modulus operator, and average operator, respectively. More recently, Chan et al. [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF] proposed a new deep learning architecture named principal components analysis network (PCANet) whose convolutional layer, nonlinear layer, pooling layer are constructed by principal components analysis (PCA) filters, binary hashing, and block-wise histograms, respectively. Although PCANet uses the most basic and simple operations, it is quite on par with and often better than the state-of-the-art techniques of feature selection for most image classification tasks, including face images, hand-written digits, texture images, and object images. Meanwhile, Chan et al. [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF] proposed the linear discriminant analysis network (LDANet) as a variant of PCANet. The construction method of PCANet attracts the attentions of many researchers. Gan et al. [START_REF] Gan | A deep graph embedding network model for face recognition[END_REF] proposed a deep graph embedding network (GENet) for face recognition and Feng et al. [START_REF] Feng | DLANet: A manifold-learning-based discriminative feature learning network for scene classification[END_REF] presented a discriminative locality alignment network (DLANet) for scene classification. Jia et al. [START_REF] Jia | 2DPCANet: Dayside aurora classification based on deep learning[END_REF] proposed two-dimensional PCANet (2DPCANet) for dayside aurora classification. Qin et al. [START_REF] Qin | DeepFish: Accurate underwater live fish recognition with a deep architecture[END_REF] combined PCANet and spatial pyramid pooling for underwater live fish recognition. Zhao et al. [START_REF] Zhao | Multi-level modified finite radon transform network for image upsampling[END_REF] proposed multi-level modified finite radon transform network (MMFRTN) for image upsampling and Lei et al. [START_REF] Lei | Learning stacked image descriptor for face recognition[END_REF] proposed stacked image descriptor for face recognition. Li et al. [START_REF] Li | Feature learning based on SAE-PCA network for human gesture recognition in RGBD images[END_REF] proposed SAE-PCA network for human gesture recognition in RGBD (Red, Green, Blue, Depth) images. Zeng et al. [START_REF] Zeng | Tensor object classification via multilinear discriminant analysis network[END_REF] proposed a multilinear discriminant analysis network (MLDANet) for tensor object classification.

On the other hand, color image processing also attracts the attention of many researchers since color images can provide much larger amount of information of the real-world objects when compared with gray images [START_REF] Van De Sande | Evaluating color descriptors for object and scene recognition[END_REF][START_REF] Zhang | Simple yet effective color principal and discriminant feature extraction for representing and recognizing color images[END_REF][START_REF] Khemchandani | Color image classification and retrieval through ternary decision structure based multi-category TWSVM[END_REF][START_REF] Zheng | Image segmentation by generalized hierarchical fuzzy C-means algorithm[END_REF][START_REF] Wen | A rapid learning algorithm for vehicle classification[END_REF]. The classical representation of color images usually combines the values of Red (R), Green (G), Blue (B) -channels into one vector. Under this simple representation, the relationships between each color channel pixel of color image are destroyed, and the dimension of images is three times of that of gray-scale images. As a result, it is crucial to seek a way to represent color images more properly by taking the spatial relationships between R, G, B channels into consideration. Quaternion [START_REF] Hamilton | On quaternions[END_REF][START_REF] Zhang | Quaternions and matrices of quaternions[END_REF] is shown as powerful mathematical tool for representing color images [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF][START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF][START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF][START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF][START_REF] Guo | Multifocus color image fusion based on quaternion curvelet transform[END_REF][START_REF] Chen | Color image analysis by quaternion-type moments[END_REF][START_REF] Guo | Quaternion Fourier-Mellin moments for color images[END_REF][START_REF] Gai | New banknote defect detection algorithm using quaternion wavelet transform[END_REF] and also rotation operation [START_REF] Kuipers | Quaternions And Rotation Sequences: A Primer With Applications To Orbits, Aerospace And Virtual Reality[END_REF][START_REF] Liwicki | Full-angle quaternions for robustly matching vectors of 3D rotations[END_REF]. In this way, a color image can be represented as a compact quaternion, which preserves intrinsic structure of color image and the spatial relationships between R, G, B channels. Bihan and Sangwine [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF] as well as Pei et al. [START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF] proposed a new quaternion feature extraction method for color image, called quaternion principal components analysis (QPCA), which is able to extract more robust and informative features from color image than classical PCA. These studies emphasize that QPCA is superior to PCA in color image representation.

At the present stage, most of the above feature extraction algorithms (SIFT, LBP, PCANet, LDANet, etc.) were only designed for gray image classification and they will face many issues when simply be applied to color image classification, such as degrading performance, dimension curse and so on [START_REF] Shakhnarovich | Face recognition in subspaces[END_REF]. Thus, similar researches that are suitable for color image classification need to be studied.

In this paper, to address the above issues, we propose a new color image feature extraction algorithm, namely, quaternion principal component analysis network (QPCANet) which is a quaternion deep learning architecture and extends the principals of PCANet approach from real domain to quaternion domain. The input of QPCANet utilize the quaternion expression of color image, where the values of R, G, B, channels are put into three imaginary parts of quaternion, respectively. In the processing stage of QPCANet, we first make use of quaternion convolutional layer with QPCA filter to extract higher level semantic quaternion feature of color image. We then use quaternion binary hashing operation to construct the nonlinear layer. For the last stage, i.e. pooling layer, the block histogram operation for four components of quaternion feature, which is extracted from the above layer, are processed to generate the feature vector for classification. The performance of QPCANet feature on various classifiers, for example, k-nearest neighbor (KNN), large scale linear support vector machine (large scale SVM) [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF][START_REF] Gu | Incremental Support Vector Learning for Ordinal Regression[END_REF][START_REF] Gu | A Robust Regularization Path Algorithm for ν-Support Vector Classification[END_REF], are also studied in detail in this paper. Moreover, QPCANet are then evaluated and compared with PCANet on various color image databases for many classification tasks, including face recognition, object recognition, texture classification and land use classification.

The contributions of the paper are as follows. First, we propose a novel feature extraction algorithm, namely, QPCANet, which shows better performance on color image classification than that of RGB PCANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF], RGB LDANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF], RGB SIFT [START_REF] Van De Sande | Evaluating color descriptors for object and scene recognition[END_REF], and RGB LBP [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]; Second, we show that, despite its simplicity and generality of quaternion architecture, one layer QPCANet achieve better result than two layer PCANet. Last but not least, to our best knowledge, we first extend simple deep learning network to quaternion domain. Crucially, this work gives a totally new perspective of deep learning neural networks on quaternion domain. We empirically show that even a low-level hierarchy convolutional networks in quaternion domain may outperform high-level convolutional networks in real domain.

The rest of the paper is organized as follows. The quaternion algebra and the QPCA are briefly introduced in Section 2. Section 3 describes the architecture of QPCANet in details and also shows how it works. Section 4 introduces similar methods proposed by existing literature and their parameter setting in experiments. The parameter model of QPCANet is also discussed in this section. Classification performances of QPCANet and PCANet are evaluated and compared on various color image datasets in Section 5. Section 6 concludes the paper.

Quaternion algebra and QPCA

In this section, we briefly review the quaternion algebra and QPCA.

Quaternion algebra

A quaternion number x is a hypercomplex number, which consists of one real part and three imaginary parts:

( ) ( ) ( ) ( ) x S x I x J x K x      i j k  (1)
where  denotes the quaternion number field, ( ), ( ), ( ), ( ) S x I x J x K x   (  denotes the real number field) and i , j and k are three imaginary units obeying the following rules:

2 2 2 1, , ,           i j k ijk ij ji k jk kj i ki ik j (2)
, which shows that the quaternion multiplication is non-commutative. We call x a pure quaternion when ( ) 0 S x  .

The conjugate and l 2 norm of a quaternion x is defined as, respectively,

* ( ) ( ) ( ) ( ) x S x I x J x K x     i j k (3) 2 2 2 2 ( ) ( ) ( ) ( ) x S x I x J x K x     (4)
If 1 x  , we call x a unit quaternion. For a complete review of the properties of quaternion, please refer to [START_REF] Zhang | Quaternions and matrices of quaternions[END_REF].

QPCA

A color image can be represented as a pure quaternion matrix

m n   Q 
whose elements are as follows [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF]:

( , ) ( , ) ( , ) ( , ) , 1 ,1 R G B s t st st st s m t n        Q Q Q Q i j k (5) Where ( , ) R s t Q , ( , ) G s t Q
, and ( , ) B s t Q are represented as red, green, and blue components of the pixel at position s, t), respectively. Suppose we have a set of quaternion images { , 1,..., }

m n i i N    Q 
. We vectorize all quaternion images and denote them as { , 1,..., }

mn i i N   Q  , which are then concatenated to obtain 1 2 [ , , , ] mn N N    Q Q Q Q   (6)
Then, the covariance matrix is given by

1 H N    C QQ (7)
where

[ ] E   Q Q Q = , [ ]
E Q is a matrix whose columns are centered from Q and the superscript H is the conjugate transposition operator. Then, the QPCA of color images is given by: H  C WΩW [START_REF] Bengio | Representation learning: a review and new perspectives[END_REF] Note that C is a quaternionic Hermitian matrix, Ω is a real diagonal matrix [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF] and each column of

mn mn   W 
is an eigenvector of C according to quaternion eigenvalue decomposition [START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF]. Then quaternion feature of color images can be written as:

H f  Q W Q (9)
Each column of f Q is a quaternion feature corresponds to each color image.

One may note that QPCA is invariant to spatial rotation of color image shown in [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF].

However, if you concatenate the values of R, G, B channels, and then use the ordinary PCA, the process is not invariant to spatial rotation of color images. That is why QPCA is superior to the ordinary PCA when dealing with color images.

Quaternion principal component analysis network (QPCANet)

The architecture of the proposed QPCANet is depicted in Fig. 1. In this section, we analyze the structure of two-layer QPCANet (QPCANet-2) for color image classification, and we show how to build a multi-layer QPCANet.

The first QPCA layer

Suppose that we have

N quaternion images { , 1,..., } m n i i N    Q 
and their corresponding labels for training. For simplicity, we assume that the patch size is k1  k2, k1 and k2 are positive odd integers. We collect all (mk11)  (nk21) quaternion patches around each pixel of the ith quaternion image i Q . Then each quaternion patch is centered by subtracting its mean. Thus we get zero-mean quaternion patches for the ith quaternion image. We note these quaternion features

[ , , , ]

i i i i m k n k        q
q q q , and each column of i q belongs to 1 2 k k  . Repeating the above process, we can get all quaternion patches of N input patterns for training. By constructing the same matrix for all quaternion patches and by putting them together, we obtain

1 2 1 2 ( 1 ) ( 1 ) 1 2 [ , , , ] k k N m k n k N         q q q q  (10)
The covariance matrix of q is computed as:

1 1 2 ( 1 ) ( 1 )
H N m k n k     
qq C [START_REF] Olshausen | Emergence of Simple-Cell Receptive Field Properties by Learning a Sparse Code for Natural Images[END_REF] The above matrix admits a quaternion eigenvalue decomposition:

1 1 1 H  C W ΩW (12)
where

1 2 1 2 1 k k k k   W 
is a unitary matrix that contains the eigenvectors of 1 C and

1 2 1 2 k k k k   Ω 
is a diagonal matrix with eigenvalues on its diagonal. The values on the diagonal of Ω are arranged in decreasing magnitude order and the corresponding eigenvectors, i.e. the principal components of q , are arranged accordingly in 1 W . The bigger the eigenvalue is, the more important the quaternion principal component will be. These principal component vectors are also called QPCA filter bank.

Let i L be the desired number of QPCA filters in the i th layer. We choose then the first 1

L eigenvectors of 1 W to form a new matrix 1 2 1 k k L   W 
where each of its columns is seen as a

filter 1 l W : 1 2 1 , 1,2, , k k l l L    W   (13)
These filters form the filter bank of quaternion image { , 1,..., }

i i N  Q
. The filter bank captures the main variations of all quaternion patches.

Then the first QPCA layer executes a convolution using the QPCA filter bank W 1 to get a set of quaternion feature maps

1 { , 1,..., } l i l L  F of i Q as follows: 1 1 , 1,2, , , 1,2,
,

l i i l i N l L     F Q W   (14)
where can also be used as a quaternion input pattern. By repeating the above process, one can expect to derive high-level features.

1 { , 1,..., } l i l L  F is the lth feature map of i Q .The

The second QPCA layer

Repeating the same process as the first stage, we collect all the quaternion patches which have the same size as the first layer from all the quaternion feature maps of all images. We denote these quaternion patch as a patch set

, 1 1 2 { , 1,..., ; 1,...,( 1)( 1)} l a i l L a m k n k       r ,
L L L L L m n i i i i i i   G G G G G G     as follows: , 2 1 2
{ , 1,..., ; 1,..., }

l l i i l L L     G F W    (15)
One or more additional stages can be stacked as above. But in our experiments, we found that QPCANet-2 is enough to get a good accuracy in color image classification and thus a deep architecture is not necessarily required.

Quaternion feature maps coding

Each quaternion feature maps extracted by the second QPCA layer should be binarized, weighted, and then summed in order to reduce the complexity of quaternion feature maps and make the extracted feature more separable. Thus, each quaternion feature

l l l l l i i i i i S I J K     G G G G G i j k      map , , , , , ( ) ( ) ( ) ( ) ,
is binarized by applying the Heaviside step function H) to its four parts (the value of H) is one for positive entries and zero otherwise). The binarized quaternion feature map is denoted by

, , , , , ( ) ( ) ( ) ( )
l l l l l i i i i i S I J K     G G G G G           i j k,
and the resulting maps are weighted to form a new quaternion pattern:

        2 1 , 1 2 L l l l l l l i i i i i i S I J K         T G T T T T     i j k . ( 16
) Note that the pixel values of         , , , l l l l i i i i S I J K T T T T are integers belonging to the interval 2 0, 2 1 L      .
Next, we discuss how to code these binarized quaternion feature maps

l i T . Since l i T has four parts         , , , l l l l i i i i S I J K T T T T ,
      , , l l l i i i I J K T T
T , we obtain four vectors, which are then concatenated to derive the feature vector of the quaternion image Q i , which can be written as follows:

                2 1 4 (2 )
hist , hist , hist , hist

L L B l l l l i i i i i B S B I B J B K        f T T T T  (12)
Since the feature of quaternion image is computed by histogram from a large amount of blocks which picked from several feature maps, the feature vector of i Q is sparse and very long. Thus we feed these feature vectors to classifies, for example, KNN, SVM [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF][START_REF] Gu | Incremental Support Vector Learning for Ordinal Regression[END_REF][START_REF] Gu | A Robust Regularization Path Algorithm for ν-Support Vector Classification[END_REF], for classification.

Chan et al. [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF] showed that large scale linear SVM [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF] can obtain higher accuracy than other classifiers when dealing with the sparse and long feature vector, and they also suggested that, when using PCANet, non-overlapping blocks are suitable for face images, whereas the overlapping blocks are appropriate for hand-written digits, textures, and object images. We found that the recommendation is also suitable for QPCANet.

Experimental Setting and Parameters model selecting

Introduction of the experimental datasets and the preprocessing

In this subsection, we introduce the datasets and the preprocessing used in our experiments.

UC Merced Land Use dataset. UC Merced land use dataset [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF] is manually extracted from large images from the United States Geological Survey (USGS) National Map Urban Area Imagery collection for various urban areas around the country. The dataset, composed of 256 × 256 pixels RGB images with pixel resolution of one foot, is manually partitioned into 21 classes.

Each class contains 100 images. Fig. 2 shows 21 images, each of which represents a distinct class.

Some representative images from the dataset are also shown in Fig. 3, from which we can see that the dataset contains abundant rotation information due to a different point of view of aerial photography. Thus, the land objects appear in one class are mainly the same, the only difference is the angles of view. Obviously, one may see that there is a large intra-class variability, which will increase the difficulty of color image classification. Every images are scaled to 32 × 32 pixels to reduce the computational cost in our experiment.

Georgia Tech face database. Georgia Tech face database [START_REF] Nefian | Georgia tech face database[END_REF] Caltech-101 object dataset. Caltech-101 dataset [START_REF] Li | Learning generative visual models from few training examples: An incremental bayesian approach tested on 101 object categories[END_REF] contains a total of 9146 images including color and gray-level images, belonging to 101 distinct objects, including faces, watches, ants, pianos, etc. In our experiment, we only consider 8733 color images in Caltech-101 dataset and discard all gray-level images. All color images were resized to be 32 × 32 pixels in order to reduce the computational time.

Introduction of our compared methods

In this paper, similar with PCANet, we only compare the results of single descriptor but not compare the methods that are composed of many different descriptors. Specifically, we compare our proposed QPCANet method with original gray PCANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF], and its two variants, that is, RGB PCANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF], and RGB LDANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF], and also with some widely used feature extractors, including RGB SIFT [START_REF] Van De Sande | Evaluating color descriptors for object and scene recognition[END_REF] and RGB LBP [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. The compared methods are described as follows.

With respect to RGB SIFT, we first perform conventional SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] to R, G, B channels separately to obtain three feature vectors which are then concentrated together to get one feature vector. RGB SIFT shows better performance than conventional SIFT when dealing with color images. The conventional SIFT descriptor is extracted from a 16×16 block and the grid spacing is set to 6. Following the method proposed by Yang et al. [START_REF] Yang | Linear spatial pyramid matching using sparse coding for image classification[END_REF], we computer sparse coding for each SIFT feature vector of all images and then histogram by three pyramid levels. For RGB LBP, conventional LBP are separately performed with respect to R, G, B channels to obtain three vectors, which are then combined to form one feature vector. The number of filter size in LBP is set to 8 and the radius is set to 1. When dealing with RGB PCANet, we use the same method as [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF] to preprocess color image, i.e. we gather three color channels together to form a long vector representation of color image and then regard it as a gray image to process it. Because RGB LDANet has the same architecture as RGB PCANet, so we directly use the same parameter model and preprocessing method. The postfixes -1 -2 appeared in later sections denote the number of layers in network architecture (e.g. RGB LDANet-1, RGB LDANet-2).

The impact of classifiers and various parameters

In this subsection, we focus on the problem of classifier selection and parameter models setting for QPCANet by using UC Merced land use dataset.

The impact of classifiers. We first discuss which classifier is more suitable for the feature extracted by the proposed QPCANet. We test three kinds of classifier, including KNN with Euclidean distance, KNN with cosine similarity, and large scale linear SVM. The parameters of two-layer QPCANet (QPCANet-2) and two-layer RGB PCANet (RGB PCANet-2) are set as follows: the patch size, the block size, the overlapping ratio, and the number of filters are 3×3, 8×8, 0.6, and L 1 = L 2 = 8, respectively. The recognition accuracy of QPCANet-2 and RGB PCANet by using three classifiers is shown in Table 1. One can see that the performance of RGB PCANet-2 is slightly better than that of QPCANet-2 by using KNN classifier, however, the performance of both QPCANet-2 and RGB PCANet-2 is not very good in KNN classifier in general.

The reason of experimental results is that KNN with cosine similarity outperform that with Euclidean distance in high dimensional vector space, especially for sparse entries [START_REF] Levent | Finding clusters of different sizes, shapes, and densities in noisy, high dimensional data[END_REF][START_REF] Strehl | Impact of similarity measures on web-page clustering[END_REF]. One may note that a feature vector generated by QPCANet is a combination of histograms computed statistically from blocks, thereby possessing lots of sparse entries. Considering Euclidean distance of two vectors are computed from the difference of each entry within both vectors, this sort of measure cannot be avoided to assign same weight to every entry of vector. The sparse entries in a high dimensional vector should be assigned smaller weight because they do not contain abundant information as much as other entries possessed normal value. Thus, the capacity of classification will decay when apply to high dimensional feature vector. On the contrary, cosine similarity pays more attention to the estimation of difference between two angles of two feature vectors, it means that cosine similarity typically focuses on principal components of feature vector to discriminate two vectors, rather than trivial entries.

The performance of QPNCANet-2 by using large scale linear SVM is much better than that of RGB PCANet. Therefore, a proper classifier significantly promotes the performance of QPCANet-2. Why large scale linear SVM is more suitable for QPCANet-2 feature and RGB PCANet-2 feature than KNN? Because the length of RGB PCANet-2 feature is three times of that of gray PCANet-2 and the length of QPCANet-2 feature is four times of that of gray PCANet-2.

Both QPCANet-2 and RGB PCANet-2 will obtain a very long but very sparse feature vector from color image, which lead to very worse performance by using KNN classifier but good performance by using large scale linear SVM. Therefore, we use large scale linear SVM as classifier in the following experiments.

The impact of the number of filters L i . We then discuss the effect of the number of filters L i on the recognition accuracy of QPCANet, whose patch size, block size, and overlapping ratio are set to 3×3, 8×8, and 0.6, respectively. Then we use a greedy algorithm to find the optimal number of QPCA filters in every stages of QPCANet. Fig. 5 shows the result of the recognition accuracy of one-staged QPCANet (QPCANet-1) in relation with the number of filters L 1 , whose values vary from 2 to 9. One can see that QPCANet-1 achieves the best result in L1 = 9. But we still choose L1 = 8 when taking the octave filters into consideration as that of gray PCANet [START_REF] Chan | PCANet: A simple deep learning baseline for image classification?[END_REF]. We then draw in Fig. 4 the result of the recognition accuracy of two-staged QPCANet (QPCANet-2) in relation with the number of filters L 2 , whose values vary from 2 to 9, meanwhile, L 1 is set to 8. It can be seen from Fig. 4 that QPCANet-2 achieves a good performance for L 2 = 8, however, when L 2 ≥ 8, the recognition rate begins decreasing. Considering the simplicity of applications of QPCANet, we should make the parameter model as simple as possible. Thus, although some fine-tuned L 1 and L2 values could lead to performance improvement, we decided to set L1 = L2 = 8.

The impact of the patch size and the block size. We then study the effect of the patch size and the block size on the recognition accuracy of QPCANet, whose number of filters is set to L 1 = L 2 = 8 and overlapping ratio is set to 0.6, respectively. Similar to gray PCANet, the patch size is set to 3×3, 5×5, 7×7. We set the block size to 4×4, 8×8, 16×16, which are the submultiple divisors of the original image size. The results are shown in Table 2. One can see that the best recognition rate is obtained when the patch size is 3×3 and the block size is 8×8. We can also note that the accuracy in UC Merced land Use dataset declines slowly with the increase of the patch size, which means that appropriate patch size and block size should be chosen in QPCANet. Generally speaking, the smaller the patch size is, the finer the scale of the image feature is. Furthermore, modest block size can capture enough energy to represent image in the coding layer of QPCANet.

We then study the effect of the patch size and the block size on the recognition accuracy of QPCANet, whose number of filters is set to L1 = L2 = 8 and overlapping ratio is set to 0.6, respectively. Similar to gray PCANet, the patch size is set to 3×3, 5×5, 7×7. We set the block size to 4×4, 8×8, 16×16, which are the submultiple divisors of the original image size.

We then research the effect of the overlapping ratio on the recognition accuracy of QPCANet, whose patch size, block size, and number of filters are set to 3×3, 8×8, and

L 1 = L 2 = 8,
respectively. The overlapping ratio varies from 0.1 to 0.9. Fig. 6 demonstrates the recognition rate in this experiments. One can see that the recognition rates of different overlapping ratio fluctuate at 78%. We choose overlapping ratio 0.6 to get the best performance in UC Merced Land Use.

The parameter model is concluded in Table 3.

Experimental Results and Discussion for different image datasets

Testing on UC Merced land use database

For the UC Merced land use database [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF], we randomly selected 80 training images per class, and the remaining ones were used for testing. The QPCANet is trained with the number of filters

L 1 = L 2 = 8
, the block size is set to 8 × 8 and the overlapping ratio is fixed to 0.6. We changed the patch size k1 × k2 from 3 × 3 to 7 × 7. The same parameters are used for RGB PCANet-2 and Gray PCANet. The results are reported in Table 4. Note that with the increase of patch size, the recognition accuracy of all networks decreases. The best performance for various patch sizes is highlighted in bold. In UC Merced Land use dataset, one-staged QPCANet (QPCANet-1) outperforms all the other one-staged networks. It is even better than some two-staged networks.

Compared to RGB PCANet-2, the proposed QPCANet-2 increases the recognition accuracy by more than 6% in the case of patch size 3 × 3. The performance of Gray PCANet is worse than that of RGB PCANet and QPCANet in all cases.

The best performance of one-staged RGB LDANet (RGB LDANet-1), two-staged RGB LDANet (RGB LDANet-2), RGB SIFT, and RGB LBP are also given in Table 5. It turns out that QPCANet performs much better than RGB PCANet, gray PCANet, RGB LDANet, RGB SIFT and RGB LBP in UC Merced Land Use database. The reason may be due to the nature of UC Merced land use database which contains a large amount of rotation type images. Quaternion representation gives a concise and all-in-one representation for rotation, and is known for its successful application in computer graphics, computer vision, and orbital mechanics of satellites [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF][START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF][START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF][START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF][START_REF] Guo | Multifocus color image fusion based on quaternion curvelet transform[END_REF][START_REF] Chen | Color image analysis by quaternion-type moments[END_REF][START_REF] Guo | Quaternion Fourier-Mellin moments for color images[END_REF][START_REF] Gai | New banknote defect detection algorithm using quaternion wavelet transform[END_REF][START_REF] Kuipers | Quaternions And Rotation Sequences: A Primer With Applications To Orbits, Aerospace And Virtual Reality[END_REF][START_REF] Liwicki | Full-angle quaternions for robustly matching vectors of 3D rotations[END_REF], etc. We suspect that the quaternion model of color image gives more rotation invariance to QPCANet when compared to RGB PCANet. In the next section, we further experimentally verify this conjecture.

features that are less sensitive to the rotation of color images. The performance of LDANet, RGB LBP, and RGB SIFT are not as good as QPCANet and RGB PCANet.

Texture discrimination in CURet Dataset

In our experiments, CURet database is randomly split into a training set and a testing set, with

23 training images for each class. The QPCANet is trained with the number of filters

L 1 = L 2 = 8,
the overlapping ratio is set to 0.5, and the block size is fixed to 8 × 8. The RGB PCANet and Gray PCANet use the same parameters as that of QPCANet.

The average recognition rates of QPCANet and PCANet over 5 different random splits are given in Table 10 and the performance of other methods whose parameter model shown in section 4.3 are given in Table 11. We can see that QPCANet outperforms other networks whether the stage is one or two. When the patch size is set to 7 × 7, the performance of QPCANet-1 approximately equals that of RGB PCANet-2 and outperform that of Gray PCANet-2. We also notice that RGB PCANet-1 outperforms Gray PCANet-1. The performance of RGB PCANet is much higher than that of Gray PCANet due to the use of color information and boundary information of texture dataset at the same time. The performance of RGB LDANet is neither good nor bad. RGB SIFT perform not well in texture database due to inconspicuous key points.

RGB LBP outperform RGB SIFT slightly on this texture database.

Performance in Caltech-101

In Caltech-101 dataset, we randomly select 15 color images as training images per class and the remaining images are used for testing. With respect to the parameters setting of QPCANet, the overlapping ratio of block is set to 0.5 to preserve appropriate feature invariance for objects in Caltech-101 database.

For fair comparisons, appropriate parameters of RGB PCANet and gray PCANet are required.

For both RGB PCANet and Gray PCANet, we enlarged patch size from 3×3 to 7×7. The block size was fixed to 7×7 and overlapping ratio is set to 0.5. The average accuracy over 5 drawing of the training set is listed in Table 12 and the best accuracy rates of different methods whose use the showing parameter model is given in Table 13. We see that the best performance (63.78%) is achieved by QPCANet-2 when patch size is 3 × 3. With the increase of the patch size the accuracy decreases. This illustrates that too big patch size is not appropriate for color object recognition with QPCANet. We think there are at least two reasons: (1) generally speaking, the smaller patch size is, the finer scale of feature is. When the patch size increases, the invariance increases, but it comes with a partial loss of information [START_REF] Bruna | Invariant scattering convolution networks[END_REF], that is, a good choice of the patch size is the balances between the invariance enhancement and information loss. What's more important for the performance of QPCANet is that it should choose the small patch size to capture more details from color images, however, we do not worry too much about invariance since QPCA has enough rotation invariance for color images [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF][START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF]; (2) the convolutional layer of QPCANet is constructed by QPCA, which can capture enough energy to represent and classify color images even with very small patch sizes [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Bihan | Quaternion principal component analysis of color images[END_REF][START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF].

The performance improvement from QPCANet-1 to QPCANet-2 is not as large as that of RGB PCANet. But we note that the accuracy of QPCANet-1 is much better than other methods. In fact, the behavior of QPCANet-1 is very close to that of Gray PCANet-2 and RGB PCANet-2. This implies that we can achieve an acceptable result by just apply a simple QPCANet-1. Thus, QPCANet with only one layer still has competitive performance when comparing with other methods.

By comparing Gray PCANet with RGB PCANet, a strange phenomenon draws our attention.

Why the accuracy of RGB PCANet-1 is lower than that of Gray PCANet-1, but RGB PCANet-2 outperforms Gray PCANet-2? First, let us consider some pictures extracted from Caltech-101 and reported in Fig. 8. We note that every horizontal pair of objects has the same main colors, but the objects are not belonging the same class. Thus, we cannot do a right classification only with color information. For color image classification, color information and boundary information seem equally important. To some extent, color information is harmful to color image classification when an algorithm has not enough representation ability. That is, color is a burden for some classification algorithms. It seems that RGB PCANet-1 is such an algorithm, which is not able to utilize simultaneously color information and boundary information. For Gray PCANet-1, we changed color images to gray-scale ones and thus only the boundary information, which is effectively utilized by Gray PCANet-1, is preserved. Therefore, it is not surprising that Gray PCANet-1 outperforms RGB PCANet-1.

Surprisingly, a significant improvement of performance happens when upgrade RGB PCANet-1 to RGB PCANet-2 by adding one additional layer. RGB PCANet-2 seems to have enough ability to deal with both color information and boundary information. When compared with RGB PCANet-1, the representation ability of RGB PCANet-2 is significantly higher. RGB PCANet-2 also performs better than Gray PCANet-2 by utilizing color information of objects.

Gray PCANet-2 does not perform significantly better than Gray PCANet-1 in this database. We can conclude that the accuracy of color object recognition is difficult to improve by only utilizing boundary information.

Why QPCANet outperforms the others methods (RGB PCANet, Gray PCANet, RGB LDANet, RGB SIFT and RGB LBP) in color image classification? We think that the reason is the utilization of quaternion model of color image. The quaternion representation of color image preserves the underlying spatial structures and relationships between R, G, B channels. QPCA is shown to be more suitable for color images representation, and allows enhancing the robustness of color images features [START_REF] Bihan | Quaternion principal component analysis of color images[END_REF][START_REF] Pei | Quaternion matrix singular value decomposition and its applications for color image processing[END_REF]. The shortcoming of RGB PCANet is that it neglects the structure and unity of color images.

To conceptualize the learned QPCANet filters of patch size 3 × 3, we draw them in Fig. All the filters depicted in Fig. 9 and 10 are scaled within [0, 1]. Some learned filters exhibit similar shapes however; their numerical values are not the same. We notice that some similar QPCA filters occur several times in the different parts of quaternion. Such a filter redundancy may improve the intra-class invariance of feature maps. For the first three rows in RGB PCA filters, we notice that the last seven RGB PCA filters of each row have the same shape although their numerical values are not the same. It means that the similar filtering operations are performed in R, G, B channels.

Conclusion

In

2 1 L

 21 So we need coding each part of l i T separately. Let us first consider ( ) l i S T for example. We divide ( ) l iS T into B blocks. For each block, we simply compute the histogram (with 2  bins) of the decimal values in each block. After this coding process, we concatenate the histograms of the B blocks into one vector, that is, the same coding operation on another three parts

  contains images of 50 people and for each individual, 15 color images were collected at the Center for Signal and Image Processing at Georgia Institute of Technology. All images in the database are 640 × 480 pixels and the average size of faces in these images is 150 × 150 pixels. Most of the images were taken in two different sessions to take into account the variations in illumination conditions, facial expression, and appearance. In our experiment, the face images are cropped with the dimension 128 × 160 to guarantee classification accuracy. Some examples are shown in Fig. 4. To reduce the computational complexity, we resize all cropped images to 64 × 64 pixels. CURet texture database. The CURet texture dataset [59] contains images of 61 materials that broadly span the range of different surfaces that we commonly see in our environment. Every class is composed of 92 images of the size of 200 × 200 pixels. All images are linearly scaled to 32 × 32 pixels for experimental simplicity.

  9. In this figure, the first four rows correspond to the real part and three imaginary parts of the first stage QPCA filters, and the others correspond to the four parts of the second stage QPCA filters. The learned RGB PCANet filters of patch size 3 × 3 are shown in Fig. 10. The first three rows represent the one stage RGB PCA filters where each row corresponds to R, G, B channels, respectively. The last row represents the PCA filters of the second stage.

Fig. 1 .

 1 Fig. 1. Architecture of two-stage QPCANet

Fig. 2 .Fig. 3 .Fig. 4 . 3 Fig. 5 .

 23435 Fig. 2. Some representative examples of classes in the UC Merced land use database

Fig. 6 .Fig. 7 .

 67 Fig.6. The overlapping ratio of QPCANet changes from 0.1 to 0.9

Fig. 8 .

 8 Fig. 8. Some pictures have the same main colors but are not belonging to the same category of objects

  where ,

	2 { , 1,..., } 2 L  W  	, which are then convoluted by l i F to obtain L1  L2 quaternion maps
	11	, ,	1	2	,	21	, ,	2	2	,	1	1	, ,	1 2
															l a r is the i
	ath quaternion patch of l i F . The set	, { , 1,..., ; 1 l a i l L a  r		1,...,(	1 m k  	1)(	2   n k	1)}	are then sent
	to QPCA algorithm to get the filter bank	2 W in the second layer.
	For the filter bank
															10

2 W , we only pick the first 2 L column to form a new filter bank

Table 9

 9 The accuracy (%) of different methods on Rotated Georgia Tech face database dataset

	Method	Recognition rates
	RGB SIFT	45.83
	RGB LBP	42.22
	QPCANet-1	84.94
	QPCANet-2	79.17
	RGB PCANet-1	80.56
	RGB PCANet-2	74.17
	Gray PCANet-1	69.44
	Gray PCANet-2	63.61
	RGB LDANet-1	86.94
	RGB LDANet-2	64.72

Table 10

 10 The accuracy (%) of texture classification on CURet dataset

	Patch Size	QPCANet-1 QPCANet-2	RGB PCANet-1	RGB PCANet-2	Gray PCANet-1	Gray PCANet-2
	3 × 3	96.94	98.66	92.28	96.36	90.71	93.70
	5 × 5	97.81	98.57	95.60	98.08	95.60	96.56
	7 × 7	98.05	98.40	96.29	98.03	96.03	96.20

Table 11

 11 The accuracy (%) of texture classification on CURet dataset

	Method	Recognition rates
	RGB SIFT	60.91
	RGB LBP	74.25
	QPCANet-1	98.05
	QPCANet-2	98.66
	RGB PCANet-1	96.29
	RGB PCANet-2	98.08
	Gray PCANet-1	96.03
	Gray PCANet-2	96.56
	RGB LDANet-1	94.61
	RGB LDANet-2	95.58

Acknowledgement

This work was supported by the National Natural Science Foundation of China (No. 61201344, 61271312, 61401085, 61572258, 11301074), and by the Project Sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, by the Qing Lan Project and the '333' project (No. BRA2015288), and by the Open Fund of China-USA Computer Science Research Center (KJR16026). The authors are also grateful to the anonymous reviewers for their constructive comments and suggestions to greatly improve the quality of this work and the clarity of the presentation.

Face recognition on Georgia Tech face database

For the Georgia Tech face database, we randomly pick up 10 images in each class as training images, and the others are used for testing images. We find that one-layer networks provide excellent results. Therefore, two-layer networks are not considered here. The number of filters L1 in networks is set to 8 and the non-overlapping block is of size 8 × 8. The average accuracy of classifications over 10 times experiments are listed in Table 6.

We also compare the best performance of different methods in Table 7. We can see from the table that QPCANet performs better than RGB PCANet, gray PCANet, RGB LDANet, RGB SIFT and RGB LBP in Georgia Tech face database.

Next, we experimentally verify the hypotheses which we find in the experiment of UC Merced Land use database: QPCANet outperforms RGB PCANet when many rotated color images are contained in the database.

Firstly, we choose only one cropped front face image from 20 individuals. These images are zero padded to form a square image of the size 220 × 220 pixels. A representative padded image is shown in Fig. 7(a). The zero-padding process ensures that the rotation operation will not induce any border effect. For each individual (i.e. class), we rotate the padded front image from 0° to 360° by a step of 10° to obtain 36 images for each class. Some rotated images are shown in Fig. 7(b). We randomly select 18 images from each class to train networks and the remaining ones are used in the testing procedure. All images are then linearly scaled to 64 × 64 pixels. For QPCANet and PCANet, the number of filters, and the (non-overlapping) block size are set to L 1 = L 2 = 8, and 8×8, respectively. The recognition rates of QPCANet, RGB PCANet, and gray PCANet in terms of different parameter models are shown in Table 8.

The results of best performance of different methods shown in section 4.3 are listed in Table 9.

Surprisingly, the accuracy of all networks decreases with the increase of stage. When the patch size is 5 × 5, the recognition accuracy of QPCANet-1 reaches 84.94%, which outperforms all the other methods. We can conclude that, compared with RGB PCANet, QPCANet allows to extract


Benoît Gaüzère

Sébastien Bougleux

Luc Brun

Approximating Graph Edit Distance using GNCCP

The graph edit distance (GED) is a flexible and widely used dissimilarity measure between graphs. Computing the GED between two graphs can be performed by solving a quadratic assignment problem (QAP). However, the problem is NP complete hence forbidding the computation of the optimal GED on large graphs. To tackle this drawback, recent heuristics are based on a linear approximation of the initial QAP formulation. In this paper, we propose a method providing a better local minimum of the QAP formulation than our previous proposition based on IPFP. We adapt a convex concave regularization scheme initially designed for graph matching which allows to reach better local minimum and avoids the need of an initialization step. Several experiments demonstrate that our method outperforms previous methods in terms of accuracy, with a time still much lower than the computation of a GED.

Introduction

Graphs provide a flexible framework to represent data including relationships between elements. In addition, graphs come with an underlying powerful theory which allows to infer a lot of information from this representation. However, graph's space is not an euclidean space. This last point avoids the use of classic machine learning methods mainly designed to operate in euclidean spaces. Several approaches aim to bridge the gap between graph space and euclidean spaces in order to combine machine learning and graphs. A first historical approach consists in embedding graphs onto euclidean spaces by computing a set of descriptors describing graphs. Even if this method is straightforward and can be easily controlled by choosing the information to keep, the loss of structural information induced by the euclidean embedding may constitute a major drawback for some applications. An extension of this approach is based on the well known kernel trick. Kernel methods implicitly define a scalar product between embedding of graphs in some Hilbert spaces, without requiring an explicit formulation of the Reproducing Kernel Hilbert Space (RKHS) associated to the kernel. The use of this trick allows to combine graphs and powerful kernel methods such as Kernel Ridge Regression or Support Vector Machines. However, constraints on the design of graph kernels often complexify the consideration of fine similarities or dissimilarities between graphs.

Another strategy consists in operating directly in the graph space. One of the most used measure is the graph edit distance (GED) [START_REF] Bunke | Inexact graph matching for structural pattern recognition[END_REF][START_REF] Sanfeliu | A distance measure between attributed relational graphs for pattern recognition[END_REF]. The GED of two graphs may be understood as the minimal amount of distortion required to transform a source graph into a target one. This distortion is encoded by an edit path, defined as a sequence of edit operations which includes nodes and edges substitutions, removals and insertions. Depending on the context, each edit operation e included in an edit path γ is associated to a non-negative cost c(e). The sum of all edit operation costs included within γ defines the cost A(γ) associated to this edit path. The minimal cost A(γ ) among all edit paths Γ (G 1 , G 2 ) defines the GED between G1 and G2:

GED(G 1 , G2) = A(γ ) = min γ∈Γ (G1,G2) A(γ) (1)
with A(γ) = e∈γ c(e). The edit path γ corresponds to an optimal edit path. GED has been widely used by the structural pattern recognition community [START_REF] Serratosa | Speeding up fast bipartite graph matching through a new cost matrix[END_REF][START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF][START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Riesen | A novel software toolkit for graph edit distance computation[END_REF] despite the fact that such distance comes along with several drawbacks. First of all, computing the GED of two graphs requires to find a path having a minimal cost among all possible paths, which is a NP-complete problem [START_REF] Lin | Hardness of approximating graph transformation problem[END_REF].

Computing an exact GED is generally done using A algorithm. In practice, due to its high complexity, the computation of an exact GED is intractable for graphs having more than 10 nodes [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF]. Such a limitation restricts applications of GED on real datasets, hence motivating the graph community to focus on heuristics providing suboptimal solutions of Eq. 1.

In this paper, we propose a method to compute an accurate GED approximation. In Section 2, we first show the close relationship between GED and nodes' mappings. This relationship leads to the formulation of the graph edit distance as a quadratic assignment problem (QAP). Then, Section 3 first reviews the method used in [START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF] to find a local minimum of the QAP associated to GED. Then, we propose a more accurate and reliable optimizer of our QAP formulation. Section 4 shows the effectiveness of our proposal on chemoinformatics problems.

GED as a quadratic assignment problem

Exact GED computation is based on a tree search algorithm which finds an optimal edit path among possible ones. The resulting algorithm follows the intuition given by the formal definition of the GED. However, this approach has an exponential complexity with the number of nodes.

Another formulation of the GED is based on its relationship with nodes' mapping. First, let us consider two sets of nodes

V 1 = {v 1 , . . . , v n } and V 2 = {u 1 , . . . , u m } of two graphs G 1 and G 2 , with n = |V 1 | and m = |V 2 |. The substi- tution of node v i ∈ V 1 to node u j ∈ V 2 can
be encoded by mapping v i to u j . Insertions/removals of nodes can not be encoded by a node to node mapping since a removed node will no longer appear in the target graph. Therefore, we augment the two sets V 1 and V 2 by adding enough null elements to encode the removal or insertion of any node:

V 1 = V 1 ∪ { 1 , . . . , m } and V 2 = V 2 ∪ { 1 , . . . , n }. Note that |V 1 | = |V 2 | = n + m.
Insertion of a node u j ∈ V 2 can now be represented by mapping an element to u j . In the same way, a mapping v i → encodes the removal of node v i . It has be shown [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF][START_REF] Neuhaus | A quadratic programming approach to the graph edit distance problem[END_REF] that, under mild assumptions on edit paths, each mapping between the two sets V 1 and V 2 of two graphs G 1 and G 2 corresponds to one and only one edit path transforming G 1 into G 2 . Considering this bijective relationship, finding the optimal edit path relies on finding an optimal mapping between the two sets V 1 and V 2 which minimizes the total mapping cost of nodes and edges.

Let us consider a bijective function φ :

V 1 → V 2 .
The mapping cost induced by this mapping can be defined as a sum of two terms:

S(V 1 , V 2 , φ) = L v (V 1 , V 2 , φ) + Q e (V 1 , V 2 , φ) (2)
The first term L v (V 1 , V 2 , φ) encodes the cost induced by nodes' mappings and the second term Q e (V 1 , V 2 , φ) the cost induced by edges' mappings. Let us consider the matrix X ∈ Π, with Π representing the set of binary doubly stochastic matrices: X encodes a mapping function Φ iff X i,j = 1 with φ(i) = j and X i,j = 0 otherwise. X corresponds then to a mapping or permutation matrix. The cost induced by nodes mapping can then be defined as:

L v (V 1 , V 2 , φ) = n i=1 c(v i → φ(v i )) + m i=1 c( i → φ( i )) = n+m i=1 n+m j=1 C i,j X i,j (3)
with

C ∈ R (n+m)×(n+m) +
encoding costs associated to edit operations on nodes. The bipartite GED, proposed by [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF], can be understood as an approximation of the problem expressed in Eq. 2. In bipartite approach, only the linear term

L v (V 1 , V 2 , φ) is considered. Finding the optimal mapping w.r.t. L v (V 1 , V 2 , φ)
corresponds to a Linear Sum Assignment Problem (LSAP) defined as:

X = argmin X∈Π vec(C) T vec(X) (4)
where vec(X) encodes the vectorization of matrix X, i.e. the concatenation of its rows into one single vector. For a sake of clarity, vec(X) will be denoted x.

Eq. 4 can be resolved using well known algorithms such as Hungarian [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF] in O((n + m) 3 ). The optimal mapping encoded by X encodes a set of nodes edit operations. Using simple graphs, edges' edit operations are inferred from the nodes' mapping to complete the edit path. Finally, the sum of costs associated to this edit path is taken as the approximation of GED. Note that since this edit path may not be optimal, the computed GED may be an overestimation of the exact one.

Each entry C i,j of matrix C encodes the cost of the edit operation induced by mapping the i-th element of V 1 to the j-th element of V 2 . Table 1 summarizes the general term C i,j . Computing an optimal mapping using only C as defined in Table 1 will only take into account information about nodes, ignoring all the structural information encoded by edges. However, since the quality of the GED

∈ V1 ∈ V2 substitution of vi by uj c(vi → uj) ∈ V1 ∈ V2 removal of vi c(vi → ) ∈ V1 ∈ V2 insertion of uj c( → uj)) ∈ V1 ∈ V2 none 0
approximation depends directly from the computed mapping, costs encoded in matrix C may also estimate costs induced by edit operations on edges. To include this information, methods based on bipartite graph matching [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF][START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF][START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF] define a cost matrix C augmented with the costs induced by the mapping of nodes' neighborhoods. The differences between different methods are mainly based on the size of the radius considered around each node. Results presented by [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF] have shown that, as expected, the quality of the approximation increases as long as we take into account a larger radius around each node, and thus more structural information. However, we also observe an asymptotic gain for radius greater than 3, hence showing the limit of this method. Obviously, computational times also increase as we increase the radius. Cost induced by edges' edit operations are inferred from edges' mapping. However, mapping φ is only defined between two extended sets of nodes and an edge is defined by a couple of nodes. Therefore, edge mappings can be directly deduced from the mapping of their incident nodes. Let e = (v i , v j ) ∈ E 1 , its mapped edge in V 2 corresponds to e = (φ(v i ), φ(v j )). If e ∈ E 2 , then the mapping corresponds to an edge substitution, else if e / ∈ E 2 , then the edge has been removed. Reciprocally, if e = (u i , u j ) ∈ E 2 and e = (φ -1 (u i ), φ -1 (u j )) / ∈ E 1 , then this edge has been inserted. The cost associated to edges' operations is encoded by a sum of three terms, one for each kind of edit operation:

Q e (V 1 , V 2 , φ) = (i,j)∈E1,φ(i)=k, φ(j)=l,(k,l)∈E2 c((i, j) → (k, l)) substitutions + (i,j)∈E1,φ(i)=k, φ(j)=l,(k,l) ∈E2 c((i, j) → ) removals + (k,l)∈E2,i=φ -1 (k), j=φ -1 (l),(i,j) ∈E1 c( → (k, l)) insertions (5)
Following the quadratic assignment formulation given by [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF], the term Q e (V 1 , V 2 , φ) can be expressed as a quadratic term depending on a matrix D ∈ R (n+m) 2 ×(n+m) 2 encoding edges' mapping costs :

Q e (V 1 , V 2 , φ) = x T Dx (6) Table 2. General term of matrix D (i, j) (k, l) edit operation D ik,jl ∈ E1 ∈ E2 substitution of (i, j) by (k, l) c((i, j) → (k, l)) ∈ E1 ∈ E2 removal of (i, j) c((i, j) → ) ∈ E1 ∈ E2 insertion of (k, l) into E1 c( → (k, l)) ∈ E1 ∈ E2 none 0
Table 2 summarizes all mapping costs encoded within matrix D. Note that the computation of x T Dx has, in a naive computation, a O((n + m) 4 ) computational and memory complexity. However, the vector x ∈ {0, 1} (n+m) 2 is sparse since it encodes a mapping and thus only (n + m) terms differs from 0, which reduces the computational complexity if we only process them (Alg. 2, lines 1 and 3). Therefore, with a proper approach, the computation of (x T Dx) has thus a O((n + m) 2 ) complexity. Note that this complexity is given for x encoding a permutation matrix. If more than n + m elements of x are different from 0, the associated complexity will be higher. We also avoid the storage of D since we can efficiently compute each term trough a function d (Alg.2, line 7) which computes D ij,kl according to Table 2. The overall memory complexity can thus be reduced to the storage of n + m mappings. Plugging Eq. 3 and Eq. 6 into Eq. 2, the cost of transforming G 1 into G 2 including both nodes and edge costs according to a mapping φ can then be written as:

S(V 1 , V 2 , φ) = Q e (V 1 , V 2 , φ) + L v (V 1 , V 2 , φ) = x T Dx + c T x (7)
with D and c = vec(C) being fully dependant on the two input graphs G 1 and G 2 and x encoding the mapping φ. Given a pair of graphs, Eq. 7 can be rewritten as only depending on x as S(x) = x T Dx + c T x. Given the Alg. 2, each edge is processed twice, i.e. once as (i, j) and once as (j, i). In case of undirected graphs, it will thus count each edit operation twice. We have x T Dx = 1 2 x T (D + D T )x and to handle both directed and undirected Algorithm 1 Computation of x T Dx(φ)

1: for all j, l | x jl = 0 do // O(n + m) if x ∈ Π 2: (x T D) jl = 0 3: for all i, k | x ik = 0 do // O(n + m) if x ∈ Π 4: δij ← (i, j) ∈ E1 5: δ kl ← (k, l) ∈ E2 6: (x T D) jl ← (x T D) jl + d(δij, δ kl , c(• → •) // see Table 2 7:
end for 8: end for 9: return x k+1 graphs, we introduce a matrix ∆ : [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF]. Hence, ∆ is symmetric in both cases. Including the linear part into the quadratic one leads to:

S(x) = 1 2 x T ∆x + c T x with ∆ = D if G 1 and G 2 are undirected D + D T else (8) Note that D is symmetric if both G 1 and G 2 are undirected
S(x) = x T ∆x, with ∆ = 1 2 ∆ + diag(c) (9)
Then, computing the GED of two graphs leads to minimize the quadratic function given in Eq. 9 under the constraint that x is a permutation matrix. This problem corresponds thus to the QAP:

GED(G 1 , G 2 ) = min x∈Π x T ∆x (10)
Note that quadratic formulation of GED have already been proposed: [START_REF] Justice | A binary linear programming formulation of the graph edit distance[END_REF] proposes a binary quadratic program transformed into a binary linear program but restricted to undirected and unlabeled edges. Slightly after, [START_REF] Neuhaus | A quadratic programming approach to the graph edit distance problem[END_REF] proposed another quadratic program but where optimization is only focused on node substitutions.

Resolution of the QAP

As stated in Introduction, computing an exact GED is intractable for most applications. Reducing the complexity relies thus on finding the best possible mapping, i.e. the one having the lowest mapping cost, in an acceptable computational time. Rather than finding an exact solution of a linear approximation of the problem, another strategy consists in optimizing the QAP associated to GED (Eq. 10). However, since ∆ is neither positive nor negative definite, finding a global minimum of Eq.10 is, conversely to LSAP, a NP-complete problem. Since this problem is NP-complete, most of algorithms find thus approximate solutions by relaxing the original problem to the set of doubly stochastic matrices.

As we resolve a QAP and not a generic quadratic problem, the solution should be a permutation matrix encoding a mapping. Since most of classic approaches used to resolve quadratic problems do not take this constraint into account, the final projection of the continuous solution to a permutation matrix may alter the approximation of the GED. Interger-Projected Fixed Point (IPFP) method [START_REF] Leordeanu | An integer projected fixed point method for graph matching and map inference[END_REF] has originally been proposed to resolve graph matching problems formalized as a maximization of a QAP. Conversely to classical methods used to resolve quadratic problems, IPFP allows to approximate the optimal solution by a gradient descent approach taking into account the nature of the solution. The algorithm mainly iterates over two steps: (i) compute a gradient direction which maximizes a linear approximation of the QAP in discrete domain and (ii) maximize the QAP in the continuous domain between the previous solution and the one found at step (i).

Algorithm 2 IPFP(x 0 , ∆, ε)

1: k = 0 2: repeat 3: b ← argmin b∈Π (x k T ∆) 4: α ← (x k T ∆)b -2S(x k ) + c T x k 5: β ← S(b ) + S(x k ) -(x k T ∆)b + c T x k 6:
t ← -α/2β 7:

x k+1 ← x k + t (b -x k ) 8: k ← k + 1 9: until ||x k -x k-1 || < ε 10: return argmax b∈Π x T k b
IPFP has been adapted from its original definition to GED computation [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF][START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF] by considering a minimization rather than a maximization and by considering the matrix ∆. The algorithm corresponding to this adaptation is shown in Alg. 3. Line 3 corresponds to the computation of the gradient direction. This step corresponds to a LSAP and can be efficiently computed using an Hungarian algorithm such as in bipartite GED framework. Then, a line search is performed to minimize the quadratic objective function in continuous domain. This optimization relies on the computation of a step size t (line 6) which can be computed analytically. Finally, the solution at convergence is projected onto the set of mapping matrices (line 10). Note that the stochastic matrix x k is projected by line 10 onto its closest permutation matrix. However, the distance betwenn x k and its projection can not be bounded a priori. In [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF][START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF], we reported the results obtained by using this approach and, as expected, it reaches a better accuracy than LSAP based methods. However, due to the non convexity of the objective function, the accuracy of the approximation is strongly dependant on the initialization x 0 .

Graduated NonConvexity and Concativity Procedure (GNCCP) [START_REF] Liu | GNCCP-Graduated NonConvexity and Concavity Procedure[END_REF][START_REF] Zaslavskiy | A path following algorithm for the graph matching problem[END_REF] is a path following algorithm which consists in approximating the solution of a QAP by using a convex-concave relaxation. This approach brings several advantages over IPFP. First, no initial mapping is required. This may avoid the accuracy variations induced by the initialization step observed with IPFP [START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF]. Second, this algorithm converges towards a mapping matrix [START_REF] Liu | GNCCP-Graduated NonConvexity and Concavity Procedure[END_REF]. This second property allows to avoid the projection step required at the end of IPFP algorithm which may alter the accuracy of the approximation. GNCCP algorithm is based on a weighted sum of a convex and a concave relaxation of Eq. 9:

S ζ (x) = (1 -|ζ|)S(x) + ζx T x (11)
For ζ = 1, S ζ (x) is equal to x T x which corresponds to a fully convex objective function. Conversely, for ζ = -1, S ζ (x) = -x T x which now defines a concave function. GNCCP algorithm starts with ζ = 1, and since the problem is convex, the initialization does not alter the result and will not influence the quality of the approximation as in [START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF]. Then, the algorithm, detailed in Alg. 3, smoothly Algorithm 3 GNCCP(x 0 , c, ∆, k max )

1: ζ = 1, d = 0.1, x = 0 2: while ζ > -1 & x / ∈ A do 3:
// Resolution of relaxed QAP according to ζ using Frank-Wolfe like algorithm 4:

x ← argmin x∈Π (1 -|ζ|)( 1 2 x T ∆x) + ζx T x 5:

ζ ← ζ -d 6: end while 7: return x k+1 interpolates convex and concave relaxations by passing from ζ = 1 to ζ = -1 with steps of size d until convergence is reached by having a mapping matrix or ζ = -1. Note that the minimum of the concave relaxation is a mapping matrix [START_REF] Zhou | Factorized graph matching[END_REF][START_REF] Zaslavskiy | A path following algorithm for the graph matching problem[END_REF]. Hence, conversely to IPFP, we do not need to perform a final projection step which may alter the approximation (end of Alg. 3).

For a given ζ, each iteration of GNCCP minimizes the quadratic functional defined in Eq. 11 using a Frank-Wolfe like algorithm. In this paper, we adapt IPFP algorithm used in [START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF] to perform this step and we remove last projection step (Alg. 3, line 10). However, since we optimize the relaxed objective function (Eq. 11), we have to update the linear subproblem and line search step. In our modified version of Alg. 3, the gradient direction b has now to minimize S ζ (x):

∂S ζ (x) ∂x = (1 -|ζ|) ∂S(x) ∂x + ζ2x T = (1 -|ζ|)(x T ∆) + ζ2x T (12)
Then, line 3 of Alg. 3 is updated with:

b ← argmin b∈Π [(1 -|ζ|)(x T ∆) + 2ζx T ]b (13)
Given the gradient direction b , line search step has also to be updated to minimize S ζ (x) rather than S(x). After some calculus, lines 4 and 5 are updated with:

α ζ ← [(1 -|ζ|)(x k T ∆) + 2ζx k T ]b -2S ζ (x k ) + (1 -|ζ|)c T x k (14)
β ζ ← S ζ (b ) + S ζ (x k ) -[(1 -|ζ|)(x T ∆) + 2ζx k T ]b + (1 -|ζ|)c T x k (15)
Note that some of the terms of Eq.

Experiments

Our method based on GNCCP is tested on 3 real world chemical datasets 4 . Similarly to our previous papers using these datasets [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF][START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF][START_REF] Bougleux | Graph edit distance as a quadratic program[END_REF], edge and node substitutions costs have been set to 1 and 3 for insertions and deletions. To avoid bias due to arbitrary order of nodes in data files, each graph adjacency matrix is randomly permuted. We compare our method to others methods computing GED. First, we used an exact GED computed using A but only available for small graphs (line 1). Second, we show results obtained by LSAP based methods using different information around each node: incident edges and adjacent nodes [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF] (line 2), bag of random walks up to 3 edges (line 3) [START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF], subgraphs up to a radius 3 (line 4) [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF]. We also show results obtained by IPFP approach, detailed in Section 3, with two different initializations: a random one (line 5) and one using [START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF] (line 6). Finally, results obtained by our method with d = 0.1 are displayed (line 7). First of all, Table 3 shows clearly that our method allows to reach the best accuracy. Considering the two first datasets Alkane and Acyclic for which we have a ground truth computed using A (line 1), average error of approximation is divided by over 2 for acyclic and by over 4 for alkane, hence showing a very good approximation of GED. These improvements allow to reduce the relative error to about 10%. Concerning PAH dataset, the quality of the approximation must be deduced from average distance since we can not compute the exact GED. Since the edit distance is always overestimated, the lowest edit distance may correspond to a better approximation. The quality of our approximation is also validated on PAH dataset since the average edit distance is reduced by about 20%. Despite the good approximation offered by GNCCP, we also observe an important increase of computationnal times. This phenomenon is not surprising since GNCCP iterates over IPFP. For d = 0.1, we may perform about 20 iterations in the worst case (Alg. 3). In practice, we observe an average of 11 iterations before reaching convergence. Nonetheless, computational times increase by a factor 40. This increase is due to the computation of quadratic term x T ∆x. Indeed, for ζ far from 0, the matrix x may contains a lot of entries different from 0. Therefore, the complexity stated from Alg. 2 no longer holds.

Conclusion

Considering the relationship between graph edit distance and nodes' mapping, we present in this paper a quadratic assignment problem formalizing the graph edit distance with substitutions, insertions and removals of nodes and edges. Following a previous optimization scheme based on a Frank Wolfe algorithm, we adapt a convex concave relaxation framework to the resolution of the QAP associated to edit distance. This approach does not require any initialization step and converges towards a mapping matrix. These two properties allow to compute a more reliable and robust approximation of the graph edit distance.

The experiments conducted on real world chemical datasets valid our hypothesis.

  14 and 15 have been already computed: [(1 -|ζ|)(x T ∆) + 2ζx k T ]b corresponds to the optimal cost computed by LSAP resolver on line 3 of Alg. 3, c T x k to the linear term included into ∆ and computed for the linear subproblem and S ζ (x k ) corresponds to the score of objective function computed at previous iteration. Finally, step size t is computed in the same way as in Alg. 3 with t = -α ζ /2β ζ (line 6).

Table 1 .

 1 Edit operations and costs encoded in matrix C

	i	j	edit operation	cost Ci,j

Table 3 .

 3 Accuracy and complexity scores. d is the average edit distance, e the average error and t the average computational time.

	Algorithm	d	Alkane e t	d	Acyclic e t	d	PAH t
	A *	15 -1.29	17 -6.02	-	-
	LSAP [11]	35 18	10 -3 35 18	10 -3 138	10 -3
	LSAP Random Walks [4] 33 18	10 -3 31 14	10 -2 120	10 -2
	LSAP K-graphs[3]	26 11 2.27	28 9 0.73	129 2.01
	IP F P Random init	22.6 7.1 0.007 23.4 6.1 0.006 63 0.04
	IP F P Init [4]	20.5 5 0.006 20.7 3.4 0.005 52.5 0.037
	GN CCP	16.7 1.2 0.46	18.8 1.5 0.33	41.8 6.24

Datasets are available at https://iapr-tc15.greyc.fr/links.html