Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. from __future__ import division import warnings from typing import Tuple import numpy as np import scipy.sparse as sp from scipy import linalg as LA from cvxpy.atoms.affine.wraps import psd_wrap from cvxpy.atoms.atom import Atom from cvxpy.expressions.expression import Expression from cvxpy.interface.matrix_utilities import is_sparse from cvxpy.utilities.linalg import sparse_cholesky class CvxPyDomainError ( Exception ): class QuadForm ( Atom ): _allow_complex = True def __init__ ( self , x , P ) -> None : """Atom representing :math:`x^T P x`.""" super ( QuadForm , self ) . __init__ ( x , P ) def numeric ( self , values ): prod = values [ 1 ] . dot ( values [ 0 ]) if self . args [ 0 ] . is_complex (): quad = np . dot ( np . conj ( values [ 0 ]) . T , prod ) else : quad = np . dot ( np . transpose ( values [ 0 ]), prod ) return np . real ( quad ) def validate_arguments ( self ) -> None : super ( QuadForm , self ) . validate_arguments () n = self . args [ 1 ] . shape [ 0 ] if self . args [ 1 ] . shape [ 1 ] != n or self . args [ 0 ] . shape not in [( n , 1 ), ( n ,)]: raise ValueError ( "Invalid dimensions for arguments." ) if not self . args [ 1 ] . is_hermitian (): raise ValueError ( "Quadratic form matrices must be symmetric/Hermitian." ) def sign_from_args ( self ) -> Tuple [ bool , bool ]: """Returns sign (is positive, is negative) of the expression. return ( self . is_atom_convex (), self . is_atom_concave ()) def is_atom_convex ( self ) -> bool : """Is the atom convex? P = self . args [ 1 ] return P . is_constant () and P . is_psd () def is_atom_concave ( self ) -> bool : """Is the atom concave? P = self . args [ 1 ] return P . is_constant () and P . is_nsd () def is_atom_log_log_convex ( self ) -> bool : """Is the atom log-log convex? return True def is_atom_log_log_concave ( self ) -> bool : """Is the atom log-log concave? return False def is_incr ( self , idx ) -> bool : """Is the composition non-decreasing in argument idx? return ( self . args [ 0 ] . is_nonneg () and self . args [ 1 ] . is_nonneg ()) or \ ( self . args [ 0 ] . is_nonpos () and self . args [ 1 ] . is_nonneg ()) def is_decr ( self , idx ) -> bool : """Is the composition non-increasing in argument idx? return ( self . args [ 0 ] . is_nonneg () and self . args [ 1 ] . is_nonpos ()) or \ ( self . args [ 0 ] . is_nonpos () and self . args [ 1 ] . is_nonpos ()) def is_quadratic ( self ) -> bool : """Is the atom quadratic? return True def has_quadratic_term ( self ) -> bool : """Always a quadratic term. return True def is_pwl ( self ) -> bool : """Is the atom piecewise linear? return False def name ( self ) -> str : return " %s ( %s , %s )" % ( self . __class__ . __name__ , self . args [ 0 ], self . args [ 1 ]) def _grad ( self , values ): x = np . array ( values [ 0 ]) P = np . array ( values [ 1 ]) D = ( P + np . conj ( P . T )) @ x return [ sp . csc_matrix ( D . ravel ( order = "F" )) . T ] def shape_from_args ( self ) -> Tuple [ int , ... ]: return tuple () class SymbolicQuadForm ( Atom ): Symbolic form of QuadForm when quadratic matrix is not known (yet). def __init__ ( self , x , P , expr ) -> None : self . original_expression = expr super ( SymbolicQuadForm , self ) . __init__ ( x , P ) self . P = self . args [ 1 ] def get_data ( self ): return [ self . original_expression ] def _grad ( self , values ): raise NotImplementedError () def is_atom_concave ( self ) -> bool : return self . original_expression . is_atom_concave () def is_atom_convex ( self ) -> bool : return self . original_expression . is_atom_convex () def is_decr ( self , idx ) -> bool : return self . original_expression . is_decr ( idx ) def is_incr ( self , idx ) -> bool : return self . original_expression . is_incr ( idx ) def shape_from_args ( self ) -> Tuple [ int , ... ]: return self . original_expression . shape_from_args () def sign_from_args ( self ) -> Tuple [ bool , bool ]: return self . original_expression . sign_from_args () def is_quadratic ( self ) -> bool : return True def decomp_quad ( P , cond = None , rcond = None , lower = True , check_finite : bool = True ): Compute a matrix decomposition. Compute sgn, scale, M such that P = sgn * scale * dot(M, M.T). The strategy of determination of eigenvalue negligibility follows the pinvh contributions from the scikit-learn project to scipy. Parameters ---------- P : matrix or ndarray A real symmetric positive or negative (semi)definite input matrix cond, rcond : float, optional Cutoff for small eigenvalues. Singular values smaller than rcond * largest_eigenvalue are considered negligible. If None or -1, suitable machine precision is used (default). lower : bool, optional Whether the array data is taken from the lower or upper triangle of P. The default is to take it from the lower triangle. check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. The default is True; disabling may give a performance gain but may result in problems (crashes, non-termination) if the inputs contain infinities or NaNs. Returns ------- scale : float induced matrix 2-norm of P M1, M2 : 2d ndarray A rectangular ndarray such that P = scale * (dot(M1, M1.T) - dot(M2, M2.T)) if is_sparse ( P ): # TODO: consider using QDLDL instead, if available. try : sign , L , p = sparse_cholesky ( P ) if sign > 0 : return 1.0 , L [ p , :], np . empty (( 0 , 0 )) else : return 1.0 , np . empty (( 0 , 0 )), L [:, p ] except ( ValueError , ModuleNotFoundError ): P = np . array ( P . todense ()) # make dense (needs to happen for eigh). w , V = LA . eigh ( P , lower = lower , check_finite = check_finite ) if rcond is not None : cond = rcond if cond in ( None , - 1 ): t = V . dtype . char . lower () factor = { 'f' : 1e3 , 'd' : 1e6 } cond = factor [ t ] * np . finfo ( t ) . eps scale = max ( np . absolute ( w )) if scale == 0 : w_scaled = w else : w_scaled = w / scale maskp = w_scaled > cond maskn = w_scaled < - cond # TODO: allow indefinite quad_form if np . any ( maskp ) and np . any ( maskn ): warnings . warn ( "Forming a nonconvex expression quad_form(x, indefinite)." ) M1 = V [:, maskp ] * np . sqrt ( w_scaled [ maskp ]) M2 = V [:, maskn ] * np . sqrt ( - w_scaled [ maskn ]) return scale , M1 , M2