Set of solutions to the least-squares problem via QR decomposition

The set mathbf{S} of solutions to the least-squares problem

 min_x : |Ax - y|_2^2,

where A in mathbf{R}^{m times n}, and y in mathbf{R}^m are given, can be expressed in terms of the full QR decomposition of A:

 AP = QR , ;; R = left( begin{array}{cc} R_{11} & R_{12}  0 & 0 end{array} right),

where R_{11} in mathbf{R}^{r times r} is upper triangular and invertible, R_{12} in mathbf{R}^{r times (n-r)}, P is a permutation matrix, and Q is m times m and orthogonal.

Precisely we have mathbf{S} = { x_0 + Nz ::: z in mathbf{R}^{n-r} }, with N a matrix whose columns span the nullspace of A:

 x_0 = P left(begin{array}{c} R_{11}^{-1}y_1  0 end{array}right), ;; N = left(begin{array}{c} R_{11}^{-1}R_{12}  I end{array}right) in mathbf{R}^{n times (n-r)}.

Proof: Since Q and P are orthogonal, we have, with bar{x} := P^Tx, bar{y} := Q^Ty:

 Ax-y=APP^Tx-y = QR bar{x} - y = Q(R bar{x} - Q^Ty) = Q(R bar{x} - bar{y})

Exploiting the fact that Q leaves Euclidean norms invariant, we express the original least-squares problem in the equivalent form:

 min_{bar{x}} : |R bar{x} - bar{y}|_2.

Once the above is solved, and bar{x} is found, we recover the original variable x with x = Pbar{x}.

Now let us decompose bar{x} and bar{y} in a manner consistent with the block structure of R: bar{x} = (x_1,x_2), bar{y} = (y_1,y_2), with x_1,y_1 two r-vectors. Then

 R bar{x} - bar{y} = left( begin{array}{c} R_{11}x_1+R_{12}x_2-y_1  -y_2 end{array} right),

which leads to the following expression for the objective function:

 |R bar{x} - bar{y}|_2^2 = |R_{11}x_1+R_{12}x_2-y_1|_2^2+|y_1|_2^2.

The optimal choice for the variables x_1,x_2 is to make the first term zero, which is achievable with

 x_1 = R_{11}^{-1}(y_1-R_{12}x_2),

where x_2 is free, and describes the ambiguity in the solution. The optimal residual is |y_1|_2^2.

We are essentially done: with x_2 = z, we can write

 P^Tx = bar{x} = left(begin{array}{c} x_1  x_2 end{array}right), left(begin{array}{c} R_{11}^{-1}y_1  0 end{array}right) +  left(begin{array}{c} R_{11}^{-1}R_{12}  I end{array}right)z,

that is: x = Pbar{x} = x_0 + Nz, with

 x_0 = P left(begin{array}{c} R_{11}^{-1}y_1  0 end{array}right), ;; N = Pleft(begin{array}{c} R_{11}^{-1}R_{12}  I end{array}right) in mathbf{R}^{n times (n-r)}.