=Paper=
{{Paper
|id=Vol-3754/paper15
|storemode=property
|title=A stable computation of multivariarte apporximate GCD based on SVD and lifting technique
|pdfUrl=https://ceur-ws.org/Vol-3754/paper15.pdf
|volume=Vol-3754
|authors=Masaru Sanuki
|dblpUrl=https://dblp.org/rec/conf/sycss/Sanuki24
}}
==A stable computation of multivariarte apporximate GCD based on SVD and lifting technique==
<pdf width="1500px">https://ceur-ws.org/Vol-3754/paper15.pdf</pdf>
<pre>
                         A Stable Computation of Multivariarte Apporximate GCD
                         Based on SVD and Lifting Technique
                         Masaru Sanuki1,*,†
                         1
                             Institute of Medicine, University of Tsukuba, Ten-noudai 1-1-1, Tsukuba-shi, Ibaraki 305-8575, Japan


                                        Abstract
                                        For univariate polynomials, the approximate GCD can be obtained by computing the null space of the subresultant
                                        matrix of given polynomials. In this study, for multivariate polynomials, we propose a method for computing null
                                        space of the subresultant matrix within polynomials stably and efficiently, which is based on the SVD (singular
                                        value decomposition) and lifting techniques. Therefore, we show the multivariate approximate GCD can be also
                                        computed by using subresultant matrix. In addition, we describe an ill-conditioned case (initial factors have
                                        approximate common factor) and solve them.

                                        Keywords
                                        Approximate GCD, Lifting technique, Ill-conditioned cases




                         1. Preliminaries
                         Let 𝐹 (𝑥, 𝑡) and 𝐺(𝑥, 𝑡) be multivariate polynomials in F[𝑥, 𝑡1 , . . . , 𝑡ℓ ] = F[𝑥, 𝑡] (𝑥 is the main variable
                         and 𝑡 = (𝑡1 , . . . , 𝑡ℓ ) are sub-variables), and be expressed as
                                                                    ˜ (𝑥, 𝑡)𝐶(𝑥, 𝑡) + ∆𝐹 = 𝑓𝑚 (𝑡)𝑥𝑚 + . . . + 𝑓0 (𝑡),
                                                         𝐹 (𝑥, 𝑡) = 𝐹
                                                         𝐺(𝑥, 𝑡) = 𝐺˜ (𝑥, 𝑡)𝐶(𝑥, 𝑡) + ∆𝐺 = 𝑔𝑛 (𝑡)𝑥𝑛 + . . . + 𝑔0 (𝑡).

                         Here, 𝐹 ˜,𝐺˜ , 𝐶, ∆𝐹 , ∆𝐺 are polynomials in F[𝑥, 𝑡], and when ||∆𝐹 || ≪ ||𝐹 || and ||∆𝐺 || ≪ ||𝐺||, 𝐶 is
                         called an approximate factor of 𝐹 and 𝐺. In particular, the approximate factor of maximum degree is
                         called approximate GCD, which is denoted by gcd(𝐹, 𝐺).
                            Various algorithm exist for the approximate GCD of univariate polynomials. However, there are
                         few stable all-purpose methods for a large number of variables in a multivariate case. Numerical-
                         based methods are stable but significantly less efficient, so we have tried to improve efficiency by
                         combining lifting methods [4]. In this study, we challenge the stable computation of the null space of
                         the subresultant within polynomial entries.
                            First, we review the method for the subresultant matrix for multivariate polynomials. For the resultant
                         within polynomials, we propose a QRGCD-like method over truncated power-series polynomials, it
                         is efficient [5]. For the null space of the subresultant matrix, Gao et al. and Zeng-Dayton proposed
                         SVD-based methods for numeric matrices at the same conference [2, 7], where the SVD is the singular
                         value decomposition for matrix. These matrices are sparse and the size are also huge extremely although
                         the degree of given polynomial is not large. Lifting techniques is known for solving of equation modulo
                         an ideal 𝐼 and lifting them to solution modulo 𝐼 2 , 𝐼 3 , . . . in order to get the ideal adic completion.
                         Here 𝐼 is an ideal as 𝐼 = ⟨𝑡1 − 𝑠1 , . . . , 𝑡ℓ − 𝑠ℓ ⟩ with (𝑠1 , . . . , 𝑠ℓ ) ∈ Fℓ (in this paper, (𝑠1 , . . . , 𝑠ℓ ) is the
                         origin). For multivariate GCD computation, the EZ-GCD method is well-known lifting method based
                         on Hensel’s lemma, however, its approximate computation will be unstable when initial factors have an
                         approximate common factor [8].
                            In this paper, we propose a stable multivariate approximate GCD computation, which is based on the
                         SVD and lifting techniques. It is able to compute the approximate GCD even though initial factors have
                         approximate common factors.

                          SCSS 2024: 10th International Symposium on Symbolic Computation in Software Science, August 28–30, 2024, Tokyo, Japan
                         *
                           Corresponding author.
                          $ sanuki@md.tsukuba.ac.jp (M. Sanuki)
                                        © 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).


CEUR
                  ceur-ws.org
Workshop      ISSN 1613-0073
Proceedings




                                                                                                               99
2. Framework of algorithm
In this paper, we discuss non-singular case only, i.e., 𝐹 (𝑥, 0) · 𝐺(𝑥, 0) ̸= 0 and 𝑓𝑚 (0) · 𝑔𝑛 (0) ̸= 0. For
non-singular case, every polynomial 𝑃 (𝑥, 𝑡) is transform to 𝑃 (𝑥, 𝑇, 𝑡) = 𝑃 (𝑥, 𝑇 𝑡1 , . . . , 𝑇 𝑡ℓ ), with 𝑇
is the total-degree variable. Every polynomial 𝑃 (𝑥, 𝑇, 𝑡) is represented as the sum of homogeneous
polynomials w.r.t. the total-degree variable 𝑇 ;

                  𝑃 (𝑥, 𝑇, 𝑡) = 𝑃 (0) (𝑥) + 𝑇 · 𝛿𝑃 (1) (𝑥, 𝑡) + . . . + 𝑇 𝑤 · 𝛿𝑃 (𝑤) (𝑥, 𝑡) + . . . ,
               𝑃 (𝑤) (𝑥, 𝑇, 𝑡) = 𝑃 (0) (𝑥) + 𝑇 · 𝛿𝑃 (1) (𝑥, 𝑡) + . . . + 𝑇 𝑤 · 𝛿𝑃 (𝑤) (𝑥, 𝑡).

   In non-singular case, the following two conditions exist: deg(gcd(𝐹, 𝐺)) ≤ deg(gcd(𝐹 (0) , 𝐺(0) ))
and gcd(𝐹, 𝐺)|gcd(𝐹 (0) , 𝐺(0) ). In such situations, lifting algorithms can be applied. The proposed
algorithm is discussed in the next section.

2.1. Computing cofactors of 𝐹 and 𝐺 via lifting method
Let 𝒮𝑖 (𝐹, 𝐺) = 𝒮𝑖 ∈ F[𝑡, 𝑇 ](𝑚+𝑛−2𝑖)×(𝑚+𝑛−2𝑖) be an 𝑖th-subresultant matrix of 𝐹 and 𝐺 w.r.t. 𝑥, and
be represented as
                            ⎛                                                   ⎞
                                 𝑓𝑚                       𝑔𝑛
                                  ..    ..                 ..     ..
                                   .       .                .        .
                            ⎜                                                   ⎟
                            ⎜                                                   ⎟
                            ⎜                                                   ⎟
                   𝒮𝑖 = ⎜ 𝑓𝑚−𝑛−𝑘
                            ⎜                  𝑓𝑚      𝑔𝑛−𝑚−𝑘            𝑔𝑛     ⎟
                                   ..   ..      ..          ..    ..      ..
                                                                                ⎟
                                    .      .     .           .       .     .
                            ⎜                                                   ⎟
                            ⎝                                                   ⎠
                                             𝑓𝑚−𝑛−𝑘                    𝑔𝑛−𝑚−𝑘
                                 (0)               (1)                        (𝑤)
                            = 𝒮𝑖       + 𝑇 · 𝛿𝒮𝑖         + . . . + 𝑇 𝑤 · 𝛿𝒮𝑖        + ...,
         (0)          (0)                                        (𝑤)
where 𝒮𝑖 = 𝛿𝒮𝑖 ∈ F(𝑚+𝑛−2𝑖)×(𝑚+𝑛−2𝑖) and 𝛿𝒮𝑖 ∈ F[𝑡](𝑚+𝑛−2𝑖)×(𝑚+𝑛−2𝑖) for 𝑤 ≥ 1.
  When 𝑘 = deg(gcd(𝐹, 𝐺)) = deg(gcd(𝐹 (0) , 𝐺(0) )), it is well-known as the null space of 𝒮𝑘−1
corresponds to 𝐺
               ˜ and 𝐹
                     ˜ , and rank(𝒮𝑘−1 ) = 𝐾 − 1 where 𝐾 = 𝑚 − (𝑘 − 1) + 𝑛 − (𝑘 − 1).

computation of cofactors for univariate part: SVD
                                                                                             (0)
Cofactors of 𝐹 (0) and 𝐺(0) can be obtained from the null space of 𝒮𝑘−1 . In this paper, we compute the
                  (0)                                                   (0)
null space of 𝒮𝑘−1 using the SVD [1]. Using the SVD of 𝒮𝑘−1 , we obtain the following decomposition:

                                                                                                       𝑣 𝑇1
                                                                       ⎛                           ⎞⎛       ⎞
                                                                           𝜎1
                      (0)                                                           ..             ⎟ ⎜ .. ⎟
                  𝒮𝑘−1 = 𝑈 Σ𝑉 𝑇 =
                                             (︀                     )︀ ⎜
                                                  𝑢1 · · · 𝑢𝐾          ⎝                 .         ⎠⎝ . ⎠,
                                                                                              𝜎𝐾        𝑣 𝑇𝐾

where 𝐾 = 𝑚 − (𝑘 − 1) + 𝑛 − (𝑘 − 1), 𝑈 and 𝑉 are orthogonal matrices, and 𝜎𝑖 are singular vectors
                                                                       (0)
with 𝜎1 ≥ 𝜎2 ≥ . . . ≥ 𝜎𝐾−1 ≫ 𝜎𝐾 ≥ 0, respectively. Then, 𝑣 𝐾 ∈ Ker(𝒮𝑘−1 ), and it is one of the
                 (0
solutions of 𝒮𝑘−1 𝑧 = 0 is 𝑧 = 𝑟 (0) = 𝑣 𝐾 ;

                                                   (0)
                                         ⎛                  ⎞
                                                  𝑔˜𝑛−𝑘
                                     ⎜         ⎟
                                     ⎜         ⎟
                                     ⎜    (0)  ⎟
                                     ⎜ 𝑔˜0
                                               ⎟ , with ||𝑣 𝐾 ||2 = 1.
                                               ⎟
                                𝑣𝐾 = ⎜
                                     ⎜    (0)
                                        ˜
                                       −𝑓
                                               ⎟
                                     ⎜    𝑚−𝑘 ⎟
                                     ⎜         ⎟
                                     ⎝         ⎠
                                           (0)
                                         ˜
                                       −𝑓 0




                                                              100
Computation of cofactors for multivariate part: lifting method
Suppose we have 𝑟 (𝑤) = 𝑟 (𝑤−1) +𝛿𝑟 (𝑤) . Here, 𝛿𝑟 (𝑤) is a vector generated by homogenious polynomials
with total-degree 𝑤 w.r.t. 𝑇 . Then, 𝛿𝑟 (𝑤+1) is generated as follows. Note that the following consists.
                               𝒮𝑘−1 𝑟 (𝑤+1) ≡ 0         (mod 𝑇 𝑤+2 )
                                                       𝑤+1
                               (0)
                                                       ∑︁ (𝑗)
                              𝒮𝑘−1 𝛿𝑟 (𝑤+1)        = −     𝒮𝑘−1 𝛿𝑟 (𝑤−𝑗) = 𝛿𝑝(𝑤) .
                                                         𝑗=1

Now, 𝛿𝑟 (𝑤+1) and 𝛿𝑝(𝑤+1) are transformed bases from 𝑒1 , . . . , 𝑒𝐾 to 𝑣 1 , . . . , 𝑣 𝐾 and 𝑢1 , . . . , 𝑢𝐾 ,
respectively, as follows:
         ⎛ (𝑤) ⎞                                           ⎛ (𝑤) ⎞
             𝛿𝑧1                                                 𝛿𝑝1
   (𝑤)   ⎜ .. ⎟           (𝑤)           (𝑤)        (𝑤)     ⎜ .. ⎟                   (𝑤)                  (𝑤)
𝛿𝑧     = ⎝ . ⎠ = 𝛿𝑧ˆ1 𝑣 1 + . . . + 𝛿𝑧ˆ𝐾 𝑣 𝐾 , 𝛿𝑝      = ⎝ . ⎠ = 𝛿𝑝               ˆ1 𝑢1 + . . . + 𝛿𝑝   ˆ𝑁 𝑢𝐾 .
               (𝑤)                                                                (𝑤)
            𝛿𝑧𝐾                                                                𝛿𝑝𝐾
                     (𝑤+1)           (𝑤+1)
Then, we obtain 𝛿𝑧ˆ𝑖         = 𝛿𝑝
                                ˆ𝑖           /𝜎𝑖   (𝑖 = 1, . . . , 𝐾 − 1). Therefore, 𝛿𝑟 (𝑤+1) is constructed, as
follows.
                              (1)                       (1)
              𝛿𝑟 (𝑤+1) = 𝛿𝑝                     ˆ𝐾−1 /𝜎𝐾−1 𝑣 𝐾−1 + F[𝑇, 𝑡]𝑤+1 · 𝑣 𝐾 ,
                          ˆ1 /𝜎1 𝑣 1 + . . . + 𝛿𝑝
where F[𝑇, 𝑡]𝑤+1 is homogeneous polynomial set with total-degree 𝑤 + 1 w.r.t. 𝑇 , and we have the
following as a candidate solution.
                        𝑤                              𝑤
                               (𝑤)                              (𝑤)
                       ∑︁                             ∑︁
        𝑟 (𝑤) = 𝑡𝐾 +         𝛿𝑝
                              ˆ1 /𝜎1 𝑣 1 + . . . +             ˆ𝐾−1 /𝜎𝐾−1 𝑣 𝐾−1 + F[𝑇, 𝑡][1,𝑤] · 𝑣 𝐾 ,
                                                              𝛿𝑝
                       𝑗=1                             𝑗=1

where F[𝑇, 𝑡][1,𝑤] = ∪𝑤
                      𝑗=1 F[𝑇, 𝑡]𝑗 . To compute the approximate GCD of 𝐹 and 𝐺, we need to determine
𝑞 (𝑤) = 𝛿𝑞 + . . . + 𝛿𝑞(𝑤) ∈ F[𝑇, 𝑡][1,𝑤] s.t.
          (1)

                               𝑤                                𝑤
                                      (𝑤)                              (𝑤)
                              ∑︁                               ∑︁
          𝑟 (𝑤) (𝑞) = 𝑡𝐾 +          𝛿𝑝
                                     ˆ1 /𝜎1 𝑣 1 + . . . +             ˆ𝐾−1 /𝜎𝐾−1 𝑣 𝐾−1 + 𝑞 (𝑤) · 𝑣 𝐾 ,
                                                                     𝛿𝑝
                              𝑗=1                              𝑗=1

  To determine the approximate GCD, we must determine 𝑞 (𝑖) . The following example shows one
approach to determining each undetermined element 𝛿𝑞 (𝑖) for 1 ≤ 𝑖 ≤ 𝑤..

Example1
Polynomials 𝐹 (𝑥, 𝑡1 , 𝑡2 ) and 𝐺(𝑥, 𝑡1 , 𝑡2 ) having an approximate GCD 𝐶(𝑥, 𝑡1 , 𝑡2 ) = 𝑥3 + (1 + 𝑡2 −
2𝑡1 + 𝑡1 2 )𝑥 + 3 are expressed as
              𝐹 (𝑥, 𝑡1 , 𝑡2 ) = (𝑥3 + (𝑡2 2 + 𝑡1 + 𝑡1 − 2)𝑥2 − 1) × 𝐶(𝑥, 𝑡1 , 𝑡2 ) + 𝑀𝑒𝑝𝑠 ,
              𝐺(𝑥, 𝑡1 , 𝑡2 ) = (𝑥3 + (2𝑡2 2 − 𝑡1 + 3)𝑥2 − 1) × 𝐶(𝑥, 𝑡1 , 𝑡2 ) + 𝑀𝑒𝑝𝑠 ,
where 𝑀𝑒𝑝𝑠 is the machine epsilon.
                                                                          (0)
  In this example, 𝑘 = 2 is already known (𝐾 = 8). Then, one solution of 𝒮1 𝑧 = 0 is 𝑧 = 𝑟 (0) = 𝑣 8 ;
                                           −0.242535625036333
                                     ⎛                               ⎞
                                     ⎜     −0.727606875108999        ⎟
                                     ⎜ −2.24840273230668 × 10−15 ⎟
                                     ⎜                               ⎟
                                     ⎜                               ⎟
                        (0)
                                     ⎜      0.242535625036333        ⎟
                       𝑟 = 𝑣8 = ⎜    ⎜                               ⎟.
                                     ⎜      0.242535625036333        ⎟
                                                                     ⎟
                                     ⎜
                                     ⎜     −0.485071250072665        ⎟
                                                                     ⎟
                                     ⎝ −1.32375311946987 × 10−15 ⎠
                                           −0.242535625036333




                                                          101
                                     (1)                   (1)
A candidate of 𝛿𝑟 (1) |𝛿𝑞(1) =0 = 𝛿𝑝                     ˆ7 /𝜎7 𝑣 7 + 𝛿𝑞 (1) × 𝑣 8 is
                                   ˆ1 /𝜎1 𝑣 1 + . . . + 𝛿𝑝
                                                                                           ⎛        (1)     ⎞
                                                                                                  𝛿𝑔˜𝑛−𝑘
                 0.0713340073 · · · 𝑡1 + 0.0285336029 · · · 𝑡2
         ⎛                                                               ⎞
                                                                                                    (1)
                                                                                                𝛿𝑔˜𝑛−𝑘−1
                                                                                        ⎜             ⎟
                                                                                        ⎜             ⎟
         ⎜      −0.0285336029 · · · 𝑡1 + 0.0856008088 · · · 𝑡2       ⎟                  ⎜             ..
                                                                                                      ⎟
         ⎜ 3.65419500 · · · × 10−15 𝑡1 − 4.96824803 · · · × 10−15 𝑡2 ⎟                                 .
         ⎜                                                           ⎟                  ⎜             ⎟
                                                                                        ⎜             ⎟
         ⎜                                                           ⎟                  ⎜       (1)   ⎟
                −0.0713340073 · · · 𝑡1 − 0.0285336029 · · · 𝑡2                              𝛿𝑔˜0
𝛿𝑧 (1) = ⎜                                                           ⎟ + 𝛿𝑞 (1) · 𝑣 8 = ⎜
         ⎜                                                           ⎟                  ⎜             ⎟
                                                                                                      ⎟.
                −0.0713340073 · · · 𝑡1 − 0.0285336029 · · · 𝑡2                                   (1)
                                                                                        ⎜ −𝛿𝑓𝑚−𝑘 ⎟
         ⎜                                                           ⎟                  ⎜             ⎟
         ⎜                                                           ⎟
         ⎜
         ⎜      −0.0998676103 · · · 𝑡1 − 0.185468419 · · · 𝑡2        ⎟
                                                                     ⎟
                                                                                        ⎜     (1)
                                                                                        ⎜ −𝛿𝑓𝑚−𝑘−1
                                                                                                      ⎟
                                                                                                      ⎟
         ⎝ 4.77577504 · · · × 10−16 𝑡1 − 1.10469359 · · · × 10−15 𝑡2 ⎠
                                                                                               ..
                                                                                        ⎜             ⎟
                                                                                                .
                                                                                        ⎜             ⎟
                0.0713340073 · · · 𝑡1 + 0.0285336029 · · · 𝑡2                           ⎝             ⎠
                                                                                                  (1)
                                                                                           −𝛿𝑓0

Generally, it is difficult to determine 𝛿𝑞 (1) properly.
   However, assuming that cofactors are also not dense or the approximate GCD is monic, sev-
                                                                                             (1)
eral coefficients will be zero. In this example, assume the 1st element is zero, 𝛿𝑞 (1) is 𝛿𝑞1 =
(0.0713340073636269 𝑡1 + 0.0285336029454512 𝑡2 )/0.242535625036333 and 𝛿𝑟 (1) becomes

                                           0
                     ⎛                                            ⎞
                     ⎜
                     ⎜            0.242535625036331𝑡1             ⎟
                                                                  ⎟
                     ⎜
                     ⎜                     0                      ⎟
                                                                  ⎟
                     ⎜
                     ⎜                     0                      ⎟
                                                                  ⎟.
                     ⎜
                     ⎜                     0                      ⎟
                                                                  ⎟
                     ⎜ −0.242535625036331𝑡1 − 0.242535625036334𝑡2 ⎟
                     ⎜                                            ⎟
                     ⎝                     0                      ⎠
                                           0

It is unlikely that many factors will be close to zero simultaneously, and this can only happen if the result
is correct. Unlike the EZ-GCD method, it is more efficient because it can extract each undetermined
coefficient at each lifting step. So that, “check zeros” is very efficiency.
    If the coefficients are dense, lc(lc(𝐹 ), lc(𝐺)) or lc(lc(gcd(𝐹, 𝐺))) should be calculated in advance so
that the elements can be determined uniquely.

2.2. Computing approximate GCD
After obtaining cofactors, the approximate GCD is computed by solving and 𝐹                                =
(𝑓𝑚 , 𝑓𝑚−1 , . . . , 𝑓0 )𝑇 ∈ F[𝑡]𝑚+1 .


                           𝑓˜ 𝑚−𝑘
                         ⎛                       ⎞
                               ..
                                                    ⎛      ⎞ ⎛      ⎞
                                      ..               𝑐𝑘        𝑓𝑚
                                .        .
                         ⎜                        ⎟
                                                  ⎟ ⎜ 𝑐𝑘−1 ⎟
                                                           ⎟ ⎜ 𝑓𝑚−1 ⎟
                         ⎜                        ⎟⎜           ⎜    ⎟
                         ⎜ 𝑓˜              ˜
                         ⎜
                                 0         𝑓 𝑚−𝑘 ⎟      .    =
                                                  ⎟⎝ . ⎠ ⎝ . ⎟    .
                                                        .         . ⎠
                                                    ⎜      ⎟   ⎜
                                      ..       ..
                         ⎜
                                         .      .
                         ⎜                        ⎟
                                                       𝑐0        𝑓0
                         ⎝                        ⎠
                                             𝑓˜
                                             0

This linear equation is solve as following step.
              (0)
   1. Solve 𝒞𝑚+1,𝑘+1 (𝐹 ˜ ) · 𝑐(0) = 𝐹 (0) . Actually, we utilize the SVD as in the former case.
                            (0)       ˜ ) · 𝛿𝑐(𝑤) = 𝛿𝐹 (𝑤) − ∑︀ 𝐶 (𝑖)
   2. Lifting step: solve 𝒞𝑚+1,𝑘+1 (𝐹                                          ˜
                                                                   𝑖 𝑚+1,𝑘+1 (𝐹 ) × 𝛿𝑐
                                                                                       (𝑤−𝑖) . This step can

      also be solved using SVD.
   3. Return 𝑐𝑘 𝑥𝑘 + · · · + 𝑐0 as an approximate GCD.




                                                     102
3. Solve in ill-conditioned cases
In this section, we demonstrate that our method is stable for ill-conditioned cases [8, 5]. On the other
word, we deals with cases where the initial factor is an approximate common factor. In this case, the
EZ-GCD method is unstable since large cancellation errors occur [6].

Example 2 (initial factors have approximate common factor)
Compute the approximate GCD of 𝐹 and 𝐺, where both polynomials are monic.

       𝐹 (𝑥, 𝑡1 , 𝑡2 ) = (𝑥3 + (𝑡2 2 + 𝑡1 + 𝑡2 − 2)𝑥2 − 1)(𝑥 − 1.0003 + 2𝑡2 − 𝑡1 2 )𝐶 + 𝑀𝑒𝑝𝑠 ,
       𝐺(𝑥, 𝑡1 , 𝑡2 ) = (𝑥3 + (2𝑡2 2 − 𝑡1 + 3)𝑥2 − 1)(𝑥 − 1.0005 + 𝑡1 + 𝑡2 + 𝑡1 𝑡2 )𝐶 + 𝑀𝑒𝑝𝑠 ,
       𝐶(𝑥, 𝑡1 , 𝑡2 ) = 𝑥3 + (1 + 𝑡2 − 2𝑡1 + 𝑡1 2 )𝑥 + 3.

Initial factors 𝐹 (0) and 𝐺(0) have an approximate common factor (𝑥 − 1.0002) with tolerance 𝑂(10−5 ).
                     (0)
Sigular values of 𝒮3 (𝐹, 𝐺) are 19.8 > 18.3 > 14.5 > 12.8 > 8.2 > 4.4 > 1.1 > 0.6 ≫ 1.5×10−5 ≫
1.1×10−16 . Because give polynomials are monic, the leading coefficient of cofactors and the approximate
GCD are also monic, respectively.
   When 𝑤 = 1, Adjusting the 1st element of 𝛿𝑟 (1) by 𝑧 𝐾 only, we obtained the following.

                                             0.
                   ⎛                                                        ⎞
               ⎜
               ⎜    −7.25814180424500 × 10−8 𝑡1 + 0.176741414005228𝑡2       ⎟
                                                                            ⎟
               ⎜        0.707053518826616𝑡1 + 0.530224242015704𝑡2           ⎟
               ⎜ −6.96526170074208 × 10−14 𝑡1 + 6.29774010718620 × 10−14 𝑡2 ⎟
               ⎜                                                            ⎟
               ⎜                                                            ⎟
           (1)
               ⎜       −0.176741268890038𝑡1 − 0.176741414005196𝑡2           ⎟
         𝛿𝑟 = ⎜⎜
                                      −15                           −16
                                                                            ⎟
               ⎜ 1.42941214420489 × 10 𝑡1 + 6.38378239159465 × 10 𝑡2 ⎟
                                                                            ⎟
               ⎜
               ⎜       −0.176741268889989𝑡1 − 0.530224096948041𝑡2           ⎟
                                                                            ⎟
               ⎜
               ⎜        0.176794218725546𝑡 1 +  0.883759874841642𝑡2
                                                                            ⎟
                                                                            ⎟
                                       −14
               ⎝ −1.18932641512970 × 10 𝑡 − 7.57727214306669 × 10 𝑡 ⎠−15
                                                    1                                 2
                         −7.25814328778052 × 10−8 𝑡1 + 0.353482755476628𝑡2

                      ˜ (1) 𝐺(1) − 𝐹
The perturbation is ||𝐹            ˜ (1) 𝐺(1) || ≈ 𝑂(10−8 ) ≈ 𝜎𝐾 /𝜎𝐾−1 . On the other hand, by adjusting
𝑧 𝐾−1 , we obtained the following, It can be confirmed that the solution is not accurate; ||𝐹˜ (1) 𝐺(1) −
˜ (1) 𝐺(1) || ≈ 𝜎𝐾 .
𝐹

                                             0.
              ⎛                                                          ⎞
              ⎜
              ⎜       −0.1988511825793107𝑡1 − 0.14357803747786974𝑡2      ⎟
                                                                         ⎟
              ⎜
              ⎜       0.11055165805122452𝑡 1 −  0.43065120320683287𝑡2    ⎟
                                                                         ⎟
              ⎜
              ⎜    0.0002976768351589665𝑡 1  + 0.00047951294108034004𝑡 2 ⎟
                                                                         ⎟
              ⎜      0.022318043342197766𝑡1 + 0.14391342019466513𝑡2      ⎟
                                        −5
              ⎜                                                          ⎟
              ⎜ −0.7183155936049679 × 10 𝑡1 − 0.000011570991832604571𝑡2 ⎟
              ⎜                                                          ⎟
              ⎜
              ⎜       0.02212670015656562𝑡 1 −  0.20987748805445672𝑡2
                                                                         ⎟
                                                                         ⎟
              ⎜
              ⎜       −0.22096310056080598𝑡1 + 0.243032215144183𝑡2       ⎟
                                                                         ⎟
              ⎝   0.00007010876634662433𝑡1 + 0.00011293475597320968𝑡2    ⎠
                     −0.19880976332099576𝑡1 + 0.03323002423516758𝑡2

  When 𝑤 = 2, by lifting step and adjusting the 1st element of 𝛿𝑟 (2) by 𝑧 𝐾 , we obtained the following




                                                  103
vector.
                                              0.
⎛                                                                                                 ⎞
⎜     0.353120013467623𝑡2 2 + 0.177466845429488𝑡1 𝑡2 − 0.000362979823316206𝑡1 2                   ⎟
⎜
⎜                                2
       −0.354747432749536𝑡2 + 0.355659122269873𝑡1 𝑡2 − 0.177830208344029𝑡1               2        ⎟
                                                                                                  ⎟
                       −13   2                         −12                              −12   2
⎜                                                                                                 ⎟
⎜ 8.84819995 · · · × 10 𝑡2 − 4.59634934 · · · × 10 𝑡1 𝑡2 + 1.03052288 · · · × 10 𝑡1 ⎟
                                2 − 0.177466845422696𝑡 𝑡 + 0.000362979818887450𝑡 2
⎜                                                                                                 ⎟
⎜
⎜   0.000362669479650586𝑡     2                          1  2                              1      ⎟
                                                                                                  ⎟
⎜ −9.08967345 · · · × 10−13 𝑡2 2 − 3.63698827 · · · × 10−12 𝑡1 𝑡2 − 1.36436695 · · · × 10−12 𝑡1 2 ⎟
⎜                                                                                                 ⎟
⎜
⎜    −0.176378671991595𝑡2 2 − 0.000725503949587463𝑡1 𝑡2 + 0.177104321289445𝑡1 2                   ⎟
                                                                                                  ⎟
⎜                                2
       −0.177413730546595𝑡2 − 0.352031674994445𝑡1 𝑡2 − 0.354208569999507𝑡1               2        ⎟
⎜                                                                                                 ⎟
⎝ −9.15635622 · · · × 10−13 𝑡 2 + 2.71640175 · · · × 10−12 𝑡 𝑡 − 1.68760838 · · · × 10−13 𝑡 2 ⎠
                                2                               1 2                             1
     −0.000362669478412958𝑡2 2 + 0.000725503952765407𝑡1 𝑡2 − 0.177104321291319𝑡1 2

                                                                                 (0)
Adjusting only 𝑣 𝐾 is not accurate. Therefore, adjusting 𝑣 𝐾 and 𝑣 𝐾−1 in ker 𝒮𝑘 we have the following,
                                                                                ˜ (2) 𝐺(2) − 𝐹
and it obtains the expected solution one . In this case, perturbation becomes ||𝐹            ˜ (2) 𝐺(2) || ≈
𝜎𝐾 , it is better.
  The SVD is stable even if the matrix is irregular. Thus, the SVD of 𝒮 (0) is also stable even if initial
factors have an approximate common factor. On the other hand, a lifting method using the Bezout
matrix is unstable since initial matrix is assumed to be regular [4]. Hence, our method is more stable
and efficient than existing methods.


References
 [1] R. Corless, P. Gianni, B. Trager and S. Watt, The singular value decomposition for polynomial systems,
     Proc. of ISSAC’95, ACM Press, 1995, 195–207.
 [2] S. Gao, E. Kaltofen, J. P. May, Z. Yang and L. Zhi, Approximate factorization of multivariate
     polynomials via differential equations, Proc. of ISSAC’04, ACM Press, 2004, 167–174.
 [3] M. Ochi, M-T. Noda and T. Sasaki, Approximate greatest common divisor of multivariate polynomials
     and its application to ill-conditioned systems of algebraic equations, J. Inform. Proces., 14 (1991),
     292–300.
 [4] M. Sanuki. Computing multivariate approximate GCD based on Barnett’s theorem, Proc. of Symbolic-
     Numeric Computation 2009 (SNC 2009), 2009, 149–157.
 [5] M. Sanuki and T. Sasaki, Computing approximate GCDs in ill-conditioned cases, International
     Workshop of Symbolic-Numeric Computation 2007 (SNC2007), ACM Press, 2007, 170-179, 25–27.
 [6] T. Sasaki and S. Yamaguchi, An analysis of cancellation error in multivariate Hensel construction
     with floating-point number arithmetic, Proc. of ISSAC’98, ACM Press, 1998, 1–8.
 [7] Z. Zeng and B. H. Dayton, The approximate GCD of inexact polynomials part II: A multivariate
     algorithm, Proc. of ISSAC’04, ACM Press, 2004, 320–327.
 [8] L. Zhi and M-T. Noda, Approximate GCD of Multivariate Polynomials, Proc. of ASCM2000, World
     Scientific, 2000, 9–18.




                                                    104

</pre>