|
| 1 | +# Incremental QR Least Squares |
| 2 | + |
| 3 | +## Overview & Motivation |
| 4 | + |
| 5 | +Identification and calibration runs produce long logs — motor current against acceleration, joint |
| 6 | +torque against motion — that a microcontroller cannot store. Accumulating the normal equations fits |
| 7 | +in memory but squares the condition number, which single precision cannot afford. Real regressors |
| 8 | +are also often rank-deficient: a parameter never shows up in the data, or only in a fixed |
| 9 | +combination with another, and a plain solver then returns garbage. |
| 10 | + |
| 11 | +Incremental QR least squares rotates each arriving row into a small triangular summary and forgets |
| 12 | +it: memory stays fixed however long the log runs and the data matrix is never squared. The final |
| 13 | +solve decomposes the summary into singular directions, inverts the strong ones, drops the weak ones, |
| 14 | +and reports which parameter combinations the data could actually determine. |
| 15 | + |
| 16 | +## Mathematical Theory |
| 17 | + |
| 18 | +### Streaming QR Factor |
| 19 | + |
| 20 | +Rows $(\mathbf{y}_k^T, \tau_k)$ with $\mathbf{y}_k \in \mathbb{R}^P$ arrive one at a time, $m$ |
| 21 | +unbounded. The task is $\min_{\boldsymbol{\pi}} \|Y\boldsymbol{\pi} - \boldsymbol{\tau}\|_2$ over the |
| 22 | +stacked rows, taking the minimum-norm minimizer. Write |
| 23 | +$Y = Q \begin{bmatrix} R \\ 0 \end{bmatrix}$ with $Q$ orthogonal and $R$ upper triangular |
| 24 | +$P \times P$, and split $Q^T\boldsymbol{\tau} = (\mathbf{d}, \boldsymbol{\tau}_\perp)$. Orthogonal maps |
| 25 | +preserve norms, so |
| 26 | + |
| 27 | +$$\|Y\boldsymbol{\pi} - \boldsymbol{\tau}\|^2 = \|R\boldsymbol{\pi} - \mathbf{d}\|^2 + \rho, |
| 28 | +\qquad \rho = \|\boldsymbol{\tau}_\perp\|^2$$ |
| 29 | + |
| 30 | +The problem lives in $R$, $\mathbf{d}$ and $\rho$, the part of the right-hand side no |
| 31 | +$\boldsymbol{\pi}$ can reach: $P^2 + P + 2$ scalars with the row count. $Q$ is never formed. |
| 32 | + |
| 33 | +**Row update.** Place the new row beneath $[\,R \mid \mathbf{d}\,]$. For $j = 1, \dots, P$, skip the |
| 34 | +step if $y_j = 0$; otherwise rotate row $j$ of $R$ and the incoming row in their common plane: |
| 35 | + |
| 36 | +$$\begin{bmatrix} c & s \\ -s & c \end{bmatrix} |
| 37 | +\begin{bmatrix} R_{j,j:P} & d_j \\ y_{j:P} & \tau \end{bmatrix}, |
| 38 | +\qquad c = \frac{R_{jj}}{\ell},\; s = \frac{y_j}{\ell},\; \ell = \sqrt{R_{jj}^2 + y_j^2}$$ |
| 39 | + |
| 40 | +This zeroes $y_j$, leaves $R_{jj} = \ell \ge 0$ and touches only columns $j$ to $P$. After the sweep |
| 41 | +the incoming row is zero in every column, so its leftover right-hand side $\tilde{\tau}$ reads |
| 42 | +$0 \cdot \boldsymbol{\pi} = \tilde{\tau}$ and adds $\tilde{\tau}^2$ to $\rho$ whatever |
| 43 | +$\boldsymbol{\pi}$ is. A zero pivot degenerates to a swap and a zero pair to the identity. |
| 44 | + |
| 45 | +**Why not normal equations.** $Y^TY = R^TR$, yet the singular values of $R$ are those of $Y$: the |
| 46 | +factor works at $\kappa(Y)$, the normal equations at $\kappa(Y)^2$. With single-precision |
| 47 | +$\epsilon \approx 1.2 \times 10^{-7}$ they fail near $\kappa(Y) \approx \epsilon^{-1/2} \approx 3 |
| 48 | +\times 10^3$; the factor lasts to about $\epsilon^{-1}$. |
| 49 | + |
| 50 | +### Rank-Revealing Solve |
| 51 | + |
| 52 | +Take $R = U\Sigma V^T$ with $\sigma_0 \ge \dots \ge \sigma_{P-1} \ge 0$. In the coordinates |
| 53 | +$\mathbf{z} = V^T\boldsymbol{\pi}$, |
| 54 | + |
| 55 | +$$\|R\boldsymbol{\pi} - \mathbf{d}\|^2 = \sum_{i=0}^{P-1} \left(\sigma_i z_i - \mathbf{u}_i^T\mathbf{d}\right)^2$$ |
| 56 | + |
| 57 | +The rank is $r = \#\{\, i : \sigma_i > \text{tol} \cdot \sigma_0 \,\}$. For $i < r$ the term is zeroed by |
| 58 | +$z_i = \mathbf{u}_i^T\mathbf{d} / \sigma_i$; for $i \ge r$ it cannot be removed and $z_i$ is free, so |
| 59 | +the minimum-norm choice is $z_i = 0$. The discarded terms are charged to the residual: |
| 60 | + |
| 61 | +$$\hat{\boldsymbol{\pi}} = \sum_{i < r} \mathbf{v}_i \, \frac{\mathbf{u}_i^T\mathbf{d}}{\sigma_i}, |
| 62 | +\qquad |
| 63 | +\text{RMS} = \sqrt{\frac{\rho + \sum_{i \ge r} (\mathbf{u}_i^T\mathbf{d})^2}{m}}$$ |
| 64 | + |
| 65 | +The first $r$ rows of $V^T$ are the combinations $\mathbf{v}_i^T\boldsymbol{\pi}$ the data determine; |
| 66 | +the rest are the unidentifiable directions. Moving $\boldsymbol{\pi}$ along them leaves the fit |
| 67 | +unchanged, and $\hat{\boldsymbol{\pi}}$ is orthogonal to them. |
| 68 | + |
| 69 | +Scaling a row and its right-hand side by $\sqrt{w_k}$ gives weighted least squares. Every step is |
| 70 | +orthogonal, so row order changes the result only by rounding. With no rows there is no result; a |
| 71 | +reset clears the state. |
| 72 | + |
| 73 | +## Complexity Analysis |
| 74 | + |
| 75 | +| Operation | Time | Space | Notes | |
| 76 | +|------------------------|-----------------|----------------|------------------------------------------------------| |
| 77 | +| Row, sparse (best) | $O(P)$ | $P^2 + P + 2$ | Zero entries skipped; one short rotation per nonzero | |
| 78 | +| Row, dense (worst) | $O(P^2)$ | $P^2 + P + 2$ | $P$ rotations on $P - j + 1$ pairs, $P$ square roots | |
| 79 | +| Solve | $O(P^3)$ | $O(P^2)$ extra | One $P \times P$ SVD, independent of $m$ | |
| 80 | +| Whole stream, $m$ rows | $O(mP^2 + P^3)$ | $P^2 + P + 2$ | Memory does not grow with $m$ | |
| 81 | + |
| 82 | +At $P = 60$ the state is about 14 KB in single precision. Accumulation can run at the sample rate; |
| 83 | +the solve is an offline or commissioning step. No heap, no recursion. |
| 84 | + |
| 85 | +## Step-by-Step Walkthrough |
| 86 | + |
| 87 | +**Full rank.** Rows $(1, x)$ at $x = (0, 1, 2, 3)$, $\boldsymbol{\tau} = (2.1, 4.9, 8.1, 10.9)$: |
| 88 | + |
| 89 | +| Rows | $R_{00}$ | $R_{01}$ | $R_{11}$ | $d_0$ | $d_1$ | $\rho$ | |
| 90 | +|------|----------|----------|----------|--------|--------|--------| |
| 91 | +| 1 | 1 | 0 | 0 | 2.1 | 0 | 0 | |
| 92 | +| 2 | 1.4142 | 0.7071 | 0.7071 | 4.9497 | 1.9799 | 0 | |
| 93 | +| 3 | 1.7321 | 1.7321 | 1.4142 | 8.7180 | 4.2426 | 0.0267 | |
| 94 | +| 4 | 2 | 3 | 2.2361 | 13 | 6.6188 | 0.0320 | |
| 95 | + |
| 96 | +Row 1 meets a zero pivot and is adopted; row 2 rotates by $45^\circ$ and fits two points exactly; |
| 97 | +from row 3 on, $\rho$ grows. With $\sigma \approx (4.1000, 1.0908)$, $r = 2$ and back-substitution |
| 98 | +gives $\hat{\pi}_1 = d_1 / R_{11} = 2.96$, $\hat{\pi}_0 = (13 - 3 \cdot 2.96)/2 = 2.06$, with |
| 99 | +RMS $= \sqrt{0.032 / 4} = 0.089443$. |
| 100 | + |
| 101 | +**Rank-deficient.** Rows $(1, x, x)$ with $\boldsymbol{\tau} = 2 + 3x$: two equal columns leave the |
| 102 | +third row of $R$ zero, $\rho = 0$, and $\sigma = (5.540482,\, 1.141517,\, 0)$, so $r = 2$. Only |
| 103 | +$\pi_1 + \pi_2 = 3$ is determined; the minimum-norm split is $\hat{\boldsymbol{\pi}} = (2, 1.5, 1.5)$. |
| 104 | +The last row of $V^T$ is the unidentifiable direction $(0, 1, -1)/\sqrt2$. |
| 105 | + |
| 106 | +Verified against double-precision normal equations on random 3-parameter problems of 5 to 205 rows: |
| 107 | +parameters agree to about $3 \times 10^{-6}$. |
| 108 | + |
| 109 | +## Pitfalls & Edge Cases |
| 110 | + |
| 111 | +- **Column scaling.** The rank test is relative to $\sigma_0$, so a column in large units hides the |
| 112 | + others. Scale columns to comparable magnitude before accumulating. |
| 113 | +- **Tolerance.** Too small admits noise-level directions and the solution grows like $1/\sigma_i$; too |
| 114 | + large discards real ones, visible as a higher RMS. About $10^{-5}$ is a single-precision start. |
| 115 | +- **RMS is not a noise estimate.** It divides by $m$, not $m - r$; scale its square by $m / (m - r)$ |
| 116 | + for an unbiased variance. |
| 117 | +- **Few or repeated rows.** Fewer rows than parameters is valid, with rank at most $m$; rows that |
| 118 | + repeat earlier ones leave the rank unchanged. |
| 119 | +- **Long runs.** $\rho$ only grows; once it is about $2^{24}$ times one term, new terms vanish and the |
| 120 | + RMS turns optimistic. One NaN row poisons the state until a reset. |
| 121 | + |
| 122 | +## Variants & Generalizations |
| 123 | + |
| 124 | +- **Ridge.** Feed $P$ extra rows $\sqrt{\lambda}\,\mathbf{e}_j$ with zero right-hand side. Replacing |
| 125 | + $1/\sigma_i$ by $\sigma_i / (\sigma_i^2 + \lambda)$ filters smoothly instead of cutting off. |
| 126 | +- **Forgetting.** Scaling $R$, $\mathbf{d}$ by $\sqrt{\lambda}$ and $\rho$ by $\lambda$ before each row |
| 127 | + gives exponentially weighted square-root recursive least squares. |
| 128 | +- **Merging logs.** Rotating the rows of one summary $[\,R \mid \mathbf{d}\,]$ into another, and adding |
| 129 | + the $\rho$ values and counts, combines separately accumulated segments. |
| 130 | + |
| 131 | +## Applications |
| 132 | + |
| 133 | +- Inertial-parameter identification of robot arms, rank-deficient by construction |
| 134 | +- Motor and drive parameter identification from long logs |
| 135 | +- Friction identification: if the axis only ever moves one way, the sign column duplicates the offset |
| 136 | + column and the unidentifiable direction says so |
| 137 | +- Sensor calibration from data streams too long to store |
| 138 | + |
| 139 | +## Connections to Other Algorithms |
| 140 | + |
| 141 | +- [Givens Rotation](../math/GivensRotation.md): the plane rotation that zeroes one entry per step. |
| 142 | +- [Singular Value Decomposition](../solvers/SingularValueDecomposition.md): reveals rank and directions. |
| 143 | +- [QR Decomposition](../solvers/QrDecomposition.md): the same factor, built row by row, no $Q$ kept. |
| 144 | +- [Recursive Least Squares](RecursiveLeastSquares.md): forgetting and a $P \times P$ covariance; here |
| 145 | + all rows weigh equally and rank is revealed. |
| 146 | +- [Linear Regression](LinearRegression.md): batch normal equations that store the samples. |
| 147 | +- [Total Least Squares](TotalLeastSquares.md): also SVD-based, but for noisy regressors, not rank. |
| 148 | +- [Friction Compensation](../controllers/FrictionCompensation.md): its Coulomb and viscous parameters |
| 149 | + are linear and identifiable this way. |
| 150 | + |
| 151 | +## References & Further Reading |
| 152 | + |
| 153 | +- G. H. Golub, C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins University Press, 2013 — |
| 154 | + §5.1.8–5.3 (Givens QR, updating) and §5.5 (rank-deficient least squares). |
| 155 | +- Å. Björck, *Numerical Methods for Least Squares Problems*, SIAM, 1996. |
0 commit comments