Least Squares with Blocks / Updates

Question

I have a continuous-time system that I want to fit via least squares. I just send $N$ digital samples $x[n]$ through the system and receive (via analog signal chain, ADC etc) $N$ digital samples $y[n]$. My system can be described as

$$ \mathbf{y} = \mathbf{X} \mathbf{c} $$

where $\mathbf{X}$ is a matrix composed of $x[n]$. Hence I solve for the coefficients via:

$$ \hat{\bf c} = (\mathbf{X}^{H} \mathbf{X})^{-1} \mathbf{X}^{H} \mathbf{y} $$

which works wonderful and it's straight forward and easy. Note - $(\cdot)^H$ denotes the Hermitian transpose

Now $\mathbf{y}$ are heavily contaminated by noise so I need many samples. But on each receive call I can only receive a certain maximum number of samples. Even if I would receive an infinite amount - calculating the pseudo inverse becomes very expensive.

What is the easiest way and direct generalization of the above least squares set up so that I can iteratively receive blocks and refine the estimate after each block? The method should match a hypothetical, full least squares solution as closely as possible.

There are a couple of concepts which may be related but I am not sure how they relate to each other:

1. I could first receive, say 1000 batches of repetitive $\mathbf{y}$ and average them and solve one least squares. Is this equivalent to solving a larger least squares system with all samples? Disadvantage: I need to send repetitive data.
2. I am not looking for LMS or similar sampled-based methods. Ideally I would just use Least Squares ... if the number of measurements would not be as excessive
3. What about RLS oder IRLS? To my understanding, they are also sample based (not block based) methods.
4. I have heard about gradient descent + Newton type:

$$ \mathbf{c}_{n+1} = \mathbf{c}_{n} = \mu \left(\frac{\partial^2 J}{\partial \mathbf{c}^{*} \partial \mathbf{c}}\right)^{-1} \frac{\partial J}{\partial \mathbf{c}^{*}} $$

which, in my linear model would expand to:

$$ \mathbf{c}_{n+1} = \mathbf{c}_{n} = \mu (\mathbf{X}^{H}\mathbf{X})^{-1} \mathbf{X} \mathbf{e} $$

where $\mathbf{e}$ is the error vector between actual measurements and model with the current estimate $\mathbf{c}_n$. It requires that the unknowns are smaller than the number of samples in each iteration. This looks very close to what I am looking at ... but I am not sure if it is the right thing and how it relates to the other methods.

Answer 1

RLS is derived using the Woodbury Identity for rank-1 update.

I would go back and re-derive the RLS approach in light of Woodbury and generalize it to the block update case. I vaguely recall doing this before, but I can't recall when/what-for at the moment.

Here's the identity:

$$ (A+UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1} $$

[Aside: I will try to use upper case for matrices, lower case with underbar for column vectors, and simply lower case for scalars - though I don't think we will have any scalars. I also use superscript $H$ to denote Hermitian transpose]

Now, consider your data vector and matrix at time 'k' and the resulting LS solution:

$$ \underline{c}_k = arg\min_{c_k}||\underline{y}_k - X_k\underline{c}_k||_2^2 $$

Then, at time 'k+L' - after another block of 'L' samples, we have

$$ \underline{c}_{k+L} = arg\min_{\underline{c}_{k+L}}||\underline{y}_{k+L} - X_{k+L}\underline{c}_{k+L}||_2^2 $$

Next, define some terms for the block update

$$ \underline{y}_{k+L} = \left[\begin{array}{l}\underline{y}_k\\ \underline{\tilde{y}}_{k+L}\end{array}\right] $$ $$ {X}_{k+L} = \left[\begin{array}{l}X_k\\ {\tilde{X}}_{k+L}\end{array}\right] $$

Now, solve for $\underline{c}_{k+L}$ using the Woodbury identity (in line 3), state, and the solution from time $k$:

\begin{eqnarray} \underline{c}_{k+L} &=& (X^H_{k+L}X_{k+L})^{-1} X^H_{k+L}\underline{y}_{k+L}\\ &=&(X^H_{k}X_{k}+\tilde{X}^H_{k+L}\tilde{X}_{k+L})^{-1} (X^H_{k}\underline{y}_{k} + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L})\\ &=&[(X^H_{k}X_{k})^{-1} + (X^H_{k}X_{k})^{-1}\tilde{X}^H_{k+L}(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}\tilde{X}_{k+L}(X^H_{k}X_{k})^{-1}][X^H_k\underline{y}_k + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L}]\\ &=&\underline{c}_k+W^H_{k+L}\underline{\tilde{y}}_{k+L} + [W^H_{k+L}(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}W_{k+L}](\underline{v}_k+\tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L}) \end{eqnarray}

where

$$ W_{k+L} = \tilde{X}_{k+L}(X^H_{k}X_{k})^{-1} $$ $$ \underline{v}_k = X^H_k\underline{y}_k $$ $$ \underline{v}_{k+L} = \underline{v}_k + \tilde{X}^H_{k+L}\underline{\tilde{y}}_{k+L} $$

I would keep state for the latter intermediates and for $(X^H_{k}X_{k})^{-1}$. The expensive part here is usually the calculation of $(I+\tilde{X}_{k+L}\tilde{X}^H_{k+L})^{-1}$, since it is a rank $L>1$ update.

Apologies if I have typos above. I was deriving as I typed, and eventually I couldn't see the wysiwyg output anymore. ;) Hopefully, you can see the jist of the approach, nonetheless. There certainly could be further intermediate and state optimization. I didn't go back to RLS to remind me of how it was wrung to minimize operations, but it probably goes somewhere along these lines.

Finally, given the above block update approach, the question for you is - which is preferred: $L$ rank 1 updates (RLS) or 1 rank $L$ ("block RLS" as above) update? This depends on your situation and environment. A block approach may be better if the rank-1 updates would be done in an interpreted loop (e.g. Matlab), and one may be better than the other for accumulated precision errors.

Answer 2

If you have a good guess at the inverse of $A$ called $G_0$ to begin with, the approach you might be thinking of looks like this:

$$ P_n = G_n A = I + E_n $$

$$ ( I - E_n ) G_n A = ( I - E_n )( I + E_n ) = I - E_n^2 \approx I $$

Note that:

$$ ( I - E_n ) = ( 2 I - P_n ) $$

Therefore:

$$ G_{n+1} = ( 2 I - P_n ) G_n $$

This will rapidly converge to a better inverse and requires way fewer calculations. I'm not sure what it is called.

Ced

Answer 3

You can check out Block-based LMS. Without further information, I think it's hard to come up with a more robust suggestion.

Hint: Kalman filter is indeed optimal (and is suggested by others), but only for Gaussian systems (which is not stated in the question).