"""
Dual Problem Interface for Quadratically Constrained Quadratic Programming (QCQP).
This module provides interfaces for solving QCQP problems with shared projection
constraints using dual optimization methods. It includes both sparse and dense
implementations optimized for different matrix structures.
"""
__all__ = ["SparseSharedProjQCQP", "DenseSharedProjQCQP"]
import copy
from collections import OrderedDict
from typing import Any, cast
import numpy as np
import scipy.linalg as la
import scipy.sparse as sp
import sksparse.cholmod
from numpy.typing import ArrayLike
from dolphindes.types import ComplexArray, FloatNDArray
from ._base_qcqp import _SharedProjQCQP
[docs]
class SparseSharedProjQCQP(_SharedProjQCQP):
"""Sparse QCQP with projector-structured and optional general quadratic constraints.
This class specializes the generic _SharedProjQCQP implementation for the case
where ALL quadratic matrices (A0, A1, A2, each projector block, and any B_j)
are sparse. It supports:
1. A family of "shared" (projection-structured) constraints defined through
sparse projector matrices with shared sparsity structure.
2. An optional list of additional (general) quadratic equality constraints
parameterized by matrices B_j, vectors s_2j, and scalars c_2j.
Primal problem (maximization form):
maximize_x - x^† A0 x + 2 Re(x^† s0) + c0
subject to Re( - x^† A1 P_j A2 x + 2 x^† A2^† P_j^† s1 ) = 0 (shared)
Re( - x^† A2^† B_j A2 x + 2 x^† A2^† s_2j + c_2j ) = 0 (general)
Dual feasibility relies (heuristically) on at least one projector direction
such that A0 + λ A1 P_j A2 becomes PSD for sufficiently large λ. This is not
programmatically verified; users are responsible for supplying a suitable
projector set (the second multiplier is chosen for this role).
Parameters
----------
A0 : sp.csc_array
Objective quadratic matrix (Hermitian expected).
s0 : ArrayLike
Objective linear vector (complex allowed).
c0 : float
Objective constant term.
A1 : sp.csc_array
Left quadratic factor in projector constraints.
A2 : sp.csc_array
Right quadratic factor used in both projector and general constraints.
s1 : ArrayLike
Linear term coupled with projector constraints.
Plist : list[ArrayLike]
list of 2D arrays representing the projector matrices P_j.
B_j : list[sp.csc_array] | None
(Optional) list of general constraint middle matrices (between A2^† and A2).
s_2j : list[ArrayLike] | None
(Optional) list of linear term vectors for general constraints.
c_2j : ArrayLike | None
(Optional) array of constant terms for general constraints.
verbose : int, default 0
Verbosity level (≥1 prints preprocessing info).
Attributes
----------
A0, A1, A2 : scipy.sparse.csc_array
Stored sparse matrices
B_j : list[scipy.sparse.csc_array]
General constraint matrices (empty if none supplied).
s0, s1, s_2j : ComplexArray, ComplexArray, list[ComplexArray]
Complex vectors for objective / constraints.
c0 : float
Objective constant.
c_2j : FloatNDArray
Real constants for general constraints (length matches B_j).
Proj : :class:`dolphindes.util.Projectors`
dolphindes Projectors object representing all projector matrices P_j.
n_gen_constr : int
Number of general constraints (len(B_j)).
precomputed_As : list[sp.csc_array]
Symmetrized matrices [Sym(A1 P_j A2)] for projectors followed by
[Sym(A2^† B_j A2)] for general constraints (if any).
Fs : ComplexArray
Columns are A2^† P_j^† s1 (projector-only part used in derivatives).
Acho : sksparse.cholmod.Factor | None
Most recent numeric CHOLMOD factor (``None`` until the first
factorization). The symbolic analysis is held separately in
``_Acho_symbolic`` and reused across factorizations.
current_dual : float | None
Cached optimal dual value after solve_current_dual_problem().
current_lags : FloatNDArray | None
Cached Lagrange multipliers (projector first, then general).
current_grad : FloatNDArray | None
Gradient of dual at current_lags (if computed).
current_hess : FloatNDArray | None
Hessian of dual at current_lags (only when no general constraints).
current_xstar : ComplexArray | None
Primal maximizer associated with current_lags.
use_precomp : bool
Whether precomputation of constraint matrices/vectors is enabled.
verbose : int
Stored verbosity level.
Performance Notes
-----------------
- CHOLMOD symbolic analysis is performed once (via _initialize_Acho) based on
an example A(lags); subsequent factorizations reuse the sparsity pattern.
- Precomputation accelerates repeated evaluations for moderate constraint counts.
See Also
--------
DenseSharedProjQCQP : Dense analogue using LAPACK factorization.
_SharedProjQCQP : Base abstract class with core logic.
"""
# Symbolic CHOLMOD analysis (sparsity/fill pattern), computed once by
# _initialize_Acho during construction and reused by every numeric
# factorization. Declared here so its type is visible to static checkers.
_Acho_symbolic: "sksparse.cholmod.Factor"
def __repr__(self) -> str:
"""Return a concise string summary (size and projector count)."""
return (
f"SparseSharedProjQCQP of size {self.A0.shape[0]}^2 with "
f"{self.n_proj_constr} projectors."
)
def __deepcopy__(self, memo: dict[int, Any]) -> "SparseSharedProjQCQP":
"""Copy this instance."""
# custom __deepcopy__ because the CHOLMOD factorizations (the symbolic
# analysis and the cached numeric factors) are not pickle-able.
skip = {"Acho", "_Acho_symbolic", "_factor_cache"}
new_QCQP = SparseSharedProjQCQP.__new__(SparseSharedProjQCQP)
for name, value in self.__dict__.items():
if name not in skip:
setattr(new_QCQP, name, copy.deepcopy(value, memo))
new_QCQP._factor_cache = OrderedDict() # start with an empty factor cache
new_QCQP._initialize_Acho() # Recompute the symbolic Cholesky analysis.
# TODO: update Acho with current_lags if applicable
return new_QCQP
[docs]
def compute_precomputed_values(self) -> None:
"""Precompute constraint data then initialize symbolic factorization."""
super().compute_precomputed_values()
self._initialize_Acho()
def _initialize_Acho(self) -> sksparse.cholmod.Factor:
"""
Symbolically analyze sparsity/fill pattern for Cholesky factorization.
The symbolic analysis is stored in ``self._Acho_symbolic`` and reused by
every subsequent numeric factorization (which only needs the shared
sparsity pattern). ``self.Acho`` holds the most recent numeric factor and
is reset here since no numeric factorization has been computed yet.
Returns
-------
sksparse.cholmod.Factor
Analyzed (symbolic) factorization handle (``self._Acho_symbolic``).
"""
random_lags = np.random.rand(self.n_proj_constr + self.n_gen_constr)
A = self._get_total_A(random_lags)
assert sp.issparse(A)
A = sp.csc_array(A)
if self.verbose > 1:
print(
f"analyzing A of format and shape {type(A)}, {A.shape} "
f"and # of nonzero elements '{A.nnz}"
)
self._Acho_symbolic = sksparse.cholmod.analyze(A)
self.Acho = None
return self._Acho_symbolic
def _factorize(self, A: sp.csc_array) -> sksparse.cholmod.Factor:
"""
Return a fresh numeric Cholesky factor of A, reusing the symbolic analysis.
``Factor.cholesky`` (unlike ``cholesky_inplace``) returns a new factor
object, so distinct factorizations can coexist in the factor cache.
Parameters
----------
A : sp.csc_array
Current Hermitian (PSD / PD) system matrix.
Returns
-------
sksparse.cholmod.Factor
Numeric factorization of A.
"""
assert self._Acho_symbolic is not None
return self._Acho_symbolic.cholesky(sp.csc_array(A))
def _Acho_solve(self, b: ComplexArray) -> ComplexArray:
"""
Solve A x = b using the current CHOLMOD factorization.
Parameters
----------
b : ComplexArray
Right-hand side vector (or multiple RHS as columns).
Returns
-------
ComplexArray
Solution x = A^{-1} b.
"""
assert self.Acho is not None
return cast(ComplexArray, self.Acho.solve_A(b))
[docs]
def is_dual_feasible(self, lags: FloatNDArray) -> bool:
"""
Check PSD feasibility of A(lags) via attempted Cholesky factorization.
Parameters
----------
lags : FloatNDArray
Full Lagrange multiplier vector.
Returns
-------
bool
True if factorization succeeds (A is PSD), False otherwise.
"""
assert self._Acho_symbolic is not None
key = self._lags_key(lags)
if key in self._factor_cache:
# A(lags) was already factorized successfully -> positive definite
# (see the cache invariant in _store_factor). Count the hit as a use
# so the entry is not evicted before a genuinely older one.
self._factor_cache.move_to_end(key)
return True
A = self._get_total_A(lags)
assert sp.issparse(A)
A = sp.csc_array(A)
try:
factor = self._factorize(A)
# Accessing the factor forces the numeric decomposition to run, so a
# non-PD matrix raises here rather than lazily at solve time.
factor.L()
except sksparse.cholmod.CholmodNotPositiveDefiniteError:
return False
# Feasible: cache the factor so the imminent get_dual at the same lags
# (line-search objective evaluation) reuses it instead of refactorizing.
self._store_factor(key, A, factor)
return True
[docs]
class DenseSharedProjQCQP(_SharedProjQCQP):
"""Dense QCQP with projector-structured constraints.
Dense analogue of SparseSharedProjQCQP; uses scipy.linalg for
Cholesky factorization. Inherits full problem specification from
_SharedProjQCQP.
Parameters
----------
A0 : ArrayLike
Objective quadratic term.
s0 : ArrayLike
Objective linear term.
c0 : float
Objective constant.
A1 : ArrayLike
Left quadratic factor in projector constraints.
s1 : ArrayLike
Projector constraint linear term.
Plist : list[ArrayLike]
list of 2D arrays representing the projector matrices P_j.
A2 : ArrayLike | None, default None
Right quadratic factor (defaults to identity if None).
verbose : int, default 0
Verbosity level.
Notes
-----
- All quadratic matrices must be dense (or convertible) for this class.
- General constraints can be supplied via the base constructor if extended.
"""
def __init__(
self,
A0: ArrayLike,
s0: ArrayLike,
c0: float,
A1: ArrayLike,
s1: ArrayLike,
Plist: list[ArrayLike],
Pstruct: ArrayLike | None = None,
A2: ArrayLike | None = None,
B_j: list[ArrayLike] | None = None,
s_2j: list[ArrayLike] | None = None,
c_2j: ArrayLike | None = None,
verbose: int = 0,
):
if A2 is None:
n = int(np.asarray(s0).shape[0])
A2 = sp.eye_array(n, format="csc")
super().__init__(
A0,
s0,
c0,
A1,
A2,
s1,
Plist,
Pstruct,
B_j=B_j,
s_2j=s_2j,
c_2j=c_2j,
verbose=verbose,
)
def __repr__(self) -> str:
"""Return a concise string summary (size and projector count)."""
return (
f"DenseSharedProjQCQP of size {self.A0.shape[0]}^2 with "
f"{self.n_proj_constr} projectors."
)
def _factorize(self, A: ArrayLike) -> Any:
"""
Return a dense Cholesky factorization (cho_factor tuple) of A.
Parameters
----------
A : ArrayLike
Hermitian positive (semi)definite matrix to factor.
Returns
-------
Any
The ``scipy.linalg.cho_factor`` result ``(c, lower)``.
"""
return la.cho_factor(A)
def _Acho_solve(self, b: ComplexArray) -> ComplexArray:
"""
Solve A x = b using stored dense Cholesky factorization.
Parameters
----------
b : ComplexArray
Right-hand side vector (or stacked RHS matrix).
Returns
-------
ComplexArray
Solution x = A^{-1} b.
Notes
-----
``check_finite`` is disabled because it rescans the whole n^2 factor on
every call, which nearly doubles the cost of a solve. The factor comes from
``_factorize``, where ``cho_factor`` does validate, so non-finite values
are still caught once per factorization rather than once per solve -- and
this is called hundreds of times per factorization.
"""
return cast(
ComplexArray, la.cho_solve(self.Acho, b, check_finite=False)
)
[docs]
def is_dual_feasible(self, lags: FloatNDArray) -> bool:
"""
Check PSD feasibility of A(lags) via dense Cholesky attempt.
Parameters
----------
lags : FloatNDArray
Full Lagrange multiplier vector.
Returns
-------
bool
True if Cholesky succeeds, False otherwise.
"""
key = self._lags_key(lags)
if key in self._factor_cache:
# A(lags) was already factorized successfully -> positive definite
# (see the cache invariant in _store_factor). Count the hit as a use
# so the entry is not evicted before a genuinely older one.
self._factor_cache.move_to_end(key)
return True
A = self._get_total_A(lags)
try:
factor = self._factorize(A)
except la.LinAlgError:
return False
# Feasible: cache the factor so the imminent get_dual at the same lags
# (line-search objective evaluation) reuses it instead of refactorizing.
self._store_factor(key, A, factor)
return True