Iterative Row Sampling

Mu Li, Gary L. Miller, Richard Peng

Introduction

Overview

Finding such a B\mathbf{B} is equivalent to reducing the size of a regression problem involving A\mathbf{A} since:

This means finding a shorter (1±ϵ)(1\pm\epsilon) approximation to the n×(d+1)n\times(d+1) matrix [A,b][\mathbf{A},\mathbf{b}], and solving a regression problem on this approximation gives a solution within 1+O(ϵ)1+O(\epsilon) of the minimum.

The main framework of our algorithm is iterative in nature and relies on the two-way connection between row sampling and estimation of sampling probabilities. A crude approximation to A\mathbf{A}, A′\mathbf{A}^{\prime} allows us to compute equally crude approximations of sampling probabilities, while such probabilities in turn lead to higher quality approximations. The computation of these sampling probabilities can in turn be sped up using a high quality approximation of A′\mathbf{A}^{\prime}. Our algorithm is based on observing that as long as A′\mathbf{A}^{\prime} has smaller size, we have made enough progress for an iterative algorithm. A single step in this algorithm consists of computing a small but crude approximation A′\mathbf{A}^{\prime}, finding a higher quality approximation to A′\mathbf{A}^{\prime}, and using this approximation to find estimates of sampling probabilities of the rows of A\mathbf{A}. This leads to a tail-recursive process that can also be viewed as an iterative one where the calls generate a sequence of gradually shrinking matrices, and sampling probabilities are propagated back up the sequence. An example of such a sequence is given in Figure 1.

We will term the creation of the coarse approximation as reduction, and the computation of the more accurate approximation based on it recovery. As in the figure, we will label the matrices A\mathbf{A} that we generate, as well as their approximations using the indices (l)(l).

reduction: creates a smaller version of A(l)\mathbf{A}(l), A(l+1)\mathbf{A}(l+1) with fewer rows either by a projection or a coarser row sampling process. Equilvaent to moving rightwards in the diagram.

recovery: finds a small, high quality approximation of A(l)\mathbf{A}(l), B(l)\mathbf{B}(l) using information obtained from A(l)\mathbf{A}(l), A(l+1)\mathbf{A}(l+1), and B(l+1)\mathbf{B}(l+1). This is done by estimating leverage scores p(l)\mathbf{p}(l) and is equivalent to moving leftwards in the diagram.

Preliminaries

This sampling process can be formalized in several ways, leading to similar results both theoretically and experimentally [IW12]. We will treat it as a blackbox \textscSample(A,p)\textsc{Sample}(\mathbf{A},\mathbf{p}) that takes a set of probabilities over the rows of A\mathbf{A} and samples them accordingly. It keeps row ii or A\mathbf{A} with probability min⁡{1,pi}\min\{1,p_{i}\}, and rescales it appropriately so the expected value of this row is preserved. The two key properties of \textscSample(A,p)\textsc{Sample}(\mathbf{A},\mathbf{p}) that we will use repeatedly are:

It returns B\mathbf{B} with at most O(∣p∣1)O(|\mathbf{p}|_{1}) rows.

Its running time can be bounded by O(n+∣p∣1log⁡n)O(n+|\mathbf{p}|_{1}\log{n}).

The convergence of sampling relies on matrix Chernoff bounds, which can be viewed as generalizations of single variate concentration bounds. Necessary conditions on the probabilities can be formalized in several ways, with the most common being statistical leverage scores. Although these values have been studied in statistics, their use in algorithms is more recent. To our knowledge, their first use in a limited row sampling setting was in spectral sparsification of graphs [SS08]. The most general definition of pp-norm leverage scores is based on the row norms of a basis of the column space of A\mathbf{A}. However, significant simpliciations are possible when p=2p=2, and this alternate view is crucial for in our algorithm. As a result, we will state the relevant convergence results for Sample separately in Sections 4 and 5.

They show that statistical leverage scores are closely associated the probabilities needed for row sampling, and give algorithms that efficiently approximate these values. We will also formalize an observation implicit in previous results that both the sampling and estimation algorithms are very robust. The high error-tolerance of these algorithms makes them ideal as core routines to build iterative algorithms upon.

where ai\mathbf{a}_{i} is the ii-th row of A\mathbf{A}.

To our knowledge, the first near tight bounds for row sampling using statistical leverage scores were given in [AW02], and various extensions and simplifications were made since [RV07, Ver09, Har11, AT11]. They can be stated as follows:

The importance of statistical leverage scores can be reflected in the following fact, which implies that we can obtain B\mathbf{B} with O(dlog⁡d)O(d\log{d}) rows.

(see e.g. [SS08]) Given n×dn\times d matrix A\mathbf{A}, and let τ\boldsymbol{\tau} be the leverage score w.r.t. A\mathbf{A}. Assume A\mathbf{A} has rank rr, then

Although it is tempting to directly obtain high quality approximations of the leverage scores, their computation also requires a high quality approximation of ATA\mathbf{A}^{T}\mathbf{A}, leading us back to the original problem of row sampling. Our way around this issue relies on the robustness of concentration bounds such as Lemma 4.1. Sampling using even crude estimates on leverage scores can lead to high quality approximations [DDH+09, DM10, DMMS11, DMIMW12, AT11]. Therefore, we will not approximate ATA\mathbf{A}^{T}\mathbf{A} directly, and instead obtain a sequence of gradually better approximations. The need to compute sampling probabilities using crude approximations leads us to define a generalization of statistical leverage scores.

If B1\mathbf{B}_{1} and B2\mathbf{B}_{2} satisfies:

Then for any vector x\mathbf{x} we have:

This representation leads to faster algorithms for estimating stretch using the Johnson-Lindenstrauss transform. This tool is used in a variety of settings from estimating effective resistances [SS08] to more generally leverage scores [DMIMW12]. We will use the following randomized projection theorem:

2 Reductions and Recovery

Our reduction and recovery processes are based on projecting A\mathbf{A} to one with fewer rows, and moving the estimates on the projection back to the original matrix. Our key operation is to combine every RR rows into kk rows, where RR and kk are set to dθd^{\theta} and O(c/θ)O(c/\theta) respectively. By padding A\mathbf{A} with additional rows of zeros, we may assume that the number of rows is divisible by RR. We will use nb=n/Rn_{b}=n/R to denote the number of blocks, and use the notation ⋅(b)\cdot_{(b)} to index into the bbth block. Our key step is then a (R,k)(R,k)-reduction of the rows:

A (R,k)(R,k)-reduction of A\mathbf{A} describes the following procedure:

For each block A(b)\mathbf{A}_{(b)}, pick U(b)\mathbf{U}_{(b)} to be a k×Rk\times R random Gaussian matrix with entries picked independently from N(0,1)\mathcal{N}(0,1) and compute \mathbf{A}\textrm{\tiny\downarrow}_{(b)}=\mathbf{U}_{(b)}\mathbf{A}_{(b)}.

Concatenate the blocks \mathbf{A}\textrm{\tiny\downarrow}_{(b)} together vertically to form \mathbf{A}\textrm{\tiny\downarrow}.

We first show that projections preserve the stretch of blocks w.r.t. A\mathbf{A}. This can be done by bounding the effect of U(b)\mathbf{U}_{(b)} on the norm of each column of A(b)(ATA)†12\mathbf{A}_{(b)}(\mathbf{A}^{T}\mathbf{A})^{{\dagger}\frac{1}{2}}. It follows directly from properties of the Johnson-Lindenstrauss projections described in Lemma 4.5, and we’ll give its proof in Appendix B.

Assume R=dθ≥e2R=d^{\theta}\geq e^{2} for some constant θ\theta and let \mathbf{A}\textrm{\tiny\downarrow} be a (R,k)(R,k)-projection of A\mathbf{A}. For any constant c>0c>0 there exists a constant k=O(c/θ)k=O(c/\theta) such that

holds for all block b=1,…,nbb=1,\ldots,n_{b} with probability at least 1−d−c1-d^{-c}.

However, generalized stretches w.r.t. A\mathbf{A} and \mathbf{A}\textrm{\tiny\downarrow} are evaluated under the norms given by the inverses of these matrices, (ATA)+(\mathbf{A}^{T}\mathbf{A})^{+} and (\mathbf{A}\textrm{\tiny\downarrow}^{T}\mathbf{A}\textrm{\tiny\downarrow})^{+}. As a result, we need to bound the operator bound between these two pseudoinverses, which we obtain using the following lemma.

Let C\mathbf{C} and D\mathbf{D} be symmetric positive semi-definite matrices and let P\mathbf{P} be the orthogonal projection operator onto the range space of C\mathbf{C}. Then:

This is straightforward when both C\mathbf{C} and D\mathbf{D} are full rank, or share the same null space. However, as pseudo-inverses do not act on the null space, it is crucial that we’re only considering vectors of the form ai′\mathbf{a}^{\prime}_{i}. This Lemma is proven in Appendix B. Combining it with bounds in the other direction allows us to bound the distortion caused by switching reference from A\mathbf{A} to \mathbf{A}\textrm{\tiny\downarrow}.

For any constant cc, there exists a constant c′c^{\prime}, such that with probability at most 1−d−c1-d^{-c}, we have for each row ii of \mathbf{A}\textrm{\tiny\downarrow}, denoted by \mathbf{a}\textrm{\tiny\downarrow}_{i}, satisfies

Proof Denote by A(b)\mathbf{A}_{(b)} the bb-th block of A\mathbf{A} and \mathbf{A}\textrm{\tiny\downarrow}_{(b)} the corresponding block in \mathbf{A}\textrm{\tiny\downarrow}, by Lemma 4.9,

Since each U(b)\mathbf{U}_{(b)} consists of k×Rk\times R independent random variables chosen from N(0,1)\mathcal{N}(0,1), ∥U(b)∥F2\|\mathbf{U}_{(b)}\|_{F}^{2} is distributed as N(0,kR)\mathcal{N}(0,kR). This gives:

Further note that \mathbf{a}\textrm{\tiny\downarrow}_{i} is completely contained within the range space of \mathbf{A}\textrm{\tiny\downarrow}^{T}\mathbf{A}\textrm{\tiny\downarrow}. Therefore for all ii, \mathbf{P}\mathbf{a}\textrm{\tiny\downarrow}_{i}=\mathbf{a}\textrm{\tiny\downarrow}_{i} and:

For any constant cc, there exists a setting of constants such that for any R=dθR=d^{\theta}, we have with probability at least 1−d−c1-d^{-c}

3 Iterative Algorithm

Assume A(l)\mathbf{A}(l) and B(l)\mathbf{B}(l) satisfy the following condition

Then for any constant cc, there is a setting of the constants such that

holds with probability at least 1−d−c1-d^{-c}.

Proof The given condition implies that 23B(l)\sqrt{\frac{2}{3}}\mathbf{B}(l) satisfies the condition needed for Lemma 4.6 with κ=3\kappa=3. Let the constants c′=c+log⁡d2c^{\prime}=c+\log_{d}2, then with probability at least 1−d−c′1-d^{-c^{\prime}} we have:

Since A(l)\mathbf{A}(l) is a projection of A(l−1)\mathbf{A}(l-1), we can index corresponding blocks in them. Apply Corollary 4.12, then with probability at least 1−d−c′1-d^{-c^{\prime}}, we have

holds for all blocks bb with probability at least 1−2d−c′1-2d^{-c^{\prime}}, which is equal to 1−d−c1-d^{-c} by the definition of c′c^{\prime}.

and B\mathbf{B} has O(dlog⁡dϵ−2){O}(d\log{d}\epsilon^{-2}) rows, each being a scaled copy of some row of A\mathbf{A},

Proof We first show correctness via. induction backwards on ll. Define c′=c+1c^{\prime}=c+1, we show that B(l)\mathbf{B}(l) has O(dR4log⁡dϵ(l)−2)O(dR^{4}\log{d}\epsilon(l)^{-2}) rows and satisfies

with probability at least 1−3(L−l)dc′1-3(L-l)d^{c^{\prime}} for each ll.

As k=O(c/θ)k=O({c}/{\theta}) is a constant, one (R,k)(R,k)-projection decreases the number of rows by a factor of O(R)O(R). After L=log⁡R(n/d)L=\log_{R}(n/d) projections, we get that A(L)\mathbf{A}(L) has O(d)O(d) rows. Therefore the base case where l=Ll=L follows from B(L)=A(L)\mathbf{B}(L)=\mathbf{A}(L).

For the inductive step, we assume that the inductive hypothesis holds for l≥1l\geq 1 and try to show it for l−1l-1. As ϵ(l)\epsilon(l) was set to 1/21/2, we have:

This allows us to invoke Lemma 4.13, which combined with Lemma 4.1 gives that with probability 1−d−c′1-d^{-c^{\prime}}, the inductive hypothesis also holds for l−1l-1.

We will use this Fact with one of pp or qq being 22, in which case it gives:

If 1≤p≤21\leq p\leq 2, ∥x∥2≤∥x∥p≤d12−1p∥x∥2\left\|\mathbf{x}\right\|_{2}\leq\left\|\mathbf{x}\right\|_{p}\leq d^{\frac{1}{2}-\frac{1}{p}}\left\|\mathbf{x}\right\|_{2}

If 2≤p2\leq p, d1p−12∥x∥2≤∥x∥p≤∥x∥2d^{\frac{1}{p}-\frac{1}{2}}\left\|\mathbf{x}\right\|_{2}\leq\left\|\mathbf{x}\right\|_{p}\leq\left\|\mathbf{x}\right\|_{2}

Let A\mathbf{A} be an n×dn\times d matrix of rank rr, p∈[1,∞]p\in[1,\infty] and qq be its dual norm such that 1p+1q=1\frac{1}{p}+\frac{1}{q}=1. Then an n×rn\times r matrix U\mathbf{U} is an (α,β,p)(\alpha,\beta,p)-well-conditioned basis for the column space of A\mathbf{A} if the columns of U\mathbf{U} span the column space of A\mathbf{A} and:

∣∣∣U∣∣∣p≤α{\left|\kern-0.6458pt\left|\kern-0.6458pt\left|\mathbf{U}\right|\kern-0.6458pt\right|\kern-0.6458pt\right|}_{p}\leq\alpha.

where Ui∗\mathbf{U}_{i*} is the ii-th row of U\mathbf{U} and cpc_{p} is a constant depending only on pp. Then with probability at least 1−d−c1-d^{-c}, \textscSample(A,p,ϵ)\textsc{Sample}(\mathbf{A},p,\epsilon) returns B\mathbf{B} satisfying

The estimation of ∥(AQ)i∗∥pp\left\|(\mathbf{A}\mathbf{Q})_{i*}\right\|_{p}^{p} and sampling can then be done in a way similar to Section 4. When 1≤p≤21\leq p\leq 2, we can compute O(1)O(1) approximations using pp-stable distributions [Ind06] in a way analogous to Section 4.2.1. of Clarkson et al. [CDMI+12]. When 2≤p2\leq p, we will use the 22-norm as a surrogate at the cost of more rows. As all of our calls to Sample will be using probabilities estimated via. the same matrix, we will the estimation of pp-norm leverage scores and sampling as a single blackbox.

For any constant cc, there exist an algorithm \textscEstimateAndSampleP(A,C,α,β,R,ϵ)\textsc{EstimateAndSampleP}(\mathbf{A},\mathbf{C},\alpha,\beta,R,\epsilon) that given a A\mathbf{A}, C\mathbf{C} such that AC\mathbf{A}\mathbf{C} is a (α,β,p)(\alpha,\beta,p)-well-conditioned basis for A\mathbf{A}, returns a matrix B\mathbf{B} with probability at least 1−d−c1-d^{-c} such that:

Proof We start by show that UTU\mathbf{U}^{T}\mathbf{U} is close to the identity matrix as an operator. Also, since CTC\mathbf{C}^{T}\mathbf{C} is a full rank matrix and CT(CCT)†C\mathbf{C}^{T}(\mathbf{C}\mathbf{C}^{T})^{{\dagger}}\mathbf{C} is a projection operator onto the column space of C\mathbf{C}, we have CT(CCT)†C=I.\mathbf{C}^{T}(\mathbf{C}\mathbf{C}^{T})^{{\dagger}}\mathbf{C}=\mathbf{I}. Taking pseudoinverses of the given condition on C\mathbf{C} gives:

Substituting it into U=AC\mathbf{U}=\mathbf{A}\mathbf{C} then gives:

This allows us to infer that ∣∣∣U∣∣∣2≤2d{\left|\kern-0.6458pt\left|\kern-0.6458pt\left|\mathbf{U}\right|\kern-0.6458pt\right|\kern-0.6458pt\right|}_{2}\leq\sqrt{2d}, and for any vector x\mathbf{x}, ∥x∥2≤2∥Ux∥2\left\|\mathbf{x}\right\|_{2}\leq\sqrt{2}\left\|\mathbf{U}\mathbf{x}\right\|_{2}.

Next we find values of α\alpha and β\beta that meet the requirements of a well-conditioned basis given in Definition 5.2. Let qq be the dual norm for pp which satisfies 1p+1q=1\frac{1}{p}+\frac{1}{q}=1.

First consider the case where 1≤p≤21\leq p\leq 2. We can view all entries of the matrix U\mathbf{U} as a vector of length nr≤ndnr\leq nd vector. Apply Fact 5.1 gives:

Which gives ∣∣∣U∣∣∣2≤2d12{\left|\kern-0.6458pt\left|\kern-0.6458pt\left|\mathbf{U}\right|\kern-0.6458pt\right|\kern-0.6458pt\right|}_{2}\leq\sqrt{2}d^{\frac{1}{2}} and therefore α=2(nd)1p−12d12\alpha=\sqrt{2}(nd)^{\frac{1}{p}-\frac{1}{2}}d^{\frac{1}{2}}. For the second part, given any vector x\mathbf{x}, we have

Now we consider the case where p≥2p\geq 2 similarly.

Combining the bounds from these two cases on pp gives that UU is a (α,β,p)−(\alpha,\beta,p)-well-conditioned basis, where

It can be checked that in both cases the stated bound on αβ\alpha\beta holds. ■\blacksquare

And the number of rows in B\mathbf{B} can be bounded by:

Proof It suffices to verify both conditions of Definition 5.2 holds for U=AC\mathbf{U}=\mathbf{A}\mathbf{C}. For ∣∣∣U∣∣∣p{\left|\kern-0.6458pt\left|\kern-0.6458pt\left|\mathbf{U}\right|\kern-0.6458pt\right|\kern-0.6458pt\right|}_{p}, we can treat its ppth power as a summation over the columns of U\mathbf{U} and get:

Proof of Lemma 5.6: By the guarantees of RowCombineL2 given in Theorem 4.14, we can set its constants so that with probability at least 1−d−c′1-d^{-c^{\prime}} we have:

And the number of rows in B\mathbf{B} can be bounded by

To bound the runtime, we first show that if B\mathbf{B} has nbn_{b} rows, then in the next iteration the number of rows in B\mathbf{B} can be bounded by (nb/n∗)cpn∗(n_{b}/n*)^{c_{p}}n* where cp<1c_{p}<1 is a constant based on pp. When 1≤p≤21\leq p\leq 2, Lemma 5.6 gives that there exist some constant c0c_{0} such that the new row count can be bounded by:

So n∗=c02pd4p+O(θ)log⁡2pdn^{*}=c_{0}^{\frac{2}{p}}d^{\frac{4}{p}+O(\theta)}\log^{\frac{2}{p}}{d} and cp=1−p2c_{p}=1-\frac{p}{2} suffices. Similarly for the case where 2≤p2\leq p, it can be checked that

Gives that the new row count can be bounded by (nb/n∗)cpn∗(n_{b}/n^{*})^{c_{p}}n^{*} where cp=p2−1c_{p}=\frac{p}{2}-1. As we can set θ\theta to any arbitrary constant, the constant in front of its exponent can also be removed.

2 Fewer Rows by Iterating Again

A closer look at the proof of Lemma 5.8 shows that a significant increase in the number of rows comes from dividing by the 1−∣1−p2∣1-|1-\frac{p}{2}| term in the exponent of nn. As a result, the row count can be further reduced if the leverage scores are computed via. a p′p^{\prime}-norm approximation where qq is between 22 and pp. For simplicity we only show this improvement for the case where 1≤p≤p′≤21\leq p\leq p^{\prime}\leq 2.

We will start by proving a generalization of Lemma 5.5.

This leads to a result analogous to Lemma 5.6. Solving for the fixed point of this process allows us to prove Theorem 2.2.

Proof of Theorem 2.2: Similar to the proof of Lemma 5.8, in each iteration we reduce the number of rows from nn to O(n1−pp′(d4p′log⁡2p′d)pp′−p2Rpd2log⁡2pp′+1d)O\left(n^{1-\frac{p}{p^{\prime}}}\left(d^{\frac{4}{p^{\prime}}}\log^{\frac{2}{p^{\prime}}}{d}\right)^{\frac{p}{p^{\prime}}-\frac{p}{2}}R^{p}d^{2}\log^{\frac{2p}{p^{\prime}}+1}{d}\right). Ignoring terms in RR and log⁡d\log{d}, we have that the number of rows converges doubly exponentially towards:

This is minimized when p′=2p^{\prime}=\sqrt{2}, giving 4p′+2p′−2=42−2\frac{4}{p^{\prime}}+2p^{\prime}-2=4\sqrt{2}-2. As we can set R=dO(θ)R=d^{O(\theta)}, the number of rows in B\mathbf{B} can be bounded by O(d42−2+θ)O(d^{4\sqrt{2}-2+\theta}) for any constant θ\theta. ■\blacksquare

References

Appendix A Properties and Estimation of Stretch and Leverage Scores

We now give proofs for estimating leverage scores and row sampling that we stated in Sections 4 and 5.

Proof of Fact 4.4: Note that both stretch and the Frobenius norm acts on the rows independently. Therefore it suffices to prove this when A′\mathbf{A}^{\prime} has a single row, aka. A′=a\mathbf{A}^{\prime}=\mathbf{a}. In this case the cyclic property of trace gives:

Proof of Lemma 4.3: The condition given implies that the null spaces of B1TB1\mathbf{B}_{1}^{T}\mathbf{B}_{1} and B2TB2\mathbf{B}_{2}^{T}\mathbf{B}_{2} are identical, giving:

Applying this to the vector x\mathbf{x} gives:

A.2 Estimation of Generalized Stretch

Based on this fact, we can estimate these scores using randomized projections in a way that’s by now standard [SS08, DMIMW12]. Pseudocode of our estimation algorithm is shown in Algorithm 4, while the error analysis is nearly identical to the ones given in Section 4 of [SS08] and Section 3.2. of [DMIMW12].

We remark that (BTB)†12(\mathbf{B}^{T}\mathbf{B})^{{\dagger}\frac{1}{2}} can be replaced by any matrix whose product with its transpose equals to BTB\mathbf{B}^{T}\mathbf{B}. nee candidate for this is B(BTB)†\mathbf{B}(\mathbf{B}^{T}\mathbf{B})^{{\dagger}}, and using it would avoid computing the 1/21/2 power of a matrix. However, from a theoretical point of view both of these operations take O(mdω−1)O(md^{\omega-1}) time, and we omit this extra step for simplicity.

where the last inequality is due to the assumption that R≤e2R\leq e^{2}. By a suitable choice of constants in k=O(log⁡Rndc)=O(log⁡Rd)k=O(\log_{R}{nd^{c}})=O(\log_{R}d) this can be made the above probability less than n−1d−cn^{-1}d^{-c}, taking a union bound over the nn rows gives Part 1.

A.3 Estimation of pp-Norm Leverage Scores

We will estimate the values of ∥Ui∗∥p\left\|\mathbf{U}_{i*}\right\|_{p} using similar dimensionality reduction theorems. Specifically, we utilize a result on pp-stable distributions first shown by Indyk [Ind06].

Note that the result from [Ind06] was only stated in terms of obtaining 1±ϵ1\pm\epsilon approximations, which leads to a factor of log⁡d\log{d} on the leading term. However, these bounds can be obtained analogously by the fact that pp-stable distributions have bounded derivative.

The guarantees on the output then follows from computing Lemma 5.3, while the total running time follows from the cost of evaluating A(CΠT)\mathbf{A}(\mathbf{C}\Pi^{T}). ■\blacksquare

Appendix B Deferred Proofs from Section 4

Where the last inequality is by 1−R−1−ln⁡R≤−12ln⁡R1-R^{-1}-\ln{R}\leq-\frac{1}{2}\ln{R} with the assumption that R≥e2R\geq e^{2}. If we kk to 4(c+1)θ−1log⁡dn4(c+1)\theta^{-1}\log_{d}n and substitute R=dθR=d^{\theta}, we get:

Proof of Lemma 4.10: Consider an orthonormal basis for the range space of C\mathbf{C}, v1…vrank(C)\mathbf{v}_{1}\ldots\mathbf{v}_{\textrm{rank}(\mathbf{C})}. Since C+D⪰C\mathbf{C}+\mathbf{D}\succeq\mathbf{C}, this basis can be extended to an orthonormal basis to the range space of C+D\mathbf{C}+\mathbf{D} by adding vrank(C)+1…vrank(C+D)\mathbf{v}_{\textrm{rank}(\mathbf{C})+1}\ldots\mathbf{v}_{\textrm{rank}(\mathbf{C}+\mathbf{D})}. It suffices to prove the claim under this basis system. Here C\mathbf{C} and D\mathbf{D} can be rewritten as by proper rotation:

where C11\mathbf{C}_{11} and D22\mathbf{D}_{22} are strictly positive definite. Furthermore, since D\mathbf{D} is positive semi-definite we have that D11−D12D22−1D12T\mathbf{D}_{11}-\mathbf{D}_{12}\mathbf{D}_{22}^{-1}\mathbf{D}_{12}^{T} is also positive semi-definite. For any vector x\mathbf{x}, PCx\mathbf{P}_{\mathbf{C}}\mathbf{x} gives a vector that’s non-zero only in the first rank(C)\textrm{rank}(\mathbf{C}) entries. Let this part be x1\mathbf{x}_{1}. Then evaluating (C+D)†PCx=[y1;y2](\mathbf{C}+\mathbf{D})^{{\dagger}}\mathbf{P}_{\mathbf{C}}\mathbf{x}=[\mathbf{y}_{1};\mathbf{y}_{2}] becomes solving the following system:

The second equation gives y2=−D22−1D12Ty1\mathbf{y}_{2}=-\mathbf{D}_{22}^{-1}\mathbf{D}_{12}^{T}\mathbf{y}_{1}. Substituting it into the first one gives:

Note that this is the same as taking the partial Cholesky factorization onto the range space of PC\mathbf{P}_{\mathbf{C}}. Combining things gives:

Since both D11−D12D22−1D12T\mathbf{D}_{11}-\mathbf{D}_{12}\mathbf{D}_{22}^{-1}\mathbf{D}_{12}^{T} and D11\mathbf{D}_{11} are positive definite, we have C11⪯C+D11−D12D22−1D12T\mathbf{C}_{11}\preceq\mathbf{C}+\mathbf{D}_{11}-\mathbf{D}_{12}\mathbf{D}_{22}^{-1}\mathbf{D}_{12}^{T} and therefore:

holds for every x\mathbf{x}. ■\blacksquare