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}$ has (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)

so 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$ . This helps because running the RSVD on $(\vec{A}\vec{A}^\T)^q\vec{A}$ scales the sketch’s component along the $i$ -th singular direction by $\sigma_i^{2q+1}$ : shrinking the tail relative to the top makes the sketch, and hence $\vec{Q}$ , line up more closely with the leading singular directions. As illustrated in the figure below, the small singular values become smaller relative to the large ones. Drag the slider to raise the singular values to the power $2q+1$ and observe that the tail is damped substantially relative to the leading singular values as $q$ increases.

import numpy as np
import pandas as pd
import altair as alt

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

# Pre-compute the damped spectrum across the whole slider range so that the
# browser (Vega) can filter by q client-side without re-running Python.
q_max = 8
poly_data = pd.DataFrame([
    {'i': int(i), 'sigma': float(Λ[i] ** (2 * q + 1)), 'q': int(q), 'series': 'damped'}
    for q in range(q_max + 1) for i in range(d)
])
ref_data = pd.DataFrame({'i': np.arange(d), 'sigma': Λ, 'series': 'original'})

q_slider = alt.binding_range(min=0, max=q_max, step=1, name='q: ')
q_param = alt.param(name='q_param', value=1, bind=q_slider)

color = alt.Color('series:N',
                  scale=alt.Scale(domain=['original', 'damped'], range=['#b0b0b0', '#1e3668']),
                  legend=alt.Legend(title=None, orient='top-right'))

# Original singular values (q = 0), always shown for reference
ref = alt.Chart(ref_data).mark_circle(size=55, opacity=0.9).encode(
    x=alt.X('i:Q', title='index i', scale=alt.Scale(domain=[-0.5, d - 0.5]),
            axis=alt.Axis(tickCount=5, grid=False)),
    y=alt.Y('sigma:Q', title='singular value', scale=alt.Scale(domain=[0, 1.03]),
            axis=alt.Axis(grid=False)),
    color=color,
)

# Damped singular values for the value of q chosen by the slider
active_pts = alt.Chart(poly_data).transform_filter(alt.datum.q == q_param).mark_point(
    filled=True, size=60).encode(x='i:Q', y='sigma:Q', color=color)

chart = (ref + active_pts).add_params(q_param).properties(
    width=480, height=300,
    title={'text': {'expr':
        '"singular values raised to the power 2q+1   (q = " + q_param + ")"'}}
).configure_view(strokeWidth=0)

chart

Loading...

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$ .

5.2 Subspace Iteration