Power Iteration

Recall the singular value decomposition (SVD) of a matrix A = UDV^{T} = \sum_{i=1}^{r}\sigma_i \boldsymbol{u^{(i)}}\boldsymbol{v^{(i)^{T}}}.

We showed that V_{k} = span\{\boldsymbol{v^{(1)}}, ..., \boldsymbol{v^{(k)}}\} is the best-fit subspace for A, so

V_{k} = \arg\min_{V: dim(V) \leq k}\sum_{i=1}^{m}d(A_{(i)}, V)^2

We have that A_{(i)}V_{k}V_{k}^{T} is the projection of A_{(i)} to V_{k}, and so AV_{k}V_{k}^{T} is the projection of A to V_{k}.

Let A_k be defined as AV_{k}V_{k}^{T} = \sum_{i=1}^{k}\sigma_i \boldsymbol{u^{(i)}}\boldsymbol{v^{(i)^{T}}}.

Theorem 1.
\|A-A_{k}\|_{F} = \min_{B: rank(B) \leq k}\|A-B\|_{F}
AND
\|A-A_{k}\|_{2} = \max_{\|\boldsymbol{v}\| = 1}\|(A-A_{k}) \cdot \boldsymbol{v}\| = \sigma_{k+1} = \min_{B: rank(B) \leq k}\|A-B\|_{2}

Proof.
We proved part (1) last lecture. We will prove part (2) here. Note that

A-A_{k} = \sum_{i=k+1}^{r}\boldsymbol{u^{(i)}}\boldsymbol{v^{(i)^{T}}}
and
\|A-A_{k}\|_{2} = \sigma_{k+1}.

We need to show that for any B: rank(B)=k,
B=WZ^{T} where W, Z have \leq k columns and \|A-B\|_{2} \geq \sigma_{k+1}.

Since we know that Z has only k columns, \exists \boldsymbol{v} \in span\{\boldsymbol{v^{(1)}},...,\boldsymbol{v^{(k+1)}}\} such that
Z^{T}\boldsymbol{v} = 0,
\boldsymbol{v} = \sum_{i=1}^{k+1}\alpha_{i}\boldsymbol{v^{(i)}}, and
\|\boldsymbol{v}\|=1.

(A-B) \cdot \boldsymbol{v} = A \cdot \boldsymbol{v} = \sum_{i=1}^{k+1}\sigma_{i}\alpha_{i}\boldsymbol{u^{(i)}}
\|(A-B) \cdot \boldsymbol{v}\|^{2} = \sum_{i=1}^{k+1}\sigma_{i}^{2}\alpha_{i}^{2} \geq \sigma_{k+1}^{2}. \square

We’ve seen that the SVD is useful in more than one way. Now we will look at how to compute it, a central question in numerical analysis. A simple method, power iteration, is presented here:

Set \boldsymbol{x} \leftarrow e_i for i=1,2,\ldots,n
Repeat t times:{
\boldsymbol{x} \leftarrow A^{T}A\boldsymbol{x}
\boldsymbol{x} \leftarrow \frac{\boldsymbol{x}}{\|\boldsymbol{x}\|}
}

Theorem 2.
Assume (x\cdot v^{(1)})^2 \ge \frac{1}{n}. After t = \frac{\log n}{2\varepsilon} iterations, \boldsymbol{x} satisfies \|A\boldsymbol{x}\| \geq (1-\varepsilon) \sigma_{1}.

Proof.
\boldsymbol{x} = \sum_{i=1}^{r}\alpha_{i}\boldsymbol{v^{(i)}}
A^{T}A = \sum_{i=1}^{r}\sigma_{i}^{2}\boldsymbol{v^{(i)}}\boldsymbol{v^{(i)^{T}}}
(A^{T}A)^l = \sum_{i=1}^{r}\sigma_{i}^{2l}\boldsymbol{v^{(i)}}\boldsymbol{v^{(i)^{T}}}

\frac{\boldsymbol{x}^{T}(A^{T}A)^{l}A^{T}A(A^{T}A)^{l}\boldsymbol{x}}{\boldsymbol{x}^{T}(A^{T}A)^{l}(A^{T}A)^{l}\boldsymbol{x}} = \frac{\sum_{i=1}^{r}\sigma_{i}^{4l+2}\alpha_{i}^{2}}{\sum_{i=1}^{r}\sigma_{i}^{4l}\alpha_{i}^{2}}

Since \sigma_{1} \geq \sigma_{2} \geq ..., the weight on higher \sigma_{i} is increasing.

\sum \sigma_{i}^{4l}\alpha_{i}^{2} = \sum \sigma_{i}^{4l}(\alpha_{i}^{2})^{\frac{2l}{2l+1}}(\alpha_{i}^{2})^{\frac{1}{2l+1}}
\leq (\sum \sigma_{i}^{4l+2}\alpha_{i}^{2})^{\frac{2l}{2l+1}}(\sum \alpha_{i}^{2})^{\frac{1}{2l+1}}
= (\sum \sigma_{i}^{4l+2}\alpha_{i}^{2})^{\frac{2l}{2l+1}}

\therefore \frac{\boldsymbol{x^{l^{T}}}A^{T}A\boldsymbol{x^{l}}}{\boldsymbol{x^{l^{T}}}\boldsymbol{x^{l}}} \geq (\sum \sigma_{i}^{4l+2}\alpha_{i}^{2})^{\frac{1}{2l+1}} \geq \alpha_{1}^{\frac{2}{2l+1}}\sigma_{1}^{2}.

Therefore, after l iterations,
\frac{\|A\boldsymbol{x^{l}}\|^{2}}{\|\boldsymbol{x^{l}}\|^{2}} \geq (\frac{1}{n})^{\frac{1}{2l+1}}\sigma_{1}^{2}
= e^{-\frac{\ln n}{2l+1}}\sigma_{1}^{2}

Setting 2l + 1 = \frac{\ln n}{\varepsilon} \geq e^{-\varepsilon}\sigma_{1}^{2} \geq (1-\varepsilon)\sigma_{1}^{2}. \square

We note that the choice of the initial x in the algorithm, one of the axis vectors e_i, must satisfy the assumption that (e_i \cdot v^{(1)})^2 \ge \frac{1}{n}. Instead of trying all n basis vectors, we could also start with a random unit vector.

Claim.
\mathbb{E}(\alpha_{1}^{2}) = \frac{1}{n}

Proof.
\sum x_{i}^{2} = 1 symmetric.

So, \mathbb{E}(x_{1}^{2}) = \mathbb{E}(x_{2}^{2}) = ... = \mathbb{E}(x_{n}^{2}).

\implies n\mathbb{E}(x_{1}^{2}) = 1 \implies \mathbb{E}(x_{1}^{2}) = \alpha_{1}^{2} = \frac{1}{n}. \square

Moreover, \alpha_1^2 \ge \frac{1}{n} with some constant probability. So a few runs of the algorithm will suffice with high probability.

 

Advertisements