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}$).

• @Quantum_Oli One possibility is to try to rewrite the algorithm at the and of this paper: math.iupui.edu/~ccowen/Downloads/33NumRange.pdf into Mathematica. But I expect that something already has to exist in Mathematica and here someone may know. – Twi Apr 20 '17 at 8:56
• Somehow I have the feeling that this question needs to be explained more precisely. Otherwise you either won't get answers, or you will get answers which will not be useful to you. – Szabolcs Apr 20 '17 at 9:32
• @Szabolcs Maybe I should add that it would be nice to see where the boundary of the numerical range is located and what is its shape. The "random sampling" solution which you suggest below produces somewhat blurred figures and does not give a very good intuition. – Twi Apr 20 '17 at 10:18
• @Szabolcs , Twi, Perhaps ConvexHullMesh@ReIm@Table... – Michael E2 Apr 20 '17 at 11:21
• Wouldn't $x$ be a point in an n-dimensional space, and you are trying to plot a surface in that space in the new space "mapped" by A? The set of all $x$ defines the surface of the hypersphere of radius =1? Not sure how you plot it for $n>3$. – MikeY Apr 20 '17 at 13:37

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. • What matrix did you use for your example? – J. M. will be back soon Apr 20 '17 at 13:05
• I have included my choice in the answer. Thanks for pointing this out. – kalvotom Apr 20 '17 at 14:21

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] /;
Internal`EffectivePrecision[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]]}] 