A New Deflation Method For Verifying the Isolated Singular Zeros of Polynomial Systems

TL;DR

This paper introduces a novel deflation method for refining and verifying isolated singular zeros in polynomial systems, utilizing linear combinations and derivatives to improve accuracy and efficiency, especially for high multiplicity zeros.

Contribution

It presents the first deflation approach based on linear combinations of polynomials, with strategies to reduce system size and enhance numerical stability.

Findings

Effective for high multiplicity zeros in large systems

Works well for non-polynomial systems with isolated singular zeros

Provides strategies to accelerate the deflation process

Abstract

In this paper, we develop a new deflation technique for refining or verifying the isolated singular zeros of polynomial systems. Starting from a polynomial system with an isolated singular zero, by computing the derivatives of the input polynomials directly or the linear combinations of the related polynomials, we construct a new system, which can be used to refine or verify the isolated singular zero of the input system. In order to preserve the accuracy in numerical computation as much as possible, new variables are introduced to represent the coefficients of the linear combinations of the related polynomials. To our knowledge, it is the first time that considering the deflation problem of polynomial systems from the perspective of the linear combination. Some acceleration strategies are proposed to reduce the scale of the final system. We also give some further analysis of the tolerances we use, which can help us have a better understanding of our method.The experiments show that our method is effective and efficient. Especially, it works well for zeros with high multiplicities of large systems. It also works for isolated singular zeros of non-polynomial systems.

Tables

${𝐩}_i$	x	y
${\widetilde{𝐩}}_1$	$0.000006851$	$-0.000004368$
${\widetilde{𝐩}}_2$	$0.0000000000000$	$-0.0000021841948$
${\widetilde{𝐩}}_3$	$-0.0000000000000$	$-0.0000010920974$
${\widetilde{𝐩}}_4$	$-0.0000000000000$	$-0.0000005460487$
${\widetilde{𝐩}}_5$	$-0.0000000000000$	$-0.0000002730243$
${\widetilde{𝐩}}_6$	$-0.0000000000000$	$-0.0000001365122$
${\widetilde{𝐩}}_7$	$-0.0000000000000$	$-0.0000000682561$
${\widetilde{𝐩}}_8$	$0.0000000000000$	$-0.0000000341280$
${\widetilde{𝐩}}_9$	$-0.0000000000000$	$-0.0000000170640$
${\widetilde{𝐩}}_{10}$	$-0.0000000000000$	$-0.0000000853201$

${𝐩}_i$	x	y	$\alpha$
${\widetilde{𝐩}}_1$	0.000006851	$-0.000004368$	$-0.091484814324$
${\widetilde{𝐩}}_2$	0.0000002081968	0.0000001993959	$-0.0500009466172$
${\widetilde{𝐩}}_3$	0.0000000003977	$-0.0000000000002$	$-0.0500000007560$
${\widetilde{𝐩}}_4$	0.0000000000000	0.0000000000000	$-0.0500000000000$

Table 3. Table 1: Comparison of VDSS and viss for simple systems

System	$var$	$mul$	Verified acc		Max err	times		Final size
			VDSS	viss		VDSS	viss	VDSS	viss
DZ1	4	131	true	e-322	0	0.3066	0.3337	4	16
DZ2	3	16	e-14	e-14	0	0.2989	0.7343	3	24
DZ3	2	4	e-14	e-15	e-14	0.8780	1.0093	3	10
cbms1	3	11	true	e-322	0	0.1851	0.1107	3	6
cbms2	3	8	true	e-322	0	0.2546	0.1271	3	6
mth191	3	4	e-14	e-14	e-32	0.3118	0.1221	4	6
KSS	10	638	e-14	e-14	0	8.2295	0.3036	19	20
RuGr09	2	4	e-323	e-14	0	0.1567	0.4955	2	8
LZ	100	3	e-320	e-14	0	2.0197	13.3068	100	300
Ojika1	2	3	e-14	e-14	0	0.7636	0.3447	5	6
Ojika2	3	2	e-14	e-14	e-16	0.3936	0.2942	5	6
Ojika3	3	2	e-14	e-14	0	0.3967	0.3427	4	6
Ojika4	3	3	e-14	e-14	0	0.1851	1.0621	3	9
Decker2	3	4	e-323	e-14	0	0.1752	0.4650	3	8
Caprasse	4	4	e-14	e-14	e-31	2.0180	0.5126	6	8
Cyclic9	9	4	e-14	e-14	e-15	5.9266	3.6878	12	18

Table 4. Table 2: Comparison of VDSS and viss for large systems

System	$var$	$mul$	Verified acc		Max err	times		Final size
			VDSS	viss		VDSS	viss	VDSS	viss
LZ2000	2000	3	e-319	–	0	448.07	–	2000	–
simple1	5	9	e-14	e-14	0	0.29	8.20	5	45
addvar2	4	12	e-14	e-13	e-14	11.10	250.67	6	32
reduce3	4	24	e-14	e-14	e-11	12.21	317.50	7	12
unre3	4	36	e-15	e-14	e-13	4.08	360.32	4	32
unre5	8	576	e-14	e-14	e-13	24.26	229.83	8	64
big1	20	512	e-14	e-15	e-12	29.92	1724.09	20	160
big2	20	8192	e-14	e-14	e-12	40.90	1751.61	20	160
big3	30	196608	e-15	e-14	e-15	155.18	425.51	30	240
rankone2	6	32	e-15	e-15	e-15	6.8693	1.5199	11	12
rankone3	6	96	e-15	e-14	e-14	12.44	136.54	11	48
breadth2	5	$2^5$	e-322	–	0	0.20	–	5	–
large3	100	$3^{100}$	e-323	e-319	0	187.88	647.86	100	400
large6	500	$4^{100}$	e-321	e-34	0	905.00	3262.78	500	2000
large8	500	$4^{300}$	e-321	–	0	1745.85	–	500	–
nonpoly3	3	64	e-322	e-14	0	0.19	6.62	3	36

Table 5. Table 3: Comparison of VDSS and other methods for four examples

	Method A		Method B		Method C		Method D		Method E
	$Poly$	$Var$	$Poly$	$Var$	$Poly$	$Var$	$Poly$	$Var$	$Poly$	$Var$
1	16	4	22	4	22	4	16	4	8	4
2	24	11	11	3	12	3	12	3	5	3
3	32	17	6	2	6	2	6	2	4	3
4	96	41	54	3	54	3	22	3	5	3

Equations

$${\mathbf{d}}_{\mathbf{x}}^{\mathbf{\gamma }}(f)=\frac{1}{{\gamma }_1!\cdots {\gamma }_n!}\cdot \frac{{\partial }^{|\mathbf{\gamma }|}f}{\partial x_1^{{\gamma }_1}\cdots \partial x_n^{{\gamma }_n}},\forall f\in \mathbb{C}[\mathbf{x}],$$

$${\bm{\mathrm{J}}}({\mathbf{F}})=\left(\begin{array}[]{ccc}\frac{\partial f_{1}}{\partial x_{1}}&\ldots&\frac{\partial f_{1}}{\partial x_{n}}\\ \vdots&\ddots&\vdots\\ \frac{\partial f_{n}}{\partial x_{1}}&\ldots&\frac{\partial f_{n}}{\partial x_{n}}\\ \end{array}\right).$$

$$\exists \epsilon >0:\{\mathbf{y}\in {\mathbb{C}}^n:\Vert \mathbf{y}-\mathbf{p}\Vert <\epsilon \}\cap {\mathbf{F}}^{-1}(0)=\{\mathbf{p}\},$$

$$\mathrm{rank}(\mathrm{J}(\mathbf{F})(\mathbf{p}))<n.$$

$$f(\mathbf{x})=f(\mathbf{p})+\sum _{j=1}^n\frac{\partial f(\mathbf{p})}{\partial x_j}(x_j-p_j)+\sum _{1\le i,j\le n}\frac{{\partial }^{2}f(\mathbf{p})}{\partial x_i\partial x_j}(x_i-p_i)(x_j-p_j)+\dots .$$

$$\frac{\partial f(\mathbf{p})}{\partial x_j}=0,\forall 1\le j\le n.$$

$$\left|\frac{\partial f(\tilde{\mathbf{p}})}{\partial x_j}\right|<\theta ,\forall 1\le j\le n.$$

$$f=\sum _{d=1}^{deg(f)}{f}_d.$$

$$f(\mathbf{x})=x_1+3x_3+4x_4-x_1^2+x_3^2-x_4^2-x_2^3.$$

$$|f(\mathbf{p})|=0,|{\mathrm{J}}_2(f)(\mathbf{p})|=0,|{\mathrm{J}}_i(f)(\mathbf{p})|\ne 0,i=1,3,4,$$

$$\begin{array}{rl}f(\mathbf{x})=& 0.003002001+0.998(x_1-0.001)-3\cdot 10^{-6}(x_2+0.001)+3.004(x_3-0.002)\\ & +4.002(x_4+0.001)-(x_1-0.001{)}^2+3\cdot 10^{-3}(x_2+0.001{)}^2+(x_3-0.002{)}^2\\ & -(x_4+0.001{)}^2-(x_2+0.001{)}^3.\end{array}$$

$$g=\sum _{i=1}^n\sum _j{a}_{i,j}\frac{{\partial }^{|{\mathbf{\gamma }}_{i,j}|}{f}_i}{\partial {\mathbf{x}}^{{\mathbf{\gamma }}_{i,j}}},$$

$$\mathrm{rank}(\mathrm{J}(g_1,\dots ,g_s)(\mathbf{p}))=s.$$

$$\mathrm{rank}(\mathrm{J}(g_1,\dots ,g_s,g_{s+1})(\mathbf{p}))=s+1.$$

$$\dim \mathbb{V}(g_1,\dots ,g_s,f_{s+1})>\dim \mathbb{V}(f_1,\dots ,f_s,f_{s+1}).$$

$$\mathrm{rank}(\mathrm{J}(f_1,\dots ,f_s)(\mathbf{p}))=s.$$

$$f_i=\sum _{k=1}^n{a}_{ik}{x}_k+{\mathrm{T}}_i,$$

$$L:\left\{\begin{array}[]{l}y_{i}=\sum\limits_{k=1}^{n}a_{ik}x_{k},1\leq i\leq s\\ y_{i}=x_{i},\ i=s+1,\ldots,n.\end{array}\right.\ \$$

$$\{\begin{array}{rl}{F}_i& =y_i+L^{-1}({\mathrm{T}}_i),i=1,\dots ,s,\\ F_{s+i}& =\sum _{j=1}^s{b}_{i,j}{y}_j+L^{-1}({\mathrm{T}}_{s+i}),i=1,\dots ,k.\end{array}$$

$$\dim \mathbb{V}(F_1,\dots ,F_{s+k})\le n-s-1.$$

$$\begin{array}{rl}{F}_{s+1}^{\prime }& =\frac{{\partial }^{d-1}{F}_l}{\partial y_{s+1}^{d_{s+1}-1}{y}_{s+2}^{d_{s+2}}\cdots y_n^{d_n}}=\sum _{i=1}^n{\gamma }_i{y}_i+{\mathrm{T}}_l^{\prime },d=\sum _{j=s+1}^n{d}_j.\end{array}$$

$$\mathrm{rank}(\mathrm{J}(F_1,\dots ,F_s,F_{s+1}^{\prime })({\mathbf{p}}^{\prime }))=s+1.$$

$$f_i^{\prime }=L(F_i)=f_i(i=1,\dots ,s),f_{s+1}^{\prime }=L(F_{s+1}^{\prime })\text{with}\mathrm{rank}(\mathrm{J}(f_1^{\prime },\dots ,f_s^{\prime },f_{s+1}^{\prime })(\mathbf{p}))=s+1.$$

$$\mathrm{rank}(\mathrm{J}(f_1,\dots ,f_r)(\mathbf{p}))=\mathrm{rank}(\mathrm{J}(\mathbf{F})(\mathbf{p}))=r,0\le r\le n.$$

$${\mathbf{F}}_1=\{f_1^{(1)},\dots ,f_r^{(1)},f_{r+1}^{(1)}\},$$

$${\mathbf{F}}_1(\mathbf{p})=0\text{and}\mathrm{rank}(\mathrm{J}(f_1^{(1)},\dots ,f_r^{(1)},f_{r+1}^{(1)})(\mathbf{p}))=r+1.$$

$${\mathbf{F}}_2=\{f_1^{(2)},\dots ,f_{r+1}^{(2)},f_{r+2}^{(2)}\},$$

$${\mathbf{F}}_2(\mathbf{p})=0\text{and}\mathrm{rank}(\mathrm{J}(f_1^{(2)},\dots ,f_{r+1}^{(2)},f_{r+2}^{(2)})(\mathbf{p}))=r+2.$$

$${\mathbf{F}}_{n-r}=\{f_1^{(n-r)},f_2^{(n-r)},\dots ,f_n^{(n-r)}\},$$

$${\mathbf{F}}_{n-r}(\mathbf{p})=0\text{and}\mathrm{rank}(\mathrm{J}(f_1^{(n-r)},f_2^{(n-r)},\dots ,f_n^{(n-r)})(\mathbf{p}))=n.$$

Full text

1 Introduction

Solving polynomial systems with singular zeros is always a challenge in algebraic and geometric computation. For an isolated simple zero of a polynomial system, the classical Newton’s method is widely used and quadratic convergent. However, for singular zeros of a polynomial system, Newton’s method is not fit for the original system directly because it converges slowly or even doesn’t converge in a bad situation. What’s more, it is an ill-posed problem to compute an isolated singular zero of a polynomial system or a nonlinear system, since a small perturbation of coefficients may transform an isolated singular zero into a cluster of simple zeros.

Therefore, finding methods to keep the quadratic convergence of Newton’s method for singular zeros is a way to handle this problem. Given a polynomial system with an isolated singular zero, we can construct a new system owing the same singular zero as an isolated simple one. Based on this idea, in recent years, there are many symbolic or symbolic-numerical methods coming up to deal with this problem. The basic idea is the deflation techniques [1, 3, 4, 6, 7, 8, 10, 25], which usually have two basic strategies: adding new equations only or both new equations and new variables to the original system.

Deflation for an isolated singular solution originated from the ideas of Ojika [19, 20, 21]. T. Ojika et al. present a deflation algorithm for determining the multiple zeros for a system of nonlinear equations. Through triangulating the Jacobian matrix of the original system at an approximate zero, new equations, which comes from the minors of the Jacobian matrix, are introduced to the original system to reduce the multiplicity until they get a system which is regular at the singular zero.

In [5], Giusti and Yakoubsohn propose a construction, which is based on two operations: deflating and kerneling, to determine a regular system without adding new variables. In the deflating, all the partial derivatives of the polynomials, which are zero at the multiple zero, are introduced to replace the corresponding polynomials. The kerneling operation consists of adding the polynomials given by the nonzero numerators of the coefficients of the Schur complement of the Jacobian matrix of the original system to the original system.

In [8], Hauenstein and Wampler define a strong deflation by only adding new equations coming from the one order differential of the Jacobian matrix of the original system to the original system. Different from [5], at each iteration step, both the number and the degree of the added equations are reduced.

In [9], Mourrain et al. give a method which uses a single linear differential form defined from the Jacobian matrix of the input system, and defines the deflated system by applying this differential form to the original system.

These above methods do introduce new equations and finally get a new system owing the isolated singular zero of the original system as a simple zero. In order to get the new polynomials, one needs to compute the determinant of some polynomial matrices. Thus the degree of the polynomials in the new system may be very high.

In the following, denote $n$ as the number of both the variables and the equations in the original system, and $\mu$ as the multiplicity of the isolated singular zero of the original system.

In [18], Yamamoto introduces new equations and new variables to the original system simultaneously. New variables are used to bring some perturbations of the original system and the Jacobian matrix of the original system, which produce new equations.

In [12, 13, 24], Leykin et al. present an effective modification of Newton’s method to restore quadratic convergence for isolated singular solutions of polynomial systems. Different from [18], new variables are only introduced to the Jacobian matrix of the original system, which produce new equations. Meanwhile, they also prove that the number of deflation stages is bounded by $\mu$.

In [2], Dayton and Zeng modify the method in [12] and further prove that the number of deflation steps is bounded by the depth of the dual space. For the special case of breathe one, they also propose a modified deflation method, which is based on duality analysis, to reduce the final size $2^{\mu-1}n\times 2^{\mu-1}n$ of deflated system in [12] to $\mu n\times\mu n$.

In [23], by introducing a smoothing parameter to the original system and $n-1$ new variables to the Jacobian matrix of the original system, which produces new equations, Rump and Graillat consider the case of the double zero of the original system. In [15], based on the parameterized multiplicity structure, Li and Zhi generalize the algorithm in [23] to deflate the breath-one isolated singular zero of the original system. Their final deflated regular system is of size $\mu n\times\mu n$.

In [17], based on the given multiplicity structure of the original system, which depends on the accuracy of the given approximate multiple zero, Mantzaflaris and Mourrain give a method to find a (small) perturbed system of the original system and then first compute a deflated system in one deflation step. The size of the final deflated system is equal to $\mu n\times\mu n$.

In [16], by lifting the independent perturbations in the first-order differential system appearing in [18] back to the original system, Li and Zhi modify the method in [18] and also prove that the modified deflation technique terminates after a finite number of steps bounded by the depth of the dual space. The size of the final modified regularized system is bounded by $2^{\mu-1}n\times 2^{\mu-1}n$.

In [9], by introducing some variables to represent the coefficients of the dual basis, Mourrain et al. give a method to deflate the original system and determine the multiplicity structure simultaneously. They also show that the number of variables and equations in this method is bounded by $n+n\mu(\mu-1)/2$ and $n\mu+n(n-1)(\mu-1)(\mu-2)/4$. However, one point worth noting is that this method needs to know the monomial basis of the original system first.

These methods introduce new variables and new equations to the original system simultaneously. By repeatedly using these deflation constructions, they will get an augmented system finally, which has an isolated simple zero, whose partial projection corresponds to the isolated singular zero of the original system.

Main contribution. In this paper, given a polynomial system ${\mathbf{F}}\subset{\mathbb{C}}[\mathbf{x}]$ with an isolated singular zero $\mathbf{p}$, by computing the derivatives of the input polynomials directly or the linear combinations of the related polynomials, we propose a new deflation method to construct a final deflated system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})$, which has an isolated simple zero $(\mathbf{p},\hat{\bm{\alpha}})$, whose projection corresponds to the isolated singular zero $\mathbf{p}$ of the input system. New variables $\bm{\alpha}$ are introduced to represent the coefficients of the linear combinations of the related polynomials to ensure the accuracy of the numerical implementation. Moreover, we also prove that the size of our deflation system depends on the depth or the multiplicity of $\mathbf{p}$.

Compared to the previous methods, our method has the following differences:

For the input system ${\mathbf{F}}$, we can, if needed, compute the derivatives of every $f_{i}$ to get the needed polynomials, which are regular at $\mathbf{p}$ at the beginning. Then, we put all these polynomials together to construct a system ${\mathbf{F}}_{0}$, such that the rank $r$ of its Jacobian matrix at $\mathbf{p}$ is maximal. In some cases, we have $r=n$, which means that we need not introduce new variables. 2. 2.

We compute the derivatives of the linear combinations of the related polynomials to get some polynomials which are regular at $\mathbf{p}$. Here we introduce new variables to represent the coefficients of the linear combination. 3. 3.

Considering that we know only the approximate zero $\tilde{\mathbf{p}}$ in actual computations, we use a tolerance $\theta$ to judge if a polynomial is $\theta$-regular or $\theta$-singular at $\tilde{\mathbf{p}}$ and another tolerance $\varepsilon$ to judge the numerical rank of the Jacobian matrix. As long as the tolerance $\theta$ is chosen properly, we will get the same judgement in numerical case as in the exact case. Thus, our deflation system usually has the same isolated zero as the input system. Inspired by the work [12] of Leykin et al., we also give some further analysis on the tolerances $\theta$ and $\varepsilon$, which tells us that our final system is a perturbed system with a bounded perturbation in the worst case. To make our final system as accurate as possible, we also analyse the case that the tolerance $\theta$ is not introduced.

Thanks to the above acceleration strategies, the size of the final system in our actual computations is much less than that we give in theory. Furthermore, we implement our method in Matlab. The experiments show that our method is effective and efficient, especially for large systems with singular zeros of high multiplicities. Besides, for the non-polynomial systems, our method is also applicable.

The paper is organized as below. We introduce some notations and preliminaries in the next section. In Section 3, we give a new deflation idea to construct a deflated square system from the input system with an isolated singular zero. An effective version of our method is given in Section 4. Some numerical experiment results are given to demonstrate the performance of our algorithm in Section 5 and at last, we draw a conclusion in Section 6.

2 Notations and Preliminaries

Let $\mathbb{C}$ be the complex field and $\mathbb{C}[\mathbf{x}]=\mathbb{C}[x_{1},\ldots,x_{n}]$ be the polynomial ring. Denote ${\mathbf{F}}=\{f_{1},f_{2},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ as a polynomial system and $\deg(f_{i})$ as the degree of the polynomial $f_{i}$. Similarly, $\deg({\mathbf{F}})=\max\limits_{f_{i}\in{\mathbf{F}}}{\deg(f_{i})}$. Let $\mathbf{p}=(p_{1},\ldots,p_{n})\in\mathbb{C}^{n}$. ${\mathbf{F}}(\mathbf{p})={\mathbf{0}}$ denotes that $\mathbf{p}$ is a zero of ${\mathbf{F}}(\mathbf{x})=\mathbf{0}$.

Let ${\mathbb{V}}({\mathbf{F}})\subset{\mathbb{C}}^{n}$ denote the variety defined by ${\mathbf{F}}$ and $\dim{\mathbb{V}}({\mathbf{F}})$ denote the dimension of ${\mathbb{V}}({\mathbf{F}})$.

Let $\mathbf{d}^{\bm{\gamma}}_{\mathbf{x}}:\mathbb{C}[\mathbf{x}]\rightarrow\mathbb{C}[\mathbf{x}]$ denote the differential functional defined by

[TABLE]

where ${\bm{\gamma}}=(\gamma_{1},\ldots,\gamma_{n})\in\mathbb{N}^{n}$ with $\mathbb{N}=\{0,1,2,\ldots\}$ and $|\bm{\gamma}|=\sum\limits_{i=1}^{n}\gamma_{i}$.

Denote ${\mathrm{rank}}(A)$ as the rank of a matrix $A$. Denote ${\bm{\mathrm{J}}}({\mathbf{F}})$ as the Jacobian matrix of ${\mathbf{F}}$. That is,

[TABLE]

For a polynomial $f\in\mathbb{C}[\mathbf{x}]$, let ${\bm{\mathrm{J}}}(f)$ denote $(\frac{\partial f}{\partial x_{1}},\frac{\partial f}{\partial x_{2}},\ldots,\frac{\partial f}{\partial x_{n}})$ and ${\bm{\mathrm{J}}}_{i}(f)=\frac{\partial f}{\partial x_{i}}$. Let $\bm{\mathrm{J}}({\mathbf{F}})(\mathbf{p})$ denote the value of a function matrix ${\bm{\mathrm{J}}}({\mathbf{F}})$ at a point $\mathbf{p}$, similarly for ${\bm{\mathrm{J}}}(f)(\mathbf{p})$.

Definition 1**.**

An isolated zero of ${\mathbf{F}}(\mathbf{x})={\mathbf{0}}$ is a point $\mathbf{p}\in\mathbb{C}^{n}$ which satisfies:

[TABLE]

where ${\mathbf{F}}^{-1}(\mathbf{0})\triangleq\{\mathbf{p}\in{\mathbb{C}}^{n}:{\mathbf{F}}(\mathbf{p})={\mathbf{0}}\}$.

Definition 2**.**

We call an isolated zero $\mathbf{p}\in{\mathbb{C}}^{n}$ of ${\mathbf{F}}(\mathbf{x})={\mathbf{0}}$ an isolated singular zero if and only if

[TABLE]

Otherwise, $\mathbf{p}$ is an isolated regular( simple) zero of ${\mathbf{F}}(\mathbf{x})={\mathbf{0}}$.

The Taylor series expansion (Taylor expansion for short) of $f\in{\mathbb{C}}[\mathbf{x}]$ at $\mathbf{p}=(p_{1},\ldots,p_{n})\in{\mathbb{C}}^{n}$ is

[TABLE]

Definition 3**.**

Let ${\mathbf{p}}\in{\mathbb{C}}^{n}$ and $f(\mathbf{p})=0$. We say $f\in{\mathbb{C}}[\mathbf{x}]$ is singular at $\mathbf{p}$ if

[TABLE]

Otherwise, we say $f$ is regular at $\mathbf{p}$.

Definition 4**.**

Let $f\in{\mathbb{C}}[\mathbf{x}]$, $\tilde{\mathbf{p}}\in{\mathbb{C}}^{n}$ and a tolerance $\theta>0$, s.t. $|f(\tilde{\mathbf{p}})|<\theta$. We say $f$ is $\theta$-singular* at $\tilde{\mathbf{p}}$ if*

[TABLE]

Otherwise, we say $f$ is $\theta$-regular* at $\tilde{\mathbf{p}}$.*

Lemma 5**.**

*Let $f\in\mathbb{C}[\mathbf{x}]\setminus{\mathbb{C}}$, s.t. $f(\mathbf{p})=0$. Then there exists at least a $\bm{\gamma}\in\mathbb{N}^{n}$, s.t. $\mathbf{d}^{\bm{\gamma}}_{\mathbf{x}}(f)$ is regular at $\mathbf{p}$. *

Proof.

Without loss of generality, we assume $\mathbf{p}={\mathbf{0}}$. Then $f$ can be rewritten as a sum of homogeneous polynomials as

[TABLE]

Since $f\not\equiv 0$, there exists at least a $\bm{\gamma}^{\prime}\in{\mathbb{N}}^{n}$ such that $\mathbf{d}^{\bm{\gamma}^{\prime}}_{\mathbf{x}}(f)(\mathbf{p})\neq 0$. Thus there exists at least a $\bm{\gamma}\in{\mathbb{N}}^{n}$ such that $\mathbf{d}^{\bm{\gamma}}_{\mathbf{x}}(f)$ is regular at $\mathbf{p}$. ∎

Now, we give an example to explain this and illustrate Definitions 3, 4 and Lemma 5.

Example 1**.**

Let $f=x_{1}+3\,x_{3}+4\,x_{4}-x_{1}^{2}+x_{3}^{2}-x_{4}^{2}-x_{2}^{3}$. For the exact point $\mathbf{p}=(0,0,0,0)$, we have the Taylor expansion of $f$ at $\mathbf{p}$ is:

[TABLE]

Since

[TABLE]

we know that $f$ is regular at $\mathbf{p}$. So is $\mathbf{d}^{(0,2,0,0)}_{\mathbf{x}}(f)=-3\,x_{2}$.

Similarly, for the approximate point $\tilde{\mathbf{p}}=(0.001,-0.001,0.002,-0.001)$ and a tolerance $\theta=0.01$, we have:

[TABLE]

%beginequation Since $|f(\tilde{\mathbf{p}})|=0.003002001<\theta$, $|\frac{\partial f}{\partial x_{2}}(\tilde{\mathbf{p}})|=3\cdot 10^{-6}<\theta$, $|\frac{\partial f}{\partial x_{1}}(\tilde{\mathbf{p}})|=0.998>\theta$, $|\frac{\partial f}{\partial x_{3}}(\tilde{\mathbf{p}})|=3.004>\theta$, $|\frac{\partial f}{\partial x_{4}}(\tilde{\mathbf{p}})|=4.002>\theta$, thus $f$ is $\theta$-regular at $\tilde{\mathbf{p}}$.

From this example, it’s easy to see that when compared with the exact case, the approximate zero $\tilde{\mathbf{p}}$ brings a small perturbation in the coefficients of the Taylor expansion of $f$ at $\tilde{\mathbf{p}}$. However, once given a proper $\theta$, we could acquire the same judging result as the exact case. For the above example, $f$ is regular at $\mathbf{p}$ and it is also $\theta$-regular at $\tilde{\mathbf{p}}$.

Definition 6**.**

Denote the operation set $\Delta\triangleq\{{+},\ {\cdot},\ {\partial}\}$, where $``\ {+}\ "$ denotes the sum of two polynomials, $``\ {\cdot}\ "$ denotes scalar multiplication and $``\ {\partial}\ "$ the differential of a polynomial. Given a polynomial system ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ and $\mathbf{p}\in{\mathbb{C}}^{n}$ such that ${\mathbf{F}}(\mathbf{p})=\mathbf{0}$, we define a polynomial set $\Delta_{\mathbf{p}}({\mathbf{F}})$, which satisfies:

1. $(1)$

${\mathbf{F}}\subset\Delta_{\mathbf{p}}({\mathbf{F}})$; 2. $(2)$

$\{a\,h|h\in\Delta_{\mathbf{p}}({\mathbf{F}}),a\in{\mathbb{C}}\backslash\{0\}\}\subset\Delta_{\mathbf{p}}({\mathbf{F}})$; 3. $(3)$

$\{h_{1}+h_{2}|h_{1}(\mathbf{p})+h_{2}(\mathbf{p})=0,h_{1},h_{2}\in\Delta_{\mathbf{p}}({\mathbf{F}})\}\subset\Delta_{\mathbf{p}}({\mathbf{F}})$; 4. $(4)$

$\{\frac{\partial h}{\partial x_{i}}|\frac{\partial h}{\partial x_{i}}(\mathbf{p})=0,i\in\{1,\ldots,n\},h\in\Delta_{\mathbf{p}}({\mathbf{F}})\}\subset\Delta_{\mathbf{p}}({\mathbf{F}})$.

Especially, for one polynomial $f\in{\mathbb{C}}[\mathbf{x}]$, we have the corresponding set $\Delta_{\mathbf{p}}(f)$.

The following lemma shows the relationship between the polynomials in $\Delta_{\mathbf{p}}({\mathbf{F}})$ and the polynomials in ${\mathbf{F}}$.

Lemma 7**.**

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ and $\mathbf{p}\in{\mathbb{C}}^{n}$, s.t. ${\mathbf{F}}(\mathbf{p})=\mathbf{0}$. $\forall g\in\Delta_{\mathbf{p}}({\mathbf{F}})$, we have

[TABLE]

where $a_{i,j}\in{\mathbb{C}}$ and $\bm{\gamma}_{i,j}\in{\mathbb{N}}^{n}$.

Proof.

The proof is obvious. ∎

We illustrate Definition 6 and Lemma 7 by the following example.

Example 2**.**

Let ${\mathbf{F}}=\{f_{1}=(x+y)^{2}+x^{3},f_{2}=x+y+y^{3}\}$. $\mathbf{p}=(0,0)$ is an isolated zero of ${\mathbf{F}}=0$. Let $h_{1}=\frac{\partial f_{1}}{\partial x}=2\,(x+y)+3\,x^{2}$, $h_{2}=\frac{\partial f_{1}}{\partial y}=2\,(x+y)$, $h_{3}=h_{1}-2\ f_{2}=3\,x^{2}-2\,y^{3}$, $h_{4}=\frac{\partial h_{3}}{\partial x}=6\,x$, $h_{5}=\frac{\partial^{2}h_{3}}{\partial y^{2}}=-12\,y$. It is clear that $h_{i}\in\Delta_{\mathbf{p}}({\mathbf{F}}),i=1,\ldots,5$ and $h_{i}$ has the form as (2).

3 Computing a deflated square system

Given a polynomial system with a multiple zero, Newton-type method usually is not used directly on the input system since it converges slowly or even doesn’t converge. Thus, deflation techniques are developed to transform the input system into another deflated system, which is regular at some zero whose certain projection is the given multiple zero.

In this section, given a polynomial system ${\mathbf{F}}\subset{\mathbb{C}}[\mathbf{x}]$ with an isolated singular zero $\mathbf{p}\in{\mathbb{C}}^{n}$, by employing some differential operations on the input polynomials or on the linear combinations of the related polynomials, we propose a new method to construct a new square system ${\mathbf{F}}^{\prime}\subset{\mathbb{C}}[\mathbf{x}]$, which satisfies that $\mathbf{p}$ is a simple zero of ${\mathbf{F}}^{\prime}={\mathbf{0}}$. We also prove the existence of ${\mathbf{F}}^{\prime}$ and show some properties of it.

First, let’s see a simple example to explain our idea.

Example 3**.**

Let ${\mathbf{F}}=\{f_{1}=x-y+x^{2},f_{2}=x-y+y^{2}\}$ with a 3-fold isolated zero $\mathbf{p}=(0,0)$. Obviously, $f_{1}$ and $f_{2}$ are already regular at $\mathbf{p}$. However, it is easy to find that the terms with degree one of $f_{1}$ and $f_{2}$ are linear dependent. Using $f_{2}-f_{1}$ to eliminate these terms of degree one, we get the polynomial $h=y^{2}-x^{2}$ and two new polynomials $\frac{\partial h}{\partial x}=-2\,x$, $\frac{\partial h}{\partial y}=2\,y$, which are both regular at $\mathbf{p}$. Selecting the two polynomials $f_{1}$ and $\frac{\partial h}{\partial y}$, we get a new square system ${\mathbf{F}}^{\prime}=\{x-y+x^{2},2\,y\}$, which has a regular zero $\mathbf{p}=(0,0)$. Moreover, it’s a system without perturbation.

Based on the idea in the above simple example, now we show our technique to construct a deflated square system below.

Assume that we have got the polynomials $g_{1},\ldots,g_{s}$, which are regular at $\mathbf{p}$, from the input polynomials $f_{1},\ldots,f_{s}$ such that

[TABLE]

Given one more polynomial $f_{s+1}$, we want to compute another polynomial $g_{s+1}$, s.t.

[TABLE]

Using only $g_{1},\ldots,g_{s}$ and $f_{s+1}$, we may not get the suitable $g_{s+1}$ if

[TABLE]

The input polynomials are needed in this case. Thus, we use $\{g_{1},\ldots,g_{s}\}\cup\{f_{1},\ldots,f_{s+1}\}$ to compute $g_{s+1}$. We will show how to compute $g_{s+1}$ in the following lemma.

Lemma 8**.**

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{s},f_{s+1},\ldots,f_{s+k}\}\subset{\mathbb{C}}[x_{1},\ldots,x_{n}](k\geq 1)$ and $\mathbf{p}\in{\mathbb{C}}^{n}$, s.t. ${\mathbf{F}}(\mathbf{p})={\mathbf{0}}$ and ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p}))=s$. Assume $\dim{\mathbb{V}}({\mathbf{F}})\leq n-s-1$ and $\deg({\mathbf{F}})=m(m>1)$. Then we can get a polynomial system ${\mathbf{F}}^{\prime}=\{f^{\prime}_{1},\ldots,f^{\prime}_{s},f^{\prime}_{s+1}\}$, which satisfies:

${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}}^{\prime})(\mathbf{p}))=s+1$, and $f_{j}^{\prime}\in\Delta_{\mathbf{p}}({\mathbf{F}})(1\leq j\leq s+1)$; 2. 2.

$\deg({\mathbf{F}}^{\prime})\leq m$.

Proof.

Without loss of generality, we assume that $\mathbf{p}$ is the origin and

[TABLE]

In the following, we consider the case of $s>0$, since if $s=0$, we can use the operators $\partial$ on $f_{i}(1\leq i\leq s+k)$ to get some polynomials, which are regular at $\mathbf{p}$.

To construct a polynomial system ${\mathbf{F}}^{\prime}$, s.t. ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}}^{\prime})(\mathbf{p}))=s+1$, we consider the rest polynomials $\{f_{s+1},\ldots,f_{s+k}\}$. Our proof is constructive.

First, $f_{i}(i=1,\ldots,s)$ has the form :

[TABLE]

where $\mathrm{T}_{i}\in{\mathbb{C}}[\mathbf{x}]$ and $\deg(\mathrm{T}_{i})=0$ or $\deg(\mathrm{T}_{i})\geq 2$. It’s easy to know that the row vector $\mathbf{a}_{i}=(a_{i1},\ldots,a_{in})(1\leq i\leq s)$ of the Jacobian matrix of $(f_{1},\ldots,f_{s})$ at $\mathbf{p}$ are linear independent since (3) holds.

Therefore, we can consider the following linear coordinate transformation $L$:

[TABLE]

With a realignment of the sequence of the variables $\{x_{1},\ldots,x_{n}\}$, we can always have the first $s$ columns of the coefficient matrix of $L$ being linear independent. Then $L$ is invertible. Denote the inverse of $L$ as $L^{-1}$. Let $\mathbf{p}^{\prime}=L(\mathbf{p})$ and $F_{i}=L^{-1}(f_{i})\in\mathbb{C}[{y}_{1},\ldots,{y}_{n}]$. We have:

[TABLE]

Since $\dim\mathbb{V}({\mathbf{F}})\leq n-s-1$ and $L^{-1}$ is invertible, it is obvious that

[TABLE]

Therefore, noticing that the terms with degree one of all $F_{i}(i=1,\ldots,s+k)$ in (4) contain only $s$ variables, there must be at least one of $\{L^{-1}(\mathrm{T}_{i}),i=1,\ldots,s+k\}$ containing at least one term, which has the form of $y_{s+1}^{d_{s+1}}y_{s+2}^{d_{s+2}}\cdots y_{n}^{d_{n}}$, such that $\sum\limits_{j=s+1}^{n}d_{j}>1$.

It’s easy to prove the claim. Suppose all $L^{-1}(\mathrm{T}_{i})(1\leq i\leq s+k)$ contain no terms of the form of $y_{s+1}^{d_{s+1}}y_{s+2}^{d_{s+2}}\cdots y_{n}^{d_{n}}$. Then, all the terms of $F_{i}(1\leq i\leq s+k)$ have the form of $y_{1}^{d_{1}}\cdots y_{s}^{d_{s}}y_{s+1}^{d_{s+1}}\cdots y_{n}^{d_{n}}$, $\sum\limits_{j=1}^{s}d_{j}>0$. In this case, the system $\{F_{1},\ldots,F_{s+k}\}$ vanishes on $\{y_{1}=0,\ldots,y_{s}=0\}$. Thus, we can verify easily that $\dim{\mathbb{V}}(F_{1},\ldots,F_{s+k})=n-s$, which contradicts with $\dim{\mathbb{V}}(F_{1},\ldots,F_{s+k})\leq n-s-1$. Thus, the claim is true.

Without loss of generality, we suppose that $L^{-1}(\mathrm{T}_{l})(l\in\{1,\ldots,s+k\})$ has the term with the form of $y_{s+1}^{d_{s+1}}y_{s+2}^{d_{s+2}}\cdots y_{n}^{d_{n}}$ and take the variable $y_{s+1}$ for example, i.e. $d_{s+1}\neq 0$. Further, we ask for the term with a lowest degree among all this kind of terms and denote the lowest degree as ${d}$. Then, we have:

[TABLE]

It is easy to see that $\gamma_{s+1}\neq 0,\ \deg(F_{s+1}^{\prime})<\deg(F_{l})$.

Thus, we have a new system $\{F_{1},\ldots,F_{s},F^{\prime}_{s+1}\}$. It’s easy to check that

[TABLE]

Finally, after doing the transformation $L$ on $F_{i}(1\leq i\leq s)$ and $F^{\prime}_{s+1}$, we have the new system ${\mathbf{F}}^{\prime}=\{f^{\prime}_{1},\ldots,f^{\prime}_{s+1}\}$, where

[TABLE]

By the definition of $\Delta_{\mathbf{p}}({\mathbf{F}})$ (see Definition 6), we can find that $f_{i}^{\prime}\in\Delta_{\mathbf{p}}({\mathbf{F}})(1\leq i\leq s+1)$. Therefore, we finished the first part of the proof.

From Lemma 7 and (5), it is easy to know that the maximal degree of $f^{\prime}_{i}(i=1,\ldots,s+1)$ is no larger than $m$. That is, $\deg({\mathbf{F}}^{\prime})\leq m$. Thus, we complete the proof. ∎

Now, we consider constructing a square system, which is regular at an isolated singular zero of the input system.

Theorem 9**.**

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{N}\}\subset{\mathbb{C}}[\mathbf{x}](N\geq n)$ be a polynomial system. $\mathbf{p}\in{\mathbb{C}}^{n}$ an isolated singular zero of ${\mathbf{F}}={\mathbf{0}}$ and $\deg({\mathbf{F}})=m$. Then there exists a square polynomial system ${\mathbf{F}}^{\prime}=\{f^{\prime}_{1},\ldots,f^{\prime}_{n}\}\subset\Delta_{\mathbf{p}}({\mathbf{F}})$, s.t.

$\mathbf{p}$* is an isolated regular zero of ${\mathbf{F}}^{\prime}={\mathbf{0}}$;* 2. 2.

$\deg({\mathbf{F}}^{\prime})\leq m$.

Proof.

Without loss of generality, assume that $\mathbf{p}$ is the origin. In the following, we will construct a square system by the polynomials in $\Delta_{\mathbf{p}}({\mathbf{F}})$.

First, we can choose a system ${\mathbf{F}}_{0}$ from ${\mathbf{F}}$, denoted as ${\mathbf{F}}_{0}=\{f_{1},\ldots,f_{r}\}$, whose Jacobian matrix at $\mathbf{p}$ has a maximal rank, s.t.

[TABLE]

If $r=n$, we finish the proof. Noticing that when $r=0$, we need only considering at least one of the polynomials in $f_{1},\ldots,f_{N}$ and can always get at least one polynomial, which is regular at $\mathbf{p}$ by Lemma 5. Thus, in the following, we consider the case of $1\leq r<n$.

First, considering the system $\{f_{1},\ldots,f_{r},f_{r+1},\ldots,f_{N}\}$, by Lemma 8, we can get a system

[TABLE]

s.t.

[TABLE]

Using the technique in Lemma 8, when considering the system ${\mathbf{F}}\cup\{f_{1}^{(1)},\ldots,f_{r}^{(1)},f_{r+1}^{(1)}\}$, we can get a system

[TABLE]

s.t.

[TABLE]

Repeat this process $n-r$ times and finally, we get a square system

[TABLE]

s.t.

[TABLE]

Thus, our final square system

[TABLE]

By Lemma 8, it is obvious that the maximal degree of $f^{\prime}_{i}(1\leq i\leq n)$ is no larger than $m$. That is, $\deg({\mathbf{F}}^{\prime})\leq m$. ∎

Remarks**.**

1. In the above construction process, we repeat $n-r$ times to get the deflated system ${\mathbf{F}}^{\prime}$. If considering all the variables simultaneously, we get more than one eligible polynomial each time in (5). Thus, the number of times in actual computation is less than $n-r$.

2. In the beginning of our construction, we also can compute all the related polynomials of all the input polynomials, which are regular at $\mathbf{p}$. Then, we choose a system from these polynomials, whose Jacobian matrix at $\mathbf{p}$ has a maximal rank. That’s to say that we make $r$ as big as possible to reduce our repeating steps.

Theorem 9 tells us that given a polynomial system ${\mathbf{F}}$ with an isolated singular zero $\mathbf{p}$, we can construct a new square system ${\mathbf{F}}^{\prime}$, which is regular at $\mathbf{p}$ and moreover, the degree of the polynomials in ${\mathbf{F}}^{\prime}$ does not increase. We give an example to illustrate our method.

Example 4**.**

(DZ2 [2]) Let ${\mathbf{F}}=\{f_{1}=x_{1}^{4},f_{2}=x_{1}^{2}x_{2}+x_{2}^{4},f_{3}=x_{3}+x_{3}^{2}-7\,x_{1}^{3}-8\,x_{1}^{2}\}$, which has a 16-fold zero $\mathbf{p}=(0,0,-1)$. The maximal degree of $f_{1},f_{2},f_{3}$ is 4. First, by the Taylor expansions of $f_{1},f_{2},f_{3}$ at $\mathbf{p}$, we have:

[TABLE]

It’s easy to find that only $f_{3}$ is regular at $\mathbf{p}$. Since $s={\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p}))=1$ and $\dim{\mathbb{V}}(f_{3},f_{2})=1$, we consider the system $\{f_{3},f_{2}\}$ directly. By Lemma 8, we have a system

[TABLE]

which satisfies ${\mathrm{rank}}({\bm{\mathrm{J}}}(f_{1}^{(1)},f_{2}^{(1)})(\mathbf{p}))=2$.

Next, we continue to consider the system $\{f_{1}^{(1)},f_{2}^{(1)}\}\cup{\mathbf{F}}$. Since $\dim{\mathbb{V}}(f_{1}^{(1)},f_{2}^{(1)},{\mathbf{F}})=0$, by Lemma 8, we have a system

[TABLE]

which satisfies ${\mathrm{rank}}({\bm{\mathrm{J}}}(f_{1}^{(2)},f_{2}^{(2)},f_{3}^{(2)})(\mathbf{p}))=3$.

Thus, we acquire the final square system ${\mathbf{F}}^{\prime}=\{f_{3},x_{2},4\,x_{1}\}$. It’s easy to check that $\mathbf{p}$ is a simple zero of ${\mathbf{F}}^{\prime}={\mathbf{0}}$ and the degree of every polynomial in ${\mathbf{F}}^{\prime}$ is no more than 4.

In this example, we repeat $n-s=2$ times to acquire the final square system ${\mathbf{F}}^{\prime}$. In fact, as what we say in Remark 2 of Theorem 9, computing twice is not necessary. Noticing that when computing $f_{2}^{(1)}=\mathbf{d}^{(2,0,0)}_{\mathbf{x}}(f_{2})=x_{2}$, we also can get $\mathbf{d}^{(1,1,0)}_{\mathbf{x}}(f_{2})=2\,x_{1}$. They are both regular at $\mathbf{p}$. It is easy to check that

[TABLE]

Thus, we obtain another square system ${\mathbf{F}}^{\prime}=\{f_{3},x_{2},2\,x_{1}\}$.

4 An effective version of our deflation method

In the section, by introducing some new variables to represent the coefficients of the linear combinations, we give an effective version of our deflation method. The deflated system produced by our deflation method has a simple zero, whose partial projection corresponds to the isolated singular zero of the input system. Furthermore, we also analyze the influences of the given tolerances $\theta$ and $\varepsilon$ to our method and show how to adjust their values to get a deflated system as exact as possible.

4.1 Parametric deflation system

Given a polynomial system ${\mathbf{F}}$ with an isolated singular zero $\mathbf{p}$, by employing some differential operations on the input polynomials directly or on the linear combinations of the related polynomials, we give a method to construct a new polynomial system ${\mathbf{F}}^{\prime}$ in Section 3, which satisfies that $\mathbf{p}$ is a simple zero of ${\mathbf{F}}^{\prime}={\mathbf{0}}$.

However, in practice, we can just get an approximate zero $\tilde{\mathbf{p}}$. As what we say in Example 1, the inexact value of $\tilde{\mathbf{p}}$ usually brings perturbations in the coefficients when doing the Taylor series expansions of the input polynomials at $\tilde{\mathbf{p}}$. Therefore, we can not do exact computations when adding two or more polynomials together. The inexact computations would produce a perturbed system of ${\mathbf{F}}^{\prime}$, which will lead to a bad final deflation result. We show an example to illustrate this case.

Example 5**.**

Continue with Example 3. Given an approximate zero

[TABLE]

Using the method in Theorem 9, we have $\tilde{h}=f_{2}+\tilde{\alpha}f_{1}$. By solving a Least Square problem, we can get $\tilde{\alpha}=-0.9984909264232$. Finally, we get an inexact system

[TABLE]

Obviously, we can not get a good result by the system $\widetilde{{\mathbf{F}}}^{\prime}$.

With a simple analysis, we can find that we couldn’t get an exact coefficient $\alpha$ of the linear combination of the polynomials with an approximate zero.

In the following, by introducing some new variables to represent the coefficients of the linear combinations, we give an effective version of our deflation method. Finally, the effective version of our deflation method will usually produce an exact deflated system, which has a simple zero, whose partial projection corresponds to the isolated singular zero of the input system. Furthermore, we also provide the size bound of our method. To our knowledge, it is the first time that considering the deflation of the polynomial system from the perspective of linear combination.

Similarly, before giving our theoretical results, we also show our main idea with a simple example first.

Example 6**.**

Still consider Example 3. Once given an approximate zero of the input system: $\tilde{\mathbf{p}}=(0.0006721,0.0008381)$, by Example 5, we know the coefficient $\tilde{\alpha}$ is inexact. Now we introduce a new variable $\alpha_{1}$. Let $h=f_{2}+\alpha_{1}f_{1}$ and compute

[TABLE]

Similar as in Example 5, we have $\tilde{\bm{\alpha}}_{1}=-0.9984909264232$. Given a tolerance $\varepsilon=0.05$, we have

[TABLE]

Do once again this process and introduce two new variables $\alpha_{2},\alpha_{3}$. Let

[TABLE]

and compute

[TABLE]

By solving another Least Square problem, we get the approximate values:

[TABLE]

Then, we have

[TABLE]

Thus, we get a polynomial system

[TABLE]

whose Jacobian matrix at $(\tilde{\mathbf{p}},\tilde{\alpha}_{1},\tilde{\alpha}_{2},\tilde{\alpha}_{3})$ has a full rank under the tolerance $\varepsilon$. In fact, we can find that $(0,0,-1,2,1)$ is a simple zero of $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})={\mathbf{0}}$ and the partial projection $(0,0)$ of $(0,0,-1,2,1)$ corresponds to the isolated singular zero $\mathbf{p}$ of the input system ${\mathbf{F}}$.

Given a polynomial system with an isolated zero, we have the following lemma.

Lemma 10**.**

[12]** Let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ be a polynomial system. $\mathbf{p}\in{\mathbb{C}}^{n}$ is an isolated singular zero of ${\mathbf{F}}={\mathbf{0}}$. $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{n})\in{\mathbb{C}}^{n}$ is a nonzero row vector, which satisfies ${\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p})\bm{\lambda}^{\mathrm{T}}={\mathbf{0}}$. For the new system

[TABLE]

we have the multiplicity of $\mathbf{p}$ in $\{{\mathbf{F}},{\mathbf{G}}\}={\mathbf{0}}$ is lower than the multiplicity of $\mathbf{p}$ in ${\mathbf{F}}={\mathbf{0}}$.

Remarks**.**

In Remark 2.1 of [8], the authors mentioned that deflation could also be constructed using the left null space. That is, we can replace ${\mathbf{G}}$ by the following system

[TABLE]

where $\bm{\lambda}\,{\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p})={\mathbf{0}}$. Furthermore, we have the following lemma.

Lemma 11**.**

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ be a polynomial system. $\mathbf{p}\in{\mathbb{C}}^{n}$ be an isolated singular zero of ${\mathbf{F}}={\mathbf{0}}$. Assume ${\mathrm{rank}}({\bm{\mathrm{J}}}(f_{1},\ldots,f_{s})(\mathbf{p}))={\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p}))=s$. Consider the augmented system

[TABLE]

where

[TABLE]

Then, we have:

there exists a unique $\hat{\bm{\alpha}}\in{\mathbb{C}}^{s}$ such that the system ${\mathbf{G}}$ has an isolated zero at $(\mathbf{p},\hat{\bm{\alpha}})$. 2. 2.

the multiplicity of ${\mathbf{G}}$ at $(\mathbf{p},\hat{\bm{\alpha}})$ is lower than that of ${\mathbf{F}}$ at $\mathbf{p}$.

Proof.

Let

[TABLE]

Denote the matrix $A=(a_{ij}),i=1,\ldots,s,j=1,\ldots,n$ and the row vector $b=(a_{s+1,1},\ldots,a_{s+1,n})$.

On one hand, when we fix $\mathbf{x}=\mathbf{p}$, the system

[TABLE]

is a linear system with respect to the variables $\alpha_{1},\ldots,\alpha_{s}$. Furthermore, it is easy to check that $\hat{\bm{\alpha}}$, which is determined by $AA^{\mathrm{T}}\hat{\bm{\alpha}}=-Ab^{\mathrm{T}}$, is the unique zero of ${\mathbf{H}}(\mathbf{p},\bm{\alpha})={\mathbf{0}}$. That is, there exists a unique $\hat{\bm{\alpha}}$ such that the system ${\mathbf{G}}$ has an isolated zero at $(\mathbf{p},\hat{\bm{\alpha}})$.

On the other hand, with the row operations, we could reduce the system ${\mathbf{G}}$ to the system

[TABLE]

where $l_{i}(\mathbf{x})$ are rational expressions and $\hat{\alpha}_{i}=l_{i}(\mathbf{p})$. Thus, considering the multiplicity of ${\mathbf{G}}$ at $(\mathbf{p},\hat{\bm{\alpha}})$ is equivalent to considering the multiplicity of ${\mathbf{G}}(\mathbf{x},\hat{\bm{\alpha}})$ at $\mathbf{p}$. Note that $(\hat{\alpha}_{1},\ldots,\hat{\alpha}_{s},1,0,\ldots,0)\,{\bm{\mathrm{J}}}({\mathbf{F}})(\mathbf{p})={\mathbf{0}}$. By Lemma 10 and (6), we know the second part holds. Thus, we finished the proof. ∎

In the above lemma, we construct $n$ new polynomials $h_{1},\ldots,h_{n}$. In fact, we can get them from the following way. Note that

[TABLE]

We know easily that ${\bm{\mathrm{J}}}(f_{s+1})(\mathbf{p})$ and ${\bm{\mathrm{J}}}(f_{1})(\mathbf{p}),\ldots,{\bm{\mathrm{J}}}(f_{s})(\mathbf{p})$ are linearly dependent. Thus, we can do the linear combination between $f_{s+1}$ and $f_{1},\ldots,f_{s}$ to eliminate this linear relationship. Let

[TABLE]

where new variables $\alpha_{i}$ are used to represent the coefficients of the linear combination. Compute all the derivatives of $g$ with respect to the variables $x_{1},\ldots,x_{n}$ and we get

[TABLE]

Thus, the above lemma tells us that after doing the linear combination of polynomials between $f_{s+1}$ and $f_{1},\ldots,f_{s}$, we get an augmented system ${\mathbf{G}}$, which satisfies that the multiplicity of ${\mathbf{G}}$ at $(\mathbf{p},\hat{\bm{\alpha}})$ is lower than that of ${\mathbf{F}}$ at $\mathbf{p}$. By repeating using the linear combination between polynomials in the original system and its related derivatives, we can construct a final deflated system, which processes an isolated simple zero. Denote $\mu$ be the multiplicity of ${\mathbf{F}}$ at $\mathbf{p}$. We do this repetitive process at most $\mu$ times.

Further, based on Lemma 11, we have the following theorem.

Theorem 12**.**

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ be a polynomial system. $\mathbf{p}\in{\mathbb{C}}^{n}$ be an isolated singular zero of ${\mathbf{F}}={\mathbf{0}}$. Denote $m=\deg({\mathbf{F}})$. Then there exists a square polynomial system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})=\{g_{1},\ldots,g_{t}\}\subset{\mathbb{C}}[\mathbf{x},\bm{\alpha}]$, s.t.

$({\mathbf{p}},\hat{\bm{\alpha}})\in{\mathbb{C}}^{t}$* is an isolated simple zero of $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})={\mathbf{0}}$;* 2. 2.

$t$* is bounded by $2^{\mu}\,n$, where $\mu$ is the multiplicity of $\mathbf{p}$ in ${\mathbf{F}}$;* 3. 3.

$\deg(\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha}))\leq m$.

Next, based on Lemma 11 and Theorem 12, we give an effective algorithm to compute a deflated square system from the input system with an approximate isolated singular zero below. It is an effective version of Lemma 8. $\theta$ is a tolerance to detect the regularity of the polynomials and we will talk about it in next subsection. $\varepsilon$ is another tolerance to judge the numerical rank of the Jacobian matrix at an approximate zero and we also talk about it in next subsection.

Remarks**.**

1. The termination and correctness of the algorithm is guaranteed by Lemma 11 and Theorem 12.

2. In the above algorithm, we compute polynomials of every $f_{i}$, which are regular at $\mathbf{p}$ at the beginning. Then, we put all these polynomials together to compute a system ${\mathbf{F}}_{0}$, such that the rank of its Jacobian matrix at $\mathbf{p}$ is maximal. This operation can make $r$ as big as possible. In some cases, we have $r=n$, which means we need not introduce new variables, such as Example 8. The aim of this preprocessing operation can speed up our algorithm.

Now, we give two examples to illustrate Algorithm 1.

Example 7**.**

Consider a polynomial system ${\mathbf{F}}=\{f_{1}=-\frac{9}{4}+\frac{3}{2}\,x_{1}+2\,x_{2}+3\,x_{3}+4\,x_{4}-\frac{1}{4}\,x_{1}^{2},f_{2}=x_{1}-2\,x_{2}-2\,x_{3}-4\,x_{4}+2\,x_{1}x_{2}+3\,x_{1}x_{3}+4\,x_{1}x_{4},f_{3}=8-4\,x_{1}-8\,x_{4}+2\,x_{4}^{2}+4\,x_{1}x_{4}-x_{1}x_{4}^{2},f_{4}=-3+3\,x_{1}+2\,x_{2}+4\,x_{3}+4\,x_{4}\}$. Given an approximate singular zero

[TABLE]

of ${\mathbf{F}}={\mathbf{0}}$ and the tolerance $\varepsilon=0.005$.

First, we have the Taylor expansion of $f_{3}$ at $\tilde{\mathbf{p}}$:

[TABLE]

[TABLE]

Consider the tolerance $\theta=0.05$. Since

[TABLE]

we get a polynomial

[TABLE]

which is $\theta$-regular at $\tilde{\mathbf{p}}$. Similarly, by the Taylor expansion of $f_{1},f_{2},f_{4}$ at $\tilde{\mathbf{p}}$, we have that $f_{1},f_{2},f_{4}$ are all $\theta$-regular at $\tilde{\mathbf{p}}$.

Thus, by Algorithm 1, we have ${\mathbf{G}}=\{f_{1},f_{2},-8+4\,x_{1}+4\,x_{4}-2\,x_{1}x_{4},f_{4}\}$. Compute

[TABLE]

We can choose ${\mathbf{H}}_{1}=\{h_{1}=f_{1},h_{2}=f_{2},h_{3}=-8+4\,x_{1}+4\,x_{4}-2\,x_{1}x_{4}\}$ from ${\mathbf{H}}={\mathbf{G}}\cup{\mathbf{F}}$. To $h_{4}=f_{4}\in{\mathbf{H}}\setminus{\mathbf{H}}_{1}$, let

[TABLE]

First, by solving a Least Square problem:

[TABLE]

we get an approximate value:

[TABLE]

Then, compute

[TABLE]

and we get the polynomial set

[TABLE]

which satisfies

[TABLE]

Thus, we get the final square system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})=\mathbf{H}_{1}$ and the point $\tilde{\bm{\alpha}}=(\tilde{\alpha}_{1},\tilde{\alpha}_{2},\tilde{\alpha}_{3})=(-1.000006509,-0.9997557989,0.000106178711)$.

In this example, given the input polynomial system ${\mathbf{F}}$ with an approximate singular zero $\tilde{\mathbf{p}}$, we can get a final square system by Algorithm 1 with only one step. In fact, $\alpha_{3}$ is not necessary to be introduced in this example by noticing that we can acquire a needed square system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})$ by using $F=f_{4}+\alpha_{1}f_{1}+\alpha_{2}f_{2}$. We give another example to illustrate the case that we do not introduce new variables.

Example 8**.**

(DZ2) Continue with Example 4. Given an approximate isolated singular zero

[TABLE]

and a tolerance $\varepsilon=0.005$, we use the Taylor series to expand $f_{i}(i=1,2,3)$ at $\tilde{\mathbf{p}}$ and compare all the coefficients with a tolerance $\theta=\varepsilon$. For $f_{1}$, we have

[TABLE]

It is obvious that only the absolute value of the coefficient of $(x_{1}-\tilde{p}_{1})^{4}$ is bigger than $\theta$. Therefore, compute $\mathbf{d}^{(3,0,0)}_{(x_{1},x_{2},x_{3})}(f_{1})=4x$, which is $\theta$-regular at $\tilde{\mathbf{p}}$. Similarly, for $f_{2},f_{3}$, we have the corresponding polynomial(s): $\{2x_{1},x_{2}\}\ \text{and}\ f_{3}.$ Thus, we have ${\mathbf{G}}=\{4\,x_{1},2\,x_{1},x_{2},f_{3}\}$. It is easy to check that

[TABLE]

Thus, we get the needed square system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x})={\mathbf{G}}=\{4x_{1},x_{2},f_{3}\}$.

In the above two examples, we assume that we have a right judgement on the tolerances $\theta$ and $\varepsilon$. In fact, the choice of the tolerances $\theta$ and $\varepsilon$ is important to our algorithm. Next, we give some analysis of them.

4.2 The analysis of $\theta$ and $\varepsilon$

As what we say in Example 1, $\theta$ is an important parameter in deciding if a polynomial is $\theta$-regular at $\tilde{\mathbf{p}}$. The other important parameter involved in our actual computation is $\varepsilon$, which is used to judge the numerical rank of the Jacobian matrix. Therefore, in this section, we will give some analysis about the parameters $\theta$ and $\varepsilon$.

First, we point out that $\theta$ is related to the absolute values of the coefficients of the Taylor expansion of the polynomial at its approximate zero.

For example, given a polynomial $f=x^{2}+10000\,y^{2}$ with an approximate zero $\tilde{\mathbf{p}}=(\tilde{p}_{1},\tilde{p}_{2})=(0.0006851,-0.0004368)$, we have the Taylor expansion of $f$ at $\tilde{\mathbf{p}}$:

[TABLE]

Given $\theta=0.5$, we have

[TABLE]

Thus, we draw the conclusion that $f$ is $\theta$-regular at $\tilde{\mathbf{p}}$. However, considering that the lowest degree of $f$ is 2, we know that $f$ is singular at the exact zero $\mathbf{p}=(0,0)$ actually, which is a different result from the case of $\tilde{\mathbf{p}}$. That means $\theta$ is not chosen properly. The main reason is that the coefficient of $f$ has a great fluctuation or the accuracy of $\tilde{\mathbf{p}}$ is not high enough. If given another approximate zero $\tilde{\mathbf{q}}=(\tilde{q}_{1},\tilde{q}_{2})=(0.000006851,-0.000004368)$ with higher precision, we have:

[TABLE]

By this time, using the same $\theta=0.5$, we have $f$ is $\theta$-singular at $\tilde{\mathbf{q}}$, which is the same judgement as the exact case of $\mathbf{p}$.

In actual computation, to deal with this case, we give one solution: For a nonzero polynomial $f\in{\mathbb{C}}[\mathbf{x}]$, let $\Gamma_{f}$ be a set of the absolute values of all the coefficients of $f$. We denote the maximal and minimal ones inside $\Gamma_{f}$ as $M=\max(\Gamma_{f})$ and $m=\min(\Gamma_{f})$ respectively. If $m/M\leq 10^{-a}$, we regard that the coefficients of $f$ fluctuate a lot and take $\epsilon=(m+M)/2M$; Else, we take $\theta=(m+M)/(2M\times 10^{a})$, where $a\in\mathbb{N}$ is related to the precision of the given approximate zero $\tilde{\mathbf{p}}$. For example, if the accuracy of the given approximate zero $\tilde{\mathbf{p}}$ has three significant digits, we can take $a=3$. Of course, we can overcome this problem thoroughly by refining the approximate zero to a higher precision with the input system if the Jacobian matrix of the system at $\tilde{\mathbf{p}}$ is numerically nonsingular.

In summary, the reason for the above situation is that we judge a polynomial, which is singular at the exact zero $\mathbf{p}$, as a polynomial being $\theta$-regular at the approximate zero $\tilde{\mathbf{p}}$.

The other situation is that a polynomial, which is regular at the exact zero $\mathbf{p}$, may be judged as a polynomial being $\theta$-singular at the approximate zero $\tilde{\mathbf{p}}$.

For example, consider the polynomial $f=\frac{1}{20}x+x^{2}+10000y^{2}$ with the approximate zero $\tilde{\mathbf{q}}=(\tilde{q}_{1},\tilde{q}_{2})=(0.000006851,-0.000004368)$. We have:

[TABLE]

Still use $\theta=0.5$ and we get the judgement that $f$ is $\theta$-singular at $\tilde{\mathbf{q}}$. In fact, $f$ is regular at $\mathbf{p}=(0,0)$. One way to deal with this case is that we can take a smaller $\theta$. When we take $\theta=0.05$, we will acquire the appropriate result.

From the above analysis about the tolerance $\theta$, we know that the choice of $\theta$ is crucial to our method. We give a further theoretical analysis about the tolerance $\theta$ below. Here, we assume that the judgement of the other tolerance $\varepsilon$, which is used to decide the numerical rank of the Jacobian matrix at the approximate zero, is correct.

Let $\theta$ be a tolerance. Assume that we have computed an intermediate system ${\mathbf{H}}=\{h_{1},\ldots,h_{s}\}\subset{\mathbb{C}}[\mathbf{x}^{\prime}]$. Denote $\mathbf{x}^{\prime}=(\mathbf{x},\bm{\alpha})$. Assume that $\mathbf{p}$ is an isolated singular zero of the original system. The exact value of $\bm{\alpha}$ related to the coefficients of linear combinations is $\hat{\bm{\alpha}}$. Denote $\mathbf{p}^{\prime}=(\mathbf{p},\hat{\bm{\alpha}})$. Let $\tilde{\mathbf{p}}^{\prime}$ be an approximate zero of ${\mathbf{H}}$ related to $\mathbf{p}^{\prime}$ such that

[TABLE]

Next, we consider one more polynomial $h\in{\mathbb{C}}[\mathbf{x}^{\prime}]$. If $h$, which is regular at $\mathbf{p}^{\prime}$, is judged as being $\theta$-singular at $\tilde{\mathbf{p}}^{\prime}$, we may get a perturbed system finally. Specifically, compute the Taylor expansion of $h$ at $\tilde{\mathbf{p}}^{\prime}$:

[TABLE]

Since $h$ is $\theta$-singular at $\tilde{\mathbf{p}}^{\prime}$, we know that $|h(\tilde{\mathbf{p}}^{\prime})|<\theta$ and all $|\frac{\partial h(\tilde{\mathbf{p}}^{\prime})}{\partial x_{j}}|<\theta$. Thus, we compute

[TABLE]

If there exists some $j$ such that

[TABLE]

we may derive a perturbed system in the end, where $\frac{\partial h}{\partial x_{j}}$ has and only has one perturbed term $\frac{\partial h(\mathbf{p}^{\prime})}{\partial x_{j}}$ compared to the polynomial $\frac{\partial h}{\partial x_{j}}-\frac{\partial h(\mathbf{p}^{\prime})}{\partial x_{j}}$ which vanishes at $\mathbf{p}^{\prime}$.

For other cases, if

[TABLE]

it is clear that $\frac{\partial h}{\partial x_{j}}$ vanishes at $\mathbf{p}^{\prime}$. Thus it is exact. If

[TABLE]

according to our constructive method, we should do the linear combination

[TABLE]

and compute its derivatives. Thus the perturbed term $\frac{\partial h(\mathbf{p}^{\prime})}{\partial x_{j}}$ disappears. We will get an exact polynomial which vanishes at $\mathbf{p}^{\prime}$ in the end. Notice that if $h_{i}$’s have perturbed terms, which are constants $h_{i}(\mathbf{p}^{\prime})$. We know that if we compute the derivatives of $f$, these terms will disappear. Thus whether $h_{i}$’s have perturbed terms or not, the polynomials in the final deflated system derived by the linear combinations vanish at the exact zero $\mathbf{p}^{\prime}$.

Now we consider the case that $h$ is regarded as $\theta$-regular at $\tilde{\mathbf{p}}^{\prime}$ while it is singular at $\mathbf{p}^{\prime}$. If

[TABLE]

we will do the linear combination of $h$ and $h_{1},\ldots,h_{s}$ and compute its derivatives. It is obvious that this operation has no influence on our result. Usually the case

[TABLE]

will not happen. It is related to the numerical computation of the rank of the Jacobian matrix of $({\mathbf{H}},h)$ at $\tilde{\mathbf{p}}^{\prime}$.

As a summary of the foregoing analysis, we have:

Let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}\subset{\mathbb{C}}[\mathbf{x}]$ be a polynomial system. $\tilde{\mathbf{p}}\in{\mathbb{C}}^{n}$ is an approximate zero of ${\mathbf{F}}={\mathbf{0}}$ and $\theta$ is a tolerance. According to our method, we acquire a final system $\widetilde{{\mathbf{F}}}^{\prime}\subset{\mathbb{C}}[\mathbf{x},\bm{\alpha}]$. During we compute the final system $\widetilde{{\mathbf{F}}}^{\prime}$,

if we judge a polynomial, which is singular at the exact zero $\mathbf{p}$, as being $\theta$-regular at $\tilde{\mathbf{p}}$, the final system $\widetilde{{\mathbf{F}}}^{\prime}$ is accurate. 2. 2.

if we judge a polynomial, which is regular at the exact zero $\mathbf{p}$, as being $\theta$-singular at $\tilde{\mathbf{p}}$, the final system $\widetilde{{\mathbf{F}}}^{\prime}=\widetilde{{\mathbf{F}}}+\bm{\vartheta}$, is a perturbed system, where $\widetilde{{\mathbf{F}}}$ is an accurate system and $\bm{\vartheta}$ is the perturbed term, which satisfies $\max\limits_{i}|\vartheta_{i}|<\theta$.

In actual computation, to make our method as accurate as possible, we give an adaptive adjustment step at the end of our algorithm. To be specific, assume that the initial tolerance $\theta=\theta_{1}$. After the refining steps, denote the refined zero as $\bar{\mathbf{p}}$. We compute the Taylor expansions of all the related polynomials in computing the system $\widetilde{{\mathbf{F}}}^{\prime}$ at $\bar{\mathbf{p}}$, including all the input polynomials. We denote the maximal absolute value of both the coefficients of the polynomials, which are judged as $\theta_{1}$-singular at $\tilde{\mathbf{p}}$ and the polynomials, which are judged as $\theta_{1}$-regular at $\tilde{\mathbf{p}}$, as $\theta_{2}$. It is also the term named “Max err” in Tables 1 and 2 in the next section.

It is easy to imagine that $\theta_{2}\leq\theta_{1}$ usually. If $\theta_{2}$ has a very higher precision than $\theta_{1}$, such as $\theta_{1}=10^{-2}$ and $\theta_{2}=10^{-13}$, we are sure that our conclusion is exact. If $\theta_{2}>\theta_{1}$ or $\theta_{2}$ still has a bad accuracy, such as $\theta_{1}=10^{-2}$ and $\theta_{2}=10^{-1}$ or $\theta_{2}=10^{-4}$, we will take a smaller $\theta<\min\{\theta_{1},\theta_{2}\}$ and repeat our method again.

After repeating our method several times, if $\theta_{2}$ is still bad, we will merely get a perturbed system.

Now, we give two examples to explain the above analysis.

Example 9**.**

Given a polynomial system ${\mathbf{F}}=\{f_{1}=x+x^{2}+10000y^{2},f_{2}=x^{2}+10000y^{2}\}$ with an approximate zero

[TABLE]

Consider the tolerances $\varepsilon=0.05$ and $\theta=0.5$. By the Taylor expansions of $f_{i}$ at $\tilde{\mathbf{p}}$, we know that $f_{1},f_{2}$ are both $\theta$-regular at $\tilde{\mathbf{p}}$.

Next, according to Algorithm 1, we compute

[TABLE]

Thus, we can use Newton’s method to refine $\tilde{\mathbf{p}}$ to a higher accuracy and get

[TABLE]

At this time, it’s easy to check that $f_{1}$ is $\theta$-regular at $\tilde{\mathbf{p}}^{\prime}$ and $f_{2}$ is $\theta$-singular at $\tilde{\mathbf{p}}^{\prime}$. Therefore, for $f_{2}$, we have

[TABLE]

which are both $\theta$-regular at $\tilde{\mathbf{p}}^{\prime}$. Furthermore,

[TABLE]

Thus, we get the final system $\widetilde{{\mathbf{F}}}^{\prime}=\{f_{1},\ 20000y\}$. After applying Newton’s method, we get the refined zero $\bar{\mathbf{p}}=(\bar{p}_{1},\bar{p}_{2})=10^{-16}\cdot(0.53016,0)$.

At last, we check if our chosen $\theta$ is proper. We compute the Taylor expansion of all the polynomials, which is judged as $\theta$-singular at $\tilde{\mathbf{p}}$, at the refined zero $\bar{\mathbf{p}}$ and get:

[TABLE]

Thus, we have

[TABLE]

which means that our final system $\widetilde{{\mathbf{F}}}^{\prime}$ is more accurate than before.

Example 10**.**

Consider the system ${\mathbf{F}}=\{f_{1}=x+x^{2}+2xy+10000y^{2},f_{2}=\frac{1}{20}x+x^{2}+2xy+10000y^{2}\}$ with an approximate zero

[TABLE]

Let the tolerances $\varepsilon=0.05$ and $\theta=0.5$. Similarly, by the Taylor expansions of $f_{i}$ at $\tilde{\mathbf{p}}$, we know that $f_{1}$ is $\theta$-regular at $\tilde{\mathbf{p}}$ and $f_{2}$ is $\theta$-singular at $\tilde{\mathbf{p}}$. Therefore, we have

[TABLE]

Compute

[TABLE]

Thus, we get the final system

[TABLE]

Obviously, $\widetilde{{\mathbf{F}}}_{1}^{\prime}$ is a perturbed system and $\vartheta_{2}=\frac{1}{20}$ is the perturbed term, which satisfies $|{\vartheta}_{2}|<\theta$. It’s easy to imagine that with $\widetilde{{\mathbf{F}}}_{1}^{\prime}$, we could not get a good result. The main reason is that $\theta=0.5$ is too big, which leads to a wrong judgement on whether $f_{2}$ is $\theta$-regular at $\tilde{\mathbf{p}}$.

If given another smaller tolerance $\theta^{\prime}=0.05$, we will get a right judgement that $f_{2}$ is $\theta^{\prime}$-regular at $\tilde{\mathbf{p}}$. Thus, we consider the linear combination of $f_{1}$ and $f_{2}$. Let $f=f_{2}+\alpha f_{1}$ and compute

[TABLE]

where $\alpha$ is a new variable and its initial value $\tilde{\alpha}=-0.050076986$. Compute

[TABLE]

Thus, we get the final system $\widetilde{{\mathbf{F}}}^{\prime}=\{f_{1},g_{1},g_{2}\}$. Similarly, we consider applying Newton’s method on the final system $\widetilde{{\mathbf{F}}}^{\prime}$ and get the refined zero:

[TABLE]

Then, we check the coefficients of the terms with degree one of the Taylor expansion of $f$ at $\bar{\mathbf{p}}$ and get

[TABLE]

Thus, we are sure that our final system $\widetilde{{\mathbf{F}}}^{\prime}$ is accurate. Here, “0” is not exact zero but means in Matlab machine accuracy.

From the above two examples, we can see that once given an appropriate tolerance $\theta$, we can make sure that our final system is accurate. Otherwise, what we acquired is just a perturbed system, such as the system $\widetilde{{\mathbf{F}}}^{\prime}_{1}$ in Example 10.

Next, we continue analyzing the other tolerance $\varepsilon$, which is used to judge the numerical rank of a matrix. That is, we determine the numerical rank by comparing the absolute values of the singular values of the Jacobian matrix at approximate zero with the tolerance $\varepsilon$. Specifically, assume that we have computed an intermediate system ${\mathbf{H}}=\{h_{1},\ldots,h_{s}\}\subset{\mathbb{C}}[\mathbf{x}^{\prime}]$. Denote $\mathbf{x}^{\prime}=(\mathbf{x},\bm{\alpha})$. Assume that $\mathbf{p}$ is an isolated singular zero of the original system. The exact value of $\bm{\alpha}$ related to the coefficients of linear combinations is $\hat{\bm{\alpha}}$. Denote $\mathbf{p}^{\prime}=(\mathbf{p},\hat{\bm{\alpha}})\in{\mathbb{C}}^{t}$. Let $\tilde{\mathbf{p}}^{\prime}\in{\mathbb{C}}^{t}$ be an approximate zero of ${\mathbf{H}}$ related to $\mathbf{p}^{\prime}$ such that

[TABLE]

Next, we consider one more polynomial $h_{s+1}\in{\mathbb{C}}[\mathbf{x}^{\prime}]$. Given the tolerance $\theta$, we can compute a polynomial $h$ from $h_{s+1}$, which is $\theta$-regular at $\tilde{\mathbf{p}}^{\prime}$. Denote

[TABLE]

For simplicity, we denote the deflated system as ${\mathbf{H}}^{\prime}$, which comes from $\{{\mathbf{H}},\,h_{s+1}\}$ after one step deflation, and its corresponding exact zero as $\mathbf{q}$, whose partial projection is $\mathbf{p}^{\prime}$.

According to the above analysis of $\theta$, for $h_{s+1}$, we have the following cases:

if $\theta$ is chosen properly, that is, we judge $h_{s+1}$, which is regular or singular at $\mathbf{p}^{\prime}$, as being $\theta$-regular or $\theta$-singular at $\tilde{\mathbf{p}}^{\prime}$ respectively, we know that $h$ is regular at $\mathbf{p}^{\prime}$. Thus, we have:

1. (a).

if $r_{2}=r_{1}$, we, of course, get an exact system ${\mathbf{H}}^{\prime}$. That is, ${\mathbf{H}}^{\prime}(\mathbf{q})={\mathbf{0}}$. 2. (b).

if $r_{2}<r_{1}$, according to our algorithm, we consider do the linear combination:

[TABLE]

and compute all the derivatives of $g$ with respect to all variables: $g_{i}=\frac{\partial g}{\partial x_{i}^{\prime}},\ i=1,\ldots,t$. Correspondingly, ${\mathbf{H}}^{\prime}=\{h_{1},\ldots,h_{s},g_{1},\ldots,g_{t}\}$. Note that ${\bm{\mathrm{J}}}(h)(\mathbf{p}^{\prime})$ and ${\bm{\mathrm{J}}}(h_{1})(\mathbf{p}^{\prime}),\ldots,{\bm{\mathrm{J}}}(h_{s})(\mathbf{p}^{\prime})$ are actually linear independent, which means that the equations $(\alpha_{1},\ldots,\alpha_{s},1)\,{\bm{\mathrm{J}}}({\mathbf{H}})(\mathbf{p}^{\prime})={\mathbf{0}}$ has no solution. Thus, although we can give the initial value $\tilde{\alpha}_{j}$ of $\alpha_{j}$ by solving a Least Squares problem, the linear independent will bring us some inexact polynomials $g_{i}$, which means $g_{i}(\mathbf{q}^{\prime})\neq 0$. Further, we may get a perturbed system ${\mathbf{H}}^{\prime}$. That is, ${\mathbf{H}}^{\prime}(\mathbf{q})\neq{\mathbf{0}}$. 3. (c).

if $r_{1}<r_{2}$, we will add $h$ to the system ${\mathbf{H}}$ directly and get an exact system ${\mathbf{H}}^{\prime}=\{h_{1},\ldots,h_{s},h\}$. 2. 2.

if $\theta$ is chosen too big, that is, we judge $h_{s+1}$, which is regular at $\mathbf{p}^{\prime}$, as being $\theta$-singular at $\tilde{\mathbf{p}}^{\prime}$, we may get a perturbed polynomial $h$, which means that $h(\mathbf{p}^{\prime})\neq 0$ and $h-h(\mathbf{p}^{\prime})$ is regular at $\mathbf{p}^{\prime}$. Thus, we have:

1. (a).

if $r_{2}=r_{1}$, only the choice of $\theta$ affects our final result. Thus, we may get a perturbed system ${\mathbf{H}}^{\prime}$ in this case. 2. (b).

if $r_{2}<r_{1}$, with a similar discussion with the case of 1(b), we get a perturbed system ${\mathbf{H}}^{\prime}$. 3. (c).

if $r_{1}<r_{2}$, we add $h$ to $\{h_{1},\ldots,h_{s}\}$ directly and get a perturbed system ${\mathbf{H}}^{\prime}$. 3. 3.

if $\theta$ is chosen too small, that is, we judge $h_{s+1}$, which is singular at $\mathbf{p}^{\prime}$, as being $\theta$-regular at $\tilde{\mathbf{p}}^{\prime}$, we know that $h=h_{s+1}$ is singular at $\mathbf{p}^{\prime}$. Thus, we have:

1. (a)

if $r_{2}=r_{1}$, only the choice of $\theta$ affects our result. Thus, the system ${\mathbf{H}}^{\prime}$ is exact in this case. 2. (b)

if $r_{2}<r_{1}$, similar to the case of 1(b), we consider do the linear combination:

[TABLE]

One difference from 1(b). is that $h$ is singular at $\mathbf{p}^{\prime}$. Thus, all the $g_{i}=\frac{\partial g}{\partial x_{i}^{\prime}}$ are exactly vanish at $\mathbf{q}$. So, the system ${\mathbf{H}}^{\prime}$ is also exact in this case. 3. (c)

if $r_{1}<r_{2}$, we add $h$ to $\{h_{1},\ldots,h_{s}\}$ directly and get an exact system ${\mathbf{H}}^{\prime}$.

With the above analysis, we know that the choice of the tolerances $\theta$ and $\varepsilon$ has an influence on our final deflated system: an exact deflated system or a perturbed system. To judge which case a final deflated system belongs to, we use the following judgment method:

Denote the input system as ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}$, the final deflated system $\widetilde{{\mathbf{F}}}^{\prime}$. Noticing that we use Newton’s method to refine the system $\widetilde{{\mathbf{F}}}^{\prime}$, thus, we denote Newton’s iteration sequence as $\{\tilde{\mathbf{p}}_{l},\,l\geq 1\}$ and the final certified zero $\tilde{\mathbf{p}}^{\prime}$.

• •

First, we check if Newton’s iteration sequence $\{\tilde{\mathbf{p}}_{l},\,l\geq 1\}$ is quadratic convergence. If not, we claim that our deflated system is a perturbed system.

• •

If it is, we compute

[TABLE]

• •

Next, we give a tolerance $\theta^{\prime}$, which is usually a very small value, and compare the magnitude of $\theta^{\prime}$ and $\Delta$. If $\Delta<\theta^{\prime}$, we regard the final deflated system $\widetilde{{\mathbf{F}}}^{\prime}$ as an exact system; otherwise, $\widetilde{{\mathbf{F}}}^{\prime}$ is a perturbed system.

Of course, for the exact case, we are done. For the perturbed case, we hope to make our final deflated system as accurate as possible by adjusting the values of $\theta$ and $\varepsilon$. However, we still do not have a good idea on how to distinguish the effect of the two tolerances $\theta$ and $\varepsilon$ on the final system. Fortunately, noting that the tolerance $\theta$ is used to accelerate our algorithm and is not necessary, therefore, according to the remark of Lemma 10, we can use the deflation construction 6 to compute the final system. In this case, we just need to consider the tolerance $\varepsilon$, which is used to judge the numerical rank of the Jacobian matrix. That’s to say, even if the first two steps of Algorithm 1 are removed, our deflation construction process can still work well. Considering the possible judgment, our final system can also be a perturbed system. Next, we give a possible modified method to overcome this case.

let ${\mathbf{F}}=\{f_{1},\ldots,f_{n}\}$ be the input system, $\tilde{\mathbf{p}}\in{\mathbb{C}}^{n}$ be the initial approximate zero. Let $\varepsilon$ and $\theta^{\prime}$ be the given tolerances.

• •

First, assume ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\tilde{\mathbf{p}}),\varepsilon)=n$. We apply Newton’s method on the system ${\mathbf{F}}$ and get the refined zero $\tilde{\mathbf{p}}^{\prime}$. Next, we check if Newton’s iteration sequence is quadratic convergence. If it is, we continue comparing the magnitude of $\theta^{\prime}$ and $\Delta$. If $\Delta<\theta^{\prime}$, we regard ${\mathbf{F}}$ as a system with an isolated simple zero; otherwise, we know ${\mathbf{F}}$ is a system with a multiple zero and ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\tilde{\mathbf{p}}),\varepsilon)<n$.

• •

Assume ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\tilde{\mathbf{p}}),\varepsilon)=n-1$. After using the deflation construction in Algorithm 1 once(from step 3 to step 10), we get a deflation system $\widetilde{{\mathbf{F}}}_{1}$ and an approximate zero $\tilde{\mathbf{p}}_{1}$. Then, we consider all the possibilities of ${\mathrm{rank}}({\bm{\mathrm{J}}}(\widetilde{{\mathbf{F}}}_{1})(\tilde{\mathbf{p}}_{1}),\varepsilon)$. For every case, we go on our deflation construction in Algorithm 1 and use our mentioned judging method to check which case the final deflated system belongs to. As long as the final deflated system is judged to be an exact system, we will stop our deflation process; Otherwise, we know ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\tilde{\mathbf{p}}),\varepsilon)<n-1$.

• •

Assume ${\mathrm{rank}}({\bm{\mathrm{J}}}({\mathbf{F}})(\tilde{\mathbf{p}}),\varepsilon)=n-2$. We consider as the case of $n-1$.

About the above judgment process, we have two things to say:

The above judgement process must terminate in finite steps considering that our deflation construction terminates in finite steps.

• 2.

In the above judgement process, we traverse all the possibilities of the rank of the Jacobian matrix. Thus, there must be at least one case that we get an exact deflated system.

Now, we give an example to illustrate our idea.

Example 11**.**

Continue with Example 10. Here, we only use the tolerance $\varepsilon=0.05$ to judge the numerical rank. First, we compute

[TABLE]

Thus, we consider using Newton’s method on the system ${\mathbf{F}}$ directly. Given an iterative error $10^{-8}$, we have the following Newton’s iteration sequence:

We can check easily that Newton’s iteration sequence $\{\tilde{\mathbf{p}}_{j},\ j=1,\ldots,10\}$ is linear convergence. According to our judging criteria, we know that

[TABLE]

Next, according to our construction process in Algorithm 1, we let

[TABLE]

and compute

[TABLE]

We have $\tilde{\alpha}=-0.091484814324$.

Next, let $\widetilde{{\mathbf{F}}}_{1}=\{f_{1},g_{1},g_{2}\}$ and $\tilde{\mathbf{p}}_{1}=(\tilde{\mathbf{p}},\tilde{\alpha})$. By our given revised method above, we continue considering all the possibilities of ${\mathrm{rank}}({\bm{\mathrm{J}}}(\widetilde{{\mathbf{F}}}_{1})(\tilde{\mathbf{p}}_{1}),\varepsilon)$. For example, we consider the case of

[TABLE]

Similarly, given the iterative error $10^{-8}$, by using Newton’s method on the system $\widetilde{{\mathbf{F}}}_{1}$, we get the following iteration sequence:

It is easy to check that the iteration sequence $\{\tilde{\mathbf{p}}_{j},j=1,2,3,4\}$ is quadratic convergence. Furthermore, given a tolerance $\theta^{\prime}=10^{-12}$, we can compute

[TABLE]

and verify that $\Delta<\theta^{\prime}$. Thus, we regard the final deflated system $\widetilde{{\mathbf{F}}}^{\prime}=\widetilde{{\mathbf{F}}}_{1}$ as an exact system. At the same time, we stop our deflation process.

Until now, we have finish all the discussions about the tolerances $\theta$ and $\varepsilon$. Once given a polynomial system ${\mathbf{F}}\subset{\mathbb{C}}[\mathbf{x}]$ with an isolated singular zero, we use Algorithm 1 to compute a new system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})$, which has a simple zero. What’s more, according to the analysis of the tolerances $\theta$ and $\varepsilon$, our final system $\widetilde{{\mathbf{F}}}^{\prime}$ is an accurate system usually. For the perturbed case, we also give one ergodic way to adjust our final result as accurate as possible. Thus, we can use the final system $\widetilde{{\mathbf{F}}}^{\prime}$ to certify the isolated zeros of the input system.

In the following, by using the algorithm verifynlss in INTLAB[22], we give an example to explain how we certify the isolated singular zero of the input system.

Example 12**.**

Continue with Example 7. Applying the algorithm verifynlss, we get the system:

[TABLE]

*and two verified inclusions *

[TABLE]

and

[TABLE]

*For the deflated system $\widetilde{{\mathbf{F}}}(\mathbf{x},\bm{\alpha})$, we affirm that there is a unique isolated simple zero $(\hat{\mathbf{x}},\hat{\bm{\alpha}})\in(\mathbf{X},\bm{\mathcal{A}})$, such that $\widetilde{{\mathbf{F}}}(\hat{\mathbf{x}},\hat{\bm{\alpha}})={\mathbf{0}}$. What’s more, the projection $\hat{\mathbf{x}}$ of $(\hat{\mathbf{x}},\hat{\bm{\alpha}})$ corresponds to the isolated singular zero of the input system ${\mathbf{F}}$. That’s to say, we certified the isolated singular zeros of the original system. *

5 Experiments and results

We implement our method in Matlab of Algorithm 1. The code and some examples can be found in http://www.mmrc.iss.ac.cn/~jcheng/VDSS. In this section, we show the results of the experiment and the comparison of our method with some other methods. We do the experiments in Matlab R2012b with INTLAB-V5.5 on a computer with Windows 7, Intel $i7$ processor and $8GB$ memory.

in [16], by modifying the method proposed by Yamamoto [18], they give a deflation method to compute a regular and square augmented system. which can be used to prove the existence of an isolated singular solution of a slightly perturbed system. Moreover, by applying INTLAB function verifynlss[22], they also give an algorithm $\mathbf{viss}$ to compute verified error bounds. However, noticing that their method is essentially a deflation method. Thus, we also implement our algorithm based on INTLAB function verifynlss.

In Table 1, we compare our algorithm $\mathbf{VDSS}$ with the algorithm $\mathbf{viss}$. These examples are relatively simple and small scale, which can be found in [2, 16]. We also list them in http://www.mmrc.iss.ac.cn/~jcheng/VDSS/fun.m. We denote $var$ the number of polynomials, $mul$ the multiplicity and Verified acc the final verified accuracy, which is measured by the breath of the verified inclusion $\mathbf{X}$. And Max err is $\delta_{2}$ as mentioned in Section 4.2. We use a same initial accuracy $10^{-4}$ for all the examples. “true” means we get two same endpoints of the verified inclusion. When the term for Max err is “0”, it does not mean Max err is exactly zero and only shows in Matlab machine precision.

From Table 1, we can see that our algorithm is effective. On one hand, the verified accuracy of our method is never worse than $\mathbf{viss}$ for all these examples. On the other hand, thanks to our acceleration strategies, our practical size and computing time are smaller than those of $\mathbf{viss}$ in most cases.

We also compare our method with viss for large-scale polynomial systems. All the examples in Table 2 can be found in http://www.mmrc.iss.ac.cn/~jcheng/VDSS/example.m. The example LZ2000 can be found in [14]. The example nonpoly3 is a non-polynomial nonlinear system. We construct the other examples as below: First, we produce some polynomials randomly to form a zero-dimensional system $\{f_{1},\ldots,f_{n}\}$, which has a simple zero $\mathbf{p}$ and $\deg(f_{i})\geq 2$ usually. The final systems have the form: ${\mathbf{F}}=\{f_{i}^{d_{i}}+g_{i},1\leq i\leq n,g_{i}\in\{f_{1}^{d^{\prime}_{1}},\ldots,f_{n}^{d^{\prime}_{n}},0\},d_{i}\geq 1,d^{\prime}_{i}\geq 1,1\leq i\leq n\}$. The new systems are always dense polynomial systems. The examples named simple1, reduce3, big1, big2, big3, large3, large6, large8 are of the form that $g_{i}=0(1\leq i\leq n)$; The examples named addvar3, unre3, unre5, rankone2, rankone3 are of the form that $g_{i}$ are not all zeros. The ranks of the Jacobian matrices of the examples rankone2, rankone3 at the zeros both are one. In Table 2, “-” means there is no results with the code.

From Tables 2, we can see that for the examples with more variables and high multiplicity, our method has a better result regardless of the verified accuracy, computing time or the final scale.

We also test the example in [11] with the form: $\{x_{1}^{3}-x_{1}^{2}-x_{2}^{2},x_{2}^{3}+x_{2}^{2}-x_{3},\ldots,x_{n-1}^{3}+x_{n-1}^{2}-x_{n},x_{n}^{2}\}$. The example named breath2 in Table 2 has this form for $n=5$. The method in [9] can compute this example for $n=6$ and it takes 659.59 seconds with the final size for 321 variables and 819 polynomials. We test the cases for $n=6,n=1000$ and $n=2000$ with our code, it takes 0.228965 seconds, 165.274439 seconds and 1036.773847 seconds respectively without introducing new variables.

For our method, although we introduce new variables, the size of our final deflated system is small in experiments. And further, we also compare our method with the other four deflation methods [9] on the following four small systems.

1. 1:

$\{x_{1}^{4}-x_{2}x_{3}x_{4},\ x_{2}^{4}-x_{1}x_{3}x_{4},\ x_{3}^{4}-x_{1}x_{2}x_{4},\ x_{4}^{4}-x_{1}x_{2}x_{3}\}$ at $(0,0,0,0)$ with $\mu=131$; 2. 2:

$\{x^{4},\ x^{2}y+y^{4},z+z^{2}-7x^{3}-8x^{2}\}$ at $(0,0,-1)$ with $\mu=16$; 3. 3:

$\{14x+33y-3\sqrt{5}x^{2}-12\sqrt{5}xy-12\sqrt{5}y^{2}-6\sqrt{5}+x^{3}+6x^{2}y+12xy^{2}+8y^{3}+\sqrt{7},41x-18y-\sqrt{5}+8x^{3}-12x^{2}y+6xy^{2}-y^{3}+12\sqrt{7}xy-12\sqrt{7}x^{2}-3\sqrt{7}y^{2}-6\sqrt{7}\}$ at $\mathbf{p}\approx(1.5055,0.36528)$ with $\mu=5$; 4. 4:

$\{2x_{1}+2x_{1}^{2}+2x_{2}+2x_{2}^{2}+x_{3}^{2}-1,(x_{1}+x_{2}-x_{3}-1)^{3}-x_{1}^{3},(2x_{1}^{3}+5x_{2}^{2}+10x_{3}+5x_{3}^{2}+5)^{3}-1000x_{1}^{5}\}$ at $(0,0,-1)$ with $\mu=18$.

The result (see also in [9]) is below, where method A is in [2, 12], method B is in [8], method C is in [5], method D is in [9], method E is our method $\mathbf{VDSS}$. In Table 3, we denote $Poly$ the number of the polynomials of the final deflation system and $Var$ the number of the variables in the final deflation system. Noting that our final system does not always contain all the polynomials of the input system, therefore, we will contain the number of the different polynomials in the input system, which is not contained in the final system, into $Poly$.

In Table 3, for system 1, 2 and 4, our method matches the best of the other four methods and simultaneously has a smallest deflated system in the five methods. For system 3, although our final system has one more variable than method D, we have less polynomials.

6 Conclusion

In this paper, we develop a new deflation method for refining or verifying the isolated singular zeros of polynomial systems. Given a polynomial system ${\mathbf{F}}\subset{\mathbb{C}}[\mathbf{x}]$ with an isolated singular zero $\mathbf{p}$, by computing the derivatives of the input polynomials directly or the linear combinations of the related polynomials, we prove constructively that there exists a final deflated system $\widetilde{{\mathbf{F}}}^{\prime}(\mathbf{x},\bm{\alpha})$, which has an isolated simple zero $(\mathbf{p},\hat{\bm{\alpha}})$, whose partial projection corresponds to the isolated singular zero $\mathbf{p}$ of the input system ${\mathbf{F}}$. New variables $\bm{\alpha}$ are introduced to represent the coefficients of the linear combinations of the related polynomials to ensure the accuracy of the numerical implementation.

Compared to the previous deflation methods, on one hand, our method also has an output size depending on the depth or the multiplicity of $\mathbf{p}$ in theory. On the other hand, thanks to the acceleration strategies we proposed in the paper, the size of the final system in our actual computations is much less than that we give in theory. The results of the experiments we conduct give a very persuasive argument for this.

In order to essentially have a deeper understanding of our approach, we also give some further analysis of the tolerances $\theta$ and $\varepsilon$ we use. The results of the analysis tells us that our final system is a perturbed system with a bounded perturbation in the worst case. To make our final system as accurate as possible, we also analyse the case that the tolerance $\theta$ is not introduced.

Acknowledgement

The work is partially supported by NSFC Grants 11471327.