AN APPROACH TO STABILIZING THE FAST ARRAY RLS ADAPTIVE FILTER USING HOMOGENEOUS COORDINATES IN PROJECTIVE GEOMETRY Todd K. Moon∗ and Kevin Hencke† and Jacob H. Gunther‡
ABSTRACT The fast array RLS (FARLS) adaptive filter employs a circular and hyperbolic rotation to transform from a pre-array to a post-array. The algorithm is known to be numerically unstable for finite-precision implementations. It is believed that the difficulty lies in the numerical computation of the hyperbolic rotation. We propose the use of homogeneous coordinates in projective geometry to transform the hyperbolic rotation to a circular rotation. This rotation (being circular) is numerically stable. However, the transformation to the new circular representation in the projective coordinates is still numerically difficult. Nevertheless, this method may lead to a new approach that completely circumvents the instability. 1. INTRODUCTION Fast Array RLS (FARLS) adaptive filter algorithms have complexity linear in the filter length [1], but it is known that for finite-precision (even double-precision) implementations operating on real-world data they may go unstable. This instability introduces additional complexity and awkwardness into the adaptive processing, which must be monitored for instability onset and restarted. Two key processing steps of the FARLS algorithm are a circular rotation and a hyperbolic rotation which together drive two locations of the “pre-array” to 0. In these rotations, a point on a conic locus (circle or ellipse) is rotated along that locus until one of the point coordinates is equal to 0. It is believed that the source of the numeric instability is the hyperbolic rotation. A point that is too close to the asymptotes — too far out along the hyperbola — cannot with numerical reliability be moved along the hyperbola to its zero crossing point. This paper introduces a novel attempt to stabilize the adaptive processing algorithm by transforming the hyperbolic rotation to a circular rotation by means of projective Information Dynamics Laboratory and Electrical and Computer Engineering Dept., Utah State University, Logan, UT 84332,
[email protected] University of Maryland,
[email protected] Utah State University,
[email protected] This work was supported by the NSA REU in Coding and Communications, Summer 2009
coordinates. In general terms, it has been hoped that the extra degrees of freedom provided by projective coordinates might enable the algorithm to be stabilized. More specifically, in projective coordinates, all nonsingular conic sections are equivalent: a circle is an ellipse is is a hyperbola. Expressing the hyperbolic rotation in projective coordinates converts it into a circular rotation. This conversion actually involves two conceptual steps: expression as an ellipse, then scaling the ellipse to a circle. The results of this attempt are mixed. We have successfully made the transformation to projective coordinates and expressed the hyperbola as a circle, employing this in a modified FARLS algorithm. However, the scaling of the ellipse to a circle still experiences numerical difficulties when the ellipse is extremely eccentric. Geometrically this corresponds to precisely the same condition that the hyperbola approaches its asymptotes. The result is that (for somewhat different reasons, and with more insight gained), the projective approach still fails to reliably stabilize the FARLS algorithm. Nevertheless, the method described here introduces new perspective to the problem, and may provide a valuable stepping stone to a fully stabilized algorithm. 2. SUMMARY OF THE FARLS ALGORITHM AND ITS STABILIZATION Due to space limitations, we forgo steps of the development (see [1] for a lucid development — our presentation is transposed from that — see also [2, 3, 4]) and simply state that the FARLS algorithm begins with the formation of a prearray matrix A whose components (with dimensions indicated, where M is the length of the filter) are · ¸ A1×1 D1×(M +1) . A= B2×1 E2×(M +1) A matrix Θ is sought that transforms the pre-array A into a post-array matrix B of the form · ¸ X Y B = ΘA = . (1) 02×1 Z The matrix Θ must be 1 ⊕ S-unitary, that is, must satisfy the equation · ¸ · ¸ 1 H 1 Θ Θ= . S S
£ 0 ¤ where S = 10 −1 . Let J = 1 ⊕ S. The significant aspect of the transformation (1) is that it places zeros in the last two elements of the first column of B. Let us denote the first column of A (before the trans£ ¤T formation) as v = a b c , and call this the “priming vector.” Furthermore, let us factor Θ as Θ = Ω2 Ω1 , where Ω1 rotates circularly on the first two coordinates, and Ω2 rotates hyperbolically (due to the signature of S) on the first and third coordinates. We can write θ1,11 θ1,12 0 θ2,11 0 θ2,12 1 0 . Ω1 = θ1,21 θ1,22 0 Ω2 = 0 0 0 1 θ2,21 0 θ2,22 Then θ1 =
· θ1,11 θ1,21
θ1,12 θ1,22
¸ and
θ2 =
· θ2,11 θ2,21
θ2,12 θ2,22
¸
are unitary and S-unitary rotation matrices, respectively: θ1H θ1 = I and θ2H Sθ2 = S. That is, θ1 produces a circular rotation and θ2 produces a hyperbolic rotation. The circular rotation θ1 is such that · ¸ · ¸ d a a d θ1 = or, equivalently, Ω1 b = 0 , b 0 c c with a2 + b2 = d2 . This can be stably computed using (for example) a Givens rotation. This is followed by a hyperbolic rotation on the vector v2 = [ dc ], · ¸ · ¸ d e θ2 = θ2 v2 = c 0 such that d2 − c2 = e2 . The point (d, c) thus lies on the hyperbola x2 − y 2 = e2 , and the hyperbolic rotation amounts to moving along this hyperbola to the axis-crossing point (e, 0). Computing the hyperbolic rotation is known to be numerically unstable in finite precision arithmetic [5, 6, 7, 8, 9]. This has led to the development of rescue schemes [10, 11, 12] which restart the algorithm whenever they detect instability. Various methods of computing the hyperbolic rotation have been presented, including the orthogonal-diagonal method [13] which presents the rotation θ2 as an eigendecomposition q · ¸ · ¸ c+d 0 1 1 1 c−d q 1 −1 √1 . θ2 = √ d−c 1 1 2 −1 1 2 0 c+d (2) This is claimed to be numerically stable [1, p. 866], but experiments performed at our laboratory in float precision in C on speech data still result in instabilities.
Fig. 1. A hyperbola placed into R3 at z = 1. R3 . By representing R2 as a plane in R3 , we can generalize equations of the two-dimensional conic sections to re-create cones in R3 . We do this by the use of what are known as homogeneous coordinates, which embed an n-dimensional field Rn in an n + 1-dimensional field by appending a 1 to the end of every coordinate vector from the n-dimensional field, creating the space denoted as RP n , which may be viewed as Rn+1 . Homogeneous coordinates represent projective geometry; see, for example, [14]. In our case, for any vector x = (x, y) ∈ R2 , the corresponding point in homogeneous space is xh = (x, y, 1) ∈ RP 2 . This can be visualized as placing R2 into R3 at the plane z = 1. Figure 1 demonstrates this with a hyperpola on the plane. Also indicated on the hyperbola are the points (d, c, 1) and the desired rotation point (e, 0, 1), corresponding to the points on the original hyperbola (d, c) and (e, 0). In homogeneous coordinates, vectors equivalent up to a nonzero multiple are identified as the same point. Explicitly, if h1 = λh2 for λ 6= 0 then h1 ≡ h2 . Any two equivalent vectors in RP 2 represent the same single vector in R2 . This means that points in R2 become lines in R3 (representing RP 2 ). Rotation of points along a hyperbola in R2 can be generalized to motion on the cone in R3 of which the hyperbola is a section. This motion may be achieved, specifically, by a rotation of points about any ellipse found by cutting the cone by a plane perpendicular to its axis (excepting the trivial case of cutting the cone at its single-point tip, since the origin is not meaningful in homogeneous coordinates). What this means is that we can take our hyperbola x2 − 2 y = e2 , represent it as a cone in R3 , and do rotations on an elliptic cross-section of that cone, rather than having to do rotations on the hyperbola. The equation for the hyperbola x2 − y 2 = e2 when embedded into R3 represents the section of a cone whose axis is the x axis. Such a cone is determined in general by the equation (k1 x)2 = (k2 y)2 + (k3 z)2 for ki ∈ R
3. ROTATION IN PROJECTIVE SPACE It is well known that lines, ellipses, parabolas, and hyperbolas in R2 represent the intersection of a cone with a plane in
or (k1 x)2 − (k2 y)2 = (k3 z)2 .
Fig. 2. The hyperbola of figure 1 are generalized to a cone in projective space Thus, as we know z to be 1 for our hyperbolic cross-section of the cone, this general cone must specialize to x2 − y 2 = (ez)2 . This expresses the curve in homogeneous coordinates (there is the same degree for each term). The cone can thus be written as x2 = y 2 + (ez)2
(3)
Cross-sections of this cone taken along planes x = x0 have the form (x0 )2 = y 2 + (ez)2 or e 1 2 y) + ( z)2 , (4) x0 x0 so that the cone is not a circular cone, but rather elliptical, having elliptical cross sections. When we move into homogeneous coordinates, the addition of a new coordinate and the scaling equivalence means that all points in R2 become lines in homogeneous space, and any point along this line (except the origin) can be used to represent the original point in R2 . We can thus reverse the homogenization without loss of information: any 2-dimensional cross-section of the variety (provided that this cross-section does not pass through the origin) yields the same information geometrically as the original hyperbola. By this process, we have generalized a hyperbola into a cone, and we can then turn the cone into an ellipse by cutting along the plane x = d. This ellipse preserves all the data of the original hyperbola: The overall concept is portrayed in figure 3. The original hyperbola (blue) has two points on it (red): (d, c) and (e, 0). In homogeneous coordinates, the points become lines (red). The elliptical cross section of the cone at x = d is shown (green). There is a point-for-point equivalence between points on the hyperbola and points on the ellipse. Movement on the ellipse is equivalent to movement on the hyperbola. The intersection of the line representing (e, 0) and the ellipse represents to point to rotate to on the ellipse. The elliptical cone (4) can be converted to a circular cone by scaling each point on the ellipse by the coordinate change matrix 1 e or C = e 1 C= (5) 1 1 e 1=(
Fig. 3. The original hyperbola at z = 1 with the elliptical cross section of the cone from the cross section at x = d. Provided that e 6= 0, these matrices are interchangeable since eCx = Cex = C(ex) ≡ Cx in homogeneous coordinates. Applying C to a point in R3 redefines the axes to turn (4) into 1 1 1 = ( y)2 + ( z)2 . (6) x0 x0 This is a circle lying on the cross-sectioning plane at x = x0 , upon which a conventional circular rotation may be applied. Unfortunately, the transformation C requires knowing e2 = d2 − c2 , which may be difficult to compute. We use the rotation in (2) to determine e (but not to apply to the rest of the pre-array). This minimizes the numerical difficulty. From v2 above, we define the homogenized coordinate point · ¸ d v2 h2 = = c 1 1 and apply C to transform R3 . Let d Ch2 = k2 = c . 1 e
Now k2 lies on the surface of the circular cone. To zero the second component of k2 , we need only rotate it around the x-axis. This can be accomplished by a rotation matrix 1 θ3,11 θ3,12 (7) θ3,21 θ3,22 where θ3 is circular rotation matrix which zeros the the first component of the vector it is applied to. Putting all the pieces together, we first scale to circular, rotate, then scale back, producing the transformation on the
6. REFERENCES [1] A. H. Sayed, Fundamentals of Adaptive Filtering. Hoboken, NJ: Wiley Interscience, 2003. [2] A. Sayed and T. Kailath, “Structured matrices and fast RLS adaptive filtering,” in Prosc. 2nd. IFAC Workshop on Algorithms and Architectures for Real-Time Control, (Seoul, Korea), pp. 211–216, 1992. [3] A. Sayed and T. Kailath, “Extended Chandresekahr recursions,” IEEE Trans. Automatic Control, vol. 39, no. 3, pp. 619–623, 1994.
Fig. 4. Cone with moderate eccentricity. homogenized vector h2 as 1 d e e Θ3,1,1 Θ3,1,2 C c = λ 0 ≡ 0 Ω3 h2 = C −1 Θ3,2,1 Θ3,2,2 1 1 1 (8) The point (e, 0, 1) so produced lies on the z = 1 plane, and represents the intersection of the original hyperbola with the x axis. To apply this now the FARLS algorithm, the circularrotation transformation Ω1 is first applied across the columns of the pre-array A to obtain a mid-array A0 . Then the transformation Ω3 is applied to the homogenized coordinates obtained from the first and third rows of A0 . These are then placed back to create the post array B. 4. POTENTIAL NUMERICAL DIFFICULTIES While the rotation is now circular in the projective coordinates, the transformation Ω3 in (8) requires division by e in either C or C −1 . As the value d2 −c2 = e2 nears zero, small perturbations in d and c give rise to dramatically different values of e. Thus, while mathematically e 6= 0, numerically e may be 0 or close enough to 0 that division by e breaks the update. This can be appreciated geometrically from the cone plotted in figure 4. As e becomes smaller, the cone becomes increasingly eccentric. In this limit, the cone become planar and rotation on the cone is impossible. 5. SUMMARY AND CONCLUSION This paper describes a new approach to the hyperbolic rotation in the update rule of the FARLS algorithm. This method offers the hope of stabilization since all rotations are circular (bounded), but unfortunately still requires computation of e2 = d2 − c2 and its reciprocal. Further work along this line, making additional use of projective coordinates to produce an algorithm which rides out the unstable periods, is necessary.
[4] A. Sayed and T. Kailath, DSP Handbook, ch. Recursive leastsquares adaptive filters. CRC Press, 1998. [5] P. A. Regalia, “Numerical stability issues in fast least-squares adaptation algorithms,” Optical Engineering, vol. 41, pp. 1144–1152, 1992. [6] P. A. Regalia, “Numerical stability properties of a QR-based fast lesat squares algorithm,” Optical Engineering, vol. 41, pp. 2006–2109, 1993. [7] D. Slock, “The backward consistency concept and roundoff error propagation dynamics in RLS algorithms,” Optical Engineering, vol. 31, pp. 1153–1169, 1992. [8] J. Bunch, R. C. Leborne, and I. Proudler, “A conceptual framework for consistency, conditioning, and stability issues in signal processing,” IEEE Trans. Sig. Proc., vol. 49, no. 9, pp. 1971–1981, 2001. [9] J. Cioffi, “Limited-precision effects in adaptive filtering,” IEEE Trans. Circ. Syst., vol. 34, no. 7, pp. 821–833, 1987. [10] D. Lin, “On the digital implementation of the fast Kalman algorithm,” IEEE Trans. Acoust., Speech, and Sig. Proc., vol. 32, pp. 998–105, 1984. [11] E. Eleftheriou and D. Falconer, “Tracking properties and steady-state performance of RLS adaptive filters,” IEEE Trans. Acoust., Speech and Sig. Proc., vol. 34, pp. 1097– 1110, 1986. [12] J. Cioffi and T. Kailath, “Fast recursive least-squares transversal filters for adaptive filtering,” IEEE Trans. Acoust., Speech, and Sig. Proc., vol. 32, pp. 304–337, 1984. [13] S. Chandrasekaran and A. Sayed, “Stabilizing the generalized Schur algorithm,” SIAM J. Matrix Anal. Appl., vol. 17, pp. 152–216, 1996. [14] D. Cox, J. Little, and D. O’Shea, Using Algebraic Geometry. New York: Springer-Verlag, 1998.
