2026-09-13
I'm used to thinking of regression as a one-off calculation. You have some regressand \( Y \), some regressors \( X \), and you want to estimate a linear model \( Y = X \beta + \epsilon \). All the data you want, \( (Y, X) \), already exists, so it's just a matter of making your estimate and you're done.
However, what if you get more data and want to make updates? Now, you're in an online world where you are making updates to your estimate as data flows in.
At each point you have a new observation \( (y_n, x_n ) \) and you want to make the best estimate possible, \( \hat{\beta}_n \), given all past data.
The naive approach would be to recalculate the usual OLS formula each time. However, this might stink in cases with limited memory or compute. There's some room to get clever with accumulating parts of it, but you're always stuck with a costly inversion every time some data rolls in.
The trick to online OLS calculations is to use the Sherman-Morrison formula:
For an invertible, square matrix \( A \) and column vectors \( u, v \).
The Sherman-Morrison formula is a big deal because it lets us skip costly matrix inversion when we update \( A \) under certain conditions. As long as we have the inverse \( A^{-1} \), we can get to the updated inverse with just simpler operations.
Broadly, we want to take some matrix update \( A + u v^\top \) and update its inverse by adding some term \( B \). This simple form means we need a lot of cancelation, since the expanded multiplication already gives an \( A^{-1} A = I_n \) term.
Now, we can't invert \( A + u v^\top \) to find \( B \), since calculating the inverse is the whole trick! However, we can write this equation to look closer to the formula by eliminating \( A \) with its inverse.
\[ B (I_n + u v^\top A^{-1}) = - A^{-1} u v^\top A^{-1} \]
Now we're very close, but we still have a matrix \( u v^\top A^{-1} \) where we want a scalar \( v^\top A^{-1} u \). The last thing to note is that, since \( v^\top A^{-1} u \) is a scalar, it can be pulled out of a large series of matrix multiplications.
Conveniently, \( A^{-1} u v^\top A^{-1} \) turns this matrix multiplication into a scalar multiplication, so it's just a matter of picking \( B \propto A^{-1} u v^\top A^{-1} \). In this case:
\[ B = - \frac{1}{1 + v^\top A^{-1} u} A^{-1} u v^\top A^{-1} \]
In our case, we'll use \( u = v = x_n \). Note that the matrix we want to invert is additively separable in each piece of new data. \[ X_n^\top X_n = \sum_{i = 1}^n x_i x_i^\top \]
What's more, the term for each \( x_n \) is exactly the format in Sherman-Morrison, except simpler.
\[ X_n^\top X_n = X_{n - 1}^\top X_{n - 1} + x_n x_n^\top \]
Since we only care about the inverse, \( (X_n^\top X_n)^{-1} \), this formula will let us make updates without tracking the full dataset.
That formula is a bit of a mess, so it'll be easier to define an intermediate.
\[ M_n = \frac{(X_{n - 1}^\top X_{n - 1})^{-1} x_n}{1 + x_n^\top (X_{n - 1}^\top X_{n - 1})^{-1} x_n} \]
From here, we can do the full update with only a few simple operations from our starting point \( (X_{n - 1}^\top X_{n - 1})^{-1} \). \[ (X_n^\top X_n)^{-1} = (X_{n - 1}^\top X_{n - 1})^{-1} - M_n x_n^\top (X_{n - 1}^\top X_{n - 1})^{-1} \]
Updating \( (X^\top X)^{-1} \) isn't very useful unless we either update \( X^\top Y \) or recalculate the whole estimate somehow.
Fortunately, we can just use the same linear algebra tricks that we've already seen! We're going to split the transposed matrix multiplication into a sum to break out the latest data.
Let's start with the first term, the one with \( X_{n-1}^\top Y_{n-1} \). Remember that the Sherman-Morrison formula for \( (X_n^\top X_n)^{-1} \) has two terms, both involving \( (X_{n-1}^\top X_{n-1})^{-1} \)? Well, we're multiplying that by \( X_{n-1}^\top Y_{n-1} \), so we're poised to get back a bunch of \( \hat{\beta}_{n-1} \), our starting point.
Now for the second term, with \( x_n y_n \). We're multiplying our formula by \( x_n \), which gives us another opportunity to use that scalar drop-out trick we used in explaining the Sherman-Morrison formula to begin with. Again, things will cancel wonderfully.
So let's put that all together now.
The update to our estimate is just our helper \( M_n \) times the prediction error at the previous model.
So to sum up, we only store the current \( \hat{\beta} \) and \( (X^\top X)^{-1} \). When we get new data \( (x_n, y_n) \), we:
Let's use a basic Python implementation.
from dataclasses import dataclass
import numpy as np
@dataclass
class ShermanMorrison:
xtx_inv: np.ndarray
beta: np.ndarray
@classmethod
def calculate(cls, X: np.ndarray, Y: np.ndarray) -> "ShermanMorrison":
xtx_inv = np.linalg.inv(X.T @ X)
beta = xtx_inv @ X.T @ Y
return ShermanMorrison(xtx_inv, beta)
def update(self, x: np.ndarray, y: np.ndarray) -> "ShermanMorrison":
x = x.reshape(-1, 1)
M = self.xtx_inv @ x / (1 + x.T @ self.xtx_inv @ x)
xtx_inv = self.xtx_inv - M @ (x.T @ self.xtx_inv)
beta = self.beta + M @ (y - x.T @ self.beta)
return ShermanMorrison(xtx_inv, beta)One thing I haven't mentioned yet is that even the online algorithm needs to start with the standard OLS calculation. Sherman-Morrison only gives us updates, not the initial point. In practice, you can just start with whatever data you have, provided it is enough to identify \( \beta \).
Here, let's use it to check. The following script takes the first 50 observations as a starting point, then runs ten updates. If everything works out, then the estimated \( \hat{\beta} \) should be the same whether it was calculated normally or with this algorithm.
if __name__ == "__main__":
rng = np.random.default_rng(0)
beta = np.array([1.0, 0.5, -1.0]).reshape(-1, 1)
N = 60
# Generate some random data to use
X = rng.normal(size=(N, beta.shape[0]))
Y = X @ beta + rng.normal(size=(X.shape[0], 1))
# Initialize with 50 points of data
first_50 = ShermanMorrison.calculate(X[:50], Y[:50])
print(f"At N=50:\n\t{first_50.beta.flatten()}")
# Update for each of the remaining points, one at a time
updated = first_50
for i in range(50, N):
updated = updated.update(X[i:i+1], Y[i])
print(f"At N={N}, Sherman-Morrison:\n\t{updated.beta.flatten()}")
# Compare to the usual OLS calculation on the full sample
full_data = ShermanMorrison.calculate(X, Y)
print(f"At N={N}:\n\t{full_data.beta.flatten()}")
difference = np.abs((updated.beta - full_data.beta)).sum()
print(f"Direct vs Online:\n\t{difference}")When I run this, I get pretty similar results at N = 60. We don't recover beta exactly, but both calculations on the full dataset give approximately the same result.
At N=50:
[ 1.02340832 0.64158546 -0.62784762]
At N=60, Sherman-Morrison:
[ 0.91446094 0.64090721 -0.65613307]
At N=60:
[ 0.91446094 0.64090721 -0.65613307]
Direct vs Online:
6.661338147750939e-16Well, just 10 updates is pretty small. There's not a lot of room for error to creep in. Let's crank N up! As you bump it up to a million or so, you'll see the error increase a bit, but it's still so small (~1e-14) that it won't matter much.
The big thing you may notice from running on large N is the time! It's much faster to calculate on the full sample than Sherman-Morrison update one at a time. This algorithm is not a substitute for the usual situation where OLS works; it's only for situations where you need the updated estimate at each point in time!