1 Towards An Exact Combinatorial Algorithm for LP Decoding of Turbo Codes Michael Helmling and Stefan Ruzika Mathematical Institute University of Koblenz-Landau, Campus Koblenz 56070 Koblenz, Germany Email: {helmling, ruzika}@uni-koblenz.de Abstract—We present a novel algorithm that solves the turbo ML (LP) decoding is equivalent to solving a (LP-relaxed) code LP decoding problem in a fininte number of steps by shortest path problem with additional linear side constraints. Euclideandistanceminimizations,whichinturnrelyonrepeated Sofar,twomethodsforsolvingtheLPhavebeenproposed: shortest path computations in the trellis graph representing the General purpose LP solvers like CPLEX [7] are based on the turbocode.Previousattemptstoexploitthecombinatorialgraph structure only led to algorithms which are either of heuristic matrix representation of the LP problem. They utilize either natureordonotguaranteefiniteconvergence.Anumericalstudy the simplex method or interior point approaches [8], but do 3 showsthatouralgorithmclearlybeatstherunningtime,uptoa not exploit any structural properties of the specific problem. 1 factorof100,ofgenericcommercialLPsolversformedium-sized 0 Lagrangian relaxation in conjunction with subgradient opti- codes, especially for high SNR values. 2 mization [1], [9], on the other hand, utilizes this structure, n Index Terms—LP decoding, turbo codes, combinatorial opti- buthaspracticallimitations,mostnotablyitusuallyconverges mization a very slowly. J This paper presents a new approach to solve the LP 7 I. INTRODUCTION decoding problem exactly by an algorithm that exploits its 2 graphical substructure, thus combining the analytical power Since its introduction by Feldman et al. in 2002 [1], Linear of the LP approach with the running-time benefits of a com- ] Programming based channel decoding has gained tremendous T binatorial method which seems to be a necessary requirement interestbecauseofitsanalyticalpower—LPdecodingexhibits I for practical implementation. Our basic idea is to construct . themaximumlikelihood(ML)certificateproperty[2],andthe s an alternative polytope in the space defined by the additional c decoding behavior is completely determined by the explicitly constraints (called constraints space) and show how the LP [ described “fundamental” polytope [3]—combined with note- solution corresponds to a specific point zQ of that polytope. 1 worthy error-correcting performance and the availability of Then, we show how to computationally LfiPnd zQ by a geo- v efficient decoding algorithms. LP metric algorithm that relies on a sequence of shortest path 3 Turbo codes, invented by Berrou et al. in 1993 [4], are a computations in the trellis graphs. 6 class of concatenated convolutional codes that, together with 3 The reinterpretation of constrained optimization problems a heuristic iterative decoding algorithm, feature remarkable 6 in constraints space was first developed in the context of 1. error-correcting performance. multicriteria optimization in [10], where it is applied to mini- WhilethefirstpaperonLPdecoding[1]actuallydealtwith 0 mum spanning tree problems with a single side constraint. In turbo codes, the majority of publications in the area of LP 3 2010, Tanatmis [11] applied this theory to the turbo decoding 1 decodingnowfocusonLDPCcodes[5]whichprovidesimilar problem. His algorithm showed a drastic speedup compared : performance (cf. [6] for a recent overview). Nevertheless, v to a general purpose LP solver, however it only works for up turbo codes have some analytical advantages, most impor- i totwoconstraints,whileinreal-worldturbocodesthenumber X tantly the inherent combinatorial structure by means of the of constraints equals the information length. r trellis graph representations of the underlying convolutional a By adapting an algorithm by Wolfe [12] to compute in a encoders. ML Decoding of turbo codes is closely related polytope the point with minimum Euclidean norm, we are to shortest path and minimum network flow problems, both able to overcome these limitations and decode turbo codes being classical, well-studied topics in optimization theory for with lengths of practical interest. The algorithm is, compared which plenty efficient solution methods exist. The hardness to previous methods, advantageous not only in terms of of ML decoding is caused by additional conditions on the runningtime,butalsogivesvaluableinformationthatcanhelp path through the trellis graphs (they are termed agreeability to improve the error-correcting performance. Furthermore, constraints in [1]) posed by the turbo code’s interleaver. Thus branch-and-bound methods for integer programming-based ML decoding depend upon fast lower bound computations, The authors are with the Optimization Research Group, Department of Mathematics, University of Kaiserslautern, Erwin-Schroedinger-Strasse, mostly givenby LPrelaxations, andcan often besignificantly 67663Kaiserslautern,Germany.Email:{helmling,ruzika}@mathematik.uni- improved by dedicated methods that evaluate combinatorial kl.de properties of the LP solutions. Since our LP decoder contains WewouldliketoacknowledgetheCenterforMathematicalandComputa- tionalModelling(CM)andtheDAAD such information, it could also be considered a step towards 2 x 0 0 0 0 0 0 0 0 x∈Fk2 C eC(x) 1 π 1 1 1 1 1 1 C eC(π(x)) 1 0 2 2 2 2 2 2 Fig. 1: Turbo encoder with two convolutional encoders C , a Cb and interleaver π. 3 3 1 3 3 3 3 S S 1 2 IP-based algorithms with the potential of practical implemen- Fig. 2: Excerpt from a trellis graph with four states and tation. initial state 0. The style of an edge indicates the according information bit, while the labels refer to the single parity bit. II. BACKGROUNDANDNOTATION A. Definition of Turbo Codes Ak-dimensionalsubspaceC ofthevectorspaceFn (where input and output bit, respectively. The following description 2 F2 ={0,1}denotesthebinaryfield),iscalledan(n,k)binary ofconvolutionalcodeslooselyfollows[6,SectionV.C],albeit linear block code, where n is the block length and k the the notation slightly differs. information (or input) length. One way to define a code is by We denote a trellis by T = (V,E) with vertex set V and an appropriate encoding function eC, for which any bijective edge set E. Vertices are indexed by time step and state; i.e., linear mapping from Fk2 onto C qualifies. This paper deals vi,s denotesthevertexcorrespondingtostates∈{0,...,2d− with turbo codes [4], a special class of block codes built by 1}attimei∈{1,...,k+d+1}.Anedgeinturnisidentified interconnecting (at least) two convolutional codes (see e.g. by the time and state of its tail vertex plus its input label, so [13]). For the sake of clear notation, we focus on turbo codes e denotes the edge outgoing from v with input bit b ∈ i,s,b i,s asusedinthe3GPPLTEstandard[14]—i.e.,systematic,par- {0,1}. We call vertical “slices”, i.e., the subgraphs induced allely concatenated turbo codes with two identical terminated by the edges of a single time step, segments of the trellis. rate-1constituentencoders—despitethefactthatourapproach Formally, the segment at time i is is applicable to arbitrary turbo coding schemes. An in-depth covering of turbo code construction can be found in [15]. S =(V ,E ) i i i An(n,k)turbocodeTC =TC(C,π)isdefinedbyarate-1 where V ={v ∈V :j ∈{i,i+1}} i j,s convolutional (n ,k) code C with constraint length d and a permutation π ∈CS such that n = k +2·n . Because we and Ei ={ej,s,b ∈E :j =i}. k C consider terminated convolutional codes only (i.e., there is a Becausetheinitialandfinalstateoftheconvolutionalencoder designated terminal state of the encoder), the final d bits of arefixed,theleadingaswellasthetrailingdsegmentscontain the information sequence (also called the tail) are not free to less than 2d vertices. Figure 2 shows the first few segments choose and thus can not carry any information. Consequently, of a trellis with d=2. those bits together with the corresponding d output bits are By construction, the paths from the starting node to the considered part of the output, which yields n = k +2·d C end node in a trellis of a convolutional code C are in one-to- and a code rate slightly below 1. Let e :Fk −→FnC be the C 2 2 one correspondence with the codewords of C: Let I ⊂ E associated encoding function. Then, the encoding function of j j and O ⊂ E be those edges of S whose input label TC is defined as j j j and output label, respectively, is a 1. The correspondence eTC :Fk2 −→Fk2+2·nC between a codeword y ∈ Fk2+2·d and the according path e (x)=(x|e (x)|e (π(x))) P =(e ,...,e ) in T is given by TC C C 1 k+d where π(x) = (x ,...,x ). In other words, the code- (cid:40) π(1) π(k) e ∈I for 1≤i≤d word for an input word x is obtained by concatenating y =1⇔ k+i k+i (1) i e ∈O for d<i≤k+2·d, • a copy of x itself, i−d i−d • a copy of x encoded by C, and where the first part accounts for the d “input” tail bits that are • acopyofx,permutedbyπandencodedbyC afterwards. prepended by convention. From (1), for each e∈E an index Figure 1 shows a circuit-type visualization of this definition. set J (e) can be computed with the property that e ∈ P ⇒ C y = 1 for all j ∈ J (e). In our case, |J (e)| varies from 0 j C C B. Trellis Graphs of Convolutional Codes (for edges in S , i ≤ k, with output label 0) to 2 (for edges i A convolutional code with a specific length is represented in Si, k+1≤i≤k+d, with both input and output label 1). naturally by its trellis graph, which is obtained by unfolding The path-codeword relation can be exploited for maxi- thecode-definingfinitestatemachineinthetimedomain:Each mum likelihood (ML) decoding, if the codewords are trans- vertex of the trellis represents the state at a specific point in mitted through a memoryless binary-input output-symmetric time, while edges correspond to valid transitions between two (MBIOS) channel: Let λ ∈ Rk+2·d be the vector of LLR subsequent states and exhibit labels with the corresponding valuesof thereceived signal.Ifwe assignto eachedge e∈E 3 the cost By also transforming (2) and (4), we obtain (cid:88) (cid:88) c(e)= λ , (TC-IP) min c(e)·f j e j∈JC(e) e∈E1∪E2 s.t. f1 ∈P1 (8) path it can be shown [2] that the shortest path in T corresponds to f2 ∈P2 (9) the ML codeword. path (cid:88) (cid:88) f = f i=1,...,k(10) e e e∈I1 e∈I2 i π(i) C. Trellis Representation of Turbo Codes f ∈{0,1}, e∈E. (11) e For turbo codes, we have two isomorphic trellis graphs, D. Polyhedral Theory Background T1 and T2, according to the two component convolutional Besidescodingtheory,thispaperrequiressomebitsofpoly- encoders. Let formally T = (G1 ∪ G2,E1 ∪ E2), and by hedraltheory.Apolytopeistheconvexhullofafinitenumber P = P1 ◦ P2 denote the path that consists of P1 in of points: P = conv(v ,...,v ). It can be described either T1 and P2 in T2. Only certain paths, called agreeable, 1 n by its vertices (or extreme points), i.e., the unique minimal actually correspond to codewords; namely, an agreeable path set fulfilling this defining property, or as the intersection of a P ◦ P = (e1,...,e1 ,e2,...,e2 ) must obey the k 1 2 1 k+d 1 k+d finite number of halfspaces: P = (cid:84)m {x : aTx ≤ b }. An consistency constraints i=1 i i inequality aTx ≤ b is called valid for P if it is true for all x ∈ P. In that case, the set F = {x ∈ P : aTx = b} is e1 ∈I1 ⇔ e2 ∈I2 for i=1,...,k (2) a,b i i π(i) π(i) called the face induced by the inequality. For any r satisfying aTr ≥ b (aTr > b) we say that the inequality separates because both encoders operate on the same information word, (strongly separates) r from P. onlythat itispermuted forthesecond encoder.Consequently, ML decoding for turbo codes can be formulated as finding III. THELPRELAXATIONANDCONVENTIONAL the shortest agreeable path in T. If an agreeable path contains SOLUTIONMETHODS e1 ∈I1,itmustalsocontaine2 ∈I2 ,andthusi∈J (e) i i π(i) π(i) C ML decoding of general linear block codes is known to be for both e1 and e2 . To avoid counting the LLR value λ i π(i) i NP-hard [17]. While the computational complexity of TC- twice in the objective function, we use the modified cost IP is still open, it is widely believed that this problem is NP-hardaswell,whichwouldimplythatnopolynomial-time (cid:40) c(e)= (cid:88) λˆ with λˆ = λ2j if 1≤j ≤k, (3) algorithmcansolveTC-IPunlessP=NP1.Byrelaxing(11) j j λ otherwise. to f ∈[0,1], we get the LP relaxation (referred to as TC-LP) j e j∈JC(e) of the integer program TC-IP, which in contrast can be solved efficiently by the simplex method or interior point approaches Then,theMLdecodingproblemforturbocodescanbestated [8]. Feldman et al. [1] were the first to analyze this relaxation as the combinatorial optimization problem and attested it reasonable decoding performance. (cid:88) A general purpose LP solver, however, does not make use (TC-ML) min c(e) (4) of the combinatorial substructure contained in TC-IP via (8) e∈P=P1∪P2 and(9)andthuswastessomepotentialofsolvingtheproblem s.t. P1 is a path in T1 (5) more efficiently—while LPs are solvable in polynomial time, P2 is a path in T2 (6) theydonotscaletoowell,andthenumberofvariables(about P is agreeable 2·|V|=(k+d)·2d+2) and constraints (|V|+k) in TC-LP is very large (practical values of d range roughly from 3 to 8). The codeword variables y can be included into (TC-ML) by Notethatwithouttheconsistencyconstraints(10),wecould i the constraints solveTC-LPbysimplycomputingshortestpathsinbothtrellis graphs, which is possible in time O(k + d), even in the  (cid:80) fe for 1≤i≤k presence of negative weights, because the graphs are acyclic  2 y = JC(e)(cid:51)i (7) [19]. A popular approach for solving optimization problems i (cid:80)  fe for i>k that comprise “easy” subproblem plus some “complicating” JC(e)(cid:51)i additional constraints is to solve the Lagrangian dual [20] by (cid:80) subgradient optimization. If we define g (f) = f − where the factor 21 is analogical to (3). However, these (cid:80)e∈I2 fe, the constraints (10) can beicompactlyer∈eIwi1riteten variables are purely auxiliary in the LP and thus not needed. as π(i) ItisstraightforwardtoformulateTC-MLasanintegerlinear g (f)=0 for i=1,...,k. i program by introducing a binary flow variable f ∈{0,1} for e each e∈E1∪E2. The constraints (5) and (6) can be restated 1Note that with state-of-the-art software and prohibitive computational in terms of flow conservation and capacity constraints [16] effort, ML turbo decoding can be simulated off-line on desktop computers; which define the path polytopes P1 and P2 , respectively. see[18] path path 4 The Lagrangian relaxation with multiplier µ∈Rk is defined 2) If f represents an agreeable path in T, then D(f) is as located on the (k+1)st axis (henceforth called c-axis or A ). k c (TC-LR) z(µ)=min (cid:88) c(e)·fe+(cid:88)µk·gi(f) 3) If v is a vertex of Q and v =D(f) for some f ∈Ppath, then f is also a vertex of P . e∈E1∪E2 i=1 path s.t. f1 ∈P1 Inthesituationthatv =D(f)wewillalsowritef =D−1(v) path withthemeaningthat f isany preimageof v,whichneednot f2 ∈P2 path be unique. f ∈{0,1}, e∈E e We consider the auxiliary problem For all µ ∈ Rk, the objective value of TC-LR is smaller (TC-LP ) zQ =min v Q LP k+1 or equal to that of TC-LP. The Lagrangian dual problem s.t. v ∈Q is to find multipliers µ that maximize this objective, thus v ∈A (12) minimizing the gap to the LP solution. It can be shown c that in the optimal case both values coincide. Note that the the solution of which is the lower “piercing point” of the axis feasible region of TC-LR is the combined path polytope of A throughQ.Notethatdueto(12),kofthek+1variablesin c both T and T , so it can be solved by a shortest path 1 2 TC-LP(Q)arefixedtozero,thustheproblemisinasenseone- routineinbothtrelliseswithmodifiedcosts,andtheintegrality dimensional, the feasible region being the (one-dimensional) conditiononf isfulfilledautomatically.ApplyingLagrangian projection of Q onto A . Nevertheless, the following theroem c relaxationtoturbodecodingwasalreadyproposedbyFeldman shows that TC-LP and TC-LP are essentially equivalent. Q etal.[1]andfurtherelaboratedbyTanatmisetal.[9];thelatter Theorem 1: Let v be an optimal solution of TC-LP LP Q referencecombinestheapproachwithaheuristictotightenthe with objective value zQ and f = D−1(v ) ∈ P the LP LP LP path integrality gap between TC-LP and TC-IP. corresponding flow. Then zQ = z , the optimal objective LP LP The Lagrangian dual is typically solved by a subgradient value of TC-LP, and f is an optimal solution of TC-LP. LP algorithm that iteratively adjusts the multipliers µ, converging Proof: First we show zQ ≤ z . Let f be an optimal LP LP LP (under some mild conditions) to the optimal value [20]. solution of TC-LP with cost c(f ) = z . Then D(f ) = LP LP LP However, the convergence is often slow in practice and the (0,...,0,z ) by definition of D, since f is feasible and LP LP limitisnotguaranteedtobeeverreachedexactly.Additionally, thus g (f ) = ··· = g (f ) = 0. Hence D(f ) ∈ A ∩Q 1 LP k LP LP c the dual only informs us about the objective value; recovering with D(f ) =z , from which it follows that zQ ≤z . LP k+1 LP LP LP the actual solution of the problem requires additional work. If we assume on the other hand that zQ <z , there must LP LP In summary, subgradient algorithms suffer from three major be a v ∈ A ∩Q such that v < z . By definition of D c k+1 LP flaws.Themainresultofthispaperisanalternativealgorithm thisimpliestheexistenceofaflowf =D−1(v)withg (f)= 1 which exhibits none of these. ··· = g (f) = 0, hence a feasible one, and c(f) = v < k k+1 z , contradicting optimality of z . LP LP IV. ANEQUIVALENTPROBLEMINCONSTRAINTSSPACE While we do not have an explicit representation of Q—by meansofeitherverticesorinequalities—athand,wecaneasily Like Lagrangian dualization, our algorithm also uses a minimize linear functionals over Q: relaxedformulationofTC-IPwithmodifiedobjectivefunction Observation 1: The problem that resembles TC-LR. However, via geometric interpretation of the image of the path polytope in the “constraints space”, (LP ) minγTv Q as defined below, the exact LP solution is found in finitely s.t. v ∈Q many steps. can be solved by first computing an optimal solution f∗ of the weighted sum problem A. The Image Polytope Q Let P = P1 ×P2 be the feasible region of TC-LR. (cid:88)k path path path (TC-WS) min γ ·g (f)+γ ·c(f) We define the map i i k+1 i=1 D:Ppath →Rk+1 s.t. f ∈Ppath f (cid:55)→(g1(f),...,gk(f),c(f))T and then taking the image of f∗ under D. As noted before, (cid:80) this can be achieved within running time Ø(n). where c(f) = c(e) · f is a short hand for the e∈E1∪E2 e Note that TC-WS is closely related to TC-LR: as long as objective function value of TC-LP. For a path f, the first k γ (cid:54)= 0, we get the same problem by setting µ = γi coordinates vi,i=1,...,k, of v =D(f) tell if and how the k+1 i γk+1 in TC-LR. conditiong (f)=0isviolated,whilethelastcoordinatev i k+1 equals the cost of f. Let Q = D(P ) be the image of the path path polytope under D. The following results are immediate: B. Solving TC-LP with Nearest Point Calculations Q Lemma 1: Our algorithm solves TC-LP by a series of nearest point Q 1) Q is a polytope. computations between Q and reference points ri on A , the c 5 lastofwhichgivesafaceofQcontainingtheoptimalsolution Case 1: a(r)Ty < b(r). Then y ∈/ NF(r), but y ∈ N N N v . NF(s) by definition, which proves the claim for this case. LP For each r ∈Rk+1, we denote by Case2:a(r)Ty =b(r).Froma(r)Tr >br anda(r)Ts= N b(r) we obtain NP(r)=argmin(cid:107)v−r(cid:107) 2 v∈Q a(r)Tr >a(r)Ts the nearest point to r in Q with respect to Euclidean norm ⇒a(r)T(r−s)>0 and define ⇒a(r) (r −z )>0 k+1 k+1 k+1 a(r)=r−NP(r) ⇒r <s , (14) k+1 k+1 b(r)=a(r)T NP(r). where we have used again a(r) < 0 and the fact that k+1 The following well-known result will be used frequently r,s∈A . c below. Applying Lemma 3 to s as reference point we obtain Lemma 2: The inequality a(s) =(s−y ) <0, hence k+1 N k+1 a(r)Tv ≤b(r) (13) (y ) >s N k+1 k+1 isvalidforQandinducesafacecontainingNP(r),whichwe ⇒(y ) (s −r )>s (s −r ) by (14) N k+1 k+1 k+1 k+1 k+1 k+1 call NF(r). If r ∈/ Q, (13) strongly separates r from Q. ⇒yT(s−r)>sT(s−r) N The following theorem is the foundation of our algorithm. ⇒yTs−yTr+sTr−sTs>0 (15) Theorem 2: There exists an ε>0 such that for all r inside N N theopenlinesegment(cid:0)v ,v −(0,...,0,ε)T(cid:1)thecondition LP LP Plugging the definitions into a(r)Ty = b(r) = a(r)Ts N v ∈NF(r) yields (r − x )Ty = (r − NP(r))Ts or rTs − rTy = LP N N N NP(r)Ts−NP(r)Ty .Usingthiswecontinuefrom(15)with N holds. Our constructive proof of Theorem 2 shows how find a point ⇒yTs+NP(r)Ts−NP(r)Ty −sTs>0 N N inside the interval mentioned in the theorem. The outline is ⇒ NP(r)T(s−y )>sT(s−y ) as follows: At first, start with a reference point r ∈ A that N N c is guaranteed to be located below v . Then, we iteratively ⇒a(s)T NP(r)>a(s)Ts>b(s) LP compute NF(r) and update r to be the intersection of A c Thus, NP(r) ∈/ NF(s) = {v ∈ Q : a(s)Tv ≤ b(s)}, but with the hyperplane defining NF(r). The following lemmas NP(r)∈NF(r) by definition, so those faces must differ. show that this procedure is valid and finite. Now we show the auxiliary result that if two inequalities The first result is that the hyperplane defining NF(r) is induce the same face, then also every convex combination of always oriented “downwards”. them does. Lemma 3: Let r = (0,...,0,ρ)T with ρ < zQ and let LP Lemma 5: Let P be a polytope, x1,x2 ∈ P, and r1,r2 ∈/ a(r)Tv ≤b(r) be the inequality defined in (13). Then, P. If the inequalities a(r) <0. k+1 H1 :(r1−x1)Tx≤(r1−x1)Tx1 Proof: Assuming a(r) ≥ 0, we obtain a(r)Tv = k+1 LP a(r) zQ ≥ a(r) ρ = a(r)Tr > b(r), which contradicts and k+1 LP k+1 H2 :(r2−x2)Tx≤(r2−x2)Tx2 v ∈ Q by Lemma 2. Note that the equalities hold because LP both v and r are elements of A , the first inequality stems LP c both induce the same face F of P, then also from the assumptions on a(r) and ρ, and the second k+1 follows from Lemma 2. H¯ :(r¯−x¯)Tx≤(r¯−x¯)Tx¯ Nextweshowthatupdatingr leadstoadifferentnearestface, with r¯=λr1+(1−λ)r2, x¯=λx2+(1−λ)x2, 0≤λ≤1, unless we have arrived at the optimal solution. is valid and induces F. Lemma 4: Under the same assumptions as in Lemma 3, let Proof: We first show that H¯ is valid. For x∈P s=(0,...,0,s ) with s = b(r) be the point where k+1 k+1 a(r)k+1 theseparatinghyperplaneandA intersect.IfNF(r)=NF(s), (r¯−x¯)Tx=λ(r1−x1)Tx+(1−λ)(r2−x2)Tx c then s=vLP. ≤λ(r1−x1)Tx1+(1−λ)(r2−x2)Tx2 Proof: We use contraposition to show that s (cid:54)= v implies NF(r) (cid:54)= NF(s), so assume s (cid:54)= v . We know thLaPt =(cid:0)λ(r1−x1)+(1−λ)(r2−x2)(cid:1)T LP a(r)Tv ≤ b(r) is valid for Q and a(r)Ts = a(r)k+1sk+1 = (cid:0)λx1+(1−λ)x2(cid:1) b(r) by construction. This implies that s ∈/ Q; otherwise =(r¯−x¯)Tx¯, we would have s = v because for all ζ < s , LP k+1 a(r)T(ζe ) > b(r), so s would really be the lowest point where we have used the fact that H1 is satisfied with equality k+1 on A that is also in Q and thus optimal. for x=x2 and vice versa because of the assumptions. Since c It follows that y = NP(s) (cid:54)= s. Since y ∈ Q and we have shown that H¯ is valid, it must induce a face F¯. It N N a(r)Tv ≤b(r) is valid for Q, we have a(r)Ty ≤b(r). remains to show that F =F¯. N 6 “F ⊆ F¯”: x ∈ F fulfills both H1 and H2 with equality, An illustration of the process in two dimensions is given in so we can carry out the above calculation with a “=” in the Figure 3. second line to conclude x∈F¯. “F¯ ⊆ F”: Let x ∈ F¯ and assume x ∈/ F, which implies C. Solving the Nearest Point Problems (ri − xi)Tx < (ri − xi)Txi for i ∈ {1,2}. Then λ(r1 − x1)Tx1+(1−λ)(r2−x2)Tx2 >λ(r1−x1)Tx+(1−λ)(r2− It remains to show how to solve the nearest point problems x2)Tx = (r¯−x¯)Tx = (r¯−x¯)Tx¯ = λ(r1 −x1)Tx1 +(1− arising in the discussion above. To that end, we utilize an λ)(r2−x2)Tx2, which is a contradiction. algorithm by Wolfe [12] that finds in a polytope the point The above lemma is used to show that the part of A that lies with minimum Euclidean norm. Wolfe’s algorithm elaborates c below v dissects into intervals such that reference points on a set of vertices of the polytope that are obtained via LP within one interval yield the same face of Q. minimization of linear objective functions. In our situation, Lemma 6: If r1 < r2 < x∗ and NF(r1) = NF(r2), then this means that LPQ has to be solved repeatedly, which by d NF(r)=NF(r1) for all r ∈[r1,r2]. Observation1boilsdowntothelinear-timesolvableweighted Proof: Let vi =NP(ri) for i∈{1,2}. By Lemma 5, for sum shortest path problem TC-WS. Note that by subtracting eachλ∈(0,1)andr¯=λr1+(1−λ)r2,v¯=λv1+(1−λ)v2, r from the results of LPQ and adding r to the final result, it holds the algorithm can be used to calculate the minimum distance between Q and r also in the case r (cid:54)=0. {v ∈Q:(r¯−v¯)Tv =(r¯−v¯)Tv¯}=NF(r1), The algorithm in [12] maintains in each iteration a subset and applying the converse statement from Lemma 2 follows P of the vertex set V(Q) and a point x such that x = v¯=NP(r¯), so NF(r¯)=NF(r1) as claimed. NP(aff(P)) lies in the relative interior of conv(P), where Nowwehavealletheingredientsathandtoproveourtheorem. aff(P) is the affine hull of P. Such a set is called a corral, Proof of Theorem 2: First we show that there exists at and we denote the nearest point in aff(P) by vaff. P least one r with the desired properties. Initially P = {v0} for an arbitrary vertex v0 and x = v0. Choosesomearbitraryr0 ∈A withr0 <zQ (thusr0 ∈/ Note that then vaff = v and P is indeed a corral. Then the c k+1 LP P 0 Q). If v ∈NF(r0), we are done. Otherwise, Lemma 4 tells following is executed iteratively (we explain afterwards how LP us how to find an r1 with r1 > r0 such that NF(r1) (cid:54)= the computations are actually carried out): k+1 k+1 NF(r0).IteratingthisargumentandassumingthatvLP isnever 1) Solve p=argminv∈Q(xTv). contained in the induced face results in a sequence (ri)i with 2) If p = 0 (0 is optimal) or xTp = xTx (x is optimal), ri+1 > ri for all i. Because of Lemma 6, NF(ri+1) (cid:54)= stop.Otherwise,setP :=P∪{p}andcomputey =vaff. k+1 k+1 P NF(ri) implies NF(ri+1)(cid:54)=NF(rl) for all 0≤l<i+1, so 3) If y is in the relative interior of conv(P), P is a corral. that all NF(ri) are distinct. But since there are only finitely Set x:=y and continue at 1). manyfacesofQ,thiscannotbetrue,soeventuallytheremust 4) Determine z ∈ conv(P)∩conv{x,y} with minimum be an i∗ such that vLP ∈NF(ri∗). distance to y; z will be a boundary point of conv(P). Now let r∗ ∈ Ac be any such point whose existence we 5) Remove from P some point that is not on the smallest have just proven, vN = NP(r∗) and λ ∈ (0,1]. Let r¯ = face of conv(P) containing z, set x:=z, and continue λr∗+(1−λ)vLP and v¯=λvN +(1−λ)vLP. We use similar at 3). arguments as in the proof of Lemma 5 to show that The algorithm will eventually find a corral P such that the (r¯−v¯)Tv ≤(r¯−v¯)Tv¯ (16) nearest point of Q equals vaff. P The computations in each step are performed as follows: induces NF(r∗). 1) This matches the solution of TC-WS. For v ∈Q, 2) IfweinterchangeablyusethesymbolP forboththeset (r¯−v¯)Tv =(λ(r∗−vN)+(1−λ)(vLP−vLP))Tv ofpointsandthematrixthatcontainstheelementsofP =λ(r∗−v )Tv as columns, every v ∈ aff(P) can be characterized by N some λ ∈ R|P| such that v = Pλ and eTλ = 1. Thus, ≤λ(r∗−v )Tv N N the subproblem of determining vaff can be written as =λ(λ(r∗−v )Tv +(1−λ)(r∗−v )Tv ) P N N N LP =λ(r∗−vN)T(λvN +(1−λ)vLP) min (cid:107)Pλ(cid:107)22 =λTPTPλ =(r¯−v¯)Tv¯. s.t. eTλ=1. So the inequality is valid, and since again for v ∈ NF(r∗) It can be shown [12] that this is equivalent to solving equality holds in the third line, we know that the face F¯ the system of linear equations induced by (16) contains NF(r∗). (eeT +PTP)µ=e Nowletv ∈F¯,i.e.,(r¯−x¯)Tv =(r¯−x¯)Tx¯.Fromtheabove 1 (17) equations we conclude λ(r∗−v )Tv =λ(r∗−v )Tv , and λ= µ N N N (cid:107)µ(cid:107) because λ>0 this implies v ∈NF(r∗). 1 Becausetheaboveholdsforany0<λ≤1,wecanchoose As an efficient method to solve (17), Wolfe suggests to r¯ arbitrarily close to v on A , which completes the proof. maintainanuppertriangularmatrixRsuchthatRTR= LP c eeT +PTP. Then the solution µ can be found by first 7 A A c c Q Q ri+1 vi+1 ri ri vi vi (a) Step i: vi =D(fi) is found as nearest point to some refer- (b)Stepi+1:NotethattheinducedfaceofQhereisafacet, ence point ri−1. The intersection of the separating hyperplane while it was a 0-dimensional face in step i. with the axis A , r1, will be the reference point of the next c iteration. A A c c Q Q ri+2 =v ri+2 =vi+3 =v LP LP vi+2 ri+1 vi+2 vi+1 vi+1 (c) Step i + 2 (zoomed in): The facet induced in this step (d) Step i+3: Optimality is detected by vi+3 = ri+2. The intersectsAc atvLP,butthealgorithmcannotyetdetectthis. solution D−1(vLP) is returned. Fig. 3: Schematic execution of Algorithm 1 in image space solving RTµ¯ = e for µ¯ and then Rµ = µ¯ for µ; both 5) A point not contained in the smallest face of conv(P) can be done by a simple backward substitution. When containing z is not needed for the convex description of (cid:80) P changes, R can be updated relatively easily without z = λ v; thus it can be identified by λ =0. v∈P v v the necessity of a complete recomputation [12]. 3) y is in the relative interior of conv(P) if and only if all D. Recovering the Optimal Flow and Pseudocodeword coefficients λ in the convex representation of y satisfy i λ >0. Sofarwehaveshownhowtocomputetheoptimalobjective i (cid:80) 4) By construction x ∈ conv(P). Let x = λ v and value, but not the LP solution, i.e. the flow f ∈ P and v∈P v LP path (cid:80) (cid:80) (cid:80) y = µ v, where λ = µ = 1, but the (pseudo)codeword y. The algorithm yields its solution v v∈P v v∈P v v∈P v LP µ ≤0 for at least one p. The goal can then be restated by means of a convex combination of extreme points of Q: p as finding the minimal θ ∈ [0,1] such that z = θx+ θ t t (1−θ)y ∈conv(P). Substituting the above expressions (cid:88) (cid:88) v = λ v , λ ≥0, λ =1. yields LP i i i i (cid:88) i=1 i=1 z = (θλ +(1−θ)µ )v, θ v v During its execution the preimage paths f =D−1(v ) can be v∈P i i stored alongside with the v . Then, the LP-optimal flow f andtheconditionisthatallcoefficientsarenonnegative. i LP is obtained by summing up the paths with the same weight Thus, for all v with µ ≤0, v coefficients λ, i.e., µ θ ≥ p µ −v t p p (cid:88) f = λ D−1(v ). LP i i must hold. In summary, θ can be computed by the rule i=1 (cid:26) (cid:26) (cid:27)(cid:27) µ θ =min 1,max p : µ <0 . In order to get the corresponding pseudocodeword, a simple µ −v p p p computation based on (7) suffices. 8 For most applications, however, the values of y are of Algorithm 1 Combinatorial Turbo LP Decoder (CTLP) interest only in the case that the decoder has found a valid 1: Initialize edge cost c(f) by (3) codeword, i.e., t = 1 in the above sum. In such a case, the 2: f0 ←argmin{c(f): f ∈Ppath}. most recent solution of (TC-WS) is an agreeable path that 3: v0 ←D(f0). immediately gives the codeword. No intermediate paths have 4: r0 ←(0,...,0,v0 )T k+1 to be stored, which can save a substantial amount of space 5: i←0 and running time. 6: while vi (cid:54)=ri do (cid:13) (cid:13) 7: vi+1 ←NP(ri)=argmin (cid:13)v−ri(cid:13) v∈Q 2 E. Efficient Reference Point Updates 8: fi+1 =D−1(vi+1) 9: ai+1 ←vi+1−ri As suggested by the proof of Theorem 2, the nearest point 10: bi+1 ←a(i+1)Tvi+1 algorithm is run iteratively, and between two runs the k+1st (cid:18) (cid:19)T component of r is increased by means of the rule rk+1 = 11: ri+1 ← 0,...,0,abii++11 a(br()rk)+1. This section describes how some information from 12: i←i+1 k+1 the previous iteration can be re-used to provide an efficient 13: end while warm start for the next nearest point run. 14: return fi AssumethatiniterationithepointNP(ri)=vi+1 hasbeen found, inducing the face NF(ri) defined by a(ri)Tv ≤ b(ri) of Q. Recall that NPA internally computes the minimum l norm of Q−ri. Thus, it outputs v¯i+1 = vi+1 −ri as th2e VI. NUMERICALRESULTS convex combination of t ≤ k+1 points v¯ = v −ri, all of A. Running Time Comparison j j whicharelocatedonthecorrespondingfaceNˆF(ri)ofQ−ri: To evaluate the computational performance of our algo- rithm, we compare its running time with the commercial t t (cid:88) (cid:88) v¯i+1 =vi+1−ri = λ (v −ri)= λ v¯ general purpose LP solver CPLEX [7] which is said to be j j j j one of the most competitive implementations available. j=1 j=1 Simulations were run using LTE turbo codes with block- Inthesubsequentnearestpointcalculation,thenormofQ− ri+1 is minimized. Obviously NˆF(ri) corresponds to a face lengths 132, 228, and 396, respectively, and a three- NˆF(ri+1) of Q−ri+1, and we can initialize the algorithm dimensional turbo code with blocklength 384 (taken from [21]) with various SNR values. For each SNR value, we have with that face by simply adding ri−ri+1 to v¯ and each v¯ , j generated up to 105 noisy frames, where the computation was j =1,...,t, which yields stoppedwhen200decodingerrorsoccured.Thisshouldensure (cid:88)t sufficientsignificanceoftheaverageresultsshowninTablesI– v¯i+1+ri−ri+1 =vi+1−ri+1 = λ (v −ri+1) j j IV. j=1 TABLEI:AverageCPUtimeperinstance(in 1 sofseconds) and all vj−ri+1 are vertices of Q−ri+1. Note that ri−ri+1 100 for the (132,40) LTE turbo code is zero in all but the last component, so this update takes only t≤k+1 steps. SNR 0 1 2 3 4 5 In order to warm-start the nearest point algorithm, the CPLEX 9.1 9.5 9.6 9.6 9.6 9.8 CTLP 1.4 0.9 0.5 0.29 0.24 0.22 auxiliary matrix R has to be recomputed as well. Using its ratio 6.5 10.6 19 33 40 45 definition RTR=eeT +VTV TABLEII:AverageCPUtimeperinstance(in 1 sofseconds) we can efficiently compute R by Cholesky decomposition. 10 for the (228,72) LTE turbo code After these updates we can directly start the nearest point algorithm in Step 2. Numerical experiments have shown that SNR 0 1 2 3 4 5 thisspeedsupLPdecodingbyafactoroftwo.Inparticular,the CPLEX(×10−1) 3.1 3.4 4.2 4.6 4.7 4.7 computationtimeoftheCholeskydecompositionisnegligible. CTLP(×10−1) 0.7 0.4 0.15 0.05 0.04 0.04 ratio 4.4 8.5 28 92 118 118 V. THECOMPLETEALGORITHM Algorithm 1 formalizes the procedure developed in Sec- TABLEIII:AverageCPUtimeperinstance(in 1 sofseconds) 10 tion IV in pseudocode. The initial reference point r0 is for the (396,128) LTE turbo code generated by first minimizing c(f) on P (thus solving path SNR 0 1 2 3 4 TC-WS with γ = (0,...,0,1) and projecting the result in CPLEX 4.4 4.2 3.6 3.3 3.2 constraints space onto Ac (Line 4). Thereby we ensure that CTLP 6.3 4.1 0.6 0.09 0.08 either r0 ∈/ Q or it is located on the boundary of Q, in which ratio 0.7 1 6 37 40 case it already is the optimal LP solution. The solution of the nearest point problem and the recovery of the original flow As one can see, the benefit of using the new algorithm is are encapsulated in Lines 7 and 8. largerforhighSNRvalues.Thisbecomesmosteminentforthe 9 TABLE V: Statistical data for the (132,40) LTE turbo code; 10−1 average per-instance counts SNR 0 2 4 facedim 25.2 3.6 0.01 ) integral 0.26 0.89 0.9995 s CPLEX nce( CTLP tmriavjioarl 2021 05.133 0.464 sta mainloops 4.36 1.9 0.7 n e/i 10−2 m 10−2 UpdateR ti SP LstSq 10−3 GenSol ) s ( e 0 1 2 3 4 5 nc a 10−4 SNRb(dB) e/inst m Fig. 4: CPU time comparison for the (132,40) LTE Turbo ti 10−5 code depending on the SNR value (note the logarithmic time scale). 10−6 TABLE IV: Average CPU time per instance (in seconds) for 0 1 2 3 4 5 a (384,128) 3-D turbo code SNR (dB) b SNR 0 1 2 3 4 CPLEX 1.4 1.2 0.9 0.72 0.57 Fig. 5: Average per-instance CPU time spent on various CTLP 4.5 3.1 0.8 0.04 0.014 subroutinesofthealgorithmdecodingthe(132,40)LTEcode ratio 0.31 0.39 1.1 18 41 (SP=shortest path, LstSq=solution of least squares problems, GenSol=generation of solution in path space). 3-D code for which the dimension of Q is the highest, where the new algorithm is slower than CPLEX for SNRs up to 2. B. Numerical Stability The reason for this behavior can be explained by analyzing For larger codes, the dimension of Q becomes very large statisticalinformationaboutvariousinternalparametersofthe which leads to numerical difficulties in the nearest point algorithm when run with different SNR values: algorithm: the equation systems solved during the execution sometimeshaverank“almostzero”whichleadstodivisionby • Theaveragedimensionoftheoptimalnearestface,found very small numbers, resulting in the floating-point value NaN. in the last iteration of the algorithm, drops substantially Careful adjustment of the tolerance values for equivalence with increasing SNR. Intuitively, it is not surprising that checks help to eliminate this problem at least for the block finding a face that needs less vertices to describe can be lengths presented in this numerical study. found more efficiently. In addition, it has proven beneficial to divide all objective • In particular, the share of instances for which the LP values by 10 in advance. Intuitively, this compresses Q along solution is integral (and thus, the face dimension is zero) the c-axis, evening out the extensiveness of the polytope in increases with the SNR. the different dimensions (note that for all axes other than c, • There are some trivial instances where the initial short- the values only range from −1 to 1). est path among both trellis graphs is already a valid codeword. This occurs more often for low channel noise VII. IMPROVINGERROR-CORRECTINGPERFORMANCE and allows for extremely fast solution (no nearest point calculations have to be carried out). As discussed above, Algorithm 1 can be easily modified to • The average number of major cycles of the nearest point returnalistofpathsfi,i=1,...,t,suchthattheLPsolution algorithm performed per instance is seen to drop rapidly isaconvexcombinationofthatpaths.Eachfi canbesplitinto with increasing SNR. a paths fi1 and fi2 through trellis T1 and T2, respectively. • Likewise, the the average number of main loops (Line 6 A path in a trellis, in turn, can uniquely be extended to a of Algorithm 1) drops, reducing the required calls to codeword.Thus,wehaveatotalof2tcandidatecodewords.By CTLP. selecting among them the codeword with minimum objective function value, we obtain a heuristic decoder (Heuristic A in Table V exemplarily contains the average per-instance values the following) that always outputs a valid codeword, and has of these parameters for the (132,40) LTE code and SNRs 0, thepotentialofabettererror-correctingperformancethanpure 2, and 4. LP decoding. 10 100 Finally, the concepts presented here might be of inner- mathematical interest as well. Optimization problems that 10−1 are easy to solve in principle but have some complicating constraints are very common in mathematical optimization. Being able to efficiently solve their LP relaxation is a key 10−2 component of virtually all fast exact or approximate solution R E algorithms. F 10−3 ACKNOWLEDGEMENTS 10−4 WewouldliketoacknowledgetheGermanResearchCoun- LPDecoding cil(DFG),theGermanAcademicExchangeService(DAAD), HeuristicA 10−5 HeuristicB andtheCenterforMathematicalandComputationalModelling MLDecoding ((CM)2) of the University of Kaiserslautern for financial support. 0 1 2 3 4 5 SNR (dB) REFERENCES b [1] J. Feldman, D. R. Karger, and M. Wainwright, “Linear programming- Fig. 6: Decoding performance of the proposed heuristic en- based decoding of turbo-like codes and its relation to iterative ap- hancements on the (132,40) LTE turbo code. proaches,” in Proc. 40th Allerton Conf. Commun. Control Computing, 2002. [2] J.Feldman,“Decodingerror-correctingcodesvialinearprogramming,” Ph.D.dissertation,MassachusettsInstituteofTechnology,2003. [3] P. O. Vontobel and R. Koetter, “Graph-cover decoding and finite- A slightly better decoding performance, at the cost of once lengthanalysisofmessage-passingiterativedecodingofLDPCcodes,” more increased running time, is reached if we consider not arXiv:cs/0512078v1[cs.IT],2005. only the paths that constitute the final LP solution but rather [4] C. Berrou, A. Glavieux, and P. Thitimajshima, “Near shannon limit allintermediatesolutionsofTC-WS.Wecallthismodification error-correctingcodinganddecoding:Turbo-codes,”inIEEEInt.Conf. Commun.,May1993,pp.1064–1070. Heuristic B. [5] R. G. Gallager, “Low-density parity-check codes,” IRE Trans. Inf. Simulation results for the (132,40) LTE code are shown in Theory,vol.8,no.1,pp.21–28,Jan.1962. [6] M.Helmling,S.Ruzika,andA.Tanatmis,“Mathematicalprogramming Figure 6. As one can see, the frame error rate indeed drops decoding of binary linear codes: Theory and algorithms,” IEEE Trans. notably when using the heuristics, but for low SNR values Inf.Theory,vol.58,no.7,pp.4753–4769,Jul.2012. therestillremainsasubstantialgaptotheMLdecodingcurve. [7] “IBM ILOG CPLEX optimization studio,” Software Package, 2011, At 5dB, Heuristic B empirically reaches ML performance; version12.4. [8] A. Schrijver, Theory of linear and integer programming. John Wiley for lower SNR values it is comparable to a Log-MAP turbo &Sons,1986. decoder with 8 iterations. [9] A.Tanatmis,S.Ruzika,andF.Kienle,“ALagrangianrelaxationbased decoding algorithm for LTE turbo codes,” in Proc. Int. Symp. Turbo CodesandIterativeInf.Proc.,Brest,France,Sep.2010,pp.369–373. VIII. CONCLUSIONANDOUTLOOK [10] S. Ruzika, “On multiple objective combinatorial optimization,” Ph.D. dissertation,UniversityofKaiserslautern,2007. We have shown how the inherent combinatorial network- [11] A. Tanatmis, “Mathematical programming approaches for decoding of flow structure of turbo codes in form of the trellis graphs can binary linear codes,” Ph.D. dissertation, University of Kaiserslautern, Kaiserslautern,Germany,Aug.2010. beutilizedtoconstructahighlyefficientLPsolver,specialized [12] P. Wolfe, “Finding the nearest point in a polytope,” Math. Program., forthatclassofcodes.Thedecreaseinrunningtime,compared vol.11,pp.128–149,1976. to a general purpose solver, is dramatic, and in contrast to [13] D. J. C. MacKay, Information Theory, Inference, and Learning Algo- rithms. CambridgeUniversityPress,2003. classical approaches based on Lagrangian dualization, the [14] TS 36.212 v11.0.0: LTE E-UTRA Mutliplexing and Channel algorithm is guaranteed to terminate after a finite number of Coding, 3rd Generation Partnership Project (3GPP) Std., Oct. steps with the exact LP solution. 2012. [Online]. Available: http://www.etsi.org/deliver/etsi ts/136200 136299/136212/11.00.00 60/ts 136212v110000p.pdf It is still an open question, however, if and how the LP [15] S.LinandD.Costello,Jr.,ErrorControlCoding,2nded. UpperSaddle can be solved in a completely combinatorial manner. The River,NJ:Prentice-Hall,Inc.,2004. nearest point algorithm suggested in this paper introduces a [16] R.K.Ahuja,T.L.Magnanti,andJ.B.Orlin,NetworkFlows. Prentice- Hall,1993. numerical component, which is necessary at this time but [17] E. Berlekamp, R. J. McEliece, and H. C. A. van Tilborg, “On the rather undesirable since it can lead to numerical problems in inherent intractability of certain coding problems,” IEEE Trans. Inf. high dimension. Theory,vol.24,no.3,pp.954–972,1978. [18] A.Tanatmis,S.Ruzika,M.Punekar,andF.Kienle,“Numericalcompar- Another direction for further research is to examine the isonofIPformulationsasMLdecoders,”inIEEEInt.Conf.Commun., usefulness of our decoder as a building block of branch-and- 2010,pp.1–5. bound methods that solve the integer programming problem, [19] T.H.Cormen,C.E.Leiserson,R.L.Rivest,andC.Stein,Introduction toAlgorithms,2nded. MITPress,2001. i.e., ML decoders. Several properties of the decoder suggest [20] G. L. Nemhauser and L. A. Wolsey, Integer and Combinatorial Op- that this might be a valuable task. For instance, the shortest timization. Wiley-Interscience series in discrete mathematics and pathscanbecomputedevenfasterifaportionofthevaribales optimization,JohnWiley&Sons,1988. [21] E. Rosnes, M. Helmling, and A. Graell i Amat, “Pseudocodewords of is fixed, or the algorithm could be terminated prematurely if linear programming decoding of 3-dimensional turbo codes,” in Proc. thereferencepointexceedsaknownupperboundatthecurrent IEEEInt.Symp.Inform.Theory,St.Petersburg,Russia,Jul./Aug.2011, node of the branch-and-bound tree. pp.1643–1647.