ellalgo.oracles namespace

Submodules

ellalgo.oracles.ldlt_mgr module

LDLTMgr (LDLT Manager)

This code defines a class called LDLTMgr, which stands for LDLT Manager. The purpose of this class is to perform a special kind of matrix factorization called LDLT factorization on symmetric matrices. This factorization is useful in various mathematical and engineering applications, especially when dealing with linear systems and eigenvalue problems.

The main input for this class is a symmetric matrix, which can be provided either as a NumPy array or through a function that returns individual matrix elements. The primary output is a boolean value indicating whether the input matrix is positive definite (a special property of matrices). Additionally, the class can produce other outputs like a witness vector (if the matrix is not positive definite) and a square root of the matrix (if it is positive definite).

An important aspect of this code is its use of “lazy evaluation”. This means it doesn’t need the entire matrix upfront but can work with just a function that provides matrix elements as needed. This can be more efficient for large matrices or when the matrix is defined by a formula rather than stored values.

The class also includes methods to check if the matrix is positive definite (‘is_spd’), calculate a witness vector if it’s not (‘witness’), and compute a symmetric quadratic form (‘sym_quad’). These additional functionalities make the class versatile for various matrix-related tasks.

One of the key transformations happening in this code is the factorization itself. It’s breaking down the original matrix into simpler parts (L, D, and L transpose), which can be used to solve problems more easily. This factorization is done in a way that avoids using square roots, which can be computationally expensive.

Overall, the LDLTMgr class provides a set of tools for working with symmetric matrices, with a focus on determining their properties (like positive definiteness) and performing useful calculations efficiently. It’s designed to be flexible and can handle both standard matrices and those defined by functions, making it useful in a variety of mathematical and engineering contexts.

class ellalgo.oracles.ldlt_mgr.LDLTMgr(ndim: int)[source]

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[source]

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[source]

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[source]

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[source]

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[source]

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[source]

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[source]

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

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

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[source]

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

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

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[source]

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

Raises:

None

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

LMIOracle

This code defines a class called LMIOracle, which is designed to solve a specific type of mathematical problem known as a Linear Matrix Inequality (LMI) constraint. The purpose of this code is to determine if there exists a solution that satisfies the given constraint, and if not, to provide information about why it’s not feasible.

The LMIOracle class takes two inputs when it’s created: mat_f and mat_b. These are matrices that define the LMI constraint. mat_f is a list of numpy arrays, and mat_b is a single numpy array. These matrices represent the mathematical relationship that needs to be satisfied.

The main function in this class is called assess_feas, which takes a numpy array xc as input. This function tries to determine if xc is a feasible solution to the LMI constraint. If it is feasible, the function returns None. If it’s not feasible, it returns what’s called a “cut” - a tuple containing information about why the solution isn’t feasible.

To achieve its purpose, the code uses a technique called LDLT factorization. This is a way of breaking down a matrix into simpler parts, which can help determine if the matrix satisfies certain properties. The LDLTMgr class (which is used but not defined in this code snippet) handles this factorization.

The assess_feas function works by first creating a new matrix using the input xc and the original matrices mat_f and mat_b. It then tries to perform the LDLT factorization on this new matrix. If the factorization is successful, it means xc is a feasible solution, and the function returns None. If the factorization fails, it means xc is not feasible, and the function calculates and returns the “cut” information.

An important part of the logic is the get_elem function inside assess_feas. This function calculates each element of the new matrix based on the original matrices and the input xc. This is where the main mathematical operation of the LMI constraint is performed.

In summary, this code provides a way to check if a given solution satisfies a complex mathematical constraint, and if not, it provides information about why the solution doesn’t work. This could be useful in optimization problems where you’re trying to find a solution that satisfies certain mathematical conditions.

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

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[source]

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

Lowpass Oracle

This code implements a LowpassOracle, which is used to design a low-pass filter for signal processing. A low-pass filter allows low-frequency signals to pass through while attenuating high-frequency signals. The main purpose of this code is to help optimize the design of such a filter by providing a way to assess whether a given set of filter coefficients meets certain specifications.

The code defines a class called LowpassOracle that takes several inputs when initialized:

  1. ndim: The number of filter coefficients

  2. wpass: The end of the passband (frequencies that should pass through)

  3. wstop: The end of the stopband (frequencies that should be attenuated)

  4. lp_sq: The lower bound for the squared magnitude response in the passband

  5. up_sq: The upper bound for the squared magnitude response in the passband

  6. sp_sq: The upper bound for the squared magnitude response in the stopband

The main outputs of this code are produced by two methods: assess_feas and assess_optim. These methods take a set of filter coefficients as input and determine whether they meet the specified requirements or how close they are to meeting them.

The LowpassOracle achieves its purpose through a series of checks on the frequency response of the filter. It uses a pre-computed spectrum matrix to efficiently calculate the frequency response at different points. The code then checks if the response falls within the specified bounds for the passband and stopband.

The important logic flow in this code involves iterating through different frequency points and checking the filter’s response at each point. If any violations of the specifications are found, the code returns information about the violation, which can be used to adjust the filter coefficients.

A key data transformation happening in this code is the conversion from filter coefficients to frequency response. This is done using the pre-computed spectrum matrix, which allows for efficient calculation of the response at many frequency points.

The code also includes a helper function called create_lowpass_case, which sets up a specific instance of the LowpassOracle with predefined parameters. This function can be used to quickly create a standard test case for filter design.

Overall, this code provides a tool for iteratively designing and optimizing low-pass filters by giving feedback on how well a set of coefficients meets the desired specifications. It’s part of a larger optimization process where the coefficients would be adjusted based on the feedback from this oracle until a satisfactory filter design is achieved.

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

Bases: OracleOptim

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

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][source]

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[source]

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 Oracle

This code defines several classes that implement oracles for profit maximization problems. An oracle, in this context, is a function that helps solve optimization problems by providing information about the feasibility and optimality of potential solutions.

The main class, ProfitOracle, is designed to solve a specific type of profit maximization problem. It takes as input parameters related to production (like price, scale, and limits) and output elasticities. The goal is to find the optimal input quantities that maximize profit, given certain constraints.

The ProfitOracle class has methods to assess the feasibility of a solution (assess_feas) and to find the optimal solution (assess_optim). These methods take as input a vector y (representing input quantities in log scale) and a gamma value (representing the current best solution). They output “cuts”, which are linear constraints that help narrow down the search for the optimal solution.

The code also includes two variations of the profit oracle:

  1. ProfitRbOracle: This is a robust version of the profit oracle that can handle some uncertainty in the input parameters.

  2. ProfitQOracle: This version deals with discrete (integer) input quantities, as opposed to continuous ones.

The main logic flow in these classes involves calculating various economic functions (like Cobb-Douglas production functions) and their gradients. The code uses these calculations to determine if a given solution is feasible and to guide the search towards the optimal solution.

The output of these oracles is typically a “cut” (a linear constraint) and possibly an updated best solution (gamma). These outputs are used by an external optimization algorithm (not shown in this code) to iteratively improve the solution until the optimal one is found.

For beginners, it’s important to understand that this code is implementing mathematical optimization techniques. While the details might be complex, the basic idea is to efficiently search for the best solution to a profit maximization problem, given certain constraints and economic relationships.

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

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[source]

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][source]

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[source]

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[source]

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[source]

Gradient for y₁ ≤ k constraint.

Returns:

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

grad2(gamma: float) → ndarray[source]

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)[source]

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][source]

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])[source]

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][source]

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 Code

This code implements spectral factorization, which is a mathematical technique used in signal processing. The main purpose of this code is to compute a minimum-phase impulse response that satisfies a given auto-correlation. In simpler terms, it’s trying to find a special sequence of numbers (the impulse response) that, when processed in a certain way, matches a given pattern of relationships between data points (the auto-correlation).

The code contains two main functions: spectral_fact and inverse_spectral_fact.

The spectral_fact function takes one input: r, which is the top-half of the auto-correlation coefficients. This input should be a list or array of numbers. The function outputs h, which is the impulse response that gives the desired auto-correlation.

To achieve its purpose, the spectral_fact function follows these steps:

  1. It determines the length of the input and creates an oversampled version of it.

  2. It computes a logarithmic representation of the input in the frequency domain.

  3. It applies a mathematical operation called the Hilbert transform.

  4. It combines the results of steps 2 and 3 to create a complex representation.

  5. Finally, it converts this representation back to the time domain to get the impulse response.

The function uses several mathematical operations like Fourier transforms, logarithms, and complex number manipulations to achieve this. These operations help transform the data between different representations (time domain and frequency domain) and extract the necessary information to compute the impulse response.

The inverse_spectral_fact function does the opposite of spectral_fact. It takes the impulse response h as input and attempts to reconstruct the original auto-correlation coefficients. This function is simpler and uses a mathematical operation called convolution to compute its result.

Overall, this code provides tools for working with signal processing problems, particularly those involving auto-correlations and impulse responses. It’s useful in fields like audio processing, communications, and data analysis where understanding the relationships between data points over time is important.

Spectral Factorization Process Diagram:

```svgbob
       Auto-correlation
             |
             v
    +-------------------+
    |  Oversampling     |
    +-------------------+
             |
             v
    +-------------------+
    | Log computation   |
    |  alpha(w) = 1/2*  |
    |  ln(R(w))         |
    +-------------------+
             |
             v
    +-------------------+
    | Hilbert Transform |
    |  phi(w) = H[alpha]|
    +-------------------+
             |
             v
    +-------------------+
    | Complex Rep.      |
    |  H(exp(jTw)) =    |
    |  alpha(w) + j*phi |
    +-------------------+
             |
             v
    +-------------------+
    | Inverse FFT       |
    |  (Time Domain)    |
    +-------------------+
             |
             v
     Impulse Response
```
ellalgo.oracles.spectral_fact.inverse_spectral_fact(h: ndarray) → ndarray[source]

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[source]

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