Interpolative cur and nystrom decompositions

Input a finite real or complex matrix $A\in {\mathbb{F}}^{m\times n}$ with $m\ge n\ge 1$, an integer $k$ with $0\le k\le n$, a factor $f\ge 1$, a finite swap cap $S\in \mathbb{N}\cup \{0\}$, and an arithmetic mode: exact arithmetic or finite precision. Finite-precision inputs and the arithmetic format have the represented-dyadic meaning specified by the called QR algorithm; all exact certificate calculations below use unbounded-size rational arithmetic on those represented values. In that mode supply the called QR procedure's finite real format, rounding rule and optional certificate-pass budget $C\in {\mathbb{N}}_0$ (default $1$); forward them unchanged on the interior call. Its exact unitary $W$ is reconstructed from the final sweep's dyadic reflector record, and its score and residual bounds may be symbolic nonnegative square roots of exact rationals. The direct ID check below uses only rational candidate entries and does not require evaluating those square roots. The procedure returns an algorithm status in success, rank-deficient-input, singular-selected-block, swap-budget-exhausted, or range-failure. On a successful branch it returns the selected index list $J$ and coefficient matrix $T$; finite-precision outputs are candidates until their certification status is checked. For $1\le k<n$, it also returns the RRQR diagnostics specified in step 4, including the rank status, selected-block conditioning data, and score data. In finite precision it returns separately the underlying RRQR certificate status and certificate data, and an ID-certificate status for the returned $J,T$.

1. If $k=0$, return success, $J=\mathrm{\varnothing }$, and the unique $0\times n$ matrix $T$. Then $A(:,J)T=0$. Set the block, rank, and conditioning diagnostics to not applicable and both certificate statuses to not-required. This branch includes $A=0$.
2. If $k=n$, return success, $J=(1,\dots ,n)$, and $T=I_n$. Then $A(:,J)T=A$. Set the block, rank, and conditioning diagnostics to not applicable and both certificate statuses to not-required.
3. If $1\le k<n$, call Strong rank-revealing QR with finite swaps and explicit statuses on $(A,k,f,S)$ in the requested arithmetic mode. In exact arithmetic, propagate each nonsuccess status without returning $J,T$. In finite precision, propagate each nonsuccess status, including range-failure, without returning $J,T$; retain any partial diagnostics supplied by the called procedure. On either mode's nonsuccess branch set both certificate statuses to not-applicable and return no ID-certificate data.
4. On success, let $\mathrm{\Pi }$ be the returned column permutation. In exact arithmetic retain the full factorization $A\mathrm{\Pi }=QR$ with unitary $Q$ and partition its upper-trapezoidal factor after its first $k$ rows and columns as $R=\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right]\mathrm{.}$ In finite precision use the returned candidate factors $\widehat{Q},\widehat{R}$ and the corresponding blocks; the exact unitary transformation representation $W$ remains part of the certificate data in step 5. Let $J$ be the ordered list of original column indices in the first $k$ positions of $\mathrm{\Pi }$. Set $X=R_{11}^{-1}{R}_{12}$ in exact arithmetic, or use the returned candidate $X$ in finite precision, and define $T$ by $T\mathrm{\Pi }=[I_{k}X]$. Retain for the final output $J,T$, the full factor ($Q$ or $\widehat{Q}$), and the diagnostics $(\mathrm{\Pi },R_{11},R_{12},R_{22},X,({\omega }_i{)}_i,({\gamma }_j{)}_j,(s_{ij}{)}_{ij},{\max }_{i,j}\mathrm{\mid }{X}_{ij}\mathrm{\mid },{\max }_{i,j}{s}_{ij})$, where ${\omega }_i=\mathrm{\parallel }{e}_i^{\ast }{R}_{11}^{-1}{\mathrm{\parallel }}_2$, ${\gamma }_j=\mathrm{\parallel }{R}_{22}(:,j){\mathrm{\parallel }}_2$, and $s_{ij}=(\mathrm{\mid }{X}_{ij}{\mathrm{\mid }}^2+({\omega }_i{\gamma }_j{)}^2{)}^{1\mathrm{/}2}$. In finite precision these are candidate values; retain validated enclosures only when supplied by the certificate. Continue to step 5 before returning.
5. In exact arithmetic, the returned decomposition and diagnostics refer to the exact successful factorization; set both certificate statuses to not-required. In finite precision, retain the called procedure's RRQR certificate status and any completed data $(W,(b_{ij}),\epsilon )$, including its dyadic reflector record and exact squared quantities, unchanged and independently of the following direct ID check; fields not computed on an uncertified branch remain absent. If a candidate $X$ or other required output is nonfinite, return range-failure without $J,T$, with both certificate statuses not-applicable. Otherwise form $T$ by copying the candidate entries and inserting exact zeros and ones, so no new rounding changes $X$. For the actual returned $J,T$, compute in exact rational complex arithmetic $D=A-A(:,J)T,r^2={\sum }_{a=1}^m{\sum }_{b=1}^n\mathrm{\mid }{D}_{ab}{\mathrm{\mid }}^2,c^2={\max }_{a,b}\mathrm{\mid }{T}_{ab}{\mathrm{\mid }}^2\mathrm{.}$ Check the dimensions, distinct valid indices in $J$, the entrywise identity $T(:,J)=I_k$, and $c^2\le f^2$. Set the ID-certificate status to certified exactly when these checks pass, and otherwise to uncertified. Return the exact rational data $(D,r^2,c^2)$ and the check outcomes in either case. The residual bound is the symbolic nonnegative square root $\sqrt{r^2}$, so its reporting requires no floating-point square root or range conversion. A certified ID therefore obeys $T(:,J)=I_k,{\max }_{a,b}\mathrm{\mid }{T}_{ab}\mathrm{\mid }\le f,\mathrm{\parallel }A-A(:,J)T{\mathrm{\parallel }}_2\le \sqrt{r^2}\mathrm{.}$ This direct certificate is independent of the RRQR certificate: it can certify the actual interpolation and this explicit residual bound even when RRQR score certification is unavailable. It does not certify the RRQR singular-value comparison or replace its score and unitary-reconstruction conditions. No iterative certificate search or rounding-based equality test is used. Return success with these outputs and the two separately assigned certificate statuses.

The error theorem below concerns exact-arithmetic success. It gives no finite-precision error bound for a returned candidate; a strong singular-value-relative bound additionally requires verified RRQR score and factorization information. Step 5 supplies its separate direct bound for the actual represented output, without asserting that it equals the exact-arithmetic trailing-block bound.

Let $M\in {\mathbb{F}}^{p\times q}$ be any finite real or complex matrix, and let $M^{†}$ be its unique Moore-Penrose inverse. Then $P=M{M}^{†}$ and $Q=M^{†}M$ are self-adjoint idempotents. The range of $P$ is $\mathrm{Ran}(M)$, and the range of $Q$ is $\mathrm{Ran}(M^{†})$. Thus both are orthogonal projections. If $M$ has full row rank, then $M{M}^{†}=I_p$.

Let $A\in {\mathbb{F}}^{m\times n}$ over $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$ with $m,n\ge 1$, and let $I=(i_1,\dots ,i_k)$ and $J=(j_1,\dots ,j_k)$ be ordered lists of distinct row and column indices. For $k=0$, both lists and all selected matrices have zero size, with ordinary empty-matrix multiplication. Define $C=A(:,J)\in {\mathbb{F}}^{m\times k}$, $R=A(I,:)\in {\mathbb{F}}^{k\times n}$, and the cross matrix $W=A(I,J)\in {\mathbb{F}}^{k\times k}$, using the same orders in $I,J,W$.

The skeleton CUR core is $U_{\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{l}}=W^{†}$ and its reconstruction is ${\widehat{A}}_{\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{l}}=C{W}^{†}R$. The projected least-squares CUR core is $U_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}=C^{†}A{R}^{†}$ and its reconstruction is ${\widehat{A}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}=C(C^{†}A{R}^{†})R=(C{C}^{†})A(R^{†}R)$.

For $k=0$, both reconstructions are the zero $m\times n$ matrix.

The Moore-Penrose theorem in the dependencies guarantees that these formulas are single-valued for singular as well as nonsingular $C,R,W$. The two core constructions may yield the same matrix: for $A=[1]$ and $I=J=(1)$, both cores equal $[1]$. The projection residual bound on the A page applies only to ${\widehat{A}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}$; it is not asserted for the skeleton core $W^{†}$.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, let $A\in {\mathbb{F}}^{n\times n}$ be Hermitian positive semidefinite with $n\ge 1$, and let $S\in {\mathbb{F}}^{n\times s}$ with $s\ge 0$ be any finite sketch matrix. Put $W=S^{\ast }AS\in {\mathbb{F}}^{s\times s}$. The Nyström approximation associated with $(A,S)$ is ${\widehat{A}}_{\mathrm{n}\mathrm{y}\mathrm{s}}=AS{W}^{†}{S}^{\ast }A$.

The Moore-Penrose inverse is defined for every finite matrix, including singular $W$, and is unique by the dependencies. No rank, sampling, or probabilistic condition is part of this definition.

For every admissible input to Strong-RRQR interpolative decomposition with explicit statuses, in either exact arithmetic or its specified represented finite-precision mode, the procedure terminates with exactly one listed algorithm status. The interior branch makes at most $S$ swaps and at most $S+1$ score-test passes. In finite precision the reported input rank is the exact rank of the represented matrix, not a numerical-rank estimate for unrepresented data. Every certificate status is assigned in finite work, without a search until success.

On a finite-precision successful interior branch, the direct certificate data $(D,r^2,c^2)$ refer to the actual returned $J,T$ and satisfy

$D=A-A(:,J)T,r^2={\sum }_{a,b}\mathrm{\mid }{D}_{ab}{\mathrm{\mid }}^2,c^2={\max }_{a,b}\mathrm{\mid }{T}_{ab}{\mathrm{\mid }}^2\mathrm{.}$

An ID-certificate status certified guarantees

$T(:,J)=I_k,{\max }_{a,b}\mathrm{\mid }{T}_{ab}\mathrm{\mid }\le f,\mathrm{\parallel }A-A(:,J)T{\mathrm{\parallel }}_2\le \sqrt{r^2}\mathrm{.}$

This direct bound needs no successful RRQR certificate. It asserts no singular-value-relative error bound for an uncertified RRQR candidate.

Let $A\in {\mathbb{F}}^{m\times n}$ over $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, with $m\ge n\ge 1$, and let $0\le k\le n$, $f\ge 1$. The spectral norm is $\mathrm{\parallel }M{\mathrm{\parallel }}_2={\sup }_{\mathrm{\parallel }x{\mathrm{\parallel }}_2=1}\mathrm{\parallel }Mx{\mathrm{\parallel }}_2$, and ${\sigma }_j(M)$ is the $j$th singular value in decreasing order. The procedure Strong-RRQR interpolative decomposition with explicit statuses returns the exact zero approximation for $k=0$ and the exact reconstruction for $k=n$. For $1\le k<n$, assume $\mathrm{rank}(A)=n$ and its exact-arithmetic strong-RRQR call returns success, and write $g=(1+f^{2}k(n-k){)}^{1\mathrm{/}2}$ for the returned blocks. Then $T(:,J)=I_k$, ${\max }_{i,j}\mathrm{\mid }{T}_{ij}\mathrm{\mid }\le f$, and

$\mathrm{\parallel }A-A(:,J)T{\mathrm{\parallel }}_2=\mathrm{\parallel }{R}_{22}{\mathrm{\parallel }}_2\le g{\sigma }_{k+1}(A)\mathrm{.}$

No decomposition or error conclusion is asserted when the procedure returns a rank-loss or budget status. This theorem gives no finite-precision error bound for a returned candidate; such a bound requires verified score bounds and a reconstruction certificate for the actual returned interpolative decomposition.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, let $m,n\ge 1$, and let $A\in {\mathbb{F}}^{m\times n}$. Choose an integer $0\le k\le \min (m,n)$ and ordered lists $I$ and $J$ of $k$ distinct valid row and column indices, respectively. Let $C=A(:,J)$, $R=A(I,:)$, and $W=A(I,J)\in {\mathbb{F}}^{k\times k}$ as in Skeleton and projected least-squares CUR cores. If $\mathrm{rank}(A)=\mathrm{rank}(W)=k$, then the skeleton CUR reconstruction satisfies $C{W}^{†}R=A$. For $k=0$, use the empty-matrix convention: a rank-zero $A$ is zero and the empty product reconstructs the zero $m\times n$ matrix.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, $m,n\ge 1$, $A\in {\mathbb{F}}^{m\times n}$, and choose an integer $0\le k\le \min (m,n)$ and ordered lists $I,J$ of $k$ distinct valid row and column indices, respectively. Set $C=A(:,J)$ and $R=A(I,:)$ as in Skeleton and projected least-squares CUR cores, and define $P_C=C{C}^{†}$ and $P_R=R^{†}R$. The projected least-squares CUR approximation is ${\widehat{A}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}=C(C^{†}A{R}^{†})R=P_{C}A{P}_R$. With $\mathrm{\parallel }X{\mathrm{\parallel }}_F=({\sum }_{i,j}\mathrm{\mid }{X}_{ij}{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$,

$\mathrm{\parallel }A-{\widehat{A}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}{\mathrm{\parallel }}_F\le \mathrm{\parallel }(I_m-P_C)A{\mathrm{\parallel }}_F+\mathrm{\parallel }A(I_n-P_R){\mathrm{\parallel }}_F\mathrm{.}$

This bound is for the projected least-squares core only; it is not asserted for the skeleton core $W^{†}$.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, $n\ge 1$, $s\ge 0$, $A\in {\mathbb{F}}^{n\times n}$ be Hermitian positive semidefinite, $S\in {\mathbb{F}}^{n\times s}$, and $B=A^{1\mathrm{/}2}S$. If $P_{\mathrm{Ran}(B)}$ denotes the orthogonal projection onto the range of $B$, then

${\widehat{A}}_{\mathrm{n}\mathrm{y}\mathrm{s}}=A^{1\mathrm{/}2}{P}_{\mathrm{Ran}(B)}{A}^{1\mathrm{/}2}=A^{1\mathrm{/}2}B(B^{\ast }B{)}^{†}{B}^{\ast }{A}^{1\mathrm{/}2}=A^{1\mathrm{/}2}B{B}^{†}{A}^{1\mathrm{/}2}\mathrm{.}$

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, $n\ge 1$, $s\ge 0$, $A\in {\mathbb{F}}^{n\times n}$ be Hermitian positive semidefinite, and $S\in {\mathbb{F}}^{n\times s}$. In the Loewner order, where $B⪯C$ means $C-B$ is positive semidefinite, the Nyström approximation in The positive-semidefinite Nyström approximation satisfies $0⪯{\widehat{A}}_{\mathrm{n}\mathrm{y}\mathrm{s}}⪯A$. If the sketch captures the range in the sense that $\mathrm{Ran}(AS)=\mathrm{Ran}(A)$, then ${\widehat{A}}_{\mathrm{n}\mathrm{y}\mathrm{s}}=A$. Here rank, when used in the proof, means dimension of the range.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, $n\ge 1$, $s\ge 0$, and let $A\in {\mathbb{F}}^{n\times n}$ be Hermitian positive semidefinite with eigenvalues ${\lambda }_1\ge \cdots \ge {\lambda }_n\ge 0$. Let $S\in {\mathbb{F}}^{n\times s}$, and choose an integer $0\le k\le n$ and an orthonormal eigenbasis $U=(U_1{U}_2)$ with the first $k$ eigenvectors in $U_1$. Set ${\mathrm{\Lambda }}_1=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_k)$, ${\mathrm{\Lambda }}_2=\mathrm{diag}({\lambda }_{k+1},\dots ,{\lambda }_n)$, ${\mathrm{\Omega }}_1=U_1^{\ast }S$, and ${\mathrm{\Omega }}_2=U_2^{\ast }S$. Empty blocks have norm zero, ${\lambda }_{n+1}:=0$, and the $0\times s$ matrix has full row rank vacuously. If ${\mathrm{\Omega }}_1$ has full row rank, then, deterministically,

$\mathrm{\parallel }A-{\widehat{A}}_{\mathrm{n}\mathrm{y}\mathrm{s}}{\mathrm{\parallel }}_2\le {\lambda }_{k+1}+\mathrm{\parallel }{\mathrm{\Lambda }}_2^{1\mathrm{/}2}{\mathrm{\Omega }}_2{\mathrm{\Omega }}_1^{†}{\mathrm{\parallel }}_2^2\mathrm{.}$

Here $\mathrm{\parallel }M{\mathrm{\parallel }}_2={\sup }_{\mathrm{\parallel }x{\mathrm{\parallel }}_2=1}\mathrm{\parallel }Mx{\mathrm{\parallel }}_2$ is the induced Euclidean operator norm. No probability assumption on $S$ is made.