[6] E. D. Engeberg and S. G. Meek, “Adaptive object slip prevention for [10] A. K. Padthe, B. Drincic, J. Oh, D. D. Rizos, S. D.

os, S. D. Fassois, and D. S. Bern-

prosthetic hands through proportional-derivative shear force feedback,” stein, “Duhem modeling of friction-induced hysteresis,” IEEE Control Syst.
in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, 2008, pp. Mag., vol. 28, no. 5, pp. 90–107, 2008.
1940–1945. [11] K. Astrom and C. Canudas De Wit, “Revisiting the LuGre friction mod-
[7] E. D. Engeberg and S. G. Meek, “Adaptive sliding mode control for el,” IEEE Control Syst. Mag., vol. 28, no. 6, pp. 101–114, 2008.
prosthetic hands to simultaneously prevent slip and minimize deformation [12] F. Al-Bender, V. Lampaert, and J. Swevers, “Modeling of dry sliding
of grasped objects,” IEEE/ASME Trans. Mechatronics, vol. 18, no. 1, pp. 376– friction dynamics: From heuristic models to physically motivated models
384, 2013. and back,” Chaos: Interdiscip. J. Nonlinear Sci., vol. 14, no. 2, pp. 446–460, 2004.
[8] Z. Ding, N. Paperno, K. Prakash, and A. Behal, “An adaptive control- [13] Y. F. Liu, J. Li, Z. M. Zhang, X. H. Hu, and W. J. Zhang, “Experimental
based approach for 1-click gripping of novel objects using a robotic ma- comparison of five friction models on the same test-bed of the micro stick-
nipulator,” IEEE Trans. Control Syst. Technol., 2018. doi: 10.1109/TCST.2018 slip motion system,” Mech. Sci., vol. 6, pp. 15–28, Mar. 2015.
.2821651. [14] A. Prach, J.-J. Cabibihan, N. V. Thakor, and D. S. Bernstein, “Pareto-
[9] S. Andersson, A. Söderberg, and S. Björklund, “Friction models for slid- front analysis of a monotonic PI control law for slip suppression in a robot-
ing dry, boundary and mixed lubricated contacts,” Tribology Int., vol. 40, no. ic manipulator,” in Proc. Int. Conf. Robotics Biomimetics, 2017, pp. 2728–2733.
4, pp. 580–587, 2007. [15] J. K. Hale, Ordinary Differential Equations. New York: Wiley, 1969.

Recursive Least Squares for Real-Time Implementation

SYED ASEEM UL ISLAM and DENNIS S. BERNSTEIN

M
any estimation and control problems involve a pro- Note that (2) has the form Ax = b, where A denotes U, x
cess of the form denotes i, and b denotes Y.
In the presence of noise corrupting the data Y and U, (2)
y k = z k i, (1) may not have a solution. In this case, it is useful to replace
(2) by a least squares optimization problem of the form
where k = 0, 1, 2, f is the discrete-time step correspond-
ing to the continuous-time step size Ts, the scalar or vector k

y k ! R p is the measurement at step k, the matrix z k ! R p # n

J k (it) _ / (y i - z i it) T (y i - z i it) + (it - i 0) T R (it - i 0)
i =0
is the regressor at step k whose entries consist of current = (Y - Uit) T (Y - Uit) + (it - i 0) T R (it - i 0), (4)
and past data, and i ! R n is a column vector of n unknown
parameters. The objective is to use y k and z k to estimate the where R is a positive semidefinite (and thus, by definition,
components of i. In applications, y k and z k are corrupted symmetric) matrix, and i 0 is an initial estimate of i. Assum-
by noise, and thus (1) does not hold exactly. This motivates ing that R is chosen such that the inverse in (5) exists, the regu-
the need for the least squares estimates of i given below. larization term (it - i 0) T R (it - i 0) weights the initial estimate
The measurements y k and the data in z k are typically and ensures that J k has a unique global minimizer. In particu-
obtained from a continuous-time process and, as such, are lar, the batch least squares (BLS) minimizer of (4) is given by
available at the sample times kTs, where Ts is the sample
T -1
interval. The batch approach to this problem is to collect i opt,R = (U U + R) (U T Y + Ri 0) . (5)
a large amount of data and then apply least squares op-
timization to the collected data to compute an estimate Note that the inverse required to compute (5) is of size n # n,
of i. In particular, collecting data over the time window and thus the computational requirement of the inverse is of
i = 0, f, k, it follows from (1) that order n 3. In addition to the inverse, three matrix multipli-
cations are needed. Note also that the memory needed to
Y = Ui, (2) store U grows with k. Furthermore, if U has full column
rank, then R can be set to zero, and thus (5) becomes
where
T -1
y0 i opt,0 = (U U) U T Y. (6)
Y _ > h H , U _ > h H . (3)
z0

yk zk In the case where (2) has a solution and U has full column
rank, (6) is the unique solution of (2).
Digital Object Identifier 10.1109/MCS.2019.2900788
In many applications, computational speed and memo-
Date of publication: 17 May 2019 ry are limited. One way to alleviate these requirements is

82 IEEE CONTROL SYSTEMS MAGAZINE » JUNE 2019 1066-033X/19©2019IEEE

to recursively update an estimate of i opt,R using each addi-
tional measurement y k . A recursive algorithm of this type Summary

R
is especially convenient for real-time applications. ecursive least squares (RLS) is a technique used for
Recursive least squares (RLS) is an iterative implementa- minimizing a quadratic cost function, where the mini-
tion of BLS that significantly reduces the computational mizer is updated at each step as new data become avail-
and storage requirements of BLS. The purpose of this ar- able. RLS is more computationally efficient than batch least
ticle is to provide a statement of RLS that highlights its squares, and it is extensively used for system identification
real-time implementation along with a self-contained deri- and adaptive control. This article derives RLS and empha-
vation (see “Summary”). sizes its real-time implementation in terms of the availability
Variations of RLS have been studied for more than half of the data as well as the time needed for the computation.
a century. An early exposition is given in [1], which em-
phasizes the real-time utility of RLS relative to BLS. Ap-
plications of RLS to adaptive control are discussed in [2, P1 _ A 0-1,
pp. 41, 103]. Numerous extensions of RLS have been devel-
oped to address initialization, forgetting, and numerical it follows from Lemma 2 (see “Three Useful Lemmas”)
stability, for example, [3]–[6]. The development of RLS that with A = P 0-1, C = (1/m) I, U = z T0 , and V = z 0 that
is closest to the present article is given in [7, pp. 26–28].
P1 = 1 (P 0-1 + 1 z T0 z 0) -1
m m
RECURSIVE LEAST SQUARES
This section provides a statement and proof of the RLS al- = 1 P0 - 1 P0 z 0T (mI + z 0 P0 z T0 ) -1 z 0 P0 .
m m
gorithm. This result involves a recursive algorithm for opti-
mizing J k at each step k. The optimization of J k updates the Hence, (9) is satisfied for k = 0. Additionally, because A 0 is
estimate i k of i as measurements and data become avail- positive definite, it follows from Lemma 1 (see “Three Use-
able. As an extension of (4), the cost function (7) includes a ful Lemmas”) that the unique minimizer i 1 of J 0 is given
forgetting factor m, which provides higher weighting to by
more recent measurements and data.
-1
i 1 = -A 0 b 0
Theorem 1 = P1 z 0T y 0 + mP1 P 0-1 i 0
For all k $ 0, let z k ! R p # n and y k ! R p. Furthermore, let = P1 z T0 y 0 + P1 (P 1-1 - z T0 z 0) i 0
n
i 0 ! R , P0 ! R
n#n
be positive definite, and m ! (0, 1] . For = P1 z T0 y 0 + (I - P1 z T0 z 0) i 0
T
all k $ 0, denote the minimizer of the function = i 0 + P1 z 0 (y 0 - z 0 i 0) .

k
J k (it) _ / m k -i ( y i - z i it) T ( y i - z i it) + m k +1 (it - i 0) T P 0-1 (it - i 0) Hence, (10) is satisfied for k = 0.
(7) i =0 Now, let k $ 1. Then J k (it) can be written as

by J k (it) = it T A k it + 2b kT it + c k,

i k +1 _ argmin J k (it) . (8)

it ! R
n where
k
Then, for all k $ 0, i k +1 is given by Ak _ / m k -i z Ti z i + m k +1 P 0-1,
i =0

Pk +1 = 1 Pk - 1 Pk z kT (mI + z k Pk z Tk ) -1 z k Pk, (9)

k
m m b k _ - / m k -i z Ti y i - m k +1 P 0-1 i 0,
i =0
T
i k +1 = i k + Pk +1 z k (y k - z k i k) . (10) k
ck _ / m k -i y Ti y i + m k +1 i 0 P 0-1 i 0 .
i =0
Proof
Note that J 0 (it) can be written as Furthermore, A k and b k can be written recursively as

J 0 (it) = it T A 0 it + 2b 0T it + c 0,
A k = mA k -1 + z Tk z k,
where
b k = mb k -1 - z Tk y k .
A 0 _ z T0 z 0 + mP 0-1,
b 0 _ -z T0 y 0 - mP 0-1 i 0,
Defining
c 0 _ y T0 y 0 + mi 0 P 0-1 i 0 .

Defining Pk +1 _ A k-1,

JUNE 2019 « IEEE CONTROL SYSTEMS MAGAZINE 83

 RLS is more computationally efficient than batch least squares,
and it is extensively used for system identification
and adaptive control.

it follows from Lemma 2 (see “Three Useful Lemmas”) Note that, if m = 1, R is positive definite, and P0 _ R -1,
with A = P k-1, C = ^1 mh I, U = z Tk , and V = z k that then the RLS minimizer i k +1 given by (8) is equal to the BLS
minimizer i opt,R given by (5).
Pk +1 = [m (A k -1 + 1 z Tk z k)] -1 The notation i k +1 for the minimizer of J k emphasizes
m
the fact that i k +1, which is based on data available up to
= 1 (P k-1 + 1 z Tk z k) -1
m m step k, is not available until the update given by (9) and (10)
= 1 Pk - 1 Pk z kT (mI + z k Pk z Tk ) -1 z k Pk . is completed, which occurs at step k + 1.
m m
The derivation of RLS given by the proof of Theorem 1 is
Hence, (9) is satisfied. Furthermore, the minimizer i k +1 of an extension of the RLS derivation given in [7]. In particu-
J k is given by lar, the derivation given in [7, pp. 26–28] is based on the cost
function (7) but without the regularization term involving
-1
i k +1 =-A k b k . P0 . As can be seen in the proof of Theorem 1, this term
guarantees that A 0 is nonsingular in the first step and A k
Because A k is positive definite, it follows from Lemma 1 is nonsingular in the inductive step. In the case of BLS, the
that role of P0 is played by the matrix R.
-1
i k + 1 = -A k b k
The computational requirements of RLS are primarily
= -1 T
A ( y - mb k -1) determined by n. In particular, Pk is n # n, and thus the
k zk k

= -1 T
A ( y + mA k -1 i k) computational requirement for updating Pk given by (9)
k zk k

= A [ y + (A k - z Tk z k) i k]
-1 T is of order n 2 . Furthermore, the inverse in (9) is of size
k zk k
-1 T p # p, which, since typically p % n, is much less demand-
= A y + (I - A -k 1 z Tk z k) i k
k zk k
T ing than the inverse required by BLS. In addition, the
= i k + Pk +1 z k (y k - z k i k) .
storage requirements of RLS are of order n 2, which does
Hence, (10) is satisfied. not grow with k. Consequently, the computational and

Three Useful Lemmas

T (A + UCV) -1 = A -1 - A -1 U (C -1 + VA -1 U) -1 VA -1. 
he following result is the quadratic minimization lemma. (S4)

LEMMA 3
LEMMA 1 Let A ! R n # n be positive semidefinite, let B ! R n # m, and let
Let A ! R n # n, assume that A is positive definite, let b ! R n and C ! R n # n be positive definite. Then,
c ! R, and define f : R n " R by
[I - AB T (C + BAB T) -1 B] AB T = AB T (C + BAB T) -1 C.  (S5)
f (x) _ x T Ax + 2b T x + c.  (S1)

Then the unique minimizer of f is PROOF

Note that
x opt =-A -1 b,  (S2)

and the minimum value of f is

AB T (C + BAB T) -1 C = AB T (C + BAB T) -1 (C + BAB T - BAB T)
f (x opt) = c - b A b. 
T -1
(S3) = AB T [I - (C + BAB T) -1 BAB T]
= AB T - AB T (C + BAB T) -1 BAB T
The following result is the matrix inversion lemma [S1, p. 304]. = [I - AB T (C + BAB T) -1 B] AB T .

LEMMA 2
REFERENCE
Let A ! R n # n, U ! R n # p, C ! R p # p, and V ! R p # n, and assume [S1] D. S. Bernstein, Scalar, Vector, and Matrix Mathematics: Theory,
that A, C, and A + UCV are nonsingular. Then Facts, and Formulas. Princeton, NJ: Princeton Univ. Press, 2018.

84 IEEE CONTROL SYSTEMS MAGAZINE » JUNE 2019

memory requirements of RLS are
significantly less than those of BLS. Data Data Data Data
y0 y1 yk yk+1
ALTERNATIVE ik UPDATE φ0 φ1 φk φk+1
WITH INVERSE θ0
The following result is a variation of
Theorem 1. In this formulation, the Compute Compute t
updates of Pk and i k are reversed. P1, θ1 to Pk +1, θk +1 to Ts
0 1 k k+1
Minimize J0 (2, ..., k – 1) Minimize Jk

Theorem 2
For all k $ 0, let z k ! R p # n and y k ! R p. θ1 θk θk+1
Furthermore, let i 0 ! R n, let P0 ! R n # n
be positive definite, and let m ! (0, 1] . FIGURE 1 A real-time implementation of recursive least squares. Data at step k, which cor-
For all k $ 0, denote the minimizer responds to time t = kTs, are used to compute the minimizer i k +1 of the cost J k. Because
of the function (7) by (8). Then, for all of the time needed for the computation, the estimate i k +1 of the unknown parameter i is
not available until step k + 1.
k $ 0, i k +1 is given by

T T -1
i k +1 = i k + Pk z k (mI + z k Pk z k ) (y k - z k i k), (11) under the Dynamic Data-Driven Applications Systems
grant FA9550-16-1-0071 (http://www.1ddas.org/).
Pk +1 = 1 Pk - 1 Pk z kT (mI + z k Pk z Tk ) -1 z k Pk . (12)
m m
AUTHOR INFORMATION
Syed Aseem Ul Islam (aseemisl@umich.edu) received the
Proof B.Sc. degree in aerospace engineering from the Institute of
Using (9) to substitute Pk +1 into (10) yields Space Technology, Islamabad, Pakistan, and is currently a
Ph.D. student in flight dynamics and control at the Univer-
i k +1 = i k + 8
1 P - 1 P z T (mI + z P z T) -1 z P B z T (y - z i ) sity of Michigan, Ann Arbor. His research interests include
k k k k k k k k k k k k
m m
data-driven adaptive control for aerospace applications.
= i k + 1 6I - Pk z kT (mI + z k Pk z Tk ) -1 z k@ Pk z kT (y k - z k i k) Dennis S. Bernstein received the Sc.B. degree in ap-
m
= i k + Pk z kT (mI + z k Pk z kT) -1 (y k - z k i k), plied mathematics from Brown University, Providence,
Rhode Island, and the Ph.D. degree in control engineer-
where the last equality follows from Lemma 3 (see “Three ing from the University of Michigan, Ann Arbor, where
Useful Lemmas”). Hence, (11) holds. Finally, (12) is identical he is currently a professor in the Aerospace Engineering
to (9). Department. He is the author of Scalar, Vector, and Matrix
Mathematics (Princeton University Press, 2018). His re-
REAL-TIME IMPLEMENTATION search interests include estimation and control for aero-
OF RECURSIVE LEAST SQUARES space applications.
In many applications, it is desirable to implement RLS so
that the estimate i k is available in real time without la- REFERENCES
tency. Note that the estimate i k +1 given by (10) depends [1] A. Albert and R. W. Sittler, “A method for computing least squares esti-
mators that keep up with the data,” SIAM J. Control, vol. 3, no. 3, pp. 384–417,
on measurements available up to and including step k, 1965.
namely, y k and z k . Since time is needed to compute i k +1, [2] K. Astrom and B. Wittenmark, Adaptive Control, 2nd ed. New York: Do-
the updated estimate i k +1 is not available at time k; rather, ver, 2013.
[3] M. E. Salgado, G. C. Goodwin, and R. H. Middleton, “Modified least
it is available at the next step, namely, k + 1. Consequently, squares algorithm incorporating exponential resetting and forgetting,” Int.
the minimizer of J k is denoted by i k +1, where the subscript J. Control, vol. 47, no. 2, pp. 477–491, 1988.
k + 1 conveys the fact that the minimizer of J k is not avail- [4] S. Dasgupta and Y.-F. Huang, “Asymptotically convergent modified re-
cursive least-squares with data-dependent updating and forgetting factor
able until step k + 1. In contrast, the notation used in [7, for systems with bounded noise,” IEEE Trans. Inf. Theory, vol. 33, no. 3, pp.
p. 27] is i k . Figure 1 shows how the measurements and data 383–392, 1987.
that are available at step k are used during the time inter- [5] P. Stoica and P. Ahgren, “Exact initialization of the recursive least-
squares algorithm,” Int. J. Adapt. Control Signal Process., vol. 16, pp. 219–230,
val [kTs, (k + 1) Ts] to compute the next estimate i k +1 . Jan. 2002.
[6] L. Ljung and T. Soderstrom, Theory and Practice of Recursive Identification.
ACKNOWLEDGMENTS Cambridge, MA: MIT Press, 1985.
[7] C. Johnson, Lectures on Adaptive Parameter Estimation. Englewood Cliffs,
The authors thank Ankit Goel, Jonathan How, and Tam NJ: Prentice Hall, 1988.
Nguyen for helpful suggestions. This research was partial-

ly supported by the Air Force Office of Scientific Research

JUNE 2019 « IEEE CONTROL SYSTEMS MAGAZINE 85