I.6.5.1: Yates Augmentation

A linear least squares problem is balanced if the design matrix $A$ is orthogonal. In the balanced case the problem of minimizing $$f(x)=(b-Ax)'(b-Ax)\tag{1}$$ is easily solved, with solution $x=A'b$. Computing $x$ does not involve any matrix inversion. In the case of a balanced factorial design it simply involves computing means of rows, columns, slices, and so on.

If some elements of $b$ are missing then we can partition $A$ and $b$ into a missing and non-missing part as in $$A=\begin{bmatrix}A_1\\A_2\end{bmatrix}\qquad b=\begin{bmatrix}b_1\\b_2\end{bmatrix},$$ with $A_1$ and $b_1$ the non-missing part, and minimize $$f(x)=(b_1-A_1x)'(b_1-A_1x).$$ Now $A_1'A_1=I-A_2'A_2<I$, and the optimal $x$ can no longer be computed with a simple matrix multiplication.

If one want to avoid matrix inversion, then we can use the basic approach suggested by Yates [1933]. We define the augmentation $g$ by $$g(x,z)=(b_1-A_1x)'(b_1-A_1x)+(z-A_2x)'(z-A_2x).$$ Because $f(x)=\min_z g(x,z)$ this leads to an easy block relaxation algorithm. $$\begin{align*} z^{(k+1)}&=A_2x^{(k)},\\ x^{(k+1)}&=A_1'b_1+A_2'z^{(k+1)}. \end{align*}$$ This gives $$\begin{align*} z^{(k+1)}&=A_2A_1'b_1+A_2A_2'z^{(k)},\\ x^{(k+1)}&=A_1'b_1+A_2'A_2x^{(k)}. \end{align*}$$ As we have shown in section blockrelaxation:twoblockleastsquares this implies an iteration radius equal to the largest eigenvalue of $A_2'A_2$.

Note that Yates augmentation can be used to transform any linear least squares problem to a balanced problem, even if there are no missing data. In minimizing $f(x)$ of $(1)$ we first check if $A'A\lesssim I.$ If this is the case, there is no need to normalize. If we do not have $A'A\lesssim I$ we start by dividing $A$ by $\tau(A)=\sqrt{\mathbf{tr}\ A'A}$. This new normalized $A$, say $\tilde A$, now satisfies $\tilde A'\tilde A\lesssim I$. Then find any $G$ such that $\tilde A'\tilde A+G'G=I$, and iterate according to $$\begin{align*} z^{(k+1)}&=Gx^{(k)},\\ x^{(k+1)}&=\tilde A'b+G'z^{(k+1)}, \end{align*}$$ which amounts to $$x^{(k+1)}=\tilde A'b+(I-\tilde A'\tilde A)x^{(k)}.$$ The iteration radius is $$\kappa=1-\frac{\lambda_{\text{min}}(A'A)}{\mathbf{tr}\ A'A}.$$ We can do better if we compute $\tilde A$ by dividing $A$ by its trace norm, i.e. its largest singular value. Then $$\kappa=1-\frac{\lambda_{\text{min}}(A'A)}{\lambda_{\text{max}}(A'A)}.$$ There is R code for Yates augmentation in the file yates.R.