Re: [Help-glpk] Applying a threshold to the solution using GMPL?

[Michael Hennebry]

On Sat, 9 Nov 2013, Reginald Beardsley wrote:

> I'm solving Ax=b using an L1 norm.  In some cases I get a number of relatively 
> small values in x that I would like to suppress based on the fraction of the 
> result that they contribute so as to make the solution more sparse.

> From this description, doing what you seem want
in one swell foop would require binary variables
and might be difficult to express even with them.
The desired range is discontinous.
A two-stage process is in order.

> Unfortunately I'm unable to recognize from the examples in the distribution 
> (4.45) or the discussion in the AMPL book how to implement this using GMPL.  
> The following model file shows what I'm trying to do.  I'd like to apply a 
> constraint such that:
>
> (sum{i in I}(A[i,j,k]*X[j,k]/b[i]))/ii >= threshold

If your L1 norm was the right call before, you should use it again.

On the second stage, put an upper limit on the original L1 norm
and minimize a weighted sum of X's:

> # array limits
> param ii;
> param jj;
> param kk;
>
> # array indices
> set I := 1..ii;
> set J := 1..jj;
> set K := 1..kk;
>
> param b{I} ,>=0;
>
> param A{I,J,K} ,>=0;
>
> # solution variables
> var X{J,K} ,>=0;
>
> # slack variables
> var u{I} ,>=0;
> var v{I} ,>=0;
>
> # objective

new objective:
sum{j in J, k in K} w[j,k]*X[j,k]
where w[j,k]= sum{i in I}b[i]/(A[i,j,k]*X1opt[j,k])

> minimize error: sum{i in I}(u[i]+v[i]);

new constraint, make the new L1 norm no more than twice as bad:
sum{i in I}(u[i]+v[i]) <= obj1opt*2

> #constraints
> s.t. eq{i in I}:sum{j in J} (sum{k in K}A[i,j,k]*X[j,k])
>                            + u[i] - v[i] = b[i];
> solve;
>
> # solution output
>
> for{k in K}{
>   printf{j in J:X[j,k]>0} "Tst: %15.8g %8.6f \n"
> }
>
> end;

I suspect that the stage 1 output could be the stage 2 input.
A more complex version of roughly the same thing is to
combine the objective functions of stage 1 and stage 2.
I'm not all that sure that either would work well.
One might get a lot of cases of close, but no cigar.

Using the API might be simpler.
Select each X in order.
    Force it to zero.
    Check to see whether the L1 norm is too big.
    If it is, unforce the X.

As you might need to sort the X's, the API might be simpler.
Sorting can be done with GMPL, but you might not like it.

--
Michael   address@hidden
"On Monday, I'm gonna have to tell my kindergarten class,
whom I teach not to run with scissors,
that my fiance ran me through with a broadsword."  --  Lily