Source code for picos.expressions.exp_quadratic

# 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())