Does anybody know a simple way to implement a recursive least squares function in Python?
I want a fast way to regress out a linear drift ([1 2 ... n], where n is the number of time points up until now) from my incoming signal every time it updates. RLS is typically what is used to do this, because the computing time does not increase as the number of time points increase.
The least squares fit of a line to data t[], x[] is given by
x = xbar + (C/V)*(t-tbar)
where
xbar = Sum{ x[i]} / N
tbar = sum{ t[i]} / N
V = Sum{ (t[i]-tbar)^2 } / N
C = Sum{ (x[i]-xbar)*(t[i]-tbar) } / N
You can compute xbar,tbar,V and C incrementally like this:
Initially
N = 0
xbar = tbar = C = V = 0
Incorporating data t,x:
N += 1
f = 1.0/N
dx = x - xbar
dt = t - tbar
xbar += f*dx
tbar += f*dt
V = (1.0-f)*(V + f*dt*dt)
C = (1.0-f)*(C + f*dx*dt)
Note that until you have at least two data points V will be zero, and so there is no line. Note also that each x[] could be a vector; as long as xbar and C are also computed as vectors the same formulae work.