Outer Spherical Approximations

Localization > Back | Outer Spherical Approximations
  • Outer spherical approximation via gridding

  • Outer spherical approximation via Lagrange relaxation

  • Improved spherical approximation via affine Lagrange multipliers

Let us consider a candidate sphere mathbf{S}_0, of center x_0 and radius R. We formulate the condition that the sphere mathbf{S}_0 contains the intersection of all the spheres of center mathbf{S}_i, as follows:

 (ast) ~~~:~~~ |x-x_0|_2^2 le R_0^2 mbox{ for every } x mbox{ such that } |x-x_i|_2^2 le R_i^2, ;; i=1,ldots,m.

This conditions appears to be hard to check.

Outer spherical approximation via gridding

A first approach, which works well only in moderate dimensions (2D or 3D), simply entails gridding the boundary of the intersection. In 2D, we can parametrize the boundary explicitly, as a curve, using an angular parameter. For each angular direction theta in [0,2pi], we can easily find the point that is located on the boundary of the intersection, in the direction given by theta: we simply maximize t such that the point (tcostheta,tsin theta) is inside everyone of the spheres. (There is an explicit formula for the maximal value.)

One we have computed N points on the boundary, x^{(k)}, k=1,ldots,N, we simply solve the SOCP

 min_{x_0,R_0} : R_0 ~:~ R_0 ge |x_0 - x^{(k)}|_2, ;; k=1,ldots,N.

The gridding approach suffers from the fact that we need to grid finely to be able to guarantee that the spherical approximation we are computed indeed contains all the points on the intersection. On the other hand, too fine a grid results in a large number of constraints, which in turn adversely impacts the time needed to solve the above problem.

alt text 

Outer spherical approximation to the intersection via gridding. This provides an estimated point (the center of the outer shpere), with an estimate of the uncertainty around it. We have used only 13 points on the boundary to enforce the constraint that the outer sphere contains the intersection. As a result, our approximation is not really an outer approximation, as it does not contain all the intersection points.

alt text 

Outer spherical approximation to the intersection via gridding. Here we have used 63 points on the boundary. The approximation can now be termed an outer approximation as it (almost!) does not leave out any point in the intersection.

Outer spherical approximation via Lagrange duality

We now develop an alternate approach which relies on weak duality concepts.

A sufficient condition

A sufficient condition for the above condition (ast) to hold is that there exist a non-negative vector y in mathbf{R}_+ such that

 mbox{ for every } x ~:~ |x-x_0|_2^2 le R_0^2 + sum_{i=1}^m y_i left( |x-x_i|_2^2 - R_i^2 right).

Indeed, it is easily checked that if a point x belongs to all the spheres, then indeed the sum above is less or equal than 0 when y ge 0.

Minimizing the radius R_0 subject to the condition that the above holds for some y ge 0 will result in a conservative (larger) estimate of the amount of uncertainty.

SOCP formulation

Alternatively we can express the above condition as:

 R_0^2 ge F(x_0,y) := sum_{i=1}^m y_i R_i^2 +max_x : left(|x-x_0|_2^2 - sum_{i=1}^m y_i |x-x_i|_2^2 right) .

The problem of minimizing the radius R_0 subject to the above condition can thus be written as

 (R_0^2)^ast : = min_{x_0, y ge 0} :  F(x_0,y).

It turns out that we can express the function F(x_0,y) in closed-form, as proven here:

 F(x_0,y) = left{ begin{array}{ll}  x_0^Tx_0 + y^T z + frac{|Xy-x_0|_2^2}{S-1} & mbox{if } S:= sum_{i=1}^m y_i > 1,  x_0^Tx_0 + y^T z & mbox{if } sum_{i=1}^m y_i = 1, ;; x_0 = Xy ,  +infty & mbox{otherwise,} end{array} right.

where for notational convenience we use the matrix X := (x_1,ldots x_m) and the vector z with components z_i := R_i^2-x_i^Tx_i, i=1,ldots,m. (Remember that our measurement consistency condition holds, so that z ge 0.)

We are led to consider the sub-problem of minimizing F(x_0,y) over x_0, with y fixed. If S:=sum_{i=1}^m y_i = 1, then we must have x_0 = Xy. If S> 1, x_0 must solve

 min_{x_0} : x_0^Tx_0 + frac{|Xy-x_0|_2^2}{S - 1}.

Again, this a convex quadratic problem, with a (unique) solution obtained by taking derivatives. A minimizer is x_0(y), where

 x_0(y) := frac{1}{sum_{i=1}^m y_i}{Xy},

that is, x_0 is the weighted average, with weights given by yge 0, of the points x_i. Note that the expression remains valid when sum_i y_i = 1, since then we must have x_0 = Xy.

Plugging the above value of x_0 in the problem leads to a problem in variable y only:

 (R_0^2)^ast = min_{y} : y^Tz + frac{1}{sum_{i=1}^m y_i}|Xy|_2^2 ~:~ y ge 0, ;; sum_{i=1}^m y_i ge 1.

The above problem can be solved as an SOCP with rotated second-order cone constraints:

 min_{y,alpha} : y^Tz + alpha~:~  alpha ( sum_{i=1}^m y_i ) ge |Xy|_2^2, ;; y ge 0, ;; sum_{i=1}^m y_i ge 1.

A simpler formulation (consistent measurements)

We can obtain a simpler formulation that involves QP only. We now use our measurement consistency assumption, which translates as z ge 0. First we replace the variable y with two new variables S,p, and a new constraint. The new variables S,p are defined as

 S := sum_{i=1}^m y_i, ;; p_i = frac{y_i}{S}, ;; i=1,ldots,m,

and the new constraint is sum_{i=1}^m p_i = 1. We obtain the new problem

 (R_0^2)^ast = min_{S,p} : S left( p^Tz + |Xp|_2^2 right) ~:~ p ge 0, ;; sum_{i=1}^m p_i= 1, ;; S ge 1.

Since z ge 0, the term inside the parentheses is non-negative, and thus S = 1 is optimal. Hence the problem reduces to a QP:

 (R_0^2)^ast = min_{p} : p^Tz + |Xp|_2^2 ~:~ p ge 0, ;; sum_{i=1}^m p_i= 1 .

This provides the optimal radius. The optimal point is

 x_0^ast = frac{1}{S} Xy = Xp = sum_{i=1}^m x_i p_i.
alt text 

Outer spherical approximation via Lagrange duality. This provides an estimated point (the center of the outer sphere), with an estimate of the uncertainty around it. The approximation is very conservative.

An improved condition via affine duality

To improve our condition, we start from the equivalent condition for intersection containment:

 t-2x_0^Tx+x_0^Tx_0 le R_0^2 mbox{ for every } x,t mbox{ such that } t = x^Tx, ;; t-2x_i^Tx+x_i^Tx_i le R_i^2, ;; i=1,ldots,m.

We now consider a relaxed of the form

 mbox{ for every } x,t ~:~ t-2x_0^Tx+x_0^Tx_0 le R_0^2 + tau (t-x^Tx) + sum_{i=1}^m y_i(x,t) left( t-2x_i^Tx+x_i^Tx_i- R_i^2 right),

where tau ge 0, and the y_i's are now non-negative functions of (x,t). (Before, we chose them to have constant values.) The above is still hard to check in general, but becomes easy if we make the assumption that the functions (x,t) rightarrow y_i(x,t) are affine.

We proceed by taking the vector y(x,t) := (y_i(x,t))_{i=1}^m to be of the form

 y(x,t) = y + Y^Tx + t nu,

where y in mbox{bf R}^m, Y in mbox{bf R}^{2 times m} and nu in mbox{bf R}^m will be variables. (The case before will be recovered upon setting Y=0, nu = 0.)

The condition above becomes a single quadratic condition on (x,t):

 mbox{ for every } x,t ~:~ t-2x_0^Tx+x_0^Tx_0 le R_0^2 + tau (t-x^Tx) + sum_{i=1}^m (y + Y^Tx + t nu)_i left( t - 2x^Tx_i +|x_i|_2^2 - R_i^2 right).

The above can be equivalently expressed as a positive semi-definiteness condition on a symmetric matrix that is affine in y,Y,nu:

 left( begin{array}{ccc} nu^Tmathbf{1} & -nu^TX^T+frac{1}{2}mathbf{1}^TY^T & frac{1}{2}(mathbf{1}^Ty - nu^Tz -tau-1) -Xnu-frac{1}{2}Ymathbf{1}  &tau I - YX^T-XY^T & x_0 - Xy - frac{1}{2}Yzfrac{1}{2}(mathbf{1}^Ty - nu^Tz -tau-1) & (x_0 - Xy - frac{1}{2}Yz)^T & R_0^2 - y^Tz-eta end{array}right) succeq 0, ;; eta ge x_0^Tx_0.

It remains to express the condition that the affine functions y_i(x,t) should be non-negative on the intersection. Precisely, we seek to enforce that for every i in {1,ldots,m}, the condition

 y_i + Y_i^Tx + nu_i x^Tx ge 0 mbox{ for every } x mbox{ such that }  x^Tx-2x_j^Tx+x_j^Tx_j le R_j^2, ;; j=1,ldots,m

holds (here, Y_i stands for the i-th column of matrix Y). Now, we apply direct Lagrange duality. A sufficient condition for the above to hold is that there exist non-negative numbers N_{ij}, j=1,ldots,m, such that

 mbox{ for every } x ~:~ y_i + Y_i^Tx + nu_i x^Tx ge sum_{j=1}^m N_{ij} (R_j^2 - x^Tx+2x_j^Tx-x_j^Tx_j ).

Again, the above is equivalent to a single linear matrix inequality in the variables y,Y,nu and N:

 left( begin{array}{cc} nu_i - sum_{j=1}^m N_{ij} & frac{1}{2} Y_i - sum_{j=1}^m N_{ij}x_j  (frac{1}{2} Y_i - sum_{j=1}^m N_{ij}x_j)^T & y_i - sum_{j=1}^m N_{ij} z_j end{array}right) succeq 0.

Putting this together, we obtain the SDP in variables x_0,R_0^2,y,Y,nu,eta,N:

 min : R_0^2 ~:~ begin{array}[t]{l} left( begin{array}{ccc} nu^Tmathbf{1} & -nu^TX^T+frac{1}{2}mathbf{1}^TY^T & frac{1}{2}(mathbf{1}^Ty - nu^Tz -tau-1) -Xnu-frac{1}{2}Ymathbf{1}  &tau I - YX^T-XY^T & x_0 - Xy - frac{1}{2}Yzfrac{1}{2}(mathbf{1}^Ty - nu^Tz -tau-1) & (x_0 - Xy - frac{1}{2}Yz)^T & R_0^2 - y^Tz-eta end{array}right) succeq 0, ;; eta ge x_0^Tx_0  left( begin{array}{cc} nu_i - sum_{j=1}^m N_{ij} & frac{1}{2} Y_i - sum_{j=1}^m N_{ij}x_j  (frac{1}{2} Y_i - sum_{j=1}^m N_{ij}x_j)^T & y_i - sum_{j=1}^m N_{ij} z_j end{array}right) succeq 0, ;; i=1,ldots,m, end{array}
alt text 

Improved outer spherical approximation via affine Lagrange duality. The refined approximation is now exact, in contrast with the previous approximation. The improvement is due to the fact that our relaxation involves Lagrange dual variables that are not constant anymore, but are affine in the primal variables.