The Woodbury Matrix Identity
Updating an inverse cheaply — the trick behind Kalman filters and online regression.
On this page
The one-sentence idea: if you already know the answer to a hard problem, and the problem only changes a little, there's often a shortcut to the new answer that's much cheaper than solving the whole thing again from scratch — Woodbury is exactly that shortcut, specifically for matrix inverses.
Beginner: inverting a matrix (section 1.6) is expensive — roughly O(n³). But what if you already have the inverse of a matrix, and you only need to update it slightly (say, adding one new data point to a dataset)? Recomputing the whole inverse from scratch every time would be wasteful.
Intermediate: the Woodbury matrix identity (also called the matrix inversion lemma) gives a formula for the inverse of a matrix after a low-rank update, expressed entirely in terms of the original inverse — turning an expensive full O(n³) re-inversion into a much cheaper update.
Advanced: this is exactly the trick that makes online/streaming algorithms practical: a Kalman filter updates its covariance estimate every time a new measurement arrives, and a Gaussian process can incorporate one new observation, both without ever recomputing a full matrix inverse from scratch at every single step.
Dense as this looks, the key point is simple: if A is n×n but U, C, V represent only a small (rank-k) update, the right-hand side only ever needs to invert a small k×k matrix instead of a full n×n one.
The cleanest proof of a claimed inverse formula is to just multiply it by the original matrix and check the result is the identity. Let B = A⁻¹ − A⁻¹U(C⁻¹ + VA⁻¹U)⁻¹VA⁻¹ be the claimed inverse of (A + UCV). Multiply them together and distribute:
Group the last three terms, all of which share a common right-hand factor of (C⁻¹ + VA⁻¹U)⁻¹VA⁻¹ once UCVA⁻¹ is rewritten as UC(C⁻¹ + VA⁻¹U)(C⁻¹ + VA⁻¹U)⁻¹VA⁻¹:
The bracket is exactly CC⁻¹ − I = 0, so everything past the leading I vanishes and (A + UCV)B = I — confirming B really is the inverse, purely by algebra, with no assumption needed beyond the relevant inverses existing.
Where this is used: recursive least squares and Kalman filters (an O(n³) covariance re-inversion at every timestep would make real-time filtering infeasible), Gaussian process regression when adding one new observation, and low-rank adaptation (LoRA) style updates to large weight matrices in deep learning.
Let A = [[2, 0], [0, 2]] (so A⁻¹ = [[0.5, 0], [0, 0.5]]), and update it with u = [1, 1] so A_new = A + uuᵀ = [[3, 1], [1, 3]]. Direct inversion gives A_new⁻¹ = [[0.375, −0.125], [−0.125, 0.375]]. Sherman-Morrison instead computes A⁻¹u = [0.5, 0.5], then 1 + uᵀA⁻¹u = 1 + 0.5 + 0.5 = 2, giving A_new⁻¹ = A⁻¹ − (A⁻¹u)(A⁻¹u)ᵀ / 2 = [[0.5,0],[0,0.5]] − [[0.125,0.125],[0.125,0.125]] = [[0.375, −0.125], [−0.125, 0.375]] — the exact same answer, using only vector operations on the already-known A⁻¹, never re-inverting the full 2×2 matrix.
This special rank-1 case is called the Sherman-Morrison formula — the same idea as Woodbury, just for the simplest possible update.
#include <cmath>
#include <cstdio>
#include <vector>
#include <random>
#include <algorithm>
std::vector<double> solveLinearSystem(std::vector<std::vector<double>> A, std::vector<double> b) {
int n = static_cast<int>(b.size());
for (int col = 0; col < n; ++col) {
int pivot = col;
for (int r = col + 1; r < n; ++r)
if (std::fabs(A[r][col]) > std::fabs(A[pivot][col])) pivot = r;
std::swap(A[col], A[pivot]);
std::swap(b[col], b[pivot]);
for (int r = col + 1; r < n; ++r) {
double factor = A[r][col] / A[col][col];
for (int c = col; c < n; ++c) A[r][c] -= factor * A[col][c];
b[r] -= factor * b[col];
}
}
std::vector<double> x(n, 0.0);
for (int row = n - 1; row >= 0; --row) {
double s = b[row];
for (int c = row + 1; c < n; ++c) s -= A[row][c] * x[c];
x[row] = s / A[row][row];
}
return x;
}
std::vector<std::vector<double>> matrixInverse(const std::vector<std::vector<double>>& A) {
int n = static_cast<int>(A.size());
std::vector<std::vector<double>> inv(n, std::vector<double>(n));
for (int k = 0; k < n; ++k) {
std::vector<double> e(n, 0.0);
e[k] = 1.0;
std::vector<double> col = solveLinearSystem(A, e);
for (int i = 0; i < n; ++i) inv[i][k] = col[i];
}
return inv;
}
int main() {
std::mt19937 rng(0);
std::uniform_real_distribution<double> unit(0.0, 1.0);
const int n = 200;
std::vector<std::vector<double>> A(n, std::vector<double>(n, 0.0));
for (int i = 0; i < n; ++i) A[i][i] = 2.0;
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
A[i][j] += unit(rng) * 0.01;
auto A_inv = matrixInverse(A); // expensive, done once
std::vector<double> u(n);
for (double& x : u) x = unit(rng);
std::vector<std::vector<double>> A_new(n, std::vector<double>(n));
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
A_new[i][j] = A[i][j] + u[i] * u[j];
// Sherman-Morrison: A_inv - (A_inv u)(A_inv u)^T / (1 + u^T A_inv u)
std::vector<double> AinvU(n, 0.0);
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j) AinvU[i] += A_inv[i][j] * u[j];
double denom = 1.0;
for (int i = 0; i < n; ++i) denom += u[i] * AinvU[i];
std::vector<std::vector<double>> A_new_inv_fast(n, std::vector<double>(n));
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
A_new_inv_fast[i][j] = A_inv[i][j] - (AinvU[i] * AinvU[j]) / denom;
for (int k : {0, 5}) {
std::vector<double> e(n, 0.0);
e[k] = 1.0;
std::vector<double> xDirect = solveLinearSystem(A_new, e);
bool ok = true;
for (int i = 0; i < n; ++i)
if (std::fabs(xDirect[i] - A_new_inv_fast[i][k]) > 1e-6) ok = false;
std::printf("column %d matches: %s\n", k, ok ? "true" : "false");
}
return 0;
}- Kalman filters (robotics, GPS, finance) use exactly this identity to update state covariance estimates efficiently every time a new sensor reading arrives.
- Online/recursive least squares updates a regression model's solution as new data streams in, without recomputing the full normal-equation inverse from section 1.6 each time.
- Gaussian process libraries use Woodbury-style updates to add new training points incrementally, avoiding a full
O(n³)Cholesky refactorization (section 1.15) every time.
- Applying Woodbury when the "update" isn't actually low-rank — the whole benefit disappears if k is comparable to n; it's specifically a low-rank-update trick.
- Forgetting numerical stability still matters — repeated incremental updates can accumulate floating-point drift, and long-running online systems often periodically recompute a fresh, exact inverse to correct for it.
Going deeper
The Woodbury identity is, algebraically, a generalization of the simple scalar fact that 1/(a+bc) can be rewritten in terms of 1/a when bc is "small" relative to a — matrix inversion has a genuine analogue of this same idea, just dressed up in more notation.
At the master level: the Woodbury identity is one of the standard tools that makes Bayesian linear regression and Gaussian process regression tractable at scale — both rely on repeatedly manipulating covariance matrices under low-rank updates, and naive full re-inversion at every step would make either approach computationally infeasible for any real dataset size.