Mining Order-preserving Submatrices under Data Uncertainty: A Possible-world Approach and Efficient Approximation Methods

Before applying the dynamic programming algorithm, we first need to preprocess each row \(g=\langle [\ell _1, r_1], [\ell _2, r_2], \ldots , [\ell _m, r_m]\rangle\) to obtain a set of subintervals demarcated by \(\ell _i\) , \(r_i\) \((i=1,\ldots ,m)\) . These subintervals, denoted by \(I_1, \ldots , I_s\) , are ordered in increasing order of value. Note that \(s\le 2m\) . We also denote the interval of \(I_i\) as \([\ell (I_i), r(I_i)]\) .

Figure 4 shows the subintervals obtained by preprocessing row \(g_1\) in the data matrix in Figure 3: \(I_1=[38, 49]\) , \(I_2=[49, 50]\) , \(I_3=[49, 79]\) , \(I_4=[79, 80]\) , \(I_5=[80, 81]\) , \(I_6=[81, 110]\) , and \(I_7=[110, 115]\) .

Obviously, for each row \(g\) , given an interval element \(D[g][t]\) and a subinterval \(I_i\) , we must have: either (1) \(I_i\subseteq [\ell _t, r_t]\) or (2) \(I_i\cap [\ell _t, r_t]=\emptyset\) or \(\lbrace \ell _t\rbrace\) or \(\lbrace r_t\rbrace\) . For example, in Figure 4, \(I_2\subseteq t_1\) , \(I_2\subseteq t_2\) , \(I_2\cap t_3=\emptyset\) , and \(I_2\cap t_4=r(I_2)=\ell _4=\lbrace 50\rbrace\) .

Since each \(D[g][t]\) is assumed to be uniform, we have:

Therefore, we can obtain the constant probability density of \(D[g][t]\) on any subinterval \(I_i\) in \(O(1)\) time by checking whether \(I_i\subseteq [\ell _t, r_t]\) . Since \(f_t(x)\) is constant on \(I_i\) , we abbreviate it as \(f_t\) from now on. For example, in Figure 4, consider \(g_1\) and subinterval \(I_2\) : \(f_1=1/31\) , since \(I_2\subseteq [49, 80]\) , while \(f_4=0\) except at the boundary 50 (which is immaterial for continuous distribution). As we shall see next, Property 1 is critical to the efficiency of computing supporting probability.

Assume data matrix \(D\) is \(n\, \times \, m\) , preprocessing each row into intervals \(I_1, \ldots , I_s\) requires sorting \(s\le 2m\) values with the cost of \(O(s\log s)\) time. Overall, processing the \(n\) rows of \(D\) takes \(O(ns\log s) = O(nm\log m)\) .

Given a row \(g\) , a pattern \(P=(t_{i_1}\lt t_{i_2}\lt \cdots \lt t_{i_\ell })\) , and the subintervals \(I_1, I_2, \ldots , I_s\) by preprocessing \(g\) , we compute \(p_g(P)\) using dynamic programming (DP) as follows. We abuse the notation \(t_i\) to mean the interval \(D[g][t_i]\) , since \(g\) is given.

The DP algorithm first creates a 2D array \(A\) with \(\ell\) rows and \(s\) columns. The element \(A[j][k]\) denotes the probability of the event \(t_{i_1}\lt \cdots \lt t_{i_j}\) with \(t_{i_1}, \ldots , t_{i_j}\) located in the interval consisting of the first \(k\) subintervals, i.e., \(\bigcup _{i=1}^k I_i\) . Note that the value of \(A[\ell ][s]\) is exactly \(p_g(P)\) .

Our algorithm assumes that different column intervals of a row are independent, which is natural, since different microarray tests or RFID readers are usually independent. This assumption is used in many places below, such as the Proof of Theorem 3.1 and the last term’s product in Equation (6).

Probability Evaluation on One Subinterval. Before describing the recursive formula for computing \(A[j][k]\) , we first explain how to compute the probability of the event \(t_{i_1}\lt t_{i_2}\lt \cdots \lt t_{i_\ell }\) with \(t_{i_1}, \ldots , t_{i_ \ell }\) located in \(I_k\) . We denoted this probability by \(P_{I_k}(t_{i_1}\lt \cdots \lt t_{i_\ell })\) .

Theorem 3.1 implies the following corollary for row \(g\) :

Evaluation of \(A[j][k]\) . Now, we are ready to present the recursive formula for computing \(A[j][k]\) . Let us first consider the base case when \(k=1\) . In this case,which can be computed using Equation (4).

When \(k\gt 1\) , the event “ \(t_{i_1}\lt t_{i_2}\lt \cdots \lt t_{i_j}\) with \(t_{i_1},t_{i_2},\ldots ,t_{i_j}\in \bigcup _{i=1}^k I_i\) ” can be decomposed into the following disjoint events:

According to the above discussion, we obtain

In fact, if we define \(A[j][0]=0\) for any \(j\) , then we can even compute \(A[j][1]\) using Equation (6).

Note that computing \(A[j][k]\) involves:

According to Equation (4), computing \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) takes \(O(j-z+1)\) time. Thus, computing \(A[j][k]\) takes \(\sum _{z=1}^j O(j-z+1)=O(j^2)\) time if we compute each \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) individually.

We can actually compute \(A[j][k]\) in \(O(j)\) time as follows: From Equation (4), we can derive the following recursive formula for computing \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) :

Therefore, if we compute \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) with \(z\) from \(j\) down to 1, then each \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) can be computed from \(P_{I_k}(t_{i_{z+1}}\lt \cdots \lt t_{i_j})\) using Equation (7) in O(1) time. Thus, computing \(A[j][k]\) takes \(\sum _{z=1}^j O(1)=O(j)\) time.

Accordingly, computing \(p_g(P)\) requires computing all elements in the \(\ell \times s\) array \(A\) , which takes \(\sum _{j=1}^\ell \sum _{k=1}^s O(j)=O(\ell ^2s)=O(\ell ^2m)\) time (recall that \(s\le 2m\) ). We can further optimize the computation: If we find \(f_{t_{i_z}}=0\) when evaluating \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) , then we can terminate early, since according to Equation (4), \(P_{I_k}(t_{i_{z^{\prime }}}\lt \cdots \lt t_{i_j})=0\) for any \(z^{\prime }\le z\) .

From now on, let us call the array for dynamic programming as DP-array. We mine patterns by pattern-growth, i.e., pattern \(t_{i_1}\lt \cdots \lt t_{i_j}\lt t_{i_{j+1}}\) is checked after pattern \(t_{i_1}\lt \cdots \lt t_{i_j}\) .

We now consider how to reuse the DP-array for \(t_{i_1}\lt \cdots \lt t_{i_j}\) to compute the DP-array for \(t_{i_1}\lt \cdots \lt t_{i_j}\lt t_{i_{j+1}}\) .

In the base case, the pattern is a singleton \(t_i\) . Assume that the interval of \(D[g][t_i]\) is \([\ell _i, r_i]\) , then there is only one subinterval \(I_1=[\ell _i, r_i]\) and the DP-array is a \(1\times 1\) array with \(A[1][1]=Pr_{I_1}(t_i)=1\) .

We now consider how to incrementally compute the DP-array of a pattern grown by one more interval. Referring to the data matrix in Figure 3 again, let us focus on the computation of \(p_{g_1}(t_3\lt t_4\lt t_1)\) .

Figure 5 illustrates the evaluation process, where the array on the left is the DP-array for pattern \(t_3\lt t_4\) that is already computed, and the array on the right is the DP-array for pattern \(t_3\lt t_4\lt t_1\) that is to be computed. The split points of subintervals are also marked above the DP-arrays. Note that \(A[j][k]=Pr\lbrace (t_{i_1}\lt \cdots \lt t_{i_j})\ \wedge \ (t_{i_1},\ldots ,t_{i_j}\lt r(I_k))\rbrace\) . For example, the value of cell \(E\) in the left array is the probability of the event \(t_3\lt t_4\) with \(t_3, t_4\lt 110\) . Obviously, the value of cell \(E\) in the right array is exactly the same. In fact, we can copy the values of cells \(A\) – \(F\) in the left array directly into the corresponding cells in the right array.

However, for pattern \(t_3\lt t_4\lt t_1\) , the introduction of \(t_1\) with interval \([49, 80]\) adds two more split points. For example, \(I_2=[79, 110]\) of the left array is now split into two subintervals for the right array: \(I_3=[79, 80]\) and \(I_4=[80, 110]\) . However, the value of cell 4 (the probability of the event \(t_3\lt t_4\) with \(t_3, t_4\lt 80\) ), for example, is not computed in the left array, and therefore it has to be computed using Equation (6). In fact, the values of cells 1–9 have to be computed over the right array using Equation (6).

Algorithm 1 shows how we compute the DP-array \(A_1\) for pattern \(t_{i_1}\lt \cdots \lt t_{i_j}\) from the DP-array \(A_0\) for pattern \(t_{i_1}\lt \cdots \lt t_{i_{j-1}}\) and the interval \([\ell _{i_j}, r_{i_j}]\) for \(t_{i_j}\) . First, the new ordered list of split points, \(L_1\) , is constructed by merging \(\lbrace \ell _{i_j}, r_{i_j}\rbrace\) with the old split point list \(L_0\) in Line 2, which takes \(O(s_1)\) time ( \(s_1=|L_1|\) ). Then, Lines 5–12 compute the first \((j-1)\) rows of \(A_1\) : If the column-right-marks of the elements \(A_1[row][col]\) and \(A_0[row][pos]\) are the same (Line 8), then the value of \(A_1[row][col]\) is copied from \(A_0[row][pos]\) ; otherwise, \(A_1[row][col]\) has to be computed using the dynamic programming formula (Line 12). Finally, since \(A_0\) does not have the \(j\) th row, elements \(A_1[j][col]\) have to be computed using dynamic programming (Lines 14–15).

We now analyze the cost of this incremental computation of \(p_g(t_{i_1}\lt \cdots \lt t_{i_j})\) . Since at most two new split points \(\ell _{i_j}\) and \(r_{i_j}\) are introduced when a pattern is grown with \(t_{i_j}\) , for each of the first \((j-1)\) rows in the DP-array, there are at most two elements to compute using dynamic programming. Together with the new \(j\) th row, there are totally \(2(j-1)+(s-1)=O(j)\) elements (note that the number of split points \(s\lt 2j\) ) to compute using dynamic programming, which takes \(O(j^2)\) time. The remaining \((j-1)\cdot s\) elements are copied from the old DP-array, which takes \(O(j\cdot s)\) time. Therefore, the total time complexity is \(O(j^2)\) , quadratic to the pattern length \(j\) .

For a pattern \(P\) of length \(\ell\) , time complexity of computing \(p_g(P)\) is reduced from \(O(\ell ^2m)\) in Section 3.2 to \(O(\ell ^2)\) here.

We have been assuming that each matrix entry conforms to a uniform distribution defined over its interval, which is a proper assumption to avoid inductive bias when there are only several replicates.

When there are sufficient number of replicates, we can infer the underlying distribution (e.g., Gaussian distribution) and learn the parameters from the replicates or fit the underlying distribution by kernel density estimation. Then, we can discretize the PDF using an equi-width or equi-depth histogram. Note that our dynamic programming framework is still applicable here, though each entry may now introduce more than two split points. In fact, if we discretize each distribution into \(n_s\) ranges (i.e., \(n_s+1\) split points), then for a pattern \(P\) of length \(\ell\) , time complexity of computing \(p_g(P)\) is increased from \(O(\ell ^2)\) in Section 3.3 to \(O(n_s^2\cdot \ell ^2)\) here. The weakness of histogram-based approximation is that, the probability density within each range is treated as uniform, so we need to discretize a PDF with many small ranges to reach reasonable accuracy, but the time complexity of computing \(p_g(P)\) grows quadratically with \(n_s\) .

We consider an alternative solution that fits a PDF with a cubic spline, which is a spline constructed of piecewise third-order polynomials passing through a set of control points. This approach is generally applicable to the PDF of an arbitrary distribution. We next illustrate this approach using the Gaussian distribution that is commonly used to model noises in measurements, as well as the exponential distribution as a second example. We also provide a discussion on how to choose appropriate control points for a general distribution.

Control Points for the Gaussian Distribution. Given a set of replicates for a matrix entry, we compute the sample mean \(\mu\) and variance \(\sigma ^2\) to obtain a Gaussian distribution \(\mathcal {N}(\mu , \sigma ^2)\) ; we then fit its PDF with a cubic spline using critical points at \(\mu -4\sigma\) , \(\mu -3\sigma\) , \(\mu -2\sigma\) , \(\mu -\sigma\) , 0, \(\mu +\sigma\) , \(\mu +2\sigma\) , \(\mu +3\sigma\) , \(\mu +4\sigma\) . Figure 6 illustrates the PDF of \(\mathcal {N}(0, 1)\) as well as the cubic spline fit from the PDF’s nine critical points; we can see that the cubic spline matches tightly with the Gaussian PDF, so the probability approximation is much more accurate. Also, since 99.73% values in a Gaussian distribution lies within \([\mu -3\sigma , \mu +3\sigma ]\) , we only consider the six intervals \([\mu -3\sigma , \mu -2\sigma ]\) , \([\mu -2\sigma , \mu -\sigma ]\) , \([\mu -\sigma , 0]\) , \([0, \mu +\sigma ]\) , \([\mu +\sigma , \mu +2\sigma ]\) , \([\mu +2\sigma , \mu +3\sigma ]\) when computing \(p_g(P)\) , but we also use \(\mu \pm 4\sigma\) as critical points to ensure smooth fitting on the border of \(\mu \pm 3\sigma\) . For each matrix entry \(D[g][t_j]\) , we now have seven split points \(\mu -3\sigma\) , \(\mu -2\sigma\) , \(\mu -\sigma\) , 0, \(\mu +\sigma\) , \(\mu +2\sigma\) , \(\mu +3\sigma\) , and the PDF between any two consecutive split points is determined by a unique third-order polynomial \(p^{(j)}(x)=a^{(j)}_0+a^{(j)}_1x+a^{(j)}_2x^2+a^{(j)}_3x^3\) .

Control Points for the Exponential Distribution. Given a set of replicates for a matrix entry, \(\lbrace v_1, v_2, \ldots , v_k\rbrace\) , we compute shift parameter \(min = \min _i\lbrace v_i\rbrace\) , shifted mean \(\mu =\frac{1}{N}\sum _i (v_i - min)\) , and rate parameter \(\lambda = \frac{1}{\mu }\) to obtain an exponential distribution \(Exp(\lambda , min)\) whose PDF is given by \(\lambda e^{-\lambda (x-min)}\) \((x\ge min)\) . Figure 7 illustrates the PDF of a standard exponential distribution \(Exp(1, 0)\) as well as the cubic spline fit from the PDF’s six critical points at 0, \(\mu\) , \(2\mu\) , \(3\mu\) , \(4\mu\) , and \(5\mu\) . Note that 99% values of an exponential distribution lies within \([0, 4.61\mu ]\) , so we only consider five intervals \([0, \mu ]\) , \([\mu , 2\mu ]\) , \([2\mu , 3\mu ]\) , \([3\mu , 4\mu ]\) , and \([4\mu , 5\mu ]\) when computing \(p_g(P)\) . For each matrix entry \(D[g][t_j]\) , we now have six split points \(min\) , \(min+\mu\) , \(min+2\mu\) , \(min+3\mu\) , \(min+4\mu\) , \(min+5\mu\) , and the PDF between any two consecutive split points is determined by a unique third-order polynomial \(p^{(j)}(x)=a^{(j)}_0+a^{(j)}_1x+a^{(j)}_2x^2+a^{(j)}_3x^3\) . Guidelines on Choosing Control Points for a General Distribution. Now that we have seen how to choose good control points for Gaussian and exponential distributions, we are ready to discuss how to choose proper control points for a general distribution. We remark that in our context, the most reasonable distributions are unimodal distributions where most probability density is clustered around the ground-truth value of a matrix entry, such as the exponential family that comprises the natural sets of distributions to consider in reality. For such a general distribution, the probability density drops quickly as the value increases beyond the mean value, so we only need to focus on fitting the limited value range that accounts for the vast majority of the value possibility (e.g., 99% or more), since ignoring the rest would not compromise much accuracy. For this value range considered, we can either divide the range evenly as we do for Gaussian and exponential distributions or by selecting intervals in sequence where each interval is selected by increasing its length till the point when a goodness-of-fit metric drops below a user-defined threshold (binary search applies here for efficiency). The number of intervals or the goodness-of-fit threshold can be adjusted according to users’ accuracy needs, but there is an efficiency-accuracy tradeoff that users need to consider. However, as we have shown in Figures 6 and 7, a small number of five to six intervals is often sufficient to achieve very good fitting quality for a common distribution falling into the exponential family.

Computing \(p_g(P)\) . We can still use the same computation framework of Algorithm 1 in Section 3.3 when \(D[g][t_j]\sim \mathcal {N}(\mu , \sigma ^2)\) (or when \(D[g][t_j]\sim Exp(\lambda , min)\) ). One difference is that Line 2 now merges the split points in \(L_0\) with the seven (or six) split points of \(D[g][t_j]\) to obtain the new split point list \(L_1\) . For any two consecutive split points \(\ell ,r\) in \(L_1\) , the PDF on the interval \([\ell , r]\) is defined by a unique third-order polynomial for each of \(t_{i_1}\) , \(t_{i_2}\) , \(\ldots\) , \(t_{i_j}\) .

The other difference in Algorithm 1 is the computation of \(A_1[row][col]\) as specified in Lines 12 and 15, which uses Equations (5) or (6) where the key operation is to compute the probability \(P_{I_k}(t_{i_1}\lt t_{i_2}\lt \cdots \lt t_{i_j})\) on subinterval \(I_k=[\ell (I_k), r(I_k)]\) , with \(\ell (I_k)\) and \(r(I_k)\) being two consecutive points in \(L_1\) . We next discuss its computation assuming \(p^{(j)}(x)=a^{(j)}_0+a^{(j)}_1x+a^{(j)}_2x^2+a^{(j)}_3x^3\) :

Even though we cannot derive a closed form formula as Equation (4) to evaluate Equation (8), Equation (8) can still be computed analytically to avoid expensive numerical evaluation. Specifically, let us define an order- \(k\) polynomial \(poly(x, \mathbf {a})=\sum _{i=0}^k a_ix^i\) , which we implement as a class Polynomial that maintains an array of coefficients \(\mathbf {a}=[a_0, a_1, \ldots , a_k]\) and that supports two operations: (1) multiplying another polynomial (which generates a new polynomial with a new coefficient array) and (2) computing the following integral where \(y\) is a variable while \(\ell\) is a constant:where \(b_i\) and \(c_i\) are constants. In other words, the integral over \(poly(x, \mathbf {a})=\sum _{i=0}^k a_ix^i\) gives a new order- \((k+1)\) polynomial \(poly(y, \mathbf {b})=\sum _{i=0}^{k+1} b_iy^i\) .

Using Equation (9), we can compute Equation (8) recursively from \(x_1\) all the way till \(x_n\) as follows:

Regarding the time complexity, evaluating the polynomial multiplication in Equation (10) takes \(O([4(j-1)+1]\times [3+1])=O(j)\) time, and evaluating the integral in Equation (10) to obtain Equation (11) also takes \(O(j)\) time. So, evaluating the probability in Equation (12) takes \(O(n)\) time, so according to Equation (6), computing \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) takes \(O(j-z+1)\) time. Thus, computing \(A[j][k]\) takes \(\sum _{z=1}^j O(j-z+1)=O(j^2)\) time, as we compute each \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) individually.

Accordingly, computing \(p_g(P)\) requires computing all elements in the \(\ell \times s\) array \(A\) , which takes \(\sum _{j=1}^\ell \sum _{k=1}^s O(j^2)=O(\ell ^3s)=O(\ell ^3m)\) time, where the number of split points \(s\le 7m\) in our Gaussian approximation scheme in Figure 6, and \(s\le 6m\) in our exponential approximation scheme in Figure 7. We can further optimize the computation:If \(I_k\) is outside the \([-3\sigma , 3\sigma ]\) -interval (of the Gaussian distribution) or \([0, 5\mu ]\) -interval (of the exponential distribution) of some column value \(t_{i_z}\) when we evaluate \(P_{I_k}(t_{i_z}\lt \cdots \lt t_{i_j})\) , then we can terminate early, since we have \(P_{I_k}(t_{i_{z^{\prime }}}\lt \cdots \lt t_{i_j})=0\) for any \(z^{\prime }\le z\) .

We define the support of a pattern \(P\) in data matrix \(D=(G, T)\) as the number of rows in \(G\) that supports \(P\) . Since each row \(g\) supports \(P\) only with probability \(p_g(P)\) , the support of \(P\) is a random variable, denoted as \(X\) .

We now consider how to compute the PMF of \(X\) , using the supporting probabilities \(p_{g}(P)\) of all rows \(g\in G\) .

Naïve Method. Let \(f_i(c)\) be the PMF of the support of \(P\) in matrix \(D_i\) , which consists of the first \(i\) rows in \(D\) , where \(c=0,1,\ldots ,i\) ( \(f_i(c)=0\) for other values of \(c\) ). Then, we have the following recursive formula:

This is because \(P\) ’s support in \(D_{i+1}\) is \(c\) , iff (1) \(P\) ’s support in \(D_i\) is \((c-1)\) and \(g_{i+1}\) supports \(P\) , or (2) \(P\) ’s support in \(D_i\) is \(c\) and \(g_{i+1}\) does not support \(P\) . Note that Equation (13) holds only for \(c=1,\ldots ,i\) , and the products of probabilities are due to row independence (the same for later equations and thus omitted).

When \(c=0\) , we havesince the support of \(P\) in \(D_{i+1}\) is 0, if and only if the support of \(P\) in \(D_i\) is 0 and \(g_{i+1}\) does not support \(P\) .

When \(c=i+1\) , we havesince the support of \(P\) in \(D_{i+1}\) is \((i+1)\) , if and only if the support of \(P\) in \(D_i\) is \(i\) and \(g_{i+1}\) supports \(P\) .

Furthermore, we have the base casesince there is no row in \(D_0\) and the support of \(P\) must be 0.

For a data matrix \(D\) with \(n\) rows, the PMF of \(X\) , denoted by \(f_P(c)\) , is equal to \(f_n(c)\ (c=0,\ldots ,n)\) , since \(D=D_n\) . To compute the PMF \(f_P(c)\) , we start with \(f_0(c)\) and recursively compute \(f_{i+1}(c)\) from \(f_i(c)\) until \(f_n(c)\) is computed. According to Equations (13), (14), and (15), it takes \(O(i)\) time to compute \(f_i\) and thus \(O(n^2)\) time to compute \(f_P(c)\) .

Divide-and-conquer Algorithm. The naïve method takes time quadratic to the number of rows, which does not scale well. We now describe another algorithm for computing the support PMF \(f_P(c)\) , which achieves better time complexity by the divide-and-conquer strategy.

Given a set \(S\) of rows, let us define \(f^S(c)\) as the PMF of the support of \(P\) in row set \(S\) , then our ultimate goal is to compute \(f^D(c)\) . To compute \(f^S(c)\) , we first divide the rows in \(S\) into two sets \(S_1\) and \(S_2\) of equal size. Let us denote \(|S|\) by \(n\) , then \(S_1\) contains the first \(\lfloor \frac{n}{2}\rfloor\) rows, and \(S_2\) contains the remaining \(\lceil \frac{n}{2}\rceil\) rows.

Assume that \(f^{S_1}(c)\) and \(f^{S_2}(c)\) are already computed, then \(f^S(c)\) can be obtained by the following formula:since the event that “the support of \(P\) in \(S\) is \(c\) ” can be decomposed into the disjoint events “the support of \(P\) in \(S_1\) is \(i\) , and the support of \(P\) in \(S_2\) is \((c-i)\) ” for \(i=0,\ldots ,c\) .

Note that \(|S|=n\) , while \(|S_1|=\lfloor \frac{n}{2}\rfloor\) and \(|S_2|=\lceil \frac{n}{2}\rceil\) . In Equation (17), we define \(f^{S_1}(c)=0\) for \(c\gt \lfloor \frac{n}{2}\rfloor\) , and define \(f^{S_2}(c)=0\) for \(c\gt \lceil \frac{n}{2}\rceil\) .

According to Equation (17), \(f^S\) is the convolution of \(f^{S_1}\) and \(f^{S_2}\) . Therefore, \(f^S\) can be computed from \(f^{S_1}\) and \(f^{S_2}\) in \(O(n\log n)\) time using Fast Fourier Transform (FFT) [8].

Our divide-and-conquer algorithm for computing \(f^S\) is described as follows: \(S\) is first divided into two row sets \(S_1\) and \(S_2\) of equal size; then, PMFs \(f^{S_1}\) and \(f^{S_2}\) are computed by recursion; finally, \(f^S\) is computed as the convolution of \(f^{S_1}\) and \(f^{S_2}\) using FFT. Since each recursion step takes \(O(n\log n)\) time (due to FFT), the overall time complexity of computing \(f^D(c)\) is \(O(n\log ^2 n)\) .

The base case for recursion is when \(S=\lbrace g\rbrace\) , in which case, we directly return \(f^S\) where \(f^S(0)=1-p_g\) and \(f^S(1)=p_g\) . In reality, we find that recursion down to \(|S|=1\) does not provide the best performance. The most efficient configuration is to stop recursion when \(|S|\le 500\) and compute \(f^S\) directly using the naïve method. We adopted this implementation.

In the previous subsection, we described how to compute the PMF of the support of pattern \(P\) in data matrix \(D\) . Once the PMF is computed, we can decide whether pattern \(P\) is p-frequent using Equations (1) and (2) in Section 2.

However, this two-step approach is time-consuming, since it requires to compute the whole PMF vector \(f_P(c),\ c=0,\ldots ,n\) . In fact, to determine whether pattern \(P\) is frequent, it is not always necessary to compute \(f_P\) to the end. The theorem below states that pattern \(P\) is frequent in \(D\) , as long as it is found to be frequent in a subset of \(D\) . As a result, in our divide-and-conquer algorithm, in each recursion step (that computes \(f^S\) ), we will check the frequentness of \(P\) over \(S\) using Equations (1) and (2), and if it is found to be frequent, we terminate the frequentness checking immediately and conclude that \(P\) is frequent.

The frequentness checking operations described above is still very expensive (i.e., \(O(n\log ^2n)\) time). We now present three pruning rules for pruning infrequent patterns. These rules can be checked efficiently in \(O(n)\) time, and if any rule determines that a pattern \(P\) is infrequent, we do not need to do the expensive frequentness checking for \(P\) . The proofs of these rules can be found in Appendix B.3.

(1) Count-prune. Let \(cnt(P)=|\lbrace g\in G\,|\,p_g(P)\gt 0\rbrace |\) , then pattern \(P\) is not p-frequent if \(cnt(P)\lt \tau _{row}\) .

(2) Markov-prune. Pattern \(P\) is not p-frequent if

(3) Exponential-prune. Let \(\mu =E(X)\) and \(\delta =\frac{\tau _{row}-\mu -1}{\mu }\) . When \(\delta \gt 0\) , pattern \(P\) is not p-frequent if

To further improve the efficiency of examining pattern frequentness using probabilistic frequentness, instead of computing the exact PMF vector, we can use the Poisson or Gaussian distribution to approximate the PMF, which reduces the time complexity from \(O(n\log ^2 n)\) to \(O(n)\) . We next describe the two PMF approximation techniques using Poisson and Gaussian distributions, respectively.

Consider the matrix in Figure 3; it is obvious that \(g_2\) can never support pattern \(t_3\lt t_4\) , i.e., \(p_{g_2}(t_3\lt t_4)=0\) . We now describe how to determine whether \(p_{g}(P)=0\) efficiently without actually evaluating the supporting probability.

Given a pattern \(P=(t_{i_1}\lt t_{i_2}\lt \cdots \lt t_{i_j})\) , we define its valid interval as the interval \([\ell _P, r_P]\) such that \(p_{g}(P)\gt 0\) if and only if \(t_{i_j}\in [\ell _P, r_P]\) . We now consider how to compute the valid interval of a pattern \(P\) .

In the base case when \(P=t_{i_1}\) , its valid interval is exactly the interval \(D[g][t_{i_1}]=[\ell _{i_1}, r_{i_1}]\) .

For pattern \(P=(t_{i_1}\lt \cdots \lt t_{i_j})\) \((j\gt 1)\) , we compute its valid interval using that of pattern \(P^{\prime }=(t_{i_1}\lt \cdots \lt t_{i_{j-1}})\) . Since row \(g\) supports \(P^{\prime }\) if and only if \(t_{i_{j-1}}\in [\ell _{P^{\prime }}, r_{P^{\prime }}]\) , row \(g\) would support \(P\) if and only if \(\exists t_{i_{j-1}}\in [\ell _{P^{\prime }}, r_{P^{\prime }}]\) such that \(t_{i_{j-1}}\lt t_{i_j}\) , or equivalently, \(t_{i_j}\in [\max \lbrace \ell _{i_j}, \ell _{P^{\prime }}\rbrace , r_{i_j}]\triangleq [\ell _P, r_P]\) (where \([\ell _{i_j}, r_{i_j}]\) is the interval \(D[g][t_{i_j}]\) ). Note that \([\ell _P, r_P]=\emptyset\) if \(\ell _P\gt r_P\) , and in this case \(p_g(P)=0\) and \(g\) can be filtered.

In our mining algorithm when checking pattern \(P^{\prime }\) , for each row \(g\) , we maintain the following information: (1) its valid interval \([\ell _{P^{\prime }}, r_{P^{\prime }}]\) , (2) the ordered list of split points, and (3) the DP-array. Since our mining algorithm checks patterns by pattern-growth, \(P\) will be checked after \(P^{\prime }\) . To process \(g\) for pattern \(P\) , we will first compute \([\ell _P, r_P]=[\max \lbrace \ell _{i_j}, \ell _{P^{\prime }}\rbrace , r_{i_j}]\) . If \(\ell _P\gt r_P\) , then we drop \(g\) from further consideration, since it does not contribute to the support of \(P\) . Otherwise, we will compute the DP-array with that of \(g\) for \(P^{\prime }\) , using the algorithm of Section 3.3, to obtain \(p_g(P)\) .

Pattern Anti-monotonicity. Suppose pattern \(P^{\prime }\) is a sub-pattern of pattern \(P\) , then we have the following conclusions: (1) For any row \(g\) , \(g_p(P)\lt g_p(P^{\prime })\) . This is because in any possible world instantiation of \(g\) , \(P^{\prime }\) must be supported if \(P\) is supported. (2) When pattern frequentness is defined using expected support, if \(P^{\prime }\) is infrequent, then \(P\) must be infrequent. This is because \(\sum _{g\in G}g_p(P)\lt \sum _{g\in G}g_p(P^{\prime })\) . (3) If \(P^{\prime }\) is probabilistically infrequent, then \(P\) must be probabilistically infrequent. This is because in any possible world of \(D\) , the support of \(P^{\prime }\) is at least the support of \(P\) .

Frequentness Checking. Since our mining algorithms check patterns by pattern-growth, when checking pattern \(P=(t_{i_1}\lt \cdots \lt t_{i_j})\) , we make use of the information of the rows for computation when checking \(P^{\prime }=(t_{i_1}\lt \cdots \lt t_{i_{j-1}})\) .

Given a pattern \(P\) , let \(\mathbf {G_{P}}\) be the set of rows \(\mathbf {g\in G}\) with \(\mathbf {p_g(P)\gt 0}\) , where each row is associated with its valid interval, split point list and DP-array.

The expected-support-based frequentness checking of pattern \(P\) is described by Algorithm 2. In Line 2, we prune \(P\) using the fact that \(p_g(P)\le 1\) and thus \(\sum _{g\in G_{P}}p_g(P)\le |G_{P}|\) . Furthermore, we maintain two variables \(sum\) and \(sum_0\) . Let \(C\) be the set of rows already processed, then \(sum =\sum _{g\in C}p_g(P)\) and \(sum_0=\sum _{g\in (G_{P}-C)}p_g(P^{\prime })\) . Line 7 prunes \(P\) using the fact that \(\sum _{g\in G_{P}}p_g(P)=sum + \sum _{g\in (G_{P}-C)}p_g(P)\le sum+sum_0\) .

Algorithm 3 checks the p-frequentness of pattern \(P\) , where the pruning rules described in Section 4.3 are used: Line 2 applies Count-prune, Lines 7–8 apply Markov-prune, and Line 9 applies Exponential-prune.

Mining by Prefix-projection. Our first mining algorithm is based on the idea of prefix-projection used by sequential pattern mining [37], where the current pattern \(P=(t_{i_1}\lt \cdots \lt t_{i_j})\) is processed using the projected row set \(G_{P^{\prime }}\) of \(P^{\prime }=(t_{i_1}\lt \cdots \lt t_{i_{j-1}})\) .

The recursive algorithm, denoted as \(DFS(P, G_{P^{\prime }})\) ,¹ finds all the frequent patterns with prefix \(P\) from \(G_{P^{\prime }}\) . Specifically, \(DFS(P, G_{P^{\prime }})\) first checks whether \(P\) is frequent using \(G_{P^{\prime }}\) (Algorithm 2 or 3). If so, then it recursively calls \(DFS(t_{i_1}\lt \cdots \lt t_{i_j}\lt t, G_{P})\) for all \(t\in T-\lbrace t_{i_1},\ldots ,t_{i_j}\rbrace\) . The mining algorithm starts by calling \(DFS(\emptyset , G)\) .

For each recursion that processes \(P\) , \(G_{P}\) is maintained so the subsequent recursions for \(t_{i_1}\lt \cdots \lt t_{i_j}\lt t\) may use it. The maintenance of \(G_{P}\) incurs a space cost of \(O(n\cdot j^2)\) , since a DP-array of size \(j\times s=O(j^2)\) is maintained for each row in \(G_{P}\) , and the space can only be released after its corresponding recursion is done.

Since our algorithm works in a depth-first manner, at most \(m\) projected row sets are maintained at any time, with the total space cost \(\sum _{j=1}^{m} O(n\cdot j^2)=O(n\cdot m^3)\) .

For a pattern \(P\) of length \(\ell\) , computing \(p_g(P)\) for all rows \(g\) using the pay-as-you-go algorithm of Section 3.3 takes \(O(n\ell ^2)\) time; then frequency checking takes \(O(n)\) (respectively, \(O(n\log ^2 n)\) ) for expected support (respectively, p-frequentness with FFT computation). If \(C\) patterns are checked, then the time cost is bounded by \(C\cdot O(n\ell ^2 + n) = O(Cnm^2)\) (respectively, \(C\cdot O(n\ell ^2 + n\log ^2 n) = O(Cn(m^2+\log ^2 n)\) ). Additionally, row preprocessing of Section 3.1 takes \(O(nm\log m)\) time.

We remark that the above analysis is very loose, and \(m\) can be replaced by \(\ell _{max}\) , the length of the longest pattern checked.

The benefit of this algorithm is that, for any pattern \(P\) checked, its DP-arrays (of the rows in \(G_{P}\) ) is computed exactly once; however, it can be expensive when the column set \(T\) is large, since the recursions have to be called for all \(t\in T-\lbrace t_{i_1},\ldots ,t_{i_j}\rbrace\) and column candidate pruning is lacking.

Mining by Apriori. The Apriori algorithm works as follows: Starting with the set of all length-1 patterns, we construct length- \(j\) frequent patterns from the set of all length- \((j-1)\) frequent patterns, until there is no frequent pattern left.

We organize the set of length- \(j\) frequent patterns using a prefix-tree \(T_j\) , like the FP-tree proposed in Reference [15]. When computing length- \(j\) frequent patterns \(P=(t_{i_1}\lt \cdots \lt t_{i_j})\) , we need \(G_{P^{\prime }}\) for \(P^{\prime }=(t_{i_1}\lt \cdots \lt t_{i_{j-1}})\) .

However, the space cost is prohibitive if we maintain \(G_{P^{\prime }}\) for each frequent length- \((j-1)\) pattern \(P^{\prime }\) due to the breadth-first search order. Therefore, we choose to recompute \(G_{P^{\prime }}\) over the prefix-tree \(T_{j-1}\) in a depth-first manner (like the prefix-projection method without pattern pruning), and when recursing to the last tree-node \(t_{i_{j-1}}\) , we check whether \(P=(t_{i_1}\lt \cdots \lt t_{i_{j-1}}\lt t_{i_j})\) is frequent for all possible column candidates \(t_{i_j}\in C\) (we will discuss how to compute the candidate set \(C\) later), and if \(P\) is checked to be frequent, then it is inserted to the new prefix-tree \(T_j\) .

Unlike the prefix-projection method, we need to recompute \(G_{P}\) when checking all prefix-trees \(T_z\) ( \(z\ge j\) ). However, the computation in the same tree is still shared. For example, \(G_{t_1\lt t_2}\) is used for computing both \(G_{t_1\lt t_2\lt t3}\) and \(G_{t_1\lt t_2\lt t_4}\) .

We now consider how to determine \(t_{i_j}\) ’s candidate set \(C\) .

Thus, we choose \(C\) to include all child nodes of node \(t_{i_{j-2}}\) in \(T_{j-1}\) . Compared with the prefix-projection method, \(C\) is much smaller than \((T-\lbrace t_{i_1},\ldots ,t_{i_{j-1}}\rbrace)\) and thus the recursion has a much smaller fan-out.

The pattern anti-monotonicity also allows for the following pattern pruning:

Theorem 5.2 is efficient to check and is thus examined before the more expensive frequentness checking of Algorithms 2 or 3, which requires computing the DP-array for all rows in \(G_{P}\) .

In the Apriori algorithm, computing \(p_g(P)\) of patterns of length \(\ell\) requires re-computing the DP-array for frequent patterns of length \(\lt \ell\) ; this brings an additional factor of \(\ell _{max}\) to our prefix-projection algorithm’s time complexity, since each \(P\) is re-processed for at most \(\ell _{max}\) times. In return, the algorithm enjoys a small recursion fan-out, which, however, is not reflected in time complexity analysis.

Datasets. As Table 2 shows, five publicly available real datasets were used in our experiments, including three microarray gene expression datasets, a movie rating dataset, and an RFID user trace dataset.

Among the three biology datasets, GDS2712 and GDS2002² are microarray datasets of the baker’s yeast Saccharomyces cerevisiae from the Gene Expression Omnibus (GEO) database [3], where each matrix entry has three replicates. We also tested some other datasets such as GDS2713, GDS2715, and GDS2003, and we found that the results are similar and hence omitted to avoid redundancy in presentation. For GDS2712, the 7,826 rows shown in Table 2 are obtained from preprocessing, where we removed those unidentified genes and control probes, since they cannot be used to calculate the biological significance p-value against the database of ground-truth gene functional categories.

The other biology dataset, GAL,³ is regarding yeast galactose utilization [32], which was also used by References [12, 33]. As Table 2 shows, GAL contains 205 gene probes (rows) and 20 experimental conditions (columns) with four replicates for each entry in the matrix.

For the gene expression microarray datasets, we construct a uniform value interval for every matrix entry as \([min, max]\) , where \(min\) (respectively, \(max\) ) is the minimum (respectively, maximum) replicated value of the entry from repeated experiments. As for Gaussian value distribution, given a matrix entry with replicates, we compute their sample mean \(\mu\) and sample variance \(\sigma\) and assume that the matrix entry has a value following \(\mathcal {N}(\mu , \sigma ^2)\) . We fit its PDF with a cubic spline as we described in Section 3.4, and we only consider the six intervals \([\mu -3\sigma , \mu -2\sigma ]\) , \([\mu -2\sigma , \mu -\sigma ]\) , \([\mu -\sigma , 0]\) , \([0, \mu +\sigma ]\) , \([\mu +\sigma , \mu +2\sigma ]\) , \([\mu +2\sigma , \mu +3\sigma ]\) when computing \(p_g(P)\) . We similarly fit the PDF of exponential value distributions following our discussion in Section 3.4.

Besides the gene expression datasets, we also used a movie rating dataset MovieLens⁴ with 100,000 ratings from 943 users on 1,682 movies [16] to evaluate the algorithms. The movie rating dataset is an incomplete matrix, and we used TensorFlow to conduct factorization-based matrix completion, which provides a complete 943 \(\times\) 1,682 user rating matrix \(M\) , where \(M_{ij}\) estimates User \(i\) ’s rating towards Movie \(j\) . The loss function of factorization is given bywhere \(M^0\) is the observed rating matrix with missing values, \(M=UV\) is the estimated rating matrix, \(U\in \mathbb {R}^{943\times 10}\) , \(V\in \mathbb {R}^{10\times 100}\) , \(|X|\) takes the element-wise absolute value, reduce_sum \((X)\) sums all elements of \(X\) , and \(\mathcal {P}_{\Omega }(X)\) projects to a matrix where the \((i, j)\) -th element equals \(X_{ij}\) if matrix entry \((i, j)\) in \(M^0\) is observed, and 0 otherwise. Put simply, the goal is to minimize the difference between \(M^0\) and \(M\) for those elements that are observed in \(M^0\) . Implementation in TensorFlow is simple, with operators like tf.abs for taking the absolute value, tf.subtract for matrix subtraction, tf.gather to collect elements at observed matrix locations, and tf.reduce_sum to compute the sum of elements. We learn \(U\) and \(V\) by initializing them with truncated normal sampling, running stochastic gradient descent with a learning rate of \(10^{-3}\) and a staircase decay of rate 0.96, and running for 100,000 steps.

To consider only the popular movies, we select the top-100 movies that get the most user ratings, which generates a \(943\times 100\) submatrix \(M_s\) of \(M\) . Since each new rating \(r\) is now a low-rank approximation, we introduced uncertainty to each rating \(r\) . A rating in the original data takes its value from \(\lbrace 1,2,3,4,5\rbrace\) ; after matrix completion, a rating is a real value like \(r=3.4\) , in which case, we assign an interval \([3, 4]\) to the matrix entry, and for OPSMRM, we assign replicates \(\lbrace 3, 4\rbrace\) . Other reasonable methods to assign uncertainty to the values can also be used. Note that user ratings are intrinsically uncertain: (1) the discrete rating values \(\lbrace 1,2,3,4,5\rbrace\) are coarse-grained and preferences among the movies with the same rating are not captured; (2) different users have different bias on the rating scale, with some giving rating 5 only to movies they like while others giving rating 4 also to such movies. Such uncertainty should be captured by relaxing the strict values in \(\lbrace 1,2,3,4,5\rbrace\) .

Finally, we also prepared an RFID dataset,⁵ which contains the traces of 12 voluntary users in a building with 195 RFID antennas installed. Only 80 antennas are active, i.e., they detect the trace of at least one user for at least one time. The trace of a user consists of a series of (period, antenna) pairs, which indicates the time period during which the user is continuously detected by the antenna. When a user is detected by the same antenna at different time periods, we split the raw trace into several new traces in which the user is detected by any antenna for at most one time. We then generate a matrix having 105 rows (i.e., user traces) and 80 columns (i.e., antennas), where the \((i, j)\) -th element records the potential interval of time that User \(i\) is near Antenna \(j\) , and an OPSM thus records a group of users that visits a group of locations (captured by antennas) in the same time order. To run OPSMRM for comparison, the set of timestamps within each time period is enumerated to generate the corresponding set-valued matrix.

Evaluation Metrics. The following metrics are studied to evaluate result quality, time efficiency, and space efficiency, respectively:

(1) Result Quality: For the microarray gene expression datasets, we use the biological significance of the mined OPSMs to demonstrate the result quality of our proposed method. We adopt a widely used metric, p-value [12, 21, 33], which measures the association between OPSMs mined and the known gene functional categories. Specifically, a smaller p-value indicates a stronger association between an OPSM and gene categories, i.e., biologically more significant. Following Reference [12], we also consider four exponential-scale p-value ranges as significance levels, such as \([0,10^{{-}40})\) , \([10^{{-}40}, 10^{{-}30})\) , \([10^{{-}30}, 10^{{-}20})\) , and \([10^{{-}20}, \infty)\) (the actual ranges depend on the concrete datasets). To compare the result quality of our algorithms with those of the existing algorithms, we use the number and proportion of OPSMs mined at each significance level.

For the movie rating dataset, we define a Kendall tau score (KTS) motivated by the concept of Kendall tau distance: For a mined OPSM with movie order \(t_{1}\lt t_{2}\lt \cdots \lt t_{k}\) , we compute a KTS for each user \(g_i\) in the OPSM, which equals the fraction of all possible \(C_k^2\) movie pairs \((t_i, t_j)\) where the rating order is consistent with that in our \(943\times 100\) submatrix \(M_s\) ; the KTS of the OPSM is then computed as the average KTS over all its users.

For the RFID dataset, each mined OPSM consists of a set of users sharing a common subroute passing a subset of antennas. We define a measure called the trace matching score (TMS) as follows: For each minded OPSM with subroute \(P=(t_{1}\lt t_{2}\lt \cdots \lt t_{k})\) , we compute a TMS for each user \(g_i\) in the OPSM as follows:where \(\mathcal {T}_{g_i}\) denotes the manually annotated ground-truth trace of \(g_i\) , which is also provided by the RFID dataset. The TMS of the OPSM is then computed as the average TMS over all its users.

(2) Time Efficiency: We evaluate the time efficiency of the algorithms by considering the following two metrics:

Here, SP-Time equals TR-Time divided by the number of patterns mined.

(3) Space Efficiency: We report the peak memory consumption of each program run in our experiments.

All experiments were repeated for 10 times and the reported metrics were averaged over the 10 runs (although the results observed from different runs are actually quite stable/similar).

For the purpose of succinct presentation, we abuse the term “more than” (similarly for “less than” and “fewer than”): When we say that A is \(k\) times more than B, we mean \(A=k\cdot B\) , not \(A=(1+k)\cdot B\) . The terms “more,” “less,” and “fewer” are meant to indicate the directions of comparison.

As an overview of our findings, our ES and PF algorithms are robust and output higher-quality OPSMs than OPSMRM and POPSM, and generally find more OPSMs. Our algorithms also consistently outperform OPSMRM and POPSM in terms of SP-Time in the context of uniform value distributions. ES consumes the least amount of memory, and PF consumes more memory than ES but generally less memory than OPSMRM and POPSM. PF performs the best w.r.t. result quality and quantity. Among other findings, using Gaussian value distributions as the uncertain model generally outperforms the uniform distribution model scheme, but the running time is two orders of magnitude longer; using exponential value distributions exhibits a similar performance. Finally, our parallel algorithm for prefix-projection-based mining scales well with the number of mining threads.

The rest of this section is organized as follows: Section 7.2 reports the results on GDS2712 and GDS2002 in terms of the result quality and quantity, running time, and peak memory cost. There, we use GDS2712 to demonstrate that our algorithms are efficient even when exponential value distribution is adopted. Due to space limitation, we report the results on GAL in Appendix C. Section 7.3 then uses the GDS2712 dataset to study the effect of parameters \(\tau _{cut}\) and \(\tau _{prob}\) , and Section 7.4 compares row selection strategies to justify the use of \(\tau _{cut}\) . Subsequently, Section 7.5 reports the vertical scalability experiments of our parallel prefix-projection-based mining algorithm, with additional scalability results in Appendix D. Section 7.6 reports our results on MovieLens, and Section 7.7 replicates the MovieLens to generate datasets of different sizes to study the algorithm scalability. Finally, Section 7.8 reports the results on the RFID user trace dataset.

Size Thresholds \(\tau _{row}\) and \(\tau _{col}\) for Testing. Both datasets are used in Reference [33] for evaluation, and following Reference [33], we fix \(\tau _{cut} = 0.6\) for them. In fact, we also tested the other GDS datasets in Reference [33] such as GDS2713, GDS2715, and GDS2003, and the results are very similar and thus omitted. Since the results are similar for various \(\tau _{prob}\) values, we only show results when \(\tau _{prob} = 0.5\) to save space (see Section 7.3 for more results on the effect of \(\tau _{prob}\) ). While many \((\tau _{row}, \tau _{col})\) pairs are possible, we would like to show a limited number of typical combinations so they can fit in one figure to be readable. For GDS2712, we set \(\tau _{col} = 4\) and vary \(\tau _{row}\) among \(\lbrace 400, 500, 600, 700, 800\rbrace\) to show the effect \(\tau _{row}\) . For GDS2002, we vary \(\tau _{row}\) among \(\lbrace 200, 300, 400\rbrace\) but also vary \(\tau _{col}\) among \(\lbrace 2, 3\rbrace\) to show the effect of \(\tau _{col}\) .

Result Quality. Figure 8 presents the fraction of mined OPSMs that fall in each significance level for GDS2712. We see that our algorithms find larger fractions of high-quality OPSMs than POPSM and OPSMRM, as they have taller white bars (representing the highest significance level with p-value \(\in [0, 10^{-20})\) ). For example, when \((\tau _{row}, \tau _{col}) = (800, 4)\) , our proposed ES(U), PF(U), PF-G(U), and PF-P(U) all have \(90\%\) patterns falling into the highest significance level, while POPSM and OPSMRM have only \(50\%\) . Moreover, our proposed ES(U), PF(U), PF-G(U), and PF-P(U) discovered a significantly larger number of OPSM patterns compared with POPSM and OPSMRM as shown in Table 3. We can see that our algorithms using uniform value distributions generate consistently more OPSMs than POPSM and OPSMRM (and of higher quality as well). Furthermore, if we look at the numbers of OPSMs falling in the highest significance level, then when \((\tau _{row}, \tau _{col}) = (800, 4)\) , our proposed ES(U), PF(U), PF-G(U), and PF-P(U) have 60, 60, 55, and 60 patterns falling into the highest significance level, respectively, while POPSM and OPSMRM have only 5 and 1, respectively.

As for our algorithms using Gaussian value distributions, while Figure 8 shows that they do not always have a taller white bar than the uniform distribution variants (e.g., much taller when \((\tau _{row}, \tau _{col}) = (400, 4)\) but shorter when \((\tau _{row}, \tau _{col}) = (800, 4)\) ), the absolute OPSM counts found in each significance level dominates those of the uniform distribution variants, as we can see from Table 3, showing that using the more expensive Gaussian distribution algorithm variants does pay off in result quality.

Figure 9 presents the fraction of mined OPSMs that fall in each significance level for GDS2712, similar to Figure 8 except that we replace those algorithms adopting Gaussian value distributions with those adopting exponential value distributions. We can see that using exponential value distributions leads to much shorter white bars than using uniform distributions when \((\tau _{row}, \tau _{col}) = (800, 4)\) , even though the absolute result counts are consistently higher, as shown in Table 3.

Comparing Figure 8 with Figure 9, we can see that using Gaussian value distributions is a better assumption than using exponential value distributions, since the former has a taller white bar. The same holds true when comparing absolute result counts as shown in Table 3, where using Gaussian value distributions consistently leads to more results for every \((\tau _{row}, \tau _{col})\) setting.

For the future experiments on the other datasets, we only consider uniform and Gaussian value distributions, since using exponential value distributions leads to a lower result quality than using Gaussian value distributions, while both are much more expensive than using uniform value distributions due to the need of applying the expensive cubic spline approach.

As Table 3 shows, among our algorithms, PF performs better than ES, as well as PF-G and PF-P. This verifies that considering PMF gives better results than considering only expectation. Also note that PF-G and PF-P produce results not far behind PF, and in fact, PF-P(G) gives the same results as PF(G), and PF-P(E) gives the same results as PF(E), showing that our approximations are accurate.

As for GDS2002, Figure 10 presents the fraction of mined OPSMs that fall in each significance level, while Table 4 shows the absolute counts. We observe similar results as on GDS2712. Specifically, Figure 10 shows that our algorithms find larger fractions of high-quality OPSMs than POPSM and OPSMRM, as they have taller white bars. The absolute counts in all the significance levels are also much higher, as Table 4 shows, with the absolute counts found by the Gaussian distribution algorithm variants in each significance level dominating those of the uniform distribution variants. This shows that using the more expensive Gaussian distribution algorithm variants does pay off in result quality. The comparative performances of the various algorithms also exhibit the same observations as those on GDS2712. Time Efficiency. We hereby report the experimental results on time efficiency for GDS2712. The time results on GDS2002 are very similar and thus omitted.

Figures 11 and 12 show the total running time TR-Time and average single-pattern running time SP-Time, respectively, of our proposed methods adopting uniform value distributions, as well as the existing algorithms OPSMRM and POPSM for different size thresholds \((\tau _{row}, \tau _{col})\) . From Figure 11, we can see that POPSM runs faster than our methods, and OPSMRM has comparable total running time with our methods, but recall that both POPSM and OPSMRM mined significantly less OPSM patterns than our algorithms and their results are of poorer quality w.r.t. all different size thresholds \((\tau _{row}, \tau _{col})\) . For example, as Table 3 shows, when \((\tau _{row}, \tau _{col}) = (800, 4)\) , our ES(U) discovered 33.5 times more OPSMs than OPSMRM and 6.7 times more than POPSM.

According to Figure 12, OPSMRM has the worst SP-Time, and our methods have comparable average single-pattern running time (SP-time) with POPSM. Among our methods, PF-G and PF-P algorithms are always faster than PF, and in some cases the DFS algorithms run faster than our Apri algorithms: For example, DFS-ES runs faster than Apri-ES. This is because the GDS2712 dataset is small with merely 7 columns, so many of our pruning techniques are not effective. As we shall see later, on the GAL dataset with 20 columns and MovieLens dataset with 100 columns, Apriori-based algorithms are faster than DFS-based ones.

Figures 13 and 14 show the total running time TR-Time and average single-pattern running time SP-Time, respectively, of our proposed methods adopting Gaussian value distributions. We can see that the comparative performances of our algorithms are very similar to those in Figures 11 and 12 in the context of uniform value distributions. However, algorithms using Gaussian value distributions are two orders of magnitude more expensive.

Figures 15 and 16 show the total running time TR-Time and average single-pattern running time SP-Time, respectively, of our proposed methods adopting exponential value distributions. We can see that the performances of our algorithms are very similar to those in Figures 13 and 14 in the context of Gaussian value distributions, actually slightly more expensive but comparable. This shows that our cubic spline approach is efficient and can handle various value distributions.

Memory Efficiency. We hereby report the experimental results on memory efficiency for GDS2712. The memory cost results on GDS2002 are very similar and thus omitted.

Figure 17 shows the peak memory usage of various algorithms using uniform value distributions. We can see that our algorithms have a consistently lower peak memory consumption than OPSMRM. As an illustration, when \((\tau _{row}, \tau _{col})=(400, 4)\) , our Apri-ES algorithm consumes 6.6 times less memory than OPSMRM. Our less memory usage does not compromise the number of mined OPSMs, as Table 3 already indicates. our algorithms discovered many times more OPSMs than OPSMRM and POPSM. For example, when \((\tau _{row}, \tau _{col})=(400, 4)\) , our DFS/Apri-ES algorithms discovered 2.5 times more OPSMs than OPSMRM and 3.7 times more OPSMs than POPSM.

Our DFS/Apri-PF algorithms consume more memory than DFS/Apri-ES due to the need of processing PMF vectors (cf. Section 4.2), but they generate higher-quality OPSMs (in terms of p-value) than ES (cf. Figure 8). However, even DFS/Apri-PF algorithms have a much lower memory consumption compared with OPSMRM. While POPSM tends to consume slightly less memory than our algorithms, it outputs much fewer OPSMs.

Figure 18 shows the peak memory usage of our algorithms using Gaussian value distributions, and compared with Figure 17, we can see that their memory cost is one order of magnitude higher. This is reasonable, because each matrix entry now introduces seven split points (cf. Figure 6) rather than two as in the uniform interval case, and the recursive integral computation as specified in Equation (11) now requires the maintenance of a high-order polynomial for probability computation.

Figure 19 shows the peak memory usage of our algorithms using exponential value distributions, and compared with Figure 18, we can see that the memory cost is very similar.

We next report the experiments on the effect of \(\tau _{cut}\) and \(\tau _{prob}\) using the GDS2712 dataset. The results on the other datasets are similar and thus omitted.

Effect of Inclusion Threshold \(\tau _{cut}\) . We fix \(\tau _{row} = 400\) and \(\tau _{col}=3\) , since various methods generate comparable number of OPSMs with these parameters. We also fix \(\tau _{prob} = 0.5\) for DFS-PF and Apri-PF, but vary \(\tau _{cut}\) among \(\lbrace 0.2,0.3,0.4,0.5,0.6\rbrace\) . Figure 20 shows the fraction of OPSMs in each significance level for every algorithm with different inclusion threshold \(\tau _{cut}\) . Table 5 presents the number of OPSMs mined with different inclusion threshold \(\tau _{cut}\) , where we can see that our algorithms consistently find more OPSMs than POPSM and OPSMRM at various values of \(\tau _{cut}\) in all significance levels, and those with Gaussian value distributions are better than those with uniform value distributions. In terms of the distribution of OPSMs among various significance levels, we can see from Figure 20 our algorithms consistently have a tall white bar (with those using Gaussian distributions being taller), while OPSMRM is only competitive when \(\tau _{cut}\) is very low (however, the absolution count is much lower, as shown in Table 5).

Effect of Probability Threshold \(\tau _{prob}\) . We study the effect of \(\tau _{prob}\) used in PF by fixing \((\tau _{row} = 600, \tau _{col}=4)\) and \(\tau _{cut}=0.2\) ; and varying \(\tau _{prob}\) from 0.2 to 0.95 with a step length of 0.05.

First consider PF(U) algorithms. Table 6 shows the number of OPSMs mined at each value of \(\tau _{prob}\) , and Figure 21 shows their distribution in different significance levels. We can see that our approach is very robust, not sensitive to parameter \(\tau _{prob}\) , and has consistent performance in terms of result quality. Figure 22 shows the TR-Time of PF, PF-G, and PF-P algorithms w.r.t. diffierent \(\tau _{prob}\) . We can see that our DFS-based algorithms are faster than Apriori-based ones, which is often the case when the number of results are small. However, we tested that when the number of results reach thousands, our Apriori-based algorithms are actually faster, since the additional pruning power of the Apriori-based algorithms outweighs the pruning cost. Recall, for example, Figures 36 and 37 for GAL. So, both algorithms have their merits.

We also tested PF(G) algorithms, with Table 7 showing the number of OPSMs mined at each value of \(\tau _{prob}\) , Figure 23 showing their distribution in different significance levels, and Figure 24 showing the TR-Time of PF, PF-G, and PF-P algorithms w.r.t. different \(\tau _{prob}\) . We can observe similar results, except that the TR-Time is two orders of magnitude slower (compare Figure 24 with Figure 22), and the result quality is higher, as can been seen from the taller white bars in Figure 23 than in Figure 21, and the larger number of results in Table 7 than in Table 6.

Recall that after (Step 1) mining patterns \(P=(t_{i_1}, t_{i_2}, \ldots , t_{i_\ell })\) that are deemed significant, we (Step 2) select those rows whose supporting probabilities are at least \(\tau _{cut}\) into the OPSM for each pattern \(P\) .

We remark that our row selection strategy in Step 2 is reasonable, since \(\tau _{cut}\) directly reflects the confidence level of supporting probability, which is intuitive to end-users. An alternative approach is to select \(k\) rows whose supporting probabilities are the highest, but it is difficult to ask end-users to set a proper \(k\) . One way is to set \(k\) as \(\tau _{row}\) , but this tends to be an underestimate, since each of the \(k\) rows has \(p_g(P)\lt 100\%\) while \(\tau _{row}\) is defined with respect to each deterministic possible world.

To compare these two strategies, Figure 25 shows our previous experiments on GDS2712 assuming uniform value distributions, where algorithms ending with “(U)” are those that select rows with \(p_g\ge \tau _{cut}=0.6\) into the OPSMs, and algorithms ending with “(k)” are those that select rows with \(k=\tau _{row}\) highest values of \(p_g\) into the OPSMs. We can see that the latter row selection method consistently delivers fewer highly significant OPSMs, as the white bars are much shorter. Also, Figure 26 shows our previous experiments on GDS2712 assuming Gaussian value distributions, and we can obtain the same observation that selecting rows with \(k=\tau _{row}\) highest values of \(p_g\) is inferior to our adopted strategy of selecting rows by \(\tau _{cut}\) in terms of result quality.

Recall that Section 6 presents a parallel version of our prefix-projection-based mining algorithm. We now explore its vertical scalability by varying the number of mining threads to be 1, 2, 4, 8, 16, 32, and 64. Figure 27 shows the scalability curve for our four algorithms with uniform value distributions on GDS2712 when \((\tau _{row}, \tau _{col})=(400, 4)\) . Figure 28 shows the scalability curve for our four algorithms with uniform value distributions on GDS2002 when \((\tau _{row}, \tau _{col})=(200, 2)\) . We can see that the running time clearly reduces as we increase the number of mining threads. For example, in Figure 28, our parallel PF takes 56.27 seconds when mining with one thread, but the time reduces to 3.80 seconds when mining with 32 threads, and the time reduces to 2.89 seconds when mining with 64 threads. Appendix D reports the results for our algorithms with Gaussian value distributions, where we observe even better speedups. We also report the speedups of our parallel algorithm over the other datasets in Appendix D, as well as the peak memory usage results.

Recall from Section 7.1 that we use a completed movie rating matrix to find OPSMs where users share the same preference order over a set of movies, and that a Kendall tau score (KTS) is defined to judge how well an OPSM pattern is reflected in its users’ movie ratings. We now report the results on the movie rating dataset where we added uncertain intervals to the discrete scale ratings.

Time Efficiency. Table 8 shows the TR-Time of various algorithms with different \((\tau _{row}, \tau _{col})\) (without loss of generality, we fix \(\tau _{cut} = 0.3\) and \(\tau _{prob} = 0.6\) ), where OOM means out of memory. Note that when \(\tau _{row} = 200\) , POPSM always runs out of memory after 10 hours. OPSMRM has the shortest TR-Time but it failed to produce any OPSMs that satisfy the \(\tau _{col}\) threshold. We see that our methods are tens of times faster than POPSM, and ES methods are faster than PF, as they do not process PMF vectors. Also, our Apriori-based algorithms are faster than their DFS-based counterparts due to the effective pattern pruning. Finally, approximation algorithms PF-P and PF-G provide reasonable speedup to PF.

Effectiveness. Without loss of generality, we consider \(\tau _{cut} = 0.3\) , \(\tau _{prob} = 0.6\) , \(\tau _{row} = 300\) , and vary \(\tau _{col}=5\) and visualize the top-10 OPSMs mined by DFS/Apri-ES, DFS/Apri-PF, and POPSM with the highest KTS. We remark that in this set of experiments, our DFS-ES consumes 2,582 \(\times\) less memory than POPSM and 3,294 \(\times\) less memory than OPSMRM.

Table 9(a) lists the top-10 OPSM patterns with the highest KTS produced by ES, Table 9(b) by PF, Table 9(c) by PF-P, Table 9(d) by PF-G, and Table 9(e) by POPSM. To be space-efficient, we use abbreviations for movie names, and their meanings are listed in Table 10. We can see that the OPSMs mined by our methods generally have a higher KTS.

Memory Efficiency. Table 11 shows the peak memory usage of our algorithms, POPSM and OPSMRM w.r.t. size thresholds \((\tau _{row}, \tau _{col})\) , where OOM represents Out of Memory. When \((\tau _{row}, \tau _{col}) = (200, 4)\) or \((200, 5)\) , POPSM ran out of memory after running for 10 hours (or more accurately, 603 minutes). It is clear that our Expected Support (ES)-based mining algorithms, namely, Apri-ES and DFS-ES, have a consistently much lower peak memory consumption than all the other algorithms. As an illustration, when \((\tau _{row}, \tau _{col})=(300, 5)\) , DFS-ES consumes 2,582 times less memory than POPSM and 3,294 times less memory than OPSMRM, and our Apri-ES consumes 720 times less memory than POPSM and 918 times less memory than OPSMRM. Also, our PMF approximation-based PF algorithms also use memory comparable to (or more accurately, slightly larger than) the ES counterparts, and the memory usage is three orders of magnitude less than their exact PF counterparts. This shows that our PMF approximation allows PF-based mining algorithms to scale to much larger datasets given the same memory budget, even though the speedup ratio is not as significant.

To examine how well our algorithms scale to the number of rows and columns of a data matrix \(D\) , we duplicate the rows and columns of our \(943\times 100\) movie rating submatrix \(M_s\) to generate larger datasets for running scalability experiments.

To test row scalability, we duplicated the rows of \(M_s\) for 100, 200, 300, 400, and 500 times, and ran the various algorithms on them. Without loss of generality, we set \(\tau _{row} = 0.6 * n\) (where \(n\) is the row number) and fix \(\tau _{col} = 5\) and \(\tau _{prob} = 0.6\) . Unfortunately, OPSMRM runs out of memory even on the smallest data with \(M_s\) ’s rows duplicated for 100 times, so we cannot report its result. The results for the other algorithms are shown in Figure 29, where we observe that POPSM is much slower than our algorithms, and that Apri-ES, Apri-PF-P, and Apri-PF-G are much faster than the other algorithms. Overall, the time of all algorithms scale linearly with row number \(n\) , which matches our derived time complexity \(O(Cn(m^2+\log ^2 n))\) in Section 5.2 (i.e., almost linear to \(n\) and quadratic to \(m\) ).

To test column scalability, we duplicated the columns of \(M_s\) for 2, 4, 8, and 16 times and ran the various algorithms on them. Here, we set \(\tau _{row} = 700, \tau _{col} = 5\) , and \(\tau _{prob} = 0.6\) . The results are shown in Figure 30, where we observe that the running time of various algorithms increase quickly with column number \(m\) , which aligns with our analysis. Although OPSMRM has a similar running time to our algorithms, the memory consumption rises quickly with more columns. Also, OPSMRM runs out of memory when the columns are duplicated for four times and thus, we only plot two points for it in Figure 30. Finally, POPSM is not only slower than all our algorithms, its curve slope is also steeper.

Recall from Section 7.1 that we use an RFID user trace dataset to find OPSM patterns that reveal a group of users who travel to a sequence of locations in the same time order, and that a trace matching score (TMS) is defined to judge how well an OPSM pattern is reflected in its users’ actual travel trajectories that are manually annotated. We now report the results on the RFID user trace dataset where uniform time intervals are imposed at each visit location (i.e., antenna). We also compare our methods with POPSM and OPSMRM. However, OPSMRM always runs out of memory in this set of experiments and is thus not reported.

In this set of experiments, we use \(\tau _{prob}=0.8\) and \(\tau _{cut}=0.6\) for PF algorithms. Figure 31 shows the fraction of OPSMs in each significance level for our algorithms and POPSM with different \((\tau _{row}, \tau _{col})\) , and Table 12 presents the number of OPSMs mined. We can see that POPSM mines very few patterns compared with our algorithms, and all OPSMs are with TMS in the range of \([0.8, 0.9)\) , while our algorithms are able to find OPSMs with TMS \(\gt 0.9\) (i.e., black bars). Moreover, PF-G still achieves performance close to that of PF even when PF-P fails to get close. This verifies that PF-G is more robust than PF-P, thanks to its use of both the expected support and the standard deviation (i.e., \(\mu\) and \(\delta\) ) while FP-P only uses the expected support, as we have indicated at the end of Section 4.4.2.

Figures 32 and 33 show the total running time TR-Time and average single-pattern running time SP-Time, respectively, of our algorithms and POPSM for different size thresholds \((\tau _{row}, \tau _{col})\) . From Figure 32, we can see that POPSM runs faster than our methods, but recall that POPSM mined significantly less OPSM patterns than our algorithms and their results are of poorer quality. Also, PF-P algorithms are also fast but recall that the number of OPSMs found is much smaller than our other algorithms. According to Figure 33, POPSM has the worst SP-Time, followed by PF-P algorithms, due to their small number of OPSMs found.

Figure 34 shows the peak memory usage of our algorithms and POPSM. We can see that the memory cost of PF algorithms is much higher than the other algorithms. In this set of experiments ES and PF-G have similar result quantity and quality, showing that ES is effective in certain datasets, and PF-G is a good approximation method that reduces the cost significantly while retaining similar results as with PF.

Let us define \(x_{n+1}=r\) . Then, if we can prove that the following equation holds for any positive integer \(i\) then Theorem 3.1 is proved by setting \(i=n\) .

We prove Equation (23) by induction:

Proof. Suppose that the rows in set \(S\) is divided into two sets \(S_1\) and \(S_2\) . It is sufficient to prove that, when \(P\) is p-frequent in \(S_1\) , it is also p-frequent in \(S\) .

Let us denote the support of \(P\) in a set \(S\) by \(X^S\) , and denote the PMF (and CDF) of \(X^S\) by \(f^S(c)\) (and \(F^S(c)\) ). When \(P\) is p-frequent in \(S_1\) , according to Equation (2), we have

According to Equation (24), \(F^{S_1}(\tau _{row}-1)\le 1-\tau _{prob}\) . If we can prove \(F^S(\tau _{row}-1)\le F^{S_1}(\tau _{row}-1)\) , then we are done, since this implies \(F^S(\tau _{row}-1)\le F^{S_1}(\tau _{row}-1)\le 1-\tau _{prob}\) , or equivalently, \(Pr\lbrace X^S\ge \tau _{row}\rbrace =1-F^S(\tau _{row}-1)\ge \tau _{prob}\) .

We now prove \(F^S(\tau _{row}-1)\le F^{S_1}(\tau _{row}-1)\) . Let us denote \(\tau _{row}^{\prime }=\tau _{row}-1\) . Then, we obtain

(1) Count-Prune. Let \(cnt(P)=|\lbrace g\in G\,|\,p_g(P)\gt 0\rbrace |\) , then pattern \(P\) is not p-frequent if \(cnt(P)\lt \tau _{row}\) .

(2) Markov-Prune. Pattern \(P\) is not p-frequent if

(3) Exponential-Prune. Let \(\mu =E(X)\) and \(\delta =\frac{\tau _{row}-\mu -1}{\mu }\) . When \(\delta \gt 0\) , pattern \(P\) is not p-frequent if

The problem of top- \(k\) queries (or ranking queries) has been extensively studied over uncertain databases, in which each tuple is associated with an uncertain score, and the objective is to find top- \(k\) tuples with the largest ranking scores in the probabilistic sense.

Many semantics for top- \(k\) queries on uncertain data are proposed in the literature, and they can be categorized into 3 classes, based on the uncertain data models they use:

The problem of ranking uncertain tuples is relevant to our OPSM mining problem, since both problems consider the order relationships of uncertain numerical values. In fact, the uncertain database considered in the ranking problem is equivalent to a row of the data matrix in our OPSM mining problem. OPSMRM [33] adopts the discrete attribute-level uncertain data model, while we study OPSM mining under the continuous attribute-level uncertain data model.

We focus on the problem of mining frequent patterns (such as itemsets and sequences) over uncertain data in this subsection.

Frequent Itemset Mining. The problem of frequent itemset mining has been extensively studied in the context of uncertain data. Earlier works use expected support to measure pattern frequentness [1, 7]. However, [5, 35] find that the use of expected support may render important patterns missing, and thus, recent research focuses more on using p-frequentness [5, 25]. See the excellent survey of [26] for a more complete overview of the many important studies in this area.

In fact, the transaction database for frequent itemset mining can also be represented as a data matrix \(D\) , where each row corresponds to a transaction and each column corresponds to an item. If transaction \(g\) contains item \(t\) , then \(D[g][t]=1\) ; otherwise, \(D[g][t]=0\) .

Sequential Pattern Mining. The problem of sequential pattern mining has also been studied over uncertain data: Liu et al. [22] evaluate pattern frequentness based on expected support, while Zhao et al. [37] evaluate pattern frequentness based on p-frequentness. Yang et al. [31] study the problem of mining long sequential patterns in a noisy environment.

The problem of OPSM mining can be reduced into a sequential pattern mining problem: each row in the data matrix can be treated as a sequence (called row sequence hereafter) as Figure 1 illustrates, where the elements are the columns \(t_1, \ldots , t_m\) . However, unlike general sequential pattern mining problems, in OPSM mining, each element can appear at most once in a row sequence, and it appears exactly once if there is no missing data.

Another difference between OPSM mining and sequential pattern mining is that, the identities of the supporting row sequences are important in OPSM mining but are immaterial in sequential pattern mining. Recall that an OPSM may refer to a set of coexpressed genes, or a group of users with similar preferences.