On Low Rank Matrix Approximations with Applications to Synthesis Problem in Compressed Sensing

Introduction

We call the matrix $A$ $s$-good if whenever the true signal $x$ is $s$-sparse (i.e., has at most $s$ nonzero entries) and there is no observation errors ($\epsilon=0$), $x$ is the unique optimal solution to the optimization program $\min\{\|w\|_{1}:Aw=Ax\}.$

To the best of our knowledge, nearly the strongest verifiable sufficient condition for $A$ to be $s$-good is as follows (cf ):

(here and in what follows $\|X\|_{\infty}=\max\limits_{i,j}|X_{ij}|$, $X_{ij}$ being the elements of $X$). We address the reader to for details concerning the derivation, the link to the necessary and sufficient condition of $s$-goodness and its comparison to traditional non-verifiable sufficient conditions for $s$-goodness based on Restricted Isometry or Restricted Eigenvalue Property and a verifiable sufficient condition based on mutual incoherence.

In this paper we consider the synthesis problem of Compressed Sensing as follows:

Given $s$ and an $M\times n$ matrix $A$, extract from it an $m\times n$ submatrix $A_{m}$, certified to be $s$-good, with $m$ as small as possible.

Then we can approach the synthesis problem as follows:

Choosing $\epsilon<{1\over 2s}-\mu$ and invoking (2), we ensure that the output $A_{m}$ of the above procedure is $s$-good. This simple observation motivates our interest to the problem of approximating a given matrix by a matrix of specified (low rank) in the uniform norm.

Note that in the existing literature on low rank approximation of matrices the emphasis is on efficient construction when the approximation error is measured in the Frobenius norm (for the Frobenius norm $\|A\|_{F}=\left(\sum_{i,j}A_{ij}^{2}\right)^{1/2}$). Though the Singular Value Decomposition (SVD) gives the best rank $k$ approximation in terms of all the norms that are invariant under rotation (e.g., the Frobenius norm and the spectral norm), its computational cost may be prohibitive for applications involving large matrices. Recently, the properties of fast low rank approximations in the Frobenius norm based on the randomized sampling of rows (or columns) of the matrix (see, e.g., ) or random sampling of a few individual entries (see and references therein) has been studied extensively. Another randomized fast approximation based on the preprocessing by the Fast Fourier Transform or Fast Hadamard Transform has been studied in . Yet we do not know explicit bounds available from the previous literature which concern numerically efficient low rank approximations in the uniform norm.

In this work, we aim at developing efficient algorithms for building low rank approximation of a given matrix in the uniform norm. Specifically, we consider two types of low rank approximations:

Let $W=Y^{T}A$, where $Y$ and $A$ are known $M\times n$ matrices. We consider the approximation $W_{k}=Y_{k}A^{T}_{k}$ of $W$ such that the matrices $Y_{k}$ and $A_{k}$ of dimension $m_{k}\times n$, $m_{k}\leq k\leq M$, are composed of multiples of the rows of the matrices $Y$ and $A$ respectively. We show that a fast (essentially, of numerical complexity $O(kMn^{2})$) approximation $W_{k}$ can be constructed which satisfies

where $L(Y,A)=\sum_{i}\|y_{i}\|_{\infty}\|a_{i}\|_{\infty}$ and $y_{i}^{T},a_{i}^{T}$ denote the $i$-th rows of $Y$ and $A$ respectively. Note that for moderate values of $L(Y,A)=O(1)$ and $k<n/2$ this approximation is “quasi-optimal”, as we know (cf, e.g. [5, Proposition 4.2]) that (for certain matrices $W$) the accuracy of such an approximation cannot be better than $O(k^{-1/2})$.

where $D$ is the maximal Euclidean norm of rows of $M$ and $N$. We show that when $A$ is an $n\times n$ identity matrix the above bound is unimprovable up to a logarithmic factor.

In this paper we propose two types of construction of fast approximations: we consider the randomized construction, for which the accuracy bounds above hold in expectation (or with significant probability). We also supply “derandomized” versions of the approximation algorithms which does not require random sampling of matrices and attains the same accuracy bounds as the randomized method.

Low rank approximation in Compressed Sensing

In this section we suppose to be given $s$ and an $M\times n$ matrix $A$ and our objective is to extract from $A$ a submatrix $A_{k}$ which is composed of, at most, $k$ rows of $A$, with as small $k$ as possible, which is $s$-good. We assume that $A$ admits a “goodness certificate” $Y$. Namely, we are given an $M\times n$ matrix $Y$ such that

and we are looking for $A_{k}$ and the corresponding $Y_{k}$ such that $\|I_{n}-Y_{k}^{T}A_{k}\|<{1\over 2s}$.

The starting point of our developments is the following simple

we have $\|z\|_{\infty}-\beta\ln[2d]\leq V_{\beta}(z)\leq\|z\|_{\infty}$;

if $\beta_{1}\leq\beta_{2}$ then $V_{\beta_{1}}(z)\geq V_{\beta_{2}}(z)$;

Lemma 2.1 has the following immediate consequence:

Proof. Let $\beta\geq\beta^{\prime}$. By applying items (ii) and (iii) of the lemma we get:

When taking the expectation (first conditional to $\xi_{1},...,\xi_{k-1}$), due to ${\mathbf{E}}\{\xi_{k}|\xi_{1},...,\xi_{k-1}\}=0$ a.s., we obtain

which is (6). Now let us set $\beta^{\prime}=\beta=\sqrt{\sum_{i=1}^{k}\sigma^{2}_{i}\over 2\ln[2d]}$. Since $V_{\beta}(0)=0$ we conclude that

On the other hand, by item (i) of Lemma 2.1,

Denoting $y^{T}_{i}$ and $a^{T}_{i},\,i=1,...,M$, $i$-th rows of $Y$ and $A$, respectively, let us set

Now let $\Xi$ be random rank 1 matrix taking values $z_{i}a_{i}^{T}$ with probabilities $\pi_{i}$, and let $\Xi_{1},\Xi_{2},...$ be a sample of independent realizations of $\Xi$. Consider the random matrix

Then $W_{k}$ is, by construction, of the form $Y_{k}^{T}A_{k}$, where $A_{k}$ is a random $m_{k}\times n$ submatrix of $A$ with $m_{k}\leq k$.

As an immediate consequence of Proposition 2.1 we obtain the following statement:

In particular, the probability of the event

is $\geq 1/2$, and whenever this event takes place, we have in our disposal a matrix $Y_{k}$ and a $m_{k}\times n$ submatrix $A_{k}$ of $A$ with $m_{k}\leq k$ such that

Discussion.

2 Derandomization

Looking at the proof of Proposition 2.1, we see that the construction of $A_{k}$ and $Y_{k}$ can be derandomized. Indeed, (6) implies that

Specifically, the above bound is satisfied for every $i$ such that

and because $\pi_{i}\geq 0$ and $\sum_{i}\pi_{i}(z_{i}a_{i}^{T}-W)=0$, the latter inequality is certainly satisfied for some $i$.

Now assume that given a sequence $\beta_{0}\geq\beta_{1}\geq...$ of positive reals, we build a sequence of matrices $S_{i}$ according to the following rules:

Then for every $k\geq 1$ the matrix $U_{k}=k^{-1}S_{k}$ is of the form $Y_{k}^{T}A_{k}-W$, where $A_{k}$ is a $m_{k}\times n$ submatrix of $A$ with $m_{k}\leq k$, and

Given $S_{k}$ we solve $M$ one-dimensional convex optimization problems

If the bisection algorithm is used to find $t^{*}_{i}$, solving the problem (16) for one $i$ to the relative accuracy $\epsilon$ requires $O(n^{2}\ln\epsilon^{-1})$ elementary operations. The total numerical complexity of the step of the method is $O(Mn^{2}\ln\epsilon^{-1})$.

Given $S_{k}$, we solve $M$ convex optimization problems

Note that due to the structure of $V_{\beta}$ to solve (17) it suffices to find a solution to the system

Since the equations of the system (C.) are independent, one can use bisection to find the component $u_{j}$ of the solution. Note that due to the convexity of the left-hand side of the equation in (C.), even faster algorithm of Newton family can be used. Finding a solution to the relative accuracy $\epsilon$ to each equation then requires $O(n\ln\epsilon^{-1})$ arithmetical operations, and the total complexity of solving (17) becomes $O(Mn^{2}\ln\epsilon^{-1})$.

Note that the numerical schemes of this section should be initialized with matrices $Y$ and $W=Y^{T}A$. We can do as follows:

where $\mu$ is a certain fraction of ${1\over 2s}$. Assuming the problem is feasible for the chosen $\mu$, we get in this way the “initial point” – the matrix $W=Y^{T}A$.

3 Numerical illustration

Here we report on preliminary numerical experiments with the synthesis problem as posed in the introduction. In our experiment, $A$ is square, specifically, this is the Hadamard matrix $H_{11}$ of order 2048.

Recall that the Hadamard matrix $H_{\nu}$, $\nu=0,1,...$ is a square matrix of order $2^{\nu}$ given by the recurrence

whence $H_{\nu}$ is a symmetric matrix with entries $\pm 1$ and $H_{\nu}^{T}H_{\nu}=2^{\nu}I_{2^{\nu}}$.

The goal of the experiment was to extract from $A=H_{11}$ an $m\times 2048$ submatrix $A_{m}$ which satisfies the relation (cf. (1))

with $s=10$; under this requirement, we would like to have $m$ as small as possible. In Compressed Sensing terms, we are trying to solve the synthesis problem with $A=H_{11}$; in low rank approximation terms, we want to approximate $I_{2048}$ in the uniform norm within accuracy $<0.05$ by a rank $m$ matrix of the form $Y_{m}^{T}A_{m}$, with the rows of $A_{m}$ extracted from $H_{11}$. The advantages of the Hadamard matrix in our context is twofold:

The error bound (13) is proportional to the quantity $L$ defined in (8). By the origin of this quantity, we clearly have $\|Y^{T}A\|_{\infty}\leq L$, whence $L\geq 1-\mu>1-{1\over 2s}\geq 1/2$ by (3). On the other hand, with $A=H_{\nu}$ being an Hadamard matrix, setting $Y=2^{-n}YH_{\nu}$, so that $Y^{T}A-I_{2^{\nu}}$, we ensure the validity of (3) with $\mu=0$ and get $L=1$, that is, $\mu$ is as small as it could be, and $L$ is nearly as small as it could be.

Whenever $A_{m}$ is a submatrix of $H_{\nu}$, the optimization problem in the left hand side of (21) is easy to solve.

Item 2 deserves an explanation. Clearly, the optimization program in (21) reduces to the series of $n=2048$ LP programs

The experiment was organized as follows. As it was already mentioned, we used $\nu=11$ (that is, $n=2048$) and $s=10$ (that is, the desired uniform norm of approximating $I_{2048}$ by $Y_{m}^{T}A_{m}$ was 0.05). We compared two approximation policies:

“Active” approximation, which is obtained from algorithm A′ by the same refinement as in the previous item.

In our experiments, we ran every policy 6 times. The results were as follows: “Blind” policy ${\cal B}$: the rank of $0.05$-approximation of $W=I_{2048}$ varied from 662 to 680. “Active” policy ${\cal A}$: the rank of $0.05$-approximation of $W$ varied from 617 to 630. Note that in both algorithms the resulting matrix $A_{m}$ is built “row by row”, and the certified levels of goodness of the intermediate matrices $A^{1},A^{2},...$ are computed. In the below table we indicate, for the most successful (resulting in the smallest $m$) of the 6 runs of each algorithm, the smallest values of $k$ for which $A^{k}$ was certified to be $s$-good, $s=1,2,...,10$:

Finally, we remark that with $A$ being the Hadamard matrix $H_{\nu}$, the “no refinement” versions of our policies would terminate according to the criterion $\|I_{n}-{1\over k}A_{k}^{T}A_{k}\|_{\infty}<{1\over 2s}$, which, on a closest inspection, is nothing but a slightly spoiled version of the goodness test based on mutual incoherence The mutual incoherence test is as follows: given a $k\times n$ matrix $B=[b_{1},...,b_{n}]$ with nonzero columns, we compute the quantity $\mu(B)=\max\limits_{i\neq j}|b_{i}^{T}b_{j}|/b_{i}^{T}b_{i}$ and claim that $B$ is $s$-good for all $s$ such that $s<{1+\mu(B)\over 2\mu(B)}$. With the Hadamard $A$, the “no refinement” criterion for our scheme is nothing but $s<{1\over 2\mu(A^{k})}$.. In the experiments we are reporting, this criterion is essentially weaker that the one based on (21): for the best, over the 6 runs of the algorithms ${\cal A}$ and ${\cal B}$, 10-good submatrices $A_{m}$ of $H_{11}$ matrices we got the test based on mutual incoherence certifies the levels of goodness as low as 5 (in the case of ${\cal B}$) and 7 (in the case of ${\cal A}$).

Low rank approximation of arbitrary matrices

Given a positive integer $k$, consider the random matrix

2 The norm associated with Proposition 3.1

This relation clearly defines a norm, and one clearly has $\|A\|=\|A^{T}\|$.

The next result summarizes the basic properties of the norm we have introduced.

(i) $\|A\|_{\infty}\leq\|A\|\leq\sqrt{\min[m,n]}\|A\|_{\infty}$.

(ii) $\|A\|\leq\|A\|_{2,2}$, where $\|A\|_{2,2}$ is the usual spectral norm of $A$ (the maximal singular value).

(iii) If $A$ is symmetric positive semidefinite, then $\|A\|=\|A\|_{\infty}$.

(iv) If the Euclidean norms of all rows (or all columns) of $A$ are $\leq D$, then $\|A\|\leq D$.

3 Lower bound

When $n\geq 2k$, the $\|\cdot\|_{\infty}$ error of any approximation of the unit matrix $I_{n}$ by a matrix of rank $k$ is at least

Proof [cf. [5, Proposition 4.2]] Let $\alpha(n,k)$ be the minimal $\|\cdot\|_{\infty}$ error of approximation of $I_{n}$ by a matrix of rank $\leq k$; this function clearly is nondecreasing in $n$. Let $\nu$ be an integer such that $k<\nu\leq n$, and $A$ be an $\nu\times\nu$ matrix of rank $\leq k$ such that $\|I_{\nu}-A\|_{\infty}=\alpha:=\alpha(\nu,k)$. By variational characterization of singular values, at least $\nu-k$ singular values of $I_{\nu}-A$ are $\geq 1$, whence ${\mathop{\hbox{\rm Tr}}}([I_{\nu}-A][I_{\nu}-A]^{T})\geq\nu-k$. On the other hand, $\|I_{\nu}-A\|_{\infty}\leq\alpha$, whence ${\mathop{\hbox{\rm Tr}}}([I_{\nu}-A][I_{\nu}-A]^{T})\leq\nu^{2}\alpha^{2}$. We conclude that $\alpha^{2}\geq{\nu-k\over\nu^{2}}$ for all $\nu$ with $k<\nu\leq n$, whence $\alpha^{2}\geq{1\over 4k}$ when $n\geq 2k$. $\square$

References

Appendix A Proof of Lemma 2.1

Properties (i) and (ii) are immediate consequences of the definition of $V_{\beta}$. Observe that $V_{\beta}$ is convex and continuously differentiable with

Appendix B Problems (22) in the case of Hadamard matrix AA

These problems clearly have equal optimal values, due to

Appendix C Proof of Proposition 3.1

The reasoning to follow is completely standard. Let us fix $i$, $1\leq i\leq m$, and $j$, $1\leq j\leq n$, and let $\xi\sim{\cal N}(0,I_{d})$, $\mu=D^{-1/2}p^{T}_{i}\xi$, $\nu=D^{-1/2}q_{j}^{T}\xi$, and $\alpha=D^{-1}A_{ij}$. Then $[\mu;\nu]$ is a normal random vector with ${\mathbf{E}}\{\mu^{2}\}\leq 1$, ${\mathbf{E}}\{\nu^{2}\}\leq 1$ and ${\mathbf{E}}\{\mu\nu\}=\alpha$. We can find a normal random vector $z=[u;v]\sim{\cal N}(0,I_{2})$ such that $\mu=au$, $\nu=bu+cv$; note that $a^{2}\leq 1$, $b^{2}+c^{2}\leq 1$ and $ab={\mathbf{E}}\{\mu\nu\}=\alpha$. Note that $\mu\nu=z^{T}Bz$ with B=\left[\begin{array}[]{cc}ab&ac/2\cr ac/2&0\cr\end{array}\right]. Denoting $\lambda_{1}$, $\lambda_{2}$ the eigenvalues of $B$, we have

where the last inequality follows from $|2\gamma\lambda_{s}|\leq 1/2$, for $s=1,2$, and $-\ln(1-r)-r\leq r^{2}$ when $|r|\leq 1/2$. Using (26) we obtain,

Setting $\gamma={t\over 4k^{1/2}}$ (this results in $0<\gamma\leq 1/4$ due to $k^{1/2}\geq t$), we get

Letting $\kappa_{-}=\hbox{\rm Prob}\{k[A_{k}]_{ij}<A_{ij}k-Dk^{1/2}t\}$, we have

for all $\gamma\in(0,1/4]$, whence, same as above,

Since this relation holds true for all $i,j$, we conclude that

Appendix D Proofs for section 3.2

First claim: there exist $M,N$ such that the matrix {\cal A}=\left[\begin{array}[]{c|c}M&A\cr\hline\cr A^{T}&N\cr\end{array}\right] is positive semidefinite and has all diagonal entries, and then all entries, in $[-\|A\|,\|A\|]$. Let ${\cal A}={\cal B}{\cal B}^{T}$; then the rows in ${\cal B}$ have Euclidean norms $\leq\sqrt{\|A\|}$. Representing ${\cal B}=[P;Q]$ with $m$ rows in $P$ and $n$ rows in $Q$, the relation $[P;Q][P;Q]^{T}={\cal A}$ implies that $A=PQ^{T}$.

Second claim: Let $A=PQ^{T}$ with the Euclidean norms of rows in $P,Q$ not exceeding $\sqrt{D}$. Then 0\preceq\left[\begin{array}[]{c}P\cr Q\cr\end{array}\right]\left[\begin{array}[]{c}P\cr Q\cr\end{array}\right]^{T}=\left[\begin{array}[]{c|c}PP^{T}&A\cr\hline\cr A^{T}&QQ^{T}\cr\end{array}\right] and the diagonal entries in $M=PP^{T}$ and $N=QQ^{T}$ do not exceed $D$. $\square$

Proof of Proposition 3.3.

(i): The first inequality in (i) is evident. Let us prove the second. W.l.o.g. we can assume $\|A\|_{\infty}\leq 1$. In this case our statement reads

Assume, on the contrary, that $\hbox{\rm Opt}>D$. Since the semidefinite problem defining Opt is strictly feasible, the dual problem

By $(a)$, letting $L=\hbox{\rm Diag}\{\sqrt{\lambda_{i}}\}$, $R=\hbox{\rm Diag}\{\sqrt{\rho_{j}}\}$, we have $V=LWR$ with certain $W$, $\|W\|_{2,2}\leq 1$ ($\|\cdot\|_{2,2}$ is the usual matrix norm, the maximum singular value), thus

where the concluding $\leq$ is due to $(b)$, and $(*)$ is given by the following reasoning: w.l.o.g. we can assume that $n\leq m$. Since $W$ is of the matrix norm $\leq 1$, the columns $U_{j}$ of $U=[|W_{ij}|]_{i,j}$ satisfy $\|U_{j}\|_{2}\leq 1$, whence

The resulting inequality in (27) contradicts $(c)$; we have arrived at a desired contradiction. (i) is proved. (ii): This is evident, since \left[\begin{array}[]{c|c}\|A\|_{2,2}I_{m}&A\cr\hline\cr A^{T}&\|A\|_{2,2}I_{n}\cr\end{array}\right]\succeq 0. (iii): This is evident, since for $A\succeq 0$ we have \left[\begin{array}[]{c|c}A&A\cr\hline\cr A&A\cr\end{array}\right]\succeq 0. (iv): Since $\|A\|=\|A^{T}\|$, it suffices to consider the case when the rows of $A$ are of the norm not exceeding $D$. In this case, the result is readily given by the fact that \left[\begin{array}[]{c|c}D^{-1}AA^{T}&A\cr\hline\cr A^{T}&DI_{n}\cr\end{array}\right]\succeq 0. $\square$