# Modifying U=mxn SVD Algorithm to U=mxm Algorithm

I have painstakingly ported this Python source "svd.py" to C++. I confirm it works for the example it comes with. While testing another example (this one, from Wikipedia), the assert statement trips because $m < n$.

I find this somewhat of a curious constraint. It seems to fail mainly because the matrix $U$ (named $u$ in the source) is initialized as an $m \times n$ matrix (as opposed to the $m \times m$ matrix it is in the definition of a SVD). It's pretty simple to initialize $U$ instead to an $m \times m$ matrix (and, say, pad with zeroes), but this causes some of the code to break for some other examples.

Why would $U$ be an $m \times n$ matrix, and how can I fix the algorithm to work on all sizes of matrices?

EDIT: If the SVD is $M=U_1\Sigma_1 V_1^*$, I came up with the idea of computing the SVD of $M^*$ to get $M^*=V_1\Sigma_1^*U_1^*=U_2\Sigma_2 V_2^*$. From here, I use that substitution to get the desired $U_1$, $\Sigma_1$, $V_1$ out. This produces a result that multiplies back through correctly, but is it a proper SVD?

• Why are you porting that particular svd algorithm when there are tons of well-implemented svd's in C++ already? – Batman Aug 29 '14 at 20:32
• @Batman it needs to integrate with some existing C++ sources in a very particular way, and external libraries aren't an option. I wasn't able to find a compact enough SVD algorithm to adapt to my system. – imallett Aug 29 '14 at 20:34
• @IanMallett Since there are plenty of implementations, you can either use a library (I don't see why you can't, but I don't know your constraints), either extract the important source code and integrate it in your system. There are good implementations in C (for example from gsl), Fortran (for example from Lapack), and probably C++ as well. So I don't really see the problem. Usually, the simple way to go is just to find a Lapack compiled for your OS/compiler, and integrate it, within any language (even Excel's VBA). – Stop hurting Monica Aug 29 '14 at 21:05
• Ad EDIT: Yes, it is a correct "economy" SVD of a "wide" matrix. It depends on the application whether or not you need a "full" SVD. – Algebraic Pavel Sep 8 '14 at 14:08

1. Sounds like the routine is built for overdetermined systems. The standard theoretical trick is to work with $\mathbf{A}^{*}$ in lieu of $\mathbf{A}$. An underdetermined system is going to contend with infinite solution space.
2. Given a matrix of rank $\rho$, $$\mathbf{A} \in \mathbb{C}^{m\times n}_{\rho},$$ the singular value decomposition is \begin{align} \mathbf{A} &= \mathbf{U} \, \Sigma \, \mathbf{V}^{*} \\ % &= % U \left[ \begin{array}{cc} \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}} \end{array} \right] % Sigma \left[ \begin{array}{cccccc} \sigma_{1} & 0 & \dots & & & \dots & 0 \\ 0 & \sigma_{2} \\ \vdots && \ddots \\ & & & \sigma_{k} \\ & & & & 0 & \\ \vdots &&&&&\ddots \\ 0 & & & & & & 0 \\ \end{array} \right] % V \left[ \begin{array}{c} \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \\ \color{red}{\mathbf{V}_{\mathcal{N}}}^{*} \end{array} \right] \\ % & = % U \left[ \begin{array}{cccccccc} \color{blue}{u_{1}} & \dots & \color{blue}{u_{k}} & \color{red}{u_{k+1}} & \dots & \color{red}{u_{n}} \end{array} \right] % Sigma \left[ \begin{array}{cc} \mathbf{S}_{k\times k} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{array} \right] % V \left[ \begin{array}{c} \color{blue}{v_{1}^{*}} \\ \vdots \\ \color{blue}{v_{k}^{*}} \\ \color{red}{v_{k+1}^{*}} \\ \vdots \\ \color{red}{v_{n}^{*}} \end{array} \right] % \end{align} Your problem may correspond to $\rho = n$.
3. The least squares solution highlights the uniqueness issue. Given the solution vector $x\in\mathbb{C}^{n}$, and the data vectors $b \in \mathbb{C}^{m}$, and the restriction that $b\notin \color{red}{\mathcal{N}\left( \mathbf{A}\right)}$, the least squares solution to $\mathbf{A}x = b$ is $$x_{LS} = \color{blue}{\mathbf{A}^{\dagger} b} + \color{red}{\left( \mathbf{I}_{n} - \mathbf{A}^{\dagger}\mathbf{A} \right) y}, \quad y\in\mathbb{C}^{n}$$ If the nullspace $\color{red}{\mathcal{N}\left( \mathbf{A}^{*}\right)}$ is trivial, then $\mathbf{A}^{\dagger}\mathbf{A} =\mathbf{I}_n$ and the solution is unique.
4. So yes, we can build the Hermitian conjugate with SVD components. Using the thin form $$\mathbf{A} = \color{blue}{\mathbf{U}_{\mathcal{R}}} \mathbf{S} \, \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} % \qquad \Rightarrow \qquad % \mathbf{A}^{*} = \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \mathbf{S} \, \color{blue}{\mathbf{U}_{\mathcal{R}}}$$