Numerically stable computation of $F((ax+b)^k)$

Let $F: \mathbb{R} \rightarrow \mathbb{R}$ be a linear map. I want to evaluate an expression of the type $$F((ax+b)^k)$$ in terms of $F(x)$ for some fixed value of $x$ (I already know $F(x)^r$ for $r=1,...,k$).

$x$ is typically small ($0<x<1$), and $k$ is typically about $15$. $a$ is always positive (about $30$) and $b$ is always negative (about $-30$) so that $ax+b$ is always between $0$ and $1$. Further it is also known that the range of the map $F$ is $[0,1]$.

Using the binomial theorem to expand $(ax+b)^k$ involves very large numbers, which cause errors due to overflow (in Matlab).

Since it is already known that the argument $(ax+b)^k$ as well as its image $F((ax+b)^k)$ are always small, I am interested to know if there are methods to evaluate $F((ax+b)^k)$ in terms of $F(x)$ stably or without involving large numbers. Any help will be much appreciated.

EDIT: The actual problem I want to solve is not quite the same. I am describing it below.

I have two functions $f,g:\mathbb{R} \to \mathbb{R}$. $g$ is positive and unit normalized, i.e. $\int_{-\infty}^{\infty} g(x)dx=1$ The function $f$ is not known, but it is known that $f(x) \in [0,1]$ for all $x$. Also the quantity $\int_{-\infty}^{\infty} g(x) f(x)^r dx$ is known for $r=1,...,k$ (which is, of course, in $[0,1]$). I want to calculate the quantity $$I = \int_{-\infty}^{\infty} g(x) (a f(x) + b)^k dx$$ where $a,b \in \mathbb{R}$. It is known that $a>0,b<0$ and $a,b$ are chosen such that $af(x)+b \in [0,1]$. Hence $I \in [0,1]$. $a$ is typically about $30$, and $k$ is about $15$.

Since only $\int_{-\infty}^{\infty} g(x) f(x)^r dx$ is known I have no option but to expand $(a f(x) + b)^k = \sum_{r=0}^k \binom{k}{r} a^r f(x)^r b^{k-r}$ using the binomial theorem. However this involves the product of large numbers like $\binom{k}{r}$ and powers of $a$ and $b$. Since $I \in [0,1]$, compuing $I$ by adding and subtracting such large nos. does not seem like a good idea.

In the original question I intended the linear map $F$ to correspond to the map $f \mapsto \int_{-\infty}^{\infty} g(x) f(x) dx$, but that was clearly not a good example.

• if $x$ and $ax+b$ vary between $[0,1]$, does this imply that a and b are always choosen such that they compensate each other? But then they cannot both be about 30. If a and b are always of the same magnitude, then $a^k F((x+b/a)^k)$ should be more feasible.
– Bort
Jul 24, 2017 at 15:08
• If $F$ is linear, doesn't that mean $F(x)=\alpha x$ for some scalar $\alpha$, which you can find out explicitly? Why are you expanding $(ax+b)^k$ if you can evaluate it directly? Can you be more explicit (preferably with code)? It would be nice to have a testcase that shows directly how what you're doing fails in matlab. Jul 24, 2017 at 17:38
• @Bort You are right, I only know that $a$ is about $30$. $b$ may be smaller, but can still be large enough ($>10$) so that the powers of $b$ are substantially large. Your suggested method is indeed better, but still involves large binomial coefficients. Anyway, my original question was quite vague. I have made it more explicit now. Jul 26, 2017 at 4:38
• @Kirill Sorry, my original question was indeed not clear. It was just a (wrongly) simplified version of my actual problem. I thought the actual problem might be too narrow in scope to get sufficient attention from the community. I have edited my question now with the real problem. Hopefully it is clearer now why I need to use binomial theorem expansion. Maybe this question is better suited to Math Stackexchange? Jul 26, 2017 at 4:45

I think your problem is ill-conditioned and that this is not a numerical issue: the condition number with respect to $F_j$ is proportional to $\binom{k}{j}a^j b^{k-j}$, which with yours numbers can be huge (e.g., $30^{15}$ for $j=k$). Think for example about how the answer would change if $F_k$ was changed by $10^{-16}$%—this has nothing to do with the exact way in which you compute that number, only how you define it. So since the number you are trying to compute is so sensitive to minute changes in the input data, this is not a well-posed problem.

Any correct algorithm, however numerically accurate, will be just as unstable because it is a property of the problem as posed. In principle, this could change if you picked a different but similar problem with a different parametrization, but not as it stands.

This is probably not an answer, more like random thoughts. We know that $$\int g(x) f(x)^kdx=F_k$$. Let's assume you can invert $y=f(x)$, we get $$\int g(f^{-1}(y)){f^{-1}}^\prime(y) y^k dy =F_k$$ Hence, we can interpret $F_k$ as the moments of the function: $$g(f^{-1}(y)){f^{-1}}^\prime(y) = \sum_k \frac{y^k}{k!}F_k$$ which is $$\frac{g(x)}{f^\prime(x)} = \sum_k \frac{f(x)^k}{k!}F_k$$ .

Now, let's do the same exercise with your shifted form: $$\int g(x) (a f(x)+b)^kdx$$ with $y(x)=a f(x)+b$ and $x(y)=f^{-1}\left(\frac{y-b}{a}\right)$. Hence, $\frac{dx}{dy}={f^{-1}}^\prime\left(\frac{y-b}{a}\right)\frac{1}{a}$.

Again we get:

$$\int g\left(f^{-1}\left(\frac{y-b}{a}\right)\right){f^{-1}}^\prime\left(\frac{y-b}{a}\right)\frac{1}{a} y^k dy =\hat{F}_k$$.

With the same interpretation: $$g\left(f^{-1}\left(\frac{y-b}{a}\right)\right){f^{-1}}^\prime\left(\frac{y-b}{a}\right)\frac{1}{a}=\sum_k \frac{y^k}{k!}\hat{F}_k$$ $$\frac{g(x)}{a f^\prime(x)} = \sum_k \frac{\left(a f(x)+b\right)^k}{k!}\hat{F}_k$$

Since both series should be equal in powers of $f(x)^k$: $$a \sum_k \frac{\left(a f(x)+b\right)^k}{k!}\hat{F}_k=\sum_k \frac{f(x)^k}{k!}F_k$$ we have a relation between $\hat{F}_k$ and $F_k$.

Is this easier to compute? Is this even correct? I don't know, but maybe it's helpful.