# TELET: A Monotonic Algorithm to Design Large Dimensional Equiangular Tight
Frames for Applications in Compressed Sensing

R. Jyothi and P. Babu

###### Abstract

An Equiangular tight frame (ETF) - also known as the Welch-bound-equality
sequences - consists of a sequence of unit norm vectors whose absolute inner
product is identical and minimal. Due to this unique property, these frames
are preferred in different applications such as in constructing sensing
matrices for compressed sensing systems, robust transmission and quantum
computing. Construction of ETFs involves solving a challenging non-convex
minimax optimization problem, and only a few methods were successful in
constructing them, albeit only for smaller dimensions. In this paper, we
propose an iterative algorithm named TEchnique to devise Large dimensional
Equiangular Tight-frames (TELET-frames) based on the majorization minimization
(MM) procedure - in which we design and minimize a tight upper bound for the
ETF cost function at every iteration. Since TELET is designed using the MM
approach, it inherits useful properties of MM such as monotonicity and
guaranteed convergence to a stationary point. Subsequently, we use the derived
frames to construct optimized sensing matrix for compressed sensing systems.
In the numerical simulations, we show that the proposed algorithm can generate
complex and real frames (in the order of hundreds) with very low mutual
coherence value when compared to the state-of-the-art algorithm, with a slight
increase in computational cost. Experiments using synthetic data and real
images reveal that the optimized sensing matrix obtained through the frames
constructed by TELET performs better, in terms of image reconstruction
accuracy, than the sensing matrix constructed using state-of-the-art methods.

###### Index Terms:

Equiangular Tight Frame, Majorization Minimization, Minimax Problem,
Compressed Sensing, Sensing Matrix

##  I Introduction

Orthonormal basis is a powerful tool and has applications in diverse fields
such as in computer vision, signal and image processing [1]. Sometimes,
instead of orthonormal basis, one may prefer the use of an overcomplete
spanning set of vectors. An example where such a situation occurs is in Gabor
analysis, wherein orthonormal Gabor basis with good time-frequency
localization do not exist, while it is not difficult to find overcomplete
Gabor system with excellent time-frequency localization [2]. The overcomplete
spanning set of vectors in finite dimensions is called a _Frame_. Formally, a
frame is defined as a sequence of vectors
$\\{{\boldsymbol{x}}_{i}\\}_{i=1}^{N}$ drawn from $\mathbb{C}^{d}$ ($d\leq N$)
that satisfies the generalized Parseval condition:

$\begin{array}[]{ll}a\|{\boldsymbol{v}}\|_{2}^{2}\leq\displaystyle\sum_{n=1}^{N}|{\boldsymbol{v}}^{H}{\boldsymbol{x}}_{i}|^{2}\leq
b\|{\boldsymbol{v}}\|_{2}^{2}\>\>{\textrm{for
all}}\>\>{\boldsymbol{v}}\in\mathbb{C}^{d}\end{array}$ (1)

where $\|\cdot\|_{2}$ is the Euclidian norm, $(\cdot)^{H}$ is the Hermitian
operator, $a$ and $b$ are the lower and upper frame bounds, respectively. If
$a=b$ in (1), then the frame is called $a$-tight frame, and if $a=b=1$, it is
a Parseval frame. When the frame is $a$-tight, the condition in (1) is
equivalent to the statement that
${\boldsymbol{X}}{\boldsymbol{X}}^{H}=a\,{\mathbf{I}}_{d}$, where
${\boldsymbol{X}}=[{\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\cdots,{\boldsymbol{x}}_{N}]\in\mathbb{C}^{d\times
N}$ is the frame synthesis operator. This implies that the $d$ non-zero
singular values of ${\boldsymbol{X}}$ are equal and is equal to $\sqrt{a}$.
Hence, being a tight frame is a spectral constraint on the matrix
${\boldsymbol{X}}$, which as discussed later in this section is used by the
state-of-the-art algorithms to construct $a$-tight frames. When each
${\boldsymbol{x}}_{i}$ of an $a$-tight frame has unit norm i.e. if
$\|{\boldsymbol{x}}_{i}\|_{2}=1$, then such an $a$-tight frame is called as a
Unit Norm Tight Frame (UNTF). These frames are commonly used in the
construction of signature sequences of CDMA systems [3]. In applications such
as sparse approximation [4], robust transmission [5], and quantum computing
[6], it is required that the vectors comprising the UNTFs be maximally
uncorrelated. The maximum correlation among the pair of vectors in a frame can
be estimated using the mutual coherence metric $\mu({\boldsymbol{X}})$ which
is defined as:

$\begin{array}[]{ll}\mu({\boldsymbol{X}})=\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\rm
maximize}\>\dfrac{\left|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}\right|}{\|{\boldsymbol{x}}_{i}\|_{2}\|{\boldsymbol{x}}_{j}\|_{2}}\end{array}$
(2)

Frames with small mutual coherence are known as incoherent frames. A lower
bound (known as the Welch bound) on the minimum achievable correlation for any
arbitrary frame is given as:

$\begin{array}[]{ll}\mu({\boldsymbol{X}})\geq\sqrt{\dfrac{N-d}{d(N-1)}}\end{array}$
(3)

Unit norm tight frames which achieve equality in the Welch bound are known as
Equiangular Tight Frame (ETF). Geometrically, it corresponds to a frame in
which each pair of distinct unit norm vectors meet at the same angle
$\theta={\rm{cos}}^{-1}\left(\sqrt{\dfrac{N-d}{d(N-1)}}\right)$. Note that,
all orthonormal bases are ETFs with $\mu({\boldsymbol{X}})$ equal to zero.
Hence, ETFs can be perceived as a generalization of orthonormal bases.
However, ETFs do not exist for all values of $d$ and $N$. A necessary
condition for the existence of a real ETF (${\boldsymbol{x}}_{i}$’s taking
only real values) is given by $N\leq\dfrac{1}{2}d(d+1)$ while a complex ETF
(${\boldsymbol{x}}_{i}$’s taking complex values) can only exist when $N\leq
d^{2}$ [7]. Construction of an ETF is challenging as it involves solving a
challenging non-convex minimax optimization problem. The existing methods are
based on the alternating projection framework and they focus on the
construction of the Gram matrix
${\boldsymbol{G}}={\boldsymbol{X}}^{H}{\boldsymbol{X}}$ of the ETF rather than
directly constructing the frame itself [8, 9, 10]. These methods consists of
two steps, where in Step. I, the algorithms [8, 9, 10] finds a Gram matrix
$\bar{{\boldsymbol{G}}}^{t}$ with the desired structural constraint using the
shrinkage operator defined as:

$\begin{array}[]{ll}\bar{{\boldsymbol{G}}}^{t}(i,j)=\left\\{\begin{array}[]{ll}G^{t}(i,j)&i\neq
j,|G^{t}(i,j)|\leq\eta\\\ \eta e^{i{\rm{arg}}G^{t}(i,j)}&i\neq
j,|G^{t}(i,j)|>\eta\\\ 1&i=j\\\ \end{array}\right..\end{array}$ (4)

In Step. II, the nearest $\alpha$-tight matrix
${{\boldsymbol{G}}^{t}_{\alpha}}$ to $\bar{{\boldsymbol{G}}}^{t}$ is obtained,
according to
${{\boldsymbol{G}}^{t}_{\alpha}}=\alpha{\boldsymbol{U}}{\boldsymbol{U}}^{*}$,
where ${\boldsymbol{U}}$ is a $N\times d$ matrix whose columns are the
eigenvectors corresponding to the $d$ largest eigenvalues of
$\bar{{\boldsymbol{G}}}^{t}$. Hence at every iteration, the algorithms in [8,
9, 10] alternates between finding a Gram matrix with the desired structural
and spectral constraint. After performing alternating projection, the tight
frames are obtained from the output Gram matrix using either eigenvalue
decomposition (EVD) or rank-revealing QR factorization. The algorithms in [8,
9, 10] differs essentially in the choice of the threshold $\eta$ used in the
shrinkage operator in (4) and the tightness constraint $\alpha$. The
algorithms in [8], [9] uses the threshold parameter $\eta$ as $\sqrt{1/d}$ and
[10] computes $\eta$ as $\sqrt{{(N-d)}/{d(N-1)}}$. The tightness constraint
$\alpha$ in the case of [8] and [10] is chosen as $\sqrt{{N}/{d}}$ while [9]
computes $\alpha=\displaystyle\sum_{k=1}^{d}\lambda_{k}/{d}$ where
$\lambda_{k}$’s are the eigenvalues of $\bar{{\boldsymbol{G}}}^{t}$. One of
the drawbacks of the alternating projection based methods is that they cannot
construct ETFs for large values of $N$ and $d$ [8]. This could be because of
the high complexity of the algorithm due to the computation of EVD at every
iteration. [11]. Different from the alternating projection based methods, the
authors in [12, 13, 14, 15, 16] proposed algorithms which constructs real
equiangular frames by directly modifying the frame synthesis operator
${\boldsymbol{X}}$ and not its Gram matrix ${\boldsymbol{G}}$. The authors in
[14] and [15] proposed an iterative algorithm wherein at every iteration the
frame synthesis operator is constructed column-by-column and each update
involves solving a convex optimization problem using interior point method -
which makes the algorithm computationally expensive. We next discuss an
application of ETFs to design the optimal sensing matrix in the compressed
sensing (CS) systems.

###  I-A Application to compressed sensing

Based on the theoretical findings of [17], ETFs can be employed in the CS
framework. In the following we briefly discuss the CS model and the need for
optimal sensing matrix in CS systems. The CS framework involves reconstructing
the high-dimensional $K$-sparse signal vector
${\boldsymbol{u}}\in\mathbb{C}^{N\times 1}$ from its corresponding low-
dimensional measurement vector ${\mathbf{y}}\in\mathbb{C}^{d\times 1}$
($d<<N$) via a carefully constructed sensing matrix
${\boldsymbol{\Theta}}\in\mathbb{C}^{d\times N}$ using the relation
${\mathbf{y}}={\boldsymbol{\Theta}}{\boldsymbol{u}}$. The vector
${\boldsymbol{u}}$ is $K$-sparse, since it can be represented using $K$ non-
zero coefficients over a fixed dictionary
${\boldsymbol{\Psi}}\in\mathbb{C}^{N\times N}$ i.e.
${\boldsymbol{u}}={\boldsymbol{\Psi}}{\boldsymbol{s}}$ (where
${\boldsymbol{s}}\in\mathbb{C}^{N}$ has $K$ non-zero entries i.e.
$\|{\boldsymbol{s}}\|_{0}=|\\{i:s_{i}\neq 0\\}|=K$). Hence, the low-
dimensional measurement vector ${\mathbf{y}}$ can be expressed as
${\mathbf{y}}={\boldsymbol{\Theta}}{\boldsymbol{\Psi}}{\boldsymbol{s}}$. Now,
given the measurement vector ${\mathbf{y}}$, the dictionary
${\boldsymbol{\Psi}}$ and the sensing matrix ${\boldsymbol{\Theta}}$, the
sparse signal vector ${\boldsymbol{s}}$ can be reconstructed by solving the
following $\ell_{0}$ minimization problem:

$\begin{array}[]{ll}\underset{{\boldsymbol{s}}}{\rm
minimize}\>\|{\boldsymbol{s}}\|_{0}\>\>\>\>{\rm{s.t.}}\>\>{\mathbf{y}}={\boldsymbol{\Theta}}{\boldsymbol{\Psi}}{\boldsymbol{s}}\end{array}$
(5)

The problem in (5) is NP-hard and can be solved using greedy algorithms such
as orthogonal matching pursuit (OMP) [18, 19] or Basis Pursuit (BP) algorithm
[20]. As shown in [17], the $K$-sparse signal ${\boldsymbol{s}}$ can be
uniquely recovered from the measurements
${\mathbf{y}}={\boldsymbol{\Theta}}{\boldsymbol{\Psi}}{\boldsymbol{s}}$ as a
solution to (5) as long as:

$\begin{array}[]{ll}\|{\boldsymbol{s}}\|_{0}<\dfrac{1}{2}\left[1+\dfrac{1}{\mu({\boldsymbol{X}}_{\rm{ED}})}\right]\end{array}$
(6)

where ${\boldsymbol{X}}_{\rm{ED}}={\boldsymbol{\Theta}}{\boldsymbol{\Psi}}$ is
called the equivalent dictionary. Note that the above upper bound is dictated
by mutual coherence of $\mu({\boldsymbol{X}}_{\rm{ED}})$, smaller the value of
$\mu({\boldsymbol{X}}_{\rm{ED}})$, larger the sparsity of signal
${\boldsymbol{s}}$ which can be successfully recovered. Also, as reported in
[18], the incoherence property of ${\boldsymbol{X}}_{\rm{ED}}$ plays an
important role in the performance of the OMP and BP algorithm. Hence, given
the dictionary ${\boldsymbol{\Psi}}$, we need to find the optimal sensing
matrix ${\boldsymbol{\Theta}}$ such that ${\boldsymbol{X}}_{\rm{ED}}$ is as
close as possible to an ETF.
Elad in [21] was the first to come up with the idea of optimizing the sensing
matrix ${\boldsymbol{\Theta}}$ such that the equivalent dictionary is
incoherent. However, instead of decreasing the mutual coherence of the
equivalent dictionary, they worked on decreasing its average mutual coherence.
Since then many different approaches have appeared. Majority of the approaches
[22, 23, 24, 25] designed the sensing matrix ${\boldsymbol{\Theta}}$ such that
the Gram matrix of the equivalent dictionary ${\boldsymbol{X}}_{\rm{ED}}$ is
close to a target Gram matrix ${\boldsymbol{G}}_{{\rm{Target}}}$. The authors
in [10], obtained the optimal sensing matrix by first constructing the
incoherent dictionary ${\boldsymbol{X}}_{\rm{ED}}$ as an ETF using the
alternating projection method (as discussed earlier) and then found the
optimal sensing matrix ${\boldsymbol{\Theta}}$ over the fixed dictionary by
solving the following least squares problem:

$\begin{array}[]{ll}\underset{{\boldsymbol{\Theta}}}{\rm
minimize}\|{\boldsymbol{X}}_{\rm{ED}}-{\boldsymbol{\Theta}}{\boldsymbol{\Psi}}\|_{F}^{2}\end{array}$
(7)

Sometimes, even though the optimal sensing matrix is constructed such that the
equivalent dictionary is as close as possible to an ETF, it may not perform
well. This could happen when the sparse representation error (SRE) defined as
${\boldsymbol{E}}={\boldsymbol{U}}-{\boldsymbol{\Psi}}{\boldsymbol{S}}$, where
${\boldsymbol{U}}=[{\boldsymbol{u}}_{1},{\boldsymbol{u}}_{2},\cdots{\boldsymbol{u}}_{R}]$
is the original signal set and ${\boldsymbol{S}}$ is its corresponding sparse
coefficient matrix
${\boldsymbol{S}}=[{\boldsymbol{s}}_{1},{\boldsymbol{s}}_{2},\cdots{\boldsymbol{s}}_{R}]$
is not taken into account. Projecting SRE to the measurement domain through
the sensing matrix ${\boldsymbol{\Theta}}$ that is not properly designed can
affect the performance of the CS systems. Hence, the optimal sensing matrix
must be constructed by also taking the SRE into consideration. Huang et.al. in
[9] found the optimal sensing matrix such that equivalent dictionary is close
to a target frame and also by taking the SRE into account:

$\begin{array}[]{ll}\underset{{\boldsymbol{\Theta}},{\boldsymbol{X}}_{\rm{ED-
Target}}}{\rm minimize}\>\>\omega\|{\boldsymbol{X}}_{{\rm{ED-
Target}}}-{\boldsymbol{\Theta}}{\boldsymbol{\Psi}}\|_{F}^{2}+(1-\omega)\|{\boldsymbol{\Theta}}({\boldsymbol{U}}-{\boldsymbol{\Psi}}{\boldsymbol{S}})\|_{F}^{2}\end{array}$
(8)

where ${\boldsymbol{X}}_{\rm{ED-Target}}$ is a certain target equivalent
dictionary and $\omega$ is the trade-off factor which takes the value from
$[0,1]$.

###  I-B Contributions and outline

In this paper, we propose an algorithm named TEchnique to devise Large
dimensional Equiangular Tight-frames (TELET-frames) based on Majorization
Minimization (MM) procedure to construct ETFs. We also use TELET to construct
optimal sensing matrix for the CS systems. Since our method does not involve
computing EVD and requires only computing some matrix vector product, it is
scalable to construct medium to large ETFs. As the proposed algorithm is based
on MM, it enjoys nice properties such as monotonicity and guaranteed
convergence to a stationary point of the problem. The major contributions of
the paper are as follows:

  1. 1.

A MM based iterative algorithm TELET is proposed to construct the ETFs. The
computational complexity of TELET is dictated only by matrix vector products
and does not involve computing EVD or matrix inversion at any iteration. We
also show a computational efficient way of implementing TELET. Subsequently,
we show that the frames constructed using TELET can be used to construct
optimal sensing matrix for the CS system using the model proposed in [9].

  2. 2.

We also propose a way to accelerate the convergence of the TELET algorithm.
The monotonic convergence of the TELET algorithm to the stationary point of
the frame design problem is proved.

  3. 3.

We show through numerical simulations that when compared to the state-of-the-
art algorithms, TELET can construct frames (especially of large dimensions, of
order $N=1000$) with low mutual coherence but with higher run time. Various
numerical simulations using synthetic data and real images are also conducted
to compare the performance of the sensing matrix constructed via TELET with
the other state-of-the-art algorithms used in the CS systems.

The paper is organized as follows. In Sec. II, we formulate the ETF design
problem and briefly discuss the MM procedure. Next, in Sec. III we develop the
proposed algorithm TELET using the MM procedure. In Sec. IV, we compare the
proposed algorithm with the state-of-the-art algorithm through numerical
simulations and conclude the paper in Sec. V.

##  II ETF Problem Formulation and Preliminaries

A unit norm ETF is constructed by minimizing the maximum absolute inner
product between pairs of distinct vectors in a frame:

$\displaystyle\underset{{\boldsymbol{x}}_{1},\cdots,{\boldsymbol{x}}_{N}}{\rm
minimize}$ $\displaystyle\underset{i,j=1,2,\cdots N\>(i\neq j)}{\rm
maximize}\>\left|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}\right|$ (9)
$\displaystyle\text{subject to}\>\>\|{\boldsymbol{x}}_{i}\|_{2}=1\quad
i=1,2,\cdots N$

Suppose
${\boldsymbol{x}}=[{\boldsymbol{x}}_{1}^{T},{\boldsymbol{x}}_{2}^{T},\cdots{\boldsymbol{x}}_{N}^{T}]^{T}\in\mathbb{C}^{Nd\times
1}$ is formed by stacking $\\{{\boldsymbol{x}}_{i}\\}_{i=1}^{N}$ and
${\boldsymbol{S}}_{l}$ be a $d\times dN$ selection matrix defined as
${\boldsymbol{S}}_{l}=[{\boldsymbol{0}}_{d\times(l-1)d},{\mathbf{I}}_{d},{\boldsymbol{0}}_{d\times(N-l)d}]$,
then we can obtain ${\boldsymbol{x}}_{l}$ using the following relation:

$\begin{array}[]{ll}{\boldsymbol{x}}_{l}={\boldsymbol{S}}_{l}{\boldsymbol{x}},\>l=1,2,\cdots,N\end{array}$
(10)

Using the above relation, we rewrite the absolute inner product
$|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}|$ in (9) as:

$\begin{array}[]{ll}|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}|=|{\boldsymbol{x}}^{H}{\boldsymbol{S}}_{j}^{H}{\boldsymbol{S}}_{i}{\boldsymbol{x}}|=|{\boldsymbol{x}}^{H}{\mathbf{A}}_{ij}{\boldsymbol{x}}|\end{array}$
(11)

Then using (11), the problem in (9) becomes:

$\displaystyle\underset{\boldsymbol{\boldsymbol{\boldsymbol{x}}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\text{\text{maximum}}}$}$ $\displaystyle\bar{f}_{ij}({\boldsymbol{x}})$
(12) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}=1\quad i=1,2,\cdots
N$

where
$\bar{f}_{ij}({\boldsymbol{x}})=2\left(|{{{{\boldsymbol{x}}}}}^{H}{{{{{\mathbf{A}}}_{ij}}}}{{{{\boldsymbol{x}}}}}|^{2}\right)$.
Please note that we have squared the objective function (which will not change
the optimum) and scaled by two for our convenience. The problem in (12) can be
further rewritten as:

$\displaystyle\underset{\boldsymbol{\boldsymbol{{{\boldsymbol{Y}}}}},\boldsymbol{\boldsymbol{\boldsymbol{x}}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\text{\text{maximum}}}$}$
$\displaystyle\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}){\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}}\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}})$
(13) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}=1\quad i=1,2,\cdots
N$
$\displaystyle\boldsymbol{\boldsymbol{\boldsymbol{Y}}}=\boldsymbol{\boldsymbol{x}}\boldsymbol{\boldsymbol{x}}^{H},$

where
${\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}}=\text{vec}\left({{{\mathbf{A}}}_{ij}}\right)\text{vec}^{H}\left(\left({{{\mathbf{A}}}_{ij}}\right)^{H}\right)+\text{vec}\left(\left({{{\mathbf{A}}}_{ij}}\right)^{H}\right)\text{vec}^{H}\left({{{\mathbf{A}}}_{ij}}\right)$.
We solve the above non-convex minimax problem using the MM procedure which
will be explained in the next subsection.

###  II-A MM for the minimax problem

Consider the following optimization problem:

$\begin{array}[]{ll}\underset{{\boldsymbol{x}}\in\chi}{\rm
minimize}\>f({\boldsymbol{x}})\end{array}$ (14)

where the objective function $f({\boldsymbol{x}})$ is equal to the maximum
value of the $i^{th}$ sub-function $\bar{f}_{i}({\boldsymbol{x}})$ i.e.
$f({\boldsymbol{x}})=\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{f}_{i}({\boldsymbol{x}})$ . The MM based algorithms solves the
above problem by first constructing a surrogate function
$g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ which majorizes the objective
function $f({\boldsymbol{x}})$ at the current iterate ${\boldsymbol{x}}^{t}$
and in the next step, the surrogate function is minimized to get the next
iterate.

${\boldsymbol{x}}^{t+1}\in\underset{{\boldsymbol{x}}\in\chi}{\rm
arg\>min}\>g\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)$ (15)

Majorization combines the tangency and the upper bound condition, i.e., a
function qualifies to be a surrogate function if it satisfies the following
conditions:

$g\left({\boldsymbol{x}}^{t}|{\boldsymbol{x}}^{t}\right)=f\left({\boldsymbol{x}}^{t}\right)$
(16) $g\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)\geq
f\left({\boldsymbol{x}}\right)$ (17)

Constructing $g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ for the minimax
problem in (14) does not look obvious, but can be worked out as follows, let:

$\begin{array}[]{ll}g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})=\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{g}_{i}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})\end{array}$ (18)

where each $\bar{g}_{i}({\boldsymbol{x}})$ is a tight upper bound on
$\bar{f}_{i}({\boldsymbol{x}})$ satisfying the following conditions:

$\begin{array}[]{ll}\bar{g}_{i}({\boldsymbol{x}}^{t}|{\boldsymbol{x}}^{t})=\bar{f}_{i}({\boldsymbol{x}}^{t})\end{array}$
(19)
$\begin{array}[]{ll}\bar{g}_{i}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})\geq\bar{f}_{i}({\boldsymbol{x}})\end{array}$
(20)

We now show that the surrogate function
${g}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ defined in (18) satisfies
condition (16) and (17):
_Proof of $g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ satisfying condition
(16)_:

$\begin{array}[]{ll}g({\boldsymbol{x}}^{t}|{\boldsymbol{x}}^{t})=\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{g}_{i}({\boldsymbol{x}}^{t}|{\boldsymbol{x}}^{t})=\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{f}_{i}({\boldsymbol{x}}^{t})=f({\boldsymbol{x}}^{t})\end{array}$
(21)

where the first and second equality is by (18) and (19) and the last equality
is by the definition of $f({\boldsymbol{x}})$.
_Proof of $g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ satisfying condition
(17)_:

$\begin{array}[]{ll}g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})=\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{g}_{i}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})\geq\underset{i=1,2,\cdots,I}{\rm
maximize}\>\bar{f}_{i}({\boldsymbol{x}}^{t})=f({\boldsymbol{x}}^{t})\end{array}$
(22)

where the first equality is given by (18), second inequality is by (20) and
the last equality is given by the definition of $f({\boldsymbol{x}})$. Using
(15), (16) and (17) it can be shown that the sequence of points
$\\{{\boldsymbol{x}}^{t}\\}$ generated by the MM procedure monotonically
decrease the objective function. Summary of different techniques used to
construct the surrogate function can be found in ([26], [27]). Further, since
the proposed algorithm is developed using MM, it can be shown using Theorem 1
in ([28]) that TELET will converge to a stationary point of the problem in
(9). Due to lack of space, we do not include the proof of convergence here.

##  III Algorithm for constructing large dimensional equiangular tight frames

In this section, we develop TEchnique to devise Large dimensional Equiangular
Tight-frames construction algorithm (TELET) using the MM procedure. For better
readability, we present our MM based algorithm in two stages: majorization
function construction and minimization. We also propose a way to accelerate
the convergence of TELET algorithm using SQUAREM acceleration scheme.

###  III-A Majorization Function construction

In this subsection, we construct the surrogate function for the objective
function in (13):

$\displaystyle\underset{\boldsymbol{\boldsymbol{{{\boldsymbol{Y}}}}},\boldsymbol{\boldsymbol{\boldsymbol{x}}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\text{\text{maximum}}}$}$ $\displaystyle\bar{f}_{ij}({\boldsymbol{Y}})$
(23) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}=1\quad i=1,2,\cdots
N$
$\displaystyle\boldsymbol{\boldsymbol{\boldsymbol{Y}}}=\boldsymbol{\boldsymbol{x}}\boldsymbol{\boldsymbol{x}}^{H},$

where
$\bar{f}_{ij}({\boldsymbol{Y}})=\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}){\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}}\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}})$.
As discussed in Sec. II.A, constructing majorizing function for
$\text{$\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\text{\text{maximum}}}$}\bar{f}_{ij}({\boldsymbol{Y}})$ boils down to
finding a majorizing function one for each $\bar{f}_{ij}({\boldsymbol{Y}})$.
Thus, we can focus on $\bar{f}_{ij}({\boldsymbol{Y}})$ only. Next, we discuss
the following lemmas which will be used to find a tight upper bound for
$\bar{f}_{ij}({\boldsymbol{Y}})$.

###### Lemma III.1

Let $f({\boldsymbol{x}})\in\mathbb{C}^{md\times 1}\rightarrow\mathbb{R}$ be a
continuous twice differentiable function in ${\boldsymbol{x}}$ and if
$f({\boldsymbol{x}})$ has a bounded curvature, then there exists a matrix
${\boldsymbol{M}}$ such that
$\boldsymbol{\boldsymbol{M}}\succeq\nabla^{2}f(\boldsymbol{\boldsymbol{\boldsymbol{x}}})$
and $f({\boldsymbol{x}})$ can be upper bounded as:

$\begin{array}[]{ll}f({\boldsymbol{x}})\leq
f({\boldsymbol{x}}^{t})+{\rm{Re}}\left(\nabla
f({\boldsymbol{x}}^{t})^{H}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})\right)+\dfrac{1}{2}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})^{H}{\boldsymbol{M}}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})\end{array}$
(24)

where ${\boldsymbol{x}}={\boldsymbol{x}}^{t}$ is the value taken by
${\boldsymbol{x}}$ at the $t^{th}$ iteration. The upper bound for
$f({\boldsymbol{x}})$ is quadratic and differentiable in ${\boldsymbol{x}}$.

###### Proof:

Suppose there exists a matrix ${\boldsymbol{M}}$ such that
$\boldsymbol{\boldsymbol{M}}\succeq\nabla^{2}f(\boldsymbol{\boldsymbol{\boldsymbol{x}}})$,
then we have the following inequality by second order Taylor expansion:

$\begin{array}[]{ll}f({\boldsymbol{x}})\leq
f({\boldsymbol{x}}^{t})+{\rm{Re}}\left(\nabla
f({\boldsymbol{x}}^{t})^{H}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})\right)+\dfrac{1}{2}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})^{H}{\boldsymbol{M}}({\boldsymbol{x}}-{\boldsymbol{x}}^{t})\end{array}$
(25)

and equality is achieved at ${\boldsymbol{x}}={\boldsymbol{x}}^{t}$. ∎

###### Lemma III.2

Given any ${{\boldsymbol{x}}={\boldsymbol{x}}^{t}}$, a concave function
$f\left({\boldsymbol{x}}\right)\in\mathbb{C}^{md\times
1}\rightarrow\mathbb{R}$ can be upper bounded as:

$\begin{array}[]{ll}{f}\left({\boldsymbol{x}}\right)\leq{f}\left({\boldsymbol{x}}^{t}\right)+{\rm{Re}}\left(\left(\nabla
f\left({\boldsymbol{x}}^{t}\right)\right)^{H}\left({\boldsymbol{x}}-{\boldsymbol{x}}^{t}\right)\right)\end{array}$
(26)

The upper bound for $f\left({\boldsymbol{x}}\right)$ is linear in
${\boldsymbol{x}}$.

###### Proof:

Since $f({\boldsymbol{x}})$ is concave, linearizing it around
${{\boldsymbol{x}}^{t}}$ using the first order Taylor series gives the above
inequality. ∎

Since $\bar{f}_{ij}({\boldsymbol{Y}})$ in (23) is quadratic and differentiable
in ${\textrm{vec}}({\boldsymbol{Y}})$, we apply Lemma III.1 with
${\boldsymbol{\boldsymbol{M}}}_{ij}=\lambda_{\text{max}}({\boldsymbol{\boldsymbol{\Phi}}}_{ij})\boldsymbol{\boldsymbol{I}}_{N^{2}d^{2}}=d{\boldsymbol{I}}_{N^{2}d^{2}}$
(refer Appendix 1 for proof) and use the relation
$\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}})\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}})=N^{2}$
to get the following surrogate function:

$\begin{array}[]{ll}\tilde{g}_{ij}\left({\boldsymbol{Y}}|{\boldsymbol{Y}}^{t}\right)=-\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}^{t})\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}^{t})+2N^{2}d+2\text{Re}(\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}^{t})\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}))-2d\text{Re}(\text{vec}^{H}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}}^{t})\text{vec}(\boldsymbol{\boldsymbol{\boldsymbol{Y}}})).\end{array}$
(27)

Next, we rewrite
$\tilde{g}_{ij}\left({\boldsymbol{Y}}|{\boldsymbol{Y}}^{t}\right)$ in terms of
${\boldsymbol{x}}$ by using
${\boldsymbol{Y}}={\boldsymbol{x}}{\boldsymbol{x}}^{H}$:

$\displaystyle\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)=-2|(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}{{{\mathbf{A}}}_{ij}}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}|^{2}+2(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}{{\boldsymbol{B}}}^{t}_{ij}\boldsymbol{\boldsymbol{\boldsymbol{x}}})-2d(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}})+2N^{2}d$
(28)

where
$\boldsymbol{B}^{t}_{ij}=\boldsymbol{\boldsymbol{\boldsymbol{A}}}_{ij}((\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}\boldsymbol{\boldsymbol{\boldsymbol{A}}}_{ij}^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})+\boldsymbol{\boldsymbol{\boldsymbol{A}}}_{ij}^{H}((\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}\boldsymbol{\boldsymbol{\boldsymbol{A}}}_{ij}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})$.
Using (11), ${{\boldsymbol{B}}}^{t}_{ij}$ can be written compactly as
${{\boldsymbol{B}}}^{t}_{ij}={\mathbf{A}}_{ij}\left(\left({\boldsymbol{x}}_{j}^{t}\right)^{H}{\boldsymbol{x}}_{i}^{t}\right)+{\mathbf{A}}_{ij}^{H}\left(\left({\boldsymbol{x}}_{i}^{t}\right)^{H}{\boldsymbol{x}}_{j}^{t}\right)$.
Note that the above surrogate function
$\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)$ has a
dominant convex quadratic term in ${\boldsymbol{x}}$:
$2(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}{{\boldsymbol{B}}}^{t}_{ij}\boldsymbol{\boldsymbol{\boldsymbol{x}}})$
which makes the resulting surrogate minimization problem intractable under the
unit norm constraint. To overcome this, we first add and subtract
$\lambda_{\text{max}}({{\boldsymbol{B}}}^{t}_{ij})$ to
$\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)$ and obtain
the surrogate function:

$\begin{array}[]{ll}\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)=-2|(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}{{{\mathbf{A}}}_{ij}}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}|^{2}+2(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}\hat{{\boldsymbol{B}}}^{t}_{ij}\boldsymbol{\boldsymbol{\boldsymbol{x}}})+2\lambda_{\text{max}}({{\boldsymbol{B}}}^{t}_{ij})N-2d(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}})+2N^{2}d\end{array}$
(29)

where
$\hat{{\boldsymbol{B}}}^{t}_{ij}={{\boldsymbol{B}}}^{t}_{ij}-(\lambda_{\text{max}}({{\boldsymbol{B}}}^{t}_{ij})\boldsymbol{\boldsymbol{I}}_{Nd})={{\boldsymbol{B}}}^{t}_{ij}-|\left({{\boldsymbol{x}}_{i}^{t}}\right)^{H}{{{\boldsymbol{x}}}_{j}^{t}}|\boldsymbol{\boldsymbol{I}}_{Nd}$
(see Appendix 2 for proof of
$\lambda_{\text{max}}({{\boldsymbol{B}}}^{t}_{ij})=|\left({{\boldsymbol{x}}_{i}^{t}}\right)^{H}{{{\boldsymbol{x}}}_{j}^{t}}|$).
Next, since the surrogate function
$\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)$ is now
concave in ${\boldsymbol{x}}$, we majorize it once again using Lemma III.2:

$\begin{array}[]{ll}\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})=-2|(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}{{{\mathbf{A}}}_{ij}}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}|^{2}+2\left(-(\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}\hat{{\boldsymbol{B}}}^{t}_{ij}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}+2\text{Re}\left((\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t})^{H}{\hat{{\boldsymbol{B}}}^{t}_{ij}}\boldsymbol{\boldsymbol{\boldsymbol{x}}}\right)\right)+2\lambda_{\text{max}}({{\boldsymbol{B}}}^{t}_{ij})N-\\\
2d\left(-N^{2}+2\text{Re}\left(N\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{H}\boldsymbol{\boldsymbol{\boldsymbol{x}}}^{t}\right)\right)+2N^{2}d.\vspace{-3.5mm}\end{array}$
(30)

The above surrogate function
$\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ is linear in
${\boldsymbol{x}}$ and hence can now be minimized. Note that although we do
double majorization i.e. find surrogate of the surrogate function, since
$\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ is a tighter surrogate
for $\tilde{g}_{ij}\left({\boldsymbol{x}}|{\boldsymbol{x}}^{t}\right)$, it can
be viewed as a direct surrogate to the objective function
$\bar{f}_{ij}({\boldsymbol{x}})$ in (12). The surrogate function in (30) can
be rewritten compactly as:

$\displaystyle\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})=4\text{Re}(\boldsymbol{\boldsymbol{x}}{}^{H}{{{\boldsymbol{d}}^{t}_{ij}}})+s^{t}_{ij}$
(31)

where
$s^{t}_{ij}=-6|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j}|^{2}+4N|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j}|+4N^{2}d$
and
${\boldsymbol{d}}_{ij}^{t}=\left({\mathbf{A}}_{ij}{\boldsymbol{x}}^{t}\right)(({\boldsymbol{x}}^{t}_{j})^{H}{\boldsymbol{x}}^{t}_{i})+\left({\mathbf{A}}_{ij}^{H}{\boldsymbol{x}}^{t}\right)(({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})-\left(|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})|{\mathbf{I}}_{Nd}\right){\boldsymbol{x}}^{t}-Nd{\boldsymbol{x}}^{t}$.
The superscript $t$ in $s^{t}_{ij}$ and ${\boldsymbol{d}}_{ij}^{t}$ denotes
that the computation of these quantities depends on the previous iterate value
${\boldsymbol{x}}^{t}$. The only computational complexity involved in
computing $\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ is the matrix
vector product ${\mathbf{A}}_{ij}{\boldsymbol{x}}^{t}$ and
${\mathbf{A}}_{ij}^{H}{\boldsymbol{x}}^{t}$. Exploiting the sparse structure
of ${\mathbf{A}}_{ij}$, the matrix vector product can be carried out
efficiently as:
${\mathbf{A}}_{ij}{\boldsymbol{x}}^{t}=[{\boldsymbol{0}}^{T}_{d(j-1)\times
1},({{\boldsymbol{x}}_{i}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-j)\times
1}]^{T}$ and
${\mathbf{A}}_{ij}^{H}{\boldsymbol{x}}^{t}=[{\boldsymbol{0}}^{T}_{d(i-1)\times
1},({{\boldsymbol{x}}_{j}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-i)\times
1}]^{T}$ Then the majorizing function for $\text{$\underset{i,j=1,2,\cdots
N\>(i\neq j)}{\text{\text{maximum}}}$}\bar{f}_{ij}({\boldsymbol{x}})$ is:

$\begin{array}[]{ll}g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})=\text{$\underset{i,j=1,2,\cdots
N\>(i\neq
j)}{\text{\text{maximum}}}$}\bar{g}_{ij}({\boldsymbol{x}}|{\boldsymbol{x}}^{t})\end{array}$
(32)

Hence at any iteration, given ${\boldsymbol{x}}^{t}$, the surrogate
minimization problem becomes:

$\displaystyle\underset{\boldsymbol{\boldsymbol{\boldsymbol{x}}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\text{\text{maximum}}}$}$ $\displaystyle
4\text{Re}(\boldsymbol{\boldsymbol{x}}{}^{H}{{{\boldsymbol{d}}^{t}_{ij}}})+s^{t}_{ij}$
(33) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}=1,\,i=1,2,\cdots,N$

###  III-B Solution to the Minimization Problem

The surrogate minimization problem in (33) is non-convex in ${\boldsymbol{x}}$
because of the equality constraint. However, the equality constraint can be
relaxed with the inequality constraint as the optimal solution of the relaxed
problem would lie only on the boundary of the relaxed constraint set [29]. The
relaxed problem is:

$\displaystyle\underset{{\boldsymbol{x}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{i,j=1,2,\cdots,N(i\neq
j)}{\text{\text{maximum}}}$}$ $\displaystyle
4\text{Re}(\boldsymbol{\boldsymbol{x}}{}^{H}\boldsymbol{\boldsymbol{d}}^{t}_{ij})+s^{t}_{ij}$
(34) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}\leq
1,\,i=1,2,\cdots,N\ $

Next, we rewrite the inner maximization problem over the indices as
maximization over a simplex variable:

$\displaystyle\underset{{\boldsymbol{x}}}{\text{\text{minimize}}}$
$\displaystyle\text{$\underset{{\boldsymbol{q}}\geq
0,{\boldsymbol{1}}^{T}{\boldsymbol{q}}=1}{\text{\text{maximum}}}$}$
$\displaystyle
4\text{Re}\left(\left({{\boldsymbol{D}}}^{t}{\boldsymbol{q}}\right)^{H}{\boldsymbol{x}}\right)+{\boldsymbol{q}}^{T}{\boldsymbol{s}}^{t}$
(35) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}\leq
1,\,i=1,2,\cdots,N\ $

where ${\boldsymbol{s}}^{t}=[s^{t}_{12},\,s^{t}_{13}\cdots
s^{t}_{N(N-1)}]^{T}\in\mathbb{R}^{N(N-1)\times 1}$ and
${{\boldsymbol{D}}}^{t}=[{{\boldsymbol{d}}}^{t}_{12},{{\boldsymbol{d}}}^{t}_{13},\cdots,{{\boldsymbol{d}}}^{t}_{N(N-1)}]\in\mathbb{R}^{2Nd\times
N(N-1)}$. The objective function in the above problem is bilinear in
${\boldsymbol{q}}$ and ${\boldsymbol{x}}$ and the constraints are compact
convex sets, then by minimax theorem [30], we can swap minimax to maximin
without altering the solution:

$\displaystyle\text{$\underset{{\boldsymbol{q}}\geq
0,{\boldsymbol{1}}^{T}{\boldsymbol{q}}=1}{\text{\text{maximum}}}$}$
$\displaystyle\underset{{\boldsymbol{x}}}{\text{\text{minimize}}}$
$\displaystyle
4\text{Re}\left(\left({{\boldsymbol{D}}}^{t}{\boldsymbol{q}}\right)^{H}{\boldsymbol{x}}\right)+{\boldsymbol{q}}^{T}{\boldsymbol{s}}^{t}$
(36) subject to $\displaystyle\|{\boldsymbol{x}}_{i}\|_{2}\leq
1,\,i=1,2,\cdots,N\ $

which is equivalent to

$\displaystyle\text{$\underset{{\boldsymbol{q}}\geq
0,{\boldsymbol{1}}^{T}{\boldsymbol{q}}=1}{\text{\text{maximize}}}$}$
$\displaystyle h({\boldsymbol{q}})$ (37)

where

$\begin{array}[]{ll}h({\boldsymbol{q}})=\underset{{\boldsymbol{x}}}{{\textrm{minimize}}}\>4\text{Re}\left(\left({{\boldsymbol{D}}}^{t}{\boldsymbol{q}}\right)^{H}{\boldsymbol{x}}\right)+{\boldsymbol{q}}^{T}{\boldsymbol{s}}^{t}\\\
\>\>\>\>\>\>\>\>\>\>\>\>\>{\textrm{subject
to}}\>\|{\boldsymbol{x}}_{i}\|_{2}\leq 1,\,i=1,2,\cdots,N\end{array}$ (38)

The above problem has a closed-form solution and is given as

${\boldsymbol{x}}^{t+1}_{i}=-\dfrac{{{{\boldsymbol{a}}}^{t}_{i}}}{\|{{{\boldsymbol{a}}}^{t}_{i}}\|_{2}},\>i=1,2,\cdots,N$
(39)

where ${\boldsymbol{a}}^{t}$ is made of blocks
${\boldsymbol{a}}^{t}_{1},{\boldsymbol{a}}^{t}_{2},\cdots,{\boldsymbol{a}}^{t}_{N}$
and ${\boldsymbol{a}}^{t}={{\boldsymbol{D}}}^{t}{\boldsymbol{q}}$ where
${\boldsymbol{q}}$ is obtained by solving the problem in (37). However, the
problem in (37) does not have a closed-form solution and hence we solve it
iteratively using Mirror Descent Algorithm (MDA) whose update step is given
as:

$\displaystyle{\boldsymbol{q}}^{k+1}={\underset{{\boldsymbol{q}}\geq
0,{\boldsymbol{1}}^{T}{\boldsymbol{q}}=1}{\text{\text{arg\>max}}}}$
$\displaystyle{({\boldsymbol{h}}^{k})}^{T}{\boldsymbol{q}}+\dfrac{1}{\gamma_{k}}{\boldsymbol{B}}_{\psi}\left({\boldsymbol{q}},{\boldsymbol{q}}^{k}\right)$
(40)

where ${\boldsymbol{q}}^{k}$ represents the value taken by ${\boldsymbol{q}}$
at the $k^{th}$ iteration of the MDA algorithm (this should not be confused
with superscript $t$ which is used as the iteration index in the outer loop),
${\boldsymbol{h}}^{k}\in\partial h({\boldsymbol{q}}^{k})$ is the subgradient
of $h({\boldsymbol{q}})$, $\gamma_{k}>0$ is the step size given by
$\gamma_{k}=\dfrac{\eta}{\sqrt{k}}$ where $\eta$ is some constant and
${\boldsymbol{B}}_{\psi}\left({\boldsymbol{q}},{\boldsymbol{q}}^{k}\right)$ is
Bregman-like distance generated by $\psi$ and is given by
${\boldsymbol{B}}_{\psi}\left({\boldsymbol{q}},{\boldsymbol{q}}^{k}\right)=\psi({\boldsymbol{q}})-\psi({\boldsymbol{q}}^{k})-\nabla^{T}\psi\left({\boldsymbol{q}}^{k}\right)\left({\boldsymbol{q}}-{\boldsymbol{q}}^{k}\right)$.
Since the constraint in the maximization problem is a unit simplex, similar to
[31], we choose $\psi({\boldsymbol{q}})$ as:

$\begin{array}[]{ll}\psi({\boldsymbol{q}})=\left\\{\begin{array}[]{ll}\displaystyle\sum_{i=1}^{N(N-1)}q_{i}\textrm{log}\>q_{i}&{\boldsymbol{q}}\in\mathcal{Q}\\\
+\infty&\rm{otherwise}\\\ \end{array}\right..\end{array}$ (41)

where
$\mathcal{Q}=\\{{\boldsymbol{q}}\in\mathbf{R}^{N(N-1)}|{\boldsymbol{1}}^{T}{\boldsymbol{q}}=1,{\boldsymbol{q}}\geq
0\\}$. Hence, the update step in (40) is simplified to:

$\begin{array}[]{ll}q_{j}^{k+1}=\dfrac{q_{j}^{k}\textrm{exp}\left(-\gamma_{k}h_{j}^{k}\right)}{\displaystyle\sum_{j=1}^{N(N-1)}q_{j}^{k}\textrm{exp}\left(-\gamma_{k}h_{j}^{k}\right)}\>j=1,2,\cdots
N(N-1)\end{array}$ (42)

The subgradient ${\boldsymbol{h}}^{k}$ is given as:

$\begin{array}[]{ll}{\boldsymbol{h}}^{k}=4\text{Re}\left(({{\boldsymbol{D}}}^{t})^{H}{\mathbf{y}}^{k}\right)+{{\boldsymbol{s}}}^{t}\end{array}$
(43)

where

$\begin{array}[]{ll}{\mathbf{y}}^{k}=\underset{{\mathbf{y}}}{\rm{min}}\>4\text{Re}\left(\left(\left({{\boldsymbol{D}}}\right)^{t}{\boldsymbol{q}}^{k}\right)^{H}{\mathbf{y}}\right)\\\
\hskip 14.22636pt\>{\textrm{subject to}}\>\|{\mathbf{y}}_{i}\|_{2}\leq
1,\,i=1,2,\cdots,N\end{array}$ (44)

The above problem has closed-form solution and is given by
${\mathbf{y}}^{k}_{i}=-\dfrac{{{{\boldsymbol{a}}}^{k,t}_{i}}}{\|{{{\boldsymbol{a}}}^{k,t}_{i}}\|_{2}}$
for $i=1,2\cdots N$ and
${\boldsymbol{a}}^{k,t}={{\boldsymbol{D}}}^{t}{\boldsymbol{q}}^{k}$ (we have
used the superscript $k,t$ for the vector ${\boldsymbol{a}}$ here since it
depends on the $k^{th}$ iteration value of ${\boldsymbol{q}}$ of the MDA
algorithm and also on the $t^{th}$ iteration value of ${\boldsymbol{D}}$). The
update step in (42) is repeated until convergence and once the optimal
${\boldsymbol{q}}$ is found, ${\boldsymbol{x}}$ can be calculated using (39).
The pseudocode of the proposed algorithm is shown in Algorithm 1.

Algorithm 1: Pseudocode of the proposed algorithm
---
Input: Number of vectors in frame $N$ and the dimension of each vector $d$.
Initialize: Set $t=0$. Initialize ${{{\boldsymbol{x}}}^{0}}$
Repeat:
1) Compute the following for $i,j=1,2,\cdots,N$ and $i\neq j$
a) Compute ${\boldsymbol{d}}^{t}_{ij}=[{\boldsymbol{0}}^{T}_{d(j-1)\times
1},({{\boldsymbol{x}}_{i}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-j)\times
1}]^{T}(({\boldsymbol{x}}^{t}_{j})^{H}{\boldsymbol{x}}^{t}_{i})+[{\boldsymbol{0}}^{T}_{d(i-1)\times
1},({{\boldsymbol{x}}_{j}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-i)\times
1}]^{T}(({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})-(|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})|+Nd){\boldsymbol{x}}^{t}$.
b) Compute
$s^{t}_{ij}=-6|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j}|^{2}+4N|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j}|+4N^{2}d$.
2) Compute
${{\boldsymbol{D}}}^{t}=[{{\boldsymbol{d}}^{t}_{12}},{{\boldsymbol{d}}^{t}_{13}},\cdots,{{\boldsymbol{d}}^{t}}_{(N(N-1))}]$
and ${\boldsymbol{a}}^{t}={{\boldsymbol{D}}^{t}}{\boldsymbol{q}}$ where
${\boldsymbol{q}}$ is obtained by solving the maximization problem in (37)
using the MDA algorithm as discussed in Subsection. III-B.
3) Compute
${\boldsymbol{x}}^{t+1}_{i}=-\dfrac{{{{\boldsymbol{a}}}^{t}_{i}}}{\|{{{\boldsymbol{a}}}^{t}_{i}}\|_{2}}$
for $i=1,2,\cdots,N$.
$t\leftarrow t+1$, until convergence

Before we end this subsection, we will here discuss the computational
complexity of TELET. As shown in the algorithmic development, TELET consists
of an inner loop which implements MDA to solve the surrogate minimization
problem and an outer loop which updates each element of
${\boldsymbol{d}}^{t}_{ij}$ and $s^{t}_{ij}$. We first discuss the complexity
of MDA - which is mostly dictated by the matrix vector product
$({{\boldsymbol{D}}}^{t})^{H}{\mathbf{y}}^{k}$ and
$({{\boldsymbol{D}}}^{t}){\boldsymbol{q}}^{k}$. The columns of the matrix
${{\boldsymbol{D}}^{t}}$ denoted as ${{\boldsymbol{d}}^{t}_{ij}}$ has a
special structure which we utilize to carry out this matrix vector operation
efficiently. Note that each ${\boldsymbol{d}}^{t}_{ij}$ is made by adding a
sparse vector: $[{\boldsymbol{0}}^{T}_{d(j-1)\times
1},({{\boldsymbol{x}}_{i}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-j)\times
1}]^{T}(({\boldsymbol{x}}^{t}_{j})^{H}{\boldsymbol{x}}^{t}_{i})+[{\boldsymbol{0}}^{T}_{d(i-1)\times
1},({{\boldsymbol{x}}_{j}^{t}})^{T},{\boldsymbol{0}}^{T}_{d(N-i)\times
1}]^{T}(({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})$\- which
consists of $2d$ non-zero elements and a non-sparse vector:
$(|({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j})|+Nd){\boldsymbol{x}}^{t}$
\- which is nothing but a scalar vector product. Exploiting this structure,
the matrix vector products can be carried out efficiently with complexity of
$\mathcal{O}\left(N^{2}d\right)$ per iteration. In the case of outer loop, the
major computations are that of computing the inner products
$({\boldsymbol{x}}^{t}_{i})^{H}{\boldsymbol{x}}^{t}_{j}$ which requires only
$\mathcal{O}\left(Nd\right)$ computation per iteration. Hence, after
neglecting the lower terms, the computational complexity of TELET is equal to
$\mathcal{O}\left(N^{2}d\right)$.

###  III-C Proof of convergence

We denote the objective function of the problem in (9) as
$f_{{}_{ETF}}({\boldsymbol{x}})$. Then, from the derivation of the proposed
algorithm, we know that $f_{{}_{ETF}}({\boldsymbol{x}})$ is majorized by
$g({\boldsymbol{x}}|{\boldsymbol{x}}^{t})$ in (32) at ${\boldsymbol{x}}^{t}$
over the constraint $\|{\boldsymbol{x}}_{i}\|_{2}\leq 1$. Using (15), (21) and
(22) it can be shown that the sequence of points $\\{{\boldsymbol{x}}^{t}\\}$
generated by the MM procedure monotonically decrease the objective function:

$\begin{array}[]{ll}f_{{}_{ETF}}({\boldsymbol{x}}^{t+1})\leq
g({\boldsymbol{x}}^{t+1}|{\boldsymbol{x}}^{t})\leq
g({\boldsymbol{x}}^{t}|{\boldsymbol{x}}^{t})=f_{{}_{ETF}}({\boldsymbol{x}}^{t})\end{array}$
(45)

where the first inequality and the last equality are obtained by using (21)
and (22). The second inequality is by (15). Since
$f_{{}_{ETF}}({\boldsymbol{x}})$ in (9) is bounded below by zero and the
points $\\{{\boldsymbol{x}}^{t}\\}$ monotonically decrease the objective
function, it is guaranteed that the sequence
$\\{f_{{}_{ETF}}({\boldsymbol{x}}^{t})\\}$ will converge to a finite value.
In the following, we show that the sequence $\\{{\boldsymbol{x}}^{t}\\}$
converges to the stationary point of the problem in (9). Assume that there is
a subsequence $\tilde{{\boldsymbol{x}}}^{r_{j}}$ converging to a limit point
${{\mathbf{p}}}$. Then from (21), (22) and (45) we get:

$\begin{array}[]{ll}g(\tilde{{\boldsymbol{x}}}^{r_{j+1}}|\tilde{{\boldsymbol{x}}}^{r_{j+1}})=f_{{}_{\textrm{ETF}}}({\tilde{{\boldsymbol{x}}}^{r_{j+1}}})\leq
f_{{}_{\textrm{ETF}}}({\tilde{{\boldsymbol{x}}}^{r_{j}+1}})\\\ \leq
g(\tilde{{\boldsymbol{x}}}^{r_{j}+1}|\tilde{{\boldsymbol{x}}}^{r_{j}})\leq
g(\tilde{{\boldsymbol{x}}}|\tilde{{\boldsymbol{x}}}^{r_{j}})\end{array}$ (46)

where $g(.)$ is the surrogate function as defined in (32). Then, letting
$j\rightarrow\infty$, we get:

$\begin{array}[]{ll}g({{\mathbf{p}}}|{{\mathbf{p}}})\leq
g(\tilde{{\boldsymbol{x}}}|{{\mathbf{p}}})\end{array}$ (47)

which implies $g^{\prime}({{\mathbf{p}}}|{{\mathbf{p}}})\geq 0$. Since the
first order behavior of surrogate function is same as function
$f_{{}_{\textrm{ETF}}}({\tilde{{\boldsymbol{x}}}})$ ([28]),
$g^{\prime}({{\mathbf{p}}}|{{\mathbf{p}}})\geq 0$ implies
$f^{\prime}_{{}_{\textrm{ETF}}}({{{\mathbf{p}}}})\geq 0$. Hence,
${{\mathbf{p}}}$ is the stationary point of
$f_{{}_{\textrm{ETF}}}({\tilde{{\boldsymbol{x}}}})$ and therefore the proposed
algorithm converges to the stationary point of the ETF problem.

###  III-D SQUAREM acceleration scheme

As discussed in Sec. II, the choice of surrogate function will dictate the
convergence speed of the MM based algorithm. In our case, since the surrogate
function is constructed through double majorization of the objective function,
it may lead to slower convergence of the proposed algorithm. Squared Iterative
Method (SQUAREM) [32] is based on Cauchy-Brazilai-Brownein method, to
accelerate the convergence rate. Let $F_{\rm{{MM}}}(\cdot)$ denote the
nonlinear fixed-point iteration map of the proposed algorithm. Then the MM
updating scheme can be expressed as:

${\boldsymbol{x}}^{t+1}=F_{\rm{{MM}}}({\boldsymbol{x}}^{t})$ (48)

The pseudocode of the proposed algorithm accelerated using SQUAREM
acceleration scheme is shown in Algorithm. 2.

Algorithm 2: Proposed algorithm acceleration scheme via SQUAREM
---
Input: Number of vectors in frame $N$ and the dimension of each vector $d$.
Initialize: Set $t=0$. Initialize ${{{\boldsymbol{x}}}^{0}}$
Repeat:
1) $\bar{{\boldsymbol{x}}}^{1}=F_{MM}({\boldsymbol{x}}^{t})$ ,
$\bar{{\boldsymbol{x}}}^{2}=F_{MM}(\bar{{\boldsymbol{x}}}^{1})$
2) ${\mathbf{r}}=\bar{{\boldsymbol{x}}}^{1}-{\boldsymbol{x}}^{t}$,
${\boldsymbol{v}}=\bar{{\boldsymbol{x}}}^{2}-\bar{{\boldsymbol{x}}}^{1}-{\mathbf{r}}$
3) Compute step length
$\alpha=-\dfrac{\|{\mathbf{r}}\|_{2}}{\|{\boldsymbol{v}}\|_{2}}$
4) Compute
${\boldsymbol{x}}^{t+1}={\boldsymbol{x}}^{t}-2\alpha{\mathbf{r}}+\alpha^{2}{\boldsymbol{v}}$
and normalize the frame.
5) Modify the step length $\alpha$ if $\underset{i,j=1,2,\cdots N\>(i\neq
j)}{\rm maximize}\>\bar{f}_{ij}({\boldsymbol{x}}^{t})>\underset{i,j=1,2,\cdots
N\>(i\neq j)}{\rm maximize}\>\bar{f}_{ij}({\boldsymbol{x}}^{t+1})$ as shown in
section 6 in [32]. Update
${\boldsymbol{x}}^{t+1}={\boldsymbol{x}}^{t}-2\alpha{\mathbf{r}}+\alpha^{2}{\boldsymbol{v}}$
and normalize the frame.
$t\leftarrow t+1$, until convergence

##  IV Numerical Simulations and CS Application

In this section we first compare TELET (whose convergence speed was
accelerated using the scheme discussed in the Subsection. III.D) with the
state-of-the-art algorithms used to construct the ETFs. When compared to the
other state-of-the-art algorithm, we show that TELET is monotonic and is able
to construct new ETFs of medium to large dimensions. Next, using synthetic
data and real images we compare the performance of optimal sensing matrix
constructed using TELET with the state-of-the-art algorithms used in the CS
systems. All the algorithms were implemented in MATLAB. The simulations were
carried out in a PC with $2.40$ GHz processor with $16$ GB RAM.

###  IV-A Comparison of TELET with the state-of-the-art algorithms

1\. In this simulation, we compare the convergence plot of the proposed
algorithm and that of the state-of-the-art algorithms for real and complex
ETFs. In the case of complex ETFs, we compare our algorithm with [8, 9, 10]
and [13, 12, 16] and in the case of real ETFs, we compare our algorithm with
the algorithms in [8, 9, 10] and [12, 13]. We denote the frames created by our
method and others as ${\textrm{Frame}}_{\textrm{TELET}}$,
${\textrm{Frame}}_{\textrm{Tropp}}$ [8], ${\textrm{Frame}}_{\textrm{Xiong}}$
[9], ${\textrm{Frame}}_{\textrm{Katsaggelos}}$ [10],
${\textrm{Frame}}_{\textrm{CBGC}}$[13],
${\textrm{Frame}}_{\textrm{ICBP}}$111The code for ICBP is obtained from
https://www.nt.tuwien.ac.at/christian-doppler-laboratory/cd-download/ [12] and
${\textrm{Frame}}_{\textrm{BCASC}}$222The codes for CBGC and BCASC are
obtained from https://codeocean.com/capsule/3494920/tree/v1 [16]. The
algorithm were initialized using the following technique: Generate a large
collection of random vectors ($>>N$) ${\boldsymbol{x}}_{i}$’s whose $j^{th}$
element was generated using:

$\begin{array}[]{ll}x_{j}={\rm{exp}}(j2\pi\phi_{j})\end{array}$ (49)

where $\phi_{j}$ is randomly generated from a uniform distribution from
$[0,1]$. Remove the vectors having the largest inner product. Repeat this
process until the desired number of vectors $N$ are reached. All the
algorithms were made to run until the following criteria was met:

$|\mu(\textbf{X}^{k})-\mu_{\textrm{CB}}|<10^{-5}$ (50)

where $\mu_{\textrm{CB}}$ is the composite bound (which is a tighter bound
when compared to the Welch bound in the region $N>d^{2}$), or until certain
maximum number of iterations was met. The composite bound in the case of
complex frames is defined as:

$\mu_{\textrm{CB}}=\begin{cases}\text{if }N\leq d^{2}:&\\\
\sqrt{\dfrac{N-d}{d(N-1)}},&\\\ \text{if }d^{2}<N\leq 2(d^{2}-1):&\\\
{\rm{max}}\left(\sqrt{\dfrac{1}{d}},\sqrt{\dfrac{2N-d^{2}-d}{(d+1)(N-d)}},1-2N^{-\dfrac{1}{d-1}}\right)&\\\
\text{if }N>2(d^{2}-1):&\\\
{\rm{max}}\left(\sqrt{\dfrac{2N-d^{2}-d}{(d+1)(N-d)}},1-2N^{-\dfrac{1}{d-1}}\right)&\end{cases}$

In the case of real frames, composite bound is defined as:

$\mu_{\textrm{CB}}={\rm{max}}\left(\sqrt{\dfrac{N-d}{d(N-1)}},\sqrt{\dfrac{3N-d^{2}-2d}{(d+2)(N-d)}}\right)$

Fig. 1.a and Fig. 1.b compares the convergence plot of the frame design
algorithms for complex and real ETF, respectively. From Fig. 1.a it can be
seen that for some frame dimensions such as $3\times 6$ and $4\times 7$, the
algorithms $\textrm{Frame}_{\textrm{Kastaggelos}}$ and
$\textrm{Frame}_{\textrm{Xiong}}$ converges to a higher mutual coherence value
when compared to the other state-of-the-art algorithms. Also, from Fig. 1.a,
it can be seen that the remaining algorithms converges closer to the composite
bound, except for frame dimensions $3\times 5$, $4\times 6$ and $5\times 7$,
where ${\textrm{Frame}}_{\textrm{Tropp}}$ algorithm converges to a larger
coherence value. Also, from Fig. 1.b it can be seen that almost always
${\textrm{Frame}}_{\textrm{ICBP}}$, ${\textrm{Frame}}_{\textrm{CBGC}}$ and the
proposed algorithm converge to similar coherence value for the frame
dimensions considered in Fig. 1.b.

(a) Convergence plot: ${\rm
max}\>\left|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}\right|$ vs. iteration
for complex ETFs of different dimensions.

(b) Convergence plot: ${\rm
max}\>\left|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}\right|$ vs. iteration
for real ETFs of different dimensions.

Figure 1: Convergence plot for different values of $d$ and $N$ for real and
complex ETF.

2\. We now compare the algorithms with respect to the following metrics:
mutual coherence and run time for frames of small, medium and large
dimensions. The maximum number of iterations was set equal to $10^{4}$ in the
case of small and medium frame dimensions and $10^{3}$ in the case of large
frame dimensions (due to larger computational complexity of the methods in
comparison we have reduced the number of iterations for large frame
dimensions). Table. I \- Table. IV compares the mutual coherence value and run
time of the algorithms for frames of various dimensions, respectively. In the
case of CBGC and BCASC we report the cumulative sum of run time which was
required to solve each sub-problem (problem (6) in [13]). For the other
algorithms, we report the run time either until the condition in (50) was met
or the run time corresponding to the minimum mutual coherence value achieved
by the algorithms. In Table. I \- Table. IV we have highlighted the algorithms
which performs the best with respect to mutual coherence value and run time,
respectively. From Table. I it can be seen that when compared to the other
algorithms ${\textrm{Frame}}_{\textrm{Tropp}}$ and
${\textrm{Frame}}_{\textrm{Xiong}}$ have higher mutual coherence value for
small frame dimensions. Hence, we did not include them in the comparison for
medium and large frame dimensions. From Table. III it can be seen that
${\textrm{Frame}}_{\textrm{CBGC}}$ and ${\textrm{Frame}}_{\textrm{BCASC}}$
takes more time to converge when compared to the other algorithms. This could
be because unlike the other algorithms, ${\textrm{Frame}}_{\textrm{CBGC}}$ and
${\textrm{Frame}}_{\textrm{BCASC}}$ requires solving sub-problems for various
values of $p$ to obtain frames with minimum coherence. Hence, we did not
include them for comparison for large frame dimensions. Also, from Table. I
and Table. II it can be observed that out of the $54$ frame dimensions
considered, the proposed algorithm achieves the lowest coherence value $42$
times. Especially, for large frame dimensions, it can be observed that the
proposed algorithm almost always achieves lower mutual coherence value when
compared to the state-of-the-art algorithm. From Table. IV it can be seen that
for all the medium and large frame dimensions considered,
${\textrm{Frame}}_{\textrm{ICBP}}$ always takes lesser time when compared to
the proposed and state-of-the-art algorithms. In Fig. 2 we also show the
convergence plot for some of the large frame dimensions considered in Table.
II. From Fig. 2 it can be seen that for some of the frame dimensions such as
$25\times 800$ and $30\times 1200$, the proposed algorithm even though
achieves lower coherence value, has not converged and had to be stopped
because of the maximum iteration limit.

Table I: Comparison of Mutual coherence of complex frames of small and medium dimensions. $(d,N)$ |  $\mu_{\textrm{CB}}$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$ |  ${\textrm{Frame}}_{\textrm{CBGC}}$ |  ${\textrm{Frame}}_{\textrm{BCASC}}$ |  ${\textrm{Frame}}_{\textrm{Xiong}}$ |  ${\textrm{Frame}}_{\textrm{Tropp}}$
---|---|---|---|---|---|---|---
$(2,8)$ |  $0.75$ |  0.7941 |  0.7941 |  $0.7942$ |  $0.7950$ |  $0.8110$ |  $0.9270$
$(3,16)$ |  $0.6202$ |  0.6483 |  $0.6547$ |  $0.6491$ |  $0.6499$ |  $0.6870$ |  0.6970
(4,5) |  0.25 |  0.25 |  0.25 |  0.25 |  0.25 |  0.25 |  0.25
(4,6) |  0.3162 |  0.3273 |  0.3274 |  0.3276 |  0.3282 |  0.4425 |  0.3536
(4,7) |  0.3536 |  0.3536 |  0.3536 |  0.3536 |  0.3536 |  0.5 |  0.3536
(4,8) |  0.3780 |  0.3780 |  0.3780 |  0.3780 |  0.3780 |  0.5 |  0.3780
(4,9) |  0.3853 |  0.4019 |  0.4131 |  0.4022 |  0.4025 |  0.5 |  0.416
(4,10) |  0.4082 |  0.4110 |  0.4112 |  0.4118 |  0.4124 |  0.5 |  0.414
(4,19) |  0.5 |  0.5 |  0.5188 |  0.5 |  0.5 |  0.5330 |  0.5313
(4,20) |  0.5 |  0.5272 |  0.5371 |  0.5278 |  0.5264 |  0.5398 |  0.5932
(5,6) |  0.2 |  0.2 |  0.2 |  0.2 |  0.2 |  0.2 |  0.2
(5,7) |  0.2582 |  0.2664 |  0.2665 |  0.2667 |  0.2676 |  0.3658 |  0.3154
(5,8) |  0.2928 |  0.2953 |  0.2952 |  0.2954 |  0.2957 |  0.4472 |  0.3055
(5,9) |  0.3162 |  0.3203 |  0.3216 |  0.3205 |  0.3212 |  0.4472 |  0.3331
(5,10) |  0.3333 |  0.3333 |  0.3333 |  0.3333 |  0.3333 |  0.4472 |  0.3334
(5,16) |  0.3830 |  0.3892 |  0.3941 |  0.3895 |  0.3905 |  0.4472 |  0.4041
(5,26) |  0.4472 |  0.4486 |  0.4511 |  0.4495 |  0.4479 |  0.4564 |  0.4671
(6,37) |  0.4082 |  0.4180 |  0.4176 |  0.4147 |  0.4164 |  0.4258 |  0.4357
(20,30) |  0.1313 |  0.1315 |  0.1317 |  0.1314 |  0.1317 |  - |  -
(20,50) |  0.1750 |  0.1753 |  0.1753 |  0.1751 |  0.1754 |  - |  -
(20,100) |  0.2010 |  0.2152 |  0.2118 |  0.2109 |  - |  - |  -
(30,40) |  0.0925 |  0.0930 |  0.0930 |  0.0927 |  - |  - |  -
(30,50) |  0.1166 |  0.1168 |  0.1169 |  0.1167 |  - |  - |  -
(30,100) |  0.1535 |  0.1606 |  0.1553 |  0.1549 |  - |  - |  -
Table II: Comparison of Mutual coherence of large dimension complex frames constructed by different algorithms $(d,N)$ |  $\mu_{\textrm{CB}}$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$
---|---|---|---
(23,500) |  0.2039 |  0.2178 |  0.3257
(23,600) |  0.2163 |  0.3195 |  0.3401
(23,700) |  0.2285 |  0.3218 |  0.3632
(23,800) |  0.2371 |  0.3334 |  0.3721
(23,900) |  0.2435 |  0.3431 |  0.3895
(23,1000) |  0.2485 |  0.3511 |  0.3758
(23,1200) |  0.2558 |  0.3580 |  0.3978
(25,500) |  0.1951 |  0.2807 |  0.3047
(25,600) |  0.1960 |  0.2069 |  0.3391
(25,700) |  0.2067 |  0.3032 |  0.3265
(25,800) |  0.2171 |  0.3143 |  0.3632
(25,900) |  0.2248 |  0.3249 |  0.3619
(25,1000) |  0.2308 |  0.3292 |  0.3702
(25,1200) |  0.2393 |  0.3424 |  0.3428
(27,500) |  0.1874 |  0.2718 |  0.2913
(27,600) |  0.1882 |  0.2791 |  0.3011
(27,700) |  0.1888 |  0.2909 |  0.3327
(27,800) |  0.1975 |  0.2999 |  0.3388
(27,900) |  0.2067 |  0.3086 |  0.3358
(27,1000) |  0.2137 |  0.3144 |  0.3517
(27,1200) |  0.2237 |  0.3307 |  0.3632
(30,800) |  0.1792 |  0.2823 |  0.2443
(30,900) |  0.1796 |  0.3153 |  0.3376
(30,1000) |  0.1886 |  0.2962 |  0.2650
(30,1200) |  0.2013 |  0.3337 |  0.3561
(40,800) |  0.1542 |  0.2466 |  0.2480
(40,900) |  0.1546 |  0.2454 |  0.2596
(40,1000) |  0.1550 |  0.1619 |  0.2737
(40,1200) |  0.1555 |  0.2745 |  0.2988
(50,1000) |  0.1379 |  0.2229 |  0.2368
Table III: Comparison of run time (seconds) of the algorithms to construct complex frames of small and medium dimensions. $(d,N)$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$ |  ${\textrm{Frame}}_{\textrm{CBGC}}$ |  ${\textrm{Frame}}_{\textrm{BCASC}}$ |  ${\textrm{Frame}}_{\textrm{Xiong}}$ |  ${\textrm{Frame}}_{\textrm{Tropp}}$
---|---|---|---|---|---|---
$(2,8)$ |  $0.8408$ |  $0.16$ |  $9.68$ |  $5749.9$ |  $0.034$ |  0.014
$(3,16)$ |  13.4313 |  3.58 |  496.23 |  7017.9 |  0.049 |  0.061
(4,5) |  0.2082 |  0.0734 |  0.0080 |  0.1288 |  0.0075 |  0.0066
(4,6) |  40.8189 |  0.1344 |  6.6730 |  58.18 |  0.0025 |  0.2492
(4,7) |  1.6661 |  0.3182 |  0.054 |  1.3887 |  0.1556 |  0.0053
(4,8) |  2.8156 |  0.9512 |  11.754 |  180.116 |  0.37 |  0.23
(4,9) |  7.1192 |  1.7749 |  30.1607 |  161.78 |  0.189 |  0.3126
(4,10) |  16.1572 |  1.54 |  16.2422 |  146.28 |  0.4609 |  0.4701
(4,19) |  82.1215 |  4.3984 |  103.5152 |  2909.2 |  0.1584 |  3.5493
(4,20) |  224.2072 |  4.3790 |  223.69 |  5871 |  0.2831 |  0.0425
(5,6) |  0.3619 |  0.0839 |  0.0315 |  0.4355 |  0.0047 |  0.0040
(5,7) |  0.2017 |  0.1359 |  16.8251 |  2354.6 |  0.0265 |  0.5062
(5,8) |  21.5216 |  1.1905 |  21.3538 |  3226.7 |  0.4841 |  0.5052
(5,9) |  20.6065 |  1.5550 |  118720 |  81271 |  0.2590 |  0.4304
(5,10) |  17.6087 |  0.7675 |  48.6249 |  10.1563 |  0.3976 |  0.4967
(5,16) |  158.11 |  3.7181 |  396.65 |  3390.4 |  1.0742 |  0.0345
(5,26) |  124.627 |  6.2891 |  1927.3 |  9932.3 |  6.3617 |  0.2744
(6,37) |  548.18 |  1.7595 |  1760.1 |  31121 |  7.708 |  1.5603
(20,30) |  1024.7 |  36.54 |  5443.2 |  122950 |  - |  -
(20,50) |  1083.8 |  61.1703 |  6321.6 |  343680 |  - |  -
(20,100) |  7246.7 |  99.9807 |  35599 |  - |  - |  -
(30,40) |  1749.9 |  53.6401 |  9260.3 |  - |  - |  -
(30,50) |  2778.5 |  116.5102 |  3928 |  - |  - |  -
(30,100) |  6397.8 |  286.835 |  43210 |  - |  - |  -
Table IV: Comparison of run time (seconds) of the proposed and state-of-the-art algorithms to construct large dimensional complex frames. $(d,N)$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$
---|---|---
(23,500) |  9627.5 |  1528
(23,600) |  7053.6 |  1268.1
(23,700) |  15235 |  3498.4
(23,800) |  15307 |  2668.4
(23,900) |  22147 |  3269.4
(23,1000) |  42103 |  9255.4
(23,1200) |  48966 |  7868.4
(25,500) |  6778 |  1597.1
(25,600) |  13701 |  2832.3
(25,700) |  10573 |  2014.5
(25,800) |  20653 |  7385
(25,900) |  28069 |  9803.2
(25,1000) |  45089 |  7142.4
(25,1200) |  86220 |  11260
(27,500) |  9658.9 |  1581.8
(27,600) |  1390.7 |  6952
(27,700) |  24756 |  3900.3
(27,800) |  14244 |  2777.8
(27,900) |  28770 |  10211
(27,1000) |  37806 |  13187
(27,1200) |  80170 |  19889
(30,800) |  31733 |  5093.1
(30,900) |  47413 |  7525.5
(30,1000) |  43566 |  7661.9
(30,1200) |  89700 |  11891
(40,800) |  16350 |  3732
(40,900) |  35033 |  7978.6
(40,1000) |  31342 |  5591.0
(40,1200) |  55046 |  11676
(50,1000) |  46612 |  9186.1

|  |  |  |  |
---|---|---|---|---|---
$23\times 700$ | $23\times 800$ | $23\times 900$ | $25\times 700$ | $25\times 800$ | $25\times 900$

|  |  |  |  |
---|---|---|---|---|---
$27\times 700$ | $27\times 800$ | $27\times 900$ | $30\times 900$ | $30\times 1000$ | $30\times 1200$

Figure 2: Convergence plot: ${\rm
max}\>\left|\textbf{x}_{i}^{H}\textbf{x}_{j}\right|$ vs. iteration of complex
frames of large dimensions. Red line - converge of
${\textrm{Frame}}_{\textrm{TELET}}$; Magenta line - convergence of
${\textrm{Frame}}_{\textrm{ICBP}}$ and Blue line - composite bound.

3\. We now repeat the above experiments for real frames. All the algorithms
were initialized with the same random vectors with the initialization scheme
as discussed in the previous simulation. Similar to the previous simulation,
the maximum number of iterations was set equal to $10^{4}$ in the case of
small and medium frame dimensions and $10^{3}$ in the case of large frame
dimensions. Table. V \- Table. VIII compares the mutual coherence and run time
of the algorithms for frames of different dimensions, respectively. In Table.
V \- Table. VIII we have highlighted the algorithms which performs the best
with respect to mutual coherence value and run time. Similar to the
observations made for complex case, from Table. V it can seen that
${\textrm{Frame}}_{\textrm{Xiong}}$ and ${\textrm{Frame}}_{\textrm{Tropp}}$
have higher mutual coherence for the considered small frame dimensions and
from Table. VII it can be seen that CBGC algorithm takes more time to
converge. Also, from Table. VI it can be seen that the proposed algorithm
always achieves the lowest coherence value for the large frame dimensions when
compared to the other state-of-the-art algorithms. Moreover, from Table. V and
Table. VI it can be seen that out of the $54$ frames considered, the proposed
algorithm achieves the least mutual coherence value $43$ times. Also, from
Table. VII and Table. VIII, it can be seen that
${\textrm{Frame}}_{\textrm{ICBP}}$ takes lesser time for medium and large
frame dimensions. In Fig. 3 we also show the convergence plot of some of the
large frame dimensions considered in Table. VI. Similar to the observation
made for complex case in Fig. 2, it can be seen in Fig. 3 that for some of the
frame dimensions (such as $23\times 1000$ and $25\times 1000$), the proposed
algorithm even though achieves lower coherence value, has not converged and is
stopped due to the maximum iteration limit.
We now conclude the results of the simulations conducted for various
dimensions of complex and real frames. As can be observed from the
simulations, in the case of complex frames, the two best algorithms (among the
competing algorithms) are the proposed algorithm
${\textrm{Frame}}_{\textrm{TELET}}$ and ${\textrm{Frame}}_{\textrm{ICBP}}$
algorithm. The mutual coherence value attained by the two best algorithms are
comparable for small and medium dimension frames. However, with the increase
in the value of the frame dimensions, it can be seen that the proposed
algorithm almost always achieves lower mutual coherence value when compared to
${\textrm{Frame}}_{\textrm{ICBP}}$ algorithm. We also calculated the
percentage decrease in the mutual coherence value achieved by the proposed
algorithm with respect to coherence value achieved by the
${\textrm{Frame}}_{\textrm{ICBP}}$ using the formula:

${\textrm{Mutual coherence decrease
in}}\>\%=\left(\dfrac{\mu_{\textrm{ICBP}}-\mu_{\textrm{TELET}}}{\mu_{\textrm{ICBP}}}\right)\times
100$ (51)

where $\mu_{\textrm{ICBP}}$ and $\mu_{\textrm{TELET}}$ are the mutual
coherence values obtained by ${\textrm{Frame}}_{\textrm{ICBP}}$ and
${\textrm{Frame}}_{\textrm{TELET}}$ algorithms, respectively for the large
frame dimensions mentioned in Table. II. We found that the average percentage
decrease in the mutual coherence value achieved by the proposed algorithm with
respect to ${\textrm{Frame}}_{\textrm{ICBP}}$ to be equal to $11.58\%$. Also,
for some large frame dimensions, it was noticed that while the
${\textrm{Frame}}_{\textrm{ICBP}}$ algorithm had converged, the proposed
algorithm had not converged and still achieved lower coherence value. However,
for large dimensional frames, the average run time of
${\textrm{Frame}}_{\textrm{ICBP}}$ algorithm was found to be about $5$ times
computationally faster than the proposed algorithm. Similarly, in the case of
real frames, for small and medium dimensions, the mutual coherence value
attained by ${\textrm{Frame}}_{\textrm{TELET}}$ and
${\textrm{Frame}}_{\textrm{ICBP}}$ algorithm were comparable. However, with
the increase in the value of the frame dimensions, the
${\textrm{Frame}}_{\textrm{TELET}}$ algorithm always achieving lower coherence
value when compared to ${\textrm{Frame}}_{\textrm{ICBP}}$. For large
dimensions, the average percentage decrease in the mutual coherence value
achieved by the proposed algorithm with respect to
${\textrm{Frame}}_{\textrm{ICBP}}$ was found to be equal to $33.12\%$. Also,
similar to complex case, for some of the large dimensional frames, it was
found that while ${\textrm{Frame}}_{\textrm{ICBP}}$ algorithm had converged to
a larger coherence value, the proposed algorithm had not converged and yet
achieved a lower coherence value. However, with respect to average run time,
${\textrm{Frame}}_{\textrm{ICBP}}$ was found to be faster than
${\textrm{Frame}}_{\textrm{TELET}}$ for large frames.

Table V: Comparison of Mutual coherence of real frames constructed by different algorithms for small and medium dimensions $(d,N)$ |  $\mu_{\textrm{CB}}$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$ |  ${\textrm{Frame}}_{\textrm{CBGC}}$ |  ${\textrm{Frame}}_{\textrm{Xiong}}$ |  ${\textrm{Frame}}_{\textrm{Tropp}}$
---|---|---|---|---|---|---
$(2,8)$ |  $0.8165$ |  0.9239 |  0.9239 |  0.9239 |  $0.99$ |  $0.9602$
$(3,16)$ |  0.7125 |  0.7947 |  0.7947 |  0.7947 |  0.9799 |  0.9778
(4,5) |  0.25 |  0.25 |  0.25 |  0.25 |  0.25 |  0.25
(4,6) |  0.3162 |  0.3333 |  0.3334 |  0.3338 |  0.4998 |  0.3536
(4,7) |  0.3536 |  0.3909 |  0.3904 |  0.3905 |  0.5 |  0.4438
(4,8) |  0.3780 |  0.3780 |  0.3952 |  0.3780 |  0.5 |  0.3780
(4,9) |  0.3953 |  0.4343 |  0.4344 |  0.4344 |  0.5 |  0.5192
(4,10) |  0.4082 |  0.4343 |  0.4343 |  0.4344 |  0.5 |  0.4407
(4,19) |  0.6055 |  0.6467 |  0.6774 |  0.65529 |  0.7149 |  0.7377
(4,20) |  0.6124 |  0.6548 |  0.6546 |  0.6553 |  0.7114 |  0.9065
(5,6) |  0.2 |  0.2 |  0.2 |  0.2 |  0.2 |  0.2
(5,7) |  0.2582 |  0.2863 |  0.2862 |  0.2865 |  0.3960 |  0.3861
(5,8) |  0.2928 |  03291 |  0.3304 |  0.3338 |  0.4472 |  0.4015
(5,9) |  0.3162 |  0.3334 |  0.3334 |  0.3336 |  0.4472 |  0.4148
(5,10) |  0.3333 |  0.3333 |  0.3855 |  0.3333 |  0.4472 |  0.3333
(5,16) |  0.4109 |  0.4472 |  0.4472 |  0.4472 |  0.4472 |  0.4636
(5,26) |  0.5408 |  0.5950 |  0.6316 |  0.5910 |  0.7150 |  0.8141
(6,37) |  0.5040 |  0.5583 |  0.5345 |  0.5354 |  0.5764 |  0.8225
(20,30) |  0.1313 |  0.1406 |  0.141 |  0.1402 |  - |  -
(20,50) |  0.1750 |  0.1970 |  0.1949 |  0.1937 |  - |  -
(20,100) |  0.2010 |  0.2609 |  0.2465 |  0.2471 |  - |  -
(30,40) |  0.0925 |  0.0994 |  0.1006 |  0.0984 |  - |  -
(30,50) |  0.1166 |  0.1266 |  0.1258 |  0.1243 |  - |  -
(30,100) |  0.1535 |  0.1919 |  0.1768 |  0.1764 |  - |  -
Table VI: Comparison of Mutual coherence of large dimension real frames constructed by different algorithms $(d,N)$ |  $\mu_{\textrm{CB}}$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$
---|---|---|---
(23,500) |  0.2785 |  0.3703 |  0.5850
(23,600) |  0.2914 |  0.3869 |  0.6139
(23,700) |  0.3002 |  0.4006 |  0.6027
(23,800) |  0.3065 |  0.4449 |  0.6817
(23,900) |  0.3113 |  0.4220 |  0.6703
(23,1000) |  0.3151 |  0.4697 |  0.6785
(23,1200) |  0.3206 |  0.4476 |  0.6736
(25,500) |  0.2536 |  0.3468 |  0.4878
(25,600) |  0.2692 |  0.3641 |  0.5772
(25,700) |  0.2796 |  0.4028 |  0.6091
(25,800) |  0.2871 |  0.3900 |  0.6403
(25,900) |  0.2928 |  0.4003 |  0.6437
(25,1000) |  0.2972 |  0.4340 |  0.6799
(25,1200) |  0.3036 |  0.4486 |  0.6926
(27,500) |  0.2286 |  0.3275 |  0.4336
(27,600) |  0.2474 |  0.3440 |  0.5119
(27,700) |  0.2598 |  0.3597 |  0.5715
(27,800) |  0.2686 |  0.3704 |  0.6461
(27,900) |  0.2752 |  0.3816 |  0.6444
(27,1000) |  0.2803 |  0.3911 |  0.6480
(27,1200) |  0.2878 |  0.4070 |  0.6604
(30,800) |  0.2417 |  0.3458 |  0.5086
(30,900) |  0.25 |  0.3607 |  0.6008
(30,1000) |  0.2564 |  0.3655 |  0.6076
(30,1200) |  0.2655 |  0.3851 |  0.6549
(40,800) |  0.1542 |  0.2909 |  0.3357
(40,900) |  0.1680 |  0.3015 |  0.3773
(40,1000) |  0.1809 |  0.3119 |  0.3954
(40,1200) |  0.1985 |  0.3275 |  0.4488
(50,1000) |  0.1379 |  0.2788 |  0.3217
Table VII: Comparison of run time (seconds) of the proposed and state-of-the-art algorithms to construct real frames of small, medium and large dimensions $(d,N)$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$ |  ${\textrm{Frame}}_{\textrm{CBGC}}$ |  ${\textrm{Frame}}_{\textrm{Xiong}}$ |  ${\textrm{Frame}}_{\textrm{Tropp}}$
---|---|---|---|---|---
$(2,8)$ |  $0.1920$ |  $0.0287$ |  $1.9849$ |  $0.0009$ |  0.0004
(3,16) |  11.46 |  0.0915 |  15.4802 |  0.0159 |  0.059
(4,5) |  0.3121 |  0.0039 |  62 |  0.0031 |  0.00003
(4,6) |  0.5354 |  0.0288 |  0.1023 |  0.2349 |  0.3027
(4,7) |  6.8206 |  0.1115 |  2.1150 |  0.1972 |  0.3523
(4,8) |  0.5354 |  0.4750 |  11.9228 |  0.2812 |  0.2473
(4,9) |  1.43 |  0.0758 |  1.7694 |  0.2693 |  0.1309
(4,10) |  7.8786 |  0.079 |  1.1954 |  0.2651 |  0.176
(4,19) |  12.218 |  5.4175 |  818.0925 |  0.0324 |  0.66
(4,20) |  24.2077 |  0.4647 |  49.6893 |  0.0070 |  0.0025
(5,6) |  0.5059 |  0.0067 |  0.0078 |  0.0016 |  0.0012
(5,7) |  1.8895 |  0.32 |  0.7191 |  0.007 |  0.0003
(5,8) |  2.5685 |  0.7696 |  0.2503 |  0.2406 |  0.0122
(5,9) |  43.7984 |  0.2205 |  0.8046 |  0.4533 |  0.0005
(5,10) |  3.7073 |  1.3572 |  0.0341 |  0.2521 |  0.0025
(5,16) |  19.586 |  1.0214 |  0.8949 |  0.159 |  0.1139
(5,26) |  74.3653 |  8.9504 |  524.6116 |  0.0348 |  0.0098
(6,37) |  95.505 |  0.8462 |  47.4128 |  2.5269 |  0.0029
(20,30) |  127.97 |  6.128 |  206.79 |  - |  -
(20,50) |  56.78 |  11.7099 |  1700.7 |  - |  -
(20,100) |  1626.6 |  59.1218 |  12201 |  - |  -
(30,40) |  181.3425 |  4.4492 |  322.132 |  - |  -
(30,50) |  289.8671 |  10.9571 |  2221.2 |  - |  -
(30,100) |  2608 |  52.827 |  14339 |  - |  -
Table VIII: Comparison of run time (seconds) of the proposed and state-of-the-art algorithms to construct large dimension real frames. $(d,N)$ |  ${\textrm{Frame}}_{\textrm{TELET}}$ |  ${\textrm{Frame}}_{\textrm{ICBP}}$
---|---|---
(23,500) |  7447.6 |  936.7474
(23,600) |  12993 |  1871
(23,700) |  20205 |  2764.8
(23,800) |  18562 |  2523.5
(23,900) |  19512 |  1878.8
(23,1000) |  35523 |  4342.4
(23,1200) |  57621 |  4086.1
(25,500) |  9089.4 |  745.7
(25,600) |  16677 |  1420.4
(25,700) |  13243 |  1769.4
(25,800) |  15117 |  1202.5
(25,900) |  19670 |  1937.4
(25,1000) |  34711 |  4371.5
(25,1200) |  81627 |  7430.7
(27,500) |  7142.5 |  1029.8
(27,600) |  15357 |  1505
(27,700) |  22302 |  2247.6
(27,800) |  29570 |  4348.4
(27,900) |  25403 |  3374.8
(27,1000) |  33407 |  6421
(27,1200) |  54222 |  4083
(30,800) |  33546 |  2915.4
(30,900) |  13589 |  1816
(30,1000) |  22993 |  2778.5
(30,1200) |  57005 |  4112.9
(40,800) |  19930 |  2856.1
(40,900) |  46389 |  4809.8
(40,1000) |  34270 |  7634.8
(40,1200) |  77984 |  6932.5
(50,1000) |  26612 |  3332.7

|  |  |  |  |
---|---|---|---|---|---
$23\times 800$ | $23\times 900$ | $23\times 1000$ | $25\times 800$ | $25\times 900$ | $25\times 1000$

|  |  |  |  |
---|---|---|---|---|---
$27\times 800$ | $27\times 900$ | $27\times 1000$ | $30\times 800$ | $30\times 900$ | $30\times 1000$

Figure 3: Convergence plot: ${\rm
max}\>\left|\textbf{x}_{i}^{H}\textbf{x}_{j}\right|$ vs. iteration of real
frames of large dimensions.Red line - converge of
${\textrm{Frame}}_{\textrm{TELET}}$; Magenta line - convergence of
${\textrm{Frame}}_{\textrm{ICBP}}$ and Blue line - composite bound.

###  IV-B Compressed sensing application

As discussed in the introduction section, an optimal sensing matrix must be
constructed such that its equivalent dictionary
${\boldsymbol{X}}_{{\rm{ED}}}={\boldsymbol{\Theta}}{\boldsymbol{\Psi}}$ has
small mutual coherence and also by taking the sparse representation error into
account. Therefore, using the model proposed in [9], we solve the following
problem to construct an optimal sensing matrix:

$\begin{array}[]{ll}\underset{{\boldsymbol{\Theta}},{\boldsymbol{X}}_{\rm{ED-
Target}}}{\rm minimize}\>\>\omega\|{\boldsymbol{X}}_{\rm{ED-
Target}}-{\boldsymbol{\Theta}}{\boldsymbol{\Psi}}\|_{F}^{2}+(1-\omega)\|{\boldsymbol{\Theta}}{\boldsymbol{E}}\|_{F}^{2}\end{array}$
(52)

where $\omega$ is the trade-off factor which takes value in $[0,1]$,
${\boldsymbol{X}}_{\rm{ED-Target}}$ is a target equivalent dictionary with
properties of the ETF. The above problem is solved iteratively using the
alternating minimization scheme wherein minimization is performed with respect
to one matrix while keeping the other matrix fixed and vice-versa. At the
$t^{th}$ iteration, the sensing matrix ${\boldsymbol{\Theta}}^{t}$ is fixed
and the equivalent dictionary ${\boldsymbol{X}}_{\rm{ED-Target}}^{t+1}$ is
updated using the TELET algorithm shown in Algorithm 2. Next, using the
updated equivalent dictionary the problem in (52) can be solved - which has a
closed-form solution and can be found leveraging Theorem $2$ in [9].
We now conduct simulations to demonstrate the performance of the CS system
using the optimized sensing matrices constructed using the above described
approach and the state-of-the-art algorithms ([21, 25, 22, 9]). We denote the
CS system of the proposed algorithm as ${\rm{CS_{TELET}}}$ and the algorithms
in [21], [25], [22] and [9] as ${\rm{CS_{Elad}}}$, ${\rm{CS_{Sanei}}}$,
${\rm{CS_{Sapiro}}}$ and ${\rm{CS_{Xiong}}}$, respectively.
_1\. Experiments on synthetic data_
To evaluate the CS systems under various conditions, we first obtained a set
of lower dimensional measurement vectors $\\{{\mathbf{y}}_{t}\\}_{t=1}^{R}$
($R$ is the number of monte-carlo experiments) using the following
experimental procedure. We first generated $R$ number of vectors
$\\{{\boldsymbol{s}}_{t}\\}_{t=1}^{R}$ with $K$ non-zero entries of length
$N$. The non-zero elements of ${\boldsymbol{s}}_{t}$ were randomly positioned
and they were generated using Gaussian distribution with zero mean and unit
variance. Then, the sparse signal vector ${\boldsymbol{u}}^{*}_{t}$ was
generated as
${\boldsymbol{u}}^{*}_{t}={\boldsymbol{\Psi}}{\boldsymbol{s}}_{t}$, where the
dictionary ${\boldsymbol{\Psi}}$ of dimension $N\times N$ dictionary was
chosen as the Haar wavelet matrix. Next, we obtain the following contaminated
signals:

$\begin{array}[]{ll}{\boldsymbol{u}}_{t}={\boldsymbol{u}}^{*}_{t}+\epsilon_{t}={\boldsymbol{\Psi}}{\boldsymbol{s}}_{t}+\epsilon_{t},\>\>t=1,2\cdots
R\end{array}$ (53)

where $\\{\epsilon_{t}\\}_{t=1}^{R}$ represents the sparse representation
error and is assumed to Gaussian distributed noise with zero mean and variance
$\sigma^{2}$. The lower dimensional measurement vectors
$\\{{\mathbf{y}}_{t}\\}_{t=1}^{R}$ each of length $d$ were obtained using
${\mathbf{y}}_{t}={\boldsymbol{\Theta}}{\boldsymbol{u}}_{t}$ where
${\boldsymbol{\Theta}}$ is the optimized sensing matrix constructed using the
proposed approach ${\rm{CS_{TELET}}}$ and the other state-of-the-art
algorithms. From the obtained measurement vectors, the high dimensional signal
vector $\hat{{\boldsymbol{u}}}_{t}$ is reconstructed using
$\hat{{\boldsymbol{u}}}_{t}={\boldsymbol{\Psi}}\hat{{\boldsymbol{s}}}_{t}$
where $\hat{{\boldsymbol{s}}}_{t}$ is obtained by the Basis Pursuit algorithm
[20]. The reconstruction accuracy is quantified using the mean square error
(MSE) metric defined as:

$\begin{array}[]{ll}{\rm{MSE}(\hat{{\boldsymbol{U}}})}=\dfrac{1}{d\times
R}\|\hat{{\boldsymbol{U}}}-{\boldsymbol{U}}^{*}\|_{F}^{2}\end{array}$ (54)

where
${\boldsymbol{U}}^{*}=[{\boldsymbol{u}}_{1}^{*},{\boldsymbol{u}}_{2}^{*},\cdots,{\boldsymbol{u}}_{R}^{*}]$,
$\hat{{\boldsymbol{U}}}=[\hat{{\boldsymbol{u}}}_{1},\hat{{\boldsymbol{u}}}_{2},\cdots,\hat{{\boldsymbol{u}}}_{R}]$,
each $\hat{{\boldsymbol{u}}}_{t}$ is the reconstructed version of
${\boldsymbol{u}}_{t}^{*}$ and is equal to
$\hat{{\boldsymbol{u}}}_{t}={\boldsymbol{\Psi}}\hat{{\boldsymbol{s}}}_{t}$. We
evaluate the CS systems by varying the number of measurements $d$ and the
sparsity level $K$. In the first experiment we fixed the number of
measurements $d$ and varied sparsity level $K$ while in the second experiment
we varied sparsity level $K$ for fixed number of measurements $d$. The
parameters of the first experiment are as follows: Haar wavelet matrix of
dimension $N=32$ was chosen, the number of measurements $d=10$, the number of
monte-carlo experiments $R=50$, noise variance $\sigma^{2}=0.25$, $\omega=0.5$
[9] and varied the sparsity of the vector $K$ from $2$ to $7$. In the second
experiment, we fixed the signal sparsity $K=4$ and varied the number of
measurements $d$ from $10$ to $15$ with the same dimension of the dictionary
$N$, noise variance $\sigma^{2}$ and the number of monte-carlo experiments $R$
as that of the first experiment. Fig. 4.a and Fig. 4.b shows the MSE
($\hat{{\boldsymbol{U}}}$) vs. sparsity $K$ and MSE ($\hat{{\boldsymbol{U}}}$)
vs. measurement dimension $d$ of the proposed and the state-of-the-art
algorithms, respectively. From Fig. 4 it can be seen that the proposed
algorithm has lower MSE ($\hat{{\boldsymbol{U}}}$) for all the values of
sparsity $K$ and measurement dimension $d$.

(a) MSE ($\hat{{\boldsymbol{U}}}$) vs. Sparsity $K$ for Haar wavelet
dictionary

(b) MSE ($\hat{{\boldsymbol{U}}}$) vs. measurement dimension $d$ for Haar
wavelet dictionary

Figure 4:  MSE ($\hat{{\boldsymbol{U}}}$) vs. Sparsity $K$ and measurement
dimension $d$ for Haar wavelet dictionary using BP reconstruction algorithm

Next, we compare the two competing algorithms ${\rm{CS_{Xiong}}}$ and
${\rm{CS_{TELET}}}$ in the context of image reconstruction from lower
dimensional patches.
_2\. Experiments on image data_
Images from USC-SIPI database [33] were obtained to evaluate the performance
of ${\rm{CS_{Xiong}}}$ and ${\rm{CS_{TELET}}}$ in the context of image
reconstruction. This database is commonly used for conducting image processing
research and contains images in TIFF format. We choose three different images
of different sizes from the database - ‘House’ image of size $256\times 256$,
‘Boat’ image of size $512\times 512$ and a ‘Male’ image of size $1024\times
1024$. Each image of size $Q\times Q$ is first divided into non-overlapping
patches of size $8\times 8$ and each patch is vectorized to get
${\boldsymbol{u}}_{t}$’s of size $64\times 1$. Then from each patch, we
obtained the compressed lower dimensional vector ${\mathbf{y}}_{t}$’s of size
$d\times 1$ using the relation
${\mathbf{y}}_{t}={\boldsymbol{\Theta}}{\boldsymbol{u}}_{t}$ where
${\boldsymbol{\Theta}}$ is the optimized sensing matrix constructed using
${\rm{CS_{TELET}}}$ and ${\rm{CS_{Xiong}}}$ algorithm. For ‘House’ image we
obtained $31.25\>\%$ ($d=20$) compressed measurements from each patch. In the
case of ‘Boat’ image we took $35.9\>\%$ ($d=23$) compressed measurements from
each patch and finally for ‘Male’ image we obtained $40.62\>\%$ ($d=30$)
compressed measurements. From the obtained compressed measurements, the high
dimensional vector $\hat{{\boldsymbol{u}}}_{t}$ is reconstructed using
$\hat{{\boldsymbol{u}}}_{t}={\boldsymbol{\Psi}}\hat{{\boldsymbol{s}}}_{t}$
where we have chosen ${\boldsymbol{\Psi}}$ as the DCT matrix,
$\hat{{\boldsymbol{s}}}_{t}$ is obtained by using the basis pursuit algorithm.
Then, each reconstructed vector $\hat{{\boldsymbol{u}}}_{t}$ is reshaped to
get back a patch of size $p\times p$. Fig. 5 visually compares the original
and the reconstructed images using ${\rm{CS_{Xiong}}}$ and ${\rm{CS_{TELET}}}$
algorithm. From Fig. 5 it can be seen that ${\rm{CS_{TELET}}}$ algorithm is
able to reconstruct the last two test images with high fidelity. In the case
of ‘House’ image, even though the reconstructed image using
${\rm{CS_{TELET}}}$ looks better than the image reconstructed using
${\rm{CS_{Xiong}}}$, its fidelity is low. This could be because of the smaller
size of the ‘House’ image and also because of the higher compression from each
patch. The quantitative accuracy of the reconstructed images is evaluated
using PSNR metric as defined in [34]. image of the proposed algorithm was
found to be $25.910$ dB, $25.66$ dB and $28.244$ dB, respectively while that
of the algorithm in [9] was found to $24.457$ dB, $25.12$ dB and $27$ dB,
respectively. Hence, ${\rm{CS_{TELET}}}$ has higher PSNR for all the images,
with mean SNR of $26.604$ dB, when compared to ${\rm{CS_{Xiong}}}$ which has
mean SNR equal to $25.526$ dB.

|  |
---|---|---
1.(a) | 1\. (b) | 1\. (c)

|  |
---|---|---
2\. (a) | 2\. (b) | 2\. (c)

|  |
---|---|---
3\. (a) | 3\. (b) | 3\. (c)

Figure 5: (a) Original image (b) Reconstructed image by the proposed algorithm
(c) Reconstructed image by the algorithm in [9]

##  V Conclusion

In this paper, we proposed an iterative algorithm TELET based on the
Majorization Minimization principle to construct equiangular tight frames.
Since the proposed algorithm is developed using Majorization Minimization
procedure, it enjoys nice properties such as monotonicity and guaranteed
convergence to a stationary point of the ETF problem. We also apply TELET
algorithm to construct optimized sensing matrix for compressed sensing
systems. We show through numerical simulations that when compared to the
state-of-the-art algorithms TELET can construct complex and real frames with
very low mutual coherence value, especially for large frame dimensions.
However, with respect to run time, when compared to state-of-the-art
algorithm, it is found that TELET is computationally slower. We also compared
the performance of the optimized sensing matrix obtained using TELET algorithm
with the state-of-the-algorithms used in compressed sensing systems using both
synthetic and real images and found that the proposed algorithm performs
better in terms of reconstruction accuracy.

##  VI Appendix

1\. Proof of $\lambda_{\text{max}}({\boldsymbol{\boldsymbol{\Phi}}}_{ij})=d$

###### Proof:

We first recall that:

${\boldsymbol{\boldsymbol{\boldsymbol{\Phi}}}_{ij}}=\text{vec}\left({{{\mathbf{A}}}_{ij}}\right)\text{vec}^{H}\left(\left({{{\mathbf{A}}}_{ij}}\right)^{H}\right)+\text{vec}\left(\left({{{\mathbf{A}}}_{ij}}\right)^{H}\right)\text{vec}^{H}\left({{{\mathbf{A}}}_{ij}}\right)$
(55)

Notice that by construction of
${\mathbf{A}}_{ij}={\boldsymbol{S}}_{j}^{H}{\boldsymbol{S}}_{i}$, (where
${\boldsymbol{S}}_{l}=[{\boldsymbol{0}}_{d\times(l-1)},{\mathbf{I}}_{d},{\boldsymbol{0}}_{d\times(N-l)d}]$
is the selection matrix), it will have $d$ off-diagonal elements equal to one
and zeros everywhere else. Let
${\boldsymbol{U}}_{ij}=\text{vec}({\mathbf{A}}_{ij})\text{vec}^{H}({{\mathbf{A}}_{ij}^{H}})={\boldsymbol{q}}_{ij}{\boldsymbol{s}}_{ij}^{H}$.
Then it is easy to see that:

$\text{Tr}({\boldsymbol{U}}_{ij}{\boldsymbol{U}}_{ij}^{H})=\|{\boldsymbol{s}}_{ij}\|_{2}^{2}\|{\boldsymbol{q}}_{ij}\|_{2}^{2}=d^{2}.$
(56) $\text{Tr}({\boldsymbol{U}}_{ij})=0$ (57)

Since ${\boldsymbol{\Phi}}_{ij}$ is the sum of two rank one matrices, its rank
can be at most two. Let $\lambda_{1}$ and $\lambda_{2}$ be the two non-zero
eigen values of ${\boldsymbol{\Phi}}_{ij}$. Then using (56) and (57) we have:

$\begin{array}[]{ll}\lambda_{1}+\lambda_{2}=\rm{Tr}({\boldsymbol{U}}_{ij}+{\boldsymbol{U}}_{ij}^{H})=0\end{array}$
(58)
$\begin{array}[]{ll}\lambda_{1}^{2}+\lambda_{2}^{2}=\rm{Tr}({\boldsymbol{U}}_{ij}+{\boldsymbol{U}}_{ij}^{H})^{2}=2\rm{Tr}({{\boldsymbol{U}}_{ij}{\boldsymbol{U}}_{ij}^{H}})=2d^{2}\end{array}$
(59)

Using (58) and (59) we get:

$\begin{array}[]{ll}\lambda_{1}\lambda_{2}=\dfrac{1}{2}\left(\left(\lambda_{1}+\lambda_{2}\right)^{2}-\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)\right)=-d^{2}\end{array}$
(60)

Then, we can find the maximum eigen value by solving the following equation:

$\begin{array}[]{ll}\lambda^{2}-d^{2}=0\end{array}$ (61)

whose solution is equal to $\lambda=\pm d$, choosing its maximum gives
$\lambda_{max}({\boldsymbol{\Phi}}_{ij})=d$ ∎

2\. Proof of
$\lambda_{\text{max}}({{\boldsymbol{B}}}_{ij})=|c_{ij}|=|{\boldsymbol{x}}_{i}^{H}{\boldsymbol{x}}_{j}|$

###### Proof:

To prove the above, we first define $\bar{{\mathbf{A}}}_{ij}$ and
$\bar{{\boldsymbol{B}}}_{ij}$ matrices as follows:

$\begin{array}[]{cc}\bar{{\mathbf{A}}}_{ij}=\end{array}\left(\begin{array}[]{cc}0&c_{ij}\\\
c_{ij}^{*}&0\end{array}\right),\>\bar{{\boldsymbol{B}}}_{ij}=\left(\begin{array}[]{cc}{\boldsymbol{0}}_{d\times
d}&c_{ij}{\mathbf{I}}_{d\times d}\\\ c_{ij}^{*}{\mathbf{I}}_{d\times
d}&{\boldsymbol{0}}_{d\times d}\end{array}\right)$

One can easily see that the eigenvalues of $\bar{{\mathbf{A}}}_{ij}$ is equal
to $\lambda_{\bar{{\mathbf{A}}}_{ij}}=\pm|c_{ij}|$. Next, the eigenvalues of
$\bar{{\boldsymbol{B}}}_{ij}$ can be obtained by solving the following
equation (which we obtained by cofactor expansion):

${\lambda_{\bar{{\boldsymbol{B}}}_{ij}}}=(-{\lambda}^{2}+|c_{ij}|^{2})^{d-1}\rm{det}(\bar{{\mathbf{A}}}_{ij}-\lambda{\mathbf{I}})=0$
(62)

Hence, the maximum eigenvalue of $\bar{{\boldsymbol{B}}}_{ij}$ will be equal
to $|c_{ij}|$. Also, note that by inserting zero vector columns and rows
between the blocks of $\bar{{\boldsymbol{B}}}_{ij}$ wont change its maximum
eigenvalue. We use these observations to prove that
$\lambda_{\text{max}}({{\boldsymbol{B}}}_{ij})=|c_{ij}|$ where
${\boldsymbol{B}}_{ij}={\mathbf{A}}_{ij}(c_{ij}^{*})+{\mathbf{A}}_{ij}^{H}(c_{ij})$.
By construction, the matrix ${\boldsymbol{B}}_{ij}$ will be a block diagonal
matrix with a main block equal to $\bar{{\boldsymbol{B}}}_{ij}$ (with or
without the zero vectors inserted between its blocks) and the remaining main
blocks equal to zero square matrices. Sometimes, depending on the value of
$(i,j)$, the entire matrix ${\boldsymbol{B}}_{ij}$ will be equal to
$\bar{{\boldsymbol{B}}}_{ij}$ with zero vectors inserted between the blocks.
For the later case, as discussed earlier, the addition of zero vectors do not
change the maximum eigenvalue of $\bar{{\boldsymbol{B}}}_{ij}$, hence the
maximum eigenvalue of ${\boldsymbol{B}}_{ij}=|c_{ij}|$. For the former case,
we make use of the property that the eigenvalues of a block diagonal matrix
will be the eigenvalues of its blocks which implies the maximum eigenvalue of
${\boldsymbol{B}}_{ij}=|c_{ij}|$. ∎

## References

  * [1] G. Bachmann, L. Narici, and E. Beckenstein, “Fourier and wavelet analysis,” _Springer Science & Business Media_, 2012.
  * [2] T. Strohmer and R. Heath, “Grassmannian frames with applications to coding and communication,” _arXiv preprint math/0301135_ , 2003.
  * [3] P. Viswanath and V. Anantharam, “Optimal sequences and sum capacity of synchronous cdma systems,” _IEEE Transactions on information theory_ , vol. 45, no. 6, pp. 1984–1991, 1999.
  * [4] D. Donoho and M. Elad, “Maximal sparsity representation via $\ell_{1}$ minimization,” _Proceedings of National Academy of Sciences_ , vol. 100, pp. 2197–2202, 2003.
  * [5] M. Fickus and D. G. Mixon, “Numerically erasure-robust frames,” _arXiv preprint arXiv:1202.4525_ , 2012.
  * [6] Y. C. Eldar and G. D. Forney, “Optimal tight frames and quantum measurement,” _IEEE Transactions on Information Theory_ , vol. 48, no. 3, pp. 599–610, 2002\.
  * [7] J. A. Tropp, “Complex equiangular tight frames,” _International Society for Optics and Photonics_ , vol. 5914, p. 591401, 2005.
  * [8] J. A. Tropp, I. S. Dhillon, R. W. Heath, and T. Strohmer, “Designing structured tight frames via an alternating projection method,” _IEEE Transactions on information theory_ , vol. 51, no. 1, pp. 188–209, 2005.
  * [9] H. Bai, S. Li, and X. He, “Sensing matrix optimization based on equiangular tight frames with consideration of sparse representation error,” _IEEE Transactions on Multimedia_ , vol. 18, no. 10, pp. 2040–2053, 2016.
  * [10] E. V. Tsiligianni, L. P. Kondi, and A. K. Katsaggelos, “Construction of incoherent unit norm tight frames with application to compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 60, no. 4, pp. 2319–2330, 2014.
  * [11] K. Jaganathan and B. Hassibi, “Reconstruction of integers from pairwise distances,” _International Conference on Acoustics, Speech and Signal Processing_ , 2012.
  * [12] B. Tahir, S. Schwarz, and M. Rupp, “Constructing grassmannian frames by an iterative collision-based packing,” _IEEE Signal Processing Letters_ , vol. 26, no. 7, pp. 1056–1060, 2019.
  * [13] H. E. A. Laue and W. P. Du Plessis, “A coherence-based algorithm for optimizing rank-1 grassmannian codebooks,” _IEEE Signal Processing Letters_ , vol. 24, no. 6, pp. 823–827, 2017.
  * [14] C. Rusu and N. Gonzalez-Prelcic, “Designing incoherent frames through convex techniques for optimized compressed sensing,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 9, pp. 2334–2344, 2016.
  * [15] C. Rusu, “Design of incoherent frames via convex optimization,” _IEEE Signal Processing Letters_ , vol. 20, no. 7, pp. 673–676, 2013.
  * [16] H. Zörlein and M. Bossert, “Coherence optimization and best complex antipodal spherical codes,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 24, pp. 6606–6615, 2015.
  * [17] D. L. Donoho and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via $\ell_{1}$ minimization,” _Proceedings of the National Academy of Sciences_ , vol. 100, no. 5, pp. 2197–2202, 2003.
  * [18] J. A. Tropp, “Greed is good: Algorithmic results for sparse approximation,” _IEEE Transactions on Information theory_ , vol. 50, no. 10, pp. 2231–2242, 2004.
  * [19] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Transactions on information theory_ , vol. 53, no. 12, pp. 4655–4666, 2007.
  * [20] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” _SIAM review_ , vol. 43, no. 1, pp. 129–159, 2001.
  * [21] M. Elad, “Optimized projections for compressed sensing,” _IEEE Transactions on Signal Processing_ , vol. 55, no. 12, pp. 5695–5702, 2007.
  * [22] J. M. Duarte-Carvajalino and G. Sapiro, “Learning to sense sparse signals: Simultaneous sensing matrix and sparsifying dictionary optimization,” _IEEE Transactions on Image Processing_ , vol. 18, no. 7, pp. 1395–1408, 2009\.
  * [23] G. Li, Z. Zhu, D. Yang, L. Chang, and H. Bai, “On projection matrix optimization for compressive sensing systems,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 11, pp. 2887–2898, 2013.
  * [24] J. Xu, Y. Pi, and Z. Cao, “Optimized projection matrix for compressive sensing,” _EURASIP Journal on Advances in Signal Processing_ , vol. 2010, no. 1, p. 560349, 2010.
  * [25] V. Abolghasemi, S. Ferdowsi, and S. Sanei, “A gradient-based alternating minimization approach for optimization of the measurement matrix in compressive sensing,” _Signal Processing_ , vol. 92, no. 4, pp. 999–1009, 2012.
  * [26] D. R. Hunter and K. Lange, “A tutorial on mm algorithms,” _The American Statistician_ , vol. 58, no. 1, pp. 30–37, 2004.
  * [27] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 794–816, 2016.
  * [28] M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence analysis of block successive minimization methods for nonsmooth optimization,” _SIAM Journal on Optimization_ , vol. 23, no. 2, pp. 1126–1153, 2013.
  * [29] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [30] J. v. Neumann, “Zur theorie der gesellschaftsspiele,” _Mathematische annalen_ , vol. 100, no. 1, pp. 295–320, 1928.
  * [31] A. Beck and M. Teboulle, “Mirror descent and nonlinear projected subgradient methods for convex optimization,” _Operations Research Letters_ , vol. 31, no. 3, pp. 167–175, 2003.
  * [32] R. Varadhan and C. Roland, “Simple and globally convergent methods for accelerating the convergence of any em algorithm,” _Scandinavian Journal of Statistics_ , vol. 35, no. 2, pp. 335–353, 2008.
  * [33] U. SIPI, “The usc-sipi image database,” 2016. [Online]. Available: http://sipi.usc.edu/services/database/data-base.html
  * [34] R. Gonzalez and R. Woods, “Digital image processing: Pearson prentice hall,” _Upper Saddle River, NJ_ , vol. 1, pp. 376–376, 2008.