How can I solve on a computer a large projection problem with redundant constraints?

Question

This question is the essence of this one. After we remove all the cruft, we can recast it as follows:

Problem: Given $b \in \mathbb{R}^n$, $C\in \mathbb{R}^{n\times m}$, and $g\in \mathrm{Range}(C^\mathrm{T})\subset \mathbb{R}^m$. Find $x\in \mathbb{R}^n$ so that $\| x-b \|$ is minimized subject to the constraint $C^\mathrm{T} x = g$.

If $C$ had full column rank, this is straight forward problem, but I don't have that, and instead I know that $g\in \mathrm{Range}(C^\mathrm{T})$. In my case $m$ is large but small compared to $n$. How can I do this?

Answer 1

We have the following linear system in $\mathrm x \in \mathbb R^n$

$$\rm A x = b \tag{1}$$

where $\mathrm A \in \mathbb R^{m \times n}$ is fat (i.e., $n > m$) and rank-deficient (i.e., $r := \mbox{rank} (\mathrm A) < m$), and $\mathrm b \in \mathbb R^m$ is in the column space of $\mathrm A$, i.e., the linear system is consistent. We would like to find the solution of $(1)$ that is closest to a given $\mathrm x_0 \in \mathbb R^n$. We thus have the following (convex) quadratic program

$$\begin{array}{ll} \text{minimize} & \| \mathrm x - \mathrm x_0 \|_2^2\\ \text{subject to} & \mathrm A \mathrm x = \mathrm b\end{array}$$

Let's look for the solution of this QP using the singular value decomposition (SVD) of $\rm A$.

---

Solution via unconstrained minimization

Let the SVD of $\rm A$ be

$$\mathrm A = \mathrm U \Sigma \mathrm V^{\top} = \begin{bmatrix} \mathrm U_1 & \mathrm U_2\end{bmatrix} \begin{bmatrix} \Sigma_1 & \mathrm O\\ \mathrm O & \mathrm O\end{bmatrix} \begin{bmatrix} \mathrm V_1^{\top}\\ \mathrm V_2^{\top}\end{bmatrix}$$

Hence, $\rm A x = b$ can be written as $\rm \mathrm U \Sigma \mathrm V^{\top} x = b$, or, $\rm \Sigma \mathrm V^{\top} x = \mathrm U^{\top} b$. Let $\rm y := V^{\top} x$. We have

$$\begin{bmatrix} \Sigma_1 & \mathrm O\\ \mathrm O & \mathrm O\end{bmatrix} \begin{bmatrix} \mathrm y_1\\ \mathrm y_2\end{bmatrix} = \begin{bmatrix} \mathrm U_1^{\top} \mathrm b\\ \mathrm U_2^{\top} \mathrm b\end{bmatrix}$$

Since $\rm b$ is in the column space of $\rm A$, it is orthogonal to the left null space of $\rm A$, i.e., $\mathrm U_2^{\top} \mathrm b = 0_{m-r}$.

$$\begin{bmatrix} \Sigma_1 & \mathrm O\\ \mathrm O & \mathrm O\end{bmatrix} \begin{bmatrix} \mathrm y_1\\ \mathrm y_2\end{bmatrix} = \begin{bmatrix} \mathrm U_1^{\top} \mathrm b\\ 0_{m-r}\end{bmatrix}$$

Note that $\mathrm y_2$ is free. The solution set is the $(m - r)$-dimensional affine space

$$\mathrm y \in \left\{ \begin{bmatrix} \Sigma_1^{-1} \mathrm U_1^{\top} \mathrm b\\ \eta\end{bmatrix} : \eta \in \mathbb R^{m-r} \right\}$$

Since $\rm x = V y$,

$$\mathrm x \in \left\{ \color{blue}{\mathrm V_1 \Sigma_1^{-1} \mathrm U_1^{\top} \mathrm b + \mathrm V_2 \eta} : \eta \in \mathbb R^{m-r} \right\}$$

where

$$\mathrm x_{\text{LN}} := \color{blue}{\mathrm V_1 \Sigma_1^{-1} \mathrm U_1^{\top} \mathrm b}$$

is the least-norm solution of $\rm A x = b$, as

$$\| \mathrm x \|_2^2 = \mathrm x^{\top} \mathrm x = \mathrm y^{\top} \underbrace{\mathrm V^{\top} \mathrm V}_{= \mathrm I_n} \, \mathrm y = \mathrm y^{\top} \mathrm y = \mathrm b^{\top} \mathrm U_1 \Sigma_1^{-2} \mathrm U_1^{\top} \mathrm b + \| \eta \|_2^2$$

is minimized when $\eta = 0_{m-r}$.

We now have an unconstrained quadratic program in $\eta \in \mathbb R^{m-r}$

$$\begin{array}{ll} \text{minimize} & \| \mathrm V_2 \eta - (\mathrm x_0 - \mathrm x_{\text{LN}}) \|_2^2\end{array}$$

Since $\rm V V^\top = V_1 V_1^\top + V_2 V_2^\top = I_n$, where $\rm V_1 V_1^\top$ and $\rm V_2 V_2^\top$ are the projection matrices that project onto the row space and the (right) null space of $\rm A$, respectively, the objective function can be written as follows

$$\begin{array}{rl} \| \mathrm V_2 \eta - (\mathrm x_0 - \mathrm x_{\text{LN}}) \|_2^2 &= \| \mathrm V_2 \eta - \mathrm V_1 \mathrm V_1^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) - \mathrm V_2 \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) \|_2^2\\ &= \| \mathrm V_1 \mathrm V_1^\top ( \mathrm x_{\text{LN}} - \mathrm x_0 ) + \mathrm V_2 \left( \eta - \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) \right) \|_2^2\\ &= \left\| \begin{bmatrix} \mathrm V_1 & \mathrm V_2\end{bmatrix} \begin{bmatrix} \mathrm V_1^\top ( \mathrm x_{\text{LN}} - \mathrm x_0 )\\ \eta - \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) \end{bmatrix} \right\|_2^2\\ &= \left\| \begin{bmatrix} \mathrm V_1^\top ( \mathrm x_{\text{LN}} - \mathrm x_0 )\\ \eta - \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) \end{bmatrix} \right\|_2^2\\ &= \| \mathrm V_1^\top ( \mathrm x_{\text{LN}} - \mathrm x_0 ) \|_2^2 + \| \eta - \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) \|_2^2\end{array}$$

which is minimized at

$$\eta^* := \mathrm V_2^\top (\mathrm x_0 - \mathrm x_{\text{LN}}) = \mathrm V_2^\top \mathrm x_0 - \underbrace{\mathrm V_2^\top \mathrm V_1}_{= \mathrm O_{(n-r) \times r}} \Sigma_1^{-1} \mathrm U_1^{\top} \mathrm b = \mathrm V_2^\top \mathrm x_0$$

and, thus, the solution of the quadratic program is

$$\boxed{\mathrm x^* := \mathrm x_{\text{LN}} + \mathrm V_2 \eta^* = \color{blue}{\mathrm V_1 \Sigma_1^{-1} \mathrm U_1^{\top} \mathrm b + \mathrm V_2 \mathrm V_2^\top \mathrm x_0}}$$