1
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ Why are you porting that particular svd algorithm when there are tons of well-implemented svd's in C++ already? $\endgroup$ – Batman Aug 29 '14 at 20:32
  • $\begingroup$ @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. $\endgroup$ – imallett Aug 29 '14 at 20:34
  • $\begingroup$ @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). $\endgroup$ – Stop hurting Monica Aug 29 '14 at 21:05
  • $\begingroup$ 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. $\endgroup$ – Algebraic Pavel Sep 8 '14 at 14:08
0
$\begingroup$
  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}}} $$

$\endgroup$

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.