ellalgo.oracles namespace

Submodules

ellalgo.oracles.ldlt_mgr module

LDL^T factorization manager for symmetric matrices.

The LDLTMgr class implements a square-root-free Cholesky decomposition (A = LDL^T) for symmetric matrices. Key features:

• Lazy evaluation: Accepts a callback for element access, avoiding full matrix construction. Useful for large or implicitly-defined matrices.

• Positive definiteness check via diagonal entry monitoring.

• Witness vector computation for non-SPD matrices, certifying infeasibility.

• Quadratic form computation using the witness vector.

• Square root (Cholesky factor) extraction for SPD matrices.

This is primarily used by LMI (Linear Matrix Inequality) oracles to check positive semidefiniteness during cutting-plane iterations.

class ellalgo.oracles.ldlt_mgr.LDLTMgr(ndim: int)

Bases: object

The LDLTMgr class implements a square-root-free version of the Cholesky decomposition, known as LDLT factorization. This method decomposes a symmetric matrix A into A = LDL^T, where L is a lower triangular matrix with ones on the diagonal, D is a diagonal matrix, and L^T is the transpose of L.

This factorization is particularly useful for Linear Matrix Inequality (LMI) oracles in optimization problems. Its main advantages include:

• Numerical Stability: By avoiding the computation of square roots, the LDLT factorization can be more numerically stable than the standard Cholesky decomposition.

• Efficiency: The square-root-free nature of the algorithm can also lead to performance improvements.

• Lazy Evaluation: The implementation supports lazy evaluation, allowing it to work with matrices that are not explicitly stored in memory.

The class provides the following capabilities:

• Check if a matrix is symmetric positive-definite (SPD).

• Find a “witness” vector that certifies that a matrix is not SPD.

• Compute the Cholesky factorization (R matrix such that A = R^T R) if the matrix is SPD.

factor(get_elem: Callable[[int, int], float]) → bool

Performs LDLT factorization using lazy element access.

The factorization proceeds row by row, computing diagonal entries and off-diagonal multipliers. If any diagonal entry becomes non-positive, factorization stops early and records the failure position.

Parameters:

get_elem – Function that returns matrix element at (i,j) position

Returns:

True if matrix is positive definite (all diagonal entries positive)

Return type:

bool

The factorization stores results in _storage: - Diagonal entries (D matrix) are stored in _storage[i,i] - Off-diagonal entries (L matrix) are stored in _storage[i,j] for j < i

Examples

>>> mat = np.array([[1.0, 0.5, 0.5], [0.5, 1.25, 0.75], [0.5, 0.75, 1.5]])
>>> ldl = LDLTMgr(3)
>>> ldl.factor(lambda i, j: mat[i, j])
True

factor_with_allow_semidefinite(get_elem: Callable[[int, int], float]) → bool

Performs LDLT factorization allowing for positive semi-definite matrices.

Similar to factor() but handles zero diagonal entries (indicating semi-definiteness) by restarting factorization from the next row.

Parameters:

get_elem – Function that returns matrix element at (i,j) position

Returns:

True if matrix is positive semi-definite (no negative diagonal entries)

Return type:

bool

This version is more tolerant of zero diagonal entries than factor(), which requires strictly positive entries for positive definiteness.

Examples

>>> mat = np.array([[1.0, 0.5, 0.5], [0.5, 1.25, 0.75], [0.5, 0.75, 1.5]])
>>> ldl = LDLTMgr(3)
>>> ldl.factor_with_allow_semidefinite(lambda i, j: mat[i, j])
True

factorize(mat: ndarray) → bool

Performs LDLT factorization on a NumPy array.

This method is a convenience wrapper around the factor method. It allows you to perform the factorization directly on a NumPy array without needing to provide a custom element access function.

Parameters:

mat (np.ndarray) – The symmetric matrix to be factorized.

Returns:

True if the matrix is positive-definite, False otherwise.

Return type:

bool

is_spd() → bool

Checks if the matrix is symmetric positive definite (SPD).

Returns:

True if the matrix is SPD (pos[1] == 0), False otherwise

Return type:

bool

The check is based on whether any diagonal entry was non-positive during factorization, which would have set pos[1] to a non-zero value.

Examples

>>> mat = np.array([[1.0, 0.5, 0.5], [0.5, 1.25, 0.75], [0.5, 0.75, 1.5]])
>>> ldl = LDLTMgr(3)
>>> ldl.factorize(mat)
True
>>> ldl.is_spd()
True

pos: Tuple[int, int]

sqrt() → ndarray

Computes the upper triangular square root matrix R where A = R^T R.

Returns:

Upper triangular matrix R

Return type:

np.ndarray

Raises:

AssertionError – If matrix is not positive definite

This is essentially the Cholesky decomposition, computed from the LDLT factors without directly computing square roots until the final step.

Examples

>>> mat = np.array([[1.0, 0.5, 0.5], [0.5, 1.25, 0.75], [0.5, 0.75, 1.5]])
>>> ldl = LDLTMgr(3)
>>> ldl.factorize(mat)
True
>>> ldl.sqrt()
array([[1. , 0.5, 0.5],
       [0. , 1. , 0.5],
       [0. , 0. , 1. ]])

sym_quad(mat: ndarray) → float

Computes the quadratic form v^T M v using the witness vector.

Parameters:

mat – The matrix M to compute the quadratic form with

Returns:

The value of v^T M v where v is the witness vector

Return type:

float

Note: witness() must be called first to set up the witness vector. The computation uses only the submatrix where positive definiteness failed.

Examples

>>> mat = np.array([[1.0, 2.0, 3.0], [2.0, 3.5, 5.0], [3.0, 5.0, 6.0]])
>>> ldl = LDLTMgr(3)
>>> ldl.factorize(mat)
False
>>> ldl.pos
(0, 2)
>>> ldl.witness() # call this before sym_quad()
np.float64(0.5)
>>> ldl.wit
array([-2.,  1.,  0.])
>>> mat_b = np.array([[1.0, 0.5, 0.5], [0.5, 1.25, 0.75], [0.5, 0.75, 1.5]])
>>> ldl.sym_quad(mat_b)
np.float64(3.25)

wit: ndarray

witness() → float

Computes a witness vector proving the matrix is not positive definite.

Returns:

The negative eigenvalue (ep) showing v^T A v = -ep < 0

Return type:

float

Raises:

AssertionError – If called on a positive definite matrix

The witness vector is stored in self.wit and can be accessed after calling this method. The vector satisfies v^T A v < 0 for the failed submatrix.

Examples

>>> mat = np.array([[1.0, 2.0, 3.0], [2.0, 3.5, 5.0], [3.0, 5.0, 6.0]])
>>> ldl = LDLTMgr(3)
>>> ldl.factorize(mat)
False
>>> ldl.witness()
np.float64(0.5)

ellalgo.oracles.lmi0_oracle module

LMI oracle with zero constant term.

This module provides a feasibility oracle for Linear Matrix Inequality (LMI) constraints of the form F(x) = Σ x_k F_k ⪰ 0, where the constant matrix B is zero. It uses LDL^T factorization to check positive semidefiniteness and returns separating hyperplanes for infeasible points.

class ellalgo.oracles.lmi0_oracle.LMI0Oracle(mat_f: List[ndarray])

Bases: object

Oracle for the Linear Matrix Inequality (LMI) constraint: F(x) ⪰ 0.

This class is a specialized oracle for solving the LMI feasibility problem where the constant matrix B is zero. The constraint is of the form:

F(x) = F₁x₁ + F₂x₂ + … + Fₙxₙ ⪰ 0

where Fᵢ are symmetric matrices and x is the vector of decision variables.

The assess_feas method checks if a given solution x satisfies the LMI constraint. If it does, the method returns None. If not, it returns a separating hyperplane (a “cut”) that separates the infeasible point from the feasible set.

assess_feas(x: ndarray) → Tuple[ndarray, float] | None

Assess the feasibility of a candidate solution x.

This method checks if the given solution x satisfies the LMI constraint F(x) ⪰ 0. It does this by constructing the matrix F(x) and performing an LDLT factorization to determine if it is positive semidefinite.

Parameters:

x (np.ndarray) – The candidate solution vector.

Returns:

None if x is feasible (i.e., the LMI constraint is satisfied). Otherwise, a tuple (g, ep) representing a separating hyperplane, where g is the subgradient and ep is the measure of violation.

Return type:

Optional[Cut]

ellalgo.oracles.lmi_old_oracle module

Legacy LMI oracle with explicit matrix construction.

This module provides a feasibility oracle for Linear Matrix Inequality (LMI) constraints. Unlike LMIOracle (which uses lazy element-wise evaluation), this implementation constructs the full LMI matrix explicitly before factorization.

class ellalgo.oracles.lmi_old_oracle.LMIOldOracle(mat_f: List[ndarray], mat_b: ndarray)

Bases: OracleFeas

Oracle for Linear Matrix Inequality constraint.

This oracle solves the following feasibility problem:

find x s.t. (B − F * x) ⪰ 0

This is a legacy implementation that constructs the full LMI matrix explicitly. For better performance with large matrices, use LMIOracle which uses lazy evaluation.

Examples

>>> import numpy as np
>>> from ellalgo.oracles.lmi_old_oracle import LMIOldOracle
>>> F1 = np.array([[1.0, 0.0], [0.0, 1.0]])
>>> F2 = np.array([[0.0, 1.0], [1.0, 0.0]])
>>> B = np.array([[2.0, 0.0], [0.0, 2.0]])
>>> oracle = LMIOldOracle([F1, F2], B)
>>> result = oracle.assess_feas(np.array([0.0, 0.0]))
>>> result is None or isinstance(result, tuple)
True

assess_feas(xc: ndarray) → Tuple[ndarray, float] | None

Assess the feasibility of a candidate solution.

This method checks if the given solution satisfies the LMI constraint (B − F₁x₁ − F₂x₂ − … − Fₙxₙ) ⪰ 0 by constructing the full matrix and performing LDLT factorization.

Parameters:

xc – The candidate solution vector x.

Returns:

None if feasible, otherwise a tuple (g, ep) containing the subgradient g and the negative eigenvalue measure ep.

Examples

>>> import numpy as np
>>> from ellalgo.oracles.lmi_old_oracle import LMIOldOracle
>>> F1 = np.array([[1.0, 0.0], [0.0, 1.0]])
>>> F2 = np.array([[0.0, 1.0], [1.0, 0.0]])
>>> B = np.array([[2.0, 0.0], [0.0, 2.0]])
>>> oracle = LMIOldOracle([F1, F2], B)
>>> oracle.assess_feas(np.array([0.0, 0.0])) is None
True

ellalgo.oracles.lmi_oracle module

LMI (Linear Matrix Inequality) feasibility oracle.

The LMIOracle class implements a feasibility oracle for LMI constraints of the form:

B - (F₁x₁ + F₂x₂ + … + Fₙxₙ) ⪰ 0

where B and Fᵢ are symmetric matrices. It uses lazy element-wise matrix construction and LDL^T factorization to check positive semidefiniteness, avoiding construction of the full matrix whenever possible.

When a point is infeasible, the oracle returns a separating hyperplane (cut) derived from the LDL^T witness vector, enabling the cutting-plane algorithm to narrow the search space.

class ellalgo.oracles.lmi_oracle.LMIOracle(mat_f: List[ndarray], mat_b: ndarray)

Bases: OracleFeas

Oracle for Linear Matrix Inequality (LMI) constraints.

This class implements the OracleFeas interface for solving semidefinite feasibility problems involving Linear Matrix Inequalities (LMIs). An LMI constraint is of the form:

B - (F₁x₁ + F₂x₂ + … + Fₙxₙ) ⪰ 0

where B and Fᵢ are symmetric matrices, and x is the vector of decision variables. The notation ⪰ 0 means that the resulting matrix is required to be positive semidefinite.

The assess_feas method checks if a given solution x satisfies the LMI constraint. If it does, the method returns None. If not, it returns a separating hyperplane (a “cut”) that separates the infeasible point from the feasible set.

assess_feas(xc: ndarray) → Tuple[ndarray, float] | None

Assess the feasibility of a candidate solution xc.

This method checks if the given solution xc satisfies the LMI constraint. It does this by constructing the matrix M(xc) and performing an LDLT factorization to determine if it is positive semidefinite.

Parameters:

xc (np.ndarray) – The candidate solution vector.

Returns:

None if xc is feasible (i.e., the LMI constraint is satisfied). Otherwise, a tuple (g, ep) representing a separating hyperplane, where g is the subgradient and ep is the measure of violation.

Return type:

Optional[Cut]

ellalgo.oracles.lowpass_oracle module

FIR low-pass filter design oracle via spectral factorization.

This module implements the approach from:

S.-P. Wu, S. Boyd, and L. Vandenberghe, “FIR Filter Design via Spectral Factorization and Convex Optimization”

The LowpassOracle class formulates the FIR filter design problem as a convex optimization over the auto-correlation coefficients. It checks passband ripple and stopband attenuation constraints using a pre-computed spectrum matrix for efficient frequency response evaluation at discretized frequency points.

Key methods:

• assess_feas: Check whether filter coefficients meet passband/stopband specs.

• assess_optim: Assess optimality, returning the maximum stopband response.

Also provides create_lowpass_case() for a standard test case with typical parameters (passband 0-0.12π, stopband 0.20-π, ±0.025dB ripple).

class ellalgo.oracles.lowpass_oracle.LowpassOracle(ndim: int, wpass: float, wstop: float, lp_sq: float, up_sq: float, sp_sq: float)

Bases: OracleOptim

assess_feas(x: ndarray) → Tuple[ndarray, float | Tuple[float, float | None] | List[float]] | None

Assess whether the given filter coefficients meet the design specifications.

This method checks the frequency response at various points in three bands: 1. Passband (0 to nwpass): Checks if response is within [lp_sq, up_sq] 2. Stopband (nwstop to end): Checks if response is below sp_sq and non-negative 3. Transition band (nwpass to nwstop): Checks if response is non-negative

Uses a round-robin approach to check different frequency points on each call to distribute the computational load across multiple iterations.

Parameters:

x (Arr) – The filter coefficients (autocorrelation coefficients)

Returns:

• None if all specifications are met

• A tuple containing:

  • The gradient of the violating constraint

  • The violation amount (or tuple of lower/upper violations)

Return type:

Optional[ParallelCut]

assess_optim(xc: ndarray, gamma: float) → Tuple[Tuple[ndarray, float | Tuple[float, float | None] | List[float]], float | None]

Assess the optimality of the current filter coefficients for the stopband.

First checks feasibility using assess_feas. If feasible, returns information about the maximum response in the stopband which can be used to further optimize the filter design.

Parameters:

• xc (Arr) – The filter coefficients (autocorrelation coefficients)

• gamma (float) – The current best stopband attenuation value to beat

Returns:

A tuple containing:

• A tuple of (gradient, (lower, upper)) for the maximum stopband response

• The maximum stopband response value (or None if not feasible)

Return type:

tuple

idx1: int = 0

ellalgo.oracles.lowpass_oracle.create_lowpass_case(ndim: int = 48) → LowpassOracle

Creates a standard low-pass filter design case with typical parameters.

Sets up a LowpassOracle instance with commonly used specifications: - Passband edge at 0.12π - Stopband edge at 0.20π - Passband ripple of ±0.025 dB - Stopband attenuation of 0.125

Parameters:

ndim (int, optional) – Number of filter coefficients. Defaults to 48.

Returns:

An initialized LowpassOracle instance with standard parameters

Return type:

LowpassOracle

ellalgo.oracles.profit_oracle module

Profit maximization oracles with Cobb-Douglas production functions.

This module implements oracles for profit maximization where the production function is of the Cobb-Douglas type: q = A ⋅ y₁^α ⋅ y₂^β.

Three oracle variants:

• ProfitOracle: Standard continuous optimization over log-space variables.

• ProfitRbOracle: Robust variant handling parameter uncertainties in elasticities, prices, and production limits.

• ProfitQOracle: Discrete (integer-quantized) variant that rounds continuous solutions and manages infeasibility via retry mechanism.

All oracles return cutting planes (subgradient + violation) that the cutting-plane algorithm uses to iteratively narrow the search space.

class ellalgo.oracles.profit_oracle.ProfitOracle(params: Tuple[float, float, float], elasticities: ndarray, price_out: ndarray)

Bases: OracleOptim

Oracle for a profit maximization problem with a Cobb-Douglas production function.

This class implements the OracleOptim interface for a specific profit maximization problem. The production function is of the Cobb-Douglas type, which is widely used in economics to represent the relationship between production inputs and the amount of output.

The optimization problem is to maximize the profit, which is the difference between the revenue from selling the product and the cost of the inputs. The problem is subject to a constraint on one of the inputs.

The assess_optim method is the core of the oracle. It takes a candidate solution (a vector of input quantities) and the current best profit, and it returns a cutting plane that helps to narrow down the search for the optimal solution.

assess_feas(xc: ndarray, gamma: float) → Tuple[ndarray, float] | None

Feasibility assessment using round-robin constraint checking.

Implements: - Alternates between checking y₁ constraint (fn1) and optimality (fn2) - Returns first violated constraint found

Parameters:

• xc – Current solution point in log-space

• gamma – Current best profit estimate

Returns:

Cut (gradient, violation) if constraint violated None if all constraints satisfied

assess_optim(xc: ndarray, gamma: float) → Tuple[Tuple[ndarray, float], float | None]

Assess the optimality of a candidate solution xc.

This method is the core of the ProfitOracle. It takes a candidate solution xc and the current best profit gamma, and it returns a cutting plane that helps to narrow down the search for the optimal solution.

The method first checks if the solution is feasible. If not, it returns a feasibility cut. If the solution is feasible, it calculates the profit at xc and generates an optimality cut.

Parameters:

• xc (Arr) – The candidate solution vector (in log-space).

• gamma (float) – The current best profit.

Returns:

A tuple containing the cutting plane (g, beta) and the updated profit gamma_new. If the solution is infeasible, gamma_new is None.

Return type:

Tuple[Cut, Optional[float]]

elasticities: ndarray

fn1(x: ndarray, _: float) → float

Constraint function for y₁ ≤ k (in log-space).

Parameters:

x – Log-scale input vector [log(y₁), log(y₂)]

Returns:

x[0] - log(k) Positive values indicate constraint violation

Return type:

Constraint violation measure

fn2(x: ndarray, gamma: float) → float

Optimality condition function for profit maximization.

Computes: - Cobb-Douglas value in log-space: log(pA) + αlog(y₁) + βlog(y₂) - Variable costs: v₁y₁ + v₂y₂ - Optimality gap: log(γ + vy) - log_Cobb

Parameters:

• x – Log-scale input vector

• gamma – Current best profit estimate

Updates intermediate values used in gradient calculations

grad1(_: float) → ndarray

Gradient for y₁ ≤ k constraint.

Returns:

Gradient vector [1, 0] since ∂(x₀ - log_k)/∂x = (1, 0)

grad2(gamma: float) → ndarray

Gradient of optimality condition function.

Computes:

∇f = [v₁y₁/(γ+vy) - α, v₂y₂/(γ+vy) - β]

Parameters:

gamma – Current profit estimate used in denominator

Uses precomputed q (v₁y₁, v₂y₂) from last fn2 call

idx: int = -1

log_Cobb: float

log_k: float

log_pA: float

price_out: ndarray

q: ndarray

vy: float

class ellalgo.oracles.profit_oracle.ProfitQOracle(params: Tuple[float, float, float], elasticities: ndarray, price_out: ndarray)

Bases: OracleOptimQ

Discrete profit oracle for integer input quantities.

Solves mixed-integer version of the profit maximization problem:

max p(A y₁^α y₂^β) − v₁y₁ − v₂y₂ s.t. x₁ ≤ k, x ∈ ℕ²

Uses continuous relaxation followed by rounding to nearest integer, with mechanisms to handle infeasible integer solutions.

assess_optim_q(xc: ndarray, gamma: float, retry: bool) → Tuple[Tuple[ndarray, float], ndarray, float | None, bool]

Discrete optimization assessment with rounding mechanism.

Workflow: 1. First try continuous solution (retry=False) 2. If infeasible, return feasibility cut 3. If feasible, round to nearest integer and check optimality 4. On retry (retry=True), check rounded solution optimality

Returns:

• Cut information

• Evaluation point (continuous or rounded)

• Updated gamma value

• Retry flag for integer feasibility check

Return type:

Tuple containing

xd: ndarray

class ellalgo.oracles.profit_oracle.ProfitRbOracle(params: Tuple[float, float, float], elasticities: ndarray, price_out: ndarray, vparams: Tuple[float, float, float, float, float])

Bases: OracleOptim

Robust profit oracle handling parameter uncertainty.

Implements robust optimization version from [Aliabadi and Salahi, 2013] considering uncertainties in: - Elasticity parameters (α, β) - Price parameters (p, v) - Production limit (k)

Uses interval-based uncertainty sets for robust constraint satisfaction.

assess_optim(xc: ndarray, gamma: float) → Tuple[Tuple[ndarray, float], float | None]

Robust optimization assessment accounting for parameter uncertainties.

Adjusts elasticities based on direction of uncertainty impact: - Decreases effective α, β when y > 0 (conservative adjustment) - Increases effective α, β when y ≤ 0 (worst-case scenario)

ellalgo.oracles.spectral_fact module

Spectral factorization for minimum-phase impulse response computation.

Implements the Kolmogorov 1939 spectral factorization approach as described in A. Papoulis, “Signal Analysis” (pp. 232-233). This is used by the LowpassOracle to convert between auto-correlation coefficients and the minimum-phase impulse response of an FIR filter.

Functions:

• spectral_fact(r): Compute minimum-phase impulse response from auto-correlation.

• inverse_spectral_fact(h): Reconstruct auto-correlation from impulse response.

The spectral factorization pipeline:

auto-correlation → oversampling → log(|R(w)|) → Hilbert transform → complex log-spectrum → IFFT → impulse response

ellalgo.oracles.spectral_fact.inverse_spectral_fact(h: ndarray) → ndarray

Computes the auto-correlation sequence from the given impulse response.

Parameters:

h (numpy.ndarray) – The impulse response sequence.

Returns:

The auto-correlation sequence, where the length is the same as the input impulse response.

Return type:

numpy.ndarray

Examples

>>> h = np.array([1.0, 0.5, 0.2])
>>> r = inverse_spectral_fact(h)
>>> isinstance(r, np.ndarray)
True
>>> r.shape == (len(h),)
True

ellalgo.oracles.spectral_fact.spectral_fact(r: ndarray) → ndarray

Computes the minimum-phase impulse response satisfying a given auto-correlation.

This function implements the Kolmogorov 1939 approach to spectral factorization, as described in pp. 232-233 of “Signal Analysis” by A. Papoulis.

Parameters:

r (numpy.ndarray) – The top-half of the auto-correlation coefficients, starting from the 0th element to the end of the auto-correlation. This should be passed in as a column vector.

Returns:

The impulse response that gives the desired auto-correlation.

Return type:

numpy.ndarray

Raises:

• ValueError – If the input array is empty or contains invalid values.

• RuntimeError – If numerical errors occur during spectral factorization (e.g., log of negative numbers, FFT errors).

Examples

>>> r = np.array([1.0, 0.5, 0.2])
>>> h = spectral_fact(r.reshape(-1, 1))
>>> isinstance(h, np.ndarray)
True
>>> h.shape == (r.shape[0], r.shape[0])
True