tessl/pypi-qpsolvers

Quadratic programming solvers in Python with a unified API

—

Pending

Overview

Eval results

Files

Data Structures

Name: tessl/pypi-qpsolvers
Author: tessl

Core classes for representing quadratic programs, solutions, and active constraint sets. These provide structured data handling with validation, optimality checking, and comprehensive solution information.

Capabilities

Problem Representation

The Problem class provides a structured representation of quadratic programs with validation and metrics.

class Problem:
    """
    Data structure describing a quadratic program.
    
    The QP is formulated as:
    minimize    1/2 x^T P x + q^T x
    subject to  G x <= h
                A x = b
                lb <= x <= ub
    """
    
    P: Union[np.ndarray, spa.csc_matrix]  # Symmetric cost matrix
    q: np.ndarray                         # Cost vector
    G: Optional[Union[np.ndarray, spa.csc_matrix]] = None  # Inequality matrix
    h: Optional[np.ndarray] = None                         # Inequality vector
    A: Optional[Union[np.ndarray, spa.csc_matrix]] = None  # Equality matrix
    b: Optional[np.ndarray] = None                         # Equality vector
    lb: Optional[np.ndarray] = None                        # Lower bounds
    ub: Optional[np.ndarray] = None                        # Upper bounds
    
    def __init__(
        self,
        P: Union[np.ndarray, spa.csc_matrix],
        q: np.ndarray,
        G: Optional[Union[np.ndarray, spa.csc_matrix]] = None,
        h: Optional[np.ndarray] = None,
        A: Optional[Union[np.ndarray, spa.csc_matrix]] = None,
        b: Optional[np.ndarray] = None,
        lb: Optional[np.ndarray] = None,
        ub: Optional[np.ndarray] = None,
    ): ...
    
    def check_constraints(self) -> None:
        """
        Check constraint dimensions and consistency.
        
        Raises:
        - ProblemError: If constraints are inconsistent or malformed
        """
    
    def cond(self, active_set: ActiveSet) -> float:
        """
        Compute condition number of the problem matrix.
        
        Parameters:
        - active_set: Active constraint set for the condition number computation
        
        Returns:
        Condition number of the symmetric matrix representing problem data
        """
    
    def unpack(self) -> Tuple[
        Union[np.ndarray, spa.csc_matrix],  # P
        np.ndarray,                         # q
        Optional[Union[np.ndarray, spa.csc_matrix]],  # G
        Optional[np.ndarray],               # h
        Optional[Union[np.ndarray, spa.csc_matrix]],  # A
        Optional[np.ndarray],               # b
        Optional[np.ndarray],               # lb
        Optional[np.ndarray],               # ub
    ]:
        """
        Unpack problem matrices and vectors as a tuple.
        
        Returns:
        Tuple of (P, q, G, h, A, b, lb, ub)
        """
    
    @property
    def has_sparse(self) -> bool:
        """
        Check whether the problem has sparse matrices.
        
        Returns:
        True if at least one of P, G, or A matrices is sparse
        """
    
    @property
    def is_unconstrained(self) -> bool:
        """
        Check whether the problem has any constraint.
        
        Returns:
        True if the problem has no constraints (G, A, lb, ub all None)
        """
    
    def save(self, file: str) -> None:
        """
        Save problem to a compressed NumPy file.
        
        Parameters:
        - file: Either filename (string) or open file where data will be saved.
                If file is a string, the .npz extension will be appended if not present.
        """
    
    @staticmethod
    def load(file: str) -> 'Problem':
        """
        Load problem from a NumPy file.
        
        Parameters:
        - file: File-like object, string, or pathlib.Path to read from.
                File-like objects must support seek() and read() methods.
        
        Returns:
        Problem instance loaded from file
        """
    
    def get_cute_classification(self, interest: str) -> str:
        """
        Get the CUTE classification string of the problem.
        
        Parameters:
        - interest: Either 'A' (academic), 'M' (modeling), or 'R' (real application)
        
        Returns:
        CUTE classification string in format "QLR2-{interest}N-{nb_var}-{nb_cons}"
        
        Raises:
        - ParamError: If interest is not 'A', 'M', or 'R'
        """

Usage example:

import numpy as np
import scipy.sparse as spa
from qpsolvers import Problem

# Define QP matrices
P = np.array([[2.0, 1.0], [1.0, 2.0]])
q = np.array([1.0, -1.0])
G = np.array([[-1.0, 0.0], [0.0, -1.0]])
h = np.array([0.0, 0.0])

# Create problem instance
problem = Problem(P, q, G, h)

# Validate the problem
problem.check_constraints()

# Check condition number with active set
from qpsolvers import ActiveSet
active_set = ActiveSet()  # Empty active set
cond_num = problem.cond(active_set)
print(f"Condition number: {cond_num}")

# Unpack for solver
P_out, q_out, G_out, h_out, A_out, b_out, lb_out, ub_out = problem.unpack()

Solution Representation

The Solution class contains comprehensive solution information including primal/dual variables and optimality checks.

class Solution:
    """
    Solution returned by a QP solver for a given problem.
    
    Contains primal solution, dual multipliers, objective value, and
    optimality verification methods.
    """
    
    problem: Problem                    # The problem this solution corresponds to
    found: Optional[bool] = None        # True if solution found, False if infeasible
    obj: Optional[float] = None         # Primal objective value
    x: Optional[np.ndarray] = None      # Primal solution vector
    y: Optional[np.ndarray] = None      # Dual multipliers (equalities)
    z: Optional[np.ndarray] = None      # Dual multipliers (inequalities)
    z_box: Optional[np.ndarray] = None  # Dual multipliers (box constraints)
    extras: dict                        # Solver-specific information
    
    def __init__(
        self,
        problem: Problem,
        extras: dict = None,
        found: Optional[bool] = None,
        obj: Optional[float] = None,
        x: Optional[np.ndarray] = None,
        y: Optional[np.ndarray] = None,
        z: Optional[np.ndarray] = None,
        z_box: Optional[np.ndarray] = None,
    ): ...
    
    def is_optimal(self, eps_abs: float) -> bool:
        """
        Check whether the solution satisfies optimality conditions.
        
        Parameters:
        - eps_abs: Absolute tolerance for residuals and duality gap
        
        Returns:
        True if solution is optimal within tolerance
        """
    
    def primal_residual(self) -> float:
        """
        Compute primal feasibility residual.
        
        Returns:
        Maximum violation of primal constraints
        """
    
    def dual_residual(self) -> float:
        """
        Compute dual feasibility residual.
        
        Returns:
        Norm of dual residual (KKT stationarity condition)
        """
    
    def duality_gap(self) -> float:
        """
        Compute duality gap between primal and dual objectives.
        
        Returns:
        Absolute difference between primal and dual objectives
        """

Usage example:

from qpsolvers import solve_problem, Problem
import numpy as np

# Create and solve problem
problem = Problem(P, q, G, h, A, b)
solution = solve_problem(problem, solver="osqp")

# Check solution status
if solution.found:
    print(f"Solution found!")
    print(f"Objective value: {solution.obj}")
    print(f"Primal solution: {solution.x}")
    
    # Check optimality
    if solution.is_optimal(eps_abs=1e-6):
        print("Solution is optimal")
    
    # Examine dual multipliers
    if solution.z is not None:
        active_inequalities = np.where(solution.z > 1e-8)[0]
        print(f"Active inequality constraints: {active_inequalities}")
    
    # Check residuals
    print(f"Primal residual: {solution.primal_residual():.2e}")
    print(f"Dual residual: {solution.dual_residual():.2e}")
    print(f"Duality gap: {solution.duality_gap():.2e}")
    
    # Solver-specific information
    if solution.extras:
        print(f"Solver extras: {solution.extras}")

else:
    print("No solution found (infeasible/unbounded)")

Active Set Representation

The ActiveSet class tracks which inequality constraints are active at the solution.

class ActiveSet:
    """
    Indices of inequality constraints that are active (saturated) at the optimum.
    
    A constraint is active when it holds with equality at the solution.
    """
    
    G_indices: Sequence[int]   # Active linear inequality constraints
    lb_indices: Sequence[int]  # Active lower bound constraints
    ub_indices: Sequence[int]  # Active upper bound constraints
    
    def __init__(
        self,
        G_indices: Optional[Sequence[int]] = None,
        lb_indices: Optional[Sequence[int]] = None,
        ub_indices: Optional[Sequence[int]] = None,
    ): ...

Usage example:

from qpsolvers import ActiveSet, solve_problem, Problem
import numpy as np

# Solve a problem
problem = Problem(P, q, G, h, lb=np.zeros(2), ub=np.ones(2))
solution = solve_problem(problem, solver="osqp")

if solution.found:
    # Identify active constraints from dual multipliers
    active_G = []
    active_lb = []
    active_ub = []
    
    if solution.z is not None:
        active_G = np.where(solution.z > 1e-8)[0].tolist()
    
    if solution.z_box is not None:
        active_lb = np.where(solution.z_box < -1e-8)[0].tolist()
        active_ub = np.where(solution.z_box > 1e-8)[0].tolist()
    
    # Create active set
    active_set = ActiveSet(
        G_indices=active_G,
        lb_indices=active_lb,
        ub_indices=active_ub
    )
    
    print(f"Active linear inequalities: {active_set.G_indices}")
    print(f"Active lower bounds: {active_set.lb_indices}")
    print(f"Active upper bounds: {active_set.ub_indices}")

Advanced Usage

Problem Validation and Debugging

from qpsolvers import Problem, ProblemError
import numpy as np

def create_validated_problem(P, q, G=None, h=None, A=None, b=None, lb=None, ub=None):
    """Create a problem with comprehensive validation."""
    
    try:
        problem = Problem(P, q, G, h, A, b, lb, ub)
        problem.check_constraints()
        
        # Check positive semi-definiteness
        eigenvals = np.linalg.eigvals(P)
        if np.any(eigenvals < -1e-12):
            print("Warning: P matrix is not positive semi-definite")
        
        # Check condition number
        cond_num = problem.condition_number()
        if cond_num > 1e12:
            print(f"Warning: Problem is ill-conditioned (cond={cond_num:.2e})")
        
        return problem
        
    except ProblemError as e:
        print(f"Problem validation failed: {e}")
        return None

Solution Analysis

def analyze_solution(solution: Solution, tolerance: float = 1e-6):
    """Comprehensive solution analysis."""
    
    if not solution.found:
        print("No solution found")
        return
    
    print(f"Objective value: {solution.obj:.6f}")
    print(f"Solution vector: {solution.x}")
    
    # Optimality check
    is_opt = solution.is_optimal(tolerance)
    print(f"Is optimal (tol={tolerance}): {is_opt}")
    
    # Residual analysis
    primal_res = solution.primal_residual()
    dual_res = solution.dual_residual()
    gap = solution.duality_gap()
    
    print(f"Primal residual: {primal_res:.2e}")
    print(f"Dual residual: {dual_res:.2e}")
    print(f"Duality gap: {gap:.2e}")
    
    # Constraint analysis
    if solution.z is not None:
        active_ineq = np.where(solution.z > tolerance)[0]
        print(f"Active inequalities: {active_ineq.tolist()}")
    
    if solution.z_box is not None:
        active_lb = np.where(solution.z_box < -tolerance)[0]
        active_ub = np.where(solution.z_box > tolerance)[0]
        print(f"Active lower bounds: {active_lb.tolist()}")
        print(f"Active upper bounds: {active_ub.tolist()}")

Working with Large Problems

import scipy.sparse as spa
from qpsolvers import Problem, solve_problem

def create_sparse_problem(n: int):
    """Create a large sparse QP problem."""
    
    # Sparse cost matrix (tridiagonal)
    P = spa.diags([1, -2, 1], [-1, 0, 1], shape=(n, n), format="csc")
    P = P.T @ P  # Make positive semi-definite
    
    # Cost vector
    q = np.random.randn(n)
    
    # Sparse inequality constraints
    G = spa.random(n//2, n, density=0.1, format="csc")
    h = np.random.rand(n//2)
    
    return Problem(P, q, G, h)

# Create and solve large problem
large_problem = create_sparse_problem(1000)
solution = solve_problem(large_problem, solver="osqp")
analyze_solution(solution)

Install with Tessl CLI