|
| 1 | +# Computing roots of Chebyshev real-valued polynomials in 1d via colleague matrices on |
| 2 | +# recursively subdivided intervals, as described in Trefethen, "Approximation |
| 3 | +# Theory and Approximation Practice", chapter 18. https://epubs.siam.org/doi/10.1137/1.9781611975949 |
| 4 | +# |
| 5 | +# Some additional implementation tricks were inspired by the chebfun source |
| 6 | +# code: https://github.com/chebfun/chebfun/blob/7574c77680d7e82b79626300bf255498271a72df/%40chebtech/roots.m |
| 7 | + |
| 8 | +using LinearAlgebra |
| 9 | +export colleague_matrix, roots |
| 10 | + |
| 11 | +# colleague matrix for 1d array of Chebyshev coefficients, assuming |
| 12 | +# trailing zero coefficients have already been dropped. |
| 13 | +function _colleague_matrix(coefs::AbstractVector{<:Number}) |
| 14 | + n = length(coefs) |
| 15 | + n <= 1 && return Matrix{float(eltype(coefs))}(undef,0,0) # 0×0 case (no roots) |
| 16 | + iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero")) |
| 17 | + |
| 18 | + if n == 2 # trivial 1×1 (degree 1) case |
| 19 | + C = Matrix{float(eltype(coefs))}(undef,1,1) |
| 20 | + C[1,1] = float(-coefs[1])/coefs[2] |
| 21 | + return C |
| 22 | + else |
| 23 | + # colleague matrix, transposed to make it upper-Hessenberg (which doesn't change eigenvalues) |
| 24 | + halves = fill(one(float(eltype(coefs)))/2, n-2) |
| 25 | + C = diagm(1 => halves, -1 => halves) |
| 26 | + C[2,1] = 1 |
| 27 | + @views C[:,end] .-= coefs[1:end-1] ./ 2coefs[end] |
| 28 | + return C |
| 29 | + end |
| 30 | +end |
| 31 | + |
| 32 | + |
| 33 | +if !isdefined(LinearAlgebra, :diagview) # added in Julia 1.12 |
| 34 | + if !isdefined(LinearAlgebra, :diagind) |
| 35 | + diagind(A::Matrix) = range(1, step=size(A,1)+1, length=min(size(A)...)) |
| 36 | + end |
| 37 | + diagview(A::Matrix) = @view A[diagind(A)] |
| 38 | +end |
| 39 | + |
| 40 | +function _colleague_matrix(coefs::AbstractVector{<:Number}, lb::Real, ub::Real) |
| 41 | + C = _colleague_matrix(coefs) |
| 42 | + # scale and shift from [-1,1] to [lb,ub] via C * (ub-lb)/2 + (ub+lb)/2 * I |
| 43 | + C .*= (ub - lb)/2 |
| 44 | + diagview(C) .+= (ub + lb)/2 |
| 45 | + return C |
| 46 | +end |
| 47 | + |
| 48 | +""" |
| 49 | + colleague_matrix(c::ChebPoly{1,<:Number}; tol=5eps) |
| 50 | +
|
| 51 | +Return the "colleague matrix" whose eigenvalues are the roots of the 1d real-valued |
| 52 | +Chebyshev polynomial `c`, dropping trailing coefficients whose relative contributions |
| 53 | +are `< tol` (defaulting to 5 × floating-point `eps`). |
| 54 | +""" |
| 55 | +function colleague_matrix(c::ChebPoly{1,<:Number}; tol::Real=5*epsvals(c.coefs)) |
| 56 | + abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol |
| 57 | + n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1) |
| 58 | + return UpperHessenberg(_colleague_matrix(@view(c.coefs[1:n]), c.lb[1], c.ub[1])) |
| 59 | +end |
| 60 | + |
| 61 | +function filter_roots(c::ChebPoly{1,<:Real}, roots::AbstractVector{<:Number}) |
| 62 | + @inbounds lb, ub = c.lb[1], c.ub[1] |
| 63 | + htol = eps(float(ub - lb)) * 100 # similar to chebfun |
| 64 | + return [ clamp(real(r), lb, ub) for r in roots if abs(imag(r)) < htol && lb - htol <= real(r) <= ub + htol ] |
| 65 | +end |
| 66 | + |
| 67 | +if isdefined(LinearAlgebra.LAPACK, :hseqr!) # added in Julia 1.10 |
| 68 | + # see also LinearAlgebra.jl#1557 - optimized in-place eigenvalues for upper-Hessenberg colleague matrix |
| 69 | + function hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T<:Union{LinearAlgebra.BlasReal, LinearAlgebra.BlasComplex}} |
| 70 | + ilo, ihi, _ = LinearAlgebra.LAPACK.gebal!('S', triu!(C.data, -1)) |
| 71 | + return sort!(LinearAlgebra.LAPACK.hseqr!('E', 'N', 1, size(C,1), C.data, C.data)[3], by=reim) |
| 72 | + end |
| 73 | +end |
| 74 | +hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(triu!(C.data, -1)) |
| 75 | + |
| 76 | +""" |
| 77 | + roots(c::ChebPoly{1,<:Real}; tol=5eps, maxsize::Integer=50) |
| 78 | +
|
| 79 | +Returns a sorted array of the real roots of `c` on the interval `[lb,ub]` (the lower and |
| 80 | +upper bounds of the interpolation domain). |
| 81 | +
|
| 82 | +Uses a colleague-matrix method combined with recursive subdivision of the domain |
| 83 | +to keep the maximum matrix size `≤ maxsize`, following an algorithm described by |
| 84 | +Trefethen (*Approximation Theory and Approximation Practice*, ch. 18). The `tol` |
| 85 | +argument is a relative tolerance for dropping small polynomial coefficients, defaulting |
| 86 | +to `5eps` where `eps` is the precision of `c`. |
| 87 | +""" |
| 88 | +function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Integer=50) |
| 89 | + tol > 0 || throw(ArgumentError("tolerance $tol for truncating coefficients must be > 0")) |
| 90 | + abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol |
| 91 | + n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1) |
| 92 | + if n <= maxsize |
| 93 | + λ = hesseneigvals!(UpperHessenberg(_colleague_matrix(@view c.coefs[1:n]))) |
| 94 | + @inbounds λ .= (λ .+ 1) .* ((c.ub[1] - c.lb[1])/2) .+ c.lb[1] # scale and shift to [lb,ub] |
| 95 | + return filter_roots(c, λ) |
| 96 | + else |
| 97 | + # roughly halve the domain, constructing new Chebyshev polynomials on each half, and |
| 98 | + # call roots recursively. Following chebfun, we split at an arbitrary point 0.004849834917525 |
| 99 | + # on [-1,1] rather than at 0 to avoid introducing additional spurious roots (since 0 is |
| 100 | + # often a special point by symmetry). |
| 101 | + @inbounds split = oftype(float(c.lb[1]), 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1] |
| 102 | + |
| 103 | + # pick a fast order for the recursive DCT, should be highly composite, say 2^m |
| 104 | + order = nextpow(2, n-1) |
| 105 | + c1 = chebinterp(c, order, c.lb[1], split; tol=2tol) |
| 106 | + c2 = chebinterp(c, order, split, c.ub[1]; tol=2tol) |
| 107 | + return vcat(roots(c1; tol=2tol, maxsize=maxsize), roots(c2; tol=2tol, maxsize=maxsize)) |
| 108 | + end |
| 109 | +end |
0 commit comments