Re: Polynomial Approximation L-inf norm

Srikanth <skt@xxxxxxxxxx> wrote in message <d97a9378-655a-42a6-
8f71-ec6a38c8792c@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>...
Hi
Thanks for the info - I tried implementing the solution, however I
noticed something that I could you some help on:

I wanted to approximate 1/(x+a) by a polynomial function, when a is a
constant. I used the L-inf norm, and linprog as suggested above, by
using two extra variables, and writing it as a constrained linear
optimization. However, my domain of x varies from 0 to 10. If I take a
to be much much larger than x, I get a good approximation. But when a
is less than x, I get horrible results (and since I use a finite
length polynomial, I seem to be getting a worse result as the power of
x increases, since my last term becomes much bigger than my preceding
terms as the number of terms increases).
I asked on sci.math, and I was shown a Maple example that gives a good
approximation for a similar problem (15th order approximation). I
wonder if there is something I am fundamentally doing wrong, or some
unwarranted assumption I am making, that is causing this problem.

http://groups.google.com/group/sci.math/browse_thread/thread/ff8d75268
826fee6/7bb105c21921d08e#7bb105c21921d08e
(the link to sci.math post)

My code terminates with the following error -
Exiting: One or more of the residuals, duality gap, or total relative
error has stalled:
the primal appears to be infeasible (and the dual unbounded).
(The dual residual < TolFun=1.00e-008.)

Here is a snippet of the code I wrote... I know that some parts might
not appear coded well, but that is because I just made some quick
changes to the actual code to illustrate my point:

x=(1:.01:50)';
k=1;
a=10;
N=15; %The maximum power i.e. the order of the polynomial
b=k./(x+a);
A = bsxfun(@power,x,0:1:N);
UB=[];
LB = [-inf*ones(1,length(0:1:N)) 0 0]';
A_LP=[A -1*ones(length(b),1) zeros(length(b),1); -A zeros(length(b),1)
-1*ones(length(b),1)];
b_LP=[b;-b];
f = [zeros(1,length(0:1:N)) 1 1]';
soln=linprog(f,A_LP,b_LP,[],[],LB,UB);

Thanks a million for your help so far... I really appreciate it.
Srikanth

I'll add that you should always think carefully
about what you are doing when you try these
things.

Here, x lives in the interval [0,50], and you are
using a 15th order polynomial. So you must
use numbers as large as 50^15. This results
in numerical trash, when much of the rest of
your array is composed of smaller numbers.

Next, when you work with high order
polynomials, don't use an interval like [0,50],
or even [0,1]. Instead use [-1,1]. Lets look at
the difference.

x = [0:.01:50]';
A = bsxfun(@power,x,0:15);

Plot the columns of A. See that the last
column stands out here.

plot(A)

What is the numerical rank of A?

rank(A)
ans =
5

Yes, it is only 5. As I said, using the [0,50]
interval with high order polynomials is a killer.
Now try it with some other intervals.

x = linspace(0,1,5001)';
A = bsxfun(@power,x,0:15);
rank(A)
ans =
16

Amazingly, the matrix is now numerically
a full rank matrix, yet we know that nothing
really changed except the scaling. Look at
the condition number, to see more clearly
how good or bad it is.

cond(A)
ans =
1.4179e+11

The condition number of A is the ratio of
the largest to the smallest singular value.
When it gets large, things get bad. So
here, even though it is not quite as large
as 1/eps, it is still rather large. Can we do
better? Yes.

x = linspace(-1,1,5001)';
A = bsxfun(@power,x,0:15);

rank(A)
ans =
16

cond(A)
ans =
2.2946e+05

You can understand the differences by
plotting the columns of A. See how much
alike they look to each other when x lives
on the interval [0,1]. See that they are more
different looking when the interval changes
to [-1,1].

Next, as Bruno points out, 5001 points
is probably overkill. You might consider
other sampling strategies too, such as a
Chebychev sampling.

The one thing that you should avoid doing
is to just throw some lines of code together
and then run the whole mess without
understanding what you are doing. Go
slowly, looking at each piece, making sure
that each step is correct before you proceed.
Think about the numerical analysis, and if
you don't understand what just happened
to you, dive into to it until you do.

John


.

Relevant Pages

  • Re: accurate polynomial approximation
    ... Interpolation makes use of information ... which should result in such a high precise approximation cannot be ... What concerns the piecewise use of least squares polynomials: ...
    (sci.math.num-analysis)
  • Re: polynomial approx. of natural logarithm
    ... Peter Spellucci wrote: ... for log on there is the chebyshev polynomial approximation ... where T_k are the chebyshev polynomials of the first kind and a_k ... Why not leave that series expansion in Tchebycheff polynomials, ...
    (sci.math.num-analysis)
  • Re: Fast exponent and logarithm, given initial estimate
    ... > take a polynomial approximation. ... > rather low precision, ... > series in chebyshev polynomials of the first kind (hopefully avaliable ... recover the actual exponent. ...
    (sci.math.num-analysis)
  • Re: Questions on polynomial interpolation
    ... Characterization of the Chebyshev spline of best approximation in nonsymmetric $L\sb 1$ norm with the positive weight for a class of continuous functions. ... Chebyshev approximation by spline functions with free knots. ... Uniqueness of the Chebyshev spline of the best approximation in the mean for the class of functions having no more than one discontinuity point. ... which is much easier to do with polynomials. ...
    (sci.math.num-analysis)
  • ? pure lucky or something deep inside when estimating a matrix
    ... both on eigenvalues and eigenvectors. ... quantitatively way to measure how good this approximation is. ... Or could there be a probablity estimate to this approach? ...
    (sci.math)