layout: true class: typo, typo-selection --- count: false class: nord-dark, middle, center # 👁️ Ellipsoid Method Revisited @luk036 👨‍💻 · 2026 📅 --- ### 🥚 Basic Ellipsoid Method - An ellipsoid ${\color{lime}\mathcal{E}}({\color{green}x_c}, {\color{coral}P})$ is specified as a set $$\{{\color{green}x} \mid ({\color{green}x}-{\color{green}x_c}){\color{coral}P}^{-1}({\color{green}x}-{\color{green}x_c}) \leq 1 \},$$ where ${\color{green}x_c}$ is the center of the ellipsoid.  --- ### 🔄 Ellipsoid Method Steps .mermaid[ <pre> flowchart TD style A fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px style B fill:#e3f2fd,stroke:#1565c0,stroke-width:3px style C fill:#fff3e0,stroke:#e65100,stroke-width:3px style D fill:#f3e5f5,stroke:#7b1fa2,stroke-width:3px style E fill:#c8e6c9,stroke:#2e7d32,stroke-width:3px style F fill:#ffcdd2,stroke:#c62828,stroke-width:3px style G fill:#e3f2fd,stroke:#1565c0,stroke-width:3px style H fill:#fff3e0,stroke:#e65100,stroke-width:3px A[Initial ellipsoid E₀] --> B[Check center xc] B --> C[Query oracle at xc] C --> D{xc feasible?} D -->|Yes| E[Return xc ✓] D -->|No| F[Compute cutting plane g,h] F --> G[New minimum-volume ellipsoid] G --> H{Ellipsoid small?} H -->|Yes| E H -->|No| B </pre> ] --- ### 🐍 Python code {.allowframebreaks} ```python import numpy as np class ell: def __init__(self, val, x): '''ell = { x | (x - xc)' * P^-1 * (x - xc) <= 1 }''' n = len(x) if np.isscalar(val): self.P = val * np.identity(n) else: self.P = np.diag(val) self.xc = np.array(x) self.c1 = float(n*n)/(n*n-1.0) def update_core(self, calc_ell, g, beta):... def calc_cc(self, g):... def calc_dc(self, g, beta):... def calc_ll(self, g, beta):... ``` --- ### 🔪 Updating the ellipsoid (deep-cut) - Calculation of minimum volume ellipsoid covering: $${\color{lime}\mathcal{E}} \cap \{{\color{firebrick}z} \mid {\color{coral}g}^\mathsf{T} ({\color{firebrick}z} - {\color{green}x_c}) + {\color{coral}h} \leq 0\}$$ - Let ${\color{coral}\tilde{g}} = {\color{coral}P}\,{\color{coral}g}$, ${\color{coral}\tau} = \sqrt{{\color{coral}g}^\mathsf{T}{\color{coral}\tilde{g}} }$, ${\color{coral}\alpha} = {\color{coral}h}/{\color{coral}\tau}$. - If ${\color{coral}\alpha} > 1$, intersection is empty. - If ${\color{coral}\alpha} < -1/{\color{coral}n}$ (shallow cut), no smaller ellipsoid can be found. - Otherwise, $$\begin{array}{lll} {\color{green}x_c}^+ &=& {\color{green}x_c} - \frac{1+{\color{coral}n}{\color{coral}\alpha}}{({\color{coral}n}+1){\color{coral}\tau}} {\color{coral}\tilde{g}} \\ {\color{coral}P}^+ &=& \frac{{\color{coral}n}^2(1-{\color{coral}\alpha}^2)}{{\color{coral}n}^2 - 1}\left({\color{coral}P} - \frac{2{\color{coral}\rho}}{(1+{\color{coral}\alpha}){\color{coral}\tau}^2} {\color{coral}\tilde{g}}{\color{coral}\tilde{g}}^\mathsf{T}\right) \end{array}$$ --- ### 🐍 Python code (updating) ```python def update_core(self, calc_ell, g, beta): Pg = self.P.dot(g) tsq = g.dot(Pg) tau = np.sqrt(tsq) alpha = beta/tau status, rho, sigma, delta = calc_ell(alpha) if status != 0: return status, tau self.xc -= (rho/tau) * Pg self.P -= (sigma/tsq) * np.outer(Pg, Pg) self.P *= delta return status, tau ``` --- ### 🐍 Python code (deep cut) ```python def calc_dc(self, alpha): '''deep cut''' if alpha == 0.0: return self.calc_cc() n = len(self.xc) status, rho, sigma, delta = 0, 0.0, 0.0, 0.0 if alpha > 1.: status = 1 # no sol'n elif n*alpha < -1.: status = 3 # no effect else: rho = (1.0+n*alpha)/(n+1) sigma = 2.0*rho/(1.0+alpha) delta = self.c1*(1.0-alpha*alpha) return status, rho, sigma, delta ``` --- ### ⏸️ Parallel Cuts - Oracle returns a pair of cuts instead of just one. - The pair of cuts is given by ${\color{coral}g}$ and $({\color{coral}h}_1, {\color{coral}h}_2)$ such that: $$\begin{array}{l} {\color{coral}g}^\mathsf{T} ({\color{green}x} - {\color{green}x_c}) + {\color{coral}h}_1 \leq 0, \\ {\color{coral}g}^\mathsf{T} ({\color{green}x} - {\color{green}x_c}) + {\color{coral}h}_2 \geq 0, \end{array}$$ for all ${\color{green}x} \in \mathcal{K}$. - Only linear inequality constraint can produce such parallel cut: $$ {\color{coral}l} \leq {\color{coral}a}^\mathsf{T} {\color{green}x} + {\color{coral}b} \leq {\color{coral}u}, \qquad {\color{coral}L} \preceq F({\color{green}x}) \preceq {\color{coral}U} $$ - Usually provide faster convergence. --- ### ⏸️ Parallel Cuts  --- ### 📐 Updating the ellipsoid - Let ${\color{coral}\tilde{g}} = {\color{coral}P}\,{\color{coral}g}$, ${\color{coral}\tau} = \sqrt{{\color{coral}g}^\mathsf{T}{\color{coral}\tilde{g}} }$, ${\color{coral}\alpha}_1 = {\color{coral}h}_1/{\color{coral}\tau}$, ${\color{coral}\alpha}_2 = {\color{coral}h}_2/{\color{coral}\tau}$. - If ${\color{coral}\alpha}_2 > 1$, it reduces to deep-cut with ${\color{coral}\alpha} = {\color{coral}\alpha}_1$. - If ${\color{coral}\alpha}_1 > {\color{coral}\alpha}_2$, intersection is empty. - If ${\color{coral}\alpha}_1 {\color{coral}\alpha}_2 < -1/{\color{coral}n}$, no smaller ellipsoid can be found. Otherwise, $${\color{green}x_c}^+ = {\color{green}x_c} - \frac{{\color{coral}\rho}}{{\color{coral}\tau}} {\color{coral}\tilde{g}}, \qquad {\color{coral}P}^+ = {\color{coral}\delta}\left({\color{coral}P} - \frac{{\color{coral}\sigma}}{{\color{coral}\tau}^2} {\color{coral}\tilde{g}}{\color{coral}\tilde{g}}^\mathsf{T}\right)$$ where $$\begin{array}{lll} {\color{coral}\xi} &=& \sqrt{4(1 - {\color{coral}\alpha}_1^2)(1 - {\color{coral}\alpha}_2^2) + {\color{coral}n}^2({\color{coral}\alpha}_2^2 - {\color{coral}\alpha}_1^2)^2}, \\ {\color{coral}\sigma} &=& \frac{1}{{\color{coral}n}+1}({\color{coral}n} + \frac{2}{({\color{coral}\alpha}_1 + {\color{coral}\alpha}_2)^2}(1 - {\color{coral}\alpha}_1{\color{coral}\alpha}_2 - \frac{{\color{coral}\xi}}{2})), \\ {\color{coral}\rho} &=& \frac{1}{2}({\color{coral}\alpha}_1 + {\color{coral}\alpha}_2) {\color{coral}\sigma} \\ {\color{coral}\delta} &=& \frac{{\color{coral}n}^2}{{\color{coral}n}^2-1} (1 - \frac{1}{2}({\color{coral}\alpha}_1^2 + {\color{coral}\alpha}_2^2 - \frac{{\color{coral}\xi}}{{\color{coral}n}})) \end{array}$$ --- ### 🐍 Python code (parallel cut) ```python def calc_ll(self, alpha): '''parallel or deep cut''' if np.isscalar(alpha): return self.calc_dc(alpha) # parallel cut a0, a1 = alpha if a1 >= 1.0: return self.calc_dc(a0) n = len(self.xc) status, rho, sigma, delta = 0, 0.0, 0.0, 0.0 aprod = a0 * a1 if a0 > a1: status = 1 # no sol'n elif n*aprod < -1.0: status = 3 # no effect else: asq = alpha * alpha asum = a0 + a1 asqdiff = asq[1] - asq[0] xi = np.sqrt(4.0*(1.0-asq[0])*(1.0-asq[1]) + n*n*asqdiff*asqdiff) sigma = (n + (2.0*(1.0 + aprod - xi/2.0)/(asum*asum)))/(n+1) rho = asum * sigma/2.0 delta = self.c1*(1.0 - (asq[0] + asq[1] - xi/n)/2.0) return status, rho, sigma, delta ``` --- ### 📚 Example: FIR filter design  - The time response is: $${\color{firebrick}y}[{\color{coral}t}] = \sum_{{\color{coral}k}=0}^{{\color{coral}n}-1}{{\color{coral}h}[{\color{coral}k}]{\color{coral}u}[{\color{coral}t}-{\color{coral}k}]}$$ --- ### 📚 Example: FIR filter design (cont'd) - The frequency response: $$H({\color{coral}\omega}) = \sum_{{\color{coral}m}=0}^{{\color{coral}n}-1}{{\color{coral}h}({\color{coral}m})e^{-j{\color{coral}m}{\color{coral}\omega}} }$$ - The magnitude constraints on frequency domain are expressed as $${\color{coral}L}({\color{coral}\omega}) \leq |H({\color{coral}\omega})| \leq {\color{coral}U}({\color{coral}\omega}), \forall {\color{coral}\omega}\in(-\infty,+\infty)$$ where ${\color{coral}L}(\omega)$ and ${\color{coral}U}(\omega)$ are the lower and upper (nonnegative) bounds at frequency ${\color{coral}\omega}$ respectively. - The constraint is non-convex in general. --- ### 📚 Example: FIR filter design (cont'd) - However, via *spectral factorization*, it can transform into a convex one: $${\color{coral}L}^2({\color{coral}\omega}) \leq R({\color{coral}\omega}) \leq {\color{coral}U}^2({\color{coral}\omega}), \forall {\color{coral}\omega}\in(0,\pi)$$ where - $R({\color{coral}\omega})=\sum_{{\color{coral}i}=-1+{\color{coral}n}}^{{\color{coral}n}-1}{{\color{green}r}({\color{coral}t})e^{-j{\color{coral}\omega}{\color{coral}t}} }=|H({\color{coral}\omega})|^2$ - ${\color{green}\mathbf{r}}=({\color{green}r}(-{\color{coral}n}+1),{\color{green}r}(-{\color{coral}n}+2),...,{\color{green}r}({\color{coral}n}-1))$ are the autocorrelation coefficients. --- ### 📚 Example: FIR filter design (cont'd) - ${\color{green}\mathbf{r}}$ can be determined by ${\color{coral}\mathbf{h}}$: $${\color{green}r}({\color{coral}t}) = \sum_{{\color{coral}i}=-{\color{coral}n}+1}^{{\color{coral}n}-1}{{\color{coral}h}({\color{coral}i}){\color{coral}h}({\color{coral}i}+{\color{coral}t})}, {\color{coral}t}\in\mathbf{Z}.$$ where ${\color{coral}h}({\color{coral}t})=0$ for ${\color{coral}t} < 0$ or ${\color{coral}t} > {\color{coral}n}-1$. - The whole problem can be formulated as: $$\begin{array}{ll} \text{min} & {\color{coral}\gamma} \\ \text{s.t.} & {\color{coral}L}^2({\color{coral}\omega}) \leq R({\color{coral}\omega}) \leq {\color{coral}U}^2({\color{coral}\omega}), \; \forall {\color{coral}\omega} \in [0,\pi] \\ & R({\color{coral}\omega}) > 0, \forall {\color{coral}\omega} \in [0,\pi] \end{array}$$ --- ### 📚 Example: Maximum Likelihood estimation $$\begin{array}{ll} \min_{{\color{coral}\kappa}, {\color{coral}p}} & \log \det (\Omega({\color{coral}p}) + {\color{coral}\kappa} I) + \mathrm{Tr}((\Omega({\color{coral}p}) + {\color{coral}\kappa} I)^{-1}{\color{coral}Y}) \\ \text{s.t.} & \Omega({\color{coral}p}) \succeq 0, {\color{coral}\kappa} \geq 0 \\ \end{array}$$ 👉 Note: 1st term is concave, 2nd term is convex - However, if there is enough samples such that ${\color{coral}Y}$ is a positive definite matrix, then the function is convex within $[0, 2{\color{coral}Y}]$ --- ### 📚 Example: Maximum Likelihood estimation (cont'd) - Therefore, the following problem is convex: $$\begin{array}{ll} \min_{{\color{coral}\kappa}, {\color{coral}p}} & \log \det V({\color{coral}p}) + \mathrm{Tr}(V({\color{coral}p})^{-1}{\color{coral}Y}) \\ \text{s.t.} & \Omega({\color{coral}p}) + {\color{coral}\kappa} I = V({\color{coral}p}) \\ & 0 \preceq V({\color{coral}p}) \preceq 2{\color{coral}Y}, {\color{coral}\kappa} \geq 0 \end{array}$$ --- class: nord-light, 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? --- ### 🧩 Mixed-Integer Convex Programming Consider: $$\begin{array}{ll} \text{minimize} & f_0({\color{green}x}), \\ \text{subject to} & f_j({\color{green}x}) \leq 0, \; \forall j=1,2,\ldots \\ & {\color{green}x} \in \mathbb{D} \end{array}$$ where - $f_0({\color{green}x})$ and $f_j({\color{green}x})$ are "convex" - Some design variables are discrete. --- ### 🔮 Oracle Requirement - The oracle looks for the nearby discrete solution ${\color{green}x_d}$ of ${\color{green}x_c}$ with the cutting-plane: $${\color{coral}g}^\mathsf{T} ({\color{green}x} - {\color{green}x_d}) + {\color{coral}h} \leq 0, {\color{coral}h} \geq 0, {\color{coral}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) --- ### 📚 Example: FIR filter design  --- count: false class: nord-dark, middle, center ## Q&A 🎤