Fixed-precision problem and Fixed-rank problem in Matrix low-rank approximation

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

The following factorizaitons are often used to express a low-rank approximation of a given matrix:
$$\begin{array}{cccc}\mathbf{A}& \approx & \mathbf{B}& \mathbf{C}\\ m\times n& & m\times k& k\times n\end{array}\text{\hspace{1em}(1)}$$
When the numerical rank $k$ is much smaller than either dimension $m$ or $n$, a factorization such as (1) allows the matrix to be stored inexpensively and to be multiplied rapidly with vectors or other matrices. The factorizations can also be used for data interpretation or to solve computational problems, such as least squares.

1 Fixed-precision problem

The basic challenge in producing low-rank approximations is a primitive question that are called the fixed-precision approximation problem.
Suppose we are given a matrix $A$ and a positive error tolerance $ϵ$. We seek a matrix $A$ with $k=k(ϵ)$ orthonormal columns such that
$$\Vert A-Q{Q}^{\ast }A\Vert \le \epsilon \text{\hspace{1em}(2)}$$
where $\Vert \cdot \Vert$ denotes the $l_2$ operator norm. The range of $Q$ is a $k$-dimensional subspace that captures most of the action of $A$, and we would like $k$ to be as small as possible.

The SVD furnishes an optimal answer to the fixed-precision problem. Let ${\sigma }_i$ denote the $j$th largest singular value of $A$. For each $j\ge 0$,
$$\underset{\mathrm{rank}(\mathbf{X})\le j}{\min }\Vert \mathbf{A}-\mathbf{X}\Vert ={\sigma }_{j+1}\text{\hspace{1em}(3)}$$
One way to construct a minimizer is to choose $X=Q{Q}^{\ast }A$, where the columns of $Q$ are $k$ dominant left singular vectors of $A$. Consequently, the minimal rank $k$ where (2) holds equals the number of singular values of $A$ that exceed the tolerance $ϵ$.

To simplify the development of algorithms, it is convenient to assume that the desired rank $k$ is specified in advance. We call the resulting problem the the fixed-rank approximaiton problem.

---

2 Fixed-rank problem

Given a matrix $A$, a target rank $k$, and a oversampling parameter $p$, we seek to construct a matrix $Q$ with $k+p$ orthonormal columns such that
$$\Vert A-Q{Q}^{\ast }A\Vert \approx \underset{\mathrm{rank}(X)\le k}{\min }\Vert A-X\Vert .\text{\hspace{1em}(4)}$$
Although there exists a minimizer $Q$ that solves the fixed-rank problem for $p=0$, the opportunity to use a small number of additional columns provides a flexibility that is crucial for the effectiveness of the computational methods we discuss.

---

3 Adapte algorithms for Fixed-rank to Fixed-precision

3.1 Algorithm for fixed rank approximation

Algorithms for the fixed-rank problem can be adapted to solve the fixed-precision problem. The connection is based on the observation that we can build the basis matrix $Q$ incrementally and, at any point in the computation, we can inexpensively estimate the residual error $\Vert A-Q{Q}^{\ast }A\Vert .$

For the fixed-rank approximation problem, we can use the following random projeciton method,

The output matrix $Q$ the above algorithm is what we want for the fixed-rank problem, where the target rank of the input matrix is specified in advance.

To handle the fixed-precision problem, where the parameter is the computational tolerance, we need a scheme for estimating how well a putative basis matrix $Q$ captures the action of the matrix $A$ . To do so, we develop a probabilistic error estimator. These methods are inspired by work of Dixon [2]; our treatment follows [3]and [4].

3.2 One important lemma and some instructions

The exact approximation error is $\Vert (I-QQ\ast )A\Vert$. It is intuitively plausible that we can obtain some information about this quantity by computing $\Vert (I-QQ\ast )A\omega \Vert$, where $\omega$ is a standard Gaussian vector. This notion leads to the following method. Draw a sequuence { ω ( i ) : i = 1 , 2 , ⋯   , r } \{\omega^{(i)}:i=1,2,\cdots,r\} {ω(i):i=1,2,⋯,r} of standard Gaussian vectors, where $r$ is a small integer that balances computational cost and reliability. Then
$$\Vert \left(\mathbf{I}-\mathbf{Q}{\mathbf{Q}}^{\ast }\right)\mathbf{A}\Vert \le 10\sqrt{\frac{2}{\pi }}\underset{i=1,\dots ,r}{\max }\Vert \left(\mathbf{I}-\mathbf{Q}{\mathbf{Q}}^{\ast }\right)\mathbf{A}{\mathbf{\omega }}^{(i)}\Vert \text{\hspace{1em}(5)}$$
with probability at least $1-10^{-r}$.

By setting $B=(I-Q{Q}^{\ast })A$ and $\alpha =10$ in the following lemma, whose proof appears in [4,section 3.4]

Lemma 1
Let $B$ be a real $m\times n$ matrix. Fix a positive integer $r$ and a real number $\alpha >1$. Draw an independent family { ω ( i ) : i = 1 , 2 , ⋯   , r } \{\omega^{(i)}:i=1,2,\cdots,r\} {ω(i):i=1,2,⋯,r} of standard Gaussian vectors. Then
$$\Vert \mathbf{B}\Vert \le \alpha \sqrt{\frac{2}{\pi }}\underset{i=1,\dots ,r}{\max }\Vert \mathbf{B}{\mathbf{\omega }}^{(i)}\Vert$$
except with probability ${\alpha }^{-r}$.

As the error estimate (5) is computationally inexpensive because it requires only a samll number of matrix-vector products. Therefore, we can make a lowball guess for the numerical rank of $A$ and add more samples if the error estimate is too large.

• The asymptotic cost of the above algorithm (randomized range finder) is preserved if we double our guess for the rank at each step. For example, we can start with 32 samples, compute another 32, then another 64, etc.
• The estimate (5) is actually somewhat crude. We can obtain a better estimate at a similar computational cost by initializing a power iteration with a random vector and repeating the process several times[3].

3.3 Algorithm for fixed precison problem (based on the algorithm for fixed rank problem)

Let us recall the fixed-precision problem. Suppose $A$ is an $m\times n$ matrix and $ϵ$ is a computational tolerance. Seek an integer $l$ and an $m\times l$ orthonormal matrix $Q^{(l)}$ such that
$$\Vert \left(\mathbf{I}-{\mathbf{Q}}^{(\ell )}{\left({\mathbf{Q}}^{(\ell )}\right)}^{\ast }\right)\mathbf{A}\Vert \le \epsilon \text{\hspace{1em}(6)}$$
The size $l$ of the basis will typically be slightly large than the size $k$ of the smallest basis that achieves this error.
The basic observation behind the adaptive scheme is that we can generate the basis in step $3$ of above algorithm incrementally. Starting with an empty basis matrix $Q^{(0)}$, the following scheme generates an orthonormal matrix whose range captures the action of $A$:

• for $i=1,2,3,\cdots$
• • Draw an $n\times 1$ Gaussian random vector ${\omega }^{(i)}$, and set $y^{(i)}=A{\omega }^{(i)}$.
• • Compute ${\tilde{\mathbf{q}}}^{(i)}=\left(\mathbf{I}-{\mathbf{Q}}^{(i-1)}{\left({\mathbf{Q}}^{(i-1)}\right)}^{\ast }\right){\mathbf{y}}^{(i)}$.
• • Normalize ${\mathbf{q}}^{(i)}={\tilde{\mathbf{q}}}^{(i)}/\Vert {\tilde{\mathbf{q}}}^{(i)}\Vert$, and form ${\mathbf{Q}}^{(i)}=\left[{\mathbf{Q}}^{(i-1)}{\mathbf{q}}^{(i)}\right]$.
• end for

How do we know when we have reached a basis $Q^{(l)}$ that verifies (6)? The answer becomes apparent once we observe that the vectors ${\tilde{q}}^{(i)}$ are precisely the vectors that appear in the error bound (5). The resulting rule is that we break the loop once we observe $r$ consecutive vecotrs ${\tilde{q}}^{(i)}$ whose norms are smaller than $ϵ/(10\sqrt{2/pi})$.

A formal description of the resulting algorithm appears as the following algorithm.

• A potential complication of the method is that the vecotrs ${\tilde{q}}^{(i)}$ become samll as the basis starts to capture most of the action of $A$. In finite-precision arithmetic, their direction is extremely unreliable. To address this problem, we simply reproject the normalized vector $q^{(i)}$ onto the $\mathrm{range}{\left({\mathbf{Q}}^{(i-1)}\right)}^{\perp }$ in steps 7 and 8 of the above algorithm.
• The failure probability stated for the above algorithm is pessimistic because it is derived from a simple union bound argument. In practive, the error estimator is reliable in a range of circumstances when we take $r=10$.

---

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] J. D. Dixon, Estimating extremal eigenvalues and condition numbers of matrices, SIAM J.Numer. Anal., 20 (1983), pp. 812–814.
[3] E. Liberty, F. F. Woolfe, P.-G. Martinsson, V. Rokhlin, and M. Tygert, Randomized algorithms for the low-rank approximation of matrices, Proc. Natl. Acad. Sci. USA, 104 (2007), pp. 20167–20172.
[4] F. Woolfe, E. Liberty, V. Rokhlin, and M. Tygert, A fast randomized algorithm for the approximation of matrices, Appl. Comput. Harmon. Anal., 25 (2008), pp. 335–366.