This notebook will provide implementation and examples from the paper Time-varying Convex Optimization, Amir Ali Ahmadi and Bachir El Khadir, 2018.
In [0]:
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
In [0]:
!pip install cvxpy
!pip install sympy
In [0]:
import numpy as np
import scipy as sp
The TV-SDP framework for CVXPY for imposing constraints of the form:
$$A(t) \succeq 0 \; \forall t \in [0, 1],$$where $$A(t)$$ is a polynomial symmetric matrix, i.e. a symmetric matrix whose entries are polynomial functions of time, and $$A(t) \succeq 0$$ means that all the eigen values of the matrix $$A(t)$$ are nonnegative.
In [0]:
def _mult_poly_matrix_poly(p, mat_y):
"""Multiplies the polynomial matrix mat_y by the polynomial p entry-wise.
Args:
p: list of size d1+1 representation the polynomial sum p[i] t^i.
mat_y: (m, m, d2+1) tensor representing a polynomial
matrix Y_ij(t) = sum mat_y[i, j, k] t^k.
Returns:
(m, m, d1+d2+1) tensor representing the polynomial matrix p(t)*Y(t).
"""
mult_op = lambda q: np.convolve(p, q)
p_times_y = np.apply_along_axis(mult_op, 2, mat_y)
return p_times_y
def _make_zero(p):
"""Returns the constraints p_i == 0.
Args:
p: list of cvxpy expressions.
Returns:
A list of cvxpy constraints [pi == 0 for pi in p].
"""
return [pi == 0 for pi in p]
def _lambda(m, d, Q):
"""Returns the mxm polynomial matrix of degree d whose Gram matrix is Q.
Args:
m: size of the polynomial matrix to be returned.
d: degreen of the polynomial matrix to be returned.
Q: (m*d/2, m*d/2) gram matrix of the polynomial matrix to be returned.
Returns:
(m, m, d+1) tensor representing the polynomial whose gram matrix is Q.
i.e. $$Y_ij(t) == sum_{r, s s.t. r+s == k} Q_{y_i t^r, y_j t^s} t^k$$.
"""
d_2 = int(d / 2)
def y_i_j(i, j):
poly = list(np.zeros((d + 1, 1)))
for k in range(d_2 + 1):
for l in range(d_2 + 1):
poly[k + l] += Q[i + k * m, j + l * m]
return poly
mat_y = [[y_i_j(i, j) for j in range(m)] for i in range(m)]
mat_y = np.array(mat_y)
return mat_y
def _alpha(m, d, Q):
"""Returns t*Lambda(Q) if d odd, Lambda(Q) o.w.
Args:
m: size of the polynomial matrix to be returned.
d: degreen of the polynomial matrix to be returned.
Q: gram matrix of the polynomial matrix.
Returns:
t*Lambda(Q) if d odd, Lambda(Q) o.w.
"""
if d % 2 == 1:
w1 = np.array([0, 1]) # t
else:
w1 = np.array([1]) # 1
mat_y = _lambda(m, d + 1 - len(w1), Q)
return _mult_poly_matrix_poly(w1, mat_y)
def _beta(m, d, Q):
"""Returns (1-t)*Lambda(Q) if d odd, t(1-t)*Lambda(Q) o.w.
Args:
m: size of the polynomial matrix to be returned.
d: degreen of the polynomial matrix to be returned.
Q: gram matrix of the polynomial matrix.
Returns:
(1-t)*Lambda(Q) if d odd, t(1-t)*Lambda(Q) o.w.
"""
if d % 2 == 1:
w2 = np.array([1, -1]) # 1 - t
else:
w2 = np.array([0, 1, -1]) # t - t^2
mat_y = _lambda(m, d + 1 - len(w2), Q)
return _mult_poly_matrix_poly(w2, mat_y)
def make_poly_matrix_psd_on_0_1(mat_x):
"""Returns the constraint X(t) psd on [0, 1].
Args:
mat_x: (m, m, d+1) tensor representing a mxm polynomial matrix of degree d.
Returns:
A list of cvxpy constraints imposing that X(t) psd on [0, 1].
"""
m, m2, d = len(mat_x), len(mat_x[0]), len(mat_x[0][0]) - 1
# square matrix
assert m == m2
# build constraints: X == alpha(Q1) + beta(Q2) with Q1, Q2 >> 0
d_2 = int(d / 2)
size_Q1 = m * (d_2 + 1)
size_Q2 = m * d_2 if d % 2 == 0 else m * (d_2 + 1)
Q1 = cvxpy.Variable((size_Q1, size_Q1))
Q2 = cvxpy.Variable((size_Q2, size_Q2))
diff = mat_x - _alpha(m, d, Q1) - _beta(m, d, Q2)
diff = diff.reshape(-1)
const = _make_zero(diff)
const += [Q1 >> 0, Q2 >> 0, Q1.T == Q1, Q2.T == Q2]
return const
In [0]:
def integ_poly_0_1(p):
"""Return the integral of p(t) between 0 and 1."""
return np.array(p).dot(1 / np.linspace(1, len(p), len(p)))
def spline_regression(x, y, num_parts, deg=3, alpha=.01, smoothness=1):
"""Fits splines with `num_parts` to data `(x, y)`.
Finds a piecewise polynomial function `p` of degree `deg` with `num_parts`
pieces that minimizes the fitting error sum |y_i - p(x_i)| + alpha |p|_1.
Args:
x: [N] ndarray of input data. Must be increasing.
y: [N] ndarray, same size as `x`.
num_parts: int, Number of pieces of the piecewise polynomial function `p`.
deg: int, degree of each polynomial piece of `p`.
alpha: float, Regularizer.
smoothness: int, the desired degree of smoothness of `p`, e.g.
`smoothness==0` corresponds to a continuous `p`.
Returns:
[num_parts, deg+1] ndarray representing the piecewise polynomial `p`.
Entry (i, j) contains j^th coefficient of the i^th piece of `p`.
"""
# coefficients of the polynomial of p.
p = cvxpy.Variable((num_parts, deg + 1), name='p')
# convert to numpy format because it is easier to work with.
numpy_p = np.array([[p[i, j] for j in range(deg+1)] \
for i in range(num_parts)])
regularizer = alpha * cvxpy.norm(p, 1)
num_points_per_part = int(len(x) / num_parts)
smoothness_constraints = []
# cuttoff values
t = []
fitting_value = 0
# split the data into equal `num_parts` pieces
for i in range(num_parts):
# the part of the data that the current piece fits
sub_x = x[num_points_per_part * i:num_points_per_part * (i + 1)]
sub_y = y[num_points_per_part * i:num_points_per_part * (i + 1)]
# compute p(sub_x)
# pow_x = np.array([sub_x**k for k in range(deg + 1)])
# sub_p = polyval(sub_xnumpy_p[i, :].dot(pow_x)
sub_p = eval_poly_from_coefficients(numpy_p[i], sub_x)
# fitting value of the current part of p,
# equal to sqrt(sum |p(x_i) - y_i|^2), where the sum
# is over data (x_i, y_i) in the current piece.
fitting_value += cvxpy.norm(cvxpy.vstack(sub_p - sub_y), 1)
# glue things together by ensuring smoothness of the p at x1
if i > 0:
x1 = x[num_points_per_part * i]
# computes the derivatives p'(x1) for the left and from the right of x1
# x_deriv is the 2D matrix k!/(k-j)! x1^(k-j) indexed by (j, k)
x1_deriv = np.array(
[[np.prod(range(k - j, k)) * x1**(k - j)
for k in range(deg + 1)]
for j in range(smoothness + 1)]).T
p_deriv_left = numpy_p[i - 1].dot(x1_deriv)
p_deriv_right = numpy_p[i].dot(x1_deriv)
smoothness_constraints += [
cvxpy.vstack(p_deriv_left - p_deriv_right) == 0
]
t.append(x1)
min_loss = cvxpy.Minimize(fitting_value + regularizer)
prob = cvxpy.Problem(min_loss, smoothness_constraints)
prob.solve(verbose=False)
return _piecewise_polynomial_as_function(p.value, t)
def _piecewise_polynomial_as_function(p, t):
"""Returns the piecewise polynomial `p` as a function.
Args:
p: [N, d+1] array of coefficients of p.
t: [N] array of cuttoffs.
Returns:
The function f s.t. f(x) = p_i(x) if t[i] < x < t[i+1].
"""
def evaluate_p_at(x):
"""Returns p(x)."""
pieces = [x < t[0]] + [(x >= ti) & (x < ti_plusone) \
for ti, ti_plusone in zip(t[:-1], t[1:])] +\
[x >= t[-1]]
# pylint: disable=unused-variable
func_list = [
lambda u, pi=pi: eval_poly_from_coefficients(pi, u) for pi in p
]
return np.piecewise(x, pieces, func_list)
return evaluate_p_at
def eval_poly_from_coefficients(coefficients, x):
"""Evaluates the polynomial whose coefficients are `coefficients` at `x`."""
return coefficients.dot([x**i for i in range(len(coefficients))])