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?

  • 2
    $\begingroup$ What's the big deal? Throw it in a Quadratic Programming solver, or formulate it as a Second Order Cone Problem (SOCP) and provide it to an SOCP solver. Perhaps the solver's pre-solver will transform it to a full rank problem, but even if not, it should still be solvable. Redundant (i.e., non-full rank) linear constraints are common and routine - no biggy. $\endgroup$ – Mark L. Stone Jun 28 '17 at 16:34
  • $\begingroup$ Thanks for the suggestion, I have never heard of a Second Order Cone Problem. You are right, this is exactly just a quadratic programming question.I guess my question is, then, about how would these quadratic programming solvers work in the presence of non-full rank constraints. I am hoping for something simple that I could implement myself. $\endgroup$ – fred Jun 28 '17 at 17:33
  • 1
    $\begingroup$ Shouldn't be a problem for a well-written solver. Could be a problem for a not so well-written solver, such as you might produce yourself if you don't know what you're doing. I suggest you use an existing solver. $\endgroup$ – Mark L. Stone Jun 28 '17 at 18:15

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\}$$


$$\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}}$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.