多项式近似因式分解的可信验证

The Verification of Approximate Polynomial Factorization

摘要: 众所周知，系数有扰动多项式的因式分解是不连续的。因此，传统的多项式因式分解对于数值计算来说是一个不适定问题。本文利用区间算法，研究多项式近似因式分解的可信计算。给定一个实多项式，本文利用已有算法计算给定多项式其因式分解流形结构，设计算法输出在该因式分解流形结构中一系数为区间的因式分解。算法保证，在该区间因式分解中存在一系数为实数的因式分解，其所对应的多项式为在确定的因式分解流形结构中与给定多项式残差最小的多项式。

Abstract: It is well-known for a polynomial with perturbed coefficients, its factorization is discontinuous. Therefore, the traditional polynomial factorization is an ill-posed problem for numerical computation. This paper is to study the trusted computing of polynomial approximate factorization on the basis of interval algorithm. Given a polynomial of real coefficients, this paper uses the existing algorithm to compute the structure of factorization manifold, and provides a verification algorithm to compute a factorization with interval coefficients in the computed factorization manifold structure. The algorithm is guaranteed that there exists a real factorization within this interval factorization such that the corresponding polynomial of the real factorization is the polynomial with the minimum residual in the computed decomposition manifold structure.

1. 引言

多项式因式分解是多项式基本运算的主要研究问题之一，也是Maple，Mathematica等计算机代数系统的主要功能之一。国内外很多学者对多项式因式分解问题进行了大量的研究。Lenstra和Lovasz [1] 首先提出了关于一元多项式因式分解的时间算法。该算法使用Berlekamp算法并结合Hensel引理对有限域上的多项式进行因式分解。Kaltofen和Von zur Gathen [2] [3] 给出了多元多项式因式分解的时间算法，并对该算法的时间复杂性进行了严格的证明。目前，精确的二元多项式因式分解时间复杂性最低的算法是由Lecerf [4] 提出的。Gao [5] 在Ruppert [6] 针对多项式微分形式相关结论的启发下，提出了一种新的因式分解算法。该算法首先利用Hilbert不可约定理，将多项式由多元降为二元，并提出了二元多项式任意域上的因式分解方法。

众所周知，多项式其系数微小的摄动都有可能改变其因式分解的结构。因此，对于系数有扰动的多项式，其因式分解计算是不连续的。于是当一个多项式系数的精确度有限时，其因式分解计算就显得十分困难。随着对多项式因式分解研究的逐步深入，学者们将研究工作从符号计算扩展到数值计算。Sasaki [7] 通过矩阵运算尽可能多的得到零和关系，提出了一种有效算法改进了多项式的近似因式分解算法。Gao、Kaltofen、May、Yang和Zhi [8] [9] 基于Ruppert和Gao针对多元多项式微分形式的研究成果，利用奇异值分解，结构总体最小二乘以及高斯–牛顿算法来计算多元多项式近似因式分解。Corless、Galligo和Giesbrecht等人 [10] [11] [12] 从几何角度研究了多项式因式分解问题。他们基于参数空间中的投影和随机环上的延拓法，从一个具有扰动的多项式中重建一个近似多项式及其不可约因子，通过建立因式分解的分层复解析流形及其管状邻域，设计了多项式因式分解的数值算法。Kahan [13] 发现了对于不适定代数问题不连续解的流形上隐藏的连续性，解决了对数据微小变化敏感的多项式因式分解的计算问题。Wu和Zeng [14] 提出了基于多项式空间几何和因式分解流形分层的数值分解概念，证明了多项式数值因式分解的存在性、唯一性、李普希茨连续性和收敛性，并提出了多项式数值因式分解算法。该算法消除了传统因式分解的病态性，将数值计算中的不适定因式分解问题完全正则化。

本文利用区间算法，在文献 [14] 工作的基础上设计了在确定因式分解流形结构下，与给定多项式残差向量2-范数最小的多项式的可信计算方法。

2. 预备部分

本文分别用 $\mathbb{N}$ 、 $ℝ$ 和 ${ℝ}^{+}$ 表示非负整数集合、实数集合和正实数集合。用 $\le$ 来表示 ${\mathbb{N}}^n$ 上的乘积顺序，对于 $\alpha ={\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n\right)}^{\text{T}}$， $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_n\right)}^{\text{T}}\in {\mathbb{N}}^n$， $\alpha \le \beta$ 当且仅当 ${\alpha }_i\le {\beta }_i$， $i=1,2,\cdots ,n$。用 ${\prec }_{lex}$ 表示 ${\mathbb{N}}^n$ 上的字典序，若向量 $\beta -\alpha$ 的最左边的非零项是正的，记为 $\alpha {\prec }_{lex}{x}^{\beta }$。假设向量 $x={\left(x_1,x_2,\cdots ,x_r\right)}^{\text{T}}$ 的分

量为变量， $M\left(x\right)$ 是一个 $m\times n$ 维矩阵，该矩阵的每一个分量都是关于变量 $x_1,x_2,\cdots ,x_r$ 的函数， $\frac{\partial M\left(x\right)}{\partial x_k}$ 表示分量为 $\frac{\partial M_{i,j}\left(x\right)}{\partial x_k}$ 的矩阵，其中 $1\le i\le n$， $1\le j\le m$， $1\le k\le r$。 $I_n$ 表示 $n\times n$ 维单位矩阵， $O_{m,n}$ 表示 $m\times n$ 维零矩阵。

令 $ℝ\left[x\right]:=ℝ\left[x_1,x_2,\cdots ,x_r\right]$ 是关于变量 $x_1,x_2,\cdots ,x_r$ 的r元多项式环。令 $x^{\alpha }:=x_1^{{\alpha }_1}{x}_2^{{\alpha }_2}\cdots x_r^{{\alpha }_r}$ 为 $ℝ\left[x\right]$ 中的单项，其中 $\alpha =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_r\right)\in {\mathbb{N}}^n$。多项式 $f\in ℝ\left[x\right]$ 是单项式的有限线性组合，即f可以写成以下形式

$f={\displaystyle \underset{\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_n\right)}^{\text{T}}\in {\mathbb{N}}^n}{\sum }{c}_{\beta }{x}^{\beta }},c_{\beta }\in ℝ$

令 ${\mathrm{deg}}_{x_i}\left(f\right)$ 表示多项式f关于变量 $x_i$ 的次数，则多项式f的次数定义为 $\mathrm{deg}\left(f\right)={\left({\mathrm{deg}}_{x_1}\left(f\right),{\mathrm{deg}}_{x_2}\left(f\right),\cdots ,{\mathrm{deg}}_{x_n}\left(f\right)\right)}^{\text{T}}$。令 $\left[f\right]$ 表示 ${\displaystyle {\prod }_{i=1}^n\left({\mathrm{deg}}_{x_i}\left(f\right)+1\right)}$ 维向量，其分量是由多项式

f的系数按字典序降序排列而成的向量。文献 [8] [9] [15] 引入了多项式w关于次数 $l={\left(l_1,l_2,\cdots ,l_n\right)}^{\text{T}}$ 的卷积矩阵 $C_l\left(w\right)$，对于一个 $\mathrm{deg}\left(w\right)=l$ 的多项式w，矩阵 $C_l\left(w\right)$ 乘以 $\left[f\right]$ 得到的向量为多项式w乘以多项式f所得的多项式其系数按照字典序排列而成的向量。

令 $\mathbb{I}ℝ$ 表示所有区间 $\left\{\left[\underset{\_}{x},\bar{x}\right]:\underset{\_}{x},\bar{x}\in ℝ,\underset{\_}{x}\le \bar{x}\right\}$ 构成的集合。区间向量 $X={\left(X_1,X_2,\cdots ,X_n\right)}^{\text{T}}$ 是一个分量为区间的n维向量组，且满足条件

$X=\left\{\left(x_1,x_2,\cdots ,x_n\right)\in {ℝ}^n:x_i\in X_i,1\le i\le n\right\}.$

区间矩阵 $A\in \mathbb{I}{ℝ}^{m\times n}$ 是具有区间项的矩阵，如果区间矩阵A中的任意实矩阵是满秩的，则称区间矩阵A满秩。如果属于区间矩阵 $M\in \mathbb{I}{ℝ}^{m\times n}$ 的任意实矩阵M，对所有的 $1\le i\le n$， $1\le j\le n$ 满足条件 $M_{i,j}=M_{j,i}$，则称区间矩阵M是区间对称矩阵。如果区间对称矩阵M中的每一个对称矩阵都为正定矩阵，则称区间对称矩阵M为区间对称正定矩阵。对于一个区间对称矩阵M，INTLAB工具箱 [16] 中的isspd函数可以验证M的正定性，即命令isspd (M)返回值为1，说明区间对称矩阵M为对称正定矩阵。

定理1 (见 [17])：设 $f:{ℝ}^n\to {ℝ}^n$， $f={\left(f_1,\cdots ,f_n\right)}^{\text{T}}$， $f_1,\cdots ,f_n$ 是连续可微函数。对于 $\tilde{x}\in {ℝ}^n$ 和区间向量 $X\in \mathbb{I}{ℝ}^n$，其中 $0\in X$，令 $f^{\prime }\left(\tilde{x}+X\right)$ 表示 $f$ 在区间 $\tilde{x}+X$ 处的区间雅可比矩阵。给定 $R\in {ℝ}^{n\times n}$ 和满足条件 $f^{\prime }\left(\tilde{x}+X\right)\subseteq M$ 的区间矩阵 $M\in \mathbb{I}{ℝ}^{m\times n}$，若

$-Rf\left(\tilde{x}\right)+\left(I_n-RM\right)X\subseteq \mathrm{int}\left(X\right),$

则存在唯一的 $\hat{x}\in \tilde{x}+X$，使得 $f\left(\hat{x}\right)=0$，其中 $\mathrm{int}\left(X\right)$ 表示 $X$ 的内部。

3. 多项式近似因式分解

假设 $\gamma p_1^{k_1}{p}_2^{k_2}\cdots p_r^{k_r}$ 是多项式 $f\in ℝ\left[x\right]$ 的近似因式分解，其中 $\gamma \in ℝ$， $k_1\le k_2\le \cdots \le k_r$， $k_i\in \mathbb{N}$， $p_1,p_2,\cdots ,p_r$ 是成对互质的无平方多项式，且对于任意的 $1\le i\le r$， $\mathrm{deg}\left(p_i\right)=m_i={\left(m_{i,1},m_{i,2},\cdots ,m_{i,n}\right)}^{\text{T}}$。Wu和Zeng在文献 [14] 中提出了基于多项式空间几何和分解流形分层的数值多项式分解的概念，并且定义多项式全体近似因式分解的子集

$F_{m_1^{k_1}{m}_2^{k_2}\cdots m_r^{k_r}}=\left\{\gamma p_1^{k_1}{p}_2^{k_2}\cdots p_r^{k_r}:\gamma \in ℝ,\mathrm{deg}{p}_i\le m_i,i=1,\cdots ,r;p_1,p_2,\cdots ,p_r是成对互质的无平方多项式的\right\}.$ (1)

Wu和Zeng称子集 $F_{m_1^{k_1}{m}_2^{k_2}\cdots m_r^{k_r}}$ 中任意一个多项式的因式分解流形结构为 $m_1^{k_1}{m}_2^{k_2}\cdots m_r^{k_r}$，并引入了以下向量

$b_i={\left(b_{i,\beta }|\beta \in {\mathbb{N}}^n,\beta \le m_i\right)}^{\text{T}},1\le i\le r,$ (2)

其中向量 $b_i$ 的分量是按照关于 $\beta$ 的字典序降序排列。然后，Wu和Zeng将求解因式分解流形结构 $m_1^{k_1}{m}_2^{k_2}\cdots m_r^{k_r}$ 上的近似因式分解问题转换成计算下述超定方程组的最小二乘解

$\varphi \left(\gamma ,\left[p_1\right],\cdots ,\left[p_r\right]\right)=\left(\begin{array}{c}\left[\gamma p_1^{k_1}\cdots p_r^{k_r}-f\right]\\ b_1^{\text{T}}\left[p_1\right]-1\\ ⋮\\ b_r^{\text{T}}\left[p_r\right]-1\end{array}\right)=0.$ (3)

4. 主要结果

4.1. 数值部分

数值部分利用文献 [8] 中的数值多项式因式分解算法来计算因式分解流形结构 $m_1^{k_1}{m}_2^{k_2}\cdots m_r^{k_r}$。然后，求超定方程组(3)的最小二乘解，得到该因式分解流形结构上高精度的近似解，记作

$f\approx \bar{\gamma }{\bar{p}}_1^{k_1}{\bar{p}}_2^{k_2}\cdots {\bar{p}}_r^{k_r}.$

我们引入扰动向量

$e_i={\left(e_{i,\beta }|\beta \in {\mathbb{N}}^n,\beta \le m_i\right)}^{\text{T}},1\le i\le r,$ (4)

其中向量 $e_i$ 的分量关于 $\beta$ 的字典序降序排列。定义扰动多项式 ${\bar{p}}_i\left(e_i\right)$ 为

${\bar{p}}_i\left(e_i\right)={\bar{p}}_i+{\displaystyle \underset{\beta \in {\mathbb{N}}^n,\beta \le m_i}{\sum }{e}_{i,\beta }{x}^{\beta }},1\le i\le r.$

为了简化表达式，令 $e={\left(e_1,e_2,\cdots ,e_r\right)}^{\text{T}}$，定义非线性函数

$F\left(e\right)={\Vert \left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\Vert }_2^2,$ (5)

定义

${\bar{q}}_i\left(e\right)=\left(\frac{k_i}{{\bar{p}}_i\left(e_i\right)}\right)\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r},1\le i\le r.$

非线性函数 $F\left(e\right)$ 的梯度为

$F^{\prime }\left(e\right)=\{\begin{array}{l}{C}_{m_1}^{\text{T}}\left({\bar{q}}_1\left(e\right)\right)\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\hfill \\ C_{m_2}^{\text{T}}\left({\bar{q}}_2\left(e\right)\right)\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\hfill \\ ⋮\hfill \\ C_{m_r}^{\text{T}}\left({\bar{q}}_r\left(e\right)\right)\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\hfill \end{array},$ (6)

其中 $C_{m_i}\left({\bar{q}}_i\left(e\right)\right)$ 是多项式 ${\bar{q}}_i\left(e\right)$ 关于次数 $m_i={\left(m_{i,1},m_{i,2},\cdots ,m_{i,n}\right)}^{\text{T}}$ 的卷积矩阵。关于变量 $e_{i,\beta }$，有 $F^{\prime }\left(e\right)$ 对 $e_{i,\beta }$ 的偏导数如下

$\begin{array}{c}\frac{\partial F^{\prime }\left(e\right)}{\partial e_{i,\beta }}\text{=}\left(\begin{array}{c}\frac{\partial C_{m_1}^{\text{T}}\left({\bar{q}}_1\left(e\right)\right)}{\partial e_{i,\beta }}\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\\ \frac{\partial C_{m_2}^{\text{T}}\left({\bar{q}}_2\left(e\right)\right)}{\partial e_{i,\beta }}\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\\ ⋮\\ \frac{\partial C_{m_r}^{\text{T}}\left({\bar{q}}_r\left(e\right)\right)}{\partial e_{i,\beta }}\left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]\end{array}\right)\\ +\left(\begin{array}{c}{C}_{m_1}^{\text{T}}\left(q_1\left(e\right)\right)\\ C_{m_2}^{\text{T}}\left(q_2\left(e\right)\right)\\ ⋮\\ C_{m_r}^{\text{T}}\left(q_r\left(e\right)\right)\end{array}\right)\frac{\partial \left[\gamma {\bar{p}}_1{\left(e_1\right)}^{k_1}{\bar{p}}_2{\left(e_2\right)}^{k_2}\cdots {\bar{p}}_r{\left(e_r\right)}^{k_r}-f\right]}{\partial e_{i,\beta }}.\end{array}$ (7)

注1：数值部分的主要工作是计算下列优化问题的稳定点

$\begin{array}{l}\min F\left(e\right)\\ \text{s}\text{.t}\text{.}\{\begin{array}{l}{b}_1^{\text{T}}\left[{\bar{p}}_1\left(e_1\right)\right]-1=0\hfill \\ b_2^{\text{T}}\left[{\bar{p}}_2\left(e_2\right)\right]-1=0\hfill \\ ⋮\hfill \\ b_r^{\text{T}}\left[{\bar{p}}_r\left(e_r\right)\right]-1=0\hfill \end{array}.\end{array}$ (8)

¨

通过引入拉格朗日算子 $\lambda ={\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_r\right)}^{\text{T}}$，定义拉格朗日函数

$G\left(e,\lambda \right)=F\left(e\right)+{\displaystyle \underset{i=1}{\overset{r}{\sum }}{\lambda }_i}\left(b_i^{\text{T}}\left[{\bar{p}}_i\left(e_i\right)\right]-1\right),$ (9)

则有

$G^{\prime }\left(e,\lambda \right)=\left[\begin{array}{c}{F}^{\prime }\left(e\right)\\ \begin{array}{c}{b}_1^{\text{T}}\left[{\bar{p}}_1\left(e_1\right)\right]-1\\ b_2^{\text{T}}\left[{\bar{p}}_2\left(e_2\right)\right]-1\\ ⋮\\ b_r^{\text{T}}\left[{\bar{p}}_r\left(e_r\right)\right]-1\end{array}\end{array}\right]+\left[\begin{array}{c}{\lambda }_1{b}_1\\ {\lambda }_2{b}_2\\ \begin{array}{c}⋮\\ {\lambda }_r{b}_r\end{array}\\ \begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\end{array}\right].$ (10)

由于 $G\left(e,\lambda \right)$ 的Hession矩阵不包含变量 ${\lambda }_i$， $1\le i\le r$，为了方便起见将 $G^{″}\left(e,\lambda \right)$ 写作 $G^{″}\left(e\right)$，显然，

$G^{″}\left(e\right)=\left[\begin{array}{cccccccc}& & & & b_1& & & \\ & F^{″}\left(e\right)& & & & b_2& & \\ & & & & & & \ddots & \\ & & & & & & & b_r\\ b_1^{\text{T}}& & & & & & & \\ & b_2^{\text{T}}& & & & & & \\ & & \ddots & & & & O_{r\times r}& \\ & & & b_r^{\text{T}}& & & & \end{array}\right].$ (11)

注2：通过选择零向量作为初始值，数值部分采用数值牛顿迭代法

$\left(\begin{array}{c}\bar{e}\\ \bar{\lambda }\end{array}\right)\leftarrow \left(\begin{array}{c}\bar{e}\\ \bar{\lambda }\end{array}\right)-G^{″}{\left(\bar{e}\right)}^{-1}{G}^{\prime }\left(\bar{e},\bar{\lambda }\right)$ (12)

计算非线性系统 $G^{\prime }\left(e,\lambda \right)=0$ 的近似解 $\left(\bar{e},\bar{\lambda }\right)$。 ¨

4.2. 验证部分

验证部分分为两个阶段，第一阶段验证非线性系统 $G^{\prime }\left(e,\lambda \right)=0$ 近似解的可信误差界，第二阶段验证优化问题(8)的稳定点是否为其局部最优解。

注 3：对于满足 $\bar{e}\in \tilde{E}$ 的区间向量 $\tilde{E}$，使用区间运算和式(6) (7) (11)计算区间雅克比矩阵 $G^{″}\left(\tilde{E}\right)$。利用定理1和verifynlss函数计算区间向量 $\hat{E}$ 和 $\hat{\Lambda }$，其满足条件存在实向量 $\hat{e}\in \hat{E}$ 和实向量 $\hat{\lambda }\in \hat{\Lambda }$，使得 $G^{\prime }\left(\hat{e},\hat{\lambda }\right)=0$。 ¨

对每个 $1\le i\le r$，令 $s_i$ 表示 $e_i$ 的维数，令 $s={\displaystyle \underset{i=1}{\overset{r}{\sum }}{s}_i}$。记 $E={\left(E_1,E_2,\cdots ,E_r\right)}^{\text{T}}$，其中 $E_i\in \mathbb{I}{ℝ}^{s_t}$。向量 $E_i$， $b_i$ 和 $e_i$ 的最后一项分别用 $E_{i,s}$， $b_{i,s}$ 和表示，令，和分别表示向量，和删除最后一项所得的向量。

不失一般性，我们假设

对每一个，存在一个关于变量的线性函数，其满足条件

显然，

定义非线性函数

则有约束优化问题(8)可以转化为下面无约束优化问题

(13)

验证部分第二阶段使用isspd函数验证区间Hession矩阵的正定性。当isspd函数返回值为1时，满足条件的实向量为问题(8)的局部最优解。

5. 主要算法

算法1

输入

输出，，，或“Failure”。

步骤1 数值部分

步骤1.1 计算多项式f的因式分解流形结构；

步骤1.2 计算系统(3)的近似解；

步骤1.3 计算系统的近似解。

步骤2 验证部分

步骤2.1 初始化，，和；

步骤2.2 当进行以下操作：

令；

；

计算；

令；

若，则令

如果，则返回“Failure”并停止。

步骤2.3 计算区间Hession阵，若isspd返回值为1，则输出，和，否则，返回“Failure”。 ¨

注4：在算法1中，给定实数x，则intval(x)输出区间。给定实数a，b，其中，则输出实区间。给定，则输出包含和零向量的最小的区间向量。 ¨

通过以上分析可以得到下述定理。

定理2：假设算法1成功输出因式分解流形结构，系数向量和摄动区间向量，则存在一个实向量，它是优化问题(8)的局部最优解。

6. 算例

例1 考虑多项式，其中多项式，，。摄动多项式r与f的次数相同且r的系数是[−5, 5]上随机的一个数。应用算法1可得

因此根据定理2，存在一个实向量，它是优化问题(8)的局部最优解。

例2 [14] 考虑多项式，其中多项式，，。应用算法1可得

因此根据定理2，存在一个实向量，它是优化问题(8)的局部最优解。

例3考虑单变量多项式，其中多项式，，。摄动多项式r与f的次数相同且r的系数是[−5, 5]上随机的一个数。应用算法1可得

因此根据定理2，存在一个实向量，它是优化问题(8)的局部最优解。

基金项目

吉林省自然科学基金(批准号：11601039)。

NOTES

^*通讯作者。