extfieldsofstudy:
Computer Science
Engineering

text:
Hybrid Contact Preintegration for Visual-Inertial-Contact State Estimation Using Factor Graphs

The factor graph framework is a convenient modeling technique for robotic state estimation where states are represented as nodes, and measurements are modeled as factors. When designing a sensor fusion framework for legged robots, one often has access to visual, inertial, joint encoder, and contact sensors. While visual-inertial odometry has been studied extensively in this framework, the addition of a preintegrated contact factor for legged robots has been only recently proposed. This allowed for integration of encoder and contact measurements into existing factor graphs, however, new nodes had to be added to the graph every time contact was made or broken. In this work, to cope with the problem of switching contact frames, we propose a hybrid contact preintegration theory that allows contact information to be integrated through an arbitrary number of contact switches. The proposed hybrid modeling approach reduces the number of required variables in the nonlinear optimization problem by only requiring new states to be added alongside camera or selected keyframes. This method is evaluated using real experimental data collected from a Cassie-series robot where the trajectory of the robot produced by a motion capture system is used as a proxy for ground truth. The evaluation shows that inclusion of the proposed preintegrated hybrid contact factor alongside visual-inertial navigation systems improves estimation accuracy as well as robustness to vision failure, while its generalization makes it more accessible for legged platforms.

I. INTRODUCTION AND RELATED WORK
Long-term state estimation and mapping for legged robots requires a flexible sensor fusion framework that allows for reducing drift and correcting past estimates as the robot perceives new information. During long-term missions, odometry systems can drift substantially since the absolute position and yaw (rotation about gravity) are unobservable [1], [2], leading to an unbounded growth in the covariance of the estimate and an undesirable expansion of the search space for data association tasks. The factor graph smoothing framework [3]- [6] offers suitable machineries for building such systems in which real-time performance is achieved by exploiting the sparse structure of the Simultaneous Localization and Mapping (SLAM) problem [7], [8]. In addition, the incorporation of loop-closures [8] into the graph, upon availability, is convenient.
Legged robot perception often involves fusing kinematic odometry, Inertial Measurement Unit (IMU) data, and visual/depth measurements to infer the robot's trajectory, controller inputs such as velocity, and calibration parameters [9]- [12]. The challenge in such perception problems is the rigorous real-time performance requirements in legged The authors are with the College of Engineering, University of Michigan, Ann Arbor, MI 48109 USA {rosshart, maanigj, ganlu, bjhuang, grizzle, eustice}@umich.edu. Fig. 1: Experiments were conducted on a Cassie-series robot designed by Agility Robotics in an indoor laboratory environment. The motion capture system is used to record the robot trajectory as a proxy for ground truth data. The Cassie-series robot has 20 degrees of freedom, 10 actuators, joint encoders, an IMU, and mounted with a Multisense S7 stereo camera. robots arising from their direct and switching contact with the environment [1], [13]- [16]. Furthermore, kinematic leg odometry involves estimating relative transformations and velocity using kinematic and contact information, which can be noisy due to the encoder noise and foot slip [2], [17].
In our previous work [18], we made progress towards a long-term, legged robot state estimator by developing two novel factors that integrate the use of forward kinematics and the notion of contact between the robot and the environment into the factor graph framework. While these novel factors improved the estimator's performance, new challenges emerged due to frequent switching contacts as the robot navigates through an environment. In particular, a new node needs to be added to the graph every time contact was made or broken. This leads to an larger number of optimization variables, which increases the complexity of the graph and ultimately slows down the estimator. This problem is only exacerbated as the number of contact points increases (quadraped or hexapod robots).
In this paper, we generalize the idea of the preintegrated contact factor in [18] to a hybrid preintegrated contact factor that alleviates the effect of frequently changing contacts with the environment. We develop a novel method for preintegrating contact information though an arbitrary number of contact switches. The proposed hybrid modeling approach reduces the number of required variables in the nonlinear optimization problem by only requiring new states to be added alongside camera or selected keyframes. The present work has the following contributions: i. A generic forward kinematic factor modeled in SE (3) using the manipulator Jacobian that supports both prismatic and revolute joints; the factor incorporates noisy encoder measurements to estimate an end-effector pose at any time-step. ii. A hybrid preintegrated contact factor modeled in SE (3) that allows for an arbitrary number of contact switches between camera or selected keyframes. iii. Real-time implementation and experimental evaluation of the derived factors on a Cassie-series biped robot in an indoor laboratory environment where ground truth data was collected using a motion capture system.
The remainder of this paper is organized as follows. Section II provides the mathematical background and preliminaries. We formulate the problem using the factor graph approach in Section III. Section IV explains the forward kinematic modeling and derives the forward kinematic factor. The proposed hybrid rigid contact model and hybrid contact preintegration are derived in Section V. Experimental evaluations of the proposed methods on a 3D biped robot (shown in Fig. 1) are presented in Section VI. Finally, Section VII concludes the paper and provides future work suggestions.

II. MATHEMATICAL BACKGROUND AND PRELIMINARIES
We first review the Lie group theory corresponding to the rotation and motion groups [19], [20]. Afterwards, we discuss the optimization technique on matrix Lie groups [21] and the choice of retraction for deriving forward kinematic and hybrid preintegrated contact factors.
Matrices are capitalized in bold, such as in X, and vectors are column-wise in lower case bold type, such as in x.
We denote e 2 Σ e T Σ −1 e. The n-by-n identity matrix and the n-by-m matrix of zeros are denoted by I n and 0 n×m respectively. The vector constructed by stacking x i , ∀ i ∈ {1, . . . , n} is denoted by vec(x 1 , . . . , x n ). The covariance of a random vector is denoted by Cov[·]. Finally, we denote the base (IMU) frame of the robot by B, the world frame by W, and contact frame by C.

A. Matrix Lie Group of Rotation and Motion in R 3
The general linear group of degree n, denoted by GLn(R), is the set of all n × n nonsingular real matrices, where the group binary operation is the ordinary matrix multiplication. The three-dimensional (3D) special orthogonal group, denoted by SO(3) = {R ∈ GL 3 (R) | RR T = I 3 , det R = +1} is the rotation group on R 3 . The 3D special Euclidean group, denoted by is the group of rigid transformations on R 3 . The Lie algebra (tangent space at the identity together with Lie bracket) of SO(3), denoted by so (3), is the set of 3 × 3 skew-symmetric matrices such that for any ω vec(ω 1 , ω 2 , ω 3 ) ∈ R 3 : and (ω ∧ ) ∨ = ω. The Lie algebra of SE(3), denoted by se(3), can be identified by 4 × 4 matrices such that for any ω, v ∈ R 3 and the twist is defined as ξ vec(ω, v) ∈ R 6 : where the wedge operator ( ∧ ) for twist is overloaded. The exponential map exp : se(3) → SE(3) can be used to map a member of se(3) around a neighborhood of zero to a member of SE(3) around a neighborhood of the identity. The logarithm map is the inverse, i.e., log : SE(3) → se(3), and exp(log(H)) = H, H ∈ SE(3). Now we can define the difference between a transformation H ∈ SE(3) and its estimate with a small perturbationH ∈ SE(3) as [20], [22]: where ∧ ∈ se(3). To define the norm and covariance of the error term, we exploit the fact that se(3) is isomorphic to R 6 , i.e., ∧ → ∈ R 6 using the ∨ operator. Therefore, we can define the 6 × 6 covariance matrix conveniently as Σ = Cov[ ]. We use the following adopted simplified notations from [6]: The adjoint representation of a Lie group is a linear map that captures the non-commutative structure of the group. For SE(3), the matrix representation of the adjoint map [20] is For any H ∈ SE(3) and ξ ∈ R 6 , the adjoint map asserts

B. Retraction Map and Optimization on Manifold
Given a retraction mapping and the associated manifold, we can optimize over the manifold by iteratively lifting the cost function of our optimization problem to the tangent space, solving the re-parameterized problem, and then mapping the updated solution back to the manifold using the retraction [21]. For SE(3) we use its exponential map as the natural retraction, R H (δξ) = HExp(δξ), where δξ ∈ R 6 is the twist defined earlier in (1). Therefore, the retractions on the base and the contact poses are defined as These retractions are used to both update the state during optimization and to compute the Jacobians of the residuals.

III. PROBLEM STATEMENT AND FORMULATION
In this section, we formulate the state estimation problem using the factor graph framework where independent sensor measurements can be incorporated by introducing additional factors based on the associated measurement models. The biped robot, shown in Figure 1, is equipped with a stereo camera, an IMU mounted on the torso, joint encoders, and a series of springs used for binary contact detection. Without loss of generality, we assume the IMU and camera are collocated with the base frame of the robot.

A. State Representation
The state variables (at time t i ) include the 3D body pose, H WB (t) ∈ SE(3), and velocity, W v WB (t) ∈ R 3 , in the world frame, the 3D contact pose in the world frame, are the gyroscope and accelerometer biases, respectively. All together, the state at any time-step i is a tuple as follows: Further, it is convenient to denote the trajectory of the state variables up to time-step k by X k k i=1 T i . Foot slip is the major source of drift in kinematic leg odometry. The inclusion of the contact pose in the state tuple allows for isolating the noise at the contact point.

B. Factor Graph Formulation
Let L ij ∈ SE(3) be a perceptual loop-closure measurement relating poses at time-steps i and j (j > i) computed from an independent sensor, e.g. using a point cloud matching algorithm. The forward kinematic measurements at timestep i are denoted by F i . The IMU and contact sensors provide measurements at higher frequencies. Between any two time-steps i and j, we denote the set of all IMU and contact measurements by I ij and C ij , respectively. Let K k be the index set of time-steps (or keyframes) up to time-step k. We denote the set of all measurements up to time-step k by By assuming the measurements are conditionally independent and are corrupted by additive zero mean white Gaussian noise, the posterior probability of the full SLAM problem can be written as The Maximum-A-Posteriori (MAP) estimate of X k can be computed by solving the following optimization problem: in which due to the noise assumption mentioned earlier is equivalent to the following nonlinear least-squares problem: where r 0 and Σ 0 represent the prior over the initial state and serves to anchor the graph, r Lij , r Iij , r Fi , r Cij are the residual terms associated with the loop closure, IMU, forward kinematic, and contact measurements respectively, i.e. Fig. 2: In this paper, we refer to two separate forward kinematics functions. The pose of the current contact frame relative to the base frame is denoted by H BC . When the robot has multiple points of contact with the environment, it is possible to transfer this contact from from one frame to another. This transfer of contact is captured by the homogeneous transform H C -C + .
the error between the measured and predicted values given the state, and Σ Lij , Σ Iij , Σ Fi , Σ Cij are the corresponding covariance matrices.

IV. FORWARD KINEMATICS
Forward kinematics refers to the process of computing the relative pose transformation between two frames of a multilink system. Each individual joint displacement describes how the child link moves with respect to the parent one. This joint displacement can either be an angle (revolute joints) or a distance (prismatic joints).
Let α ∈ R N denote the vector of joint displacements for a general robot. Without loss of generality, we define a base frame on the robot, denoted B, that is assumed to be collocated with both the IMU and the camera frames. When the robot is in contact with the static environment, we can also define a contact frame, denoted C, on the robot at the point of contact. The homogeneous transformation between the base frame and the contact frame is defined by: where R BC (α) and B p BC (α) denote the relative orientation and position of the contact frame with respect to the base frame.
When there are two points in contact with the static environment, it is possible to "transfer" the contact frame from one point to the other (shown in Figure 2). Let C − denote the old contact frame and C + denote the new contact frame. Then, the homogeneous transformation between the old frame and the new frame is defined by: where R C -C + (α) and C -p C -C + (α) denote the relative orientation and position of the new contact frame with respect to the old contact frame.

A. Measurements
We assume the robot's joints are equipped with a set of joint encoders that can measure the joint displacement. These encoder measurements,α, are assumed to be corrupted with Gaussian white noise. This is an explicit measurement coming from physical sensors located on the robot.
The geometric (or manipulator) Jacobian, denoted J (α), provides a method for computing the angular and linear velocity of an end-effector given the vector of joint velocities [19]. In a similar manner, we can use the Jacobian to map incremental angle displacements to changes in the end-effector pose. Let J BC (α) denote the body manipulator Jacobian of the forward kinematics function (5). Then, the following relationship holds: where δα and δt are incremental encoder and time quantities and C ξ BC denotes the vector of angular and linear velocities of the contact frame due to δα (measured in the contact frame). We can use this relationship (8) to factor out the noise from the forward kinematics equations. Up to a first order approximation, (5) can be factored as: A similar approximation can be found for (6), albeit with a different Jacobian.
Remark 1. In general, the manipulator Jacobian can be derived as spatial or body Jacobian [19], and based on this choice, the noise can appear on the left or right side of the rotation/rigid body transformation, respectively.

B. Forward Kinematic Factor
The goal of this section is to derive a general forward kinematic factor that can be used in the factor graph framework. This will be a unary factor that relates two poses through the forward kinematics equations while accounting for encoder noise. We derive it here for relating the base and contact frames.
The orientation and position of the contact frame with respect to the world frame are given by: Substituting in the state variables at time t i (4), the forward kinematic equation (5) yields: We can now use the first order approximation (9) to factor out the encoder noise to give the following expression: Defining the zero-mean white Gaussian forward kinematic noise term, δc i J BC (α i )w α i , allows us to write out the forward kinematic measurement model: where the forward kinematic noise is characterized by δc i ∼ N (0 6×1 , Σ Fi ). The residual error is defined in the tangent space and can be written as: The covariance is computed through the following linear transformation: where J BC (α i ) is the body manipulator Jacobian evaluated at the current encoder measurements and Σ α i denotes the encoder covariance matrix at time t i .

V. HYBRID CONTACT PREINTEGRATION
A continuous hybrid dynamical system, H, can be defined with a continuous dynamics function, f (·), a discrete transition map, ∆(·), and a switching surface, S [23]. Trajectories of the hybrid dynamical system evolve according to the continuous dynamics, until the switching surface it hit. At those moments, the state gets mapped through the discrete transition map, after which the trajectory continues according to the continuous dynamics again. The general form of this system can be expressed as follows: As long as the number of contact points is greater than or equal to one, the switching contact frame dynamics can be modeled as a hybrid system: where the switching surface, S is simply modeled as the set of all times where contact is switched from one frame to another. Since the sensor measurements are coming in at discrete time-steps, we perform Euler integration from time t to t+∆t to discretize the continuous hybrid contact dynamics (18), forming the following discrete hybrid system: Physically, the continuous dynamics function, f , describes how a single contact frame moves over time while contact is maintained. When a new contact frame is detected, the new contact pose can be computed by applying the transition map, ∆(·), which describes the homogeneous transformation between the old and new contact frames.

A. Measurements
The angular and linear velocity of the contact point is an implicit measurement that is inferred through a binary contact sensor: specifically, when this sensor indicates contact, the pose of the contact frame is assumed to remain fixed with respect to the world frame, i.e. the measured velocities are zero. In order to accommodate potential contact slip, Fig. 3: In the factor graph framework, the robot's state along a discretized trajectory denoted by red circles. Each independent sensor measurement is a factor denoted by lines that constraints the state at separate time-steps. The proposed hybrid contact factor (shown on top) allows preintegration of highfrequency contact data through an arbitrary number of contact switches. In this example, there are two contact switches, where the robot moves from left-stance (L) to right-stance (R), then back to left stance. this measured twist is modeled to be corrupted with white Gaussian noise, namely Using the state variables (4), forward kinematics definition (6), encoder measurements (7), and contact measurements, the hybrid contact dynamics (19) can be written as: where w vd k , the discrete time contact noise, is computed using the sample time, Cov[w ξd (t)] = 1 ∆t Cov[w ξ (t)].

Remark 2.
H C -C + (α k ) will depend on the specific contact frames, C − and C + , at time t k . For example, the forward kinematics function to switch from left foot contact to a right foot contact will be different than the one used to switch from right foot to left foot.

B. Preintegrating Contact Pose
The goal of this section is to formulate a general hybrid preintegrated contact factor which relates the contact pose at t i to the contact pose at t j . The hybrid nature of this factor comes from the potential switching of contact that occurs naturally in legged locomotion. Following the work on IMU preintegration theory [6], [24], we preintegrate the high-frequency contact measurements to prevent unnecessary computation allowing efficient implementation of the factor.
Let S denote the sequence of all time indices associated with contact switches, so that each s i ∈ S represents a single time index where a contact switch occurs. To integrate the hybrid contact dynamics model (20) through these contact switches, we integrate the up to the next contact switch time t si , apply the transition map ∆(·), and continue the integration. For the following derivations, a concrete example (depicted in Figure 3) is used to make the derivation easier to follow. We later extend the derivation to an arbitrary number of contact switches. For the example shown in Figure 3, there are two contact switches between times t i and t j , i.e., |S| = 2. Integrating the discrete hybrid contact model (20) from t i to t j yields: Exp(−w ξd k ∆t) .
All noise terms can be shifted to the right by using (10) to factor the noise from the forward kinematics term and the adjoint relation (2) to shift the measured kinematics to the left.
After multiplying both sides by C −1 i , we arrive at a relative pose expression that is independent of states T i and T j .
where ∆C ij |S| n=1 H C -C + (α sn ) represents the hybrid preintegrated contact measurement, and Exp(−δc ij ) groups all of the noise terms together. This noise term is a product of multiple small rigid body transformations. Therefore, it can be approximated as a summation in the tangent space through iterative use of the Baker-Campbell-Hausdorff (BCH) formula [20] (while keeping only the first order terms).
The hybrid preintegrated contact noise, δc ij , is a summation of zero-mean Gaussian terms, and is therefore also zeromean and Gaussian. It is possible to generalize this noise expression to an arbitrary number of contact switches, however, it becomes much simpler to do so when looking at the iterative propagation form in the following section.

C. Iterative Propagation
It is possible to write both the preintegrated contact measurements, ∆C ij , and the preintegrated contact noise, δc ij , in iterative update forms. This allows the terms to be updated as contact and encoder measurements are coming in. In addition, this form simplifies the expressions and allows for through the use of the adjoint map of the forward kinematics transformation, H C -C + . However, due to noisy encoders, an addition noise term (computed using the manipulator Jacobian) has to be added to compute δc + ik .
the covariance to be conveniently computed. The following proposition generalizes the hybrid preintegrated contact pose and noise iterative propagation to an arbitrary number of contact switches. The proof is given in the supplementary material.

Proposition 1 (Iterative Propagation of Hybrid Contact
Process [25]). Between any two time-steps t i and t j such that j > i, starting with ∆C ii = I 4 and δc ii = 0 3×1 , the hybrid preintegrated contact measurement, ∆C ij , and noise, δc ij , for an arbitrary number of contact switches can be computed iteratively using the following hybrid systems: δH : An abstract representation of these hybrid systems is shown in Figure 4.

D. Rigid Contact Factor
The relative contact pose expression (22) can be used to define the hybrid preintegrated contact measurement model: where the preintegrated contact noise is characterized by δc ij ∼ N (0 6×1 , Σ Cij ). In the factor graph framework, this hybrid preintegrated contact model represents a binary factor that relates the contact frame pose over consecutive time steps. The residual error is defined in the tangent space, and can be written as: The covariance is computed using the hybrid contact noise model (24), starting with Σ Cii = 0 6×6 :

VI. EXPERIMENTAL RESULTS
We now present experimental evaluations of the proposed factors. In the first experiment, we compare three odometry systems composed by Visual-Inertial-Contact (VIC), Inertial-Contact (IC), and Visual-Inertial (VI) factors. Whenever the contact factor is used, it is assumed that forward kinematic factor is also available; therefore, it is not explicitly mentioned for brevity. In the second experiment, we study the effect of losing visual data for a period of time to see how the contact factor can constrain the graph in the absence of a reliable vision system. In the third experiment, we will evaluate the use of terrain factors, where loop-closures can be added to the graph through contact frame poses.

A. Experimental Setup
All experiments are done on a Cassie-series robot designed by Agility Robotics, which has 20 degrees of freedom, 10 actuators, joint encoders, and a VN-100 IMU. A Multisense S7 stereo camera, which contains a separate IMU, is mounted on the top of the Cassie robot as shown in Fig. 1. The robot also has four springs (two on each leg) that can be used as binary contact sensors by thresholding the spring deflection measurements. The Cassie robot has two computers which run MATLAB Simulink Real-Time (for the controller) and Ubuntu (for the estimator) respectively. We use the Robot Operating System (ROS) [26] with the User Datagram Protocol (UDP) to communicate sensor data between the two computers. We also integrate the time synchronization algorithm in [27] into our system to ensure all sensory data is synchronized. The motion capture system developed by Vicon is used as a proxy for ground truth trajectories. The setup consists of 17 motion capture cameras with four markers attached to the robot to track the IMU pose. The dataset contains the stereo images (20 Hz) and IMU data (750 Hz) from the Multisense S7 camera as well as the joint encoders and IMU data from the Cassie robot (at 400 Hz each).
The proposed kinematic and contact factors were implemented using the GTSAM 4.0 library [4]. We utilized the built-in preintegrated IMU factor [6] with the iSAM2 [3] incremental solver. For visual odometry, we used the semidirect visual odometry library (SVO 2.0 [28]). The Multisense camera recorded synchronized stereo images at 20 Hz and IMU measurements in about 750 Hz. SVO processes those measurements in real-time and outputs the pose of the left camera in a fixed world frame, X i , for the current time-step i. The relative transformation of the camera from time-step i to j can be obtained using ∆X ij = X −1 i X j . We selected keyframes approximately every 0.25 seconds.

B. First Experiment: Odometry Comparison
In this experiment, we had Cassie stand in place for about 15 seconds, then slowly walk forwards and backwards along the length of the lab for approximately 45 seconds. The resulting data was used to compare the odometry performance (processed off-line) of different combinations of factors. The results are shown in Figure 5. The odometry estimate from The odometry results from a 60 second walking experiment using a Cassie-series robot. The Visual-Inertial-Contact (VIC) odometry outperformed both the Inertial-Contact (IC), and the Visual-Inertial (VI) odometry. "Ground-truth" data was collected from a Vicon motion capture system. It is important to note that no loop-closures are being performed, which helps to explain the relatively poor odometry from VI. The video of this experiment is provided at https://youtu.be/WDPhdl5g2MQ. VIC, as expected, outperforms all other combinations of factors. The Cumulative Distribution Function (CDF) of the relative position error is shown in Figure 6. The relative position CDF provides a method for analyzing the drift of an odometry estimate.
From Figure 6, it can be seen that VIC has the highest fraction of data corresponding to smaller relative position errors. This means VIC has lowest drift among all odometry systems. When the robot is walking, the hard impacts cause significant camera shake which leads to motion blur in the images. In addition, the lab environment seemingly lacked numerous quality features. These effects help to explain the relatively poor VI odometry performance.

C. Second Experiment: Vision Dropout
One of the main benefits from including the forward kinematic and contact factors is that the state estimator can be more robust to failure of the vision system. In this experiment, we simulate the effects of "vision dropout" by simply ignoring SVO visual odometry data for two 10-second periods of the experimental data described in the previous section. In other words, during a "vision dropout" period, VIC odometry reduces to IC odometry, and VI odometry Fig. 7: When vision data is lost, the covariance of the robot's base pose sharply grows for VI odometry due to the lack of additional measurements to constrain the graph. In contrast, during "vision dropout" periods, the additional contact factors allows the covariance estimate from VIC to remains close to the nominal case. reduces down to inertial (I) odometry. Figure 7 shows the log determinant of the base pose covariance for VIC and VI for these "vision dropout" experiments. The larger the log determinant, the more uncertain the estimator is about the robot's base pose. During the "vision dropout" periods, uncertainty grows for VI odometry. This sharp covariance growth is due to the lack of additional sensor measurements to add into the factor graph. In contrast, the covariance growth for VIC is hardly affected over the same dropout periods.

D. Third Experiment: Terrain Factors for Loop-Closure
Another benefit of adding the proposed contact factors comes from the addition of the contact frame poses into the robot's state. With these new state variables, it becomes simple to place additional constraints that relate the contact pose to a prior map. We test this idea on the collected experiment data by recognizing that the ground was relatively flat in the laboratory. This "zero-height" elevation data serves as our prior map. Figure 8 shows how adding this trivial constraint can reduce position drift in the z-direction. This experiment Fig. 8: Since the contact pose is now part of the estimated state, it is possible to add "terrain factors" that relate this contact pose to a prior map. Adding the simple constraint that the contact frame z-translation is zero (VIC-T) improves the drift in the z-direction when compared to the nominal Visual-Inertial-Contact (VIC) case. simply serves to illustrate the potential for "terrain factors", as the state estimate could be further improved if the robot was actually mapping out the terrain; there was actually a slight slope in the lab (as shown in the motion capture data).

VII. CONCLUSION
We developed a novel method for preintegrating contact information through an arbitrary number of contact switches using a hybrid preintegrated contact factor. The proposed approach was verified through experimental evaluations using a Cassie-series robot where a motion capture system is used as a proxy for ground truth data. Our results indicate that the fusion of contact information with IMU and vision data provide a reliable odometry system for legged robot. Furthermore, we showed that the developed visual-inertialcontact odometry system is robust to occasional vision system failures.
In the future, we plan to incorporate loop closure constraints into our factor graph framework to further improve state estimation, paving the way for long-term mapping on legged robots. We also plan to further investigate the potential utility of "terrain factors" to allow the robot's state to be corrected though detected contact on an estimated map.


Sequence to Sequence Mixture Model for Diverse Machine Translation

Sequence to sequence (SEQ2SEQ) models lack diversity in their generated translations. This can be attributed to their limitations in capturing lexical and syntactic variations in parallel corpora, due to different styles, genres, topics, or ambiguity of human translation process. In this paper, we develop a novel sequence to sequence mixture (S2SMIX) model that improves both translation diversity and quality by adopting a committee of specialized translation models rather than a single translation model. Each mixture component selects its own training dataset via optimization of the marginal log-likelihood, which leads to a soft clustering of the parallel corpus. Experiments on four language pairs demonstrate the superiority of our mixture model compared to SEQ2SEQ model with the standard and diversity encouraged beam search. Our mixture model incurs negligible additional parameters and no extra computation in the decoding time.

Introduction
Neural sequence to sequence (SEQ2SEQ) models have been remarkably effective machine translation (MT) (Sutskever et al., 2014;Bahdanau et al., 2015). They have revolutionized MT by providing a unified end-to-end framework, as opposed to the traditional approaches requiring several submodels and long pipelines. The neural approach is superior or on-par with statistical MT in terms of translation quality on various MT tasks and domains e.g. (Wu et al., 2016;Hassan et al., 2018).
A well recognized issue with SEQ2SEQ models is the lack of diversity in the generated translations. This issue is mostly attributed to the decoding algorithm (Li et al., 2016), and recently to the model (Zhang et al., 2016;Schulz et al., 2018a). The former direction has attempted to design diversity encouraging decoding algorithm, particularly beam search, as it generates translations sharing the majority of their tokens except a few trailing ones. The latter direction has investigated modeling enhancements, particularly the introduction of continuous latent variables, in order to capture lexical and syntactic variations in training corpora, resulted from the inherent ambiguity of the human translation process. 1 However, improving the translation diversity and quality with SEQ2SEQ models is still an open problem, as the results of the aforementioned previous work are not fully satisfactory.
In this paper, we develop a novel sequence to sequence mixture (S2SMIX) model that improves both translation quality and diversity by adopting a committee of specialized translation models rather than a single translation model. Each mixture component selects its own training dataset via optimization of the marginal log-likelihood, which leads to a soft clustering of the parallel corpus. As such, our mixture model introduces a conditioning global discrete latent variable for each sentence, which leads to grouping together and capturing variations in the training corpus. We design the architecture of S2SMIX such that the mixture components share almost all of their parameters and computation.
We provide experiments on four translation tasks, translating from English to German/French/Vietnamese/Spanish. The experiments show that our S2SMIX model consistently outperforms strong baselines, including SEQ2SEQ model with the standard and diversity encouraged beam search, in terms of both translation diversity and quality. The benefits of our mixture model comes with negligible additional parameters and no extra computation at inference time, compared to the vanilla SEQ2SEQ model.

Attentional Sequence to Sequence
An attentional sequence to sequence (SEQ2SEQ) model (Sutskever et al., 2014;Bahdanau et al., 2015) aims to directly model the conditional distribution of an output sequence y ≡ (y 1 , . . . , y T ) given an input x, denoted P (y | x). This family of autoregressive probabilistic models decomposes the output distribution in terms of a product of distributions over individual tokens, often ordered from left to right as, where y <t ≡ (y 1 , . . . , y t−1 ) denotes a prefix of the sequence y, and θ denotes the tunable parameters of the model. Given a training dataset of input-output pairs, denoted by D ≡ {(x, y * ) d } D d=1 , the conditional log-likelihood objective, predominantly used to train SEQ2SEQ models, is expressed as, A standard implementation of the SEQ2SEQ model is composed of an encoder followed by a decoder. The encoder transforms a sequence of source tokens denoted (x 1 , . . . , x N ), into a sequence of hidden states denoted (h 1 , . . . , h N ) via a recurrent neural network (RNN). Attention provides an effective mechanism to represent a soft alignment between the tokens of the input and output sequences (Bahdanau et al., 2015), and more recently to model the dependency among the output variables (Vaswani et al., 2017).
In our model, we adopt a bidirectional RNN with LSTM units (Hochreiter and Schmidhuber, 1997). Each hidden state h n is the concatenation of the states produced by the forward and backward RNNs, h n = [h →n , h n← ]. Then, we use a two-layer RNN decoder to iteratively emit individual distributions over target tokens (y 1 , ..., y T ). At time step t, we compute the hidden representations of an output prefix y ≤t denoted s 1 t and s 2 t based on an embedding of y t denoted M[y t ], previous representations s 1 t−1 , s 2 t−1 , and a context vector c t as, where M is the embedding table, and W and W are learnable parameters. The context vector c t is computed based on the input and attention, where W h , W s , b a , and v are learnable parameters, and a t is the attention distribution over the input tokens at time step t. The decoder utilizes the attention information to decide which input tokens should influence the next output token y t+1 .

Sequence to Sequence Mixture Model
We develop a novel sequence to sequence mixture (S2SMIX) model that improves both translation quality and diversity by adopting a committee of specialized translation models rather than a single translation model. Each mixture component selects its own training dataset via optimization of the marginal log-likelihood, which leads to a soft clustering of the parallel corpus. We design the architecture of S2SMIX such that the mixture components share almost all of their parameters except a few conditioning parameters. This enables a direct comparison against a SEQ2SEQ baseline with the same number of parameters. Improving translation diversity within SEQ2SEQ models has received considerable recent attention (e.g., Vijayakumar et al. (2016); Li et al. (2016)). Given a source sentence, human translators are able to produce a set of diverse and reasonable translations. However, although beam search for SEQ2SEQ models is able to generate various candidates, the final candidates often share majority of their tokens, except a few trailing ones. The lack of diversity within beam search raises an issue for possible re-ranking systems and for scenarios where one is willing to show multiple translation candidates to the user. Prior work attempts to improve translation diversity by incorporating a diversity penalty during beam search (Vijayakumar et al., 2016;Li et al., 2016). By contrast, our S2SMIX model naturally incorporates diversity both during training and inference.
The key difference between the SEQ2SEQ and S2SMIX models lies in the formulation of the conditional probability of an output sequence y given an input x. The S2SMIX model represents P θ (y | x) by marginalizing out a discrete latent variable z ∈ {1, . . . , K}, which indicates the selection of the mixture component, i.e., where K is the number of mixture components. For simplicity and to promote diversity, we assume that the mixing coefficients follow a uniform distribution such that for all z ∈ {1, . . . , K}, For the family of S2SMIX models with uniform mixing coefficients (10), the conditional loglikelihood objective (2) can be re-expressed as: where log(1/K) terms were excluded because they offset the objective by a constant value. Such a constant has no impact on learning the parameters θ. One can easily implement the objective in (11) using automatic differentiation software such as tensorflow (Abadi et al., 2016), by adopting a LogSumExp operator to aggregate the loss of the individual mixture components. When the number of components K is large, computing the terms P θ (y * t | y * <t , x, z) for all values of z ∈ {1, . . . , K} can require a lot of GPU memory. To mitigate this issue, we will propose a memory efficient formulation in Section 3.3 inspired by the EM algorithm.

S2SMIX Architecture
We design the architecture of the S2SMIX model such that individual mixture components can share as many parameters and as much computation as possible. Accordingly, all of the mixture components share the same encoder, which requires processing the input sentence only once. We consider different ways of injecting the conditioning signal into the decoder. As depicted in Figure 1, we consider different ways of injecting the conditioning on z into our two-layer decoder. These different variants require additional lookup tables denoted M 1 , M 2 , or M b . When we incorporate the conditioning on z into the LSTM layers, each lookup table (e.g., M 1 and M 2 ) has K rows and D LSTM columns, where D LSTM denotes the number of dimensions of the LSTM states (512 in our case). We combine the state of the LSTM with the conditioning signal via simple addition. Then the LSTM update equations take the form, for i ∈ {1, 2}. We refer to the addition of the conditioning signal to the bottom and top LSTM layers of the decoder as bt and tp respectively. Note that in the bt configuration, the attention mask depends on the indicator variable z, whereas in the tp configuration that attention mask is shared across different mixture components. We also consider incorporating the conditioning signal into the softmax layer to bias the selection of individual words in each mixture component. Accordingly, the embedding table M b has K rows and D vocab entries, and the logits from (5) are added to the corresponding row of M b as, We refer to this configuration as sf and to the configuration that includes all of the conditioning signals as all.

Separate Beam Search per Component
At the inference stage, we conduct a separate beam search per mixture component. Performing beam search independently for each component encourages diversity among the translation candidates as different mixture components often prefer certain phrases and linguistic structures over each other. Letŷ z denote the result of the beam search for a mixture component z. The final output of our model, denotedŷ is computed by selecting the translation candidate with the highest probability under the corresponding mixture component, i.e., In order to accurately estimate the conditional probability of each translation candidate based on (9), one needs to evaluate each candidate using all of the mixture components. However, this process considerably increases the inference time and latency. Instead, we approximate the probability of each candidate by only considering the mixture component based on which the candidate translation has been decoded, as outlined in (14). This approximation also encourages the diversity as we emphasized in this work.
Note that we have K mixture components and a beam search of b per component. Overall, this requires processing K × b candidates. Accordingly, we compare our model with a SEQ2SEQ model using the same beam size of K × b.

Memory Efficient Formulation
In this paper, we adopt a relatively small number of mixture components (up to 16), but to encompass various clusters of linguistic content and style, one may benefit from a large number of components. Based on our experiments, the memory footprint of a S2SMIX with K components increases by about K folds, partly because the softmax layers take a big fraction of the memory. To reduce the memory requirement for training our model, inspired by prior work on EM algorithm (Neal and Hinton, 1998), we re-express the gradient of the conditional log-likelihood objective in (11) exactly as, where with uniform mixing coefficients, the posterior distribution P (z | x, y * ) takes the form, where z (y | x) = log P θ (y | x, z). Based on this formulation, one can compute the posterior distribution in a few forward passes, which require much less memory. Then, one can draw one or a few Monte Carlo (MC) samples from the posterior to obtain an unbiased estimate of the gradient in (15). As shown in algorithm 1, the training procedure is divided into two parts. For each minibatch we compute the componentspecific log-loss for different mixture components in the first stage. Then, we exponentiate and normalize the losses as in (16) to obtain the posterior distribution. Finally, we draw one sample from the posterior distribution per input-output example, and optimize the parameters according to the loss of such a component. These two stages are alternated until the model converges. We note that this algorithm follows an unbiased stochastic gradient of the marginal log likelihood.

Experiments
Dataset. To assess the effectiveness of the S2SMIX model, we conduct a set of translation experiments on TEDtalks on four language pairs: English→French (en-fr), English→German (en-de), English→Vietnamese (en-vi), and English→Spanish (en-es).
We use IWSLT14 dataset 2 for en-es, IWSLT15  Implementation details. All of the models use a one-layer bidirectional LSTM encoder and a twolayer LSTM decoder. Each LSTM layer in the encoder and decoder has a 512 dimensional hidden state. Each input word embeddings is 512 dimensional as well. We adopt the Adam optimizer (Kingma and Ba, 2014). We adopt dropout with a 0.2 dropout rate. The minibatch size is 64 sentence pairs. We train each model 15 epochs, and select the best model in terms of the perplexity on the dev set.
Diversity metrics. Having more diversity in the candidate translations is one of the major advantages of the S2SMIX model. To quantify diversity within a set {ŷ m } M m=1 of translation candidates, we propose to evaluate average pairwise BLEU between pairs of sentences according to As an alternative metric of diversity within a set {ŷ m } M m=1 of translations, we propose a metric based on the fraction of the n-grams that are 3 https://sites.google.com/site/iwsltevaluation2016 4 https://github.com/moses-smt/mosesdecoder 5 https://nlp.stanford.edu/projects/nmt where ngram(y) returns the set of unique n-grams in a sequence y. We report average div_bleu and average div_ngram across the test set for the translation candidates found by beam search. We measure and report bigram diversity in the paper and report unigram diversity in the supplementary material.

S2SMIX configuration.
We start by investigating which of the ways of injecting the conditioning signal into the S2SMIX model is most effective. As seen in Section 3, the mixture components can be built by adding component-specific vectors to the logits (sf), the top LSTM layer (tp) or the bottom LSTM layer (bt) in the decoder, or all of them (all). Figure 2 shows the BLEU score of these variants on the translation tasks across four different language pairs. We observe that adding a component-specific vector to the recurrent cells in the bottom layer of the decoder is the most effective, and results in BLEU scores superior or onpar with the other variants across the four language pairs. Therefore, we use this model variant in all experiments for the rest of the paper.
Furthermore, Table 2 shows the number of parameters in each of the variants as well as the base SEQ2SEQ model. We confirm that the mixture model variants introduce negligible number   S2SMIX vs. SEQ2SEQ. We compare our mixture model against a vanilla SEQ2SEQ model both in terms of translation quality and diversity. To be fair, we compare models with the same number of beams during inference, e.g., we compare vanilla SEQ2SEQ using a beam size of 8 with S2SMIX-4 with 4 component and a beam size of 2 per component.
As an effective regularization strategy, we adopt label smoothing to strengthen generalisation performance (Szegedy et al., 2016;Pereyra et al., 2017;Edunov et al., 2018). Unlike conventional cross-entropy loss, where the probability mass for the ground truth word y is set to 1 and q(y ) = 0 for y = y, we smooth this distribution as: where is a smoothing parameter, and V is the vocabulary size. In our experiments, is set to 0.1.  Table 4: S2SMIX with 4 components vs SEQ2SEQ endowed with the beam-diverse decoder (Li et al., 2016) with the beam size of 4.
with the corresponding row in the bottom part for a fair comparison in terms of the effective beam size. Firstly, we observe that increasing the beam size deteriorates the BLEU score for the SEQ2SEQ model. Similar observations have been made in the previous work (Tu et al., 2017). This behavior is in contrast to our S2SMIX model where increasing the beam size improves the BLEU score, except en-es, which demonstrates a decreasing trend when beam size increases from 2 to 4. Secondly, our S2SMIX models outperform their SEQ2SEQ counterparts in all settings with the same number of bins. Figure 3 shows the diversity comparison between the S2SMIX model and the vanilla SEQ2SEQ model where the number of decoding beams is 8. The diversity metrics are bigram and BLEU diversity as defined earlier in the section. Our S2SMIX models significantly dominate the SEQ2SEQ model across language pairs in terms of the diversity metrics, while keeping the translation quality high (c.f. the BLEU scores in Table 3).
We further compare against the SEQ2SEQ model endowed with the beam-diverse decoder (Li et al., 2016). This decoder penalizes sibling hypotheses generated from the same parent in the search tree, according to their ranks in each decoding step. Hence, it tends to rank high those hypotheses from different parents, hence encouraging diversity in the beam. Table 4 presents the BLEU scores as well as the diversity measures. As seen, the mixture model significantly outperforms the SEQ2SEQ endowed with the beam-diverse decoder, in terms of the diversity in the generated translations. Furthermore, the mixture model achieves up to 1.7 BLEU score improvements across three language pairs.  Table 5: BLEU scores using greedy decoding and training time based on the original log-likelihood objective and online EM coupled with gradient estimation based on a single MC sample. The training time is reported by taking the average running time of one minibatch update across a full epoch.
Large mixture models. Memory limitations of the GPU may make it difficult to increase the number of mixture components beyond a certain amount. One approach is to decrease the number of sentence pairs in a minibatch, however, this results in a substantial increase in the training time.
Another approach is to resort to MC gradient estimation as discussed in Section 3.3. The top-part of Table 5 compares the models trained by online EM vs the original log-likelihood objective, in terms of the BLEU score and the training time. As seen, the BLEU score of the EMtrained models are on-par with those trained on the log-likelihood objective. However, online EM leads to up to 35% increase in the training time for S2SMIX-4 across four different language pairs, as we first need to do a forward pass on the minibatch in order to form the lower bound on the loglikelihood training objective.
The bottom-part of Table 5 shows the effect of online EM coupled with sampling only one mixture component to form a stochastic approximation to the log-likelihood lower bound. For each minibatch, we first run a forward pass to compute the probability of each mixture component given each sentence pair in the minibatch. We then sample a mixture component for each sentence-pair to form the approximation of the log-likelihood lower bound for the minibatch, which is then optimized using back-propagation. As we increase the number of mixture components from 4 to 8, we see about 0.7 BLEU score increase for en-de; while there is no significant change in the BLEU score for en-fr, en-vi and en-es.
Increasing the number of mixture components further to 16 does not produce gains on these datasets. Time-wise, training with online EM coupled with 1-candidate sampling should be significantly faster that the vanilla online EM and the original likelihood objective in principle, as we need to perform the backpropagation only for the selected mixture component (as opposed to all mixture components). Nonetheless, the additional computation due to increasing the number of mixtures from 4 to 8 is about 26%, which increases to about 55% when going from 8 to 16 mixture components.

Qualitative Analysis
Finally, we would like to demonstrate that our S2SMIX does indeed encourage diversity and im-
prove the translation quality. As shown in Table 6, compared with SEQ2SEQ, which mistranslates the second clause, our S2SMIX is not only capable of generating a group of correct translation, but also emitting synonyms for different mixture components. We provide more examples in the supplementary material.

Related Work
Obviously, different domains aim at different readers, thus they exhibit distinctive genres compared to other domains. A well-tuned MT system cannot directly apply to new domains; otherwise, translation quality will degrade. Based on this factor, out-domain adaptation has been widely studied for MT, ranging from data selection (Li et al., 2010;Wang et al., 2017), tuning (Luong and Manning, 2015;Farajian et al., 2017) to domain tags (Chu et al., 2017). Similarly, in-domain adaptation is also a compelling direction. Normally, to train an universal MT system, the training data consist of gigantic corpora covering numerous and various domains.This training data is naturally so diverse that Mima et al. (1997) incorporated extralinguistic information to enhance translation quality. Michel and Neubig (2018) argue even without explicit signals (gender, politeness etc.), they can handle domain-specific information via annotation of speakers, and easily gain quality improvement from a larger number of domains. Our approach is considerably different from the previous work. We remove any extra annotation, and treat domain-related information as latent variables, which are learned from corpus.
Prior to our work, diverse generation has been studied in image captioning, as some of the training set are comprised of images paired with multiple reference captions. Some work puts their efforts on decoding stages, and form a group of beam search to encourage diversity (Vijayakumar et al., 2016), while others pay more attention to adversarial training (Shetty et al., 2017;. Within translation, our method is similar to Schulz et al. (2018b), where they propose a MT system armed with variational inference to account for translation variations. Like us, their diversified generation is driven by latent variables. Albeit the simplicity of our model, it is effective and able to accommodate variation or diversity. Meanwhile, we propose several diversity metrics to perform quantitative analysis.
Finally, Yang et al. (2018) proposes a mixture of softmaxes to enhance the expressiveness of language model, which demonstrate the effectiveness of our S2SMIX model under the matrix factorization framework.

Conclusions and Future Work
In this paper, we propose a sequence to sequence mixture (S2SMIX) model to improve translation diversity within neural machine translation via incorporating a set of discrete latent variables. We propose a model architecture that requires negligible additional parameters and no extra computation at inference time. In order to address prohibitive memory requirement associated with large mixture models, we augment the training procedure by computing the posterior distribution fol-lowed by Monte Carlo sampling to estimate the gradients. We observe significant gains both in terms of BLEU scores and translation diversity with a mixture of 4 components. In the future, we intend to replace the uniform mixing coefficients with learnable parameters, since different components should not necessarily make an equal contribution to a given sentence pair. Moreover, we will consider applying our S2SMIX model to other NLP problems in which diversity plays an important role.


extfieldsofstudy:
Computer Science
Mathematics

text:
Optimization of the n-dimensional sliding window inter-channel correlation algorithm for multi-core architecture

Calculating the correlation in a sliding window is a common method of statistical evaluation of the interconnect between two sets of data. And although the calculation of a single correlation coefficient is not resource-intensive and algorithmically complex, sequential computation in a large number of windows on large data sets can take quite a long time. In this case, each value in the data, falling into different windows, will be processed many times, increasing the complexity of the algorithm and the processing time. We took this fact into account and optimized the correlation calculation in the sliding window, reducing the number of operations in the overlapping area of the windows. In addition, we developed a parallel version of the optimized algorithm for the GPU architecture. Experimental studies have shown that for a 7x7 correlation window sliding in one pixel increments, we were able to accelerate the processing of an 12 MPixel image pixels on the GPU by about 60 times compared to the serial version running on the CPU. The article presents an optimized version of the algorithm, a scheme for its parallelization, as well as the results of experimental studies.

Introduction
Correlation in a sliding (or rolling) window is a very common method of statistical evaluation of synchronisation of two data sets used in image processing [1], signal analysis [2], medicine [3], geophysics [4], etc. The functions for calculating parameters in a sliding window (sum, arithmetic mean, correlation, etc.) are included in many software packages (for example, MATLAB, R, S-PLUS) and are actively used in the analysis of time series in financial statistics [5].
In this paper, we consider the Pearson linear correlation calculated in a sliding window between sets of spatial data on a regular grid, preferably between two images. The applied problem concerns the automatic search of fishing boat lights on multispectral satellite images obtained at night time [6]. The task is to filter the false detections caused by the reflection of moonlight from the clouds by calculating the correlation in the sliding window between the data in the visible and infrared channels. In the presence of clouds, the correlation coefficient should be close to -1, since the values in the visible channel increase due to the reflection of moonlight from the clouds, while the values in the infrared channel decrease due to the lower cloud temperature as compared to the ground surface.
The basic algorithm can be represented as a synchronous passage of two sets of data by a sliding window of a given size with a given step, and calculating the Pearson linear correlation between the values of the investigated sets in this window. In this case, the values in n-dimensional sliding windows are drawn into one-dimensional vectors.
To calculate the correlation in each window, the classical Pearson linear correlation formula is used: where c is Pearson linear correlation between the x and y datasets; x i and y i are the individual sample values; and are the arithmetic means for the x and y datasets; x y n is the sample size.
Calculating the correlation in a sliding window using the classical formula leads to significant overhead, as illustrated by the example of processing 2D images, shown in Figure 1. Windows 1 and 2 have a large overlap. When calculating the correlation in each window, the elements that belong to both windows will be processed twice. If the step of the sliding window is small enough, for example, one pixel, and the window size itself is large, then each element of the image will be processed as many times as the windows cover it.
In the work presented, we optimized the basic version of the algorithm for interchannel correlation in a sliding window, taking into account the modern trend to multi-core architecture, and have tested the optimized version on graphic processes.

Optimized version of the algorithm
To reduce the overhead caused by multiple processing of the same element, we used a modified correlation calculation formula that is often used to reduce the number of computational operations associated with subtracting from each arithmetic mean value: An advantage of this formula for the problem studied in this paper is the expression of the correlation value through the sums of the quantities x , y , xy , x 2 , y 2 . Such sums for one-dimensional series can be obtained by a linear algorithm for calculating the sum in the sliding window of the movements (similar to the moving average of the rollmean R function in the zoo library ). To move from a one-dimensional series to a two-dimensional image matrix, we 1 consistently apply the "sliding summation" first in one dimension, and then in the second dimension ( Figure 2). As a result of the superposition of two one-dimensional summations with a window of length k, we obtain a matrix of the same size as the original, but at each position of which there will be a sum of elements of this matrix taken in a window of size k * k centered at a given point (provided that the one-dimensional summation also puts the result in the center of the window).
Thus, the optimized algorithm is as follows: 1. Element multiplication of matrices x and y computes the values of matrices xy, x 2 , y 2 . 2. For each matrix x, y, xy, x 2 , y 2 , a summation is made with a sliding window of length k over all rows (the results are denoted by Ax, Ay, Axy, Ax 2 , Ay 2 ). 3. For each matrix Ax, Ay, Axy, Ax 2 , Ay 2 , a summation is made with a sliding window of length k over all columns; the results are denoted by Sx, Sy, Sxy, Sx 2 , Sy 2 (the sum matrix in the windows centered on this element). 4. For each element of the resulting matrix, we calculate the correlation by the formula 1.
In order to generalize the above algorithm to the n-dimensional case, it is necessary to repeat the summation with the sliding window (steps 2-3) n times, one for each measurement.
It is worth noting that optimization can be done not only with the help of a sliding window, but also with the help of cumulative sums [7]. We have tested this approach, but in our application the cumulative sums give an increase in the execution time in comparison with the algorithm proposed in this section. Theoretical calculations (which we do not show here to save space) confirm the increase in the number of operations using cumulative sums compared to the summation in a sliding window.

Parallelizing an optimized version of the algorithm
The algorithm described in the previous section can be efficiently parallelized for the case of a multi-core architecture, for example, the GPU: 1. In step 1, each thread processes elements from a small neighborhood of the source data. 2. In step 2, each thread processes one or more rows, depending on the size of the data and the degree of parallelism. 3. In step 3, each thread processes one or more columns, depending on the size of the data and the degree of parallelism. 4. In step 4, each thread processes elements from a small neighborhood of the source data. 5. For the n-dimensional case, parallelization is performed analogously.

Implementation of algorithms and testing
The algorithms described in the previous sections were implemented in software. The basic and optimized versions were implemented in the programming language MATLAB. The parallel optimized version is implemented in C ++ programming language for CUDA libraries. In the software implementation, the correlation in the sliding window between two two-dimensional data sets is calculated.
Step of the window is one cell. Both sets of data must be the same size.
The result set is the same size as the arguments. For each pair of windows, the correlation result is written to the middle cell. Data may contain missing values that are identified by a number less than or equal to -999. If at least one missing value occurs in the correlated windows, then the resulting value for the two given windows is set to -2. Values along the edges of the resulting matrix, for which the correlation windows corresponding to them exceed the boundaries of the matrix-arguments, are set equal to -2.
Testing and comparison of software implementations was carried out on two satellite images, reprojected to a regular geographic grid. In particular, we used images in the visible and infrared channels obtained by the VIIRS sensor installed on the Suomi NPP satellite. The size of the images is 4064x3072 pixels. The size of the correlation window is 7x7 pixels.
The MATLAB version was tested on a workstation with the following characteristics: • Processor: 2 x Intel Xeon E5640 (8 cores, 2.66 GHz) • RAM: 8 GB • Operating system: Windows 7 The GPU version was tested on a computer with the following characteristics: Obviously, the overall performance of a computer with a GPU is significantly higher than that of a computer with a CPU, so it would seem that the comparison is incorrect. But the essence of the experiment is not to compare the specific performance of two different architectures, but to evaluate the effect and quality threading of a single-threaded algorithm, as well as the efficiency of running a parallel algorithm on a GPU architecture (the latter is not equivalent to the first one, since the algorithm can be paralleled well, but its implementation on the GPU architecture may not display a large increase in performance due to various technical limitations, for example, the cost of sending data can significantly exceed the cost of processing them.) If parallelization is effective, then the loading of the GPU in our test will be high, and the processing speed should be significantly reduced in comparison with the CPU, so as the overall performance of the GPU is higher. If parallelization is inefficient, then either we start with a large overhead for organizing parallelization on the GPU, or only a small fraction of the GPU will be downloaded der (and we know that one GPU core has suschestenno lower performance than a single CPU core so the GPU efficiency is achieved only through massive parallelism). In both cases, we will not get a significant performance gain on any GPU-device, and may even get the opposite effect -the slowdown of the algorithm.
The testing showed the following results: • Basic version: about 42 seconds, about 1 second -data reading and about 41 seconds -the work of the correlation algorithm itself. Of course, the above calculations can not accurately reflect the time relationship, since not all operations are performed the same time, and also because we did not take into account computer operations, for example, assigning values, but nevertheless these formulas can be used to approximate the complexity of algorithms. If we substitute in the above formulas the values used in the experiment, we get: • for the basic version of the algorithm: 10 * n 2 * A * B = {n = 7} = 490 * A * B; • for the optimized version of the algorithm: 36 * A * B. And their ratio is 490/36 ≈ 13.5, which is pretty close to the ratio of the execution time of these algorithms, obtained during testing.
Analyzing the received formulas it is possible to draw a conclusion: the number of arithmetic operations in the basic version of the algorithm depends squarely on the window size, whereas in the optimized algorithm the amount of arithmetic operations does not depend on the window size at all. Similar theoretical calculations can be made for a parallel version, but they will differ greatly from the practical implementation for the following reasons: 1. We can not in the theoretical calculations take into account the overhead associated with the initialization and support of one computational flow. Moreover, different computing architectures will require different overhead. 2. It is difficult for us to take into account the degree of parallelism in a particular architecture and the particular implementation. For example, not all parallel threads can run in parallel. For example, in the case of a GPU, a particular partition of the data space into blocks can greatly affect performance. 3. It is difficult to take into account the costs of data exchange. For example, again, in the case of a GPU, the compiler can optimize the program by placing data on registers or in local memory. And this leads to an increase or decrease in the operating time of the program.
The most critical issue of accuracy arises when all values of the correlation window are equal, in which case the denominator in the correlation calculation formula goes to zero, but even a small calculation error can lead to deviation of the calculated value from zero. As a result, the divide-by-zero error will not arise, and the resulting value will be incorrect. To avoid this situation, we propose to check all the elements for equality before carrying out the basic calculations.
Nevertheless, if the extreme values are not taken into account, then the error of the algorithms is rather small: if the results of the base version are taken as a reference result, then in our test for the optimized version the maximum error is less than 10 -3 .

Conclusion
The optimized version of the algorithm for calculating inter-channel correlation in a sliding is used in VIIRS boat detection algorithm [6] to filter false detections of fishing boat lights caused by the reflection of moonlight from clouds. Parallel implementation have reduced computational time by more than 15 times for a multicore CPU architecture, and 2 times more for a CUDA GPU. Initialization time of the GPU is not taken into account, since initialization can be done in advance, when the program starts.