Say we have a matrix A∈Rn×dA \in \mathbb{R}^{n \times d} where n≫dn \gg d. If a linear map S∈Rk×nS \in \mathbb{R}^{k \times n}, where k≪nk \ll n, has the property that, for all x∈Rdx \in \mathbb{R}^d and some ϵ∈(0,1)\epsilon \in (0, 1), the following holds:

(1−ϵ)∥Ax∥2≤∥SAx∥2≤(1+ϵ)∥Ax∥2, (1 - \epsilon) \lVert Ax \rVert_2 \leq \lVert SAx \rVert_2 \leq (1+\epsilon) \lVert Ax \rVert_2,

then we call SS an ℓ2\ell_2-subspace embedding of the column space of AA, or simply an ℓ2\ell_2-subspace embedding of AA. Essentially, SS preserves norms of the column space of AA up to ϵ\epsilon amount of error. To simplify notation, I’ll write the above expression as: ∥SAx∥2=(1±ϵ)∥Ax∥2\lVert SAx \rVert_2 = (1 \pm \epsilon) \lVert Ax \rVert_2.

We can, without loss of generality, assume that AA has orthonormal columns. That is because, if that’s not the case, we can redefine A:=UA:=U where UU comes from the SVD of A=UΣV∗A=U\Sigma V^\ast. Indeed, the sets {y  ∣  y=Ax,  x∈Rd}\{ y \;|\; y = Ax,\; x \in \mathbb{R}^d \} and {y  ∣  y=Uz,  z∈Rρ}\{ y \;|\; y=Uz, \; z \in \mathbb{R}^\rho \} where ρ=rank(A)\rho=\mathop{rank}(A) are the same.

What that means is that, we can simplify the condition to: ∥SAx∥2=(1±ϵ)∥x∥2\lVert SAx \rVert_2 = (1 \pm \epsilon) \lVert x \rVert_2. We can also require that the above hold for unit vectors xx due to the linearity of the maps involved. That is in turn equivalent to ∥I−A∗S∗SA∥2≤ϵ\lVert I - A^\ast S^\ast S A \lVert_2 \leq \epsilon where ∥⋅∥2\lVert \cdot \rVert_2 is the operator norm of its argument. To see this:

∥SAx∥2=(1±ϵ)∥x∥2⇒∣x∗(A∗S∗SA−I)x∣≤ϵ∥x∥2⇒sup⁡∥x∥2=1∣x∗(A∗S∗SA−I)x∣≤ϵ⇒sup⁡∥x∥2=1∣x∗Σx∣≤ϵ⇒σ1(A∗S∗SA−I)≤ϵ \begin{aligned} \lVert S&A x \rVert_2 = (1 \pm \epsilon) \lVert x \rVert_2 \Rightarrow \\ & \big\lvert x^\ast \big( A^\ast S^\ast S A - I \big) x \big\rvert \leq \epsilon \lVert x \rVert_2 \Rightarrow \\ & \sup_{\lVert x \rVert_2 = 1} \big\lvert x^\ast \big( A^\ast S^\ast S A - I \big) x \big\rvert \leq \epsilon \Rightarrow \\ & \sup_{\lVert x \rVert_2 = 1} \big\lvert x^\ast \Sigma x \big\rvert \leq \epsilon \Rightarrow \\ & \sigma_1(A^\ast S^\ast S A - I) \leq \epsilon \end{aligned}

where we used the SVD (A∗S∗SA−I)=UΣV∗(A^\ast S^\ast SA - I)=U \Sigma V^\ast and the fact that UU and VV have orthonormal columns.

Preserve a whole subspace, not just a point cloud The unit sphere in a fixed d-dimensional column space is mapped by the same sketch S into an annulus between radii one minus epsilon and one plus epsilon. The ambient row count changes from n to k, but every vector in the fixed subspace is covered by the guarantee. Preserve a whole subspace, not just a point cloud a fixed column space y = Ax, ||y|| = 1 same S n → k rows all directions stay in this band 1 − ε 1 + ε 1 − ε ≤ ||Sy|| ≤ 1 + ε One successful S preserves every vector in the same fixed subspace.
The same successful sketch preserves every direction in this fixed column space. A unit sphere becomes a shape trapped between radii 1 − ε and 1 + ε; linearity extends the guarantee to vectors of any length. The plane is a schematic view of a low-dimensional subspace in a larger ambient space. View largerScroll sideways for the full sketch.

What properties should SS have? Well, first and foremost, we want k≪nk \ll n so that SASA has as few rows as possible. That is, we want the smallest subspace of the column space of AA that offers a (1±ϵ)(1 \pm \epsilon)-approximation of the Euclidean norm. Additionally, we want to be able to compute SASA rather quickly. That often means we prefer a sparser SS.

Definition

Suppose that we have a distribution Π\Pi over the space of all k×nk \times n linear maps, where kk is a function of n,d,ϵ,δn, d, \epsilon, \delta, with the following property: For every S∼ΠS \sim \Pi and any fixed matrix A∈Rn×dA \in \mathbb{R}^{n \times d}, with probability at least 1−δ1 - \delta, SS is a (1±ϵ)(1 \pm \epsilon) ℓ2\ell_2-subspace embedding of AA. Then we call Π\Pi an (ϵ,δ)(\epsilon, \delta) oblivious ℓ2\ell_2-subspace embedding.

As a side-note, in addition to oblivious embeddings, we can design sampling-based sketching algorithms that are called Leverage Score Sampling. I’ll discuss these in a later post.

Construction

Theorem: Let S=1kR∈Rk×nS=\frac{1}{\sqrt{k}}R \in \mathbb{R}^{k \times n} where the entries of RR are iid standard normal random variables. If k=Θ(ϵ−2(d+log⁡(1/δ)))k = \Theta\big(\epsilon^{-2} (d + \log(1/\delta)) \big), then SS is a (1±ϵ)(1 \pm \epsilon) ℓ2\ell_2-subspace embedding for any matrix A∈Rn×dA \in \mathbb{R}^{n \times d} with probability at least 1−δ1 - \delta.

Proof: Notice that, we are trying to find an ℓ2\ell_2-subspace embedding for the column space of AA. We can also focus only on unit vectors in the column space, due to linearity. The idea behind the proof is to first lay an ϵ\epsilon-net over the column space of AA with the property that, a Johnson-Lindenstrauss Transform for the ϵ\epsilon-net is a ℓ2\ell_2-subspace embedding for the entire space. Let’s get to it.

Consider the set S={y  ∣  y=Ax,  ∥y∥2=1}\mathcal{S} = \{ y \;|\; y = Ax,\; \lVert y \rVert_2 = 1 \} and its 1/21/2-net N\mathcal{N}. We can write every y∈Sy \in S as the series: y=∑i=0∞wiy = \sum_{i=0}^{\infty} w_i where wiw_i is a scaled member of N\mathcal{N} with the property that ∥wi∥≤2−i\lVert w_i \rVert \leq 2^{-i}. To see that that’s always possible, take w0:=v0∈Nw_0 := v_0 \in \mathcal{N} and write y=w0+(y−w0)y = w_0 + (y - w_0). Clearly ∥y−w0∥2≤1/2\lVert y - w_0 \rVert_2 \leq 1/2 because w0∈Nw_0 \in \mathcal{N}. Now write y=w0+w1+(y−w0−w1)y = w_0 + w_1 + (y - w_0 - w_1) where w1=∥y−w0∥2v1w_1 = \lVert y - w_0 \rVert_2 v_1 for some v1∈Nv_1 \in \mathcal{N}. Then:

∥y−w0−w1∥2=∥y−w0∥2⋅∥y−w0∥y−w0∥2−v1∥2≤∥y−w0∥22≤14, \lVert y - w_0 - w_1 \rVert_2 = \lVert y - w_0 \rVert_2 \cdot \lVert \frac{y - w_0}{\lVert y - w_0 \rVert_2} - v_1 \rVert_2 \leq \frac{\lVert y - w_0 \rVert_2}{2} \leq \frac{1}{4},

as desired. We can continue to expand the residual as above and reason by induction.

A finite net controls an infinite sphere Approximate a unit vector by a net point and repeatedly approximate the normalized residual by another net point. With a half-net, residual norms shrink by a factor of at least two and the sum of the norms of all scaled net vectors is at most two. Bilinear errors accumulate to at most four epsilon. A finite net controls an infinite sphere 1. Pick a net direction 2. Approximate the residual 3. Repeat at half the scale residual size after each step start 1 w0 ≤ 1/2 w1 ≤ 1/4 w2 ≤ 1/8 y = w0 + w1 + w2 + … each wi is a scaled direction from the same net Σ ||wi|| ≤ 1 + 1/2 + 1/4 + … = 2 Pairwise net errors of ε accumulate to at most ε(Σ ||wi||)2 ≤ 4ε.
Choose a direction from a 1/2-net, subtract it, and repeat with the normalized residual. The residual halves at each step. Consequently Σ||wᵢ|| ≤ 2, so an additive ε error for every pair of net directions contributes at most 4ε to the squared norm of any unit vector. View largerScroll sideways for the full sketch.

Now consider a JLT for N\mathcal{N}, SS, so that with probability at least 1−δ1 - \delta, we have that ⟨Sv,Sv′⟩=⟨v,v′⟩±ϵ\langle Sv, Sv^\prime \rangle = \langle v, v^\prime \rangle \pm \epsilon for all pairs of v,v′∈Nv, v^\prime \in \mathcal{N}. If we applied the same map to any y∈Sy \in \mathcal{S}, then we’d get, with probability at least 1−δ1 - \delta:

⟨Sy,Sy⟩=⟨∑Swi,∑Swi⟩=∑i⟨Swi,Swi⟩+∑i≠j⟨Swi,Swj⟩=∑i∥wi∥22+∑i≠j⟨wi,wj⟩±ϵ∑i,j∥wi∥2∥wj∥2=⟨∑wi,∑wi⟩+O(ϵ)=⟨y,y⟩±O(ϵ). \begin{aligned} \langle Sy, Sy \rangle &= \langle \sum S w_i, \sum Sw_i \rangle \\ &= \sum_i \langle Sw_i, Sw_i \rangle + \sum_{i \neq j} \langle Sw_i, Sw_j \rangle \\ &= \sum_i \lVert w_i \rVert_2^2 + \sum_{i \neq j} \langle w_i, w_j \rangle \pm \epsilon \sum_{i, j} \lVert w_i \rVert_2 \lVert w_j \rVert_2 \\ &= \langle \sum w_i, \sum w_i \rangle + \mathcal{O}(\epsilon) \\ &= \langle y, y \rangle \pm \mathcal{O}(\epsilon). \end{aligned}

So, we have established that a JLT for the 1/21/2-net over the unit sphere in the column space of AA is an ℓ2\ell_2-subspace embedding for the column space of AA. We know that ∣N∣≤5d\lvert \mathcal{N} \rvert \leq 5^d, so that k=Ω(ϵ−2log⁡(5d/δ))=Ω(ϵ−2(d+log⁡(1/δ)))k = \Omega(\epsilon^{-2} \log (5^d/\delta)) = \Omega(\epsilon^{-2} (d + \log (1/\delta))) as desired.