# RegionPlot is slow to plot a simple region

I start out with 5x5 matrix gCorrMatrix[g] depending on a parameter $$g$$ and an unknown $$W_2$$

gCorrMatrix[g_] := {{1, 0, Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g)},
{0, Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g), 0},
{Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g), 0,
-((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2)},
{0, -((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/
g^2), 0}, {-((-1 + Subscript[W, 2])/g), 0,
-((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2), 0,
-((-1 - 2*g + Subscript[W, 2] + 4*g*Subscript[W, 2] -
g^2*Subscript[W, 2]^2)/g^3)}}


I define a function of this matrix that returns an inequality

gRegion[g_] := Reduce[Eigenvalues[gCorrMatrix[g]] >= 0, Subscript[W, 2]]


that gives a range of values of $$W_2$$ for each $$g$$ where gCorrMatrix[g] is positive-definite. Evaluating gives

In[157]:= Timing[gRegion[.3]]

During evaluation of In[157]:= Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

Out[157]= {0.083676, 0.703257 <= Subscript[W, 2] <= 0.714234}


Clearly this is computationally expensive to evaluate over and over again. Mathematica just hangs when I try to plot

RegionPlot[gRegion[g], {g, -1/8, .5}, {Subscript[W, 2], .5, 2}]


This surprised me because gRegion[g] is very simple (just a line segment and independent of $$g$$). Is there a way I can tell RegionPlot to just plot each line segment of $$W_2$$ sequentially? This really shouldn't take more than a few seconds. I thought I could change this with the Method option but I can't see anything in the documentation aside from Method->"Automatic".

Redefine your matrix to explicitly depend on the values of $$g$$ and $$W_2$$ (and, while you're at it, lose the Subscript expressions that make everything a bit more complicated):

ClearAll[gCorrMatrix]
gCorrMatrix[g_, w_] :=
{
{1, 0,w, 0, -((-1 + w)/g)},
{0, w, 0, -((-1 + w)/g), 0},
{w, 0, -((-1 + w)/g), 0, -((1 - w - 2*g*w)/g^2)},
{0, -((-1 + w)/g), 0, -((1 - w - 2*g*w)/g^2), 0},
{-((-1 + w)/g), 0, -((1 - w - 2*g*w)/g^2), 0, -((-1 - 2*g + w + 4*g*w - g^2*w^2)/g^3)}
}


Then create a predicate function that checks whether the matrix above is positive definite when given numerical values of the two parameters $$g$$ and $$W_2$$:

ClearAll[gPosQ]
gPosQ[g_?NumericQ, w_?NumericQ] :=
PositiveDefiniteMatrixQ[gCorrMatrix[g, w]]


You can then use that function as a predicate in RegionPlot:

RegionPlot[
gPosQ[g, w],
{g, -1/8, 0.5}, {w, 0.5, 2},
PlotPoints -> 60, MaxRecursion -> 8
]


Adjust the values of PlotPoints and MaxRecursion to attain your preferred plotting goals.

• Nice trick using ?_NumericQ to prevent any symbolic manipulation. One comment—the region in the neighborhood of $g\approx 1/2$ has vanishing measure and hence isn't shown by RegionPlot. This isn't a defect—it's in the documentation—but it is a limitation of this type of analysis. As I increase the size of the correlation matrix (computed before) and improve the bound, RegionPlot sees less and less of the graph. It's still very good for visualizing the uncertainty in $W_2$, but extracting the solution function $W_2(g)$ needs a different approach. Commented Mar 18, 2022 at 23:13

Edit

Clear[gCorrMatrix, sol, lowsol, upsol];
gCorrMatrix[
g_] = {{1, 0, Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g)}, {0,
Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g),
0}, {Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/
g^2)}, {0, -((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2),
0}, {-((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2),
0, -((-1 - 2*g + Subscript[W, 2] + 4*g*Subscript[W, 2] -
g^2*Subscript[W, 2]^2)/g^3)}} /. Subscript[W, 2] -> w2;
sol[g_] :=
Reduce[Eigenvalues[gCorrMatrix[Rationalize[g, 0]]] >= 0, w2];
lowsol[g_] := sol[g][[1]];
upsol[g_] := sol[g][[-1]];
Plot[{lowsol[g], upsol[g]}, {g, 0, 5}, Filling -> {1 -> {{2}, Red}},
AxesOrigin -> {0, 0}]
Plot[{lowsol[g], upsol[g]}, {g, -.1, -0.005},
Filling -> {1 -> {{2}, Cyan}}, AxesOrigin -> {0, 0}]


• When g is fixed,since this is one dimensional plot,we can use NumberLinePlot instead.
gCorrMatrix[
g_] = {{1, 0, Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g)}, {0,
Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g),
0}, {Subscript[W, 2], 0, -((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/
g^2)}, {0, -((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2),
0}, {-((-1 + Subscript[W, 2])/g),
0, -((1 - Subscript[W, 2] - 2*g*Subscript[W, 2])/g^2),
0, -((-1 - 2*g + Subscript[W, 2] + 4*g*Subscript[W, 2] -
g^2*Subscript[W, 2]^2)/g^3)}} /. Subscript[W, 2] -> w2;
plot[g_] :=
NumberLinePlot[
Reduce[Eigenvalues[gCorrMatrix[Rationalize[g, 0]]] >= 0, w2], w2]
plot[.3]
plot[-.001]


• To draw the region
region = ImplicitRegion[
Eigenvalues[gCorrMatrix[g]] >= 0 // Thread, {g, w2}];
RegionPlot[RegionMember[region, {g, w2}], {g, -.2, .5}, {w2, .5, 8},
MeshFunctions -> {#1 &}, Mesh -> 5, MeshStyle -> Red]
RegionPlot[
Eigenvalues[gCorrMatrix[g]] >= 0], {g, -.2, .5}, {w2, .5, 8}]


solve[g_] :=
Reduce[Eigenvalues[gCorrMatrix[Rationalize[g, 0]]] >= 0, w2];
reg[g_] := RegionProduct[Point[{{g}}], ImplicitRegion[solve[g], {w2}]];

right = RegionUnion[Table[reg[g], {g, .01, .5, .005}]] // Region
left = RegionUnion[Table[reg[g], {g, -.1, -.01, .005}]] // Region
Show[left, right, PlotRange -> All, AspectRatio -> 1, Axes -> True]