Would it be a solution to evaluate $L$ at a set of points $x_i$ and morph $G$ with a function $G_m$ in such a way that $G_m[x]G[x]=L[x]$ for those points? That is, with $G_m[x]$ defined as the interpolated function through the set of points ${x_i,L[x_i]/G[x_i]}$.
For example, let your two functions be:
l[x_] := PDF[GammaDistribution[5, 2], x]
g[x_] := PDF[NormalDistribution[8, (3 E^4)/(16 Sqrt[2 \[Pi]])], x]
They share the same maximum, and are somewhat similar but not closely:
Plot[{l[x], g[x]}, {x, 1, 20}, Frame -> True, Axes -> False]

Then with:
gm = Interpolation[Table[l[x]/g[x], {x, 1, 21}]];
the plots are much more similar:
Plot[{l[x], gm[x] g[x]}, {x, 1, 20}, Frame -> True, Axes -> False]
NIntegrate[l[x], {x, 1, 20}]
(*
==> 0.9705751963
*)
NIntegrate[g[x], {x, 1, 20}]
(*
==> 0.9550846595
*)
NIntegrate[gm[x] g[x], {x, 1, 20}]
(*
==> 0.9703985212
*)
In this case, the result of integrating the morphed 'el cheapo' function is within 0.02% of the value of the 'expensive' function. Of course, for this to work the functions should be sufficiently smooth and G should never get too close to zero.
The set $x_i$ could be obtained with the EvaluationMonitor
option in the integration of $G$:
Last@ Reap[
NIntegrate[g[x], {x, 1, 20}, EvaluationMonitor :> Sow[x]]]
(*
==> {{1.151189079, 1.891291464, 3.335416098, 5.384541554, 7.843511075,
10.5, 13.15648893, 15.61545845, 17.6645839, 19.10870854,
19.84881092, 1.07559454, 1.445645732, 2.167708049, 3.192270777,
4.421755537, 5.75, 7.078244463, 8.307729223, 9.332291951,
10.05435427, 10.42440546, 10.57559454, 10.94564573, 11.66770805,
12.69227078, 13.92175554, 15.25, 16.57824446, 17.80772922,
18.83229195, 19.55435427, 19.92440546}}
*)
This set is for one reason or another a meaningful set and could be used as anchorpoints for the interpolating function.
Actually, now that I come to think about it a bit more, this may not be so useful in many cases as the number of evaluations of the expensive function may not be drastically reduced. We need to evaluate the sample points and this set may be as large as the set NIntegrate
might need. So, this saves time only if we can reduce the number of points used in the interpolation function. But if we could halve that set this would save a factor $2^n$ for $n$-dimensional functions.
Method
option? Maybe it can help, but I'm not sure it's enough for your case. $\endgroup$ – FJRA Mar 14 '12 at 13:56