Chebyshev approximation via iterated weighted least squares fits

Question

I have the task of finding a Chebyshev approximation for a time-series; I want to check different types of functions, e.g. polynomials, rational functions, harmonics, etc.

I know that the Remez algorithm solves that problem for polynomials but its application to other functions has its own challenges which incurs a lot of implementation overhead for testing different function classes; that contrasts weighted least squares approximation, whose application to various function classes is straight forward.

That led me to the

Questions:

• can Chebyshev approximations be calculated by iteratively calculating weighted least squares approximations by giving values with higher approximation error a higher weight in the next iteration?
• is there an optimal algorithm for calculating the weights for the next iteration from the errors of the current iteration when "optimal" means yielding the fastest rate of convergence?

This question is primarily aimed at linear least squares, i.e. $\min\limits_{\alpha_i}\max\limits_{x_j}\left|y(x_j)-\sum\alpha_if_i(x_j)\right|$
but answers related to other interpretations are also appreciated.

To put some flesh to the bones in reply to fedja's request for a more specific statement of the problem:

the concrete problem I have is to "reconstruct" a smooth function from almost regularly spaced samples, whose values were rounded to integer values; therefore least squares approximation isn't acceptable because individual errors may become arbitrarily large, but I know that they must not exceed $0.5$ in the problem at hand.

The functions I intend to use for approximation are linear combinations of damped harmonics, $f_{ij}(x):=e^{a_ix}\sin(b_jx)$ with fixed values for the $a_i$ and the $b_j$.

The task is then to determine the set of $\lbrace c_{ij}$ factors for the function's linear combination that yields the best Chebyshev approximation.

I'm aware that it is not guaranteed that the error bounds may still be violated for a specific set of functions, but that can be overcome by extending that set.

It may be assumed that the least squares approximation has a high number of sign-changes in the error, so that balancing the magnitude of positive and negative errors should be possible.

Answer 1

I tried several techniques and, as usual, the simplest one turn out to work the best, so I'm posting it here.

Abstractly, we just want to find the point $x$ in an affine $m$-dimensional space $E$ of $\mathbb R^n$ ($n\approx 500, m\approx 10$, if I interpreted your "five functions" remark right) with the least $\ell^\infty$ norm $\|x\|_\infty$. In other words, we start expanding the cube from the origin until it touches $E$. What we have is an oracle (your implementation of weighted least squares) that for every weight $w$ returns the point $x(w)$ in $E$ of the least $\ell^2(w)$ norm.

WARNING: For my algorithm to work, this oracle should be decent, i.e., it should return the point sufficiently close to the true minimizer. If the least square problem is nearly singular and the oracle returns some garbage instead of the good approximation to the true solution occasionally, it may cause the unpredictable behavior of my algorithm as well. So if you want to go to larger $m$, I suggest that you think really hard about how to implement LS in a numerically stable way.

If we had $m=n-1$ (the affine hyperplane case), the task would be trivial. Take some initial weight $w_0$ (say, uniform) and ask the oracle for $x(w_0)$. Then you'll know the point $x(w_0)$ on the ellipsoid at which it is tangent to $E$. But then you know $E$! So, all you need is to formulate the new oracle query such that the output would be the point with all coordinates equal in absolute value. The elementary geometric considerations show that it happens for the weight $w[k]=w_0[k]|x(w_0)[k]|$. You switch to that weight and, voila, the oracle outputs the Chebysev approximation.

The problem is that our $m\ll n$. The key point is that for the $\ell^\infty$ approximation only some $m+1$ coordinates really matter (the alternance theorem; note, however, that you don't have a Chebyshev system, so the alternance does not have to be classical, i.e., when you go from left to right, the signs do not need to flip every time, you can even get multiplicity 2 alternance points). The only issue is that we don't know which ones.

So, we start with approximately equally spaced $m+1$ points on the interval when we are talking about smooth functions (theoretically we can start with anything, but for smooth functions starting with close points may cause numerical stability problems in your oracle). We put equal weights there with $0$ weight everywhere else, request $x(w_0)$, modify the weight to $w[k]=|x(w_0)[k]|w_0[k]$ as above (it still has $m+1$ non-zero entries) call the oracle for the second time, and get a vector $x(w)$ (of length $n$) which has equal in absolute value entries in the weighted positions and some entries everywhere else. Let $M$ be the common absolute value of weighted entries. If all other entries are not greater than $M$ in absolute value, we are done: we have both a vector and a certificate for its optimality.

Otherwise let $|x(w)[q]|>M$ be the largest entry. Then we want to replace one of the weighted positions with $q$. Which one? The rule is that after replacement and doing the above procedure, we should get the vector whose entry in the replaced position is not greater in absolute value than the common absolute value of the now weighted positions. In other words, we solve the codimension 2 Chebyshev problem with $m$-dimensional subspace of $m+2$-dimensional space. I do it by the exhaustive choice ($m$ possibilities for the alternance positions). Note that the new candidate for the alternance will have the alternance value strictly greater than $M$, so we cannot do this replacement procedure forever. In finitely many iterations we are guaranteed to stop and get the true solution.

I should confess that I do not have a good theoretical upper bound on the number of iterations needed. However it practice on your suggested sample it works reasonably fast and experimenting with other functions in generic position corroborates this observation. The asymptote code is below. I cut some corners there because the code was experimental (so my oracle really sucks), but it still clearly shows that this approach is both viable and easy to implement. If you want to try it on Alberta Asymptote website, replace uuu<=100 with uuu<=2 (otherwise it will time out and that website does not show you the console output until the program terminates).

Feel free to ask questions if something is unclear :-)

srand(seconds());

int tot=0;

for(int uuu=1; uuu<=100; ++uuu)
{

int n=500; int m=9, md=m+1;
real[] x, w, mu;
for(int k=0;k<n;++k) {x[k]=-pi+2*pi*k/(n-1); w[k]=0; mu[k]=1;}
x=sort(x);
int[] ind; int r=quotient(n-1,md);
for(int k=0;k<md;++k) {w[r*k]=1/md; ind[k]=r*k;}


real[] g(real x)
{
real[] u={exp(-x),1,exp(x),exp(-x)*sin(pi*x),sin(pi*x),exp(x)*sin(pi*x),
exp(-x)*sin(2*pi*sqrt(2)*x),sin(2*pi*sqrt(2)*x),exp(x)*sin(2*pi*sqrt(2)*x),};
return u;
}

real[] c; for(int k=0;k<m;++k) c[k]=2*unitrand()-1;
real f(real x)
{
real s=0; for(int k=0;k<m;++k) s+=c[k]*g(x)[k];
return floor(s+0.5);
}

real ff(real x, real[][] C)
{
real s=f(x); real[] gg=g(x);
for(int i=0;i<m;++i) s-=C[i][0]*gg[i];
return s;
}
real[][] A, B, C;

int MM=2000; real[] wold; real Mold=10000000; int mm;
for(mm=1;mm<=MM;++mm)
{


void work()
{

for(int i=0;i<m;++i)
{
A[i]=new real[]; 
real ss=0; for(int q=0;q<md;++q) ss+=w[ind[q]]*g(x[ind[q]])[i]*f(x[ind[q]]);
B[i]=new real[]{ss};
for(int j=0;j<m;++j)
{
real s=0; for(int q=0;q<md;++q) s+=w[ind[q]]*g(x[ind[q]])[i]*g(x[ind[q]])[j];
A[i][j]=s;
}

}
real[][] IA=inverse(A);
C=IA*B;
for(int k=0;k<5;++k) 
{
real[][] dC=IA*(B-A*C); C+=dC; 
}
real[][] ACB=A*C-B;
real s=0; for(int q=0;q<m;++q) s+=abs((ACB)[q][0]);
if(s>1/10^11) write(s);
}






real M=0, Mx=0, Mxx=0; int qq=0;




void setM()
{
M=0; Mx=0; Mxx=0; qq=0; 
real[] fff;
for(int q=0;q<n;++q) fff[q]=abs(ff(x[q],C));
for(int q=0;q<n;++q) 
if(Mx<fff[q]) {Mx=fff[q]; qq=q;}
for(int q=0;q<md;++q) Mxx=max(Mxx,fff[ind[q]]);
}



void update()
{
for(int q=0; q<md;++q) w[ind[q]]*=abs(ff(x[ind[q]],C));
real ww=sum(w);
for(int q=0;q<n;++q) w[q]/=ww;
}

work(); update(); work(); setM();


for(int k=0;k<m;++k) write(C[k][0],c[k], abs(C[k][0]-c[k]));
write("################");
for(int q=0;q<n;++q) if(w[q]>0){write(uuu, mm, q,x[q],ff(x[q],C), w[q]);}
write("***** "+string(Mx)+"  "+ string(Mold)+" ******");
write(qq);



Mold=min(Mx,Mold);

bool flag=false; for(int s=0;s<md;++s) if(ind[s]==qq) flag=true;
if(flag) break;

int ss=0; real MMM=100000000; int qqold=qq;
w[qqold]=1/md; for(int s=0;s<md;++s) w[ind[s]]=1/md; wold=copy(w); int[] indd=copy(ind);
bool flag=true;
for(int s=0; s<md && flag; ++s)
{
w[indd[s]]=0; ind[s]=qqold; work(); update(); work(); setM();
if(abs(ff(x[indd[s]],C))<=Mxx){ss=s; flag=false;}
w=copy(wold); ind=copy(indd);
}

w=copy(wold); ind=copy(indd);
w[ind[ss]]=0; ind[ss]=qqold; 

}

bool fflag=true;
for(int q=0;q<n;++q) if(abs(ff(x[q],C))>=0.5) fflag=false;
tot+=mm;
if(fflag) write(string(uuu)+" success! "+string(tot/uuu)); else {write("failure!"); pause();} 
}