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

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
8
$\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

Improve this answer
$\endgroup$
0
5
$\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

Improve this answer
$\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

Improve this answer
$\endgroup$
2
  • $\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$ Commented Apr 20, 2017 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$ Commented Apr 20, 2017 at 15:51

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.