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

Question

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.

Answer 1

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.

Answer 2

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.