EURASIP Journal on Applied Signal Processing 2005:8, 1235–1250 c ○ 2005 Hindawi Publishing Corporation On Extended RLS Lattice Adaptive Variants: Error-Feedback, Normalized, and Array-Based Recursions

Error-feedback, normalized, and array-based recursions represent equivalent RLS lattice adaptive filters which are known to offer better numerical properties under finite-precision implementations. This is the case when the underlying data structure arises from a tapped-delay-line model for the input signal. On the other hand, in the context of a more general orthonormality-based input model, these variants have not yet been derived and their behavior under finite precision is unknown. This paper develops several lattice structures for the exponentially weighted RLS problem under orthonormality-based data structures, including error-feedback, normalized, and array-based forms. As a result, besides nonminimality of the new recursions, they present unstable modes as well as hyperbolic rotations, so that the well-known good numerical properties observed in the case of FIR models no longer exist. We verify via simulations that, compared to the standard extended lattice equations, these variants do not improve the robustness to quantization, unlike what is normally expected for FIR models.


INTRODUCTION
In a recent paper [1], a new framework for exploiting data structure in recursive-least-squares (RLS) problems has been introduced. As a result, we have shown how to derive RLS lattice recursions for more general orthonormal networks other than tapped-delay-line implementations [2]. As is well known, the original fast RLS algorithms are obtained by exploiting the shift structure property of the successive rows of the corresponding input data matrix to the adaptive algorithm. That is, consider two successive regression (row) vec- One can exploit this relation to obtain the LS solution in a fast manner. The key for extending this concept to more general structures in [1,3] was to show that, although the above equality no longer holds for general orthonormal models, it is still possible to relate the entries of {u M,N , u M,N+1 } asū where Φ M is an M × (M − 1) structured matrix induced by the underlying orthonormal model. Figure 1 illustrates such structure for which the RLS lattice algorithm of [1] was derived. They constitute what we will refer to in this paper as the a-posteriori-based lattice algorithm, since all these recursions are based on a posteriori estimation errors. Now, it is a well-understood fact that several other equivalent lattice structures exist for RLS filters that result from tapped-delayline models. These alternative implementations are known in the literature as error-feedback, array-based (also referred to as QRD lattice), and normalized lattice algorithms (see, e.g., [4,5,6,7,8]). In [9], all such variants were further extended to the special case of Laguerre-based filters, as we have explained in [1]. Although all these forms are theoretically equivalent, they tend to exhibit different performances when considered under finite-precision effects.
In this paper, we will derive all such equivalent lattice implementations for input data models based on the structure Figure 1: Transversal orthonormal structure for adaptive filtering.
of Figure 1. The use of orthonormal bases can provide several advantages. For example, in some situations, long FIR models can be replaced by shorter compact all-pass-like models as Laguerre filters (see, e.g., [10,11]). From the adaptive filtering point of view, this can represent large savings in computational complexity. The conventional IIR adaptive methods [12,13] present serious problems of stability, local minima, and slow convergence, and in this sense the use of orthogonal bases offers a stable and global solution, due to their fixed poles location. Moreover, orthonormality guarantees good numerical conditioning for the underlying estimation problem, in contrast to other equivalent system descriptions (such as the fixed-denominator model and the partial-fraction representation-see further [2]). The most important application of such structured RLS problems is in the field of line echo cancelation corresponding to long channels, whereby FIR models can be replaced by short orthonormal IIR models. Other applications include channelestimate-based equalization schemes, where the feedforward linear equalizer can be similarly replaced by an orthonormal IIR structure.
After obtaining the new algorithms we will verify the their performance through computer simulations under finite-precision arithmetic. As a result, the new forms turn out to exhibit an unstable behavior. Besides nonminimality of their corresponding algorithm states, they present unstable modes or hyperbolic rotations in their recursions, unlike the corresponding fast variants for FIR models (the latter, in contrast, is free from hyperbolic rotations and unstable modes, and present better numerical properties, despite nonminimality). As a consequence, the new variants do not show improvement in robustness to the quantization effect, compared to the standard RLS lattice recursions of [1], which remains the only reliable extended lattice structure. This discussion on the numerical effects is provided in Section 9.
However, before starting our presentation we call the attention of the reader for an important point. Our main goal in this paper is the development of the equivalent RLS recursions that are normal extensions of the FIR case, and to present some preliminary comparisons based on computer simulations. A complete analytical error analysis for each of these algorithms is not a simple task and is beyond the scope of this paper. Nevertheless, the algorithms derivation is by itself a starting point for further development and improvements on such variants, which is a subject for future research. Moreover, we will provide a brief review and discussion on the minimality and backward consistency properties in order to explain (up to a certain extent) the stability of these variants from the point of view of error propagation. This will be pursued in Section 9, while commenting on the sources of numerical errors in each case.
Notation. In this paper, A ⊕ B is the same as diag{A, B}. We also denote * as the conjugate and transpose of a vector. Since we will be dealing with order-recursive variables, we will write, for example, H M,N , the order-M data matrix up to time N. The same goes for u M,N , e M (N), and so on.

A MODIFIED RLS ALGORITHM
We first provide a brief review of the regularized least-squares problem, but with a slight modification in the definitions of the desired vector, denoted by y N , and the weighting matrix W N . Thus given a column vector y N ∈ C N+1 , and a data matrix H N ∈ C (N+1)×M , the exponentially weighted least-squares problem seeks the column vector w ∈ C M that solves where Π M is a positive regularization matrix, is a weighting matrix defined in terms of a forgetting factor λ satisfying 0 λ < 1, and t is an arbitrary scaling factor. The symbol * denotes complex conjugate transposition. Moreover, we define d N as a growing length vector whose entries are assumed to change according to the following rule: for some scalar θ. 1 Note that the regularized problem in (4) can be conveniently written as where W N = (λ N+L ⊕ λ N+L−1 ⊕ · · · ⊕ t), and where we have factored Π −1 M as for some matrix A M,L . This assumes that the incoming data has started at some point in the past, depending on the number of rows L of A M,L (see [1]). Hence, defining the extended quantities where x i,−1 represents a column of A M,L and h i,N denotes a column of H M,N , as well as we can express (4) as a pure least-squares problem: Therefore, the optimal solution of (11), denoted by w M,N , is given by where We denote the projection of y N onto the range space of Using (5) and the fact that in addition to the matrix inversion formula, it is straightforward to verify that the following (modified) RLS recursions hold: with w M,−1 = 0 M and P M,−1 = Π M . These recursions tell us how to update the weight estimate w M,N in time. The wellknown exponentially weighted RLS algorithm corresponds to the special choice θ = t = 1. The introduction of the scalars {θ, t} allows for a level of generality that is convenient for our purposes in the coming sections.

STANDARD LATTICE RECURSIONS
Algorithm 1 shows the standard extended lattice recursions that solves the RLS problem when the underlying input regression vectors arise from the orthonormal network of Figure 1. The matrix Π M as well as all the initialization variables are obtained according to an offline procedure as described in [1]. The main step in this initialization procedure is the computation of Π M , which remains unchanged for the new recursions we will present in the next sections. The reader should refer to [1] for the details of its computation. Figure 2 illustrates the structure of the mth section of this extended lattice algorithm.

ERROR-FEEDBACK LATTICE FILTERS
Observe that all the reflection coefficients defined for the aposteriori-based lattice algorithm are computed as a ratio in which the numerator and denominator are updated via separate recursions. An error-feedback form is one that replaces the individual recursions for the numerator and denominator quantities by equivalent recursions for the reflection coefficients themselves. In principle, one could derive the recursions algebraically as follows. Consider for instance and some algebra will result in a relation between κ M (N) and Figure 2: A lattice section.
We will not pursue this algebraic procedure here. Instead, we will follow the arguments used in [9] which highlights the interpretation of the reflection coefficients in terms of a leastsquares problem. This will allow us to invoke the recursions we have already established for the modified RLS problem of Section 2 and to arrive at the recursions for the reflection coefficients almost by inspection.

A priori estimation errors
One form of error-feedback algorithm is the one based on a priori, as opposed to a posteriori, estimation errors. They are defined as where now the a posteriori weight vector w Following the same arguments as in Section III of [1], it can be verified that the last entries of these errors satisfy the following order-update relations in terms of the same reflec- The above recursions are well known and they are obtained regardless of data structure. Now, recall that the a-posteriori-based algorithm still requires the recursions is referred to as the a posteriori rescue variable. As we will see in the upcoming sections, similar arguments will also lead to the quantities {β M (N), ν M (N)}, where ν M (N) will be similarly defined as the a priori rescue variable corresponding to v M (N). These in turn will allow us to obtain recursions for their corresponding reflection coefficients {kb M (N), k v M (N)}. Moreover, we will verify that ν M (N) is the actual rescue quantity used in the fixed-order fast transversal algorithm, and which is based on a priori estimation errors.

Exploiting data structure
The procedure to find a recursion for β M,N follows similarly to the one for the a posteriori error b M,N . Thus, beginning from its definition where Φ is the matrix that relates {H M+1,N−1 ,H M,N }, and using the following relations into (22) (see [1]): we obtain, after some algebra, where is the normalized gain vector, defined by the corresponding fast fixed-order recursions. Thus, defining the a priori rescue variable we haveβ In order to obtain a recursion for ν M (N), consider the order-update recursion for k M,N , that is, Taking the complex transposition of (28) and multiplying it from the left by φ M+1 we get Of course, an equivalent recursion for χ M (N) can be obtained, considering the time update for w b M,N , which can be written as Hence, multiplying (30) from the left by φ * M+1 we get Now, it only remains to find recursions for the reflection co-

Time updates for {kb
} We now obtain time relations for the reflection coefficients by exploiting the fact that these coefficients can be regarded as least-squares solutions of order one [9,14].
We begin with the reflection coefficient where, from (31) and Section 5.1 of [1], the numerator and denominator quantities satisfy Now define the angle normalized errors in terms of the square root of the conversion factor γ M (N). It then follows from the above time updates for χ M (N) and can be recognized as the inner products which are written in terms of the following vectors of angle normalized prediction errors: In this way, the defining relation (32) for κb M (N) can be written as which shows that κb M (N) can be interpreted as the solution of a first-order weighted least-squares problem, namely that of projecting the vector (a M+1 W N b M,N ) onto the vector v M+1,N−1 . This simple observation shows that κb M (N) can be readily time updated by invoking the modified RLS recursion introduced in Section 2. That is, by making the identification θ → λ, and t → −1, we have This last equation is obtained from the update forβ M (N) in (27). Similarly, the weight κ v M (N) can be expressed as and therefore, making the identification θ = λ −1 , and t → 1, we can justify the following time update: A similar approach will also lead to the time updates of Algorithm 2 shows the a-priori-based lattice recursions with error feedback. 2

Alternative recursions for the reflection coefficients {κ
} that are based on a posteriori errors can also be obtained. The resulting reflection coefficients updates possess the advantage of avoiding the multiplicative factor λ −1 in the corresponding error-feedback recursions, which represent a potential source of instability of the algorithm.
Thus consider for example the first equality of (38). It can be written as Recalling thatγ M (N) has the updatē we have that 2 Observe that the standard lattice filter obtained in [1] performs feedback of several estimation error quantities from a higher-order problem, for example, b M+1 (N − 1), into the computation of b M (N). The definition of error feedback in fast adaptive filters, however, has been referred to as the feedback of such estimation errors into the computation of the reflection coefficients themselves instead.

Initialization
For m = 0 to M, set µ is a small positive number Algorithm 2: The a-priori-based extended RLS lattice filter with error feedback.
so that we can write (41) as In a similar fashion, we can obtain the following recursion for κ v M (N) from (40): where we have used the fact that We can thus derive similar updates for the other reflection coefficients. Algorithm 3 is the resulting a-posterioribased algorithm.

NORMALIZED EXTENDED RLS LATTICE ALGORITHM
A normalized lattice algorithm is an equivalent variant that replaces each pair of cross-reflection coefficient updates, that is, {κ

Recursion for η M (N)
We start by defining the coefficient along with the normalized prediction errors This yields µ is a small positive number However, from the time-update recursion for ζb M (N) and ζ f M (N), the following relations hold: Substituting these equations into (50), we obtain the desired time-update recursion for the first reflection coefficient: This recursion is in terms of the errors {b M (N), f M (N)}. We thus need to determine order updates for these errors. Thus dividing the order-update equation for b M+1 (N) by Using the order-update relation for ζ b M (N) we also have In addition, the relation, for γ M (N), can be written as Therefore substituting (54) and (56) into (53), we obtain

Recursion for ω M (N)
In a similar vein, we introduce the normalized error and the coefficient Using the order update for ζ 1/2 M (N) and γ M (N), we can establish the following recursion: In order to simplify this equation, we need to relate ζb M (N) to ζ b M+1 (N − 1) andγ M (N) to γ M+1 (N − 1). Recalling the alternative update for ζb M (N): we get where we have defined the reflection coefficient In order to relate {γ M (N), γ M+1 (N − 1)}, we resort to the alternative relation of Algorithm 1: which can be written as Substituting (65) and (68) into (63), we obtain Thus recall that these quantities satisfy the following order updates: which lead to the following relations: Taking the square root on both sides of (72) and substituting into (70), we get (73)

Recursion for ϕ M (N)
for some {m, n, p}. The values of the resulting {m, n, p} can be determined from the equality AΘ M Θ * M A * = AA * = BB * , which gives Similarly, we can implement a J-unitary (hyperbolic) transformation Θb M that lower triangularizes the following prearray of numbers: for some {m , n , p }, which can be determined from the equality CΘb M Θb * M C * = CJC * = DD * . This gives Proceeding similarly we can derive two additional array transformations, all of which are listed in Algorithm 5. The transformations that introduce the zero entries in the postarrays at the desired locations. We illustrate the mth array lattice section that corresponds to this algorithm in Figure 3.

SIMULATION RESULTS
We have performed several simulations in Matlab in order to verify the behavior of the proposed lattice variants under finite-precision arithmetic. Under infinite precision, all lattice filters are equivalent. However, unlike what is normally expected from the corresponding standard FIR lattice filter variants, which tend to be more reliable in finite precision, we observed that the new algorithms exhibit some unstable patterns when compared with the standard lattice recursions of Algorithm 1, in the sense that at some point the algorithms diverge or possibly achieve a much higher MSE.
In the sequel, we will present some simulation results for the algorithms obtained in this paper. We have tested the algorithms over 500 runs, for an exact system identification of a 5-tap orthonormal basis. The noise floor was fixed at -50 dB.
We have observed different behaviors throughout the several scenarios set for testing. In order to characterize their performance up to a certain extent, we have selected some of these settings, which we believe to be the most relevant ones for this purpose.
Experiment 1 (comparison among all algorithms under Matlab precision). Here we have set the forgetting factor λ = 0.99, a typical value, for all the recursions. As a result, we have observed ( Figure 4) the best performance for the a posteriori error-feedback version, followed by the a priori error-feedback recursion exhibiting a slightly higher MSE.
Algorithm 5: The array-based extended RLS lattice filter.
Although the normalized recursion appears to have an even higher MSE and the QR lattice algorithm does not converge, we observed that their behavior changes when λ = 1. In order to observe their behavior more closely, we tested each algorithm separately in fixed-point arithmetic, as shown next. Experiment 2 (a priori error-feedback algorithm for different values of λ). This is shown in Figure 5. In these simulations, we have arbitrarily limited the number of fixed-point quantization steps to 16 bits. Unlike what is observed in experiment 1, this algorithm diverges at some point depending on the value of λ used.  Experiment 3 (a posteriori error-feedback algorithm for different values of λ). This is shown in Figure 6. Similarly to experiment 2, we have kept the quantization steps to 16 bits. Although the scenario of λ = 0.98 seems to be more steady, we still observed divergence when running the experiment over 10 000 data samples. Experiment 4 (performance for 24-bit quantization and for λ = 1). In experiment 1, these algorithms appeared to have a higher MSE compared to the error-feedback versions, and  does not seem to converge to the noise floor. However, we have noticed that in the case where λ = 1, they present some convergence up to a certain instant (depending on how big the number of quantization steps is), and then diverge in two different ways. This is illustrated in Figure 7, for a 24bit fixed-point quantization. For a less number of bits, the performance becomes worse.
In summary, all lattice variants become unstable, except for the standard equations of Algorithm 1, of which the  corresponding MSE curve presents unstable behavior only for λ = 1. Figure 8 illustrates this fact for λ = 1, 0.9, and 0.85, considering a 5-tap orthonormal model.

BACKWARD CONSISTENCY AND MINIMALITY ISSUES
The key structural problem for the unstable behavior observed in all the above algorithms relies on the nonminimality of the state vector in the backward prediction portion of each algorithm. In order to elaborate on this point with some depth, we will briefly review the concepts of minimality and backward consistency for FIR structures ( [15,16,17,18]) and extend these arguments to the algorithms of this paper. Error analysis in fast RLS algorithms is performed in the prediction portion of the recursions, since the flow of the information required for the optimal least-squares solution is one way to the joint process estimation section. The prediction part of the algorithm is a nonlinear mapping of the form where s i denotes the state vector that contains the variables propagated by the underlying algorithm. In finite-precision implementations, however, the actual mapping propagates the perturbed state s i , that is, where δ i is result of quantization. In this case, one's primary goal is to show exponential stability of this system, that is, to show (or not) that the influence of such perturbation should decay sufficiently fast to zero as i → ∞. Thus let S i denote the set of all state values {s i } for which the mapping (86) is exponentially stable. This set includes all state values that can be reached in exact arithmetic, as the input {u(i), u(i − 1), . . .} varies over all realizations that are persistently exciting (we will return to the persistency of excitation issue shortly, considering the general orthonormal basis studied here). Clearly, the state errors i = s i − s i will remain bounded provided that the system (86) is exponentially stable for all states s i and the perturbation δ i does not push s i outside S i . Now, in order to fully understand round-off error effect in a given algorithm, one must consider three aspects in its analysis: (1) error generation, that is, the properties of the round-off error δ i ; (2) error accumulation, that is, how the overall state errors i is affected by the intermediate errors generated at different time instants; (3) error propagation, in the sense that it is assumed that from a certain time instant, no more round-off errors are made, and the propagation of the accumulated errors from that point onward is observed.
Since it is not our intention to embark on a detailed error analysis of these algorithms, we will consider only error propagation stability, which is equivalent to exponential stability of the time recursions. (A conventional stability analysis of an algorithm is usually difficult to pursue, due to the nonlinear nature of system T . This can be accomplished, however, via local linearization and Lyapunov methods, although this requires considerable effort).
A richer approach to the stability problem relies on checking the so-called backward consistency property, that is, if a computed solution (with round-off errors) is the exact solution to a perturbed problem in exact arithmetic. In other words, denote by S i the set of all state vectors that are reachable in finite-precision arithmetic, which vary according to the implementations of the algorithm recursions (that is, the effect of word length plus rounding/truncation). Then, if S i ⊂ S i , the algorithm is said to be backward consistent in this context. Now, an algorithm is said to be nonminimal, if the elements of the state vector s i can be expressed in a reduced dimension with no loss of information, that is, loosely speaking, when the recursions that constitute the algorithm propagate redundancy. In this case, these redundant components can be expressed in terms of constraints of the form f s i = 0, ∀i and u(i) which defines a manifold. In this case, there always exist local perturbations that push s i outside this manifold and therefore out of the stability region S i . It is proved in [15,16,17,18], for fast FIR RLS recursions, that the minimal dimension of s i is 2M + 1, for all i ≥ 2M. In this case, none of the FIR counterparts of the algorithms presented in this paper is minimal, which propagate around 5M variables (this does not mean that stable propagation of the variables cannot be proved without resorting to backward consistency issues, as for example in the case of the a priori error-feedback FIR lattice filter [19]). Now, these minimal components for the FIR case can be established, once the connection with the minimal components of fast transversal filters of all least-squares orders is recognized, thus resulting in a 2M + 1 minimal dimension (see [17]). We have shown in [3], for the general orthonormal basis of this paper, that the same defining components of the minimum state vector in a fast transversal FIR algorithm also define the minimum components of an orthonormalitybased FTF, 3 therefore, using the arguments of [17], one can conclude that this holds similarly for the order-recursive algorithms of this paper, resulting in 2M + 1 minimal parameters. The nonminimal character of the above recursions can be intuitively seen, by following their derivation. It is like solving fast transversal least-squares of all orders, except that the necessity of calculating the augmented Kalman gain vectorǧ M+1,N (in order to propagate g M,N ) on which the FTF is based, is replaced by the use of the augmented scalarb M+1,N (in order to propagate b M (N)), on which the lattice recursions are based.
Clearly, the extended lattice algorithms of this paper are nonminimal; they propagate an additional redundancy that eventually leads to divergence. Consider for instance, the array-based lattice filter obtained. Besides the 5M variables propagated in the forward and backward prediction section, , v M+1 }, which are updated by Θb. Now, note that in the FIR case, even though minimality of these recursions is violated, error propagation in all variables is exponentially stable, since the prearray is scaled by √ λ. The analysis in this case is sufficient for an individual section, since any given lattice section is constituted only by lower-order ones. For the extended array lattice of this paper, however, two facts contribute for divergence. First, variables are further propagated via an unstable mode, that is, λ −1/2 . Second, it makes use of hyperbolic rotations, which are well known to be naturally unstable (unless some care is taken). We recall that in the case of FIR models, only circular, and therefore stable rotations are needed.
In the case of error-feedback algorithms, besides nonminimality, the presence of λ −1 in the recursion for the reflection coefficient κ v M (N) contributes similarly to divergence. This behavior can also be observed in the standard lattice algorithm if one attempts to propagate the minimum cost ζb M (N) via its inverse ζ v M (N), which also depends on λ −1 . The existence of such recursion is also a source of instability in fast fixed-order RLS counterparts. For the normalized lattice algorithm, because all prediction errors are normalized, these recursions become independent of the forgetting factor (except for {ζ 0 (N), ζ b 0 (N)} in the initial step). This, however, eliminates the need for recursion ζb m (N) = σ mγm (N)ζ f m(N − 1), an enforcing relation that represents a source of good numerical conditioning for the algorithm. This relation helps in reducing the redundant variables in fast RLS recursions and is one of the relations forming the manifold of S i in such recursions (as observed for the FTF algorithm in [15,16] in the case of FIR models). 4 This relation has been further extended to the general orthonormal model of [3] and has been used in the standard and it turns out that this is also the case for the lattice recursions, since we have shown that this relation holds for every order-M LS problem.
It is important to add that the fast array algorithm considered in this paper (also called a rotation-based algorithm) is one among a few other fast QR-based recursions. It computes the forward and backward prediction errors in ascending order, leading to the conventional lattice networks studied in this paper. Still, other QR variants are also possible, one of them in fact resulting in a minimal realization for FIR structures [17]; it computes the forward prediction errors in ascending order, but the backward prediction errors in descending order. This is an important case of study, whose extension to the orthonormal basis case will be pursued elsewhere.

Persistency of excitation
One must distinguish between stability due to the algorithm structure and ill-conditioning of the underlying input data matrix. Ill-conditioned data may push the states x i outside the interior of S i , undermining the exponential stability of the recursions. The question then is whether the change of basis would affect the numerical error properties in the above lattice extensions. Now, one of the main purposes of using orthonormal bases lies in the well-studied good numerical conditioning offered by such structures (where input conditioning remains unaltered). Note that this is not the case for an arbitrary model realization as a fixed denominator, of partialfraction representation as pointed out in [2]. Consider the vector of basis functions associated with the orthonormal model: Under the assumption of stationarity of the input sequence u(n), the input correlation matrix associated with these basis functions is given by B e jω B * e jω S u e jω dω, where S u (e jω ) is the power spectral density of u(n). Now the conditioning of the (linear) estimation problem associated with the orthonormal model is related to the condition number of R B , namely, k(R B ), which shows that orthonormality of the basis functions basically leaves the condition number unaltered when passing through such model. Therefore we can say that for the same condition number, the extended lattice recursions behave worse than their tapped-delay-line lattice counterparts.

CONCLUSION
In this work, we developed several lattice forms for RLS adaptive filters based on general orthonormal realizations. One technique is based on propagating the reflection coefficients in time. A second form is based on propagating a fewer number of normalized reflection coefficients. A third form is based on propagating angle-normalized quantities via unitary and J-unitary rotations. Even though the algorithms are all theoretically equivalent, they differ in computational cost and in robustness to finite-precision effects. The new algorithms, besides nonminimality, present unstable modes as well as hyperbolic rotations, so that the well-known good numerical properties observed in the class of FIR models no longer exist for the extended fast recursions derived. In this context, the standard lattice recursions of Algorithm 1 represent up to now the most numerically reliable algorithm for this class of input data structures. We remark that the development of the above recursions represent an initial step towards further refinement of these algorithms and it is not our purpose to provide all answers on these extended lattice filter variants in this presentation. Although our presentation lacks a more precise analysis of the numerical behavior of these algorithms, we believe that the arguments in Section 9 suffice as a preliminary explanation for the unstable behavior observed in all lattice variants. As future works, we will look into the numerical issues of these algorithms in detail as well as pursue a minimal fast QR realization similarly to [17] for FIR models.