<< Chapter < Page Chapter >> Page >

where x ^ is the new weighted least squares solution of [link] which is used to partially update the previous value x ( k - 1 ) using a convergence up-date factor 0 < q < 1 which gave convergence over a larger range of around 1 . 5 < p < 5 but but it was slower.

A second improvement showed that a specific up-date factor of

q = 1 p - 1

significantly increased the speed of convergence. With this particular factor, the algorithm becomes a form of Newton's method which has quadratic convergence.

A third modification applied homotopy [link] , [link] , [link] , [link] by starting with a value for p which is equal to 2 and increasing it each iteration (or each few iterations) until it reached the desired value, or, in the case of p < 2 , decrease it. This made a significant increase in both therange of p that allowed convergence and in the speed of calculations. Some of the history and details can be found applied to digital filter design in [link] , [link] .

A Matlab program that implements these ideas applied to our pseudoinverse problem with more equations than unknowns (case 2a) is:

% m-file IRLS1.m to find the optimal solution to Ax=b %  minimizing the L_p norm ||Ax-b||_p, using IRLS. %  Newton iterative update of solution, x, for  M > N. %  For 2<p<infty, use homotopy parameter K = 1.01 to 2 %  For 0<p<2, use K = approx 0.7 - 0.9 %  csb 10/20/2012 function x = IRLS1(A,b,p,K,KK) if nargin < 5, KK=10;  end; if nargin < 4, K = 2;  end; if nargin < 3, p = 10; end; pk = 2;                                      % Initial homotopy value x  = pinv(A)*b;                              % Initial L_2 solution E = []; for k = 1:KK                                 % Iterate    if p >= 2, pk = min([p, K*pk]);           % Homotopy change of p       else pk = max([p, K*pk]); end    e  = A*x - b;                             % Error vector    w  = abs(e).^((pk-2)/2);                  % Error weights for IRLS    W  = diag(w/sum(w));                      % Normalize weight matrix    WA = W*A;                                 % apply weights    x1  = (WA'*WA)\(WA'*W)*b;                 % weighted L_2 sol.    q  = 1/(pk-1);                            % Newton's parameter    if p > 2, x = q*x1 + (1-q)*x; nn=p;       % partial update for p>2       else x = x1; nn=2; end                 % no partial update for p<2    ee = norm(e,nn);   E = [E ee];            % Error at each iteration end plot(E)

This can be modified to use different p 's in different bands of equations or to use weighting only when the error exceeds a certain threshold to achieve aconstrained LS approximation [link] , [link] , [link] . Our work was originally done in the context of filter design but others have done similar things insparsity analysis [link] , [link] , [link] .

This is presented as applied to the overdetermined system (Case 2a and 2b) but can also be applied to other cases. A particularly important application of this section is tothe design of digital filters.

The underdetermined system with more unknowns than equations

If one poses the l p approximation problem in solving an underdetermined set of equations (case 3 from Chapter 3), it comes from defining the solution norm as

Get Jobilize Job Search Mobile App in your pocket Now!

Get it on Google Play Download on the App Store Now




Source:  OpenStax, Basic vector space methods in signal and systems theory. OpenStax CNX. Dec 19, 2012 Download for free at http://cnx.org/content/col10636/1.5
Google Play and the Google Play logo are trademarks of Google Inc.

Notification Switch

Would you like to follow the 'Basic vector space methods in signal and systems theory' conversation and receive update notifications?

Ask