# How to plot the field of values (numerical range) of a matrix

I found some difficulties in plotting the set $$W(\mathbb{A}):=\{\langle x,\mathbb{A}x \rangle \mid \|x\|=1\},$$ where $\mathbb{A}\in\mathbb{C}^{n,n}$ is a given complex matrix and $\langle\cdot,\cdot\rangle$ stands for the Euclidean inner product in $\mathbb{C}^{n}$. Is there any (perhaps straightforward) way to visualize $W(\mathbb{A})$ in Mathematica?

A comment: I am particularly interested in the localization of $W(\mathbb{A})$ in $\mathbb{C}$ with emphasis on its boundary and its shape. I am looking for a procedure to give a plot of $W(\mathbb{A})$ without any additional assumptions on the matrix $\mathbb{A}$. The size of $\mathbb{A}$, however, need not be too big, say $n\leq100$.

Remark: $W(\mathbb{A})$ is a convex subset of $\mathbb{C}$ located in the disc $|z|\leq\|\mathbb{A}\|$ (the spectral norm of $\mathbb{A}$).

The algorithm in the article mentioned above tries to find the boundary of the numerical range "from outside" and therefore may give more accurate results.

My naive Mathematica implementation is:

lsEigenvalue[H_] := Module[{emax, emaxspace, emin, eminspace, es},
(* Find the largest and smallest eigenvalue of a Hermitian matrix H together with the corresponding eigenspace (naive implementation) *)

es = Sort[Eigenvalues[H], Re[#1] < Re[#2] &];
emin = First[es];
emax = Last[es];
eminspace = Orthogonalize[NullSpace[H - emin*IdentityMatrix[Length[H]]]];
emaxspace = Orthogonalize[NullSpace[H - emax*IdentityMatrix[Length[H]]]];
Return[{{emin, eminspace}, {emax, emaxspace}}]
]

plotNR[A_, n_] := Module[{t = 0., td = 2 π/n, Ht, Kt, points = {}, segments = {}, emax, emaxspace, emin, eminspace, vp, vm, Q, R},
PrintTemporary[ProgressIndicator[Dynamic[t], {0, 2 π}]];
(* data for numeric range plot *)
While[t < 2 π,
Ht = (Exp[-I t]*A + Exp[I t]*ConjugateTranspose[A])/2;
{emax, emaxspace} = Last[lsEigenvalue[Ht]];

Which[(* One dimensional eigenspace *)
Length[emaxspace] == 1,

vp = First[emaxspace];
AppendTo[points, Conjugate[vp].A.vp],

(* Two or greater dimension -- almost does not happen? *)

Length[emaxspace] > 1,

Kt = (Exp[-I t]*A - Exp[I t]*ConjugateTranspose[A])/(2 I);
Q = Transpose[emaxspace];
R = ConjugateTranspose[Q].Kt.Q;

{{emin, eminspace}, {emax, emaxspace}} = lsEigenvalue[R];

vp = Q.First[emaxspace];
vm = Q.First[eminspace];

AppendTo[segments, {Conjugate[vm].A.vm, Conjugate[vp].A.vp}],

(* Fail *)
True,
Print["Error"]
];

t = t + td;
];
Return[{DeleteDuplicates[points], DeleteDuplicates[segments]}]
]


An example result for a $2\times 2$ matrix

$$\begin{pmatrix} -1 & i \\ 2 & 3i \end{pmatrix}.$$

Red dots are computed by randomly sampling vectors and computing the quadratic form, blue dots are computed by the algorithm above. The following is a somewhat naïve and inefficient Monte Carlo sampler that samples with higher probability far from a given point, and approximately uniformly according to the argument of complex numbers. This it forces the samples closer to the boundary of $W$.

We start by choosing $n$ and $A$:

n = 15;
a = RandomComplex[{-1 - I, 1 + I}, {n, n}];


If we sample uniformly from $\{x \mid \|x\|=1\}$, the result looks "blurry" because all the points are near the centre of $W$. But this at least allows us to estimate where its centre is:

xs = Normalize /@ RandomComplex[{-1 - I, 1 + I}, {100000, n}];
points = Table[Conjugate[x].a.x, {x, xs}];
center = Mean[points]
(* 0.213328 - 0.141669 I *) Choose a random starting point for the Monte Carlo sampler:

x = Normalize @ RandomComplex[{-1 - I, 1 + I}, n];
y = Conjugate[x].a.x;


The vector arghist will keep a histogram of the complex arguments. We use it to force the sampler away from "directions" it has already explored. This ensures that it will go around the boundary.

arg2ind[z_] := 1 + Floor[(π + Arg[z]) k/(2 π)]

k = 100;
arghist = ConstantArray[0., k];


The parameters may need to be tuned for each problem:

beta1 = 200; (* lager: close to the boundary *)
beta2 = 0.5; (* forcing away from already explored portions of the boundary *)
stepsize = 0.05;
result = Union @ First @ Last @ Reap @ Do[
xp =
Normalize[x + RandomComplex[stepsize {-1 - I, 1 + I}, n]];
yp = Conjugate[xp].a.xp;
If[
RandomReal[] <
Exp[beta1 (Norm[yp - center] - Norm[y - center]) +
beta2 (arghist[[arg2ind[y - center]]] -
arghist[[arg2ind[yp - center]]])],
{x, y} = {xp, yp}
];
arghist[[arg2ind[y]]] += 1;
Sow[y],
{100000}
];

Graphics[{AbsolutePointSize, Point @ ReIm[result], Red,
PointSize[Large], Point @ ReIm[center]}, Frame -> True, Axes -> True,
AxesOrigin -> {0, 0}] Since you say $W$ is convex, we can take the convex hull of these points per Michael's comment:

ConvexHullMesh @ ReIm[result] • It seems that you forgot the complex conjugation in the inner product, $\langle x,y\rangle =\sum_{i}\overline{x_{i}}y_{i}$. So instead of x.A.x, you should have Conjugate[x].A.x, etc. – Twi Apr 20 '17 at 12:26
• @Twi Fixed. The Monte Carlo method could be improved to sample uniformly (and more quickly) from $W$, but since someone has already posted an implementation of the algorithm from the paper, there is no point to it. I only fixed the mistake you pointed out without improving on the efficiency for oblong regions. – Szabolcs Apr 20 '17 at 15:51

Here is some compact code for evaluating the complex curve corresponding to the field of values of a matrix, using the method of Johnson:

fval[mat_?SquareMatrixQ, t_?NumericQ] /;
InternalEffectivePrecision[mat] < ∞ || InexactNumberQ[t] := Module[{tm, v},
tm = (# + ConjugateTranspose[#])/2 &[mat Exp[I t]];
v = Quiet[Check[First[Eigenvectors[tm, 1, Method -> {"Arnoldi",
"Criteria" -> "RealPart"}]],
MaximalBy[Transpose[Eigensystem[tm]], First][[1, -1]],
Eigenvectors::arall]];
(Conjugate[v].mat.v)/(Conjugate[v].v)]


One can then use ParametricPlot[] for the visualization:

(* Grcar matrix, https://arxiv.org/abs/1203.2390 *)
grcar[r : _Integer?Positive : 3, n_Integer?Positive] :=
SparseArray[{{j_, k_} /; j == k + 1 :> -1, {j_, k_} /; 0 <= k - j <= r :> 1}, {n, n}]

mat = grcar; eigs = Eigenvalues[N[mat]];

ParametricPlot[ReIm[fval[mat, t]], {t, 0, 2 π},
Epilog -> {AbsolutePointSize, ColorData[97, 2], Point[ReIm[eigs]]}]
` 