{'fraction_non_alphanumeric': 0.04587137100090468, 'fraction_numerical': 0.01574965046467637, 'mean_word_length': 3.875889724310777, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-Calibration and Biconvex Compressive Sensing

The design of high-precision sensing devises becomes ever more difficult and expensive. At the same time, the need for precise calibration of these devices (ranging from tiny sensors to space telescopes) manifests itself as a major roadblock in many scientific and technological endeavors. To achieve optimal performance of advanced high-performance sensors one must carefully calibrate them, which is often difficult or even impossible to do in practice. In this work we bring together three seemingly unrelated concepts, namely Self-Calibration, Compressive Sensing, and Biconvex Optimization. The idea behind self-calibration is to equip a hardware device with a smart algorithm that can compensate automatically for the lack of calibration. We show how several self-calibration problems can be treated efficiently within the framework of biconvex compressive sensing via a new method called SparseLift. More specifically, we consider a linear system of equations y = DAx, where both x and the diagonal matrix D (which models the calibration error) are unknown. By"lifting"this biconvex inverse problem we arrive at a convex optimization problem. By exploiting sparsity in the signal model, we derive explicit theoretical guarantees under which both x and D can be recovered exactly, robustly, and numerically efficiently via linear programming. Applications in array calibration and wireless communications are discussed and numerical simulations are presented, confirming and complementing our theoretical analysis.


Introduction
The design of high-precision sensing devises becomes ever more difficult and expensive. Meanwhile, the need for precise calibration of these devices (ranging from tiny sensors to space telescopes) manifests itself as a major roadblock in many scientific and technological endeavors. To achieve optimal performance of advanced high-performance sensors one must carefully calibrate them, which is often difficult or even impossible to do in practice. It is therefore necessary to build into the sensors and systems the capability of self-calibration -using the information collected by the system to simultaneously adjust the calibration parameters and perform the intended function of the system (image reconstruction, signal estimation, target detection, etc.). This is often a challenging problem requiring advanced mathematical solutions. Given the importance of this problem there has been a substantial amount of prior work in this area. Current self-calibration algorithms are based on joint estimation of the parameters of interest and the calibration parameters using standard techniques such as maximum likelihood estimation. These algorithms tend to be very computationally expensive and can therefore be used only in relatively simple situations. This has prevented their widespread applications in advanced sensors.
In this paper we take a step toward developing a new framework for self-calibration by bringing together three seemingly unrelated concepts, namely Self-Calibration, Compressive Sensing, and Biconvex Optimization. By "lifting" the problem and exploiting sparsity in the signal model and/or the calibration model we express the biconvex problem of self-calibration as a convex optimization problem, which can be solved efficiently via linear programming. This new method, called SparseLift, ensures equivalence between the biconvex problem and the convex problem at the price of a modest increase in the number of measurements compared to the already perfectly calibrated problem. SparseLift is inspired by PhaseLift, an algorithm for phase retrieval [11,8], and by the blind deconvolution framework via convex programming due to Ahmed, Recht, and Romberg [1].
Self-calibration is a field of research in its own right, and we do not attempt in this paper to develop a framework that can deal with any kind of self-calibration problems. Instead our goal is more modest-we will consider a special, but at the same time in practice very important, class of self-calibration problems.
To that end we first note that many self-calibration problems can be expressed as (either directly or after some modification such as appropriate discretization) where y represents the given measurements, A(h) is the system matrix depending on some unknown calibration parameters h and x is the signal (function, image, ...) one aims to recover and w is additive noise. As pointed out, we do not know A(h) exactly due to the lack of calibration. Our goal is to recover x and h from y. More precisely, we want to solve If A(h) depends linearly on h, then (1.2) is a bilinear inverse problem. In any case, the optimization problem (1.2) is not convex, which makes its numerical solution rather challenging.
In many applications (1.1) represents an underdetermined system even if A(h) is perfectly known. Fortunately, we can often assume x to be sparse (perhaps after expressing it in a suitable basis, such as a wavelet basis). Moreover, in some applications h is also sparse or belongs to a low-dimensional subspace. Hence this suggests considering min where f (x, h) is a convex function that promotes sparsity in x and/or h. Often it is reasonable to assume that (1.3) is biconvex, i.e., it is convex in h for fixed x and convex in x for fixed h. In other cases, one can find a close convex relaxation to the problem. We can thus frequently interpret (1.3) as a biconvex compressive sensing problem.
However, useful as (1.3) may be, it is still too difficult and too general to develop a rigorous and efficient numerical framework for it. Instead we focus on the important special case y = D(h)Ax, (1.4) where D(h) is a diagonal matrix depending on the unknown parameters h and A is known. The reader may think that we finally reached a level of simplicity that is trivial to analyze mathematically and no longer useful in practice. However, as is often the case in mathematics, a seemingly simple model can be deceptive. As a matter of fact, on the one hand a rigorous analysis of the diagonal calibration model above requires quite non-trivial mathematical tools, and on the other hand this "simple" model does arise in numerous important applications, including radar, wireless communications, geolocation, direction finding, X-ray crystallography, cryo-electron microscopy, astronomy, medical imaging, etc. We will briefly discuss two such applications in more detail below.
Before moving on to applications, we pause for a moment to discuss the issue of identifiability. It is clear that if (D, x) is a pair of solutions to (1.4), then so is (αD, 1 α x) for any non-zero α. Therefore, it is only possible to uniquely reconstruct D and x up to a scalar. Fortunately for most applications this is not an issue, and we will henceforth ignore this scalar ambiguity when we talk about recovering D and x.

Array self-calibration and direction-of-arrival estimation
Consider the problem of direction-of-arrival estimation using an array with L antennas. Assuming that n signals are impinging on the array, the output vector (using the standard narrowband model) is given by (1.5) By discretizing the array manifold a(θ k , ψ k ), where (θ, ψ) are the azimuth and elevation angles, respectively, we can arrange equation (1.5) as y = D(h)Ax + w, (1.6) where w now also includes discretization error. Here a(θ, ψ) is a completely known function of (θ, ψ) and does not depend on calibration parameters, so A is a known matrix. The calibration is captured by the matrix D. An important special case is when D(h) is a diagonal matrix whose diagonal elements represent unknown complex gains associated with each of the antennas. The calibrating issue may come from position offsets of each antenna and/or from gain discrepancies caused by changes in the temperature and humidity of the environment. If the true signal directions (θ 1 ,ψ 1 ), . . . , (θ n ,ψ n ) lie on the discretized grid of angles (θ 1 , ψ 1 ) · · · , (θ N , ψ N ) corresponding to the columns of A, then the N × 1 vector x is n-sparse, otherwise x is only approximately n-sparse. In the direction finding problem we usually have multiple snapshots of the measurements y(t), in which case y and x in (1.6) become matrices, whose columns are associated with different snapshots [18].

The Internet of Things, random access, and sporadic traffic
The Internet of Things (IoT) will connect billions of wireless devices, which is far more than the current Fourth-Generation (4G) wireless system can technically and economically accommodate. One of the many challenges of 5G, the future and yet-to-be-designed wireless communication system, will be its ability to manage the massive number of sporadic traffic generating IoT devices which are most of the time inactive, but regularly access the network for minor updates with no human interaction [34]. Sporadic traffic will dramatically increase in the 5G market and it is a common understanding among communications engineers that this traffic cannot be handled within the current 4G random access procedures. Dimensioning the channel access according to classical information and communication theory results in a severe waste of resources which does not scale towards the requirements of the IoT. This means among others that the overhead caused by the exchange of certain types of information between transmitter and receiver, such as channel estimation, assignment of data slots, etc, has to be avoided. In a nutshell, if there is a huge number of devices each with the purpose to sporadically transmit a few bits of data, we must avoid transmitting each time a signal where the part carrying the overhead information is much longer than the part carrying the actual data. The question is whether this is possible.
In mathematical terms we are dealing with the following problem (somewhat simplified, but it captures the essence). For the purpose of exposition we consider the single user/single antenna case. Let h denote the channel impulse response and x the signal to be transmitted, where x is a sparse signal, its non-zero coefficients containing the few data that need to be transmitted. We encode x by computing z = Ax, where A is a fixed precoding matrix, known to the transmitter and the receiver. However, the non-zero locations of x are not known to the receiver, and nor is h, since this is the overhead information that we want to avoid transmitting, as it will change from transmission to transmission. The received signal is given by the convolution of h and z, i.e., y = h * z. The important point here is that the receiver needs to recover x from y, but does not know h either. This type of problem is often referred to as blind deconvolution. After applying the Fourier transform on both sides and writing D = diag(ĥ) (wherê denotes the Fourier transform), we obtainŷ = D Ax, which is of the form (1.4) modulo a change of notation.

Notation and outline
Matrices are written in uppercase boldface letters (e.g. X), while vectors are written in lowercase boldface letters (e.g. x), and scalars are written in regular font or denoted by Greek letters. For a matrix X with entries X ij we denote X 1 := ij |X ij |, X * is the nuclear norm of X, i.e., the sum of its singular values; X ∞ = max ij |X ij |, X is the operator norm of X and X F is its Frobenius norm. Tr(X) denotes the trace of X. If x is a vector, x is its 2-norm and x 0 denotes the number of non-zero entries in x, andx denotes the complex conjugate of x.
The remainder of the paper is organized as follows. In Section 2 we briefly discuss our contributions in relation to the state of the art in self-calibration. Section 3 contains the derivation of our approach and the main theorems of this paper. Sections 4 and 5 are dedicated to the proofs of the main theorems. We present numerical simulations in Section 6 and conclude in Section 7 with some open problems.

Related work and our contributions
Convex optimization combined with sparsity has become a powerful tool; it provides solutions to many problems which were believed intractable in the fields of science, engineering and technology. Two remarkable examples are compressive sensing and matrix completion. Compressive sensing [16,12,17], an ingenious way to efficiently sample sparse signals, has made great contributions to signal processing and imaging science. Meanwhile, the study of how to recover a low-rank matrix from a seemingly incomplete set of linear observations has attracted a lot of attention, see e.g. [5,7,29].
Besides drawing from general ideas in compressive sensing and matrix completion, the research presented in this work is mostly influenced by papers related to PhaseLift [3,11,26] and in particular by the very inspiring paper on blind deconvolution by Ahmed, Recht, and Romberg [1]. PhaseLift provides a strategy of how to reconstruct a signal x from its quadratic measurements y i = | x, z i | 2 via "lifting" techniques. Here the key idea is to lift a vector-valued quadratic problem to a matrix-valued linear problem. Specifically, we need to find a rank-one positive semi-definite matrix X = xx * which satisfies linear measurement equations of X: y i = Tr(XZ i ) where Z i = z i z * i . Generally problems of this form are NP-hard. However, solving the following nuclear norm minimization yields exact recovery of xx * : min X * , subject to y i = Tr(XZ i ), X 0.
(2.1) This idea of "Lifting" also applies to blind deconvolution [1], which refers to the recovery of two unknown signals from their convolution. The model can be converted into the form (3.1) via applying the Fourier transform and under proper assumptions. X 0 = h 0 x T 0 can be recovered exactly by solving a nuclear norm minimization program in the following form, min X * , subject to y i = Tr(XA i ), A i is a rank-one matrix, (2.2) if the number of measurements is at least larger than the dimensions of both h and x. Compared with PhaseLift (2.1), the difference is that X can be asymmetric and there is of course no guarantee for positivity. Ahmed, Recht, and Romberg derived a very appealing framework for blind deconvolution via convex programming, and we borrow many of their ideas for our proofs.
The idea of our work is motivated by the two examples mentioned above and our results are related to and inspired by those of Ahmed, Recht, and Romberg. However, there are important differences in our setting. First, we consider the case where one of the signals is sparse. It seems unnecessary to use as many measurements as the signal length and solve (2.2) in that case. Indeed, as our results show we can take much fewer measurements in such a situation. This means we end up with a system of equations that is underdetermined even in the perfectly calibrated case (unlike [1]), which forces us to exploit sparsity to make the problem solvable. Thus [1] deals with a bilinear optimization problem, while we end up with a biconvex optimization problem. Second, the matrix A in [1] is restricted to be a Gaussian random matrix, whereas in our case A can be a random Fourier matrix as well. This additional flexibility is relevant in real applications. For instance a practical realization of the application in 5G wireless communications as discussed in Section 1.2 and Section 6.3 would never involve a Gaussian random matrix as spreading/coding matrix, but it could definitely incorporate a random Fourier matrix.
Bilinear compressive sensing problems have also been invesstigated in e.g. [14,32,24]. These papers focus on the issue of injectivity and the principle of identifiability. Except for the general setup, there is little overlap with our work. Self-calibration algorithms have been around in the literature for quite some time, see e.g. [19,30,2,27] for a small selection, with blind deconvolution being a prominent special case, cf. [25] and the references therein. While many of these algorithms are definitely useful and interesting, most of them do not provide any guarantees of recovery. Also the proposed methods are very different from ours and do not involve sparsity, convex optimization, or ideas from compressive sensing. To the best of our knowledge, our method is the first work to bring self-calibration problems into the framework of biconvex compressive sensing. Uncertainties in the sensing matrix have been analyzed for instance in [22,13], but with the focus on quantifiying the uncertainty, while our current work can be seen as a means to mitigate the uncertainty in the sensing matrix.
3 Self-calibration, biconvex optimization, and SparseLift As discussed in the introduction, we are concerned with y = D(h)Ax + w, where the diagonal matrix D and x are unknown. Clearly, without further information it is impossible to recover x (or D) from y since the number of unknowns is larger than the number of measurements. What comes to our rescue is the fact that in many applications additional knowledge about D and x is available. Of particular interest is the case where x is sparse and the diagonal entries of D belong to some known subspace.
Hence, suppose we are given can be interpreted as the unknown calibrating parameters lying in the range of a linear operator B. E.g., B may be a matrix consisting of the first k columns of the DFT (Discrete Fourier Transform) matrix which represents a scenario where the diagonal entries in D are slowly varying from entry to entry. One may easily formulate the following optimization problem to recover h and x, We can even consider adding an ℓ 1 -component in order to promote sparsity in x, However, both programs above are nonconvex. For (3.2), it is an optimization with a bilinear objective function because it is linear with respect to one vector when the other is fixed. For (3.3), the objective function is not even bilinear but biconvex [20].

SparseLift
To combat the difficulties arising from the bilinearity of the objective function in (3.2), we apply a novel idea by "lifting" the problem to its feasible set [8,11] . It leads to a problem of recovering a rank-one matrix hx T instead of finding two unknown vectors. Denote by a i the i-th column of A T and by b i the i-th column of B * . A little linear algebra yields Let X := hx T and define the linear operator A : Now the problem is to find a rank-one matrix X satisfying the linear constraint A popular method is to take advantage of the rank-one property and use nuclear norm minimization.
where X * denotes the nuclear norm of X which is the sum of singular values.
Before we proceed, we need a bit of preparations. By using the inner product defined on C k×N as X, Z := Tr(XZ * ), the adjoint operator A * (u) : We also need a matrix representation Φ : vec(X) reshapes matrix X into a long column vector in lexicographic order. For each l, b * l Xa l = (a l ⊗ b l ) T vec(X) follows from the property of Kronecker product, whereb l is the complex conjugate of b l . It naturally leads to the block form of Φ * , With all those notations in place, the following formula holds Φ also has an explicit representation of its columns. Denote φ i,j := diag( b i ) a j , i = 1, 2, · · · , k and j = 1, 2, · · · , N. b i and a j are the i-th and j-th column of B and A respectively. Both b i and a j are of the size L × 1.
Then we can write Φ into . Concerning the support of x 0 , X 0 and vec(X 0 ), we can assume without loss of generality that the support of x 0 is its first n entries and thus X 0 := h 0 x T 0 have its first n columns nonzero. By a slight abuse of notation, we write Ω for the supports of x 0 , X 0 and vec(X 0 ). Let x Ω and X Ω denote the orthogonal projections of x and X on Ω, respectively. Define A Ω and A * Ω as the restriction of A and A * onto Ω, One can naturally define Φ Ω : Even if we stipulate that the degrees of freedom in D is smaller than L by assuming that the diagonal entries of D belong to some subspace, without any further knowledge about x the self-calibration problem is not solvable, since the number of unknowns exceeds the number of measurements. Fortunately, in many cases we can assume that x is sparse. Therefore, the matrix X = hx T is not only of rank one but also sparse. A common method to recover X from y is to use a linear combination of · * and · 1 as an objective function (See [26] for an example in phase retrieval) and solve the following convex program (3.10) tries to impose both sparsity and rank-one property of a matrix. While (3.10) seems a natural way to proceed, we will choose a somewhat different optimization problem for two reasons. First, it may not be straightforward to find a proper value for λ, see also Figures 2 and 3 in Section 6. Second, in [28] it is shown that if we perform convex optimization by combining the two norms as in (3.10), then we can do no better, order-wise, than an algorithm that exploits only one of the two structures. Since for any matrix Z there holds Z 1 ≥ Z * , this suggests that it may suffice to simply try to minimize X 1 , which in turn may already yield a solution with small nuclear norm. Indeed, as well will demonstrate, ℓ 1 -minimization performs well enough to recover X 0 exactly under certain conditions although it seemingly fails to take advantage of the rank-one property. Hence we consider This idea first lifts the recovery problem of two unknown vectors to a matrix-valued problem and then exploits sparsity through ℓ 1 -minimization 1 . We will refer to this method as SparseLift. It takes advantage of the "linearization" via lifting, while maintaining the efficiency of compressive sensing and linear programming.
If the observed data are noisy, i.e., y = A(X 0 ) + w with w < η, we take the following alternative convex program min X 1 , subject to It is useful to rewrite (3.11) in matrix-vector form. Noticing that A(V ) = Φv for v = vec(V ), V ∈ C k×N and v ∈ C kN ×1 , we have the following equivalent formulation of (3.11) where v 0 := vec(X 0 ). The matrix form of (3.12) is obtained similarly, This model attempts to recover X 0 from noisy observation y, where y = Φv 0 + w and w ≤ η.

Main results
Our main finding is that SparseLift enables self-calibration if we are willing to take a few more additional measurements. Furthermore, SparseLift is robust in presence of additive noise. In our theory we consider two choices for the matrix A: in the one case A is a Gaussian random matrix; in the other case A consists of rows chosen uniformly at random with replacement from the N × N Discrete Fourier Transform (DFT) matrix F with F * F = F F * = NI N . In the latter case, we simply refer to A as a random Fourier matrix.
Theorem 3.1. Consider the model (3.1) and assume that B is an L × k tall matrix with B * B = I k , h is a k × 1 dense vector and x is an n-sparse and N × 1 vector. Then the solutionX to (3.11) is equal to A is an L × N random Fourier matrix with L < N and where γ = 2N(log(2kN) + 1) + 1.
In both cases C α > 0 is growing linearly with α and µ max is defined as the largest absolute value in the matrix √ LB, i.e., The definition of µ max seems strict but it contains a lot of useful and common examples. For example, a matrix B has µ max = 1 if each column is picked from an L × L unitary DFT matrix and B * B = I k . We expect the magnitude of each entry in B does not vary too much from entry to entry and thus the rows of B are incoherent. Theorem 3.2. If the observed y is contaminated by noise, namely, y = A(X 0 ) + w with w ≤ η, the solutionX given by (3.12) has the following error bound.
1. If A is a Gaussian random matrix, with probability of success at least 1 − O(L −α+1 ). L satisfies (3.15) and both C 0 and C 1 are both constants.

If
A is a random Fourier matrix, with probability of success at least 1 − O(L −α+1 ). L and P satisfy (3.16) and both C 0 and C 1 are both constants The solutionX obtained by SparseLift, see (3.11) or (3.12) for noisy data, is not necessarily rank-one. Therefore, we compute the singular value decomposition ofX and setĥ = √σ 1û1 andx = √σ 1v1 , whereσ 1 is the leading singular value andû 1 andv 1 are the leading left and right singular vectors, respectively. Corollary 3.3. Letσ 1û1v T 1 be the best Frobenius norm approximation toX, then there exists a scalar α 0 and a constant c 0 such that 11,9] also derived error bounds of recovering the signals from the recovered "lifted" matrix by using the leading eigenvector or leading left and right singular vectors. [11,9] focus on PhaseLift and the recovered vector fromX is unique up to a unit-normed scalar. Our corollary is more similar to Corollary 1 in [1] because both of them deal with a bilinear system instead of a quadratic system like PhaseLift. There are two signals to recover fromX in our case and thus they are unique up to a scalar α 0 . Without prior information, no theoretical estimate can be given on how large α 0 is.
The proof of Corollary 3.3 consists of considering the leading eigenvalue and eigenvector of the Hermitian matricesX * X andXX * and then using the sin-θ Theorem in [15]. The same idea is implemented in the proof of Theorem 1.2 from [11] where the detailed reasoning can be found.

Proof of Theorem 3.1
The outline of the proof of Theorem 3.1 follows a meanwhile well-established pattern in the literature of compressive sensing and low-rank matrix recovery, but as is so often the case, the difficulty lies in the technical details. We first derive sufficient conditions for exact recovery. The fulfillment of these conditions involves computing norms of certain random matrices as well as the construction of a dual certificate. The main part of the proof is dedicated to these two steps. While the overall path of the proof is similar for the case when A is a Gaussian random matrix or a random Fourier matrix, there are several individual steps that are quite different.

A sufficient condition for exact recovery
We will first give a sufficient condition for exact recovery by using SparseLift. Then we prove all the proposed assumptions are valid. This proposition can be regarded as a sufficient condition of exact recovery of a sparse vector via ℓ 1 minimization in compressive sensing (See Theorem 4.26 in [17]). Based on (4.1), we propose an alternative sufficient condition which is much easier to use. Proposition 4.2. X 0 is the unique minimizer if there exists a Y ∈ Range(A * ) satisfying the following conditions, for any H : A(H) = 0.
Proof: The proof starts with (4.1). By splitting the inner product of sgn(X 0 ) − Y , H on Ω and Ω ⊥ , the following expression can guarantee the exact recovery of l 1 minimization, By using the Hölder inequality, we arrive at the stronger condition The second equality follows from A(H Ω ) + A(H Ω ⊥ ) = 0. Substitute (4.5) into (4.4) and we only need to prove the following expression in order to achieve exact recovery,

Weak restricted isometry property
In this section, we will prove Lemma 4.3, the weak version of the restricted isometry property, which achieves the first condition in (4.2). With Lemma 4.3, one can easily derive holds true with probability at least 1 − L −α+1 where I Ω is the identity operator on Ω if 1. a l is the l-th row of an L×N Gaussian random matrix A and L ≥ C α µ 2 max kn max{log L/δ 2 , log 2 L/δ}. 2. a l is the l-th row of an L × N random Fourier matrix A and L ≥ C α µ 2 max kn log L/δ 2 . In both cases, C α grows linearly with respect to α.
In fact, it suffices to prove Lemma 4.6 and Lemma 4.7 because Lemma 4.3 is a special case of them when Q = L, P = 1.

Two ingredients of the proof
Before we move to the proof of Lemma 4.6 and Lemma 4.7, we introduce two important ingredients which will be used later. One ingredient concerns the rows {b l } L l=1 of B. There exists a disjoint partition of index set {1, 2, · · · , L} into {Γ p } P p=1 with each |Γ p | = Q and Q > Cµ 2 max k log(L) such that The existence of such a partition follows from the proof of Theorem 1.2 in [10] and is used in Section 3.2 of [1]. (4.6) also can give a simple bound for the operator norm of both T p and T −1 p .
The other ingredient concerns the non-commutative Bernstein inequality. For Gaussian and Fourier cases, the inequalities are slightly different. This following theorem, mostly due to [31], is used for proving the Fourier case of (4.19). Compute the variance parameter then for all t ≥ 0 For the Gaussian case we replace Z l ≤ R by Z l ψ 1 ≤ R where the norm · ψ 1 of a matrix is defined as The concentration inequality is different if a l is Gaussian because of the unboundedness of S. We are using the form in [23]. Theorem 4.5. For a finite sequence of independent M ×M random matrices Z l with R := max 1≤l≤Q Z l ψ 1 and S and σ 2 as defined in (4.8) and (4.9), we have the tail bound on the operator norm of S, where C 0 is an absolute constant.

A key lemma and its proof
For each subset of indices Γ p , we denote Under certain conditions, it is proven that A * p,Ω A p,Ω does not deviate from the identity much for each p. We write this result into Lemma 4.6. Lemma 4.6. For any equal partition {Γ p } P p=1 of {1, 2, · · · , L} with |Γ p | = Q such that (4.6) holds true, we have with probability at least 1 − P L −α ≥ 1 − L −α+1 if we choose Q ≥ C α kn log(L) where C α grows linearly with respect to α and 1. a l is the l-th row of an L × N Gaussian random matrix A and Q ≥ C α µ 2 max kn log 2 L. 2. a l is the l-th row of an L × N random Fourier matrix A and Q ≥ C α µ 2 max kn log L. First, we prove the following key lemma. Lemma 4.6 is a direct result from this key lemma and from (4.6). We postpone the proof of Lemma 4.6 to Section 4.2.3. Lemma 4.7. For any fixed 0 < δ ≤ 1 and equal partition

13)
with probability at least 1 − P L −α ≥ 1 − L −α+1 if 1. a l is the l-th row of an L×N Gaussian random matrix A and Q ≥ C α µ 2 max kn max{log L/δ 2 , log 2 L/δ}. 2. a l is the l-th row of an L × N random Fourier matrix A and Q ≥ C α µ 2 max kn log L/δ 2 . In both cases C α grows linearly with respect to α.
Proof: We prove the Gaussian case first; the proof of the Fourier case is similar. Let a l be a standard Gaussian random vector. a l,Ω denotes the projected vector a l on the support of x 0 with | supp(x 0 )| = n. We first rewrite the linear operator A * p,Ω A p,Ω (X) − T p X Ω into matrix form. Using the Kronecker product we define Z l as follows: Every Z l is a centered random matrix since E(a l,Ω a * l,Ω ) = I N,Ω , which follows from restricting a l a * l to the support Ω and taking its expectation. I N,Ω is an N × N matrix with I n as its n × n leading submatrix. It is easy to check that vec(A * p,Ω A p,Ω (X) − T p X Ω ) = l∈Γp Z l vec(X Ω ).
More precisely, a l,Ω in Z l needs to be replaced byā l,Ω according to the property of Kronecker product. But there is no problem here since each a l is a real random vector. For simplicity, we denote Z l (X) := b l b * l X(a l,Ω a * l,Ω − I N,Ω ). In order to use Theorem 4.5 to estimate L l=1 Z l , we need to know Z l ψ 1 and σ 2 . We first give an estimate for the operator norm of Z l .
Here a l,Ω 2 is a χ 2 random variable with n degrees of freedom. This implies that Z l ψ 1 ≤ Cµ 2 max kn L , following from Lemma 6 in [1]. Therefore there exists an absolute constant C b such that The second step is to estimate σ 2 . It suffices to consider Z * l Z l (X) because Z l is a Hermitian matrix.
The expectation of Z * l Z l (X) is following from E(a l,Ω a * l,Ω − I N,Ω ) 2 = (n + 1)I N,Ω .
Now we have estimates of R and σ 2 and can apply the concentration inequality (4.10). Z l can be simply treated as a kn × kn matrix for each l because its support is only kn × kn. There holds log √ QR σ 2 ≤ C 1 log L for some constant C 1 because Q, R and σ 2 are at most of order L. Hence, If we let Q ≥ C 0 µ 2 max kn( 5 2 + C b δ log L)(α log L + log 2kn)/δ 2 and also properly assume L > kn, then Q ≥ C α µ 2 max kn max{log L/δ 2 , log 2 L/δ} and where C α is a function which grows linearly with respect to α. Now we take the union bound over all p = 1, 2, · · · , P and have This completes the proof for the Gaussian case.
We proceed to the Fourier case. The proof is very similar and the key is to use Theorem 4.4. Now we assume a l is the l-th row of a random Fourier matrix A and thus E(a l a * l ) = I N . Here Z l = (ā l,Ωā * l,Ω − I N,Ω ) ⊗ b l b * l because a l is a complex random vector. However, it does not change the proof at all except the notation. Like the proof of the Gaussian case, we need to know R = max l∈Γp Z l and also σ 2 . Following the same calculations that yielded (4.15), (4.16) and (4.18), we have Then we apply Theorem 4.4 and obtain

Estimating the operator norm of A
The operator norm γ of A with A will be given for both the Gaussian case and the Fourier case such that A(X) F ≤ γ X F for any X with high probability. For the Gaussian case, Lemma 1 in [1] gives the answer directly. Lemma 4.8. For A defined in (3.4) with A a Gaussian random matrix and fixed α ≥ 1, with probability at least 1 − L −α .
For the Fourier case, a direct proof will be given by using the non-commutative Bernstein inequality in Theorem 4.4. Lemma 4.9. For A defined in (3.4) with A a random Fourier matrix and fixed α ≥ 1, its operator norm γ yields γ ≤ 2N(log(2kN) + 1) + 1 (4.21) with probability at least 1 − L −α if we choose L ≥ αk log(L).
Proof: First note that A * A = A 2 . Therefore it suffices to estimate the norm of A * A, The second equality follows from the fact that L l=1 b l b * l = I k . For each centered random matrix Z l := (ā lā * l − I N ) ⊗ b l b * l , we have vec(A * A(X)) = L l=1 Z l vec(X) + vec(X) and That means R = max 1≤l≤L (ā lā * . In order to estimate the variance of each Z l , we only have to compute E(Z l Z * l ) since Z l is a Hermitian matrix.
Hence the variance of L l=1 Z l is Now we apply Bernstein's inequality and for any t ≥ 0, we have If we pick t = 2N(log(2kN) + 1), then If L ≥ αµ 2 max k log(L), we have

Construction of the dual certificate
In the last section we have shown that (4.2) holds true. Now we will prove that a Y ∈ Range(A * ) satisfying (4.3) exists by explicitly constructing it for both, the Gaussian and the Fourier case.

A is a Gaussian random matrix
Exact dual certificate via least square method: It follows from Lemma 4.3 and 4.6 that Φ * Ω Φ Ω is invertible on Ω with (Φ * Ω Φ Ω ) −1 ≤ 2 if L > C 0 kn log 2 L for some scalar C 0 . We construct the exact dual certificate via (see [6]) Let Y = A * (p); it suffices to consider q and verify whether q Ω = vec(X 0 ) and q Ω ⊥ ∞ ≤ 1 2 .
Proof: Note Φ * Ω ⊥ p = ( φ * s,t p) 1≤s≤k,t>n and φ s,t = diag( b s ) a t is independent of p for t > n because p only uses a j : 1 ≤ j ≤ n. In this case, φ * s,t p = a * t diag( b * s )p where a t is an L × 1 standard Gaussian random vector. Due to the independence between diag( b * s )p and Therefore, φ * s,t p is a centered Gaussian random variable with variance σ 2 at most 2µ 2 max kn L . For any x > 0, we have Taking the union bound over all pairs (s, t) of indices ( φ * s,t ) satisfying t > n gives Pr max If we pick x = 1 2 and σ x ≤ 2 2µ 2 max kn L , Pr max Choosing L > 16αµ 2 max nk(log(L) + log(kN)) impies that the probability of {max 1≤s≤k,t>n | φ * s,t p| > 1 2 } will be less than L −α .

A is a random Fourier matrix
An inexact dual certificate Y will be constructed to satisfy the conditions given in (4.3) via the golfing scheme, which was first proposed by David Gross in [21] and has been widely used since then, see e.g. [1,17,26]. The idea is to use part of the measurements in each step and construct a sequence of random matrices Y p recursively in Range(A * ). The sequence converges to sgn(X 0 ) exponentially fast and at the same time keeps each entry of Y p small on Ω ⊥ . We initialize Y 0 := 0 and set We denote Y := Y P . Let W p := Y p,Ω − sgn(X 0 ) be the residual between sgn(X 0 ) and Y p on Ω and W 0 := Y 0 − sgn(X 0 ) = − sgn(X 0 ). W p yields the following relations due to (4.23), where A p,Ω is defined in (4.11). Following from Lemma 4.6, W p F decreases exponentially fast with respect to p with probability at least 1 − L −α+1 , The value of P can be determined to make The next step is to bound the infinity norm of Y on Ω ⊥ . We first rewrite the recursive relation (4.23); Y has the following form, Our aim is to achieve Y Ω ⊥ ∞ ≤ 1 2 . It suffices to make P Ω ⊥ is the operator which orthogonally projects a matrix onto Ω ⊥ . For each p, has the same support as W p−1 does and thus this gives L for all 1 ≤ p ≤ P with probability at least 1 − L −α+1 if L ≥ P Q and Q ≥ C α µ 2 max kn(log(kN) + log(L)). Theorem 4.11. Conditioned on Lemma 4.6 and (4.24), for fixed p, .
In particular this gives kn(log(kN) + log(L)) and taking the union bound for all p.
Proof: For any z = vec(Z) where Z : k × N has support Ω, we have (3.4) and (3.5). Therefore, where e s is a kN × 1 unit vector with its s-th entry equal to 1. Here s > kn because only the first kn entries of z are nonzero. The main idea of the proof is to estimate | e s , Φ * p Φ p z | for each s ∈ Ω ⊥ and take the union bound over all s ∈ Ω ⊥ . For any fixed s ∈ Ω ⊥ , Let z l := e s , φ l φ * l z . Notice E(z l ) = 0 because b l b * l Z and Z have the same support and thus any entry of b l b * l Z on Ω ⊥ equals to 0. Therefore each z l is a bounded and centered random variable. We can use scalar version of Bernstein inequality (a special version of Theorem 4.4) to estimate l∈Γp z l . First we estimate |z l |, There holds |e * s φ l | ≤ ā l ⊗ b l ∞ = ā l ∞ b l ∞ ≤ µmax √ L , which follows from the fact that a l is a row of DFT matrix and each entry of a l is of unit norm. Moreover, z has support Ω, The variance of l∈Γp z l is estimated as Now we apply the Bernstein inequality in Theorem 4.4 and have the estimate for l∈Γp z l , In particular, if we choose z = vec(W p−1 ) with z ≤ 2 −p+1 √ kn and t = 2 −p−1 Q L , we have We take the union bound over all s ∈ Ω ⊥ and obtain

Stability theory
Now we consider the noisy case by solving the convex program (3.14). Our goal is to prove Theorem 3.2.
The main ingredient of our proof is Theorem 4.33 in [17]. A variant of that theorem, adapted to our setting, reads as follows. Theorem 5.1. Let φ s,t be the columns of Φ ∈ C L×kN , v 0 = vec(X 0 ) ∈ C kN ×1 with support on Ω and y = Φv 0 + w with w ≤ η. For δ, β, θ, η, τ ≥ 0 and δ < 1, and that there exists a vector q = Φ * p ∈ C kN ×1 with p ∈ C L such that where C 0 and C 1 are two constants.
Actually, we have already proven in the previous section that we can choose δ = 1 2 , θ = 1 2 . The only two unknown components are τ and β. τ follows from computing the norm of p. We need to know the coherence µ of Φ to estimate β and claim that β ≤ 1. Define the coherence as the largest correlation between two different columns of the L × kN matrix Φ, where φ s,t is defined in (3.7).

Estimation of µ and β.
It is easy to obtain a bound of Φ * Ω φ s,t for 1 ≤ s ≤ k and t > n with µ, It suffices to bound µ by 1/ √ kn in order to achieve β ≤ 1. The following lemma (see Lemma A.9 in [17]) will be used in the proof. Lemma 5.2. For any matrix C, any entry is bounded by its operator norm, i.e., A is a Gaussian random matrix and L ≥ 4Cαµ 2 max max{kn log(L), √ kn log 2 L} andα = 2 log(kN) + α.

2.
A is a random Fourier matrix and L ≥ 4Cαµ 2 max kn log(L) andα = 2 log(kN) + α. Here, in both cases Cα is a function growing at most linearly with respect toα.
where γ is the operator norm of A and is certainly greater than 1. Applying Theorem 5.1 gives the desired result.

Proof of Theorem 3.2 when A is a random Fourier matrix
Proof: We have an explicit form of p from (4.26) and Y is an inexact dual certificate. The norm of p is For each A p (W p−1 ), which follows from (4.24) and A * p,Ω A p,Ω ≤ 3Q 2L in Lemma 4.6. Therefore, This implies τ = 2 √ 2P where P is chosen in (4.25). We also have δ = 1 2 , θ = 1 2 and β ≤ 1 which are the same as Gaussian case and ρ = θ + β 4 √ 2γ(1−δ) < 1. Now we have everything needed to apply Theorem 5.1. The solutionX to (3.14) satisfies with P = O(log(N)), and C 0 , C 1 are two constants.

.1 Recovery performance of SparseLift
We investigate the empirical probability of success of SparseLift over different pairs of (k, n) for fixed L and N. We fix L = 128, N = 256 and make both k and n range from 1 to 15. We choose B as the first k columns of an L × L unitary DFT matrix and A is sampled independently for each experiment. x has random supports with cardinality n and the nonzero entries of h and x yield N(0, 1). For each (k, n), we run (3.11) 10 times and compute the probability of success. We classify a recovery as a success if the relative error is less than 1%, The results are depicted in Figure 1. The phase transition plots are similar for the case when A is a Gaussian random matrix or a random Fourier matrix. In this example choosing kn < 70 gives satisfactory recovery performance via solving SparseLift.
For comparison, we illustrate the performance when we solve the convex optimization problem that combines the nuclear norm and sparsity as in (3.10), see Figures 2 and 3. The values for λ are λ = 0.1 (the best value for λ for this experiment) and λ = 10. The settings of other parameters are exactly the same as the scenario of solving SparseLift. Figure 2 shows that (3.10) with λ = 0.1 performs similar to solving (3.11). However, if λ = 10 the results are clearly worse, in particular for the Gaussian case, see Figure 3. It supports our hypothesis that the performance of recovery is sensitive to the choice of λ and SparseLift is good enough to recover calibration parameters and the signal. This is in line with the findings in [28] which imply that convex optimization with respect to a combination of the norms · 1 and · * will not perform significantly better than using only one norm, the ℓ 1 -norm in this case.

Self-calibration and direction-of-arrival estimation
We examine the array self-calibration problem discussed in Section 1.1. We consider a scenario of three fully correlated equal-power signals. Note that the case of fully correlated signals is much harder than the case of uncorrelated (or weakly correlated) signals. Standard methods like MUSIC [33] or methods utilizing a mixed ℓ 1 -ℓ 2 sparsity model as in [35] fail completely for fully correlated signals even in the absence of any calibration error. The reason is that these methods rely on the additional information obtained by taking multiple snapshots. Yet, for fully correlated signals all snapshots are essentially identical (neglecting noise) and therefore no new information is contained in multiple snapshots.
We use a circular array consisting of L = 64 omni-directional sensors. We assume that the calibration error is changing slowly across sensors, as for example in the case of sensors affected by changing temperature or humidity. There are of course many possibilities to model such a smooth change in calibration across sensors. In this simulation we assume that the variation can be modeled by a low-degree trigonometric polynomial (the periodicity inherent in trigonometric polynomials is fitting in this case, due the geometry of the sensors); in this example we have chosen a polynomial of degree 3. To be precise, the calibration error, i.e., the complex gain vector d is given by d = Bh, where B is a 64 × 4 matrix, whose columns are the first four columns of the 64 × 64 DFT matrix and h is a vector of length 4, whose entries are determined by taking a sample of a multivariate complex Gaussian random variable with zero mean and unit covariance matrix. The stepsize for the discretization of the grid of angles for the matrix A is 1 degree. The directions of arrival are set to −10 • , 5 • , and 20 • , and the SNR is 25dB. Since the signals are fully correlated we take only one snapshot.
We apply SparseLift, cf. (3.14), to recover the angles of arrival. Since the sensing matrix A does not fulfill the conditions required by Theorem 3.2 and the measurements are noisy, we cannot expect SparseLift to return a rank-one matrix. We therefore extract the vector with the angles of arrival by computing the right singular vector associated with the largest singular value ofX. As can be seen in Figure 6.2 the estimate obtained via self-calibration provides a quite accurate approximation to the true solution.  Figure 2: Phase transition plot of performance via · 1 + λ · * minimization for λ = 0.1.

5G, the Internet of Things, and SparseLift
We briefly discussed in Section 1.2 how self-calibration problems of the form analyzed in this paper are anticipated to play a role in 5G wireless communication systems in connection with the Internet of Things, sporadic traffic and random access. Let us assume we deal with a CDMA communication scheme. To accommodate as many sporadic users as possible, we consider an "overloaded" CDMA system. This means that the matrix A containing the so-called spreading sequences (aptly called spreading matrix) has many more columns than rows. We use a random Fourier matrix for that purpose. Furthermore, we assume that the transmission channel is frequency-selective and time-invariant for the time of transmission of one sequence, but it may change from transmission to transmission. The channel is unknown at the receiver, it only knows an upper bound on the delay spread (the length of the interval containing the non-zero entries of h). Moreover, the transmitted signal experiences additive Gaussian noise. The mathematician who is not familiar with all this communications jargon may rest assured that this is a standard scenario in wireless communication. In essence it translates to a linear system of equations of the form y = DAx + w, where A is a flat random Fourier matrix, D is of the form D = diag(Bh) and w is AWGN. To be specific, we choose the transmitted signal x to be of length N = 256 with n = 5 data-bearing coefficients (thus, x 0 = 5), the locations of which are chosen at random. The spreading sequences are of length L = 128 (thus the system is overloaded by a factor of 2). The delay spread is assumed to be 5 and the SNR ranges from 0dB to 80dB (the latter SNR level is admittedly unrealistic in practical scenarios).
The result is depicted in Figure 6.3; for comparison we also show the results when using a Gaussian random matrix as spreading matrix. We plot the relative error (on a dB scale) for different levels of SNR (also on a dB scale). The horizontal axis represents SNR = 10 log 10 X 0 2 F w 2 and the vertical axis represents 10 log 10 (Avg.Err.). We run 10 independent experiments for each SNR level and take the average relative error. The error curves show clearly the desirable linear behavior between SNR and mean square error with respect to the log-log scale.  Figure 3: Phase transition plot of performance via · 1 + λ · * minimization for λ = 10.

Conclusion and outlook
Self-calibration is an emerging concept that will get increasingly important as hardware and software will become more and more intertwined. Several instances of self-calibration, such as blind deconvolution, have already been around for some time. In this work we have taken a step toward building a systematic mathematical framework for self-calibration. There remain more open questions than answers at this stage. Some obvious extensions of the results in this paper will be to consider the case where both x and h are sparse. Another useful generalization is the scenario where x is replaced by a matrix, like in the directionof-arrival problem with multiple snapshots. Lifting the self-calibration problem in this case would lead to a tensor-valued (instead of a matrix-valued) underdetermined linear system [18]. The computational complexity for solving the associated relaxed convex problem will likely be too expensive for practical purposes and other methods will be needed. Extensions of the Wirtinger-Flow approach in [4] seem to be a promising avenue to derive computationally efficient algorithms for this and other self-calibration scenarios. Another important generalization is the case when there is mutual coupling between sensors, in which case D is no longer a diagonal matrix. Yet another very relevant setting consists of D depending non-linearly on h, as is often the case when the sensors are displaced relative to their assumed locations.
In conclusion, we hope that it has become evident that self-calibration is on the one hand an inevitable development in sensor design in order to overcome some serious current limitations, and on the other hand it is a rich source of challenging mathematical problems.