To find the minimum lplp approximate solution,
we propose the iterative reweighted least squared (IRLS) error algorithm
which starts with unity weighting, W=IW=I, solves for an
initial xx with Equation 15, calculates a new error from Equation 13,
which is then used to set a new weighting matrix WW
W
=
d
i
a
g
(
w
n
)
(
p
-
2
)
/
2
W
=
d
i
a
g
(
w
n
)
(
p
-
2
)
/
2
(16)to be used in the next iteration of Equation 15.
Using this, we find a new solution xx and repeat until convergence
(if it happens!). To converge, this should be a contraction map which
converges to a fixed point which is the solution. A simple Matlab program
that implements this algorithm is:
% m-file IRLS0.m to find the optimal solution to Ax=b
% minimizing the L_p norm ||Ax-b||_p, using IRLS.
% csb 11/10/2012
function x = IRLS0(A,b,p,KK)
if nargin < 4, KK=10; end;
x = pinv(A)*b; % Initial L_2 solution
E = [];
for k = 1:KK % Iterate
e = A*x - b; % Error vector
w = abs(e).^((p-2)/2); % Error weights for IRLS
W = diag(w/sum(w)); % Normalize weight matrix
WA = W*A; % apply weights
x = (WA'*WA)\(WA'*W)*b; % weighted L_2 sol.
ee = norm(e,p); E = [E ee]; % Error at each iteration
end
plot(E)
Program 1.
This core idea has been repeatedly proposed and developed in different
application areas over the past 50 years with a variety of success [10], [35], [3]. Used
in this basic form, it reliably converges for 1.5<p<31.5<p<3. In 1970 [29], a
modification was made to partially update the solution each iteration using
x
(
k
)
=
q
x
^
(
k
)
+
(
1
-
q
)
x
(
k
-
1
)
x
(
k
)
=
q
x
^
(
k
)
+
(
1
-
q
)
x
(
k
-
1
)
(17)where x^x^ is the new weighted least squares solution of Equation 15
which is used to only partially update the previous value x(k-1)x(k-1).
The first application of this partial update optimized the value for qq each
iteration to give a more robust convergence but is slowed the total algorithm
considerably.
A second improvement used a specific up-date factor of
q
=
1
p
-
1
q
=
1
p
-
1
(18)to significantly increased the speed of convergence. With this particular factor, the
algorithm becomes a form of Newton's method which has quadratic convergence [28], [7]
but initial convergence was less robust.
A third modification applied homotopy [6], [43], [41], [45], [28], [23]
by starting with a value for pp 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<2p<2, decrease it.
This improved the initial performance of the algorithm and made a significant increase in both the
range of pp that allowed convergence and in the speed of calculations. Some
of the history and details can be found applied to digital filter design in [3], [7].
A Matlab program that implements these ideas applied to our pseudoinverse
problem with more equations than unknowns 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)
Program 2.
This can be modified to use different pp's in different bands of equations or to
use weighting only when the error exceeds a certain threshold to achieve a
constrained LS approximation [3], [7], [53].
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 to
the design of digital filters.
Preformatted PDF for this module
