Fast Algorithm for computing a matrix transform used to detect trends in noisy data
AA fast algorithm for computing a matrix transform usedto detect trends in noisy data
Dan Kestner a , Glenn Ierley b , Alex Kostinski a, ∗ a Physics Department, Michigan Technological University, 1400 Townsend Dr. HoughtonMI. 49931 b Mathematics Department, Michigan Technological University, 1400 Townsend Dr.Houghton MI. 49931
Abstract
A recently discovered universal rank-based matrix method to extract trendsfrom noisy time series is described in [1] but the formula for the outputmatrix elements, implemented there as an open-access supplement MATLABcomputer code, is O ( N ), with N the matrix dimension. This can becomeprohibitively large for time series with hundreds of sample points or more.Based on recurrence relations, here we derive a much faster O ( N ) algorithmand provide code implementations in MATLAB and in open-source JULIA.In some cases one has the output matrix and needs to solve an inverse problemto obtain the input matrix. A fast algorithm and code for this companionproblem, also based on the above recurrence relations, are given. Finally, inthe narrower, but common, domains of (i) trend detection and (ii) parameterestimation of a linear trend, users require, not the individual matrix elements,but simply their accumulated mean value. For this latter case we provide ayet faster O ( N ) heuristic approximation that relies on a series of rank one ∗ Corresponding author.
E-mail address: [email protected]
Preprint submitted to Computer Physics Communications January 28, 2020 a r X i v : . [ phy s i c s . d a t a - a n ] J a n atrices. These algorithms are illustrated on a time series of high energycosmic rays with N > × . Keywords:
Time-Series; Trend; Noise; Rank; Complexity Reduction
PROGRAM SUMMARY
Program Title:
Pfromdata, QofP, mbasisandcoeffs, nonzerop, Qavgapprox, PofQ,Testing
Program Files doi: http://dx.doi.org/xx.xxxxx/xxxxx.x (to be assigned by jour-nal)
Licensing provisions:
MIT (Julia)
Programming language:
MATLAB and Julia
Nature of problem:
An order-rank data matrix and its transform to a stable formare used repeatedly to detect and/or extract trends from noisy data. An efficientyet accurate calculation of the matrix transform is therefore required.
Solution method:
We introduce and apply an analytic recursion relation, whichspeeds up the execution of the matrix transform from O ( N ) arithmetic opera-tions to O ( N ). Since this matrix transform is called often during optimization,our improvement allows for far shorter optimization times, for a given sample size.For example, a transform whose time is extrapolated to an unrealistic 75 days on aDell personal laptop computer with a 2.2 GHz quad-core AMD processor running32 bit MATLAB version R2015b on 64 bit Windows 10 ( N = 5000), now takes afraction of a second. References [1] Universal Rank-Order Transform to Extract Signals from Noisy Data, GlennIerley and Alex Kostinski, Phys. Rev. X 9 031039 (2019) . Introduction A broadly-applicable rank-based approach for detection and extraction ofgenerally non-linear trends in noisy time series has recently been introduced[1] and we shall now briefly review the mathematical essentials. The inputtime series is segmented into n t samples, with each sample having n T datapoints. A square n T × n T population matrix P is then calculated such that P j,k is the population (number) of data points with order j (position inthe sample), and rank k (position in the sample after an ascending sort)[1].Alternatively, P j,k can also be viewed as a 2D probability density function(pdf) or a histogram over the plane defined by rank and times axes. Thematrix P is illustrated below in (1). P = P , P , . . . P ,k − P ,k P ,k +1 . . . P ,n T P , P , . . . P ,k − P ,k P ,k +1 . . . P ,n T ... ... ... ... ... ... P j − , P j − , . . . P j − ,k − P j − ,k P j − ,k +1 . . . P j,n T P j, P j, . . . P j,k − P j,k P j,k +1 . . . P j,n T P j +1 , P j +1 , . . . P j +1 ,k − P j +1 ,k P j +1 ,k +1 . . . P k +1 ,n T ... ... ... ... ... ... P n T , P n T , . . . P n T ,k − P n T ,k P n T ,k +1 . . . P n T ,n T (1)To “zoom in” on the trends hidden in P , the Q -transform was introduced[1] as follows 3 j,k = (cid:32) (cid:80) jm =1 (cid:80) kn =1 P m,n + (cid:80) n T m = j +1 (cid:80) n T n = k +1 P m,n jk + ( n T − j )( n T − k ) − (cid:80) jm =1 (cid:80) n T n = k +1 P m,n + (cid:80) n T m = j +1 (cid:80) kn =1 P m,n ( n T − j ) k + j ( n T − k ) (cid:33) n T n t (2)To understand the construction, consider the division of P into quadrantsfor calculation of Q j,k as shown on the RHS of (1). Each element of Q isthe difference between the average matrix element of the combined upperleft and lower right quadrants of P , and the average matrix element of thecombined upper right and lower left quadrants, normalized by the overallaverage matrix element of P , (cid:104) P (cid:105) ≡ (cid:80) n T m =1 (cid:80) n T n =1 P m,n /n T = n t /n T .The number of operations (+ , − , × , ÷ ) required to compute Q using equa-tion (2) is O ( n T ). For large n T and repeated calls, as will often be neededin applications, the computation time can become prohibitively long. In ad-dition, setting (cid:104) Q (cid:105) = 0 where angular brackets denote average over matrixelements, functions as a trend detector when the functional form of the trendis not available and we shall illustrate it on the time series of cosmic rays in5. To that end, our purpose in this paper is four-fold:(i) present a O ( n T ) algorithm for computing the Q -transform and itsMATLAB implementation;(ii) supply an open source (Julia) implementation;(iii) present an efficient calculation of (cid:104) Q (cid:105) , where (cid:104) Q (cid:105) is the average matrixelement of Q . The departure of this (scalar) quantity from zero is used todetect presence of trend[1].(iv) provide an illustrative example from a long cosmic ray time series;4 . Derivation of the Algorithm To begin, note that equation (2) can be simplified by making use ofconstraints on P that each row and column sum to n t . Thus, the sums ofelements in the four quadrants of P, entering the numerator of (2) are notindependent. Numbering the quadrants as 1-4 beginning from the upperright, moving counter-clockwise, and calling the sums of elements in eachquadrant i as Σ ( i ) , we haveΣ (1) + Σ (2) = jn t Σ (2) + Σ (3) = kn t Σ (3) + Σ (4) = ( n T − j ) n t Σ (4) + Σ (1) = ( n T − k ) n t (3)This system of four equations in four unknowns (Σ ( i ) , i = 1 −
4) is under-determined and when recast as a 4x4 matrix equation, has a matrix rank ofthree. Thus, only one of the four Σ ( i ) is independent and we picked Σ (2) forthat purpose. Σ (1) = jn t − Σ (2) Σ (2) = Σ (2) Σ (3) = kn t − Σ (2) Σ (4) = ( n T − j − k ) n t + Σ (2) (4)This can be substituted back into equation (2), making explicit a ( j, k )index on Σ ( i ) . 5 j,k = (cid:32) ( n T − j − k ) n t + 2Σ (2) j,k jk + ( n T − j )( n T − k ) − ( j + k ) n t − (2) j,k ( n T − j ) k + j ( n T − k ) (cid:33) n T n t (5)Define D as a ( n T − × ( n T −
1) matrix, whose elements are the productof the two denominators in equation (2): D j,k = ( jk + ( n T − j )( n T − k ))( j ( n T − k ) + ( n T − j ) k ) (6)The Q matrix can be expressed compactly in terms of Σ (2) j,k and D j,k . Q j,k D j,k = 2 n T n t (cid:32) Σ (2) j,k − n t n T jk (cid:33) (7)The motivation for this is that the second quadrant sum Σ (2) j,k satisfies arecurrence relation.Σ (2) j,k = Σ (2) j − ,k + Σ (2) j,k − − Σ (2) j − ,k − + P j,k (8)Taken together with equation (7), this yields a recurrence relation for Q . Q j,k = 1 D j,k (cid:32) D j,k − Q j,k − + D j − ,k Q j − ,k − D j − ,k − Q j − ,k − + 2 n T n t ( P j,k − n t n T ) (cid:33) (9)The algorithm used to calculate Q via (9) is described in Algorithm 2and its MATLAB and Julia implementations accompany this manuscript.The full Q matrix is calculable in O ( n T ) operations, as is seen by observing6hat each element of Q can be calculated in O (1) from a small number ofneighboring Q elements and some constants, and that the total number ofelements in Q is ( n T − . Algorithm 1: An O ( n T ) implementation of the Q -transform, using re-currence. The (1,1) element is found first. Rows and columns are foundby moving rightwards or downwards from diagonal elements, which arefound from previous rows/columns. Data: n T × n T population matrix P Result: ( n T − × ( n T −
1) matrix Q First, verify that input matrix P satisfies row and column sumconstraints;Calculate (1 ,
1) element using equation (2);Evaluate the remainder of the first row of Q using recursion. Zero ornegative indices of Q in equation (9) mean Q may be replaced byzero there, as a more careful consideration of values of Q near thematrix edge confirms;Repeat for the remainder of the first column of Q ;Calculate (2,2) from surrounding three elements;Evaluate remainder of row/column 2 of Q recursively;...Calculate ( n, n ) element of Q from upper left surrounding elements;Evaluate remainder of row/column n of Q , using recursion;...Calculate ( n T − , n T −
1) element7 . Analytical Results
One key result of this paper is equation (9), just derived. This permitsan O ( n T ) method for calculating Q that is much faster than the O ( n T ) bruteforce evaluation of equation (2), especially for large n T . Another essentialresult is the transformation for P , given Q . This was obtained by rearrangingequation (9) as follows. P j,k = (cid:32) D j,k Q j,k − D j,k − Q j,k − − D j − ,k Q j − ,k + D j − ,k − Q j − ,k − (cid:33) n t n T + n t n T (10)Not only does this transformation turn out to be stably computable butalso efficiently so. In fact, it can be accomplished also in O ( n T ) operationsand is implemented in MATLAB and Julia programs in the accompanyingfiles. These results allow analyses of previously inaccessible data because ofthe prohibitively long computation times. The confirmation of the speed upof equation (9) over equation (2) directly in terms of CPU time is given inFig. 1 below.As a sample application, possible because of the computational improve-ment provided by calculating the Q -transform recursively rather than by thedirect evaluation of the double sums in equation (2), we choose n T = 5000 ≈ . , which lies just outside of the axis range shown in Fig. 1. A calculationusing the recursive result in equation (9) takes about 0.7 seconds[2]. In com-parison, using Fig. 1 to extrapolate the O ( n T ) curve out to log ( n T ) = 12 . igure 1: A comparison of calculation times for the Q-transform. A recursive approachreduces the CPU time from O ( n T ) to O ( n T ). The blue curve shows results from animplementation of the brute force approach of equation (2), while the red curve showstiming data from the optimized calculation of equation (9). The contrast in complexityorder becomes apparent for larger size matrices, log ( n T ) (cid:39)
4. Note that the calculationtime becomes intractable for n T (cid:39) ≈ shown in the figure.
4. Fast Algorithm for Calculating (cid:104) Q (cid:105) We now turn to efficient calculation of (cid:104) Q (cid:105) , the mean matrix element of Q in the special case of large n T and small n t . In such single sample (timeseries) cases, the matrix P is sparse, consisting of N T − N T zeroes and neednot be stored in memory in its entirety. Rather, only indices of the non-zero matrix elements, found by independently sorting each of the repeated n t trials, may suffice to calculate Q . Also, it was shown in [1] that theentire Q information is not required when one is concerned merely with trend9etection or parameter estimation of a linear trend. One such application(time series of cosmic ray arrivals) is discussed in the next section. In suchcases, it suffices to calculate only element-averaged metric (cid:104) Q (cid:105) ≡ n T − n T − (cid:88) j,k =1 Q j,k (11)A sufficiently large value of (cid:104) Q (cid:105) implies the presence of a trend. Here,”sufficiently large” is with reference to a fiducial value expected for pure noise,a formula for which in terms of n T and n t is given in [1]. Using the resultsof Section 3, Q , and thus (cid:104) Q (cid:105) (which has ( n T − terms), can be computedfrom P in O ( n T ) operations. In order to calculate (cid:104) Q (cid:105) , only accumulatedmean value is needed and the rest of the matrix elements need not be stored,thereby reducing complexity of the calculation. To that end, it is shown in [1]that (cid:104) Q (cid:105) can be written as the sum of elements of the Hadamard product of a( n T × n T ) matrix m and P itself (Appendix B.3 of [1]). The matrix m is notdetermined by the data, and depends solely on n t and n T , the former beingonly an inverse multiplicative constant. Once m is constructed, the meanelement of Q is easily accessed as (cid:104) Q (cid:105) = (cid:80) n T j =1 (cid:80) n T k =1 m j,k P j,k . For sparse P , n t (cid:28) nT , this means a calculation of order n T , much faster than the O ( n T )needed to calculate the matrix Q explicitly prior to averaging. As it stands, m takes O ( n T ) to construct, which can become prohibitive for large n T . Notonly this, but m is memory limited to about n T = √ × for an 8GBRAM, being of type double (8 bytes per element). Since for n t (cid:28) n T thematrix P is sparse, in this case we only need calculate the elements of m corresponding to nonzeros in P . This greatly reduces demands on memory,allowing n T , the number of data points per trial, to be as large as about 10 n t = 1 and a 8 GB RAM. Also, there is an O ( n T ) way to approximateany element of m in O (1), operations, reducing the calculation of m forsparse P to O ( n T ). For nonsparse P , when the full matrix m is needed, theapproximation scheme gives m in O ( n T ), due the number of elements needed.We also find that m may be calculated exactly, up to accumulated roundingerrors due to finite machine precision, by an O ( n T ) recursive algorithm, muchlike for Q in section 2. The advantage in this case is that m need only becalculated once, for each n T , and then it applies to any dataset of the samesize parameter n T , modulo a rescaling due to varying n t . This saves thetime needed for calculating Q itself each time. The disadvantage of usingrecursion to find m as compared with using the approximate approach is thatrecursion can only create contiguous rectangular blocks of m . In contrast,the approximation for m allows only the needed elements to be calculated,regardless of how they are spatially related in the matrix. To this we nowturn.To describe how the matrix m is approximated, which is the most efficientway to calculate (cid:104) Q (cid:105) for large n T and n t (cid:28) n T , we first note the following: therank one matrix that is an outer product of the first column of m with itself,normalized by 1 /m (1 , m . Seeing that this approximation is rank one suggests the approximationcan be improved by adding another rank one term. This turns out to beso. One simply takes a linear combination of the first two columns of m ,and since the first row and column of m are already exact, chooses thiscombination such that a new vector is obtained with vanishing first element.The outer product of this vector with itself is zero along the first row and11olumn, and thus does not alter the previous exact rank one approximationthere. If one then adds this new rank one matrix to the old, while choosing ascalar coefficient to match any of the elements in the original second columnof m , one obtains a rank two approximation of m that is exact in the first two rows and columns. This holds true for any symmetric matrix, as a littlealgebra can show. Moreover, this extends easily: a third ”basis vector” maybe obtained with vanishing first and second elements by judiciously mixingthe first three exact columns of m . Again adding the outer product of thisnew vector with itself to the rank 2 approximation, with an appropriateconstant chosen to match an arbitrary element of the exact third column of m , a rank 3 approximation is obtained that is exact in the first three rowsand columns. And so on. By generalization, beginning with r columns of m yields a rank r approximation, exact in the first r rows and columns of m . Since m always has the same form when regarded as a two dimensionalfunction, the effect of increasing the number of rows/columns n T is to bringnearby columns closer numerically. Therefore, practically, difficulties arisewith the above approximation scheme due to nearby columns of m becominglinearly dependent for large n T , but these can be circumvented by avoidingsuccessive columns, but picking columns increasingly separated with n T , soas to roughly maintain proportionate horizontal locations in the matrix m . Italso proves advantageous to mix the columns such that the zeros in the newcolumn vectors are also spaced out proportionately within the matrix and notmerely adjacent and at the beginning. Heuristically, the optimal case is whenthe zeros in the new column vectors have the same spacing as the columns.This improves the conditioning of a certain matrix that must be inverted12n this process. Also, we find that for uniform column spacing, there is anoptimal rank approximation of r = 6. In this case, the optimal column/rowzero spacings s are given empirically by s = [0 . . × n T ].Taking non-uniformly spaced columns of the matrix m yields generally muchbetter results, as found for example by using MATLAB’s Genetic Algorithmto find the optimal columns of m for a rank r = 6 expansion with row zerospacings matching the column spacings. We also choose the coefficients ofthe rank 1 matrices in order to optimize the approximation of m along itsdiagonal, via least squares (MATLAB’s backslash operator).We note that this approximation is essentially a low rank matrix approx-imation that uses low rank matrix completion, topics that both arise in datascience[3].With this approximation of m , we now show that the cost of computing m is reduced from O ( n T ) to O ( n T ) for the overhead, and O (1) operationsfor each element after that. We also show that the memory overhead is also O ( n T ), plus the cost of each element computed after that (between O ( n T )and O ( n T )).The equation for m is now derived. As shown in [1], unwrapping matrices P and Q to vectors allows to write Q = M P , where M is a matrix withdimensions ( n T − × n T . m is then defined by m = 1 T M = [111 ... M andhas dimension 1 × n T . In this paper, by comparing with equation (2) andcarefully converting between matrix indices and linear vector indices, we finda closed form solution for m rewrapped as a n T × n T matrix, and satisfying (cid:104) Q (cid:105) = (cid:80) j,k ( m ◦ P ) j,k , where ◦ denotes the Hadamard product.13 igure 2: O ( n T ) vs O ( n T ) calculations of (cid:104) Q (cid:105) ≡ average over all elements of Q , in thelimit n t << n T . The relative difference δ ≡ ( < Q > approx. − < Q > exact ) / < Q > exact isshown in the inset. The linear time includes the time to ”precompute” the needed subsetof elements of the matrix m (see the text) and take the Hadamard product with P . Thequadratic time includes the time to calculate, store, and average the full Q matrix. m j,k = n T n t ( n T − (cid:18) j − (cid:88) m =1 k − (cid:88) n =1 d (1) m,n + n T − (cid:88) m = j n T − (cid:88) n = k d (1) m,n − j − (cid:88) m =1 n T − (cid:88) n = k d (2) m,n − n T − (cid:88) m = j k − (cid:88) n =1 d (2) m,n (cid:19) (12) d (1) m,n = 1 mn + ( n T − m )( n T − n ) d (2) m,n = 1 m ( n T − n ) + ( n T − m ) n (13)Since only the non-zero entries of P contribute to the element-wise prod-uct with m , and P can have as few as O ( n T ) non-zero entries (when n t = 1),14 Q (cid:105) can be computed in as low as O ( n T ) operations, after overhead that isalso O ( n T ), leaving a grand total of O ( n T ) operations. This is in contrastto the O ( n T ) operations it would take to compute Q using the recursive al-gorithm of section 2, and then sum and average the elements of Q . Thisdifference is illustrated in figure 2.
5. Illustration on Time Series of Cosmic Rays
To illustrate the importance of the numerical acceleration for trend de-tection as just described in the previous section, we pick an example fromcosmic ray physics. The data consists of 49,223 events (only 1% of the totaldata is available to general public), in a form of a time series of arrivals withvarious energies (see Fig. 3 from data in [4]).Energy-resolved flux (spectrum) plays the central role in the field andit is universally assumed that the underlying time series are statisticallystationary. Are they? Here we ask whether the time series in Fig.3 arestationary and we use (cid:104) Q (cid:105) to test the hypothesis. Stationarity implies that (cid:104) Q (cid:105) ≈ (cid:104) Q (cid:105) = 0 via the asymptotics in equation (10) ofreference [1]. Calculation of the auto-correlation function for this cosmicray data shows that it is uncorrelated (“white”) so using (cid:104) Q (cid:105) -asymptotics isparticularly relevant.For sake of consistency, we tested a variety of data partitioning, but withthe same product n t n T . Table 1 shows the importance of the approximate (cid:104) Q (cid:105) algorithm. For n t approaching unity, dimensions of Q are ∼ × and the O ( N ) algorithm is crucial. Table 2 shows the calculations. To15 igure 3: Time series of energies of high energy cosmic rays. There are a total of 49223events between 2004 and 2019. Energy scale is 10 eV = 1EeV. The data is uncorrelated(white), and has a non Gaussian distribution. Data is taken from the Pierre Auger Obser-vatory Public Event Explorer [4], and represents 1% of all data taken by the observatory. our surprise, the (cid:104) Q (cid:105) -test consistently detects a presence of a trend beyondreasonable doubt. Specifically, (cid:104) Q (cid:105) = 0 .
06 gives the confidence limit of 19 σ (taking the case n t (cid:39)
100 for specificity). The associated linear trend is largeenough to affect the spectrum and cast doubt on the traditional power-lawanalysis as the latter implies stationarity via the Wiener-Khintchin theorem.
6. Concluding Remarks
In conclusion, we have discovered an O ( N ) calculation of a previously O ( N ) matrix transform with applications in trend detection from noisy data.This increases the efficiency of the transform, and allows access to previouslyout-of-reach data sample lengths N . For the special case of a small number16 T Time Q (Sec) Time Q sum (Sec) Time m Basis andCoefficients (Sec) Time m i,j P i,j sum (Sec)49140 - - 0.968 1.58E-0224570 - - 0.484 1.15E-029828 31.8 8.70E-02 0.200 1.14E-024095 2.85 1.53E-02 8.39E-02 7.37E-03468 1.91E-02 3.95E-04 1.05E-02 7.02E-0391 5.50E-04 2.94E-05 2.73E-03 7.69E-0312 4.99E-05 2.34E-05 1.22E-03 7.11E-03 Table 1: To investigate the range of possibilities, we compare several partitions of thetime-series of cosmic rays. By trimming the data from 49223 down to 49140 datapoints,the number of distinct integers n t , n T such that n t n T equals the total number of keptdata points is maximized, providing many partitions for study. The product n t n T is heldconstant. In this table, timing dependence of partitions is shown. For all entries shown, n t = 49140 /n T is integer. One can see that the fast approximate calculation for (cid:104) Q (cid:105) is possible for all possible partitions, while those based on full calculation of Q becomememory limited past about n T = 10 , as indicated by the dashes. The approximatecalculation stores neither P nor Q , and is able to be performed up to the maximum n T .Thus, the approximate method is the only way to probe these partitions. For n T greaterthan a few hundred, the approximate calculation, consisting of approximating basis vectorsand coefficients for m , followed by a Hadamard product with P , is much faster than directevaluation of Q and its mean element. of samples n t , we present also an O ( N ) calculation of trend detection metric (cid:104) Q (cid:105) which bypasses the need to carry out the full Q -transform. Open accesscomputer codes are provided for both of these calculations.17 T (cid:104) Q (cid:105) exact (cid:104) Q (cid:105) approximate δ (cid:104) Q (cid:105) /σ white noise Table 2: Companion to Table 1 for (cid:104) Q (cid:105) calculation values under different data partitions.Both exact and approximate calculations give results that are roughly independent of datapartition for sufficiently large n T (cid:39) . For all entries, n t = 49140 /n T . Accuracy of theapproximate calculation is δ = (cid:104) Q (cid:105) approximate / (cid:104) Q (cid:105) exact −
1, and is excellent, being betterthan 10 − where n T is small enough that the exact calculation can be performed andcompared to. (cid:104) Q (cid:105) , normalized by its standard deviation from zero for a comparable whitenoise process, is about 18, indicating the presence of a trend in the data.
7. Declaration of Interests
The authors have no competing interests to declare.
8. Funding
This work was supported by the National Science Foundation grant AGS-1639868. 18 ppendix A. Mathematical Identities for matrix m used in Soft-ware
From (12), it can be shown necessarily that m j,k = m k,j , a symmetricmatrix, and m j,nT − k +1 = − m j,k , a matrix odd under vertical or horizontalinversion. For the upper left matrix quadrant j ≤ n T − j and k ≤ n T − k itcan be shown from (12) that the following is necessary: m j,k = n T n t ( n T − (cid:18) n T − j (cid:88) m = (cid:98) nT (cid:99) +1 ψ ( k − n T ( n T − m ) n T − m ) − ψ (1 − k − mn T n T − m ) n T − m + 2 n T ( n T − k + 1) (cid:12)(cid:12)(cid:12) cos (cid:16) n T π (cid:17)(cid:12)(cid:12)(cid:12) (cid:19) (A.1)Here, ψ is the polygamma function of order 0 (e.g. MATLAB psi func-tion, Julia module SpecialFunctions’ polygamma function with zero as thefirst argument). From this, the following can be shown: m j +1 ,k = m j,k − n T n t ( n T − ψ ( k − n T ( n T − j ) n T − j ) − ψ (1 − k − jn T n T − j ) n T − j (A.2)This allows recursion down the columns of m in an exact calculation.Let the subtracted quantity from m j,k be f ( j, k ). The following is againnecessary: m j,j = m j +1 ,j +1 + f ( j, j ) + f ( j, j + 1) (A.3)This permits recurrence along the diagonal of m in an exact calculation.Finally, it also is necessary that m satisfy the following:19 j,k + m j − ,k − − m j,k − − m j − ,k = 2( d (1) j − ,k − + d (2) j − ,k − ) (A.4)Here d (1) j,k and d (2) j,k are given in equation (13) of the main text. This causes m as a whole to be calculable in O ( n T ) operations, in contrast to O ( n T ) fromequation (12) ( O ( n T ) per element). References [1] G. Ierley, A. Kostinski, Universal rank-order transform to extract signalsfrom noisy data, Phys. Rev. X 9 (2019) 031039. doi:10.1103/PhysRevX.9.031039 .[2] All calculations were performed with the first author’s personal Dell In-spiron 15 laptop, equipped with 2.2GHz AMD quadcore processor, using64 bit Windows 10. The best time scalings were provided by the slower32 bit version of MATLAB, and thus this is what is shown. For Tables 1and 2, 64 bit MATLAB was used, being the faster version.[3] L. T. Nguyen, J. Kim, B. Shim, Low-rank matrix completion: A contem-porary survey, IEEE Access 7 (2019) 94215–94237.[4] The Pierre Auger Collaboration, Pierre Auger Observatory Public Data,The Public Event Explorer, http://labdpr.cab.cnea.gov.ar/ED-en/index.phphttp://labdpr.cab.cnea.gov.ar/ED-en/index.php