EURASIP Journal on Advances in Signal Processing manuscript No.
(will be inserted by the editor)

Computable performance guarantees for compressed sensing
matrices

arXiv:1604.02769v2 [] 4 Feb 2018

Myung Cho · Kumar Vijay Mishra · Weiyu Xu

Received: date / Accepted: date

Abstract The null space condition for ℓ1 minimization in compressed sensing is a necessary and sufficient
condition on the sensing matrices under which a sparse
signal can be uniquely recovered from the observation
data via ℓ1 minimization. However, verifying the null
space condition is known to be computationally challenging. Most of the existing methods can provide only
upper and lower bounds on the proportion parameter that characterizes the null space condition. In this
paper, we propose new polynomial-time algorithms to
establish upper bounds of the proportion parameter.
We leverage on these techniques to find upper bounds
and further develop a new procedure - tree search algorithm - that is able to precisely and quickly verify
the null space condition. Numerical experiments show
that the execution speed and accuracy of the results
obtained from our methods far exceed those of the previous methods which rely on Linear Programming (LP)
relaxation and Semidefinite Programming (SDP).
Keywords compressed sensing · null space condition ·
ℓ1 minimization · performance guarantee · sensing
matrix
1 Introduction
Compressed sensing is an efficient signal processing technique to recover a sparse signal from fewer samples
than required by the Nyquist-Shannon theorem, reducing time and energy spent in sampling operation. These
advantages make compressed sensing attractive in various signal processing areas [1].
Myung Cho, Kumar Vijay Mishra, and Weiyu Xu
Dept. of ECE, University of Iowa, Iowa City, IA, 52242
E-mail:
(myung-cho,
kumarvijay-mishra,
weiyuxu)@uiowa.edu.

In compressed sensing, we are interested in recovering the sparsest vector x ∈ Rn that satisfies the underdetermined equation y = Ax. Here R is the set of real
numbers, A ∈ Rm×n , m < n, is a sensing matrix, and
y ∈ Rm is the observation or measurement data. This
is posed as an ℓ0 minimization problem:
minimize kxk0
subject to y = Ax,

(1)

where kxk0 is the number of non-zero elements in vector
x. The ℓ0 minimization is an NP-hard problem. Therefore, we often relax (1) to its closest convex approximation - the ℓ1 minimization problem:
minimize kxk1
subject to y = Ax.

(2)

It has been shown that the optimal solution of ℓ0 minimization can be obtained by solving ℓ1 minimization
under certain conditions (e.g. Restricted Isometry Property or RIP) [2,3,4,5,6]. For random sensing matrices,
these conditions hold with high probability. We note
that RIP is a sufficient condition for sparse recovery
[7].
A necessary and sufficient condition under which a
k-sparse signal x, (k ≪ n) can be uniquely obtained via
ℓ1 minimization is Null Space Condition (NSC) [8,9,3].
A matrix A satisfies NSC for a positive integer k if
||zK ||1 < ||zK ||1

(3)

holds true for all z ∈ {z : Az = 0, z 6= 0} and for
all subsets K ⊆ {1, 2, ..., n} with |K| ≤ k. Here K is
an index set, |K| is the cardinality of K, zK is the
part of the vector z over the index set K, and K is
the complement of K. NSC is related to the proportion
parameter αk defined as

2

Myung Cho et al.
αk ,

maximize

maximize

{z: Az=0, z6=0}{K: |K|≤k}

kzK k1
.
kzk1

(4)

The αk is the optimal value of the following optimization problem:
maximize

z,{K: |K|≤k}

kzK k1

subject to kzk1 ≤ 1, Az = 0,

(5)

where K is a subset of {1, 2, . . . , n} with cardinality
at most k. The matrix A satisfies NSC for a positive
integer k if and only if αk < 21 . Equivalently, NSC can
be verified by computing or estimating αk . The role of
αk is also important in recovery of an approximately
sparse signal x via ℓ1 minimization where a smaller αk
implies more robustness [8,9,10].
We are interested in computing αk and, especially,
finding the maximum k for which αk < 21 . However,
computing αk to verify NSC is extremely expensive
and was reported in [7] to be NP-hard. Due to the
challenges in computing αk , verifying NSC explicitly
for deterministic sensing matrices remains a relatively
unexamined research area. In [11,12,3,8], convex relaxations were used to establish upper or lower bounds of
αk (or other parameters related to αk ) instead of computing the exact αk . While [3,11] proposed semidefinite programming based methods, [8,12] suggested linear programming relaxations to obtain the upper and
lower bounds of αk . For both methods, computable performance guarantees on sparse signal recovery were reported via bounding αk . However, these bounds of αk
√
could only verify NSC with k = O( n), even though
theoretically k can grow linearly with n.
Our work drastically departs from these prior methods [11,12,3,8] that provide only the upper and lower
bounds. In our solution, we propose the pick-l-element
algorithms (1 ≤ l < k), which compute upper bounds
of αk in polynomial time. Subsequently, we leverage on
these algorithms to develop the Tree Search Algorithm
(TSA) - a new method to compute an exact αk by significantly reducing computational complexity of an exhaustive search method. This algorithm offers a way to
control a smooth tradeoff between complexity and accuracy of the computations. In the conference precursor
to this paper, we had introduced Sandwiching Algorithm (SWA) [13], which employs a branch-and-bound
method. Although SWA can also be used to calculate
the exact αk , it has a disadvantage of greater memory usage than TSA. On the other hand, TSA provides
memory and performance benefits for high-dimensional
matrices (e.g., up to size ∼ 6000 × 6000).
It is noteworthy that our methods are different from
RIP or the neighborly polytope framework for analyzing the sparse recovery capability of random sensing matrices. For example, prior works such as [6,22]

employ the neighborly polytope to predict theoretical
lower bounds on recoverable sparsity k for a randomly
chosen Gaussian matrix. However, our methods do not
resort to a probabilistic analysis and are applicable for
any given deterministic sensing matrix. Also, our algorithms have the strength of providing better bounds
than existing methods [11,12,3,8] for a wide range of
matrix sizes.

1.1 Main contributions
We summarize our main contributions as follows:
(i) Faster algorithms for high dimensions. We
designed the pick-l algorithm (and its optimized
version), where l is a chosen integer, to provide
upper bounds on αk . We are able to show that
when l increases, the optimized pick-l algorithm
provides tighter upper bound on αk . Numerical
experiments show that, even with l = 2 or 3, the
pick-l algorithm already provides better bound
on αk than the previous algorithms based on the
LP [8] and SDP [3]. For large sensing matrices,
the pick-1-element algorithm can be significantly
faster than the LP and SDP methods.
(ii) Novel formulations using branch-and-bound.
Based on the pick-l algorithm, we propose a branchand-bound tree search approach to compute tighter
bounds or even the exact value of αk . To the
best of our knowledge, this tree search algorithm
is the first branch-and-bound algorithm to verify NSC for ℓ1 minimization. This branch-andbound approach heavily depends on the pick-l algorithm developed in this paper. For example, the
LP [8] and SDP [3] methods cannot be directly
adapted to provide an efficient branch and bound
approach, due to their lack of subset-specific upper bounds on αk . In numerical experiments, we
demonstrated that the tree search algorithm reduced the execution time to precisely calculate αk
by around 40-8000 times, compared to the exhaustive search method.
(iii) Simultaneous upper and lower bounds. The
branch-and-bound tree search algorithm simultaneously maintains upper and lower bounds of αk
during the run-time. This approach has two benefits. Firstly, if one is interested in merely certifying
the NSC for a positive k rather than obtaining the
exact αk , then one can terminate the TSA early
to shorten the running time. This can be done as
soon as the global upper (lower) bound drops below (exceeds) 1/2 and, therefore, concluding that
the NSC for the positive k is satisfied (not sat-

Computable performance guarantees for compressed sensing matrices

isfied). Secondly, consider the case when TSA is
terminated early due to, say, constraints on running time. Then, the process still yields meaningful bounds on αk via the record of continuously
maintained upper and lower bounds.
(iv) New results on recoverable sparsity. For a
certain l < k, we can compute αl or its upper
bound by using the branch-and-bound tree search
algorithm (for example, based on the pick-1-element
algorithm). We introduce a novel result (Lemma
3), which can use αl to lower bound the recoverable sparsity k. This approach of lower bounding the recoverable sparsity k is useful when l is
too large to perform
 the pick-l algorithm directly
(which requires nl enumerations).

3

puters. Our goal is to devise algorithms that can rapidly
yield an exact value of αk . As an initial step, we develop
a method to compute an upper bound of αk in polynomial time, which is called the pick-l-element algorithm
(or simply, pick-l algorithm), where l is a chosen integer
such that 1 ≤ l < k.
Let us define the proportion parameter for a given
index set L such that |L| = l, denoted by αl,L , as
αl,L ,

maximize

{z: Az=0, z6=0}

kzL k1
.
kzk1

(6)

(6) is the partial optimization problem of (4) only considering the vector z in the null space of A for a fixed
index set L. We can obtain αl,L by solving the following
optimization problem:
maximize kzL k1
z

1.2 Notations and preliminaries
We denote the sets of real numbers, and positive integers as R and Z+ respectively. We reserve uppercase letters K and L for index sets, and lowercase letters k, l ∈ Z+ for their respective cardinalities. We
also use | · | to denote the cardinality of a set. We assume k > l ≥ 1 throughout the paper. For vectors or
scalars, we use lowercase letters, e.g., x, k, l. For a vector x ∈ Rn , we use xi for its i-th element. If we use an
index set as a subscript of a vector, it represents the
partial vector over the index set. For example, when
x ∈ Rn and K = {1, 2}, xK represents [x1 , x2 ]T . We reserve uppercase A for a sensing matrix whose dimension
is m × n. Since the number of columns of a sensing matrix A is n, the full index set we consider is {1, 2, ..., n}.
In addition,we represent nl numbers of subsets as Li ,
i = 1, ..., nl , where Li ⊂ {1, 2, ..., n}, |Li | = l. We use
the superscript * to represent an optimal solution of
an optimization problem. For instance, z ∗ and K ∗ are
the optimal solution of (5). Since we need to represent
an optimal solution for each index set Li , we use the
superscript i∗ to represent an optimal solution for an
index set Li , e.g., z i∗ . The maximum value of k such
that both αk < 21 and αk+1 ≥ 21 hold true is denoted
by the maximum recoverable sparsity kmax .
2 Pick-l-element Algorithm
Consider
 a sensing matrix with n columns. Then, there
are nk subsets K each of cardinality k. When n and
k are large, exhaustive search over these subsets to
compute αk is extremely expensive. For example, when
n = 100 and k = 10, it takes a search over 1.7310e+13
subsets to compute αk - a combinatorial task that is beyond the technological reach of common desktop com-

subject to kzk1 ≤ 1, Az = 0.

(7)

Since (7) is maximizing a convex function for a given
subset L, we cast (7) as 2l linear programming problems by considering all the possible sign patterns of
every element of zL (e.g., if l = 2 and L = {1, 2}, then,
||zL ||1 = |z1 | + |z2 | can correspond to 2l = 4 possibilities: z1 + z2 , z1 − z2 , −z1 + z2 , and −z1 − z2 ). αl,L is
equal to the maximum among the 2l objective values.
The pick-l algorithm uses αl,L ’s obtained from different index sets to compute an upper bound of αk .
Algorithm 1 shows the steps of the pick-l algorithm
in detail. The following Lemmata show that the pick-l
algorithm provides an upper bound of αk . Firstly, we
provide Lemma 1 to derive the upper bound of the proportion parameter for a fixed index set K, and then, we
show that the pick-l algorithm yields an upper bound
of αk in Lemma 2.
Lemma 1 (Cheap Upper Bound (CUB) for a given
subset K) Given a subset K, we have
CU B(αk,K ) ,

1
k−1
l−1



X

αl,Li ≥ αk,K .

(8)

{Li ⊆K, |Li |=l}

Proof Suppose that when z = z i∗ and z = z ∗ , we
achieve the optimal value of (6) for given index sets
i∗
kzL
k1
i
kz i∗ k1 and αk,K

appears k−1
l−1 times

Li and K respectively, i.e., αl,Li =
∗
kzK
k1
kz ∗ k1 .

=

Since each element of K
in
{Li ⊆ K, |Li | = l}, we obtain the following inequality:
αk,K =

≤

∗
kzK
k1
=
kz ∗ k1

1

1
k−1
l−1



X

X

∗
kzL
k1
i

{Li ⊆K, |Li |=l}

k−1
l−1 {Li ⊆K, |Li |=l}

i∗
kzL
k1
i

kz i∗ k1

kz ∗ k1

= CU B(αk,K ).

The inequality is from the optimal value of (6) for each
index set Li .
⊓
⊔

4

Myung Cho et al.

Algorithm 1: Pick-l-element algorithm, 1 ≤ l < k
for computing an upper bound of αk
1: Given a matrix A, calculate αl,L ’s for all the subsets L,
|L| = l, via (7).

2: Sort these nl different values of αl,L ’s in descending
order like (10).
3: Compute an upper bound of αk via (9).
4: If the upper bound of αk is larger than 1, then, set the
upper bound to 1. If the upper bound is less than 21 ,
then NSC for k ∈ Z+ is satisfied.

Lemma 2 The pick-l algorithm provides an upper bound
of αk , namely
αk ≤

(kl)
X

1

αl,Li ,

k−1
l−1 i=1

(9)

where αl,L1 ≥ αl,L2 ≥ · · · ≥ αl,Li ≥ · · · ≥ αl,L n . (10)
(l)

Proof Without loss ofgenerality, we assume that when
z = z i∗ , i = 1, 2, ..., nl , αl,Li ’s are obtained in descending order like (10). It is noteworthy that αk,K is defined
for a fixed K set; however, αk is the maximum value
over all the subsets with cardinality k. Suppose that
when z = z ∗ and K = K ∗ , αk is achieved in (4). From
the aforementioned definitions and similar argument as
in Lemma 1, we have:
αk = α

k,K ∗

1

≤

k−1
l−1



X

αl,Li ≤

{Li ⊆K ∗ , |Li |=l}

1
k−1
l−1



(kl)
X

αl,Li .

i=1

The steps 2 and 3 in Algorithm 1, which are sorting
αl,L ’s and computing an upper bound of αk with sorted
αl,L ’s via (9), can also be done by solving the following
optimization problem without sorting operation:
maximize
n
)
l

γi , 1≤i≤(

γi αl,Li

i=1

1

subject to 0 ≤ γi ≤
(nl)
X

γi ≤

i=1

k−1
l−1

k
.
l

, 1 ≤ i ≤(nl),

(11)


× k = kl . Therefore, for the
( )  l
optimal value, the first kl largest αl,Li ’s are chosen
1
.
with the coefficient k−1
( l−1 )
The upshot of the pick-l algorithm
 is that we can
reduce number of operations from nk enumerations to
Here, we note that

1

k−1
l−1


. For example, when n = 300, k = 20, and l = 2,
the number of operations is reduced by around 1026
times. Moreover, as n increases, the reduction rate increases. With the reduced enumerations, we can still
have non-trivial upper bounds of αk through the pickl-element algorithm. We will present the performance
of the pick-l algorithm in Section 5 showing that the
pick-l algorithm provides better upper bounds than the
previous research [3,8] even when l = 2. Furthermore,
thanks to the pick-l algorithm, we can design a new
algorithm based on a branch-and-bound search to calculate αk by using upper bounds of αk obtained from
the pick-l algorithm. It is noteworthy that the cheap
upper bound introduced in Lemma 1 can provide upper bounds on αk,K for specific subsets K, which enable
our branch-and-bound method to calculate αk or more
precise bounds on αk . However, LP relaxation method
[8] and SDP method [3] do not provide upper bounds
on αk,K for specific subsets K, which overwhelms LP
and SDP methods to be used in the branch-and-bound
method.
Since we are also interested in kmax , we introduce
the following Lemma 3 to bound the maximum recoverable sparsity kmax .
Lemma 3 The maximum recoverable sparsity kmax satisfies


1/2
− 1 ≤ kmax ,
k(αl ) , l ·
αl

The first inequality is from Lemma 1, and the last inequality is from the assumption that αl,Li ’s are sorted
in descending order.
⊓
⊔

(nl)
X

n
l

(12)

where ⌈.⌉ is the ceiling function.
Proof To prove
this lemma, we will show that when

−
1,
αk < 21 . This can be concluded from
k = l · 1/2
αl
the upper bound of αk given as follows:
αk = αk,K ∗ ≤

≤

1

X

k−1
l−1 {Li ⊆K ∗ , |Li |=l}
k
l
 αl
k−1
l−1



= αl ·

αl,Li

k
.
l

(13)


Note that there are kl terms in the summation. From
(13), if αl · kl < 12 , then αk < 12 . In other words, if
1
k < l · 1/2
αl , then αk < 2 . Since k is a positive integer,
 1/2 
when k = l · αl −1, αk < 12 . Therefore, the maximum
recoverable sparsity
kmax should be larger than or at

least equal to l · 1/2
− 1.
⊓
⊔
αl

It is noteworthy that in [8, Section 4.2.B], the authors
introduced lower bound on k based on α1 , i.e., k(α1 ).
However, in Lemma 3, we provide a more general result.
Furthermore, in Lemma 3, instead of using αl , we can
use an upper bound of αl to obtain

 the recoverable
sparsity k; namely, k(U B(αl )) =

l·

1/2
UB(αl )

−1 ≤

Computable performance guarantees for compressed sensing matrices

kmax , where U B(αl ) represents an upper bound of αl .
Since the proof follows the same track as the proof of
Lemma 3, we omit the proof.
Finally, we introduce the following proposition to
compare our algorithm to LP method [8] theoretically.
Proposition 1 For any integer k ≥ 1, let αpick1
be
k
the upper bound on αk provided by the pick-1-element
algorithm according to Lemma 2. Let αLP
be the upper
k
bound on αk provided by the LP method [8] according
to the following definition (namely Equation (4.25) in
[8] with β = ∞) :
=
αLP
k

minimize

Y =[y1 ,...,yn ]∈Rm×n



maximize ||(I − Y T A)ej ||k,1 ,
1≤j≤n

where ej is the standard basis vector with the j-th element equal to 1, and k · kk,1 stands for the sum of k
maximal magnitudes of components of a vector. Then
we have:
αpick1
≥ αLP
k .
k

(14)

For readability, we place the proof of Theorem 1 in Appendix A.
The LP method can provide tighter upper bounds
on αk than the pick-1-element algorithm, however this
comes at a cost of solving a big optimization problem
of design dimension mn. When m and n are large, the
complexity of computing αLP
k can be prohibitive (please
see Table 2).

5

In the following lemmata, we show that the optimized pick-l algorithm produces an upper bound of αk
and this bound is tighter than that of the basic pick-l algorithm introduced in (11). The last lemma establishes
that as l increases, the upper bound of αk decreases or
stays the same.
Lemma 4 The optimized pick-l algorithm provides an
upper bound of αk .
Proof The strategy to prove Lemma 4 is to show that
one feasible solution of (15) gives an upper bound of
αk . Suppose when K = K ∗ , αk is achieved, i.e., αk =
1
αk,K ∗ . For a feasible solution, let us choose γi = k−1
( l−1 )
when Li ⊆ K ∗ , and γi = 0 otherwise, which we can
easily check whether it satisfies the first and second
constraints of (15). For the third constraint, let us check
the case when b = l first. For b = l, we can choose an
arbitrary index set B such that |B| = b = l. For the
chosen B, there is only one Li such that B ⊆ Li , which
is itself, i.e., B = Li . For other chosen B’s, it is the
same. Hence, the third constraint represents
γi ≤

1
k−1
l−1

, i = 1, 2, ...,

We can tighten the upper bound of αk in the pick1
in
l algorithm by replacing the constant factor k−1
( l−1 )
(9) with optimized coefficients at the cost of additional
complexity, which we call as the optimized pick-l algorithm. This optimized pick-l algorithm is mostly useful
from a theoretical perspective. In practice, it gives improved but similar performance in calculating the upper
bound of αk to the basic pick-l algorithm described in
Section 2. As a theoretical merit of the optimized pick-l
algorithm, we can show that as l increases, the upper
bound of αk becomes smaller or stays the same.
The optimized pick-l algorithm provides an upper
bound of αk via the following optimization problem:
maximize
n
)
l

γi , 1≤i≤(

(nl)
X

X

γi ≤ 1.

(17)

n
l

), |B|=1}


Note that there are n−1
’s which have an
l−1 numbers of L
i
n−1
index set B as a subset. Among l−1 numbers of γi ’s,
only γi ’s whose corresponding Li ’s are the subsets of K ∗
1
. Since each element in Li such that Li ⊆ K ∗
are k−1
( l−1 )


n
∗
appears k−1
l−1 times in {i : Li ⊆ K , 1 ≤ i ≤ l },
the summation of γi , where the corresponding
Li ’s are

1
× k−1
=
1,
which
the subset of K ∗ , becomes k−1
l−1
( l−1 )
satisfies (17). Basically, the third constraint makes that
for an index, the summation of coefficients related to
the index is limited to 1. In the same way, for 1 < b < l,
the chosen γi is a feasible solution
P of (15). From this
1
α
feasible solution, we have k−1
( l−1 ) {i: Li ⊆K ∗ , |Li |=l} l,Li
for the optimal value, which is an upper bound of αk
as shown in (13).
⊓
⊔
Lemma 5 The optimized pick-l algorithm provides a
tighter, or at least the same, upper bound of αk than
the basic pick-l algorithm introduced in (11).

γi αl,Li

i=1

subject to γi ≥ 0, 1 ≤ i ≤(nl),
(nl)
X

(16)

For b = 1, the third constraint represents
{i: B⊆Li , 1≤i≤(

3 Optimized Pick-l Algorithm

n
.
l

Proof We will show that the optimization problem (11)
is a relaxation of (15). As in the proof of Lemma 4,
i=1

for b = l, the third constraint of (15) represents (16),
k−b
X
l−b
∀ b ∈ Z+ s.t. 1 ≤ b ≤ l,
γi ≤ k−1 , ∀ B with |B| = b
. which is involved in the first constraint of (11). Since
n
l−1
{i: B⊆Li , 1≤i≤( )}
the third constraint of (15) considers other b values such
l
γi ≤

k
,
l

(15)

6

Myung Cho et al.

that 1 ≤ b < l, (15) has more constraints than (11).
Therefore, the optimized pick-l algorithm, which is (15),
provides a tighter or at least the same upper bound than
the basic pick-l algorithm.
⊓
⊔
Lemma 6 The optimized pick-l algorithm provides a
tighter or at least the same upper bound than the optimized pick-p algorithm when l > p.
Proof We can upper bound the objective function of
(15) by using (8) as follows:
1

maximize
n
γi , 1≤i≤( )
l


l−1

p−1

(nl)
X

X

γi

i=1

γi ≤

i=1

X

=

≤
(18)
=

γi ≤

k−b
l−b
,
k−1
l−1



+

∀ b ∈ Z s.t. 1 ≤ b ≤ l,
.
∀ B with |B| = b

γj
n
p

{j: B⊆Pj , 1≤j≤(

{j: Pj ⊂Li , |Pj |=p}

k
,
l

′

X

=

αp,Pj

subject to γi ≥ 0, 1 ≤ i ≤(nl),
(nl)
X

where the second equality is obtained from the fact that
γi , which is a coefficient of αl,Li , appears pl times in
P(np) P
j=1
{i: Pj ⊂Li } γi . The final inequality is from the
second constraint of (18). The third constraint in (19)
can be deduced from the following inequality:

1
l−1
p−1



l−1 
p−1

1
l−1
p−1

k−b
p−b
,
k−1
p−1



n−b n−p
p−b
l−p
n−b
l−b



X

n
p

{j: B⊆Pj , 1≤j≤(

1



)}
X



n−b n−p
l−p
p−b

n−b
l−b

γi

)} {i: Pj ⊂Li , 1≤i≤(nl)}
X

γi
n
l

{i: B⊂Li , 1≤i≤(
k−b
l−b
,
k−1
l−1



)}

1≤b≤p

1 ≤ b ≤ p,

where the second equality
 is from the fact that for a
fixed Pj , there are n−p
Pj ⊂
l−p numbers of Li ’s, where


n
n−b
,
1
≤
Note that
in
the
objective
function
of
(18),
each
α
Li , i = 1, ..., l ; for a fixed B, there are p−b numbers
p,Pj




j ≤ np , appears n−p
l−p times. Let us define
of Pj ’s, where B ⊂ Pj , j = 1, ..., np , and n−b
numl−b

n
bers of Li ’s, where B ⊂ Li , i = 1, ..., l . Since (19) is
X
′
1
γi .
γj , l−1 
obtained from the relaxation of (18), the optimal value
p−1 {i: Pj ⊂Li , 1≤i≤(n)}
l
of (19) is larger or equal to the optimal value of (18).
(19) is just the optimized pick-p algorithm. Thus, when
We can relax (18) to the following problem, which turns
l > p, the optimized pick-l algorithm provides a tighter
out to be the same as the optimized pick-p algorithm:
or at least the same upper bound than the optimized
(np)
pick-p algorithm.
⊓
⊔
X
n
l

{i: B⊆Li , 1≤i≤(

)}

′

γj αp,Pj

maximize
′
n
γj , 1≤j≤( )
p

By using larger l in the pick-l algorithm, we can obtain
a tighter upper bound of
subject to γj ≥ 0, 1 ≤ j ≤( ),
 αk . However, for a certain
n
l,
we
need
to
enumerate
l possibilities, and this be(np)
X
comes infeasible when l is large. Moreover, when l < k,
k
′
(19)
γj ≤ ,
the pick-l algorithm only gives an upper bound of αk ,
p
j=1

instead of an exact value of αk . There is, however, a
k−b
X
′
p−b
b ∈ Z+ s.t. 1 ≤ b ≤ p, need to find tighter bounds on α , or to even find the
γj ≤ k−1 , ∀
.
k
∀ B with |B| = b

p−1
exact
value
of
α
,
when
k
is
too
large for nk enumer{j: B⊆Pj , 1≤j≤(n)}
k
p
ations of exhaustive search [14,15,16]. To this end, we
The relaxation is shown by checking the constraints.
propose a new branch-and-bound tree search algorithm
The first constraint of (19) is trivial to obtain. For the
to find tighter bounds on αk than Lemma 2 provides, or
second constraint, we can obtain the second constraint
to even find the exact αk . Our branch-and-bound tree
of (19) from the following relations:
search algorithm is enabled by the pick-l algorithms introduced in Sections 2 and 3.
n
n
j=1

′

( p)
X

′

γj =

(p )
X

j=1

j=1

n
p

1
l−1
p−1

1



l−1 
p−1

γi
n
l

{i: Pj ⊂Li , 1≤i≤(

(nl)
l X

γi
p i=1
1 lk
k
≤ l−1 
= ,
p
l
p
p−1

=

X

)}

4 Tree Search Algorithm
To find the index set K ∗ which leads to the maximum αk,K (among all possible index set K’s), the Tree
Search Algorithm (TSA) performs a best-first branchand-bound search [23] over a tree structure representing

Computable performance guarantees for compressed sensing matrices

different subsets of {1, 2, ..., n}. In its essence, for each
subset J with cardinality no bigger than k, TSA calculates an upper bound of αk,K , which is valid for any
set K (with cardinality k) such that J ⊆ K. If this upper bound is smaller than a lower bound of αk , TSA
will not further explore any of J’s supersets, leading
to reduced average-case computational complexity. For
simplicity, we will describe the TSA based on pick-1element algorithm, simply called 1-Step TSA. However,
we remark we can also extend the TSA to be based
on pick-l-element (l ≥ 2) algorithm, by calculating upper bounds of αk,K based on the results of the pick-lelement algorithm.
4.1 Tree structure
A tree node J represents an index subset of {1, ..., n}
such that |J| ≤ k. We have the following rule:
[R1] A parent node is a subset of each of its child node(s).
A node that has no child is referred to as a leaf node.
We call the cardinality of the index set corresponding
to J as J’s height. The tree structure follows the “legitimate order”, which ensures that any new index in the
child node is bigger than the indices of its parent node.
[R2] “Legitimate order” - Let P and C denote the parent
node, and the child node. Then, any index in P must
be smaller than any index in C \ P .
Fig. 1 illustrates this rule in a tree with k = 2 and
n = 3.

7

an upper bound B(J), which is no smaller than αk,K
for any child node K (with cardinality k) of node J. If
B(J) is smaller than the lower bound on αk , then the
algorithm will not explore the child nodes of the tree
node J.
In our algorithm, we calculate B(J) as
B(J ) = αj,J +

t
X

α1,{i+max(J)} ,

(20)

i=1

where j + t = k, max(J) represents the largest index
in J, and α1,{1} ≥ α1,{2} ≥ ... ≥ α1,{n} . We obtain
this descending order by permuting the columns of the
sensing matrix A in descending order of α1,{i} ’s as the
pre-computation step of TSA. For example, in Fig. 1,
for k = 2, B({1}) = α1,{1} + α1,{2} . In order to justify
that B(J) is an upper bound of αk,K for all node K
such that J ⊆ K, we provide the following lemma.
Lemma 7 Given α1,{1} ≥ α1,{2} ≥ ... ≥ α1,{n} , B(J) =
P
αj,J + ti=1 α1,{i+max(J)} , where j +t = k, and max(J)
represents the largest index in J, is an upper bound of
αk,K for all nodes K such that J ⊆ K.
Proof For any subset K such that J ⊆ K, we can write
αk,K = αj+t,{J∪T } , where j + t = k and T = K \ J.
Then, following exactly the same line of argument as in
the proof of Lemma 1, we have
αk,K ≤ αj,J + αt,T ,
P
and αt,T is no larger than tj∈T α1,{j} . Finally, since
P
α1,{i} ’s are sorted in the descending order, j∈T α1,{j} ≤
Pt
i=1 α1,{i+max(J)} . Note that, due to the legitimate order [R2], the smallest element of the index set T is no
less than 1 + max(J). In conclusion, for all nodes K
such that J ⊆ K, B(J) becomes an upper bound of
αk,K .
⊓
⊔
4.3 Best-first tree search strategy

Fig. 1: A tree structure following the legitimate order for
k = 2 and n = 3.

4.2 Basic idea of a branch-and-bound approach for
calculating αk
We use a branch-and-bound approach over the tree
structure to calculate αk . This method maintains a
lower bound on αk (how to maintain this lower bound
will be explained in Subsection 4.3). When the algorithm explores a tree node J, the algorithm calculates

TSA adopts a best-first tree search strategy for the
branch-and-bound approach. We first describe a basic
version of the best-first tree search strategy, and then
introduce two enhancements to this strategy in the next
subsection.
In its basic version, TSA starts with a tree having
only the root node, and sets the global lower bound
of αk as 0. In each iteration, TSA selects a leaf tree
node J with the largest B(J), and expands the tree by
adding the child nodes of J to the tree. For each of these
newly added child nodes, say Q, TSA then calculates
the upper bound B(Q) in (20). Note that if a newly
added child node Q has k elements, TSA will calculate

8

Myung Cho et al.

αk,Q , which is a lower bound on αk . For this k-element
Q, if the newly calculated αk,Q is bigger than the global
lower bound of αk , TSA will set the global lower bound
equal to αk,Q . TSA will terminate if a leaf tree node J
has the largest B(J) among all the leaf nodes, and that
B(J) is no bigger than the global lower bound on αk .
From standard theories of the branch-and-bound
approach, this TSA will output the exact αk . Also, in
this process, the global lower bound will keep increasing
until it is equal to an upper bound of αk (the largest
B(J) among leaf nodes).

Algorithm 2: Tree search algorithm based on the
pick-1-element algorithm (1-Step TSA)
Input: A ∈ Rm×n , k, l ← 1 ⊲ 1-Step TSA, i.e., l = 1
Output: αk
1
2

3

4
5

6
7

4.4 Two enhancements

8
9
10

We incorporate two novel features to TSA in order
to reduce the computational complexity. Firstly, when
TSA attaches a new node Q to a node J in the tree
structure, TSA computes B(Q) as (21):

11
12
13
14
15
16

B(Q) = αj,J + α1,Q\J +

t
X

17

α1,{i+max(Q)} ,

(21)

18

⊲ Pre-computation:
compute αl,{i} for i = 1, ..., n via (7)
permute columns of A in descending order of α1,{i} ’s,
so that α1,{1} ≥ ... ≥ α1,{n}
⊲ Tree expansion:
P
start with root node ∅, where B(∅) = k
i=1 α1,{i} , in
a tree structure Υ
Loop
J ← a node that has the largest B(·) among all
the leaf nodes in Υ
j ← |J |
if αj,J is not calculated then
compute αj,J via (7) and update B(J ) via (20)
expand Υ from the parent of J
⊲ See [R3]
else
if j = k then
αk ← B(J )
break
else
expand Υ from J
⊲ See [R3]
end
end
EndLoop

i=1

where j + t + 1 = k, max(Q) represents the largest
index in Q, and α1,{1} ≥ α1,{2} ≥ ... ≥ α1,{n} . Thus,
without calculating αj+1,Q (which involves higher computational complexity), we can still have B(Q) as an
upper bound of αk,K for any child node K (with cardinality k) of the node Q.
Secondly, when TSA adds a new node Q as the child
of node J in the tree structure (assuming αj,J has already been calculated), TSA does not need to add all
of J’s child nodes to the tree at the same time. Instead,
TSA only adds the node J’s unattached child node Q
with the largest B(Q) as defined in (21). Namely, the
index Q \ J is no bigger than the index Q′ \ J, where
Q′ is any unattached child of the node J. We note that
B(Q) is an upper bound on B(Q′ ) (according to (21))
for any other unattached child node Q′ of the node J.
Thus, for any child node K (of cardinality k) of node J’s
unattached child nodes, B(Q) is still an upper bound
of αk,K .
Algorithm 2 shows detailed steps of TSA, based on
the pick-1-element algorithm (namely, l = 1, 1-Step
TSA). In the description, we define “expanding the tree
from a node J” as follows:
[R3] “Expanding the tree from a node J” - Attaching
a new node Q to the node J, where B(Q) is the
largest value defined as (21) among the node J’s all
the unattached child nodes.

4.5 Advantage of the tree search algorithm
Due to the nature of the branch-and-bound approach,
we can obtain a global upper bound and a global lower
bound of αk while TSA runs. As the number of iterations increases in TSA, we can obtain tighter and
tighter upper bounds on αk , which is the largest B(·)
among leaf nodes. By using the global upper bound
of αk , we can obtain a lower bound of the recoverable
sparsity k via Lemma 3. Thus, even if the complexity
of TSA is too high to finish in a timely manner, we can
still obtain a lower bound on the recoverable sparsity k
by early terminating TSA.
We note that the methods based on LP [8] and SDP
[3] also provide upper bounds on αk . However, they are
unable to determine upper bounds of αk,K , which is
for a specific index set K. This prevents the use of LP
and SDP methods in our branch-and-bound method for
computing αk .

5 Numerical Experiments
We conducted extensive simulations to compute αk and
its upper/lower bounds using the pick-l algorithms and
TSA. In this section, we call the pick-l algorithms introduced in Section 2 and 3 as simply the (basic) pick-l
and the optimized pick-l algorithms respectively.

Computable performance guarantees for compressed sensing matrices

For same matrices, we compared our methods with
LP relaxation [8] approach and SDP method [3]. We
assessed the computational complexity in terms of execution time of the algorithms.1 In addition, we carried
out numerical experiments to demonstrate the computational complexity of TSA empirically.
For LP method in [8] and SDP method in [3], we
used the Matlab codes2 provided by the authors. Consistent with previous research, we used CVX [17] - a
package for specifying and solving convex programs for the SDP method, and MOSEK [18] - a commercial
LP solver - for the LP method. In our own algorithms,
we used MOSEK to solve (7). Also, to be consistent
with the previous research, matrices were generated
from the Matlab code provided by the authors of [3] at
http://www.di.ens.fr/~aspremon/NSPcode.html. For
valid bounds, we rounded down lower bounds on αk and
exact αk , and rounded up upper bounds on αk to the
nearest hundredth.
5.1 Performance comparison
Firstly, we considered Gaussian matrices and partial
Fourier matrices sized from n = 40 to n = 6144. We
chose n = 40 so that our results can be compared with
the simulation results in [3].
5.1.1 Low-dimensional sensing matrices
Sensing matrices with n = 40: We considered sensing matrices of row dimension m = 0.5n, 0.6n, 0.7n,
0.8n, where n = 40. For every matrix size, we randomly generated 10 different realizations of Gaussian
and partial Fourier matrices. So in total we used 80
different n = 40 sensing matrices for the numerical experiments in Tables 7 and 8. We normalized all of the
matrix columns so that they have a unit ℓ2 -norm. The
entries of Gaussian matrices were i.i.d standard Gaussian N (0, 1). The partial Fourier matrices had m rows
randomly draw from the full Fourier matrices. We compared our algorithms - pick-1-element, pick-2-element,
pick-3-element and TSA - to LP and SDP methods.
For readability, we place the numerical results for these
small sensing matrices in Appendix B.
For each matrix size and type, we increased k from
1 to 5 in unit steps. Tables 7 (a) and 8 (a) show the median values of αk . (To be consistent with the previous
research [3], in which the authors used the median value

9

of αk to compare the SDP method with the LP method,
we provided the median values obtained from 10 random realizations of sensing matrix.) From the median
value of αk , we obtained the recoverable sparsity kmax
such that αkmax < 1/2 and αkmax +1 > 1/2. In addition, we calculated the arithmetic mean of kmax ’s. For
the arithmetic mean, we obtained each kmax from each
random realization, and computed the arithmetic mean
of ten kmax ’s. Compared with LP and SDP methods, we
obtained bigger or at least the same recoverable sparsity
kmax by using pick-2, pick-3 and TSA. It is noteworthy
that we obtained the exact αk for k = 1, 2, ..., 5 by using TSA, while LP and SDP methods only provided the
exact αk for k = 1. We observed that αk < 1/2 but the
upper bound of αk > 1/2 holds true in several cases,
e.g., α5 in 32 × 40 Gaussian matrices, α4 in 28 × 40
Gaussian matrices, α3 in 24 × 40 Gaussian matrices, α3
in 20 × 40 partial Fourier matrices, and α4 in 24 × 40
partial Fourier matrices. Additionally, this can also be
established by the arithmetic mean of kmax in Tables 7
(a) and 8 (a).
To compare the computational complexity, we calculated the geometric mean of the algorithms’ execution time, to avoid biases for the average. Tables 7 (b)
and 8 (b) list the average execution time. We also ran
the Exhaustive Search Method (ESM) to find αk , and
compared its execution time with that of TSA. In calculating α5 , on average, 3-Step TSA reduced the computational time by around 86 times for 20 × 40 Gaussian
matrices, and by 94 times for 20 × 40 partial Fourier
matrices, compared to ESM. For 32 × 40 Gaussian matrix and partial Fourier matrix, the speedup compared
to the best l-Step TSA, l = 1, 2, 3, becomes around
1760 times and 182 times respectively. We observed
that when m/n = 0.5 e.g. 20 × 40 sensing matrices, in
general, the 3-step TSA provides the fastest result for
k = 5. On the other hand, for m/n = 0.8 (e.g. 32 × 40
case), the 2-Step TSA is the quickest in finding an exact αk for k = 5; however, for k > 5, the fastest l-step
TSA cannot be determined from either experiments or
theory.

Sensing matrices with n = 256: We assessed the
performance of the pick-l algorithm for sensing matrices with n = 256. We carried out numerical experiments on 128 × 256 Gaussian matrices in Fig. 2 (a) and
64 × 256 partial Fourier matrices in Fig. 2 (b). Here,
for 10 sensing matrices, we obtained the median value
of upper bounds of αk using the pick-l algorithm and
compared the result with LP relaxation method [8]. We
1
We conducted our experiments on HP Z220 CMT with
omitted SDP method [3] from this experiment due to
Intel Core i7-3770 dual core CPU @3.4GHz clock speed and
its very high computational complexity. For the pick-3
16GB DDR3 RAM, using Matlab (R2013b) on Windows 7.
2
LP method from http://www2.isye.gatech.edu/~ nemirovs/ algorithm in Fig. 2 (a), we calculated an upper bound
and SDP method from http://www.di.ens.fr/~ aspremon/NSPcode.html.
of α3 via TSA, and used this result to calculate upper

10

Myung Cho et al.

Table 2: Lower bound on k and execution time (Gaussian Matrix with n = 1024)
(a) Lower bound on k
matrix A
102 × 1024
205 × 1024
307 × 1024
410 × 1024
512 × 1024
614 × 1024
717 × 1024
819 × 1024
922 × 1024

(a) 128 × 256 Gaussian matri-(b) 64 × 256 Partial Fourier
ces
matrices

Fig. 2: Median upper bounds of αk from the pick-l algorithm
and the LP relaxation method.

Pick-1
2
4
5
7
9
12
16
21
32

k(U B(α2 )b )
3
4
6
8
10
13
17
23
36

k(α1 )
2
4
5
7
9
12
15
20
30

LPa
2
4
5
7
9
12
16
21
32

(b) Execution time (Unit: second)

Table 1: Lower bound on k and execution time (Gaussian Matrix with n = 512)

matrix A
102 × 1024
205 × 1024
307 × 1024
410 × 1024
512 × 1024
614 × 1024
717 × 1024
819 × 1024
922 × 1024

(a) Lower bound on k
matrix A
102 × 512
205 × 512
307 × 512
410 × 512

Pick-1
2
5
10
14

Pick-2
3
7
17
27

LPa
2
5
10
14
a

(b) Execution time (Unit: second)
matrix A
102 × 512
205 × 512
307 × 512
410 × 512

Pick-1
53.7
114.8
309.7
133.1

Pick-2
2.96e4
6.36e4
1.19e5
5.03e4

b

LPa
50.8
105.1
333.0
510.0

a Linear Programming [8]

bounds of αk , k = 3, 4, ..., 8 via (13). Fig. 2 (a) and
(b) demonstrate that, with an appropriate choice of l,
the upper bound of αk obtained via the pick-l algorithm can be tighter than that from the LP relaxation
method. For example, for 128 × 256 Gaussian matrices,
LP relaxation often determines the maximum recoverable sparsity as 5, while the pick-2 algorithm improves it
to 6. In the pick-3 algorithm, the maximum recoverable
sparsity is 7 (α7 = 0.49). For 64 × 256 partial Fourier
matrices, the maximum recoverable sparsity from LP
relaxation and the pick-2 algorithm are 3 and 4 respectively.
Sensing matrices with n = 512: We further conducted numerical experiments on Gaussian sensing matrices with n = 512. The simulation results in Table 1
clearly demonstrate that the pick-2 algorithm provides
larger lower bound on the recoverable sparsity k than
the LP method [8]. Especially, when Gaussian sensing
matrix is 410×512, the lower bound on k obtained from
the pick-2 algorithm is almost twice larger than that of
the LP method.

Pick-1
237
452
796
1207
1952
2150
1337
838
386

k(U B(α2 )b )
24 hours
24 hours
24 hours
24 hours
24 hours
24 hours
24 hours
24 hours
24 hours

k(α1 )
237
452
796
1207
1952
2150
1337
838
386

LPa
200
429
723
1073
1600
2217
2992
3904
4730

Linear Programming [8]
Upper bound of α2 obtained from 1-Step TSA after 24 hours’ run

5.1.2 High-dimensional sensing matrices
Sensing matrix with n ≥ 1024: We conducted numerical experiments for Gaussian sensing matrices with n
from 1024 to 6144. We show these numerical experiments in Tables 2 and 3, where we calculated the lower
bound on the recoverable sparsity k and obtained the
corresponding execution time. The SDP method [3] was
not applicable in these experiments due to its very high
computational complexity. In Table 2, we ran TSA for
1 day (24 hours) and obtained an upper bound of α2 ,
denoted by U B(α2 ). With the upper bound of α2 , we
obtained a lower bound of k, denoted by k(U B(α2 )),
via Lemma 3. Our numerical results in Tables 2 and
3 clearly show that our pick-l algorithm outperforms
the LP method in recoverable sparsity k or execution
time. We note that although our pick-1-element algorithm provides the same recoverable sparsity k as the
LP method [8] in Tables 2 and 3, the complexity of LP
method can be 10 times higher than our method on
m × n Gaussian matrices, where m is large.
For extremely large sensing matrices, e.g, 4014 ×
4096 and 6021×6144, the LP and SDP methods cannot
provide any lower bound on k due to unreasonable computational time. However, our pick-l algorithm can still
provide the lower bound on k efficiently. Table 3 shows
the lower bound on k and the execution time for these

Computable performance guarantees for compressed sensing matrices

11

Table 3: Lower bound on k and execution time (Gaussian Matrix)

Table 4: αk comparison and execution time (Gaussian
Matrix)

(a) Lower bound on k

(a) αk comparison

matrix A
512 × 2048
2007 × 2048
4014 × 4096
1024 × 6144
6021 × 6144
6134 × 6144

Pick-1
7
102
152
8
190
558

k(α1 )
6
90
139
8
174
406

a

LP
7
102
N/Ab
8
N/A
N/A

(b) Execution time (Unit: second)
matrix A
512 × 2048
2007 × 2048
4014 × 4096
1024 × 6144
6021 × 6144
6134 × 6144
a

Pick-1
7.51e3
6.71e2
9.12e3
2.18e5
3.89e4
1.37e4

k(α1 )
7.51e3
6.71e2
9.12e3
2.18e5
3.89e4
1.37e4

b
Linear Programming [8]
Not Available
c
Estimated time (15 hours for 4% calculations)
d
Estimated time (15 hours for 1% calculations)
e
Estimated time (10 hours for 1% calculations)

LP
6.63e3
7.19e4
15 daysc
1.61e5
65.5 daysd
41.7 dayse

matrix A
28 × 40
40 × 50

Algo.
Basic pick-3
Optimized pick-3
3-Step TSA
Basic pick-3
Optimized pick-3
3-Step TSA

α4
0.52
0.52
0.47
0.40
0.39
0.36

α5
0.64
0.63
0.54
0.48
0.47
0.41

α6
0.75
0.75
0.62
0.57
0.55
0.46

α7
0.86
0.85
0.67
0.65
0.62
0.51

α8
0.97
0.96
0.72-0.78
0.72
0.70
0.57-0.59

(b) Execution time (Unit: second)
matrix A

Algo.
Basic pick-3
Optimized pick-3
Basic pick-3
40 × 50
Optimized pick-3
28 × 40

α4
249.28
420.97
748.88
3.31e3

α5
249.28
410.43
748.88
3.49e3

α6
249.28
422.14
748.88
3.26e3

α7
249.28
422.41
748.88
3.26e3

α8
249.28
460.52
748.88
3.31e3

better or at least equal upper bound on αk to the basic
pick-l algorithm, with additional complexity. In spite of
the increased complexity of the optimized pick-l algorithm, it has an important theoretical merit, which is
Lemma 6.

large dimensional matrices, where our verified recoverable sparsity k can be as large as 558 for a 6134 × 6144
5.3 Complexity of tree search algorithm
sensing matrix. We obtained the estimated time for the
LP method by running the Matlab code obtained from
In this subsection, we carried out numerical experihttp://www2.isye.gatech.edu/~nemirovs/, which showsments to demonstrate the computational complexity of
the percentage of the calculation on screen.
TSA empirically on randomly chosen Gaussian sensing
matrices. Fig. 3 (a) and (b) show the distribution of
execution time and the distribution of number of nodes
in height 5 attached to the tree structure in TSA re5.2 Comparison between the optimized pick-l
spectively. For m = 0.5n, we generated 100 random realgorithm and the basic pick-l algorithm
alizations of Gaussian matrices and computed α5 using
3-Step TSA. The maximum
number
of leaf node whose
We compared the basic pick-l algorithm introduced in


n
40
cardinality
is
k
is
=
=
6.58008e5.
From Fig.
Section 2 to the optimized pick-l algorithm in Section
k
5
3 (b), we note that for 90 % of the cases, 3-Step TSA
3 on Gaussian sensing matrices 28 × 40 and 40 × 50
was terminated before 1.6 % of all the possible height-5
for l = 3 and k = 4, 5, ..., 8. Table 4 demonstrates that
nodes were attached to the tree structure.
when l = 3 and k = 4, 5, ..., 8, the optimized pick-l
We provided the execution time of TSA for differentalgorithm provided tighter upper bounds on αk than
sized randomly chosen Gaussian matrices in Fig. 4. We
the basic pick-l algorithm. This is because when l is
compared the execution time of TSA to ESM. Fig. 4 (a)
large and k > l, (15) includes more constraints, which
shows that when k = 1, 1-Step TSA provides almost
leads to the reduced size of the feasible set, than the
similar
performance to ESM. This is because 1-Step
case when k and l are small. Hence, the optimal value
TSA calculates all the α1,{i} ’s as a pre-computation,
of (15), which is the result from the optimized pick-l,
which is the same procedure as ESM. However, for k > l
can be smaller than or equal to that of (11), which is
as shown in Fig. 4 (b), (c), and (d), TSA can find αk
the basic pick-l. Additionally, we provided the exact αk
with reduced computation by using all the αl,L ’s, while
values obtained from TSA in order to check how tight
it is required to compute all the αk,K ’s in ESM. In order
the bounds obtained from the basic pick-l and the opto compute αk , we achieved a speedup of around 100
timized pick-l are. In terms of the execution time, the
times via 2-Step TSA compared to ESM for k = 3, 4.
optimized pick-l algorithm, which computes (15), was
around 1.7 and 4.4 times slower than the basic pick-l
In addition, in Fig. 5, we compared the execution
on 28 × 40 and 40 × 50 Gaussian matrix respectively.
time of TSA to ESM by varying k with n fixed on ranIn summary, the optimized pick-l algorithm provides
dom Gaussian matrices. For the best execution time of

12

Myung Cho et al.

(a)

(b)

Fig. 3: Histograms of the TSA (based on the pick-3 algorithm) to find α5 on 100 randomly chosen 20 × 40 Gaussian
sensing matrices for each method. (a) Execution time. (b)
Number of nodes in height 5.

(a) k = 1

(b) k = 2

(a) 80 × 100

(b) 160 × 200

Fig. 6: Global Lower Bound (GLB) and Global Upper
Bound (GUB) in TSA on Gaussian sensing matrices. (a)
For (k, l) = (5, 3), we obtained (GLB,GUB)=(0.27,0.28) after
167501 iterations. (b) For (k, l) = (4, 2), we obtained (GLB,
GUB)=(0.15,0.17) after 148101 iterations.

Finally, Fig. 6 gives illustrations of the values of
the global lower and upper bounds, for 80 × 100 and
160 × 200 Gaussian sensing matrices, as the number of
iterations in TSA increases. As we can see, the global
upper and lower bounds get close very quickly. This
implies that we can sometimes terminate TSA early
and still obtain tight bounds on αk .

5.4 Application to network tomography problem

(c) k = 3

(d) k = 4

Fig. 4: The execution time of TSA in log scale as a function
of n on randomly chosen m × n Gaussian matrices, where
m = n/2.

(a) n = 40

(b) n = 50

Fig. 5: The execution time of TSA in log scale as a function
of k on randomly chosen m × n Gaussian matrices, where
m = n/2.

TSA, we used different l values for TSA. For n = 40
and n = 50, 3-Step TSA reduced the execution time to
find α5 by around 100 times and 300 times respectively,
compared with ESM .

We apply our new tools introduced in this paper to verify NSC for sensing matrices in network tomography
problems [14,15,16,19,20,21]. In an undirected graph
model for the communication network, the communication delay over each link can be determined by sending
packets through probing paths that are composed of
connected links. The delay of each path is then measured by adding the delays over its links. Generally
most links are uncongested, and only a few congested
links have significant delays. It is, therefore, reasonable
to think of finding the link delays as a sparse recovery problem. This sparse problem can be expressed in
a system of linear equations y = Ax, where the vector
y ∈ Rm is the delay of m paths, the vector x ∈ Rn
is the delay vector for the n links, and A is a sensing matrix. The element Aij of A is 1, if and only
if path yi , i ∈ {1, 2, ..., m}, goes through link j,
j ∈ {1, 2, ..., n}; otherwise Aij equal to 0 (see Fig. 7).
The indices of nonzero elements in the vector x correspond to the congested links.
In our numerical experiments to verify NSC in network tomography problems, the paths for sending data
packets were generated by random walks of fixed length.
Table 5 summarizes the results of our experiments. We
note that by using TSA, one can exactly verify that a
total of k = 2 and k = 4 congested link delays can be
uniquely found by solving ℓ1 minimization problem (2)

Computable performance guarantees for compressed sensing matrices

13

Table 6: αk and execution time in a large network
model having 300 nodes and 400 edges
(a) αk values
Algo.
α1 α2 α3 α4
α5
α6
1-Step TSA 0.07 0.13 0.15 0.18
0.20
0.22-0.26a
2-Step TSA 0.07 0.13 0.15 0.18 0.20-0.23a 0.22-0.28a
(b) Execution time (Unit: second)

Fig. 7: (a) A simple example of a network tomography graph.
W , X, Y , and Z are nodes in the network, and Path1, 2,
3, and 4 are the probing paths through which the packets
are sent. (b) The sensing matrix corresponding to the graph
shown in (a). The rows and columns of the matrix represent
probing paths and edges respectively.

Table 5: αk and execution time in network tomography
problems
(a) αk values
matrix A

Algo.
1-Step TSA
33 × 66 2-Step TSA
3-Step TSA
1-Step TSA
53 × 105 2-Step TSA
3-Step TSA

a

α1
0.28
0.28
0.28
0.20
0.20
0.20

α2
0.41
0.41
0.41
0.29
0.29
0.29

α3
0.50
0.50
0.50
0.36
0.36
0.36

α4
α5
kmax
0.57
0.62
2
0.57 0.62-0.64
2
0.57
0.62
2
0.45 0.52-0.54
4
0.45 0.49-0.56
4
0.45
0.52
4

Random walk step: 20

(b) Execution time (Unit: second)
matrix A

Algo.
1-Step TSA
2-Step TSA
33 × 66
3-Step TSA
ESM
1-Step TSA
2-Step TSA
53 × 105
3-Step TSA
ESM

α1
0.74
0.74
0.74
0.64
1.31
1.31
1.31
1.28

α2
3.62
3.62
3.62
3.94
30.61
116.12
116.12
127.28

α3
28.94
43.94
1.69e3
1.63e3
608.90
143.99
7.95e3
8.70.e3

α4
404.11
541.70
1.73e3
1.4e4a
5.35e3
1.05e3
7.93e3
9.6e4a

α5
5.94e4
1 day
3.70e4
1.8e5a
1 day
1 day
1.38e4
1.9e6a

a

Exhaustive search method (Estimated execution time = average time to solve (7) (=0.02 second) for an index set × total
number of index sets)

for the randomly generated network measurement matrices 33 × 66 (12-node complete graph) and 53 × 105
(15-node complete graph) respectively. For ESM, we estimated the execution time by multiplying the unit time
to solve (7) and the total number of cases in the exhaustive search. We obtained the unit time to solve (7) by
calculating the arithmetic mean from 100 trials. For a
53 × 105 matrix, 3-Step TSA substantially reduced the
execution time to find α5 around 137 times compared
to ESM.
We further carried out numerical experiments on
even larger network model having 300 nodes and 400
edges. We created a random spanning tree for a network
model by using random walk approach [24]. At each

Algo.
α1
α2
α3
α4
α5
α6
1-Step TSA 63.37 65.70 599.96 5.49e3 8.60e4
1 day
2-Step TSA 63.37 3.46e4 3.54e4 4.03e4
1 day
1 day
ESM
73.22 3.20e4 1.59e6b 1.58e8b 1.25e10b 8.22e11b
a

Lower bound - upper bound
Exhaustive search method (Estimated Operation time = average time to solve (7) (=0.15 second) for an index set × total
number of index sets)
b

probing path, we randomly chose a node among 300
nodes as a starting point of random-walk and walked
100 times along the network connection. We obtained a
320 × 400 matrix corresponding to the network model.
We calculated αk values via l-Step TSA, where l = 1, 2.
In terms of the execution time, in Table 6, we compared TSA with ESM, where the unit time to solve (7)
was obtained by calculating the arithmetic mean from
100 trials. Especially, 1-Step TSA reduced the execution time to find α4 by around 28700 times compared
to ESM.
5.5 Discussion
In this section, we discuss the strengths and weaknesses
of our proposed algorithms, compared with earlier research [3,8].
1. Comparisons with LP and SDP. Our proposed
pick-1-element algorithm can achieve similar performance as the LP [8] and SDP methods [3]. However, our pick-1-element algorithm has the clear advantage of being more computationally efficient for
large dimensional sensing matrices. Please see Table
3, where the LP and SDP methods cannot provide
the performance bounds on recoverable sparsity k
due to high computational complexity. On the other
hand, in Table 3, our pick-1-element algorithm can
efficiently provide bounds on recoverable sparsity k.
The LP method has high computational complexity
because it has to deal with a large convex program
of design dimension mn, which leads to prohibitive
computational complexity when m and n are large
[8].
In our pick-1-element algorithm, we proposed the
novel idea of sorting α1,{i} ’s (see Lemma 2), which

14

leads to improved performance bounds on αk and
recoverable sparsity k. This sorting idea, combined
with Lemma 2, provides us with larger recoverable
sparsity bound k, than purely using α1 for bounding
recoverable k in [8, Section 4.2.B].
2. Set-specific upper bounds. Our proposed pickl-element algorithm (l ≥ 2) is novel, and can provide improved bounds on αk and recoverable sparsity k, using polynomial computational complexity
in n when l is fixed. This approach is not practical
when l is large. However, pick-2-element and pick3-element algorithm can already provide improved
performance bounds, compared with the previous
research [8,3].
The fact that we can obtain upper bounds on αk ,
based on the results of pick-l-element (l ≥ 2) algorithm, is new and non-trivial (see Lemma 2, Lemma
3 and Lemma 4). For example, if we know α5 ≤ 0.22,
we can use Lemma 3 to obtain that α11 ≤ 0.22 ×
11/5 < 0.5.
Our pick-l-element algorithm can provide set-specific
upper bound for αk,K , laying the foundation for our
branch-and-bound TSA.
3. Computational complexity of TSA. We proposed TSA to find precise values for αk with significantly reduced average-case computational complexity than ESM. The computational complexity
of TSA is dependent on n, sparsity k, and a chosen
constant l. When k, n and l are large enough, finding αk via TSA is still computationally expensive.
In the worst case, TSA has the same computational
complexity as ESM. However, our extensive simulations ranging from Fig. 3 to Fig. 5 and from Table
5 to Table 8 show that on average, TSA can greatly
reduce the computational complexity of finding αk
compared with ESM.
Moreover, since TSA maintains an upper bound and
a lower bound of αk during its iterations, one can
always early terminate TSA, and still get improved
performance bounds on αk than the LP and SDP
methods. We can use TSA to find an exact value of
αl , where l < k, and then use Lemma 3 to bound
αk .
4. Use of data structures. We used Object-Oriented
Programming (OOP) to implement TSA in Matlab
[25], because the OOP makes it easy to handle treetype structures. In OOP, we defined a class and created objects from the class to store property of each
node J, e.g., B(J), in the tree. In order to make a
connection between two tree nodes, we used doubly
linked list data structure as a part of the object.
However, in case readers would like to implement
the algorithm using alternative data structures, we

Myung Cho et al.

have provided implementation-agnostic pseudocode
of our algorithm in Algorithm 2.
5. Difference from phase transition works. There
has been extensive research on the phase transitions of various sparse recovery algorithms such as
Basis Pursuit (BP), Orthogonal Matching Pursuit
(OMP), and Approximate Message Passing (AMP)
[1]. However, our research is different from the research on phase transition in two aspects. Firstly,
our work and the previous works [8,3] are focusing
on worst-case performance guarantee (recovering all
the possible k-sparse signals), while the research
on phase transition is considering the average-case
performance guarantee for a single k-sparse signal
with fixed support and sign pattern. Secondly, the
phase transition bounds are mostly for random matrices. Hence, for a given deterministic sensing matrix, phase transition results cannot be used for that
particular matrix.

6 Conclusion
In this paper, we consider the problem of verifying the
null space condition in compressed sensing. Calculating the proportional parameter αk that characterizes
the null space condition of a sensing matrix is a nonconvex optimization problem, and also known to be NPhard in [7]. In order to verify the null space condition,
we proposed novel and simple enumeration-based algorithms, which are called the basic and optimized pick-l
algorithms, to obtain upper bounds of αk . With these
algorithms, we further designed a new algorithm called
the tree search algorithm to gain a global solution to
the non-convex optimization problem of verifying the
null space condition. Numerical experiments show that
our algorithms outperform the previously proposed algorithms [3,8] in performance as well as speed.

Declarations
Availability of data and material
All the codes used for the numerical experiments are available
at the following link:
https://sites.google.com/view/myungcho/software/nsc.

Acknowledgements
We thank Alexandre d’Aspremont from CNRS at Ecole Normale Superieure, Anatoli Juditsky from Laboratoire Jean Kuntzmann at Universite Grenoble Alpes, and Arkadi Nemirovski
from Georgia Institute of Technology for helpful discussions
and providing codes for the simulations in [3] and [8].

Computable performance guarantees for compressed sensing matrices

Funding
The work of Weiyu Xu is supported by Simons Foundation
318608, KAUST OCRF-2014-CRG-3, NSF DMS-1418737 and
NIH 1R01EB020665-01.

Competing interests
The authors declare that they have no competing interests.

Authors’ contributions
Myung Cho and Weiyu Xu designed the algorithms. Myung
Cho implemented the algorithms. Kumar Vijay Mishra checked
the implementation of the algorithms and helped to polish the
manuscript. All authors read and approved the final manuscript.

References
1. YC Eldar, G Kutyniok, Compressed Sensing: Theory and
Applications (Cambridge University Press, 2012).
2. EJ Candès, T Tao, Decoding by linear programming.
IEEE Trans. Inf. Theory. 51(12):4203–4215 (2005).
3. A d’Aspremont, L El Ghaoui, Testing the nullspace property using semidefinite programming. Math. prog. 127(1),
123-144 (2011).
4. EJ Candès, J Romberg, T Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete
frequency information. IEEE Trans. Inf. Theory. 52(2),
489-509 (2006).
5. EJ Candès, J Romberg, T Tao, Stable signal recovery
from incomplete and inaccurate measurements. Comm.
on pure and appl. math. 59(8), 1207-1223 (2006).
6. D Donoho, Neighborly polytopes and sparse solution
of underdetermined linear equations. Technical report
(Stanford University. Dept. of Statistics 2005).
7. AM Tillmann, ME Pfetsch, The computational complexity of the restricted isometry property, the nullspace
property, and related concepts in compressed sensing.
IEEE Trans. Inf. Theory. 60(2), 1248-1259 (2014).
8. A Juditsky, A Nemirovski, On verifiable sufficient conditions for sparse signal recovery via ℓ1 minimization.
Math. prog. 127(1), 57-88 (2011).
9. A Cohen, W Dahmen, R DeVore, Compressed sensing
and best k-term approximation. J. American Math. Soc.
22(1),211-231 (2009).
10. W Xu, B Hassibi, Precise stability phase transitions for
ℓ1 minimization: a unified geometric framework. IEEE
Trans. Inf. Theory. 57(10), 6894-6919 (2011).
11. K Lee, Y Bresler, Computing performance guarantees for
compressed sensing. In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (2008), pp. 5129-5132.
12. G Tang, A Nehorai, Verifiable and computable ℓ∞ performance evaluation of ℓ1 sparse signal recovery. In Proceedings of Conference on Information Sciences and Systems
(CISS) (2011), pp. 1-6.
13. M Cho, W Xu, New algorithms for verifying the null
space conditions in compressed sensing. In Proceedings of
Asilomar Conference on Signals, Systems and Computers
(2013), pp. 1038-1042.

15

14. W Xu, E Mallada, A Tang, Compressive sensing over
graphs. In Proceedings of IEEE International Conference
on Computer Communizations (INFOCOM) (2011), pp.
2087-2095.
15. MH Firooz, S Roy, Network tomography via compressed
sensing. In Proceedings of IEEE Global Telecommunications Conference (GLOBECOM) (2010), pp. 1-5.
16. MJ Coates, RD Nowak, Network tomography for internal
delay estimation. In Proceedings of IEEE International
Conference on Acoustics, Speech, and Signal Processing
(ICASSP) (2001), vol. 6, pp. 3409-3412.
17. M Grant, S Boyd, CVX: Matlab software for
disciplined convex programming, version 2.1 beta.
http://cvxr.com/cvx (2012).
18. MOSEK ApS, The MOSEK optimization toolbox
for MATLAB manual. Version 7.1 (Revision 31).
http://docs.mosek.com/7.1/toolbox/index.html
(2015).
19. Y Tsang, M Coates, RD Nowak, Network delay tomography. IEEE Trans. Signal Process. 51(8), 2125-2136
(2003).
20. Y Vardi, Network tomography: estimating sourcedestination traffic intensities from link data. J.American
Stat. Assoc. 91(433), 365-377 (1996).
21. R Castro, M Coates, G Liang, R Nowak, B Yu, Network
tomography: recent developments. Stat. Sci. 19(3), 499517 (2004).
22. DL Donoho, J Tanner, Precise Undersampling Theorems.
Proc. IEEE 98(6), 913-924 (2010).
23. G Brassard, P Bratley, Fundamentals of algorithmics
(Englewood Cliffs: Prentice Hall, 1996).
24. DB Wilson, Generating random spanning trees more
quickly than the cover time. In Proceedings of the twentyeighth annual ACM symposium on Theory of computing
(1996), pp. 296-303.
25. Mathworks,
Object-Oriented
Programming
in
MATLAB.
Retrieved
May
26,
2017,
from
https://www.mathworks.com/discovery/objectoriented-programming.html

A Proof of Proposition 1
Proof Let us denote the sum of k maximal magnitudes of
elements of x ∈ Rn as
X
|xi |.
||x||k,1 = maximize
|K|≤k

i∈K

We use i1 , i2 , ..., and ik (or j1 , j2 , ..., and jk ) to denote k
distinct integers between 1 and n. For a matrix, say A, we use
Ai,j to represent its element in the i-th row and j-th column.


T
=
minimize
maximize
||(I
−
Y
A)e
||
αLP
j
k,1
k
m×n
1≤j≤n

Y =[y1 ,...,yn ]∈R

=

≤

=

k
X

|(I − Y T A)it ,j |



|(I − Y T A)it ,jt |



minimize



minimize



minimize



k
X
maximize
||eit − AT yit ||∞

Y =[y1 ,...,yn ]∈Rm×n

Y =[y1 ,...,yn ]∈Rm×n

Y =[y1 ,...,yn ]∈Rm×n

= maximize
i1 ,i2 ,...,ik

X
k

maximize

i1 ,i2 ,...,ik ,j

maximize

t=1
k
X

i1 ,i2 ,...,ik ,
j1 ,j2 ,...,jk t=1

i1 ,i2 ,...,ik

t=1

(minimize
||eit − AT yit ||∞ )
m×1

t=1

yit ∈R



,

(22)

16

Myung Cho et al.

where we can exchange the order of “maximize” and
“minimize” in the last equality because ||eit − AT yit ||∞
only depends on yit .
Moreover, according to the equations for “αi ” between (4.29) and (4.30) in [8] (taking β there to be
∞),

Table 7: αk comparison and execution time - Gaussian
Matrix
(a) αk comparison
A

(m×n)

minimize
||eit − AT yit ||∞
m×1
yit ∈R

20 × 40


= maximize eT
it x : Ax = 0, ||x||1 ≤ 1
x

= α1,{it } .

Combining this with (22), αLP
is no bigger than
k
the upper bound calculated by Lemma 2 (based on the
pick-1-element algorithm). Namely,
≤ αpick1
.
αLP
k
k

24×40

(23)
⊓
⊔
28×40

B Sensing matrices with n = 40
Here, we provide the numerical results for small sensing matrices with n = 40 to compare our methods to LP [8] and
SDP [3] methods.

32×40

Algo.

α1

α2

α3

α4

α5

kmax d

pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LPa
SDPb
ESMc
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM

0.28
0.28
0.28
0.28
0.28
0.28
0.28
0.28
0.28
0.23
0.23
0.23
0.23
0.23
0.23
0.23
0.23
0.23
0.18
0.18
0.18
0.18
0.18
0.18
0.18
0.18
0.18
0.14
0.14
0.14
0.14
0.14
0.14
0.14
0.14
0.14

0.55
0.45
0.45
0.45
0.45
0.45
0.50
0.49
0.45
0.46
0.37
0.37
0.37
0.37
0.37
0.41
0.41
0.37
0.36
0.31
0.31
0.31
0.31
0.31
0.34
0.34
0.31
0.29
0.24
0.24
0.24
0.24
0.24
0.27
0.27
0.24

0.81
0.66
0.57
0.57
0.57
0.57
0.67
0.66
0.57
0.67
0.53
0.46
0.46
0.46
0.46
0.56
0.55
0.46
0.53
0.46
0.41
0.41
0.41
0.41
0.49
0.48
0.41
0.42
0.37
0.33
0.33
0.33
0.33
0.38
0.38
0.33

1
0.85
0.76
0.67
0.67
0.67
0.84
0.81
0.67
0.87
0.69
0.61
0.57
0.57
0.57
0.71
0.70
0.57
0.70
0.59
0.54
0.49
0.49
0.49
0.61
0.60
0.49
0.55
0.47
0.44
0.40
0.40
0.40
0.49
0.48
0.40

1
1
0.92
0.75
0.75
0.75
0.98
0.95
0.75
1
0.85
0.75
0.65
0.65
0.65
0.84
0.82
0.65
0.86
0.72
0.66
0.57
0.57
0.57
0.72
0.71
0.57
0.67
0.58
0.53
0.47
0.47
0.47
0.58
0.57
0.47

1/1.1
2/1.9
2/1.9
2/1.9
2/1.9
2/1.9
2/1.6
2/1.8
2/1.9
2/2.0
2/2.1
3/2.8
3/2.8
3/2.8
3/2.8
2/2.0
2/2.0
3/2.8
2/2.0
3/3.0
3/3.0
4/3.5
4/3.5
4/3.5
3/3.0
3/3.0
4/3.5
3/3.0
4/3.8
4/4.2
5/4.9
5/4.9
5/4.9
4/3.9
4/4.0
5/4.9

(b) Execution time
A

(m×n)

20×40

24×40

28×40

32×40

a
c

Algo.

α1

α2

α3

pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Stpe TSA
2-Stpe TSA
3-Stpe TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM

0.35
0.35
0.35
0.50
0.50
0.50
0.55
56.92
0.35
0.44
0.44
0.44
0.50
0.50
0.50
0.84
62.18
0.44
0.58
0.58
0.58
0.52
0.52
0.52
1.12
71.27
0.58
0.42
0.42
0.42
0.55
0.55
0.55
0.70
56.12
0.42

0.35
10.96
10.96
2.14
13.20
13.20
0.55
6.02e3
10.96
0.44
13.00
13.00
2.05
12.92
12.92
0.94
5.59e3
13.00
0.58
14.67
14.67
1.41
13.54
13.54
1.20
5.55e3
14.67
0.42
13.29
13.29
1.14
14.22
14.22
0.71
7.17e3
13.29

0.35
10.96
313.65
11.78
14.11
320.20
0.58
5.14e3
313.65
0.44
13.00
311.27
9.63
13.60
319.27
0.88
4.89e3
311.27
0.58
14.67
326.80
4.39
13.82
327.79
1.12
4.90e3
326.80
0.42
13.29
331.80
2.89
14.32
340.87
0.72
5.43e3
331.80

(Unit: second)
α4
α5
0.35
10.96
313.65
128.98
58.93
346.43
0.55
5.12e3
4.5e3
0.44
13.00
311.27
77.45
35.08
378.10
0.83
4.75e3
4.6e3
0.58
14.67
326.80
32.43
29.35
404.23
1.09
4.98e3
4.7e3
0.42
13.29
331.80
13.50
18.13
336.29
0.70
5.07e3
4.9e3

b
Linear Programming [8]
Semidefinite Programming [3]
d
Exhaustive Search Method
median / arithmetic mean

0.35
10.95
313.65
1.62e3
3.77e3
695.53
0.56
5.61e3
6.0e4
0.44
13.00
311.27
429.48
634.62
481.29
0.82
5.37e3
6.4e4
0.58
14.67
326.80
119.86
126.62
383.61
0.68
4.72e3
6.9e4
0.42
13.29
331.80
40.67
40.35
355.06
0.70
4.79e3
7.1e4

Computable performance guarantees for compressed sensing matrices

Table 8: αk comparison and execution time - Partial
Fourier Matrix
(a) αk comparison
A

(m×n)

20×40

24×40

28×40

32×40

Algo.

α1

α2

α3

α4

α5

kmax d

pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LPa
SDPb
ESMc
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM

0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.12
0.12
0.12
0.12
0.12
0.12
0.12
0.12
0.12
0.09
0.09
0.09
0.09
0.09
0.09
0.09
0.09
0.09

0.39
0.36
0.36
0.36
0.36
0.36
0.39
0.39
0.36
0.31
0.27
0.27
0.27
0.27
0.27
0.31
0.31
0.27
0.25
0.23
0.23
0.23
0.23
0.23
0.25
0.25
0.23
0.19
0.17
0.17
0.17
0.17
0.17
0.19
0.19
0.17

0.59
0.55
0.47
0.47
0.47
0.47
0.59
0.59
0.47
0.47
0.42
0.38
0.38
0.38
0.38
0.47
0.47
0.38
0.37
0.35
0.32
0.32
0.32
0.32
0.37
0.37
0.32
0.29
0.27
0.25
0.25
0.25
0.25
0.29
0.29
0.25

0.78
0.73
0.64
0.61
0.61
0.61
0.78
0.78
0.61
0.62
0.55
0.51
0.49
0.49
0.49
0.62
0.62
0.49
0.50
0.47
0.44
0.41
0.41
0.41
0.50
0.50
0.41
0.38
0.36
0.35
0.33
0.33
0.33
0.38
0.38
0.37

0.98
0.91
0.80
0.70
0.70
0.70
0.98
0.98
0.70
0.78
0.69
0.64
0.59
0.59
0.59
0.78
0.78
0.59
0.62
0.58
0.54
0.50
0.50
0.50
0.62
0.62
0.50
0.48
0.44
0.43
0.39
0.39
0.39
0.48
0.48
0.39

2/2.0
2/2.2
3/2.7
3/2.7
3/2.7
3/2.7
2/2.0
2/2.0
3/2.7
3/2.8
3/3.0
3/3.4
4/3.5
4/3.5
4/3.5
3/2.8
3/2.8
4/3.5
4/3.6
4/4.0
4/4.0
4/4.0
4/4.1
4/4.1
4/3.6
4/3.6
4/4.1
5/4.7
5/4.7
5/4.7
5/4.7
5/4.7
5/4.7
5/4.7
5/4.7
5/4.7

(b) Execution time
A

(m×n)

20×40

24×40

28×40

32×40

a
c

Algo.

α1

α2

α3

pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM
pick-1
pick-2
pick-3
1-Step TSA
2-Step TSA
3-Step TSA
LP
SDP
ESM

0.31
0.31
0.31
0.47
0.47
0.47
0.49
33.93
0.31
0.39
0.39
0.39
0.48
0.48
0.48
0.62
41.13
0.39
0.43
0.43
0.43
0.50
0.50
0.50
0.65
40.51
0.43
0.57
0.57
0.57
0.52
0.52
0.52
0.86
46.51
0.57

0.31
10.85
10.85
9.72
11.97
11.97
0.77
2.34e3
10.85
0.39
11.51
11.51
12.12
12.52
12.52
0.56
2.39e3
11.51
0.43
13.29
13.29
8.70
12.99
12.99
0.67
2.17e3
13.29
0.57
17.24
17.24
6.39
13.38
13.38
0.89
2.41e3
17.24

0.31
10.85
260.41
70.57
18.54
291.29
0.53
2.65e3
260.41
0.39
11.51
302.86
76.21
21.46
306.43
0.66
2.66e3
302.86
0.43
13.29
341.05
31.53
16.85
317.40
0.71
2.29e3
341.05
0.57
17.24
385.26
22.35
18.65
326.40
0.78
2.62e3
385.26

(Unit: second)
α4
α5
0.31
10.85
260.41
329.28
45.18
297.45
0.59
2.91e3
4.1e3
0.39
11.51
302.86
407.67
107.00
426.17
0.59
2.63e3
4.5e3
0.43
13.29
341.05
272.68
47.45
410.47
0.67
2.80e3
4.7e3
0.57
17.24
385.26
101.67
49.46
476.55
0.75
2.53e3
4.8e3

b
Linear Programming [8]
Semidefinite Programming [3]
d
Exhaustive Search Method
median / arithmetic mean

0.31
10.85
260.41
3.60e3
3.36e3
633.12
0.51
2.60e3
6.0e4
0.39
11.51
302.86
2.77e3
1.83e3
1.36e3
0.58
2.56e3
6.4e4
0.43
13.29
341.05
731.90
544.79
553.67
0.75
2.63e3
6.5e4
0.57
17.24
385.26
451.62
372.35
1.02e3
0.76
2.75e3
6.8e4

17

