Contents
RLS and LMS
55
3.1
RLS estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
3.2
Windowed RLS for time-varying systems . . . . . . . . . . . . . . .
58
3.2.1
Sliding window RLS . . . . . . . . . . . . . . . . . . . . . .
58
3.2.2
Exponentially weighted RLS . . . . . . . . . . . . . . . . . .
60
3.2.3
Directionally weighted RLS1 . . . . . . . . . . . . . . . . . .
62
LMS and NLMS . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
3.3.1
Steepest-descent algorithm . . . . . . . . . . . . . . . . . . .
65
3.3.2
LMS algorithm . . . . . . . . . . . . . . . . . . . . . . . . .
68
3.3.3
LMS versus RLS convergence . . . . . . . . . . . . . . . . .
71
3.3.4
LMS for time-varying set-ups . . . . . . . . . . . . . . . . .
72
3.3.5
LMS variants . . . . . . . . . . . . . . . . . . . . . . . . . .
72
3.3.6
Normalized LMS (NLMS) . . . . . . . . . . . . . . . . . . .
73
3.3
1 Optional
reading.
54
Chapter 3
RLS and LMS
In chapter 2 it was mentioned that Wiener lter theory and least squares estimation
theory form the basis for the two common adaptive ltering procedures, namely the
LMS algorithm and the RLS algorithm. These two algorithms are derived and briey
analyzed in this chapter.
We reverse the order, and focus on the RLS procedure rst, which is derived from least
squares estimation theory. The LMS algorithm may then be viewed as a pruned version of the RLS algorithm. This is in contrast to the usual approach in which the LMS
algorithm is derived from Wiener lter theory, based on steepest-descent optimization
and instantaneous estimates of the process statistics.
For RLS, only the basic principles are described. Practical algorithms are derived in
later chapters.
3.1 RLS estimation
In chapter 2, least squares estimation was explained from a batch-mode processing perspective, where a complete batch of data samples is given, and then used to compute
an optimal lter. Many of the applications given in Chapter 1, however, require realtime processing, where the signal processing device has to keep pace with sampling
devices that produce new data in each time step. Hence, what is needed is a recursive
parameter estimation procedure, where new data samples -as they are fed in- are used
to re-compute or update of all the stored information, and then to produce the corresponding output.
We rst focus on a recursive version of least squares estimation, referred to as RLS (for
recursive least squares) estimation. Only the basic idea, and corresponding formulas will be given. These formulas may be used for theoretical purposes, or so-to-speak
for innite precision implementation. For nite precision implementation, alternative so-called square-root algorithms should be used instead, which are better behaved
numerically. These will be derived in Chapter 4.
We recall the least squares criterion
L
JLS (w)
e2
k
k=1
55
CHAPTER 3. RLS AND LMS
56
L
(dk uT w)2
k
k =1
where L refers to the length of the data set. To obtain a more common notation further
on, we change notation as follows
k
JLS (w)
e2
l
l =1
k
(dl uT w)2
l
l =1
This is now viewed as the LS criterion at time k. In RLS estimation, the main idea is
to compute at each time k an optimal lter from all past data. Furthermore, the lter at
time k + 1 is computed from the lter at time k, instead of from scratch which saves
on computation. We now outline how this can be done.
Assume that we have the optimal solution at time k
k
arg{min (dl uT w)2 }
l
wLS (k)
U (k) U (k)]1 U (k)T d(k)]
l =1
{z
Xuu(k)]1
}|
{z
Xdu(k)
with
2
6
6
4
d(k) = 6
d1
d2
.
.
.
7
7
7
5
U (k ) = 6
6
6
4
uT
1
uT
2
.
.
.
3
7
7
7:
5
uT
k
dk
Computing wLS (k + 1) from scratch would come down to computing
k+1
wLS (k + 1)
arg{min (dl uT w)2 }
l
l =1
U (k + 1)T U (k + 1)]1
U (k + 1)T d(k + 1)]
Xuu(k+1)]1
Xdu(k+1)
{z
}|
{z
with
2
6
6
4
d(k + 1) = 6
d1
d2
.
.
.
d k +1
7
7
7
5
U (k + 1 ) = 6
6
6
4
uT
1
uT
2
.
.
.
3
7
7
7:
5
uT+1
k
Recall that the matrix inversion Xuu (k + 1)]1 requires (N3 ) arithmetic operations1.
This could be a lot of computation do each time some new data arrives. What is sought
is an (N2 ) scheme that computes wLS (k + 1) directly from wLS (k).
Xuu() is an N N matrix.
3.1. RLS ESTIMATION
57
Observe that
Xuu(k + 1) = Xuu (k) + uk
+1
uT+1
k
and
Xdu(k + 1) = Xdu(k) + uk
+1
dk+1 :
The matrix inversion lemma2 then states that
1
+ T
Xuu (k + 1)]1 = Xuu(k)]1 Xuu1 k uT ukX1uuukk+11Xuu+k1
uk
( )]1
( )]
+ k 1
( )]
which allows us to compute Xuu (k + 1)]1 from Xuu (k)]1 in only
steps. By making use of these equations, it is readily shown that
wLS (k + 1) = wLS (k) +
Xuu(k +{z)]1uk
1
+1
Kalman gain kk+1
ON
( 2)
(dk+1 uT+1 wLS (k))
k
} |
{z
arithmetic
a priori residual k+1
The above framed formulas constitute the so-called standard recursive least squares
algorithm. Note that wLS (k + 1) is computed from wLS (k) with only (N2 ) arithmetic
operations. This is indeed what we aimed at.
We should reiterate that the above formulas should only be used for theoretical purposes. It has been shown theoretically and experimentally that the standard recursive
least squares algorithm when it is run for a long period of time suffers from numerical round-off error build-up3. For practical computations, square-root algorithms
should be used, which are derived in Chapter 4.
The updating formula for wLS shows that wLS (k + 1) is equal to wLS (k) up to a correction term with two ingredients. The so-called Kalman gain vector kk+1 gives the
direction in which wLS is modied and the a priori residual k+1 which controls the
amount by which it is changed. The a priori residual k+1 is the prediction error for
the data at time k + 1, when the weight vector wLS (k) from the previous time step is
employed. Note that if it is zero, wLS (k) ts the data at time k + 1 and so no adaptation
is necessary, i.e. wLS (k + 1) = wLS (k). A useful alternative formula for the Kalman
gain vector which is readily derived and will be used later on is
kk+1 =
Xuu (k + 1)]1uk
+1
Xuu(k)]1
1 + uT 1 Xuu (k)]1 uk
k
+
+1
uk+1 :
Finally, the above recursive least squares algorithm obviously has to be initialized properly. We may do this by rst applying batch-mode processing to the rst portion of the
data sequence, i.e. u1 ; : : : ; uN and d1 ; : : : ; dN . A simpler procedure consists in starting
from an initial estimate
Xuu(0) = I
wLS (0) = wo
2 (A + BCD)1 = A1 A1 B(C1 + DA1 B)1 DA1 .
3 It has been shown that unstability results from the loss
off errors in the updating formula.
of positive deniteness of
Xuu(k), due to round-
CHAPTER 3. RLS AND LMS
58
where is a (small) constant. It can be veried that this changes the cost function into
JLS (w(k)) = 2 w wo
2
2+
e2
l
l =1
The rst term, in effect, penalises movement away from some nominal value wo . However with time (increasing k) this rst term is insignicant compared with the second
one and can be ignored.
3.2 Windowed RLS for time-varying systems
In this section it is shown how the RLS algorithm of the previous section can be modied so that it can deal with time-varying situations too. The RLS algorithm is supplemented with a windowing or weighting technique, which effectively places a (sliding)
window over the given data. It is assumed that the system is sufciently time-invariant
within the window, so that time averaging may be applied.
We should reiterate that the RLS algorithms given here are standard RLS type algorithms which are known to be numerically unstable. In practice, square-root versions
of these algorithms should be used (to be derived in Chapter 4).
In adaptive ltering, time-varying may refer to either nonstationarity of the stochastic
processes involved, or to true time-variation in the system (plant, channel, etc.) itself.
In an acoustic echo cancellation problem, for example, the rst form of time-variation
would correspond to the nonstationarity of the speaker signals; the latter form of timevariation refers to changes in the acoustic channel (say, whenever an object is moved).
In practice, the second form is really the only relevant form of time-variation as far
as the LS/RLS procedures are concerned. By virtue of the deterministic cost function
use here, the nonstationarity of the signal statistics is unimportant and indeed effectively ignored in the LS/RLS setting. This is not the case howvere in the stochastic setting, where nonstationarity strictly disallows any form of time averaging (even if, for
example, in the acoustic echo example, the system/channel itself happens to be timeinvariant). The stochastic setting is therefore sometimes less appropriate.
3.2.1 Sliding window RLS
For time-varying systems, the basic assumption is that the system is sufcienly timeinvariant within a short time segment, so that time averaging in this segment is indeed a
sensible thing to do. This is, admittedly, a weak piece of reasoning, but given only one
sequence of data samples (which may be viewed as one realization of the stochastic
processes involved), an obvious and natural way out (for the time being).
Restricting the averaging to a suitable length window, the least squares criterion at time
k may be modied into
k
JLS (w(k))
e2
l
l =kL+1
where L is the length of the data window over which e2 is averaged. The corresponding
l
RLS algorithm is referred to as sliding window RLS, or RLS with nite memory, and is
easily derived from the standard RLS algorithm.
3.2. WINDOWED RLS FOR TIME-VARYING SYSTEMS
59
Assume that we have the optimal solution at time k
k
wLS (k)
arg{min
w
T
2
(dl ul w) }
l =kL+1
U (k) U (k)]1 U (k)T d(k)]
T
{z
}|
Xuu(k)]1
{z
Xdu(k)
with
2
6
6
4
d(k) = 6
dkL+1
dkL+2
.
.
.
7
7
7
5
U (k) = 6
6
6
4
uT +1
kL
uT +2
kL
.
.
.
3
7
7
7:
5
uT
k
dk
The aim is again to compute wLS (k + 1) with only
ON
( 2)
operations, where
k +1
wLS (k + 1)
=
=
T
2
(dl ul w) }
k=kL+2
U (k + 1)T U (k + 1)]1 U (k + 1)T d(k + 1)]
arg{min
w
{z
}|
Xuu(k+1)]1
{z
Xdu(k+1)
and
2
6
6
4
d(k + 1) = 6
dkL+2
dkL+3
.
.
.
7
7
7
5
U (k + 1 ) = 6
uT +2
kL
6 uT
6 kL+3
.
.
4
.
3
7
7
7:
5
uT+1
k
d k +1
It is observed that
Xuu(k + 1) = Xuu (k) ukL
|
{z
+1
uT +1 +uk+1 uT+1
kL
k
}
k|k
Xuu +1
and
Xdu(k + 1) = Xdu(k) u{z
kL
|
+1
dkL+1 +uk+1 dk+1 :
}
k|k
Xdu+1
Note that the update formulas are similar to the standard RLS procedure except that,
along with the term that is to be added (update), there is a term that has to be subtracted
(downdate). Not surprisingly, the complete time update procedure turns out to consist
of an update with (uk+1 ; dk+1 ) following a downdate with (ukL+1 ; dkL+1 ).
The formulas for the downdate are similar to those for the update (standard RLS). One
can easily prove that
k|k
Xuu
+1
]1 =
k 1
+1 T +
Xuu(k)]1 + Xuu1uT ukLXuuukL11uXuu+1k
k
kL
( )]1
( )]
kL+1
( )]
CHAPTER 3. RLS AND LMS
60
k|k
wLS +1 = wLS (k)
k|k
Xuu
+1
]1 ukL+1 (dkL+1 uT +1 wLS (k))
kL
The formulas for the update are those of the standard RLS algorithm (with appropriate
notational changes)
k|k
Xuu (k + 1)]1 = Xuu
+1
k|k
wLS (k + 1) = wLS +1 +
]1
k|k
k|k
Xuu +1]1 uk+1 uT+1 Xuu +1]1
k
k|k+1 1
T
1+uk+1 Xuu ] uk+1
Xuu(k + 1)]1uk
+1
(dk+1 uT+1 wk|k+1 ) :
LS
k
The above four formulas constitute a sliding window RLS algorithm.
Remark 2.10. It should be noted rightaway that downdating (unlike updating) is often a dangerous step from a numerical point of view. The minus sign in the
denominators is indicative of this: due to round-off errors the net effect of updating with (uk+1 ; dk+1 ) at time k + 1 and downdating with (uk+1 ; dk+1 ) at time
k + L + 1, is not going to be zero, so one can expect at least a linear error buildup. In conclusion, the above method, as it stands (even in square-root form), is
not really useful. Remedies will be treated later.
Remark 2.11. In a stochastic analysis (see section 2.3.3) the formula for the MSE in
the case of sliding window RLS, is obviously given as
JMSE (wLS (k)) = JMSE (wWF (k)) + JMSE (wWF (k))
|
{z
N
:
L
}
Jex (k)
where L is the window length. However, in a time-varying environment, this
only gives the excess MSE due to weight vector noise (as a result of noise in
the estimated process statistics). In addition, there will be excess MSE due to
weight vector lag, which results from the fact that the algorithm computes some
average optimal weight vector over the time window, which will be different from
the instantaneous optimal weight vector (wWF (k) = Xuu (k)1 Xdu (k)). This leads
to
JMSE (wLS (k)) = JMSE (wWF (k)) + JMSE (wWF (k))
|
{z
noise
Jex (k)
N
+ f (L)
L |{z}
}
lag
Jex (k)
lag
where Jex is expected to be proportional to L, but is obviously difcult to quanlag
tify (and so we will not make an attempt). As Jex (k) is proportional to L and
noise (k) is inversely proportional to L, the optimal choice for L should be such
Jex
that both terms are equal. However, nding such an optimal L is not an easy
task.
3.2.2 Exponentially weighted RLS
By far the most popular alternative to nite memory RLS is exponentially weighted
RLS. The least squares criterion at time k is modied to
k
JLS (w(k))
2(kl)e2
l
l =1
3.2. WINDOWED RLS FOR TIME-VARYING SYSTEMS
61
where 0 < < 1 is weighting factor or forget factor. The memory is innite, but old
values of e2 are given a small weight, and so do not play a signicant role. Roughly
l
1
speaking, only the latest 1 data points are given large weights and this is a measure
of the memory of the algorithm.
One can easily verify that
k
arg{min 2(kl ) (dl uT w)2 }
l
wLS (k)
U (k) U (k)]1 U (k)T d(k)]
l =1
{z
Xuu(k)]1
}|
{z
Xdu(k)
with
2
6
6
4
d(k) = 6
k1 d1
k2 d2
.
.
.
7
7
7
5
U (k) = 6
k1 uT
1
6 k2 uT
6
2
.
.
4
.
0 d k
3
7
7
7:
5
0 uT
k
which leads to
Xuu(k + 1) = 2 Xuu(k) + uk
+1
uT+1
k
and
Xdu(k + 1) = 2Xdu (k) + uk
+1
d k +1 :
and hence
k|k
Xuu
+1
]1 =
1
2
Xuu(k)]1
k|k
Xuu (k + 1)]1 = Xuu
+1
wLS (k + 1) = wLS (k) +
]1
k|k
k|k
Xuu +1]1 uk+1 uT+1 Xuu +1]1
k
k|k+1 1
T
1+uk+1 Xuu ] uk+1
Xuu(k + 1)]1uk
+1
(dk+1 uT+1 wLS (k))
k
The three framed formulas constitute the exponentially weighted RLS algorithm. It is
seen that the exponential weighting introduces only one simple additional step, i.e. the
weighting of the covariance matrix.
Remark 2.12. Exponential weighting not only leads to an extremely simple RLS algorithm, but it is also crucial in view of numerical stability. Even though the
standard RLS forms are known to be numerically unstable, the corresponding
square-root forms (to be derived Chapter 4), are indeed numerically stable and
owe their stability to the weighting scheme. Roughly speaking, exponential weighting indeed not only wipes out old data, but also old errors.
For this reason, some form of exponential weighting is usually also included in
CHAPTER 3. RLS AND LMS
62
the sliding window RLS algorithms4 , as well as in the directionally weighted RLS
scheme (see next section).
Remark 2.13. In a stochastic analysis (see section 2.3.1), one proves that (for small
values of 1 ) the formula for the MSE is given as
JMSE (wLS (k)) = JMSE (wWF (k)) + JMSE (wWF (k))
|
{z
2
N (1 )
)
+ f(
2 } | 1 }
{z
noise
Jex (k)
lag
Jex (k)
lag
2
where L has been replaced by the effective window length 1 . Again Jex is
2
expected to be proportional to 1 , but difcult to quantify (and so we will not
make an attempt), and hence choosing an optimal is not an easy task.
3.2.3 Directionally weighted RLS5
The exponentially weighted RLS may give rise to numerical trouble when the incoming
data are not persistently exciting. This refers to the situation where the input signal
does not contain the information we are trying to nd. In the extreme case where the
data are all-zero, e.g. from time k on, it is seen that all information is gradually wiped
out, while no new information is added, so that X1 blows up. An appealing alternauu
tive strategy consists in exponentionally discounting old data only when there is new
incoming information (with vector data this can correspond to certain subspaces or
directions). This prevents the matrix X1 from blowing up when the excitation is not
uu
uniform over the parameter space (persistent).
We will not try and give all the theoretical background of the resulting so-called directionally weighted RLS algorithm, only indicate how the formulas come about. The
algorithm is as follows and differs from the exponentially weighted RLS algorithm only
in the rst formula :
k|k
Xuu
+1
]1 =
T X
k 1
Xuu(k)]1 12 Xuu uT Xuk ukk 1uuk k
2
u
uu
( )]1
( )]
k|k
Xuu (k + 1)]1 = Xuu
+1
wLS (k + 1) = wLS (k) +
]1
( )]
k|k
k|k
Xuu +1]1 uk+1 uT+1 Xuu +1]1
k
k|k+1 1
T
1+uk+1 Xuu ] uk+1
Xuu(k + 1)]1uk
+1
(dk+1 uT+1 wLS (k))
k
The rst formula may be explained as follows. First, it is seen (see second formula) that
k|k
Xuu (k + 1)]1 is derived from the intermediate Xuu +1 ]1 by a rank-one update, which
means that the corresponding (multidimensional) quadratic cost function is affected
4 The Us are are then exponentially weighted, but still truncated matrices (i.e. with always L rows). The
basic updating formula is derived from
Xuu(k + 1) = 2 Xuu(k) 2L ukL+1 uT +1 +uk+1uT+1
kL
k
|
{z
X+
k|k 1
uu
5 Optional
reading.
3.3. LMS AND NLMS
63
k|k
only in one direction, namely Xuu +1 ]1 uk+1 . In order to apply the weighting only to
that particular direction, one can calculate the dyad corresponding to that direction,
remove it from Xuu (k + 1), weigh the dyad and then include it again, i.e.
rank = N
}|
Xuu(k + 1)
k+1|k+2
=
uu
uk+1 uT+1
uk+1 uT+1
k
k
2
+ T
Xuu(k + 1)]1uk+1 uk+1 Xuu(k + 1)]1uk+1
uT+1
k
{z
rank = N 1
The underbraced matrix represents a positive semidenite (rank N 1) quadratic form,
k|k
which is zero in the Xuu +1 ]1 uk+1 direction, i.e.
uT+1
k
k|k
Xuu
+1
T
]
uk
Xuu (k + 1) uT Xuu+k1 uk1+11uk+1
k+1
T
( + )]
k|k
Xuu
+1
uk+1 = 0
which means that the relevant dyad is indeed extracted. The dyad is then multiplied by
2 and added again. By changing the time indices, one obtains
k|k
Xuu
+1
]=
Xuu (k)] (1 2)
uT
k
uk uT
k
:
Xuu (k)]1 uk
k|k
By applying the matrix inversion lemma one obtains the Xuu +1 -formula of the given
algorithm.
Remark 2.14. Some additional exponential weighting is sometimes included, leading
to
k|k
Xuu
+1
1
] =
1
2
Xuu(k)]1 1 2
Xuu(k)]1 uk uT Xuu (k)]1
k
uT Xuu (k)]1 uk
k
where is a second forget factor (mostly 0 < < < 1).
3.3 LMS and NLMS
O
O
We have seen that for recursive processing the RLS algorithm achieves an (N) reduction in the number of arithmetic operations per time update, which is done by invoking
the matrix inversion lemma. The computational complexity is reduced to (N2 ) operations per time update. In this section, we focus on the family of least mean squares
(LMS) algorithms, which have an (N) computational complexity, i.e. achieve another (N) gain over RLS. The price to be paid is that these cheaper algorithms give
only approximate solutions at each time instant.
The original LMS algorithm is a simple one line updating algorithm, which may be
derived from the RLS algorithm by ignoring the covariance matrix update and setting
Xuu (k + 1)]1 I in the updating formula for w, leading to
wLMS (k + 1) = wLMS (k) + uk+1 (dk+1 uT+1 wLMS (k))
k
Here is the step-size parameter, the signicance of which will be explained later when
we will also derive some bounds on its value.
CHAPTER 3. RLS AND LMS
64
input signal
u[k]
u[k-1]
u[k-2]
u[k-3]
w0[k-1]
output signal
w1[k-1]
w2[k-1]
w3[k-1]
desired signal
e[k]
d[k]
w0[k]
w1[k]
w2[k]
w3[k]
a
b
w
-
a-bw
a
b w
w+ab
Figure 3.1: LMS signal ow graph
Note that this one line algorithm is remarkably simple, requiring not more than roughly
2N multiply-add operations per time update. A signal ow graph is straightforwardly
derived and shown in Figure 3.1 (see also Figure 1.21). In contrast, the convergence/stability analysis of the LMS algorithm is, however, remarkably complex (see below).
A popular alternative to the standard LMS algorithm that we will mention rightaway
is normalized LMS (NLMS), which may be derived in a similar fashion by setting
Xuu (k)]1 I, leading to6
wNLMS (k + 1) = wNLMS (k) +
1 + uT+1uk+1
k
uk+1 (dk+1 uT+1 wNLMS (k)):
k
A more general form of NLMS is
wNLMS (k + 1) = wNLMS (k) + 2 +uT
k+1 uk+1
uk+1 (dk+1 uT+1 wNLMS (k))
k
where is a constant. Clearly, NLMS may be viewed as LMS with a time-varying step
size = +uT u .
k+1 k+1
As already mentioned the convergence/stability analysis of the LMS algorithm turns
out to be a remarkably difcult mathematical task. A completely general analysis seems
to be lacking even. The LMS algorithm is usually analyzed in the (Wiener lter type)
stochastic framework, where simplifying independence assumptions can be made.
Again, such assumptions are virtually always false in practice, however the analyses
do provide results which are found to be in agreement with experiments and computer
simulations, and hence seem to provide reliable design guidelines.
We will not make an attempt to reproduce such LMS analyses, only indicate some major results. As already mentioned, the analysis is mostly done in a purely stochastic
6 By
employing the alternative formula for the Kalman gain vector.
3.3. LMS AND NLMS
65
framework. It is therefore instructive to rst indicate how LMS may be derived from
an iterative steepest-descent algorithm, applied to the Wiener ltering problem. This
is explained rst.
3.3.1 Steepest-descent algorithm
Recall that the MSE criterion is given as
JMSE (w) = {e2 }
k
and that the solution may be obtained by solving the Wiener-Hopf equations
E {u{zu } w
|
}
T
k k
Xuu
WF =
E {u d }
k k
| {z }
Xdu
(so we are back to the given statistics problem for a while). Instead of solving the set
of equations by means of a direct method (such as, e.g., Gauss elimination), one can
apply an iterative procedure. The method of steepest descent begins with picking an
initial value w(0), and then iterates as follows 7 .
w(n + 1) = w(n) +
JMSE (w)
w
i
w=w(n)
i.e.
w(n + 1) = w(n) + (Xdu Xuu w(n))
So the next guess w(n + 1) is equal to the present guess w(n) plus a correction in
the steepest-descent direction, which is known to be opposite to that of the gradient
vector in w(n). Again, is the step-size parameter, that tells how far one moves in the
steepest-descent direction. Moving too far in that direction might actually overshoot
the minimum and result in instability. Bounds on will be derived further on.
Before moving on to LMS, we give a few properties of this steepest-descent algorithm
in the following remarks. We should reiterate that we are back to the given statistics
case for while, where we assume that the process statistics (i.e. the Xuu and Xdu ) are
given, and that n is the iteration index (rather than a time index) in the iterative search
uu
for wWF = X1 Xdu .
Remark 2.15. With Xuu wWF
w(n + 1) wWF ]
Xdu , the above iteration is readily transformed into
=
=
X
X
(I uu ) w(n) wWF ]
n+1
(I uu )
w(0) wWF ]
which means that the weight-error vector gets multiplied by (I Xuu ) in each
iteration. The aim obviously is to have the weight error reduced in each step,
which is only the case when the eigenvalues of (I Xuu ) are strictly smaller
than 1 in absolute value. For analysis purposes it is convenient to make use of
7 The
factor
1
2
is used merely for convenience, i.e. to avoid a factor 2 in the second line.
CHAPTER 3. RLS AND LMS
66
the eigenvalue decomposition of Xuu . As Xuu is symmetric and assumed to be
positive denite, one has
Xuu = Quu uuQT
uu
where Quu is an orthogonal matrix (QT Quu = I) with the eigenvectors in its
uu
columns, and uu is a diagonal matrix with the corresponding eigenvalues uu =
diag{i }. This leads to
w(n + 1) wWF ]
=
=
Quu (I uu)n+1 QT w(0) wWF ]
uu
Quu diag{1 i}n+1 QT w(0) wWF ]
uu
It follows that for stability, or convergence of the steepest-descent algorithm, one
should have
1 < 1 i < 1
i:
With this, one will have w() = wWF . Since the eigenvalues i are all real and
positive8 , the necessary and sufcient condition for convergence is
0<<
2
max
where max is the largest eigenvalue of Xuu . This is illustrated in Figure 3.2,
where 1 i is plotted for various choices of . As will become clear in the
sequel, the optimal choice for convergence is
opt =
2
max +
min
such that
(1 optmin ) = (1 optmax ) =
max min
max + min
This is also illustrated in Figure 3.2. From the last formula it already follows
then whenever max
min , convergence is going to be very slow ( will be
small and hence many iterations will be required).
In the absence of knowledge of the eigenvalues of Xuu , one often can use an (admittedly sometimes conservative) upper bound for max , namely
N
max < i = trace{Xuu }
i=1
In the FIR case one has trace{Xuu } = Nxuu (0) = N {u2 }, leading to
k
0<<
8 This
2
NE {u2 }
k
follows from the fact that uu is symmetric and assumed positive denite.
3.3. LMS AND NLMS
67
1
=0
=1/ max
max
opt
min
=2/ max
-1
Figure 3.2: Step-size parameters for the steepest-descent algorithm
Remark 2.16. It is instructive to investigate the transient behavior of both the weight
vectors w(n) and the corresponding MSE JMSE (w(n)). First
QT w(n) wWF ]
uu
|
{z
w(n)
diag{1 i}n QT w(0) wWF ]
uu
|
{z
w(0)
where w(n) represents the weight-error vector expanded in terms of the basis of
eigenvectors (i.e. columns of Quu ). It is seen that the weight error in the direction
of the i-th eigenvector gets multiplied by 1 i in each iteration. The error
reduction in this i-th direction is slow if 1 i 1. This typically happens when
Xuu is ill-conditioned, i.e. when max min. The large max then results in
a small so that 1 min 1.
This effect is illustrated in Figure 3.3, where the contours of JMSE (w) are plotted
for a 2-dimensional problem, i.e. with wT = w1 w2 (see also Figure 2.6).
From
JMSE (w)
=
=
JMSE (wWF ) + (w wWF )T Xuu (w wWF )
T
JMSE (wWF ) + w uu w
it is seen that the contours (loci of constant JMSE (w)) are ellipsoids (ellipses in
2-D) with principal axes dened by the eigenvectors of Xuu . Note that, as dened
above, the major axis corresponds to the q2 direction. The gradient directions
are orthogonal to the contours. For a sufciently small , (see Figure 3.3.a)
convergence is observed in both directions with faster convergence in the q1 direction (eigenvector corresponding to the largest eigenvalue 1 ). As a result
it is seen that the error in the q1 -direction is reduced much faster than the error in
the q2 -direction, for small . The effect gets worse, when the contours are more
stretched, i.e. when 1 or 2 . If is now increased, fewer iterations are
required to reach the minimum of the error surface (see Figure 3.3.b). However
the weight vector estimates no longer appraoch the optimim smoothly but show a
tendancy to oscillate about it. If is too big, the iteration directly overshoots the
minimum in the q1 -direction. This is illustrated in Figure 3.3.c and 3.3.d, where
|11| = 1 and |11| > 1 respectively. Note that because 2 < 1 , the weight
vector in the q2 direction continues to converge smoothly to the optimum.
&
CHAPTER 3. RLS AND LMS
68
mu = 0.01
mu = 0.045
40
40
30
30
w2
50
w2
50
20
20
10
10
10
20
30
40
50
10
20
w1
30
40
50
40
50
w1
mu = 0.0469 = 2/lambda
mu = 0.05
40
40
30
30
w2
50
w2
50
20
20
10
10
10
20
30
40
50
10
20
w1
30
w1
Figure 3.3: Steepest-descent algorithm
Finally, by substituting
w(n) = diag{1 i}n w(0)
in the formula for JMSE (w), one obtains
JMSE (w(n)) = JMSE (wWF ) + N 1 i (1 i)2n wi (0)2
i=
where wi (0) is the initial error in the direction of the i-th eigenvector. The curve
obtained by plotting JMSE (w(n)) versus the iteration number n, is called the learning curve. It is seen that it consists of a sum of N exponentials. All these exponentials decay with different time constants. The worst (slowest decay) components are again those that have |1 i| 1.
3.3.2 LMS algorithm
We now return to the given data case, where Xuu and Xdu are unknown, and derive
a given data version of the steepest-descent algorithm, i.e. the LMS algorithm. We
recall the steepest descent iteration
w(n + 1)
w(n) + ( {uk dk } {uk uT }w(n)):
k
The rst (subtile) modication is that the iteration index n is changed into the time
index k, which means that the steepest-descent algorithm is turned into a recursive al-
3.3. LMS AND NLMS
69
gorithm with one weight-update (one iteration) per time step. So we have
w(k + 1)
w(k) + ( {uk dk } {uk uT }w(k)):
k
The second and last modication is that of leaving out the expectation operators. In
effect this means that we substitute an instantaneous estimate for the gradient, sometimes referred to as the stochastic gradient or the noisy gradient. This leads to
wLMS (k + 1) = wLMS (k) + uk+1(dk+1 uT+1 wLMS (k))
k
which is indeed the LMS formula we gave at the beginning of the section.
Substituting the stochastic gradient for the true gradient has a tangible impact on the
convergence properties. Also, it makes the theoretical convergence analysis dramatically more complex. As already mentioned, we will not make an attempt here to present
an LMS analysis, rather only indicate some of the major effects of the substitution. The
two issues we will be interested in are the effect in the allowable step-size , and the
resulting so-called excess MSE or mismatch.
1. First, it is instructive to rewrite the LMS formula as
wLMS (k + 1) wWF ]
T
(I uk+1 uk+1 )
wLMS (k) wWF ]
T
+ uk+1 (dk+1 uk+1 wWF ):
One can now take the expected value of both sides of the equation. The last term
then vanishes because of the orthogonality principle. This leads to
E {w
LMS (k + 1) wWF }
E { I u
T
k+1 uk+1 )
wLMS (k) wWF ]}
The right-hand side term is problematic. However it turns out one can make9
some independence assumptions that leads to
E {w
LMS (k + 1) wWF }
=
=
E {u
(I
{wLMS (k) wWF }
uu) {wLMS (k) wWF }:
(I X
T
k+1 uk+1 })
From this, one concludes that the expected behavior for LMS is the steepest
descent behavior. One therefore deduces that the LMS algorithm is convergent
in the mean, i.e. that the mean of the wLMS (k) converges to wWF as L ,
provided that the step-size parameter statises the steepest-descent bound
0<<
2
max
2. The true difference between LMS and the steepest-descent algorithm is that LMS
works with instantaneous estimates of the gradient directions, while the steepestdescent algorithm computes true gradient directions. As a result, even when the
LMS algorithm reaches wWF at some point, it will continue to bounce around,
because the estimated gradient directions are never going to be exactly equal to
zero. This in turn results in
JMSE (wLMS ()) > JMSE (wWF ):
9 The
details are omitted here.
CHAPTER 3. RLS AND LMS
70
In Figure 2.6, for example, the average point in the {w1 ; w2 }-plane will be wWF ,
but the average point on the error-performance surface will be above the JMSE (wWF )level. The excess MSE Jex () and mismatch
are then dened as follows
JMSE (wLMS ())
=
=
JMSE (wWF ) + Jex ()
JMSE (wWF ) (1 + ):
One proves that (under various technical assumptions) the LMS algorithm is
convergent in the mean square, i.e. JMSE (wLMS (k)) converges to a nite steadystate value if, and only if, the step-size parameter also satises
N
2 i
<
1:
i=1
For small values of (compared to
0<<
2
max )
this reduces to
2
N 1 i
i=
This bound is clearly much tighter than the bound for convergence in the mean.
In the FIR case one has N 1 i = trace{Xuu } = Nxuu (0) = N {u2 }, leading
i=
k
to
0<<
2
NE {u2 }
k
It is seen that the upper bound for is inversely proportional to the lter length
N.
Finally, it is found that if the above conditions are satised
Jex () = JMSE (wWF )
i
N 1 2i
i=
i
1 N 1 2
i=
{z
For small values of this reduces to
Jex () JMSE (wWF )
N
i
2 i
=1
| {z }
M
For the FIR case that is
Jex () JMSE (wWF ) N {u2 } :
2 {z k }
|
The mismatch
is often expressed as a percentage. A mismatch of 10% is
ordinarily considered to be satisfactory, meaning that the LMS algorithm produces an MSE that is (upon convergence) 10% greater than the minimum MSE,
i.e JMSE (wWF ). Hence, a choice for could be
<
0:2
N 1 i
i=
3.3. LMS AND NLMS
71
N=10
N=100
1.5
1.5
0.5
0.5
0
0
50
100
0
0
500
1000
Figure 3.4: LMS and RLS convergence
Obviously, the choice of must be based on a compromise between fast convergence and small excess mean-square error. A (much) smaller jeopardizes fast
convergence even further.
3.3.3 LMS versus RLS convergence
It has been noted that the upper bound for is inversely proportional to the lter length
N. For large N, where compared to RLS the LMS algorithm indeed achieves a large
reduction in the computational complexity (cf. (N) instead of (N2 ) operations per
update), one may expect poor convergence behavior. Furthermore, it is usual that the
larger the value of N, the larger the eigenvalue spread of the covariance matrix will be
which again means slower convergence.
This is illustrated in Figure 3.4. We consider two cases, namely N = 10 and N = 100.
Recall that the excess MSE for RLS (in the unwindowed/time-invariant case) is proportional to N , and so if 10% excess MSE is considered to be satisfactory, an RLS alk
gorithm needs k = 10 N data samples to achieve this. This is indicated in Figure 3.4,
where the excess MSE is plotted as a function of the number of samples k. The solid
line curves are the RLS curves, with 0.1 excess MSE at k = 10 N.
Let us now assume that LMS converges to the same excess MSE level, for the same
number of data points k. This allows us to check what the initial weight-vector error
should be in order to achieve this. Recall that the expected convergence behavior for
LMS is the steepest-descent behavior, so that the MSE converges (to the 10% excess
MSE level) exponentially. The smallest eigenvalue of the correlation matrix denes
the slowest decay, corresponding to (see the JMSE (w) -formula for the steepest-descent
algorithm)
(1 min)
2k
First, if (in the ideal case) all eigenvalues are equal, we have
=
0 :2
N 1 i
i=
0:2
N min
and so
(1 min)
2k
= (1
0:2 2k
) :
N
CHAPTER 3. RLS AND LMS
72
Such curves, passing through MSE=0.1 at k = 10 N are indicated in Figure 3.4, by
the upper dotted lines. It is seen that (both for N = 10 and N = 100) the initial MSE
(for L = 0) can be quite large, so that there are apparently no convergence problems.
However, if (in a more realistic setting) we assume that the eigenvalue spread is linear, i.e Ni = i min, we have
=
0 :2
N 1 i
i=
0:4
N2
min
and so
(1 min)
2k
= (1
0:4 2k
) :
N2
Such curves, passing through MSE=0.1 at k = 10 N are indicated in Figure 3.4, by the
lower dotted lines. It is seen that (for N = 10, and even more for N = 100) the initial
MSE (for k = 0) has to be close to the the nal excess MSE, so that here we denitely
have convergence problems. That is to say the LMS algorithm will take a long time
to converge with respect to the RLS algorithm (unless we are lucky and guess a good
value for the initial weight vector wLMS (0)).
3.3.4 LMS for time-varying set-ups
The LMS algorithm works with an instantaneous estimate of the gradient, which at time
k is computed from data at time k only. So in a sense, LMS works with a rectangular
window of length L = 1. It can do so because only gradient information is needed,
and no correlation/cross-correlation information (unlike in RLS). Hence the LMS algorithm may be applied to time-varying set-ups immediately.
To allow fast time-varying systems to be tracked properly, the step-size parameter
should be chosen sufciently large. With too small a , the LMS solution will lag behind. As a consequence, the lag error has a mean-square value that decreases with
an increase in . The total MSE is
JMSE (wLMS (k)) = JMSE (wWF (k)) + JMSE (wWF (k))
|
{z
noise
Jex (k)
1
N
i + f ( )
2 i
=1
}
| {z }
lag
Jex (k)
The choice of must be based on a compromise between fast tracking (large ) and
noise
small Jex (small ). For rapidly time-varying set-ups, even with the largest possible
value for , the lag error may dominate. LMS is then inappropriate and so one must
rely on more complex algorithms (such as RLS type algorithms).
3.3.5 LMS variants
A number of closely related variants to LMS exist. Normalized LMS has been mentioned already, and is also treated in the next section.
3.3. LMS AND NLMS
73
In Block LMS, the gradient vectors are averaged over several (K) iterations prior to
making adjustments of the lter coefcients. The lter coefcients are updated only
once every K iterations
wLMS ((k + 1) K) = wLMS (k K) +
K1
u
(d
uT
wLMS (k K)):
K i kK+i kK+i kK+i
=0
The averaging operation reduces the noise in the estimate of the gradient vector, and
hence the excess mean-square error. Block processing also allows the computation
of the gradients (as well as the prediction errors) in the frequency domain, which (for
large K) leads to additional computational saving (through the use of FFT techniques).
This is referred to as frequency domain LMS. Frequency domain LMS and the related
class of transform domain LMS techniques will be treated later.
To further reduce the number of multiplications in the LMS algorithm, the signed-error
algorithm has been suggested. The updating formula is
wLMS (k + 1) = wLMS (k) + uk+1 sign(dk+1 uT+1 wLMS (k)):
k
This may be viewed as a stochastic gradient algorithm that attempts to minimize the
least absolute value of the error, i.e. J (w) = {|ek |}.
An alternative way to reduce the numerical complexity is to apply the signum function
to the vector uk+1 in the updating formula (signed-regressor LMS)
wLMS (k + 1) = wLMS (k) + sign(uk+1 ) (dk+1 uT+1 wLMS (k))
k
or to the complete gradient estimate (sign-sign LMS)
wLMS (k + 1) = wLMS (k) + sign(uk+1 ) sign(dk+1 uT+1 wLMS (k)):
k
Such algorithms are simple to implement since extracting the sign of a number is easy
to do in hardware (cf multiplication). Howevere they are not very well understood, and
sometimes even show divergence for suitably pathological input signals.
3.3.6 Normalized LMS (NLMS)
From the analysis of the LMS algorithm, we know that the stable range for the step
size is inversely proportional to the sum of the eigenvalues of Xuu and hence the
power of the input signals uk . To end up with an algorithm that is robust with respect
to power variations, one may introduce a normalization of the step-size by means of a
quantity that behaves roughly as the power. One such quantity is uT uk . This results in
k
the already cited normalized LMS (NLMS) algorithm, the updating formula of which
is repeated here for convenience
wNLMS (k + 1) = wNLMS (k) + 2 +uT
k+1 uk+1
uk+1 (dk+1 uT+1 wNLMS (k))
k
The is a (small) constant that prevents overow in the event of uT uk = 0.
k
Sufce it to say that in an NLMS convergence/stability analysis it often assumed that
2 + uT+1 uk+1
k
Xuu }
2 + trace{
constant
CHAPTER 3. RLS AND LMS
74
so that the LMS fomulas can be used. From this it is seen the NLMS is convergent in
the mean as well as convergent in the mean square if
0<<2 :
This may be viewed as a normalized version of the -bounds for LMS. A 10% excess
MSE is obtained with (roughly)
< 0:2 :
Remark 2.17. We have derived the NLMS algorithm as a natural modication to the
LMS algorithm. It is interesting to not, however, that NLMS may be derived as
the solution to a specic optimization problem, which provides further insight in
the method. Given w(k), we may want to compute a w(k + 1) that minimizes
J(w(k + 1)) = 2 w(k + 1) w(k)
T
2
2
2 + (dk+1 uk+1 w(k + 1)) :
i.e. compute w(k + 1) in order to best t the current data without changing its
value too much from the previous time (i.e. w(k)). One can prove10 that the so
lution to this problem is indeed given by the NLMS formula, albeit with = 1.
Note that for 0, the optimization problem is effectively turned into a constrained optimization problem
J(w(k + 1)) = w(k + 1) w(k)
2
2
subject to
dk+1 = uT+1 w(k + 1):
k
So in the NLMS algorithm with = 0 and = 1, the w(k + 1) is the weight vector
that ts the last observations exactly, while it is as close to w(k) as possible.
Remark 2.19. 11 Interestingly, it turns out that the NMLS may also be derived from
the directionally weighted RLS algorithm, when a time varying forget factor is
used. This is a third way of deriving the NLMS algorithm, and is explained only
briey here.
Let us rst assume that in the directionally weighted RLS algorithm (section 2.3.3)
k|k
Xuu
+1
= I:
With this, one obtains
u
uT+1
k
k+1 uk+1
1
Xuu (k + 1)]1 = I 1 +kuT
+
w(k + 1) = w(k) + 1+uT 1 u
k+1 k+1
uk+1 (dk+1 uT+1 w(k))
k
10 Compare the above cost function with the one given at the end of section 3.1 (RLS initialization). The
NLMS formula then follows directly from a single, properly initialized, RLS iteration.
11 Optional reading.
3.3. LMS AND NLMS
75
The last formula is indeed the NLMS formula with = 1 and = 1. One can
verify that by setting
2+1 =
k
1
1 + uT+1uk+1
k
the next intermediate covariance matrix is also going to be equal to the identity
matrx, i.e.
k
Xuu 1|k
+
+2
=I
which closes the loop. In conclusion, NLMS (with = = 1) may be viewed as
a directionally weighted RLS algorithm, where in each time update the weighting
factor is chosen such that the covariance matrix basically remains unchanged
always. The end result is then that the covariance matrix updating formulas may
indeed be left out, so that the complexity of the algorithm is reduced by a factor
(N ).