sci.math.num-analysis #25318 (1 + 88 more) (1)+-[1]
From: Alan Horwitz \-[1]
[1] Re: efficient Chebyshev approximation?
Reto Koradi wrote:
>
> I'm trying to solve a problem that, according to the books I have handy,
> is called "discrete Chebyshev approximation":
>
> r = A * x - b max(|r_i|) = min!
>
> Where A is a given matrix of size m * n, b is a given vector of size n,
> and I'm looking for the vector x (length m). This is similar to the
> more common linear least squares problem, only that the maximum norm
> is to be minimized.
>
> I solved this by transforming it to a LP problem:
>
> sum(k=1..m)(A_ik * x_k) + eps >= b_i (i=1..n)
> sum(k=1..m)(A_ik * x_k) - eps <= b_i (i=1..n)
> min: eps
>
> This basically works, but it takes MUCH longer than I expected, I'd like
> it to be a few orders of magnitude faster...
>
> I also tried taking the solution of the least squares problem as a
> start value, and throwing the whole thing into a generic minimizer
> (Powell's method). Without much success, the problem seems to have
> tons of local minima.
>
> So is there a better way? n can easily be a few 1000, and m a 100.
> The matrix can get somewhat sparse for the really big problems, but not
> extremely so.
> --
> Reto Koradi (kor@mol.biol.ethz.ch, http://www.mol.biol.ethz.ch/~kor)
> You might want to look at "Approximation Theory" by Cheney. Chapter 2 has nume
rous algorithms
for the Maximum norm solution of inconsistent linear equations.
End of article 25318 (of 25395) -- what next? [npq]
sci.math.num-analysis #25361 (0 + 88 more) (1)+-(1)
From: pmcallis@capmark.funb.com (Pat McAllister) \-[1]
[1] Re: efficient Chebyshev approximation?
There is probably not an easier way to solve the problem than your original LP a
pproach, although I might suggest trying a more sophisticated nonlinear optimize
r. Your problems with Powell's method don't result from local minima -- the obj
ective function is convex, so there are no local minima except the global minimu
m -- but from the non-differentiability of the objective function. It is piecew
ise linear, and your generic optimization code may have trouble trying to fit "d
erivatives" and second derivati
ves whenever the path crosses the boundary between piecewise linear segments. Y
ou should at least try a code that uses analytical first derivatives (e.g., quas
i-Newton), which should give much better performance. Also, how good is your LP
code? I would think that a modern LP code would solve this problem in a matter
of a few minutes for the problem sizes you are talking about. Unless you have
a real time application, I would think that should be good enough.
--
Pat McAllister (Pat.McAllister@funb.com)
End of article 25361 (of 25395) -- what next? [npq]
sci.math.num-analysis #25376 (3 + 88 more) [1]
From: Reto Koradi
[1] Chebyshev approximation (results so far)
Some days ago I posted a query for an efficient solution for a discrete
Chebyshev approximation problem (I've been told that it's called
L-infinity problem).
The best thing I have so far is the CHEB function from netlib. I converted
it from FORTRAN to C (by hand!), if anybody wants my version just ask.
The results are very good, it's around 300 times faster than the solution
I had before (generic Simplex algorithm). As an example, a problem with
1125 equations and 8 unknowns is solved in 0.02 seconds on a DEC Alpha.
I have one problem with it: it can give negative values for the
coefficients, I would prefer to have them all >= 0.0. The generic
Simplex algorithm did that automatically. It would also be useful
to have an algorithm where you can fix the LOWER bound for the error,
and only minimize the UPPER bound.
It was also suggested to use a Remez algorithm. I tried implementing
one, but it doesn't work well. I ordered references from the library now,
and might look into it some more. But I doubt that it will be much
faster than what I have now.
Thanks to John Gerard Malecki, Bob Craig, Lee Killough and Brian ???
for their help! (Even) better solutions are still welcome.
--
Reto Koradi (kor@mol.biol.ethz.ch, http://www.mol.biol.ethz.ch/~kor)
End of article 25376 (of 25395) -- what next? [npq]
sci.math.num-analysis #25374 (2 + 88 more) ( )--[1]
From: jdadson@ix.netcom.com (Jive Dadson )
Newsgroups: sci.stat.math,sci.math.num-analysis
[1] C++ demo of Chebyshev
Date: 5 Jan 1996 04:58:32 GMT
Organization: Netcom
Lines: 135
Distribution: inet
Message-ID: <4cib5o$aeu@cloner2.ix.netcom.com>
References: <1996Jan3.222733.22531@threetek.dialix.oz.au>
NNTP-Posting-Host: ix-sc7-08.ix.netcom.com
X-NETCOM-Date: Thu Jan 04 8:58:32 PM PST 1996
Xref: math.ohio-state.edu sci.stat.math:8689 sci.math.num-analysis:25374
Here's one reason why I like C++ for numerical stuff.
Many numerical algorithms require a function as a parameter.
Furthermore, the output of many of these algorithms is really a
function. For example, lets say you have a function F, and you
want to create a Chebyshev approximation of it. You probably call the
"buildChebyshev" routine with F as one of the parameters. The
routine returns to you the Chebyshev approximation CF. At least
that's what you want. You must settle for it returning an array of
reals. To evaluate CF, you can't just call CF like any other function.
Instead you must call an "evaluateChebyshev" routine with the array of
reals as an argument. You even have to supply also the length of
the array, and the bounds in which the approximation is valid. If you
get those messed up, you lose. (If you are lucky, the program crashes.)
So far I have only described minor inconveniences. But here
is the really tragic part: You can not pass the function CF to other
routines that take functions as parameters! For example, maybe you
already have a tricky, well-tested routine that performs some
kind of convolution of two functions. You can not pass CF to it as
one of those parameters. It wants a function, not an array of reals
and an "a" and a "b". That is a bummer of catastrophic proportions.
In C++, it is very easy to do just what you want to do.
Look at this bit of code:
Cheby cExp(exp, -5, 5, 50); // Create a Cheby approx of exp
out << cExp(3) << endl; // Print the Cheby approximation
// of exp(3)
That declares a function called "cExp" that is the Chebyshev
approximation of the exponential function "exp", valid in the interval
-5, 5, and based on 50 coefficients. Notice that I call it just like
any regular function.
Here is how the magic was performed. King's X. I like to make the
following declaration:
typedef double real;
That means today I'm using double precision for reals. Tommorrow, maybe
I'll change my mind.
Okay, time in. Now the neat stuff. I have an "abstract base class"
called "RealToReal".
class RealToReal {
public:
virtual ~RealToReal() {;}
virtual real operator() (real r) const = 0;
};
There isn't any more to it. That's it. Five lines of code. In C++,
classes are lattices, like Dana Scott described back in the sixties.
The "virtual" means that this is the top, or most abstract level in a
lattice. The class RealToReal can have many sub-classes. The above
declaration says that all subclasses of RealToReal will have an
operator "()" that takes a real as a parameter, and returns a real. The
"const" means that evaluating the function does not change its
behavior. (A random number generator would not quality as a
RealToReal.) What we do now is create "concrete" classes of functions
that map the reals to the reals. One such class of functions is
Chebyshev approximations. Here is the complete declaration for class
Cheby:
class Cheby: public RealToReal {
public:
Cheby(const RealToReal &func,real a_,real b_,int n_);
Cheby(real(*func)(real), real a_, real b_, int n_);
real operator() (real r) const;
private:
Vec coeffs;
real a;
real b;
int n;
void init(RealToReal &func);
};
The first line declares that it "is a" RealToReal. ("Is a" pronounced
as one word, is part of the object-oriented lexicon.) It means in this
case that any routine that is declared to take a RealToReal as a
parameter can take a Cheby and run with it. The routine need never have
heard of a Cheby. In fact it might have been written and compiled
before Cheby was invented.
The "private" part of class Cheby keeps all the relevant
info tucked safely away together. You can not mess up and pass the
wrong "a" and "b" with an array of coefficients "coeffs" to an
"evaluateCheby" routine. They live together. You can't even
accidentally re-use them. They are "private" and safe. And besides,
there IS no routine "evaluateCheby". There's just Cheby.
The third line declares that you can create ("construct") a Cheby
using any RealToReal. For example, if you have a kernel system that
maps reals to reals, you can construct the ChebyChev approximation of
it with a few keystrokes, provided only that the kernel system "is a"
RealToReal. What? You wrote your kernel system before you knew it
should be a "RealToReal"? No problem. In a few keystrokes you can
create a class that "is a" RealToReal and "is a" kernel system also,
and use that!
A few minutes ago, I wrote the complete class Cheby. It took
considerably less time to do (cribbing from "Numerical Recipes in C")
than it has taken to write this letter. About 15 minutes total. It is
90 lines long. (This letter is pushing 130.) To test it, I wrote a
little program that creates the Chebyshev approximation of the
exponential function in the interval -5, 5. Then I created the
Chebyshev approximation of THAT in the interval -4.5, 4.5. (Not
something you would generally want to do, I guess.) I wrote out x,
exp(x) pairs to plot and convince myself that it worked. Sure did.
First time. Here is the complete test-program:
int main()
{
ofstream out("c:\\stats\\out.dat"); // a file to plot with gnuplot
Cheby cExp(exp, -5, 5, 50); // Create Cheby approximation of exp
Cheby c2Exp(cExp, -4.5, 4.5, 50); // Create an approximation of THAT.
// Evaluate the second Chebyshev approximation a lot of times ...
real delta = .1;
for(real x = -4; x<4; x+= delta) {
out << x << " " << c2Exp(x) << endl; // write x exp(x) to output
}
return 0;
}
I think that is totally cool.
End of article 25374 (of 25395) -- what next? [npq]