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.
