Error bounds in signal compression represented by truncated Moore-Penrose biorthogonal bases using von Neumann wavelets

Question

I was reading and trying to reproduce the results in the arXiv preprint of Periodic Gabor Functions with Biorthogonal Exchange: A Highly Accurate and Efficient Method for Signal Compression by Asaf Shimshovitz et al., and a few questions arose. In brief, the paper outlines a new (as of 2012) family of wavelet-like expansions (closely related to the von Neumann lattice functions) which are time-frequency localized, and I have questions regarding how badly the method can fail when they're potentially used in signal compression or quantum-chemical calculations.

Most of the details of the paper can be stripped away to give a simple linear algebra procedure, as follows. Given a vector $\mathbf{v}\in\mathbb{C}^n$ and a not necessarily linearly-independent basis $\mathbf{G}=\{\mathbf{g}_1,\mathbf{g}_2,...,\mathbf{g}_n\}$ of column vectors spanning $\mathbb{C}^n$, the following identity holds: $$\mathbf{B}\mathbf{G}^\dagger \mathbf{v}=\mathbf{v}$$ where $\mathbf{B}=(\mathbf{G}^\dagger)^{-1}$. In Dirac notation, this trivial statement reads as $$\sum_{j=1}^n|\mathbf{b}_j\rangle\langle\mathbf{g}_j|\mathbf{v}\rangle=\sum_{j=1}^n\beta_j|\mathbf{b}_j\rangle=|\mathbf{v}\rangle$$ which can be interpreted as saying that given any vector $\mathbf{v}$ and basis $\mathbf{G}$, the expansion coefficients $\beta_j$ of $\mathbf{v}$ in the basis $\mathbf{B}$ biorthogonal to $\mathbf{G}$ are simply given by the inner products $\langle\mathbf{g}_j|\mathbf{v}\rangle$.

With a clever choice of $\mathbf{G}$, one often finds that the majority of the inner products $\langle\mathbf{g}_j|\mathbf{v}\rangle$ nearly vanish, and one is naively tempted to approximate $\mathbf{v}$ as $$|\mathbf{v}\rangle\approx\sum_{\substack{j=1 \\ \langle\mathbf{g}_j|\mathbf{v}\rangle\geq\delta}}^n|\mathbf{b}_j\rangle\langle\mathbf{g}_j|\mathbf{v}\rangle$$ where $\delta$ is some compression threshold. In matrix notation, this reads $$\mathbf{v}\approx\mathbf{B}_{tr}(\mathbf{G}_{tr})^\dagger\mathbf{v}$$ where $\mathbf{G}_{tr}=\{\mathbf{g}_j|\langle\mathbf{g}_j|\mathbf{v}\rangle\geq\delta\}$ is the basis $\mathbf{G}$ truncated to only include the basis columns which have considerable overlap with $\mathbf{v}$, and $\mathbf{B}_{tr}$ is the corresponding truncated version of $\mathbf{B}$.

This approximation is decent but not optimal; the optimal approximation is actually $$\mathbf{v}\approx\mathbf{B}_{tr}'(\mathbf{G}_{tr})^\dagger\mathbf{v}$$ where $\mathbf{B}_{tr}'=(\mathbf{G}_{tr}^\dagger)^+$ where $+$ denotes the Moore-Penrose pseudoinverse; a proof of optimality is simple and relies on basic facts about pseudoinverses, so I'll skip it. The use of $\mathbf{B}_{tr}'$ instead of $\mathbf{B}_{tr}$ is referred to as the Porat (Genossar) correction in the paper.

Simply put, my question is:

• Are there any reliable bounds on the error inherent in this approximation, given knowledge of the coefficients $\langle\mathbf{g}_j|\mathbf{v}\rangle$ and the compression threshold $\delta$? For example, if we choose $\delta$ so that $\frac{\sum_{\langle\mathbf{g}_j|\mathbf{v}\rangle\geq\delta}|\langle\mathbf{g}_j|\mathbf{v}\rangle|^2}{\sum_{j=1}^n|\langle\mathbf{g}_j|\mathbf{v}\rangle|^2}=0.99$, can we say that the truncation-compressed vector is 99% close to the original? Or is it much worse than that? Or can we say nothing? If we can say nothing, is there any further info which would allow us to say something useful?

The paper gives no statements regarding error bounds, which I assume the authors had set aside for later to figure out; this is probably largely due to the fact that the original motivation for the von Neumann wavelet procedure was for quantum mechanical calculations, rather than signal compression applications. Nevertheless, knowing bounds on how reliable the technique is seems like an essential facet of any compression method. For now, let's assume that the basis vectors $\mathbf{g}_j$ of $\mathbf{G}$ are normalized to one ($\mathbf{g}_j^\dagger\mathbf{g}_j=1$). My first thoughts on attacking this were that the error can be quantified by the expression $$\epsilon:=\frac{|\mathbf{v}-\mathbf{B}_{tr}'\mathbf{G}_{tr}^\dagger\mathbf{v}|^2}{|\mathbf{v}|^2}=\frac{|(I-\mathbf{A}^+\mathbf{A})\mathbf{v}|^2}{|\mathbf{v}|^2}=\frac{|P_{R(\mathbf{G}_{tr})}\mathbf{v}|^2}{|\mathbf{v}|^2}$$ where $\mathbf{A}=\mathbf{G}_{tr}^\dagger$ and $P_{R(\mathbf{G}_{tr})}$ is the projector operator onto the subspace spanned by the columns of $\mathbf{G}_{tr}$. But I'm stuck here, and not sure where to go (I'm not particularly good at math). Does anyone have any helpful clues or intuition?

As an unrelated but interesting bit of visual imagery, here is a picture of the signal $\frac{k^2}{67108864}+\sin \left(\frac{k^3}{64000000}-\frac{k^2}{4900}+\frac{k}{4}\right) \cos \left(\frac{k}{3}\right)$ for $1\leq k\leq 16384$ in the von Neumann domain (brightness corresponds to absolute value of coefficient and color corresponds to phase angle in complex plane, with red being 1, green being $i$, cyan being -1, and purple being $-i$, with the time axis on bottom and the frequency axis on the side):

Answer 1

[I am part of this research group, so I know this issue in-depth]

There is absolutely no approximation in the method. I'll start at the very beginning:

Start with some finite Hilbert basis, $\mathcal{H}$. In Asaf's paper this is a Fourier grid - the space of band-limited cyclic functions. The Fourier grid prescribes a rectangular area in phase-space, which may be spanned by the spectral basis or the DVR/pseudo-spectal basis (the latter is similar to the periodic-sinc functions, but not identical to it). Once this is done, we can forget about the continuous functions, and switch to a discrete representation, using the coefficients of the pseudo-spectral basis (i.e. the sampling values at the Fourier grid points).

Let us say we have $N$ such Fourier grid points. Any set of $N$ linearly independent vectors, i.e. any non-singular matrix $G$ can serve as a (generally non-orthogonal) basis to the Fourier grid. Meaning any function within the Hilbert space of the Fourier grid, $\mathcal{H}$, can be accurately spanned by vectors from $G$.

From here we proceed as you described, defining a bi-orthogonal partner to $G$, $B:=G^{-\dagger}$. By construction $B$ also spans $\mathcal{H}$.

Given any vector $v\in\mathcal{H}$, we can expand it exactly in the $B$ basis with coefficients $G^\dagger v$. If we choose $G$ to be composed of well-localized functions, then there is a good chance the coefficient vector will be sparse, allowing efficient computation. The von Neumann grid is one such choice of a basis $G$.

Where does the issue of approximation come in? Only when $v$ is not in $\mathcal{H}$. Then it has to be projected into the space, with loss of information.

Note that once you've described $v$ by the sampling values at the Fourier grid points, you have already "projected" it into the space (this is not a proper projection, which must be performed by integration with the DVR functions).