<< Chapter < Page | Chapter >> Page > |
where is the new weighted least squares solution of [link] which is used to partially update the previous value using a convergence up-date factor which gave convergence over a larger range of around but but it was slower.
A second improvement showed that a specific up-date factor of
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 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 , decrease it. This made a significant increase in both therange of 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 '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.
If one poses the approximation problem in solving an underdetermined set of equations (case 3 from Chapter 3), it comes from defining the solution norm as
Notification Switch
Would you like to follow the 'Basic vector space methods in signal and systems theory' conversation and receive update notifications?