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]