MINT: Tensor Decomposition on Stacked Recurrence Matrices for Time Series Data Mining

# MINT: Tensor Decomposition on Stacked Recurrence Matrices for Time Series Data Mining
Source: [https://arxiv.org/html/2608.04157](https://arxiv.org/html/2608.04157)
Kaamil Kaka∗,2Equal contribution\.University of Texas at Dallas, Correspondence:[kmk230008@utdallas\.edu](https://arxiv.org/html/2608.04157v1/mailto:[email protected])2Neuralix AIEvangelos E\. Papalexakis3Zachary Zimmerman33University of California, Riverside, Department of Computer ScienceVikram Jayaram2

###### Abstract

Recurrence plots are a time series data mining primitive applied to a variety of domains \(e\.g\. star light curves, sound waveforms, CCT telemetry\)\. This work proposes tensorized self\-similarity matrices as a primitive for univariate time series datasets \(N×nN\\times n\) ofNNtime series of lengthnnwith a subsequence window of lengthmm, and whose tensor\-based nature is naturally extensible to multivariate datasets\. The proposed method to compute this primitive computes dot plots of sizeN×\(n−m\+1\)×\(n−m\+1\)N\\times\(n\-m\+1\)\\times\(n\-m\+1\)from these datasets, where the subsequent tensor is mined using tensor decomposition methods to mine for co\-clustered patterns\. We demonstrate our results in mass rapid transit, electricity demand, wind turbine, and car traffic data, finding the MINT pipeline effectively co\-clusters cross\-sensor patterns in highly regular datasets containing motifs at regular intervals\.

## 1Introduction

When the same discordant phenomenon appears across dozens of sensors simultaneously, do we haveNNanomalies or one? Which sensors share not just the time stamp of an event, but itsshape? Which subsequences of a sensor’s history participate in a cross\-series pattern? These questions are notpredictionproblems, butstructure discoveryproblems, requiring interpretable decomposition, not learned representations\.

The Matrix Profile\[[28](https://arxiv.org/html/2608.04157#bib.bib28)\]is a primitive for motif and discord mining, compressing the full self\-similarity structure of a time series into an efficient nearest\-neighbor distance vector\. This compression is lossy; the Mplot\[[22](https://arxiv.org/html/2608.04157#bib.bib22)\], which retains the full pairwise self\-similarity distances between all subsequences, contains richer structure than the Matrix Profile alone\. Prior work demonstrates that Mplots support rich visual analysis of individual time series, enabling practitioners to identify motifs, discords, and phase relationships by eye\. However, Mplots alone are not sufficient as an analytical method\. Visual inspection does not scale with dataset size, and with tens or hundreds of series, there is no accompanying method for extracting cross\-series structure or summarizing temporal alignment across the collection\. By retaining the full Mplot for each series and stacking them into a tensor, we consider the entire dataset simultaneously\. The resulting tensor exhibits low\-rank structure recoverable by interpretable tensor decomposition methods \(e\.g\. CPD, Tucker, etc\.\), corresponding to co\-clustering information that has no natural representation in per\-series analysis\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/x1.png)Figure 1\.1:NMF applied to the synthetic dataset\.Top\)Representative series from each phase group over the first 30 days, gray shading marks the first sawtooth\.Bottom Left\)The temporal factorWTW^\{T\}\(time×\\timescomponents\)\.Bottom Right\)The series\-mode factorHTH^\{T\}, columns sorted by ground\-truth phase group, revealing a block structure consistent withT1T\_\{1\},T2T\_\{2\},T3T\_\{3\}\.To illustrate what is lost when the self\-similarity representation is abandoned in favor of scalable analytical methods, consider applying NMF directly to a raw data matrix of 60 synthetic time series, each a sinusoidal signal with additive noise and periodic sawtooth anomalies \(Figure[1\.1](https://arxiv.org/html/2608.04157#S1.F1)\)\. The series belong to three latent groups,T1T\_\{1\},T2T\_\{2\}, andT3T\_\{3\}, differing only in the phase offset at which anomalies begin: 0, 10, and 20 days respectively\. This is a controlled setting where the ground truth is known, making it possible to evaluate what a method recovers and what it does not\. NMF recovers the sensor\-mode grouping structure inHTH^\{T\}: when series are sorted by ground\-truth group membership, a clear block pattern emerges\. However,WW, operating in raw time space, does not surface when each group’s characteristic pattern occurs: a 10\-day phase difference spans only 240 samples out of 8,761, indistinguishable at figure scale without additional post\-hoc analysis\.

We now construct and decompose a tensorized self\-similarity matrix on the same dataset\. Figure[1\.2](https://arxiv.org/html/2608.04157#S1.F2)shows the resulting factor matricesAA,BB, andCC\. The structure is legible and internally consistent across all three modes\. InAA, the three phase groups emerge as distinct blocks when sorted by ground\-truth membership\. This sensor\-mode picture is confirmed by the temporal factors: inBBandCC, the bright bands align with the analytically computed phase\-offset positions forT1T\_\{1\},T2T\_\{2\}, andT3T\_\{3\}; these dashed lines are drawn after the fact from ground truth as a validation step, not as inputs to the decomposition\. All three matrices independently arrive at the same component\-to\-phase assignment, providing cross\-modal validation that the decomposition is recovering genuine periodic structure rather than an artifact of the representation\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/x2.png)Figure 1\.2:MINT applied to the same dataset\.Top\)The sensor\-mode factorAA, rows sorted by phase group, recovering the same three\-group block structure as NMF’sHH\.Bottom\)The temporal factorsBBandCC, operating in Mplot self\-similarity space\. Dashed lines mark analytically computed phase\-offset positions forT1T\_\{1\}\(red\),T2T\_\{2\}\(green\), andT3T\_\{3\}\(blue\)\. Bright bands inBBandCCalign with these positions, directly encoding when each group’s sawtooth pattern recurs\.The contribution of this work is not a new decomposition method; rather, it proposes a new object to decompose, the tensor of self\-similarity matrices\. We introduce tensorized self\-similarity matrices as a new time series data mining primitive, and propose MINT as a first method instantiating this primitive under Euclidean distance\. We show this representation concentrates cross\-series temporal structure in a form that standard tensor decomposition can extract, enabling analysis that was previously limited to what a practitioner could see by eye in individual Mplots\.

## 2Motivation and Related Work

##### Recurrence plots and Mplots\.

The essential insight of encoding a time series’ self\-similarity data in a matrix has a rich history of iterative development and application in many fields\. In dynamical systems analysis,recurrence plots\(RPs\) utilize pairwise state similarities to succinctly characterize a nonlinear time series’ dynamics\[[6](https://arxiv.org/html/2608.04157#bib.bib6)\]\. Since Eckmann et al\.\[[4](https://arxiv.org/html/2608.04157#bib.bib4)\]introduced the concept in 1987, extensions such as the cross recurrence plot \(CRP\)\[[16](https://arxiv.org/html/2608.04157#bib.bib16)\]and joint recurrence plot \(JRP\)\[[20](https://arxiv.org/html/2608.04157#bib.bib20)\]have generalized its scope to bivariate and multivariate time series, respectively\. RP\-based methods are tailored to analyzing and comparing recurrence behaviors among time series and have been used effectively in time series classification and anomaly detection\[[24](https://arxiv.org/html/2608.04157#bib.bib24),[3](https://arxiv.org/html/2608.04157#bib.bib3)\]\. Nevertheless, they are sensitive to thresholding and choice of embedding, and thus do not alone provide a single, continuous representation of both intra\-series and inter\-series self\-similarity\.

Recently, Kang et al\. proposed the RP\-GCN methodology\[[9](https://arxiv.org/html/2608.04157#bib.bib9)\], which exploits the topological recurrence properties of the RP for time series classification\. Given a multivariate time series, an RP is computed for each univariate time series\. The RPs are then stacked along a new axis to produce a33\-tensor, whose graph representation is finally processed by a graph convolutional network \(GCN\)\. While RP\-GCN captures the recurrence structure of a time series dataset, it only uses it as an intermediate representation to produce a global classification result, not as a mineable object for structure discovery\. It does not perform subsequence\-granularity co\-clustering, nor can its final result be easily mined to provide such information\.

Mplots\[[22](https://arxiv.org/html/2608.04157#bib.bib22)\]leverage the Matrix Profile as a recurrence plot, finding patterns in star light curves, animal sound waveforms, CCT telemetry in mice, and pulsus paradoxus \(an exaggerated drop in blood pressure\)\. A more formal definition is given in[section˜3](https://arxiv.org/html/2608.04157#S3), but an Mplot is roughly a matrix whose entries are the Euclidean distances between every pair of subsequences\. Intuitively, Mplots measure the self\-similarity of a time series in Euclidean Distance sliding window fashion; they cluster correlated subsequences in a single series, thus capturing key structural information\. Moreover, they are interpretable, enabling visual observation of time series trends, motifs, and discords\. For our objectives, the Mplot is thus the richest object in this line of work, but it has never been treated as a mineabledataset\-level\(rather thanseries\-level\) primitive, either for univariate or multivariate datasets\.

##### Tensor methods for time series\.

Tensor algebra algorithms historically developed in tandem with applications in psychometrics, chemometrics, and signal processing, but in the past two decades they have borne substantial fruit in machine learning and data mining\[[23](https://arxiv.org/html/2608.04157#bib.bib23)\]\. Prior work on tensor methods for time series largely treats tensors as theobserveddata structure \(e\.g\., annn\-tensor time series whose entries may be stacked to form a\(n\+1\)\(n\+1\)\-tensor\), modeling their dynamics for forecasting or dimensionality reduction through tensor decomposition \(e\.g\. CP or Tucker\) or tensor autoregressive models\[[1](https://arxiv.org/html/2608.04157#bib.bib1),[2](https://arxiv.org/html/2608.04157#bib.bib2)\]\. Cross\-series self\-similarity structure is absent from these formulations entirely, hindering their ability to capture subtle co\-clustering information at the subsequence level\.

##### Time series clustering\.

One of the fundamental tasks in time series analysis isclustering, partitioning data samples into groups in such a way that in\-group similarity is maximized and inter\-group similarity is minimized\. Individual clusters share domain\-relevant features and thus constitute potentialpatternsrevealing the latent structure of the time series\[[18](https://arxiv.org/html/2608.04157#bib.bib18)\]\.

Whole time\-series clustering \(WSC\) learns clusters of individual series within a collection of distinct time series\. Paparrizos et al\.\[[18](https://arxiv.org/html/2608.04157#bib.bib18)\]propose a comprehensive taxonomy of this space incorporating both classical and deep learning strategies, classifying clustering techniques into \(i\) distance\-based, \(ii\) distribution\-based, \(iii\) subsequence\-based, and \(iv\) representation\-learning\-based methods\. Subsequence\-based WSC methods in particular extract representative sequences to differentiate among time series clusters\. While they may leverage similarity data at the subsequence level—for example, the Matrix Profile in MPNCMI clustering\[[12](https://arxiv.org/html/2608.04157#bib.bib12)\]—they ultimately associate only a single label to each time series\. Moreover, they do not natively detect time\-aligned correlations between subsequences in different series\.

Subsequence\-based WSC is distinguished from subsequence clustering, which clusters time series subsequences obtained by sliding windows rather than the individual time series themselves\[[18](https://arxiv.org/html/2608.04157#bib.bib18)\]\. It is important to distinguish subsequence clustering from subsequenceco\-clustering\. Keogh and Lin\[[10](https://arxiv.org/html/2608.04157#bib.bib10)\]demonstrated that, for virtually all datasets, the results obtained by subsequence clustering are essentially random\. This “meaninglessness” is an inextricable result of the way subsequences are sampled by overlapping sliding windows, rather than features of the input data itself\. Subsequenceco\-clustering, on the other hand, clusters patterns along several modes \(in both time and sensor space, in our case\), rather than just in the time domain\. We leave formal verification for future work, but we propose that the tensorized self\-similarity matrix \(and the MINT pipeline\) constitutes an existence proof that this setting forces clustering methods to find genuine co\-occurrence\.

##### Conclusion\.

None of these works, individually or combined, produces an interpretable, mineable representation of cross\-series subsequence self\-similarity structure\.

## 3Definitions and Methodology

Atime seriesT=\{t1,t2,…,tn\}T=\\\{t\_\{1\},t\_\{2\},\.\.\.,t\_\{n\}\\\}is a sequence of real numbers, and asubsequenceT​\[i,j\]T\[i,j\]of lengthmmis a contiguous subset of values fromTTstarting at indexiiand ending atjj, wherej=i\+m−1j=i\+m\-1\. Adistance profileD​PA​Bj,mDP\_\{AB\}^\{j,m\}is the vector of distances between each subsequence in a reference time seriesTAT\_\{A\}and a query subsequenceTB​\[j,j\+m−1\]T\_\{B\}\[j,j\+m\-1\]\. We show examples of these primitives in Figure[3\.1](https://arxiv.org/html/2608.04157#S3.F1)\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/x3.png)Figure 3\.1:Top\)A time seriesTof Mass Rapid Transit usage\.Middle\)A subsequenceQfromTBottom\)Adistance profilecomputed using Mueen’s Algorithm for Similarity Search\[[17](https://arxiv.org/html/2608.04157#bib.bib17)\]ofQagainstT\.In this work, we consider Self\-Join Mplots, which leverage a special case of an AB\-Join Matrix Profile, whereTAT\_\{A\}andTBT\_\{B\}are the same time series\[[22](https://arxiv.org/html/2608.04157#bib.bib22)\]\.

###### Definition 3\.1\.

AnMplotis the visualized matrix of distance profiles\. Each row of this matrix isD​PA​Bj,mDP^\{j,m\}\_\{AB\}\.

There are resources available on how to read dot/recurrence plots\[[15](https://arxiv.org/html/2608.04157#bib.bib15),[19](https://arxiv.org/html/2608.04157#bib.bib19),[14](https://arxiv.org/html/2608.04157#bib.bib14),[22](https://arxiv.org/html/2608.04157#bib.bib22)\], but to keep this work self\-contained, we reiterate a few key points, discussing “light values” as having low similarity and “dark values” having high similarity:

- •The Mplot is a subsequence self\-similarity matrix; consequently, a visible central diagonal has a distance of zero\.
- •The appearance of evenly spaced diagonals parallel to the central diagonal indicates periodicity or cycles in the time series\.
- •Mplots can result in a type of pattern where the slope changes, either “up” or “down”, indicating patterns are either slowing down or speeding up\.
- •A streak that curves in both directions suggests a pair of subsequences that match after one locally “warps” in order to match the other\. This is a useful characteristic that allows us to spot patterns in a warping\-invariant manner\.

The key novelty in this work is to treat Mplots as decomposable matrices and to extend this perspective to decomposing constructed stacked Mplots using classical tensor decomposition methods\.

###### Definition 3\.2\.

Ana×ba\\times bmatrix is a function of two integer indices; that is, a functionM:\[1\.\.a\]×\[1\.\.b\]→ℝM:\[1\.\.a\]\\times\[1\.\.b\]\\to\\mathbb\{R\}\. We view a matrix as a “2\-dimensional array” whereM​\(i,j\)M\(i,j\)denotes the entry at theiith row andjjth column\.

Given ana×ba\\times bmatrixMM,

1. 1\.Itsentrywiseℓ1\\ell\_\{1\}normis the sum of the absolute values of its entries: ‖M‖1=∑i=1a∑j=1b\|M​\(i,j\)\|\.\|\|M\|\|\_\{1\}=\\sum\_\{i=1\}^\{a\}\\sum\_\{j=1\}^\{b\}\|M\(i,j\)\|\.
2. 2\.Itsnuclear normis the sum of its singular values: ‖M‖∗=∑i=1min⁡\(a,b\)σi​\(M\)\.\|\|M\|\|\_\{\*\}=\\sum\_\{i=1\}^\{\\min\(a,b\)\}\\sigma\_\{i\}\(M\)\.

###### Definition 3\.3\.

Atensorof ordernnis a function ofnninteger indices; that is, a functionT:\[1\.\.k1\]×…×\[1\.\.kn\]→ℝT:\[1\.\.k\_\{1\}\]\\times\.\.\.\\times\[1\.\.k\_\{n\}\]\\to\\mathbb\{R\}\. The positive integerklk\_\{l\}is the size of the tensor along thellth mode or axis\.

Informally, a tensor of ordernnis an “nn\-dimensional matrix\.” We say a tensorTTisrank\-11if it is the outer product ofnnvectors\.

T=𝐚1⊗…⊗𝐚n\.T=\\mathbf\{a\}\_\{1\}\\otimes\.\.\.\\otimes\\mathbf\{a\}\_\{n\}\.
On a similar line, we say thatTTisrank\-RRif it can be expressed using a minimum ofRRouter products ofnnvectors\[[23](https://arxiv.org/html/2608.04157#bib.bib23)\]\.

T=∑i=1R𝐚1,i⊗…⊗𝐚n,i\.T=\\sum\_\{i=1\}^\{R\}\\mathbf\{a\}\_\{1,i\}\\otimes\.\.\.\\otimes\\mathbf\{a\}\_\{n,i\}\.
Our analysis primarily concerns tensors constructed by stacking Mplots along one axis\. Our use of the termtensorrefers to a third\-order tensor \(3\-tensor\) unless otherwise specified\. For a given time series dataset𝒟=\{T1,T2,…,TN\}\\mathcal\{D\}=\\\{T\_\{1\},T\_\{2\},\.\.\.,T\_\{N\}\\\}, the resulting tensors constructed from stacked Mplots are of sizeN×\(n−m\+1\)×\(n−m\+1\)N\\times\(n\-m\+1\)\\times\(n\-m\+1\), wherennis the length of the time series in the dataset,mmthe subsequence length parameterizing the Mplots, andNNthe number of time series in𝒟\\mathcal\{D\}\. The methods used to analyze the objects defined above are:

Nonnegative Robust PCA\. We adopt a Robust PCA formulation derived from the Stable PCP of Zhou et al\.\[[29](https://arxiv.org/html/2608.04157#bib.bib29)\]\. Given ann×nn\\times nmatrixMM, we find a nonnegative low rank matrixLLand a sparse matrixSSsuch that

whereZZis a bounded entry\-wise noise term\. Ifλ=1n\\lambda=\\frac\{1\}\{\\sqrt\{n\}\}is the standard RPCA sparsity parameter, then we consider the problem

minL≥0,S​‖L‖∗\+λ​‖S‖1\+μ2​‖M−L−S‖F2\.\\min\_\{L\\geq 0,S\}\|\|L\|\|\_\{\*\}\+\\lambda\|\|S\|\|\_\{1\}\+\\frac\{\\mu\}\{2\}\|\|M\-L\-S\|\|^\{2\}\_\{F\}\.We obtain a heuristic estimate\(L^,S^\)\(\\hat\{L\},\\hat\{S\}\)by alternating minimization\. Initially,L^=S^=0\\hat\{L\}=\\hat\{S\}=0\. In each iteration,L^\\hat\{L\}is computed by singular\-value thresholding ofM−S^M\-\\hat\{S\}at level1μ\\frac\{1\}\{\\mu\}followed by a nonnegativity clamp, whileS^\\hat\{S\}is computed by soft\-thresholding ofM−L^M\-\\hat\{L\}at levelλμ\\frac\{\\lambda\}\{\\mu\}\. We iterate until the relative iterate change inL^\\hat\{L\}falls belowϵ=10−5\\epsilon=10^\{\-5\}or after a maximum of50005000iterations\.

AssumeZZhas i\.i\.d\.N​\(0,σ2\)N\(0,\\sigma^\{2\}\)entries\. With penalty12​μ​‖M−L−S‖F2\\frac\{1\}\{2\\mu\}\|\|M\-L\-S\|\|\_\{F\}^\{2\}, Zhou et al\. prescribeμ=2​n​σ\\mu=\\sqrt\{2n\}\\sigma\. Our penalty isμ2​‖M−L−S‖F2\\frac\{\\mu\}\{2\}\|\|M\-L\-S\|\|^\{2\}\_\{F\}, so the equivalent choice isμ=1/\(2​n​σ\)\\mu=1/\(\\sqrt\{2n\}\\sigma\)\. We inherit the noise\-calibrated choice ofμ\\mufrom the original formulation\[[29](https://arxiv.org/html/2608.04157#bib.bib29)\]\. MINT adds a nonnegativity constraint onLL, and we adoptμ\\mu\(which tunes the residual term to the noise levelσ\\sigma\) as a reasonable heuristic for our form of the objective\.111Zhou et al\. establish stable recovery guarantees for Stable PCP under the same conditions as PCP\. We leave analogous guarantees for the nonnegative\-constrained variant and our projected iteration to future work\.

It remains to estimateσ\\sigmagivenMM\. We accomplish this through a two\-step algorithm\.

1. 1\.We recover a coarse low\-rank approximationMτM\_\{\\tau\}by hard\-thresholding singular values at the optimal thresholdτ=2\.858⋅ym​e​d\\tau=2\.858\\cdot y\_\{med\}, whereym​e​dy\_\{med\}is the median empirical singular value\[[5](https://arxiv.org/html/2608.04157#bib.bib5)\]\. We then estimate the entry\-wise noise by the residualZτ=M−MτZ\_\{\\tau\}=M\-M\_\{\\tau\}\.
2. 2\.BecauseZτZ\_\{\\tau\}is effectively symmetric, each off\-diagonal residual appears twice; we compute the median absolute deviation \(MAD\) over its strict upper triangular values to count each value once\. Then we estimate σ^=MAD0\.6745\.\\hat\{\\sigma\}=\\frac\{\\text\{MAD\}\}\{0\.6745\}\.This parallels the formulation in Sato and Ono, Section IV\-I\[[21](https://arxiv.org/html/2608.04157#bib.bib21)\]\.

Canonical polyadic decomposition \(CPD\)\. Given an order\-nntensorTTand a rankRR, CPD returnsnnmatrices

A1=\[𝐚1,1,…,𝐚1,R\],A\_\{1\}=\[\\mathbf\{a\}\_\{1,1\},\.\.\.,\\mathbf\{a\}\_\{1,R\}\],An=\[𝐚n,1,…,𝐚n,R\],A\_\{n\}=\[\\mathbf\{a\}\_\{n,1\},\.\.\.,\\mathbf\{a\}\_\{n,R\}\],such that

T=∑i=1R𝐚1,i⊗…⊗𝐚n,i\.T=\\sum\_\{i=1\}^\{R\}\\mathbf\{a\}\_\{1,i\}\\otimes\.\.\.\\otimes\\mathbf\{a\}\_\{n,i\}\.In other words, thekkth matrix gives thekkth terms of each outer product in the outer product decomposition\[[23](https://arxiv.org/html/2608.04157#bib.bib23)\]\. We interpret each of theRRcolumns in the matrices as representing a uniquecomponentorpattern\. If we stack Mplots along the first axis of a33\-tensor, then we interpretAk​\(i,j\)A\_\{k\}\(i,j\)as the expression of componentjjin Mplotii\(ifkkequals11\) or in the distance profile of subsequenceii\(ifkkequals22or33\)\. We introduce the tensorized self\-similarity matrix as a new time series data mining primitive, instantiated by MINT \(MplotsINtoTensor\) as an initial method of computation\. MINT is described by the pipeline shown in Figure[3\.2](https://arxiv.org/html/2608.04157#S3.F2):

1. 1\.Data loading and preprocessing\.Given a dataset𝒟\\mathcal\{D\}, we extract time seriesT1,T2,…,TN∈𝒟T\_\{1\},T\_\{2\},\.\.\.,T\_\{N\}\\in\\mathcal\{D\}, aligned on a common axis with nearly uniform intervals\. While implementations of the Matrix Profile can handle missing values in the input data, MINT requires NaNs to be preprocessed\.
2. 2\.Mplot Computation\.We use the SCAMP library\[[30](https://arxiv.org/html/2608.04157#bib.bib30)\]to compute a set of Mplots for allTi∈𝒟T\_\{i\}\\in\\mathcal\{D\}\. It appears at first glance that the time complexity of Mplot construction isO​\(m​n2\)O\(mn^\{2\}\), wheremmis the subsequence length parameterizing the Mplot computation\. However, advances in GPU implementations\[[30](https://arxiv.org/html/2608.04157#bib.bib30),[22](https://arxiv.org/html/2608.04157#bib.bib22)\]render this process more reflective ofO​\(n2\)O\(n^\{2\}\)in practice\. We limit Mplots to ak×kk\\times kwindow, wherek<nk<nis a predefinedpooling parameter; in the SCAMP library, pooling is accomplished usingmax\(\)for a similarity matrix andmin\(\)for a distance matrix\. Prior work in hyperspectral imaging\[[8](https://arxiv.org/html/2608.04157#bib.bib8)\]shows that lossy compression of sensor data can retain classification\-relevant structure; we observe the analogous property for min\-pooled distance matrices\.
3. 3\.Tensor Construction\.For each MplotMiM\_\{i\}, we use Nonnegative Robust PCA to decompose it into a low\-rank matrixLiL\_\{i\}and a sparse matrixSiS\_\{i\}\. Using exact SVD, this can be accomplished in approximatelyO​\(niter​k3\)O\(n\_\{\\text\{iter\}\}k^\{3\}\)time, whereniter≤5000n\_\{\\text\{iter\}\}\\leq 5000is the number of iterations\. In practice, Robust PCA is the primary bottleneck of the proposed pipeline, particularly with large Mplots\. The low\-rank matrices are mean\-centered, then stacked to form a 3\-tensorTTof dimensionsN×k×kN\\times k\\times k, which is then magnitude\-normalized\.
4. 4\.CPD Decomposition\.We use TensorLy’s\[[11](https://arxiv.org/html/2608.04157#bib.bib11)\]CANDECOMP/PARAFAC method to perform CPD onTTinO​\(niter​k2​N\)O\(n\_\{\\text\{iter\}\}k^\{2\}N\)time \(wherenitern\_\{\\text\{iter\}\}is again the number of iterations\), thereby producing three matrices A=\[𝐚1,…,𝐚R\],B=\[𝐛1,…,𝐛R\],C=\[𝐜1,…,𝐜R\],A=\[\\mathbf\{a\}\_\{1\},\.\.\.,\\mathbf\{a\}\_\{R\}\],B=\[\\mathbf\{b\}\_\{1\},\.\.\.,\\mathbf\{b\}\_\{R\}\],C=\[\\mathbf\{c\}\_\{1\},\.\.\.,\\mathbf\{c\}\_\{R\}\],of dimensionsN×RN\\times R,k×Rk\\times R, andk×Rk\\times Rrespectively\. The rankR∈\[2\.\.8\]R\\in\[2\.\.8\]is estimated by locating the elbow point of the reconstruction error curve using the Kneedle algorithm\. If no elbow is found, a flag is raised and the procedure defaults to22components\. The process returns threefactor matrices, which are then further analyzed\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/figures/mint_image.png)Figure 3\.2:Overview of the MINT pipeline\. Analysis of peaks in individual loadings enables easy visual identification of cross\-sensor phenomena, such as Lunar New Year festivities\.This setting is distinct from multivariate time series analysis, where thedddimensions are semantically coupled by construction\. In our setting, univariate series are independent sensors whose co\-clustering structure islatent, recoverable without labels or assumptions, rather than given\. The reader will also appreciate that this primitive is grounded in tensor mining, and naturally extends to multivariate time series datasets\.

It remains to justify mean\-centering low\-rank matrices in step \(3\)\. Each series’ Mplot \(hence low\-rank matrix\) has a varying amount of self\-similarity, from highly regular or periodic to noisy\. A series’ baseline self\-similarity level, considered as a mean\-determined offset of the similarity matrix, might describe it in isolation, but it will not contribute information about when patterns recur or which series also share them\. If not subtracted, the offset carries a large Frobenius norm \(“energy”\) contribution, and CPD will spend a component ranking sensors by overall self\-similarity\. With a limited rank range, mean\-centering of low\-rank matrices enables us to surface genuine co\-clustering information rather than per\-slice global similarity backgrounds\.

## 4Datasets

We present our results through a set of five case studies on diverse datasets from various domains with different sampling intervals\. TheMass Rapid Transit \(MRT\)dataset\[[27](https://arxiv.org/html/2608.04157#bib.bib27)\]contains entry and exit ridership data from 108 Taipei stations sampled hourly during operating hours over a total length of 517 days between November 1, 2015 and March 31, 2017\. Sample subsequences from this dataset are shown in Figure[3\.2](https://arxiv.org/html/2608.04157#S3.F2)and Figure[3\.1](https://arxiv.org/html/2608.04157#S3.F1)\. Sample time series from theSpatio\-Temporal Traffic \(LargeST\)dataset\[[13](https://arxiv.org/html/2608.04157#bib.bib13)\]are shown in Figure[5\.1](https://arxiv.org/html/2608.04157#S5.F1)\. This dataset consists of traffic data from a total of 8,600 sensors in California, USA over five years, at a base sampling rate of 5 minutes\. Prior to analysis, we extracted data from the year 2019 and aggregated them into 15\-minute means\. TheElectricity Load Diagrams dataset\[[25](https://arxiv.org/html/2608.04157#bib.bib25)\]contains the quarter\-hourly electricity consumption \(in kW\) of 370 clients from 2011 to 2014\. We used only 2013, a year in which most clients had complete records, to minimize long initial zero sequences, and converted the data to units of kWh before analysis\.

Wind Farm A of theCARE to Comparedataset\[[7](https://arxiv.org/html/2608.04157#bib.bib7)\]is based on data from the EDP open platform222https://www\.edp\.com/en/innovation/open\-data/dataand consists of five wind turbines of an onshore wind farm in Portugal\. The wind farm data contains2222datasets with a total of8686features\. For a given turbine, various sensor data traces are available, including but not limited to the ambient temperature, wind direction, temperature recordings from various turbine components \(nacelle controller, gearbox, hub controller, generator, etc\.\), voltage, and current\. 54 unique sensor traces are represented by derived statistical features such as the1010\-minute average, max, min, and standard deviation\.

Testing was performed on Dataset 0, comprising one year of training data and a prediction window containing a labeled anomaly window preceding a generator bearing failure, surrounded by padding\. Before processing, only the1010\-minute average of each sensor trace is taken, which results in5454considered sensors\. Sample subsequences from Dataset 0 are shown in Figure[6\.1](https://arxiv.org/html/2608.04157#S6.F1)\.

TheElectrical Power Demand \(OPSD\)dataset333https://open\-power\-system\-data\.org/is the recorded electrical demand of 17 European countries at hourly sampling intervals\. Each country has at least one time series, all of which span dates from December 31, 2014 to September 30, 2020, or 50,401 time steps\.

### 4\.1Case Study: OPSD

We provide an existence proof for MINT’s ability to surface real structure which no single\-series method finds\. We run the pipeline on the OPSD dataset\. Due to a small number of sensors \(18\), we impute NaN\-filled sequences with random noise\. From the method, we obtain three factor matricesAA,BB, andCC, which are then visually analyzed for patterns\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/figures/opsd_figure.jpg)Figure 4\.1:\(1\) Transformed factor matrixCCfor the OPSD dataset\. \(2\) Analysis of Components 2 and 3 indicates that they correspond to distinctive Christmas and seasonal electricity demand patterns, respectively\.Figure[4\.1](https://arxiv.org/html/2608.04157#S4.F1)depicts this analysis applied on factor matrixCC\(in particular, its entrywise absolute value\|C\|\|C\|\)\. The vertical axis corresponds to time, with each pixel signifying a week\-length interval\. For each component, we identify highly expressed positions \(indicesiifor which\|C​\(i,j\)\|\|C\(i,j\)\|is a high value\) and map them back to dates in the original data\.

From this process, we conclude that the MINT components compartmentalize distinct cross\-sensor temporal patterns\. For example, two out of the five highest expressed positions in Component 2 correspond to Christmas 2015 and 2019, while two more are in the 2014 and 2016 holiday seasons\. Component 3, on the other hand, corresponds to seasonal electricity patterns\. Its signal follows an annual schedule with high peaks in winter, a slightly smaller peak in summer, and troughs in the intervening fall and spring\. In particular, even though Component 3 troughs during the holiday season, it still attains a background “winter season” level that exceeds fall and spring\. We expect short\-scale holiday phenomena and long\-scale seasonal periodicities to correspond to distinct electricity demand signatures, and this is reflected by MINT’s component segregation\.

The tensorial perspective is crucial: while a single\-series method \(such as the Matrix Profile\) applied on each sensor’s time series would detect distinct electricity demand patterns each winter, it would be difficult to distinguish intensities and isolate contributing factors \(e\.g\. cold temperatures versus holiday festivities\) without additional contextual information\. In contrast, the highlighted components in the MINT tensor generalize information from diverse time series\. Thus, MINT not only identified but alsoseparatedChristmas events and seasonal patterns acrossalltime series in the dataset, regardless of other confounding phenomena that might occur in each time series\.

## 5Experimental Setup

This work’s primary novelty lies in identifying co\-clustering phenomena in time series, capturing both seasonality and shapelets simultaneously inO​\(N⋅n2\+N⋅k3\)O\(N\\cdot n^\{2\}\+N\\cdot k^\{3\}\)time, wherennis the time series length,kkis the side length of the pooled Mplots, andNNis the number of series transformed into Mplots to construct the tensor\. Deep clustering approaches for time series typically rely on learned representations and iterative optimization\[[18](https://arxiv.org/html/2608.04157#bib.bib18)\], resulting in unpredictable runtime and resource requirements that stand in contrast to the analyzable computational complexity of MINT\. We test our claims to co\-clustering through robustness testing by Noise Injection\.

### 5\.1Single\-Trial Procedure

Let us describe a univariate time seriesTTwith a set ofNNsensors, where the notationTiT\_\{i\}indicates theiith sensor\. We chooseb<N/2b<N/2and a set𝒲\\mathcal\{W\}of disjoint windows of lengthmmfromTT\. We construct𝒲\\mathcal\{W\}by dividingTiT\_\{i\}into windows of lengthmmand randomly selecting\|𝒲\|\|\\mathcal\{W\}\|of these windows\. For each traceTi∈TT\_\{i\}\\in T, we randomly assign one of the following transformationsθx\\theta\_\{x\}to produceTi^\\hat\{T\_\{i\}\}:

- •bbtime seriesTi∈TT\_\{i\}\\in TareCompletely Random\(θ1​\(Ti\)\\theta\_\{1\}\(T\_\{i\}\)\):Ti^\\hat\{T\_\{i\}\}is reassigned to be random noise\. We measure the trace mean and standard deviation to construct the signal of random normal noise\.
- •bbtime series areMostly Random\(θ2​\(Ti\)\\theta\_\{2\}\(T\_\{i\}\)\):Ti^\\hat\{T\_\{i\}\}is transformed to be identical to its normal time series only in the\|𝒲\|\|\\mathcal\{W\}\|disjoint windows, and is otherwise random noise, as described above\.
- •N−2​bN\-2btime series areMostly Normal\(θ3​\(Ti\)\\theta\_\{3\}\(T\_\{i\}\)\):Ti^\\hat\{T\_\{i\}\}is identical to its normal time series everywhere except the\|𝒲\|\|\\mathcal\{W\}\|disjoint windows, which are replaced with random noise\.

In Figure[5\.1](https://arxiv.org/html/2608.04157#S5.F1), we show the original time seriesTiT\_\{i\}and its three transformed series:θ1​\(Ti\)\\theta\_\{1\}\(T\_\{i\}\),θ2​\(Ti\)\\theta\_\{2\}\(T\_\{i\}\),θ3​\(Ti\)\\theta\_\{3\}\(T\_\{i\}\)\. After computing the Mplots for allT^=\{T1^,…,TN^\}\\hat\{T\}=\\\{\\hat\{T\_\{1\}\},\.\.\.,\\hat\{T\_\{N\}\}\\\}, we decompose the tensor using CPD\. For each resulting componentϕ\\phiin the decomposition, we obtain a set of top\-sssensorsS​\(ϕ\)S\(\\phi\)from a distinctiveness pipeline\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/x4.png)Figure 5\.1:Four subsequences from four time series from the LargeST dataset, surrounding one disjoint windowtt\. The Completely Randomθ1​\(Ti\)\\theta\_\{1\}\(T\_\{i\}\), Mostly Randomθ2​\(Ti\)\\theta\_\{2\}\(T\_\{i\}\), and Mostly Normalθ3​\(Ti\)\\theta\_\{3\}\(T\_\{i\}\)transforms on the originalTiT\_\{i\}are shown\. In Mostly Random subsequences, we seettcontains the original data, and the remainder of the series is random noise\. In the Mostly Normal series,ttis replaced with random noise, and the remainder of the series is the original data\.#### 5\.1\.1Distinctiveness Pipeline

Given theiith componentϕ\\phi, we expect the sensors \(or multidimensional traces\)TjT\_\{j\}for which\|A​\(j,i\)\|\|A\(j,i\)\|is highest would exhibit the highest “expression” ofϕ\\phi\. Similarly, among subsequences within the sensors, those indexed byjjwith the largest\|B​\(j,i\)\|\|B\(j,i\)\|and\|C​\(j,i\)\|\|C\(j,i\)\|values are expected to most strongly expressϕ\\phi\. Combining these observations, we conclude that the top\-sslengthmmsubsequences of the top\-sssensors represent “archetypal patterns” that best capture the structural features encoded byϕ\\phi\. We identify these sequences for allϕ\\phithrough thisDistinctiveness Pipeline:

1. 1\.Use CPD to find the factor matricesAA,BB, andCC\.
2. 2\.In componentϕ\\phiofAA, find the top\-ssindices and find the corresponding “sensor” names and store in a setS​\(ϕ\)S\(\\phi\)\.
3. 3\.In componentϕ\\phiofBB\(orCC\), find the top\-ssindices and map them to their corresponding time series intervals\. To account for pooling effects, we pad these intervals minimally: if an interval\[a,b\]\[a,b\]is mapped and each pixel on the Mplot corresponds tockc\_\{k\}positions, we extend the interval byckc\_\{k\}on both sides, yielding\[a−ck,b\+ck\]\[a\-c\_\{k\},b\+c\_\{k\}\]\. Store these adjusted intervals in a setI​\(ϕ\)I\(\\phi\)\.
4. 4\.For each interval\[i,j\]∈I​\(ϕ\)\[i,j\]\\in I\(\\phi\)and each sensorς∈S​\(ϕ\)\\varsigma\\in S\(\\phi\), consider the subsequenceTς​\[i,j\]T\_\{\\varsigma\}\[i,j\]\.

The choice ofBBorCCis arbitrary because \(up to numerical error\) they are effectively equivalent; for our experiments, we fix the factor matrix used asCC\.

### 5\.2Validation Criteria

UsingS​\(ϕ\)S\(\\phi\)andI​\(ϕ\)I\(\\phi\), we consider four properties to demonstrate the utility of MINT to find co\-clustering information\.

1. 1\.𝐏𝟏\\mathbf\{P\_\{1\}\}: There exists aϕ\\phisuch that for allI∈I​\(ϕ\)I\\in I\(\\phi\)there existsI′∈𝒲I^\{\\prime\}\\in\\mathcal\{W\}for whichI′∩I≠∅I^\{\\prime\}\\cap I\\neq\\emptyset\.
2. 2\.𝐏𝟐\\mathbf\{P\_\{2\}\}: Each component strongly reflects sensors transformed by exactly one of the transformationsθ1\\theta\_\{1\},θ2\\theta\_\{2\},θ3\\theta\_\{3\}\. That is, for eachϕ\\phi, there existsx∈\{1,2,3\}x\\in\\\{1,2,3\\\}such thatθx\\theta\_\{x\}correspondstoϕ\\phi, i\.e\. more than half of the elements ofS​\(ϕ\)S\(\\phi\)are transformed by the sameθx\\theta\_\{x\}\.
3. 3\.𝐏𝟑\\mathbf\{P\_\{3\}\}: There exists at most one componentϕ\\phis\.t\.S​\(ϕ\)S\(\\phi\)containsθ1\\theta\_\{1\}\(Completely Random\) sensors\.
4. 4\.𝐏𝟒\\mathbf\{P\_\{4\}\}:θ3\\theta\_\{3\}\(Mostly Normal\)correspondsto asignificant majorityof the components\. We fix the threshold for this criterion at⌊0\.7​R⌋\\lfloor 0\.7R\\rfloor, whereRRis the number of components\. This choice permits small decompositions \(e\.g\.,R=2R=2\) to remain valid even if only one component is classified asMostly Normal\.

Each of the properties reflects a desirable aspect of the cross\-mode co\-clustering behavior that our pipeline is designed to exhibit\.𝐏1\\mathbf\{P\}\_\{1\}measures the ability of MINT to account for strong regularities in similar positions across series\.𝐏2\\mathbf\{P\}\_\{2\}ensures similarly transformed sensors are grouped in similar categories, and𝐏3\\mathbf\{P\}\_\{3\}and𝐏4\\mathbf\{P\}\_\{4\}formalize the assumption that the dominant sources of variation in the transformed dataset reflect the original data \(as opposed to noise imputation\)\.

### 5\.3Hypothesis Testing

In the random\-selection model, the top\-sssensors and the left endpoints of top\-ssintervals are selected uniformly and without replacement from\[1\.\.N\]\[1\.\.N\]and\[1\.\.n\]\[1\.\.n\]respectively, independently for each of theRRcomponents\. Each property is a deterministic function of the resulting category count, so its chance pass rate can be bounded above by a hypergeometric or binomial tail probability\. With the experimental constants defined in Section[5\.4](https://arxiv.org/html/2608.04157#S5.SS4), chance pass rates are bounded in Appendix[B](https://arxiv.org/html/2608.04157#S2a)and characterized below\.

p1chance\\displaystyle p\_\{1\}^\{\\text\{chance\}\}≤0\.00584,\\displaystyle\\leq 0\.00584,p2chance\\displaystyle p\_\{2\}^\{\\text\{chance\}\}≤0\.176,\\displaystyle\\leq 0\.176,p3chance\\displaystyle p\_\{3\}^\{\\text\{chance\}\}≤0\.325,\\displaystyle\\leq 0\.325,p4chance\\displaystyle p\_\{4\}^\{\\text\{chance\}\}≤0\.570\.\\displaystyle\\leq 0\.570\.Therefore, to differentiate genuine co\-clustering robustness from that expected under random selection, we must examine theproportionof incidence for each property\. To this end, the procedure defined in Section[5\.1](https://arxiv.org/html/2608.04157#S5.SS1)is run on the same time seriesTTfor 50 trials\. Across all trials, we compute the proportion of trials for which each of𝐏1,𝐏2,𝐏3,and​𝐏4\\mathbf\{P\}\_\{1\},\\mathbf\{P\}\_\{2\},\\mathbf\{P\}\_\{3\},\\text\{and \}\\mathbf\{P\}\_\{4\}is satisfied \(thus obtainingp^1,p^2,p^3,and​p^4\\hat\{p\}\_\{1\},\\hat\{p\}\_\{2\},\\hat\{p\}\_\{3\},\\text\{ and \}\\hat\{p\}\_\{4\}respectively\)\. To ensure that a chance process lies inside the null region, we define the hypotheses

\(H0,i\):\\displaystyle\(H\_\{0,i\}\):pi≤0\.7,\\displaystyle~p\_\{i\}\\leq 0\.7,\(H1,i\):\\displaystyle\(H\_\{1,i\}\):pi\>0\.7,\\displaystyle~p\_\{i\}\>0\.7,fori=1,2,3,4\.\\displaystyle\\text\{for $i=1,2,3,4$\}\.We reject or fail to reject each null hypothesis using a one\-tailed z\-test for proportions with a Bonferroni\-corrected significance ofα=0\.0125\\alpha=0\.0125\.

### 5\.4Experiment Details

We performed the experiment pipeline on four datasets: \(1\) Taipei MRT entrance data, \(2\) 2019 LargeST Traffic data, \(3\) 2013 Electricity Load Diagrams data, and \(4\) CARE to Compare wind turbine data\. In the four experiments, all sensor columns were used with the exception of those resulting in Mplots containing NaNs\. This resulted in yields of 108/108 sensors for Taipei MRT, 8303/8600 sensors for LargeST, 332/370 sensors for Electricity Load Diagrams, and 43/54 sensors for CARE to Compare\. In experiments \(1\) and \(4\), each trial was tested on all NaN\-free sensors, while in experiments \(2\) and \(3\), each trial was tested on a sample of100100randomly selected NaN\-free sensors\.

Experiments were conducted on a Windows 10 virtual machine with1616vCPUs and6464GiB of RAM\. The trials were conducted in parallel and independently on separate deterministic seeds\. Trials that yielded NaNs in Mplots, when identified, were discarded altogether and re\-tried on a new seed\. In total,5050valid trials were obtained for each dataset\.

For the four datasets, we set the subsequence lengthsm1=146m\_\{1\}=146,m2=m3=95m\_\{2\}=m\_\{3\}=95, andm4=140m\_\{4\}=140, respectively\. These selections correspond to practical use cases in identifying week\-length, day\-length, and day\-length motifs respectively\. Following\[[22](https://arxiv.org/html/2608.04157#bib.bib22)\], we select subsequence lengthsmmto be slightly less than the length of the motifs of interest to help ensure these motifs are visible in the resulting Mplot\. For most datasets and practical applications, Mplots are generally forgiving with respect to the exact value ofmm\. We set the pooling parameterk=⌈3​nm⌉k=\\lceil\\frac\{3n\}\{m\}\\rceil, so that each Mplot pixel corresponds tock≈m3c\_\{k\}\\approx\\frac\{m\}\{3\}positions\. We define the co\-clustering window set𝒲\\mathcal\{W\}such that it covers 10% of each time series\. We chooses=6s=6as the minimum number of top\-sssensors for which our four properties remain nontrivial\. Finally, we setb=⌈0\.25​N⌉b=\\lceil 0\.25N\\rceilas the number of sensors transformed by each ofθ1\\theta\_\{1\}andθ2\\theta\_\{2\}\.

## 6Results and Discussion

Table[6\.1](https://arxiv.org/html/2608.04157#S6.T1)lists the proportionsp^i\\hat\{p\}\_\{i\}fori=1,2,3,4i=1,2,3,4obtained after conducting5050trials on the four datasets tested\.

Taipei MRTLargeSTElectricityLoad DiagramsCARE to Compare\(Wind turbines\)p^1\\hat\{p\}\_\{1\}1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)0\.92\(<0\.001<0\.001\)p^2\\hat\{p\}\_\{2\}1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)p^3\\hat\{p\}\_\{3\}1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)p^4\\hat\{p\}\_\{4\}1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)1\.00\(<0\.001<0\.001\)

Table 6\.1:Estimated proportionsp^1\\hat\{p\}\_\{1\}–p^4\\hat\{p\}\_\{4\}with p\-values \(in parentheses\) for hypothesis testsH0,1H\_\{0,1\}–H0,4H\_\{0,4\}\.For all four experiments, the null hypothesesH0,iH\_\{0,i\}\(fori=1,2,3,4i=1,2,3,4\) were rejected at the0\.01250\.0125significance level\. We thus conclude that for the tested datasets, each of𝐏1\\mathbf\{P\}\_\{1\},𝐏2\\mathbf\{P\}\_\{2\},𝐏3\\mathbf\{P\}\_\{3\}, and𝐏4\\mathbf\{P\}\_\{4\}occurs inmost\(that is, more than 70 percent\) of the tensor decompositions produced by the MINT pipeline\. Since the random\-selection model for each property falls inside the null region, this suggests that MINT meaningfully co\-clusters time\-aligned cross\-series phenomena\.

![Refer to caption](https://arxiv.org/html/2608.04157v1/x5.png)Figure 6\.1:Last week of data for four time series from the CARE to Compare dataset\. Observe the time\-aligned degenerate flatlines\.Observe that properties𝐏2,𝐏3,and​𝐏4\\mathbf\{P\}\_\{2\},\\mathbf\{P\}\_\{3\},\\text\{and \}\\mathbf\{P\}\_\{4\}were satisfied on all5050trials for the four datasets\. However, only92%92\\%of trials satisfied𝐏1\\mathbf\{P\}\_\{1\}in the CARE to Compare wind turbine data, as compared to100%100\\%satisfaction in the other three datasets\.

We speculate that there may be two underlying reasons for the CARE dataset’s unique behavior\.

1. 1\.Unlike the other three datasets, CARE to Compare hasheterogeneoussensors: distinct sensors \(e\.g\. ambient temperature versus wind direction\) measure different quantities\. Therefore, aligned windows inθ3\\theta\_\{3\}\-transformed sensor traces are less globally similar and thus less likely to intersect the “archetypal patterns” \(i\.e\. top\-ssintervals\) of a given component\.
2. 2\.Degenerate flatlines predominate within the post\-anomaly tails of specific sensors, as illustrated in Figure[6\.1](https://arxiv.org/html/2608.04157#S6.F1)\. These highly time\-aligned cross\-series patterns, if not disturbed by noise, may become top\-ssintervals or even dominate a high\-energy CPD component, decreasing the incidence of top\-ssintervals in the firstRRcomponents intersecting𝒲\\mathcal\{W\}\-windows\.

While we leave formal verification of \(1\) and \(2\) and their precise linkage to incidence rates of𝐏1\\mathbf\{P\}\_\{1\}to future work, the results of all four hypotheses retaining their significance lend nontrivial support to the tensorized self\-similarity matrix primitive as instantiated in the MINT pipeline, even in heterogeneous datasets\.

## 7Conclusions and Future Work

We introduce the tensorized self\-similarity matrix as a novel time series data mining primitive, instantiated through MINT as an initial method of computation\. Where neural approaches implicitly aggregate cross\-series information through learned representations, this primitive makes cross\-series structure explicitly recoverable and interpretable\. Each component of the decomposition directly indexes the sensors, subsequences, and time intervals that participate in a shared pattern\. The low\-rank structure recovered by tensor decomposition corresponds to interpretable cross\-series phenomena, as demonstrated across MRT, wind turbine, traffic, and electricity datasets spanning multiple domains and sampling regimes\. Future work includes extension to alternative self\-similarity measures \(e\.g\. dynamic time warping, phase\-space embeddings\), and additional application to multivariate time series datasets to further characterize the primitive’s behavior across diverse analytical settings\.

## Acknowledgments\.

The work at UCR was partially supported by the National Science Foundation under CAREER grant IIS \#2046086\. The views and conclusions in this paper are those of the authors and should not be interpreted as representing any funding agencies\.

In the final two weeks leading to submission, GPT\-5\.4 and Claude Fable 5 were used to debug LaTeX compilation issues, and audit a lengthy first draft of the Introduction, for grammar corrections in the Related Works\. To audit, prompts asked for sentence\-by\-sentence draft evaluation, to determine where content might be redundant\. No text was generated a priori\. The authors referenced the generated audit information to,by their own efforts, rewrite sentences in a condensed manner to meet the submission page limits\.

## References

- \[1\]S\. Bolivar, S\.\-C\. Huang, and R\. Chen,Analysis of tensor time series, Annual Review of Statistics and Its Application, 13 \(2026\), pp\. 369–398,[doi:10\.1146/annurev\-statistics\-042424\-063308](http://dx.doi.org/10.1146/annurev-statistics-042424-063308)\.
- \[2\]J\. Chang, J\. He, L\. Yang, and Q\. Yao,Modelling matrix time series via a tensor CP\-decomposition, Journal of the Royal Statistical Society Series B: Statistical Methodology, 85 \(2023\), pp\. 127–148,[doi:10\.1093/jrsssb/qkac011](http://dx.doi.org/10.1093/jrsssb/qkac011),[https://doi\.org/10\.1093/jrsssb/qkac011](https://doi.org/10.1093/jrsssb/qkac011)\.
- \[3\]Y\. Chen, S\. Su, and H\. Yang,Convolutional neural network analysis of recurrence plots for anomaly detection, International Journal of Bifurcation and Chaos, 30 \(2020\), p\. 2050002,[doi:10\.1142/S0218127420500029](http://dx.doi.org/10.1142/S0218127420500029),[https://doi\.org/10\.1142/S0218127420500029](https://doi.org/10.1142/S0218127420500029)\.
- \[4\]J\.\-P\. Eckmann, S\. O\. Kamphorst, and D\. Ruelle,Recurrence plots of dynamical systems, Europhysics Letters, 4 \(1987\), p\. 973,[doi:10\.1209/0295\-5075/4/9/004](http://dx.doi.org/10.1209/0295-5075/4/9/004),[https://doi\.org/10\.1209/0295\-5075/4/9/004](https://doi.org/10.1209/0295-5075/4/9/004)\.
- \[5\]M\. Gavish and D\. L\. Donoho,The optimal hard threshold for singular values is4/34/\\sqrt\{3\}, IEEE Transactions on Information Theory, 60 \(2014\), pp\. 5040–5053,[doi:10\.1109/TIT\.2014\.2323359](http://dx.doi.org/10.1109/TIT.2014.2323359)\.
- \[6\]B\. Goswami,A brief introduction to nonlinear time series analysis and recurrence plots, Vibration, 2 \(2019\), pp\. 332–368\.
- \[7\]C\. Gück, C\. M\. A\. Roelofs, and S\. Faulstich,CARE to compare: A real\-world benchmark dataset for early fault detection in wind turbine data, Data, 9 \(2024\),[doi:10\.3390/data9120138](http://dx.doi.org/10.3390/data9120138),[https://www\.mdpi\.com/2306\-5729/9/12/138](https://www.mdpi.com/2306-5729/9/12/138)\.
- \[8\]V\. Jayaram, B\. E\. Usevitch, and O\. M\. Kosheleva,Detection from hyperspectral images compressed using rate distortion and optimization techniques under JPEG2000 part 2, in Proceedings of the 11th IEEE Digital Signal Processing Workshop \(DSP’04\), Taos Ski Valley, New Mexico, USA, Aug\. 2004, pp\. 111–114\.
- \[9\]H\. Kang, T\.\-H\. Lee, and J\. Lee,A graph convolutional network for time series classification using recurrence plots, Applied Intelligence, 55 \(2025\), p\. 972,[doi:10\.1007/s10489\-025\-06841\-3](http://dx.doi.org/10.1007/s10489-025-06841-3),[https://doi\.org/10\.1007/s10489\-025\-06841\-3](https://doi.org/10.1007/s10489-025-06841-3)\.
- \[10\]E\. Keogh and J\. Lin,Clustering of time\-series subsequences is meaningless: Implications for previous and future research, Knowledge and Information Systems, 8 \(2005\), pp\. 154–177,[doi:10\.1007/s10115\-004\-0172\-7](http://dx.doi.org/10.1007/s10115-004-0172-7)\.
- \[11\]J\. Kossaifi, Y\. Panagakis, A\. Anandkumar, and M\. Pantic,TensorLy: Tensor learning in Python, Journal of Machine Learning Research \(JMLR\), 20 \(2019\)\.
- \[12\]H\. Li and M\. Chen,Time series clustering based on normal cloud model and complex network, Applied Soft Computing, 148 \(2023\), p\. 110876,[doi:https://doi\.org/10\.1016/j\.asoc\.2023\.110876](http://dx.doi.org/https://doi.org/10.1016/j.asoc.2023.110876),[https://www\.sciencedirect\.com/science/article/pii/S1568494623008943](https://www.sciencedirect.com/science/article/pii/S1568494623008943)\.
- \[13\]X\. Liu, Y\. Xia, Y\. Liang, J\. Hu, Y\. Wang, L\. Bai, C\. Huang, Z\. Liu, B\. Hooi, and R\. Zimmermann,LargeST: A benchmark dataset for large\-scale traffic forecasting, in Advances in Neural Information Processing Systems \(NeurIPS\) 2023 — Datasets & Benchmarks Track, 2023\.
- \[14\]T\. March, S\. Chapman, and R\. Dendy,Recurrence plot statistics and the effect of embedding, Physica D: Nonlinear Phenomena, 200 \(2005\), pp\. 171–184\.
- \[15\]N\. Marwan, M\. Carmen Romano, M\. Thiel, and J\. Kurths,Recurrence plots for the analysis of complex systems, Physics Reports, 438 \(2007\), pp\. 237–329,[doi:https://doi\.org/10\.1016/j\.physrep\.2006\.11\.001](http://dx.doi.org/https://doi.org/10.1016/j.physrep.2006.11.001),[https://www\.sciencedirect\.com/science/article/pii/S0370157306004066](https://www.sciencedirect.com/science/article/pii/S0370157306004066)\.
- \[16\]N\. Marwan and J\. Kurths,Nonlinear analysis of bivariate data with cross recurrence plots, Physics Letters A, 302 \(2002\), pp\. 299–307,[doi:https://doi\.org/10\.1016/S0375\-9601\(02\)01170\-2](http://dx.doi.org/https://doi.org/10.1016/S0375-9601(02)01170-2),[https://www\.sciencedirect\.com/science/article/pii/S0375960102011702](https://www.sciencedirect.com/science/article/pii/S0375960102011702)\.
- \[17\]A\. Mueen, S\. Zhong, Y\. Zhu, M\. Yeh, K\. Kamgar, K\. Viswanathan, C\. Gupta, and E\. Keogh,The fastest similarity search algorithm for time series subsequences under Euclidean distance, August 2022\.[http://www\.cs\.unm\.edu/˜mueen/FastestSimilaritySearch\.html](http://www.cs.unm.edu/~mueen/FastestSimilaritySearch.html)\.
- \[18\]J\. Paparrizos, F\. Yang, and H\. Li,Bridging the gap: A decade review of time\-series clustering methods\., arXiv preprint arXiv:2412\.20582, \(2024\)\.
- \[19\]T\. Rakthanmanon, B\. Campana, A\. Mueen, G\. Batista, B\. Westover, Q\. Zhu, J\. Zakaria, and E\. Keogh,Addressing big data time series: Mining trillions of time series subsequences under dynamic time warping, ACM Transactions on Knowledge Discovery from Data \(TKDD\), 7 \(2013\), pp\. 1–31\.
- \[20\]M\. C\. Romano, M\. Thiel, J\. Kurths, and W\. von Bloh,Multivariate recurrence plots, Physics Letters A, 330 \(2004\), pp\. 214–223,[https://api\.semanticscholar\.org/CorpusID:5746162](https://api.semanticscholar.org/CorpusID:5746162)\.
- \[21\]K\. Sato and S\. Ono,Joint background\-anomaly\-noise decomposition for robust hyperspectral anomaly detection via constrained convex optimization, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 19 \(2026\), pp\. 20217–20235,[doi:10\.1109/JSTARS\.2026\.3698930](http://dx.doi.org/10.1109/JSTARS.2026.3698930)\.
- \[22\]M\. Shahcheraghi, R\. Mercer, J\. M\. D\. A\. Rodrigues, A\. Der, H\. F\. S\. Gamboa, Z\. Zimmerman, K\. Mauck, and E\. J\. Keogh,Introducing Mplots: scaling time series recurrence plots to massive datasets, J\. Big Data, 11 \(2024\), p\. 96,[https://doi\.org/10\.1186/s40537\-024\-00954\-1](https://doi.org/10.1186/s40537-024-00954-1)\.
- \[23\]N\. D\. Sidiropoulos, L\. De Lathauwer, X\. Fu, K\. Huang, E\. E\. Papalexakis, and C\. Faloutsos,Tensor decomposition for signal processing and machine learning, IEEE Transactions on Signal Processing, 65 \(2017\), pp\. 3551–3582,[https://ieeexplore\.ieee\.org/abstract/document/7891546/](https://ieeexplore.ieee.org/abstract/document/7891546/)\.
- \[24\]D\. F\. Silva, V\. M\. D\. Souza, and G\. E\. Batista,Time series classification using compression distance of recurrence plots, in 2013 IEEE 13th International Conference on Data Mining, 2013, pp\. 687–696,[doi:10\.1109/ICDM\.2013\.128](http://dx.doi.org/10.1109/ICDM.2013.128)\.
- \[25\]A\. Trindade,ElectricityLoadDiagrams20112014\.[https://archive\.ics\.uci\.edu/dataset/321/electricityloaddiagrams20112014](https://archive.ics.uci.edu/dataset/321/electricityloaddiagrams20112014)\.UCI Machine Learning Repository\.
- \[26\]W\. Uhlmann,Vergleich der hypergeometrischen mit der Binomial\-Verteilung, Metrika, 10 \(1966\), pp\. 145–158,[doi:10\.1007/BF02613425](http://dx.doi.org/10.1007/BF02613425)\.
- \[27\]C\.\-C\. M\. Yeh, Y\. Zhu, H\. A\. Dau, A\. Darvishzadeh, M\. Noskov, and E\. Keogh,Online amnestic DTW to allow real\-time golden batch monitoring, in Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD ’19, New York, NY, USA, 2019, Association for Computing Machinery, p\. 2604–2612,[https://doi\.org/10\.1145/3292500\.3330650](https://doi.org/10.1145/3292500.3330650)\.
- \[28\]C\. M\. Yeh, Y\. Zhu, L\. Ulanova, N\. Begum, Y\. Ding, H\. A\. Dau, D\. F\. Silva, A\. Mueen, and E\. J\. Keogh,Matrix Profile I: all pairs similarity joins for time series: A unifying view that includes motifs, discords and shapelets, in IEEE 16th International Conference on Data Mining, ICDM 2016, December 12\-15, 2016, Barcelona, Spain, F\. Bonchi, J\. Domingo\-Ferrer, R\. Baeza\-Yates, Z\. Zhou, and X\. Wu, eds\., IEEE Computer Society, 2016, pp\. 1317–1322,[https://doi\.org/10\.1109/ICDM\.2016\.0179](https://doi.org/10.1109/ICDM.2016.0179)\.
- \[29\]Z\. Zhou, X\. Li, J\. Wright, E\. Candès, and Y\. Ma,Stable principal component pursuit, in 2010 IEEE International Symposium on Information Theory, ISIT 2010 \- Proceedings, IEEE International Symposium on Information Theory \- Proceedings, 2010, pp\. 1518–1522,[doi:10\.1109/ISIT\.2010\.5513535](http://dx.doi.org/10.1109/ISIT.2010.5513535)\.2010 IEEE International Symposium on Information Theory, ISIT 2010 ; Conference date: 13\-06\-2010 Through 18\-06\-2010\.
- \[30\]Z\. Zimmerman, K\. Kamgar, N\. S\. Senobari, B\. Crites, G\. J\. Funning, P\. Brisk, and E\. J\. Keogh,Matrix Profile XIV: scaling time series motif discovery with GPUs to break a quintillion pairwise comparisons a day and beyond, in Proceedings of the ACM Symposium on Cloud Computing, SoCC 2019, Santa Cruz, CA, USA, November 20\-23, 2019, ACM, 2019, pp\. 74–86,[https://doi\.org/10\.1145/3357223\.3362721](https://doi.org/10.1145/3357223.3362721)\.

## Appendices

## ACode

## BBounding Chance Pass Rates

For a time seriesTTofNNlength\-nnsensor traces, letmmdenote subsequence length andckc\_\{k\}the Mplot pixel width\. Suppose the random\-selection model holds forTT\. In the single\-trial procedure, a set of chosen windows𝒲\\mathcal\{W\}is constructed with coveragef=m​\|𝒲\|nf=\\frac\{m\|\\mathcal\{W\}\|\}\{n\}, while fractions ofpb,pb,and​1−2​pbp\_\{b\},p\_\{b\},\\text\{and \}1\-2p\_\{b\}of the time series are transformed byθ1,θ2,and​θ3\\theta\_\{1\},\\theta\_\{2\},\\text\{ and \}\\theta\_\{3\}respectively\. From each of theR∈\[2\.\.8\]R\\in\[2\.\.8\]components of the tensor decomposition, the top\-sssensors and top\-ssintervals are obtained\. The experimental constants in Section[5\.4](https://arxiv.org/html/2608.04157#S5.SS4)fixck≈m/3c\_\{k\}\\approx m/3,f≤0\.10f\\leq 0\.10,pb=0\.25p\_\{b\}=0\.25, ands=6s=6\. For clarity we assumepb​Np\_\{b\}Nis an integer so thatb=pb​Nb=p\_\{b\}Nexactly\.

We now establish bounds on the four chance pass rates\.

###### Proposition B\.1\.

In the random\-selection model,

p1chance≤0\.00584\.p\_\{1\}^\{\\text\{chance\}\}\\leq 0\.00584\.

###### Proof 1\.

Prior to tensor decomposition,w=n⋅fmw=\\frac\{n\\cdot f\}\{m\}𝒲\\mathcal\{W\}\-windows are planted\. An interval of lengthl=m\+2​ckl=m\+2c\_\{k\}intersects a𝒲\\mathcal\{W\}\-window\[i\.\.i\+m\)\[i\.\.i\+m\)if its left endpoint lies in\(i−l\.\.i\+m\)\(i\-l\.\.i\+m\), so it intersects one of thewwwindows with probability at most

w⋅\(m\+l\)n=n​fm⋅2​m\+2​ckn=2​f⋅\(1\+ckm\)\.w\\cdot\\frac\{\(m\+l\)\}\{n\}=\\frac\{nf\}\{m\}\\cdot\\frac\{2m\+2c\_\{k\}\}\{n\}=2f\\cdot\(1\+\\frac\{c\_\{k\}\}\{m\}\)\.For a given component, all top\-ssintervals intersect𝒲\\mathcal\{W\}\-windows with probability at most\(2​f⋅\(1\+ckm\)\)s\(2f\\cdot\(1\+\\frac\{c\_\{k\}\}\{m\}\)\)^\{s\}\. Taking a union bound overRRcomponents, we yield

p1chance≤R​\(2​f⋅\(1\+ckm\)\)s\.p\_\{1\}^\{\\text\{chance\}\}\\leq R\\left\(2f\\cdot\(1\+\\frac\{c\_\{k\}\}\{m\}\)\\right\)^\{s\}\.Withs=6s=6,f≤0\.10f\\leq 0\.10,ckm≤1/2\\frac\{c\_\{k\}\}\{m\}\\leq 1/2, andR≤8R\\leq 8we obtain the desired inequality\.

###### Proposition B\.2\.

In the random\-selection model,

p2chance≤0\.176\.p\_\{2\}^\{\\text\{chance\}\}\\leq 0\.176\.

###### Proof 2\.

LetSϕ,iS\_\{\\phi,i\}denote the number of top\-sssensors in componentϕ\\phitransformed byθi\\theta\_\{i\}\. Then under the random\-selection model

Sϕ,i∼\{Hypergeometric​\(N,pb​N,s\),i=1,2Hypergeometric​\(N,\(1−2​pb\)​N,s\),i=3\.S\_\{\\phi,i\}\\sim\\begin\{cases\}\\text\{Hypergeometric\}\(N,p\_\{b\}N,s\),&i=1,2\\\\ \\text\{Hypergeometric\}\(N,\(1\-2p\_\{b\}\)N,s\),&i=3\\end\{cases\}\.We also define the auxiliary variable

Bϕ,i∼\{Binomial​\(s,pb\),i=1,2Binomial​\(s,1−2​pb\),i=3\.B\_\{\\phi,i\}\\sim\\begin\{cases\}\\text\{Binomial\}\(s,p\_\{b\}\),&i=1,2\\\\ \\text\{Binomial\}\(s,1\-2p\_\{b\}\),&i=3\\end\{cases\}\.Using independence of components and observing that only one ofθ1,θ2,θ3\\theta\_\{1\},\\theta\_\{2\},\\theta\_\{3\}can transform more than half ofϕ\\phi’s top\-sssensors, we find that

p2chance=Pr⁡\(∀ϕ∈\[1\.\.R\],∃i∈\{1,2,3\},Sϕ,i\>s2\)p\_\{2\}^\{\\text\{chance\}\}=\\Pr\(\\forall\_\{\\phi\\in\[1\.\.R\]\}\\,,\\exists\_\{i\\in\\\{1,2,3\\\}\}\\,,S\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)=∏ϕ=1R∑i=13Pr⁡\(Sϕ,i\>s2\)\.=\\prod\_\{\\phi=1\}^\{R\}\\sum\_\{i=1\}^\{3\}\\Pr\(S\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)\.By Satz 5 of Uhlmann\[[26](https://arxiv.org/html/2608.04157#bib.bib26)\],

∀i∈\{1,2,3\},Pr⁡\(Sϕ,i\>c\)<Pr⁡\(Bϕ,i\>c\)\\forall\_\{i\\in\\\{1,2,3\\\}\}\\,,\\Pr\(S\_\{\\phi,i\}\>c\)<\\Pr\(B\_\{\\phi,i\}\>c\)when

pb,1−2​pb≤cs−1−c\(s−1\)​\(N\+1\)=cs−1​NN\+1\.p\_\{b\},1\-2p\_\{b\}\\leq\\frac\{c\}\{s\-1\}\-\\frac\{c\}\{\(s\-1\)\(N\+1\)\}=\\frac\{c\}\{s\-1\}\\frac\{N\}\{N\+1\}\.These conditions are satisfied forc=s2c=\\frac\{s\}\{2\}in the random\-selection model with the provided experiment constants \(pb=0\.25p\_\{b\}=0\.25ands≤Ns\\leq Nan even number suffice\), so we can upper\-bound by the binomial tail probability\.

∏ϕ=1R∑i=13Pr⁡\(Sϕ,i\>s2\)≤∏ϕ=1R∑i=13Pr⁡\(Bϕ,i\>s2\)\\prod\_\{\\phi=1\}^\{R\}\\sum\_\{i=1\}^\{3\}\\Pr\(S\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)\\leq\\prod\_\{\\phi=1\}^\{R\}\\sum\_\{i=1\}^\{3\}\\Pr\(B\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)=\(∑i=13Pr⁡\(Bϕ,i\>s2\)\)R\.=\\left\(\\sum\_\{i=1\}^\{3\}\\Pr\(B\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)\\right\)^\{R\}\.Withs=6s=6,pb=0\.25p\_\{b\}=0\.25, andR≥2R\\geq 2, we obtain the desired inequality\.

###### Proposition B\.3\.

In the random\-selection model,

p3chance≤0\.325\.p\_\{3\}^\{\\text\{chance\}\}\\leq 0\.325\.

###### Proof 3\.

LetEiE\_\{i\}denote the number of components where there exists a top\-sssensor transformed byθi\\theta\_\{i\}\. Notice that

Ei∼Binomial​\(R,Pr⁡\(Sϕ,i≥1\)\)\.E\_\{i\}\\sim\\text\{Binomial\}\(R,\\Pr\(S\_\{\\phi,i\}\\geq 1\)\)\.Hence,

p3chance=Pr⁡\(E1≤1\)=\(1−Pr⁡\(Sϕ,1≥1\)\)Rp\_\{3\}^\{\\text\{chance\}\}=\\Pr\(E\_\{1\}\\leq 1\)=\(1\-\\Pr\(S\_\{\\phi,1\}\\geq 1\)\)^\{R\}\+R⋅Pr⁡\(Sϕ,1≥1\)​\(1−Pr⁡\(Sϕ,1≥1\)\)R−1,\+R\\cdot\\Pr\(S\_\{\\phi,1\}\\geq 1\)\(1\-\\Pr\(S\_\{\\phi,1\}\\geq 1\)\)^\{R\-1\},which increases asPr⁡\(Sϕ,1≥1\)\\Pr\(S\_\{\\phi,1\}\\geq 1\)decreases\. Moreover

Pr⁡\(Sϕ,1≥1\)=1−Pr⁡\(Sϕ,1=0\)\\Pr\(S\_\{\\phi,1\}\\geq 1\)=1\-\\Pr\(S\_\{\\phi,1\}=0\)=1−\(\(1−pb\)​Ns\)\(Ns\)=1−∏i=0s−1\(1−pb\)​N−iN−i\.=1\-\\frac\{\\binom\{\(1\-p\_\{b\}\)N\}\{s\}\}\{\\binom\{N\}\{s\}\}=1\-\\prod\_\{i=0\}^\{s\-1\}\\frac\{\(1\-p\_\{b\}\)N\-i\}\{N\-i\}\.Every term withi\>0i\>0of the product

∏i=0s−1\(1−pb\)​N−iN−i\\prod\_\{i=0\}^\{s\-1\}\\frac\{\(1\-p\_\{b\}\)N\-i\}\{N\-i\}increases inNN, soPr⁡\(Sϕ,1≥1\)\\Pr\(S\_\{\\phi,1\}\\geq 1\)decreases inNN\. Thus, asNNincreases,p3chancep\_\{3\}^\{\\text\{chance\}\}increases\. In fact,

limN→∞Pr⁡\(Sϕ,1≥1\)=Pr⁡\(Bϕ,1≥1\)=1−\(1−pb\)s,\\lim\_\{N\\to\\infty\}\\Pr\(S\_\{\\phi,1\}\\geq 1\)=\\Pr\(B\_\{\\phi,1\}\\geq 1\)=1\-\(1\-p\_\{b\}\)^\{s\},⟹Pr⁡\(Sϕ,1≥1\)≥1−\(1−pb\)s\.\\implies\\Pr\(S\_\{\\phi,1\}\\geq 1\)\\geq 1\-\(1\-p\_\{b\}\)^\{s\}\.Therefore,

p3chance≤\(1−pb\)s​R\+R⋅\(1−\(1−pb\)s\)​\(1−pb\)s​\(R−1\)\.p\_\{3\}^\{\\text\{chance\}\}\\leq\(1\-p\_\{b\}\)^\{sR\}\+R\\cdot\(1\-\(1\-p\_\{b\}\)^\{s\}\)\(1\-p\_\{b\}\)^\{s\(R\-1\)\}\.Withpb=0\.25p\_\{b\}=0\.25,s=6s=6, andR≥2R\\geq 2, we obtain the desired inequality\.

RankRRUpper bound forp4chancep\_\{4\}^\{\\text\{chance\}\}20\.57030\.27440\.42650\.22660\.11170\.19080\.0991Table B\.1:Upper bounds forp4chancep\_\{4\}^\{\\text\{chance\}\}as a function of the selected rankRR\.###### Proposition B\.4\.

In the random\-selection model,

p4chance≤0\.570\.p\_\{4\}^\{\\text\{chance\}\}\\leq 0\.570\.

###### Proof 4\.

LetCiC\_\{i\}denote the number of components thatθi\\theta\_\{i\}corresponds to\. Notice

Ci∼Binomial​\(R,Pr⁡\(Sϕ,i\>s2\)\)\.C\_\{i\}\\sim\\text\{Binomial\}\(R,\\Pr\(S\_\{\\phi,i\}\>\\frac\{s\}\{2\}\)\)\.Letps:=Pr⁡\(Sϕ,3\>s2\)p\_\{s\}:=\\Pr\(S\_\{\\phi,3\}\>\\frac\{s\}\{2\}\)\. Then

p4chance=Pr⁡\(C3≥⌊0\.7​R⌋\)=∑r=⌊0\.7​R⌋RPr⁡\(C3=r\)p\_\{4\}^\{\\text\{chance\}\}=\\Pr\(C\_\{3\}\\geq\\lfloor 0\.7R\\rfloor\)=\\sum\_\{r=\\lfloor 0\.7R\\rfloor\}^\{R\}\\Pr\(C\_\{3\}=r\)=∑r=⌊0\.7​R⌋R\(Rr\)​\(ps\)r​\(1−ps\)R−r\.=\\sum\_\{r=\\lfloor 0\.7R\\rfloor\}^\{R\}\\binom\{R\}\{r\}\\left\(p\_\{s\}\\right\)^\{r\}\\left\(1\-p\_\{s\}\\right\)^\{R\-r\}\.Moreover,

∂∂ps​p4chance=\(R⌊0\.7​R⌋\)​⌊0\.7​R⌋​ps⌊0\.7​R⌋−1​\(1−ps\)R−⌊0\.7​R⌋,\\frac\{\\partial\}\{\\partial p\_\{s\}\}p\_\{4\}^\{\\text\{chance\}\}=\\binom\{R\}\{\\lfloor 0\.7R\\rfloor\}\\lfloor 0\.7R\\rfloor p\_\{s\}^\{\\lfloor 0\.7R\\rfloor\-1\}\(1\-p\_\{s\}\)^\{R\-\\lfloor 0\.7R\\rfloor\},which is positive, sop4chancep\_\{4\}^\{\\text\{chance\}\}increases withpsp\_\{s\}\. Now, by the Uhlmann bound

ps≤Pr⁡\(Bϕ,3\>s2\),p\_\{s\}\\leq\\Pr\(B\_\{\\phi,3\}\>\\frac\{s\}\{2\}\),so withs=6s=6and1−2​pb=0\.51\-2p\_\{b\}=0\.5, we can upper\-boundp4chancep\_\{4\}^\{\\text\{chance\}\}for2≤R≤82\\leq R\\leq 8\. The calculated values are listed in Table[B\.1](https://arxiv.org/html/2608.04157#S2.T1)\. Taking the maximum upper bound over all ranks thus yields the desired inequality\.