layout: true class: typo, typo-selection --- count: false class: nord-dark, middle, center # 🥥 The Ellipsoid Method and Its Amazing Oracles (Part 1) @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 ] --- class: nord-light, middle, center ## 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. --- ## 🎬 Introduction --- ### Common Perspective of Ellipsoid Method - It is widely believed to be inefficient in practice for large-scale problems. - Convergent rate is slow, even when using deep cuts. - Cannot exploit sparsity. - It has since then supplanted by the interior-point methods. - Used only as a theoretical tool to prove polynomial-time solvability of some combinatorial optimization problems. --- #### 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. --- class: middle, center ## 🥥 Cutting-plane Method Revisited --- ### Convex Set - Let $\mathcal{K} \subseteq \mathbb{R}^n$ be a convex set 🥚. - Consider the feasibility problem: - Find a point ${\color{darkmagenta}x}^* \in \mathbb{R}^n$ in $\mathcal{K}$, - or determine that $\mathcal{K}$ is empty (i.e., there is no feasible solution)  .mermaid[ <pre> graph LR A[Convex Set 𝒦] --> B[Search Space S] B --> C[Cutting Plane] C --> D[Smaller S⁺] D --> E[Repeat Until Found] style A fill:#4caf50 style B fill:#ff9800 style C fill:#f44336 style D fill:#2196f3 style E fill:#9c27b0 </pre> ] --- #### 🔮 Separation Oracle .column-2.column-norule[ - When a separation oracle $\Omega$ is *queried* at $x_0$, it either - asserts that $x_0 \in \mathcal{K}$, or - returns a separating hyperplane between $x_0$ and $\mathcal{K}$: $$ g^\mathsf{T} ({\color{darkmagenta}x} - x_0) + \beta \le 0, \\ \beta \geq 0, g \neq 0, \; \forall {\color{darkmagenta}x} \in \mathcal{K} $$  ] .mermaid[ <pre> graph TD A[Query Point x₀] --> B{Separation Oracle Ω} B -->|"x₀ ∈ 𝒦"| C["✅ Feasible"] B -->|"x₀ ∉ 𝒦"| D[Return Hyperplane] D --> E["gᵀ(x - x₀) + β ≤ 0"] E --> F["Cut off half-space"] style A fill:#ff9800 style B fill:#f44336 style C fill:#4caf50 style D fill:#ff9800 style E fill:#2196f3 style F fill:#9c27b0 </pre> ] --- #### 🔮 Separation Oracle (cont'd) - $(g, \beta)$ is called a *cutting-plane*, or cut, because it eliminates the half-space $\{ {\color{darkmagenta}x} \mid g^\mathsf{T} ({\color{darkmagenta}x} - x_0) + \beta > 0\}$ from our search. - 🌓 If \(\beta=0\) (\(x_0\) is on the boundary of halfspace that is cut), the cutting-plane is called *central cut*. - 🌒 If \(\beta>0\) (\(x_0\) lies in the interior of halfspace that is cut), the cutting-plane is called *deep cut*. - 🌔 If \(\beta<0\) (\(x_0\) lies in the exterior of halfspace that is cut), the cutting-plane is called *shallow cut*. .mermaid[ <pre> graph LR subgraph "Deep Cut β > 0" DC["x₀ inside → more removed"] end subgraph "Central Cut β = 0" CC["x₀ on boundary → half removed"] end subgraph "Shallow Cut β < 0" SC["x₀ outside → less removed"] end style DC fill:#f44336 style CC fill:#ff9800 style SC fill:#2196f3 </pre> ] --- #### Subgradient - $\mathcal{K}$ is usually given by a set of inequalities $f_j({\color{darkmagenta}x}) \le 0$ or $f_j({\color{darkmagenta}x}) < 0$ for $j = 1 \cdots m$, where $f_j({\color{darkmagenta}x})$ is a convex function. - A vector $g \equiv \partial f(x_0)$ is called a subgradient of a convex function $f$ at $x_0$ if $f(z) \geq f(x_0) + g^\mathsf{T} (z - x_0)$. - Hence, the cut $(g, \beta)$ is given by $(\partial f(x_0), f(x_0))$ Remarks: - If $f({\color{darkmagenta}x})$ is differentiable, we can simply take $\partial f(x_0) = \nabla f(x_0)$ --- #### 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⁻¹/²ⁿ at each step, ensuring convergence. .mermaid[ <pre> graph TD A[Initial Ellipsoid E₀] --> B[Query Oracle at Center xₖ] B --> C{Feasible?} C -->|Yes| D[Found Solution ✅] C -->|No| E[Get Cutting Plane] E --> F[Update Ellipsoid] F --> G[Volume Shrinks] G --> B style A fill:#ff9800 style B fill:#2196f3 style C fill:#f44336 style D fill:#4caf50 style E fill:#ff9800 style F fill:#9c27b0 style G fill:#2196f3 </pre> ] --- #### Key components of Cutting-plane method - A cutting plane oracle $\Omega$ - A search space $\mathcal{S}$ initially large enough to cover $\mathcal{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 \}$ --- #### Outline of Cutting-plane method - Given initial $\mathcal{S}$ known to contain $\mathcal{K}$. - Repeat - Choose a point $x_0$ in $\mathcal{S}$ - Query the cutting-plane oracle at $x_0$ - If $x_0 \in \mathcal{K}$, quit - Otherwise, update $\mathcal{S}$ to a smaller set that covers: $$\mathcal{S}^+ = \mathcal{S} \cap \{z \mid g^\mathsf{T} (z - x_0) + \beta \le 0\}$$ - If $\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 ``` --- #### From Feasibility to Optimization $$ \begin{array}{ll} \text{minimize} & f_0({\color{darkmagenta}x}), \\ \text{subject to} & {\color{darkmagenta}x} \in \mathcal{K} \end{array} $$ - The optimization problem is treated as a feasibility problem with an additional constraint $f_0({\color{darkmagenta}x}) \le {\color{coral}\gamma}$. - $f_0({\color{darkmagenta}x})$ could be a convex or a *quasiconvex function*. - ${\color{coral}\gamma}$ is also called the *best-so-far* value of $f_0({\color{darkmagenta}x})$. --- #### Convex Optimization Problem - Consider the following general form: $$ \begin{array}{ll} \text{minimize} & {\color{coral}\gamma}, \\ \text{subject to} & \Phi({\color{darkmagenta}x}, {\color{coral}\gamma}) \le 0, \\ & {\color{darkmagenta}x} \in \mathcal{K}, \end{array} $$ where $\mathcal{K}'_{\color{coral}\gamma} = \{ {\color{darkmagenta}x} \mid \Phi({\color{darkmagenta}x}, {\color{coral}\gamma}) \le 0\}$ is the ${\color{coral}\gamma}$-sublevel set of $\{ {\color{darkmagenta}x} \mid f_0({\color{darkmagenta}x}) \le {\color{coral}\gamma}\}$. - 👉 Note: $\mathcal{K'}_{\color{coral}\gamma} \subseteq \mathcal{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}$. --- #### Shrinking - Another possible way is, to update the best-so-far ${\color{coral}\gamma}$ whenever a feasible solution ${\color{darkmagenta}x}'$ is found by solving the equation: $$\Phi({\color{darkmagenta}x}', {\color{coral}\gamma}_\text{new}) = 0 \, .$$ - If the equation is difficuit to solve but ${\color{coral}\gamma}$ is also convex w.r.t. $\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{K}_{\color{coral}\gamma}$. - Repeat - Choose a point $x_0$ in $\mathcal{S}$ - Query the separation oracle at $x_0$ - If $x_0 \in \mathcal{K}_{\color{coral}\gamma}$, update ${\color{coral}\gamma}$ such that $\Phi(x_0, {\color{coral}\gamma}) = 0$. - Update $\mathcal{S}$ to a smaller set that covers: $$\mathcal{S}^+ = \mathcal{S} \cap \{z \mid g^\mathsf{T} (z - x_0) + \beta \le 0\} $$ - If $\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]. $$ \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}$$ --- #### Profit maximization in Convex Form - By taking the logarithm of each variable: - ${\color{darkmagenta}x}_1 = \log y_1$, ${\color{darkmagenta}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{darkmagenta}x}_1} + v_2 e^{ {\color{darkmagenta}x}_2}) - (\alpha {\color{darkmagenta}x}_1 + \beta {\color{darkmagenta}x}_2) \le \log(pA) \\ & {\color{darkmagenta}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) ``` --- #### Area of Applications - 🛡️ Robust convex optimization - oracle technique: affine arithmetic - Semidefinite programming - oracle technique: Cholesky or $LDL^\mathsf{T}$ factorization - Parametric network potential problem - oracle technique: negative cycle detection --- class: middle, center ## 🛡️ Robust Convex Optimization --- ### 🛡️ Robust Optimization Formulation - Consider: $$ \begin{array}{ll} \text{minimize} & \sup_{q \in \mathbb Q} f_0({\color{darkmagenta}x},q), \\ \text{subject to} & f_j({\color{darkmagenta}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} & f_0({\color{darkmagenta}x},q) \leq {\color{coral}\gamma} \\ & f_j({\color{darkmagenta}x},q) \leq 0, \; \forall q \in {\mathbb Q}, \; j = 1,2,\cdots,m. \end{array} $$ --- #### 📚 Example - Profit Maximization Problem (convex) $$ \begin{array}{ll} \text{max} & {\color{coral}\gamma} \\ \text{s.t.} & \log({\color{coral}\gamma} + \hat{v}_1 e^{ {\color{darkmagenta}x}_1} + \hat{v}_2 e^{ {\color{darkmagenta}x}_2}) - (\hat{\alpha} {\color{darkmagenta}x}_1 + \hat{\beta} {\color{darkmagenta}x}_2) \le \log(\hat{p}\,A) \\ & {\color{darkmagenta}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{darkmagenta}x}_1 > 0$, $\alpha = \bar{\alpha} - e_1$, else $\alpha = \bar{\alpha} + e_1$ - if ${\color{darkmagenta}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 $f_j(x_0, q) > 0$ for some $j$ and $q = q_0$, then - the cut $(g, \beta)$ = $(\partial f_j(x_0, q_0), f_j(x_0, q_0))$ - If $f_0(x_0, q) \geq {\color{coral}\gamma}$ for some $q = q_0$, then - the cut $(g, \beta)$ = $(\partial f_0(x_0, q_0), f_0(x_0, q_0) - {\color{coral}\gamma})$ - Otherwise, $x_0$ is feasible, then - Let $q_{\max} = \argmax_{q \in \mathbb Q} f_0(x_0, q)$. - ${\color{coral}\gamma} := f_0(x_0, q_{\max})$. - The cut $(g, \beta)$ = $(\partial f_0(x_0, q_{\max}), 0)$ Remark: for more complicated problems, affine arithmetic could be used [@liu2007robust]. --- class: middle, center ## Matrix Inequalities --- ### Problems With Matrix Inequalities Consider the following problem: $$ \begin{array}{ll} \text{find} & {\color{darkmagenta}x}, \\ \text{subject to} & F({\color{darkmagenta}x}) \succ 0, \end{array} $$ - $F({\color{darkmagenta}x})$: a matrix-valued function - $A \succ 0$ denotes $A$ is positive semidefinite. --- #### 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{darkmagenta}x}, \\ \text{subject to} & v^\mathsf{T} F({\color{darkmagenta}x}) v > 0, \; \forall v \in \mathbb{R}^N - 0^N \end{array}$$ - Consider $v^\mathsf{T} F({\color{darkmagenta}x}) v$ is concave for all $v \in \mathbb{R}^N$ w. r. t. ${\color{darkmagenta}x}$, then the above problem is a convex programming. - Reduce to _semidefinite programming_ if $F({\color{darkmagenta}x})$ is linear w.r.t. ${\color{darkmagenta}x}$, i.e., $F({\color{darkmagenta}x}) = F_0 + {\color{darkmagenta}x}_1 F_1 + \cdots + {\color{darkmagenta}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]. $$ --- #### Naïve implementation - Then, start with $a_{11} = d_{11}$, the basic algorithm of LDLT factorization is: $$ \begin{array}{l} 1~ \text{for}~i=1:n \\ 2~ \quad \text{for}~ j=1:i-1 \\ 3~ \qquad s = a_{ij} - \sum_{k=1}^{j-1} d_{kk} l_{ik} l_{jk} \\ 4~ \qquad l_{ij} = s / d_{jj} \\ 5~ \quad \text{end} \\ 6~ \quad d_{ii} = a_{ii} - \sum_{k=1}^{j-1} d_{kk} l_{ik} l_{jk} \\ 7~ \text{end} \end{array} $$ - Invoke $p^3$ FLOP's, where $p$ is the place the algorithm stops. --- #### Storage representation First, we 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]. $$ --- #### Improved implementation - Then, start with $a_{11} = t_{11}$, the improved 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, a vector $v$ can be constructed to certify that $v^T A v \leq 0$. - Let $L_{1:p}$ denote the partial sub-matrix $L(1:p, 1:p)$ where $p$ is the row of failure. - Then $v = [L_{1: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 $A(1:p, 1:p) \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 $(g, \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{darkmagenta}x}) = A_0 + {\color{darkmagenta}x}_1 A_1 + \cdots + {\color{darkmagenta}x}_n A_n$ - Problem $\min_x \| A({\color{darkmagenta}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{darkmagenta}x}) \\ A^\mathsf{T}({\color{darkmagenta}x}) & {\color{coral}\gamma}\,I \end{array} \right) \succ 0, \end{array} $$ - Binary search on ${\color{coral}\gamma}$ can be used for this problem. --- #### 📚 Example - Estimation of Correlation Function $$ \begin{array}{ll} \min_{ {\color{darkmagenta}\kappa, p} } & \| \Sigma({\color{darkmagenta}p}) + {\color{darkmagenta}\kappa} I - Y \| \\ \text{s. t.} & \Sigma({\color{darkmagenta}p}) \succ 0, {\color{darkmagenta}\kappa} \geq 0 \; .\\ \end{array} $$ - Let $\rho(h) = \sum_i^n {\color{darkmagenta}p}_i \Psi_i(h)$, where - $p_i$'s are the unknown coefficients to be fitted - $\Psi_i$'s are a family of basis functions. - The covariance matrix $\Sigma({\color{darkmagenta}p})$ can be recast as: $$\Sigma({\color{darkmagenta}p}) = {\color{darkmagenta}p}_1 F_1 + \cdots + {\color{darkmagenta}p}_n F_n$$ where $\{F_k\}_{i,j} =\Psi_k( \| s_j - s_i \|_2)$ --- #### 🧪 Experimental Result .pull-left[  : Data Sample (kern=0.5) ] .pull-right[  : Least Square Result ] --- #### 🧪 Experimental Result II .pull-left[  : Data Sample (kern=1.0) ] .pull-right[  : Least Square Result ] --- #### 🧪 Experimental Result III .pull-left[  : Data Sample (kern=2.0) ] .pull-right[  : Least Square Result ] --- class: middle, center ## Multi-parameter Network Problem --- ### 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{darkmagenta}x}, {\color{red}u} \\ \text{subject to} & {\color{red}u_j} - {\color{red}u_i} \le h_{ij}({\color{darkmagenta}x}), \; \forall (i, j) \in {\color{lime}E} , \end{array} $$ - $h_{ij}({\color{darkmagenta}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{darkmagenta}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{darkmagenta}x} \\ \text{subject to} & w_k({\color{darkmagenta}x}) \ge 0, \forall C_k \in \mathcal{C} , \end{array} $$ - $C_k$ is a cycle of $G$ - $w_k({\color{darkmagenta}x}) = \sum_{ (i,j)\in C_k} h_{ij}({\color{darkmagenta}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 $C_k$ under $x_0$, then - the cut $(g, \beta)$ = $(-\partial w_k(x_0), -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{darkmagenta}\pi'} - {\color{darkmagenta}\psi'} \le {\color{coral}\gamma} \\ & {\color{red}u_i'} - {\color{red}u_j'} \le {\color{darkmagenta}\pi'} - a_{ij}', \; \forall a_{ij} \neq 0 \,, \\ & {\color{red}u_j'} - {\color{red}u_i'} \le a_{ij}' - {\color{darkmagenta}\psi'}, \; \forall a_{ij} \neq 0 \,, \\ \text{variables} & {\color{darkmagenta}\pi'}, {\color{darkmagenta}\psi'}, {\color{red}u'} \, . \end{array} $$ where $k'$ denotes $\log( | k | )$ and $x = ({\color{darkmagenta}\pi'}, {\color{darkmagenta}\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{darkmagenta}T_\text{CP} } / {\color{darkmagenta}\beta} \\ \text{subject to} & {\color{red}u_i} - {\color{red}u_j} \le {\color{darkmagenta}T_\text{CP} } - F_{ij}^{-1}({\color{darkmagenta}\beta}), & \forall (i,j) \in {\color{lime}E}_s \,,\\ & {\color{red}u_j} - {\color{red}u_i} \le F_{ij}^{-1}(1 - {\color{darkmagenta}\beta}), & \forall (j,i) \in {\color{lime}E}_h \,, \\ & {\color{darkmagenta}T_\text{CP} } \ge 0, \, 0 \le {\color{darkmagenta}\beta} \le 1 \, , \\ \text{variables} &{\color{darkmagenta}T_\text{CP} }, {\color{darkmagenta}\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{darkmagenta}\beta}$, say ${\color{darkmagenta}\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{darkmagenta}T_\text{CP} } - {\color{darkmagenta}\beta} {\color{coral}\gamma} \le 0\\ & {\color{red}u_i} - {\color{red}u_j} \le {\color{darkmagenta}T_\text{CP} } - F_{ij}^{-1}({\color{darkmagenta}\beta}), & \forall (i,j) \in {\color{lime}E}_s \,,\\ & {\color{red}u_j} - {\color{red}u_i} \le F_{ij}^{-1}(1 - {\color{darkmagenta}\beta}), & \forall (j,i) \in {\color{lime}E}_h \,, \\ & {\color{darkmagenta}T_\text{CP} } \ge 0, \, 0.8 \le {\color{darkmagenta}\beta} \le 1 \, , \\ \text{variables} & {\color{darkmagenta}T_\text{CP} }, {\color{darkmagenta}\beta}, {\color{red}u}. \end{array} $$ --- class: nord-dark, middle, center .column-2.column-norule[ ## Q&A 🎤  ]