Lecture 10: Recursive Least Squares Estimation

1 1. Lecture 10: Recursive Least Squares Estimation Overview Recursive Least Squares Estimation ;. The exponentially weighted Least Squares Recursive -in-time solution Initialization of the algorithm Recursion for MSE criterion Examples: Noise canceller, Channel equalization, Echo cancellation Reference : Chapter 9 from S. Haykin- Adaptive Filtering Theory - Prentice Hall, 2002. Lecture 10 2. Recursive Least Squares Estimation Problem statement Given the set of input samples {u(1), u(2), .. , u(N )} and the set of desired response {d(1), d(2).}

2 , d(N )}. In the family of linear filters computing their output according to M. X. y(n) = wk u(n k), n = 0, 1, 2, .. (1). k=0. Find recursively in time the parameters {w0 (n), w1 (n), .. , wM 1 (n)} such as to minimize the sum of error Squares n X n X M. X 1. E(n) = E(w0 (n), w1 (n), .. , wM 1 (n)) = (n, i)[e(i)2 ] = (n, i)[d(i) wk (n)u(i k)]2. i=i1 i=i1 k=0. where the error signal is M. X 1. e(i) = d(i) y(i) = d(i) wk (n)u(i k). k=0. and the forgetting factor or weighting factor reduces the influence of old data 0 < (n, i) 1, i = 1, 2.

3 , n usually taking the form (0 < < 1). (n, i) = n i , i = 1, 2, .. , n 2. Lecture 10 3. The exponentially weighted Least Squares solution Writing the criterion with an exponential forgetting factor n X n X M. X 1. E(n) = E(w0 (n), w1 (n), .. , wM 1 (n)) = n i [e(i)2 ] = n i [d(i) wk (n)u(i k)]2. i=i1 i=i1 k=0. Make the following variable changes: . u0 (i) = n i u(i); d0 (i) = n i d(i) (2). Then the criterion rewrites n X M. X 1 n X M. X 1. E(n) = n i [d(i) wk (n)u(i k)]2 = [d0 (i) wk (n)u0 (i k)]2. i=i1 k=0 i=i1 k=0.

4 Which is the standard LS criterion, in the new variables u0 (i), d0 (i). The LS solution can be obtained as n X n n n 0 0 T 1 X 0 0 X n i T 1 X. w(n) = ( u (i)u (i) ) u (i)d (i) = ( u(i)u(i) ) n i u(i)d(i) = [ (n)] 1 (n) =. i=i1 i=i1 i=i1 i=i1. where we will denote (making use of Pre-windowing assumption, that data before i = 1 is zero). n X. (n) = n i u(i)u(i)T. i=1. n X. (n) = n i u(i)d(i). i=1. Lecture 10 4. Recursive in time solution We want to find a Recursive in time way to compute w(n) = [ (n)] 1 (n).

5 Using the information already available at time n 1, w(n 1) = [ (n 1)] 1 (n 1). We therefore will rewrite the variables (n) and (n) as functions of (n 1) and (n 1). n X n 1. X. (n) = n i u(i)u(i)T = n 1 i u(i)u(i)T + u(n)u(n)T = (n 1) + u(n)u(n)T. i=1 i=1. Xn n 1. X n 1 i (n) = n i u(i)d(i) = u(i)d(i) + u(n)d(n) = (n 1) + u(n)d(n). i=1 i=1. The matrix inversion formula If A and B are M M positive definite matrices, D is a N N matrix, and C is a M N matrix which are related by A = B 1 + CD 1 C T. then A 1 = B BC(D + C T BC) 1 C T B.

6 Proof Exercise. Lecture 10 5. Derivation of the algorithm Applying the matrix inversion formula to (n) = (n 1) + u(n)u(n)T. we obtain 1 1 1 2 1 (n 1)u(n)uT (n) 1 (n 1). (n) = (n 1) . 1 + 1 uT (n) 1 (n 1)u(n). Denoting P (n) = 1 (n). and 1 P (n 1)u(n). k(n) = = Exercise = P (n)u(n). 1 + 1 uT (n)P (n 1)u(n). we obtain P (n) = 1 P (n 1) 1 k(n)uT (n)P (n 1). We are now able to derive the main time-update equation, that of w(n). w(n) = [ (n)] 1 (n) = P (n) (n) = P (n)( (n 1) + u(n)d(n)) = P (n)( (n 1)w(n 1) + u(n)d(n)).

7 = P (n)(( (n) u(n)u(n)T )w(n 1) + u(n)d(n)) = w(n 1) P (n)u(n)u(n)T w(n 1) + P (n)u(n)d(n). = w(n 1) + P (n)u(n)(d(n) u(n)T w(n 1)) = w(n 1) + P (n)u(n) (n) = w(n 1) + k(n) (n). Lecture 10 6. where (n) = d(n) u(n)T w(n 1). is the innovation process (apriori errors). Now we can collect all necessary equations to form the RLS. algorithm: 1 P (n 1)u(n). k(n) =. 1 + 1 uT (n)P (n 1)u(n). (n) = d(n) u(n)T w(n 1). w(n) = w(n 1) + k(n) (n). P (n) = 1 P (n 1) 1 k(n)uT (n)P (n 1). Lecture 10 7. Initialization of RLS algorithm In RLS algorithm there are two variables involved in the recursions (those with time index n 1): w(n 1), Pn 1.

8 We must provide initial values for these variables in order to start the recursions : w(0). If we have some apriori information about the parameters w this information will be used to initialize the algorithm. Otherwise, the typical initialization is w(0) = 0. P0. 1. Recalling the significance of P (n). n X. P (n) = 1 (n) = [ n i u(i)u(i)T ] 1. i=i1. the exact initialization of the recursions uses a small initial segment of the data u(i1 ), u(i1 + 1), .. , u(0). to compute 0. X. P (0) = 1 (0) = [ i u(i)u(i)T ] 1.

9 I=i1. However, it is not a simple matter to select the length of data required for ensuring invertibility of (0) ! Lecture 10 8. 2. The approximate initialization is commonly used, it doesn't require matrix inversion: P (0) = I. There is an intuitive explanation of this initialization. The significance P (n) = 1 (n) const. E(w(n) w)(w(n). T. w). can be proven. Thus, P (n) is proportional to the covariance matrix of the parameters w(n). Since our knowledge of these parameters at n = 0 is very vague, a very high covariance matrix of the parameters is to be expected, and thus we must assign a high value to.

10 The recommended value for is > 100 u2. For large data length, the initial values assigned at n = 0 are not important, since they are forgotten due to exponential forgetting factor . Lecture 10 9. Summary of the RLS algorithm Given data u(1), u(2), u(3), .. , u(N ) and d(1), d(2), d(3), .. , d(N ). 1. Initialize w(0) = 0, P0 = I. 2. For each time instant, n = 1, .. , N , Compute = uT (n)P (n 1) M 2 f lops = + u M f lops T. k(n) = M f lops (n) = d(n) wT (n 1)u(n) M f lops w(n) = w(n 1) + k(n) (n) M f lops 0. P = k(n) M 2 f lops P (n) = 1 (P (n 1) P 0 ) M 2 f lops We used flop as an abbreviation for one addition (subtraction) + one multiplication (floating point operations).