MOSEK ApS

Optimization problem over nonnegative trigonometric polynomials.

We consider the trigonometric polynomial

$$H(\omega) = x_0 + 2 \sum_{k=1}^n [ \Re(x_k) \cos(\omega k) + \Im(x_k) \sin(\omega k) ],$$

where $H(\omega)$ is a real valued function paramtrized by the complex vector $x\in {\bf C}^{n+1}$, and where $\Im(x_0)=0$.

The example shows how to construct a non-negative polynomial $H(\omega)\geq 0, \: \forall \omega$ that satisfies,

$$1 - \delta \leq H(\omega) \leq 1 + \delta, \quad \forall \omega \in [0, \omega_p]$$

while minimizing $\sup_{\omega\in [\omega_s,\pi]} H(\omega)$ over the variables $x$.

In the signal processing literature such a trigonometric polynomial is known as (the squared amplitude respons of) a Chebyshev lowpass filter.

A squared amplitude response $H(\omega)$ is always symmetric around $0$, so $\Im(x_k)=0$, and we consider only

$$H(\omega) = x_0 + 2 \sum_{k=1}^n x_k \cos(\omega k) $$

over the real vector $x\in {\bf R}^{n+1}$. The concepts in the example are readily applied to the case with $x\in {\bf C}^{n+1}$, however.

References:

  1. "Squared Functional Systems and Optimization Problems", Y. Nesterov, 2004.

  2. "Convex Optimization of Non-negative Polynomials: Structured Algorithms and Applications", Ph.D thesis, Y. Hachez, 2003.

In [13]:
import mosek
from   mosek.fusion import *
from   math import cos, pi
import numpy as np

Nonnegativity everywhere

Nesterov proved in [1] that $H(\omega) \geq 0, \: \forall \omega$ if and only if $$x_i = \langle T_i^{n+1}, X \rangle, \quad X \in {\mathcal H}^{n+1}_+,$$ where ${\mathcal H}_+$ is the cone of Hermitian semidefinite matrices and $T_i$ is a Toeplitz matrix $$[T_i]_{kl} = \left \{ \begin{array}{ll}1, & k-l=i\\0 & \text{otherwise}.\end{array} \right .$$ For example, for $n=2$ we have $$ T_0 = \left[\begin{array}{ccc} 1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1 \end{array} \right], \quad T_1 = \left[\begin{array}{ccc} 0 & 0 & 0\\1 & 0 & 0\\0 & 1 & 0 \end{array} \right], \quad T_2 = \left[\begin{array}{ccc} 0 & 0 & 0\\0 & 0 & 0\\1 & 0 & 0 \end{array} \right]. $$ In our case we have $\Im(x_i)=0$, i.e., $X\in {\mathcal S}^{n+1}_+$ is a real symmetric semidefinite matrix.

We define the cone on nonnegative trigonometric polynomials as $$ K^n_{0,\pi} = \{ x\in \mathbf{R} \times \mathbf{C}^n \mid x_i = \langle X, T_i\rangle, \: i=0,\dots,n, \: X\in\mathcal{H}_+^{n+1} \}. $$

In [23]:
def T_dot_X(n, i, X, a=1.0):
    print n,i
    return Expr.dot(Matrix.diag(n, a, -i), X)
In [24]:
def trigpoly_0_pi(M, x):
    '''x[i] == <T(n+1,i),X>'''
    n = int(x.size()-1)
    X = M.variable("X", Domain.inPSDCone(n+1))
    
    for i in range(n+1):
        M.constraint(Expr.sub(T_dot_X(n+1,i,X), x.index(i)), Domain.equalsTo(0.0))

Note that we have dropped the imaginary part of $X$.

Nonnegativity on $[0,a]$

Similarly, $H(\omega)$ is nonnegative on $[0,a]$ if and only if

$$x_i = \langle X_1, T_i^{n+1} \rangle + \langle X_2, T_{i+1}^n \rangle + \langle X_2, T_{i-1}^n \rangle - 2 \cos(a)\langle X_2, T_{i}^n \rangle, \quad X_1 \in {\mathcal H}^{n+1}_+, \: X_2 \in {\mathcal H}^n_+, $$ or equivalently $$ K^n_{0,a} = \{ x\in \mathbf{R}\times \mathbf{C}^n \mid x_i = \langle X_1, T_i^{n+1} \rangle + \langle X_2 , T_{i+1}^n \rangle + \langle X_2 , T_{i-1}^n \rangle - 2\cos(\alpha)\langle X_2 , T_i^n\rangle, \: X_1\in \mathcal{H}_+^{n+1}, X_2\in\mathcal{H}_+^n \}. $$

In [25]:
def trigpoly_0_a(M, x, a):
    '''x[i] == <T(n+1,i),X1> + <T(n,i+1),X2> + <T(n,i-1),X2> - 2*cos(a)*<T(n,i),X2>'''
    n = int(x.size()-1)
    X1 = M.variable(Domain.inPSDCone(n+1))
    X2 = M.variable(Domain.inPSDCone(n))

    for i in range(n+1):
        M.constraint(Expr.sub(Expr.add([ T_dot_X(n+1,i,X1), 
                                         T_dot_X(n,i+1,X2), 
                                         T_dot_X(n,i-1,X2), 
                                         T_dot_X(n,i,X2,-2*cos(a))]),
                              x.index(i)), Domain.equalsTo(0.0))

Note that we have dropped the imaginary part of $X_1$ and $X_2$.

Nonnegativity on $[a,\pi]$

Finally, $H(\omega)$ is nonnegative on $[a,\pi]$ if and only if

$$x_i = \langle X_1, T_i^{n+1} \rangle - \langle X_2, T_{i+1}^n \rangle - \langle X_2, T_{i-1}^n \rangle + 2 \cos(a)\langle X_2, T_{i}^n \rangle, \quad X_1 \in {\mathcal S}^{n+1}_+, \: X_2 \in {\mathcal S}^n_+, $$ or equivalently $$ K^n_{a,\pi} = \{ x\in \mathbf{R}\times \mathbf{C}^n \mid x_i = \langle X_1, T_i^{n+1} \rangle - \langle X_2 , T_{i+1}^n \rangle - \langle X_2 , T_{i-1}^n \rangle + 2\cos(\alpha)\langle X_2 , T_i^n\rangle, \: X_1\in \mathcal{H}_+^{n+1}, X_2\in\mathcal{H}_+^n \}. $$

In [26]:
def trigpoly_a_pi(M, x, a):
    '''x[i] == <T(n+1,i),X1> - <T(n,i+1),X2> - <T(n,i-1),X2> + 2*cos(a)*<T(n,i),X2>'''
    n = int(x.size()-1)
    X1 = M.variable(Domain.inPSDCone(n+1))
    X2 = M.variable(Domain.inPSDCone(n))

    for i in range(n+1):
        M.constraint(Expr.add([ T_dot_X(n+1,i,X1), 
                                T_dot_X(n,i+1,X2,-1.0), 
                                T_dot_X(n,i-1,X2,-1.0), 
                                T_dot_X(n,i,X2,2*cos(a)),
                                Expr.neg(x.index(i))    ]),
                     Domain.equalsTo(0.0))

Note that we have dropped the imaginary part of $X_1$ and $X_2$.

An epigraph formulation

The epigraph $H(\omega) \leq t$ can now be characterized simply as $-u \in K^n_{[a,b]}, \: u=(x_0-t, \, x_{1:n}).$

In [27]:
def epigraph(M, x, t, a, b):
    '''Models 0 <= H(w) <= t, for all w in [a, b], where
         H(w) = x0 + 2*x1*cos(w) + 2*x2*cos(2*w) + ... + 2*xn*cos(n*w)'''
    n  = int(x.size()-1)    
    u = M.variable(n+1, Domain.unbounded())
    M.constraint(Expr.sub(t,Expr.add(x.index(0), u.index(0))), Domain.equalsTo(0.0)) 
    M.constraint(Expr.add(x.slice(1,n+1), u.slice(1,n+1)), Domain.equalsTo(0.0))
    
    if a==0.0 and b==pi:
        trigpoly_0_pi(M, u)
    elif a==0.0 and b<pi:
        trigpoly_0_a(M, u, b)
    elif a<pi and b==pi:
        trigpoly_a_pi(M, u, a)
    else:
        raise ValueError("invalid interval.")

A hypograph formulation

Similarly, the hypograph $H(\omega) \geq t$ can be characterized as $u \in K^n_{[a,b]}, \: u=(x_0-t, \, x_{1:n}).$

In [28]:
def hypograph(M, x, t, a, b):
    '''Models 0 <= t <= H(w), for all w in [a, b], where
         H(w) = x0 + 2*x1*cos(w) + 2*x2*cos(2*w) + ... + 2*xn*cos(n*w)'''

    n  = int(x.size()-1)    
    u0 = M.variable(1, Domain.unbounded())    
    M.constraint(Expr.sub(t,Expr.sub(x.index(0), u0)), Domain.equalsTo(0.0))
    u = Var.vstack( u0, x.slice(1, n+1) )

    if a==0.0 and b==pi:
        trigpoly_0_pi(M, u)
    elif a==0.0 and b<pi:
        trigpoly_0_a(M, u,  b)
    elif a<pi and b==pi:
        trigpoly_a_pi(M, u, a)
    else:
        raise ValueError("invalid interval.")

Putting it together

In [29]:
n = 10
M = Model("trigpoly")
x = M.variable("x", n+1, Domain.unbounded())

Global nonnegativity

In [30]:
# H(w) >= 0
trigpoly_0_pi(M, x)
11 0
11 1
11 2
11 3
11 4
11 5
11 6
11 7
11 8
11 9
11 10

Passband specifications

In [31]:
wp = pi/4
delta = 0.05

# H(w) <= (1+delta), w in [0, wp]    
epigraph(M, x, 1.0+delta, 0.0, wp)

# (1-delta) <= H(w), w in [0, wp]
hypograph(M, x, 1.0-delta, 0.0, wp)
11 0
10 1
10 -1
10 0
11 1
10 2
10 0
10 1
11 2
10 3
10 1
10 2
11 3
10 4
10 2
10 3
11 4
10 5
10 3
10 4
11 5
10 6
10 4
10 5
11 6
10 7
10 5
10 6
11 7
10 8
10 6
10 7
11 8
10 9
10 7
10 8
11 9
10 10
---------------------------------------------------------------------------
DimensionError                            Traceback (most recent call last)
<ipython-input-31-752dc0f9c3f1> in <module>()
      3 
      4 # H(w) <= (1+delta), w in [0, wp]
----> 5 epigraph(M, x, 1.0+delta, 0.0, wp)
      6 
      7 # (1-delta) <= H(w), w in [0, wp]

<ipython-input-27-fa020733500e> in epigraph(M, x, t, a, b)
     10         trigpoly_0_pi(M, u)
     11     elif a==0.0 and b<pi:
---> 12         trigpoly_0_a(M, u, b)
     13     elif a<pi and b==pi:
     14         trigpoly_a_pi(M, u, a)

<ipython-input-25-149899b26b29> in trigpoly_0_a(M, x, a)
      7     for i in range(n+1):
      8         M.constraint(Expr.sub(Expr.add([ T_dot_X(n+1,i,X1), 
----> 9                                          T_dot_X(n,i+1,X2),
     10                                          T_dot_X(n,i-1,X2),
     11                                          T_dot_X(n,i,X2,-2*cos(a))]),

<ipython-input-23-4f9ef9fb3b2c> in T_dot_X(n, i, X, a)
      1 def T_dot_X(n, i, X, a=1.0):
      2     print n,i
----> 3     return Expr.dot(Matrix.diag(n, a, -i), X)

/home/aszek/.local/lib/python2.7/site-packages/mosek/fusion/impl/_implementation.pyc in diag(*args)
  35423       return mosek_fusion_Matrix._diag_alt_ID(*args)
  35424     elif mosek_fusion_Matrix._match_diag_IDI(*args): # int32,double,int32
> 35425       return mosek_fusion_Matrix._diag_IDI(*args)
  35426     elif mosek_fusion_Matrix._match_alt_diag_IDI(*args): # int32,double,int32
  35427       return mosek_fusion_Matrix._diag_alt_IDI(*args)

/home/aszek/.local/lib/python2.7/site-packages/mosek/fusion/impl/_implementation.pyc in _diag_IDI(_0, _1, _2)
  35666     return (mosek.fusion.Matrix._diag__3DI(mosek.fusion.Utils.Tools._makevector_DI(_1,(_0 + _2)),_2))
  35667    else:
> 35668     raise mosek_fusion_DimensionError._ctor_S("Diagonal index out of bounds")
  35669   @staticmethod
  35670   def _match_diag_ID(*args):

DimensionError: Diagonal index out of bounds

Stopband specifications

In [12]:
ws = wp + pi/8

# H(w) < t, w in [ws, pi]
t = M.variable("t", 1, Domain.greaterThan(0.0))
epigraph(M, x, t, ws, pi)
---------------------------------------------------------------------------
DimensionError                            Traceback (most recent call last)
<ipython-input-12-321bdaca7120> in <module>()
      3 # H(w) < t, w in [ws, pi]
      4 t = M.variable("t", 1, Domain.greaterThan(0.0))
----> 5 epigraph(M, x, t, ws, pi)

<ipython-input-7-fa020733500e> in epigraph(M, x, t, a, b)
     12         trigpoly_0_a(M, u, b)
     13     elif a<pi and b==pi:
---> 14         trigpoly_a_pi(M, u, a)
     15     else:
     16         raise ValueError("invalid interval.")

<ipython-input-6-8bfb08a61e8e> in trigpoly_a_pi(M, x, a)
      7     for i in range(n+1):
      8         M.constraint(Expr.add([ T_dot_X(n+1,i,X1), 
----> 9                                 T_dot_X(n,i+1,X2,-1.0),
     10                                 T_dot_X(n,i-1,X2,-1.0),
     11                                 T_dot_X(n,i,X2,2*cos(a)),

<ipython-input-3-c4425f4dd099> in T_dot_X(n, i, X, a)
      1 def T_dot_X(n, i, X, a=1.0):
----> 2     return Expr.dot(Matrix.diag(n, a, -i), X)

/home/aszek/.local/lib/python2.7/site-packages/mosek/fusion/impl/_implementation.pyc in diag(*args)
  35423       return mosek_fusion_Matrix._diag_alt_ID(*args)
  35424     elif mosek_fusion_Matrix._match_diag_IDI(*args): # int32,double,int32
> 35425       return mosek_fusion_Matrix._diag_IDI(*args)
  35426     elif mosek_fusion_Matrix._match_alt_diag_IDI(*args): # int32,double,int32
  35427       return mosek_fusion_Matrix._diag_alt_IDI(*args)

/home/aszek/.local/lib/python2.7/site-packages/mosek/fusion/impl/_implementation.pyc in _diag_IDI(_0, _1, _2)
  35666     return (mosek.fusion.Matrix._diag__3DI(mosek.fusion.Utils.Tools._makevector_DI(_1,(_0 + _2)),_2))
  35667    else:
> 35668     raise mosek_fusion_DimensionError._ctor_S("Diagonal index out of bounds")
  35669   @staticmethod
  35670   def _match_diag_ID(*args):

DimensionError: Diagonal index out of bounds

Setting the objective

In [12]:
M.objective(ObjectiveSense.Minimize, t)
In [13]:
#M.setLogHandler(sys.stdout)
In [14]:
M.solve()
In [15]:
x.level()
Out[15]:
array([ 0.33580937,  0.25254317,  0.1388759 ,  0.02044169, -0.04797007,
       -0.05156956, -0.01627207,  0.01837092,  0.02948826,  0.02271378,
       -0.00952669])
In [16]:
t.level()
Out[16]:
array([ 0.0727316])

Plotting the amplitude response

In [17]:
def H(x, w): return x[0] + 2*sum([x[i]*cos(i*w) for i in range(1,len(x))])
w  = np.linspace(0, pi, 100)
T = [ H(x.level(), wi) for wi in w ]
In [18]:
import bokeh
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import FixedTicker, BoxAnnotation
In [19]:
output_notebook()
Loading BokehJS ...
In [20]:
p = figure(x_axis_label='w', x_range=(0,pi), y_axis_label='H(w)', y_range=(0,1 + 2*delta))
p.line(w, T, line_width=2)
p.yaxis[0].ticker=FixedTicker(ticks=[0, t.level()[0], 1-delta, 1+delta])
p.add_layout(BoxAnnotation(top=t.level()[0], left=ws, fill_alpha=0.1, fill_color='red'))
p.add_layout(BoxAnnotation(top=1+delta, bottom=1-delta, right=wp, fill_alpha=0.1, fill_color='red'))
show(p)
Out[20]:

<Bokeh Notebook handle for In[20]>

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License. The MOSEK logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of MOSEK or the Fusion API are not guaranteed. For more information contact our support.