 # A multidimensional analog to the Burrows-Wheeler transform

We show how to perform multidimensional pattern matching over an n-dimensional grid of text spanning a total of s characters with nength, an analog to the Burrows-Wheeler transform. Nength exploits a Fourier duality between two kinds of grid products to map a search problem that naively takes O(s^2) arithmetic operations to an equivalent problem that takes O(s s) arithmetic operations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Prefix

Suppose is an -dimensional rectangular grid of numbers. Let have size in dimension for . The total number of entries of is then given by . Denote the set of integers as and the set of positive integers as . Let mean and such that gives the same positive remainder as . Write entry of as or equivalently , and take for and . All grids discussed in this paper are special cases of .

Call some grid the pattern and another grid the text. Define a binary operator of multiplication over grids with the same dimensions where, for example, is given by

 [M]v0v1…vn−1=∑{wj}[P]w0w1…wn−1[T]v0−w0,v1−w1,…,vn−1−wn−1 for∑{wj}:=s0−1∑w0=0s1−1∑w1=0…sn−1−1∑wn−1=0 and v0,v1…,vn−1∈Z. (1)

Constrain so each of its nonzero entries is an integer from the alphabet , where is the size of . Suppose has nonzero entries . Constrain such that for , for some , and if and only if . This makes every nonzero entry of a distinct power of . Consider the query grid formed by replacing every given nonzero entry of with some for , and specify a query value as

 q=r−1∑w=0qw|Ω|rwqw∈Ω. (2)

Observe that if and only if every nonzero entry of takes the same value at the same position in a reversed rotation of given by

 [R]t0t1…tn−1=[T]v0−t0,v1−t1,…vn−1−tn−1t0,…,tn−1∈Z. (3)

So imagine writing the text by reversing a different grid via

 [T]t0t1…tn−1=[Tr]−t0,−t1,…,−tn−1t0,…,tn−1∈Z. (4)

Obtaining effectively searches for patterns matching the support of . Each entry of uniquely determines the of some query , so all entries of with the same value are at positions in that mirror the positions of all matches in to the same query. Call the search product operator, and call any grid formed by multiplying grids with the search product operator a search product. In the worst case, computing takes arithmetic operations: scalar products are summed to obtain each of entries. This paper shows how to compute with arithmetic operations.

## Ii Infix

The essential lesson of the Burrows-Wheeler transform (BWT) [burrows1994block] is that a redundant representation of a patterned object can collapse to an economical representation that simplifies pattern matching over that object [ferragina2000opportunistic, ferragina2005indexing]. Inspired by this lesson, rewrite redundantly as an -level circulant matrix . Think of as a circulant matrix whose every entry is itself a circulant matrix, whose every entry is itself a circulant matrix, and so on up to a depth of . Refer to entry of the block of the block of the block of the block matrix composed of a total of entries as . In analogy to how is indexed, take for , , and  . Then define as follows:

 [˜G](α0,β0),(α1,β1),…,(αn−1,βn−1)=[G]t0t1…tn−1 (5) such that βw−αw≡swtw for w∈{0,1,…,n−1} and α0,α1,…αn−1,β0,β1,…βn−1,t0,t1,…,tn−1∈Z.

Now consider a search product , where the are grids. Construct the corresponding -level circulant representations of these grids using the prescription of (5). The matrix product is computed as

 [˜S](α0,β0),(α1,β1),…,(αn−1,βn−1) =∑{γjk}([˜G0](α0,γ00),(α1,γ10),…,(αn−1,γn−1,0) ×[˜G1](γ00,γ01),(γ10,γ11),…,(γn−1,0,γn−1,1) ×…×[˜Gm−1](γ0,m−2,β0),(γ1,m−2,β1),…,(γn−1,m−2,βn−1)) =∑{γjk}([G0]γ00−α0,γ10−α1,…,γn−1,0−αn−1 ×[G1]γ01−γ00,γ11−γ10,…,γn−1,1−γn−1,0 ×…×[Gm−1]β0−γ0,m−2,β1−γ1,m−2,…,βn−1−γn−1,m−2) =∑{wjk}([G0]w00w10…wn−1,0[G1]w01−w00,w11−w10,…,wn−1,1−wn−1,0 ×…×[Gm−1]β0−α0−w0,m−2,β1−α1−w1,m−2,…,βn−1−αn−1−wn−1,m−2) =[M]β0−α0,β1−α1,…,βn−1−αn−1 for ∑{γjk}:=m−2∑ℓ=0⎛⎝s0−1∑γ0ℓ=0s1−1∑γ1ℓ=0…sn−1−1∑γn−1,ℓ=0⎞⎠, ∑{wjk}:=m−2∑ℓ=0⎛⎝s0−1∑w0ℓ=0s1−1∑w1ℓ=0…sn−1−1∑wn−1,ℓ=0⎞⎠, and α0,α1,…αn−1,β0,β1,…βn−1∈Z. (6)

Above, a change of variables is performed to substitute the dummy indices for the dummy indices . From (I) and (II), recovers the content of the search product , and in particular recovers the content of the search product .

A given -level circulant matrix representing the grid is diagonalized by an

-dimensional discrete Fourier transform matrix

defined via

 [F](α0,β0),(α1,β1),…,(αn−1,βn−1)=n−1∏w=0ραwβwsw√sw. (7)

Above, the th primitive root of unity is given by for and , and . The diagonal representation of is composed of at most nonzero entries, which may be arranged in a grid whose dimensions are the same as those of . Let correspond to the grids in the same way corresponds to . Because multiplication of diagonal matrices is equivalent to a Hadamard (entrywise) product, any search product of the is dual to a Hadamard product of the . Computing such a Hadamard product and transforming back to an -level circulant representation recovers the corresponding search product.

Call the duality between search products of the and Hadamard products of the nength duality. Call a nength of , and say that is obtained by nengthening .

Consider the Hadamard product , where and are the nengths of and , respectively. This product requires at most arithmetic operations, considerably fewer than the operations required by naive computation of the search product . Obtaining a nength takes arithmetic operations via a fast Fourier transform [cooley1965algorithm], and so does transforming back to a circulant representation. So nength maps a naively search operation to an operation.

## Iii Suffix

Reflect on the string case, where , to see how nength is an analog to the BWT. Recall that a string is encoded as integers in , and write all rotations of the string as the single-level circulant matrix . But instead of pursuing a lexical sort of the rows of this matrix and peeling off the BWT from the last column of the result, diagonalize with a discrete Fourier transform matrix and peel off the diagonal of the result to form a nength. A lexical sort makes sense only in one dimension, and its replacement with a discrete Fourier transform is apropos: the BWT’s sort maps periodic behavior in a string to runs of the same character in a transformed string, while a Fourier transform maps periodic behavior in a signal to peaks in Fourier space.

It has been 25 years since the BWT was published [burrows1994block], and this paper is our paean to it. We are excited to see where extensions of the BWT will be in another 25 years.

## Acknowledgements

We thank Aleksey Cherman, Ben Langmead, Craig T. Parsons, and Vanessa Vanzieleghem for their helpful feedback on early drafts of this manuscript. Abhi Nellore thanks Samia M. Hamed and Paul T. Spellman for helpful discussions. He also thanks Rachel A. Ward and the Oden Institute for Computational Engineering and Sciences at The University of Texas at Austin for hosting him as a visiting scholar August-October 2016, when he began work on extensions of the BWT.