# coding: utf-8
# ------------------------------------------------------------------------------
# Copyright (C) 2019 Maximilian Stahlberg
# Based on the original picos.expressions module by Guillaume Sagnol.
#
# This file is part of PICOS.
#
# PICOS is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PICOS is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------------
"""Implements :class:`QuadraticExpression`."""
import operator
import sys
from collections import OrderedDict as odict
from collections import namedtuple
import cvxopt
import cvxopt.cholmod
import numpy
from .. import glyphs
from ..apidoc import api_end, api_start
from ..caching import (borrow_cache, cached_property,
cached_selfinverse_unary_operator,
cached_unary_operator)
from ..constraints import (ConicQuadraticConstraint, Constraint,
ConvexQuadraticConstraint,
NonconvexQuadraticConstraint)
from .data import convert_operands, cvx2np, should_be_sparse
from .exp_affine import AffineExpression, ComplexAffineExpression, Constant
from .expression import Expression, refine_operands, validate_prediction
from .variables import BaseVariable
_API_START = api_start(globals())
# -------------------------------
PSD_PERTURBATIONS = 3
r"""Maximum number of singular quadratic form perturbations.
PICOS uses NumPy either or CHOLMOD to compute a Cholesky decomposition of a
positive semidefinite quadratic form :math:`Q`. Both libraries require
:math:`Q` to be nonsingular, which is not a requirement on PICOS' end.
If either library rejects :math:`Q`, then PICOS provides a sequence of
:data:`PSD_PERTURBATIONS` perturbed quadratic forms :math:`Q + \epsilon I` for
increasing :math:`\epsilon` until the perturbed matrix is found positive
definite or until the largest :math:`\epsilon` was tested unsuccessfully.
If this is zero, PICOS will only decompose quadratic forms that are nonsingular.
If this is one, then only the largest epsilon is tested.
"""
[docs]class QuadraticExpression(Expression):
"""A scalar quadratic expression of the form :math:`x^TQx + a^Tx + b`."""
# --------------------------------------------------------------------------
# Initialization and factory methods.
# --------------------------------------------------------------------------
[docs] def __init__(
self, string, quadraticPart={}, affinePart=AffineExpression.zero(),
scalarFactors=None, copyDecomposition=None):
r"""Initialize a scalar quadratic expression.
This constructor is meant for internal use. As a user, you will most
likely want to build expressions starting from
:mod:`~picos.expressions.variables` or a :func:`~picos.Constant`.
:param str string: A symbolic string description.
:param quadraticPart: A :class:`dict` mapping PICOS variable pairs to
CVXOPT matrices. Each entry :math:`(x, y) \mapsto A_{xy}` represents
a quadratic form
:math:`\operatorname{vec_x}(x)^T A_{xy} \operatorname{vec_y}(y)`
where :math:`\operatorname{vec_x}` and :math:`\operatorname{vec_y}`
refer to the isometric real
:mod:`~.picos.expressions.vectorizations` that are used by PICOS
internally to store the variables :math:`x` and :math:`y`,
respectively. The quadratic part :math:`Q` of the expression is then
given as :math:`Q = \sum_{x, y} A_{xy}`.
:param affinePart: The affine part :math:`a^Tx + b` of the expression.
:type affinePart: ~picos.expressions.AffineExpression
:param scalarFactors: A pair :math:`(u, v)` with both :math:`u` and
:math:`v` scalar real affine expressions representing a known
:math:`x^TQx + a^Tx + b = uv` factorization of the expression.
:type factorization:
tuple(~picos.expressions.AffineExpression)
:param copyDecomposition: Another quadratic expression with equal
quadratic part whose quadratic part decomposition shall be copied.
:type copyDecomposition:
~picos.expressions.QuadraticExpression
"""
# Ensure correct usage within PICOS.
assert isinstance(quadraticPart, dict) and all(
isinstance(xy, (tuple, list)) and len(xy) == 2
and isinstance(xy[0], BaseVariable)
and isinstance(xy[1], BaseVariable)
and isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
and A.size == (xy[0].dim, xy[1].dim)
for xy, A in quadraticPart.items())
assert isinstance(affinePart, ComplexAffineExpression) \
and len(affinePart) == 1
assert not copyDecomposition \
or isinstance(copyDecomposition, QuadraticExpression)
assert not scalarFactors or (
isinstance(scalarFactors, (tuple, list)) and len(scalarFactors) == 2
and all(isinstance(x, AffineExpression) for x in scalarFactors)
and all(len(x) == 1 for x in scalarFactors))
Expression.__init__(self, "Quadratic Expression", string)
# Check affine part.
affinePart = affinePart.refined
if not isinstance(affinePart, AffineExpression):
raise NotImplementedError(
"PICOS does not support complex-valued quadratic expressions "
"and the affine part of {} is not real.".format(self._symbStr))
elif len(affinePart) != 1:
raise TypeError("Affine part of {} is not scalar."
.format(self._symbStr))
# Store quadratic part.
self._quads = {}
for xy, A in quadraticPart.items():
# Ignore coefficients that are already zero.
if not A:
continue
if xy[0] is xy[1]:
# Store the unique symmetric form for two reasons:
# - Hermitian forms concerning the original variables are
# supplied as quadratic forms concerning the isometric real
# vectorizations of the variables. These matrices are in
# general still complex, but their symmetric form is not.
# - Coefficients that cancel out are set to numeric zeros and
# can be converted to structural zeros by _sparse_quads.
A = 0.5*(A + A.T)
if hash(xy[1]) < hash(xy[0]):
# Store only one coefficient per unordered pair.
xy, A = (xy[1], xy[0]), A.T
if xy in self._quads:
self._quads[xy] = self._quads[xy] + A
else:
self._quads[xy] = A
# Remove coefficients of zero.
self._quads = {xy: A for xy, A in self._quads.items() if A}
# Check if coefficients are real.
for (x, y), A in self._quads.items():
if A.typecode == "z":
if any(A.imag()):
raise NotImplementedError(
"PICOS does not support complex-valued quadratic "
"expressions and the quadratic part of {} for variables"
" {} and {} is not real or equivalent to a real form."
.format(self._symbStr, x.string, y.string))
self._quads[xy] = A.real()
# Store affine part.
self._aff = affinePart
# Store a known factorization into real scalars.
if scalarFactors:
a, b = scalarFactors
if a.equals(b):
self._sf = (a, a) # Speedup future equality checks.
else:
self._sf = (a, b)
else:
self._sf = None
# Copy a decomposition from another quadratic expression.
if copyDecomposition is not None:
borrow_cache(self, copyDecomposition, ("_Q_and_x", "L", "quadroot"))
# --------------------------------------------------------------------------
# Abstract method implementations and method overridings, except _predict.
# --------------------------------------------------------------------------
@cached_unary_operator
def _get_refined(self):
if not self._quads:
return self._aff.refined.renamed(self._symbStr)
else:
return self
Subtype = namedtuple("Subtype", (
"convex", "concave", "qterms", "rank", "haslin", "factors"))
@cached_unary_operator
def _get_subtype(self):
convex = self.convex
concave = False if convex and self._quads else self.concave
qterms = self.num_quad_terms
haslin = not self._aff.constant
try:
if convex:
rank = self.rank
elif concave:
rank = (-self).rank
else:
rank = None
except ValueError: # No quadratic part.
rank = 0
factors = len(set(self.scalar_factors))
return self.Subtype(convex, concave, qterms, rank, haslin, factors)
def _get_value(self):
value = self._aff._get_value()
for xy, A in self._quads.items():
x, y = xy
xVal = x._get_internal_value()
yVal = y._get_internal_value()
value = value + xVal.T * A * yVal
return value
@cached_unary_operator
def _get_variables(self):
return self._aff.variables.union(
var for vars in self._quads for var in vars)
def _is_convex(self):
try:
self.L
except ValueError:
return True # Expression is affine.
except ArithmeticError:
return False
else:
return True
def _is_concave(self):
try:
(-self).L
except ValueError:
return True # Expression is affine.
except ArithmeticError:
return False
else:
return True
def _replace_variables(self, var_map):
name_map = {old.name: new.name for old, new in var_map.items()}
string = self.string
if isinstance(string, glyphs.GlStr):
string = string.reglyphed(name_map)
quads = {(var_map[var1], var_map[var2]): coef
for (var1, var2), coef in self._quads.items()}
if not self._sf:
sf = None
elif self._sf[0] is self._sf[1]:
sf = (self._sf[0]._replace_variables(var_map),)*2
else:
sf = tuple(f._replace_variables(var_map) for f in self._sf)
return QuadraticExpression(
string, quads, self.aff._replace_variables(var_map), sf)
# --------------------------------------------------------------------------
# Python special method implementations, except constraint-creating ones.
# --------------------------------------------------------------------------
def __len__(self):
# Faster version that overrides Expression.__len__.
return 1
@convert_operands(sameShape=True)
def __add__(self, other):
"""Denote addition from the right hand side."""
if isinstance(other, QuadraticExpression):
string = glyphs.clever_add(self.string, other.string)
quadPart = {}
varPairs = set(self._quads.keys()).union(other._quads.keys())
for vars in varPairs:
if vars in self._quads:
if vars in other._quads:
quadPart[vars] = self._quads[vars] + other._quads[vars]
else:
quadPart[vars] = self._quads[vars]
else:
quadPart[vars] = other._quads[vars]
affPart = self._aff + other._aff
return QuadraticExpression(string, quadPart, affPart)
elif isinstance(other, AffineExpression):
string = glyphs.clever_add(self.string, other.string)
return QuadraticExpression(
string, self._quads, self._aff + other, copyDecomposition=self)
else:
return NotImplemented
@convert_operands(sameShape=True)
def __radd__(self, other):
"""Denote addition from the left hand side."""
if isinstance(other, QuadraticExpression):
assert False, "Impossible case."
elif isinstance(other, AffineExpression):
result = self.__add__(other)
# NOTE: __add__ always creates a fresh expression.
result._symbStr = glyphs.clever_add(other.string, self.string)
return result
else:
return NotImplemented
@convert_operands(sameShape=True)
def __sub__(self, other):
"""Denote substraction from the right hand side."""
if isinstance(other, QuadraticExpression):
string = glyphs.clever_sub(self.string, other.string)
quadPart = {}
varPairs = set(self._quads.keys()).union(other._quads.keys())
for vars in varPairs:
if vars in self._quads:
if vars in other._quads:
quadPart[vars] = self._quads[vars] - other._quads[vars]
else:
quadPart[vars] = self._quads[vars]
else:
quadPart[vars] = -other._quads[vars]
affPart = self._aff - other._aff
return QuadraticExpression(string, quadPart, affPart)
elif isinstance(other, AffineExpression):
string = glyphs.clever_sub(self.string, other.string)
return QuadraticExpression(
string, self._quads, self._aff - other, copyDecomposition=self)
else:
return NotImplemented
@convert_operands(sameShape=True)
def __rsub__(self, other):
"""Denote substraction with self on the right hand side."""
if isinstance(other, QuadraticExpression):
assert False, "Impossible case."
elif isinstance(other, AffineExpression):
result = (-self).__add__(other)
# NOTE: __add__ always creates a fresh expression.
result._symbStr = glyphs.clever_sub(other.string, self.string)
return result
else:
return NotImplemented
@cached_selfinverse_unary_operator
def __neg__(self):
string = glyphs.clever_neg(self.string)
quadPart = {xy: -A for xy, A in self._quads.items()}
affPart = -self._aff
factors = (self._sf[0], -self._sf[1]) if self._sf else None
return QuadraticExpression(string, quadPart, affPart, factors)
@convert_operands(scalarRHS=True)
def __mul__(self, other):
"""Denote scaling from the right hand side."""
if isinstance(other, AffineExpression):
if not other.constant:
raise NotImplementedError("You may only multiply a PICOS "
"quadratic expression with a constant term.")
factor = other.value
string = glyphs.clever_mul(self.string, other.string)
quadPart = {xy: factor*A for xy, A in self._quads.items()}
affPart = self._aff*other
if self._sf:
r = factor**0.5
if self._sf[0] is self._sf[1]:
factors = (r*self._sf[0],)*2
else:
factors = (r*self._sf[0], r*self._sf[1])
else:
factors = None
return QuadraticExpression(string, quadPart, affPart, factors)
else:
return NotImplemented
@convert_operands(scalarRHS=True)
def __rmul__(self, other):
"""Denote scaling from the left hand side."""
if isinstance(other, AffineExpression):
result = self.__mul__(other)
# NOTE: __mul__ always creates a fresh expression.
result._symbStr = glyphs.clever_mul(other.string, self.string)
return result
else:
return NotImplemented
@convert_operands(scalarRHS=True)
def __truediv__(self, other):
"""Denote division by a constant scalar."""
if isinstance(other, AffineExpression):
if not other.constant:
raise NotImplementedError("You may only divide a PICOS "
"quadratic expression by a constant term.")
result = self.__mul__(1.0 / other.value)
# NOTE: __mul__ always creates a fresh expression.
result._symbStr = glyphs.div(self.string, other.string)
return result
else:
return NotImplemented
def __div__(self, other):
"""Denote division by a scalar constant."""
# Python 2 uses this instead of __truediv__.
return self.__truediv__(other)
# --------------------------------------------------------------------------
# Decomposition related methods.
# --------------------------------------------------------------------------
@cached_property
def _sparse_quads(self):
"""Quadratic coefficients (re-)cast to sparse matrices."""
return {xy: cvxopt.sparse(M) for xy, M in self._quads.items()}
def _Q_and_x_or_A_and_y(self, augmentedMatrix):
"""Either :attr:`_Q_and_x` or :attr:`_A_and_y`."""
# Find occurences of the scalar entries of the variable's isometric real
# vectorization in the quadratic part.
nonzero = {}
for (x, y), A in self._sparse_quads.items():
nonzero.setdefault(x, set())
nonzero.setdefault(y, set())
nonzero[x].update(A.I)
nonzero[y].update(A.J)
# For the augmented matrix, find additional linear occurences.
if augmentedMatrix:
for x, a in self._aff._sparse_coefs.items():
nonzero.setdefault(x, set())
nonzero[x].update(a.I)
nonzero = odict((var, sorted(nonzero[var]))
for var in sorted(nonzero, key=lambda v: v.id))
# Compute each (multidimensional) variable's offset in Q/A.
offset, offsets = 0, {}
for var, nz in nonzero.items():
offsets[var] = offset
offset += len(nz)
# Compute Q.
V, I, J = [], [], []
for (x, y), A in self._sparse_quads.items():
B = A[nonzero[x], nonzero[y]]
V.extend(B.V)
I.extend(offsets[x] + i for i in B.I)
J.extend(offsets[y] + j for j in B.J)
# Turn Q into A if requested.
if augmentedMatrix:
for x, a in self._aff._sparse_coefs.items():
b = a[nonzero[x]]
V.extend(b.V)
I.extend(offsets[x] + i for i in b.I)
J.extend([offset]*len(b.J))
V.append(self._aff._const[0])
I.append(offset)
J.append(offset)
offset += 1
# Finalize Q/A.
Q = cvxopt.spmatrix(V, I, J, (offset,)*2, tc="d")
Q = 0.5*(Q + Q.T)
if not Q:
if augmentedMatrix:
raise ValueError(
"The expression {} is zero.".format(self._symbStr))
else:
raise ValueError(
"The quadratic part of {} is zero.".format(self._symbStr))
# Compute x.
x = None
for var, nz in nonzero.items():
n, r = len(nz), var.dim
F = cvxopt.spmatrix([1]*n, range(n), nz, (n, r))
y = AffineExpression(
glyphs.Fn("F")(var.string), n, coefficients={var: F})
x = y if x is None else x // y
# Turn x into y if requested.
if augmentedMatrix:
x = AffineExpression.from_constant(1) if x is None else x // 1
return Q, x
@cached_property
def _Q_and_x(self):
"""Pair containing :attr:`Q` and :attr:`x`."""
return self._Q_and_x_or_A_and_y(augmentedMatrix=False)
@cached_property
def _A_and_y(self):
"""Pair containing :attr:`R` and :attr:`y`."""
return self._Q_and_x_or_A_and_y(augmentedMatrix=True)
@property
def Q(self):
"""The coefficient matrix :math:`Q` of the expression, condensed.
This equals the :math:`Q` of the quadratic expression
:math:`x^TQx + a^Tx + b` but with zero rows and columns removed.
The vector :attr:`x` is a condensed version of :math:`x` that refers to
this matrix, so that ``q.x.T*q.Q*q.x`` equals the quadratic part
:math:`x^TQx` of the expression ``q``.
:raises ValueError: When the quadratic part is zero.
"""
return self._Q_and_x[0]
@property
def x(self):
"""The stacked variable vector :math:`x` of the expression, condensed.
This equals the :math:`x` of the quadratic expression
:math:`x^TQx + a^Tx + b` but entries corresponding to zero rows and
columns in :attr:`Q` are removed, so that ``q.x.T*q.Q*q.x`` equals the
:math:`x^TQx` part of the expression ``q``.
:raises ValueError: When the quadratic part is zero.
"""
return self._Q_and_x[1]
@property
def A(self):
r"""An affine-augmented quadratic coefficient matrix, condensed.
For a quadratic expression :math:`x^TQx + a^Tx + b`, this is
:math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}`
but with zero rows and columns removed.
The vector :attr:`y` is a condensed version of :math:`x` that refers to
this matrix, so that ``q.y.T*q.A*q.y`` equals the expression ``q``.
:raises ValueError: When the expression is zero.
"""
return self._A_and_y[0]
@property
def y(self):
"""See :attr:`A`.
:raises ValueError: When the expression is zero.
"""
return self._A_and_y[1]
def _L_or_M(self, augmentedMatrix):
"""Either :attr:`L` or :attr:`M`."""
Q = self.A if augmentedMatrix else self.Q
n = Q.size[0]
# Define a sequence of small numbers for perturbation.
# TODO: These numbers are empirical; find worst case values.
spread = [abs(v) for v in Q.V if v]
minEps = 1 * min(spread) * sys.float_info.epsilon
maxEps = n * max(spread) * sys.float_info.epsilon
epsilons = set([0.0])
for i in range(PSD_PERTURBATIONS):
epsilons.add(minEps + i * (maxEps - minEps) / PSD_PERTURBATIONS)
epsilons = sorted(epsilons)
# Check if Q is really sparse.
sparse = should_be_sparse(Q.size, len(Q))
# Attempt to find L.
L = None
if sparse: # Use CHOLMOD.
F = None
for eps in epsilons:
P = Q + cvxopt.spdiag([eps]*n) if eps else Q # Perturbed Q.
if F is None:
F = cvxopt.cholmod.symbolic(P)
try:
cvxopt.cholmod.numeric(P, F)
except ArithmeticError:
pass
else:
L = cvxopt.cholmod.getfactor(F)
break
else: # Use NumPy.
Q = cvx2np(Q)
for eps in epsilons:
P = Q + numpy.diag([eps]*n) if eps else Q # Perturbed Q.
try:
L = numpy.linalg.cholesky(P)
except numpy.linalg.LinAlgError:
pass
else:
break
if L is None:
if not PSD_PERTURBATIONS:
raise ArithmeticError("The expression {} has an (augmented) "
"quadratic form that is either not positive semidefinite "
"or singular. Try enabling PSD_PERTURBATIONS.")
elif augmentedMatrix:
raise ArithmeticError("The expression {} is not a squared norm "
"as its {} representation is not positive semidefinite."
.format(self._symbStr, glyphs.matrix(glyphs.vertcat(
glyphs.horicat("Q", glyphs.div("a", 2)),
glyphs.horicat(glyphs.div(glyphs.transp("a"), 2), "b")
))))
else:
raise ArithmeticError("The expression {} is not convex as its "
"quadratic part is not positive semidefinite."
.format(self._symbStr))
# Remove near-zeros in the order of n·sqrt(maxEps).
cutoff = 2 * n * maxEps**0.5
if sparse:
VIJ = list(t for t in zip(L.V, L.I, L.J) if abs(t[0]) > cutoff)
if VIJ: # Don't zero the matrix.
V, I, J = zip(*VIJ)
L = cvxopt.spmatrix(V, I, J, L.size)
else:
mask = abs(L) <= cutoff
if not numpy.all(mask): # Don't zero the matrix.
L[mask] = 0
L = cvxopt.sparse(cvxopt.matrix(L))
return L
@cached_property
def L(self):
"""The :math:`L` of an :math:`LL^T` Cholesky decomposition of :attr:`Q`.
:returns: A CVXOPT lower triangular sparse matrix.
:raises ValueError: When the quadratic part is zero.
:raises ArithmeticError: When the expression is not convex, that is when
the matrix :attr:`Q` is not numerically positive semidefinite.
"""
return self._L_or_M(augmentedMatrix=False)
@cached_property
def M(self):
"""The :math:`M` of an :math:`MM^T` Cholesky decomposition of :attr:`A`.
:returns: A CVXOPT lower triangular sparse matrix.
:raises ValueError: When the expression is zero.
:raises ArithmeticError: When the expression can not be written as a
squared norm, that is when the extended matrix :attr:`A` is not
numerically positive semidefinite.
"""
return self._L_or_M(augmentedMatrix=True)
@cached_property
def quadroot(self):
r"""Affine expression whose squared norm equals :math:`x^TQx`.
For a convex quadratic expression ``q``, this is equal to the vector
``q.L.T*q.x`` with zero rows removed.
:Construction:
Let :math:`x^TQx` be the quadratic part of the expression with :math:`Q`
positive semidefinite and :math:`Q = LL^T` a Cholesky decomposition.
Then,
.. math::
x^TQx
&= x^TLL^Tx \\
&= (L^Tx)^TL^Tx \\
&= \langle L^Tx, L^Tx \rangle \\
&= \lVert L^Tx \rVert^2.
Note that removing zero rows from :math:`L^Tx` does not affect the norm.
:raises ValueError: When the quadratic part is zero.
:raises ArithmeticError: When the expression is not convex, that is when
the quadratic part is not numerically positive semidefinite.
"""
L = self.L
n = L.size[0]
nz = sorted(set(L.J))
nnz = len(nz)
# F removes rows corresponding to zero columns in L.
F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n))
result = (F*L.T)*self.x
# NOTE: AffineExpression.__mul__ always creates a fresh expression.
result._symbStr = glyphs.Fn("quadroot({})")(self.string)
return result
@cached_property
def fullroot(self):
r"""Affine expression whose squared norm equals the expression.
For a convex quadratic expression ``q``, this is equal to the vector
``q.M.T*q.y`` with zero rows removed.
:Construction:
For a quadratic expression :math:`x^TQx + a^Tx + b` with
:math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}`
positive semidefinite, let :math:`A = MM^T` be a Cholesky decomposition.
Let further :math:`y = \begin{bmatrix}x^T & 1\end{bmatrix}^T`. Then,
.. math::
x^TQx + a^Tx + b
&= y^TAy \\
&= y^TMM^Ty \\
&= (M^Ty)^TM^Ty \\
&= \langle M^Ty, M^Ty \rangle \\
&= \lVert M^Ty \rVert^2.
Note that removing zero rows from :math:`M^Ty` does not affect the norm.
:raises ValueError: When the expression is zero.
:raises ArithmeticError: When the expression is not convex, that is when
the quadratic part is not numerically positive semidefinite.
"""
M = self.M
n = M.size[0]
nz = sorted(set(M.J))
nnz = len(nz)
# F removes rows corresponding to zero columns in M.
F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n))
result = (F*M.T)*self.y
# NOTE: AffineExpression.__mul__ always creates a fresh expression.
result._symbStr = glyphs.Fn("fullroot({})")(self.string)
return result
@property
def rank(self):
"""The length of the vector :attr:`quadroot`.
Up to numerical considerations, this is the rank of the (convex)
quadratic coefficient matrix :math:`Q` of the expression.
:raises ArithmeticError: When the expression is not convex, that is when
the quadratic part is not numerically positive semidefinite.
"""
try:
return len(self.quadroot)
except ValueError:
return 0 # No quadratic part.
@cached_property
def num_quad_terms(self):
"""The number of terms in the simplified quadratic form."""
num_terms = 0
for A in self._sparse_quads.values():
for i, j in zip(A.I, A.J):
if i >= j:
num_terms += 1
return num_terms
@property
def is_squared_norm(self):
r"""Whether the expression can be written as a squared norm.
If this is :obj:`True`, then the there is a coefficient vector
:math:`c` such that
.. math::
\lVert c^T \begin{bmatrix} x \\ 1 \end{bmatrix} \rVert^2
= x^TQx + a^Tx + b.
If the expression is also nonzero, then :attr:`fullroot` is
:math:`c^T \begin{bmatrix} x^T & 1 \end{bmatrix}^T` with zero entries
removed.
"""
try:
self.M
except ValueError:
assert self.is0
return True
except ArithmeticError:
return False
else:
return True
# --------------------------------------------------------------------------
# Properties and functions that describe the expression.
# --------------------------------------------------------------------------
@property
def quadratic_forms(self):
"""The quadratic forms as sparse matrices."""
return dict(self._sparse_quads)
@property
def is0(self):
"""Whether the quadratic expression is zero."""
return not self._quads and self._aff.is0
@property
def scalar_factors(self):
r"""Decomposition into scalar real affine expressions.
If the expression is known to be equal to :math:`ab` for scalar real
affine expressions :math:`a` and :math:`b`, returns the pair ``(a, b)``.
Otherwise, returns the empty tuple ``()``.
Note that if :math:`a = b`, then also ``a is b`` in the returned tuple.
"""
if self._sf:
return tuple(self._sf)
else:
return ()
# --------------------------------------------------------------------------
# Methods and properties that return modified copies.
# --------------------------------------------------------------------------
@cached_property
def quad(self):
"""Quadratic part :math:`x^TQx` of the quadratic expression."""
factors = self._sf if self._sf and self._aff.is0 else None
return QuadraticExpression(glyphs.quadpart(self._symbStr), self._quads,
scalarFactors=factors, copyDecomposition=self)
@cached_property
def aff(self):
"""Affine part :math:`a^Tx + b` of the quadratic expression."""
return self._aff.renamed(glyphs.affpart(self._symbStr))
@cached_property
def lin(self):
"""Linear part :math:`a^Tx` of the quadratic expression."""
return AffineExpression(glyphs.linpart(self._symbStr),
coefficients=self._aff._coefs)
@cached_property
def cst(self):
"""Constant part :math:`b` of the quadratic expression."""
return AffineExpression(glyphs.cstpart(self._symbStr),
constant=self._aff._const)
# --------------------------------------------------------------------------
# Constraint-creating operators, and _predict.
# --------------------------------------------------------------------------
@classmethod
def _predict(cls, subtype, relation, other):
assert isinstance(subtype, cls.Subtype)
if issubclass(other.clstype, QuadraticExpression):
raise NotImplementedError("Constraint outcome prediction not "
"supported when comparing two quadratic expressions.")
if relation == operator.__le__:
if issubclass(other.clstype, AffineExpression) \
and other.subtype.dim == 1:
if subtype.convex:
return ConvexQuadraticConstraint.make_type(
qterms=subtype.qterms, rank=subtype.rank,
haslin=(subtype.haslin or not other.subtype.constant))
else:
return NonconvexQuadraticConstraint.make_type(
qterms=subtype.qterms)
elif relation == operator.__ge__:
if issubclass(other.clstype, AffineExpression) \
and other.subtype.dim == 1:
if subtype.concave:
return ConvexQuadraticConstraint.make_type(
qterms=subtype.qterms, rank=subtype.rank,
haslin=(subtype.haslin or not other.subtype.constant))
elif subtype.factors and other.subtype.constant \
and other.subtype.nonneg:
# See __ge__ for the variants below.
# Variant 1:
return ConicQuadraticConstraint.make_type(
qterms=subtype.qterms, conic_argdim=1,
rotated=(subtype.factors == 2))
# Variant 2:
# return RSOCConstraint.make_type(argdim=1)
# Variant 3:
# if subtype.factors == 1:
# return AffineConstraint.make_type(dim=1, eq=False)
# else:
# return RSOCConstraint.make_type(argdim=1)
else:
return NonconvexQuadraticConstraint.make_type(
qterms=subtype.qterms)
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __le__(self, other):
LE = Constraint.LE
if isinstance(other, QuadraticExpression):
le0 = self - other # Unpredictable outcome.
if le0.convex:
return ConvexQuadraticConstraint(self, LE, other)
elif self.is_squared_norm and other.scalar_factors:
return ConicQuadraticConstraint(self, LE, other)
else:
return NonconvexQuadraticConstraint(self, LE, other)
elif isinstance(other, AffineExpression):
if self.convex:
return ConvexQuadraticConstraint(self, LE, other)
else:
return NonconvexQuadraticConstraint(self, LE, other)
else:
return NotImplemented
@convert_operands(sameShape=True)
@validate_prediction
@refine_operands()
def __ge__(self, other):
GE = Constraint.GE
if isinstance(other, QuadraticExpression):
le0 = other - self # Unpredictable outcome.
if le0.convex:
return ConvexQuadraticConstraint(self, GE, other)
elif other.is_squared_norm and self.scalar_factors:
return ConicQuadraticConstraint(self, GE, other)
else:
return NonconvexQuadraticConstraint(self, GE, other)
elif isinstance(other, AffineExpression):
if self.concave:
return ConvexQuadraticConstraint(self, GE, other)
elif self._sf and other.constant and other.value >= 0:
# Handle the case of c <= x*x and c <= x*y with c pos. constant.
root = Constant(glyphs.sqrt(other.string), other.value**0.5)
# Variant 1: Return a ConicQuadraticConstraint that can be
# reformulated to an (R)SOCConstraint. This is consistent
# with all other quadratic expression comparison outcomes.
other = QuadraticExpression(other.string, affinePart=other,
scalarFactors=(root, root)) # Unrefined QE.
return ConicQuadraticConstraint(self, GE, other)
# Variant 2: Return an RSOCConstraint. This is somewhat
# consistent with Norm returning an SOCConstraint and it is
# also legacy behavior. However this may not be what the
# user wants, because it silently forces x, y >= 0!
# return RSOCConstraint(root, self._sf[0], self._sf[1],
# customString = glyphs.le(other.string, self.string))
# Variant 3: Return either an RSOCConstraint or an affine
# constraint that equals an SOCConstraint. This is faster
# than Variant 2 but even more confusing, because the
# constraint would transform from e.g. 1 <= x*x to
# sqrt(1) <= x without any mention of cones.
# if self._sf[0] is self._sf[1]:
# # NOTE: The following is not allowed by SOCConstraint
# # since this is effectively an affine constraint.
# # return SOCConstraint(root, self._sf[0])
#
# return root <= self._sf[0]
# else:
# return RSOCConstraint(root, *self._sf)
else:
return NonconvexQuadraticConstraint(self, GE, other)
else:
return NotImplemented
# --------------------------------------
__all__ = api_end(_API_START, globals())