4
$\begingroup$

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

$\endgroup$
7
$\begingroup$

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.

enter image description here

$\endgroup$
4
$\begingroup$

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 *)

Mathematica graphics

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[1], Point @ ReIm[result], Red, 
  PointSize[Large], Point @ ReIm[center]}, Frame -> True, Axes -> True, 
 AxesOrigin -> {0, 0}]

Mathematica graphics

Since you say $W$ is convex, we can take the convex hull of these points per Michael's comment:

ConvexHullMesh @ ReIm[result]

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Twi Apr 20 '17 at 12:26
  • 1
    $\begingroup$ @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. $\endgroup$ – Szabolcs Apr 20 '17 at 15:51
3
$\begingroup$

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[32]; eigs = Eigenvalues[N[mat]];

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

field of values of a Grcar matrix

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.