为了更好的数学公式和代码块体验，建议在 plweb 阅读本文． ### 背景 回顾<discussion=68a1a46d3734f520583f97cb>《基于半正定规划的多项式平方和分拆》</discussion>（以下称“前文”），我们将多项式的 SOS 转换为含参矩阵半正定： $$\boldsymbol Q(\boldsymbol y)=\boldsymbol Q_0+\sum_{i=1}^ty_i\boldsymbol Q_i\succeq0.$$ 为了找出符合条件的 $\boldsymbol y$，我们使用半正定规划（SDP）求解器计算 $$\min\ f(y_1,\dots,y_n)\quad\text{s.t.}\ \boldsymbol Q_0+\sum_{i=1}^ty_i\boldsymbol Q_i\succeq0.$$ 这里 $f$ 可以取任何函数，因为我们只需要一组可行解，例如前文中取 $\operatorname{tr}\boldsymbol Q$，也有很多时候直接取常函数． 但是这个基础的算法有明显的缺陷，例如前文中如果改目标函数为常函数， ``` sdpsol=SemidefiniteOptimization[0,VectorGreaterEqual[{newmat,0},{"SemidefiniteCone",Length[monos]}],matvars]; (* 输出 {x[2,2]->1.99999,x[3,2]->2.50607*10^-6,x[3,3]->1.71319,x[4,2]->-0.856597,...} *) ``` 这里并非每个值都接近一个简洁的有理数，如果随意近似（或直接代入浮点值）到原矩阵，求特征值得 ``` Eigenvalues[newmat/.sdpsol] (* 输出 {10.0128,10.0128,8.27956,1.68632,1.68631,1.50618,0.000401756,0.000401755,2.31893*10^-7,7.84312*10^-8,7.72458*10^-8,1.99659*10^-8,-1.56603*10^-8,-1.23339*10^-8,-1.19924*10^-8} *) ``` 有若干负值，因此得到的不是真正的半正定矩阵．本文将介绍一种优化的算法，提高找出精确有理解的成功率． ### 问题 首先分析出现小负特征值的原因．注意矩阵 $\boldsymbol Q$ 满足 $\boldsymbol x^\top\boldsymbol Q\boldsymbol x=f(a,b,c)$，其中 $\boldsymbol x$ 是幂次减半的所有单项式构成的列向量、$f(a,b,c)$ 是待证不等式．当取等条件存在，即 $f(a_0,b_0,c_0)=0$ 能成立时，相应就有 $\boldsymbol x_0^\top\boldsymbol Q\boldsymbol x_0=0$．根据定义，$\boldsymbol Q$ 不可能是正定矩阵，所以无论如何取 $\boldsymbol y$，至多能让 $\boldsymbol Q$ 为半正定矩阵． 利用线性代数中半正定锥（Positive semidefinite cone）的概念，一个形象的理解方式是， - 所有的 $\boldsymbol y$ 位于一个较高维的空间， - 使得 $\boldsymbol Q(\boldsymbol y)$ 半正定的 $\boldsymbol y$ 只构成空间中的一个较低维的面 $S$． SDP 求解器能给出这个面上一个坐标的值，然而一旦我们进行截断或简单有理近似，这个点就会不受控制地向某个方向发生微小偏移，进而脱离 $S$． 为了解决问题，我们应当主动寻找半正定矩阵 $\boldsymbol Q$．一种基本的方法是——基于取等条件列方程求解． ### 方法 熟知以下结论， >若半正定矩阵 $\boldsymbol M$ 和 向量 $\boldsymbol v$ 满足 $\boldsymbol v^\top\boldsymbol M\boldsymbol v=0$，则 $\boldsymbol M\boldsymbol v=\boldsymbol0$． 取 $\boldsymbol M=\boldsymbol Q$，$\boldsymbol v$ 为单项式向量 $\boldsymbol x$ 代入取等条件的值，则有关于 $y_1,\dots,y_t$ 的线性方程组 $\boldsymbol Q(\boldsymbol y)\boldsymbol x_{\text{eq}}=0$．若取等条件足够多则该方法能直接解出 $\boldsymbol y$，最后验证 $\boldsymbol Q(\boldsymbol y)\succeq0$ 即可． **例题**（AoPS）证明：对任意实数 $a$，$b$，$c$ 有 $$16 a^4-48 a^3 b-48 a^3 c+61 a^2 b^2-12 a^2 b c+\\168 a^2 c^2+12 a b^3-96 a b^2 c+72 a b c^2-144 a c^3\\{}+4 b^4-24 b^3 c+57 b^2 c^2-36 b c^3+36 c^4\ge0.$$ 取 $\boldsymbol x=\left(a^2,b^2,c^2,a b,b c,a c\right)^\top$，用前文的方法得 Gram 矩阵 $\boldsymbol Q=$ $$ \begin{pmatrix} 16 & m_{1,2} & m_{1,3} & -24 & m_{1,5} & -24 \\ m_{1,2} & 4 & m_{2,3} & 6 & -12 & m_{2,6} \\ m_{1,3} & m_{2,3} & 36 & m_{3,4} & -18 & -72 \\ -24 & 6 & m_{3,4} & 61-2 m_{1,2} & -m_{2,6}-48 & -m_{1,5}-6 \\ m_{1,5} & -12 & -18 & -m_{2,6}-48 & 57-2 m_{2,3} & 36-m_{3,4} \\ -24 & m_{2,6} & -72 & -m_{1,5}-6 & 36-m_{3,4} & 168-2 m_{1,3} \\ \end{pmatrix} $$ 这个矩阵记为 `newM`．求原不等式的取等条件，可以利用齐次性，不妨设 $a=1$，运行 ``` sols=Solve@Reduce[bds==0&&a==1,Reals]; (*{{a->1,b->1/2,c->1/4},{a->1,b->1/2,c->2},{a->1,b->4,c->2}}*) ``` 解方程 ``` params=Solve@Thread[Expand@Flatten[x.newM/.sols]==0] (*代码解析：x.newM/.sols 会将 sols 每组解分别代入，返回多重列表； 用 Thread 构建列表每项均为 0 的方程组*) ``` 结果是 $$\begin{array}c\displaystyle\left\{m_{1,2}\to -\frac{16}{7},m_{1,3}\to \frac{48}{7},m_{1,5}\to \frac{120}{7},\right.\\\displaystyle\left.m_{2,3}\to \frac{24}{7},m_{2,6}\to -\frac{12}{7},m_{3,4}\to \frac{144}{7}\right\}.\end{array}$$ 最后代入得 $$\boldsymbol Q= \begin{pmatrix} 16 & -\frac{16}{7} & \frac{48}{7} & -24 & \frac{120}{7} & -24 \\ -\frac{16}{7} & 4 & \frac{24}{7} & 6 & -12 & -\frac{12}{7} \\ \frac{48}{7} & \frac{24}{7} & 36 & \frac{144}{7} & -18 & -72 \\ -24 & 6 & \frac{144}{7} & \frac{459}{7} & -\frac{324}{7} & -\frac{162}{7} \\ \frac{120}{7} & -12 & -18 & -\frac{324}{7} & \frac{351}{7} & \frac{108}{7} \\ -24 & -\frac{12}{7} & -72 & -\frac{162}{7} & \frac{108}{7} & \frac{1080}{7} \\ \end{pmatrix}. $$ 它的特征值有三个**精确**的 $0$，另外三个约等于 $203$，$108$，$14.5$，所以 $\boldsymbol Q$ 半正定．作 LDL 分解，得 $\boldsymbol Q = \boldsymbol L\boldsymbol D\boldsymbol L^\top$，其中 $$\begin{aligned}\boldsymbol L &= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{7} & 1 & 0 & 0 & 0 & 0 \\ \frac{3}{7} & \frac{6}{5} & 1 & 0 & 0 & 0 \\ -\frac{3}{2} & \frac{7}{10} & 1 & 1 & 0 & 0 \\ \frac{15}{14} & -\frac{13}{5} & -\frac{1}{2} & 0 & 1 & 0 \\ -\frac{3}{2} & -\frac{7}{5} & -2 & 0 & 0 & 1 \end{pmatrix},\\ \boldsymbol D &= \text{diag}\left(16, \frac{180}{49}, \frac{972}{35}, 0, 0, 0\right).\end{aligned}$$ 所以 $$\boldsymbol L\sqrt{\boldsymbol D}=\begin{pmatrix} 4 & 0 & 0 \\ -\frac{4}{7} & \frac{6\sqrt{5}}{7} & 0 \\ \frac{12}{7} & \frac{36\sqrt{5}}{35} & \frac{18\sqrt{105}}{35} \\ -6 & \frac{3\sqrt{5}}{5} & \frac{18\sqrt{105}}{35} \\ \frac{30}{7} & -\frac{78\sqrt{5}}{35} & -\frac{9\sqrt{105}}{35} \\ -6 & -\frac{6\sqrt{5}}{5} & -\frac{36\sqrt{105}}{35} \end{pmatrix}.$$ 最终结果是 $f_1^2+f_2^2+f_3^2$，其中 $$\begin{aligned}f_1^2 &= \frac{243}{35} (2 a-c)^2 (b-2 c)^2, \\[4pt] f_2^2 &=\frac{9}{245} (b-2 c)^2 (7 a+10 b-6 c)^2, \\[4pt] f_3^2 &= \frac{4}{49} \left(14 a^2-21 a b-21 a c-2 b^2+15 b c+6 c^2\right)^2. \end{aligned}$$