Help regarding simple, linear minimization problem

Dear all,

I am relatively new to MATLAB, so first of all I would like to thank you in advance for your active participation and support in this forum. I have been coding in R and Python and I switched to MATLAB last week, so I'm confident with the coding, though sometimes I have some hard time to cope with understanding some specific packages like the optimization one. I would really appreciate your help in my problem which I state in what follows:

I have two variables, called ' μ' and ' σ' and I would like to solve the following minimization problem. So, in these variables (μ,σ) for a unit (say 'io') I would like to find what is the minimum value δ in order for the equation (1) to apply, so that given the variables α,β which can vary (with α,β>0 and a+b=1) α*μ io+β*σ io + δ to be equal or larger than α*μ i-β*σ i (for every i~=io).

min δ α*μ io+β*σ io + δ >= α*μ i-β*σ i (for every i~=io) s.t α,β>0 α+β=1

I know that this is possible and I assume that it is very easy, but looking at MATLAB's documentation regarding the linprog command, I was lost with how to place these equations and the constraints e.g. A, b, Aeq beq etc. in that form.

Understandably, I would like to run this problem for every i, but I can do that in a loop afterwards. It's the problem formulation that I find hard to do in a recognizable way for linprog.

I would like to thank you all very much in advance once more and I'm looking forward to your comments!

All the best, M.

Answer 1

[John D'Errico]

Sorry, but your question is not very clear. Partly because you are being sloppy in your notation. You have a,b,alpha,beta. Sometimes you use a, sometimes alpha. Are a and alpha the same? I assume so, the same for b, beta.

Are u and sigma vectors? I'd probably assume they are vectors, both of length n.

Now it sounds like for some given index j, you want to find alpha,beta,delta, that satisfy a given set of inequalities. So there are three unknowns, alpha,beta,delta? One linear constraint? So start with

Aeq = [1 1 0];
beq = 1;

Your objective f is the vector

f = [0 0 1];

So when linprog takes the dot product of f with the unknowns [alpha,beta,delta], it gets delta.

Lower bounds?

lb = [0 0 0];

Upper bounds?

ub = [1 1 inf];

Now, assume that u and sigma are column vectors. Given i0, this should get you close:

i = setdiff(1:n,i0);
A = [u(i) - u(i0),sigma(i0) - sigma(i),ones(n-1,1)];
b = zeros(n-1,1);

abd = linprog([0 0 1],A,b,Aeq,beq,lb,ub);

Without any vectors u and sigma as an example, its hard to help you more.

[Menelaos]

Dear John,

First of all thank you very much for your prompt and helpful response, and please accept my apologies for not being specific enough.

The two variables, μ and σ are 10x1 column vectors. Alpha (α), beta (β) and delta (δ) are essentially varying numbers (1x1 each) just to satisfy these inequalities, where α>0, β>0, δ is unbounded (any value -inf to inf) and α+β=1.

My question is what is the minimum value δ to satisfy the following model:

α*μ_io - β*σ_io + δ >= α*m_i - β*σ_i

where μ_io is let's say is the μ(1,1) in the first loop and thus μ_i are μ(2,1) to μ(10,1) and similarly for σ_io and σ_i.

To give you a better picture, those are sub-samples of the actual data I want to optimize. So for instance, in the first case, I want to take the μ and σ of Germany, so io is Germany and I want to solve this problem:

α*μ_Germany - β*σ_Germany + δ >= α*μ_i - β*σ_i for every i

where i are all other countries except from German (your script works great in splitting io from i as you ve written it) and α,β>0 and α+β=1 (no bound on δ whatsoever)

I'm not sure whether it is clear enough now? By running your script it seems to work, but I get no feasible solution, whereas running the same problem in excel using Open Solver I get a feasible solution. Maybe it does not try to find the minimum value of δ to satisfy this problem? δ should be technically unbounded, but I expect it to be between 0 and 10 for example

Thanks again for your interest and support in advance!

All the best, M.

[John D'Errico]

Um, if alpha > 0 AND beta> 0, AND alpha + beta == 1, then I think we can apply upper bounds of 1 to both alpha and beta. Will you give me that much credit?

If you got a value of 10 for delta, then any solver you used was clinically insane. Actually, it was just not trying to minimize delta.

Just think. You take the product of numbers alpha and u. Then the product of beta times sigma. alpha and beta are ABSOLUTELY assured to lie between 0 and 1, and we know that mu and sigma are also always small numbers.

So the term (α*μ_io - β*σ_io) is small. It will almost always be less than 1 here. Depending on the values of alpha and beta, it will be small. It might even be negative. But 1 will be a conservative and certain upper bound on that term. Will you give me that too?

We can say the same thing for the right hand side, thus (α*m_i-β*σ_i). For each value of i, we know this to be less than 1. Again, will you give me that?

How small can those same terms be? Again, the same arguments suggest they will never be less than -1, for any values of alpha, beta.

Therefore, delta is NEVER large. 10 is silly. In fact, given my arguments above, 2 is a very conservative upper bound.

Next, was it really necessary to paste in a picture of your blasted data???? You could not paste in the numbers themselves, so I could just copy them into MATLAB? Instead, you force me to type them in by hand...

[Menelaos]

Dear John,

First of all, yes, we can apply upper bounds of 1 to both alpha and beta since their sum is 1.

Second, yes a delta value of 10 is too large, it was purely an example for other type of data I use, but in these specific data, delta can never be negative, moreover the solutions I got for some countries I tried in Excel's Solver are:

e.g. deltaGermany=0.031742969 / deltaFrance=0.032023 etc. and only deltaNetherlands, Austria and Norway are equal to zero (or extremely close to), but no negative delta (Which does not make sense mathematically).

Finally, I pasted the picture since I did not expect you to go through the data, but rather to give you a visualised example of how they look. Since you asked, here are the actual data

mu 0.44688 0.28922 0.51666 0.27797 0.51894 0.61086 0.47107 0.53304 0.70692 0.62844

sigma 0.044611 0.044702 0.01382 0.092989 0.078143 0.047995 0.012679 0.07094 0.072996 0.063957

Thanks again for your help!

All the best, M.

[John D'Errico]

Ok. so having now typed in your vectors of data...

mu = [0.44688; .28922; .51666; .27797; .51894; .61086 ;.47107; .53304 ; .70692; .62844];
sigma = [.044611; .044702; .01382; .092989; .078143; .047995; .012679; .07094; .072996; .063957];

They should be close, if not exact.

For everything to revolve around Germany, we have

i0 = 1;
i = setdiff(1:n,i0);

flhs = @(alpha) alpha*mu(i0) - (1-alpha)*sigma(i0);
frhs = @(alpha) max(alpha*mu(i) - (1-alpha)*sigma(i));

I've created two functions that I can plot.

flhs is the term (α*μ_io - β*σ_io). frhs is the term max(α*μ_i - β*σ_i).

Both need be a function only of alpha, since we know that beta=1-alpha.

Note that I used max in the right hand side, because we need only the unique value of delta that satisfies ALL such terms, and does so as a very minimum. So we just need to find the largest value from all of those right hand sides.

ezplot(flhs,[0,1,-0.2,1])
hold on
ezplot(frhs,[0,1,-0.2,1])
Warning: Function failed to evaluate on array inputs; vectorizing the function may speed up its evaluation and avoid the need to loop over array elements. 
> In ezplot>ezplot1 (line 488)
  In ezplot (line 144) 
grid on

Ignore the warning there. I was too lazy to fix it.

The line in blue is the left hand side term, green the right hand side.

The value of delta necessary to solve this is the smallest possible value of delta that we need to raise ANY part of the blue curve above the green curve. That looks to be roughly 0.05 here. It will happen on the left end of the curve, near alpha=0, so beta=1. I'm not sure if the lines are parallel down there, so I can't be sure exactly what value of alpha makes this happen.

Given that, now, lets return to what I wrote in my answer. It looks like I made a mistake when I wrote down A. I was missing a minus sign on one term.

Aeq = [1 1 0];
beq = 1;
f = [0 0 1];
lb = [0 0 0];
ub = [1 1 inf];

i0 = 1;
i = setdiff(1:n,i0);
A = [u(i) - u(i0),sigma(i0) - sigma(i),-ones(n-1,1)];
b = zeros(n-1,1);
abd = linprog([0 0 1],A,b,Aeq,beq,lb,ub)
abd =
  -1.5744e-16
            1
     0.031932

I was pretty close just from the plot. alpha=0, beta=1, delta=0.031932.

Just loop over i0 above to get delta for each country, saving the result in a vector. Thus, for netherlands, i0 is 3.

i0 = 3;
i = setdiff(1:n,i0);
A = [u(i) - u(i0),sigma(i0) - sigma(i),-ones(n-1,1)];
b = zeros(n-1,1);
abd = linprog([0 0 1],A,b,Aeq,beq,lb,ub)

Optimal solution found.

abd =
            0
            1
     0.001141

And delta is quite small. For austria, i0 would be 7, so I get

abd =
     0
     1
     0

Part of me wonders if it should be easily proven that alpha is always 0 at the optimum, given your data. I'm not sure it is worth the effort, but it might be of interest to you.

[John D'Errico]

deltaopt = zeros(10,1);
alphaopt = zeros(10,1);

for i0 = 1:10
  i = setdiff(1:n,i0);
  A = [u(i) - u(i0),sigma(i0) - sigma(i),-ones(n-1,1)];
  b = zeros(n-1,1);
  abd = linprog([0 0 1],A,b,Aeq,beq,lb,ub);
  deltaopt(i0) = abd(3);
  alphaopt(i0) = abd(1);
end

[(1:10)',alphaopt,deltaopt]
ans =
          1  -1.5744e-16     0.031932
          2            0     0.032023
          3            0     0.001141
          4            1            0
          5  -2.2681e-16     0.065464
          6     0.081644     0.007302
          7            0            0
          8   1.9025e-16     0.058261
          9  -1.1493e-16     0.060317
         10    1.615e-16     0.051278

So alpha was not always 0 at the optimum.

[Menelaos]

Dear John,

First of all let me thank you once more for your contribution in this! I can tell you that the script works for sure, In fact, I just tested it in the whole sample (instead of this small sub-sample) and for many years over subsequent loops accordingly and it works miracles; thus, kudos is the least I can say here!

In regard to your observation about alpha, that was my initial thought but it is not always zero at the optimum, though it is extremely close (look at the zeros after the decimal).

I would like to thank you very much once more for your help, especially as a newcomer in Matlab I couldn't ask for a better treatment.

All the best,

M.