18
$\begingroup$

Can anyone point me in the direction of an online (recursive) algorithm for Tikhonov Regularisation (regularised least squares)?

In an offline setting, I would calculate $\hat\beta=(X^TX+λI)^{−1}X^TY$ using my original data set where $λ$ is found using n-fold cross validation. A new $y$ value can be predicted for a given $x$ using $y=x^T\hat\beta$.

In an online setting I continually draw new data points. How can I update $\hat\beta$ when I draw new additional data samples without doing a full recalculation on the whole data set (original + new)?

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ Your Tikhonov-regularized least squares is perhaps more commonly called Levenberg-Marquardt in statistical circles, even when applied to pure-linear problems (as here). There's a paper about online Levenberg Marquardt here. I don't know if that's any help. $\endgroup$
    – Glen_b
    Commented Jan 10, 2014 at 1:46

5 Answers 5

Reset to default
17
$\begingroup$

$\hat\beta_n=(XX^T+λI)^{−1} \sum\limits_{i=0}^{n-1} x_iy_i$

Let $M_n^{-1} = (XX^T+λI)^{−1}$, then

$\hat\beta_{n+1}=M_{n+1}^{−1} (\sum\limits_{i=0}^{n-1} x_iy_i + x_ny_n)$ , and

$M_{n+1} - M_n = x_nx_n^T$, we can get

$\hat\beta_{n+1}=\hat\beta_{n}+M_{n+1}^{−1} x_n(y_n - x_n^T\hat\beta_{n})$

According to Woodbury formula, we have

$M_{n+1}^{-1} = M_{n}^{-1} - \frac{M_{n}^{-1}x_nx_n^TM_{n}^{-1}}{(1+x_n^TM_n^{-1}x_n)}$

As a result,

$\hat\beta_{n+1}=\hat\beta_{n}+\frac{M_{n}^{−1}}{1 + x_n^TM_n^{-1}x_n} x_n(y_n - x_n^T\hat\beta_{n})$

Polyak averaging indicates you can use $\eta_n = n^{-\alpha}$ to approximate $\frac{M_{n}^{−1}}{1 + x_n^TM_n^{-1}x_n}$ with $\alpha$ ranges from $0.5$ to $1$. You may try in your case to select the best $\alpha$ for your recursion.

I think it also works if you apply a batch gradient algorithm:

$\hat\beta_{n+1}=\hat\beta_{n}+\frac{\eta_n}{n} \sum\limits_{i=0}^{n-1}x_i(y_i - x_i^T\hat\beta_{n})$

Improve this answer
$\endgroup$
3
  • $\begingroup$ What if I update my regressor each time with batch samples of new data, where each successive batch is drawn from a slightly different distribution? i.e. non IID. In this case I would want the regressor to take into account the new data, but not affect its predictions in the locality of the old data (previous batches)? Can you point me to any literature you may feel useful? $\endgroup$
    – rnoodle
    Commented Jan 15, 2014 at 19:18
  • $\begingroup$ Good question, but sorry currently I cannot tell how much would it affect your model if you are still using the batch gradient formula in the answer, or approximating by applying the matrix form directly: eta^(-alpha)*X(Y-X'beta_n) where X, Y are your new batch samples $\endgroup$
    – lennon310
    Commented Jan 15, 2014 at 19:34
  • $\begingroup$ hi, it seems that the regularization coefficient does not be involved in the recursive update formula? or does it only matter in the initialization of the inverse of M matrix? $\endgroup$
    – Matics
    Commented Jul 30, 2019 at 14:21
6
$\begingroup$

A point that no one has addressed so far is that it generally doesn't make sense to keep the regularization parameter $\lambda$ constant as data points are added. The reason for this is that $\| X \beta -y \|^{2}$ will typically grow linearly with the number of data points, while the regularization term $\| \lambda\beta \|^{2}$ won't.

Improve this answer
$\endgroup$
4
  • $\begingroup$ That's an interesting point. But exactly why does it "not make sense"? Keeping $\lambda$ constant surely is mathematically valid, so "not make sense" has to be understood in some kind of statistical context. But what context? What goes wrong? Would there be some kind of easy fix, such as replacing the sums of squares with mean squares? $\endgroup$
    – whuber ♦
    Commented Aug 6, 2015 at 16:59
  • $\begingroup$ Replacing the sum of squares with a scaled version (e.g. the mean squared error) would make sense, but simply using recursive least squares won't accomplish that. $\endgroup$
    – Brian Borchers
    Commented Aug 6, 2015 at 19:09
  • $\begingroup$ As for what would go wrong, depending on your choice of $\lambda$, you'd get a very underregularized solution with a large number of data points or a very overregularized solution with a small number of data points. $\endgroup$
    – Brian Borchers
    Commented Aug 6, 2015 at 19:10
  • 1
    $\begingroup$ One would suspect that, but if $\lambda$ is tuned initially after receiving $n$ data points and then more data points are added, whether the resulting solutions with more data points and the same $\lambda$ are over- or under-regularized would depend on those new datapoints. This can be analyzed by assuming the datapoints act like an iid sample from a multivariate distribution, in which case it appears $\lambda$ should be set to $N/n$ at stage $N$. This would change the updating formulas, but in such a regular and simple way that efficient computation might still be possible. (+1) $\endgroup$
    – whuber ♦
    Commented Aug 6, 2015 at 19:16
3
$\begingroup$

Perhaps something like Stochastic gradient descent could work here. Compute $\hat{\beta}$ using your equation above on the initial dataset, that will be your starting estimate. For each new data point you can perform one step of gradient descent to update your parameter estimate.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I have since realise that SGD (perhaps minibatch) is the way to go for online problems like this i.e. updating function approximations. $\endgroup$
    – rnoodle
    Commented Sep 11, 2015 at 10:22
2
$\begingroup$

Here is an alternative (and less complex) approach compared to using the Woodbury formula. Note that $X^TX$ and $X^Ty$ can be written as sums. Since we are calculating things online and don't want the sum to blow up, we can alternatively use means ($X^TX/n$ and $X^Ty/n$).

If you write $X$ and $y$ as :

$$ X = \begin{pmatrix} x_1^T \\ \vdots \\ x_n^T \end{pmatrix}, \quad y = \begin{pmatrix} y_1 \\ \vdots \\ y_n \end{pmatrix}, $$

we can write the online updates to $X^TX/n$ and $X^Ty/n$ (calculated up to the $t$-th row) as:

$$ A_t = \left(1 - \frac{1}{t}\right) A_{t-1} + \frac{1}{t}x_t x_t^T, $$

$$ b_t = \left(1 - \frac{1}{t}\right) b_{t-1} + \frac{1}{t}x_t y_t. $$

Your online estimate of $\beta$ then becomes

$$\hat\beta_t = (A_t + \lambda I)^{-1}b_t.$$

Note that this also helps with the interpretation of $\lambda$ remaining constant as you add observations!

This procedure is how https://github.com/joshday/OnlineStats.jl computes online estimates of linear/ridge regression.

Improve this answer
$\endgroup$
1
$\begingroup$

In linear regression, one possibility is updating the QR decomposition of $X$ directly, as explained here. I guess that, unless you want to re-estimate $\lambda$ after each new datapoint has been added, something very similar can be done with ridge regression.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.