Subspace Iteration - Randomized Numerical Linear Algebra with Examples

The Randomized SVD (RSVD) produces an approximation

\widehat{\vec{A}} = \vec{Q}\vec{Q}^\T \vec{A}, \quad \vec{Q} = \Call{orth}(\vec{A}\vec{\Omega}),

(5.4)

where $\vec{\Omega}\sim \operatorname{Gaussian}(n,b)$ . Ideally $\vec{Q}$ is well-aligned with the dominant subspace of $\vec{A}$ . However, the approximation can be poor if the singular values of $\vec{A}$ decay slowly or if the spectral gap is small.

One way to mitigate this, is to damp down the tail of $\vec{A}$ relative to the leading singular values. In particular, observe that if $\vec{A}$ as (thin) SVD $\vec{A} = \vec{U}\vec{\Sigma}\vec{V}^\T$ ,

(\vec{A}\vec{A}^\T)^q\vec{A} = \vec{U} \vec{\Sigma}^{2q+1} \vec{V}^\T.

(5.5)

The singular values of $(\vec{A}\vec{A}^\T)^q\vec{A}$ are the singular values of $\vec{A}$ raised to the power $2q+1$ . Thus, as illustrated in Figure 5.1 the small singular values become smaller relative to the large ones.

Four subplots showing singular value decay for different powers q, demonstrating how small singular values are damped relative to large ones as q increases — Figure 5.1:Observe that the singular values tail of $(\vec{A}\vec{A}^\T)^q\vec{A}$ is damped substantially relative to that of $\vec{A}$ .

This leads to the Randomized Subspace Iteration (RSI) approximation:

\widehat{\vec{A}} = \vec{Q}\vec{Q}^\T \vec{A}, \quad \vec{Q} = \Call{orth}((\vec{A}\vec{A}^\T)^q\vec{A}\vec{\Omega}).

(5.6)

Observe that $(\vec{A}\vec{A}^\T)^q\vec{A}\vec{\Omega}$ can be computed by sequential products with $\vec{A}$ and $\vec{A}^\T$ . In particular, we never need to form the (potentially large) matrix $\vec{A}\vec{A}^\T$ explicitly.

Note that for numerical stability reasons, it is often recommended to re-orthogonalize $\vec{Y}$ after each multiplication by $\vec{A}$ or $\vec{A}^\T$ .

import numpy as np
import matplotlib.pyplot as plt

d = 20
Λ = np.sort(np.linspace(0,.9,d) + 1e-1*np.random.rand(d))[::-1]
Λ[0] = 1


fig,axs = plt.subplots(1,4,figsize=(8,4),sharey=True)

qs = [1,2,4]
axs[0].plot(np.arange(d),Λ,ls='None',marker='o',label='original singular values')
axs[0].set_title('original ($q=0$)')
axs[0].set_xlabel('index $i$')

for i,q in enumerate(qs):
    ax = axs[i+1]

    ax.plot(np.arange(d),Λ**(2*q+1),ls='None',marker='o',label=r'$\sigma_i(A)^{2q+1}$')

    ax.set_title(f'$q={q}$')

    ax.set_xlabel('index $i$')

axs[0].set_ylabel('singular value')
plt.savefig('SI.svg')

Randomized Numerical Linear Algebra with Examples

Randomized SVD

Randomized Numerical Linear Algebra with Examples

Block Krylov Iteration

5.2Subspace Iteration