Computing the SVD factorization on C++ (using the proof of the existence of the SVD factorization)

I am doing a C++ program that computes the SVD factorization of a real matrix A without using any known library of algebra that contains the implementation. In addition, QR descomposition is not allowed in any step (because one of my goals is to compare QR and SVD descomposition in a least squares problem).

I know there are a lot of methods, with well tecniques, but my method is simple, far from being the fastest and better one, but it is based on a proof of the theorem about SVD factorization that I have done at University. I know that it is not the best way (I accept ideas and proposals if you want to tell me something better). The proof of the theorem (and the diagram of my algorithm) is:

Teorem:

Let $A \in M_{m,n} (\mathbb{R})$, with $m \geq n$. Then exists the $A$ = $U \cdot \Sigma \cdot V^T $, with these properties:

$U\in M_{m,n} (\mathbb{R})$ satisfies $U \cdot U^T = I_n$ (orthogonal matrix)

$\Sigma \in M_{n,n} (\mathbb{R})$ is a diagonal matrix with $\Sigma =$ diag$(\sigma_1, \sigma_2,..., \sigma_n)$

$V\in M_{n,n} (\mathbb{R})$ satisfies $V \cdot V^T = I_n$ (orthogonal matrix)

Proof:

Construct $S = A^T \cdot A$, so $S \in M_{n,n} (\mathbb{R})$. We have that $S$ is a symmetric and positive semi-definite matrix (you can check that as an exercise, but if you want to know why it is true ask me and then I'll write the proof of that).

According to that, we can find an orthonormal basis of $\mathbb{R}^n$ of eigenvectors, $\mathbb{R}^n = [u_1,u_2,...,u_n]$, with $<u_i,u_j>=d_{ij}$, with $d_{ij}=0$ if $j\neq i$ and $d_{ij}=1$ if $j > = i$. Let be $\lambda_1, \lambda_2, ..., \lambda_n$ the respective eigenvalues. Let be $1 \leq r \leq n$, with $\lambda_1, \lambda_2,..., > \lambda_r \geq 0$, and $\lambda_{r+1}, \lambda_{r+2}, \lambda_{n} = 0$ (remember that $S$ is a positive semi-definite matrix, so eigenvalues can't be negative).

Then, $\sigma_i = \sqrt{\lambda_i}$, so our $\Sigma = $diag$(\sigma_1, > \sigma_2,..., \sigma_n)$ it is now constructed. In addition, we will say that $V=(v_1,v_2,...,v_n)$, and our matrix $V$ is now constructed too. At the end we will see that all these definitions make sense and the descomposition works. Finally, let's find the matrix U.

We introduce the vectors $u_i = \frac{1}{\sigma_i} \cdot A \cdot v_i$, for $1 \leq i \leq r$. These are unitary and two-by-two orthonormal vectors (you can check that as an exercise, if you want to know the proof tell me). Now we choose $u_{r+1}, u_{r+2},..., u_{n} \in > \mathbb{R}^m$ that makes {$u_1, u_2,...,u_r,u_{r+1},u_{r+2},...,u_n$} an orthonomal basis of $\mathbb{R}^m$. And we finally have our last matrix $U = (u_1, u_2, ... , u_n)$. Let's check if the factorization is correct.

It is easy to see that $A \cdot V = U \cdot \Sigma$ : we have that $ A > \cdot v_j = \sigma_j \cdot \frac{1}{\sigma_j} \cdot A \cdot v_j = > \sigma_j \cdot u_j$ for $1 \leq j \leq r$ and $A \cdot v_j = 0$ for $r+1 \leq j \leq n$, so the equality it's true. Then, as the last step, we see that $V$ is an orthogonal matrix (because the column vectors make an orthonormal basis), so $V^{-1} = V^T$ and $A \cdot V > \cdot V^T = U \cdot \Sigma \cdot V^T \iff A = U \cdot \Sigma \cdot > V^T$, as we wanted to proof $\blacksquare$

First, I want to use a method like the power method to find the eigenvalues and eigenvectors of $S$, knowing that $S$ is a symmetric positive semi-definite matrix. The problem of using the power method is that it does not use the properties of $S$, and it finds only one eigenvector for one eigenvalue, and I want to find all the eigenvectors linked with a given eigenvalue... in fact, I want to know a basis of the eigenspace corresponding to a given eigenvalue.

First question: Do you know a good algorithm (not too difficult) to find all the eigenvectors for each eigenvalue? (not only one, if more than one eigenvector have the same eigenvalue). Does this method give me the orthonormal basis of eigenvectors? I can't use the QR algorithm (I currently saw an algorithm to find the eigenspace of an eigenvalue using QR factorization).

If I have this algorithm, it's very easy to find $V$, $\Sigma$ and the $u_1,u_2,...,u_r$ of the proof. The last step is to find $u_{r+1},...,u_{n}$ vectors that make {$u_1,...,u_n$} an orthonormal basis.

Second (and last) question: Do you know a good algorithm to, given some two-by-two orthonormal vectors, complete them to a basis of the corresponding vectorial space? I don't know if it's useful to find linearly independent vectors and use the Gram-Schmidt algorithm to make an orthonormal basis: in fact, if I do this I'm using a half of a part of the QR implementation (that can be done with Gram-Schmidt), so if it is an alternative to this I will be grateful to read it.

Lots of thanks and sorry for my English.

From Spain, have a nice day.

Solution 1:

First of all, please read a standard textbook on matrix computations. The best is the "Matrix Computations" by Gene Golub and Charles Van Loan [GVL]. [BTW, the personalised car number plate of (late) Professor Golub was DR SVD]. Once you have learned what is in GVL, proceed to the LAPACK project algorithms. see E. Anderson et al. "LAPACK Users' Guide", SIAM, 1999. There is an updated online version at www.netlib.org/lapack/lug/.

I assume that "$.$" means matrix multiplication (linear algebraist usually do not use the dot for multiplication; so I avoid it).

$UU^t$ is not equal to $I_n$ unless $m=n$. However, $U^t U = I_n$.

Except in special circumstances, $A A^t$ (or $A^tA$) is not formed by numerical analyst in floating-point arithmetic. See Golub and Van Loan.

It is difficult to summarise all the methods available for computing the SVD. However, after reading GVL, please google SVD together with "inverse iterations" and "power iterations" (usually to find few of the singular vectors once the singular values are known), "QR algorithms", "qd algorithms", "divide and conquer".

The existence of the SVD was first proved in the 19th century with more formal results coming in the early to the middle 20th century. Unfortunately, there is no need to reinvent the wheel (it has been done several times already).