3
$\begingroup$

Suppose I have a system:

$$ Ax = b $$

where $A$ is a $m$ by $n$ matrix which is less than full rank (neither full column nor row rank). In my particular case $m<n$.

I'd like a combination of a minimum norm solution (for elements of $x$ which are not determined and a least squares solution for those that are determined or over determined).

I know I can approximate this with Tikhonov regularization:

$$ \mathop{\text{minimize}}\limits_x \|Ax-b\|^2 + \alpha\|x\|^2. $$

But this is not ideal as the regularizer tugs on all elements of $x$ not just those left undetermined. If I make $\alpha$ too large then my $x$ at determined or over determined elements will not be optimal (w.r.t. $\|Ax-b\|^2$). If $\alpha$ is too small, I imagine I run into conditioning issues.

Is there no way to formulate this correctly using some sort null space/QR decomposition?

Update: See accepted answer, but for completeness here's a full recipe: \begin{align} Ax &= b\\ AE &= QR &\text{rank-revealing QR decomposition of A}\\ AE &= Q\left(\begin{array}{c}I\\0\end{array}\right)\left(\begin{array}{cc}R_1&R_2\end{array}\right) & \text{only keep $r$ (rank) many rows of $R$}\\ \left(\begin{array}{c}R_1^T\\R_2^T\end{array}\right)F &= VT_\text{full} &\text{rank-revealing QR decomposition of $\left(\begin{array}{c}R_1^T\\R_2^T\end{array}\right)$} \\ \left(\begin{array}{c}R_1^T\\R_2^T\end{array}\right)F &= V\left(\begin{array}{c}I\\0\end{array}\right)T & \text{only keep $r$ many rows of $T$}\\ AE &= Q\left(\begin{array}{c}I\\0\end{array}\right)FT^T\left(\begin{array}{cc}I&0\end{array}\right)V^T \\ x &= A^+b\\ x &= EV\left(\begin{array}{c}I\\0\end{array}\right)T^{-T}F^T\left(\begin{array}{cc}I&0\end{array}\right)Q^Tb \end{align}

$\endgroup$
  • $\begingroup$ Preemptive strike: Here's an example of an under and over determined system: $$\left(\begin{array}{ccc}1&0&0\\1&0&0\\0&1&0\end{array}\right)\left(\begin{array}x_1\\x_2\\x_3\end{array}\right) = \left(\begin{array}{c}4\\5\\6\end{array}\right)$$ $\endgroup$ – Alec Jacobson May 6 '14 at 22:12
  • $\begingroup$ Does SVD not give you what you want? $\endgroup$ – Igor Khavkine May 7 '14 at 8:54
  • $\begingroup$ My matrices are very large and sparse, so complete SVD is not really an option (I assume). $\endgroup$ – Alec Jacobson May 7 '14 at 14:32
  • $\begingroup$ The SVD will give the answer, but it's not always the minimal amount of work to get a desired datum. They are implementations for sparse cases. stats.stackexchange.com/questions/41259/… I expect that using one of them is probably better than trying to implement one's own minimalist approach (with QR rather than full SVD). $\endgroup$ – Conder May 8 '14 at 4:06
  • $\begingroup$ I'm happy to implement my own approach with QR if it avoids a dense matrix of singular vectors... $\endgroup$ – Alec Jacobson May 8 '14 at 12:00
3
$\begingroup$

The Moore-Penrose pseudoinverse is probably what you're looking for. The pseudoinverse solution $A^+b$ is the smallest norm $x$ such that $\|Ax-b\|_2$ is minimized. It can be computed using QR decomposition although you have to use rank-revealing QR when $A$ does not have full column rank. A more complete explanation is given here.

$\endgroup$
  • $\begingroup$ So if I have $x = A^+b$ and $A^+=(A^TA)^{-1}A^T$ and the reduced QR decomposition of $A$ is $A=Q_1R_1$, then I'll want to solve using $R_1^TR_1x = A^Tb$ apply back-,forward- substitution with $R_1$? But if $A$ doesn't have full column rank, won't $R_1$ be rank deficient? Or what do you mean by "reduced QR"? $\endgroup$ – Alec Jacobson May 6 '14 at 22:55
  • $\begingroup$ @mangledorf Sorry, I think the right term is rank-revealing QR decomposition. This presentation explains the method and its application to this problem. $\endgroup$ – Alex Becker May 6 '14 at 23:23
  • $\begingroup$ I think you do want the pseudoinverse. An alternative way to compute, take random (orthogonal) right multiplicators $T$ until from QR you decompose $AT=QR$ where $R$ is of the special form $\Bigl({U|B\atop Z}\Bigr)$ where $Z$ is numericaly zero (possibly empty) and $U$ is invertible and upper triangular. When $B$ is empty (full column rank), you are done (take $U^{-1}$, suitably zero-extended). Else you can recurse (only once) to compute the pseudo inverse of $(U|B)$ by taking the transpose before and after, and if the same parameters are used for zero-detection, this has full column rank. $\endgroup$ – Conder May 7 '14 at 4:25
  • $\begingroup$ @Conder can you elaborate? Is this the same rank-revealing QR that AlexBecker writes of? How can $B$ be empty and full column rank? What does suitably zero-extended mean? Does $B$ participate in the final solve? $\endgroup$ – Alec Jacobson May 8 '14 at 2:03

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.