(*Auxiliary function to compute all principal minors of a matrix*)
princminors[X_]:=(
minors={};
For[k=1,k<=Length[X],k++,
ktuples=Subsets[Range[Length[X]],{k}];
For[i=1,i<=Length[ktuples],i++,
minors=AppendTo[minors,Det[X[[ktuples[[i]],ktuples[[i]]]]]];
];
];
Return[minors];
);

(*Ex1*)
M=ArrayFlatten[{{{{1+x-y,x-1},{x-1,1+y}},0},{0,{{1+y,x,1},{x,x,y},{1,y,1}}}}];
M//MatrixForm
region=And@@Thread[princminors[M]>=0];(*Lists that all principal minors are nonnegative*)
regionPlot=RegionPlot[region,{x,-1,3},{y,-1,3},PlotPoints->80](*Plots region where all principal minors are nonnegative*)

(*Minors from first block*)
minors1=Thread[princminors[({{1+x-y, -1+x}, {-1+x, 1+y}})]>=0]
RegionPlot[#,{x,-5,5},{y,-5,5},PlotPoints->80,PlotStyle->Opacity[.1]]&/@minors1
plotminors1=Show[RegionPlot[#,{x,-5,5},{y,-5,5},PlotPoints->80]&/@minors1]

(*Apparently, it suffices to have the 3 x-(x^2)+x y-(y^2)≥0 *)
region1=RegionPlot[3 x-(x^2)+x y-(y^2)≥0,{x,-5,5},{y,-5,5},PlotPoints->80,PlotStyle->{Red,Opacity[.7]},BoundaryStyle->{Thick, Red}];
Show[plotminors1,region1]

(*Use quantifier elimination to confirm the above, i.e., to show that the red ellipse is indeed contained in the other regions:*)
CylindricalDecomposition[ForAll[{x,y},Implies[3 x-(x^2)+x y-(y^2)≥0,1+x-y≥0&&1+y≥0]],{x,y}]

(*For example, note that the converse inclusion does not hold:*)
CylindricalDecomposition[ForAll[{x,y},Implies[1+x-y≥0&&1+y≥0,3 x-(x^2)+x y-(y^2)≥0]],{x,y}]

(*Minors from second block*)
minors2=Thread[princminors[({{1+y, x, 1}, {x, x, y}, {1, y, 1}})]>=0]
RegionPlot[#,{x,-2,3},{y,-2,3},PlotPoints->80,PlotStyle->Opacity[.1]]&/@minors2
plotminors2=Show[RegionPlot[#,{x,-2,3},{y,-2,3},PlotPoints->80]&/@minors2]

(*Apparently, it suffices to have x≥0,y≥0,-(x^2)+3 x y-(y^2)-(y^3)≥0 : *)
region2=RegionPlot[x≥0&&y≥0&&-(x^2)+3 x y-(y^2)-(y^3)≥0,{x,-5,5},{y,-5,5},PlotPoints->80,PlotStyle->{Red,Opacity[.7]},BoundaryStyle->{Thick, Red}];
Show[plotminors2,region2]

(*Moreover, it looks like the region2 is contained inside region1!*)
(*Apparently, it suffices to have the 3 x-(x^2)+x y-(y^2)≥0 *)
region1=RegionPlot[3 x-(x^2)+x y-(y^2)≥0,{x,-5,5},{y,-5,5},PlotPoints->80,PlotStyle->{Red,Opacity[.2]},BoundaryStyle->{Thick, Red}];
Show[region1,region2]

(*So, apparently, the spectrahedron can be described as a semialgebraic set as follows:*)
regionReduced=x≥0&&y≥0&&-(x^2)+3 x y-(y^2)-(y^3)≥0 ;

(*Note the difference in complexity between the semialgebraic descriptions: *)

regionReduced
region

(*We now check algorithmically with Quantifier Elimination that the above semialgebraic sets are indeed identical*)
CylindricalDecomposition[ForAll[{x,y},Implies[regionReduced,region]],{x,y}]
CylindricalDecomposition[ForAll[{x,y},Implies[region,regionReduced]],{x,y}] (*trivial*)

(*With the above, we answer parts a) and b) of the problem.*)

(*An alternative approach would be to follow the proof of Proposition 1 in Lecture 21*)
coeffPoly=CoefficientList[(-1)^5Det[-t IdentityMatrix[5]- M],t];
And@@Thread[coeffPoly>=0]
Reduce[And@@Thread[coeffPoly>=0]]
Reduce[regionReduced](*coincides!*)
Reduce[region](*coincides!*)

(*Moving on to part c) about minimizing and maximizing a certain target function*)
Minimize[{x-y,regionReduced},{x,y}]
Maximize[{x-y,regionReduced},{x,y}]

cx=1;cy=-1;
lines=Plot[Table[c/cy-cx/cy x,{c,-3,3,1/3}],{x,-1,3},PlotStyle->Black];
min=Plot[c/cy-cx/cy x/.c->-5/27,{x,-1,3},PlotStyle->Red];
max=Plot[c/cy-cx/cy x/.c->1,{x,-1,3},PlotStyle->Red];
Show[regionPlot,lines,min,max]

(*First approach: Lagrange multipliers*)
Grad[x-y,{x,y}]==λ Grad[-(x^2)+3 x y-(y^2)-(y^3),{x,y}]

(*thus, 1/λ = (-2 x+3 y)==-(3 x-2 y-3 (y^2)), and hence *)

sol=Solve[(-2 x+3 y)==-(3 x-2 y-3 (y^2))&&-(x^2)+3 x y-(y^2)-(y^3)==0]

x-y/.sol

(*clearly, min and max are respectively the second and third entries above.*)

(*Second approach: solve for x and substitute it in to get a function of 1 variable*)
xsol=Solve[-(x^2)+3 x y-(y^2)-(y^3)==0,x]

(*Substitute above solutions in target function and find critical points:*)
ysol1=Solve[D[x-y/.xsol[[1]],y]==0];
ysol2=Solve[D[x-y/.xsol[[2]],y]==0];
{x-y,{x,y}}/.xsol[[1]]/.ysol1
{x-y,{x,y}}/.xsol[[2]]/.ysol2

(*Clearly, min is -5/27 and max is 1.*)