TL;DR
This paper introduces a novel semidefinite programming relaxation for finding the nearest rank deficient matrix within an affine space, providing provable guarantees and improved noise robustness over existing methods.
Contribution
The authors propose the first convex SDP relaxation for the nearest structured rank deficient matrix problem with theoretical guarantees in low noise conditions.
Findings
The SDP relaxation always finds the global minimizer in low noise regimes.
The method outperforms existing techniques in noise tolerance.
Applications include system identification, GCD approximation, and camera resectioning.
Abstract
Given an affine space of matrices and a matrix , consider the problem of computing the closest rank deficient matrix to on with respect to the Frobenius norm. This is a nonconvex problem with several applications in control theory, computer algebra, and computer vision. We introduce a novel semidefinite programming (SDP) relaxation, and prove that it always gives the global minimizer of the nonconvex problem in the low noise regime, i.e., when is close to be rank deficient. Our SDP is the first convex relaxation for this problem with provable guarantees. We evaluate the performance of our SDP relaxation in examples from system identification, approximate GCD, triangulation, and camera resectioning. Our relaxation reliably obtains the global minimizer under non-adversarial noise, and its noise tolerance is significantly…
| group | I | II | II | II |
|---|---|---|---|---|
| multiplier | ||||
| equation | ||||
| indices |
| SDP succeeds | LMA and SDP succeed | |||||||||||||||
| 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
| 3 | 100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
99
% |
99
% |
95
% |
89
% |
83
% |
80
% |
76
% |
74
% |
71
% |
70
% |
| 4 | 100
% |
100
% |
100
% |
97
% |
92
% |
85
% |
79
% |
93
% |
85
% |
79
% |
73
% |
66
% |
58
% |
52
% |
||
| 5 | 100
% |
100
% |
99
% |
95
% |
87
% |
77
% |
93
% |
84
% |
78
% |
71
% |
62
% |
53
% |
||||
| 6 | 100
% |
100
% |
98
% |
91
% |
81
% |
92
% |
81
% |
75
% |
66
% |
58
% |
||||||
| 7 | 100
% |
100
% |
97
% |
89
% |
90
% |
79
% |
73
% |
63
% |
||||||||
| noise | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 |
|---|---|---|---|---|---|---|
| SDP | 100
% |
100
% |
100
% |
99
% |
98
% |
95
% |
| LMA | 100
% |
100
% |
99
% |
86
% |
67
% |
54
% |
| BFGS | 100
% |
100
% |
99
% |
86
% |
67
% |
54
% |
| 38
% |
76
% |
|||||||||||
| noise | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 |
| SDP | 100
% |
100
% |
100
% |
98
% |
94
% |
92
% |
100
% |
93
% |
80
% |
82
% |
82
% |
84
% |
| LMA | 100
% |
11
% |
0
% |
0
% |
0
% |
0
% |
100
% |
1
% |
0
% |
0
% |
0
% |
0
% |
| BFGS | 100
% |
100
% |
92
% |
68
% |
52
% |
45
% |
100
% |
38
% |
24
% |
20
% |
16
% |
17
% |
| noise | 0.0 | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 |
|---|---|---|---|---|---|---|
| SDP | 100
% |
100
% |
96
% |
91
% |
90
% |
90
% |
| LMA | 100
% |
100
% |
88
% |
71
% |
60
% |
60
% |
| BFGS | 100
% |
100
% |
88
% |
72
% |
60
% |
60
% |
| cameras on a sphere | cameras on a line | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| noise | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | |
| our SDP | 100
% |
100
% |
100
% |
100
% |
100
% |
99
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
|
| AAT [1] | 100
% |
93
% |
87
% |
81
% |
74
% |
68
% |
100
% |
0
% |
0
% |
0
% |
0
% |
0
% |
|
| our SDP | 100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
100
% |
99
% |
|
| AAT [1] | 100
% |
100
% |
99
% |
96
% |
88
% |
79
% |
100
% |
0
% |
0
% |
0
% |
0
% |
0
% |
|
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A convex relaxation to compute the
nearest structured rank deficient matrix
Diego Cifuentes
Department of Mathematics, Massachusetts Institute of Technology, Cambridge MA 02139, USA
Abstract.
Given an affine space of matrices and a matrix , consider the problem of computing the closest rank deficient matrix to on with respect to the Frobenius norm. This is a nonconvex problem with several applications in control theory, computer algebra, and computer vision. We introduce a novel semidefinite programming (SDP) relaxation, and prove that it always gives the global minimizer of the nonconvex problem in the low noise regime, i.e., when is close to be rank deficient. Our SDP is the first convex relaxation for this problem with provable guarantees. We evaluate the performance of our SDP relaxation in examples from system identification, approximate GCD, triangulation, and camera resectioning. Our relaxation reliably obtains the global minimizer under non-adversarial noise, and its noise tolerance is significantly better than state of the art methods.
Key words and phrases:
Convex relaxation, Low rank approximation, Semidefinite programming, Total least squares, Affine structure
1. Introduction
Let be positive integers with . Given an affine map and a vector , we consider the nearest point problem
[TABLE]
where we use the -norm in . More generally, we may consider the weighted -(semi)norm , where is positive (semi)definite. Denoting the affine space and the matrix , we may equivalently write (1) as
[TABLE]
where we use the (possibly weighted) Frobenius norm in .
We may further study a variant of (1) with complex variables, but this can be reduced to the real case. Indeed, a complex matrix is rank deficient if and only if the real matrix is rank deficient. So we restrict our attention to the real case.
Problem (1) is sometimes known as structured total least squares (STLS) [30, 13, 37, 24, 25]. The name comes from the total least squares method for the linear regression problem , given by minimizing the Frobenius norm subject to being rank deficient. By further imposing an affine constraint on the matrix, we arrive to (1).
There are several applications of (1) in control, systems theory, and statistics [30, 27, 28, 9, 13, 22]. Lemmerling’s thesis [24] gives a long list of concrete applications, including noisy deconvolution, linear prediction, direction of arrival estimation, information retrieval, approximate realization, and model reduction. The approximate GCD problem [22] from computer algebra is another instance of (1). There are also many applications in computer vision. Specifically, in the area of multi-view geometry one often has to minimize the following fractional programming problem:
[TABLE]
for some given vectors . After clearing denominators, this takes the form of (1) with . In particular, the triangulation, camera resectioning, and homography estimation problems can be phrased as above [20].
Problem (1) is nonconvex and computationally hard. It is usually solved in practice with local optimization methods. Although these methods are typically fast, there are no guarantees of convergence to the global optimum. An alternative approach to solve nonconvex problems is to construct a convex relaxation. This idea has been applied in several other areas, e.g., [41, 7, 36, 18, 12], often leading to global optimality results. We follow this approach in this paper and propose a novel convex relaxation for (1).
An important advantage of convex relaxations over local optimization methods is certifiability. Even if a local optimization method successfully finds the global minimum, there is no way for the user to verify this. On the other hand, convex relaxations may leverage convex duality to provide certificates of global optimality. Guaranteed performance is a critical feature in certain applications; its role in robotics is discussed in [41, Appx.A]. Convex relaxations might also be combined with local optimization methods, serving for a posteriori verification of candidate solutions.
Our convex relaxation, presented in (13), relies on semidefinite programming (SDP), and can be solved in polynomial time with interior points methods [4, §6]. The optimal value of our relaxation is always a lower bound on the global minimum of (1). We say that the relaxation is exact (or tight) if the optimal solution of the SDP is a rank-one matrix. In such a case the optimal value of the SDP agrees with the optimal value of (1) and we can recover the global minimizer of (1) from the SDP solution.
We evaluated the performance of our SDP relaxation in the following applications: approximate realization (control), approximate GCD (computer algebra), triangulation and camera resectioning (computer vision). Our experiments reveal that our SDP is often exact (and certifiably optimal) under non-adversarial noise. Furthermore, the SDP is significantly more resilient to noise than state of the art methods. Preliminary experimental results suggest that the SDP is also robust against missing data (accommodated by using a degenerate weight matrix ).
Exactness under low noise
Our main technical contribution is to show that our SDP relaxation of (1) is always exact when is sufficiently close to the set . In statistical estimation problems the point often represents a noisy sample from , so the distance from to corresponds to the magnitude of the noise. We hence refer to this setting as the low noise regime. This gives a theoretical explanation for the practical success of the relaxation in our experiments.
More precisely, we show that if satisfies a mild geometric condition, then the SDP relaxation is exact as . The geometric condition we use is non-tangentiality, as introduced in [2]. Two manifolds meet non-tangentially at a point if
[TABLE]
where denotes the tangent space at . Informally, this condition means that the two manifolds have a positive angle in directions perpendicular to . As shown in [2], non-tangentiality is a weaker condition than the more common notion of transversality (which means that is the whole space).
Assumption 1**.**
Let be the manifold of matrices with rank deficiency one. Assume that meets the affine space non-tangentially at the point .
Theorem 1.1**.**
Let such that Assumption 1 holds. Then the SDP relaxation (13) correctly recovers the global minimizer of (1) whenever is close enough to .
We will show in Section 4 that the role of Assumption 1 is to ensure that the intersection is a regularly defined smooth manifold. Due to the pivotal role of regularity (constraint qualifications) in optimization, the non-tangentiality condition becomes essential for any (local or global) optimization method for (1). For instance, Assumption 1 is needed to get linear convergence in alternating projections [2] and to prove quadratic convergence of Newton’s method [38].
Our proof of Theorem 1.1 relies on the theory of stability of SDP relaxations recently developed in [10]. More generally, the proof can be extended to the case of weighted -norms defined by a positive definite weight matrix.
Related work
Problem (1) has appeared in several works, often by the name STLS; see [30, 13, 37, 24, 25] and the references therein. All previous methods rely on local optimization. Due to nonconvexity, there are no guarantees to get the global minimum.
Recent work has shifted into a generalization of (1) known as structured low rank approximation (SLRA), given by for some given . Hence, (1) corresponds to the case . The SLRA problem has been studied extensively [9, 38, 32], notably by Markovsky [28, 27]. As before, most practical algorithms rely on local optimization, see e.g., [28, §3].
Convex optimization methods have been applied to SLRA. A common heuristic method consists in replacing the rank constraint by the convex constraint , where is the nuclear norm and is a tunable parameter. This leads to a tractable SDP, see [15]. The rationale behind this heuristic is that the nuclear norm is the convex envelope of the rank function. But there are no theoretical guarantees for the SLRA problem. The resulting matrix is not necessarily optimal, may not satisfy the rank constraint, and optimality certificates cannot be recovered.
Another convex relaxation for SLRA was proposed in [23], see also [8]. It relies on a characterization of the convex envelope for , where denotes the indicator function. This gives a better approximation of the SLRA problem than the nuclear norm, as it takes into account the -objective and the rank (but ignores the affine structure ). This relaxation allows to derive optimality certificates. However, no provable guarantees or complexity estimates are known.
To the best of our knowledge, our proposed SDP is the first convex relaxation for problem (1) with provable guarantees. It is also the first polynomial-time method for (1) that provides global optimality certificates. Special instances from computer vision have received more attention [21, 1], as we will elaborate in Section 6.3.
A closely related problem to SLRA is the rank minimization problem (RMP), given by for some parameter . RMP has been widely studied in the compressing sensing community. Several theoretical guarantees are known for both convex [35] and nonconvex methods [19]. Both cases critically rely on the restricted isometry property (RIP), which gives strong limitations on the affine structure . Note that RIP does not make sense in the context of SLRA, since it implies that there is a unique vector for which .
Since (1) is a polynomial optimization problem, the sum-of-squares (SOS) method [5] provides a hierarchy of SDP relaxations with increasing power and complexity. Our proposed relaxation involves a positive semidefinite (PSD) matrix of size , and lies in between the first and the second level of the SOS hierarchy, whose sizes are and . Though SOS relaxations often perform well in practice, there are very few theoretical results explaining this behavior. This paper contributes in this direction, by showing that our relaxation is always exact in the low noise regime.
Structure of the paper
In Section 2 we review semidefinite programming relaxations of quadratically constrained quadratic programs. In Section 3 we derive our SDP relaxation for (1). In Section 4 we elaborate on the connection between non-tangentiality, smoothness, and constraint qualifications. In Section 5 we prove Theorem 1.1. We conclude in Section 6 with experimental results.
Notation
We use lowercase boldface for vectors, uppercase italics for matrices, and uppercase calligraphic font for sets. All vectors are column vectors. The symbol “” denotes vertical concatenation of vectors or matrices. In particular, if and then . We use the -norm for vectors and the Frobenius norm for matrices unless stated otherwise. The set is denoted as .
2. Shor relaxation of a QCQP
A quadratically constrained quadratic program (QCQP) has a well studied SDP relaxation. We review this relaxation in this section.
Let denote the space of real symmetric matrices of size . Given matrices and a vector , consider the homogeneous QCQP:
[TABLE]
By homogeneous we mean that the cost function and constraints do not involve linear terms in . We can rephrase the problem in terms of the matrix . Observe that is positive semidefinite (PSD), denoted , and has rank one. Also note that , where denotes the trace inner product. Then the QCQP is equivalent to
[TABLE]
The Shor SDP relaxation of (P) is obtained by dropping the rank-one constraint:
[TABLE]
The conic dual of the above problem is the following SDP:
[TABLE]
It can be shown that (D) is also the Lagrangian dual of (P). The optimal values of the above problems always satisfy . In particular, for any feasible dual vector we have that is a lower bound for .
The SDP relaxation is exact if and the minimizer of (P*) has rank one. In this case, the minimizer of (P) can be obtained by factorizing , and the corresponding dual solution gives an optimality certificate for .
We proceed to the inhomogeneous case. Given , , , consider
[TABLE]
We can homogenize the problem by introducing an auxiliary variable . Let be obtained by prepending to . We claim that (4) is equivalent to the following homogeneous QCQP:
[TABLE]
where and . Note that by the first constraint. We may further assume that , as the problem is invariant under sign changes. Plugging in in (5) leads to (4), as we claimed.
3. Derivation of our SDP relaxation
In this section we present our SDP relaxation for (1). Our strategy consists in phrasing (1) as a QCQP, and considering the associated Shor relaxation.
3.1. The naive attempt
For convenience, we first center the problem at zero with the change of variables . Denoting , problem (1) becomes
[TABLE]
We can rewrite the above as a QCQP by introducing a variable representing a unit vector in the left kernel of . This leads to the inhomogeneous QCQP:
[TABLE]
The above is known as the kernel representation [27]. Although we can derive an SDP relaxation from (6), it does not reveal any information about the original problem.
Lemma 3.1**.**
The optimal value of the Shor relaxation of (6) is zero for any .
Proof.
The homogenized problem (see (5)) has variables , , . We may partition the PSD matrix of the Shor relaxation accordingly:
[TABLE]
We omitted the lower triangular blocks, as is symmetric. The Shor relaxation is
[TABLE]
where is the constraint matrix associated to the -th entry of . We may partition the matrix in blocks as above. Observe that all its diagonal blocks are zero since the product is bilinear. Let be the matrix with blocks
[TABLE]
Note that is feasible and . So the optimal value of the SDP is zero. ∎
3.2. Lifting the problem
Our first QCQP formulation, given in (6), did not succeed. We will derive another equivalent QCQP in a higher dimensional space, which gives rise to a better SDP relaxation. Let be obtained by prepending a one to , and consider the Kronecker products
[TABLE]
We will restate problem (6) in terms of the vector .
Since is an affine map, it can be written as
[TABLE]
where only the first matrix depends on . Let and such that
[TABLE]
so the horizontal concatenation of equals the vertical concatenation of . Observe that and . So we can rewrite (6) as
[TABLE]
Example 3.2**.**
Let , , and let , . Then
[TABLE]
Problem (8) is not a QCQP as it has linear equations and a Kronecker product constraint. We may convert the linear equations into quadratics in a simple way. Let be the canonical basis of . Then
[TABLE]
Hence, (8) is equivalent to
[TABLE]
The Kronecker product condition means that is the vectorization of an matrix of rank one (). Recall that a matrix is rank-one if and only if all its minors vanish. If we think of as a matrix, we may index its entries by pairs with , , so that . Then
[TABLE]
Injecting the above equations in (9) lead to a QCQP formulation.
3.3. Block symmetry
The quadratic equations in (10) have a nice interpretation in terms of the matrix . They say that this matrix has repeated entries: for . The pattern of repeated entries has a very special structure. We say that a matrix is -block symmetric if when divided into blocks we have that all blocks are symmetric ( must divide ). Then
[TABLE]
Example 3.3**.**
Let , , . We may index the entries of as . Reshaping as a matrix we get a single minor constraint: . The matrix can be divided into four blocks:
[TABLE]
The minor constraint ensures that is symmetric, and hence .
By the above discussion, we can rephrase (9) as
[TABLE]
The above is a QCQP, so we can derive an associated SDP relaxation. But before that we will make one last transformation that will simplify the analysis in the next sections.
There is a simple way turn an arbitrary matrix into block symmetric form. We first symmetrize by taking the average , and then symmetrize each of its blocks. We let be the result of this block symmetrization.
Since , it follows that . By substituting in (11) we obtain our final QCQP formulation:
[TABLE]
Example 3.4**.**
Retake the case from Example 3.2. There is a matrix for each and . The first and the last are:
[TABLE]
3.4. SDP relaxation
Problem (12) is a homogeneous QCQP. Its Shor relaxation is:
[TABLE]
The above pair of SDPs are our proposed convex relaxation for (1). Both SDPs have the same optimal value since the dual problem is strictly feasible (Slater’s condition). Hence, the following relations hold:
[TABLE]
Let be the optimal solutions of the primal/dual SDPs. The relaxation is exact if of (13a) has rank one. In this case the optimal values from above are all equal, the minimizer of (1) can be recovered from , and the dual variables give an optimality certificate for . Theorem 1.1 states that the relaxation is always exact in the low noise regime.
Remark 3.5* (Weighted seminorms).*
Consider the variant of (1) that uses the weighted -(semi)norm , with . Our SDP relaxation can be easily adapted to this case. In (13a) we simply replace the objective function by . Theorem 1.1 remains valid for weighted -norms (i.e., ).
Remark 3.6* (Missing data).*
Sometimes in practice the data vector has missing entries [28, §5.1], e.g., due to data corruption. We can model this situation by considering a diagonal weight matrix, where is zero if is missing and one otherwise. The proof of Theorem 1.1 is no longer valid for seminorms. Nonetheless, we will see in Example 6.2 that the SDP relaxation might still behave well in practice.
4. The role of non-tangentiality
In the previous section we rephrased problem (1) as the QCQP (12), which involves variables , with , . Let be the list of all quadratic constraints in (12) and let denote the feasible set:
[TABLE]
Let be a fixed parameter satisfying Assumption 1. This is the only parameter we consider in this section, so we will omit the -subindices to simplify the notation. Let be the optimal solution of the QCQP. Since then the optimal value must be zero, and hence has the form . We will show in this section that is a smooth manifold nearby . The non-tangentiality condition is the main ingredient behind this.
4.1. Abadie constraint qualification
We introduce here a more refined notion of smoothness, which takes into account not only the geometry, but also the algebra.
Let continuously differentiable and . The Abadie constraint qualification (ACQ) holds at a point if the set is a smooth manifold nearby and the tangent space at is given by .
Example 4.1**.**
Let The set consists of the multiples of . Since and then ACQ does not hold at . Nonetheless, there is an alternative description of for which ACQ holds, namely .
The above example illustrates that the ACQ property depends on the algebraic representation of the set . Informally, ACQ states that the local behavior of is well described by the linearization of . This condition (and similar constraint qualifications) play a pivotal role in optimization [3, §5.1]. For instance, constraint qualifications are needed to ensure the existence of Lagrange multipliers at a local minimum, or to guarantee convergence for local optimization methods.
We may also define the ACQ property for parametric sets. Let and . We say that ACQ holds at if is a smooth manifold nearby and . Note that this is equivalent to saying that the implicit set satisfies ACQ at .
Assume now that is a polynomial map. In such a case is a real algebraic variety (the zero set of polynomial equations). The ACQ property is called regularity or nonsingularity in algebraic geometry [14, §16.6]. Recall that the ideal generated by consists of all sums where each is a polynomial. Changing the defining equations preserves the ACQ property whenever the generated ideal is unchanged. More generally, the following lemma allows us to understand how behaves under a polynomial map.
Lemma 4.2**.**
Let , , be polynomial maps. Assume that the ideals generated by and satisfy . Let with , and let . Then .
Proof.
Let be the variables in . Since the entries of lie in the ideal , there is an matrix with polynomial entries such that
[TABLE]
By taking the gradient of the above equation and then evaluating at , we get
[TABLE]
So if then , as wanted. ∎
4.2. Non-tangentiality and ACQ
Recall the non-tangentiality condition from (3). This is a geometric condition, independent of the algebraic description. As shown next, non-tangentiality arises naturally when establishing ACQ for an intersection.
Lemma 4.3**.**
Let , , . Let such that ACQ holds for both at . Then ACQ holds for at if and only if meet non-tangentially at .
Proof.
As satisfies ACQ for , then . The lemma follows from:
[TABLE]
Remark 4.4*.*
The lemma also holds if we use parametric descriptions for either .
Example 4.5**.**
The set from Example 4.1 is the intersection of and . The determinantal variety satisfies ACQ away from the origin, and the linear space satisfies ACQ everywhere. As ACQ does not hold for , then they do not intersect non-tangentially.
4.3. ACQ for our problem
Our main focus is in the variety from (14). We goal is to prove the following lemma.
Lemma 4.6**.**
Under Assumption 1, we have that satisfies ACQ at .
We first analyze a simpler case. Consider the feasible set of the QCQP (6):
[TABLE]
Lemma 4.3 allows us to show that ACQ holds for at .
Lemma 4.7**.**
Under Assumption 1, we have that satisfies ACQ at .
Proof.
Consider the sets
[TABLE]
The set is the intersection of and , using a kernel description for and a parametrization for . We claim that both satisfy ACQ everywhere. This is straightforward for , as it is the image of an affine map. We proceed to . The Jacobian matrix of the defining equations is of size . The Jacobian has full row rank on since . By the implicit function theorem, is a manifold of codimension . So the codimension of the tangent space is equal to the rank of the Jacobian. It follows that ACQ holds everywhere on .
We proceed to show that , meet non-tangentially at , where . This would conclude the proof because of Lemma 4.3. By Assumption 1, the manifold of rank matrices meets non-tangentially at , so that Observe that each has exactly two unit vectors in its left kernel. In a neighborhood of , we can find a smooth map , where is a unit vector in the left kernel of , such that . The map gives a local diffeomorphism of and (nearby ). Moreover, also gives a local diffeomorphism of and . Hence is a smooth manifold nearby and
[TABLE]
It follows that meet non-tangentially at , as wanted. ∎
We are ready to prove Lemma 4.6. The intuitive idea is to show that the varieties and are isomorphic, so any ACQ point of corresponds to an ACQ point of .
Proof of Lemma 4.6.
Note that is the image of under the map
[TABLE]
Observe that the map is a closed embedding (it is variant of the projective Segre embedding). In particular, the left inverse of is the map
[TABLE]
where is obtained by reshaping as a matrix. By Lemma 4.7, the ACQ property is satisfied for . In particular, is a manifold nearby . The restriction of to this manifold is also an embedding. It follows that a is manifold nearby , and the differential gives an isomorphism .
Recall that are defined by the polynomial equations . It is easy to see that the polynomials in lie in the ideal . By Lemma 4.2, the differential , which is injective, gives an inclusion . Let
[TABLE]
which agrees with when . One can check that the polynomials in lie in the ideal . By Lemma 4.2, we also have an inclusion , and hence . Using that ACQ holds for , we get
[TABLE]
We conclude that ACQ holds for . ∎
5. Exactness under low noise
In this section we prove our main result, Theorem 1.1. To do so, we rely on the general framework to analyze the stability of SDP relaxations from [10].
5.1. Stability of SDP relaxations
Let be a parametric family of quadratic equations in variables , , , parametrized by a vector . Consider the family of QCQPs:
[TABLE]
We assume that there is a parameter such that the optimal value of () is zero. It can be shown that for such a parameter the Shor SDP relaxation is always exact. The next theorem, proved in [10], gives conditions under which the Shor relaxation of () is stable nearby , in the sense that it is also exact for any sufficiently close to .
Theorem 5.1** (SDP stability).**
Consider the parametric family (). Let be a parameter such that the minimizer of () has the form . Assume that the following conditions hold:
- •
(constraint qualification)* ACQ holds for at .*
- •
(smoothness)* The set is a smooth manifold nearby .*
- •
(not a branch point)* The right kernel of is trivial.*
- •
(restricted Slater)* There exists such that , where*
[TABLE]
Then the Shor relaxation of () is exact whenever is close enough to .
Proof.
This is a special instance of [10, Thm. 5.1] in which the branch point condition is specialized to the case of nearest point problems (see [10, Ex. 5.12]). ∎
Remark 5.2*.*
The above theorem remains valid for weighted -norms.
The QCQP (12) is a special case of (), so it suffices to verify that the conditions from Theorem 5.1 are satisfied. The last condition (restricted Slater) is particularly challenging, as discussed in [10]. The linear space from (16) consists of the Lagrange multipliers at . So we need to prove the existence of a Lagrange multiplier vector that satisfies a certain PSD condition, but the theorem does not mention how to find it. We address this issue by introducing a new condition which is easier to check and implies restricted Slater.
We proceed to explain this new condition. Assume that some terms in the summation of vanish. Let be a partition of the equations into two groups, where the second group corresponds to the vanishing terms, i.e.,
[TABLE]
We may divide the multipliers accordingly: with , . Let be the restriction of to , i.e.,
[TABLE]
Note that for because of (17). We can similarly partition the Jacobian in the form , with , .
The next lemma provides a sufficient condition for restricted Slater. We use the notation for the symmetric part of a matrix .
Lemma 5.3**.**
Let as above. Consider the linear spaces
[TABLE]
and let be the orthogonal complement of . Assume that and . Then the restricted Slater condition holds.
Proof.
We view as a linear map with domain . Let be the projection . Note that by construction. Let be the adjoint operators of . We claim that
[TABLE]
By taking the orthogonal complement of the first equation, we conclude that
[TABLE]
The above relation implies the restricted Slater condition, Hence, it suffices to show (19).
Since , the first equation in (19) is a consequence of the following identity:
[TABLE]
This identity follows by applying the rank-nullity theorem () to the linear maps , , and .
It remains to show that Let be the inclusion and be the projection , so that . Notice that since is surjective. Also note that since is defined by the equation , and hence . As , the identity (20) gives
[TABLE]
Observe that . It follows that if and only if there is some such that and , as wanted. ∎
We proceed to prove Theorem 1.1. In our case, the constraints are given in (14). The variables are , , . We need to verify that the four conditions from Theorem 5.1 are satisfied. We already showed in Lemma 4.6 that ACQ holds. We proceed to analyze the remaining three conditions.
5.2. Smoothness
Let be the variety from (14) and let be the variety from (15). Consider the sets
[TABLE]
We need to show that is a smooth manifold nearby . It is easier to first analyze . For the nominal parameter , we showed in Lemma 4.7 that satisfies ACQ at . Hence, is a smooth manifold nearby . Since , then is an affine translation of . More precisely, , where . It follows that is a smooth manifold nearby .
Note that is the image of under the map , . In fact, this map gives an isomorphism of varieties, as showed in the proof of Lemma 4.6. Therefore, is the image of under the map . It follows that is a smooth manifold nearby .
5.3. Branch point
From now on we will only focus on the nominal parameter , and hence we will no longer write the -subindices. In order to show that is not a branch point, we will first provide an explicit formula for the Jacobian .
We will identify with the tensor product . In particular, use the following bases for , , and :
[TABLE]
Recall from (7) that the map can be written in terms of matrices . For , let denote the -th columns of . Observe that . In terms of the bases from above, we have
[TABLE]
The next lemma gives a formula for (some of) the rows of .
Lemma 5.4**.**
Let , and consider the constraints from (14). Then
[TABLE]
and for any and we have
[TABLE]
where stands for the symmetric part.
Proof.
Since , then . To compute we make use of the following identity:
[TABLE]
The above identity can be checked by straightforward manipulation. It follows that
[TABLE]
Reducing the above expression we get the formula in the lemma. ∎
In order to make the following proofs more explicit, we will introduce a new assumption, which is weaker than Assumption 1.
Assumption 2**.**
The vectors as in (21) span .
Remark 5.5* (Assum1 Assum2).*
Recall that . Under Assumption 1, this matrix has rank . Since , Assumption 2 holds.
We now verify the branch point condition.
Lemma 5.6**.**
Under Assumption 2, the right kernel of is trivial.
Proof.
Let be in the right kernel. Then and for all . By Lemma 5.4, we have and for all . Since , some is nonzero. Then is orthogonal to each of , so it must be zero. ∎
5.4. Restricted Slater
In order to use Lemma 5.3, we partition the equations into the two groups shown in Table 1. Group corresponds to the first column and group to the remaining three columns. The table also gives the multipliers associated to each constraint. The complete multiplier vector consists of all .
Let us see that group satisfies (17). It can be checked that for the equations in the second and third column. As for the fourth column, let and let us show that . Since , then . Note that all the equations, except for , are homogeneous, so they satisfy . It follows that , and hence .
We proceed to compute the linear map from (18). To do so, let
[TABLE]
Group consists of equations: for and . It follows from Lemma 5.4 that . We may view the multiplier vector as an element of . Then the linear map satisfies
[TABLE]
We now compute the linear space .
Lemma 5.7**.**
Under Assumption 2, then .
Proof.
Let , . We will first show that
[TABLE]
The second containment is clear. We proceed to prove the first one. Given , we can write it in the form
[TABLE]
It suffices to show that and . Since , then and also for with . Lemma 5.4 gives formulas for these gradients. The equation says that . Assume by contradiction that . The equations imply
[TABLE]
Since then some is nonzero. Also note that some since . Similarly, some since . This is a contradiction. Hence, (23) holds.
It remains to show that , , where . The rows of are , where , . An explicit formula is given in Lemma 5.4. Let , and let us see that . The -th entry of is
[TABLE]
where we used that . Then,
[TABLE]
Consider now . The -th entry of is
[TABLE]
and thus ∎
We are ready to show the restricted Slater condition.
Lemma 5.8**.**
Under Assumption 2, then the conditions from Lemma 5.3 are satisfied. Hence, the restricted Slater condition holds.
Proof.
The formula for from (22) implies that
[TABLE]
The first equation implies that . Since by Lemma 5.7, then the second equation implies that . So the assumptions of Lemma 5.3 hold. ∎
6. Applications
We implemented our SDP relaxation using CVX [16]. The code is available in the personal website of the author: http://www.mit.edu/~diegcif/. Although not the main focus of this paper, our implementation allows problems with complex matrices, weighted norms, and missing entries. We proceed to illustrate the performance of our relaxation in various applications. The code needed to replicate the experiments is also available in the personal website of the author. We used the SDP solver Mosek with the default parameters for all experiments.
6.1. Hankel structure
Consider the problem of computing the nearest rank deficient Hankel matrix. This structure is prevalent in applications from systems theory and control, due to the fact that a superposition of exponential signals corresponds to an infinite Hankel matrix of rank . By suitably choosing a finite submatrix we obtain a rank deficient Hankel matrix. Some concrete applications are approximate realization, system identification, noisy deconvolution, and stochastic realization [25, 28, 27].
For our first experiment we consider the Hankel approximation problem in an abstract setting (without referencing any application), and illustrate the “global” behavior of our method. By global we mean that the data vector is sampled from the uniform distribution on the unit sphere. The left of Table 2 illustrates the performance of our SDP relaxation for small values of . The numbers in the table indicate the percentage of experiments for which SDP solved the problem exactly (returning an optimality certificate). For each we used 2000 uniformly random vectors . Observe that the SDP was exact for all instances with .
For comparison, the right of Table 2 illustrates the performance of local optimization using the software slra [29]. We use the Levenberg-Marquardt algorithm (LMA) as suggested in the manual, using the initialization scheme from [39] (gives better results than the default initialization). As expected, the accuracy of local optimization is lower than that of SDP. Moreover, even when the local method is exact, it is hard to certify that it indeed converged to the global optimal. In particular, we can only verify that the local method succeeded for the instances for which SDP also succeeds.
The fraction of instances for which SDP is exact should diminish as . Nevertheless, SDP always behaves well in the low noise regime, as we showed in Theorem 1.1. The following example illustrates this behavior.
Example 6.1** (Approximate realization).**
Consider the discrete linear time invariant (LTI) system with transfer function . Let be the impulse response of the system, and let be the Hankel matrix with the first entries of . The rank of is two (the order of the system). Assume now that the signal is corrupted by a Gaussian noise, and hence is full rank. The approximate realization problem consists in finding the nearest vector such that has rank two. Table 3 shows the performance of SDP for and for different standard deviations of the noise. The numbers are the percentage of experiments which are solved exactly over repetitions. Observe that SDP solves all instances below a certain noise level, as predicted by Theorem 1.1. For comparison we also show the performance of the local optimization methods LMA and BFGS. We use again the initialization from [39]. As before, we can only verify that the local methods are exact for instances for which SDP is exact.
Example 6.2** (Missing data).**
Retake the approximate realization problem from above, and assume that some entries of the vector are missing. As mentioned in Remark 3.6, we might apply our SDP by choosing a degenerate weight matrix. Table 4 shows the performance of the SDP for , considering two scenarios for the missing entries. In the left of Table 4 we delete for all , i.e., for . Our relaxation is almost as accurate as when all entries are observed (Table 3). In the right of Table 4 we delete most of the entries, only keeping when . This setting is quite complicated, due to the long gaps in between observed entries. Nonetheless, our relaxation still behaves remarkably well. For comparison, we also show the performance of LMA and BFGS. For the initialization we first complete the data using a cubic spline interpolation (linear interpolation did not work well) and then use the algorithm from [39]
6.2. Approximate GCD
Let and let be univariate polynomials of degrees , The degree- approximate GCD problem is:
[TABLE]
where denotes the -norm of the coefficient vector. The approximate GCD problem can be restated as in (1). Indeed, it is known that if and only if is rank deficient, where is the Sylvester matrix which is filled with the coefficients of ; see e.g., [22].
Example 6.3**.**
Consider the polynomials
[TABLE]
We normalize the coefficients of and then corrupt them with Gaussian noise. Table 5 illustrates the performance of SDP in computing the degree- approximate GCD. The numbers indicate the percentage of successful experiments over repetitions. Note that SDP solves all instances below a certain noise level. For comparison, we also show the performance of LMA and BFGS using the default initialization of slra.
6.3. Multi-view geometry
The pinhole camera model [17, §6] from computer vision represents 3D points by vectors , 2D points (images) by vectors , and cameras by matrices . The image of the 3D point under the camera is the following 2D point:
[TABLE]
Because of the above equation, problems from multi-view geometry often involve ratios of linear terms. The fractional programming problem from (2) captures several of these problems [20], including the triangulation problem, the camera resection problem, and the homography estimation problem.
Example 6.4** (Triangulation).**
Given cameras and noisy images of an unknown 3D point , the triangulation problem is
[TABLE]
This is a special case of (2) with and . Aholt, Agarwal, and Thomas proposed an SDP relaxation for this problem in [1], and they showed that it is exact under low noise when the camera centers are not coplanar. The SDP from [1] is smaller than ours: its PSD matrix is , as opposed to for ours. But, as we will see, ours is more precise, succeeding in almost all instances we considered.
Table 6 compares both SDP relaxations, indicating the percentage of experiments that are solved exactly (the solution is rank-one) over repetitions. The left of Table 6 considers the synthetic data set described in [1, 31]: the cameras are placed uniformly at random on the sphere of radius two pointing to the origin, and the points are generated uniformly at random inside the unit cube. For the right of Table 6 the cameras are placed uniformly on the line segment —, and the points are as before. In both cases the image size is approximately units. The second configuration is very problematic for the SDP from [1], since the camera centers are coplanar. On the other hand, our SDP (13) behaves equally well.
Example 6.5** (Resectioning).**
Given points and noisy images under an unknown projective camera , the resectioning problem is
[TABLE]
This is a special case of (2) with and . Table 7 illustrates the performance of our SDP (13) on the synthetic data set from [31]: the cameras are randomly placed on the sphere of radius two, and the points inside the unit cube. Our SDP is exact for most instances below a noise level of . For comparison, we refer to [31, Fig 2], where it is observed that their method cannot certify global optimality for any instance above a noise level of .
Kahl and Henrion proposed an SDP relaxation for this problem in [21], though without any guarantees. Table 7 also shows the performance of their SDP (the Schur formulation of order ). We use the heuristic method explained in [21] of including a small multiple of the trace in the objective function (). This promotes a rank-one solution, at the expense of a loss in optimality. As shown in Table 7, the solution still had higher rank in most of the experiments performed, even with this heuristic.
Remark 6.6*.*
The homography estimation problem is also a special instance of (2), so we may use our SDP relaxation. Another relaxation was proposed in [21], but with no guarantees. A related problem is the estimation of the essential matrix of two views, for which an SDP relaxation that is exact under low noise was recently proposed [43].
7. Future work
Theory
Our SDP always recovers the global minimizer of (1) in the low noise regime, as shown in Theorem 1.1. More generally, the theorem applies to weighted -norms of the form with . We can also use our relaxation for seminorms defined by a weight matrix , but our proof does not generalize. Nonetheless, the experiments in Table 4 suggest that Theorem 1.1 might still hold in some cases. Investigating this behavior is an interesting open problem.
Our experiments reveal that our SDP relaxation is quite resilient to noise, outperforming state of the art methods. It would be desirable to establish quantitative lower bounds on the noise tolerance. The stability theory developed in [10] includes lower bounds, but preliminary calculations indicate they are overly pessimistic. Deriving practical lower bounds on the noise tolerance is left for future work.
It is natural to ask whether our SDP relaxation can be extended to the general SLRA problem i.e., when satisfies a rank constraint. We may derive an SDP relaxation for SLRA following similar steps as in Section 3 (start with the kernel formulation, rewrite in terms of the Kronecker product, and derive the Shor relaxation). But this relaxation is never exact when . It is an interesting problem to derive a convex relaxation for SLRA with provable guarantees.
Computation
SDPs can be solved in polynomial time with interior point methods, but their high memory requirements limit their use in large scale problems. Scalability is a pervasive issue phased by SDPs, but it has been overcome in several other applications, e.g., [41, 36, 42, 12, 26]. Developing an scalable solver for our SDP relaxation is a very important open problem.
There are several known techniques for improving efficiency, see the survey paper [26]. Taking advantage of problem structure is critical (e.g., sparsity [40], symmetry [34], geometry [11]). One often also relies on first order methods instead of interior point methods. Two concrete examples are the Burer-Monteiro method [6] and the sketching technique from [42]. The PSD matrix in our relaxation has size , which can be problematic when both are large. This is due to our Kronecker product formulation (this corresponds to the second level of the SOS hierarchy). Some ways to tackle this issue are proposed in [12, 33].
We point out that our SDP relaxation might also be used in conjunction with local optimization methods in order to provide a posteriori verification of candidate solutions. This approach has been explored in e.g. [41, 18], and is often significantly more efficient than solving the SDP by itself.
Acknowledgments
The author thanks Sameer Agarwal, Pablo Parrilo, Bernd Sturmfels, Rekha Thomas, and André Uschmajew for helpful discussions and comments. This work was partially done while the author was in the Max Planck Institute in Leipzig.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] C. Aholt, S. Agarwal, and R. Thomas. A QCQP approach to triangulation. In ECCV (1) , volume 7572 of Lect. Notes Comput. Sci. , pages 654–667. Springer, 2012.
- 2[2] F. Andersson and M. Carlsson. Alternating projections on nontangential manifolds. Constr. Approximation , 38(3):489–525, 2013.
- 3[3] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty. Nonlinear programming: theory and algorithms . John Wiley & Sons, 2013.
- 4[4] A. Ben-Tal and A. Nemirovski. Lectures on modern convex optimization: analysis, algorithms, and engineering applications , volume 2. SIAM, 2001.
- 5[5] G. Blekherman, P. Parrilo, and R. Thomas. Semidefinite optimization and convex algebraic geometry , volume 13. MOS-SIAM Series Optim., 2013.
- 6[6] S. Burer and R. D. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Math. Program. , 95(2):329–357, 2003.
- 7[7] E. Candès and B. Recht. Exact matrix completion via convex optimization. Found. Comput. Math. , 9(6):717, 2009.
- 8[8] M. Carlsson. On convex envelopes and regularization of non-convex functionals without moving global minima. J. Optim. Theory Appl. , 183(1):66–84, 2019.
