A Feature Clustering Approach Based on Histogram of Oriented Optical Flow and Superpixels

Visual feature clustering is one of the cost-effective approaches to segment objects in videos. However, the assumptions made for developing the existing algorithms prevent them from being used in situations like segmenting an unknown number of static and moving objects under heavy camera movements. This paper addresses the problem by introducing a clustering approach based on superpixels and short-term Histogram of Oriented Optical Flow (HOOF). Salient Dither Pattern Feature (SDPF) is used as the visual feature to track the flow and Simple Linear Iterative Clustering (SLIC) is used for obtaining the superpixels. This new clustering approach is based on merging superpixels by comparing short term local HOOF and a color cue to form high-level semantic segments. The new approach was compared with one of the latest feature clustering approaches based on K-Means in eight-dimensional space and the results revealed that the new approach is better by means of consistency, completeness, and spatial accuracy. Further, the new approach completely solved the problem of not knowing the number of objects in a scene.

I. INTRODUCTION
Segmentation of semantic objects is a crucial first step in video understanding applications such as annotating and indexing of videos for video search and retrieval purposes. High-level object recognition typically builds upon this step. As such, finding robust solutions for segmenting of objects in videos has become an active research problem. Standard approaches to the problem include motion structure analysis and appearance-based methods. Several recent studies have shown a significant improvement on using a flow of visual features for object segmentation [1]- [4].
Earlier studies such as in [5] showed that motion and 3D world structure can be used to segmentation of objects in videos. The semantic segmentations solely based on motionderived 3D structures tend to neglect the objects with smaller regions. The study in [5] has shown that appearance-based methods can be used to detect some of the objects which are neglected by the motion and structure cues. Recent studies [1], [2] have showed that long-term flow analysis can help to accurately refine object regions. However, practically some of the objects may appear or disappear within a small fraction of the video due to the camera motion hence long term analysis may not feasible in some cases. The fast segmentation approach in [3] works reliably when satisfying certain conditions such as the object should move differently from its surrounding background in a good fraction of the video and at least a part of the object should move somewhere in the video. Such assumptions prevent the use of the algorithm for the majority of the videos due to uncertain behaviors of the salient objects. A previous study [4] has showed that a simple clustering technique can be used to cluster salient dither pattern (SDPF) features [6] by using temporal information and dither patterns, for semantic segmentation. However, it requires a human interaction for estimating the number of clusters (K) by inspecting a short sample of the video. Further, the algorithm tends to neglect some objects due to the inconsistent estimation error of K along the video. Further, the K should be re-estimated in cases of applying the algorithm to totally new set of videos.
This paper presents an improved unsupervised visual feature clustering approach for segmenting static and moving objects in videos with camera motion, by using short-term flow information and an appearance-based cue. Particularly the approach refines a set of estimated segments to form highlevel object segments by using color based cue and local motion information of each of the estimated segments. It has been shown that superpixel can be used for estimating video segments as an initial step of segmenting high-level objects [3], [7]. The comparative study mentioned in [8] has revealed that the Simple Linear Iterative Clustering (SLIC) superpixel has the state-of-the-art performance by means of time efficiency, memory efficiency and segmentation reliability. In this new approach, SLIC superpixel is used to estimate low-level segments in a weighted Lab color space. SDPF is used as the point feature due to its invariant properties [6] and the optical flow is calculated using the pyramidal implementation of Lucas-Kanade optical flow [9]. The segment refining is done by using a short term histogram of oriented optical flow (HOOF) with Radial Basis Function (RBF) kernel on Bhattacharya distance [10] and mean colors of superpixels. A refining process is employed as a bottom-up approach by considering only the neighborhood relationships, starting from the superpixel extraction to the clustering of features correspond to the objects. The superiority of the new approach is demonstrated by comparing with the results in [4].
The rest of this paper is organized as follows. In section II, related studies are briefly discussed with respect to the feature clustering. Section III describes the new approach to feature clustering. Section IV describes the experimental setup for comparing the new approach with the previous works. Section V, the experimental results are discussed with followed by conclusion Section VI.

II. RELATED STUDIES
In this section SLIC superpixels, SDPF feature and HOOF are described with respect to feature clustering.

A. SLIC Superpixel
Superpixels capture redundancy in an image and greatly reduce the complexity of subsequent image processing tasks. They have proved increasingly useful for semantic segmentation [3]. SLIC superpixels are obtained by performing a local clustering of pixels in the 5-D space defined by the L, a, b values of the CIELAB color space and x,y pixel coordinates. In order to self-contain the paper, SLIC algorithm will be explained briefly here onward.
At the beginning of the algorithm, it places a given number of approximated superpixels at nearly equal grid intervals but with avoiding pixels that have higher gradient such that the pixels which are not located on edges. This step ensures the approximation of superpixels with minimum noise. Then a localized k-means clustering is employed by using the summation of the Euclidean distance of colors in CIELAB color space and spatial distance which is weighted by a spatial proximity emphasizing factor [8]. SLIC shows a strong perceptual homogeneity within the superpixels and O(n) time complexity hence it is used in this new approach.

B. SDPF Points
This study mainly focused on clustering SDPF feature points with respect to high-level semantic object segments. SDPF is a point like low dimensional feature which is based on the concept in Compacted Dither Pattern Code (CDPC), which the dithering reduces both color and spatial data density while preserving the visual impression [11]. The adaption of the concept in CDPC by including both color and spatial information ultimately creates a point-like feature that extracted based on its saliency over a set of local dither patterns. The detailed steps of extracting SDPF features and its strong invariant properties are described in [4], [6]. The ultimate goal of the feature clustering is to recognize multiple objects by using the low dimensional SDPF descriptor hence SDPF points is used to initiate the optical flow estimation. Besides the main goal of this work, the geometrical invariant properties of SDPF is used to overcome the heavy geometrical transformation of object segments caused by the ego-motion.

C. HOOF
The concept of HOOF feature was introduced in [10] for human activity recognition. The results of the same study have shown a time series of HOOF can characterize a high-level object motion. HOOF can be built by binning each flow vector according to its primary flow angle and weighted according to its magnitude. Practically the flow estimation of feature points is very noisy especially when the magnitude of flow is very low. Theoretically HOOF feature is robust to the flow noises due to the flow angles are weighted by the flow magnitudes. As described in [10], HOOF does not reside in Euclidean space hence similarity matching should be done by using a proper matching function from the space of histogram to highdimensional Euclidean space.
Although this study is not related to activity recognition, the ability of HOOF to characterize specific motion is taken as an advantage in this approach. This approach argues that neighbors of low-level segments that belong to the same semantic object have similar HOOF features over a short period of time.

III. PROPOSED CLUSTERING APPROACH
As the first step of the clustering, SDPF points are extracted from the first frame. Then calculate the flow using the pyramidal implementation of Lucas-Kanade optical flow estimation approach with taking the SDPF points as the references for the tracking. The optical flow is estimated for a predefined number of frames (T f ) while recording the flow angle and magnitude of each tracking point for each frame. The T f was set to 3 in the experiment, in order to assess the performance of the clustering approach with a shallower history of motion.
Once the flow estimation process reaches to the T f , the clustering process is initiated by using SLIC superpixel algorithm. The original SLIC algorithm requires the input image in CIELAB color space due to two colors in CIELAB space with smaller Euclidean distance are perceptually similar. The L component in CIELAB color space closely matches the human perception of brightness [12]. Therefore, the component L is very sensitive to shadows and different amount of light intensities over the video frame. The impact of this nature of L component causes the SLIC algorithm to adhere some pixels as an edge although the pixels are actually placed on boundaries of shadows but not actual object boundaries. Therefore the weighted Lab color such that L'ab is obtained using (1).
The value of m is kept at 0.5 in the experiment because the increment of the value of α caused the SLIC to detect many false boundaries.
SLIC algorithm further requires two other parameters, the number of clusters (K) and minimum distance for color proximity (N c ). In all the experiments conducted in this study, K and N c were set to 50 and 10 respectively, based on the frame size of the videos in the dataset. A color cue is used in addition to the motion cue in a later step of the process of refining the superpixels. This color cue is obtained by averaging the colors of pixels in each superpixel segments. In the refining process, this average color similarity between two superpixel regions is used as one of the criteria for the refining. The similarity of the colors is obtained by using the Euclidean distance as defined in (2).
Where L', a and b are the weighted Lab components of average color calculated for a superpixel. The next step is to obtain HOOF for the flow of SDPF points in each superpixel. Therefore, a time series of HOOF features is calculated for each superpixel using the feature points located within the region of the superpixel. A flow vector, v= x,y T with the angle θ= tan -1 y x is in the range -π<θ≤π. By mapping the range of the θ to 0<θ≤2π the corresponding bin b can be calculated by (3).
where the B is the number of bins. In the experiments, the value of B is set to 30 as it did on its original work [10]. The contribution of the vector to the HOOF is calculated by (4).
where the is the value of b th bin in the histogram h. The HOOF for each superpixel is calculated by considering the local flow vectors correspond to the each of the superpixels. The time series of the HOOF can be calculated by using the recorded flow data within the feature points. At the end of the calculation each of the superpixel will have a set of HOOF features which's the size is equal to the predefined length of time series. As described in [10] HOOF cannot be treated as a simple vector in Euclidean space hence Euclidean-based distance measurements cannot be applied to compare two HOOF features. Hence, the similarity of histograms is measured by calculating a kernel based on Riemannian metric on a hypersphere which is created by the square root representation of the histograms. The kernel is defined as (5).
where h 1 and h 2 are the two histograms. In this study, the above kernel was implemented as mentioned in [10] using RBF kernel on Bhattachcharya distance between histograms, as shown in (6) Since it is needed to compare the similarity between two time series of HOOF, first, the RBF kernel is calculated for each pair of histograms created in the same instance of time from neighboring superpixels. Then the magnitude of the set of RBF kernel values is calculated to measure the similarity between whole time series of HOOF by (7).
Where D h is the similarity measure between the two time series of HOOF and k RBF t=1 is the RBF kernel of two HOOFs calculated for two superpixels at the time step 1 and so on.

A. Refining Process of Superpixels
The regions of high-level semantic objects are formed by refining the superpixels. The refining process is involved in following 6 steps. and T c are predefined thresholds. 3. In case the merging criteria are satisfied, check whether, either of the current superpixel or the neighbor superpixel has been given a cluster number. 4. If either of superpixels has a cluster number, then the cluster number is shared by both superpixels and merge the time series of HOOF and average colors of the two superpixels. 5. If neither of superpixels has a cluster number, then create a new cluster number by incrementing the number of existing clusters, assign to both superpixels and merge the time series of HOOF and average colors of the two superpixels. 6. Then go to the next superpixel, calculate D h and D c with all its 8 neighbors and repeat from the step 2, until there is no superpixel remains without a cluster number. The major problem of merging superpixels is to select merging criteria which accurately quantify the similarity between two superpixels, which represents the likelihood of being the two superpixels in a single semantic object. Due to the high color complexity of typical objects, the merging criteria should be not limited to appearance based cue. Therefore in this approach, the merging criteria have been designed to quantify the similarity by at least one of the color or motion cues at any given instance.
In order to find the two thresholds T h and T c , initially the maximum possible value of T h such that 3 was assigned to T h and checked the correctness of cluster merging manually, by varying T c . The observations suggested that T c can be varied effectively within a very narrow range such that [0, 30] due to the weighting process of CIELAB colors. In order to balance the color proximity measurement, T c was set to 15 and check the merging correctness manually, by varying T h . The results suggested that 1 and 15 are good enough for T h and T c respectively for most of the time. Hence, the same values were used for all the experiments conducted in this study.

IV. EXPERIMENTAL SETUP
The performance of the new approach is evaluated using the CamVid dataset [13]. The dataset contains high-resolution video with a frame rate of 30Hz, captured from the perspective of a driving automobile. The video contains both the moving and static objects. The performance evaluation requires measuring the correspondence between th semantic objects in a scene hence spatial ce the clusters were used to find the correspond of a cluster is laid within the boundary of an considered that the particular cluster belo object. The assessment was done for 100 frames comparing with the ground truth vide

A. Evaluation Criteria
One of the objectives of this stud improvement of the new clustering approach previous clustering approach described in [ same evaluation criteria are used as in the order to self-contain the paper, the performa used in both of the studies are briefly descr An objective comparison for each frame spatial accuracy measure S er [14], as defined V. EXPERIMENTAL SETUP The performance evaluation results of t approach are presented comparing with the the previous study, the completeness measur accuracy of three clustering approach in 8-d (flow angle, flow magnitude, four colors in patterns and x,y coordinates) using K-mean Expectation Maximization (EM-8D) and d (DBScan) were evaluated. In this paper, analysis of the worst approach in the prev used the DBScan is omitted in order visualization of the results of the new appro the new approach which is based on Ref using HOOF and Color cue will be referred shows the obtained results for RSHC with Means-8D and EM-8D. Fig. 1 demonstrates that the new approa small but significant improvement over the p Precisely RSHC achieved over 0.87 of aver for the whole video where as K-Means achieved around 0.86 and 0.81 of avera respectively. Although the improveme completeness is not very significant, Fig. 1 e achieved a significant improvement by mea over the frames. The consistency of the K-M 8D are lower because the less robustness he clusters to the entroid of each of dence. If a centroid n object, then it is ongs to the same visually different eo.
dy to assess the h compared to the [4]. Therefore, the previous work. In ance measurement ribed here onward.
is done using a in (8).
into the n th object, ly clustered to the umber of semantic clusive disjunction pl is defined in (9) the two sets.
the new clustering e results in [4]. In rement and spatial dimensional space n the salient dither ns (K-Means-8D), density based scan the performance ious study, which to simplify the oach. Here onward fining Superpixels d as RSHC. Fig. 1 comparing the Kach has obtained a previous methods. instance of clustering with K-Means-8D using the estimated number of clusters based on the number of feature points, and a sample instance of clustering with RSHC are shown in Fig. 2 in order to demonstrate the quality of the performance of the new approach. The Fig. 2 clearly shows that RSHC is less affected by the lack of information i.e. number of objects in the scene whereas it tends to heavy over-clustering in the K-Means based approach.
Besides the completeness and consistency, RSHC also shows a significant improvement on spatial accuracy error compared to the other two methods. The average spatial accuracy error of RSHC method is around 0.27 where the K-Means-8D and EM-8D are 0.31 and 0.42 respectively. We believe that the improvement in spatial accuracy error is gained by the ability of accurate adherence to object boundaries by SLIC algorithm with the weighted CIELAB colors.

VI. CONCLUSION
This paper presents a new visual feature clustering approach on clustering SDPF feature points to semantic objects by refining SLIC superpixels using both short-term time series of HOOF and color similarity measurements. The study has proved that the proposed scheme works consistently over different visual scenes in a video and can perform better than the methods which apply global clustering with indirect estimation of the amount of clusters. The performance evident that time series of HOOF can be used to characterize both static and moving objects. However, since HOOF does not include the spatial distribution of flow angle and magnitude within a superpixel, the issue how to incorporate motion boundaries for further improving the spatial accuracy of clustering, remains. The findings of this study will be used for detecting and recognizing multiple semantic objects using SDPF feature descriptor.


extfieldsofstudy:
Physics
Computer Science

text:
Algorithms for synthesizing management solutions based on OLAP-technologies

OLAP technologies are a convenient means of analyzing large amounts of information. An attempt was made in their work to improve the synthesis of optimal management decisions. The developed algorithms allow forecasting the needs and accepted management decisions on the main types of the enterprise resources. Their advantage is the efficiency, based on the simplicity of quadratic functions and differential equations of only the first order. At the same time, the optimal redistribution of resources between different types of products from the assortment of the enterprise is carried out, and the optimal allocation of allocated resources in time. The proposed solutions can be placed on additional specially entered coordinates of the hypercube representing the data warehouse.

Introduction
OLAP-technologies are a convenient tool for analyzing large amounts of information. Convenience and speed of analysis are due to the structuring of data based on the hypercube. And the analysis itself is carried out naturally with the aim of synthesizing managerial decisions. However, algorithms for synthesizing such solutions are often absent in OLAP-technologies. At the same time it is logical, having large volumes of actual data, to synthesize simultaneously alternative managerial decisions and offer them for final choice to the decision-maker.

Theory
It is impossible to accept any managerial decision without having the information necessary for this, usually quantitative. This requires the creation of data warehouses. They are assigned to the process of collecting, screening and preliminary processing of data in order to provide the resulting information to users for statistical analysis, and often creating analytical reports [1]. Ralph Kimball, one of the authors of the concept of data warehousing, described the data warehouse as "a place where people can access their data" [2,3,4].
The concept of OLAP was described in 1993 by Edgar Codd, a well-known database researcher and author of the relational data model [5]. In 1995, based on the requirements set by Codd, the socalled Fast Analysis of Shared Multidimensional Information (FASMI) was formulated. The relevance of using OLAP systems is to provide the ability to work with data in terms of the subject area without the knowledge of the information storage architecture [6]. The main goal of the development of OLAP technologies is to provide the user with structured information tools to explain current data and build forecasts [7].
In various OLAP products, there are two main options for organizing data [8] in multidimensional cubes: a hypercubic model and a polycubic model. Such representation would be useful in model problems [9]. In this case, the functions of sampling, processing and presentation of data can be assigned to a special OLAP machine [10].
Hypercube is initially multidimensional and its dimension is not limited in any way. For example, by postponing the volumes of products produced by its coordinates, it is possible to obtain assortment space [11]. The cube as a matter of fact represents a "photos" set of system's indicators of the enterprise condition consistently in time. This allows one to apply the step-by-step method of dynamic programming Bellman to optimize the redistribution of enterprise resources, and then, the proposed solutions can be placed on additional specially entered hypercube coordinates.
To formulate the problem of optimal redistribution of resources, it is necessary to know the efficiency functions of their investments in the production of a particular product from the assortment. In the first approximation, one can confine oneselves to a quadratic function: (1) Here is the value of the resource directed to the production of the i-th product, of τ, v are the constant coefficients chosen for the approximation.
Let us compose the Bellman function [12]: (2) Where N is the total number of products, i is the serial number of the product type, is the cumulative variables, M is the large number in the penalty function, and Z is the total resource of the enterprise: (3) After solving the problem that determines the allocated resources, it is logical to develop an algorithm for supplying them to control objects, for example, determined by the Euler-Lagrange method.
Functionals of optimal inventory management are set as follows: , The results of solving this problem can also be represented on additional specially introduced hypercube coordinates

Model and Method
Most often, production tasks are associated with the replenishment of stocks of various resources. This process can be described by a simple equation: (5) where is the expense of the i-th resource.
Then the Lagrangian can be written in the form: , (6) where is the i-th Lagrange multiplier.
To find the optimal control, it is necessary to compose the Euler equations: Substituting the Lagrangian, let us obtain the system of equations: Transforming the first two equations, let us obtain: (9) Substituting the expression of the variables from the first two equations in the third equation (8), let us obtain: .
(10) The solution of the corresponding homogeneous equation is expressed in terms of the hyperbolic functions: . (11) A particular solution must be sought in the form of the right-hand side of equation (10): .
(13) The solution of this equation, depending on the sign of the fraction [13]: (14) Here Now the algorithm of actions is clearthere is a need to analyze the data: if function is approximated by a harmonic, then the third solution fits us, if the monotone exponent, then the first one is suitable, in all other cases let us choose the second solution and so all the constants (A, , ) approximated function .
Substituting the solution obtained into the general solution and into the first equation in system (9), let us obtain: . (16) Using the initial values of the resource reserve: (17) and the condition of limited resources ( , let us take from equation (3)), one obtains: . (18)

Data and methods
Data for the test case were taken from the milk processing plant [14]. The products of this enterprise have been in demand for a long time among residents of Orenburg and the region. Natural ingredients and high quality of production are the distinctive features of all dairy products of the creamery. In addition, this production is distinguished by a large assortment -up to 10 positions and short shelf life of products from 3 to 10 days. This circumstance exacerbates the requirements for the speed and accuracy of planning production volumes.

The results and discussion
The developed algorithm includes the optimal distribution of the enterprise's production resources, primarily financial by product, approximating the resource consumption for the previous period and forecasting this expenditure for the next period in accordance with the derived formulas (15,16). In Figure 1, a sales histogram is shown, and in Figure 2    The type of the predictive function is easily chosen by minimizing the error in the approximation of real data.

Conclusion
Thus, the developed algorithms for synthesizing management decisions based on OLAP-technologies are a convenient tool of the system analyst, which allows forecasting the needs and accepted management decisions on the main types of the enterprise resources. Their advantage is the efficiency, based on the simplicity of quadratic functions and differential equations of only the first order.


extfieldsofstudy:
Computer Science
Mathematics

text:
A robust parallel algorithm for combinatorial compressed sensing

In previous work two of the authors have shown that a vector $x \in \mathbb{R}^n$ with at most $k<n$ nonzeros can be recovered from an expander sketch $Ax$ in $\mathcal{O}(\mathrm{nnz}(A)\log k)$ operations via the Parallel-$\ell_0$ decoding algorithm, where $\mathrm{nnz}(A)$ denotes the number of nonzero entries in $A \in \mathbb{R}^{m \times n}$. In this paper we present the Robust-$\ell_0$ decoding algorithm, which robustifies Parallel-$\ell_0$ when the sketch $Ax$ is corrupted by additive noise. This robustness is achieved by approximating the asymptotic posterior distribution of values in the sketch given its corrupted measurements. We provide analytic expressions that approximate these posteriors under the assumptions that the nonzero entries in the signal and the noise are drawn from continuous distributions. Numerical experiments presented show that Robust-$\ell_0$ is superior to existing greedy and combinatorial compressed sensing algorithms in the presence of small to moderate signal-to-noise ratios in the setting of Gaussian signals and Gaussian additive noise.

I. INTRODUCTION
C OMPRESSED sensing is a well studied method by which a sparse or compressible vector can be acquired by a number of measurements proportional to the number of its dominant entries [2], [3]. To fix notation, let χ n k be the set of vectors in R n that have at most k nonzero entries, let x ∈ χ n k and let A ∈ R m×n be a matrix with m < n. We will refer to A as the measurement matrix, x as the signal and y = Ax as the measurements. The goal of compressed sensing is to recover the sparsest, most parsimonious, x ∈ R n from the measurements y and the matrix A. Letting · 0 denote the number of non-zeros in x, the problem of finding x can be written as min x∈R n x 0 s.t. Ax = y.
Many algorithms have been developed to solve this problem or equivalent formulations and there are good theoretical results on when and how fast recovery of a signal is possible given certain types of measurement matrix A and signal x. These algorithms can be broadly categorized into convex optimization based algorithms like those implemented in [4]- [7] and greedy algorithms [8]- [14], and were designed and analysed for the setting of dense sensing matrices; e.g. independent RMS acknowledges the support of CONACyT This work was supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1 Furthermore, it is based on work partially supported by the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with PA Consulting Group (sub-)Gaussian entries or randomly subsampled Fourier matrices. For a more detailed introduction to compressed sensing see [15].
Here we extend an algorithm proposed in [1] which can be used to recover exactly a sparse signal from its expander sketch (see Section II for details). Specifically, [1] proposed Parallel-0 (Algorithm 2), for noiseless combinatorial compressed sensing which is guaranteed to converge in O(nnz(A) log(k)) where the sensing matrix A is an expander matrix (Definition I.1) and the signal x ∈ χ n k is dissociated in the sense of Definition I.2 or the signal is drawn independently of A. For alternative combinatorial compressed sensing algorithms see, for example, [16]- [21]. We borrow notation from combinatorics and use the shorthands  (≤k) . We denote by E m×n k,ε,d the set of (k, ε, d)expander matrices of dimension m × n.
Definition I.2 (Dissociated signals [1]). A signal x ∈ R n is said to be dissociated if j∈S1 An example of (almost surely) dissociated signals are those drawn from a continuous distribution. It is shown in [1] that if y = Ax is an expander sketch and x ∈ χ n k is dissociated, then there exists a subset T ⊂ [n] such that, for each j ∈ T , |{i ∈ [m] : y i = x j }| is bounded below by a positive constant depending on d and ε. This guarantees that if |{i ∈ supp(a j ) : y i = y }| > d/2 then x j = y . Parallel-0 (Algorithm 2) implements this observation by lettingx = 0 and estimating the decrease in y 0 when performing the updatex j ←x j + y . We denote for j ∈ [n] its neighbours by N (j) := {i ∈ [m] : |A ij | > 0}; to estimate the decrease in y 0 , Parallel-0 computes n e ← |{ ∈ N (j) : We extend their approach to the additive noise signal model of y = y + η with y = Ax and η ∈ R m by replacing (1)-(2) with scores estimating the distribution of n e and n z , e.g. (8)- (9). That is, we follow a Bayesian approach to the computation of these scores and estimate: 1) The probability of y i = 0 given that we observeŷ i .
Among our contributions are series approximations for (3)-(4) when the signals and measurements are generated according to the generating model given in Definition I.3 and illustrated in Figure 1. In what follows, we let D(R) be the set of distributions supported on R. If µ ∈ D(R), we write z ∼ µ to denote that z was drawn according to the distribution µ. We also use the notation v i i.i.d.
∼ µ to denote that each v i is drawn independently at random from µ. Finally, we use U (S) to denote the uniform distribution over a set S. Definition I.3 (Generating model GM(n, m, k, d, µ, ν)). Let n, m, k, d ∈ N be such that k < m < n and d m. Let µ, ν ∈ D(R). Then, the problem (A,ŷ) is drawn from the model GM(n, m, k, d, µ, ν) if A ∈ {0, 1} m×n andŷ ∈ R m are such that 1) each column of A has a support drawn uniformly at random from [m] (d) ; 2) supp(x) is drawn uniformly at random from [n] (k) ; , j ∈ supp(x); 6) y = Ax andŷ = y + η. We write (A,ŷ) ∼ GM(n, m, k, d, µ, ν) to denote problem instances drawn from this signal model.

Remark.
It is important to note that a matrix A ∈ R m×n generated under the model presented in Definition I.3 and Figure 1, is a (k, ε, d)-expander matrix with high probability, see [22], [15,Theorem 13.6].
Moreover, the generating model in Definition I.3 also allows us to define robust estimates for (3)-(4) for general noise and signal distributions and to any degree of accuracy under the assumption that these probability measures are available. From there we can define noisy analogues to the values n e and n z in (1)-(2) used in Parallel-0 but which are robust to additive noise. The contributions of this work are two-fold: (i) to present principled ways to compute (3)-(4) in the case where the nonzeros in η and x are drawn from continuous probability distributions; (ii) to provide a variation of Parallel-0 that is robust to noise. While other similar generating models can be considered using the techniques presented here, for ease of exposition and clarity, we restrict our description to this model.
Equations (5) and (6) allows us to quantify the uncertainty associated with computing the score for Parallel-0 under the presence of additive noise. Note that equations (5) and (6) can be easily adapted to alternative generative models, such as where the expected density of nonzeros per row varies, but for expository clarity we restrict our discussion to this somewhat generic model. It will be discussed in Section III-C that equations (5) and (6) should not be used directly, but instead be scaled by considering normalised functionsp e and p z defined as, In the most general case n e and n z can be written as the sum of individual scores q e and q z as follows A confidence threshold t > 0 needs to be given for some variants of our algorithm, so to simplify the exposition we include this parameter in the scores q e (· | t) and q z (· | t) regardless of whether it is used or not. A summary of the score functions are given in Table I.

A. Outline of the manuscript
The structure of this paper is as follows: Section II comprises a review of combinatorial compressed sensing and the Parallel-0 decoding algorithm extended here. Robust-0 decoding and the associated scores (8)- (9) are presented in Section III. In Section IV we present numerical experiments which demonstrate Robust-0 to perform superior to a number of leading greedy and combinatorial compressed sensing algorithms.
II. BACKGROUND: COMBINATORIAL COMPRESSED SENSING AND 0 -DECODING As mentioned in the previous section, the branch of combinatorial compressed sensing measures x ∈ χ n k with the adjacency matrix of an expander graph. These matrices are of very low complexity in terms of generation and storage, and also promise faster encoding and decoding than their dense counterparts, see Theorem II.2. In this section we review the basic elements of expander graphs and combinatorial compressed sensing.
There have been various algorithms proposed to reconstruct a sparse vector x from measurements y = Ax when A is an expander matrix, see [17], [18], [20] and [19]; the work presented here starts with the Parallel-0 algorithm and to improve this algorithm by making it robust to noisy measurements, results in Robust-0 (Algorithm 1). The key observation for the Parallel-0 algorithm is given by the following Lemma.
This means that at each iteration, if the residual r is nonzero, i.e. if we have not yet found the correct x, then there is Algorithm 1: Robust-0 Data: A ∈ E m×n k,ε,d ; y = Ax ∈ R m for x ∈ χ n k ; α ∈ (1, d]; µ, ν ∈ D(R); c ∈ (0, 1) Result:x ∈ R n s.t.x ≈ x Estimatep z =p z (d, k, m, n, µ, ν) as in (7); Estimatep e =p e (d, k, m, n, µ, ν) as in (7); if quantised then Use quantised scores given in Table I. else Use continuous scores given in Table I.
Recomputep e =p e (k 0 , m, n, µ, ν); end end end a set of entries in x that we can change so that we reduce the number of non-zeros in r by at least |T |(1 − 2ε)d.
Theorem II.2 (Convergence of Algorithm 2 [1]). Let A ∈ E m×n k,ε,d and let ε ≤ 1 4 , and x ∈ χ n k be dissociated. Then, To put this result into context and show its applicability, we recall the remark after Definition I.3, stating that random matrices as considered in this work are indeed expander matrices with high probability.
We furthermore emphasize the fact that the algorithm is designed in a way that allows for massively parallel implementations.

III. MAIN CONTRIBUTIONS: 0 -DECODING FOR NOISY MEASUREMENTS
We now consider the case where the measurements y are subject to additive noise, i.e. instead of y = Ax, we measurê y = y + η where η is a realization of a random variable with Parallel-0 is not able to cope with additive noise, as it needs to make decisions whether a value inŷ is zero and whether two values inŷ are equal to each other. While for very small noise levels we could consider two values as equal if they are within a certain number of standard deviations, for larger noise levels the decision becomes more challenging. A discussed in Section I we need to know p z (ŷ i ) and p e (ŷ i1 −ŷ i2 ) which correspond, respectively, to the probability of y i = 0 given that we observeŷ i and the probability of y i1 = y i2 given that we observeŷ i1 −ŷ i2 . The functions p z and p e depend on the parameters fed into the generative model in Definition I.3 and in particular on the distribution of y andŷ. Hence, since y is a vector of sparse inner products we need to understand the limiting behaviour of sparse sums.
We include Definition III.1 in order to remind the reader of some actions on measures used in this manuscript as so as to be relatively self contained.
we define the following measures.
1) The q-convolution, 2) The negative measure, 3) The symmetrized measure, 4) The measure associated with the difference of two random variables, Lemma III.2 (Limiting distribution for sparse sums of random variables). Let p ∈ (0, 1), let µ ∈ D(R) and let µ q ∈ D(R) be its the q-fold convolution. For each n ≥ 1, let be such that, Then, as n → ∞ it holds that s n Proof. Let ψ sn (t) be the characteristic function of s n . Let Taking the limit n → ∞ we see that Letting w q = q j=1 x j it holds by the independence of {x 1 , . . . , x q } that w q ∼ µ q and Now, consider a random variable z distributed according to The characteristic function of z is given by Therefore (15) equals (12). By Lévy's continuity Theroem pointwise convergence of the characteristic functions implies weak convergence of the random variables (cf. [23,Theorem 15.23]) and hence the statement follows.
Proof. To show (16), let i ∈ [m] and By our assumptions on A and x, where for two events E 1 and E 2 we let E 1 ∧E 2 be the conjunction of the events. Note that if Invoking Lemma III.2 with p n = p = dρ we obtain that y i → y * i as n → ∞ where By the independence of y * i and η i , sinceŷ i = y i + η i , the distribution ofŷ * i is given bŷ Equation (16) follows from (18) and the distributivity of the convolution operator.

A. Explicit formulas for the centred Gaussian case
We further elucidate (5)-(6) from Theorem I.4 in the case when µ and ν are Gaussian with mean zero, which are the distributions considered in Section IV. To this end, let µ = N (0, σ 2 s ) and ν = N (0, σ 2 n ). We observe that in this casē µ = µ, µ q =μ q = N (0, qσ 2 s ) andν = N (0, 2σ 2 n ). Hence, denoting by ϕ(· | σ 2 ) the probability density function of a centred Gaussian random variable with variance σ 2 , we obtain 1) Estimating the tails in the Gaussian case: We approximate the infinite sum in in the denominators of (26) and (27) by doing an approximation to the tail of this summation.
Lemma III.4 (Sums of centred density functions). Let µ i be the density function of a random varible with mean zero and variance σ 2 i and let α i > 0 be such that i α i = 1. Then, µ = i α i µ i is the density function of a random variable with mean zero and variance i α i σ 2 i .
To simplify notation let σ 2 q = qσ 2 s + σ n for q ∈ N ∪ {0}. Consider the series in the denominator of (22) Let, Note that S b /R z satisfies the conditions of Lemma III.4 so it corresponds to the density function of a centred random variable with variance A similar argument shows that .
Overall the analytical expressions fit the empirical probabilities very well, indicating that the approximations we made in the calculations above are justified. However, we observe that as ρ and especially σ increase, both functions drop significantly. This means that for large values of these parameters, the noise eventually dominates and it is difficult to decided whether values are zero or equal.

C. Scaled probabilitiesp z andp e
We mentioned in Section I that our algorithms do not implement the functions p z and p e exactly, but a scaled version of these. As observed in Figure 2, the value of max s p e (s) and max s p z (s) varies significantly as σ and ρ change. Algorithm 1 evaluates whether a given score is large or not by implementing a sweeping parameter t that is set to one at the beginning of the algorithm and decreased by a constant c after every iteration. In order to use a fixed initial t we consider the scaled probabilitiesp e andp z in (7); otherwise the initial value of t would depend on σ and ρ.

D. Adaptive k
Algorithm 1 can optionally account for the sparsity of the current estimate via the flag adaptive_k. In the noiseless model of r = Ax, an update of the formx ←x + ωe j with Parallel-0 guarantees that j ∈ supp(x) so that at the next iteration the problem with residual r − a j ω and (k − 1)sparse signal is considered. The adaptive_k flag updates the sparsity prior inx after every update in hope of having more reliable estimates of p e and p z . We will see in the numerical experiments that under the data generating model that we tested this strategy does not bring substantial benefits. We don't rule out the posssibility that there are other signal and noise distributions for which this flag becomes especially useful, but we leave that as future work.

IV. NUMERICAL EXPERIMENTS
In this section we present numerical experiments which validate the efficacy of Robust-0 decoding. In particular, we contrast Robust-0 with other state-of-the-art greedy algorithms for compressed sensing in terms of their ability to recover the measured signal for varying problem sizes (k, m, n) as well as their computational complexity. To facilitate reproducibility we begin by describing the stopping conditions and measures used to denote successful recovery in the presence of noise in Section IV-A, along with how the parameter c is varied in Section IV-B. We then present the main numerical results in Section IV-C where the algorithms phase transitions and runtime are presented, along with Sections IV-D and IV-E which show further details on Robust-0 decoding's performance as a function of noise variance and for extreme subsampling respectively.

A. Stopping conditions
We are interested in the signal model y = Ax + η where η ∼ N (0, σ 2 I m×m ). Ifx is an approximation to x the residual is r = y − Ax. Note that ifx = x, then and we should not seek reductions in the residual below η 1 which would result in fitting to the additive noise. We further account for the variance of η 1 and denote the algorithm to have successfully recovered x ifx satisfies for some C 1 ≥ 0. We should be aware that the right hand side of (30) might be greater than 1 for some choices of k, m and σ. When this happens, the stopping condition (30) becomes invalid since we expectx to have captured a proportion of the 1 -energy of x 1 . Hence, if the right hand side of (30) is greater than 1 10 we clip the upper bound at this value and use the stopping condition For the numerical experiments conducted in this section we consider nonzero entries in x drawn as x i ∼ N (0, 1), and noise η i ∼ N (0, σ 2 ) for which E [ x 1 ] = k 2 π and E [ η 1 ] = mσ 2 π and moreover Var [ η 1 ] = mσ 2 1 − 2 π , see e.g. [24].

B. Selection of parameter c in Algorithm 1
The sweeping parameter t in Algorithm 1 is initialised at 1 and updated by decreasing it by a constant value c. We observed in our experiments that the quality of the phase transitions are sensitive on the parameter c especially for low δ and ρ. We do not provide a way to fine-tune c, but we run our phase transitions with the following choices: 1) If the algorithm is quantised,  2) If the algorithm is continuous, The values in (32) and (33) were chosen heuristically for ν and µ Gaussian.

C. Phase transitions and runtime
We benchmark the variants of Robust-0 against other greedy algorithms via their phase-transitions and runtime. The user can supply two binary flags, adaptive_k and quantised which yield four different variants of Robust-0 . We assigned a unique label to each of these variants as described in Table II.
The phase transition of a compressed-sensing algorithm [25] is the largest value of k/m for which the algorithm is typically able recovery all k sparse vectors with sparsity less than k for a fixed m/n. Hence, for a fixed value of δ = m/n the phase transition of an algorithm is the largest value ρ * (δ) for which the algorithm converges for all ρ(δ) < ρ * (δ). The value ρ * (m/n) often converges to a fixed value as n → ∞, so phase transitions often partition the δ × ρ space into two regions: One in which the algorithm converges with high probability and another in which the algorithm doesn't converge with high probability. We benchmark Robust-0 against the algorithms presented in [18], [19], [26]. Specifically, our tests include the following algorithms, {Robust-0, Robust-0-adaptive, Robust-0-quantised, Robust-0-adaptive-quantised, SSMP, SMP, CGIHT}.
In the deterministic case Parallel-0 was compared against a range of combinatorial compressed sensing algorithms in [1]; out of those we have selected SSMP and SMP as these perform best and are similar in nature to Robust-0 . We also compare with CGIHT, as this algorithm was shown to be the fastest among the greedy algorithms compared in [13], [26]. Figures 3a, 3e and 3i show the phase transition curves for these algorithms with σ = 10 −3 , 10 −2 , 10 −1 respectively. The curves were computed by setting n = 2 18 , d = 7, and using the stopping condition r 1 ≤ E [ η 1 ] = mσ 2 π and a success condition (30) with C 1 = 1. The testing is done at m = δ p n for For each δ p , we set ρ = 0.01 and generate 10 synthetic problems to apply the algorithms to, with a problem generated as in GM with the given parameters and µ and ν being normal Gaussian. If at least one such problem was recovered successfully, the sparsity ratio ρ is increased by 0.01 and the experiment is repeated. Following the testing framework in [27], the recovery data is fitted using a logistic function and finally the 50% recovery transition function is computed and presented in the phase transition plots contained herein. Figures 3b, 3f and 3j show a selection map for these algorithms. Namely, these plots indicate which algorithm requires the least computational time 1 at each point in the δ × ρ space where the algorithm converges. Finally, Figures 3c, 3g and 3k show the total time for convergence in milliseconds for the fastest algorithm at each point in the δ×ρ space. The ratio of the time for the second fastest algorithm over the time for the fastest algorithm is given in Figures 3d, 3h, and 3l.
We can see from Figure 3 that CGIHT [26] dominates the upper region of the phase transition space, while the Robust-0 algorithms only converge for ρ 0.3 which is consistent with the observed phase transitions for Parallel-0 [1]. In terms of speed, Robust-0 seems to be the most competitive for σ ∈ {10 −3 , 10 −2 } and ρ 0.2. However, for large noise, σ = 10 −1 , CGIHT becomes the fastest algorithm of all. We remark that while CGIHT performs very well in our numerical tests, the current theory developed for it does not hold in the setting considered here, as it requires zero-mean columns in A. Figure 4 shows the widths for the Robust-0 algorithms. The widths measure how sharp the phase transition of an algorithm is; namely, how thin the boundary between the region of recovery with high-probability and the region of recovery where combinatorial search is needed. It has been shown that the widths of a compressed sensing algorithm tend to zero as n → ∞ when decoding with linear programming [29], and we usually expect the same behaviour for other algorithms [26]. Figure 4 show that the widths for the Robust-0 algorithms indeed decrease with n and with δ. The observed smoothness of the phase transition widths signal also suggest that the stopping conditions of the algorithm are consistent for the problem under consideration.

D. Dependence on noise variance, σ
We now investigate the extent to which the phase transitions of the algorithm decrease as we increase the noise level σ. To do this, we consider δ ∈ {0.01, 0.02, 0.1, 0.2} and for each value of δ we compute we define the grid Then at each value of σ we let ρ = 0.01 and draw ten problem instances from GM with signal distribution N (0, 1) and noise distribution N (0, σ 2 ). If at least one of the problems was recovered successfully, then we set ρ ← ρ+0.01 and repeat the experiment. We do this procedure for each of the algorithms considered and record the largest ρ having at least 50% success rate. The results are shown in Figure 5. In order to show where the clipping in (31) becomes active, the figures also show a TOL-curve which partitions the space into the region where the right hand side of (31) equals 1 10 (bottom region) and the region where it equals the right hand side of (30) (top region). We can appreciate from Figure 5 that for small δ both Robust-0 and SSMP have the best recovery capabilities, with Robust-0 being preferable for noise levels σ 10 −2 and SSMP being better suited for larger noise levels. For larger δ, CGIHT is preferable except for very low noise levels.
E. Phase transitions for extreme subsampling, δ 1 The numerical experiments of Parallel-0 in [1] showed flat phase transitions; that is, it was observed that ρ * (δ) remained approximately 0.3 as δ → 0 provided n was sufficiently large. While Robust-0 does not exhibit precisely the same behaviour in the presence of noise, we do observe that ρ * (δ) remains nontrivial even for δ as small as 10 −3 , again provided n is sufficiently large. We provide numerical evidence for this in Figure 6. For each (δ, σ) ∈ {0.001, 0.01} × {0.001, 0.01} we let ρ = 0.01 and solve ten problems drawn from GM with nonzero distribution N (0, 1) and noise distribution N (0, σ 2 ). If at least one problem instance converges, we average the run-time of the problems that converged and repeat the process with ρ ← ρ + 0.01. We plot the timing at each ρ for parameter values for which at which at least 50% of the problems were successfully recovered under the criteria (31).
It can be seen in Figure 6a-6d that the phase transition either remains nearly unchanged or increases as n increases from 2 22 to 2 24 . In particular, Figure 6a shows that for σ = 0.001 and δ = 0.001, the phase transitions for all variants of the Robust-0 algorithms in fact increase from ρ ≈ 0.08 to ρ ≈ 0.1 as n increases from 2 22 to 2 24 . Additionally, contrasting Figures  6a and 6b or 6c and 6d shows that a ten-fold increase in σ has the expected effect of reducing the phase transition and increasing the computational time. Figures 6a and 6b show the phase transition for Robust-0 and Robust-0 -adaptive remain significant even for δ as small as 10 −3 . Figures 6c and 6d show results for the same set of experiments, but for δ = 0.01 which corresponds to a ten-fold increase in δ over the value used in Figures 6a and 6b. For σ = 0.001 the phase transition of Robust-0 reaches ρ ≈ 0.17, while for σ = 0.01 the phase transitions drop to ρ ≈ 0.12 and there is an increase in the computational time.

V. CONCLUSIONS
We have shown that the decoding framework presented in [1] can be extended to the case where the measurements are corrupted by additive noise. This framework is extended by deriving the posterior distribution of an entry in the residual being zero or being equal to another residual entry given the corrupted measurements. This Bayesian approach to decoding was implemented in Robust-0 and its four variants. We show that the resulting algorithms inherits some desirable properties from Parallel-0 like high phase transitions for low δ and large n and low-latency. However, these qualities are weakened by the corruption of the measurements. Our numerical experiments show that Robust-0 should be considered in cases of moderate noise and ρ 0.3.


extfieldsofstudy:
Computer Science
Medicine

text:
MpBsmi: A new algorithm for the recognition of continuous biological sequence pattern based on index structure

A significant approach for the discovery of biological regulatory rules of genes, protein and their inheritance relationships is the extraction of meaningful patterns from biological sequence data. The existing algorithms of sequence pattern discovery, like MSPM and FBSB, suffice their low efficiency and accuracy. In order to deal with this issue, this paper presents a new algorithm for biological sequence pattern mining abbreviated MpBsmi based on the data index structure. The MpBsmi algorithm employs a sequence position table abbreviated ST and a sequence database index structure named DB-Index for data storing, mining and pattern expansion. The ST and DB-Index of single items are firstly obtained through scanning sequence database once. Then a new algorithm for fast support counting is developed to mine the table ST to identify the frequent single items. Based on a connection strategy, the frequent patterns are expanded and the expanded table ST is updated by scanning the DB-Index. The fast support counting algorithm is used for obtaining the frequent expansion patterns. Finally, a new pruning technique is developed for extended pattern to avoid the generation of unnecessarily large number of candidate patterns. The experiments results on multiple classical protein sequences from the Pfam database validate the performance of the proposed algorithm including the accuracy, stability and scalability. It is showed that the proposed algorithm has achieved the better space efficiency, stability and scalability comparing with MSPM, FBSB which are the two main algorithms for biological sequence mining.

Introduction
Biological sequence is an important component of bioinformatics data, generally including three categories: DNA sequence, RNA sequence and protein sequence [1]. Since the human genome project was completed in 2003, we have seen an explosive growth in bioinformatics data. By April 2017, there are 200877884 sequences in the GenBank [2] database. Since the release in 1982, the base number in GenBank has doubled by about every 18 months. Sequential pattern mining is an important method used for discovering frequent patterns and association rules in the data mining field [3], providing an effective way to find important rules of biological sequences. Biological sequence patterns can be used to predict human diseases and provide evidences for artificial nucleotides, artificial proteins and so on. Biological sequence pattern mining has become a hot research direction in recent years [4].
Frequent pattern mining is an important part of sequential pattern mining. Frequent pattern mining algorithms [5] contain two main categories, the Apriori based algorithm and the FP-growth based algorithm. The Apriori algorithm proposed by Agrawal et al [6] is the most commonly used in the association rule discovery. The frequent pattern mining has appeared in a variety of applications such as sequence pattern mining [7], structural mining [8], classification association mining [9], frequent pattern based on clustering [10] and so on. A number of sequential pattern mining algorithms have been proposed during the past years, such as incremental mining [11], top-k sequence pattern mining [12], maximum sequence pattern mining [13], constraint sequence pattern [14], weighted sequence pattern mining [15], closed sequence pattern mining [16][17]. To improve the efficiency and accuracy Oza K et al. [18] proposed an algorithm for regular expression constraint, weight constraint and length constraint to solve the problems of user interests, optimization for support threshold and accuracy. Xue F et al. [19] improved the PrefixSpan algorithm [20], and then proposed PrefixSpan-x to reduce unnecessary memory usage. Kemmar A et al. [21] proposed a top-k sequence pattern algorithm based on prefix projection and global constraints. This global constraint can take into account the quantity, item relation and regular expression and so on easily.
The classical sequence pattern algorithm is the footstone of biological sequence mining. The current sequence pattern algorithms are mainly improved from efficiency and precision. With the improvement of computer hardware performance, the algorithms with high memory utilization and changing time with space emerge. Lin et al. [22] proposed a fast sequence pattern algorithm named MEMISP based on memory index. This algorithm only needs scanning the sequence database once. An index table structure is proposed to record the positions of subsequence in the sequence and obtain the longer sequence pattern gradually in a recursive way. Its efficiency is higher than PrefixSpan and GSP algorithm. Then Ren et al. [23] proposed a closed sequential pattern mining algorithm MIWCSpan based on the memory index and item weights, considering the time interval of items and avoiding scanning database multiple times. Zeng et al. [24] proposed an algorithm based on the time interval weight and memory index, which further indicates the importance of time interval weight and improves the utilization of memory. Ren et al. [25] proposed an algorithm named MIFSPM, which uses the memory index structure of the frequent pattern tree. Each node in the tree stores the index table instead of storing the frequent pattern. This algorithm uses two kinds of support threshold constraint and sets the support threshold according to users' interaction.
In recent years, extensive research has been conducted on the biological sequence pattern mining focusing on the improvement of efficiency and precision. Yun [26] proposed an algorithm based on the prefix projection named BioPM, which can effectively mine the consensus sequences in the protein sequences. Chen et al. [27] proposed an algorithm MSPM based on the prefix tree and pattern expansion method. This algorithm abandons a large number of non-relevant candidate patterns, avoiding generating a large number of candidate projection sequences and many short candidate patterns which cost a lot of memory and time in dealing with complex biological sequences. Wang et al. [28] proposed an algorithm FBSB based on the bitmap storage structure and a fast sorted list. It avoids generating candidate patterns.
Parallel computing technology has also been studied for the sequence data mining [29]. Jiang et al. [30] used the Spark framework to mine sequence patterns by dealing with the uncertainty of DNA data. M Klein et al. [31] applied the Hadoop and Spark frameworks in their proposed biospark framework to deal with large data sets of biological sequences. Talouki MS et al. [32] conducted an algorithm based on a parallel prefix tree in handling protein data over several computers on a LAN network. The algorithm adopts the constraint of dynamic task assignment and selection sampling technology to avoid the machine idling and improve the precision.
The current biological sequence mining algorithms were generally focused on improvement of the classical sequence pattern algorithms or use the parallel computing. The existing sequence pattern algorithms suffer from the low efficiency, high memory utilization and low precision. To address these issues, this paper presents a biological sequence pattern mining algorithm MpBsmi based on the data index technology. This new algorithm has the following features: (1)A sequence position table ST and a sequence database index structure DB-Index are employed to improve the efficiency in scanning database and pattern support counting.
(2)A method is proposed to do the sequence pattern expansion easily through the position table and data index.
(3)A fast support counting method getting continuous sequence patterns is put forward.
(4)A novel method to prune the extended pattern is presented, which can filter out the subsequences whose support count is less than the threshold by index database, reducing the generation of a lot of candidate patterns.
Remaining of the paper is organized as follows. In Section 2, problems are defined, and the model is built. Section 3 presents our proposed algorithms including the position table, fast support counting algorithm and the updating algorithm of the position table and data index. Section 4 conducts the experiments and evaluate the performance. The discussions are made in Section 5. The paper is concluded in Section 6.

Definition and model
Biological sequences typically include DNA, RNA and protein sequences. They are different from general transaction sequences. DNA is a long chain polymer with four types of deoxy nucleotide: adenine (dAMP), thymine (dTMP), cytosine(dCMP), and guanine cytosine (dGMP). RNA is the ribonucleic acid which is a chain molecule formed by condensing two phosphate ester by DNA. RNA is composed of phosphoric acid, ribose and base. The bases of RNA are mainly in four kinds, i.e. adenine (A), guanine (G), cytosine (C) and uracil (U), where uracil (U) replaces thymine (T) in DNA. Proteins are formed by RNA translation and are large molecules containing 20 amino acid residues. In this chapter, formal concepts of biological sequence are given followed by the definitions and algorithm model. Parent sequence and child sequence. Given any two sequences BS 1 = <S 11 S 12 S 13 Á Á ÁS 1m > and BS 2 = <S 21 S 22 S 23 Á Á ÁS 2n >, and the conditions of BS 1 .type = BS 2 .type, m, n2Z + and m n are satisfied. If S 11 = S 21 , S 12 = S 22 ,Á Á Á,S 1M = S 2M , then BS 1 is the child sequence of BS 2 and BS 2 is the father sequence of BS 1 . That is to say BS 2 contains BS 1 , denoted as BS 1 BS 2 .
Example 1: Assuming a set ω = {A, C, D, E, F, G, H, I, K, L} of 10 different letters, the sequence BS 1 = <ACDEFGHIKL> can represent a protein sequence. If BS 2 = <ACDEFG>, then BS 2 BS 1 . Table 1 shows four sequences in the protein family indexed by PF00106 in the sequence database Pfam with the version 31.0 [33]. Length of biological sequence. Given a biological sequence BS 1 , the number of items contained in BS 1 is called the length of BS 1 which is written as |BS 1 |.
Biological sequence database. The database of biological sequence is a warehouse of biological sequences. It can be represented as BS_db = {BS 1 , BS 2 ,Á Á Á,BS n }, where BS_db[n] = BS n and n2Z + .
Length of BS_db. The length of BS_db is the number of sequences it contains, expressed as |BS_db|.

Definitions
Sequence position table and index database are two main data structures employed in this paper. Sequence position table records the locations of all the subsequences that appear in the sequences. The index database records all the indexes of the subsequences. Definition 1. Continuous biological sequence CBS. Given two biological sequences BS 1 = <S 11 S 12 S 13 Á Á ÁS 1m > and BS 2 = <S 21 S 22 S 23 Á Á ÁS 2n >, if S 21 = S 1y , S 2(1+1) = S 1(y+1) ,Á Á Á, S 2n = S 1(y+d) , and 1 y,y+d m,y,d,n,m2Z + , then BS 2 is called a continuous biological sequence CBS of BS 1 .
Definition 2. Subsequence. Given BS 1 = <S 11 S 12 S 13 Á Á ÁS 1m > and BS 2 = <S 21 S 22 S 23 Á Á ÁS 2n > are two biological sequences, if there exists S 2x = S 1y , S 2(x+1) = S 1(y+1) ,Á Á Á,S 2m = S 1m , and 1 x,1 y,x m n,y m n ,x,y, n, m2Z + , then BS 1 is called a subsequence of BS 2 , expressed as BS 1 Table 2 shows an example of a sequence database. From the table, we can get the items a, b, c contained in the sequences. In the first sequence, <abc> is a continuous biological sequence CBS and a subsequence of the first sequence. A part of ST for the single items is shown in Table 3. As can be obtained from Table 3 The times the index numbers containing subsequence a appear in DB-Index(a) is called the length of the sequence database index of subsequence a, expressed as |DB-Index(a)|. If every length of the sequences in the DB-Index is k, the DB-Index is also called k-DB-Index.
Definition 5. Biological sequence pattern BSP. Given a minimum support threshold minsup, a sequence database BS_db, a sequence BS 1 and BS 1 is a CBS, the support of BS 1 can be expressed as Support (BS 1 ) = |DB-Index(BS 1 )|/|BS_db|. If Support (BS 1 )!minsup, BS 1 is called a biological sequence pattern. |DB-Index(BS 1 )| is also known as the support count of the sequence BS 1 , expressed as Support (BS 1 ).

Algorithm model
The proposed algorithm to mine the pattern in the biological sequences has four steps. (1) scan the sequence database once to construct the position table 1-ST and database index 1-DB-Index of single items; then the fast support counting algorithm is used to get frequent sequence 1-BSP. (2)in order to expand (k+1)-BSP from k-BSP(k = 1, k2Z + ), we propose a sequence connection strategy named as BS-Ext to obtain (k+1)-BSP by connecting k-BSP.
(3)the index table ST and database index DB-Index will be updated to (k+1)-ST and (k+1)-DB-Index. (4) The (k+1)-BSP is obtained by the support counting algorithm. The recursive processes will be terminated until (k+1)-DB-Index is empty. These processes are shown in Fig  1. Definition 6. Sequence extension strategy BS-Ext. There are two k-BSP s BSP 1 and BSP 2 . If BSP 1 and BSP 2 can be connected, two conditions must be satisfied as below.

Build ST and DB-index
The subsequences can be obtained from scanning database. Constructing ST and DB-Index of a subsequence BS can be seen in algorithm 1.
If SP BS_i .size>0,then (7) DB-index(BS).add(i);; (8) End if (9) End if (10) i++; (11) End for End As can be seen from the algorithm 1, line (1) initializes the position table and the index database. Lines from (2) to (11) are the building processes, in which line (2) is a process of scanning databases once, line (5) finds all positions of BS appearing in BS 1 . If the BS appears in BS 1 , the index number of BS 1 in the database will be added to the DB-index in line (7).
Using Algorithm 1, we can first find the position table for all single items in the sequence database. And then the sequence extension strategy BS-Ext and ST and DB-Index updating algorithm will be used to obtain an expanded sequence.

Fast support counting and updating ST and DB-index
In this algorithm, the speed of computing support of sequence is a key factor to improve the efficiency of a sequence pattern mining algorithm. The pattern expansion technology is used to get the longer patterns. A new method of filtering out discontinuous patterns to prune and reduce unnecessary candidate patterns is described.
The algorithm is abbreviated as BSP-Distinguish. Algorithms 2. Fast pattern recognition algorithm Distinguish (BS,minsup). As can be seen from Algorithm 2, line (1) initializes support count, line (2) scans position table ST(BS), the non-empty SP(BS) is counted, and the support count is the number of nonempty SP(BS). By comparing with minsup, the BSP is obtained inline (8).
When we use the fast support counting to get 1-BSP from 1-ST, it is important to use the sequence extension BS-Ext to obtain the longer sequences. Details of this can be seen in Algorithm 3.

13) End for End
In the Algorithm 3, line (1) and line (2) traverse the k-BSP. According to line (3), the same database indexes of BSP 1 and BSP 2 are obtained. The number of the same indexes determines whether BSP 1 and BSP 2 satisfy the sequence extension strategy. In this step, the number of candidate patterns are reduced. Line (5) scans the database index DB-index. Line (6) determines whether the positions of BSP 1 and BSP 2 are adjacent to be connected. If BSP 1 and BSP 2 satisfy the extension conditions, Build-index-ST-DB(DB-index, BS_Ext) algorithm will be executed to acquire (k+1)-ST and (k+1)-DB-index after their connection. It is worth to notice that only database indexes are scanned in the process rather than the entire databases, which significantly improves the efficiency of the algorithm.

Biological sequence pattern mining
The biological sequence pattern mining algorithm Mpbsmi contains building position table ST and database index DB-index, the algorithm of frequent pattern recognition and updating ST and DB-index. The specific processes can be seen in Algorithm 4. SingleItemList.add(item) (5) End for (6) End for (7) For each item item in SingleItemList, Do In the Algorithm 4, line (1) to line (6) is the process of scanning the database to obtain all the single items. The lines from (7) to (10) get the 1-ST and1-DB-index, and then 1-BSP is found. Lines from (11) to (14) update the ST and DB-index by the sequence extension, and obtain longer frequent patterns. The algorithm is explained further in below example 4.
Example 4. In order to mine the biological sequence patterns from the biological sequence databases in Table 1, we assume minsup = 50%. The steps are as follows.

Experimental set-up and data sources
The experiment was conducted on the platform running Windows 10 operating system with 8GB RAM, CPU of E5200@2.5GHz. The proposed algorithm has been compared with MSPM [27] and FBSB [28]. All the algorithms are written in java. The experimental data of protein sequences come from protein database Pfam 31.0 [33]. The datasets contains 12 protein families, where 3000 sequences are selected and tested as shown in Table 4. The data of the algorithm scalability for data size are shown in Table 5. The data of the algorithm scalability for data length are shown in Table 6. The experiment of the relationship between the support threshold and the time cost was carried out by group experiment. The 3000 sequences were divided into 30 groups in average, in which the average length of sequence is 200.175. The support threshold changes from 5% to 40%. The results are the average values of 30 groups. The experimental results are shown in Fig 2. Although as depicted in Fig 2, the time costs of the algorithm Mpbsmi, FBSB and MSPM all decrease with the increase of the support threshold, the time consumed by the proposed algorithm Mpbsmi is significantly lower than FBSB and MSPM algorithms. The execution time of FBSB algorithm is around 10 times that of Mpbsmi algorithm. The main reason is that the Mpbsmi algorithm only processes sequence database index but the FBSB algorithm scans the entire sequence database. Meanwhile the time consumption of MSPM algorithm is significantly higher than Mpbsmi and FBSB algorithm, which is because MSPM algorithm produces a large number of basic candidate pattern s, sorting and combining the pattern is time-consuming. Mpbsmi reduces a large number of candidate patterns by sequence extension strategy BS-Ext, which reduces time cost. Experiment 2. Algorithm scalability analysis on data size of the relationship between data size and time consumption.

Performance evaluation of algorithm
The data used in the experiment is shown as Table 5. Every experimental group contains one protein family and selects 100, 200, 300, 400, 500 and 600 sequences respectively. The support threshold is set to 40%. The results are the averaged values of three groups. The experimental results are shown in Fig 3. As illustrated in Fig 3, the time consumption of algorithm MSPM rises linearly, and the time cost lines of the algorithm FBSB and Mpbsmi rise slowly. At the same time, the efficiency of Mpbsmi algorithm is noticeably higher than that of MSPM and FBSB, which indicates that Mpbsmi algorithm has better scalability of data size. Experiment 3. Algorithm scalability analysis on data length of the relationship between average sequence length and time consumption.
The experimental data are four protein families in Table 6. The average length of sequences changes from 100 to 400. Every protein family data are divided into two groups and the groups select 100 and 200 sequences respectively. The experimental results use the average values of every two groups. The support threshold is set to 40%, the experimental results are shown in Fig 4. As can be seen from Fig 4, with the average length of the sequence increasing, time consumption of the algorithm continues to increase. At the same time, the time costs of Mpbsmi are significantly lower than that of FBSB and MSPM. It is shown that the Mpbsmi algorithm has a good scalability in respect to data length. Experiment 4. The mining results of algorithm. In the case of a fixed support threshold of 40% with the same as Experiment 2, the result of the biological sequence patterns obtained are shown in Table 7: the first column shows the support threshold, the (k+1) th column and the k th column are the data set size and the corresponding number of pattern is mined, wherein 1<k<7 and k2Z + . As shown in the table, with the increment of the data set, the number of sequence pattern varies between 124 and 173, which reflects that the total single items of the biological sequences is constant with long sequence. Under the different data sets, if the support threshold is constant, the number of sequence patterns is relatively stable. When the support threshold is changed from 5% to 40% with the same data used in Experiment 1. The mining results are shown in Table 8: as Biological sequence pattern mining algorithm based on data index the support threshold increases, the number of sequence patterns is gradually reduced. In this experiment, we obtain all the continuous patterns. Because the algorithm can mine all successive sequence patterns, the algorithm's accuracy is 100% with no errors. Experiment 5. The memory and CPU usage analysis.
In the case of a fixed support threshold of 40% with the same experimental data used in Experiment 2, the memories and CPU used by the algorithm to process the data are shown in Table 9. The results are the average of three times running data. As one can see from Table 9, the peak memory usage by Mpbsmi algorithm is closer to that of the FBSB algorithm. Similarly, Mpbsmi and FBSB algorithms have the similar CPU peak values. The memory peak of MSPM algorithm is relatively higher than that of the other two algorithms but CPU peak is lower. It can be seen that Mpsmi and FBSB have lower memory and higher CPU utilization than MSPM algorithm, which improves the efficiency of algorithm and reduces the time consumption by using high CPU utilization.
The experimental results are used to calculate the P-value between the algorithm Mpbsmi and the FBSB, the algorithm Mpbsmi and the MSPM. The smaller the P-value, the more reason to think that the difference between things is a statistical indicator. As can be seen from Table 10, there are significant statistical differences in the efficiency difference of algorithm execution.

Discussion
In this paper, a new biological sequence pattern mining algorithm, the Mpbsmi, based on database index technique is presented. Compared to the two recently reported algorithms MSPM and FBSB, our proposed Mpbsmi algorithm uses a position table abbreviated ST and sequence index structure DB-Index with specially designed sequence pattern expansion and fast support counting method. Whilst the MSPM is an algorithm based on prefix tree and pattern extension, and FBSB is an algorithm based on bitmap structure and a quick sort list. By using index technology, our algorithm has shown significant performance improvement.  (1)Traditional frequent pattern mining algorithms The traditional frequent pattern algorithms [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] mentioned earlier in this paper are transaction based. They are designed for short sequences process with limited capability of processing long sequences. The biological sequences are characterized as ordered, continuous and elements repetitive. These features have constraint the efficiency of traditional algorithms. In terms of the internal data structures used by these algorithms, tables or tree structure are mostly used. The efficiency of query on large scale database is typically low. Literatures [23][24][25] have taken advantage of using the memory index as internal data structure, however they are vulnerable to memory overflow when dealing with longer sequences without sequence pattern restricting.
(2) Biological sequence pattern mining algorithm. The algorithms [26][27][28] mentioned in the paper are biological sequence pattern mining algorithms. The BioPM [26] algorithm is based on the prefix projection technology, MSPM [27] algorithm is based on prefix tree. They are both developed based on traditional sequential pattern technologies. The MSPM has shown higher efficiency than that of BioPM. FBSB [28] utilizes the bitmap as data structure and use this structure to calculate the support level. It has shown that FBSB algorithm has over performed MSPM algorithm in terms of efficiency.
(3) The algorithms based on the parallel computing technology. Another co-existing biological sequence mining algorithm stream is to take advantage of parallel computing. Literatures [29][30][31][32] used the MapReduce and Spark framework based on distributed statistical approach and are suitable for cloud computing environment running on multiple computers.
The proposed algorithm in this paper has proved its merit to obtain continuous biological sequence patterns efficiently. By further improvement, it can be used to acquire all the biological sequence patterns. There is also a plan to adapt the algorithm in cloud computing environment.

Conclusion
This paper proposes a new algorithm Mpbsmi based on data index technology for improving the efficiency of biological sequence pattern mining. Based on the index technique, a sequence position table ST is proposed to record the position information of sequences and subsequences. At the same time, the database index of sequences and subsequences is established as DB-index. The ST and DB-index are scanned only when the sequence pattern is extended. Through the sequence extension strategy and the position table, the algorithm achieves fast patterns mining and filters out a large number of invalid candidate patterns to get the entire continuous sequence patterns. Finally, the experimental results show that the proposed algorithm is superior to the existing biological sequence pattern mining algorithm such as MSPM and FBSB in terms of efficiency and scalability.