Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convert Remez algorithm to return approximants in the basis of Chebyshev polynomials of the first kind #1

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

MikaelSlevinsky
Copy link

This PR converts the package to return approximants expressed in the basis of Chebyshev polynomials of the first kind.

Advantage is better conditioning of coefficients, disadvantage is 2x slower:

Originally:

julia> using Remez

julia> @time (nc,nd,e,x) = ratfn_minimax(abs,(-1,1),34,0);
  1.261918 seconds (10.98 M allocations: 547.756 MB, 38.82% gc time)

julia> nc
35-element Array{BigFloat,1}:
  8.235814033332272170568640236994329004443663541011174921747740472739665096420163e-03
 -2.587572296130360355942489481188429445098352205569183983293097797866424183074742e-58
  1.611288605676779648712417775437950503009937547194325526952200953507314584151585e+01
  2.47069171364733380793393574903583046615551087769848485704993429067273922958783e-55 
 -8.367796860172666650270465938199765037611010689231961740560263292815387886425113e+02
 -3.353129308110700085831877900401131839260303212019805532570534617292601599738473e-53
  2.622311890996421251900283604442420603196853693432324588388626363526585076526494e+04
  1.627530022171991446982661031140207074661816507429831923222083790902496153281047e-51
 -4.756846159341440840837470455793773568141443596727679350868030019854067541386654e+05
 -3.897776109495838983513658467136516092991367596287613542219401618776343467021422e-50
                                                                                     
  9.516531782839182388388275072517402506588456392038091789532220763648784560310062e+09
  9.871633235489553377769823576451870960793634997493947111041707129150164439358101e-46
 -6.000007233385477361317902034666560168567252128686927787451801573586072239393333e+09
 -5.073625011478670491635177570270825602981991023373280885854161087447809834939939e-46
  2.549917906028978686625846833614301917665671458005505901728790309794388508891611e+09
  1.551883931785050348901004619517014917861756526076990203746742014917539796524952e-46
 -6.545295873622896870449536758741708784858064510797846197395667905393264798089009e+08
 -2.137706802999336541202926724220286802006787165249685929657853112297996137285194e-47
  7.663957060558928841372747542360797500156797900021980785523496953371201280677198e+07

With this PR:

julia> using Remez

julia> @time (nc,e,x) = remez(abs,(-1,1),34);
  2.373292 seconds (21.64 M allocations: 1.049 GB, 41.76% gc time)

julia> nc
35-element Array{BigFloat,1}:
  6.365687432172624392896774405103717782733407546999075753986094131865272413150967e-01
  8.883901568689366754394004870226839276238147621886246327034421941626848678868969e-61
  4.245154556711678131631435763533756712967409845533888803529741348638846147245724e-01
  1.836320010268516899106112733086749328458901200875359785914688828858248991243676e-60
 -8.498556224344942784857640784192156073283544528252016704622401690673614655477616e-02
  2.01368994886521494313089730124203621800738485013096419230004661789586027284053e-60 
  3.648230009526547320458040503702185990901153549021977006526384540597710387713097e-02
  2.542810171994900850123847492798942723517116822954415546921869731667242597313589e-60
 -2.031575332430068551662742183340169677631310954072101537124015341582179563289861e-02
  1.395308824764837687283901292811903584825985136828604292919782935462095773360525e-60
                                                                                     
  2.038813979407698469749337436373587874576235698178934451655926428853316080642584e-03
  2.61210184277369791668736985964943304440557406361194664666901525577854137377456e-60 
 -1.790335599474931762336064511890385211334680008923276863913551162249146954889688e-03
  1.235887224233634734317185249287094680905996300837366514504818570136503891530279e-60
  1.594706418245314196752447003851669415267394971914508827350765187302577269535231e-03
  1.395003188721555001925106938758806785087261250074939079760052218897694590536917e-60
 -1.440400823383003253615524048685604227929669235793991454776523280279664992133678e-03
  1.339983535433911429773307042323612635913205918609562377534051817069020030740909e-60
  8.922020276727770503345815863414664636548606096792478105053722127159319793853546e-03

List of changes:

  • ratfn_minimax => remez
  • poly_eval => clenshaw
  • ratfn_eval => ratfn_evalTn
  • ratfn_leastsquares now uses \ (QR) rather than creating normal
    equations
  • new method added to remez for strictly polynomial return.

…hev polynomials of the first kind

- ratfn_minimax => remez
- poly_eval => clenshaw
- ratfn_eval => ratfn_evalTn
- ratfn_leastsquares now uses \ (QR) rather than creating normal
equations
- new method added to remez for strictly polynomial return.
@simonbyrne
Copy link
Owner

Thanks for the PR.

It would be great to have this functionality, I don't want to get rid of the existing functionality: theregular monomial basis is much more convenient for special functions (which is what I used it for).

Could we do this via dispatch somehow, e.g. have a Chebyshev type?

@MikaelSlevinsky
Copy link
Author

Yes, that sounds reasonable.

The monomial basis is also a well-conditioned basis for best approximation on the unit circle (in the complex plane). From this perspective, it would also be good to have them both.

Perhaps just a keyword/arg flag such as basis=:monomial, or basis=:Chebyshev would be better for this simple dispatching? Then another flag such as domain=:interval or domain=:circle could also be added.

@simonbyrne
Copy link
Owner

Then another flag such as domain=:interval or domain=:circle could also be added.

What would the domain arg do?

@MikaelSlevinsky
Copy link
Author

The sample points for a domain=:interval would be the Chebyshev nodes.

The sample points for a domain=:circle would be equi-spaced nodes on the circle.

This would effectively allow for periodic best approximants.

@simonbyrne
Copy link
Owner

Ah right. I guess it would be nice to have this working via dispatch, but keyword arguments are probably fine for now if you think that's easier.

@MikaelSlevinsky
Copy link
Author

ApproxFun has types for domains, such as Interval, PeriodicInterval, and Circle, and types for function spaces (bases), such as Chebyshev and Taylor.

If ApproxFun were a dependency, then these types would be readily available and avoid code duplication. On the other hand, ApproxFun is under considerable development 😉.

@simonbyrne
Copy link
Owner

That makes sense: given the significant overlap, I don't think that's an unreasonable dependency.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants