Converting one matrix factorization to another

(Note: The following content is mainly excerpted from the reference [1] listed at the bottom.)

---

Suppose that we have obtained a partial decomposition of a matrix $A$ by some means:
$$\Vert A-CB\Vert \le ϵ,\text{\hspace{1em}(1)}$$
where $B$ and $C$ have rank $k$.

---

1 Preliminaries (Some definitions)

1.1 Full Decomposition and Partial Decomposition

It is possible to compute the full QR decomposition for the full SVD of an $m\times n$ matrix to double-precision accuracy with O ( m n min ⁡ { m , n } ) O(mn\min\{m,n\}) O(mnmin{m,n}) flops. (That is, full decomposition is an accurate decomposition. ) The full decompostion has the following formula,
$$A=CB\text{\hspace{1em}(2)}$$

Suppose that an $m\times n$ matrix has numerical rank $k$, where $k$ is substantially smaller than $m$ and $n$. In this case, it is possible to produce a structured low-rank decomposition that approximates the matrix well. (That is, partial decomposition is an approximate decomposition.) The partial decomposition has the following formula,
$$A=CB+E\text{\hspace{1em}(3)}$$
where $\Vert E\Vert \le ϵ$.

1.2 The Interpolative Decomposition (ID)

The final factorization identifies a collection of $k$ columns from a rank-$k$ matrix $A$ that span the range of $A$. To be precise, we can compute an index set $J=[j_1,\cdots ,j_k]$ such that
$$A=A(:,j)X,\text{\hspace{1em}(4)}$$
where $X$ is a $k\times n$ matrix that satisfies $X(:,J)=I_k$. Furthermore, no entry of $X$ has magnitude larger than two. In other words, this decomposition express each column of $A$ using a linear combination of $k$ fixed columns with bounded coefficients. Stable and efficient algorithms for computing the ID appear in the papers [2,3].

It is also possible to compute a two-sided ID,
$$A=WA(J^{\prime },J)X,\text{\hspace{1em}(5)}$$
where $J^{\prime }$ is an index set identifying $k$ of the rows of $A$, and $W$ is an $m\times k$ matrix that satisfies $W(J^{\prime },:)=I_k$ and whose entries are all bounded by two.

---

2 From given factorization to QR

Given a partial decompostion of matrix $A$ in (1), the following steps construct a partial QR factorization:

• Compute a QR factorization of $C$ so that $C=Q_1{R}_1$;
• From the product $D=R_{1}B$, and compute a QR factorization: $D=Q_{2}R$;
• From the product $Q=Q_1{Q}_2$.

The result is an orthonormal matrix $Q$ and a weekly upper-triangular matrix $R$ such that $\Vert A-QR\Vert \le ϵ$.

---

3 From given factorization to SVD

Given a partial decompostion of matrix $A$ in (1), the following steps construct a partial SVD factorization:

• Compute a QR factorization of $C$ so that $C=Q_1{R}_1$;
• Form the product $D=R_{1}B$, and compute an SVD: $D=U_2\Sigma V^{\ast }$;
• Form the product $U=Q_1{U}_2$.

The result is a diagonal matrix $\Sigma$ and orthonormal matrices $U$ and $V$ such that $\Vert A-U\Sigma V^{\ast }\Vert \le ϵ.$

---

4 From given factorization to ID

Given a partial decompostion of matrix $A$ in (1), the following one-step process can get the ID factorization:

• Compute $J$ and $X$ such that $B=B(:,J)X$.

Then $A\approx A(:,J)X$, but the approximation error may deteriorate from the initial estimate. That is, the upper bound of $\Vert A-A(:,J)X\Vert$ is larger than $ϵ$. In other words, $\Vert A-A(:,J)X\Vert \le (1+\alpha )ϵ$.

---

5 Postprocessing via Row Extraction

Given a matrix $Q$ such that the follwoing inequlity holds,
$$\Vert A-Q{Q}^{\ast }A\Vert \le ϵ\text{\hspace{1em}(6)}$$
we can obtain a rank-$k$ factorization
$$A\approx XB,\text{\hspace{1em}(7)}$$
where $B$ is a $k\times n$ matrix consisting of $k$ rows extracted from $A$.

The approximation (7) can be produced without computing any matrix-matrix products, which makes this approach to postprocessing very fast.

The drawback comes because the error $\Vert A-XB\Vert$ is usually larger than the initial error |A-QQ^{*}A|, especially when the dimensions of $A$ are large.

To obtain the factorization (7), we simply construct the ID of the matrix $Q$:
$$Q=XQ(J,:).\text{\hspace{1em}(8)}$$
The index set $J$ marks $k$ rows of $Q$ that span the row space of $Q$, and $X$ is an $m\times k$ matrix whose entries are bounded in magnitude by two and contains the $k\times k$ identity as a submatrix: $X(J,:)=I_k$. Combining (8) and (6), we reach
$$A\approx Q{Q}^{\ast }A=XQ(J,:)Q^{\ast }A,\text{\hspace{1em}(9)}$$
Since $X(J,:)=I_k$, equation (9) implies that $A(J,:)\approx Q(J,:)Q^{\ast }A$. Therefore, (7) follows when we put $B=A(J,:)$.

---

6 SVD via Row Extraction

Provided with the factorization (7), we can obtain any standard factorization (QR, SVD et.al). The following algorithm illustrates an SVD calculation. This procedure requires $O(k^2(m+n))$ flops.

The following lemma guarantees the accuracy of the computed factors.

Lemma 5.1 [1]
Let $A$ be an $m\times n$ matrix and let $Q$ be an $m\times k$ matrix that satisfy (6). Suppose that $U$, $\Sigma$, and $V$ are the matrices constructed by the above algorithm. Then
$$\Vert \mathbf{A}-\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\ast }\Vert \le [1+\sqrt{1+4k(n-k)}]\epsilon .\text{\hspace{1em}(10)}$$

Proof: The factors $U,\Sigma ,V$ constructed by the algorithm satisfy
$$\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\ast }=\mathbf{U}\mathbf{\Sigma }{\tilde{\mathbf{V}}}^{\ast }{\mathbf{W}}^{\ast }=\mathbf{Z}{\mathbf{W}}^{\ast }=\mathbf{X}{\mathbf{R}}^{\ast }{\mathbf{W}}^{\ast }=\mathbf{X}{\mathbf{A}}_{(J,:)}$$
define the approximation
$$\hat{A}=Q{Q}^{\ast }A.\text{\hspace{1em}(11)}$$
Since $\hat{A}=XQ(J,:)Q^{\ast }A$ and since $X(J,:)=I_k$, it must be that $\hat{A}(J,:)=Q(J,:)Q^{\ast }A$. Consequently,
$$\hat{A}=X\hat{A}(J,:).$$
We have the chain of relations
$$\begin{array}{ll}\Vert A-U\Sigma V^{\ast }\Vert & =\Vert A-XA(J,:)\Vert \\ & =\Vert (A-X\hat{A}(J,:))+(X\hat{A}(J,:)-XA(J,:))\Vert \\ & \le \Vert A-\hat{A}\Vert +\Vert X\hat{A}(J,:)-XA(J,:)\Vert \\ & \le \Vert A-\hat{A}\Vert +\Vert X\Vert \Vert A(J,:)-\hat{A}(J,:)\Vert .\end{array}\text{\hspace{1em}(12)}$$
Inequality (6) ensures that $\Vert A-\hat{A}\Vert \le ϵ$. Since $A(J,:)-\hat{A}(J,:)$ is a submatrix of $A-\hat{A}$, we must also have $\Vert A(J,:)-\hat{A}(J,:)\Vert \le ϵ$. Thus, (12) reduces to
$$\Vert A-U\Sigma V^{\ast }\Vert \le (1+\Vert X\Vert )ϵ.\text{\hspace{1em}(13)}$$
The bound (10) follows from (13) after we observe that $X$ contains a $k\times k$ identity matrix and that the entries of the remaining $(n-k)\times k$ submatrix are bounded in magnitude by two.

Note: Recall that the orthonormal matrix $Q$ is contruct by orthonormalizing the columns of the sample matrix $Y$($Y=A\Omega$). With finite-precision arithmetric, it is preferable to adapt the above algorithm to strat directly from the matrix $Y$. To be precise, we modify step 1 to compute $X$ and $J$ so that $Y=XY(J,:)$.

---

References:
[1] N.Halko, P.G.Martinsson, J.A.Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 2(53):217-288, 2011.
[2] H. Cheng, Z. Gimbutas, P. G. Martinsson, and V. Rokhlin, On the compression of low rank matrices, SIAM J. Sci. Comput., 26 (2005), pp. 1389–1404.
[3] M. Gu and S. C. Eisenstat, Efficient algorithms for computing a strong rank-revealing QR factorization, SIAM J. Sci. Comput., 17 (1996), pp. 848–869.