Least Squares with Examples in Signal Processing1 x

1 Least Squares with Examples in Signal Processing1 The sum of Squares ' of x is denoted by kxk22 , , Ivan Selesnick March 7, 2013 NYU-Poly kxk22 =. X. |x(n)|2 = xT x. (3). These notes address (approximate) solutions to linear equations by Least n Squares . We deal with the easy' case wherein the system matrix is full The energy' of a vector x refers to kxk22 . rank. If the system matrix is rank deficient, then other methods are In these notes, it is assumed that all vectors and matrices are real-valued. needed, , QR decomposition, singular value decomposition, or the In the complex-valued case, the conjugate transpose should be used in place pseudo-inverse, [2, 3]. of the transpose, etc. In these notes, Least Squares is illustrated by applying it to several basic problems in Signal processing: 2 Overdetermined equations 1. Linear prediction Consider the system of linear equations 2.

2 Smoothing y = Hx. 3. Deconvolution If there is no solution to this system of equations, then the system is 4. System identification overdetermined'. This frequently happens when H is a tall' matrix (more rows than columns) with linearly independent columns. 5. Estimating missing data In this case, it is common to seek a solution x minimizing the energy of For the use of Least Squares in filter design, see [1]. the error: 1 Notation J(x) = ky Hxk22 . We denote vectors in lower-case bold, , Expanding J(x) gives . x1. J(x) = (y Hx)T (y Hx) (4). x2 .. x= .. (1) = yT y yT Hx xT HT y + xT HT Hx (5).. T T T T. xN = y y 2y Hx + x H Hx. (6). We denote matrices in upper-case bold. The transpose of a vector or matrix Note that each of the four terms in (5) are scalars. Note also that the scalar in indicated by a superscript T , , xT is the transpose of x.

3 XT HT y is the transpose of the scalar yT Hx, and hence xT HT y = yT Hx. The notation kxk2 refers to the Euclidian length of the vector x, , Taking the derivative (see Appendix A), gives . q kxk2 = |x1 |2 + |x2 |2 + + |xN |2 . (2) J(x) = 2HT y + 2HT Hx x 1 For feedback/corrections, email Setting the derivative to zero, Matlab software to reproduce the Examples in these notes is available on the web or from the author.. Last edit: November 26, 2013 J(x) = 0 = HT Hx = HT y x 1. Let us assume that HT H is invertible. Then the solution is given by Set the derivatives to zero to get: x = (HT H) 1 HT y. x=. 1 T. H (11). 2. This is the Least Squares ' solution. y = Hx (12). 1. min ky Hxk22 = T. x = (H H) T. H y (7). x Plugging (11) into (12) gives In some situations, it is desirable to minimize the weighted square error, 1. y= HHT . , n wn rn2 where r is the residual, or error, r = y Hx, and wn are P.

4 2. positive weights. This corresponds to minimizing kW1/2 (y Hx)k22 where Let us assume HHT is invertible. Then W is the diagonal matrix, [W]n,n = wn . Using (7) gives = 2(HHT ) 1 y. (13). min kW1/2 (y Hx)k22 = x = (HT WH) 1 HT Wy (8). x where we have used the fact that W is symmetric. Plugging (13) into (11) gives the Least Squares ' solution: 3 Underdetermined equations x = HT (HHT ) 1 y. Consider the system of linear equations We can verify that x in this formula does in fact satisfy y = Hx by plugging y = Hx. in: Hx = H HT (HHT ) 1 y = (HHT )(HHT ) 1 y = y . If there are many solutions, then the system is underdetermined'. This X. frequently happens when H is a wide' matrix (more columns than rows). with linearly independent rows. So, In this case, it is common to seek a solution x with minimum norm. That min kxk22 y = Hx = x = HT (HHT ) 1 y.

5 (14). is, we would like to solve the optimization problem x min kxk22 (9) In some situations, it is desirable to minimize the weighted energy, , x P 2. n wn xn , where wn are positive weights. This corresponds to minimizing such that y = Hx. (10). kW1/2 xk22 where W is the diagonal matrix, [W]n,n = wn . The derivation Minimization with constraints can be done with Lagrange multipliers. So, of the solution is similar, and gives define the Lagrangian: 1. min kW1/2 xk22 y = Hx = x = W 1 HT HW 1 HT y L(x, ) = kxk22 T. + (y Hx) x (15). Take the derivatives of the Lagrangian: This solution is also derived below, see (25). L(x) = 2x HT . x 4 Regularization In the overdetermined case, we minimized ky Hxk22 . In the underde- L(x) = y Hx termined case, we minimized kxk22 . Another approach is to minimize the 2. weighted sum: c1 ky Hxk22 + c2 kxk22.

6 The solution x depends on the Setting the derivative to zero, ratio c2 /c1 , not on c1 and c2 individually.. A common approach to obtain an inexact solution to a linear system is J(x) = 0 = HT Hx + AT Ax = HT y x to minimize the objective function: = (HT H + AT A) x = HT y J(x) = ky Hxk22 + kxk22 (16). So the solution is given by where > 0. Taking the derivative, we get x = (HT H + AT A) 1 HT y . J(x) = 2HT (Hx y) + 2 x x So, Setting the derivative to zero, min ky Hxk22 + kAxk22. x (18). J(x) = 0 = HT Hx + x = HT y = x = (HT H + AT A) 1 HT y x = (HT H + I) x = HT y Note that if A is the identity matrix, then equation (18) becomes (17). So the solution is given by 6 Constrained Least Squares T 1 T. x = (H H + I) H y Constrained Least Squares refers to the problem of finding a Least Squares solution that exactly satisfies additional constraints.

7 If the additional So, constraints are a set of linear equations, then the solution is obtained as min ky Hxk22 + kxk22 = x = (HT H + I) 1 HT y (17) follows. x The constrained Least Squares problem is of the form: This is referred to as diagonal loading' because a constant, , is added to the diagonal elements of HT H. The approach also avoids the problem min ky Hxk22 (19). x of rank deficiency because HT H + I is invertible even if HT H is not. In such that Cx = b (20). addition, the solution (17) can be used in both cases: when H is tall and when H is wide. Define the Lagrangian, 5 Weighted regularization L(x, ) = ky Hxk22 + T (Cx b). A more general form of the regularized objective function (16) is: The derivatives are: J(x) = ky Hxk22 + kAxk22.. L(x) = 2HT (Hx y) + CT . where > 0. Taking the derivative, we get x . J(x) = 2HT (Hx y) + 2 AT Ax L(x) = Cx b x.

8 3. Setting the derivatives to zero, Special cases Simpler forms of (23) are frequently useful. For example , if H = I and L(x) = 0 = x = (HT H) 1 (HT y ) (21). x b = 0 in (23), then we get min ky xk22 Cx = 0. L(x) = 0 = Cx = b (22) x 1 (24). = x = y CT CCT Cy Multiplying (21) on the left by C gives Cx, which from (22) is b, so we have If y = 0 in (23), then we get C(HT H) 1 (HT y ) = b min kHxk22 Cx = b x 1 (25). or, expanding, = x = (HT H) 1 CT C(HT H) 1 CT b C(HT H) 1 HT y (HT H) 1 CT = b. If y = 0 and H = I in (23), then we get Solving for gives 1. min kxk22 Cx = b = x = CT CCT b (26). x T 1 T 1 T 1 T.. = 2 C(H H) C C(H H) H y b which is the same as (14). Plugging into (21) gives 7 Note 1 . x = (HT H) 1 HT y CT C(HT H) 1 CT C(HT H) 1 HT y b The expressions above involve matrix inverses. For example , (7) involves Let us verify that x in this formula does in fact satisfy Cx = b, (HT H) 1.

9 However, it must be emphasized that finding the Least square solution does not require computing the inverse of HT H even though the Cx = C(HT H) 1 HT y inverse appears in the formula. Instead, x in (7) should be obtained, in CT C(HT H) 1 CT. 1. C(HT H) 1 HT y b practice, by solving the system Ax = b where A = HT H and b = HT y. The most direct way to solve a linear system of equations is by Gaussian = C(HT H) 1 HT y elimination. Gaussian elimination is much faster than computing the inverse 1. C(HT H) 1 CT C(HT H) 1 CT C(HT H) 1 HT y b . of the matrix A. = C(HT H) 1 HT y C(HT H) 1 HT y b . 8 Examples =b X. Polynomial approximation So, An important example of Least Squares is fitting a low-order polynomial to data. Suppose the N -point data is of the form (ti , yi ) for 1 i N . The min ky Hxk22 Cx = b = goal is to find a polynomial that approximates the data by minimizing the x 1 energy of the residual.

10 X = (HT H) 1 HT y CT C(HT H) 1 CT C(HT H) 1 HT y b X. E= (yi p(ti ))2. (23). i 4. where p(t) is a polynomial, , p(t) = a0 + a1 t + a2 t2 . The problem can be viewed as solving the overdetermined system of equa- tions, 2.. Data y1 1 t1 t21 . y2 1 t2 t22 a0 .. 1.. a1 , .. 0 .. a 2. yN 1 tN t2N. 1. which we denote as y Ha. The energy of the residual, E, is written as 2. 0 1 2. E = ky Hak22 . 2. Polynomial approximation (degree = 2) From (7), the Least Squares solution is given by a = (HT H) 1 HT y. An 1 example is illustrated in Fig. 1. 0. Linear prediction 1. One approach to predict future values of a time-series is based on linear prediction, , 2. 0 1 2. y(n) a1 y(n 1) + a2 y(n 2) + a3 y(n 3). (27). 2. Polynomial approximation (degree = 4). If past data y(n) is available, then the problem of finding ai can be solved 1. using Least Squares .