Efficient decoding of interleaved subspace and Gabidulin codes beyond their unique decoding radius using Gröbner bases

An interpolation-based decoding scheme for L-interleaved subspace codes is presented. The scheme can be used as a (not necessarily polynomial-time) list decoder as well as a polynomial-time probabilistic unique decoder. Both interpretations allow to decode interleaved subspace codes beyond 
half the minimum subspace distance. Both schemes can decode γ insertions 
and δ deletions up to γ + Lδ ≤ L(nt − k), where nt is the dimension of the 
transmitted subspace and k is the number of data symbols from the field Fqm. 
Further, a complementary decoding approach is presented which corrects γ insertions and δ deletions up to Lγ +δ ≤ L(nt −k). Both schemes use properties 
of minimal Gr¨obner bases for the interpolation module that allow predicting 
the worst-case list size right after the interpolation step. An efficient procedure for constructing the required minimal Gr¨obner basis using the general 
K¨otter interpolation is presented. A computationally- and memory-efficient 
root-finding algorithm for the probabilistic unique decoder is proposed. The 
overall complexity of the decoding algorithm is at most O(L2n2 r) operations in 
F 
qm where nr is the dimension of the received subspace and L is the interleaving order. The analysis as well as the efficient algorithms can also be applied 
for accelerating the decoding of interleaved Gabidulin codes.


Introduction
Subspace codes have been proposed as a tool for error correction in noncoherent networks where the network topology and the in-network linear combinations performed by the intermediate nodes are not known by the transmitter and receiver [18,33]. Several code constructions, upper bounds on the size, and properties of such codes were thoroughly investigated in [2,[10][11][12]14,17,18,[33][34][35]41]. Subspace codes with efficient decoding algorithms include the Reed-Solomon (RS)-like code construction by Kötter and Kschischang [18] (KK) codes and the lifted maximum rank distance (MRD) codes (in particular based on Gabidulin codes, cf. [9,13,29]) by Silva, Kötter and Kschischang [33]. For decoding lifted MRD codes, the linear operations performed by the network must be reversed before starting the decoding process. This step is called reduction and allows to apply known error-erasure decoding schemes for MRD codes to lifted MRD codes, but with the additional computational cost of the reduction, see [33].
An interleaved MRD code consists of L matrices, which are codewords of L MRD codes (usually Gabidulin codes). One benefit of using lifted interleaved MRD codes in network coding is a reduced overhead for large interleaving orders. In [22,30,39] it was shown that probabilistic unique decoding as well as (not necessarily polynomialtime) list decoding of interleaved Gabidulin codes beyond half the minimum rank distance is possible. The main problem of list decoding subspace and MRD codes is that the list size might be exponential in the code length, see [28,38]. The approaches of Mahdavifar and Vardy [23,24] and Guruswami and Xing [15] provide list decoding schemes for subcodes and modifications of KK and MRD codes. Further, Trautmann, Silberstein and Rosenthal presented an approach for list decoding lifted Gabidulin codes [37]; here the complexity grows exponentially in the dimension of the subspace. Notice that the worst-case list size of many families of MRD codes grows exponentially with the code length, see [28], which therefore rules out polynomial-time list decoding for these code families.
In this paper, we present an interpolation-based decoding scheme for L-interleaved KK codes without the need for reduction at the receiver. This approach can be applied as a (not necessarily polynomial-time) list decoder or as a polynomialtime probabilistic unique decoder. The main contribution in this paper is that we show how a Gröbner basis for the interpolation module can increase the performance of the decoder while reducing the computational complexity. An efficient method for constructing minimal Gröbner bases is presented. Using an adapted version of the general linearized polynomial interpolation from [42], the computational complexity of the interpolation step is reduced from O(Ln 3 r ) to O(L 2 n r D) ≤ O(L 2 n 2 r ) operations in F q m , where n r is the dimension of the received space and D denotes the maximum degree of the constructed polynomial. The interleaving order L is a constant integer which we assume to be small compared to n r . Further, we propose a computationally-and memory-efficient root-finding algorithm for the unique decoder. Given a minimal Gröbner basis for the interpolation module, the algorithm reconstructs the message polynomials in O(L 2 k 2 ) operations in F q m , where k is the number of data symbols from F q m . This reduces the computational complexity for the root-finding step compared to O(L 3 k 2 ) for the recursive Gaussian elimination from [39].
Using the proposed algorithms, the total number of operations in F q m for the decoder is upper bounded by O(L 2 n 2 r ). The results are applied to solve the interpolation problem and root-finding system for interleaved Gabidulin codes efficiently [39].
This paper is structured as follows. In Section 2, we give basic definitions and describe notations. Section 3 explains the principle of our interpolation-based decoding algorithm, including calculating the maximum number of tolerated insertions and deletions, and clarifying how we generalize and improve principles from [18]. In Section 4, we outline how our ideas apply to list and unique decoding. Sections 5 and 6 provide efficient interpolation and root-finding algorithms, and Section 7 describes the entire decoding procedure. In Section 8 we compare the proposed decoding algorithm with known decoding schemes for interleaved subspace codes with respect to the error-correction performance and the computational complexity. Finally, Section 9 concludes this paper.
Parts of this work were presented at the 52nd Annual Allerton Conference on Communication, Control, and Computing, Oct. 2014, Monticello, IL, USA [7].

Preliminaries
2.1. Finite fields and subspaces. Let q be a power of a prime, and denote by F q the finite field of order q and by F q m its extension field of degree m. F N q denotes the vector space of dimension N over F q and P q (N ) is the set of all subspaces of F N q . The set of all -dimensional subspaces of F N q is called the Grassmannian and is denoted by G q (N, ).
Matrices and vectors are denoted by bold uppercase and lowercase letters, respectively, such as A and a. We index vectors and matrices beginning from zero, and [1, n] denotes the set of integers {1, 2, . . . , n}. The rank over F q and the kernel of a matrix A ∈ F m×n q are denoted by rk(A) and ker(A), respectively. For any fixed basis of F q m over F q , there is a bijective mapping between any vector a ∈ F n q m and a matrix A ∈ F m×n q . We often switch between these two representations. The row space of a matrix A ∈ F m×n q over F q is denoted by A q . For a column vector a ∈ F n q m we also use the notation a q to denote the F q -linear row space A q where A ∈ F m×n q is the representation of a over F q m . The sum U +V of two subspaces U, V in P q (N ) is the smallest subspace containing the union of U and V. The sum U + V of two subspaces U, V ∈ P q (N ) is a direct sum (or internal direct sum) U ⊕ V if U and V intersect trivially, i.e. if we have U ∩ V = {0} (see [16]). For a subspace V ∈ G q (N, ) the orthogonal subspace V ⊥ defined as The subspace distance between U, V in P q (N ) is where the last line follows since dim( A subspace code C s is a non-empty subset of P q (N ), and has minimum subspace distance d s when all subspaces in the code have distance at least d s from each other and there is at least one pair of subspaces with distance exactly d s . If every codeword is a subspace of dimension n t (t for "transmitted"), the code is called a constant-dimension code. In the following we only consider constant-dimension codes, i.e. subsets of G q (N, n t ). The code rate of a constant-dimension subspace code C s ⊆ G q (N, n t ) is defined as For a constant-dimension subspace code C s ⊆ G q (N, n t ) with minimum distance d s (C s ) the complementary code 1 is defined as the set of all orthogonal subspaces (see [18]) [18]) and code rate As channel model we use the operator channel from [18]. Such a channel has input and output alphabet P q (N ). The output U is related to the input V of dim(V) = n t by is not important for the performance of the code and can be chosen to be uniform (see [18]). The dimension of the received subspace U is thus n r = n t − δ + γ and we call δ the number of deletions and γ the number of insertions. For constantdimension codes, the subspace distance between the transmitted subspace V and the received subspace U is

2.2.
Linearized polynomials. Let a [i] def = a q i denote the i-th q-power of an element a ∈ F q m for any integer i. A nonzero polynomial of the form p(x) = d i=0 p i x [i] with p i ∈ F q m , p d = 0, is called a linearized polynomial of q-degree deg q (p(x)) = d, see [21,26]. For all a, b ∈ F q and x 1 , x 2 ∈ F q m , we have p(ax 1 + bx 2 ) = ap(x 1 ) + bp(x 2 ). Given two linearized polynomials p (1) (x) and p (2) (x) of q-degree d 1 and d 2 , their non-commutative composition p (1) ) is a linearized polynomial of q-degree d 1 +d 2 . The set of all linearized polynomials with coefficients from F q m forms a non-commutative ring L q m [x] with identity under addition + and composition ⊗. We denote the set of all linearized polynomials from L q m[x] with q-degree less than k by L q m[x] <k . For a column vector a = (a 0 a 1 . . . a n−1 ) T ∈ F n q m and p(x) ∈ L q m[x] we define p(a) def = (p(a 0 ) p(a 1 ) . . . p(a n−1 )) T . The set of all multivariate linearized polynomials of the form For brevity, we also denote the (L + 1)-variate linearized polynomial Q(x, y 1 , . . . , y L ) by Q. Further, for a vector w = (w 0 , w 1 , . . . , w L ), we denote the w-weighted q-degree deg w of a multivariate linearized polynomial Q ∈ L q m[x, y 1 , . . . , y L ] as deg w (Q) def = max w 0 + deg q (Q 0 (x)), w 1 + deg q (Q 1 (y 1 )), . . . , w L + deg q (Q L (y L )) .
The q-degree is the logarithm to the base q of the usual degree and thus the weights in deg w are additive instead of multiplicative as in the classical case.
Definition 2.1 (w-Weighted Linearized Monomial Ordering). Given a vector w = (w 0 w 1 . . . w L ) the w-weighted total order ≺ w on linearized monomials is defined for all j, j ∈ [0, L] and some ≥ 0 as The w-weighted linearized monomial ordering is also called w-weighted term over position ordering [1] since first the w-weighted degree of the term is considered and the position j is considered only if two monomials have the same w-weighted degree. To obtain the linearized monomial order ≺ w for polynomials from L q m[x, y 1 , . . . , y L ] we substitute y 0 = x.
For w = (0 k − 1 . . . k − 1) the total order ≺ w on linearized monomials is . The set of all polynomials from L q m[x, y 1 , . . . , y L ] with weighted degree less than D is denoted by L q m[x, y 1 , . . . , y L ] <D .
We identify the leading term LT(Q) of any multivariate polynomial Q as the maximum monomial under ≺ w normalized by its coefficient. The index of a univariate polynomial is defined as Ind(l(x) ⊗ y j ) = j for l(x) ∈ L q m[x]. We denote the leading position of a multivariate polynomial Q ∈ L q m[x, y 1 , . . . , y L ] by LP(Q) = Ind(LT(Q)). For a set Q ⊆ L q m[x, y 1 , . . . , y L ] we denote the set of all leading positions of the elements in Q by LP(Q) The Moore matrix of the vector a = (a 0 a 1 . . . a n−1 ) ∈ F n q m is defined as The rank of M r (a) is min{r, n} if the elements a 0 , . . . a n−1 are linearly independent over F q , see [21].

2.3.
Lifted rank-metric codes. The minimum rank distance d r of a linear code C ⊆ F n q m is defined as where X, Y are the matrix representations of x, y ∈ C. The Singleton-like bound on the minimum rank distance states that d r ≤ n − k + 1 when m ≥ n, see [9,13,29]. Codes which attain this bound are called MRD codes. A special class of MRD codes are Gabidulin codes [9,13,29], which are the analogs of Reed-Solomon codes in rank metric. Interleaved Gabidulin codes are a horizontal or vertical concatenation of L Gabidulin codes, for details see [22,30,39]. If the code rate of all component codes is the same, we call the code a homogeneous L-interleaved Gabidulin code.
Lifting an (interleaved) MRD code means that we append an identity matrix to each transposed code matrix and consider the row space of these composed matrices as a subspace code. The subspace distance of this constant-dimension code is twice the rank distance of the (interleaved) MRD code [33]. The linear combinations of the received packets at the intermediate nodes in a network coding scenario are therefore tracked by the leading identity matrix and can be inverted by computing the reduced row echelon form of the received matrix, called reduction.

2.4.
Interleaved Kötter-Kschischang subspace codes. Motivated by lifted interleaved Gabidulin codes [20,31,39] and the construction of Kötter and Kschischang [18], we define interleaved KK subspace codes as follows. T with n t ≤ m be a vector containing F q -linearly independent code locators from F q m . For fixed integers k (1) , . . . , k (L) ≤ n t , an interleaved KK subspace code ISub[L, a; n t , k (1) , . . . , k (L) ] of dimension n t and interleaving order L is defined as If k (j) = k, ∀j ∈ [1, L] we call the code a homogeneous interleaved KK subspace code and denote it by ISub[L, a; n t , k]. The dimension of the ambient vector space is N = n t + Lm and the code rate of this construction is For increasing interleaving order L, the rate loss caused by the appended code locators decreases since n t Lm. For L = 1 the codes from Definition 2.2 are equivalent to KK codes [18]. The subspace codes over F q Lm with code locators from F q m in [15] are homogeneous L-interleaved KK subspace codes over F q m (see [6]). Proof. Let V and V be two codewords generated by f (1) (x), . . . , f (L) (x) and g (1) (x), . . . , g (L) (x) with q-degrees less than k (1) , . . . , k (L) . Since dim(V) = dim(V ) = n t we have Thus the minimum distance is achieved for the maximum dimension of the intersection space V ∩ V . The dimension of V ∩ V is minimal when the evaluation polynomials of maximum q-degree, say f (j) (x) and g (j) (x), are distinct and all other evaluation polynomials are identical. Suppose dim(V ∩ V ) = r, i.e., f (j) (x) and g (j) (x) agree on r linearly independent points. Since two distinct linearized polynomials f (j) (x) and g (j) (x) of q-degree at most k (j) − 1 can agree on at most k (j) − 1 linearly independent points (see [18,Lemma 11]), it follows that r ≤ max{k (j) } − 1.
Hence we have In order to show that equality holds, let f (j) (x) and g (j) (x) be distinct but they agree on exactly k (j) − 1 points, where j is chosen such that k (j) = max i k (i) .
All other evaluation polynomials are chosen to be equal. Such a choice of f (j) (x) and g (j) (x) clearly exists since each polynomial consists of k (j) coefficients. In this case, we get that the subspace distance between the two codewords generated by f (1) (x), . . . , f (L) (x) and g (1)

Interpolation-based decoding
In [39] an interpolation-based decoding scheme for interleaved Gabidulin codes that corrects rank errors and row/column erasures was presented. The problem of decoding constant-dimension subspace codes from insertions and deletions can be mapped to a generalized rank-metric decoding problem [32,33] by computing a canonical form of the basis for the received subspace, called reduction. The interpolation-based scheme proposed in [39] considers lifted interleaved Gabidulin codes and requires the reduction of the received basis (Gaussian elimination) prior to the decoding procedure.
In the following, we propose an interpolation-based decoding approach for interleaved KK subspace codes that generalizes the Kötter-Kschischang approach [18] and is related to interpolation-based decoding of interleaved Gabidulin codes [39]. The decoder takes any basis of the received subspace without the need for the reduction which reduces the complexity at the receiver side. The decoding scheme has an improved decoding performance for insertions compared to [39]. The decoding principle consists of an interpolation step and a root-finding step which are described and analyzed in the following. For the interpolation step, we must solve the following interpolation problem.
Problem 1 (Interpolation Problem). Given the integers n r , D and k (1) , . . . , k (L) , find a nonzero (L + 1)-variate linearized polynomial from L q m[x, y 1 , . . . , y L ] of the form (12) Q (x, y 1 , . . . , y L ) = Q 0 (x) + Q 1 (y 1 ) + · · · + Q L (y L ), that satisfies the following conditions: Denote the coefficients of (12) by In the following letk = L j=1 k (j) . We can find the coefficients of Q (x, y 1 , . . . , y L ) by solving a linear system of equations where R is an n r × (D(L + 1) −k + L) matrix (see (8) for the definition of M r (a)): and Lemma 3.1. If the degree constraint is chosen as (12) fulfilling the interpolation constraints in Problem 1 exists.
Proof. The number of linearly independent equations is at most n r which must be less than the number of unknowns in order to guarantee that there is a nonzero solution. Hence we have (16) and the choice of D in (15) ensures the existence of a nonzero solution to Problem 1.
The degree constraint D can be computed at the receiver, since n r ,k and L are known at the receiver. A polynomial Q ∈ L q m[x, y 1 , . . . , y L ] fulfills the degree constraints in Problem Given a solution Q (x, y 1 , . . . , y L ) of Problem 1 we define a univariate polynomial Lemma 3.2 (Roots of Interpolation Polynomial). The linearized polynomial P (x) in (17) has at least n t − δ linearly independent roots in F q m .
where ξ 0 , . . . , ξ nt−δ−1 are elements from F q m that are linearly independent over F q . The linearized polynomial Q (x, y 1 , . . . , y L ) vanishes on all elements in the received subspace U. Since U ∩ V is a subspace of the transmitted space V we have The following theorem shows that the message polynomials f (1) (x), . . . , f (L) (x) are roots of Q (x, y 1 , . . . , y L ) under certain constraints.
Proof. By Lemma 3.2, the linearized polynomial P (x) has at least n t − δ linearly independent roots in F q m . If we choose (16) and (20) we have We say that the received subspace is decodable if (18) holds and thus call (18) the decoding region. In other words, the decoding region defines the maximum number of γ insertions and δ deletions where the decoder is able to decode. From (18) we see that interleaving makes the code more resilient against insertions while the performance for decoding deletions remains the same as in [18].
Intuitively, it makes sense that a subspace of dimension n t with n t < N/2 can tolerate less deletions than insertions. Since its dimension is n t , we can have at most n t deletions whereas we have at most N − n t > N/2 insertions. For L-interleaved KK subspace codes with n t ≈ m, we have N = (L + 1)n t and thus can correct L times more insertions than deletions.
A similar behavior can be observed for the subspace code construction in [15,24] and the folded subspace codes in [5]. Also, a similar degree of asymmetry can be seen between correctable insertions and deletions in the average list size of interleaved KK subspace codes (see (32)).
For a homogeneous interleaved code C s = ISub[L, a; n t , k] the decoding region is Figure 1 shows the decoding region of a homogeneous 4-interleaved subspace code compared to a Kötter-Kschischang code (L = 1) for d s = 8 (i.e. n t − k = 7). The figure illustrates the increased resilience against insertions due to interleaving.
we can express the normalized decoding radius τ = γ+Lδ nt as (22) τ For Lm n t the decoding radius approaches τ ≈ L (1 − R). For L = 1 the decoding region (18) is the same as in [18,Eq. 11] and improves upon [18] for higher L.  The interpolation-based decoder in [39] uses a basis of the F q m -linear solution space of (13) to solve the root-finding step, i.e., to recover the coefficients of the polynomials f (1) (19).
In our approach we use a so called minimal Gröbner basis for the interpolation module to solve the interpolation system (13) and the root-finding system (19) efficiently.

Gröbner bases
As in [42], V can be partitioned as V = j S j , where S j = {Q ∈ V : LP(Q) = j}. Hence, the subset S j contains all Q ∈ V with leading position j.
Define a set of n r functionals F i : The kernel of F i is denoted by K i and we define Notice, that all elements in K nr satisfy the first interpolation condition in Problem 1. In particular, all elements from K nr of w-weighted degree less than D are a solution for Problem 1.
Since F i are linear functionals F i (p+q) = F i (p)+F i (q) = 0 holds for any p, q ∈ K nr and i ∈ [1, n r ]. Thus K nr is closed under addition and the left symbolic product with any element from L q m[x].
In the following we call K nr the interpolation module. Properties of (left) L q m[x]modules and submodules have been studied in [19]. An important observation is, that K nr is a free L q m[x]-submodule that can be described by a basis as linearized polynomials are a special case of skew polynomials (see [27,Theorem 14]). In this paper we consider a special class of bases, called Gröbner bases, for the free L q m[x]submodule K nr . Definition 3.5 (Gröbner Basis [8,19]  The proposition follows directly from Definition 3.6 since two polynomials g i , g j ∈ G that have the same leading position would contradict that either Proof. By Lemma 3.1 there exists a polynomial satisfying the constraints in Problem 1 if D fulfills (15). Since G is a minimal Gröbner basis, the minimal (w.r.t. the weighted degree) polynomial among all polynomials satisfying the interpolation constraints is in G <D and we have |G <D | ≥ 1. Since G is an ordered minimal Gröbner basis the leading positions of g 0 , g 1 , . . . , g L are distinct [19,Proposition 7] and hence |G | ≤ L+1. We now show, that g 0 is never a solution to the interpolation problem. Let g 0 = Q 0 (y 0 ) + Q 1 (y 1 ) + · · · + Q L (y L ) be nonzero with deg w (g 0 ) < D.
Since LP(g 0 ) = 0 (i.e. leading term in y 0 ) we have (24). Thus we must have deg w (g 0 ) ≥ D and g 0 is never contained in G <D . Hence, |G <D | ≤ L.
The proof of Lemma 3.7 shows that a polynomial with leading position 0 (i.e. leading term in y 0 ) is never a solution to the interpolation problem and thus never contained in G <D .
Lemma 3.8. Let G be a minimal Gröbner basis for K nr under ≺ w . Then all solutions Q ∈ L q m[x, y 1 , . . . , y L ] <D of Problem 1 are of the form Proof. Let the set J = LP(G <D ) contain the leading positions of the elements in G <D and let J = LP(G \ G <D ). Since G <D is a subset of the minimal Gröbner basis G the leading positions of the elements in G <D are all distinct implying that the sets J and J are disjoint. We now show that there exists no element in the module generated by G \ G <D that has w-weighted degree less than D. Suppose there exists aq in the module generated by G \ G <D with deg w (q) < D. If LP(q) ∈ J then we must have thatq is an element in the module generated by G <D . If LP(q) ∈ J then we have thatq is in the module generated by is not contained in the module generated by LT(G ). This contradicts the assumption that G is a minimal Gröbner basis for K n since by definition we must have that LT(G ) and LT(K n ) span the same module. Hence all elements in L q m[x, y 1 , . . . , y L ] <D must be L q m[x]-linear combinations of the elements in G <D .
Thus every polynomial in the F q m -linear solution space of (13) can be constructed by an L q m[x]-linear combination of the polynomials in G <D which implies that d I ≥ |G <D |.
Clearly, the set K nr ∩ L q m[x, y 1 , . . . , y L ] <D is not a submodule of V , since an L q m[x]-linear combination of the form may result in a polynomial of w-weighted degree larger than D − 1.
3.3. Root-finding step using minimal Gröbner bases. The task of the rootfinding step is to find all polynomials f (j) (x) with deg q (f (j) (x)) < k (j) , ∀j ∈ [1, L], such that (19) holds. Similar to [39], this is done by solving a linear system of equations in F q m for the coefficients of f (1) Instead of using a basis for the F q m -linear solution space of (13) we use a subset of a minimal Gröbner basis of the L q m[x]-submodule K nr to solve the root-finding step. Let G be an ordered minimal Gröbner basis for K nr and define the set G <D ⊂ G as in (23). List decoding of rank-metric codes using Gröbner bases was also considered in [36]; however only for the non-interleaved case (L = 1).
We now set up the root-finding matrix using r = |G <D | polynomials Q ( ) (x, y 1 , . . . , y L ) ∈ G <D , ∈ [1, r]. Denote the coefficients of the polynomials Q ( ) ∈ G <D by We assume that f The root-finding matrix can be set up as and the roots can be found by solving the system of equations (27) Q Solving the root-finding system (27) recursively requires at most O(L 3 k 2 ) operations in F q m (see [39]). In Section 6 we present an efficient root-finding algorithm that solves (27) in O(L 2 k 2 ) operations in F q m .
The following lemma lower bounds the rank of the root-finding matrix (26) by using the properties of the minimal Gröbner basis polynomials.
for all j ∈ J and thus conclude that rk(Q) ≥ rk(Q) ≥ j∈J k (j) where the first inequality follows from the fact that we only considered a submatrix of Q.
Equipped with Lemma 3.9 we can upper bound the dimension of the affine solution space of (27).  has at most q m( j∈J k (j) ) solutions.
Proof. The dimension of the affine solution space of (27) over F q m is L j=1 k (j) − rk(Q). Let the set J = LP(G <D ) contain the leading positions of the polynomials in G <D . By Lemma 3.9 we have rank Q ≥ j∈J k (j) and thus we can have at most For the homogeneous case k (j) = k, ∀j ∈ [1, L] we have at most q m(kL−k·|G <D |) = q mk(L−|G <D |) solutions.

Application to list and unique decoding
4.1. List decoding approach. In general, the root-finding matrix Q (26) does not always have full rank. In this case, we obtain a list of roots of (19), i.e., a list of possible (interleaved) message polynomials. This decoder is not a polynomial-time list decoder but it provides with quadratic complexity the basis of the list.
An important benefit of using a Gröbner basis for the root-finding step is that Theorem 3.10 shows that the actual list size can be upper bounded if G <D is known. Hence, the actual list size can be upper bounded right after interpolating a minimal Gröbner basis for K nr without checking the rank of the root-finding matrix (26). This feature is useful if the list decoder is operated such that it declares a decoding failure if the actual list size is too large (i.e. exceeds a defined threshold).
The worst-case list size for the decoder is upper bounded as follows.

Lemma 4.1 (Maximum List Size).
Let Q be as defined in (26). Then the number of solutions of the root-finding system (27) is at most Proof. From Lemma 3.7 we know, that |G <D | ≥ 1. By Theorem 3.10 the rootfinding system (27) has at most q m( L j=1 k (j) − j∈J k (j) ) solutions which is maximized for J = arg min j∈ [1,L] For the homogeneous case k (j) = k, ∀j ∈ [1, L] this gives q m(Lk−k) = q mk(L−1) . Although the root-finding system (27) can have up to q m( L j=1 k (j) −min j∈[1,L] k (j) ) solutions, not all solutions necessarily give codewords in subspace distance up to the decoding region. Thus the proposed decoding scheme is not necessarily a polynomial time decoder since in the worst case the list size is exponential. However, this decoding scheme computes in polynomial time a basis for the list of all codewords in the decoding region (18).

Advances in Mathematics of Communications Volume 12, No. 4 (2018), 773-804
Note that if we choose n t ≈ m in (30), we observe an asymmetry between insertions and deletions of degree L, i.e., deletions affect the average list size of the code L times more than insertions.
Consider a homogeneous interleaved KK subspace code C s = ISub[L, a; n t , k]. To decode interleaved subspace codes with a probabilistic unique decoder we require the average list size to be close to one. This is fulfilled if the exponent in (30) becomes negative, i.e., if we have For n t ≈ m we get Note, that the decoding region of the list decoder (18) shows a similar order of asymmetry between insertions and deletions as (32), i.e. approximately L times more insertions γ than deletions δ can be tolerated. From (32) we also see that a good list decoder for interleaved KK subspace codes should return on average a list of size close to one if γ and δ satisfy (32). This observation motivates to consider the list decoder as a probabilistic unique decoder.

4.2.
Probabilistic unique decoder. This decoding principle can also be used as a probabilistic unique decoder. We obtain a unique solution to the root-finding problem (27) if the rank of Q is full. If Q does not have full rank, we declare a decoding failure and call the occurred error non-correctable. In the following, we show that under certain conditions the probability of a non-correctable error is very small. Proof. If we set up Q using G <D then by Lemma 3.9 the root-finding matrix Q has full rank if |G <D | = L. By Lemma 3.7 we always have 1 ≤ |G <D | ≤ L. Let d I be the dimension of the solution space of (13). Now assume that |G <D | < L and suppose that there exist d I ≥ L polynomialsQ (1) , . . . ,Q (d I ) ∈ K nr with deg w (Q (j) ) < D for all j ∈ [1, d I ] such that Since the root-finding system set up with G <D cannot have a unique solution if |G <D | < L there must exist a polynomialQ (j) for some j ∈ [1, d I ] that is not of the form (34) which contradicts that G <D is a minimal Gröbner basis for K nr ∩ L q m[x, y 1 , . . . , y L ] <D (see Lemma 3.8). Lemma 4.3 shows that if the root-finding system has a unique solution, a minimal Gröbner basis for K nr suffices to find this unique solution. Lemma 4.3 also implies that we may get a unique solution only if for the matrix R from (13) it holds that We now derive a degree constraint D u ≥ D which ensures that d I ≥ L. Then the dimension d I of the solution space of (13) satisfies d I ≥ L.
Proof. The dimension d I of the solution space of (13) is where rk(R) ≤ n r . In order to satisfy (35) we require Hence, the choice of D u in (36) ensures that d I ≥ L.
Lemma 4.4 provides the degree constraint D u for the unique decoder to obtain the set G <D for D = D u .
Since Lemma 4.3 holds in both directions we know that at least one of the L polynomials g 1 , . . . , g L ∈ G violates the degree constraint D u if there is no unique solution, i.e. |G <D | < L for D = D u . Thus, we declare a decoding failure if at least one polynomial of g 1 , . . . , g L ∈ G has w-weighted degree at least D u .
Compared to the matrix-based approach in [39] the detection of decoding failures is more efficient since computing the rank of a d I × L matrix needs more operations than determining the degree of a linearized polynomial.
To recover the message polynomials f (1) (x), . . . , f (L) (x), (19) has to be satisfied for all sets of L polynomials that are a solution of Problem 1 with D = D u . This is true if D u is less than or equal to the number of linearly independent interpolation points, i.e. if  Notice that nr+k L+1 ≥ nr+k L+1 and we therefore use the lower bound in (37) as it gives the largest decoding region.
For a homogeneous interleaved KK subspace code C s = ISub[L, a; n t , k] the probabilistic unique decoding region (37) is Compared to the decoding region of the list decoder (18), the ability for probabilistic unique decoding (37) comes with a penalty of (L − 1) on the decoding region.
The decoding region (38) shows a similar asymmetry between insertions and deletions as the decoding region that we derived from the average list size of homogeneous interleaved KK subspace codes (32). Figure 2 shows the improvement of the decoding region for the probabilistic unique decoder (38) upon Kötter-Kschischang codes (L = 1) due to interleaving. Compared to the decoding region of the list decoder in Figure 1 we see that the number of correctable insertions of the probabilistic unique decoder is decreased by (L − 1).
We now upper bound the decoding failure probability using results from [39].
Proof. Assume w.l.o.g. that the first n t − δ dimensions (x i , r , correspond to the basis of the non-corrupted subspace U ∩ V. Since the elements x 0 , . . . , x nt−δ−1 are linearly independent over F q and D u ≤ n t − δ, the rank of the (n t − δ) × D u q-Vandermonde matrix R = M Du ((x 0 , . . . , x nt−δ−1 )) T equals rk R ) = D u . The rank of R depends on rk(R ) and dim(E) = γ and we have rk(R) ≤ D u + γ. The dimension of the solution space of the interpolation system thus satisfies If the number of deletions is unknown at the receiver, we have to choose D u according to (36) in order to satisfy (20). In this case d I ≥ L nr+k L+1 + 1 −k − γ (see (40)). Substituting this in (39) gives the upper bound on the failure probability In a scenario where the number of deletions δ is known at the receiver we can choose D u = n t − δ. The dimension of the solution space of (27) then satisfies d I ≥ L(n t − δ + 1) −k − γ which is larger than in the case where δ is unknown at the receiver. For the unique decoder the maximum number of correctable insertions γ is γ max = L(n t − δ) −k. For the failure probability we get Notice that the failure probability in case δ is known at the receiver is less than or equal to the case if δ is unknown.
The following theorem summarizes the properties of the unique decoder. 1.6 · 10 6 0 0 Table 1. Simulation results of the probabilistic-unique decoder with parameters m = n t = 6, k = 2 and L = 2.

4.3.
Simulation results. In a simulation with 10 6 transmissions over an operator channel with δ = 0, γ = 5 and code parameters q = 2, m = 8, n t = 7, k (1) = k (2) = 4, L = 2, we observed a fraction of 1.5 · 10 −5 non-correctable errors (upper bound 6.1 · 10 −5 ). We simulated a code with parameters m = n t = 6, k = 2 and L = 2 over an operator channel with γ = 2, . . . , 8 and δ = 0, 1. The results are shown in Table 1 and illustrated in Figure 3.  At the receiver we calculate the dual space U of the received space U ⊥ , which is shown to be a codeword from C s which is corrupted by δ insertions and γ deletions, i.e., an insertion into U ⊥ becomes a deletion in U and vice versa. We can then perform decoding with the known interleaved decoding algorithms in the code C s .
Lemma 4.8. We have Thus, H N −nt−δ (V ⊥ ) ⊥ can be decomposed as a direct sum of V and some space of dimension δ, which is denoted by E δ .
Then, we obtain the following statement for the error-correction capability of C ⊥ s . Theorem 4.9. Given the output of the operator channel then we can reconstruct V ⊥ with probability at least Proof. We obtain: where the third line follows by Lemma 4.8 and the fourth step follows since E δ ∩E ⊥ γ = E δ . From the operator channel we have Thus the received space U corresponds to a codeword V ∈ C s that is corrupted by at most γ deletions and δ insertions. Due to (38), we can reconstruct V if Lγ + δ ≤ L(d s (C s ) − 2)/2 with high probability. Once we have V, we can simply calculate the dual space and get V ⊥ .
The asymmetry between correcting deletions and insertions therefore flips compared to (38). The code C s is more resilient against insertions while the complementary code C ⊥ s is more resilient against deletions. Notice that the list decoder can be adapted to the complementary code in a similar way.

Efficient interpolation of minimal Gröbner basis polynomials
In this section, we show how the generalization [42] of the multivariate Kötter interpolation [40] and the row reduction algorithm from [27] can be used to construct an ordered minimal Gröbner basis for the left L q m[x]-submodule K nr efficiently.

General linearized Kötter interpolation.
A single polynomial which is a solution to Problem 1 can be constructed efficiently by the general linearized polynomial interpolation algorithm by Xie, Yan and Suter [42], requiring at most O(L 2 n r D) operations in F q m . Our main contribution is to show how this algorithm can be used to construct a minimal Gröbner basis for K nr under ≺ w efficiently. Algorithm 1 in [42] iteratively constructs L + 1 polynomials g i,j in each step i that are minimal in T i,j w.r.t. the w-weighted degree for j ∈ [0, L]. The output of the algorithm in [42] is one polynomial Q * of least weighted degree among all polynomials in K nr .
Hence we modify the algorithm so that it outputs g nr,0 , . . . , g nr,L and denote this adapted version by InterpolateBasis(x T , r (1)T , . . . , r (L)T ). The pseudo-code of InterpolateBasis(·) is given in Algorithm 1. Proof. By [42, Lemma 2] each g nr,j = min{K nr ∩ S j } = min{T nr,j } with respect to the weighted degree deg w implying that LP(g nr,j ) = j. Since K nr = T nr,0 ∪ T nr,1 ∪ · · · ∪ T nr,L , the polynomials g nr,0 , . . . , g nr,L are a basis for the L q m[x]-submodule K nr . Suppose there exists a p ∈ G such that LT(p) divides LT(g nr,j ) for q ∈ G , p = g nr,j . Then g nr,j can be reduced modulo p resulting in a polynomial with leading position j and w-weighted degree less than deg w (g nr,j ). This contradicts the minimality of g nr,j in T nr,j .
Since we only modify the output of [42, Algorithm 1] the construction of the minimal Gröbner basis does not require more operations than constructing one polynomial. Thus the number of required operations in F q m is on the order of O(L 2 n r D) (see [42]). Since D ≤ n r we have O(L 2 n r D) ≤ O(L 2 n 2 r ). This interpretation of [42] can be directly applied to the decoding approach for interleaved Gabidulin codes in [39]. Fig. 4 shows the number of multiplications needed for efficient interpolation and compares this to the number of multiplications needed to solve (14) by Gaussian elimination. The figure shows that our interpolation-based decoding substantially reduces complexity for large n t .

Efficient root-finding
In [39] it was shown that the root-finding system (27) can be solved recursively with at most O(L 3 k 2 ) operations in F q m . In this section, we present an efficient rootfinding algorithm for the probabilistic unique decoder that finds all roots f (j) (x), j ∈ Algorithm 1: Adapted Linearized Polynomial Interpolation by [42] InterpolateBasis(x T , r (1)T , . . . , r (L)T ) Input : A basis for the received subspace (x i , r Output: An ordered minimal Gröbner basis G = {g nr,0 , . . . , g nr,L } for K nr under ≺ w 1 Initialize: g 0,0 = x, g 0,j = y j , j ∈ [1, L]  [1, L] using an ordered minimal Gröbner basis G <D with |G <D | = L. The complexity of this algorithm is O(L 2 k 2 ) operations in F q m . 6.1. An efficient root-finding algorithm. Let G = {g 0 , g 1 , . . . , g L } be an ordered minimal Gröbner basis under ≺ w for the left L q m[x]-submodule K nr . Suppose that all polynomials in G except g 0 (see Lemma 3.7) have weighted degree less than D u , i.e. |G <D | = L. Since G is a minimal Gröbner basis the polynomials in G <D have distinct leading positions. This property imposes a triangular structure on the linear root-finding system (27) and allows to solve it efficiently. In particular, we can solve the upper-triangular subsystem (28) for the case where |G <D | = L since we know from Lemma 3.9 that this system has full rank. Thus, using an ordered minimal Gröbner basis allows to reduce the complexity of the root-finding step.
If |G <D | = L, Algorithm 2 recursively determines all unique roots f (j) (x), j ∈ [1, L]. Algorithm 2 is based on the symbolic right division of univariate linearized polynomials and solves the root-finding system recursively.
Algorithm 2: Unique root-finding of y 1 , . . . , y L -minimal polynomials 0 (x)) and e = deg q (Q (j) j (y j )). The polynomial Q (j) (x, y 1 , . . . , y L ) is y j -minimal for any j ∈ [1, L] and thus the q-degrees fulfill b ≤ e + k (j) − 1 and deg q (Q (j) j (y j )) < e for j > j. This implies that the coefficients are q (j) j ,e = 0 for j > j. In the first iteration j = i = 1 we must solve ,e = 0 for ∈ [2, L] the calculation reduces to can uniquely be determined as which corresponds to line 5 in Algorithm 2 and it is easy to check that 1,e · (−q which reduces the q-degree of Q (1) 0 (x) by one since we enforced (44). In the next iteration i = 1, j = 2 we must solve due to the y 2 -minimality of Q (2) (x, y 1 , . . . , y s ). Hence, we can directly calculate the coefficient f (2) k (2) −1 using the updated polynomial Q (2) 0 (x) from the previous step. The algorithm computes the unique monomial f which is possible since the right hand side is computed in step i, j −1 of Algorithm 2.
6.2. Complexity analysis. Suppose we use a normal basis (which always exists). The calculation of q-powers then corresponds to a cyclic shift over the base field and its complexity can be neglected. Hence, multiplications dominate the complexity. Proof. Each step of Algorithm 2 provides L compositions of a (univariate) linearized polynomial of q-degree at most D−k with one monomial. If we consider the inversion of an element from F q m as a multiplication we have L(D − k + 2) multiplications per step. In total we need kL 2 (D − k + 2) ≤ k 2 L 2 multiplications in F q m .
A comparison of this complexity to the recursive Gaussian elimination (GE) from [39] is illustrated in Figure 5. Besides the reduction in computational cost, Algorithm 2 requires less memory than the recursive GE. An upper bound on the memory requirement is given by the following lemma. Proof. The algorithm must store L (L + 1)-variate polynomials with each at most D(L + 1) − L(k − 1) coefficients in F q m , as well as the Lk coefficients in F q m of the L message polynomials. Hence, in total we must store L 2 (D − k + 1) + L(D + k) elements in F q m .
The memory requirements of Algorithm 2 and the recursive GE are illustrated in Figure 6.

Efficient decoding algorithm
We now summarize the efficient decoding procedure for list and probabilistic unique decoding of interleaved KK subspace and Gabidulin codes. containing all tuples of message polynomials f (1) (x), . . . , f (L) (x) that correspond to codewords in subspace distance at most γ + δ from the received space can be obtained requiring at most O(L 2 n 2 r ) operations in F q m . Proof. The decoding region follows from Theorem 3.3 and the maximum list size follows from Lemma 4.1. The required minimal Gröbner basis for the interpolation module K nr can be constructed using Algorithm 1 in at most O(L 2 n 2 r ) operations in F q m . Then the unique tuple f (1) (x), . . . , f (L) (x) containing the linearized message polynomials f (j) (x) ∈ L q m[x] <k (j) for all j ∈ [1, L] that satisfy (19) can be found with probability at least The decoding procedure uses the efficient interpolation and root-finding algorithm for list decoding of KK codes and is summarized in Algorithm 3. By Theorem 3.10 the list size can be upper bounded from the set of leading positions J = LP(G \ G <D ). Thus Algorithm 3 can be modified so that it returns a list of bounded size by computing the list size using J right after the construction of G <D or declare a decoding failure if the list size exceeds a fixed threshold. In order to . . , f (L) (x) that correspond to solutions of the root-finding system (27) efficiently decode interleaved Gabidulin codes of length n, elementary dimensions k (j) and interleaving order L as defined in [22,30,39] we set n t = n r = n. Let g = {g 0 , . . . , g n−1 } ⊂ F q m with n ≤ m denote the linearly independent code locators of the interleaved Gabidulin code and denote by y (j) , j ∈ [1, L] the elementary received words. Then, Algorithm 3 called with (g T , y (1)T , . . . , y (L)T ) can decode errors of rank t up to t < L(n−k+1)

L+1
(see [39]). The complete procedure for the probabilistic unique decoder for interleaved KK codes is given in Algorithm 4. To decode an interleaved Gabidulin code, this procedure must be called with (g T , y (1)T , . . . , y (L)T ). As in [39], the decoder finds a unique solution with high probability for t ≤ L(n−k) L+1 . The computational complexity of Algorithm 3 and Algorithm 4 is dominated by the interpolation and the root-finding step. Using Algorithm 1 the complexity of constructing the minimal Gröbner basis for the interpolation module is in the order of O(L 2 n 2 r ) operations in F q m . The root-finding step using Algorithm 2 requires at most O(L 2 k 2 ) operations in F q m . Thus the overall computational complexity of the list and unique decoding algorithm is at most O(L 2 n 2 r ) operations in F q m where k = max j∈ [1,L] {k (j) }. Similar to previous work [18], we count operations in F q m , i.e., this includes multiplications, additions, and inversions.

Comparison to known decoding schemes
It has been shown that decoding KK subspace codes can be transformed to a generalized rank-metric decoding problem [33]. The main idea of this approach is to use information from the channel to decode a codeword of an interleaved Gabidulin Output: "decoding failure"

Decoding scheme
Decoding region Op. in F q m Li-Sidorenko-Silva [20,31] (L + 1)t + Lκ + Lρ ≤ L(n t − k) O(Ln 2 t ) Wachter-Zeh-Zeh [39] (L + 1)t + Lκ + Lρ ≤ L(n t − k) O(L 3 n 2 t ) Guruswami-Xing [15] (L + 1)t + κ + Lρ ≤ L(n t − k) O(L 6 n 2 t ) Bartz-Meier-Sidorenko [4] (L + 1)t + κ + Lρ ≤ L(n t − k) O(L 3 n 3 t ) This contribution (L + 1)t + κ + Lρ ≤ L(n t − k) O(L 4 n 2 t ) Table 2. Comparison of different decoding schemes for interleaved KK subspace codes. code that is corrupted by ρ row erasures, κ column erasures and t full rank errors. A full rank error t is a deletion and an insertion at the same location (i.e. same code locator α i ). Column erasures correspond to deletions at locations where no insertion occurred and row erasures correspond to insertions at positions where no deletion occurred. For comparison we only consider the homogeneous case k (j) = k for all j ∈ [1, L] since the scheme in [15] only can construct homogeneous interleaved KK subspace codes. In terms of row/column erasures and full rank errors we can express γ = κ + t and δ = ρ + t. The decoding region for the probabilistic unique decoder from Section 4 is then κ + t + L(ρ + t) ≤ L(n t − k)
The decoding regions and computational decoding complexity for known unique decoding schemes are compared in Table 2.
Recall that the interleaving order L depends on the code rate R (see (10)). Further, the impact of L on the decoding performance is quite large in practice. Thus we include L in the O-notation of the computational complexity analysis. For a detailed discussion whether to include L in the O-notation we refer to [25,Chapter 1.2]. Table 2 shows that the rank-metric based decoding schemes [15,20,31] can correct more full errors with increasing L while the performance for insertions and deletions remains the same as for non-interleaved codes [18,33].
The interpolation-based decoding schemes from Section 4 and [15] can correct the same number of full errors as the rank-metric based schemes [15,20,31] but can correct L-times more insertions.
The codes in [15] over the field F q Lm with code locators from F q m are homogeneous L-interleaved KK subspace codes over F q m [6]. In this scheme all operations are performed in the larger field F q Lm which requires higher computational complexity (≤ O(L 4 n 2 r ) ≈ O(L 6 n 2 t ) operations in F q m ). The syndrome-based scheme in [4] achieves the same decoding region as the proposed interpolation-based decoder from Section 4 requiring at most O(L 3 n 3 t ) operations in F q m . In the worst case we have γ ≤ Ln t and thus the computational complexity of the interpolation-based scheme is in the order of O(L 2 n 2 r ) < O(L 4 n 2 t ) operations in F q m .
The proposed interpolation-based decoding scheme therefore achieves the so far best-known decoding region with significantly lower complexity than [15] and lower computational complexity than [4] for n t > L. The upper bound on the failure probability is roughly 4q −m for all schemes.

Conclusion
An interpolation-based decoding scheme for interleaved KK subspace codes has been presented. We have shown that interleaved subspace codes can be made more resilient against insertions as compared to the approach from [18]. Our principle can be used as a (not necessarily polynomial-time) list decoder as well as a probabilistic unique decoder. In both cases, the procedure consists of interpolating a minimal Gröbner basis for the interpolation module followed by a root-finding step.
Two procedures for efficiently constructing required minimal Gröbner basis polynomials were proposed. The procedures are based on the general linearized Kötter interpolation and row reduction in modules. Further, a computationally-and memory-efficient root-finding algorithm for the unique decoder was presented, which exploits the structure of the minimal Gröbner basis.
The results and algorithms can also be used to accelerate interpolation-based decoding for interleaved Gabidulin codes from [39].