凸优化问题笔记（一）最小体积椭球包络问题（MVEE）数学推导

凸优化问题笔记（一）最小体积椭球包络问题（MVEE）数学推导

​ 最小体积椭球包络问题（Minimum Volume Enclosing Ellipsoid, MVEE），又称Löwner-John 椭球问题，是计算几何与凸优化领域的经典问题，目标是在 n 维欧氏空间中，寻找一个包含给定有限点集的体积最小的椭球。以下是完整的数学推导过程。

1.问题定义与椭球参数化

1.1 最小体积椭球包络问题定义

给定一个有\(m\)个点的点云集合\(X=\{x_1,x_2,...,x_m\},x_i\in{\mathbb{R}^n}\), 寻找一个椭球\(\mathcal{\varepsilon}\)，使得满足如下条件：

1. 所有点在椭球内：\(x_i\in{\varepsilon},\forall{i}=1,...,m\);
2. 椭球体积\(vol(\mathcal{E})\)最小；

1.2 椭球的数学表示

椭球有两种常用的参数化表示，对应不同的优化形式：

参数化形式	数学表达式	适用场景
中心 - 形状矩阵形式	\(\mathcal{E}=\{x\in\mathbb{R}^n|(X-c)^T\mathbf{M}^{-1}(X-c)\leq{1}\}\)	一般情况，c 为椭球中心向量，\(\mathbf{M}\in{S_{++}^{n}}\) 为正定形状矩阵
线性变换形式	\(\mathcal{E}=\{x\in\mathbb{R}^{n}|\space||\mathbf{A}x+b||_2\leq{1}\}\)	等价表示，\(\mathbf{A}\in\mathbf{S}_{++}^{n}\)，\(b=-\mathbf{A}c\)

1.3 椭球体积公式

\(n\) 维椭球体积与形状矩阵\(\mathbf{M}\) 或线性变换矩阵\(\mathbf{A}\) 的行列式直接相关：

\[vol(\mathcal{E})=\mathcal{k}_n\cdot\sqrt{det(\mathbf{M})}=\frac{\mathcal{k}_n}{det{(\mathbf{A})}} \tag{1} \]

其中\(k_n=\frac{\pi^{n/2}}{\Gamma(1+n/2)}\) 是\(n\)维单位球体的体积，\(\Gamma(\cdot)\) 为伽马函数。由于\(k_n\) 是常数，最小化体积等价于：

• 最小化\(\sqrt{det{(\mathbf{M})}}\) (中心-形状矩阵的形状)；
• 最大化\(det{(\mathbf{A})}\) (线性变换形状)；

为了凸优化求解，通常对目标函数取对数，转换为：

\[min{-\log(det(\mathbf{A}))} 或\space min{\frac{1}{2}\log(det(\mathbf{M}^{-1}))} \tag{2} \]

对数函数保持单调性，且将乘积优化转化为求和优化，避免数值下溢。

2.凸优化问题转换

2.1 标准形式建立

采用线性变换参数化，\(MVEE\) 问题可转换为以下凸优化问题：

\[\begin{aligned} & argmax_{\mathbf{A}\in{S_{++}^{n}},b\in\mathbb{R}^n} log(det(\mathbf{A})) \\ &s.t. \space ||\mathbf{A}x_i+b||_2\leq{1}, i=1,...,m \end{aligned} \tag{3} \]

等价于（将最大化转为最小化）：

\[\begin{aligned} & argmin_{\mathbf{A}\in{S_{++}^{n}},b\in\mathbb{R}^n} -log(det(\mathbf{A})) \\ &s.t. \space ||\mathbf{A}x_i+b||_2\leq{1}, i=1,...,m \end{aligned} \tag{4} \]

2.2 凸性证明

• 目标函数：\(-log(det{(\mathbf{A})})\) 是凸函数（因为\(log({det(\mathbf{A})})\)是凹函数，负号后变为凸）；
• 约束条件：每个\(||\mathbf{A}x_i+b||_2\leq1\) 是凸约束（2-范数的水平集是凸集）；
• 可行域：凸函数的交际依然是凸集；

因此，\(MVEE\) 是凸优化问题，具有唯一全局最优解。

3.拉格朗日乘数法与 KKT 条件推导

3.1 拉格朗日函数构造

引入非负拉格朗日乘数 \(\lambda_i\geq0\) (对应每个点的约束)，构造拉格朗日函数：

\[\mathcal{L}(\mathbf{A},\mathbf{b},\lambda)=-log(det{(\mathbf{A})})+\sum_{i=1}^{m}\lambda_{i}(||\mathbf{A}x_i+b||_2^2-1) \tag{5} \]

展开其中得二次范数相展开：

\[\mathcal{L}(\mathbf{A},\mathbf{b},\lambda)=-log(det(\mathbf{A}))+\sum_{i=1}^{m}\lambda_i((\mathbf{A}x_i+\mathbf{b})^T(\mathbf{A}x_i+\mathbf{b})-1) \tag{6} \]

3.2 一阶最优化条件(拉格朗日梯度为零)

对线性变换矩阵\(\mathbf{A}\)和\(\mathbf{b}\)求导并令其为零，得到\(KKT\)条件得核心部分。

（1）对中心参数\(\mathbf{b}\) 的梯度

\[\frac{\part{\mathcal{L}}}{\part{\mathbf{b}}}=2\sum_{i=1}^{m}\lambda_i(\mathbf{A}x_i+\mathbf{b})=0 \tag{7} \]

整理得到：

\[\mathbf{b}=-\mathbf{A}\frac{\sum_{i=1}^{m}\lambda_ix_i}{\sum_{i=1}^{m}\lambda_i} \tag{8} \]

这表明最优椭球中心点是给定的加权平均，权重为拉格朗日乘数\(\lambda_i\)。

(2) 对形状矩阵\(\mathbf{A}\) 的梯度：

利用矩阵导数公式\(\frac{\part{log (det(\mathbf{A}))}}{\part{\mathbf{A}}}=\mathbf{A}^{-T}=\mathbf{A}^{-1}\), 对\(\mathbf{A}\) 求导：

\[\frac{\part{\mathcal{L}}}{\part{\mathbf{A}}}=-\mathbf{A}^{-1}+2\sum_{i=1}^m\lambda_i(\mathbf{A}x_i+\mathbf{b})x_i^{T}=0 \tag{9} \]

上式整理得：

\[\mathbf{A}^{-1}=2\sum_{i=1}^{m}\lambda_i(\mathbf{A}x_i+\mathbf{b})x_i^{T} \tag{10} \]

3.3 互补松弛条件

\(\mathbf{KKT}\) 理论要求有效约束满足：

\[\lambda_i(||\mathbf{A}x_i+\mathbf{b}||-1)=0,\space \forall i=1,...,m \tag{11} \]

该式表明：

• 若 \(\lambda_i>0\)，则点\(x_i\)位于椭球边界上(\(||\mathbf{A}x_i+\mathbf{b}||_2=1\))。
• 若\(\lambda_i<0\)， 则点\(x_i\)位于椭球边界上(\(||\mathbf{A}x_i+\mathbf{b}||_2<1\))。

3.4 最优解的关键性质

通过式（8）中的中心条件\(\mathbf{b}=-\mathbf{A}\mathbf{c}\) ,其中 \(\mathbf{c}=-\frac{\sum\lambda_i x_i}{\sum\lambda_i}\) 为椭球中心，代入形状矩阵条件：

\[\mathbf{A}^{-1}=2\sum_{i=1}^{m}\lambda_i\mathbf{A}(x_i-\mathbf{c})x_i^{T} \tag{12} \]

在式（12）中，两边左乘以\(\mathbf{A}\)：

\[\mathbf{I}=2\sum_{i=1}^{m}\lambda_i(x_i-\mathbf{c})x_i^{T}\mathbf{A} \tag{13} \]

再利用\((x_i-\mathbf{c})=-\mathbf{A}^{-1}(Ax_i+\mathbf{b})\)，最终得到简化的最优性条件：

\[\sum_{i=1}^{m}\lambda_i(x_i-\mathbf{c})(x_i-\mathbf{c})^{T}=\frac{1}{2}\mathbf{A}^{-2} \tag{14} \]

4.对偶问题与解析解

4.1 对偶问题推导

​ 通过拉格朗日对偶变换，可将\(MVEE\)问题转换为对偶问题。固定\(\lambda_i\geq{0}\)，对\(\mathbf{A},\mathbf{b}\) 最大化 \(\mathcal{L}(\mathbf{A},\mathbf{b},\lambda)\)，再对\(\lambda\)最小化结果。利用中心条件\(\mathbf{b}=-\mathbf{A}\mathbf{c}\) 代入拉格朗日函数，对偶问题最终为：

\[\begin{aligned} &\min_{\lambda\in{\mathbb{R}^{m}_{+}}} \log\left(\det{\left(\sum_{i=1}^{m}\lambda_i(x_i-\mathbf{c}(\lambda))(x_i-\mathbf{c}(\lambda))^T\right)}\right)+n\log(2)\\ &s.t. \space \sum_{i=1}^{m}\lambda_i=1 \end{aligned} \tag{15} \]

其中\(\mathbf{c}(\lambda)=\sum_{i=1}^{m}\lambda_ix_i\) 是加权中心，约束\(\sum_{i=1}^{m}\lambda_i=1\) 由对偶变量归一化引入。

4.2特殊情况：固定中心的MVEE

当椭球中心固定为原点(\(\mathbf{c}=\mathbf{0}\))，问题化简为如下形式：

\[\begin{aligned} &\min_{\mathbf{A}\in\mathbf{S}^{n}_{++}} -\log(\det(\mathbf{A})) \\ &s.t.||\mathbf{A}x_i||_2\leq 1, i=1,2,3...,m \end{aligned} \tag{16} \]

此问题的对偶问题具有更简洁的形式：

\[\max_{\lambda\in\mathbb{R}_{+}^{m}} \log \det(\sum_{i=1}^m \lambda_ix_ix_i^T)-n\log(\sum_{i=1}^m \lambda_i) \tag{17} \]

最优解：

\[\mathbf{A}^{-2}=\sum_{i=1}^{m}\lambda_i x_i x_i^T \tag{18} \]

这表明固定中心的\(MVEE\) 形状矩阵是给定点的外积的加权和的逆平方根。

5.半定规划（SDP）形式与数值求解

5.1 SDP形式转化

\(MVEE\) 可转换为标准半定规划问题，便于数值求解。引入变量\(\mathbf{X}=\mathbf{A}^{-1}\)（正定矩阵），问题变为：

\[\begin{aligned} &\min_{\mathbf{X}\in{S_{++}^{n}},c\in\mathbb{R}^{n}} \log \det{(X)}\\ &s.t. (x_i-\mathbf{c})^T\mathbf{X}^{-1}(x_i-\mathbf{c})\leq{1} i=1,2,...,m \end{aligned} \tag{19} \]

利用\(Schur\)补性质，约束可改写为线性矩阵不等式\((LMI)\):

\[\left[ \begin{matrix} \mathbf{X}, &x_i-\mathbf{c} \\ (x_i-\mathbf{c})^T,&1 \end{matrix}\right] \mathcal{\succeq} 0 \space \forall{i=1,2,3,...,m} \tag{20} \]

这是凸优化中可高效求解的SDP问题。因此\(MVEE\) 可转换为极小化目标的SDP问题（内点法更擅长极小值化）：

\[\begin{aligned} &\min_{\mathbf{X}\succeq{0},c\in\mathbb{R}^n} \log(\det(\mathbf{X})) \\ & s.t. F_i(\mathbf{X},c)=\left[\begin{matrix}\mathbf{X},&x_i-\mathbf{c}\\ (x_i-\mathbf{c})^T,&1\end{matrix} \right] \succeq0, \space i=1,2,3,...,m \end{aligned} \tag{21} \]

这就是内点法可直接求解的凸SDP问题（目标函数凸，约束是凸锥约束）。

5.2 数值算法—内点法目标函数构建

​ 内点法的本质是用障碍函数“惩罚”约束违法，将带约束优化转换为一系列无约束优化问题，通过逐步减小障碍参数逼近最优解。以下MVCEE凸优化的SDP内点算法推导如下：

Step One：障碍函数的引入

对上述的SDP的锥约束 \(F_i(\mathbf{X},\mathbf{c})\succeq 0\)，引入对数障碍函数：

\[\phi(\mathbf{X},\mathbf{c})=-\sum_{i=1}^{m}\log \det(F_i(\mathbf{X},\mathbf{c})) \tag{22} \]

该函数的性质：

• 当 \(F_i(\mathbf{X},\mathbf{c})\succ{0}\) (严格满足约束)时，\(\phi(\mathbf{X},\mathbf{c})<+ \infty\);
• 当 \(F_i(\mathbf{X},\mathbf{c})\) 趋近于奇异（接近约束边界）时，\(\phi(\mathbf{X},\mathbf{c})\rightarrow +\infty\);
• 当\(\phi(\mathbf{X},\mathbf{c})\) 为凸函数（因 \(\log \det(\cdot)\)是凹函数，符号后凸，凸函数之和依然为凸）。

Step two：障碍问题与中心路径

将障碍函数融入目标函数，得到障碍问题（参数 \(\mu>0\) 为障碍参数）,得到新的目标函数：

\[\min_{\mathbf{X}\succ0,\mathbf{c}} f_{\mu}(\mathbf{X},\mathbf{c})= \log \det(\mathbf{X})-\mu\sum_{i=1}^{m}\log \det{F_i(\mathbf{X},\mathbf{c})} \tag{23} \]

其中：

• \(\mu\)越大，障碍函数的“惩罚越强”，解趋近于靠近可行域中心；
• \(\mu\rightarrow{0}\)时，障碍惩罚消失，解趋近于原SDP的最优解。

对每个\(\mu>0\)，障碍问题的最优解\((\mathbf{X}(\mu),\mathbf{c}(\mu))\) 构成曲线称为中心路径。内点法的核心是“追踪中心路径”，逐步减少\(\mu\) 直到收敛。

5.3 数值算法— 最优性条件推导（梯度+Hessian矩阵）

内点法通过牛顿法求解每个\(\mu\) 对应的障碍问题（无约束凸优化），需要先推导目标函数\(f_{\mu}(\mathbf{X},\mathbf{c})\) 的梯度和Hessian矩阵（牛顿法的核心）。

Step three：符号简化

令 \(d_i=x_i-\mathbf{c}\) 则 \(F_i(\mathbf{X},\mathbf{c})=\left[\begin{matrix}\mathbf{X},d_i \\ d_i^T,1\end{matrix}\right]\)，先给出两个关键矩阵求导公式（对称矩阵求导）：

1. \(\frac{\part{\log\det{\mathbf{A}}}}{\part{\mathbf{A}}}=\mathbf{A}^{-1}\)

2. 对 \(F_i=\left[\begin{matrix}\mathbf{X} ,d_i \\ d_i^T,1 \end{matrix}\right]\), 其逆矩阵可分块表述为：

   \[F_i^{-1}=\left[\begin{matrix}\mathbf{X}^{-1}+\mathbf{X}^{-1}d_i(1-d_i^T\mathbf{X}^{-1}d_i^{T})^{-1}d_i^T\mathbf{X}^{-1},&-\mathbf{X}^{-1}d_i(1-d_i\mathbf{X}^{-1}d_i)^{-1} \\ -(1-d_i^T\mathbf{X}^{-1}d_i)^{-1}d_i^{T}\mathbf{X}^{-1},&(1-d_i^T\mathbf{X}^{-1}d_i)^{-1}\end{matrix} \right] \tag{24} \]

简记为 \(F_i^{-1}=\left[\begin{matrix} Q_i,q_i\\ q_i^T, r_i\end{matrix}\right]\)(\(Q_i\in\mathbb{R}^n,q_i\in\mathbb{R}^n,r_i\in{\mathbb{R}}\))，简化后续推导。

Step Four : 对\(\mathbf{c}\) 的梯度（\(\nabla_{c}f_{\mu}\)）

\[\frac{\part{f_{\mu}}}{\part{\mathbf{c}}}=-\mu\sum_{i=1}^{m}\frac{\part\log\det{F_i}}{\part\mathbf{c}} \tag{25} \]

由矩阵求导规则，可以推导得到惩罚障碍项对\(\mathbf{c}\)梯度，因\(d_i=x_i-\mathbf{c}\)，对\(\mathbf{c}\) 求导等于对\(d_i\)求导取负，因此：

\[\frac{\part{F_{i}}}{\part{\mathbf{c}}}=2q_i \tag{26} \]

由此可得：

\[\nabla_{\mathbf{c}}f_{\mu}=-2\mu\sum_{i=1}^{m}q_i \tag{27} \]

障碍惩罚问题的最优条件要求为0，即为

\[\sum_{i=1}^{m}q_i=0 \tag{28} \]

Step five： 对\(\mathbf{X}\)的梯度\(\nabla_{\mathbf{X}}f_{\mu}\)

\[\frac{\part{f_{\mu}}}{\part{\mathbf{X}}}=\frac{\part{\log\det{X}}}{\part{X}}-\mu\sum_{i=1}^{m}\frac{\part\log\det{F_i}}{\part{\mathbf{X}}}=\mathbf{X}^{-1}-\mu\sum_{i=1}^{m}Q_i \tag{29} \]

最优化条件要求梯度为0，即

\[\mathbf{X}^{-1}-\mu\sum_{i=1}^{m}Q_i=0 \tag{30} \]

公式（28），公式（30）是内点求解障碍惩罚问题的核心最优性条件，需要通过牛顿迭代法迭代满足。

其Hessian矩阵的表达式如下：

5.4 内点法完整算法步骤

内地法分为“中心路径追踪”和"牛顿迭代求解障碍问题"两个核心环节，以下是标准化的算法步骤：

步骤1：初始化

• 初始参数：
  • 设置障碍参数初始值\(\mu_{0}=10\)(可调整，越大初始解越靠近可行域中心)；
  • 设置系数\(\tau=0.1\)（每次迭代缩小障碍参数的比例，通常取\(0.1-0.9\)）;
  • 收敛精度\(\epsilon=10^{-8}\) (终止条件)；
  • 牛顿迭代步长 \(\alpha\in(0,1]\) (线收索确定，保证正定)；
• 初始解：
  • 设置椭球初始中心点 \(c_0=\frac{1}{m}\sum_{i=1}^{m}x_i\) (点集的均值)；
  • 设置形状矩阵\(X_0=\mathbf{I}_n\) (单位矩阵，保证正定)；
  • 验证初始解满足\(F_i(\mathbf{X}_0,\mathbf{c}_0)\succ{0}\)，若不满足则微调（如\(X_0=tI_n,t>1\)）。

步骤2：中心路径追踪（外层循环）

当\(\mu_k>\epsilon\)时，执行以下操作：

1. 求解当前\(\mu\)的障碍问题（内层牛顿迭代）：

   对当前\((X_k,c_k)\)，迭代求解障碍问题\(\min f_{\mu_k}(\mathbf{X},\mathbf{c})\) ，直到梯度满足 \(||\nabla{f_{\mu_k}}||<\epsilon\)；

   • 计算梯度\(\nabla_c{f_{\mu_k}}\) 和 \(\nabla_{\mathbf{X}}f_{\mu_k}\);

   • 计算Hessian矩阵\(H=\left[\begin{matrix}H_{XX},H_{Xc}\\ H_{cX},H_{cc}\end{matrix}\right]\) (Hessian 是梯度的倒数，对称正定)；

   • 计算牛顿方向\(\Delta=-H^{-1}\cdot \nabla{f}_{\mu_k}\)(分解为\(\Delta_{X}\)和\(\Delta_{c}\))

   • 线性搜索确定步长\(\alpha\): 找到最大的\(\alpha\in(0,1]\),使得\(\mathbf{X}_k+\alpha\Delta_{X}\succ{0}\),且

     \[F_i(X_k+\alpha\Delta_{X})\succ{0} \tag{31} \]

   • 跟新解：\(\mathbf{X}_{k+1}=\mathbf{X}_k+\alpha \Delta_{X}\), \(c_{k+1}=c_k+\alpha\Delta_{c}\);

   • 重复直到梯度收敛。

   2.缩小障碍参数：\(\mu_{k+1}=\tau\cdot\mu_{k}\);

   3.更新当前解：\((\mathbf{X}_k,\mathbf{c}_k)=(\mathbf{X}_{k+1},\mathbf{c}_{k+1})\);

步骤3：终止条件

当\(\mu_{k}<\epsilon\) 且梯度 \(||\nabla{f_{\mu_k}}||<\epsilon\) 时，停止迭代，输出最终解：

• 逆形状矩阵 \(\mathbf{P}^{*}=(\mathbf{X}_k)^{-1}\);
• 椭球中心 \(\mathbf{c}^{*}=\mathbf{c}_k\);
• 最终椭球：\((x-c^{*})^{T}P^{*}(x-c^{*})\leq1\);

步骤4：结果转换（可选）

若需要形状矩阵\(M\)(对应之前的\(MVEE\) 定义)，则 \(\mathbf{M}^{*}=(P^{*})^{-1}=\mathbf{X}_k\)，体积为\(Vol(\mathcal{E})=k_n\cdot \sqrt{\det{\mathbf{M}^{*}}}\)。

5.5 关键实现细节

1.矩阵求逆的数值稳定性：

直接求逆\(\mathbf{F}_i^{-1}\) 易数值不稳定，可改用\(Cholesky\) 分解（因\(F_i\succ0\),Cholesky 分解更高效且稳定）。

2.线搜索的必要性：

牛顿方向是"局部最优步长"，需线搜索保证解依然在可行域内（正定），避免迭代发散。

3.Hessian矩阵的简化

高维场景下直接计算Hessian矩阵逆成本高，可用共轭梯度法求解牛顿方向$H\Delta=-\nabla{f} $，无需显式求逆。