Dynamic Spatial Verification for Large-Scale Object-Level Image Retrieval

Images from social media can reflect diverse viewpoints, heated arguments, and expressions of creativity, adding new complexity to retrieval tasks. Researchers working onContent-Based Image Retrieval (CBIR) have traditionally tuned their algorithms to match filtered results with user search intent. However, we are now bombarded with composite images of unknown origin, authenticity, and even meaning. With such uncertainty, users may not have an initial idea of what the results of a search query should look like. For instance, hidden people, spliced objects, and subtly altered scenes can be difficult for a user to detect initially in a meme image, but may contribute significantly to its composition. We propose a new approach for spatial verification that aims at modeling object-level regions dynamically clustering keypoints in a 2D Hough space, which are then used to accurately weight small contributing objects within the results, without the need for costly object detection steps. We call this method Objects in Scene to Objects in Scene (OS2OS) score, and it is optimized for fast matrix operations on CPUs. OS2OS performs comparably to state-of-the-art methods in classic CBIR problems, on the Oxford5K, Paris 6K, and Google-Landmarks datasets, without the need for bounding boxes. It also succeeds in emerging retrieval tasks such as image composite matching in the NIST MFC2018 dataset and meme-style composite imagery fromReddit.

Introduction
Image retrieval traditionally starts with the knowledge of what a user is looking for at query time and ends with a set of relevant results that match those initial expectations. This search regime must be transformed, however, in the age of social media. Nowadays, complex composite images like Internet memes are pervasive in the pools of content that are shared by users on social media [4]. Some of these i) ii) Figure 1. An example of a meme-style image from the subreddit /r/photoshopbattles. Internet memes are humorous messages that are spread on social media, often conforming to a set genre with a distinct style. In this paper, we provide spatial verification to find object-level correspondences between images, which assist in the retrieval of the donor images (highlighted in yellow) that contribute to composites like the one above.
images are of significant cultural value [37], while others represent odious extremist propaganda [10]. Both of these cases are important and timely topics of study, where the retrieval of the source content in reference image collections is paramount. The identification of near duplicate content shared between artistic images has received attention recently, and it can be solved using learned features that are both discriminative and invariant [35]. Further, the evolution of a meme can be traced by finding all of the related images, including images that contributed small donor objects ( Fig. 1), and referencing associated timestamps [24,5]. Similarly, one can debunk misleading or forged images being used for disinformation purposes by verifying identified source material [6]. In this paper we propose a new solution for spatial verification in image retrieval that supports these tasks by modeling object-level regions with the goal of retrieving the sources of small regions in large datasets. One of the major hurdles of image retrieval is the semantic gap [39,43], where a user's understanding of the content transcends the low-level representation and subsequent recall capabilities of a content-based search engine. Ideally, a good algorithm should return results that satisfy the user's higher-level intent, not just images with similar low-level features. However, retrieval approaches developed around this idea may not be sufficient for parsing composite imagery and retrieving relevant results, especially if there are aspects of an image that are not immediately apparent to the user at query time. In addition to the context of image creation and use, an image's extent is also shifted as it can be reused in composite images later on. And finally, given the free-form nature of composite imagery, an effective retrieval approach is expected to operate over a webscale database to maximize potential matching candidates.
With the above aspects in mind, we aim at adapting image retrieval to different contexts, such as composite images, especially those containing small spliced objects. For example, in Fig. 1, we see a composite image created from many smaller donor images. This type of image poses a significant challenge for existing image retrieval algorithms, because if the host (i.e., background) and donor images are matched globally, the latter group would not be highly scored since the content shared with the composite is very small. Nonetheless, tasks like meme analysis and disinformation debunking require the retrieval of each meaningful piece of content in an image under scrutiny. Because the image in question is a conglomeration from many sources, we can assume that the goal of a retrieval system should be to return instances of all images contributing to the composite.
To do so, we propose a new method of spatial verification that allows object-level instance scoring of retrieved results, without the need for costly object detection steps. We devise a feature-agnostic algorithm that utilizes a geometrically consistent voting measure inspired by Hough-based voting techniques, with the major difference of not requiring object regions of interest in the query to be known ahead of time, a manner to make the solution more appropriate to the current reality of unspecified retrieval context. As we show through experiments, the proposed method quickly and accurately localizes and ranks rigid objects contained within the query image to objects contained within a large image database. We call this method the Objects in Scene to Objects in Scene (OS2OS) score.
Through the rest of this paper, we will look at approaches that are related to the proposed OS2OS score, detail our methodology, and then perform experiments utilizing a number of relevant image retrieval datasets and features, including classic handcrafted and contemporary deep feature representations. In summary, the contributions are: • A new perspective on the problem of image retrieval, which aims at addressing the deficiencies in existing problem formulations for retrieval in cases of complex composites and other manipulated images.
• A new method called OS2OS score for spatial verification of matching objects between images, including very small objects that are important for understanding memes and other emerging Internet media. • A series of experiments showing the viability of the approach on the Oxford 5K [29], Paris 6K [31], Google-Landmarks [27], and NIST MFC2018 [25] datasets, as well as meme-style imagery from Reddit [24]. • A new experimental protocol for the Reddit Photoshop Battles dataset [33], preparing it to be used for benchmarking potential solutions to retrieve the donors of composite images.

Related Work
The typical CBIR solutions rely on multi-level representation of the images, to reduce the semantic gap between the pixel values and the system user's retrieval intent. In the lowest levels, typical methods use local features (a.k.a., keypoints) to obtain n-dimensional descriptions of the image content, ranging from handcrafted representations, such as SIFT [23] and SURF [3], to representations learned via neural networks, such as LIFT [22] and DELF [27].
In the subsequent levels, these local features are then used to index and compare, within the n-dimensional space they constitute (a.k.a., the feature space) and through euclidean distance or similar method, pairs of image localities (a.k.a., image patches). State-of-the-art large-scale indexing solutions comprise methods such as Optimized Product Quantization (OPQ) for approximate nearest neighbor (ANN) search [12], with Inverted File Indices (IVF) [20].
As proposed by Lowe [23], two features and their respective image patches are probably a match (i.e., they depict the same object or scene region in different configurations), if one is the nearest neighbor of the other within the feature space. Depending on the nature of the CBIR application, local features may be indexed in ways that use only the feature space and ignore or underutilize the scale, orientation, or (x, y) positions of the features within the images they belong to. That is the case for bag-of-features [38] and similar approaches [28,19], which are mostly useful for tasks such as retrieving semantically similar images. In the case of retrieving near-duplicates or finding spliced objects across a dataset, though, these techniques are not enough. Hence the need for additional spatial verification steps, such as the ones proposed in this work, to ensure that the local features being matched present a geometric coherence in either their scale, orientation, or (x, y) positions. Local-feature spatial verification methods can be organized into two categories, which are described below.
Hypothesis-oriented spatial verification. Methods in this category start with a set of spatial transformation hypotheses of one image towards the other (e.g., affine transformation). As originally proposed in the RANSAC [11] algorithm, these hypotheses are iteratively generated for random samples of feature matches, and are evaluated according to the overall number of matches that are inliers (i.e., matches that comply with the hypothesis). Aiming to make the process more accurate and deterministic, Philbin et al. introduced the Fast Spatial Measure (FSM) algorithm [29], which generates one hypothesis for every single feature match. Although very accurate, the major drawback of techniques from this category is the large runtime they demand, which is a quadratic function of the number of features (as we show through experiments in Fig. 4).
Hough transform-based spatial verification. Methods from this category start with Hough transforms [8,2] and the computation of histograms for their parameters, where each bin quantifies the number of agreeing feature matches. Lowe [23] proposed the adoption of four-parameter Hough transforms, computing bins with respect to the product of (i) x and (ii) y coordinates, (iii) scale, and (iv) orientation of the features. The largest histogram bin is then chosen to select a better and potentially reduced set of feature matches, prior to applying RANSAC on them. Aiming to speed up up the overall process, Avrithis and Tolias later introduced the Hough Pyramid Matching (HPM) strategy [1], which employs a hierarchical voting strategy to recursively split pairs of feature matches into bins in a top-down fashion (from coarser to finer correspondences), as a way to evaluate the pairwise affinities of matches without enumerating all pairs. Similar to HPM, Li et al. also suggested the evaluation of pairs of feature matches to develop their Pairwise Geometric Matching (PGM) strategy [21]. However, they proposed to use coarser (and therefore faster to verify) two-parameter Hough transforms, relying only on the orientation and scale change between the compared images. This indeed guarantees a significant speed-up in the spatial verification process, as we show through experiments reported in Fig. 4.
Another approach worth mentioning is the one proposed by Shen et al., namely Spatially Constrained Similarity Measure (SCSM) [36]. Although this one also makes use of four-parameter Hough transform quantization to enumerate candidate spatial transformations, its difference relies on the method to select the best transformation. By demanding the establishment of a bounding box over the features of the query object before performing the retrieval task, the algorithm uses the box center to measure the quality of the candidate transformations. Thus, the best transformation is the one that, after its application, best preserves the spatial relations among the feature positions and the box center.
Lastly, Schönberger et al. introduced the Vote and Verify (VaV) [34] method, where HPM is used to compute votes across the Hough transform hierarchy, prior to verifying transformation hypothesis for only the most voted bins.
Putting the proposed method in context with prior work. Our method belongs to the latter category of spatial verification techniques and is agnostic to the chosen local features and feature indexing approach. When compared to the literature, the novelty of this work comprises, besides the unique combination of strategies for the task at hand: • Inspired by SCSM [36], the use of match-set-wise centroids to evaluate the quality of the feature matches, without the need for selecting bounding boxes of interest in the query ahead of retrieval time. These centroids are calculated in a novel way (see Eq. 2). • Inspired by PGM [21], the use of a coarse, fast, but still accurate two-parameter Hough transform quantization, which contrary to PGM, relies on the x and y position coordinates of the matched features, instead of scale and orientation. • A novel image retrieval score (OS2OS score, see Eq. 11), which is based on two complementary measures of object-level image matching quality (namely, object centrality, and angle coherence), and allows for spatial verification while ranking images.
In the following, we detail each step of the proposed solution, as well as provide discussion on the reasons and advantages of adopting each one of the above novel aspects.

Objects in Scene to Objects in Scene (OS2OS) Score
Our proposed spatial verification method for image retrieval can be explained in five steps.
Step 1: Local feature affinity. Images are described through local feature vectors and their respective geometric data, namely (x, y) location, scale, and orientation angle. Let Q be the set of all features extracted from the query, as depicted in Fig. 2 (i), and D be the set of all features extracted from a target image database, as depicted in Fig. 2 (ii). For a particular query feature q i ∈ Q, with i = {1 . . . |Q|}, we compute the K nearest neighbors d j ∈ D in the feature space only, ignoring the images they come from. As a result, q i participates in K matches A score S for a given match m ij that can express the affinity between q i and d j is calculated via the L 2 -distance and the rank-adaptive scoring outlined in [18]: where φ is a fixed rank position of reference (usually K/2).
Step 2: Image-pairwise centroid calculation. Take an image P from the database that shares content with the query. There might be a set of feature matches between them, whose incident locations (x, y) onto the query should give a rough indication of the shared object's location within Q. Therefore, if we calculate the center of these P -wise match locations on the query, we find a point that is generally near a potential object of interest, serving as an estimation of its centroid c (see Fig. 2 (iv)). Let M be the set of feature matches m k shared between P and the query, with k = {1 . . . |M |}, and let Q m ⊆ Q be the set that contains only the query features that have a match to P . To obtain c, we apply the straightforward solution of using the Euclidean center of the query features q k ∈ Q m . Nevertheless, aiming to consider the quality of the matches while computing c, and to deal with spurious matches, we also weight the added features according to their affinity scores S(m k ) (see Eq. 1) associated with the features of P : where L(q k ) is the (x, y) location of the k-th feature q k ∈ Q m , and m k ∈ M is the respective k-th feature match between the query and P . The motivation for this is that the more similar two matched features are (i.e., the higher their value of S(.)), the more they contribute to the position of c. This strategy reduces the Hough voting noise problem described in [34], as shown via supplemental ablation experiments (see Supp. Mat.), where we report the decrease in performance due to the absence of centroid computation.
Step 3: Centroid-relative feature projection. Given the centroid c representing the query feature locations q k ∈ Q m , and their respective matched features p k ∈ P , we can estimate the translation, rotation and scaling transformation from the space of image P towards the query space, for each match m k = (q k , p k ), with k = {1 . . . |M |}: where θ(.) and σ(.) respectively provide the angle and the scale associated with the location L(.) of either q k or p k features. These transformations describe where each query feature expects the object region to be, similar to an "Rtable" of the general Hough transform [2], with the novelty that we only consider the (x, y) feature coordinates. The advantage of doing so is twofold: (i) the Hough space is twodimensional (given x and y) instead of four-dimensional, making computations faster; (ii) the Hough space maps directly to the query space, making the localization of shared objects straightforward. By applying each transformation to the (x, y) location of its respective matched feature L(p k ), we subsequently build a voting space V k for each p k ∈ P : The process of computing the voting space is shown in Fig. 3. Observe, through item (iii), that all p k features contribute to the centroid, except for the spurious p 5 one.
Step 4: Density-based feature clustering. We calculate a distinct centroid c for every image P from the database that presents matches with the query. This allows us to transform all matched feature vectors, regardless of the image they come from. After all matched feature locations have been transformed to their respective two-dimensional Hough vote space, as depicted in Fig. 2 (vi), a densitybased clustering determines which sets of feature matches are structurally consistent with the query. Instead of using a computationally intensive algorithm such as DBSCAN [9] for clustering millions of points, we apply a linear-time Figure 3. OS2OS voting space computation. (i) Translation vectors are calculated relative to a computed centroid c in Q. (ii) Translations are applied to matching features in P . (iii) Vectors are rotated by the angle differences a k to account for object rotation. Notice that the p5 spurious match votes incorrectly and therefore will not contribute to the score. (iv) Votes are binned using (x, y) coordinates only. (v) A density map is calculated, which conveniently coincides with the query space, providing visualization of matched object regions.
grid-based quantization to each V k , to find the approximate clusters of highest vote density. Because we assume, for simplicity sake, cluster morphologies to be roughly circular, we employ a square sliding window to quantize and bin V k values. The size ws of the quantization window varies according to the resolution w × h of image P : where b is a scaling factor, and 0 prevents ws from becoming too large. The magnitudes of projection vectors ||T k || are proportional to the resolution of image P . As ||T k || increases, small errors in the scale normalization and rotation transforms (Eq. 4) are amplified, resulting in lower density clustering. Unlike HPM [1], which accounts for this phenomenon by utilizing a hierarchy of window sizes for assigning votes to clusters, we employ a single window size determined by a sub-linear mapping of the image's resolution (Eq. 6). This is accomplished through : it constrains the maximum allowable vote cluster area to grow sub-linearly with the resolution of the image. The values of b and are empirically determined via an ablation study, included in the Supp. Mat.
Step 5: OS2OS filtering and scoring. By relying on the value of ws computed for an image P , the respective votes V k are quantized into cluster bins (see Fig. 3 (iv)): We consider that all V k values sharing the same bin(.) value belong to the same matched object. Let O be the set of transformed features from database image P , which participate in matches between P and the query, and whose respective votes happened to belong to the same bin according to Eq. 6. The meaning of this is that the O features likely belong to a unique object shared by the images. Subsequently, we enforce a strict one-to-one matching constraint within each O. Then, to express their affinity of features within the filtered O, we propose two novel main scoring mechanisms inspired by [41], each with a particular purpose. The first, called the Centrality Score (CS), measures the centrality of the features and is calculated as the sum of location differences between the elements o k ∈ O and their average (a.k.a., central) elementō, with k = {1 . . . |O|}: where L(.) is, again, the (x, y) location of the given feature, and pdf (.) is a function that augments the available data through a probability density function (see Fig. 3 (v)): In addition, CS is normalized by the cardinality of O, as a way to balance clusters of small objects containing few features with respect to large objects that contain many. The second mechanism is the Angle coherence Score (AS), which aims at measuring the uniformity of angles within the features of O. Consider the feature-wise difference of angles a k expressed in Eq. 4, and let A be the set of difference of angles a k computed for each feature of O. Features from keypoints belonging to a single rigid object are expected to present similar values of a k , while erroneously aggregated unrelated features are expected to present more diverse results. For that reason, we rely on the inverse of the standard deviation of A, stdv(A), to compute the angle uniformity (shifted to avoid division by zero): Finally, the OS2OS score is given by CS and AS: which lays between 0 and 1 and rewards only clusters in which each point has been transformed similarly. The logarithmic scalar log |O| penalizes clusters with very low numbers of votes, pushing their score towards zero.

OS2OS Score Evaluation
Datasets. Oxford 5k and Paris 6k. Oxford 5k [29] and Paris 6k [30] are smaller popular datasets for image retrieval performance evaluation. These two datasets contain hundreds of true positive matches per query, rather than similar sized datasets that contain only four or five [26,17].
Google-Landmarks. As a benchmark for DELF features [27], Google released the 2018 Google-Landmarks dataset [14] and subsequent evaluation protocol [13]. This dataset contains 1, 098, 461 images with 117, 703 queries in the testing protocol. While groundtruth data for the test set has not yet been released, 1, 212, 281 weak labels for the training set are available. The training set contains a total of 14, 951 unique landmarks, with each landmark having an average of 80 instances.
NIST Media Forensics Challenge 2018 (MFC2018). As part of the yearly Media Forensics Challenge run by NIST, the MFC2018 dataset [25] was constructed specifically for the task of finding related manipulated images in a forensics context. This is a large dataset, 3.1TB in size, containing 1, 031, 080 images and 3, 300 queries. Many of these queries are composites, with groundtruth results provided in the MFC2018 testing protocol. The dataset also provides groundtruth as to whether a database image contributes a majority of its content (e.g., Fig. 1 i), or only a particular object (e.g., Fig. 1 ii) to the query.
Reddit. The Reddit dataset [7] introduced by Moreira et al. [24] is collected from the Reddit Photoshop Battles [33] subreddit. Each case provides the original image and all subsequent manipulated versions of the original. This dataset contains 51, 245 images from 185 different Photoshop battles. To utilize this dataset for the OS2OS task, we generate a query set of one image chosen at random from each of the 185 cases. The rest of the images from the same battle thread are considered as the groundtruth for relevant images to these queries. We will make this new groundtruth available upon the publication of this paper.
Features. For each dataset, we report performance using both handcrafted SURF [3] and learned DELF [27] features. While other deep image representations have been proposed for image retrieval [16,42,15,32], some deep local feature  descriptors such as LIFT [22] are too slow (see Table 1) for practical use, as also observed in [27]. Other fast global descriptors such as MobileNet [16] are not applicable to the localization of multiple objects. Our region-wise matching approach requires local descriptors to match coherently within a certain spatial location, which cannot happen with a single global MobileNet descriptor.
SURF keypoints are detected in a distributed modality (DSURF), as proposed in [24]. Keypoint extraction of DSURF features automatically provides the location, scale, and rotation data needed for computing the OS2OS score. For all datasets, we extract a maximum of 5, 000 64dimensional DSURF features per image, along with their respective geometrical data. DELF features are used in similar way to DSURF. Because the DELF algorithm provides only feature scale information, we use the SURF keypoint angle algorithm [3] as an extension to provide feature angles for the DELF geometric data. We performed experiments using the default parameters and model to produce 40-dimensional local features, and cap the attention model at a maximum of 1, 000 features per image.
Indexing. For experiments on Oxford 5k and Paris 6k datasets, we keep a consistent indexing backbone among spatial verification methods. To index features, we use Optimized Product Quantization (OPQ) for Approximate Nearest Neighbor (ANN) search, with Inverted File Indices and Asymmetric Distance Computation (IVFADC) [20]. Following the vocabulary hold-out protocol suggested in [40], we utilize a randomized 1% subset of Google Landmark images to train the quantization table used for Oxford 5k, Paris 6k, and Reddit datasets. Due to the large volume of images available, a randomized 5% hold-out of the Google Landmarks dataset is used to train its own tables. These tables are obtained via OPQ matrix computation [12], and used for the IVFADC centroids. OPQ training runs for 25 epochs for 2 18

Experimental Results
Here we describe and analyze the results for a series of timing experiments, as well as the experiments for each of the five datasets described in Sec. 4.

Timing and Complexity
Feature extraction timings. Aiming to focus on largescale retrieval, we performed average timing experiments over a subset of 100 images from the Reddit dataset to analyze the tractability of different feature extraction methods on CPUs. Table 1 shows these results, justifying our selection of DSURF and DELF as candidate features for further experiments. MobileNet stands for the MobileNet-v2 architecture [16], whose features were extracted from its global pool layer. Although fast, local spatial information is not native to MobileNet, and is therefore incompatible with the task of spatial verification.
Spatial Verification Timings. We also performed timing experiments against other spatial verification and ranking methods for image retrieval. Here we compared the OS2OS score against HPM [1], PGM [21], FSM [29], SCSM [36], plain RANSAC [11], and VaV [34]. We extracted 5, 000 SURF features from a query image and varied the K retrieved nearest neighbors for each query feature from 0 to 400. For these experiments we averaged timings across 10 query images from the Google-Landmarks dataset. Standard deviations were calculated, but were too small to be plotted. While RANSAC is known to provide spatial verification for an arbitrary number of features in quadratic time [21], HPM and PGM are proven to be linear [1,21]. As can be seen in Fig. 4, OS2OS scoring is faster than HPM, PGM, and VaV, ranking nearly 400 images with 1.8 million features in only 10 seconds. All methods used Precision-recall curves for the Google-Landmarks dataset. As it can be observed, the usage of spatial verification through the OS2OS score boosts both SURF-and DELF-featurebased image instance retrieval. the same 2.7GHz single-core environment.

Image Retrieval
Oxford 5k and Paris 6k. The small-scale experiments performed on the Oxford 5k and Paris 6k datasets are meant to show that the OS2OS scoring algorithm, while designed for object-level spatial verification, is general enough to provide benefits for typical image retrieval. We compare our approach against HPM [1], PGM [21], VaV [34], and plain usage of SURF and DELF features (without spatial verification). Mean Average Precision (MAP) scores are reported in Table 2. We find that the OS2OS score significantly improves both DSURF and DELF features, suggesting that the provided spatial constraints work satisfactory for global geometric verification for instance retrieval. Additionally, we see a much larger performance improvement for DSURF, which suggests the OS2OS score helps mitigate erroneous matches from the bursty nature of SURF. Overall, the OS2OS approach is comparable to other approaches in the literature. Of noteworthy significance is the fact that the OS2OS algorithm required no bounding boxes, still performing comparably to other spatial methods. We additionally examined scores from RANSAC [11], SCSM [36], and FSM [29] for Oxford 5k, but the results were inferior to the better performing approaches in Table 2.
Google-Landmarks. This experiment shows that the OS2OS score generalizes to the instance retrieval task, as well as showcases the algorithm's scalability in the presence of many distractors. Because annotations for the index or test sets have not been released at the time of this writing, we performed our study utilizing the training set, which contains over 1 million images and provides landmark annotations. We selected 1, 000 landmarks at random for retrieval, sampling one query per landmark. Finally, we utilized the modified precision and recall measures described  in [27] to report results. Fig. 5 shows that DSURF augmented with the OS2OS score improves significantly, while also providing minor improvements to DELF.
NIST MFC2018. Unlike the experiments described thus far, the MFC2018 dataset is comprised of manipulated images. Query images from the dataset's retrieval protocol may or may not contain regions from multiple sources within the image database. We utilized the groundtruth relationship graphs to determine which images donate small objects to their respective queries. Using this data we can generate recall curves exclusively for donor image retrieval (namely, donor recall). We report both total recall scores and donor recall scores in Fig. 6. While the OS2OS score shows good boosts in total recall, we see that donor recall is more significantly improved, indicating that the OS2OS score is capable of balancing geometrically coherent matches of image regions from small donors with global matches from backgrounds.
Reddit. The Reddit provenance dataset has proven to be a difficult challenge [24]. In Table 3, we see a significant increase in retrieval performance (nearly 10% for SURF and nearly 5% for DELF) across the board, with vastly superior results when compared to VaV [34]. The near-baseline performance of VaV suggests that global spatial verification methods are not fully adequate to solve the OS2OS problem. Further, Fig. 7 shows qualitative retrieval results using OS2OS scoring, along with the object vote maps from each retrieved match. These results highlight the OS2OS score's ability to localize and appropriately weight small objects from a database image to small objects within the query, without the need for bounding boxes.

Parameter Ablation
To examine the contribution that different parameters from Sec. 3 provide for both total and donor recall metrics (see Fig. 6), we perform an ablation study using the MFC2018 dataset. We remove and re-introduce the CS score (Eq. 8), AS score (Eq. 10), logarithmic scaling (Eq. 11), and centroid calculation (Eq. 2), to empirically show that each component plays an important role in increasing the performance of the algorithm. As a consequence, the best configuration is indeed the one that uses all of these features. We also vary the bin size ws in Eq. 6, to illustrate its effect on filtering and scoring performance. These additional results are included in the Supp. Mat., due to space constraints.

Conclusions
Retrieving images that share small regions in a complex composite scene is a challenge for most image retrieval approaches. In this paper, we proposed an inexpensive scoring technique that is based on the better utilization of pretrained feature extractors and indexing techniques to yield geometrically consistent localized scores. An advantage of the technique is that it is learning-free, and works as an addon to existing feature description methods. It provides a way to include spatial verification while performing matching in a large feature space. It is optimized to be efficiently executed on CPUs without the need for high-end GPUs.
In experiments, the proposed approach improved recall for the retrieval of images for difficult emerging problems such as cultural analytics and image forensics. As with most computer vision algorithms, the proposed approach still struggles on challenging datasets obtained from the web with different styles of content correspondence. Although performance improves with the proposed scoring, retrieval for specific applications such as tracking memes on social media is worthy of additional research consideration -especially as such content grows in popularity.


BayesFit: A Tool for Modeling Psychophysical Data Using Bayesian Inference

BayesFit is a module for Python that allows users to fit models to psychophysical data using Bayesian inference. The module aims to make it easier to develop probabilistic models for psychophysical data in Python by providing users with a simple API that streamlines the process of defining psychophysical models, obtaining fits, extracting outputs, and visualizing fitted models. Our software implementation uses numerical integration as the primary tool to fit models, which avoids the complications that arise in using Markov Chain Monte Carlo (MCMC) methods [1]. The source code for BayesFit is available at https:// github.com/slugocm/bayesfit and API documentation at http://www.slugocm.ca/bayesfit/. This module is extensible, and many of the functions primarily rely on Numpy [2] and therefore can be reused as newer versions of Python are developed to ensure researchers always have a tool available to ease the process of fitting models to psychophysical data.

Introduction
Fitting statistical models to behavioural data is an important part of analyzing task performance within the field of Psychology [3]. Within the domain of psychophysics, a fundamental statistical model is the psychometric function, which relates performance -usually response accuracy or latency -to the stimulus, and is used to quantify stimulus detection or discrimination [3]. In experiments that measure response accuracy in a stimulus detection task, the psychometric function indicates how the probability of a correct detection response is related to stimulus intensity. Typically, correct responses are made more frequently as the stimulus intensity increases [3,4].
Although fitting a psychometric function to data is a common component of the analysis of behavioural data, very few tools have been designed to help researchers streamline this fitting process [1,3,4,5], and only one of these has been adapted for use in Python [1]. Furthermore, the only tool available for Python does not support versions 3.x or greater. We therefore created a module supporting newer versions of Python (3.x) that implements a Bayesian algorithm to fit psychometric functions to behavioural data with as few steps as possible. The resulting software, BayesFit, accomplishes this task by taking advantage of a simple API that provides a convenient method for users to build probabilistic models using numerical integration in Python. Our module also makes it easy for users to fit models to large numbers of datasets while constraining the parameters used across model fits.

Usage Example 1: Fitting Models to Psychophysical Data using MLE
The full API documentation for BayesFit can be found at: http://www.slugocm.ca/bayesfit/.
To demonstrate the functionality of this software, below we provide the results from a simulated twointerval, forced-choice (2-IFC) detection experiment in which each experimental trial consists of two successive, temporal intervals, and the observer must determine whether the target stimulus was presented in the first or second interval. Stimulus intensity is varied across trials, and response accuracy varies from near chance levels at low intensity to nearly perfect accuracy at high stimulus intensity (see Figure 1). Detection threshold is estimated by fitting a curve to the response accuracy data and then determining the stimulus intensity that produces an intermediate level of response accuracy (e.g., 75% correct).
BayesFit expects that data provided by the user are organized into a m-row by 3-column Numpy array, where the 1st column corresponds to the stimulus intensities used, the 2nd column the number of correct responses made by an observer at each stimulus intensity, and the 3rd column the number of trials run at each intensity level. Data in columns 2 and 3 should be organized such that values along each row correspond to the level of intensity specified along the same row in column 1 (see Figure 2).
Once data are appropriately organized, we are ready to use BayesFit to perform the model fitting procedure. We first import the BayesFit module into our workspace via: import bayesfit as bf Next, we specify the output variables and options that we want to use for our fitting procedure. BayesFit uses only one main function called fitmodel to streamline the process of fitting models to psychophysical data. If no priors are provided to the function, the fitting procedure is equivalent to performing maximum likelihood estimation (MLE). A basic example in using the fitmodel function to fit a model using MLE would be as follows: Although we have listed all options that are available to the user in fitting models using BayesFit, the only mandatory argument to provide the fitmodel function with is data. However, if additional options are not specified, default values for these options will be assigned. Complete API documentation describing each of these options, as well as their default assignments, are located at: http://www.slugocm.ca/bayesfit/. Following the fitting procedure, the fitmodel function outputs two variables that include: 1. Metrics -a dictionary containing parameter estimates from the fitted model (s), metrics about the goodnessof-fit of each parameter, arrays corresponding to the marginal distributions, the likelihood surface (posterior if priors provided), as well as a numerical approximation of the threshold from the fitted function (at level specified in options). 2. Options -the options used to fit model (s) to data.
Using these outputted variables, we can generate a plot of the fitted model using the plot_psyfcn function of the BayesFit module. For example, to visualize a psychometric function (Cumulative Normal sigmoid) fit using the BayesFit module to the data shown in Figure 1, we would specify the following snippet of code: bf.plot_psyfcn (data, options, metrics) As can be seen in the example above, only a few lines of code are needed to generate parameter estimates for our model and quickly visualize the results (see Figure 3). Although our example above was performed using only a single dataset, the BayesFit module also includes a method for fitting models to batch datasets so that the same model definitions can be used across all models fit. Once again, we encourage the reader to visit the API    , we can pass these arguments along to the fitmodel function of BayesFit during our fitting procedure. The fundamental difference in passing priors to BayesFit is that we are now using Bayesian inference methods rather than MLE in estimating parameters of the psychometric function, and thus are computing the posterior distribution. BayesFit provides users with 5 different distributions that can be used to specify priors for each parameter. These include: Therefore, if the user knew that certain values were more likely compared to others for the parameters controlling the scale and slope of the psychometric function, then the user could pass those arguments along to the fitmodel function of BayesFit as follows: From this example, it can also be seen that choosing prior distributions, and the values that define them, can be difficult. Therefore, it is best to use priors only when you are confident that the probability of certain values for parameters to-be-estimated are more likely than others based on theoretical or empirical grounds. Following the fitting procedure, the user can now visualize the priors used during the fitting procedure using the plot_priors function from BayesFit (see

Usage Example 3: Batch Fitting Models to Datasets
Occasions might arise where the user would like to automatically, rather than manually, fit multiple psychometric models across several datasets. The BayesFit module makes this task easy by requiring only two small changes to the procedure used to fit a single model. The first change is that the batch input argument for the fitmodel function must be set to True. The second change is that each dataset must be stored in a single dictionary object, with each dataset occupying a separate key (see Figure 7).
Here is a simple example of using the fitmodel function to perform batch fitting on multiple datasets:

# Store each dataset in a separate key under a # dictionary object
# Use fitmodel to fit functions across multiple # datasets metrics, options = bf.fitmodel (data = data, batch = True ) Note that only user-specified values for fixed parameters are used when performing batch fitting. Furthermore, to  prevent memory issues, only select metrics will be saved for each psychometric model fit.

Performance evaluation using simulated data
In order to assess the accuracy of parameter estimates using BayesFit, we generated 100 synthetic datasets from simulated observers in a psychophysical experiment. The datasets were generated to mimic the response of observers performing a 2-interval forced choice task with with eight equally-spaced stimulus intensities (range: 0.1-0.9). The true probability for an observer of correctly responding to a given level of intensity, x, was modelled as: where x is the intensity of the stimulus, α and β determine the scale and shape of the Weibull sigmoid function, γ is the guess rate, and λ the lapse rate. With the exception of the parameter γ which was fixed at 0.5 (i.e., 1/nafc), the parameters that define an observers› sensitivity at different stimulus intensities were drawn from the distributions: (2) These parameters were used to generate 100 binomial random responses at each stimulus intensity, for a total of 800 responses. In other words, 100 responses were simulated for each observer at each stimulus intensity. Data from each simulated observer were fit using the default settings for BayesFit, with the exception that a Weibull function was specified for the sigmoid type.
Results regarding the accuracy of parameter estimates for α, β, and λ are shown in Figures 8 and 9. Overall, Figure 7: Example of how multiple datasets should be combined into a single dictionary object before being passed as an argument to BayesFit for batch fitting. parameter estimates using BayesFit with default settings are extremely good. Estimates for α, a parameter that is typically used to represent the threshold of an observer, is very accurate. Estimates for parameters β and λ are also quite close to the true parameters used to generate data for a given observer, but are more variable in their accuracy to true values compared to parameter α. These simulations demonstrate that BayesFit can be trusted for use in estimating the parameters of models used to fit psychophsyical data. It should be noted that, when possible, greater care should be taken in choosing the options used to fit models to psychophysical data. However, these simulations serve as a worst case scenario where the user fails to specify prior knowledge that may help the fitting procedure, and yet, BayesFit still performs extremely well.

Comparison in accuracy of estimating thresholds to Psignifit 4.0
To evaluate the performance of BayesFit in comparison to Psignifit 4.0, the only other module the authors are aware of for fitting models to psychophysical data in Python (albeit 2.x), we make use of the simulated data used above for judging the accuracy of BayesFit in parameter estimation. Because BayesFit uses definitions for psychometric functions according to those specified in [6], and Psignifit 4.0 uses parameterizations of the psychometric functions according to [7], a direct comparison between parameter estimates obtained using Psignifit 4.0 and BayesFit cannot be readily made.
However, we can assess the accuracy of threshold estimates obtained from the model fit to a simulated observer's dataset using each module. Therefore, we used Psignifit 4.0 to fit models to simulated data using default settings except for specifying the sigmoid type as a Weibull function, and also defining threshold at 75% proportion correct response. Figure 10 show boxplots of the distributions of 75% threshold estimates obtained using BayesFit and Psignifit 4.0 compared to the distribution of true threshold values. Results show that threshold estimates obtained using BayesFit are remarkably similar to those obtained using Psignifit 4.0, and both distributions accurately resemble the target distribution for threshold defined at 75% correct response. These results suggest that both BayesFit and Psignifit 4.0 perform similarly in estimating thresholds for data obtained from a simulated psychophysical experiment.

Implementation and architecture
BayesFit was written in Python [8], and makes extensive use of the Numpy module [2] in constructing functions for BayesFit. The core function of BayesFit called fitmodel computes posterior estimates using numerical integration, and serves as the main function used in fitting models. Besides this core function, BayesFit uses several utility functions to compute parameter estimates. Although most of these functions require only Numpy, probability distributions currently are defined using a subset of functions from the module SciPy [9]. BayesFit also contains several plotting functions dedicated to visualizing results and analyses which make extensive use of the module Matplotlib [10]. Figure 11 shows a schematic of the architecture for BayesFit, along with displaying the module dependencies required by a given function.

Quality control
Both unit and functional testing has been performed on BayesFit. Travis CI (https://travis-ci.org/) is used to perform regular unit tests on BayesFit to meet strict quality control standards. A coverage report can also be found on the GitHub repository for BayesFit via the coverage badge.

Issue tracker
Although BayesFit undergoes extensive testing to ensure the software is functioning properly, if any issues arise during use, please let us know under the issue tracker for BayesFit on Github: https://github.com/SlugocM/bayesfit/issues.

Operating system
BayesFit should function on any software capable of running Python 3.5 or greater. Current operating systems that BayesFit has been tested on include macOS Sierra, Windows 10, and Ubuntu 16.04.

Programming language
Python (version 3.5 or greater.)

Additional system requirements
The minimum system requirements needed to run Python 3.5 or greater.

(3) Reuse potential
The need to estimate parameters of models fit to psychophysical data exists in many different domains within Psychology. Therefore, the potential for reuse of this software to be high. The core code of BayesFit has been written using only the Numpy module, and therefore many functions in BayesFit can be used in other projects that want to use numerical integration and Bayesian inference in modelling psychophysical data. Some highlighted functionality includes: 1. Checking user arguments for fitting psychometric functions. Figure 11: Schematic of the architecture for BayesFit, also displaying module dependencies for each function, whether directly or using a function that also depended upon use of a specific module.


A Comparison on Fine-grained Pre-trained Embeddings for the WMT19Chinese-English News Translation Task

This paper describes our submission to the WMT 2019 Chinese-English (zh-en) news translation shared task. Our systems are based on RNN architectures with pre-trained embeddings which utilize character and sub-character information. We compare models with these different granularity levels using different evaluating metics. We find that a finer granularity embeddings can help the model according to character level evaluation and that the pre-trained embeddings can also be beneficial for model performance marginally when the training data is limited.

Introduction
Neural Machine Translation (NMT) systems are mostly based on an encoder-decoder architecture with attention. Given a sentence x in source language, the model predicts a corresponding output sentence y in target language, which maximizes the conditional probability p(y|x). The attentionbased Recurrent Neural Network (RNN) version of this architecture has been a very popular approach to NMT (Bahdanau et al., 2015;Luong et al., 2015). Despite the success of these models, they still suffer from problems such as outof-vocabulary (OOV) words, i.e. words that have not been seen at training. To alleviate the OOV problem, we follow the methods used in word representation and segment words into smaller units. In some morphorlogically rich languages such as Chinese, a word can be divided into characters and then the characters can be further divided into smaller components called glyphs. Both character and glyph might contain semantic information and therefore utilizing such information might help alleviate the OOV problem.
Based on the RNN attention-based model (Bahdanau et al., 2015), we experiment with different granularity levels on the WMT19 Chinese-English (zh-en) news translation shared task. This paper describes our submitted systems with embeddings pre-trained on monolingual corpora. The two submitted systems use pre-trained embeddings enhanced by character and sub-character information respectively. The preprocessing methods include Chinese word segmentation, tokenization, data filtering based on rules and Byte Pair Encoding (BPE). Our baseline model is based on RNNSearch (Bahdanau et al., 2015) operating on word level and we use Long Short-Term Memory (LSTM) (Hochreiter and Schmidhuber, 1997) as encoder and decoder. For character level word embeddings, we use the Character-Enhanced Word Embedding (CWE) proposed by Chen et al. (2015). For the sub-character level embeddings, we use the Joint Learning Word Embedding (JWE) proposed by Yu et al. (2017). We use various metrics, namely BLEU (Papineni et al., 2002), METEOR (Denkowski and Lavie, 2011), TER (Snover et al., 2006) and CharacTER (Wang et al., 2016) for evaluation.
When compared with our baseline model, the models with pre-trained sub-character level embeddings on monolingual corpus show better performance, achieving an increase of +0.53 BLEU score with the sub-character level embeddings. We ran additional experiments on the character and subcharacter level pre-trained embeddings and found that the use of these embeddings can benefit the model when the training corpus size is limited.
This paper is structured as follows: Section 2 introduces the related work including the model architecture and pre-trained embeddings used in our experiment. In Section 3, data selection and preprocessing methods are described. Section 4 introduces the model architectures and hyperparameter settings. Section 5 shows the evaluation results on models with different granularity levels. Section 6 shows additional experiments to better understand our models.

Related Work
NMT has been an important task in Natural Language Processing.
A translation system aims to find the corresponding target sentence y = {y 1 , y 2 , ..., y m } given a sentence x = {x 1 , x 2 , ..., x n } in source language, in a probabilistic manner, represented as max y P (y|x). Most NMT models are based on the sequence-tosequence approach, and the RNN-based architecture  with attention (Bahdanau et al., 2015) is a popular version of such an approach. The attention mechanism functions as a dynamic calculation of the context vector. At each decoding step, a probability distribution is calculated based on the current decoder hidden state and all encoder hidden states. This distribution is defined as the attention score, representing the importance of each input token at current decoding time step. The context vector is calculated as a weighted average of all encoder hidden state vector, where the attention score is the weight. With the introduction of attention, the model does not need to rely on a single context vector to represent the whole sentence and thus can better handle long sentences.
In recent years model architectures based on convolutional neural networks (Gehring et al., 2017) and transformers (Vaswani et al., 2017) have shown competitive or better performance than RNN-based architectures. In addition, strategies such as back translation (Sennrich et al., 2016a), reranking (Neubig et al., 2015) and model ensembling have led to improvements in translation quality. In our experiments, we only experiment with RNN architectures and focus on the effect of using character and sub-character level embeddings and only use ensembling for comparison purposes.
We use the CWE model proposed by Chen et al. (2015) and the JWE model proposed by Yu et al. (2017) for pre-trained embeddings training. Both models are based on the word2vec proposed by Mikolov et al. (2013). Based on Continuous-Bagof-Word (CBOW), the CWE model construct a new word representation by summing the word embeddings with character embeddings (see Eq 1). Chen et al. also proposed a multi-prototype character embeddings where characters are tagged with additional factors, such as position and con-text cluster, for character disambiguation.
where w j is the word embeddings and c k is the embeddings of the k-th character in x j . ⊕ is the composition operator (either addition or concatenation). The JWE model proposed by Yu et al. (2017) is also based on CBOW and it utilizes character and sub-character level information. They construct a dictionary that maps each Chinese character to its sub-character components. As Figure  1 shows, words together with the characters and sub-character components within the context window are all used to predict the target word. The additional semantic information provided by character and subcharacters are shown to improve over word representation, especially in addressing outof-vocabulary words.  (Yu et al., 2017). w i−1 and w i+1 are context words. c i−1 and c i+1 represent characters in context words. s i−1 and s i+1 represent sub-characters of context characters and s i is the sub-character of target word w i .

Data and Preprocessing
We use all the parallel data provided by WMT for the zh-en translation task, including the News Commentary v14, UN Parallel Corpus V1.0 and the CWMT corpora. In addition, the Common Crawl Corpus from WMT is used as monolingual data to pre-train the embeddings. We use the newsdev2018 and newsdev2017 as validation set and the newstest2019 as our test data. We tokenize English sentences with the Moses tokenizer (Koehn et al., 2007). On the Chinese side we use Jieba for Chinese word segmentation. 1 The data preprocessing consists of filtering sentences to be added to the parallel training corpus by rules and by alignment score. Following the preprocessing criteria from submissions in previous years (Xu and Carpuat, 2018;Stahlberg et al., 2018;Haddow et al., 2018), we filter the training data based on the following criteria: • The length of sentences in both languages must be between 4 and 50.
• The maximum length ratio of sentence pairs is 1.3.
• Chinese sentences with no Chinese character are filtered out.
• English sentences with no English character are filtered out.
• Same source and target sentences are removed.
• Sentences should not contains HTML tags.
• Sentence pairs with alignment score above -65 are removed. 2 The fast_align toolkit 3 is used to calculate the alignment score for the parallel data. After the filtering, 10.38M sentence pairs are used as training data. We apply Byte-pair Encoding (BPE) (Sennrich et al., 2016b) with 30,000 merge operations on the English sentences. For Chinese sentences, we segment them into different granularity levels, including words, subwords via BPE and characters. In the character level setting, only Chinese words are separated and each character is treated as a single token. The training texts for models with pre-trained embeddings is the same as baseline, which use words as basic units.

Baseline
The baseline model is based on the bidirectional RNN architecture with attention (Bahdanau et al., 1 https://github.com/fxsjy/jieba 2 We tried different filter strategies and found this criterion gives a better performance than others. 3 https://github.com/clab/fast_align 2015). Our models are built with OpenNMT-py (Klein et al., 2017). We follow the hyperparameter setting of Deep RNN from Xu and Carpuat (2018) and use a four-layer LSTM for both the encoder and decoder. The embeddings and hidden layer size are limited to 512. We use the Adam optimizer (Kingma and Ba, 2015) with initial learning rate of 0.0005. We apply label smoothing (Szegedy et al., 2016) and dropout (Srivastava et al., 2014) of 0.1 to avoid overfitting. We use the multi-layer perception (mlp) attention as in (Bahdanau et al., 2015). The batch size is 4096 tokens per batch and the models are selected based on best performance on the validation set. All our models are trained on a GTX 1080Ti GPU.

Pre-trained Embeddings
We apply pre-trained embeddings to the two submitted systems. The character level and subcharacter level pre-trained embeddings are trained with CWE (Chen et al., 2015) and JWE (Yu et al., 2017) respectively. We trained the embeddings on the Common Crawl Corpus provided by WMT19 and fine-tuned them on the task data when training the RNN. The preprocessing for monolingual data includes Chinese word segmentation and removal of non-Chinese characters. Apart from the pretrained embeddings, the hyperparameters of the two submissions are the same as in the baseline system.

Result and Analysis
We use the CharacTER.py 4 script for Charac-TER score calculation and multeval 5 (Clark et al., 2011) to calculate BLEU, METEOR and TER scores. The evaluation results for models on word, subword and character level are presented in Table 1.
The model with BPE applied on both source and target languages (bpe2bpe) achieves higher score than other single models, with an increase of +1.18 BLEU score over the baseline system. The two models (baseline+cwe, baseline+jwe) utilizing character and sub-character information are based on pre-trained embeddings with CWE and JWE as described in Section 2. We use the source training text for the pre-trained embeddings to prevent the introduction of noise. As we can see from the BLEU scores, the model with JWE pre-4 https://github.com/rwth-i6/CharacTER 5 https://github.com/jhclark/multeval trained embeddings shows similar performance to the baseline system while the model with CWE embeddings on character level shows a marginal decrease. The METEOR and TER score presents similar trends to BLEU, whereas from the evaluation of CharacTER scores the introduction of pre-trained embeddings on both character and subcharacter levels shows better performance than the baseline.
It can also be seen from the comparison on BPE-based models that the model with CWE embeddings performs slightly worse than the bpe2bpe model, which operates on BPE on both source and target languages. The results according to CharacTER show that finer granularity embeddings can benefit the model in character level evaluations. The char2bpe model shows the worst performance according to BLEU scores, whereas the CharacTER score of this model is higher than that of other word level models. Finally, when we ensemble the baseline and four models with JWE embeddings pre-trained on different iterations, the BLEU score shows an increase of +1.26 BLEU over the baseline.
The two models with stars (apprentice-c and apprentice-g) are our official shared task submissions, with the first one operating on character level and the second, on glyph (sub-character) level. The apprentice-c model uses the CWE pretrained embeddings while the apprentice-g uses JWE embeddings. For the first, we train the pre-trained embeddings on the monolingual data (Common Crawl) and then fine-tune it on filtered parallel data during the training of RNN models. Note that we did not use back-translation to aug-ment the training data and due to time limit we apply a relatively larger learning rate than previous work to boost training speed, therefore our systems achieve relatively lower score than the previous work (Xu and Carpuat, 2018). The CWEbased model shows a better BLEU score than the baseline model. The lower performance for the apprentice-g model might have resulted from insufficient training epochs for the JWE embeddings. Due to time restrictions, we did not submit the system with the best word embeddings. In the additional experiments after the task deadline, we fine-tuned the models on the best word embeddings version and achieve a higher BLEU score of 17.43 for the apprentice-g(best) model. The Char-acTER score for the fine-tuned model is lower than other models except the two with BPE. Generally, the sub-character level models perform better than the word level and character level models.
6 Additional Experiments

Evaluating Embeddings
We have tried additional experiments to evaluate the effect of character and subcharacter level pretrained embeddings. Table 2 presents the model performance with respect to the embeddings performance in traditional word similarity and analogy tasks. We use the wordsim-240 and wordsim-297 dataset and the analogy dataset from Chen et al. (2015) for word similarity and analogy evaluation respectively. We use the evaluation script in JWE 6 for both evaluations.
From Table 2, we can see that among all models with JWE pre-trained embeddings, the one with  Table 2: Comparison of model performance and word embeddings performance. The evaluation on wordsim-240 and wordsim-297 test set shows Spearman correlation between the pre-trained embedding and human judgements. The performance on analogy indicates accuracy on analogy reasoning in "a:b::c:?" format. The number after the embeddings type represents number of training iterations.
10 iterations performs the best. When the embeddings are trained over 20 iterations, the BLEU score starts to decrease. The same pattern can be found on the CWE-based models. However, the model with 5-iteration embeddings achieves the highest BLEU score among all CWE-based models. From the embeddings performance on the analogy task, excluding the cwe5 model, we find that the embeddings performance correlates with BLEU scores. When comparing the CWEbased models with the JWE-based models, we see that on both translation quality and word embeddings evaluations, the model on finer granularity performs best.

Effect of Corpus Size
Another experiment was done to compare the effect of pre-trained embeddings on different corpora sizes. We train the word embeddings with best iteration setting and train the RNN model on different corpora sizes. Smaller corpora are created by taking 25% and 50% of the original corpus.  It can be seen from Table 3 that with smaller parallel training corpora the introduction of the pre-trained word embeddings has a more marked positive influence. When the dataset is reduced to half, all the three models show a decrease in BLEU score. However, the gap between the baseline and the cwe-based model is smaller. When the dataset is further limited to 25%, both models with pretrained embeddings perform better than the baseline, whose score does not change. Although it seems that the pre-trained embeddings, even with sub-character level semantic information involved, could only benefit marginally on the whole training data, the introduction of extra semantic information might play a more important role when the parallel training resources are limited. Here we measure the performance of models with varying sentence lengths, as shown in Figure 2. The test set is seperated into 8 subsets based on the sentence lengths and models are evaluated on each subset, the x-axis in Figure 2 represents sentence length intervals. We see that the two models with embeddings trained on a larger monolingual corpus perform better than the other models in medium-length sentences (between 30 and 50). The apprentice-c model, which uses CWE embeddings operating on character level, greatly outperforms the other models on short sentences with length less than 10. Since the sentence length is short, the tokens in the sentence are mostly composed of one or two characters, thus the model with character-based embeddings has an advantage. Regarding the two models with embeddings trained without extra monolingual data, both models show good performance on medium length sentences but perform poorly on long sentences. The introduction of pre-trained embeddings can increase the models' preference to generate shorter sentences, resulting in the model achieving lower BLEU score on long sentences.

Analysis of Model Perplexity
In order to understand the effect of pre-trained embedding on target language model, we calculate the model perplexity on the test data with models on different corpus size. The result is represented in Table 4. The model with JWE pretrained embeddings performs better on all corpus sizes, having a lower perplexity, though the difference is marginal. Similar result as the BLEU evaluation shows that the pre-trained embeddings benefit model performance on smaller corpus sizes.

Model
Perplexity

Transformer Models
Besides the RNN model, we also experimented with pre-trained embeddings and the transformer architecture. We follow the hyperparameter setting from Vaswani et al. (2017), limiting the embeddings to 512 dimensions. We compare the transformer models with and without pre-trained embeddings. The results are presented in Table 5.
From the evaluation results on BLEU and Charac-TER, the transformer models without pre-trained embeddings show better performance. We find it interesting that the embedding pre-trained with CWE decrease the performance severely, leading to a reduction of -3.85 BLEU score from the model without it. The introduction of finer granularity embeddings might not benefit the transformer performance. We hypothesize that the pretrained embedding enhanced by character and subcharacter infomation might conflict with the fixed positional encoding used in transformer.

Conclusion
This paper describes our NMT models with pretrained embeddings operating on character and sub-character levels. We participated in the WMT19 zh-en news translation shared task and submitted two systems with embeddings trained on monolingual corpus. We experimented with the effect of using fine-grained pre-trained embeddings and showed the potential benefit of using them. In additional experiments, we find that using pre-trained embeddings can better benefit the translation models when the parallel training data is limited.


Hyperspectral Data Augmentation

Data augmentation is a popular technique which helps improve generalization capabilities of deep neural networks. It plays a pivotal role in remote-sensing scenarios in which the amount of high-quality ground truth data is limited, and acquiring new examples is costly or impossible. This is a common problem in hyperspectral imaging, where manual annotation of image data is difficult, expensive, and prone to human bias. In this letter, we propose online data augmentation of hyperspectral data which is executed during the inference rather than before the training of deep networks. This is in contrast to all other state-of-the-art hyperspectral augmentation algorithms which increase the size (and representativeness) of training sets. Additionally, we introduce a new principal component analysis based augmentation. The experiments revealed that our data augmentation algorithms improve generalization of deep networks, work in real-time, and the online approach can be effectively combined with offline techniques to enhance the classification accuracy.

I. INTRODUCTION
Hyperspectral satellite imaging (HSI) captures a wide spectrum (commonly more than a hundred of bands) of light per pixel (forming an array of reflectance values). Such detailed information is being exploited by the remote sensing, pattern recognition, and machine learning communities in the context of accurate HSI classification (elaborating a class label of an HSI pixel) and segmentation (finding the boundaries of objects) in many fields [1]. Although the segmentation techniques include conventional machine learning algorithms (both unsupervised [2] and supervised [1], [3]), deep learning based techniques became the main stream [4]- [11]. They encompass deep belief networks [4], [7], recurrent neural networks [8], and convolutional neural networks (CNNs) [5], [6], [9]- [11].
Deep neural networks discover the underlying data representation, hence they do not require feature engineering and can potentially capture features which are unknown for humans. However, to train such large-capacity learners (and to avoid overfitting), we need huge amount of ground truth data. Acquiring such training sets is often expensive, timeconsuming, and human-dependent. These problems are very important real-world obstacles in training well-generalizing models (and validating such learners) faced by the remote sensing community-they are manifested by a very small number of ground truth benchmark HSI sets (there are approx. 10 widely-used benchmarks, with the Salinas Valley, Pavia University, and Indian Pines scenes being the most popular).
To combat the problem of limited, non-representative, and imbalanced training sets, data augmentation can be employed.
It is a process of synthesizing new examples following the original data distribution. Since such enhanced training sets can improve generalization of the learners, data augmentation may be perceived as implicit regularization. In computer vision tasks, data augmentation often involves simple affine (e.g., rotation or scaling) and elastic transforms of the image data [12]. These techniques, albeit applicable to HSI, do not benefit from all available information to model useful data.

A. Related literature
The literature on HSI data augmentation is fairly limited (not to mention, only one of the deep-learning HSI segmentation methods discussed earlier in this section used augmentation-simple mirroring of training samples-for improved classification [10]). In [13], the authors calculated perspectral-band standard deviation (for each class) in the training set. The augmented samples are later drawn from a zero-mean multivariate normal distribution N (0, αΣ), where Σ is a diagonal matrix with the standard deviations (for all classes) along its main diagonal, and α is a hyper-parameter (scale factor) of this technique. Albeit its simplicity, this augmentation was shown to be able to help improve generalization.
Li et al. utilized both spectral and spatial information to synthesize new samples in their pixel-block augmentation [14]. Two data-generation approaches: (i) Gaussian smoothing filtering alongside (ii) label-based augmentation were exploited in [15]. The latter technique resembles weak-labeling [16], and builds on an assumption that neighboring HSI pixels should share the same class (the label of a pixel propagates to its neighbors, and these generated examples are inserted into the training set). Thus, it may introduce wrongly-labeled samples.
Generative adversarial networks (GANs) have already attracted research attention in the context of data augmentation due to their ability of introducing invariance of models with respect to affine and appearance variations. GANs model an unknown data distribution based on the provided samples, and they are composed of a generator and discriminator. A generator should generate data samples which follow the underlying data distribution and are indistinguishable from the original data by the discriminator (hence, they compete with each other). In a recent work [17], Audebert et al. applied GAN conditioning to ensure that the synthesized HSI examples (from random distribution) belong to the specified class. Overall, all of the state-of-the-art HSI augmentation methods are aimed at increasing the size and representativeness of training sets which are later fed to train the deep learners.

B. Contribution
In this letter, we propose a novel online augmentation technique for hyperspectral data (Section II-A)-instead of arXiv:1903.05580v1 [cs.CV] 13 Mar 2019 synthesizing samples and adding them to the training set (hence increasing its size which adversely affects the training time of deep learners), we generate new examples during the inference. These examples (both original and artificial) are classified using a deep net trained over the original set, and we apply the voting scheme to elaborate the final label. To our knowledge, such online augmentation has not been exploited in HSI analysis so far (test-time augmentation was utilized in medical imaging [18], where the authors applied affine transforms and noise injection into brain-tumor images for better segmentation). Also, we introduce principal component analysis (PCA) based augmentation (Section II-B) which may be used both offline (before the training) and online. This PCAbased augmentation simulates data variability, yet follows the original data distribution (which GANs are intended to learn [17], but they are not applicable at test-time).
Our rigorous experiments performed over HSI benchmarks revealed that the online approach is very flexible-different augmentation techniques can be exploited in this setting (Section III). The results obtained using a spectral CNN indicated that the test-time augmentation significantly improves abilities of the models when compared with those trained using the original sets, and augmented using other algorithms (also, we compared our CNN with a spectral-spatial CNN from the literature whose capacity is much larger [11]). Our online approach does not sacrifice the inference speed and allows for real-time classification. We showed that the proposed PCA augmentation is extremely fast, and ensures the highest-quality generalization of the deep models for all data-split scenarios. Finally, we demonstrated that the offline and online augmentations can be effectively combined for better classification.

II. PROPOSED HYPERSPECTRAL DATA AUGMENTATION A. Online Hyperspectral Data Augmentation
Our online (test-time) data augmentation involves synthesizing A artificial samples x 1 , x 2 , . . . , x A for each incoming example x during the inference. We traverse the neighborhood of the original example and try to mitigate potential inputdependent uncertainty of the deep model. In contrast to the offline augmentation techniques, the test-time augmentation does not cause increasing the training time of the network, and it does not require defining the number of synthesized samples beforehand (also, for the majority of specific augmentation algorithms, the operation time of a trained learner would not be significantly affected since the inference is fast). We build upon the theory of ensemble learning, where elaborating a combined classifier (encompassing several weak learners) delivers high-quality generalization (it is an efficient regularizer). Here, by creating A artificial data points, we implicitly form a homogeneous ensemble of deep models (trained over the same training set T ). The final class label is elaborated using majority voting (with equal weights) over all (A + 1) samples (for low ensemble confidence, when two or more classes receive the same number of votes, we perform soft votingwe average all class probabilities, and the final class label corresponds to the class with the highest average probability).
The proposed online HSI augmentation may be considered to be a meta-algorithm, in which a specific data augmentation method is applied to synthesize samples on the fly. Although we exploited the noise injection based approach [13], and our principal component analysis based technique (see Section II-B) in this work, we anticipate that other augmentations which are aimed at modifying an existent sample can be straightforwardly utilized here. Finally, the online augmentation may be coupled with any offline technique (Section III).

B. Principal component analysis based data augmentation
In this work, we propose a new augmentation method based on PCA. Let us consider a training set T of N HSI pixels t i , where i = 1, 2, . . . , N , and each t i is b-dimensional (b denotes the number of bands in this HSI). PCA extracts a set of b (b b) projection directions (vectors) by maximizing the projected variance of a given b-dimensional dataset-the first principal component (P C 1 ) accounts for as much of the data variability as possible, and so forth. First, we center the data at the origin (hence, we subtract the average samplē is the matrix with the non-increasingly ordered eigenvalues along its main diagonal, and is the matrix with b corresponding eigenvectors (principal components) as columns. Finally, b principal components form an orthogonal base, and each sample t can be projected onto a new feature space: t = Φ T t. Importantly, each sample t can be projected back to its original space: The first step of our PCA-based data augmentation involves transforming all training samples using PCA (trained over T ). Afterwards, the first principal component 1 P C 1 (of each sample) is multiplied by a random value α drawn from a uniform distribution U(α min , α max ), where α min and α max are the hyper-parameters of our method (α is drawn independently for all original examples). This process is visualized in Fig. 1-we can observe that the synthesized examples (Fig. 1c) preserve the original data distribution (Fig. 1b) projected onto a reduced feature space, and preserve inter-class relationships. Finally, these samples are projected back onto the original space (using all principal components to ensure correct mapping), and they are added to the augmented T (if executed offline). This PCAbased augmentation can be applied in both offline and online settings (in both cases, PCA is trained over the original T ).

III. EXPERIMENTS
The experimental objective was to verify the impact of data augmentation on the deep model generalization. For online augmentation, we applied our PCA augmentation (PCA-ON), and noise injection (Noise-ON) [13], whereas for the offline setting, we used our PCA-based method (PCA), generative adversarial nets (GAN) [17], and noise injection (Noise) [13]. GAN cannot be used online, since it does not modify an incoming example, but rather synthesizes samples which follow an approximated distribution. Finally, we combined online and offline augmentation in PCA/PCA-ON (PCA augmentation is used to both augment the set beforehand, and to generate new samples at test time), and GAN/PCA-ON. For each offline technique, we at most doubled the number of original samples (unless that number would exceed the most numerous classin such case, we augment only by the missing difference). For online augmentation, we synthesize A = 4 samples, and for PCA and PCA-ON, we had α min = 0.9 and α max = 1.1.
We exploit our shallow (thus resource-frugal) 1D spectral CNN (Fig. 2) coded in Python 3.6 [19]. Larger-capacity CNNs require longer training and infer slower, hence are less likely to be deployed for Earth observation, especially on board of a satellite. The training (ADAM, learn. rate of 10 −4 , β 1 = 0.9, and β 2 = 0.999) stops, if after 15 epochs the validation set V (random subset of T ) accuracy plateaus. We train and validate the deep models using: (1) balanced T sets with random pixels (B), (2) imbalanced T sets with random pixels (IB), and (3) our patched sets (P) [19] (for fair comparison, the numbers of pixels in T and Ψ for B and IB are close to those reported in [11]). We also report the results obtained using a spectral-spatial CNN (3D-CNN) [11], trained over the original T (3D-CNN-in contrast to our CNN-does suffer from the training-test information leak problem, and the 3D-CNN results over B and IB are over-optimistic [19]). The results (over the test sets Ψ) obtained for the Salinas Valley, Indian Pines, and Pavia University datasets are gathered in Table I. Introducing augmented samples (in both online and offline settings) helped boost generalization abilities of the deep models in the majority of cases (even up to more than 8% of OA for GAN, PCA, and PCA/PCA-ON in Salinas, B; only the Noise offline augmentation deteriorated both OA and AA for P). Interestingly, exploring the local neighborhood randomly (Noise and Noise-ON) can notably deteriorate OA and AA. It usually occurs for under-represented classes (e.g., C7 in Pavia) since their examples lay close to other-class examples in the discovered feature space (therefore, they can be easily "confused" with each other). This problem is addressed by the data-distribution analysis in our PCA-based augmentations. Coupling offline and online augmentation (PCA/PCA-ON and GAN/PCA-ON) gave consistent high-quality results over all sets and all training-test splits, and dealt well with the HSI imbalance (in P, we did not ensure that examples of all classes are included in the original T , thus P is very challenging [19]).
To verify the statistical significance of the results (and see if the differences in the average per-class accuracy are important in the statistical sense), we executed two-tailed Wilcoxon's tests for each dataset split (B, IB, and P) over per-class AA for all HSI. The results reported in Table II show that applying HSI data augmentation is beneficial in most cases and delivers significant improvements in accuracy. GAN did equally well as e.g., PCA, PCA-ON, and our combined PCA/PCA-ON and GAN/PCA-ON for B, as Noise-ON, PCA-ON, and GAN/PCA-ON for IB, and as Noise-ON and PCA-ON for P. It indicates that employing time-consuming and complex deep-learning engines for data augmentation not necessarily brings larger improvements in the performance of the deep models.
Our combined approaches (PCA/PCA-ON and GAN/PCA-ON) were stable and consistently ensured high-quality generalization (as shown in Table I) of the deep models over all splits. This stability is also manifested in Table III, where we summarize the results across all sets (although PCA gave the best accuracy for B, the differences between PCA and PCA/PCA-ON and GAN/PCA-ON are not statistically significant). We can appreciate that our PCA-based augmentation (offline, online, or combined) allowed us to obtain the best generalizationvery intuitive PCA-based data-distribution analysis for synthesizing samples outperformed or worked on par with GAN in the case of difficult (small and imbalanced) sets. Finally, our CNN surpassed the accuracy elaborated using a significantly larger 3D-CNN from the literature (with a bigger capacity) for P (note that the results obtained using 3D-CNN for B and IB are over-optimistic due to the intrinsic training-test information leak problem, hence they cannot be considered reliable [11]).    To gain better insights into the augmentation performance (and its potential overhead imposed on the deep models in terms of training and/or test times), we collected the average execution times of the most important steps of the investigated methods in Table IV. It can be observed that training of GANs is very time-consuming (it was executed using NVIDIA GeForce GTX 1060), and is of orders of magnitude higher than pre-processing in other offline techniques (PCA and Noise). Although all offline augmentations affect the training time of deep networks, these differences are not dramatic. Finally, the online augmentation allowed us to classify test pixels in realtime (note that we report the inference time in ms). IV. CONCLUSIONS In this letter, we introduced a new online HSI data augmentation approach which synthesizes examples at test time. It is in contrast to other state-of-the-art hyperspectral data augmentation techniques that work offline (i.e., before the deep-network training to increase the training set cardinality and representativeness). Our experimental study, performed over three HSI benchmark sets (with different training-test data splits) and coupled with statistical tests revealed that our online augmentation is very flexible (different augmentations can be applied here), improves the generalization abilities of deep neural networks, and works in real-time. Also, we showed that combining online and offline augmentation leads to consistently well-performing models. Finally, we proposed a principal component analysis based augmentation which operates extremely fast, synthesizes high-quality data, outperforms other augmentations for small and imbalanced sets, and is applicable in online and offline settings.