layout: true class: typo, typo-selection --- count: false class: nord-dark, middle, center # 👁️ The Ellipsoid Method and Its Amazing Oracles @luk036 👨‍💻 · 2026 📅 --- class: middle .right[ > When you have eliminated the impossible, whatever remains, however > improbable, must be the truth. 🕵️‍♂️ Sir Arthur Conan Doyle, stated by Sherlock Holmes 🔍 ] --- ## Abstract 📝 The ellipsoid method is a powerful optimization technique 💪 that offers distinct advantages over interior-point methods, as it does not require the evaluation of all constraint functions. This makes it the optimal choice 🏆 for convex problems with numerous or even infinite constraints ∞. The method employs an ellipsoid as a search space 🪐 and replies on a separation oracle 🔮 to provide cutting planes ✂️ for updating it. It is worth noting that the significance of the separation oracle is often overlooked. This article evaluates the utility of the ellipsoid method in three distinct applications: robust convex optimization 🛡️, semidefinite programming 📊, and parametric network optimization 🌐. The effectiveness of separation oracles is assessed for each application. Furthermore, this article addresses the implementation issues associated with the ellipsoid method, including the utilization of parallel cuts ⚡ for updating the ellipsoid. In certain cases, the use of parallel cuts has been observed to reduce computation time ⏱️, as evidenced in the context of FIR filter design 🎛️. The article also considers discrete optimization 🧩, demonstrating how the ellipsoid method can be applied to problems involving quantized discrete design variables. The additional effort in oracle implementation is limited to locating the nearest discrete solutions. --- class: nord-light, middle, center ## 🎬 Introduction --- ### Historical Significance 🏛️ - 💡 Early Development Preliminary versions introduced by Naum Z. Shor, with significant contributions from Arkadi Nemirovski and David B. Yudin in 1972 through their approximation algorithm for real convex minimization. - 🏅 Khachiyan's Breakthrough Leonid Khachiyan's work in the late 1970s proved the polynomial-time solvability of linear programming using the ellipsoid method, a landmark result in computational complexity theory. - 💥 Theoretical Impact Demonstrated that linear programming belongs to the class P of problems solvable in polynomial time ⏳, with complexity depending on variables and coefficient size, but not on the number of constraints. --- ### Performance Considerations ⚖️ - Theoretical Strengths 💪 The ellipsoid method offers polynomial runtime complexity for convex optimization problems under certain assumptions. It's considered "R-polynomial" with operations bounded by a polynomial in problem size multiplied by the logarithm of volume factors. The number of steps required is generally $O(n^2 \log(1/ε))$, where $n$ is the problem dimension and $ε$ is the desired accuracy. - Practical Limitations 🚧 Despite theoretical power, the method often faces criticism for practical performance compared to interior-point methods. The iteration count grows quadratically with dimension, becoming a bottleneck for high-dimensional problems. The method can also be susceptible to numerical instability in practice, limiting its application for general large-scale convex optimization. --- ### But...🤔 - The ellipsoid method works very differently compared with the interior point methods. - It only requires a *separation oracle* 🔮 that provides a cutting plane ✂️. - Although the ellipsoid method cannot take advantage of the sparsity of the problem, the separation oracle is capable of take advantage of certain structural types. --- ### Consider the ellipsoid method when...🤔 - The number of design variables is moderate, e.g. ECO flow, analog circuit sizing, parametric problems - The number of constraints is large, or even infinite ∞ - Oracle can be implemented effectively. ⚡ --- ### Power of Separation Oracles ⚡ - 🔮 Oracle Definition A separation oracle is an algorithmic procedure that, given a point x, either confirms membership in a convex set Q or returns a hyperplane that separates x from Q. This approach doesn't require explicit enumeration of all constraints. - Handling Infinite Constraints ∞ The true strength of the ellipsoid method lies in its ability to work with convex sets defined by enormous or infinite constraint collections, as long as an efficient separation oracle exists. - Exploiting Problem Structure 🏗️ Effective oracles must be computationally efficient and exploit the specific structure of the optimization problem, as they're called at each iteration of the algorithm. --- class: nord-light, middle, center ### 🥥 Cutting-plane Method Revisited --- ### Convex Set 🥚 .pull-left70[ - Let $\mathcal{\color{steelblue}K} \subseteq \mathbb{R}^n$ be a convex set. - Consider the feasibility problem: - Find a point ${\color{darkblack}x}^* \in \mathbb{R}^n$ in $\mathcal{\color{steelblue}K}$, - or determine that $\mathcal{\color{steelblue}K}$ is empty (i.e., there is no feasible solution) ] .pull-right30[  ] --- ### 🔮 Separation Oracle .pull-left[ - When a separation oracle $\Omega$ is *queried* at $x_0$, it either - asserts that $x_0 \in \mathcal{\color{steelblue}K}$, or - returns a hyperplane between $x_0$ and $\mathcal{\color{steelblue}K}$: $$ {\color{darkcyan} g^\mathsf{T} } ({\color{darkblack}x} - x_0) + {\color{green} \beta} \le 0, \\ {\color{green} \beta} \geq 0, {\color{darkcyan} g} \neq 0, \; \forall {\color{darkblack}x} \in \mathcal{\color{steelblue}K} $$ ] .pull-right[  ] --- ### 🔮 Separation Oracle (cont'd) - $({\color{darkcyan} g}, {\color{green} \beta})$ is called a *cutting-plane*, or cut ✂️, because it eliminates the half-space $\{ {\color{darkblack}x} \mid {\color{darkcyan} g^\mathsf{T} } ({\color{darkblack}x} - x_0) + {\color{green} \beta} > 0\}$ from our search. - 🌓 If ${\color{green} \beta}=0$ ($x_0$ is on the boundary of halfspace that is cut), the cutting-plane is called *central cut*. - 🌒 If ${\color{green} \beta}>0$ ($x_0$ lies in the interior of halfspace that is cut), the cutting-plane is called *deep cut*. - 🌔 If ${\color{green} \beta}<0$ ($x_0$ lies in the exterior of halfspace that is cut), the cutting-plane is called *shallow cut*. --- ### Subgradient - $\mathcal{\color{steelblue}K}$ is usually given by a set of inequalities ${\color{purple}f}_j({\color{darkblack}x}) \le 0$ or ${\color{purple}f}_j({\color{darkblack}x}) < 0$ for $j = 1 \cdots m$, where ${\color{purple}f}_j({\color{darkblack}x})$ is a convex function. - A vector ${\color{darkcyan} g} \equiv \partial {\color{purple}f}(x_0)$ is called a subgradient of a convex function $\color{purple}f$ at $x_0$ if ${\color{purple}f}(z) \geq {\color{purple}f}(x_0) + {\color{darkcyan} g^\mathsf{T} } (z - x_0)$. - Hence, the cut $({\color{darkcyan} g}, {\color{green} \beta})$ is given by $(\partial {\color{purple}f}(x_0), {\color{purple}f}(x_0))$ Remarks: 💡 - If ${\color{purple}f}({\color{darkblack}x})$ is differentiable, we can simply take $\partial {\color{purple}f}(x_0) = \nabla {\color{purple}f}(x_0)$ --- ### Key components of Cutting-plane method 🗝️ - A cutting plane oracle $\Omega$ 🔮 - A search space $\mathcal{S}$ initially large enough to cover $\mathcal{\color{steelblue}K}$, e.g. - Polyhedron $\mathcal{P}$ = $\{z \mid C z \preceq d \}$ 🏗️ - Interval $\mathcal{I}$ = $[l, u]$ (for one-dimensional problem) 📏 - Ellipsoid $\mathcal{E}$ = $\{z \mid (z-x_c)P^{-1}(z-x_c) \le 1 \}$ 🥚 --- ### Algorithmic Framework ⚙️ - Initialization 🏁 Begin with an initial ellipsoid that contains the feasible region, defined by its center x⁽⁰⁾ and a positive definite matrix P₍₀₎ that determines the shape and orientation of the ellipsoid. - Oracle Consultation 🔮 At each iteration, present the current ellipsoid center to the separation oracle, which either confirms feasibility or provides a separating hyperplane that cuts through the ellipsoid. - Ellipsoid Update 🔄 Construct a new, smaller ellipsoid that contains the half-ellipsoid created by the cut. The volume shrinks by a factor of approximately $e^{-1/2n}$ at each step, ensuring convergence. --- ### Outline of Cutting-plane method 📝 - Given initial $\mathcal{S}$ known to contain $\mathcal{\color{steelblue}K}$. - Repeat 🔁 - Choose a point $x_0$ in ${\color{red}\mathcal{S} }$ - Query the cutting-plane oracle at $x_0$ - If $x_0 \in \mathcal{\color{steelblue}K}$, quit ✅ - Otherwise, update ${\color{red}\mathcal{S} }$ to a smaller set that covers: $${\color{violet}\mathcal{S}^+} = {\color{red}\mathcal{S} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - x_0) + {\color{green} \beta} \le 0\}$$ - If ${\color{violet}\mathcal{S}^+} = \emptyset$ or it is small enough, quit. --- ### Corresponding Python code 🐍 ```python def cutting_plane_feas(omega, space, options=Options()): for niter in range(options.max_iters): cut = omega.assess_feas(space.xc()) # query the oracle if cut is None: # feasible sol'n obtained return space.xc(), niter status = space.update_deep_cut(cut) # update space if status != CutStatus.Success or space.tsq() < options.tol: return None, niter return None, options.max_iters ``` --- ### Corresponding Sequence Diagram 📈 .pull-left70[ .mermaid[ <pre> sequenceDiagram participant CuttingPlane participant SearchSpace participant OracleFeas loop Cutting Plane Iterations CuttingPlane->>SearchSpace: request xc() SearchSpace-->>CuttingPlane: return xc CuttingPlane->>OracleFeas: assess_feas(xc) alt Feasible (cut is None) OracleFeas-->>CuttingPlane: return None CuttingPlane-->>CuttingPlane: Return solution (xc, niter) else Infeasible OracleFeas-->>CuttingPlane: return cut CuttingPlane->>SearchSpace: update_bias_cut(cut) alt Update successful SearchSpace-->>CuttingPlane: continue else Update failed SearchSpace-->>CuttingPlane: status CuttingPlane-->>CuttingPlane: Return (None, niter) end end end </pre> ] ] .pull-right30[ ] --- ### From Feasibility to Optimization 📈 $$ \begin{array}{ll} \text{minimize} & {\color{purple}f}_0({\color{darkblack}x}), \\ \text{subject to} & {\color{darkblack}x} \in \mathcal{\color{steelblue}K} \end{array} $$ - The optimization problem is treated as a feasibility problem with an additional constraint ${\color{purple}f}_0({\color{darkblack}x}) \le {\color{coral}\gamma}$. - ${\color{purple}f}_0({\color{darkblack}x})$ could be a convex or a *quasiconvex function*. - ${\color{coral}\gamma}$ is also called the *best-so-far* value of ${\color{purple}f}_0({\color{darkblack}x})$. --- ### Convex Optimization Problem - Consider the following general form: $$ \begin{array}{ll} \text{minimize} & {\color{coral}\gamma}, \\ \text{subject to} & {\color{purple}\Phi}({\color{darkblack}x}, {\color{coral}\gamma}) \le 0, \\ & {\color{darkblack}x} \in \mathcal{\color{steelblue}K}, \end{array} $$ where $\mathcal{\color{steelblue}K}'_{\color{coral}\gamma} = \{ {\color{darkblack}x} \mid {\color{purple}\Phi}({\color{darkblack}x}, {\color{coral}\gamma}) \le 0\}$ is the ${\color{coral}\gamma}$-sublevel set of $\{ {\color{darkblack}x} \mid {\color{purple}f}_0({\color{darkblack}x}) \le {\color{coral}\gamma}\}$. - 👉 Note: $\mathcal{\color{steelblue}K}'_{\color{coral}\gamma} \subseteq \mathcal{\color{steelblue}K}'_\epsilon$ if and only if ${\color{coral}\gamma} \le \epsilon$ (monotonicity) - One easy way to solve the optimization problem is to apply the binary search 🔍 on ${\color{coral}\gamma}$. --- ### Binary search ```python def bsearch(omega, intrvl, options=Options()): """Binary search with feasibility oracle. Args: omega: Binary search oracle implementing assess_bs() intrvl: (lower, upper) bound tuple options: Control parameters Returns: (Best g found, iterations used) """ lower, upper = intrvl T = type(upper) # Preserve numerical type (int/float) for niter in range(options.max_iters): tau = (upper - lower) / 2 if tau < options.tolerance: # Convergence check return upper, niter gamma = T(lower + tau) if omega.assess_bs(gamma): # Feasible upper = gamma else: # Infeasible -> move lower bound up lower = gamma return upper, options.max_iters ``` --- ### Binary Search Adaptor ```python class BSearchAdaptor(OracleBS): def __init__(self, omega, space, options=Options()): self.omega = omega # Gamma-sensitive feasibility oracle self.space = space # Search space for subproblems self.options = options # Subproblem solver parameters @property def x_best(self): """Current best feasible solution candidate.""" return self.space.xc() def assess_bs(self, gamma) -> bool: """Test feasibility for given g value.""" space = copy.deepcopy(self.space) # Isolate subproblem space self.omega.update(gamma) # Set current g level x, _ = cutting_plane_feas(self.omega, space, self.options) if x is not None: # Feasible solution found self.space.set_xc(x) # Update main search space return True return False ``` --- ### Shrinking - Another possible way is, to update the best-so-far ${\color{coral}\gamma}$ whenever a feasible solution ${\color{darkblack}x}'$ is found by solving the equation: $${\color{purple}\Phi}({\color{darkblack}x}', {\color{coral}\gamma}_\text{new}) = 0 \, .$$ - If the equation is difficuit to solve but ${\color{coral}\gamma}$ is also convex w.r.t. ${\color{purple}\Phi}$, then we may create a new varaible, say $z$ and let $z \le {\color{coral}\gamma}$. --- ### Outline of Cutting-plane method (Optim) - Given initial $\mathcal{S}$ known to contain $\mathcal{\color{steelblue}K}_{\color{coral}\gamma}$. - Repeat 🔁 - Choose a point $x_0$ in ${\color{red}\mathcal{S} }$ - Query the separation oracle at $x_0$ - If $x_0 \in \mathcal{\color{steelblue}K}_{\color{coral}\gamma}$, update ${\color{coral}\gamma}$ such that ${\color{purple}\Phi}(x_0, {\color{coral}\gamma}) = 0$. - Update ${\color{red}\mathcal{S} }$ to a smaller set that covers: $${\color{violet}\mathcal{S}^+} = {\color{red}\mathcal{S} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - x_0) + {\color{green} \beta} \le 0\} $$ - If ${\color{violet}\mathcal{S}^+} = \emptyset$ or it is small enough, quit. --- ### Corresponding Python code 🐍 ```python def cutting_plane_optim(omega, S, gamma, options=Options()): x_best = None for niter in range(options.max_iters): cut, gamma1 = omega.assess_optim(space.xc(), gamma) if gamma1 is not None: # better \gamma obtained gamma = gamma1 x_best = copy.copy(space.xc()) status = space.update_central_cut(cut) else: status = space.update_deep_cut(cut) if status != CutStatus.Success or space.tsq() < options.tol: return x_best, target, niter return x_best, gamma, options.max_iters ``` --- ### 📚 Example - Profit Maximization Problem 💰 This example is taken from [@Aliabadi2013Robust]. Consider the following *short-run* profit maximization problem: $$ \begin{array}{ll} \text{maximize} & p(A y_1^\alpha y_2^\beta) - v_1 y_1 - v_2 y_2 \\ \text{subject to}& y_1 \le k. \end{array} $$ - $p(A y_1^\alpha y_2^\beta)$ : Cobb-Douglas production function - $p$: the market price per unit - $A$: the scale of production - $\alpha, \beta$: the output elasticities - $y$: input quantity - $v$: output price - $k$: a given constant that restricts the quantity of $y_1$ --- ### 📚 Example - Profit maximization (cont'd) 📈 - The formulation is not in the convex form. - Rewrite the problem in the following form: $$ \begin{array}{ll} \text{maximize} & {\color{coral}\gamma} \\ \text{subject to} & {\color{coral}\gamma} + v_1 y_1 + v_2 y_2 \le p A y_1^{\alpha} y_2^{\beta}\\ & y_1 \le k. \end{array} $$ - Some readers may recognize that the problem can also be written in a geometric program by introducing one additional variable [@Aliabadi2013Robust]. --- ### Profit maximization in Convex Form - By taking the logarithm of each variable: - ${\color{darkblack}x}_1 = \log y_1$, ${\color{darkblack}x}_2 = \log y_2$. - We have the problem in a convex form: $$ \begin{array}{ll} \text{max} & {\color{coral}\gamma} \\ \text{s.t.} & \log({\color{coral}\gamma} + v_1 e^{ {\color{darkblack}x}_1} + v_2 e^{ {\color{darkblack}x}_2}) - (\alpha {\color{darkblack}x}_1 + \beta {\color{darkblack}x}_2) \le \log(pA) \\ & {\color{darkblack}x}_1 \le \log k. \end{array} $$ --- ### Corresponding Python code 🐍 ```python class ProfitOracle: def __init__(self, params, elasticities, price_out): unit_price, scale, limit = params self.log_pA = math.log(unit_price * scale) self.log_k = math.log(limit) self.price_out = price_out self.el = elasticities def assess_optim(self, x, gamma): if (fj := x[0] - self.log_k) > 0.0: # constraint return (np.array([1.0, 0.0]), fj), None log_Cobb = self.log_pA + self.el.dot(x) q = self.price_out * np.exp(x) qsum = q[0] + q[1] if (fj := math.log(gamma + qsum) - log_Cobb) > 0.0: return (q / (gamma + qsum) - self.el, fj), None Cobb = np.exp(log_Cobb) # shrinking return (q / Cobb - self.el, 0.0), Cobb - qsum ``` --- ### Main program 🖥️ ```python import numpy as np from ellalgo.cutting_plane import cutting_plane_optim from ellalgo.ell import Ell from ellalgo.oracles.profit_oracle import ProfitOracle p, A, k = 20.0, 40.0, 30.5 params = p, A, k alpha, beta = 0.1, 0.4 v1, v2 = 10.0, 35.0 el = np.array([alpha, beta]) v = np.array([v1, v2]) r = np.array([100.0, 100.0]) # initial ellipsoid (sphere) ellip = Ell(r, np.array([0.0, 0.0])) omega = ProfitOracle(params, el, v) xbest, gamma, num_iters = cutting_plane_optim(omega, ellip, 0.0) ``` --- class: nord-light, middle, center ## Area of Applications 🌍 --- ### 🛡️ Robust Optimization Applications 🏗️ - 🛡️ Robust Solutions Solutions that remain feasible under parameter uncertainty - ∞ Infinite Constraints Handling constraints that must hold for all uncertainty realizations - Various Uncertainty Sets Polyhedral, ellipsoidal, and interval uncertainty --- ### Network and Matrix Inequality 🌐 - 🖧 Parametric Network Optimization Finding optimal flows or structures in networks with variable parameters. Separation oracles can detect negative cycles using algorithms like Bellman-Ford or Floyd-Warshall, generating cutting planes that exclude infeasible parameter settings. - 𝄜 Matrix Inequality Optimization with symmetric matrix variables constrained to be positive semidefinite. Oracles leverage Cholesky factorization to check positive definiteness and generate cuts when matrices fail this property. - 🦥 Lazy Evaluation Efficiency enhancement where the oracle stops computation as soon as a violation is detected, avoiding unnecessary work and constructing cutting planes with minimal effort. --- class: middle, center ### 🛡️ Robust Convex Optimization 🏗️ --- ### Robust Optimization 🛡️ Robust optimization deals with parameters that belong to uncertainty sets, aiming to find solutions that remain feasible for all possible realizations. The ellipsoid method excels here because separation oracles can efficiently check violations and generate appropriate cutting planes without explicitly enumerating the infinite constraints. --- ### Robust Optimization Formulation 📝 - Consider: $$ \begin{array}{ll} \text{minimize} & \sup_{q \in \mathbb Q} {\color{purple}f}_0({\color{darkblack}x},q), \\ \text{subject to} & {\color{purple}f}_j({\color{darkblack}x},q) \leq 0, \; \forall q \in {\mathbb Q}, \; j = 1,2,\cdots,m, \end{array} $$ where $q$ represents a set of varying parameters. - The problem can be reformulated as: $$ \begin{array}{ll} \text{minimize} & {\color{coral}\gamma}, \\ \text{subject to} & {\color{purple}f}_0({\color{darkblack}x},q) \leq {\color{coral}\gamma} \\ & {\color{purple}f}_j({\color{darkblack}x},q) \leq 0, \; \forall q \in {\mathbb Q}, \; j = 1,2,\cdots,m. \end{array} $$ --- ### 📚 Example - Profit Maximization Problem (convex) 💰 Let us revisit the profit maximization problem. The model parameters are subject to uncertainty over a given interval. Let us now consider the case in which the parameters $\alpha$, $\beta$, $p$, $v_1$, $v_2$, and $k$ are subject to interval uncertainties, as outlined in [@Aliabadi2013Robust]: $$ \begin{array}{ll} \text{max} & {\color{coral}\gamma} \\ \text{s.t.} & \log({\color{coral}\gamma} + \hat{v}_1 e^{ {\color{darkblack}x}_1} + \hat{v}_2 e^{ {\color{darkblack}x}_2}) - (\hat{\alpha} {\color{darkblack}x}_1 + \hat{\beta} {\color{darkblack}x}_2) \le \log(\hat{p}\,A) \\ & {\color{darkblack}x}_1 \le \log \hat{k} , \end{array} $$ - Now assume that: - $\hat{\alpha}$ and $\hat{\beta}$ vary $\bar{\alpha} \pm e_1$ and $\bar{\beta} \pm e_2$ respectively. - $\hat{p}$, $\hat{k}$, $\hat{v}_1$, and $\hat{v}_2$ all vary $\pm e_3$. --- ### 📚 Example - Profit Maximization Problem (oracle) 🔮 By detail analysis, the worst case happens when: - $p = \bar{p} - e_3$, $k = \bar{k} - e_3$ - $v_1 = \bar{v}_1 + e_3$, $v_2 = \bar{v}_2 + e_3$, - if ${\color{darkblack}x}_1 > 0$, $\alpha = \bar{\alpha} - e_1$, else $\alpha = \bar{\alpha} + e_1$ - if ${\color{darkblack}x}_2 > 0$, $\beta = \bar{\beta} - e_2$, else $\beta = \bar{\beta} + e_2$ --- ### Corresponding Python code 🐍 ```python class ProfitRbOracle(OracleOptim): def __init__(self, params, elasticities, price_out, vparams): e1, e2, e3, e4, e5 = vparams self.elasticities = elasticities self.e = [e1, e2] unit_price, scale, limit = params params_rb = unit_price - e3, scale, limit - e4 self.omega = ProfitOracle(params_rb, elasticities, price_out + np.array([e5, e5])) def assess_optim(self, x, gamma): el_rb = copy.copy(self.elasticities) for i in [0, 1]: el_rb[i] += -self.e[i] if x[i] > 0.0 else self.e[i] self.omega.el = el_rb return self.omega.assess_optim(x, gamma) ``` --- ### 🔮 Oracle in Robust Optimization Formulation 🏗️ - The oracle only needs to determine: - If ${\color{purple}f}_j(x_0, q) > 0$ for some $j$ and $q = q_0$, then - the cut $({\color{darkcyan} g}, {\color{green} \beta})$ = $(\partial {\color{purple}f}_j(x_0, q_0), {\color{purple}f}_j(x_0, q_0))$ - If ${\color{purple}f}_0(x_0, q) \geq {\color{coral}\gamma}$ for some $q = q_0$, then - the cut $({\color{darkcyan} g}, {\color{green} \beta})$ = $(\partial {\color{purple}f}_0(x_0, q_0), {\color{purple}f}_0(x_0, q_0) - {\color{coral}\gamma})$ - Otherwise, $x_0$ is feasible, then - Let $q_{\max} = \argmax_{q \in \mathbb Q} {\color{purple}f}_0(x_0, q)$. - ${\color{coral}\gamma} := {\color{purple}f}_0(x_0, q_{\max})$. - The cut $({\color{darkcyan} g}, {\color{green} \beta})$ = $(\partial {\color{purple}f}_0(x_0, q_{\max}), 0)$ Remark: for more complicated problems, affine arithmetic could be used [@liu2007robust]. --- class: middle, center ## Matrix Inequalities --- ### Matrix Inequalities - Cholesky Factorization 🏗️ The ellipsoid method's separation oracle can leverage the Cholesky factorization to efficiently check the positive definiteness of symmetric matrices, a key step in solving semidefinite programming problems. - 🦥 Lazy evaluation ⏱️ By using lazy evaluation techniques, the separation oracle can construct cutting planes in O(p^3) time, where p is the row at which the Cholesky factorization fails, rather than the full O(m^3) time, making the method highly efficient. - Matrix Norm Minimization The ellipsoid method can be applied to matrix norm minimization problems, where the goal is to find the matrix with the smallest norm that satisfies certain constraints. --- ### Problems With Matrix Inequalities Consider the following problem: $$ \begin{array}{ll} \text{find} & {\color{darkblack}x}, \\ \text{subject to} & F({\color{darkblack}x}) \succ 0, \end{array} $$ - $F({\color{darkblack}x})$: a matrix-valued function - $A \succ 0$ denotes $A$ is positive semidefinite. Important application in Control theory 🎛️. --- ### Problems With Matrix Inequalities - Recall that a matrix $A$ is positive semidefinite if and only if $v^\mathsf{T} A v > 0$ for all $v \in \mathbb{R}^N - 0^N$. - The problem can be transformed into: $$\begin{array}{ll} \text{find} & {\color{darkblack}x}, \\ \text{subject to} & v^\mathsf{T} F({\color{darkblack}x}) v > 0, \; \forall v \in \mathbb{R}^N - 0^N \end{array}$$ - Consider $v^\mathsf{T} F({\color{darkblack}x}) v$ is concave for all $v \in \mathbb{R}^N$ w. r. t. ${\color{darkblack}x}$, then the above problem is a convex programming. - Reduce to _semidefinite programming_ if $F({\color{darkblack}x})$ is linear w.r.t. ${\color{darkblack}x}$, i.e., $F({\color{darkblack}x}) = F_0 + {\color{darkblack}x}_1 F_1 + \cdots + {\color{darkblack}x}_n F_n$ --- ### LDLT factorization 🏗️ - The LDLT factorization of a symmetric positive definite matrix $A$ is the factorization $A = L D L^T$, where $L$ is lower triangular with unit diagonal elements and $D$ is a diagonal matrix. - For example, $$ \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 2 \\ 1 & 1 & 3 & 1 \\ 1 & 2 & 1 & 4 \end{array}\right] = \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 \end{array}\right] \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 2 \end{array}\right] \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}\right]. $$ --- ### Storage representation 💾 First, pack the solution and the intermediate storage on a single matrix $T$ such that: $$ t_{ij} = \begin{cases} d_{ii} & \mathrm{if}\ i = j, \\ l_{ij} & \mathrm{if}\ i > j, \\ d_{ii} l_{ji} & \mathrm{if}\ j < i. \end{cases} $$ - For example, $$ T = \left[\begin{array}{cccc} d_{11} & d_{11} l_{21} & d_{11} l_{31} & d_{11} l_{41} \\ l_{21} & d_{22} & d_{22} l_{32} & d_{22} l_{42} \\ l_{31} & l_{32} & d_{33} & d_{33} l_{43} \\ l_{41} & l_{42} & l_{43} & d_{44} \end{array}\right]. $$ --- ### Implementation 💻 - Then, start with $a_{11} = t_{11}$, the implementation of LDLT factorization is: $$ \begin{array}{l} 1~ \text{for}~i=1:n \\ 2~ \quad \text{for}~ j=1:i-1 \\ 3~ \qquad t_{ji} = a_{ij} - \sum_{k=1}^{j-1} t_{ik} t_{jk} \\ 4~ \qquad t_{ij} = t_{ji} / t_{jj} \\ 5~ \quad \text{end} \\ 6~ \quad t_{ii} = a_{ii} - \sum_{k=1}^{i-1} t_{ik} t_{ki} \\ 7~ \text{end} \end{array} $$ - Invoke $\frac{p^3}{2}$ FLOP's (same as Cholesky factorization's), where $p$ is the place the algorithm stops. --- ### Witness of indefiniteness 👁️ - In the case of failure, we can construct a vector $v$ to certify that $v^T A v \leq 0$. - Let $L_{p,p}$ denote the partial sub-matrix of $L$ where $p$ is the row of failure. - Then $v = [L_{p,p}^{-T} e_p, 0, \cdots, 0]^T$, where $e_p = [0, \cdots, 0, 1]^T \in \mathbb{R}^p$ - Start with $v = e_p$, the basic algorithm is: $$ \begin{array}{l} 1~ \text{for}~i = p - 1~\text{downto}~1 \\ 2~ \quad \text{for}~ k = i~\text{to}~p \\ 3~ \qquad v_i = v_i - t_{k,i} v_k \\ 4~ \quad \text{end} \\ 5~ \text{end} \end{array} $$ --- ### Oracle in Matrix Inequalities 🔮 The oracle only needs to: - Perform a _row-based_ LDLT factorization such that $F(x_0) = L D L^\mathsf{T}$. - Let $A_{p,p}$ denotes a submatrix of $A \in \mathbb{R}^{p\times p}$. - If the process fails at row $p$, - there exists a vector $e_p = (0, 0, \cdots, 0, 1)^\mathsf{T} \in \mathbb{R}^p$, such that - $v = R_{p,p}^{-1} e_p$, and - $v^\mathsf{T} F_{p,p}(x_0) v \leq 0$. - The cut $({\color{darkcyan} g}, {\color{green} \beta})$ = $(-v^\mathsf{T} \partial F_{p,p}(x_0) v, -v^\mathsf{T} F_{p,p}(x_0) v)$ --- ### 🦥 Lazy evaluation ⏱️ - Don't construct the full matrix at each iteration! - Only O($p^3$) per iteration, independent of $N$! --- ### Python Implementation 🐍 ```python class LMIOracle: def __init__(self, F, B): self.F = F self.F0 = B self.Q = LDLTMgr(len(B)) def assess_feas(self, x: Arr) -> Optional[Cut]: def get_elem(i, j): return self.F0[i, j] - sum( Fk[i, j] * xk for Fk, xk in zip(self.F, x)) if self.Q.factor(get_elem): return None ep = self.Q.witness() g = np.array([self.Q.sym_quad(Fk) for Fk in self.F]) return g, ep ``` --- ### Google Benchmark 📊 Comparison ```terminal 2: ---------------------------------------------------------- 2: Benchmark Time CPU Iterations 2: ---------------------------------------------------------- 2: BM_LMI_Lazy 131235 ns 131245 ns 4447 2: BM_LMI_old 196694 ns 196708 ns 3548 2/4 Test #2: Bench_BM_lmi ..................... Passed 2.57 sec ``` --- ### 📚 Example - Matrix Norm Minimization - Let $A({\color{darkblack}x}) = A_0 + {\color{darkblack}x}_1 A_1 + \cdots + {\color{darkblack}x}_n A_n$ - Problem $\min_x \| A({\color{darkblack}x}) \|$ can be reformulated as $$ \begin{array}{ll} \text{minimize} & {\color{coral}\gamma}, \\ \text{subject to} & \left( \begin{array}{cc} {\color{coral}\gamma}\,I & A({\color{darkblack}x}) \\ A^\mathsf{T}({\color{darkblack}x}) & {\color{coral}\gamma}\,I \end{array} \right) \succ 0, \end{array} $$ - Binary search 🔍 on ${\color{coral}\gamma}$ can be used for this problem. --- class: middle, center ### Multi-parameter Network Problem --- ### Parametric Network Optimization - Negative Cycle Detection 🔍 The key to solving parametric network optimization problems with the ellipsoid method is the ability to efficiently detect negative cycles in the network. Algorithms like Tarjan's can exploit the network's structure to quickly identify cutting planes. - Exploiting Locality 🏘️ By taking advantage of the network's locality and other properties, the separation oracle can construct cutting planes in a highly efficient manner, further enhancing the performance of the ellipsoid method. - Diverse Applications 🌍 The ellipsoid method extends the single-parameter network problems to multi-parameter, makes it a valuable tool in many engineering and operations research domains. --- ### 🌐 Parametric Network Problem Given a network represented by a directed graph $G = ({\color{salmon}V}, {\color{lime}E})$. Consider: $$ \begin{array}{ll} \text{find} & {\color{darkblack}x}, {\color{red}u} \\ \text{subject to} & {\color{red}u_j} - {\color{red}u_i} \le {\color{purple}h}_{ij}({\color{darkblack}x}), \; \forall (i, j) \in {\color{lime}E} , \end{array} $$ - ${\color{purple}h}_{ij}({\color{darkblack}x})$ is the concave function of edge $(i,j)$, - Assume: network is large, but the number of parameters is small. --- ### 🔄 Network Potential Problem (cont'd) Given ${\color{darkblack}x}$, the problem has a feasible solution if and only if $G$ contains no negative cycle. Let $\mathcal{C}$ be a set of all cycles of $G$. $$ \begin{array}{ll} \text{find} & {\color{darkblack}x} \\ \text{subject to} & {\color{purple}w}_k({\color{darkblack}x}) \ge 0, \forall {\color{lime}C}_k \in \mathcal{C} , \end{array} $$ - ${\color{lime}C}_k$ is a cycle of $G$ - ${\color{purple}w}_k({\color{darkblack}x}) = \sum_{ (i,j)\in {\color{lime}C}_k} {\color{purple}h}_{ij}({\color{darkblack}x})$. --- ### ⚠️ Negative Cycle Finding There are lots of methods to detect negative cycles in a weighted graph [@cherkassky1999negative], in which Tarjan's algorithm [@Tarjan1981negcycle] is one of the fastest algorithms in practice [@alg:dasdan_mcr; @cherkassky1999negative]. --- ### 🧙 Oracle in Network Potential Problem - The oracle only needs to determine: - 🟢 If there exists a negative cycle ${\color{lime}C}_k$ under $x_0$, then - the cut $({\color{darkcyan} g}, {\color{green} \beta})$ = $(-\partial {\color{purple}w}_k(x_0), -{\color{purple}w}_k(x_0))$ - 🔵 Otherwise, the shortest path solution gives the value of ${\color{red}u}$. --- ### 🐍 Python Code ```python class NetworkOracle: def __init__(self, G, u, h): self._G = G self._u = u self._h = h self._S = NegCycleFinder(G) def update(self, gamma): self._h.update(gamma) def assess_feas(self, x) -> Optional[Cut]: def get_weight(e): return self._h.eval(e, x) for Ci in self._S.find_neg_cycle(self._u, get_weight): f = -sum(self._h.eval(e, x) for e in Ci) g = -sum(self._h.grad(e, x) for e in Ci) return g, f # use the first Ci only return None ``` --- ### 📚 Example - Optimal Matrix Scaling [@orlin1985computing] - Given a sparse matrix $A = [a_{ij}] \in \mathbb{R}^{N\times N}$. - Find another matrix $B = U A U^{-1}$ where $U$ is a nonnegative diagonal matrix, such that the ratio of any two elements of $B$ in absolute value is as close to 1 as possible. - Let $U = \mathrm{diag}([u_1, u_2, \ldots, u_N])$. Under the min-max-ratio criterion, the problem can be formulated as: $$ \begin{array}{ll} \text{minimize} & \pi/\psi \\ \text{subject to} & \psi \le u_i |a_{ij}| u_j^{-1} \le \pi, \; \forall a_{ij} \neq 0 , \\ & \pi, \psi, u, \, \text{positive} \\ \text{variables} & \pi, \psi, u \, . \end{array} $$ --- ### 📚 Optimal Matrix Scaling (cont'd) By taking the logarithms of variables, the above problem can be transformed into: $$ \begin{array}{ll} \text{minimize} & {\color{coral}\gamma} \\ \text{subject to} & {\color{darkblack}\pi'} - {\color{darkblack}\psi'} \le {\color{coral}\gamma} \\ & {\color{red}u_i'} - {\color{red}u_j'} \le {\color{darkblack}\pi'} - a_{ij}', \; \forall a_{ij} \neq 0 \,, \\ & {\color{red}u_j'} - {\color{red}u_i'} \le a_{ij}' - {\color{darkblack}\psi'}, \; \forall a_{ij} \neq 0 \,, \\ \text{variables} & {\color{darkblack}\pi'}, {\color{darkblack}\psi'}, {\color{red}u'} \, . \end{array} $$ where $k'$ denotes $\log( | k | )$ and $x = ({\color{darkblack}\pi'}, {\color{darkblack}\psi'} )^\mathsf{T}$. --- ```python class OptScalingOracle: class Ratio: def __init__(self, G, get_cost): self._G = G self._get_cost = get_cost def eval(self, e, x: Arr) -> float: u, v = e cost = self._get_cost(e) return x[0] - cost if u < v else cost - x[1] def grad(self, e, x: Arr) -> Arr: u, v = e return np.array([1.0, 0.0] if u < v else [0.0, -1.0]) def __init__(self, G, u, get_cost): self._network = NetworkOracle(G, u, self.Ratio(G, get_cost)) def assess_optim(self, x: Arr, gamma: float): s = x[0] - x[1] g = np.array([1.0, -1.0]) if (fj := s - gamma) >= 0.0: return (g, fj), None if (cut := self._network.assess_feas(x)): return cut, None return (g, 0.0), s ``` --- ### 📚 Example - clock period & yield-driven co-optimization $$ \begin{array}{cll} \text{minimize} &{\color{darkblack}T_\text{CP} } / {\color{darkblack}\beta} \\ \text{subject to} & {\color{red}u_i} - {\color{red}u_j} \le {\color{darkblack}T_\text{CP} } - F_{ij}^{-1}({\color{darkblack}\beta}), & \forall (i,j) \in {\color{lime}E}_s \,,\\ & {\color{red}u_j} - {\color{red}u_i} \le F_{ij}^{-1}(1 - {\color{darkblack}\beta}), & \forall (j,i) \in {\color{lime}E}_h \,, \\ & {\color{darkblack}T_\text{CP} } \ge 0, \, 0 \le {\color{darkblack}\beta} \le 1 \, , \\ \text{variables} &{\color{darkblack}T_\text{CP} }, {\color{darkblack}\beta}, {\color{red}u}. \end{array} $$ - 👉 Note that $F_{ij}^{-1}(x)$ is not concave in general in $[0, 1]$. - Fortunately, we are most likely interested in optimizing circuits for high yield rather than the low one in practice. - Therefore, by imposing an additional constraint to ${\color{darkblack}\beta}$, say ${\color{darkblack}\beta} \geq 0.8$, the problem becomes convex. --- ### 📚 Example - clock period & yield-driven co-optimization The problem can be reformulated as: $$ \begin{array}{cll} \text{minimize} & {\color{coral}\gamma} \\ \text{subject to} & {\color{darkblack}T_\text{CP} } - {\color{darkblack}\beta} {\color{coral}\gamma} \le 0\\ & {\color{red}u_i} - {\color{red}u_j} \le {\color{darkblack}T_\text{CP} } - F_{ij}^{-1}({\color{darkblack}\beta}), & \forall (i,j) \in {\color{lime}E}_s \,,\\ & {\color{red}u_j} - {\color{red}u_i} \le F_{ij}^{-1}(1 - {\color{darkblack}\beta}), & \forall (j,i) \in {\color{lime}E}_h \,, \\ & {\color{darkblack}T_\text{CP} } \ge 0, \, 0.8 \le {\color{darkblack}\beta} \le 1 \, , \\ \text{variables} & {\color{darkblack}T_\text{CP} }, {\color{darkblack}\beta}, {\color{red}u}. \end{array} $$ --- ### 📚 Example - clock period & slack co-optimization The problem can be formulated as: $$ \begin{array}{cll} \text{minimize} & {\color{coral}\gamma} \\ \text{subject to} & {\color{darkblack}T_\text{CP} } - {\color{darkblack}\beta} {\color{coral}\gamma} \le 0\\ & {\color{red}u_i} - {\color{red}u_j} \le {\color{darkblack}T_\text{CP} } - D_{ij} - {\color{darkblack}\beta}, & \forall (i,j) \in {\color{lime}E}_s \,,\\ & {\color{red}u_j} - {\color{red}u_i} \le d_{ij} - {\color{darkblack}\beta}, & \forall (j,i) \in {\color{lime}E}_h \,, \\ & {\color{darkblack}T_\text{CP} } \ge 0, \\ \text{variables} & {\color{darkblack}T_\text{CP} }, {\color{darkblack}\beta}, {\color{red}u}. \end{array} $$ --- ### Example TCG .mermaid[ <pre> graph LR A -- "7" --> E(("v2")) A -. "6" .-> E A(("v0")) -- "6" --> B(("v3")) A -. "6" .-> B D(("v1")) -- "9" --> E D -. "6" .-> E E -- "6" --> B E -. "6" .-> B B -- "8" --> C(("v4")) B -. "8" .-> C C -- "3" --> D C -. "3" .-> D %% Caption: Timing Constraint Graph (TCG). </pre> ] --- ### NetOptim Example ```python class MyOracle(OracleOptim[np.ndarray]): def __init__(self, dist: Dict[str, int], digraph): self.dist = dist self.finder = NegCycleFinder(digraph) def assess_optim(self, xc, gamma): TCP, beta = xc if TCP < 0.0: return (np.array([-1.0, 0.0]), -TCP), None if (fj := TCP - beta * gamma) > 0: return (np.array([1.0, -gamma]), fj), None def calc_weight(e): return (TCP - e["delay"] - beta if e["type"] == "s" else e["delay"] - beta ) for cycle in self.finder.howard(self.dist, calc_weight): f = -sum(calc_weight(e) for e in cycle) g = np.array([-sum(1.0 if e["type"] == "s" else 0.0 for e in cycle), len(cycle) ]) return (g, f), None # use the first cycle only return (np.array([1.0, -1.0]), 0.0), TCP / beta ``` --- ### NetOptim Example ```python digraph: Dict[str, Dict[str, Dict[str, Any]]] = { "v0": {"v3": {"type": "s", "delay": 6}, "v2": {"type": "s", "delay": 7} }, "v1": {"v2": {"type": "s", "delay": 9}, "v4": {"type": "h", "delay": 3} }, "v2": { "v0": {"type": "h", "delay": 6}, "v1": {"type": "h", "delay": 6}, "v3": {"type": "s", "delay": 6} }, "v3": { "v4": {"type": "s", "delay": 8}, "v0": {"type": "h", "delay": 6}, "v2": {"type": "h", "delay": 6} }, "v4": {"v1": {"type": "s", "delay": 3}, "v3": {"type": "h", "delay": 8} }, } xinit = np.array([7.5, 1.0]) ellip = Ell(100, xinit) dist = {"v0": 0, "v1": 0, "v2": 0, "v3": 0, "v4": 0} omega = MyOracle(dist, digraph) xbest, _, _ = cutting_plane_optim(omega, ellip, float("inf")) assert xbest is not None ``` --- class: nord-light, middle, center ## 💻 Implementation Details of Ellipsoid Method --- ### 🟡 Basic Ellipsoid Method - An ellipsoid $\mathcal{E}(x_c, P)$ is specified as a set $$\{x \mid (x-x_c)P^{-1}(x-x_c) \le 1 \},$$ where $x_c$ is the center of the ellipsoid.  --- ### 🔁 Updating the ellipsoid Calculation of minimum volume ellipsoid ${\color{violet} \mathcal{E}^+}$ covering: $$ {\color{red} \mathcal{E} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) + {\color{green} \beta} \le 0 \}. $$  --- ### 🌒 Deep-cut - Let $\tilde{g} = {\color{red} P}\,{\color{darkcyan} g}$, ${\color{brown} \tau^2} = {\color{darkcyan} g^\mathsf{T} } {\color{red} P} {\color{darkcyan} g}$. - If ${\color{blue} \tau} + n \cdot {\color{green} \beta} \leq 0$ (shallow cut 🌔), no smaller ellipsoid can be found. - If ${\color{green} \beta} > {\color{blue} \tau}$, intersection is empty. - Otherwise, $$ {\color{violet} x_c^+} = {\color{red} x_c} - \frac{\rho}{\color{brown} \tau^2 } \tilde{g}, \quad {\color{violet} P^+} = \delta\cdot\left({\color{red} P} - \frac{\sigma}{\color{brown} \tau^2} \tilde{g}\tilde{g}^\mathsf{T}\right). $$ where $$ \rho = \frac{ {\color{blue} \tau} + n \cdot {\color{green} \beta} }{n+1}, \quad \sigma = \frac{2\rho}{ {\color{blue} \tau} + {\color{green} \beta} }, \quad \delta = \frac{n^2({\color{blue} \tau} + {\color{green} \beta})({\color{blue} \tau} - {\color{green} \beta})}{(n^2 - 1){\color{brown} \tau^2} }. $$ - Example: - If $n = 4$, ${\color{blue} \tau} = 0.1$, and ${\color{green} \beta} = 0.05$, - then $\rho = 0.06$, $\sigma = 0.8$, and $\delta = 0.8$. --- ### Updating the ellipsoid (cont'd) - Even better, split ${\color{red} P}$ into two variables ${\color{red} \kappa} \cdot {\color{red} Q}$ - Let $\tilde{g} = {\color{red} Q} \cdot {\color{darkcyan} g}$, $\omega = {\color{darkcyan} g^\mathsf{T} }\tilde{g}$, ${\color{brown} \tau^2} = {\color{red} \kappa}\cdot\omega$. $$ {\color{violet} x_c^+} = {\color{red} x_c} - \frac{\rho}{\omega} \tilde{g}, \quad {\color{violet} Q^+} = {\color{red} Q} - \frac{\sigma}{\omega} \tilde{g}\tilde{g}^\mathsf{T}, \quad {\color{violet} \kappa^+} = \delta\cdot{\color{red} \kappa}. $$ - Reduce $n^2$ multiplications per iteration. - 👉 Note: - The determinant of $Q$ decreases monotonically. - The range of $\delta$ is $(0, \frac{n^2}{n^2 - 1})$. --- ### 🌓 Central Cut Calculation of minimum volume ellipsoid ${\color{violet} \mathcal{E}^+}$ covering: $$ {\color{red} \mathcal{E} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) \le 0 \}. $$  --- ### 🌓 Central Cut - A Special case of deep cut when ${\color{green} \beta} = 0$ - Deserve a separate implement because it is much simplier. - Update the ellipsoid with: $$ \rho = \frac{\color{blue} \tau}{n+1}, \qquad \sigma = \frac{2}{n+1}, \qquad \delta = \frac{n^2}{n^2 - 1}. $$ - Example: - If $n = 4$ and ${\color{blue} \tau} = 0.1$, - then $\rho = 0.02$, $\sigma = 0.4$, and $\delta = 16/15$. --- class: middle, center ### ⏸️ Parallel Cuts --- ### ⏸️ Parallel Cuts - Oracle returns a pair of cuts instead of just one. - The pair of cuts is given by ${\color{darkcyan} g}$ and $({\color{green} \beta_0}, {\color{blue} \beta_1})$ such that: $$ \begin{array}{l} {\color{darkcyan} g^\mathsf{T} } (x - {\color{red} x_c}) + {\color{green} \beta_0} \leq 0, \\ {\color{darkcyan} g^\mathsf{T} } (x - {\color{red} x_c}) + {\color{blue} \beta_1} \geq 0, \end{array} $$ for all $x \in \mathcal{\color{steelblue}K}$. - Only linear inequality constraint can produce such parallel cut: $$ l \le a^\mathsf{T} x + b \le u, \quad L \preceq F(x) \preceq U. $$ - Usually provide faster convergence. --- ### ⏸️ Parallel Cuts Calculation of minimum volume ellipsoid ${\color{violet} \mathcal{E}^+}$ covering: $$ {\color{red} \mathcal{E} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) + {\color{green} \beta_0} \le 0 \land {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) + {\color{blue} \beta_1} \ge 0 \}. $$  --- ### Updating the ellipsoid (old) - If ${\color{blue} \beta_1} \lt {\color{green} \beta_0}$, intersection is empty. - If ${\color{brown} \tau^2} + n {\color{green} \beta_0} {\color{blue} \beta_1} \leq 0$, no smaller ellipsoid can be found. - If ${\color{blue} \beta_1^2} \geq {\color{brown} \tau^2}$, it reduces to deep-cut with ${\color{green} \beta} = {\color{green} \beta_0}$ - Otherwise, update the ellipsoid with: $$ \begin{array}{lll} \sigma &=& (n + (2{\color{brown} \tau^2} + 2{\color{green} \beta_0} {\color{blue} \beta_1} - \xi)/{\color{red}(\beta_0 + \beta_1)^2} ) / (n + 1), \\ \rho &=& \sigma({\color{green} \beta_0} + {\color{blue} \beta_1}) / 2, \\ \delta &=& (n^2/2(n^2-1)) (\zeta_0 + \zeta_1 + \xi/n) / {\color{brown} \tau^2} . \end{array} $$ where $$ \begin{array}{lll} \zeta_0 &=& {\color{brown} \tau^2} - {\color{green} \beta_0^2} \\ \zeta_1 &=& {\color{brown} \tau^2} - {\color{blue} \beta_1^2} \\ \xi &=& \sqrt{4\zeta_0\zeta_1 + n^2({\color{blue} \beta_1^2} - {\color{green} \beta_0^2})^2}, \\ \end{array} $$ - Example: - If $n = 4$, ${\color{brown} \tau^2} = 0.01$, ${\color{green} \beta_0} = 0.01$, and ${\color{blue} \beta_1} = 0.04$, - then $\rho = 0.0232$, $\sigma = 0.928$, and $\delta = 1.232$. --- ### Updating the ellipsoid (old) - Break down when ${\color{green} \beta_0} + {\color{blue} \beta_1} = 0$  --- ### Updating the ellipsoid (new) - If ${\color{blue} \beta_1} \lt {\color{green} \beta_0}$, intersection is empty. - If ${\color{brown} \tau^2} + n {\color{green} \beta_0} {\color{blue} \beta_1} \leq 0$, no smaller ellipsoid can be found. - If ${\color{blue} \beta_1^2} \geq {\color{brown} \tau^2}$, it reduces to deep-cut with ${\color{green} \beta} = {\color{green} \beta_0}$ - Otherwise, update the ellipsoid with: $$ \sigma = \frac{\eta}{k}, \quad \rho = \sigma \bar{\beta}, \quad \delta = 1 + \frac{\eta}{ {\color{brown} \tau^2}(k - \eta)} (\bar{\beta}^2 \sigma - {\color{green} \beta_0}{\color{blue} \beta_1}). $$ where $$ \begin{array}{lll} \bar{\beta} &=& ({\color{green} \beta_0} + {\color{blue} \beta_1}) / 2 \\ \eta &=& {\color{brown} \tau^2} + n {\color{green} \beta_0} {\color{blue} \beta_1} \\ h &=& ({\color{brown} \tau^2} + {\color{green} \beta_0}{\color{blue} \beta_1})/2 + n \bar{\beta}^2, \\ k &=& h + \sqrt{h^2 - (n + 1) \eta \bar{\beta}^2}, \\ \end{array} $$ - Example: - If $n = 4$, ${\color{brown} \tau^2} = 1.0$, ${\color{green} \beta_0} = -0.25$, and ${\color{blue} \beta_1} = 0.25$, - then $\rho = 0$, $\sigma = 0.8$, and $\delta = 1.25$. --- ### Parallel Central Cuts Calculation of minimum volume ellipsoid ${\color{violet} \mathcal{E}^+}$ covering: $$ {\color{red} \mathcal{E} } \cap \{z \mid {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) \le 0 \land {\color{darkcyan} g^\mathsf{T} } (z - {\color{red} x_c}) + {\color{blue} \beta_1} \ge 0 \}. $$  --- ### Updating the ellipsoid - If ${\color{blue} \beta_1} \lt {\color{green} 0}$, intersection is empty. - If ${\color{blue} \beta_1^2} \geq {\color{brown} \tau^2}$, it reduces to central-cut - Otherwise, update the ellipsoid with: $$ \rho = \frac{\color{blue} \beta_1}{r + 1}, \quad \sigma = \frac{2}{r + 1}, \quad \delta = \frac{r}{r - 1/n}. $$ where $$ \begin{array}{lll} \alpha^2 &=& {\color{blue} \beta_1^2} / {\color{brown} \tau^2} \\ h &=& n \alpha^2 / 2 \\ r &=& h + \sqrt{h^2 + 1 - \alpha^2} \\ \end{array} $$ - Example: - If $n = 4$, ${\color{blue} \beta_1} = 1.0$, and ${\color{brown} \tau^2} = 4.0$, - then $\rho = 0.4$, $\sigma = 0.8$, and $\delta = 1.2$. --- ### Normal flows .mermaid[ <pre> graph LR S{central cut?} -- Yes --> A{parallel?} S -- No --> B{parallel?} A -- Yes --> L{β1 < 0} A -- No --> D[central cut] L -- No --> C{β1 >= τ} C -- No --> E[parallel central] C -- Yes --> D B -- No --> F{τ < β0} B -- Yes --> G{β1 < β0} F -- Yes --> H([NoSoln]) L -- Yes --> H([NoSoln]) F -- No --> I[deep cut] G -- Yes --> H G -- No --> J{β1 >= τ} J -- Yes --> I J -- No --> K[parallel deep] </pre> ] --- ### 📚 Example - FIR filter design  - The time response is: $$y[t] = \sum_{k=0}^{n-1}{h[k]u[t-k]}. $$ --- ### 📚 Example - FIR filter design (cont'd) - The frequency response: $$H(\omega)~=~\sum_{m=0}^{n-1}{h(m)e^{-jm\omega} }. $$ - The magnitude constraints on frequency domain are expressed as $$L(\omega)~\leq~|H(\omega)|~\leq~U(\omega),~\forall~\omega\in(-\infty,+\infty). $$ where $L(\omega)$ and $U(\omega)$ are the lower and upper (nonnegative) bounds at frequency $\omega$ respectively. - The constraint is non-convex in general. --- ### 📚 Example - FIR filter design (II) - However, via *spectral factorization* [@goodman1997spectral], it can transform into a convex one [@wu1999fir]: $$L^2(\omega)~\leq~R(\omega)~\leq~U^2(\omega),~\forall~\omega\in(0,\pi), $$ where - $R(\omega)=\sum_{i=-1+n}^{n-1}{ {\color{darkblack}r}(t)e^{-j{\omega}t} }=|H(\omega)|^2$ - $\mathbf{ {\color{darkblack}r} }=({\color{darkblack}r}(-n+1),{\color{darkblack}r}(-n+2),...,{\color{darkblack}r}(n-1))$ are the autocorrelation coefficients. --- ### 📚 Example - FIR filter design (III) - $\mathbf{ {\color{darkblack}r} }$ can be determined by $\mathbf{h}$: $${\color{darkblack}r}(t)~=~\sum_{i=-n+1}^{n-1}{h(i)h(i+t)},~t\in\mathbf{Z}, $$ where $h(t)=0$ for ${\color{coral}\gamma} < 0$ or ${\color{coral}\gamma} > n - 1$. - The whole problem can be formulated as: $$ \begin{array}{ll} \text{min} & {\color{coral}\gamma} \\ \text{s.t.} & L^2(\omega) \le R(\omega) \le U^2(\omega), \; \forall \omega \in [0,\pi] \\ & R(\omega) > 0, \forall \omega \in [0,\pi] \end{array} $$ --- ### 🧪 Experiment  --- ### 📊 Google Benchmark Result ```terminal 3: ------------------------------------------------------------------ 3: Benchmark Time CPU Iterations 3: ------------------------------------------------------------------ 3: BM_Lowpass_single_cut 627743505 ns 621639313 ns 1 3: BM_Lowpass_parallel_cut 30497546 ns 30469134 ns 24 3/4 Test #3: Bench_BM_lowpass ................. Passed 1.72 sec ``` --- class: middle, center ## 🧩 Discrete Optimization --- ### Why Discrete Convex Programming - Many engineering problems can be formulated as a convex/geometric programming, e.g. digital circuit sizing - Yet in an ASIC design, often there is only a limited set of choices from the cell library. In other words, some design variables are discrete. - The discrete version can be formulated as a _Mixed-Integer Convex programming_ (MICP) by mapping the design variables to integers. --- ### What's Wrong w/ Existing Methods? - Mostly based on relaxation. - Then use the relaxed solution as a lower bound and use the branch--and--bound method for the discrete optimal solution. - 👉 Note: the branch-and-bound method does not utilize the convexity of the problem. - What if I can only evaluate constraints on discrete data? Workaround: convex fitting? --- ### Applying the Ellipsoid Method - Handling Discrete Variables The cutting-plane method, of which the ellipsoid method is a part, can also be applied to discrete optimization problems where some design variables are restricted to discrete forms. The key is to use a separation oracle that can locate the nearest discrete solutions. - Overcoming Limitations Unlike traditional methods that rely on relaxation and branch-and-bound, the ellipsoid method can exploit the convexity of the problem, even in the presence of discrete constraints, leading to more efficient solutions. - Applications in Engineering Discrete optimization problems arise in many engineering domains, such as the design of multiplierless FIR filters, where the coefficients must be represented using a limited set of discrete values. The ellipsoid method offers a promising approach for tackling these challenges. --- ### Mixed-Integer Convex Programming Consider: $$ \begin{array}{ll} \text{minimize} & {\color{purple}f}_0(x), \\ \text{subject to} & {\color{purple}f}_j(x) \le 0, \; \forall j=1,2,\ldots \\ & x \in \mathbb{D} \end{array} $$ where - ${\color{purple}f}_0(x)$ and ${\color{purple}f}_j(x)$ are "convex" - Some design variables are discrete. --- ### 🔮 Oracle Requirement - The oracle looks for the nearby discrete solution ${\color{magenta} x_d}$ of ${\color{orange} x_c}$ with the cutting-plane: $${\color{darkcyan} g^\mathsf{T} } (x - {\color{magenta}x_d}) + {\color{green} \beta} \le 0, {\color{green} \beta} \ge 0, {\color{darkcyan} g} \neq 0$$ - 👉 Note: the cut may be a shallow cut 🌔. - Suggestion: use different cuts as possible for each iteration (e.g. round-robin the evaluation of constraints) --- ### 🧩 Discrete Cut  --- ### Discrete flows .mermaid[ <pre> graph LR B{parallel?} -- No --> F{τ < β0} B -- Yes --> G{β1 < β0} F -- Yes --> H([NoSoln]) F -- No --> A{τ + n β0 <= 0} A -- Yes --> C([NoEffect]) A -- No --> I[bias cut] G -- Yes --> H([NoSoln]) G -- No --> D{τ^2 + n β0 β1 <= 0} D -- Yes --> C([NoEffect]) D -- No --> J{β1 >= τ} J -- Yes --> I J -- No --> K[parallel bias] </pre> ] --- ### 📚 Multiplier-less FIR filter design (nnz=3)  --- ### Complementary Role in Optimization - Ellipsoid Method Strengths Efficiently handles massive constraint sets through separation oracles. Particularly valuable for robust optimization, parametric problems, and certain classes of semidefinite and discrete optimization. - Interior-Point Method Strengths Better practical performance for well-structured convex problems with manageable constraint counts. Works by traversing the interior of the feasible region rather than shrinking enclosing volumes. --- ### Companion 👫, Not Competitor 🤼‍♂️ - Complementary Strengths While the ellipsoid method may be perceived as slower than interior-point methods for solving convex problems, it offers distinct advantages, such as the ability to handle problems with a large or infinite number of constraints. - Ongoing Improvements Techniques like parallel cuts and efficient implementations have helped to improve the performance of the ellipsoid method, making it a valuable tool in the optimization landscape. - A Collaborative Approach Rather than viewing the ellipsoid method as a competitor to other optimization techniques, it should be seen as a companion, with each method offering unique strengths that can be leveraged to solve a wide range of optimization problems effectively. --- ### Conclusion: The Enduring Value of Amazing Oracles The ellipsoid method remains a cornerstone of modern optimization, with its true versatility lying in the elegant framework that leverages separation oracles. These "amazing oracles" enable tackling previously intractable challenges across robust optimization, network analysis, matrix inequalities, and discrete optimization. By understanding the strengths of the ellipsoid method and the crucial role of its oracles, researchers and practitioners can effectively leverage its capabilities to address complex optimization problems that would otherwise remain beyond reach. --- class: nord-dark, middle, center ## Q&A 🎤