inverse of erfcx (scaled complementary error function)

could anyone help me to find inverse of erfcx (scaled complementary error function)?

Answer 1

[John D'Errico]

As I recall, erfcx solves the problem that erfc goes to zero roughly like exp(x^2) for large x. For large x, erfcx is roughly

1/(x*sqrt(pi))

As we can see, this is true.

ezplot(@(x) 1./(x*sqrt(pi)),[1,10])
hold on
ezplot(@erfcx,[1,10])

And for negative x, erfcx blows up, and does so extremely fast. So you can only seriously be interested in inverting values of erfcx less than or equal to 1, where x will be non-negative.

So to invert y=erfcx(x) for a given value of y, a trivial solution is to use fzero as long as y lies in the interval [0.25,1].

For y less than say 0.25, I'd just construct a simple Newton iteration scheme, which will be even faster than fzero here.

y0 = erfcx(3)
y0 =
          0.17900115118139

>> x = 1/(y0*sqrt(pi))
x =
          3.15187684450162

>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99324120463763

>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99998665644266

>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99999999994798

>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
                         3

>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
     3

Quadratic convergence can be so nice. In fact, a quick test shows this even works quite well for larger values of y, so you need not even bother with fzero over part of the range.

In fact, if you are really lazy, don't even bother with a convergence test. Just run the Newton loop always for about a half dozen iterations for any value of y no larger then 1 and you should get full double precision accuracy.

Even better, all of this is easily vectorized.

function x = erfcxinv(y)
  % erfcx inverse, for y no larger than 1.
  x = 1./(y*sqrt(pi));
  for n = 1:6
    x = x - (erfcx(x) - y)./(2*x.*erfcx(x) - 2/sqrt(pi));
  end
end

That should work pretty cleanly. No checks to ensure that y is not greater an 1, in which case the 6 Newton iterations will probably not be sufficient, and one would want to have a tolerance test for convergence.

erfcxinv([1 .5 .25 .125])
ans =
     -1.59466173914173e-17         0.769079771061314          2.05154342412484           4.4052266527517
>> erfcx(ans)
ans =
     1  0.5  0.25  0.125

[John D'Errico]

PLEASE do NOT add an answer every time you want to make a comment. I've moved your answer into a comment. Just click on the link that lets you make a comment.

Thanks John for a quick response. Correct me if I am wrong, I take your function works only for the values of Y [0 1]. Since I need the inverse of Y greater than 1 and used your function which gave me different result than the one I got it using my alternative rookie method and checked it in wolfram alpha and was close. My code has a lot of bugs. so looking for alternative code. I would appreciate if you could help me finding inverse of Y greater than 1.

Here is what i did
%%inverse of erfcx

   erfcx_inv=erfcxinv(1/Bi);% your solution, Bi=0.25.
  erfcx_inv=-2.5923

my rookie alternative

x=-3:0.001:3;

v=round(erfcx(x),2,'significant');% round up upto 2 significant decimals

Index=find(v==(1/Bi));% find index that is equal to  a value

Index_mean=round(mean(Index)); % give me more than one so used mean

erfcx_inv=-3+Index_mean*0.001;
erfcx_inv= -0.8940 (-0.8952= wolfram alpha)

[John D'Errico]

You could have said over what domain you needed to solve this. Nobody ever does. Whine. Whine. :)

As I explicitly stated, the iterative scheme I wrote is assured to yield essentially full precision for y <= 1. And the reason one typically uses erfcx is because erfc tends to underflow very rapidly for large x. So an inverse there will be highly inaccurate.

For large y though, the only thing that needs be done differently is to change the starting estimate. I let it go one more Newton iteration just to get the last bit or so, but in my quick tests, it was already done before that point.

function x = erfcxinv(y)
  % erfcx inverse, for y no larger than 1.

    % for y <= 1, use the large x approximation for a starting value.
    k = (y <= 1);
    x = zeros(size(y));
    x(k) = 1./(y(k)*sqrt(pi));

    % for y > 1, use exp(x^2) as a very rough approximation
    % to erfcx
    x(~k) = -sqrt(log(y(~k)));

    for n = 1:7
      x = x - (erfcx(x) - y)./(2*x.*erfcx(x) - 2/sqrt(pi));
    end
  end

I suppose I should have used a tolerance, and iterated only over the elements that changed in the last iteration. But it converges so well because of the decent starting values and quadratic convergence, who cares? The extra overhead necessary for a test and keeping track of which ones did not converge yet will take more time than a simple extra iteration or so.