Re: [Help-glpk] Non-sparse problems and other questions

[Erik de Castro Lopo]

On Fri, 31 May 2002 12:41:49 +0400
"Andrew Makhorin" <address@hidden> wrote:

> 1000 rows and 3000 columns with dense matrix this is cool :+)

I thought that might be pushing the envelope a little :-).

This amount of data is required to generate a symmetric FIR filter with
2000 coefficients (ie one side alone is 1000 coefficients). I'm aiming
to be able to design filters of 20000 coefficients, just to push the 
envelope a little further ;-}.

<snip a good explanation>

> I would suggest you to use new api routines introduced in glpk 3.1.

OK, I'll do that. I have been using the version of GLPK in Debian Woody
(3.06) but I'll remove that and try to track the latest version.

> >    gm_scaling: max / min = 1.476e+18
> >    gm_scaling: max / min = 1.476e+18

<snip>

> I'd like to draw your attention on the first two lines. They mean that
> max|a[i,j]|/min|a[i,j]| = 1.5e+18, where a[i,j] is non-zero constraint
> coefficients in your problem, i.e. the constraint matrix is extremely
> badly scaled, due to that the lp solver routines are unable to provide
> any appropriate accuracy during computations. That's why you see the
> message "Numerical instablity".
> 
> I'm not an expert in FIR filters design, but I think there is something
> wrong in the problem formulation. 

The problem is that the large majority of the a[i,j] are generated by

         a[i,j] = cos ( some_function (i, j) * PI);

so although the max|a[i,j]| is <= 1.0, min|a[i,j]| could well be as small
as 1.5e-18 or worse.

> I believe there should be another approach to formulate the problem in 
> some other way. 

I've got a feeling that poor scaling is unavoidable part of this problem.
If you are interested, I am implementing the FIR design procedure explained
in this paper :

    http://www.music.Princeton.edu/classes/meteor_article.ps

The paper descibes how to use LP to design FIR filters in about 2 pages. 
Source code to an FIR filter design program using this method is also 
available but was written in pascal (there was also a version that had been
automatically converted to C using p2c) and used a very poor LP solver.

> Or the methoda you
> use and that working well for small sized problems can't be directly
> generalized for problems of such high sizes. 

That is usually the case with the FIR filter design problem. The classic
algorithm for this is McClellan-Parks-Rabiner algoithm (also known as
"Remez exchange") which is unable to design filters longer than about 300
coefficients due to numerical instability problems. 

Using GLPK I have designed filters of 1000 coefficients with nothing more
that a few "Numerial instability" warnings :-). The resulting filters were
correct and optimal. In fact, using glp_call_rsm1 I did not even get any
warnings.

Anyway, I'll have a look a version 3.1

Thanks for your help.

Regards,
Erik
-- 
+-----------------------------------------------------------+
  Erik de Castro Lopo  address@hidden (Yes it's valid)
+-----------------------------------------------------------+
     Customer:  "I'm running Windows '95.
Tech. Support:  "I see."
     Customer:  "My computer isn't working now."
Tech. Support:  "Yes, you already said that."