From 089d03b5fb88185726d8c20d540d6dca2a36e11f Mon Sep 17 00:00:00 2001 From: john verzani Date: Fri, 15 Jul 2022 15:06:47 -0400 Subject: [PATCH] Ngcd cleanup (#430) * refactor `ngcd` code * clean up doc strings * formatting and minor adjustments * use `muladd` where possible * version bump --- Project.toml | 2 +- src/Polynomials.jl | 1 - src/polynomials/multroot.jl | 109 ++-- src/polynomials/ngcd.jl | 837 ++++++++++++++++++----------- src/polynomials/pi_n_polynomial.jl | 5 +- src/polynomials/standard-basis.jl | 4 +- test/StandardBasis.jl | 20 +- 7 files changed, 598 insertions(+), 380 deletions(-) diff --git a/Project.toml b/Project.toml index a6b43055..0f47521e 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "Polynomials" uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" license = "MIT" author = "JuliaMath" -version = "3.1.4" +version = "3.1.5" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/Polynomials.jl b/src/Polynomials.jl index 897db5ae..79d82bea 100644 --- a/src/Polynomials.jl +++ b/src/Polynomials.jl @@ -23,7 +23,6 @@ include("polynomials/pi_n_polynomial.jl") include("polynomials/factored_polynomial.jl") include("polynomials/ngcd.jl") include("polynomials/multroot.jl") - include("polynomials/ChebyshevT.jl") # Rational functions diff --git a/src/polynomials/multroot.jl b/src/polynomials/multroot.jl index 3e45ee80..13cd2642 100644 --- a/src/polynomials/multroot.jl +++ b/src/polynomials/multroot.jl @@ -55,10 +55,11 @@ multiplicity structure: * `θ`: the singular value threshold, set to `1e-8`. This is used by `Polynomials.ngcd` to assess if a matrix is rank deficient by - comparing the smallest singular value to `θ ⋅ ||p||₂`. + comparing the smallest singular value to `θ ⋅ ‖p‖₂`. * `ρ`: the initial residual tolerance, set to `1e-10`. This is passed - to `Polynomials.ngcd`, the GCD finding algorithm as a relative tolerance. + to `Polynomials.ngcd`, the GCD finding algorithm as a measure for + the absolute tolerance multiplied by `‖p‖₂`. * `ϕ`: A scale factor, set to `100`. As the `ngcd` algorithm is called recursively, this allows the residual tolerance to scale up to match @@ -69,11 +70,11 @@ Returns a named tuple with * `values`: the identified roots * `multiplicities`: the corresponding multiplicities * `κ`: the estimated condition number -* `ϵ`: the backward error, `||p̃ - p||_w`. +* `ϵ`: the backward error, `‖p̃ - p‖_w`. If the wrong multiplicity structure is identified in step 1, then typically either `κ` or `ϵ` will be large. The estimated forward -error, `||z̃ - z||₂`, is bounded up to higher order terms by `κ ⋅ ϵ`. +error, `‖z̃ - z‖₂`, is bounded up to higher order terms by `κ ⋅ ϵ`. This will be displayed if `verbose=true` is specified. For polynomials of degree 20 or higher, it is often the case the `l` @@ -100,62 +101,85 @@ function multroot(p::Polynomials.StandardBasisPolynomial{T}; verbose=false, end -# The multiplicity structure, `l`, gives rise to a pejorative manifold `Πₗ = {Gₗ(z) | z∈ Cᵐ}`. -function pejorative_manifold(p::Polynomials.StandardBasisPolynomial{T}; - ρ = 1e-10, # initial residual tolerance - θ = 1e-8, # zero singular-value threshold - ϕ = 100) where {T} +# The multiplicity structure, `l`, gives rise to a pejorative manifold +# `Πₗ = {Gₗ(z) | z ∈ Cᵐ}`. This first finds the approximate roots, then +# finds the multiplicity structure +function pejorative_manifold(p::Polynomials.StandardBasisPolynomial{T,X}; + θ = 1e-8, # zero singular-value threshold + ρ = 1e-10, # initial residual tolerance + ϕ = 100 # residual growth factor + ) where {T,X} + + S = float(T) + u = Polynomials.PnPolynomial{S,X}(S.(coeffs(p))) + + nu₂ = norm(u, 2) + θ′, ρ′ = θ * nu₂, ρ * nu₂ + + u, v, w, ρⱼ, κ = Polynomials.ngcd(u, derivative(u); + satol = θ′, srtol = zero(real(T)), + atol = ρ′, rtol = zero(real(T)) + ) + ρⱼ /= nu₂ - zT = zero(float(real(T))) - u, v, w, θ′, κ = Polynomials.ngcd(p, derivative(p), - atol=ρ*norm(p), satol = θ*norm(p), - rtol = zT, srtol = zT) zs = roots(v) + ls = pejorative_manifold_multiplicities(u, + zs, ρⱼ, + θ, ρ, ϕ) + zs, ls +end + + + +function pejorative_manifold_multiplicities(u::Polynomials.PnPolynomial{T}, + zs, ρⱼ, + θ, ρ, ϕ) where {T} nrts = length(zs) ls = ones(Int, nrts) while !Polynomials.isconstant(u) - normp = 1 + norm(u, 2) - ρ′ = max(ρ*normp, ϕ * θ′) # paper uses just latter - u, v, w, θ′, κ = Polynomials.ngcd(u, derivative(u), - atol=ρ′, satol = θ*normp, - rtol = zT, srtol = zT, - minⱼ = degree(u) - nrts # truncate search + nu₂ = norm(u,2) + θ = θ * nu₂ + ρ = max(ρ, ϕ * ρⱼ) + ρ′ = ρ * nu₂ + u, v, w, ρⱼ, κ = Polynomials.ngcd(u, derivative(u), + satol = θ, srtol = zero(real(T)), + atol = ρ′, rtol = zero(real(T)), + minⱼ = degree(u) - nrts ) - - rts = roots(v) - nrts = length(rts) - - for z in rts - tmp, ind = findmin(abs.(zs .- z)) - ls[ind] = ls[ind] + 1 + ρⱼ /= nu₂ + nrts = degree(v) + for z ∈ roots(v) + _, ind = findmin(abs.(zs .- z)) + ls[ind] += 1 end end - zs, ls # estimate for roots, multiplicities + ls end + """ pejorative_root(p, zs, ls; kwargs...) Find a *pejorative* *root* for `p` given multiplicity structure `ls` and initial guess `zs`. The pejorative manifold for a multplicity structure `l` is denoted `{Gₗ(z) | z ∈ Cᵐ}`. A pejorative -root is a least squares minimizer of `F(z) = W ⋅ [Gₗ(z) - a]`. Here `a ~ (p_{n-1}, p_{n-2}, …, p_1, p_0) / p_n` and `W` are weights given by `min(1, 1/|aᵢ|)`. When `l` is the mathematically correct structure, then `F` will be `0` for a pejorative root. If `l` is not correct, then the backward error `||p̃ - p||_w` is typically large, where `p̃ = Π (x-z̃ᵢ)ˡⁱ`. +root is a least squares minimizer of `F(z) = W ⋅ [Gₗ(z) - a]`. Here `a ~ (p_{n-1}, p_{n-2}, …, p_1, p_0) / p_n` and `W` are weights given by `min(1, 1/|aᵢ|)`. When `l` is the mathematically correct structure, then `F` will be `0` for a pejorative root. If `l` is not correct, then the backward error `‖p̃ - p‖_w` is typically large, where `p̃ = Π (x-z̃ᵢ)ˡⁱ`. This follows Algorithm 1 of [Zeng](https://www.ams.org/journals/mcom/2005-74-250/S0025-5718-04-01692-8/S0025-5718-04-01692-8.pdf) """ function pejorative_root(p::Polynomials.StandardBasisPolynomial, - zs::Vector{S}, ls::Vector{Int}; kwargs...) where {S} + zs::Vector{S}, ls; kwargs...) where {S} ps = (reverse ∘ coeffs)(p) pejorative_root(ps, zs, ls; kwargs...) end # p is [1, a_{n-1}, a_{n-2}, ..., a_1, a_0] -function pejorative_root(p, zs::Vector{S}, ls::Vector{Int}; +function pejorative_root(p, zs::Vector{S}, ls; τ = eps(float(real(S)))^(2/3), maxsteps = 100, # linear or quadratic verbose=false @@ -170,7 +194,7 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int}; a = p[2:end]./p[1] # a ~ (p[n-1], p[n-2], ..., p[0])/p[n] W = Diagonal([min(1, 1/abs(aᵢ)) for aᵢ in a]) J = zeros(S, m, n) - G = zeros(S, 1 + m) + G = zeros(S, 1 + m) Δₖ = zeros(S, n) zₖs = copy(zs) # updated @@ -205,10 +229,10 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int}; if cvg return zₖs else - @info (""" + @info """ The multiplicity count may be in error: the initial guess for the roots failed to converge to a pejorative root. -""") +""" return(zₖs) end @@ -220,13 +244,11 @@ function evalG!(G, zs::Vector{T}, ls::Vector) where {T} G .= zero(T) G[1] = one(T) - ctr = 1 for (z,l) in zip(zs, ls) for _ in 1:l - for j in length(G)-1:-1:1#ctr + for j in length(G)-1:-1:1 G[1 + j] -= z * G[j] end - ctr += 1 end end @@ -241,6 +263,7 @@ function evalJ!(J, zs::Vector{T}, ls::Vector) where {T} n, m = sum(ls), length(zs) evalG!(view(J, 1:1+n-m, 1), zs, ls .- 1) + for (j′, lⱼ) in enumerate(reverse(ls)) for i in 1+n-m:-1:1 J[i, m - j′ + 1] = -lⱼ * J[i, 1] @@ -260,19 +283,21 @@ end # Defn (3.5) condition number of z with respect to the multiplicity l and weight W cond_zl(p::AbstractPolynomial, zs::Vector{S}, ls) where {S} = cond_zl(reverse(coeffs(p)), zs, ls) + function cond_zl(p, zs::Vector{S}, ls) where {S} J = zeros(S, sum(ls), length(zs)) - W = diagm(0 => [min(1, 1/abs(aᵢ)) for aᵢ in p[2:end]]) evalJ!(J, zs, ls) - F = qr(W*J) - _, σ, _ = Polynomials.NGCD.smallest_singular_value(F.R) + W = diagm(0 => [min(1, 1/abs(aᵢ)) for aᵢ in p[2:end]]) + + σ = Polynomials.NGCD.smallest_singular_value(W*J) 1 / σ end backward_error(p::AbstractPolynomial, zs::Vector{S}, ls) where {S} = backward_error(reverse(coeffs(p)), zs, ls) + function backward_error(p, z̃s::Vector{S}, ls) where {S} - # || ã - a ||w + # ‖ ã - a ‖w G = zeros(S, 1 + sum(ls)) evalG!(G, z̃s, ls) a = p[2:end]./p[1] @@ -288,8 +313,8 @@ end function show_stats(κ, ϵ) println("") println("Condition number κ(z; l) = ", κ) - println("Backward error ||ã - a||w = ", ϵ) - println("Estimated forward root error ||z̃ - z||₂ = ", κ * ϵ) + println("Backward error ‖ã - a‖w = ", ϵ) + println("Estimated forward root error ‖z̃ - z‖₂ = ", κ * ϵ) println("") end diff --git a/src/polynomials/ngcd.jl b/src/polynomials/ngcd.jl index 723fdefe..ec48949a 100644 --- a/src/polynomials/ngcd.jl +++ b/src/polynomials/ngcd.jl @@ -1,23 +1,24 @@ """ - ngcd(p,q, [k]; kwargs...) + ngcd(p, q, [k]; kwargs...) Find numerical GCD of polynomials `p` and `q`. Refer to [`NGCD.ngcd(p,q)`](@ref) for details. +The main entry point for this function is `gcd(p, q, method=:numerical)`, but `ngcd` outputs the gcd factorization -- `u, v, w` with `u*v ≈ p` and `u*w ≈ q` -- along with `Θ`, an estimate on how close `p`,`q` is to a gcd factorization of degree `k` and `κ` the GCD condition number. In the case `degree(p) ≫ degree(q)`, a heuristic is employed to first call one step of the Euclidean gcd approach, and then call `ngcd` with relaxed tolerances. """ function ngcd(p::P, q::Q, - args...; kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X}, + args...; + kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X}, S,Y,Q<:StandardBasisPolynomial{S,Y}} - if (degree(q) > degree(p)) u,w,v,Θ,κ = ngcd(q,p,args...;kwargs...) return (u=u,v=v,w=w, Θ=Θ, κ=κ) end if degree(p) > 5*(1+degree(q)) a,b = divrem(p,q) - return ngcd(q,b, args...; λ=100, kwargs...) + return ngcd(q, b, args...; λ=100, kwargs...) end # easy cases @@ -48,12 +49,14 @@ function ngcd(p::P, q::Q, q′ = PnPolynomial{R,X}(qs[nz:end]) out = NGCD.ngcd(p′, q′, args...; kwargs...) + ## convert to original polynomial type 𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X} u,v,w = convert.(𝑷, (out.u,out.v,out.w)) if nz > 1 u *= variable(u)^(nz-1) end - (u=u,v=v,w=w, Θ=out.Θ, κ = out.κ) + + (u = u, v = v, w = w, Θ = out.Θ, κ = out.κ) end @@ -64,6 +67,26 @@ Use `ngcd` to identify the square-free part of the polynomial `p`. """ square_free(p) = ngcd(p, derivative(p)).v +""" + rank_reveal(A; atol, rtol) + +This is the high rank-revealing algorithm of [Lee, Li, and Zeng](http://archive.ymsc.tsinghua.edu.cn/pacm_download/285/8743-RankRev_paper.pdf) DOI: DOI 10.1007/s11075-017-0328-7. +""" +function rank_reveal(A::AbstractMatrix{T}; kwargs...) where {T <: AbstractFloat} + m, n = size(A) + @assert m ≥ n + q = qr(A) + + R = UpperTriangular(q.R) + w = Vector{eltype(A)}(undef, n) + λ = norm(A, Inf) + τ = norm(R, Inf) + NGCD.rank_reveal!(R, w, λ, τ; kwargs...) + +end + + + ## ---- the work is done in this module module NGCD @@ -71,74 +94,125 @@ using Polynomials, LinearAlgebra import Polynomials: PnPolynomial, constructorof """ - ngcd(ps::PnPolynomial{T,X}, qs::PnPolynomial{T,X}, [k::Int]; scale::Bool=false, atol=eps(T), rtol=eps(T), satol=atol, srtol=rtol) - -Return `u, v, w, Θ, κ` where ``u⋅v ≈ ps`` and ``u⋅w ≈ qs`` (polynomial multiplication); `Θ` (`\\Theta[tab]`) is the residual error (``‖ [u⋅v,u⋅w] - [ps,qs] ‖``); and `κ` (`\\kappa[tab]`) is the numerical gcd condition number estimate. When `scale=true`, ``u⋅v ≈ ps/‖ps‖`` and ``u⋅w ≈ qs/‖qs‖`` + ngcd(p::PnPolynomial{T,X}, q::PnPolynomial{T,X}, [k::Int]; + atol = eps(real(T))^(5/6), # residual over Πₖ + rtol = eps(real(T)), + satol = atol, + srtol = rtol, + λ=one(real(T)), + scale::Bool=false + ) + +Computes numerical GCD of polynomials `p` and `q`. + +Returns ``u, v, w, Θ, κ`` where ``u⋅v ≈ p`` and ``u⋅w ≈ q`` (polynomial +multiplication); ``Θ`` (`\\Theta[tab]`) is the residual error (``‖ +[u⋅v,u⋅w] - [p,q] ‖``); and ``κ`` (`\\kappa[tab]`) is the numerical gcd +condition number estimate. When `scale=true`, ``u⋅v ≈ ps/‖ps‖₂`` and +``u⋅w ≈ qs/‖qs‖₂``. The numerical GCD problem is defined in [1] (5.4). Let ``(p,q)`` be a -polynomial pair with degree m,n. Let Ρmn be set of all such pairs. Any -given pair of polynomials has an exact greatest common divisor, ``u``, of -degree ``k``, defined up to constant factors. Let ``Ρᵏmn`` be the manifold of -all such ``(p,q)`` pairs with exact gcd of degree ``k``. A given pair ``(p,q)`` with exact gcd of degree ``j`` will -have some distance ``Θᵏ`` from ``Pᵏ``. For a given threshold ``ϵ > 0`` a numerical GCD -of ``(p,q)`` within ``ϵ`` is an exact GCD of a pair ``(p̂,q̂)`` in ``Ρᵏ`` with - -``‖ (p,q) - (p̂,q̂) ‖ <= Θᵏ``, where ``k`` is the largest value for which ``Θᵏ < ϵ``. +polynomial pair with degree ``m``, ``n``. Let ``Ρₘₙ`` be set of all +such pairs. Any given pair of polynomials has an exact greatest common +divisor, ``u``, of degree ``k``, defined up to constant factors. Let +``Ρᵏₘₙ`` be the manifold of all such ``(p,q)`` pairs with exact gcd of +degree ``k``. A given pair ``(p,q)`` with exact gcd of degree ``j`` +will have some distance ``Θᵏ`` from ``Pᵏ``. For a given threshold +``ϵ>0`` a numerical GCD of ``(p,q)`` within ``ϵ`` is an exact GCD of a +pair ``(p̂,q̂)`` in ``Ρᵏ`` with + +``‖ (p,q) - (p̂,q̂) ‖ ≤ Θᵏ``, where ``k`` is the largest value for +which ``Θᵏ < ϵ``. (In the ``ϵ → 0`` limit, this would be the exact GCD.) Suppose ``(p,q)`` is an ``ϵ`` pertubation from ``(p̂,q̂)`` where ``(p̂,q̂)`` has an exact gcd of degree ``k``, then ``degree(gcdₑ(p,q)) = k``; as ``ϵ → 0``, ``gcdₑ(p,q) → gcd(p̂, q̂)``, and -``limsup_{(p,q)→(p̂,q̂)} inf{ ‖ (u,v,w) - (û,v̂,ŵ) ‖} / ‖ (p,q) - (p̂,q̂) ‖ < κₑ(p,q)``. +``\\limsup_{(p,q)→(p̂,q̂)} \\inf{ ‖ (u,v,w) - (û,v̂,ŵ) ‖} / ‖ (p,q) - (p̂,q̂) ‖ < κₑ(p,q)``. ``κ`` is called the numerical GCD condition number. -The Zeng algorithm proposes a degree for ``u`` and *if* a triple ``(u,v,w)`` with ``u`` of degree ``k`` and ``(u⋅v, u⋅w)`` in ``Ρᵏmn`` can be found satisfying ``‖ (u⋅v, u⋅w) - (p,q) ‖ < ϵ`` then ``(u,v,w)`` is returned; otherwise the proposed degree is reduced and the process repeats. If not terminated, at degree ``0`` a constant gcd is returned. +The Zeng algorithm proposes a degree for ``u`` and then *if* a triple +``(u,v,w)`` with ``u`` of degree ``k`` and ``(u⋅v, u⋅w)`` in ``Ρᵏₘₙ`` +can be found satisfying ``‖ (u⋅v, u⋅w) - (p,q) ‖ < ϵ`` then +``(u,v,w)`` is returned; otherwise the proposed degree is reduced and +the process repeats. If not terminated, at degree ``0`` a constant gcd +is returned. -The initial proposed degree is the first ``j``, ``j=min(m,n):-1:1``, where ``Sⱼ`` is believed to have a singular value of ``0`` (``Sⱼ`` being related to the Sylvester matrix of `ps` and `qs`). The verification of the proposed degree is done using a Gauss-Newton iteration scheme holding the degree of ``u`` constant. +The initial proposed degree is the first ``j``, `j=min(m,n):-1:1`, +where ``Sⱼ`` is believed to have a singular value of ``0`` (``Sⱼ`` +being related to the Sylvester matrix of `p` and `q`). The +verification of the proposed degree is done using a Gauss-Newton +iteration scheme holding the degree of ``u`` constant. ## Scaling: -If `scale=true` the gcd of ``p/‖p‖`` and ``q/‖q‖`` is identified. Scaling can reduce the condition numbers significantly. +If `scale=true` the gcd of ``p/‖p‖₂`` and ``q/‖q₂`` is identified. When the +polynomials have large norms, scaling -- or using a relative tolerance +-- can be necessary to numerically identify the degree of the gcd. ## Tolerances: There are two places where tolerances are utilized: -* in the identification of the rank of ``Sⱼ`` a value ``σ₁ = ‖Rx‖`` is identified. To test if this is zero a tolerance of `max(satol, ‖R‖ₒₒ ⋅ srtol)` is used. +* For a given `k`, the algorithm refines values `u,v,w`. The value `Θᵏ` is estimated by the difference between ``(u ⋅ v, u ⋅ w)`` and ``(p,q)``. A tolerance of `ρ` is used to test if this is smaller than specified. The arguments `atol` and `rtol` are used to compute `ϵ=max(atol, (‖(p,q)‖₂)*κ*rtol)` -* to test if ``(u ⋅ v, u ⋅ w) ≈ (p,q)`` a tolerance of `max(atol, λ⋅rtol)` is used with `λ` chosen to be ``(‖(p,q)‖⋅n⋅m)⋅κ′⋅‖A‖ₒₒ`` to reflect the scale of ``p`` and ``q`` and an estimate for the condition number of ``A`` (a Jacobian involved in the solution). +* The value `ϵ` is also used to determine if the Sylvester matrix for a given `j`, `Sⱼ`, is singular. The theory has ``ϵ`` the same a above, but we this implementation uses `ρ = max(satol, ‖(p,q)‖₂*srtol)`, which seems to track the scaling that is needed due to floating point approximations. The theory states that if `Θᵏ < ϵ` then `σ₋₁ < ϵ √(m - j + 1)`. - -This seems to work well for a reasaonable range of polynomials, however there can be issues: when the degree of ``p`` is much larger than the degree of ``q``, these choices can fail; when a higher rank is proposed, then too large a tolerance for `rtol` or `atol` can lead to a false verification; when a tolerance for `atol` or `rtol` is too strict, then a degree may not be verified. - -!!! note: - These tolerances are adjusted from those proposed in [1]. +The default choice for `ϵ` works reasonably well for a range of polynomials, but scaling or some other choice of `ϵ` is needed for some cases. ## Specified degree: -When `k` is specified, a value for ``(u,v,w)`` is identified with ``degree(u)=k``. No tolerances are utilized in computing ``Θᵏ``. - +When `k` is specified, a value for ``(u, v, w)`` is identified with ``degree(u)=k``. No tolerances are utilized in computing ``Θᵏ``. Output: -The function outputs a named tuple with names (`u`, `v`, `w`, `Θ`, `κ`). The components `u`,`v`,`w` estimate the gcd and give the divisors. The value `Θ` estimates ``Θᵏ`` and `κ` estimates the numerical condition number. +The function outputs a named tuple with names (`u`, `v`, `w`, `Θ`, `κ`). The components `u`,`v`,`w` estimate the gcd and give the divisors. The value `Θ` is the residual error and `κ` estimates the numerical condition number. Example: +```jldoctest ngcd +julia> using Polynomials + +julia> x = variable(Polynomial{Float64}) +Polynomial(1.0*x) + +julia> p = (x+10)*(x^9 + x^8/3 + 1); + +julia> q = (x+10)*(x^9 + x^8/7 - 6/7); + +julia> degree(gcd(p, q)) +0 + +julia> degree(gcd(p, q, method=:numerical)) # u a degree 1 polynomial +1 ``` -using Polynomials -x = variable(Polynomial{Float64}) -p = (x+10)*(x^9 + x^8/3 + 1) -q = (x+10)*(x^9 + x^8/7 - 6/7) -gcd(p,q) # u a constant -gcd(p,q, method=:numerical) # u a degree 1 polynomial -Polynomials.NGCD.ngcd(coeffs(p), coeffs(q), verbose=true) # to see some computations + +This example perturbs `q` more than the default tolerance, so `atol` is set. We can see that the residual error found is on the order of `1e-9`: + +```jldoctest ngcd +julia> p = (x-1)*(x-2)*(x-3); + +julia> q = (x-2)*(x-3)*(x-4) + (1e-8)*x^2; + +julia> out = Polynomials.ngcd(p, q, atol=1e-8); + +julia> degree(out.u) +2 + +julia> round(log10(out.Θ)) +-9.0 + +julia> out = Polynomials.ngcd(p, q); + +julia> degree(out.u) +0 ``` + Reference: [1] The Numerical Greatest Common Divisor of Univariate Polynomials @@ -151,338 +225,272 @@ Note: Based on work by Andreas Varga """ function ngcd(p::PnPolynomial{T,X}, q::PnPolynomial{T,X}; - scale::Bool=false, - atol = eps(real(T)), + atol = eps(real(T))^(5/6), # residual over Θᵏ rtol = Base.rtoldefault(real(T)), - satol = eps(real(T))^(5/6), - srtol = eps(real(T)), - verbose=false, - minⱼ = -1, - λ = 1 + satol = atol, # singular tolerance + srtol = eps(real(T)), + scale::Bool=false, + λ::Real = one(real(T)), # not used + minⱼ = -1 ) where {T, X} m,n = length(p)-1, length(q)-1 + (m == 1 || n == 0) && return trivial_gcd(p, q) + @assert m >= n - ## --- begin + # scale + np₂, nq₂ = norm(p,2), norm(q,2) if scale - p ./= norm(p) - q ./= norm(q) + p ./= np₂ + q ./= nq₂ + np₂ = nq₂ = one(T) end - atol *= λ - rtol *= λ + npq₂ = sqrt(np₂^2 + nq₂^2) - # storage - A0 = zeros(T, m+1, 2) - A0[:,1] = coeffs(p) - A0[end-n:end,2] = coeffs(q) - # pre-allocate storage for Sylvester Matrices, S₁, S₂... - Q = zeros(T, m + n, m + n) + # pre-allocate storage + Q = zeros(T, m + n, m + n) # for QR decomposition of Sylvester matrices R = zeros(T, m + n, m + n) - Sₓ = hcat(convmtx(p,1), convmtx(q, m-n+1)) + uv = copy(p) # storage for uᵢ * vᵢ + uw = copy(q) # storage for uᵢ * wᵢ + x = Vector{T}(undef, m + n) # used to find σ₋₁ - uv = copy(p) - uw = copy(q) + # j is degree of proposed gcd j ≤ n ≤ m + j = n # We count down Sn, S_{n-1}, ..., S₂, S₁ + Sₓ = SylvesterMatrix(p, q, j) # initial Sylvester matrix [Cₙ₋ⱼ₊₁(p), Cₘ₋ⱼ₊₁(q)] - local x::Vector{T} + A0 = zeros(T, m+1, 2) # storage for use with extend_QR! + A0[:,1] = coeffs(p) + A0[end-n:end,2] = coeffs(q) - F = qr(Sₓ) nr, nc = size(Sₓ) # m+1, m-n+2 + F = qr(Sₓ) Q[1:nr, 1:nr] .= F.Q R[1:nc, 1:nc] .= F.R - j = n # We count down Sn, S_{n-1}, ..., S₂, S₁ + # tolerances + atol, satol, rtol, srtol = λ*atol, λ*satol, λ*rtol, λ*srtol + ρ = max(satol, npq₂ * srtol) while true - - V = view(R, 1:nc, 1:nc) - flag, σ, x = smallest_singular_value(V, satol * sqrt(1 + m - j), srtol) - verbose && println("------ degree $j ----- σ₁: $σ --- $flag") - - if (flag == :iszero || flag == :ispossible) - u, v, w = initial_uvw(Val(flag), j, p, q, x) - flag, ρ₁, σ₂, ρ = refine_uvw!(u,v,w, p, q, uv, uw, atol, rtol) - - verbose && println(" --- Θᵏ: $ρ₁ --- $flag (ρ=$(ρ))") - - if flag == :convergence - return (u=u, v=v, w=w, Θ=ρ₁, κ=σ₂) # (u,v,w) verified + V = UpperTriangular(view(R, 1:nc, 1:nc)) + xx = view(x, 1:nc) + σ₋₁ = smallest_singular_value!(xx, V, ρ * sqrt(m - j + 1)) + #@show j, σ₋₁, ρ * sqrt(m - j + 1), npq₂ + + # Lemma 7.1: If (p,q) is w/in ϵ of P^k_{mn} then σ₋₁ < ϵ√(m-j+1) + if σ₋₁ ≤ ρ * sqrt(m - j + 1) + # candidate for degree; refine u₀, vₒ, w₀ to see if ρ < ϵ + if iszero(σ₋₁) + # determinant is 0 + u, v, w = initial_uvw(Val(:iszero), j, p, q, xx) + else + u, v, w = initial_uvw(Val(:ispossible), j, p, q, xx) + end + ϵₖ, κ = refine_uvw!(u, v, w, p, q, uv, uw) + # we have limsup Θᵏ / ‖(p,q) - (p̃,q̃)‖ = κ, so + # ‖Θᵏ‖ ≤ κ ⋅ ‖(p,q)‖ ⋅ ϵ seems a reasonable heuristic. + # Too tight a tolerance and the right degree will be missed; too + # lax, and larger degrees will be accepted. We are using + # `√eps()` for `rtol`, but that may be too lax and is subject to + # change. + ϵ = max(atol, npq₂ * κ * rtol) + #@show ϵₖ, ϵ, κ + if ϵₖ ≤ ϵ + #@show :success, σ₋₁, ϵₖ + return (u=u, v=v, w=w, Θ=ϵₖ, κ=κ) end + #@show :failure, j end # reduce possible degree of u and try again with Sⱼ₋₁ # unless we hit specified minimum, in which case return it + # minⱼ = -1 if j == minⱼ - u, v, w = initial_uvw(Val(:ispossible), j, p, q, x) - flag, ρ₁, σ₂, ρ = refine_uvw!(u,v,w, p, q, uv, uw, atol, rtol) - return (u=u, v=v, w=w, Θ=ρ₁, κ=σ₂) + u, v, w = initial_uvw(Val(:ispossible), j, p, q, xx) + ϵₖ, κ = refine_uvw!(u, v ,w, p, q, uv, uw) + return (u=u, v=v, w=w, Θ=ϵₖ, κ=κ) end + # Try again with a smaller j j -= 1 nr += 1 nc += 2 nc > nr && break - extend_QR!(Q,R, nr, nc, A0) # before Q⋅R = Sⱼ, now Q⋅R = Sⱼ₋₁ + extend_QR!(Q, R, nr, nc, A0) # before Q⋅R = Sⱼ, now Q⋅R = Sⱼ₋₁ + end + return trivial_gcd(p, q) - end +end - # u is a constant - verbose && println("------ GCD is constant ------") +# fixed `k` +function ngcd(p::P, q::P, k::Int) where {T, X, P <: PnPolynomial{T,X}} - u, v, w = initial_uvw(Val(:constant), j, p, q, x) - flag, ρ₁, κ, ρ = refine_uvw!(u,v,w, p, q, uv, uw, atol, rtol) - return (u=u, v=v, w=w, Θ=ρ₁, κ=κ) + m, n = length(p)-1, length(q)-1 + Sₓ = SylvesterMatrix(p,q,k) + F = qr(Sₓ) + R = UpperTriangular(F.R) + x = zeros(T, size(Sₓ, 2)) + np = norm(p) + σ₋₁ = smallest_singular_value!(x, R) + w, v = P(x[1:(n-k+1)]), P(-x[(n-k+2):end]) + u = solve_u(v, w, p, q, k) + ρₖ, κ = refine_uvw!(u, v, w, p, q, u*v, u*w) + return (u=u, v=v, w=w, Θ=ρₖ, κ=κ) +end +function trivial_gcd(p::P, q) where {T, X, P <: PnPolynomial{T, X}} + u, v, w = one(P), p, q + return (u=u, v=v, w=w, Θ=zero(T), κ=NaN) end -# fix the degree, k -function ngcd(p::P, - q::P, - k::Int; - kwargs... - ) where {T <: AbstractFloat,X, P <: PnPolynomial{T,X}} - m::Int, n::Int = length(p)-1, length(q)-1 +# A = QR solve by iteration for largest eigenvalue of A^-1 +# modifies w in place - # u,v,w = initial_uvw(Val(:iszero), k, ps, qs, nothing) - Sⱼ = [convmtx(p, n-k+1) convmtx(q, m-k+1)] - F = qr(Sⱼ) - flag, σ, x = smallest_singular_value(F.R, eps(T) * sqrt(1 + m - k), eps(T)) - u,v,w = initial_uvw(Val(:k), flag, k, p, q, x) - flag, ρ₁, κ, ρ = refine_uvw!(u,v,w, copy(p), copy(q), copy(p), copy(q), - T(Inf), T(Inf)) # no tolerances - return (u=u, v=v, w=w, Θ=ρ₁, κ=κ) +# https://arxiv.org/abs/2103.04196 +# Find smallest singular value +# stop when value is less then θ or Δ=sⱼ - s₋₁ is small +function smallest_singular_value!(w, R::UpperTriangular{T}, + θ, + ϵₘ = eps(real(T)) + ) where {T} -end + # Cant' handle singular matrices + iszero(det(R)) && return zero(T) + nRₒₒ = norm(R, Inf) + MAXSTEPS = 50 + sⱼ = sⱼ₋₁ = typemax(real(T)) -## ----- + w .= one(T) + w ./= norm(w) -# return guess at smallest singular value and right sinuglar value, x -# for an upper triangular matrix, V -function smallest_singular_value(V::AbstractArray{T,2}, - atol=eps(real(T)), - rtol=zero(real(T))) where {T} + j = 1 + while true + sⱼ₋₁ = sⱼ + sⱼ = smallest_singular_value_one_step!(w, R) - R = UpperTriangular(V) - k = size(R)[1]/2 - if iszero(det(R)) - return (:iszero, zero(T), T[]) + if sⱼ ≤ θ || abs(sⱼ - sⱼ₋₁) ≤ max(1, sⱼ) * ϵₘ * nRₒₒ + break + end + + j += 1 + j >= MAXSTEPS && break end - m,n = size(R) + return sⱼ + +end + +# no tolerance; stop when improvment stops +function smallest_singular_value(A) + R = UpperTriangular(qr(A).R) + w = Vector{eltype(R)}(undef, size(R, 2)) + smallest_singular_value!(w, R) +end + +function smallest_singular_value!(w, R::UpperTriangular{T}) where {T} + iszero(det(R)) && return zero(T) + + MAXSTEPS = 50 - # we are testing if ‖Ax‖ ≈ 0 - # If x is a perfect 0, but x is approximated by x' due to round off - # then ‖A(x-x')‖ <= ‖A‖⋅‖x - x'‖ so we use ‖A‖ as scale factor - δ = max(atol, norm(R,Inf) * rtol) + sⱼ = typemax(real(T)) - x = ones(T, n) - y = zeros(T, m) - σ₀ = σ₁ = Inf*one(real(T)) - steps, minᵢ = 1, 5 + w .= one(T) + w ./= norm(w) + j = 1 + + w̃ = copy(w) while true - y .= R' \ x # use iteration to converge on singular value - x .= R \ y - x ./= norm(x,2) - σ₁ = norm(R * x, 2) + s′ = smallest_singular_value_one_step!(w̃, R) - if (steps <= 50) && (steps <= minᵢ || σ₁ < 0.05*σ₀) # decreasing, keep going - σ₀ = σ₁ - else - break - end - steps += 1 - end + s′ > sⱼ && break + + copy!(w, w̃) + sⱼ = s′ + j += 1 + j > MAXSTEPS && break - if σ₁ < δ - return (:ispossible, σ₁, x) - else - return (:constant, σ₁, x) end + return sⱼ + +end + + +# modify w, return sⱼ after one step +# uses R from QR factorization +function smallest_singular_value_one_step!(w, R) + ldiv!(R', w) + w ./= norm(w,2) + ldiv!(R, w) + sⱼ = 1/norm(w, 2) + w .*= sⱼ + return sⱼ end -## -------------------------------------------------- -## Refine u,v,w +# solve for u from [v,w] \ [p,q] +function solve_u(v::P, w, p, q, j) where {T, X, P<:PnPolynomial{T,X}} + A = [convmtx(v, j+1); convmtx(w, j+1)] + b = vcat(coeffs(p), coeffs(q)) + u = A \ b + return P(u) +end ## Find u₀,v₀,w₀ from right singular vector function initial_uvw(::Val{:ispossible}, j, p::P, q::Q, x) where {T,X, P<:PnPolynomial{T,X}, Q<:PnPolynomial{T,X}} # Sk*[w;-v] = 0, so pick out v,w after applying permutation - m,n = length(p)-1, length(q)-1 + m, n = length(p)-1, length(q)-1 vᵢ = vcat(2:m-n+2, m-n+4:2:length(x)) wᵢ = m-n+3 > length(x) ? [1] : vcat(1, (m-n+3):2:length(x)) - # v = 𝑷{m-j}(-x[vᵢ]) + v = P(-x[vᵢ]) w = P(x[wᵢ]) + # p194 3.9 C_k(v) u = p or Ck(w) u = q; this uses 10.2 - u = solve_u(v,w,p,q,j) - return u,v,w + u = solve_u(v, w, p, q, j) + return u, v, w end +# find u₀, v₀. w₀ when R is singular. function initial_uvw(::Val{:iszero}, j, p::P, q::Q, x) where {T,X, P<:PnPolynomial{T,X}, Q<:PnPolynomial{T,X}} m,n = length(p)-1, length(q)-1 - S = [convmtx(p, n-j+1) convmtx(q, m-j+1)] - + S = SylvesterMatrix(p,q,j) F = qr(S) R = UpperTriangular(F.R) if iszero(det(R)) - x = eigvals(R)[:,1] + x .= eigvals(R)[:,1] else - x = ones(T, size(R,2)) - x .= R \ (R' \ (x/norm(x))) - x ./= norm(x) + x .= one(T) + smallest_singular_value!(x, R) end - w = P(x[1:n-j+1]) + w = P(x[1:n-j+1]) # ordering of S is not interlaced v = P(-x[(n-j+2):end]) u = solve_u(v,w,p,q,j) return u,v,w end -function initial_uvw(::Val{:constant}, j, p::P, q, x) where {T,X,P<:PnPolynomial{T,X}} - u = one(P) - w = q - v = p - u,v,w -end - -function initial_uvw(::Val{:k}, flag, k, p::P, q, x) where {T,X,P<:PnPolynomial{T,X}} - flag == :iszero && return initial_uvw(Val(flag), k, p, q, nothing) - n = length(q)-1 - w, v = P(x[1:n-k+1]), P(-x[n-k+2:end]) - u = solve_u(v,w,p,q,k) - return (u,v,w) -end - - - -# find estimate for σ₂, used in a condition number (κ = 1/σ) -function σ₂(J) - F = qr(J) - flag, σ, x = smallest_singular_value(F.R) - σ -end - -## attempt to refine u,v,w -## check that [u * v; u * w] ≈ [p; q] -function refine_uvw!(u::U, v::V, w::W, p, q, uv, uw, atol, rtol) where {T,X, - U<:PnPolynomial{T,X}, - V<:PnPolynomial{T,X}, - W<:PnPolynomial{T,X}} - - m,n,l = length(u)-1, length(v)-1, length(w)-1 - - mul!(uv, u, v) - mul!(uw, u, w) - - ρ₀, ρ₁ = one(T), residual_error(p,q,uv,uw) - - # storage - b = zeros(T, (m+n) + (m+l) + 3) # degree(p) + degree(q) + 3 = 1 + length(p) + length(q)) - Δf = zeros(T, m + n + l + 3) - steps = 0 - - h, β = u, norm(u)^2 - minᵢ, Maxᵢ = 5, 20 - κ = NaN - A=zeros(T, JF_size(u, v, w)...) - JF!(A, h, u, v, w) - Fmp!(b, dot(h,u) - β, p, q, uv, uw) - - Δvᵢ = 1:(n+1) - Δwᵢ = (n+1+1):(n+1+l+1) - Δuᵢ = (n+1+l+1+1):length(Δf) - - while ρ₁ > 0.0 - - # update A, b, then solve A\b - qrsolve!(Δf, A, b) - - # m + n = degree(p) - # m + l = degree(q) - # b has length degree(p)+degree(q) + 3 - Δv = view(Δf, Δvᵢ) - Δw = view(Δf, Δwᵢ) - Δu = view(Δf, Δuᵢ) - - u .-= Δu - v .-= Δv - w .-= Δw - - mul!(uv, u, v) - mul!(uv, u, w) - - ρ₀, ρ′ = ρ₁, residual_error(p, q, uv, uw) - - # don't worry about first few, but aftewards each step must be productive - # though we can have really bad first steps, which we cap - if (steps <= Maxᵢ) && (steps <= minᵢ || ρ′ < 0.95 * ρ₀) && ( ρ′ < 100*ρ₁ ) - ρ₁ = ρ′ - steps += 1 - else - break - end - - # update A,b for next iteration - JF!(A, h, u, v, w) - Fmp!(b, dot(h,u) - β, p, q, uv, uw) - - end - - - # this is a heuristic - # sensitivity is Δu / Δp <= ‖ A+ ‖ = κ - # we use an estimate for ‖(p,q)‖ error times ‖A⁺‖⋅‖A‖ₒₒ - κ = 1/σ₂(A) # ≈ ‖A⁺‖ - λ = norm((norm(p), norm(q))) * (m * n) * min(1, κ) * norm(A, Inf) - ρ = max(atol, rtol * λ) - - if ρ₁ <= ρ - return :convergence, ρ₁, κ, ρ - else - return :no_convergence, ρ₁, κ, ρ - end - -end - -## ---- QR factorization - -function qrsolve!(y::Vector{T}, A, b) where {T} - y .= A \ b -end - -# # Fast least-squares solver for full column rank Hessenberg-like matrices -# # By Andreas Varga -function qrsolve!(y::Vector{T}, A, b) where {T <: Float64} - Base.require_one_based_indexing(A) - m, n = size(A) - m < n && error("Column dimension exceeds row dimension") - _, τ = LinearAlgebra.LAPACK.geqrf!(A) - T <: Complex ? tran = 'C' : tran = 'T' - LinearAlgebra.LAPACK.ormqr!('L', tran, A, τ, view(b,:,1:1)) - y .= UpperTriangular(triu(A[1:n,:]))\b[1:n] -end # extend QR to next size # Q gets a 1 in nc,nc, 0s should be elswhere -function extend_QR!(Q,R, nr, nc, A0) - +function extend_QR!(Q, R, nr, nc, A0) #old Q is m x m, old R is n x n; we add to these - n = nc-2 + n = nc - 2 m = nr - 1 k,l = size(A0) @@ -511,10 +519,90 @@ function extend_QR!(Q,R, nr, nc, A0) end +## attempt to refine u,v,w +## return residual error, ρ, estimate for 1/σ_2, κ +function refine_uvw!(u::P, v::P, w::P, + p, q, uv, uw) where {T,X, + P<:PnPolynomial{T,X}} + + mul!(uv, u, v) + mul!(uw, u, w) + ρ₁ = residual_error(p, q, uv, uw) + iszero(ρ₁) && return (ρ₁, NaN) + # storage + h, β = u, dot(u,u) # h = constant * u₀ is used + A = JF(h, u, v, w) + Δfβ = Fmp(dot(h, u) - β, p, q, uv, uw) + Δz = ones(T, length(u) + length(v) + length(w)) + n = size(A, 2) + R = UpperTriangular(Matrix{T}(undef, n, n)) + R′ = copy(R) + ũ, ṽ, w̃ = copy(u), copy(v), copy(w) + + steps = 0 + #@show steps, ρ₁ + minᵢ, Maxᵢ = 3, 20 -## Jacobian F(u,v,w) = [p,p'] is J(u,v,w) -function JF_size(u, v, w) + while ρ₁ > 0.0 + steps += 1 + refine_uvw_step!(ũ, ṽ, w̃, + Δz, A, Δfβ, R) + mul!(uv, ũ, ṽ) + mul!(uw, ũ, w̃) + ρ′ = residual_error(p, q, uv, uw) + #@show steps, ρ′ + # don't worry as much about first few, + # but afterwards each step must be productive + # terminate when no longer decreasing + #if steps < minᵢ || (steps ≤ Maxᵢ && ρ′ < 0.95*ρ₁) + if ρ′ < ρ₁ || (steps ≤ minᵢ && ρ′ ≤ 1.1*ρ₁) + ρ₁ = ρ′ + copy!(R′, R) + copy!(u.coeffs, ũ.coeffs) + copy!(v.coeffs, ṽ.coeffs) + copy!(w.coeffs, w̃.coeffs) + steps ≥ Maxᵢ && break + # update A,b for next iteration + JF!(A, h, u, v, w) + Fmp!(Δfβ, dot(h, u) - β, p, q, uv, uw) + else + steps == 1 && copy!(R′, R) + break + end + end + smallest_singular_value_one_step!(Δz, R′) # two steps, not one + σ₂ = smallest_singular_value_one_step!(Δz, R′) + κ = 1/σ₂ + return ρ₁, κ + +end + +# update u,v,w, uv, uw +# computes update step of zₖ₊₁ = zₖ - Δz; Δz = J(zₖ)⁺(fₕ(u,v,w) - [β;p;q]) +function refine_uvw_step!(u, v, w, + Δz, J⁺, Δfβ, R) + + qrsolve!(Δz, J⁺, Δfβ, R) # Δz .= J⁺ \ Δfβ + + m,n,l = length(u)-1, length(v)-1, length(w)-1 + Δuᵢ = 1:(m+1) + Δvᵢ = (m+1+1):(m+1 + n+1) + Δwᵢ = (m + 1 + n + 1 + 1):(m + n + l + 3) + + u .-= view(Δz, Δuᵢ) + v .-= view(Δz, Δvᵢ) + w .-= view(Δz, Δwᵢ) + +end + +## Jacobian of F(u,v,w) = [p, p'] is J(u,v,w) +## [h 0 0; +## Cₖ(v) Cₘ₋ₖ(u) 0; +## Cₖ(w) 0 Cₙ₋ₖ(u)] +function JF(h, u::P, v, w) where {T,X,P<:AbstractPolynomial{T,X}} + + # compute size needed to store m, k, j = length(u)-1, length(v)-1, length(w)-1 n, l = m + k, m + j @@ -525,70 +613,105 @@ function JF_size(u, v, w) ci, cj = ai, fj ei, ej = di, bj - (1 + ai + di, aj + bj + cj) -end + m, n = 1 + ai + di, aj + bj + cj -# Jacobian matrix -function JF(u::Vector{U}, v::Vector{V}, w::Vector{W}) where {U,V,W} - R = promote_type(U,V, W) - M = zeros(R, JF_size(u, v, w)...) - JF!(M, u, v, w) - M + A = zeros(T, m, n) + JF!(A, h, u, v, w) + A end +# compute Jacobian of fₕ function JF!(M, h, u::P, v, w) where {T,X,P<:AbstractPolynomial{T,X}} - du, dv, dw = length(u)-1, length(v)-1, length(w)-1 - m, n = du + dv, du + dw - - # JF_size should return these - r11,c11 = convmtx_size(u, dv+1) - r13,c13 = convmtx_size(v, du+1) - r22,c22 = convmtx_size(u, dw+1) - r23,c23 = convmtx_size(w, du+1) - - J11 = view(M, 1:r11, 1:c11) - J13 = view(M, 1:r13, c11 + c22 .+ (1:c23)) - J22 = view(M, r11 .+ (1:r22), c11 .+ (1:c22)) - J23 = view(M, r13 .+ (1:r23), (c11 + c22) .+ (1:c23)) - convmtx!(J11, u, dv+1) - convmtx!(J13, v, du+1) - convmtx!(J22, u, dw+1) - convmtx!(J23, w, du+1) - M[end, end-du:end] = coeffs(h)' + k = length(u) - 1 + dᵥ, dᵥᵥ = length(v) - 1, length(w) - 1 + m = k + dᵥ # degree(p) + n = k + dᵥᵥ # degree(q) + + M[1, 1:length(h)] .= h' + + r11, c11 = convmtx_size(v, k+1) # Ck(v) size + J11 = view(M, 2:(1+r11), 1:c11) + convmtx!(J11, v, k+1) + + r12, c12 = convmtx_size(u, m-k+1) + J12 = view(M, 2:(1+r11), (c11+1):(c11+1+c12)) + convmtx!(J12, u, m-k+1) + + r21, c21 = convmtx_size(w, k+1) + J21 = view(M, (1 + r11 + 1):(1 + r11 + r21), 1:c21) + convmtx!(J21, w, k+1) + + r23, c23 = convmtx_size(u, n - k+1) + J23 = view(M, (1 + r11 + 1):(1 + r11 + r23), + (c21 + c12 + 1):(c21 + c12 + c23)) + convmtx!(J23, u, n-k+1) return nothing end -## compute F(u,v,w) - [p, p'] = [u*v, u*w] - [p, p'] -function Fmp!(b, γ, p, q, uv, uw) - b[end] = γ - for i in 1:1+length(p)-1 - j = i - b[i] = uv[j] - p[j] +# create storage for fₕ(u,v,w) - [β; p; q]; fill in +function Fmp(Δ, p::PnPolynomial{T,X}, q, p̃, q̃) where {T,X} + b = zeros(T, 1 + length(p) + length(q)) + Fmp!(b, Δ, p, q, p̃, q̃) + b +end + +## fₕ - [β; p; q] +## Δ = h⋅uᵢ - β +function Fmp!(b, Δ, p, q, uv, uw) + b[1] = Δ + for (i, pᵢ) ∈ pairs(p) + b[2 + i] = uv[i] - pᵢ end - for i in 1+length(p):length(b)-1 - j = i - length(p) - b[i] = uw[j] - q[j] + for (i, qᵢ) ∈ pairs(q) + b[2 + length(p) + i] = uw[i] - qᵢ end return nothing end - -function residual_error(p::P,q,uv,uw) where {T,X,P<:AbstractPolynomial{T,X}} +# find ||(p,q) - (p̃, q̃)|| treating (,) as vector concatenation of coefficients +function residual_error(p::P, q, p̃, q̃) where {T,X,P<:AbstractPolynomial{T,X}} tot = zero(real(T)) - for (pᵢ, uvᵢ) in zip(p,uv) - tot += norm(pᵢ-uvᵢ)^2 + for (pᵢ, p̃ᵢ) in zip(p, p̃) + tot += abs2(pᵢ - p̃ᵢ) end - for (qᵢ, uwᵢ) in zip(q, uw) - tot += norm(qᵢ-uwᵢ)^2 + for (qᵢ, q̃ᵢ) in zip(q, q̃) + tot += abs2(qᵢ - q̃ᵢ) end sqrt(tot) end +## ---- utils +## ---- QR factorization +function qrsolve!(y::Vector{T}, A, b) where {T <: Float64} + n = size(A, 2) + R = UpperTriangular(Matrix{Float64}(undef, n, n)) + qrsolve!(y, A, b, R) +end +function qrsolve!(y::Vector{T}, A, b, R) where {T} + q = qr(A) + R .= UpperTriangular(q.R) + y .= q \ b +end + +# Fast least-squares solver for full column rank Hessenberg-like matrices +# By Andreas Varga +function qrsolve!(y::Vector{T}, A, b, R) where {T <: Float64} + # half the time of ldiv!(y, qr(A), b) + Base.require_one_based_indexing(A) + m, n = size(A) + m < n && error("Column dimension exceeds row dimension") + _, τ = LinearAlgebra.LAPACK.geqrf!(A) + R .= UpperTriangular(triu(A[1:n,1:n])) + + tran = T <: Complex ? 'C' : 'T' + LinearAlgebra.LAPACK.ormqr!('L', tran, A, τ, view(b,:,1:1)) + R = UpperTriangular(triu(A[1:n,:])) + ldiv!(y, R, view(b,1:n)) +end -## ---- utils """ convmtx(v, n::Int) @@ -631,13 +754,81 @@ function convmtx(v::AbstractPolynomial{T}, n::Int) where {T} C end +function SylvesterMatrix(p, q, j) + m, n = length(p)-1, length(q) - 1 + Sₓ = hcat(convmtx(p, n - j + 1 ), convmtx(q, m - j + 1)) +end -# solve for u from [v,w] \ [p,q] -function solve_u(v::P,w,p,q, k) where {T,X,P<:PnPolynomial{T,X}} - A = [convmtx(v,k+1); convmtx(w, k+1)] - b = vcat(coeffs(p), coeffs(q)) - u = A \ b - return P(u) +## ---- + +# non allocating numeric-rank reveal +function rank_reveal!(R::LinearAlgebra.UpperTriangular{T}, w, + nAₒₒ, nRₒₒ; + atol = Base.rtoldefault(real(T)), + rtol = eps(real(T))) where {T <: AbstractFloat} + + n = size(R, 2) + + d = prod(R[i,i] for i ∈ 1:n) + iszero(d) && return(0, zero(T)) # Cant' handle singular matrices + + + MAXSTEPS = 50 + + θ = max(atol, nAₒₒ * rtol) + + ϵₘ = nAₒₒ * eps(real(T)) + r = n + sⱼ = sⱼ₋₁ = typemax(real(T)) + + for k ∈ 1:n + sⱼ = smallest_singular_value!(w, R, θ, ϵₘ) + + if sⱼ > θ + break + end + + r -= 1 + + # W = [w W] + # can eliminate allocations here! + # RR = [τ*w'; R] + # for i ∈ 1:n + # G,_ = givens(RR, i, i+1, i) + # lmul!(G, RR) + # end + # R = UpperTriangular(RR[1:end-1, :]) + + τ = nRₒₒ + a,b = τ*w[1], R[1,1] + g,d = givens(a, b, 1, 2) + R[1,1] = d + for j in 2:n + a, b = τ*w[j], R[1, j] + R[1,j] = conj(g.c)*a + conj(g.s)*b + w[j] = R[2, j] + R[2,j] = -conj(g.s)*a + conj(g.c)*b + end + + for i ∈ 2:(n-1) + a, b = R[i, i], w[i] + g, d = givens(a, b, 1, 2) + R[i,i] = d + for j ∈ (i+1):n + a, b = R[i,j], w[j] + R[i,j] = conj(g.c)*a + conj(g.s)*b + w[j] = R[i+1, j] + R[i+1, j] = -conj(g.s)*a + conj(g.c)*b + end + end + + # n + g, d = givens(R[n,n], w[n], n-1,n) + R[n,n] = d + + end + r, sⱼ end + end diff --git a/src/polynomials/pi_n_polynomial.jl b/src/polynomials/pi_n_polynomial.jl index cc98827d..370045d7 100644 --- a/src/polynomials/pi_n_polynomial.jl +++ b/src/polynomials/pi_n_polynomial.jl @@ -21,7 +21,8 @@ struct PnPolynomial{T,X} <: StandardBasisPolynomial{T, X} end end - +PnPolynomial{T, X}(coeffs::Tuple) where {T, X} = + PnPolynomial{T,X}(T[pᵢ for pᵢ ∈ coeffs]) Polynomials.@register PnPolynomial @@ -46,7 +47,7 @@ function LinearAlgebra.mul!(pq, p::PnPolynomial{T,X}, q) where {T,X} for i ∈ 0:m for j ∈ 0:n k = i + j - @inbounds pq.coeffs[1+k] += p.coeffs[1+i] * q.coeffs[1+j] + @inbounds pq.coeffs[1+k] = muladd(p.coeffs[1+i], q.coeffs[1+j], pq.coeffs[1+k]) end end nothing diff --git a/src/polynomials/standard-basis.jl b/src/polynomials/standard-basis.jl index 1a66b6a6..8bcd2e58 100644 --- a/src/polynomials/standard-basis.jl +++ b/src/polynomials/standard-basis.jl @@ -90,13 +90,13 @@ end ## put here, not with type defintion, in case reuse is possible ## `conv` can be used with matrix entries, unlike `fastconv` function conv(p::Vector{T}, q::Vector{S}) where {T,S} - as = [p[1]*q[1]] + as = [p[1]*q[1]] z = 0 * as[1] n,m = length(p)-1, length(q)-1 for i ∈ 1:n+m Σ = z for j ∈ max(0, i-m):min(i,n) - Σ += p[1+j]*q[1 + i-j] + Σ = muladd(p[1+j], q[1 + i-j], Σ) end push!(as, Σ) end diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl index 86786af9..decbbc2c 100644 --- a/test/StandardBasis.jl +++ b/test/StandardBasis.jl @@ -862,7 +862,7 @@ end n = 4 q = p^n out = Polynomials.Multroot.multroot(q) - @test out.κ * out.ϵ > sqrt(eps()) # large forward error, l misidentified + @test (out.multiplicities == n*ls) || (out.κ * out.ϵ > sqrt(eps())) # large forward error, l misidentified # with right manifold it does yield a small forward error zs′ = Polynomials.Multroot.pejorative_root(q, rts .+ 1e-4*rand(3), n*ls) @test prod(Polynomials.Multroot.stats(q, zs′, n*ls)) < sqrt(eps()) @@ -1013,7 +1013,8 @@ end p = P([1,2,3], :x) A = [1 p; p^2 p^3] @test !issymmetric(A) - @test issymmetric(A*transpose(A)) + U = A * A' + @test U[1,2] ≈ U[2,1] # issymmetric with some allowed error for FactoredPolynomial diagm(0 => [1, p^3], 1=>[p^2], -1=>[p]) end @@ -1289,9 +1290,9 @@ end d = P([0.5490673726445683, 0.15991109487875477]); @test degree(gcd(a*d,b*d)) == 0 @test degree(gcd(a*d, b*d, atol=sqrt(eps()))) > 0 - @test degree(gcd(a*d,b*d, method=:noda_sasaki)) == degree(d) - @test_skip degree(gcd(a*d,b*d, method=:numerical)) == degree(d) # issues on some architectures - l,m,n = (5,5,5) # realiable, though for larger l,m,n only **usually** correct + @test degree(gcd(a*d,b*d, method=:noda_sasaki)) == degree(d) + @test_skip degree(gcd(a*d,b*d, method=:numerical)) == degree(d) # issues on some architectures (had test_skip) + l,m,n = (5,5,5) # sensitive to choice of `rtol` in ngcd u,v,w = fromroots.(rand.((l,m,n))) @test degree(gcd(u*v, u*w, method=:numerical)) == degree(u) @@ -1320,7 +1321,7 @@ end W(n) = prod( (x-r1*alpha(j,n))^2 + r1^2*beta(j,n)^2 for j in (n+1):2n) @testset for n in 2:2:20 p = U(n) * V(n); q = U(n) * W(n) - @test degree(gcd(p,q, method=:numerical)) == degree(U(n)) + @test degree(gcd(p,q; method=:numerical)) == degree(U(n)) end # Test 5 of Zeng @@ -1331,15 +1332,16 @@ end p = prod((x-i)^j for (i,j) in enumerate(ms)) dp = derivative(p) - @test degree(gcd(p,dp, method=:numerical)) == sum(max.(ms .- 1, 0)) + @test degree(gcd(p,dp; method=:numerical)) == sum(max.(ms .- 1, 0)) end # fussy pair x = variable(P{Float64}) - @testset for n in (2,5,10,20,50, 100) + @testset for n in (2,5,10,20,25,50, 100) p = (x-1)^n * (x-2)^n * (x-3) q = (x-1) * (x-2) * (x-4) - @test degree(gcd(p,q, method=:numerical)) == 2 + a = Polynomials.ngcd(p, q) + a.κ < 100 && @test degree(a.u) == 2 end # check for fixed k