How do I compare the 2-dimensional "faces" of a 3-dimensional matrix?

I'm solving a problem which has several solutions, but each iteration of my algorithm will output a single solution as a 2-dimensional (2 by 3) matrix. I wrote a for-loop into the algorithm so that it outputs a set of solutions as a matrix which I've called "sols", so "sols(2, 3, 100)", say, will be the output when I cycle through 100 iterations and so will display the 100 solutions as various 'faces' of the matrix.

The problem is, there are only several solutions to the problem, not all the faces of the matrix "sols" are distinct. My algorithm will output this three dimensional matrix containing the 100 solutions, but not all of 100 solutions will be different. What I'd like to do is write something into the code to then take the matrix sols(2,3,100) and create another output which only displays the distinct 2-dimension matrices (faces) of the 3-dimensional matrix.

I need the algorithm to either output the distinct 2-dimensional matrices or just the first rows of each 2-dimensional matrix, if that's any easier. A command like "unique(sols) came to mind, but of course that just outputs the distinct elements in each matrix, which is meaningless to me as I'd require either the distinct top or bottom rows of each matrix or the distinct 2-dimensional (2 by 3) matrices.

How could I do this?

.

Here's my code below, though you'd require YALMIP to run it:

sdpvar z1 z2 x(2,2,'full') t(2,2,'full') g(2,2,2,2,'full') sols(2,2,100,'full');

vc1 = [0 <= x]; vc2 = [0 <= t];  vc3 = sum(x) == 1; vc4 = sum(t) == 1;
F = [vc1, vc2, vc3, vc4];

g(:,:,1,1) = [3/4,3/4; 1,0]; g(:,:,2,1) = [3/4,1; 3/4,0];
g(:,:,1,2) = [3/8,3/8; 1,0]; g(:,:,2,2) = [3/8,1; 3/8,0];

eq1 = z1 - x(:,1).'*g(:,:,1,1)*x(:,2) <= 0;
eq2 = z1 - x(:,1).'*g(:,:,1,2)*x(:,2) <= 0;
eq3 = z2 - x(:,1).'*g(:,:,2,1)*x(:,2) <= 0;
eq4 = z2 - x(:,1).'*g(:,:,2,2)*x(:,2) <= 0;

eq5 = [[1,0]*g(:,:,1,1)*x(:,2), [1,0]*g(:,:,1,2)*x(:,2)]*t(:,1) - z1 <= 0;
eq6 = [[0,1]*g(:,:,1,1)*x(:,2), [0,1]*g(:,:,1,2)*x(:,2)]*t(:,1) - z1 <= 0;
eq7 = [x(:,1).'*g(:,:,2,1)*[1;0], x(:,1).'*g(:,:,2,2)*[1;0]]*t(:,2) - z2 <= 0;
eq8 = [x(:,1).'*g(:,:,2,1)*[0;1], x(:,1).'*g(:,:,2,2)*[0;1]]*t(:,2) - z2 <= 0;

G = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8];
objective = 0;
sol = optimize(F+G, objective);

r = recover(depends([F+G]));

for i=1:100

randomobjective = randn(length(r),1)'*r;

sol = optimize([F+G, objective >= value(objective)], randomobjective);

sols(:,:,i)=double(x);

end

sols

[Johan Löfberg]

Your claim that the set is polyhedral is not obvious after performing some eliminations. Use the simplex-structure of x and t to remove redundant parametrization

Reduced = replace(replace(F+G,x(2,:),1-x(1,:)),t(2,:),1-t(1,:))

This set still contains bilinear stuff (although less than in the original model).

>> sdisplay(sdpvar(Reduced(3:6)))

ans = 

      '-z1+0.7500*x(1,1)+x(1,2)-x(1,1)*x(1,2)'
      '-z1+0.3750*x(1,1)+x(1,2)-x(1,1)*x(1,2)'
      '-z2+x(1,1)+0.7500*x(1,2)-x(1,1)*x(1,2)'
      '-z2+x(1,1)+0.3750*x(1,2)-x(1,1)*x(1,2)'

Answer 1

[Robert Rouse]

I've resolved this issue now - by redefining the sols matrix and changing the last few lines to produce a 2-dimensional matrix I can use the unique(sols,'rows') command to output the unique rows of the solutions, which is all I require as the other set of 2nd rows of the original 2-dimensional matrices are determined by the corresponding 1st rows.

    sols(i,:)=double(x(1,:));
    end
    C = unique(sols,'rows')

Thanks.

[Johan Löfberg]

But in no sense are you theoretically guaranteed to have computed points on the boundary of the feasible set, as it is nonconvex and fmincon just as well can get stuck at basically any point

If it was convex, and you wanted to do ray-shooting to generate points on the feasible set, you should remove the objective >= value(objective), which was something I added when I though you actually had an objective, and wanted to generate new points with the same objecvtive

[Robert Rouse]

I have tested the points the algorithm outputs and they are in the equilibrium I require. Perhaps it's just happy coincidence. On removing that obj >= v(obj) it doesn't seem to affect any of the solutions at the moment, but again, as this example problem is fairly basic, perhaps it's just happy coincidence.

Thanks for all this help!

Answer 2

[Johan Löfberg]

Since your initial solve is simply a feasibility problem with the objective 0, what you are doing in the for-loop is that you are generating some weird vertices of the feasible set, where every new vertex has a larger inner product with the random direction, than the previous combination. This looks really weird. (in your mail to me, I didn't note that objective was 0)

If your goal really is to generate all optimal solutions (for some undefined objective here), you are pretty much doomed using this approach, as the feasible set is nonconvex and a ray-shooting procedure will not be applicable (unless you for some reason now that the set actually is convex, and that the nonlinear solver you use actually will be able to solve the seemingly nonconvex problem to global optimality)

If you are trying to generate vertices of the feasible set, you are also doomed using this approach (for the same reason)

Perhaps you could look at the vertices of the convex hull instead, but most likely you will have to analyze it analytically

[Robert Rouse]

In fact I'm looking for equilibria in a game theory problem known as the Public Good Game, or Free Rider Problem - only in this instance I'm using a formulation where there's uncertainty in the payoff matrix, but these uncertain payoff matrices belong to a bounded polyhedral set. The payoff matrix depends on a value $c$, and the algorithm above takes the (in this case 2) extreme values of the bounded polyhedral set and calculates an equivalent form of Nash equilibria using them. The formulation's described in M. Aghassi & D. Bertsimas "Robust Game Theory" [2006].

That's why the for loop you wrote works so well for me, and outputs multiple robust equilibria.

The only problem is I'd now like a way to compare all these 2x3 matrices representing each equilibrium point, as they're currently displayed in a 2x3x"n" matrix and looking at them one by one and picking out the unique ones isn't as efficient as I'd like.

[Johan Löfberg]

It is not a polyhedral set. If it was polyhedral, then YALMIP would select an LP solver, but you get fmincon, as it is composed of bilinear constraints. At least, it does not have a polyhedral representation, so a nonlinear solver might fail miserably in a ray-shooting prcedure

[Johan Löfberg]

BTW, the definition of g and sol as sdpvars is completely redundant (you don't use them as such)

[Robert Rouse]

Well, the procedure does output all the answers I require and that's all I need at this point.

I'm afraid I'm completely new to MATLAB and of course YALMIP, I started using it to input this problem and as I have very little understanding of it all I simply listed everything I'd be using on the sdpvar row as a memory aid. I can appreciate that both g and sols aren't necessary additions - thanks for pointing it out!

That said, I do use the sols vector after applying the algorithm to multiple problems, some of which use different dimensions to others. Using sdpvar gets round the "dimensional mismatching"

Answer 3

[Robert Rouse]

It does seem to fall down when considering playing the game with 3 players instead of two. When it outputs the solutions not all of them are feasible, and MATALB produces several warning messages about the matrix being close to singular. Here's an example code with the three player public good game, for interest's sake:

sdpvar z1 z2 z3 x(2,3,'full') t(2,3,'full') g(2,2,2,3,2,'full') sols(2,3,10,'full');
%g takes the form (2x2x2 payoff, to player #, l-value)

vc1 = [0 <= x]; vc2 = [0 <= t];  vc3 = sum(x) == 1; vc4 = sum(t) == 1;
F = [vc1, vc2, vc3, vc4];

g(:,:,1,1,1) = [3/4,3/4; 1,1];     g(:,:,2,1,1) = [3/4,3/4; 1,0]; 
g(:,:,1,2,1) = [3/4,1; 3/4,1];     g(:,:,2,2,1) = [3/4,1; 3/4,0]; 
g(:,:,1,3,1) = [3/4,3/4; 3/4,3/4]; g(:,:,2,3,1) = [1,1; 1,0]; 
g(:,:,1,1,2) = [3/4,3/4; 1,1];     g(:,:,2,1,2) = [3/4,3/4; 1,0]; 
g(:,:,1,2,2) = [3/8,1; 3/8,1];     g(:,:,2,2,2) = [3/8,1; 3/8,0]; 
g(:,:,1,3,2) = [3/8,3/8; 3/8,3/8]; g(:,:,2,3,2) = [1,1; 1,0]; 

eq1 = z1 - [x(:,1).'*g(:,:,1,1,1)*x(:,2), x(:,1).'*g(:,:,2,1,1)*x(:,2)]*x(:,3) <= 0;
eq2 = z1 - [x(:,1).'*g(:,:,1,1,2)*x(:,2), x(:,1).'*g(:,:,2,1,2)*x(:,2)]*x(:,3) <= 0;
eq3 = z2 - [x(:,1).'*g(:,:,1,2,1)*x(:,2), x(:,1).'*g(:,:,2,2,1)*x(:,2)]*x(:,3) <= 0;
eq4 = z2 - [x(:,1).'*g(:,:,1,2,2)*x(:,2), x(:,1).'*g(:,:,2,2,2)*x(:,2)]*x(:,3) <= 0;
eq5 = z3 - [x(:,1).'*g(:,:,1,3,1)*x(:,2), x(:,1).'*g(:,:,2,3,1)*x(:,2)]*x(:,3) <= 0;
eq6 = z3 - [x(:,1).'*g(:,:,1,3,2)*x(:,2), x(:,1).'*g(:,:,2,3,2)*x(:,2)]*x(:,3) <= 0;

eq7 = [[[1,0]*g(:,:,1,1,1)*x(:,2), [1,0]*g(:,:,2,1,1)*x(:,2)]*x(:,3), [[1,0]*g(:,:,1,1,2)*x(:,2), [1,0]*g(:,:,2,1,2)*x(:,2)]*x(:,3)]*t(:,1) - z1 <= 0;
eq8 = [[[0,1]*g(:,:,1,1,1)*x(:,2), [0,1]*g(:,:,2,1,1)*x(:,2)]*x(:,3), [[0,1]*g(:,:,1,1,2)*x(:,2), [0,1]*g(:,:,2,1,2)*x(:,2)]*x(:,3)]*t(:,1) - z1 <= 0;
eq9 = [[x(:,1).'*g(:,:,1,2,1)*[1;0], x(:,1).'*g(:,:,2,2,1)*[1;0]]*x(:,3), [x(:,1).'*g(:,:,1,2,2)*[1;0], x(:,1).'*g(:,:,2,2,2)*[1;0]]*x(:,3)]*t(:,2) - z2 <= 0;
eq10 = [[x(:,1).'*g(:,:,1,2,1)*[0;1], x(:,1).'*g(:,:,2,2,1)*[0;1]]*x(:,3), [x(:,1).'*g(:,:,1,2,2)*[0;1], x(:,1).'*g(:,:,2,2,2)*[0;1]]*x(:,3)]*t(:,2) - z2 <= 0;
eq11 = [[x(:,1).'*g(:,:,1,3,1)*x(:,2), x(:,1).'*g(:,:,2,3,1)*x(:,2)]*[1;0], [x(:,1).'*g(:,:,1,3,2)*x(:,2), x(:,1).'*g(:,:,2,3,2)*x(:,2)]*[1;0]]*t(:,3) - z3 <= 0;
eq12 = [[x(:,1).'*g(:,:,1,3,1)*x(:,2), x(:,1).'*g(:,:,2,3,1)*x(:,2)]*[0;1], [x(:,1).'*g(:,:,1,3,2)*x(:,2), x(:,1).'*g(:,:,2,3,2)*x(:,2)]*[0;1]]*t(:,3) - z3 <= 0;

G = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12];
objective = 0;
sol = optimize(F+G, objective);

r = recover(depends([F+G]));
for i=1:10
randomobjective = randn(length(r),1)'*r;
sol = optimize(F+G, randomobjective);
sols(:,:,i)=double(x);
end
sols

[Johan Löfberg]

It looks like you are trying to look at the convex hull of two polytopes or something. Is that what you are doing?

[Johan Löfberg]

Oh, it's cubic now. Even worse now then, and a nonlinear programming approach seems even further away from being a sensible way to enumerate the feasible set.