Typical $l_1$-recovery limit of sparse vectors represented by concatenations of random orthogonal matrices

Introduction

The recovery problem of sparse vectors from a linear underdetermined set of equations has recently attracted attention in various fields of science and technology due to its many applications, for example, in linear regression , communication , , , multimedia , , , and compressive sampling (CS) , . In such a sparse representation problem, we have the following underdetermined set of linear equations

Recent results in the parallel problem of CS, where $\boldsymbol{D}$ acts as a sensing matrix, reveal that the typical conditions for perfect $l_{1}$-recovery are universal for all random sensing matrices that belong to the rotationally invariant matrix ensembles . The standard setup, where the entries of the sensing matrix are independent standard Gaussian, is an example that belongs to this ensemble. It is also known that the conditions required for perfect recovery do not in general depend on the details of the marginal distribution related to the non-zero elements. On the other hand, we know that correlations in the sensing matrix can degrade the performance of $l_{1}$-recovery . This suggests intuitively that using a sample matrix of the rotationally invariant ensembles as $\boldsymbol{D}$ is preferred in the recovery problem when we expect to encounter a variety of dense signals $\boldsymbol{y}$. However, the set of matrix ensembles whose $l_{1}$-recovery performance are known is still limited, and further investigation is needed to assess whether the choice of $\boldsymbol{D}$ is indeed so straightforward.

The purpose of the present study is to fulfill this demand. Specifically, we examine the typical $l_{1}$-recovery performance of the matrices constructed by concatenating several randomly chosen orthonormal bases. Such construction has attracted considerable attention due to ease of implementation and theoretical elegance , , for designing sparsity inducing over-complete dictionaries for natural signals . For a practical engineering scheme, audio coding (music source coding) uses a dictionary formed by concatenating several modified discrete cosine transforms with different parameters.

By using the replica method in conjunction with the development of an integral formula for handling random orthogonal matrices, we show that the dictionary consisting of concatenated orthogonal matrices is also preferred in terms of the performance of $l_{1}$-recovery. More precisely, the matrices can result in better $l_{1}$-recovery performance than that of the rotationally invariant matrices when the density of non-zero entries of $\boldsymbol{x}$ is not uniform among the orthogonal matrix modules, while the performance is the same between the two types of matrices for the uniform densities. This surprising result further promotes the use of the concatenated orthogonal matrices in practical applications.

This paper is organized as follows. In the next section, we explain the problem setting that we investigated. In Section 3, which is the main part of this paper, we discuss the development of a methodology for evaluating the recovery performance of the concatenated orthogonal matrices on the basis of the replica method and an integral formula concerning the random orthogonal matrices. In Section 4, we explain the significance of the methodology through application to two distinctive examples, the validity of which is also justified by extensive numerical experiments. The final section is devoted to a summary.

Problem Setting

With full knowledge of $\boldsymbol{O}_{t}$ $(t=1,2,\ldots,T)$ and $\boldsymbol{y}$, the $l_{1}$-recovery is performed by solving the constrained minimization problem

For theoretically evaluating the $l_{1}$-recovery performance, we assume that the entries of $\boldsymbol{x}^{0}_{t}$, $x_{it}^{0}$ are distributed independently according to a block-dependent sparse distribution

where $0\leq\rho_{t}\leq 1$ means the density of the non-zero entries of the $t$-th block of the same size $M$ in $\boldsymbol{x}^{0}$ and $f_{t}(x)$ is a distribution whose second moment about the origin is finite, which is assumed as unity for simplicity. Intuitively, as the compression rate $\alpha=M/N=T^{-1}$ decreases, the overall density $\rho=T^{-1}\sum_{t=1}^{T}\rho_{t}$ up to which (9) can successfully recover a typical sample of the original vector $\boldsymbol{x}^{0}$ becomes smaller. However, precise performance may depend on the profile of $\rho_{t}$. The above setting allows us to quantitatively examine how such block dependence of the non-zero density affects the critical relation between $\alpha$ and $\rho$ for typically successful $l_{1}$-recovery of $\boldsymbol{x}^{0}$.

Statistical mechanics approach

and $Z^{-1}(\beta;\{\boldsymbol{x}_{t}^{0}\},\{\boldsymbol{O}_{t}\})\equiv\int\prod_{t=1}^{T}d\boldsymbol{x}_{t}\delta\left(\sum_{t=1}^{T}\boldsymbol{O}_{t}(\boldsymbol{x}_{t}^{0}-\boldsymbol{x}_{t})\right)\exp\left(-\beta\sum_{t=1}^{T}\|\boldsymbol{x}_{t}\|_{1}\right)$, constitutes the basis of our analysis. Equations (11) and (13) mean that $\hat{\boldsymbol{x}}$ can be identified with the average of the state variable $\boldsymbol{x}$ for the Gibbs-Boltzmann distribution (13) in the vanishing temperature limit $\beta\to\infty$. However, as (13) depends on $\{\boldsymbol{x}_{t}^{0}\}$ and $\{\boldsymbol{O}_{t}\}$, further averaging with respect to the generation of these external random variables is necessary for evaluating the typical properties of the $l_{1}$-recovery. Evaluation of such “double averages” can be carried out systematically using the replica method .

2 Integral formula for handling random orthogonal matrices

Let us assume that $M$-dimensional vectors $\boldsymbol{u}_{t}=\boldsymbol{x}_{t}^{0}-\boldsymbol{x}_{t}$ $(t=1,2,\ldots,T)$ are characterized by their norms as $v_{1}=M^{-1}|\boldsymbol{u}_{1}|^{2},v_{2}=M^{-1}|\boldsymbol{u}_{2}|^{2},\ldots,v_{T}=M^{-1}|\boldsymbol{u}_{T}|^{2}$, where $|\boldsymbol{u}|$ denotes the standard Euclidean norm of the vector $\boldsymbol{u}$. For these vectors, we define the function

where $\left[f(X)\right]_{X}$ generally denotes the average of $f(X)$ with respect to $X$, and ${\cal D}\boldsymbol{O}$ denotes the Haar measure of the $M\times M$ orthogonal matrices. Our claim is that by explicitly using $v_{1},v_{2},\ldots,v_{T}$, (15) can be expressed as

where $\mathop{\rm extr}_{X}\left\{f(X)\right\}$ generally denotes the extremization of function $f(X)$ with respect to $X$. Expression (17) is derived from the fact that for fixed $\boldsymbol{u}_{t}$, $\boldsymbol{O}_{t}\boldsymbol{u}_{t}$ moves uniformly on the surface of the $M$-dimensional hypersphere of radius $\sqrt{Mv_{t}}$ when $\boldsymbol{O}_{t}$ varies according to the uniform distribution of the orthogonal matrices; therefore, ${\cal D}\boldsymbol{O}_{t}$ in (15) can be replaced with a spherical measure of $M$-dimensional vector $d\boldsymbol{u}_{t}\delta(|\boldsymbol{u}_{t}|^{2}-Mv_{t})$. For details, see A.

Function $F(\{v_{t}\})$ physically represents a characteristic exponent of the probability that $M$-dimensional vectors $\boldsymbol{u}_{t}$ $(t=1,2,\ldots,T)$ form a closed loop satisfying $\sum_{t=1}^{T}\boldsymbol{u}_{t}=\boldsymbol{0}$ when they are independently and isotropically sampled under the norm constraints of $|\boldsymbol{u}_{t}|^{2}=Mv_{t}$. For small $T$, the loop condition strongly restricts the region of $\{v_{t}\}$ to which $F(\{v_{t}\})$ is well defined. In concrete terms, $F(\{v_{t}\})$ diverges to minus infinity unless an equality

are satisfied for $T=2$ and $3$, respectively. In particular, the constraint of (18) requires us to deal with the system of $T=2$ in a manner different from that of $T\geq 3$ except for the case of $\rho_{1}=\rho_{2}$, as discussed in the analysis of section 3.4.

3 Replica method

Now, we are ready to apply the replica method for analyzing the typical property of the $l_{1}$-recovery (11). For this, we evaluate the $n$-th moment of the partition function using the identity

Let us consider averaging (21) with respect to $\{\boldsymbol{O}_{t}\}$ and define for each fixed set of $\{\boldsymbol{x}_{t}^{a}\}$:

where $\boldsymbol{u}_{t}^{a}=\boldsymbol{x}_{t}^{0}-\boldsymbol{x}_{t}^{a}$. When $\{\boldsymbol{x}_{t}^{a}\}$ is placed in the configuration of the RS solution, the expression

holds for each $t$, where $\rm{T}$ stands for matrix transpose, $R_{t}=Q_{t}-2m_{t}+\rho_{t}$, and $r_{t}=q_{t}-2m_{t}+\rho_{t}$. $\boldsymbol{E}=[\boldsymbol{e}_{1}\ \boldsymbol{e}_{2}\ \ldots\ \boldsymbol{e}_{n}]$ denotes an $n\times n$ orthogonal matrix composed of the vector $\boldsymbol{e}_{1}=(n^{-1/2},n^{-1/2},\ldots,n^{-1/2})^{\rm T}$ and an orthonormal set of $n-1$ vectors $\boldsymbol{e}_{2},\boldsymbol{e}_{3},\ldots,\boldsymbol{e}_{n}$ that are orthogonal to $\boldsymbol{e}_{1}$. This indicates that $\left[\boldsymbol{u}_{t}^{1}\ \boldsymbol{u}_{t}^{2}\ \ldots\ \boldsymbol{u}_{t}^{n}\right]$ may be expressed as

3.2 Entropic part

On the other hand, inserting identities $1=\int MdQ_{t}\delta(|\boldsymbol{x}_{t}^{a}|^{2}-MQ_{t})$, $1=\int Mdq_{t}\delta(\boldsymbol{x}_{t}^{a}\cdot\boldsymbol{x}_{t}^{b}-Mq_{t})$ and $1=\int Mdm_{t}\delta(\boldsymbol{x}_{t}^{0}\cdot\boldsymbol{x}_{t}^{a}-Mm_{t})$ $(a,b=1,2,\ldots,n;t=1,2,\ldots,T)$ into (21) and taking an average concerning $\{\boldsymbol{x}_{t}^{0}\}$, in conjunction with integration with respect to dynamical variables $\{\boldsymbol{x}_{t}^{a}\}$, result in the expression

for a fixed set of $\{Q_{t},q_{t},m_{t}\}$. The conjugate variable $\hat{Q}_{t}^{a}$ is introduced for expressing a delta function as $\delta(|\boldsymbol{x}_{t}^{a}|^{2}-MQ_{t})=(4\pi\sqrt{-1})^{-1}\int_{-\sqrt{-1}\infty}^{+\sqrt{-1}\infty}d\hat{Q}_{t}^{a}\exp\left(-(\hat{Q}_{t}^{a}/2)\left(|\boldsymbol{x}_{t}^{a}|^{2}-MQ_{t}\right)\right)$, and similarly for $\hat{q}_{t}^{ab}$ and $\hat{m}_{t}^{a}$. Notation $d\hat{\boldsymbol{Q}}$ stands for an integral measure $\prod_{t=1}\left(\prod_{a=1}d\hat{Q}_{t}^{a}\prod_{a>b}d\hat{q}_{t}^{ab}\prod_{a=1}d\hat{m}_{t}^{a}\right)$, and the functions on the right hand side are defined as

where the average for $x_{t}^{0}$ is taken according to (10). The expression of (37) indicates that its rescaled logarithm is accurately evaluated using the saddle point method with respect to the conjugate variables in the large system limit $M\to\infty$. In addition, the replica symmetry guarantees that the relevant saddle point is of the RS form as $\hat{Q}_{t}^{a}=\hat{Q}_{t}$, $\hat{q}_{t}^{ab}=\hat{q}_{t}$, and $\hat{m}_{t}^{a}=\hat{m}_{t}$. As a consequence, the evaluation yields

3.3 Free energy and saddle point equations

and rescaled variables are introduced as $\beta(Q_{t}-q_{t})\to\chi_{t}$, $\hat{q}_{t}/\beta^{2}\to\hat{\chi}_{t}$, $(\hat{Q}_{t}+\hat{q}_{t})/\beta\to\hat{Q}_{t}$, and $\hat{m}_{t}/\beta\to\hat{m}_{t}$ to properly describe the relevant solution in the limit of $\beta\to\infty$. Extremization is to be performed with respect to $\{Q_{t},\chi_{t},m_{t},\hat{Q}_{t},\hat{\chi}_{t},\hat{m}_{t}\}$.

Similar to earlier studies, at the extremum characterized by a set of the saddle point equations

$Q_{t}$ and $m_{t}$ $(t=1,2,\ldots,T)$ physically denote the macroscopic averages of the recovered vector $\hat{\boldsymbol{x}}_{t}$ as $M^{-1}\left[|\hat{\boldsymbol{x}}_{t}|^{2}\right]_{\{\boldsymbol{x}_{k}^{0}\},\{\boldsymbol{O}_{k}\}}$ and $M^{-1}\left[\boldsymbol{x}_{t}^{0}\cdot\hat{\boldsymbol{x}}_{t}\right]_{\{\boldsymbol{x}_{k}^{0}\},\{\boldsymbol{O}_{k}\}}$, respectively. Here, $R_{t}=\Lambda_{t}^{-1}/\left(\sum_{t=1}^{T}\Lambda_{t}^{-1}\right)$ is provided by the extremum solution of (17) for $\{v_{t}\}=\{\chi_{t}\}$, and

For $T=2$, a Lagrange multiplier should be exceptionally introduced in (48) for enforcing $\chi_{1}=\chi_{2}$.

The success of the $l_{1}$-recovery is characterized by the condition in which $Q_{t}=m_{t}=\rho_{t}$ is satisfied at the extremum for $\forall{t}=1,2,\ldots,T$. Therefore, one can evaluate the critical relation between $\alpha=1/T$ and $\rho=T^{-1}\sum_{t=1}^{T}\rho_{t}$ by examining the thermodynamic stability of the success solution $Q_{t}=m_{t}=\rho_{t}$ $(t=1,2,\ldots,T)$ for the saddle point equations (47)–(52).

We assume $T\geq 3$ for a while, since an exceptional treatment is required for $T=2$. For obtaining the success solution, it is necessary that $\chi_{t}=0$ and $\hat{Q}_{t}=\hat{m}_{t}=+\infty$ hold. Expanding (50)–(52) under the assumption of $\hat{Q}_{t}=\hat{m}_{t}\gg 1$ yields

and we used $[(x_{t}^{0})^{2}]_{x_{t}^{0}}=\rho_{t}$, which is derived from the assumption (10). These result in the expression

holds, where $\{\Lambda_{t}\}$ is the solution of $T$ coupled equations

where we denote $\boldsymbol{J}=(J_{mn})=\left((1-2R_{m})\Lambda_{m}^{-2}\delta_{mn}\right)$ and $\boldsymbol{K}=(K_{mn})=\left((R_{m}R_{n})(\Lambda_{m}\Lambda_{n})^{-1}\right).$ This expression makes it possible to evaluate $(\partial\Lambda_{m}/\partial\chi_{n})=(\partial\chi_{m}/\partial\Lambda_{n})^{-1}$ as

where $\boldsymbol{I}=(\delta_{mn})$, $\boldsymbol{T}=(\Lambda_{m}(1-2R_{m})^{-1/2}\delta_{mn})$ and $\boldsymbol{M}=(R_{m}R_{n}(1-2R_{m})^{-1/2}(1-2R_{n})^{-1/2})$. Inserting (61), (62), and (66) into (48) yields a set of equations to determine $\{\hat{\chi}_{t}\}$ for a given set of non-zero densities $\{\rho_{t}\}$

where we set $\boldsymbol{S}=(\delta_{mn}/(R_{m}\sqrt{1-2R_{m}}))$ and used the relation $\boldsymbol{T}/\hat{Q}_{t}=(1-R_{t})\boldsymbol{S}$, which is obtained from (47) and (63).

Equations (48) and (59) indicate that the critical condition for making $\chi_{t}=0$ stable is expressed as

For characterizing the critical condition for the success of the $l_{1}$-recovery, let us suppose that the set of non-zero densities $\{\rho_{t}\}$ is provided as a function of a single parameter $\mu$ as $\{\rho_{t}(\mu)\}$. For $T\geq 3$, the critical situation is specified by an appropriate set of $2T+1$ variables of $\mu$, $\{\hat{\chi}_{t}\}$, and $\{R_{t}\}$. These are provided by $2T+1$ conditions of (68)–(70).

On the other hand, the critical condition for $T=2$ is provided differently from (68)–(70) because the constraint of $\chi_{1}=\chi_{2}$ for keeping $F(\{\chi_{1},\chi_{2}\})$ well defined requires $R_{1}=R_{2}=1/2$ for any pair of $\rho_{1}$ and $\rho_{2}$. Explicitly, the condition is provided by the following four coupled equations:

where $\eta$ is a Lagrange parameter for enforcing $\chi_{1}=\chi_{2}$. These determine four variables of $\hat{\chi}_{1},\hat{\chi}_{2},\eta,\mu$ at the critical condition.

Equations (68)–(70) for $T\geq 3$ and (71)–(74) for $T=2$ constitute the main result of this paper. For $T=2$, an equivalent result can also be obtained in a slightly different manner .

5 Validity of RS evaluation

The above calculation can be generalized to arbitrary levels of replica symmetry breaking (RSB). For example, under one-step RSB (1RSB) ansatz, where $n$ replicas of each block $t$ are classified into $n/x$ groups of an identical size $x$ and their overlaps are assumed to be $\boldsymbol{x}_{t}^{a}\cdot\boldsymbol{x}_{t}^{b}/M=q_{1t}$ if $a$ and $b$ belong to the same group and $q_{0t}(\leq q_{1t})$ otherwise, (34) and (42) are modified as

respectively, where $\Phi=-(\hat{Q}_{t}+\hat{q}_{1t})x^{2}/2+(\sqrt{\hat{q}_{1t}-\hat{q}_{0t}}z_{1}+\sqrt{\hat{q}_{0t}}z_{0}+\hat{m}_{t}x_{t}^{0})x-\beta|x|$.

In the 1RSB framework, the RS solution is regarded as a special solution of $\Delta_{t}=q_{1t}-q_{0t}=0$ and $\hat{\Delta}_{t}=\hat{q}_{1t}-\hat{q}_{0t}=0$ ($t=1,2,\ldots,T$). In the limit of $n\to 0$ and $\beta\to\infty$ keeping $\chi=\beta(Q_{t}-q_{1t})\sim O(1)$, the critical condition that a solution of $\Delta_{t}>0$ and $\hat{\Delta}_{t}>0$ bifurcates from the RS solution, which corresponds to the de Almeida-Thouless (AT) condition of the current system, is expressed as

where $\boldsymbol{H}=(2(\partial^{2}/\partial\chi_{t}\partial\chi_{s})F(\{\chi_{k}\}))$ and $\hat{\boldsymbol{H}}=(\delta_{ts}\int Dz[(\partial X(\sqrt{\hat{\chi_{t}}}z+\hat{m}_{t}x_{t}^{0};\hat{Q}_{t})/\partial(\sqrt{\hat{\chi_{t}}}z))^{2}]_{x_{t}^{0}})$.

holds for the success solution at the critical condition of the $l_{1}$-recovery. Equation (80) yields an eigenvalue of unity whose eigenvector is given as $\boldsymbol{v}_{1}\propto((1-R_{t})^{-1})$, which makes (79) hold. Similarly, (79) is also satisfied at the critical condition of the $l_{1}$-recovery of $T=2$. These validate our RS evaluation in terms of the local stability analysis, although further justification with other schemes, such as comparison with numerical experiments, is necessary for examining possibilities that the RS solution becomes thermodynamically irrelevant due to discontinuous phase transitions.

Case studies

Let us examine the significance of the developed methodology by applying it to two representative examples.

We consider the uniform density case of $\rho_{t}=\rho$ ($t=1,2,\ldots,T$) as the first example, where the uniformity allows us to solve (47)–(52) setting all variables to be independent of $t$ as $\chi_{t}=\chi$. In particular, setting $\Lambda_{t}=\Lambda$ simplifies the expressions of (47)–(49) providing $R_{t}=1/T$ and $2(\partial^{2}/\partial\chi_{t}\partial\chi_{s})F(\{\chi_{k}\})=\chi^{-2}$. This makes it unnecessary to deal with the saddle point problems of $T=2$ in an exceptional manner. As a consequence, the critical condition of the $l_{1}$-recovery is expressed compactly by using a pair of equations as

for both $T\geq 3$ and $T=2$. By setting $\alpha=1/T$, these provide a critical condition identical to that obtained for the rotationally invariance matrix ensembles in earlier studies . This indicates that for vectors of the uniform non-zero density, the $l_{1}$-recovery performance of the concatenated orthogonal matrices is identical to that of the standard setup provided by the matrix of independent standard Gaussian entries.

However, this is not the case when the non-zero density is not uniform. As a distinctive example, we examined the case of localized density, which is characterized by setting $\rho_{1}=T\rho$ and $\rho_{t}=0$ for $t=2,3,\ldots,T$. Such an assumption is plausible in handling various kinds of real world signals at least as a first approximation; due to the intrinsic nature of real world signals, one can expect that the density of the sparsest expression is localized in a certain block when the concatenation of identity and randomly chosen Fourier/wavelet bases, which is a representative example of the $T$-concatenation of orthonormal bases, is employed. Table 1 and Figure 1 show critical values of the total non-zero density $\rho=T^{-1}\sum_{t=1}^{T}\rho_{t}$ given the compression rate $\alpha=1/T$ for the uniform and localized density cases. These show that the concatenated matrices always result in better $l_{1}$-recovery performance for vectors of the localized densities, and the significance increases as $T$ becomes smaller while matrices of rotationally invariant ensembles result in identical performance as long as $\rho$ is unchanged. It is noteworthy that the performance gain is obtained without utilizing the knowledge of the profile of $\rho_{t}$ in the recovery stage. This indicates that, in addition to their ease of implementation and theoretical elegance, the concatenated orthogonal matrices are preferred for practical use in terms of their high recovery performance for vectors of non-uniform non-zero densities.

2 Numerical justification

Summary

Promising future studies include performance evaluation in the case of noisy situations and development of approximate recovery algorithms suitable for the concatenated orthogonal matrices .

Appendix A Derivation of (17)

We insert the Fourier expressions of $\delta$-function

Evaluating this by means of the saddle point method with respect to $\Lambda_{t}$ $(t=1,2,\ldots,T)$ results in

Similarly, the denominator of (84) is evaluated as

Substituting (84), (91), and (92) into (15) leads to the expression of (17).

References

References