Coverage for picos/expressions/exp_quadratic.py : 82.11%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# ------------------------------------------------------------------------------
2# Copyright (C) 2019 Maximilian Stahlberg
3# Based on the original picos.expressions module by Guillaume Sagnol.
4#
5# This file is part of PICOS.
6#
7# PICOS is free software: you can redistribute it and/or modify it under the
8# terms of the GNU General Public License as published by the Free Software
9# Foundation, either version 3 of the License, or (at your option) any later
10# version.
11#
12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
15#
16# You should have received a copy of the GNU General Public License along with
17# this program. If not, see <http://www.gnu.org/licenses/>.
18# ------------------------------------------------------------------------------
20"""Implements :class:`QuadraticExpression`."""
22import operator
23import sys
24from collections import OrderedDict as odict
25from collections import namedtuple
26from types import MappingProxyType
28import cvxopt
29import cvxopt.cholmod
30import numpy
32from .. import glyphs
33from ..apidoc import api_end, api_start
34from ..caching import (borrow_cache, cached_property,
35 cached_selfinverse_unary_operator,
36 cached_unary_operator)
37from ..constraints import (ConicQuadraticConstraint, Constraint,
38 ConvexQuadraticConstraint,
39 NonconvexQuadraticConstraint)
40from ..legacy import deprecated
41from .data import convert_operands, cvx2np, should_be_sparse
42from .exp_affine import AffineExpression, ComplexAffineExpression, Constant
43from .expression import Expression, refine_operands, validate_prediction
44from .variables import BaseVariable
46_API_START = api_start(globals())
47# -------------------------------
50PSD_PERTURBATIONS = 3
51r"""Maximum number of singular quadratic form perturbations.
53PICOS uses NumPy either or CHOLMOD to compute a Cholesky decomposition of a
54positive semidefinite quadratic form :math:`Q`. Both libraries require
55:math:`Q` to be nonsingular, which is not a requirement on PICOS' end.
56If either library rejects :math:`Q`, then PICOS provides a sequence of
57:data:`PSD_PERTURBATIONS` perturbed quadratic forms :math:`Q + \epsilon I` for
58increasing :math:`\epsilon` until the perturbed matrix is found positive
59definite or until the largest :math:`\epsilon` was tested unsuccessfully.
61If this is zero, PICOS will only decompose quadratic forms that are nonsingular.
62If this is one, then only the largest epsilon is tested.
63"""
66class QuadraticExpression(Expression):
67 """A scalar quadratic expression of the form :math:`x^TQx + a^Tx + b`."""
69 # --------------------------------------------------------------------------
70 # Initialization and factory methods.
71 # --------------------------------------------------------------------------
73 def __init__(
74 self, string, quadraticPart={}, affinePart=AffineExpression.zero(),
75 scalarFactors=None, copyDecomposition=None):
76 r"""Initialize a scalar quadratic expression.
78 This constructor is meant for internal use. As a user, you will most
79 likely want to build expressions starting from
80 :mod:`~picos.expressions.variables` or a :func:`~picos.Constant`.
82 :param str string: A symbolic string description.
83 :param quadraticPart: A :class:`dict` mapping PICOS variable pairs to
84 CVXOPT matrices. Each entry :math:`(x, y) \mapsto A_{xy}` represents
85 a quadratic form
86 :math:`\operatorname{vec_x}(x)^T A_{xy} \operatorname{vec_y}(y)`
87 where :math:`\operatorname{vec_x}` and :math:`\operatorname{vec_y}`
88 refer to the isometric real
89 :mod:`~.picos.expressions.vectorizations` that are used by PICOS
90 internally to store the variables :math:`x` and :math:`y`,
91 respectively. The quadratic part :math:`Q` of the expression is then
92 given as :math:`Q = \sum_{x, y} A_{xy}`.
93 :param affinePart: The affine part :math:`a^Tx + b` of the expression.
94 :type affinePart: ~picos.expressions.AffineExpression
95 :param scalarFactors: A pair :math:`(u, v)` with both :math:`u` and
96 :math:`v` scalar real affine expressions representing a known
97 :math:`x^TQx + a^Tx + b = uv` factorization of the expression.
98 :type factorization:
99 tuple(~picos.expressions.AffineExpression)
100 :param copyDecomposition: Another quadratic expression with equal
101 quadratic part whose quadratic part decomposition shall be copied.
102 :type copyDecomposition:
103 ~picos.expressions.QuadraticExpression
104 """
105 # Ensure correct usage within PICOS.
106 assert isinstance(quadraticPart, dict) and all(
107 isinstance(xy, (tuple, list)) and len(xy) == 2
108 and isinstance(xy[0], BaseVariable)
109 and isinstance(xy[1], BaseVariable)
110 and isinstance(A, (cvxopt.matrix, cvxopt.spmatrix))
111 and A.size == (xy[0].dim, xy[1].dim)
112 for xy, A in quadraticPart.items())
113 assert isinstance(affinePart, ComplexAffineExpression) \
114 and len(affinePart) == 1
115 assert not copyDecomposition \
116 or isinstance(copyDecomposition, QuadraticExpression)
117 assert not scalarFactors or (
118 isinstance(scalarFactors, (tuple, list)) and len(scalarFactors) == 2
119 and all(isinstance(x, AffineExpression) for x in scalarFactors)
120 and all(len(x) == 1 for x in scalarFactors))
122 Expression.__init__(self, "Quadratic Expression", string)
124 # Check affine part.
125 affinePart = affinePart.refined
126 if not isinstance(affinePart, AffineExpression):
127 raise NotImplementedError(
128 "PICOS does not support complex-valued quadratic expressions "
129 "and the affine part of {} is not real.".format(self._symbStr))
130 elif len(affinePart) != 1:
131 raise TypeError("Affine part of {} is not scalar."
132 .format(self._symbStr))
134 # Store quadratic part.
135 self._quads = {}
136 for xy, A in quadraticPart.items():
137 # Ignore coefficients that are already zero.
138 if not A:
139 continue
141 if xy[0] is xy[1]:
142 # Store the unique symmetric form for two reasons:
143 # - Hermitian forms concerning the original variables are
144 # supplied as quadratic forms concerning the isometric real
145 # vectorizations of the variables. These matrices are in
146 # general still complex, but their symmetric form is not.
147 # - Coefficients that cancel out are set to numeric zeros and
148 # can be converted to structural zeros by _sparse_quads.
149 A = 0.5*(A + A.T)
150 if hash(xy[1]) < hash(xy[0]):
151 # Store only one coefficient per unordered pair.
152 xy, A = (xy[1], xy[0]), A.T
154 if xy in self._quads:
155 self._quads[xy] = self._quads[xy] + A
156 else:
157 self._quads[xy] = A
159 # Remove coefficients of zero.
160 self._quads = {xy: A for xy, A in self._quads.items() if A}
162 # Check if coefficients are real.
163 for (x, y), A in self._quads.items():
164 if A.typecode == "z":
165 if any(A.imag()):
166 raise NotImplementedError(
167 "PICOS does not support complex-valued quadratic "
168 "expressions and the quadratic part of {} for variables"
169 " {} and {} is not real or equivalent to a real form."
170 .format(self._symbStr, x.string, y.string))
172 self._quads[xy] = A.real()
174 # Store affine part.
175 self._aff = affinePart
177 # Store a known factorization into real scalars.
178 if scalarFactors:
179 a, b = scalarFactors
180 if a.equals(b):
181 self._sf = (a, a) # Speedup future equality checks.
182 else:
183 self._sf = (a, b)
184 else:
185 self._sf = None
187 # Copy a decomposition from another quadratic expression.
188 if copyDecomposition is not None:
189 borrow_cache(self, copyDecomposition, ("_Q_and_x", "L", "quadroot"))
191 # --------------------------------------------------------------------------
192 # Abstract method implementations and method overridings, except _predict.
193 # --------------------------------------------------------------------------
195 @cached_unary_operator
196 def _get_refined(self):
197 if not self._quads:
198 return self._aff.refined.renamed(self._symbStr)
199 else:
200 return self
202 Subtype = namedtuple("Subtype", (
203 "convex", "concave", "qterms", "rank", "haslin", "factors"))
205 @cached_unary_operator
206 def _get_subtype(self):
207 convex = self.convex
208 concave = False if convex and self._quads else self.concave
209 qterms = self.num_quad_terms
210 haslin = not self._aff.constant
212 try:
213 if convex:
214 rank = self.rank
215 elif concave:
216 rank = (-self).rank
217 else:
218 rank = None
219 except ValueError: # No quadratic part.
220 rank = 0
222 factors = len(set(self._sf)) if self._sf else 0
224 return self.Subtype(convex, concave, qterms, rank, haslin, factors)
226 def _get_value(self):
227 value = self._aff._get_value()
229 for xy, A in self._quads.items():
230 x, y = xy
231 xVal = x._get_internal_value()
232 yVal = y._get_internal_value()
233 value = value + xVal.T * A * yVal
235 return value
237 @cached_unary_operator
238 def _get_mutables(self):
239 return self._aff._get_mutables().union(
240 var for vars in self._quads for var in vars)
242 def _is_convex(self):
243 try:
244 self.L
245 except ValueError:
246 return True # Expression is affine.
247 except ArithmeticError:
248 return False
249 else:
250 return True
252 def _is_concave(self):
253 try:
254 (-self).L
255 except ValueError:
256 return True # Expression is affine.
257 except ArithmeticError:
258 return False
259 else:
260 return True
262 def _replace_mutables(self, mapping):
263 name_map = {old.name: new.name for old, new in mapping.items()}
265 string = self.string
266 if isinstance(string, glyphs.GlStr):
267 string = string.reglyphed(name_map)
269 quads = {(mapping[var1], mapping[var2]): coef
270 for (var1, var2), coef in self._quads.items()}
272 if not self._sf:
273 sf = None
274 elif self._sf[0] is self._sf[1]:
275 sf = (self._sf[0]._replace_mutables(mapping),)*2
276 else:
277 sf = tuple(f._replace_mutables(mapping) for f in self._sf)
279 return QuadraticExpression(
280 string, quads, self.aff._replace_mutables(mapping), sf)
282 def _freeze_mutables(self, freeze):
283 # TODO: Allow freezing of quadratic expressions.
284 raise NotImplementedError(
285 "Partially freezing quadratic expressions is not yet supported.")
287 # --------------------------------------------------------------------------
288 # Python special method implementations, except constraint-creating ones.
289 # --------------------------------------------------------------------------
291 def __len__(self):
292 # Faster version that overrides Expression.__len__.
293 return 1
295 @convert_operands(sameShape=True)
296 def __add__(self, other):
297 """Denote addition from the right hand side."""
298 if isinstance(other, QuadraticExpression):
299 string = glyphs.clever_add(self.string, other.string)
301 quadPart = {}
302 varPairs = set(self._quads.keys()).union(other._quads.keys())
303 for vars in varPairs:
304 if vars in self._quads:
305 if vars in other._quads:
306 quadPart[vars] = self._quads[vars] + other._quads[vars]
307 else:
308 quadPart[vars] = self._quads[vars]
309 else:
310 quadPart[vars] = other._quads[vars]
312 affPart = self._aff + other._aff
314 return QuadraticExpression(string, quadPart, affPart)
315 elif isinstance(other, AffineExpression):
316 string = glyphs.clever_add(self.string, other.string)
318 return QuadraticExpression(
319 string, self._quads, self._aff + other, copyDecomposition=self)
320 else:
321 return NotImplemented
323 @convert_operands(sameShape=True)
324 def __radd__(self, other):
325 """Denote addition from the left hand side."""
326 if isinstance(other, QuadraticExpression):
327 assert False, "Impossible case."
328 elif isinstance(other, AffineExpression):
329 result = self.__add__(other)
330 # NOTE: __add__ always creates a fresh expression.
331 result._symbStr = glyphs.clever_add(other.string, self.string)
332 return result
333 else:
334 return NotImplemented
336 @convert_operands(sameShape=True)
337 def __sub__(self, other):
338 """Denote substraction from the right hand side."""
339 if isinstance(other, QuadraticExpression):
340 string = glyphs.clever_sub(self.string, other.string)
342 quadPart = {}
343 varPairs = set(self._quads.keys()).union(other._quads.keys())
344 for vars in varPairs:
345 if vars in self._quads:
346 if vars in other._quads:
347 quadPart[vars] = self._quads[vars] - other._quads[vars]
348 else:
349 quadPart[vars] = self._quads[vars]
350 else:
351 quadPart[vars] = -other._quads[vars]
353 affPart = self._aff - other._aff
355 return QuadraticExpression(string, quadPart, affPart)
356 elif isinstance(other, AffineExpression):
357 string = glyphs.clever_sub(self.string, other.string)
359 return QuadraticExpression(
360 string, self._quads, self._aff - other, copyDecomposition=self)
361 else:
362 return NotImplemented
364 @convert_operands(sameShape=True)
365 def __rsub__(self, other):
366 """Denote substraction with self on the right hand side."""
367 if isinstance(other, QuadraticExpression):
368 assert False, "Impossible case."
369 elif isinstance(other, AffineExpression):
370 result = (-self).__add__(other)
371 # NOTE: __add__ always creates a fresh expression.
372 result._symbStr = glyphs.clever_sub(other.string, self.string)
373 return result
374 else:
375 return NotImplemented
377 @cached_selfinverse_unary_operator
378 def __neg__(self):
379 string = glyphs.clever_neg(self.string)
380 quadPart = {xy: -A for xy, A in self._quads.items()}
381 affPart = -self._aff
382 factors = (self._sf[0], -self._sf[1]) if self._sf else None
384 return QuadraticExpression(string, quadPart, affPart, factors)
386 @convert_operands(scalarRHS=True)
387 def __mul__(self, other):
388 """Denote scaling from the right hand side."""
389 if isinstance(other, AffineExpression):
390 if not other.constant:
391 raise NotImplementedError("You may only multiply a PICOS "
392 "quadratic expression with a constant term.")
394 factor = other.value
395 string = glyphs.clever_mul(self.string, other.string)
396 quadPart = {xy: factor*A for xy, A in self._quads.items()}
397 affPart = self._aff*other
399 if self._sf:
400 r = factor**0.5
401 if self._sf[0] is self._sf[1]:
402 factors = (r*self._sf[0],)*2
403 else:
404 factors = (r*self._sf[0], r*self._sf[1])
405 else:
406 factors = None
408 return QuadraticExpression(string, quadPart, affPart, factors)
409 else:
410 return NotImplemented
412 @convert_operands(scalarRHS=True)
413 def __rmul__(self, other):
414 """Denote scaling from the left hand side."""
415 if isinstance(other, AffineExpression):
416 result = self.__mul__(other)
417 # NOTE: __mul__ always creates a fresh expression.
418 result._symbStr = glyphs.clever_mul(other.string, self.string)
419 return result
420 else:
421 return NotImplemented
423 @convert_operands(scalarRHS=True)
424 def __truediv__(self, other):
425 """Denote division by a constant scalar."""
426 if isinstance(other, AffineExpression):
427 if not other.constant:
428 raise NotImplementedError("You may only divide a PICOS "
429 "quadratic expression by a constant term.")
431 result = self.__mul__(1.0 / other.value)
432 # NOTE: __mul__ always creates a fresh expression.
433 result._symbStr = glyphs.div(self.string, other.string)
434 return result
435 else:
436 return NotImplemented
438 # --------------------------------------------------------------------------
439 # Decomposition related methods.
440 # --------------------------------------------------------------------------
442 @cached_property
443 def _sparse_quads(self):
444 """Quadratic coefficients (re-)cast to sparse matrices."""
445 return {xy: cvxopt.sparse(M) for xy, M in self._quads.items()}
447 def _Q_and_x_or_A_and_y(self, augmentedMatrix):
448 """Either :attr:`_Q_and_x` or :attr:`_A_and_y`."""
449 # Find occurences of the scalar entries of the variable's isometric real
450 # vectorization in the quadratic part.
451 nonzero = {}
452 for (x, y), A in self._sparse_quads.items():
453 nonzero.setdefault(x, set())
454 nonzero.setdefault(y, set())
456 nonzero[x].update(A.I)
457 nonzero[y].update(A.J)
459 # For the augmented matrix, find additional linear occurences.
460 if augmentedMatrix:
461 for x, a in self._aff._sparse_linear_coefs.items():
462 nonzero.setdefault(x, set())
463 nonzero[x].update(a.I)
465 nonzero = odict((var, sorted(nonzero[var]))
466 for var in sorted(nonzero, key=lambda v: v.id))
468 # Compute each (multidimensional) variable's offset in Q/A.
469 offset, offsets = 0, {}
470 for var, nz in nonzero.items():
471 offsets[var] = offset
472 offset += len(nz)
474 # Compute Q.
475 V, I, J = [], [], []
476 for (x, y), A in self._sparse_quads.items():
477 B = A[nonzero[x], nonzero[y]]
479 V.extend(B.V)
480 I.extend(offsets[x] + i for i in B.I)
481 J.extend(offsets[y] + j for j in B.J)
483 # Turn Q into A if requested.
484 if augmentedMatrix:
485 for x, a in self._aff._sparse_linear_coefs.items():
486 b = a[nonzero[x]]
488 V.extend(b.V)
489 I.extend(offsets[x] + i for i in b.I)
490 J.extend([offset]*len(b.J))
492 V.append(self._aff._constant_coef[0])
493 I.append(offset)
494 J.append(offset)
496 offset += 1
498 # Finalize Q/A.
499 Q = cvxopt.spmatrix(V, I, J, (offset,)*2, tc="d")
500 Q = 0.5*(Q + Q.T)
502 if not Q:
503 if augmentedMatrix:
504 raise ValueError(
505 "The expression {} is zero.".format(self._symbStr))
506 else:
507 raise ValueError(
508 "The quadratic part of {} is zero.".format(self._symbStr))
510 # Compute x.
511 x = None
512 for var, nz in nonzero.items():
513 n, r = len(nz), var.dim
514 F = cvxopt.spmatrix([1]*n, range(n), nz, (n, r))
515 y = AffineExpression(
516 glyphs.Fn("F")(var.string), n, coefficients={var: F})
517 x = y if x is None else x // y
519 # Turn x into y if requested.
520 if augmentedMatrix:
521 x = AffineExpression.from_constant(1) if x is None else x // 1
523 return Q, x
525 @cached_property
526 def _Q_and_x(self):
527 """Pair containing :attr:`Q` and :attr:`x`."""
528 return self._Q_and_x_or_A_and_y(augmentedMatrix=False)
530 @cached_property
531 def _A_and_y(self):
532 """Pair containing :attr:`R` and :attr:`y`."""
533 return self._Q_and_x_or_A_and_y(augmentedMatrix=True)
535 @property
536 def Q(self):
537 """The coefficient matrix :math:`Q` of the expression, condensed.
539 This equals the :math:`Q` of the quadratic expression
540 :math:`x^TQx + a^Tx + b` but with zero rows and columns removed.
542 The vector :attr:`x` is a condensed version of :math:`x` that refers to
543 this matrix, so that ``q.x.T*q.Q*q.x`` equals the quadratic part
544 :math:`x^TQx` of the expression ``q``.
546 :raises ValueError: When the quadratic part is zero.
547 """
548 return self._Q_and_x[0]
550 @property
551 def x(self):
552 """The stacked variable vector :math:`x` of the expression, condensed.
554 This equals the :math:`x` of the quadratic expression
555 :math:`x^TQx + a^Tx + b` but entries corresponding to zero rows and
556 columns in :attr:`Q` are removed, so that ``q.x.T*q.Q*q.x`` equals the
557 :math:`x^TQx` part of the expression ``q``.
559 :raises ValueError: When the quadratic part is zero.
560 """
561 return self._Q_and_x[1]
563 @property
564 def A(self):
565 r"""An affine-augmented quadratic coefficient matrix, condensed.
567 For a quadratic expression :math:`x^TQx + a^Tx + b`, this is
568 :math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}`
569 but with zero rows and columns removed.
571 The vector :attr:`y` is a condensed version of :math:`x` that refers to
572 this matrix, so that ``q.y.T*q.A*q.y`` equals the expression ``q``.
574 :raises ValueError: When the expression is zero.
575 """
576 return self._A_and_y[0]
578 @property
579 def y(self):
580 """See :attr:`A`.
582 :raises ValueError: When the expression is zero.
583 """
584 return self._A_and_y[1]
586 def _L_or_M(self, augmentedMatrix):
587 """Either :attr:`L` or :attr:`M`."""
588 Q = self.A if augmentedMatrix else self.Q
589 n = Q.size[0]
591 # Define a sequence of small numbers for perturbation.
592 # TODO: These numbers are empirical; find worst case values.
593 spread = [abs(v) for v in Q.V if v]
594 minEps = 1 * min(spread) * sys.float_info.epsilon
595 maxEps = n * max(spread) * sys.float_info.epsilon
596 epsilons = set([0.0])
597 for i in range(PSD_PERTURBATIONS):
598 epsilons.add(minEps + i * (maxEps - minEps) / PSD_PERTURBATIONS)
599 epsilons = sorted(epsilons)
601 # Check if Q is really sparse.
602 sparse = should_be_sparse(Q.size, len(Q))
604 # Attempt to find L.
605 L = None
606 if sparse: # Use CHOLMOD.
607 F = None
608 for eps in epsilons:
609 P = Q + cvxopt.spdiag([eps]*n) if eps else Q # Perturbed Q.
611 if F is None:
612 F = cvxopt.cholmod.symbolic(P)
614 try:
615 cvxopt.cholmod.numeric(P, F)
616 except ArithmeticError:
617 pass
618 else:
619 L = cvxopt.cholmod.getfactor(F)
620 break
621 else: # Use NumPy.
622 Q = cvx2np(Q)
623 for eps in epsilons:
624 P = Q + numpy.diag([eps]*n) if eps else Q # Perturbed Q.
626 try:
627 L = numpy.linalg.cholesky(P)
628 except numpy.linalg.LinAlgError:
629 pass
630 else:
631 break
633 if L is None:
634 if not PSD_PERTURBATIONS:
635 raise ArithmeticError("The expression {} has an (augmented) "
636 "quadratic form that is either not positive semidefinite "
637 "or singular. Try enabling PSD_PERTURBATIONS.")
638 elif augmentedMatrix:
639 raise ArithmeticError("The expression {} is not a squared norm "
640 "as its {} representation is not positive semidefinite."
641 .format(self._symbStr, glyphs.matrix(glyphs.vertcat(
642 glyphs.horicat("Q", glyphs.div("a", 2)),
643 glyphs.horicat(glyphs.div(glyphs.transp("a"), 2), "b")
644 ))))
645 else:
646 raise ArithmeticError("The expression {} is not convex as its "
647 "quadratic part is not positive semidefinite."
648 .format(self._symbStr))
650 # Remove near-zeros in the order of n·sqrt(maxEps).
651 cutoff = 2 * n * maxEps**0.5
652 if sparse:
653 VIJ = list(t for t in zip(L.V, L.I, L.J) if abs(t[0]) > cutoff)
654 if VIJ: # Don't zero the matrix.
655 V, I, J = zip(*VIJ)
656 L = cvxopt.spmatrix(V, I, J, L.size)
657 else:
658 mask = abs(L) <= cutoff
659 if not numpy.all(mask): # Don't zero the matrix.
660 L[mask] = 0
661 L = cvxopt.sparse(cvxopt.matrix(L))
663 return L
665 @cached_property
666 def L(self):
667 """The :math:`L` of an :math:`LL^T` Cholesky decomposition of :attr:`Q`.
669 :returns: A CVXOPT lower triangular sparse matrix.
671 :raises ValueError: When the quadratic part is zero.
672 :raises ArithmeticError: When the expression is not convex, that is when
673 the matrix :attr:`Q` is not numerically positive semidefinite.
674 """
675 return self._L_or_M(augmentedMatrix=False)
677 @cached_property
678 def M(self):
679 """The :math:`M` of an :math:`MM^T` Cholesky decomposition of :attr:`A`.
681 :returns: A CVXOPT lower triangular sparse matrix.
683 :raises ValueError: When the expression is zero.
684 :raises ArithmeticError: When the expression can not be written as a
685 squared norm, that is when the extended matrix :attr:`A` is not
686 numerically positive semidefinite.
687 """
688 return self._L_or_M(augmentedMatrix=True)
690 @cached_property
691 def quadroot(self):
692 r"""Affine expression whose squared norm equals :math:`x^TQx`.
694 For a convex quadratic expression ``q``, this is equal to the vector
695 ``q.L.T*q.x`` with zero rows removed.
697 :Construction:
699 Let :math:`x^TQx` be the quadratic part of the expression with :math:`Q`
700 positive semidefinite and :math:`Q = LL^T` a Cholesky decomposition.
701 Then,
703 .. math::
704 x^TQx
705 &= x^TLL^Tx \\
706 &= (L^Tx)^TL^Tx \\
707 &= \langle L^Tx, L^Tx \rangle \\
708 &= \lVert L^Tx \rVert^2.
710 Note that removing zero rows from :math:`L^Tx` does not affect the norm.
712 :raises ValueError: When the quadratic part is zero.
713 :raises ArithmeticError: When the expression is not convex, that is when
714 the quadratic part is not numerically positive semidefinite.
715 """
716 L = self.L
717 n = L.size[0]
718 nz = sorted(set(L.J))
719 nnz = len(nz)
721 # F removes rows corresponding to zero columns in L.
722 F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n))
724 result = (F*L.T)*self.x
725 # NOTE: AffineExpression.__mul__ always creates a fresh expression.
726 result._symbStr = glyphs.Fn("quadroot({})")(self.string)
727 return result
729 @cached_property
730 def fullroot(self):
731 r"""Affine expression whose squared norm equals the expression.
733 For a convex quadratic expression ``q``, this is equal to the vector
734 ``q.M.T*q.y`` with zero rows removed.
736 :Construction:
738 For a quadratic expression :math:`x^TQx + a^Tx + b` with
739 :math:`A = \begin{bmatrix}Q&\frac{a}{2}\\\frac{a^T}{2}&b\end{bmatrix}`
740 positive semidefinite, let :math:`A = MM^T` be a Cholesky decomposition.
741 Let further :math:`y = \begin{bmatrix}x^T & 1\end{bmatrix}^T`. Then,
743 .. math::
744 x^TQx + a^Tx + b
745 &= y^TAy \\
746 &= y^TMM^Ty \\
747 &= (M^Ty)^TM^Ty \\
748 &= \langle M^Ty, M^Ty \rangle \\
749 &= \lVert M^Ty \rVert^2.
751 Note that removing zero rows from :math:`M^Ty` does not affect the norm.
753 :raises ValueError: When the expression is zero.
754 :raises ArithmeticError: When the expression is not convex, that is when
755 the quadratic part is not numerically positive semidefinite.
756 """
757 M = self.M
758 n = M.size[0]
759 nz = sorted(set(M.J))
760 nnz = len(nz)
762 # F removes rows corresponding to zero columns in M.
763 F = cvxopt.spmatrix([1.0]*nnz, range(nnz), nz, (nnz, n))
765 result = (F*M.T)*self.y
766 # NOTE: AffineExpression.__mul__ always creates a fresh expression.
767 result._symbStr = glyphs.Fn("fullroot({})")(self.string)
768 return result
770 @property
771 def rank(self):
772 """The length of the vector :attr:`quadroot`.
774 Up to numerical considerations, this is the rank of the (convex)
775 quadratic coefficient matrix :math:`Q` of the expression.
777 :raises ArithmeticError: When the expression is not convex, that is when
778 the quadratic part is not numerically positive semidefinite.
779 """
780 try:
781 return len(self.quadroot)
782 except ValueError:
783 return 0 # No quadratic part.
785 @cached_property
786 def num_quad_terms(self):
787 """The number of terms in the simplified quadratic form."""
788 num_terms = 0
789 for A in self._sparse_quads.values():
790 for i, j in zip(A.I, A.J):
791 if i >= j:
792 num_terms += 1
793 return num_terms
795 @property
796 def is_squared_norm(self):
797 r"""Whether the expression can be written as a squared norm.
799 If this is :obj:`True`, then the there is a coefficient vector
800 :math:`c` such that
802 .. math::
803 \lVert c^T \begin{bmatrix} x \\ 1 \end{bmatrix} \rVert^2
804 = x^TQx + a^Tx + b.
806 If the expression is also nonzero, then :attr:`fullroot` is
807 :math:`c^T \begin{bmatrix} x^T & 1 \end{bmatrix}^T` with zero entries
808 removed.
809 """
810 try:
811 self.M
812 except ValueError:
813 assert self.is0
814 return True
815 except ArithmeticError:
816 return False
817 else:
818 return True
820 # --------------------------------------------------------------------------
821 # Properties and functions that describe the expression.
822 # --------------------------------------------------------------------------
824 @property
825 @deprecated("2.2", "This property will be removed in a future release.")
826 def quadratic_forms(self):
827 """The quadratic forms as a map from variable pairs to sparse matrices.
829 .. warning::
831 Do not modify the returned matrices.
832 """
833 return MappingProxyType(self._sparse_quads)
835 @property
836 def is0(self):
837 """Whether the quadratic expression is zero."""
838 return not self._quads and self._aff.is0
840 @property
841 def scalar_factors(self):
842 r"""Decomposition into scalar real affine expressions.
844 If the expression is known to be equal to :math:`ab` for scalar real
845 affine expressions :math:`a` and :math:`b`, this is the pair ``(a, b)``.
846 Otherwise, this is :obj:`None`.
848 Note that if :math:`a = b`, then also ``a is b`` in the returned tuple.
849 """
850 return self._sf
852 # --------------------------------------------------------------------------
853 # Methods and properties that return modified copies.
854 # --------------------------------------------------------------------------
856 @cached_property
857 def quad(self):
858 """Quadratic part :math:`x^TQx` of the quadratic expression."""
859 factors = self._sf if self._sf and self._aff.is0 else None
861 return QuadraticExpression(glyphs.quadpart(self._symbStr), self._quads,
862 scalarFactors=factors, copyDecomposition=self)
864 @cached_property
865 def aff(self):
866 """Affine part :math:`a^Tx + b` of the quadratic expression."""
867 return self._aff.renamed(glyphs.affpart(self._symbStr))
869 @cached_property
870 def lin(self):
871 """Linear part :math:`a^Tx` of the quadratic expression."""
872 return self._aff.lin.renamed(glyphs.linpart(self._symbStr))
874 @cached_property
875 def cst(self):
876 """Constant part :math:`b` of the quadratic expression."""
877 return self._aff.cst.renamed(glyphs.cstpart(self._symbStr))
879 # --------------------------------------------------------------------------
880 # Constraint-creating operators, and _predict.
881 # --------------------------------------------------------------------------
883 @classmethod
884 def _predict(cls, subtype, relation, other):
885 assert isinstance(subtype, cls.Subtype)
887 if issubclass(other.clstype, QuadraticExpression):
888 raise NotImplementedError("Constraint outcome prediction not "
889 "supported when comparing two quadratic expressions.")
891 if relation == operator.__le__:
892 if issubclass(other.clstype, AffineExpression) \
893 and other.subtype.dim == 1:
894 if subtype.convex:
895 return ConvexQuadraticConstraint.make_type(
896 qterms=subtype.qterms, rank=subtype.rank,
897 haslin=(subtype.haslin or not other.subtype.constant))
898 else:
899 return NonconvexQuadraticConstraint.make_type(
900 qterms=subtype.qterms)
901 elif relation == operator.__ge__:
902 if issubclass(other.clstype, AffineExpression) \
903 and other.subtype.dim == 1:
904 if subtype.concave:
905 return ConvexQuadraticConstraint.make_type(
906 qterms=subtype.qterms, rank=subtype.rank,
907 haslin=(subtype.haslin or not other.subtype.constant))
908 elif subtype.factors and other.subtype.constant \
909 and other.subtype.nonneg:
910 # See __ge__ for the variants below.
912 # Variant 1:
913 return ConicQuadraticConstraint.make_type(
914 qterms=subtype.qterms, conic_argdim=1,
915 rotated=(subtype.factors == 2))
917 # Variant 2:
918 # return RSOCConstraint.make_type(argdim=1)
920 # Variant 3:
921 # if subtype.factors == 1:
922 # return AffineConstraint.make_type(dim=1, eq=False)
923 # else:
924 # return RSOCConstraint.make_type(argdim=1)
925 else:
926 return NonconvexQuadraticConstraint.make_type(
927 qterms=subtype.qterms)
929 return NotImplemented
931 @convert_operands(sameShape=True)
932 @validate_prediction
933 @refine_operands()
934 def __le__(self, other):
935 LE = Constraint.LE
937 if isinstance(other, QuadraticExpression):
938 le0 = self - other # Unpredictable outcome.
939 if le0.convex:
940 return ConvexQuadraticConstraint(self, LE, other)
941 elif self.is_squared_norm and other.scalar_factors:
942 return ConicQuadraticConstraint(self, LE, other)
943 else:
944 return NonconvexQuadraticConstraint(self, LE, other)
945 elif isinstance(other, AffineExpression):
946 if self.convex:
947 return ConvexQuadraticConstraint(self, LE, other)
948 else:
949 return NonconvexQuadraticConstraint(self, LE, other)
950 else:
951 return NotImplemented
953 @convert_operands(sameShape=True)
954 @validate_prediction
955 @refine_operands()
956 def __ge__(self, other):
957 GE = Constraint.GE
959 if isinstance(other, QuadraticExpression):
960 le0 = other - self # Unpredictable outcome.
961 if le0.convex:
962 return ConvexQuadraticConstraint(self, GE, other)
963 elif other.is_squared_norm and self.scalar_factors:
964 return ConicQuadraticConstraint(self, GE, other)
965 else:
966 return NonconvexQuadraticConstraint(self, GE, other)
967 elif isinstance(other, AffineExpression):
968 if self.concave:
969 return ConvexQuadraticConstraint(self, GE, other)
970 elif self._sf and other.constant and other.value >= 0:
971 # Handle the case of c <= x*x and c <= x*y with c pos. constant.
972 root = Constant(glyphs.sqrt(other.string), other.value**0.5)
974 # Variant 1: Return a ConicQuadraticConstraint that can be
975 # reformulated to an (R)SOCConstraint. This is consistent
976 # with all other quadratic expression comparison outcomes.
977 other = QuadraticExpression(other.string, affinePart=other,
978 scalarFactors=(root, root)) # Unrefined QE.
979 return ConicQuadraticConstraint(self, GE, other)
981 # Variant 2: Return an RSOCConstraint. This is somewhat
982 # consistent with Norm returning an SOCConstraint and it is
983 # also legacy behavior. However this may not be what the
984 # user wants, because it silently forces x, y >= 0!
985 # return RSOCConstraint(root, self._sf[0], self._sf[1],
986 # customString = glyphs.le(other.string, self.string))
988 # Variant 3: Return either an RSOCConstraint or an affine
989 # constraint that equals an SOCConstraint. This is faster
990 # than Variant 2 but even more confusing, because the
991 # constraint would transform from e.g. 1 <= x*x to
992 # sqrt(1) <= x without any mention of cones.
993 # if self._sf[0] is self._sf[1]:
994 # # NOTE: The following is not allowed by SOCConstraint
995 # # since this is effectively an affine constraint.
996 # # return SOCConstraint(root, self._sf[0])
997 #
998 # return root <= self._sf[0]
999 # else:
1000 # return RSOCConstraint(root, *self._sf)
1001 else:
1002 return NonconvexQuadraticConstraint(self, GE, other)
1003 else:
1004 return NotImplemented
1007# --------------------------------------
1008__all__ = api_end(_API_START, globals())