Overview Motivation and Introduction Solving CMPs A heuristic - - PowerPoint PPT Presentation

Solving complex-valued ℓ0 minimization problems with constant modulus constraints

Frederic Matter Joint work with Tobias Fischer (Fraunhofer ITWM Kaiserslautern), Ganapati Hegde (TU Darmstadt), Marius Pesavento (TU Darmstadt), Marc Pfetsch (TU Darmstadt) and Andreas Tillmann (RWTH Aachen)

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 1

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 2

Introduction

We want to solve the following problem

Constant Modulus Program (CMP)

min

x∈CN x0

• s. t. b − Ax2 ≤

√ δ, |xn| ∈ {0, c}, ∀n ∈ [N],

with A ∈ CK×N, b ∈ CK , δ ∈ R, c ∈ R, [N] = {1, ... , n} and x0 denoting the number of nonzero entries in x.

◮ In the problem we could also allow additional constraints, such as linear

(in-)equalities.

◮ Throughout this talk, we assume w.l.o.g. that c = 1.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 3

Motivation

◮ Generalization of compressed sensing: search sparse solutions of an

• ptimization problem ⇒ use ℓ0-term in objective function.

◮ Note that the constant modulus constraint is nonconvex. ◮ We want to solve this problem exactly. ◮ constant modulus often appears in signal processing (e.g., phase retrieval).

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 4

MINLP–Reformulation

◮ auxiliary binary variables y = [y1, y2, ... , yN]T ∈ {0, 1}N. ◮ wn = Re[xn] and zn = Im[xn] as real and imaginary part of xn, respectively. ◮ w = [w1, w2, ... , wN]T and z = [z1, z2, ... , zN]T. ◮ AT = [a1, ... , aK ] ∈ CN×K .

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 5

MINLP–Reformulation

◮ auxiliary binary variables y = [y1, y2, ... , yN]T ∈ {0, 1}N. ◮ wn = Re[xn] and zn = Im[xn] as real and imaginary part of xn, respectively. ◮ w = [w1, w2, ... , wN]T and z = [z1, z2, ... , zN]T. ◮ AT = [a1, ... , aK ] ∈ CN×K .

min

w,z∈RN,y N

• n=1

yn s.t.

K

• k=1
• Re[bk] −
• Re[ak]T w − Im[ak]T z

2

+

• Im[bk] −
• Re[ak]T z + Im[ak]T w

2 ≤ δ,

w2

n + z2 n ≤ yn,

∀ n ∈ [N],

w2

n + z2 n ≥ yn,

∀ n ∈ [N],

yn ∈ {0, 1},

∀ n ∈ [N].

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 5

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 6

Solving the resulting MINLP

◮ Apply spatial branch-and-cut framework. ◮ Error bound (ℓ2-norm) constraint is convex, even in the complex case, as it

can be recast as a second-order cone constraint.

second-order cone: Ck =

• [x1, ... , xk, xk+1]T ∈ Rk+1 : [x1, ... , xk]T2 ≤ xk+1
• .

second-order cone constraint: Aix + bi2 ≤ cT

i x + di,

Ai ∈ R(ni −1)×n, bi ∈ Rni −1, ci ∈ Rn, di ∈ R.

◮ Upper bound parts of modulus constraints are convex quadratic constraints. ◮ Lower bound parts of modulus constraints are nonconvex quadratic

constraints.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 7

Solving the resulting MINLP

◮ Apply spatial branch-and-cut framework. ◮ Error bound (ℓ2-norm) constraint is convex, even in the complex case, as it

can be recast as a second-order cone constraint.

second-order cone: Ck =

• [x1, ... , xk, xk+1]T ∈ Rk+1 : [x1, ... , xk]T2 ≤ xk+1
• .

second-order cone constraint: Aix + bi2 ≤ cT

i x + di,

Ai ∈ R(ni −1)×n, bi ∈ Rni −1, ci ∈ Rn, di ∈ R.

◮ Upper bound parts of modulus constraints are convex quadratic constraints. ◮ Lower bound parts of modulus constraints are nonconvex quadratic

constraints.

⇒ exploit structure of modulus constraints to develop specialized branching and

propagation (and some very easy separation) to enforce modulus constraints.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 7

The LP relaxation

Another problem: poor linear relaxation.

⇒ Use additional linear inequalities that strengthen the LP relaxation.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 8

The LP relaxation

Another problem: poor linear relaxation.

⇒ Use additional linear inequalities that strengthen the LP relaxation.

zn wn zn ≤ yn zn ≥ −yn wn ≤ yn wn ≥ −yn wn + zn ≤ √ 2 yn wn − zn ≤ √ 2 yn −wn + zn ≤ √ 2 yn −wn − zn ≤ √ 2 yn

−yn ≤ wn ≤ yn, −yn ≤ zn ≤ yn,

wn + zn ≤

√

2 yn, wn − zn ≤

√

2 yn,

−wn + zn ≤ √

2 yn,

−wn − zn ≤ √

2 yn.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 8

Enforcing Modulus Constraints

If the solution ( ˆ w, ˆ z, ˆ y) of the LP relaxation violates w2

n + z2 n ≥ yn for some n ∈ [N]:

zn wn ( ˆ wn, ˆ zn)

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 9

Enforcing Modulus Constraints

If the solution ( ˆ w, ˆ z, ˆ y) of the LP relaxation violates w2

n + z2 n ≥ yn for some n ∈ [N]: ◮ If binary variable ˆ

yn fixed to zero: set ˆ wn, ˆ zn to zero as well.

zn wn ( ˆ wn, ˆ zn)

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 9

Enforcing Modulus Constraints

If the solution ( ˆ w, ˆ z, ˆ y) of the LP relaxation violates w2

n + z2 n ≥ yn for some n ∈ [N]: ◮ If binary variable ˆ

yn fixed to zero: set ˆ wn, ˆ zn to zero as well.

◮ If bounds of the continuous variables wn and zn are not yet restricted to one of

the four orthants: create four branching nodes, one for each orthant.

zn wn ( ˆ wn, ˆ zn)

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 9

Enforcing Modulus Constraints

If the solution ( ˆ w, ˆ z, ˆ y) of the LP relaxation violates w2

n + z2 n ≥ yn for some n ∈ [N]: ◮ If binary variable ˆ

yn fixed to zero: set ˆ wn, ˆ zn to zero as well.

◮ If bounds of the continuous variables wn and zn are not yet restricted to one of

the four orthants: create four branching nodes, one for each orthant.

◮ Otherwise:

• 1. Propagation,
• 2. Separation,
• 3. Branching.

zn wn ( ˆ wn, ˆ zn)

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 9

Propagation

If yn is fixed to one: l1, u1, l2, u2 current lower and upper bounds of wn and zn. l′

1 = max{l1, f(u2)},

u′

1 = min{u1, f(l2)},

l′

2 = max{l2, f(u1)},

u′

2 = min{u2, f(l1)},

where f(x) =

√

1 − x2. w2

n + z2 n ≥ yn needs to be

fulfilled. zn wn

u2 l2 u1 l1 f(u2) f(l2) f(u1) f(l1)

(l1, f(l1)) (u1, f(u1)) (f(l2), l2) (f(u2), u2)

l′

1 = f(u2)

u′

1 = u1

u′

2 = u2

l′

2 = f(u1) March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 10

Propagation

If yn is fixed to one: l1, u1, l2, u2 current lower and upper bounds of wn and zn. l′

1 = max{l1, f(u2)},

u′

1 = min{u1, f(l2)},

l′

2 = max{l2, f(u1)},

u′

2 = min{u2, f(l1)},

where f(x) =

√

1 − x2. w2

n + z2 n ≥ yn needs to be

fulfilled. zn wn

u2 l2 u1 l1 f(u2) f(l2) f(u1) f(l1)

(l1, f(l1)) (u1, f(u1)) (f(l2), l2) (f(u2), u2)

l′

1 = f(u2)

u′

1 = u1

u′

2 = u2

l′

2 = f(u1)

⇒ Use l′

1, u′ 1, l′ 2 and u′ 2 as new lower and upper bounds of wn and zn.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 10

Separation and Branching

If ˆ wn + ˆ zn < ˆ yn, add the cut wn + zn ≥ yn to the LP relaxation.

zn wn

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 11

Separation and Branching

If ˆ wn + ˆ zn < ˆ yn, add the cut wn + zn ≥ yn to the LP relaxation. If ˆ wn + ˆ zn ≥ ˆ yn, create two branching nodes with the inequalities

− ˆ

zn

′ − 1

ˆ wn

′

· wn + zn = 1,

ˆ zn

′

ˆ wn

′ − 1 · wn − zn =

ˆ zn

′

ˆ wn

′ − 1,

where ( ˆ wn, ˆ zn) is the current LP solution, and ˆ wn

′ = ˆ wn

√

ˆ wn2+ ˆ zn2 , ˆ

zn

′ = ˆ zn

√

ˆ wn2+ ˆ zn2 .

zn wn zn wn zn wn

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 11

Specialized Approach – Some remarks

◮ Enforcement of binary variables, second-order cone (SOC) constraint and the

(convex) upper bound part of the quadratic modulus constraints is prioritized

• ver the (nonconvex) lower bound part of the modulus constraints.

◮ A modulus constraint for enforcing is selected by a “most infeasible” rule:

Enforce a modulus constraint w2

¯ n + z2 ¯ n ≥ y¯ n, ¯

n ∈ [N], with largest violation

ρ(n) = ˆ

yn − ( ˆ w2

n + ˆ

z2

n).

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 12

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 13

A heuristic

Recall the problem

min x0 s. t. b − Ax2 ≤

√ δ, |xn| ∈ {0, 1}, ∀n ∈ [N].

Basic idea: fix sparsity level M, try iteratively to create an M-sparse solution by guessing and updating. If not successful, increase M by one and restart.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 14

A heuristic

Recall the problem

min x0 s. t. b − Ax2 ≤

√ δ, |xn| ∈ {0, 1}, ∀n ∈ [N].

Basic idea: fix sparsity level M, try iteratively to create an M-sparse solution by guessing and updating. If not successful, increase M by one and restart.

◮ guess solution: draw at random M indices and randomly choose their values

• n the unit circle, set the remaining values to zero.

◮ update solution: For two randomly chosen indices, one with non-zero value,

• ne with zero value, recompute their values with a least squares problem

where all but these two indices are fixed. Keep the smaller of the two new values, possibly swap the zero and non-zero index.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 14

A heuristic

Recall the problem

min x0 s. t. b − Ax2 ≤

√ δ, |xn| ∈ {0, 1}, ∀n ∈ [N].

◮ guess 1000 starting solutions, update each 1000 times. ◮ If desired error bound

√ δ is reached: break.

◮ Use solution as initial primal bound in the branch-and-cut tree.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 15

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 16

Application example: Joint Antenna Selection and Phase-only Beamforming

◮ Sensor array with N antenna elements, M phase shifters, M ≪ N. ◮ x ∈ CN is the transmit signal vector, xn = |xn|eiψn. ◮ channel matrix H = [h1, h2, ... , hK ] ∈ CN×K . ◮ Analog beamformers a = [a1, a2, ... , aM]T ∈ CM, assume |am| = c for m ∈ [M]. ◮ if nth antenna element connected to mth phase shifter: xn = am, else: xn = 0.

s

a1 a2 aM

1 2 M

Switches

x1 x2 xN

Phase shifters Antenna array user 1 user K H s s

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 17

Application example: Joint Antenna Selection and Phase-only Beamforming

◮ Sensor array with N antenna elements, M phase shifters, M ≪ N. ◮ x ∈ CN is the transmit signal vector, xn = |xn|eiψn. ◮ channel matrix H = [h1, h2, ... , hK ] ∈ CN×K . ◮ Analog beamformers a = [a1, a2, ... , aM]T ∈ CM, assume |am| = c for m ∈ [M]. ◮ if nth antenna element connected to mth phase shifter: xn = am, else: xn = 0. ◮ K single antenna users that need to

be served.

◮ desired output at users:

s = [s1, s2, ... , sK ]T.

◮ actual output at users:

ˆ s = HTx + n.

s

a1 a2 aM

1 2 M

Switches

x1 x2 xN

Phase shifters Antenna array user 1 user K H s s

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 17

Application example: Joint Antenna Selection and Phase-only Beamforming

Our goal

◮ Design the optimal phase

values of the ABF a.

◮ Assign appropriate antenna

elements to the phase shifters to serve users in the network.

◮ Minimize the number of active

antennas.

◮ Keep the root-mean-square

error between the desired and actual output at the user below

√ δ.

s

a1 a2 aM

1 2 M

Switches

x1 x2 xN

Phase shifters Antenna array user 1 user K H s s

Mathematical model

min

x∈CN x0

• s. t.
• s − HT x
• 2 ≤

√ δ, |xn| ∈ {0, c}, ∀n ∈ [N].

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 18

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 19

Implementation details

We have implemented the following parts:

◮ A reader that sets up the problem. It builds the error bound constraint as

either general quadratic constraint,

• r as SOC constraint: auxiliary variables that represent Ax − b are needed.

SOC-representation in SCIP:

• γ + n

i=1(αi (xi + βi))2 ≤ αn+1 (xn+1 + βn+1)

◮ A probdata structure to save A, b and δ. ◮ A constraint handler for the branching, propagation and the separation

methods for constant modulus constraints.

◮ The constraint handler has a propagation and enforcement method, the

branching and separation is done in the enforcement method.

◮ A heuristic plugin.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 20

Some computational results

First testset: 48 instances, N ∈ {16, 32, 48, 64}, K ∈ {2, 3, 4}, δ2 ∈ {0.1, 0.2} Second testset: 81 instances, N ∈ {32, 48, 64}, K ∈ {4, 5, 6}, δ2 ∈ {0.1, 0.2, 0.3}

First testset: SCIP 4.0.1 SCIP 5.0.1 Setting #opt #nodes time (s) #opt #nodes time (s) default 48 3309 36.4 46 2787 43.1 default + SOC 48 2575 26.0 48 2478 29.1 default + heur 48 396 21.5 47 448 21.9 default + SOC + heur 48 323 21.1 48 364 18.8 mod handling 48 3283 28.3 47 3398 35.5 mod handling + SOC 48 2762 18.4 48 2313 19.7 mod handling + heur 48 427 19.3 48 507 18.0 mod handling + SOC + heur 48 317 17.9 48 361 15.2 Second testset: SCIP 4.0.1 SCIP 5.0.1 Setting #opt #nodes time (s) #opt #nodes time (s) default 47 126710 1058.8 35 38884 1232.2 default + SOC 47 149482 812.0 40 60565 984.2 default + heur 56 29122 444.4 45 15310 590.2 default + SOC + heur 54 33296 455.9 51 18000 472.0 mod handling 43 154338 1070.4 34 85619 1257.9 mod handling + SOC 55 167485 666.0 49 138182 746.2 mod handling + heur 55 35912 433.6 54 28203 461.2 mod handling + SOC + heur 59 40334 360.2 58 36938 306.9 March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 21

Overview

Motivation and Introduction Solving CMPs A heuristic Application Implementation Conclusion and Outline

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 22

Summary and further research

Conclusion:

◮ Specialized methods for handling modulus constraints result in speed-up. ◮ A good primal bound has a (very) high impact on the solution process, as

significantly less nodes are needed.

◮ Explicitly building the convex error bound term as SOC greatly fastens the

solution process.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 23

Summary and further research

Conclusion:

◮ Specialized methods for handling modulus constraints result in speed-up. ◮ A good primal bound has a (very) high impact on the solution process, as

significantly less nodes are needed.

◮ Explicitly building the convex error bound term as SOC greatly fastens the

solution process. Outlook:

◮ Let the heuristic also run in deeper nodes of the branch-and-cut tree and

initialize it with current LP solution.

◮ Closely connected: vary the two parameters for the number of guesses and

updates.

◮ Apply this framework to phase retrieval: For a measured modulus |x| of a

complex signal x = |x|eiψ find the phase ψ satisfying a set of constraints.

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 23

Thanks for your attention!

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 24

References

[WSA’18] T. Fischer, G. Hegde, F. Matter, M. Pesavento, M. E. Pfetsch, A. M. Tillmann, “Joint Antenna Selection and Phase-only Beamforming using Mixed-Integer Nonlinear Programming”, in 22nd International ITG Workshop on Smart Antennas (WSA), Bochum, Germany, March 2018. [accepted]

March 07, 2018 | SCIP Workshop, Aachen | Frederic Matter | 25