Numerical evaluation of Fourier transform of Daubechies scaling funct… #921

NAThompson · 2023-01-17T03:29:02Z

…ions.

mborland

General pedantic things. I'll read the paper this week so I can be more constructive.

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

NAThompson · 2023-01-18T04:26:37Z

@jeremy-murphy : Your polynomial division algorithm is amazing; super fast and it looks like every recovered digit is correct.

jeremy-murphy · 2023-01-18T08:49:25Z

@jeremy-murphy : Your polynomial division algorithm is amazing; super fast and it looks like every recovered digit is correct.

Are you serious? I actually wouldn't have expected it.

Credit for the algorithm goes to Knuth, I just implemented it.

jeremy-murphy · 2023-01-18T09:58:02Z

Anyway, if I stop being grumpy for a moment: that's great news! :)

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

NAThompson · 2023-01-21T19:19:02Z

@jeremy-murphy , @jzmaddock , @mborland , @ckormanyos : Obviously don't want to burden you guys with a problem you think is uninteresting, but I figured I'd at least try to plea for help. The accuracy here without argument promotion just seems unacceptable-for instance without arg promotion it appears that up to 8 mantissa bits are incorrect in float precision. Arg promotion leads to extreme additional computational expense.

Maybe this is a "no solutions, only tradeoffs" situation, but I'd figure I'd at least see if there's anything obvious we can do.

jzmaddock · 2023-01-21T20:07:34Z

I'll try and take a look.

NAThompson · 2023-01-21T22:05:06Z

This is an ULP plot of the 6-vanishing moment implementation; you can see it's all over the place. The red line over the top indicates that worse values have been clipped.

mborland · 2023-01-21T22:07:48Z

This is an ULP plot of the 6-vanishing moment implementation; you can see it's all over the place. The red line over the top indicates that worse values have been clipped.

Is this with or without promotion in place?

NAThompson · 2023-01-21T22:09:33Z

Is this with or without promotion in place?

No arg promotion-actually with arg promotion this would be tautologically 0 everywhere.

NAThompson · 2023-01-21T22:23:24Z

As you might expect, we do better with fewer moments, but still not up to snuff:

NAThompson · 2023-01-21T22:26:55Z

Condition number of function evaluation revealed to be fairly bad by putting the ULP envelope around it:

jeremy-murphy · 2023-01-25T03:22:44Z

@jeremy-murphy , @jzmaddock , @mborland , @ckormanyos : Obviously don't want to burden you guys with a problem you think is uninteresting, but I figured I'd at least try to plea for help. The accuracy here without argument promotion just seems unacceptable-for instance without arg promotion it appears that up to 8 mantissa bits are incorrect in float precision. Arg promotion leads to extreme additional computational expense.

Maybe this is a "no solutions, only tradeoffs" situation, but I'd figure I'd at least see if there's anything obvious we can do.

On the contrary, it sounds very interesting, and I wish I had time to help, but I don't.

Maybe only certain operations need the extra accuracy afforded by arg promotion?

jeremy-murphy · 2023-01-26T11:48:22Z

I haven't found the literature that defines the ratio for choosing an algorithm for polynomial evaluation, but I did find this reference, which may help you:

S. Graillat, P. Langlois, and N. Louvet. Algorithms for accurate, vali-
dated and fast computations with polynomials. Japan Journal of Indus-
trial and Applied Mathematics, 26(2):215–231, 2009.

The textbook I got that reference from (Handbook of Floating-Point Arithmetic, 2nd ed.) suggests to take the "Graillat–Langlois–Louvet error-free transformation" and make it a compensated summation using the same principle as Kahan's summation. I can give you more details if you're interested.

NAThompson · 2023-01-26T16:03:25Z

I can give you more details if you're interested.

Yes, I'm interested. This particular algorithm has some trouble with compensated summation, because it's actually an infinite product where each term is a small polynomial. I think I might actually need some sort of "compensated product", if that even makes sense in floating point.

NAThompson · 2023-02-04T21:49:06Z

Well the color scheme is absolutely grotesque, but it appears to be correct:

NAThompson · 2023-06-11T20:56:30Z

@jzmaddock , @mborland : I've finally gotten this working; mind giving it a look?

…ions.

example/calculate_fourier_transform_daubechies_constants.cpp

example/fourier_transform_daubechies_ulp_plot.cpp

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

mborland · 2023-06-12T07:18:38Z

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

+// we'd immediately just have to divide through by 1/sqrt(2).
+// These numbers agree with Table 6.2, but are generated via example/calculate_fourier_transform_daubechies_constants.cpp
+template <typename Real, unsigned N> constexpr const std::array<Real, N> ft_daubechies_scaling_polynomial_coefficients() {
+  static_assert(N >= 1 && N <= 10, "Scaling function only implemented for 1-10 vanishing moments.");


Is there a way to compute the constants outside of this range on an as-needed basis like calculate_constants?

Not really, especially since I want it in a constexpr std::array . . .

The calculation is rather involved.. .

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

test/fourier_transform_daubechies_scaling_test.cpp

reporting/performance/fourier_transform_daubechies_scaling_performance.cpp

test/fourier_transform_daubechies_scaling_test.cpp

test/math_unit_test.hpp

Co-authored-by: Matt Borland <matt@mattborland.com>

…ies_scaling.hpp Co-authored-by: Matt Borland <matt@mattborland.com>

…elet.

Co-authored-by: Matt Borland <matt@mattborland.com>

NAThompson force-pushed the ft_daubechies branch 5 times, most recently from 40909f0 to c4b0791 Compare January 17, 2023 04:19

mborland reviewed Jan 17, 2023

View reviewed changes

NAThompson force-pushed the ft_daubechies branch 2 times, most recently from 5c65458 to 4adbc27 Compare January 18, 2023 04:20

mborland reviewed Jan 18, 2023

View reviewed changes

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp Outdated Show resolved Hide resolved

NAThompson force-pushed the ft_daubechies branch 2 times, most recently from 702b407 to bbea451 Compare January 19, 2023 04:00

NAThompson marked this pull request as draft January 19, 2023 04:12

NAThompson force-pushed the ft_daubechies branch from bbea451 to e0373ac Compare January 21, 2023 18:58

NAThompson force-pushed the ft_daubechies branch from e0373ac to d793bbc Compare January 21, 2023 22:05

NAThompson force-pushed the ft_daubechies branch from d793bbc to 94845c8 Compare January 21, 2023 22:24

NAThompson mentioned this pull request Jan 23, 2023

Use Estrin's Scheme for polynomial evaluation #924

Closed

NAThompson force-pushed the ft_daubechies branch from 94845c8 to cbcba28 Compare January 26, 2023 05:51

NAThompson force-pushed the ft_daubechies branch 8 times, most recently from f25ec93 to 6db459e Compare February 6, 2023 05:19

NAThompson force-pushed the ft_daubechies branch from 6db459e to 96d761e Compare February 12, 2023 23:06

NAThompson force-pushed the ft_daubechies branch from 96d761e to c6d1c51 Compare March 14, 2023 04:19

NAThompson force-pushed the ft_daubechies branch 3 times, most recently from 101e952 to 841795b Compare June 11, 2023 20:54

NAThompson marked this pull request as ready for review June 11, 2023 20:55

NAThompson force-pushed the ft_daubechies branch from 841795b to 6323bed Compare June 11, 2023 23:33

Numerical evaluation of Fourier transform of Daubechies scaling funct…

5de1353

…ions.

NAThompson force-pushed the ft_daubechies branch from 6323bed to 5de1353 Compare June 12, 2023 03:53

mborland reviewed Jun 12, 2023

View reviewed changes

NAThompson and others added 8 commits June 12, 2023 08:32

Update example/calculate_fourier_transform_daubechies_constants.cpp

a4140ee

Co-authored-by: Matt Borland <matt@mattborland.com>

Update example/fourier_transform_daubechies_ulp_plot.cpp

3b995fc

Co-authored-by: Matt Borland <matt@mattborland.com>

Update include/boost/math/special_functions/fourier_transform_daubech…

774fa2b

…ies_scaling.hpp Co-authored-by: Matt Borland <matt@mattborland.com>

Update include/boost/math/special_functions/fourier_transform_daubech…

083d45e

…ies_scaling.hpp Co-authored-by: Matt Borland <matt@mattborland.com>

Rename include file to reflect it implements both the scaling and wav…

412474d

…elet.

Add performance to docs.

99c60f1

Update test/math_unit_test.hpp

c2b0d1f

Co-authored-by: Matt Borland <matt@mattborland.com>

Add boost-no-inspect to files with non-ASCII characters.

97bd84c

NAThompson merged commit 7887d43 into develop Jun 13, 2023

NAThompson deleted the ft_daubechies branch June 13, 2023 15:05

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Numerical evaluation of Fourier transform of Daubechies scaling funct… #921

Numerical evaluation of Fourier transform of Daubechies scaling funct… #921

NAThompson commented Jan 17, 2023

mborland left a comment

NAThompson commented Jan 18, 2023

jeremy-murphy commented Jan 18, 2023

jeremy-murphy commented Jan 18, 2023

NAThompson commented Jan 21, 2023

jzmaddock commented Jan 21, 2023

NAThompson commented Jan 21, 2023

mborland commented Jan 21, 2023

NAThompson commented Jan 21, 2023

NAThompson commented Jan 21, 2023

NAThompson commented Jan 21, 2023

jeremy-murphy commented Jan 25, 2023

jeremy-murphy commented Jan 26, 2023

NAThompson commented Jan 26, 2023

NAThompson commented Feb 4, 2023 •

edited

Loading

NAThompson commented Jun 11, 2023

mborland Jun 12, 2023

NAThompson Jun 12, 2023

Numerical evaluation of Fourier transform of Daubechies scaling funct… #921

Numerical evaluation of Fourier transform of Daubechies scaling funct… #921

Conversation

NAThompson commented Jan 17, 2023

mborland left a comment

Choose a reason for hiding this comment

NAThompson commented Jan 18, 2023

jeremy-murphy commented Jan 18, 2023

jeremy-murphy commented Jan 18, 2023

NAThompson commented Jan 21, 2023

jzmaddock commented Jan 21, 2023

NAThompson commented Jan 21, 2023

mborland commented Jan 21, 2023

NAThompson commented Jan 21, 2023

NAThompson commented Jan 21, 2023

NAThompson commented Jan 21, 2023

jeremy-murphy commented Jan 25, 2023

jeremy-murphy commented Jan 26, 2023

NAThompson commented Jan 26, 2023

NAThompson commented Feb 4, 2023 • edited Loading

NAThompson commented Jun 11, 2023

mborland Jun 12, 2023

Choose a reason for hiding this comment

NAThompson Jun 12, 2023

Choose a reason for hiding this comment

NAThompson commented Feb 4, 2023 •

edited

Loading