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

3 次查看（过去 30 天）
Robert Rouse
Robert Rouse 2015-4-7
评论： Johan Löfberg 2015-4-9
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
  1 个评论
Johan Löfberg
Johan Löfberg 2015-4-9
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)'

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Robert Rouse
Robert Rouse 2015-4-9
编辑：Robert Rouse 2015-4-9
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.
  2 个评论
Johan Löfberg
Johan Löfberg 2015-4-9
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
Robert Rouse 2015-4-9
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!

请先登录，再进行评论。

更多回答（2 个）

Johan Löfberg
Johan Löfberg 2015-4-9
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
  4 个评论
显示 2更早的评论
Johan Löfberg
Johan Löfberg 2015-4-9
BTW, the definition of g and sol as sdpvars is completely redundant (you don't use them as such)
Robert Rouse
Robert Rouse 2015-4-9
编辑：Robert Rouse 2015-4-9
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"

请先登录，再进行评论。

Robert Rouse
Robert Rouse 2015-4-9
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
  2 个评论
Johan Löfberg
Johan Löfberg 2015-4-9
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
Johan Löfberg 2015-4-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Rubik's Cube 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by