text:
\begin{document}
\date{~}
\thispagestyle{empty}
\maketitle

\thispagestyle{empty}

\vspace*{-16mm}
\begin{center} { \bf Abstract  } \\[2mm]
\begin{minipage}{140mm}
This paper  describes and develops a fast and accurate  path following algorithm that computes the field of values boundary curve $\partial F(A)$ for every conceivable complex or real square matrix $A$. It relies on the matrix flow decomposition algorithm that finds a proper block-diagonal flow representation for the associated hermitean matrix flow ${\cal F}_A(t) = \cos(t) H + \sin(t) K$ under unitary similarity if that is possible. Here ${\cal F}_A(t)$  is the 1-parameter-varying linear combination of the real and skew part matrices $H = (A+A^*)/2$ and $K = (A-A^*)/(2i)$ of $A$. For indecomposable matrix flows, ${\cal F}_A(t)$ has just one block and the ZNN based field of values algorithm works with ${\cal F}_A(t)$ directly. For decomposing  flows ${\cal F}_A(t)$, the algorithm decomposes the given  matrix $A$ unitarily  into block-diagonal form $U^*AU =  \text { diag} (A_j)$  with $j > 1$ diagonal blocks $A_j$ whose individual  sizes add up to the  size of  $A$. It  then computes the field of values boundaries separately  for each diagonal block $A_j$ using  the  path following ZNN eigenvalue method.  The convex hull of all sub-fields of values boundary points $\partial F(A_j)$ finally determines the field of values boundary curve correctly for  decomposing matrices $A$. The algorithm removes standard restrictions for path following FoV methods that generally cannot deal with decomposing matrices $A$ due to possible eigencurve crossings of ${\cal F}_A(t)$. Tests and numerical comparisons are included. Our ZNN based method is coded for sequential and parallel computations and both versions run very accurately and   fast when compared with Johnson's  Francis QR eigenvalue  and Bendixson rectangle based method that computes complete eigenanalyses of ${\cal F}_A(t_k)$ for every chosen $t_k  \in {[} 0,2\pi{]}$ more slowly.
\end{minipage}\\[-1mm]
\end{center}
\thispagestyle{empty}

\noindent{\bf Keywords :}  field of values, matrix flow, single parameter varying matrices, time-varying matrix flow,\\
\hspace*{17.5mm}  decomposable matrix, numerical algorithm,  block-diagonal matrix, unitary similarity\\[-3mm]

\noindent{\bf AMS Classifications :} 15A99, 15B99,  65F99, 65F45, 15A21\\[-7mm]

\pagestyle{myheadings}
\thispagestyle{plain}
\markboth{Frank Uhlig}{Decomposing and General Matrix FoVs }

\section{Field of Values Computations for  Real or Complex Square Matrices; Preliminaries}\vspace*{-1mm}

This paper develops an algorithm to construct the {\em Field of Values} (FoV) of a  real or complex square matrix $A$ via any chosen path following method, even when $A$ is unitarily block-diagonalizable into $j > 1$ diagonal blocks $A_j$. Constructing the FoV  of properly  decomposable matrices $A$ was previously only possible by using complete matrix eigenvalue finders such as Francis' QR algorithm. Here we solve this problem in three steps by \\[0.5mm]
  (a) checking for unitary diagonalizability of a given matrix $A$ via ${\cal F}_A(t)$, see \cite {FUDecomp},  and unitarily block-diagonalizing  decomposable matrices into $j >1$ diagonal blocks, then\\[0.5mm]
   (b) using a ZNN based path following method to compute discrete FoV boundary eigendata  for each  diagonal block $A_j$ separately, and finally by \\[0.5mm]
   (c) implementing a convex hull algorithm to plot the FoV boundary curve of $A$ itself.\\[2mm]
The  {\em field of values} $F(A)$, also known as the {\em numerical range}  of a matrix $A$ is defined as
$$  F(A)\ = \ \{ x^*Ax \mid x \in \CC^n, \|x\|_2 = 1\}  \subset \CC \ .$$
Its genesis, original purpose  or inventors are unknown. In 1918/1919 Hausdorff and Toeplitz showed independently that $F(A)$ is  convex for any square matrix $A$. Some well known properties of the field of values $F(A)$ are:\\
The field of values $F(A)$ is compact for any $A$; the field of values  of a 2 by 2 matrix is an ellipse; $F(A)$ contains all eigenvalues of $A$; $F(c A) = cF(A)$ for all scalars $c \in \CC$; $F(U^*AU) = F(A)$ for any unitary matrix $U$; for normal matrices $A$, $F(A)$  is the convex hull of $A$'s eigenvalues;
for block-diagonal matrices $A=\bp A_1 & O\\O & A_2\ep$, $F(A) = \text{ convhull }  ( F(A_1), F(A_2))$.\\[1mm]
  The matrix field of values is treated in many textbooks, such as in \cite{HJ94}. The maximal and minimal distant points of the boundary $\partial F(A)$ of the field of values from the origin $0 \in \CC$ play an important role in matrix norm estimations and in control theory, respectively. The former is called the {\em numerical radius} of the matrix and the latter is the {\em Crawford number} of $A$. Generalized numerical ranges are important notions in the study of quantum computations and elsewhere.\\[2mm]
 Building on Bendixson rectangles \cite{Be} from 1902 that contain the field of values of a matrix, Johnson  \cite{J78} in 1978 established an eigendata based method to find  discrete boundary curve points on $\partial F(A)$ of the field of values  for a matrix $A$ via repeated hermitean eigendata computations. Johnson's method involves the  matrices\\[-2mm]
  $${H=  (A+A^*)/2 = H^*} \ \ \text{ and } \ \
 {K=  (A - A^*)/(2i) =  K^*}  \in \CC_{n,n} $$
   with $A = H + iK$ and the  associated 1-parameter-varying   hermitean matrix flow\\
 \begin{equation}
 \label{Aoft} {\cal F}_A(t) = \cos(t)  H +   \sin(t)   K  = ({\cal F}_A(t))^*  \ \ \ \text{ for  angles } \ \ \ 0 \leq t \leq 2\pi \ .
\end{equation}
More specifically, the normalized  eigenvectors $x(t)$ and $y(t)$ for the largest and smallest eigenvalues of each ${\cal F}_A(t)$ determine two $\partial F(A)$ { boundary points} via the quadratic form evaluations  $x(t)^*Ax(t) \text{ and } y(t)^*Ay(t) \in \CC$  and two $\partial F(A)$ tangents. One standard approach to approximate  the boundary curve $\partial F(A)$ of a given $n$ by $n$  matrix $A$ with $n$ up to around $n=10,000$  is to use the Francis implicit multi-shift QR algorithm  for many discrete parameter values $0 \leq t_i \leq 2\pi$ and find the extreme eigendata of each ${\cal F}_A(t_i)$. Then evaluate  the quadratic forms $x(t_i)^*Ax(t_i) \text{ and } y(t_i)^*Ay(t_i)$ for the respective 'extreme' eigenvectors $x(t_i)$ and $y(t_i)$ of ${\cal F}_A(t_i)$ in order to  find accurate  interpolation points on the FoV boundary curve   $\partial F(A)$. For huge dimensions $n \gg 10,000$ other ways can be used to find the extreme eigenvalues and associated eigenvectors iteratively such as Matlab's {\tt eigs}  m-file.\\[1mm]
For relatively small dimensioned matrices $A$, the global hermitean QR eigenvalue approach computes complete eigendata  sets for  each ${\cal F}_A(t_i)$ reliably. Yet  we need know all eigenvectors of ${\cal F}_A(t_i)$; the  extreme ones suffice to  plot $\partial F(A)$. Francis' implicit multi-shift QR  and other complete matrix eigensolvers offer a foolproof way to solve the FoV problem for all matrices $A$ of modest size, decomposable or not. For large dimensions $n$, with $n$ near $10,000$, global eigen methods become rather  expensive when compared with  simpler and much faster path following methods. But until now, path following methods alone could not handle the FoV problem for unitarily decomposable matrices,  see e.g, \cite[Algorithm 6.1 and p. 1733 - 1743]{LM18} for a detailed analysis. How to determine whether a given complex matrix $A$ can be block-decomposed by a fixed unitary matrix similarity  has recently been solved constructively in \cite{FUDecomp}. Attempts at unitary block-decompositions of static entry matrices and general matrix flows go back at least 90 years to the beginnings of quantum physics and quantum chemistry when von Neumann and Wigner \cite{NW29}  established two (in-)decomposability criteria for 1-parameter hermitean matrix flows from their eigencurve behavior. Namely, if for some $0 \leq t < 2\pi$ two eigencurves for ${\cal F}_A(t)$  cross each other, then their respective eigendata is associated with  separate diagonal blocks of the hermitean matrix flow  ${\cal F}_A(t)$ -- and thus by \cite{FUDecomp} also of $A$. If on the other hand two eigencurves seem to be attracted to each other and almost touch but don't actually cross but veer off like the two branches of a hyperbola do, then the respective eigencurves belong to the same unitarily indecomposable block of $A$, see \cite{NW29} and \cite{FUDecomp} for further details.\\[2mm]
Some path following methods for 1-parameter or time-varying matrix problems like ours track the solution of a derived differential equation by using numerical integrators. To do so for the matrix FoV problem, Loisel and Maxwell \cite[sections 4 and 6]{LM18}  for example,  differentiate the Bendixson/Johnson  eigenvalue equation
$$ {\cal F}_A(t) u(t) = \lambda(t)u(t) \ $$
with respect to $t$.
When additionally  requiring  that the eigenvectors $u(t)$ be normalized, they obtain a 1-parameter matrix and vector differential equation in $\CC^{n+1}$ for the matrix FoV problem. Their  ODE path following method uses the eigendata of ${\cal F}_A(0)$ as the initial value and then proceeds with the Dormand--Prince RK5(4)7M numerical integrator \cite{DP80} of fifth order accuracy $O(h^5)$. This  adaptive method can be implemented inside Matlab as {\tt ode45}. It is quite versatile, accurate  and well suited  for stiff differential equations. The computed $\partial F(A)$ points are then smoothed through  fourth degree Hermite interpolation in \cite{LM18}. Unfortunately the codes and error threshold settings etc used by Loisel and Maxwell are no longer available. For unitarily indecomposable matrices $A$, the path following results  in \cite{LM18} are as accurate as those from the Francis QR matrix eigenvalue algorithm but much faster than a complete QR {\tt eig} implementation.  The accuracy and speed of computing  indecomposable matrix FoVs was subsequently bested  again in \cite{FUZNNFoV}  by using {\em Zhang Neural Network}  (ZNN) based eigen-computations as the  path following method. Zhang Neural Networks   are specially designed  for time-varying matrix problems. They proceed in a totally different  way than classical integrators. They start from a time-varying matrix model and its global error function $E(t)$ and demand exponential decay for $E$ in time $t$. The postulated error differential equation $\dot E(t) = -\eta E(t)$ with $\eta > 0$ is then rewritten to become derivative free. And the resulting equation  is solved at each discrete time step $t_k$ for the future time step $t_{k+1}$ by one linear equations solve and using a convergent look-ahead finite difference formula that relates the solution $x(t)$ in the future at $t = t_{k+1}$ to earlier known systems data. The theoretical and practical construction  of discretized ZNN methods for various time-varying matrix problems is the subject of \cite{FUZNNsurvey}.\\[2mm]
Today  general path following methods  have become suitable for solving general time-varying matrix problems  due to our ability to decipher dense matrices and matrix flows that are unitarily block-diagonalizable,  both theoretically and computation-wise, as detailed in \cite{FUDecomp}. The unitary block-diagonalizing algorithm for decomposable static matrices $A$ of \cite{FUDecomp} is easily described. Form the Bendixson/Johnson hermitean matrix flow ${\cal F}_A(t)$. Diagonalize one flow matrix ${\cal F}_A(t_a)$ unitarily with $V$ by using Matlab's {\tt eig} function for example. Then check whether any other flow matrix ${\cal F}_A(t_b) \neq {\cal F}_A(t_a)$ exhibits a logical zero-nonzero block-diagonal structure with more than 1 diagonal block. If  so, then $A$ itself will be block-diagonalized unitarily by the same unitary similarity $V^* .. V$ -- after we have re-sorted the columns of $V$  into invariant eigenvector groups for ${\cal F}_A(t_b)$. This decomposability check  is elementary, accurate and fast. It works universally for any real or complex square matrix $A$. To check on the unitary block-diagonalizability of a dense matrix $A \in \CC_{n,n}$  with $n = 250$ for example  takes around 0.05 seconds for both decomposable and indecomposable matrices $A$. The new matrix decomposition algorithm \cite{FUDecomp} does not depend on the eigen or Jordan structure of  $A$, nor on any other conceivable property or defect of $A$. It derives solely from the properties of the associated hermitean matrix flow ${\cal F}_A(t)$.  And irrespective of  whether $A$ is found to be  unitarily block-decomposable or not by \cite{FUDecomp}, the field of values boundary curve $\partial F(A)$ data can now be computed by using any ODE path follower or the fast ZNN method \cite{FUZNNFoV} for  $A$ as a whole or for each  diagonal block  $A_j$ of a uniform block-diagonalization of $A$ separately  if $A$ can be properly decomposed. To plot the overall FoV boundary curve  in the decomposable case, we use Matlab's {\tt convhull} algorithm on the totality of all computed diagonal block FoV boundary points. This method is  new and much quicker than using a global eigensolver for either case, i.e.,  independently of whether $A$ decomposes or not.\\[1mm]
This paper combines our recent knowledge of how to  block-diagonalize both static matrices and matrix flows by unitary similarity  \cite{FUDecomp} (which had been unrealized for almost a century) with the recent introduction  of a path following ODE method \cite{LM18} to compute the field of values of indecomposable matrices.  We develop  and test one algorithm that can now find the field of values boundary for decomposable --  unavailable before -- and indecomposable matrices quickly using a path-follower that is powered by  Zhang Neural Networks  for the single-parameter-varying matrix FoV problem with high accuracy and speed.\\[2mm]
Section 2  illustrates the new computational landscape for unitarily invariant matrix problems and shows that decomposable matrices do not lend themselves well for unitarily invariant matrix and  'divide and conquer'  computations.  Section 3  explains the complexity gains of path following methods for decomposing matrices and it describes some aspects of the ZNN method \cite{FUZNNsurvey} very briefly, followed by  numerical test results for our combined FoV algorithm.  Section 4  looks back in history and  forward.

\section{Path Following Methods for  Field of Values Computations; a Warning Illustration}

In this section we describe the ill results of trying speedy path following methods to compute the field of values of a  dense matrix $A$ that is unitarily block-diagonalizable.\\
 For this purpose we construct a dense block-diagonalizable 15 by 15 matrix $A$ beginning with  the Matlab command  {\tt   B = blkdiag(1i*randn(10),-randn(5)-(3-2i)*eye(5));} . $B$ is a 15 by 15 non-normal block-diagonal matrix. Then we  create a dense 15 by 15 unitary random entry matrix $Q$ via {\tt [Q,R]= qr(randn(15));} and  form the test matrix  {\tt A = Q'*A*Q;} . Clearly the resulting matrix A is non-normal and its unitary block-diagonalizability is  hidden from view. Next one might choose any path following method to plot a discrete sample of FoV boundary points from a partial eigenanalysis of ${\cal F}_A(t_i)$ with $ 0 \leq t_i \leq 2\pi$ of $A$ and using quadratic form evaluations with $A$ as described earlier. Which path following method we choose, be it an initial value ordinary differential equations solver such as recently done in \cite{LM18}, or a ZNN based one, see \cite{FUZNNFoV} e.g., or any other method is irrelevant here. Look at the perturbing FoV  boundary curve drawing for our dense test matrix $A$  in Figure 1. \\[-5.2mm]
\begin{center}
\includegraphics[width=95mm]{F1Blk105folly.png} \\[-1mm]
Figure  1: Path following FoV output of a unitarily block-diagonalizable matrix $A$, showing the both extreme eigenvalue induced $\partial F(A)$ points for $0\leq  t \leq \pi$ of ${\cal F}_A(t)$, moving counterclockwise  separately and resulting in two partial and disjoint curves
\end{center}
\vspace*{-0.5mm}
Disjoint FoV mis-presentations such as depicted in Figure 1 or similar grotesque incomplete FoV representations are standard with  path following methods when applied unwittingly to  unitarily block-decomposable matrices. By design, these methods follow the initially chosen extreme ${\cal F}_A(t)$ eigenvalue curve faithfully and reflect the problems with computing the FoV boundary naively for decomposable matrices.\\[1mm]
Next we plot the two individual block FoV boundary curves completely  for our $A$  that are hinted at in Figure 1 after  computing the underlying unitary block decomposition of $A$ via the block-diagonalizing algorithm of \cite{FUDecomp}.\\[-6mm]
\begin{center}
\includegraphics[width=95mm]{F2Blk105two.png} \\[-1mm]
Figure  2: Path following FoV output of our dense test matrix $A$, plotting the FoV boundary of each diagonal block  for $0 \leq t < 2 \pi$ counterclockwise
\end{center}
\begin{center}
\includegraphics[width=96mm]{F3Blk105chull.png} \\[-1mm]
Figure  3: Displaying the FoV boundary curve of $A$ completely, found via a path following method from the individual FoV point data of $A$'s diagonal blocks in Figure 2 and then  extracting  their convex hull
\end{center}

\vspace*{-2mm}
For almost a century mathematicians and physicists have studied the eigenvalue curves of matrix flows $A(t)$. Hermitean matrix flows were studied by Hund \cite{FH1927},  by von Neumann and Wigner \cite{NW29}, by Johnson \cite{J78}, by Dieci et al \cite{DE99,DF01,DP03,DPP13} and others, see \cite[section 5]{FUDecomp}.
More specifically,  coalescing eigencurves and eigencurve crossing points were studied and computed in \cite{DPP13} using a Newton method based local optimization algorithm.\\
Next we display the eigencurves of our chosen block-decomposable test matrix $A$ and then use the block-size assignment algorithm of \cite{FUCoalesc} to find $A$'s hidden block-structure to locate the lead changes in the extreme eigencurves of ${\cal F}_A(t)$.\\[-8mm]
\begin{center}
\includegraphics[width=95mm]{F4Blk105ecurva.png} \\[-1mm]
Figure  4: Eigencurves for $A$ and $0 \leq t \leq 2\pi$
 \end{center}

\vspace*{-1.5mm}
In Figure 4, note that one eigencurve group  crosses other eigencurves freely, while another eigencurve group veers off hyperbolically from some of the eigencurves. This coalescing and avoidance behavior was first discovered and interpreted by von Neuman  and Wigner  in 1929 for hermitean matrix flows, see \cite{NW29}. Eigencurves that veer off from each other and do not cross are associated with the same diagonal block for a unitarily similar block representation of the flow and those that cross are associated with different diagonal blocks of a  unitarily similar block representation. The simple geometric coalescing eigencurve algorithm of \cite{FUCoalesc} shows that the visually dense matrix $A$ is unitarily decomposable into two diagonal blocks of dimensions 10 and 5 as shown  in Figure 5 where 2 colors suffice to depict the respective eigencurve groups of our 15 by 15 dense test matrix $A$.\\[-4.5mm]
\begin{center}
\includegraphics[width=92mm]{F5Blk105ecurvblka.png} \\[-1.5mm]
Figure  5: Blocked eigencurves for $A$ and $0 \leq t \leq 2\pi$
 \end{center}
  It is easy to assess the 'lead changes' among the 'maximal' eigenvalue curves of ${\cal F}_A(t)$ for our test matrix $A$ in Figures 4 and 5.  The  first lead change  occurs just a little before $\pi/2$ or at around $ 80 ^o$ when rotating  the Bendixson rectangle or the positive real axis in Figure 3 counterclockwise around the origin. The  lead among the maximal eigenvalue curves reverts back in Figure 5 to the previous eigencurve leader at an angle nearly halfway between  $\pi$ and $3/2 \ \pi$ or at around $ 180 ^o + 40 ^o = 220 ^o$ degrees. Both of these  values corroborate well with the two straight line $\partial F(A)$ parts of Figure 3 if we realize that the eigencurve crossing  angles above match the angles of the perpendiculars to the two straight line segments of the FoV boundary curve for $A$ in its upper region and  lower left region  in Figure 3.\\[-5mm]

 \section{The ZNN Path Following Method for General Indecomposable and Decomposable Matrix Fields of Values}

 \subsection{Complexity gains for path following methods}

We first assess the theoretical CPU time gains from using  path following 'divide and conquer' matrix algorithms of complexity $O(n^3)$ when applied to properly decomposable $n$ by $n$ matrices $A$ and matrix flows $A(t)$.\\
Once we  know of a unitary block-diagonalization into smaller blocks of sizes $m_i < n$ with $\sum_{i} m_i = n$  for an  $n$ by $n$ matrix $A$ or matrix flow $A(t)$,  we can speed up  any $O(n^3)$ unitarily invariant computing process to the sum of several $O(m_i^3)$  processes with $m_i < n$ for the diagonal block dimensions $m_i$ of a unitarily similar block-diagonal  representation of $A$ or $A(t)$.\\
If $A$ or $A(t)$ cannot be unitarily block reduced, then  such numerical computations  will still require the full dense matrix $O(n^3)$ effort.\\
At the other extreme, if $A$ can be diagonalized by unitary similarity into 1 by 1 diagonal block form, i.e., if $A$ is normal, then any  unitarily invariant  $O(n^3)$ process  on $A$ takes only $O(n)$ effort after -- for example -- an $O(n^3)$ QR based eigenvalue diagonalization of $A$ was obtained, to compute the desired result.  What happens in between these extremes? \\[1mm]
For decomposable $n$ by $n$ matrices and any $O(n^3)$ numerical process we look at the maximal block size $1 \leq m \leq n$ of a given matrix $A$ once its block sizes have been computed via {\tt  decHKflowFoV.m} in \cite{FUDecompMatrFoV}.\\
If $A$ has 1 by 1 diagonal blocks except for its one $m$-dimensional block, then the full $O(n^3)$ process becomes one of smaller complexity $O(m^3) = (m/n)^3 O(n^3)$ or $O(m^3) = \alpha^3 O(n^3)$ for $\alpha = m/n$ with $0< \alpha < 1$ as  the remaining one-dimensional blocks require almost no work. This is the best possible scenario when the largest block-diagonal dimension that can be unitarily achieved for $A_{n,n}$  is $m < n$. \\
The worst operations savings scenario happens when the block-diagonalization of $A$ has {\tt floor}($n/m$) blocks of the same maximal size $m$ and  only one additional smaller block, needed  so that the individual block dimensions $m_i$ add up to $n$. Until $m = n/2$ there can be at most {\tt floor}($n/m$) maximal size $m$ blocks. Thus for $m \leq n/2$ the total operations cost of dealing with {\tt floor}($n/m$) maximal size $m$ diagonal blocks is around
 $$\text{{\tt floor}}(n/m) O(m^3) \leq n/m \cdot O(m^3) = n/m \cdot (m/n)^3 O(n^3) = (m/n)^2 O(n^3) = \al^2 O(n^3)$$
  operations for $\al = m/n \leq 1/2$ or $m \leq n/2$. If $m > n/2$ is the size of the largest diagonal block for $A$, then in the 'worst case' there would be one additional indecomposable block of size $n-m < m$ resulting in $O(m^3) + O((n-m)^3)$ necessary operations. Below is a graph of the extreme bounds for the operations cost for a matrix that decomposes unitarily into diagonal blocks of maximal size $m = \alpha \cdot n$. This graph describes the CPU time gain situation in terms of $0 < \alpha = m/n \leq 1$.\\[-5mm]
\begin{center}
\includegraphics[width=100mm]{F6CPUtimegaina.png} \\[-1.5mm]
Figure  6
\end{center}

\vspace*{-2mm}
The largest gap in possible FoV operations counts for decomposable matrices with maximal block dimension $m$ and unitarily invariant computations  occurs at $\alpha = m/n = 0.5$ or when $m = n/2$, see  Figure 6. The theoretical CPU time savings at $ m = n/2$  lie between  87.5 \%  and 75 \% compared to the full $O(n^3)$ effort. Sizable runtime savings of at least 50\% occur  for all $m < 0.7 n$ as can be read off  Figure 6.\\
 Recall the 15 by 15 block-diagonalizable test matrix $A$ of Section 2 whose largest diagonal block had size 10 by 10. Here $\al = m/n = 10/15 = 0.667$. Figure 6 indicates that for any $n$ and $\al = 0.667$ the CPU time savings of block-wise computations instead of full $n$ by $n$ matrix computations would be around 67 \%.

\subsection{Comparing of  path following methods with the global Francis QR eigen algorithm for finding the FoV of indecomposable  and decomposable matrices}

How frequent are unitarily decomposable matrices?  15 non-hermitean matrices in  Matlab's test matrix gallery have difficult to compute eigenstructures and eigenvalues.\\[0.5mm]
 The {\tt 'hanowa'} matrix in this gallery is non-hermitean and normal for all dimensions. Thus it is  unitarily diagonalizable and  benefits from our block-diagonalization algorithm. We note that five out of these 15 matrices or one third are unitarily block-diagonalizable.\\[0.5mm]
  The {\tt 'redheff'} matrix block-diagonalizes for every $n$. Its maximal block size $m$ is usually near its dimension $n$ and there are generally a few additional blocks  of relatively small dimensions such as  2 by 2 or 1 by 1. The FoVs of  Redheffer's small diagonal blocks seem to be  contained within the FoV of its largest diagonal block as depicted below for $n = 27$.\\[-6mm]
\begin{center}
\includegraphics[width=106mm]{F7Redheff27.png} \\[-1mm] Figure  7
\end{center}
Matlab's {\tt 'clement'} matrices of variable sizes $n$ by $n$ always block-diagonalize into two almost equal sized blocks.\\[-5mm]
\begin{center}
\includegraphics[width=100mm]{F8Clement17.png} \\[-1mm] Figure 8
\end{center}
And their two diagonal block FoVs  are contained within each other.\\[0.5mm]
Matlab's $n$ by $n$ gallery matrix {\tt 'binomial'}  unitarily block-diagonalizes into  maximally 4 by 4 blocks when $n$ is even and into maximally 2 by 2 diagonal blocks when $n$ is odd, with additional 1 by 1 blocks so that the block dimensions add to $n$.  \\[-5.5mm]
\begin{center}
\includegraphics[width=100mm]{F9binom16.png} \\[-1mm]
Figure 9
\end{center}
 Again the  FoVs of the diagonal blocks are nested within each other.\\[0.5mm]
 The Matlab gallery matrix $A$ = {\tt 'invol'}  is not normal and it achieves uniform maximal size 2 by 2 block-diagonalizations for all dimensions $n$, with one extra 1 by 1 diagonal block for odd $n$. Its diagonal block FoV curves seem circular except for the smallest ones. They all seem to be centered at zero and the radii of the larger ones grow rapidly as $n$ increases. When $n = 15$ for example, the maximal radius of its largest 2 by 2 diagonal block FoV curve  measures around $6\times 10^{10}$ from its center at zero. Its second largest FoV  curve is also nearly circular with center near zero and it has the approximate radius of $3 \times 10^8$.\\[1mm]
 More than a quarter of relevant Matlab gallery  test matrices  are non-normal and three of them or 20 \% can be unitarily  block-diagonalized with maximal block dimensions $m = \al \cdot n$ with $\al < 0.25$. This indicates that their $\partial F(A)$ points can be computed via any path following algorithm and even block QR eigen computations in 10 \% or less of the CPU time of a  complete $n$ by $n$  Francis QR at $O(n^3)$ cost according to the speed gain graph in Figure 6. \\[1mm]
 All of our computations and plots in this paper were obtained using the Matlab m-files for unitary  block-diagonali-\\
 zation in \cite{FUMatrixDecomp} and for the field of values evaluations in  \cite{FUDecompMatrFoV}.\\[2mm]
Our ZNN based path following FoV algorithm \cite{FUZNNFoV}  searches for a unitary block-diagonalization $U^*AU$ of a given dense matrix $A\in \CC_{n,n}$ and  the related hermitean matrix flow ${\cal F}_A(t) = \cos(t)H + \sin(t)K$ with $H = (A+A^*)/2$ and $K=(A-A^*)/(2i)$ if that can be achieved. Unitary block decompositions of  matrix flows were introduced in \cite{FUDecomp}.   If  a proper block-diagonalization $U^*AU$ of $A$ with $U^*U = I_n$ has been found, $A$'s field of values boundary curve is the convex hull of the field of values of the individual diagonal blocks $\tilde A_j$  of $U^*AU$. Note  that $F(U^*AU) = F(A)$ for all unitary matrices $U$ and all $A \in \CC_{n,n}$.\\[1mm]
Ours is an unusual approach here. Recall that Francis' backward stable QR algorithm and Matlab's {\tt eig} m-file 'diagonalize' every complex matrix  $A$ by finding a non-singular 'eigenvector matrix' $V$ and a diagonal 'eigenvalue matrix' $D$ so that $\tilde A \cdot V = V \cdot D$ with $\tilde A \approx A$, i.e., they solve an adjacent eigenproblem  with $\tilde A$ in a backward stable way for $A$.  Instead we find a unitary $U$ and a proper block-diagonal matrix representation $\hat A$ for any complex matrix $A$ precisely, except for rounding errors, so that $A \cdot U = U \cdot \hat A$ if that is possible. In order to find the FoV of  a decomposable static matrix $A$ efficiently via 'divide and conquer' block methods, we rely on our knowledge of matrix flow decompositions for the related hermitean flow ${\cal F}_A(t)$ of $A$. Thus parameter-varying matrix methods help us  to plot fields of values of static matrices $A$ that block-decompose under unitary similarities. In fact  matrix decompositions for parameter-varying matrix flows now enable us for the first time to use fast eigencurve path following methods such as ZNN without concerns of eigencurve crossings. Remember that eigencurve crossings do not occur with indecomposable hermitean matrix flows due to Hund and von Neumann and Wigner \cite{FH1927, NW29}. Consequently there are no eigencurve crossings and thus no eigenvalue 'lead changes' within ${\cal F}_{\hat A_j}(t)$ for any indecomposable diagonal block $\hat A_j$ of $\hat A = U^*AU$.\\[1mm]
The FoV  codes    in \cite{FUDecompMatrFoV} first decide on the unitary decomposability of $A$. If $A$ is unitarily decomposable into  $r$ diagonal blocks $\hat A_j$ for $j = 1,...,r > 1$,  each diagonal block $\hat A_j$  is checked first for normalcy. If normal, then the FoV of $\hat A_j$  is the convex hull of $\hat A_j$'s eigenvalues.  In this case we add the eigenvalues of the normal block $\hat A_j$ to the depository of potential boundary FoV points of $A$. If $\hat A_j$ is not normal, our algorithm computes discrete  FoV boundary points of $\hat A_j$ via ZNN using  a convergent look-ahead finite difference formula of type {\tt k\_s = 4\_5}. It then adds the computed FoV boundary points for $\hat A_j$ to the depository of potential  FoV boundary points. Upon termination for all  individual blocks, i.e.,  when we know  all sub-FoVs, our algorithm uses the convex hull algorithm of  Matlab to plot the convex hull of the depository $\partial  F(A)$ point list in $\CC$.  Throughout this paper we use a convergent look-ahead   finite difference formula of type {\tt k\_s = 4\_5}  that uses $k+s = 9$ start-up values and has the truncation error order of $O(\tau^{k+2}) = O(\tau^6)$  (see \cite{FUDiscrForm19}) and  the sampling gap $\tau = 0.0002$ while  varying the exponential decay constant $50 \leq \eta \leq 240$ to achieve 15 accurate leading digits for nearly $ 32,000$ discrete FoV boundary points for $A$ that we compute.\\[2mm]
Our run-time experiments below are limited to $n = 250$  and $n = 1000$. For $n = 250$ we limit the maximal diagonal block size $m$ in our test matrices $A_{n,n}$ to $m = 10, 40, 80, 120, 160, 180 <n$ and $m = n =  250$ in turn. For $m =10$ and $n = 250$ for example, we  start with  a random entry, block-diagonal matrix $B_{250,250}$ composed of 15 blocks of size $m = 10$, 10 blocks of size $m = 8$, and 5 blocks of size $m = 4$. Then we obscure the block structure of $B$ and form the dense test matrix  $A = U^*BU$ with the identical FoV as $B$ by using a random entry 250 by 250 unitary similarity $U$ on $B$. Thereafter the actual FoV finding algorithm starts from the dense matrix $A$. First it  retrieves the original block structure of $B$ from $A$ via one hermitean eigenanalysis of  ${\cal F}_A(t_a)$ and a logical 0-1 pattern analysis for a second hermitean flow matrix ${\cal F}_A(t_b) \neq {\cal F}_A(t_a)$  under unitary similarity $V^*{\cal F}_A(t_b)V$ for the diagonalizing unitary eigenvector matrix $V$ of ${\cal F}_A(t_a)$.
This method has been established  in \cite{FUDecomp} and it is  extremely fast,  taking around 0.05 sec of CPU time for general complex matrices $A_{250,250}$ and around 0.9 sec for 1000 by 1000 matrices, both  indecomposable or decomposable ones.  Once we have discovered the  hidden block structure of our test matrices $A$ numerically,  we adapt  our  ZNN eigencode from \cite{FUZNNFoV} to evaluate the eigendata and $\partial F(\hat A_j)$ points of each diagonal block $\hat A_j$ of the block-diagonal realization of $A$ separately. We do this sequentially  for each $j$.
Finally the built-in convex hull  function of Matlab selects the FoV boundary points for $A$ itself.\\[1mm]
Figure 10 shows  30 individual block FoV curves for  our dense but decomposable example matrix $A_{250.250}$ as specified above. To plot the FoV boundary graph for this $A$, our method uses 31,426 points  out of 942,810 computed potential FoV boundary curve points for $A$'s 30 individual sub-blocks. The whole process takes 20.7  seconds of CPU time.
Figure 11 shows the block-diagonalization zero-nonzero Matlab  {\tt spy} pattern for ${\cal F}_A(t)$ that is identical to that of the unitarily block-diagonalized version of  $A$ and similar to the block structure  of the original $B_{250,250}$. Using the global Johnson type Matlab {\tt eig} $n$ by $n$ matrix QR based method  takes between 85 and 89 sec for any  maximal diagonal block size $m$ when $n = 250$. The global QR and {\tt eig} based process {\tt wber3FoV2.m} in \cite{FUDecompMatrFoV} computes 32,427 $\partial F(A)$ points and results in a FoV boundary curve whose point coordinates  agree everywhere in their leading 15 digits with those that we have computed via ZNN. Figure 12 finally depicts the computed $\partial F(A)$ curve for $A$ and  our $ m = 10$  example. \\[-5mm]
\enlargethispage{10mm}
\begin{center}
\includegraphics[width=140mm]{F10allblkFoVs250_30.png} \\[-1mm]
Figure  10:  Showing all 30 diagonal block FoV boundaries of the decomposable  $A_{250,250}$ with their starting points $\cdot$
\end{center}\enlargethispage{40mm}

\vspace*{-4.5mm}
\begin{center}
\includegraphics[width=85mm]{F11diag01patt250_30a.png} \\[-1mm]
Figure  11: Showing the block-diagonal structure of ${\cal F}_A(t_b)$ under unitary similarity with $V(t_a)$, and thus of $A$\\
\end{center}

\begin{center}
\includegraphics[width=90mm]{F12FoV250_30all.png} \\[0mm]
Figure  12: FoV boundary curve of our decomposable test matrix  $A_{250,250}$ with maximal block dimension 10,\\ as the convex hull of the partial FoV curve points of Figure 10
\end{center}

The average CPU run times to draw a $\partial F(A)$ curve such as depicted in Figure 12 above  for a given 250 by 250 dense and unitarily decomposable matrix $A$ with increasing maximal diagonal block dimensions $m$   are listed in Table 1.

\begin{center}\hspace*{-1mm}
\begin{tabular}{c||c|c|c|c|c|c|c}
$\mathbf{A_{n,n}}$ {\bf with}  $\mathbf{n = 250}$&\multicolumn{7}{c} {Max block size $\mathbf m$ for dense unitarily block decomposable  $A_{250,250}$} \\[1mm] \hline \hline
Method  for $\partial F(A)$&10&40&80&120&160&180&250\\ \hline
individual block path following ZNN &20.7  sec & 26  sec& 25  sec&  35 sec&  41 sec&  43 sec &   55 sec\\
 &&\multicolumn{5}{c|}{\emph{ZNN speed-up factors compared to {\tt eig}} :} &\\
 &{\bf 4.3} $\times$&{\bf 3.3} $\times$& {\bf 3.5} $\times$&{\bf 2.5} $\times$&{\bf 2.1} $\times$&{\bf 2} $\times$&{\bf 1.6} $\times$\\ \hline
Johnson QR {\tt eig} based& 89 sec  &86 sec & 88 sec &89 sec &87 sec &85 sec & 87 sec
\end{tabular}\\[2mm]
Table 1
\end{center}

\vspace*{-2mm}
The actual run times for our ZNN path following method and FoV computations with fixed $n$ and $m$ may vary by up to around 10 \% up or down, depending on the distribution of block sizes below the maximal block dimension $m$.\\
Both our path following ZNN method and the Francis complete  QR  eigendata method can be easily parallelized in Matlab by changing their outermost  'for' loop into a 'parfor' loop that initiates parallel processing. The sequential 'for' codes  in \cite{FUDecompMatrFoV} end with  {\tt ...for.m} and  and their parallelized 'parfor' codes with {\tt ...parfor.m}, respectively. Running both parallelized code versions on the same problem in parallel mode on a 2019 MacBook Pro with a 2.4 GHz 4  core i5 processor and 16 GB of memory  speeds up each of the computations by a factor of around three and the speed advantages of individual block path following methods over complete eigendata methods such as {\tt eig} stay about the same.\\[1mm]
Our next data table deals with run time comparisons for constructing  the FoV of general square matrices $A_{n,n}$ of sizes $n =250$ and $n = 1000$ that are either dense (top data rows) or in the lower rows allow a unitary block decomposition into six diagonal blocks of the same proportions for either dimension $n$.\\[-3mm]

\begin{center}\hspace*{-3mm}
\begin{tabular}{c|r|c||c|r|c|c}
$\mathbf{A_{n,n}}$&\multicolumn{2}{c||}{$\mathbf{n = 250}, \ m = 250$}&$\mathbf{A_{n.n}}$&\multicolumn{3}{c}{$\hspace*{-5mm}\mathbf{n = 1000},\ m = 1000$}\\
{\bf dense}&CPU times&\emph{speed-up}&{\bf dense}&CPU time &\emph{speed-up} &$k^x$ ops count\\
&&\emph{vs} {\tt eig}&&&\emph{vs} {\tt eig}& with $k = 4$
\\[1mm] \hline \hline\\[-5mm] &&&&&\\[-2mm]
Johnson QR {\tt eig}& 87 sec&   &Johnson QR {\tt eig}&3014 sec&    &$k^{2.6}$\hspace*{1.8mm}\\
single ZNN&{ 55 sec}&{  1.6}  $\times$&single ZNN&{ 1407 sec}& { 2.1}  $\times$ &${k^{2.34}}$\\[1mm] \hline \hline \hline \\[-4mm]
& \multicolumn{2}{c}{ } & &\multicolumn{3}{c}{}\\[-3mm]
$\mathbf{A_{n,n}}$&\multicolumn{2}{c||}{$\mathbf{n = 250}, \ m = 100$}&$\mathbf{A_{n.n}}$&\multicolumn{3}{c}{\hspace*{-5mm}$\mathbf{n = 1000}, \ m = 400$}\\
{\bf decomposing} into&CPU times&\emph{speed-up}&{\bf decomposing} into&CPU time&\emph{speed-up}& $k^x$ ops count\\
{[}100,80,30,30,8,2{]}&&\emph{vs} {\tt eig}&{[}400,320,120,120,32,8{]}&&\emph{vs} {\tt eig}& with $k = 4$ \\[1mm] \hline \hline\\[-5mm]
&&&&&\\[-2mm]
Johnson QR {\tt eig}& 87 sec&&Johnson QR {\tt eig}&3014 sec&     &$k^{2.6}$\hspace*{1.8mm}\\
indiv. block ZNN&26.2 sec&3.3 $\times$&indiv. block ZNN&{ 296 sec}& { 10.2}  $\times$ &${k^{1.75}}$\\
\end{tabular}\\[1.5mm]
Table 2
\end{center}

\vspace*{-2mm}
For real or complex,  dense and unitarily indecomposable matrices $A_{n,n}$ with maximal block dimension $m = n$ and large dimensions $n \gg 25$ our ZNN based path following method is faster than Matlab's Francis implicit multi-shift QR eigendata finder {\tt eig} for plotting FoVs. At $n = 25$ both methods run about equally fast at 1.78 seconds. For $n = 10$ Matlab's fully compiled {\tt eig} needs 0.2 seconds and our path finding method computes $A$'s FoV boundary curve in 0.8 seconds.
Recall from \cite{FUZNNFoV} that ZNN methods reduce matrix eigen computations to solving one linear systems plus a linear vector recursion at each angle or time step at much lower $O(n^3)$ cost than QR which relies on matrix factorizations with their intrinsically  higher $O(n^3)$ costs.\\[0.5mm]
 Progressing from $A_{250,250}$ to $A_{1000,1000}$ with $O(n^3)$ complexity methods for finding FoVs, the data in the right columns  of Table 2 should increase the  operations count  by a factor around the \emph{dimension increase factor $k$ cubed} where $k = 1000/250 = 4$ here, i.e., by around $k^3$ or $4^3$. In fact going to $k=4$ times larger  test matrices  $A$ the operations count increases for the QR based FoV method nominally by less, with an $O(k^{2.6})$ complexity increase in the top right Table 2 entry instead of $O(k^3)$. Note that the operations complexity  increase for single ZNN and FoV plotting is even lower at $O(k^{2.34})$ for indecomposable matrices when going from $n = 250$ to $n = 1000$.  For decomposable matrices $A_{n,n}$ with comparable block decomposition sizes, going from $n=250$ to $n =1000$ improves the ops count even more for our ZNN path finding method from an expected value of $4^3$ to the exceedingly low $4^{1.75}$ ops count as displayed by the bottom right entry of Table 2. \\[2mm]
A dense 'artificial' matrix $A_{52,52}$  with complicated  block structure  results in Matlab from the command \\
{\tt B = blkdiag(-2*eye(2),gallery('forsythe',6), gallery('jordbloc',8,1-1i),...\\
\hspace*{26mm}zeros(3),gallery('hanowa',8),wilkinson(12),hilb(9),...\\
\hspace*{26mm}(gallery('jordbloc',4,1+1i))');} \\
after a dense random entry unitary similarity as $A = Q^* \cdot B \cdot Q$.
 Matlab's {\tt eig} and our ZNN based path following method (after a block-diagonalization of the dense matrix $A$) plot the $\partial F(A)$  identically as the quadrilateral shown in Figure  13. The respective run times are 2.9 seconds for the {\tt eig} based Johnson QR method and around 2.07 seconds for the individual block  ZNN  method. Note that in the block-diagonalization of the dense test matrix  $A$ only three of its 37 diagonal blocks have dimensions larger than 1 -- one coming from 'forsythe' and two  from the  'jordbloc's of $B$. Otherwise $A$ is  unitarily diagonalizable. The eigenvalues of 'hanowa' lie on the vertical line through $-1 \in \CC$ and the eigenvalues of all other blocks in $A$ are real, denoted by dots on the real axis of Figure 13. The two Jordan blocks of $B$ with the shared eigenvalue $1 - i \in \CC$ generate  concentric circles as sub-fields of values boundary curves for $A$ as is well known. Their radii depend on their dimensions, see Figure 14. \\
 For comparison's sake, the convex hull computation from almost 1 million ZNN computed FoV points for  the decomposing $A_{250,250}$ matrix with $m = 10$ of Table 1 took about 0.08 seconds while the same from around 63,000 FoV points for the unitarily decomposable matrix $A_{52,52}$ takes around 0.05 seconds.

\begin{center}
\includegraphics[width=76mm]{F13nonsenseA52FoV.png} \\[0mm]
Figure 13 : FoV boundary curve of the decomposing test matrix  $A_{52,52}$
\end{center}

\vspace*{-3mm}
\begin{center}
\includegraphics[width=88mm]{F14B-d52b.png} \\[0mm]
Figure 14 : Partial FoV boundary curves of the decomposable test matrix  $A_{52,52}$ depicted under {\tt FOVZNN4\_5aeigshortFoV3for.m} \ in \cite{FUDecompMatrFoV}
\end{center}

A rather debilitating problem arose sometimes unpredictably with computational inaccuracies when replacing $H = H^*$ and $K=K^*$ by the  unitarily similar matrices $V^*HV$ and $V^*KV$ for the diagonalizing matrix $V$ of ${\cal F}_A(t_a)$ that was computed according to the algorithm of \cite{FUDecomp}. Very small un-hermitean  entry errors in $V^*HV$ and $V^*KV$  by quantities of machine constant \emph{eps} size would occasionally result in Matlab's {\tt eig} function  interpreting a partial sub-block $H_j$ or $K_j$ of the block-diagonal matrices $V^*HV$ or $V^*KV$ as non-hermitean and then compute the necessary start-up eigendata  of ${\cal F}_{A_j}(t_k) = \cos(t_k)H_j + \sin(t_k)K_j$ with  nearly correct real parts  but tainted by minuscule imaginary parts. Unfortunately however, in this case Matlab's {\tt eig} function's automatic return of ordered  eigenvalues for hermitean matrices  was lost and some computed partial sub FoVs became wild crisscrosses. This made them  useless as  starting values for path following iterations and for finding the final  convex hull FoV for a decomposing matrix $A$. To insure that  all theoretically hermitean diagonal blocks $H_j$ and  $K_j$ are recognized and treated correctly as hermitean in Matlab's {\tt eig} function in the start-up phase for ZNN we  set $\hat H_j = (H_j + H_j^*)/2 = \hat H_j^*$ and $\hat K_j = (K_j + K_j^*)/2 =\hat K_j^*$ before every partial sub-block FoV evaluation and worked with ${\cal \hat F}_{A_j}(t_k) = \cos(t_k)\hat H_j + \sin(t_k)\hat K_j$ instead. \\[-6mm]

\section{Outlook}\vspace*{-1mm}

Charlie Johnson \cite{J78} may or may not have realized four decades ago that his use of rotating Bendixson  rectangles \cite{Be}  for the field of values problem of a fixed entry static matrix $A$ had  changed the FoV problem fundamentally into a parameter-varying hermitean matrix flow problem. The Francis implicit multi-shift QR  method {\tt eig}  is generally considered the most accurate and fastest algorithm for any eigendata computation of square matrices of dimensions not exceeding  $n = 10,000$ by much. Here it evaluates the extreme eigenvalues  of the  hermitean  flow matrices ${\cal F}_A(t)$ of $A_{n,n}$ and records their associated eigenvectors' $x$ actions $x^*Ax \in \CC$   as discrete $ F(A)$ boundary curve points for $A$. Time- or parameter-varying matrix flows were not part of the Linear Algebra canon nor standard  knowledge then.\\[1mm]
 Even today,  time-varying matrix flow problems have almost exclusively been studied and used in  real-time engineering projects and their computational treatments are just beginning to appear in preprints by numerical analysts and matrix theoreticians. Currently there is only a very thin knowledge base for these problems in  western numerical analysis circles. Yet there are well over 400 papers in engineering journals on Zhang Neural Network (ZNN)  methods for time-varying matrix problems that deal with a wide range  of   on-chip and real-time applications in industry. There is a handful of books on this new, yet mostly theoretically unstudied matrix computational subject, see \cite{FUZNNsurvey} e.g..\\[2mm]
Loisel and Maxwell \cite[sect. 5, 6.2, 7.1, 7.2, 8.3; p. 1733 - 1743]{LM18} have searched for ways to solve the FoV problem for decomposing matrices and found certain seemingly natural limitations of IVP ODE path following FoV solvers.
Their observations in \cite[last sentence on p. 1743]{LM18} on this issue end in  \\
"\emph {... We were also not able to completely analyze nonnormal matrices of type 3 (when $\lambda_{\text{max}}(t)$ can be nonsimple). These limitations to our analysis are to be expected: even for the problem of computing eigenvalues, state-of-the-art eigenvalue solvers fail for some matrices.}"\\[1mm]
  Yes, our current state-of-the art  matrix eigensolvers such as Francis' multi-shift implicit QR generally  fail for some static matrices $A$ due to seemimgly  inherent limitations for matrices with  non-trivial Jordan structures or repeated eigenvalues. Likewise  naive path-following FoV methods must fail  for decomposable  matrices $A$ because of potential eigencurve crossings of the associated hermitean matrix flow  ${\cal F}_A(t)$.\\[1mm]
   But such 'limitations'  need not be 'expected' for the FoV problem at all as we were told; they are not germane for these problems  -- if we can  enlarge our horizon and understandings. Simple parameter-varying matrix flow ideas have been developed while   studying time-varying matrix problems, see \cite{FUCoalesc} and \cite{FUDecomp} for example. These ideas have now  helped  us to resolve the well 'documented', but in fact  needless FoV  path finder method limitations in a radical new way. We  have shown  how to compute general static matrix FoVs of every square matrix $A \in \CC_{n,n} \text{ or } \RR^n $, decomposable or not,  -- and how to do so accurately and efficiently -- with  the new and fast ZNN  path-following method when it is  combined with the ultra-fast matrix block-decomposition method of \cite{FUDecomp}. Besides, the matrix and matrix flow unitary decomposition results of \cite{FUDecomp}  that helped us do so here can easily be adapted as  'preconditioner' for other path finding methods and other unitarily invariant matrix problems such as  SVD computations, to least squares problems, to finding the numerical radius or Crawford number of a matrix,  and so forth.\\[1mm]
  Maybe some of the many currently open questions regarding  time-varying matrix flows, see e.g.  \cite{FUDiscrForm19}, \cite{FUDiscrForm19}, \cite{FUZNNsurvey}, can be solved in turn by using  deep or simple insights and understandings from  classic fixed entry static matrix analysis, from its concepts, and numerical methods. \\
   I do hope so.


TITLE: Perturbation theory for the generalized eigenvalue problem
QUESTION [9 upvotes]: Is there a standard reference for the perturbation theory of the generalized eigenvalue problem?
More specifically, I would like to get a systematic expansion for the problem

$(A_0 + \epsilon A_1)v = \lambda B v$

in terms of the solutions to

$A_0 v = \lambda B v$

where $A_0$, $A_1$ and $B$ are known, $n\times n$ Hermitian matrices,  $\epsilon >0 $ is a "small" parameter, and $\lambda\in \mathbb{C}$ and $v\in\mathbb{C}^n$ are the unknowns. When $B$ is positive-definite, the usual (Rayleigh-Schroedinger) perturbation series can easily be generalized by using $B$ to define an alternative inner product. However, in my problem, $B$ is not positive-definite (although it is still nonsingular).

REPLY [6 votes]: In addition to those already suggested by Federico, a good reference is Chapter 2 in Kato's book Perturbation Theory for Linear Operators, especially if you want to tackle the question from the higher perspective of functional analysis. Another reference, addressing the same problem from the alternative point of view of the calculus of variations (with an eye on computations), is I. Babuška and J. Osborn, Eigenvalue Problems, in Handbook of Numerical Analysis (Part 1), Vol. II, ed. by P.G. Ciarlet and J.L. Lions, pp. 641-787 (1991).


TITLE: Point of Borel Sets
QUESTION [0 upvotes]: One definition of a Borel set is as follows:
If $S$ is a topological space (where
it supports the concepts of open set): , then the Borel $\sigma-$ algebra on $S$ , denoted by $\mathbb{B}(S)$
, is the $\sigma-$ algebra generated by all the open sets.
I hope my question is not too broad- what is the advantage of working with a Borel Sigma algebra, as opposed to a generic sigma algebra?

REPLY [0 votes]: For example: suppose that a $\sigma $-algebra contains a topology (this is the case of Borel). Then every real valued continuous function is measurable.


TITLE: Generating functions (Counting?)
QUESTION [2 upvotes]: Determine how many ways Brian, Katie, and Charlie can split a $\$50$ dinner bill such that Brian and Katie each pay an odd number of dollars and Charlie pays at least $\$5$.

How can I use generating functions to solve this problem? Or should I use a counting approach?

REPLY [2 votes]: The generating function for the number of ways that the three can pay a bill of $n$ dollars is
$$
\begin{align}
&(x+x^3+x^5+\dots)(x+x^3+x^5+\dots)(x^5+x^6+x^7+\dots)\\[3pt]
&=x^7\frac1{(1-x^2)^2}\frac1{1-x}\\
&=x^7\frac1{(1-x)^3}\frac1{(1+x)^2}\\
&=x^7\left[\frac14\frac1{(1-x)^3}+\frac14\frac1{(1-x)^2}+\frac3{16}\frac1{1-x}+\frac18\frac1{(1+x)^2}+\frac3{16}\frac1{1+x}\right]\\
&=\sum_{k=0}^\infty\small\left[\frac14(-1)^k\binom{-3}{k}+\frac14(-1)^k\binom{-2}{k}+\frac3{16}(-1)^k\binom{-1}{k}+\frac18\binom{-2}{k}+\frac3{16}\binom{-1}{k}\right]x^{k+7}\\
&=\sum_{k=0}^\infty\small\left[\frac14\binom{k+2}{2}+\frac14\binom{k+1}{1}+\frac3{16}\binom{k}{0}+\frac18(-1)^k\binom{k+1}{1}+\frac3{16}(-1)^k\binom{k}{0}\right]x^{k+7}\\
&=\sum_{k=0}^\infty\frac1{16}\left[(2k^2+10k+11)+(-1)^k(2k+5)\vphantom{\frac1{16}}\right]x^{k+7}\\
&=\sum_{k=0}^\infty\left[\binom{\lfloor k/2\rfloor}{2}+2\binom{\lfloor k/2\rfloor}{1}+\binom{\lfloor k/2\rfloor}{0}\right]x^{k+7}\\
&=\sum_{k=0}^\infty\binom{\lfloor k/2\rfloor+2}{2}x^{k+7}\tag{1}
\end{align}
$$
Plugging $k=43$ into $(1)$ yields the coefficient of $x^{50}$ to be $253$.


TITLE: Find $x$-coordinates of all horizontal tangents on the graph of $f(x) = \sin^2x + \cos x$
QUESTION [1 upvotes]: The function is $f(x) = \sin^2x + \cos x$.
I found the derivative which was $f'(x) = 2x\cos^2x - \sin x$.
I think what you do next is make $f'(x) = 0$ so it becomes:
$0 = 2x\cos^2x - \sin x$.
After this, I'm not sure how to proceed.

REPLY [0 votes]: $\sin^2 x$ means $(\sin x)^2$, not $\sin (x^2)$
so $(\sin^2 x)' =2\sin x\cos x$


TITLE: Given $abc=1$ and $0< c \leq b \leq1\leq a$, prove that $8(a+b+c)^2\le9(1+a^2)(1+b^2)(1+c^2)$
QUESTION [1 upvotes]: I can't make progress with proving this inequality.
I have tried opening the brackets and using $abc=1$ in order to obtain the following:
$$a^2+b^2+ \frac 1{a^2b^2}+18+9a^2b^2+\frac 9{a^2}+\frac9{b^2}\ge 16 \left (ab+\frac 1a+\frac 1b \right)$$

REPLY [0 votes]: We begin by observing that if $a=b=c=1$, then the two expressions are equal to each other.
Expand everything. Use the equation $abc=1$ to get rid of compound terms. Gather terms with each variable separately. Use the preceding observation to factor each nominator.
$$\begin{align} & 9(1+a^2)(1+b^2)(1+c^2)-8(a+b+c)^2 \\ &=9\bigl(1+a^2+b^2+c^2+a^2b^2+b^2c^2+c^2a^2+a^2b^2c^2\bigr)-8\bigl(a^2+b^2+b^2+2ab+2bc+2ca\bigr) \\ &= 9\Bigl(1+a^2+b^2+c^2+{1\over c^2}+{1\over a^2}+{1\over b^2}+1\Bigr)-8\Bigl(a^2+b^2+c^2-{2\over c}-{2\over a}-{2\over b}\Bigr) \\ & ={a^4+6a^2-16a+9 \over a^2}+{b^4+6b^2-16b+9 \over b^2}+{c^4+6c^2-16c+9 \over c^2} \\ & ={(a-1)^2(a^2+2a+9) \over a^2}+{(b-1)^2(b^2+2b+9) \over b^2}+{(c-1)^2(c^2+2c+9) \over c^2}\end{align}$$
It is now obvious that this last expression is always $\geq 0$, and that it is $=0$ only in the case when $a=b=c=1$.


TITLE: Is graph automorphism Karp-reducible to graph isomorphism under hidden subgroup representation?
QUESTION [2 upvotes]: The classical representations of the graph automorphism problem is Karp-reducible to the classical representation of the graph isomorphism problem. The sketch of proof for this can be written as follows based on the detailed proof given in Section 1.5 of The Graph Isomorphism Problem - Its Structural Complexity.

Prove that the graph automorphism problem is Turing-reducible to the graph isomorphism problem.
Prove that the graph automorphism problem has a polynomial time computable or-function.
Prove that the graph isomorphism problem has both polynomial time computable and- and or-functions.
Combining the previous results, it may be argued that the graph automorphism problem is Karp-reducible to the graph isomorphism problem.

My question is whether similar reduciblity can be inferred for the problems under hidden subgroup representations (HSP). Let me define those representations.
Graph isomorphism as a hidden subgroup problem
Let the $2 n$ vertex graph $\Gamma = \Gamma_1 \sqcup \Gamma_2$ be the disjoint union of the two graphs $\Gamma_1$ and $\Gamma_2$ such that $Aut \left(\Gamma_1\right) = Aut \left(\Gamma_2\right) = \left\{e\right\}$ (a limited version). A map $\varphi : S_{2n} \to \text{Mat}\left(\mathbb{C}, N \right)$ from the group $S_{2n}$ is said to have hidden subgroup structure if there exists a subgroup $H_\varphi$ of $S_{2n}$, called a hidden subgroup, an injection $\ell_\varphi : S_{2n}/H \to \text{Mat}\left(\mathbb{C}, N \right)$, called a hidden injection, such that the diagram

is a commutative diagram, where $S_{2n}/H_{\varphi}$ denotes the collection of right cosets of $H_\varphi$ in $S_{2n}$, and where $\nu : S_{2n}/H_\varphi$ is the natural map of $S_{2n}$ onto $S_{2n}/H_\varphi$. We refer to the group $S_{2n}$ as the ambient group and to the set $\text{Mat}\left(\mathbb{C}, N \right)$ as the target set.
The hidden subgroup version of the graph isomorphism problem is to determine a hidden subgroup $H$ of $S_{2n}$ with the promise that $H$ is either trivial or $|H| = 2$.
Graph automorphism as a hidden subgroup problem
For a graph $\Gamma$ with $n$ vertices, a map $\varphi : S_{n} \to \text{Mat}\left(\mathbb{C}, N \right)$ from the group $S_{n}$ is said to have hidden subgroup structure if there exists a subgroup $\text{Aut}\left(\Gamma\right)$ of $S_{n}$, called a hidden subgroup, an injection $\ell_\varphi : S_{n}/\text{Aut}\left(\Gamma\right) \to \text{Mat}\left(\mathbb{C}, N \right)$, called a hidden injection, such that for each $g \in \text{Aut}\left(\Gamma\right)$, $g \left(\Gamma\right) = \Gamma$ and, the diagram

is commutative, where $S_{n}/\text{Aut}\left(\Gamma\right)$ denotes the collection of right cosets of $\text{Aut}\left(\Gamma\right)$ in $S_{n}$, and where $\nu : S_{n}/\text{Aut}\left(\Gamma\right)$ is the natural map of $S_{n}$ onto $S_{n}/\text{Aut}\left(\Gamma\right)$. We refer to the group $S_{n}$ as the ambient group and to the set $\text{Mat}\left(\mathbb{C}, N \right)$ as the target set.
The hidden subgroup version of the graph automorphism problem is to determine a hidden subgroup $\text{Aut}\left(\Gamma\right)$ of $S_{n}$ with the promise that $\text{Aut}\left(\Gamma\right)$ is either of trivial or non-trivial order depending on the type of $\Gamma$.
My questions:
I think it can be trivially shown that the hidden subgroup representation of the graph automorphism problem is Turing-reducible to the graph isomorphism problem by giving two input graphs as the original graph and the candidate automorphism of the original graph.
If I am correct, the remaining question is whether there is a Karp-reduction for the hidden subgroup representation. If there is, how it can be shown.

REPLY [2 votes]: I think your focus on the rigid case of GI limits you too much. Instead phrase (non-rigid) GI as an HSP in the same way, but now the goal is to determine the size of the hidden subgroup, or a generating set. The difference between the isomorphic and non-isomorphic cases will be a factor of 2 in the order of the hidden subgroup. If you phrase the problem as finding generators of the hidden subgroup, then the question is just whether any generator switches $\Gamma_1$ and $\Gamma_2$.
Now, an instance of "HSP-GA" corresponding to a graph $\Gamma$ is given by the function from $S_n \to M_n(\mathbb{C})$ defined by $f(\pi) = A(\pi(\Gamma))$ where $A(\cdot)$ denotes the adjacency matrix. In particular, $f(1) = A(\Gamma)$. Then apply the usual Karp reduction from GA to GI to get a pair of graphs $\Gamma_1, \Gamma_2$, and consider the HSP-GI instance, of the type described in the preceding paragraph, corresponding to the pair $\Gamma_1, \Gamma_2$ (that is, the disjoint union $\Gamma_1 \cup \Gamma_2$).


TITLE: What is the intuition behind the definition of the kernel of a homomorphism
QUESTION [7 upvotes]: I was starting to study some algebra (groups and homomorphisms in particular) and came across the definition of the kernel (for a group-homomorphism $f:G \rightarrow G'$):
$$\ker(f) = \{ x \in G \mid f(x) = e' \}$$
where $e'$ is the identity. I understand the definition but fail to appreciate its importance or what it means intuitively or the motivation to define such a set/object. So far I understand that it's the set of all input elements that once processed by our homomorphism, yield the identity under the other group $(G', e')$, however, I fail to really appreciate why mathematicians would bother to define (and give it a name) to such a set.

What is the motivation to define such as set" If it has a name, it must have some importance.
Why is it so special to map to the identity element of a set? What structure does it reveal about the relationship of both sets? What is the intuition? Why would mapping to the identity be useful?
What is the conceptual understanding of it beyond regurgitating/repeating the definition?

REPLY [9 votes]: There are many many reasons why kernels are important to group theory, but here's just one way of appreciating the kernel in a fairly isolated context.
If we zoom out a bit, any set-function $f: A \to B$ (here $A$ and $B$ are simply sets) naturally partitions $A$ into equivalence classes, and for $a \in A$, the equivalence class of $a$ is given by
$$[a] = \{a' \in A : f(a') = f(a)\};$$
the set of all elements of $A$ that get mapped to the same thing as $a$. The same logic applies if $f : G \to H$ isn't just a set function, but a homomorphism of groups.
With the equivalence class notation, the kernel of $f$ is simply the equivalence class of the identity $1_G$ of $G$,
$$\ker f = [e_G],$$
since any homomorphism $f: G \to H$ always sends the identity $e_G$ of $G$ to the identity $e_H$ of $H$. What can we say about arbitrary $g, g' \in G$ such that $f(g) = f(g')$? That is, what can we say about the equivalence class $[g]$ for any $g \in G$?

Claim: For a homomorphism $f: G \to H$ and $g \in G$, we have $f(g) = f(g')$ if and only if there exists some $k \in \ker f$ such that $gk = g'$; that is, $g$ and $g'$ differ by a multiple of something in the kernel of $f$. In particular, $[g] = \{gk: k \in \ker f\}$, and has size $|\ker f|$.
($\Longrightarrow$) Supposing $f(g) = f(g')$, note that there exists a unique $g^* = g^{-1}g' \in G$ such that $gg^* = g'$. Then
$$f(g) = f(g') = f(gg^*) = f(g)f(g^*),$$
and left-multiplication by $f(g)^{-1}$ shows that $e_H = f(g^*)$, hence $g^* \in \ker f$.
($\Longleftarrow$) Homework.

If you've ever heard homomorphisms described as functions that "respect" the group operation(s), the size of the kernel is a measure of just how "respectful" a given homomorphism is! A large kernel means that more of the structure of the group $G$ is "ignored" when transported to the group $H$.
Edit:
For "respectful", imagine two situations, considering $S_3$, the symmetric group of degree $3$. There's a sign homomorphism $\operatorname{sgn} : S_3 \to \{-1, 1\} = C_2$ to the multiplicative group $C_2$ sending each permutation to its sign. Its kernel is the alternating group $A_3 = \{1, (123), (132)\}$ of "even" permutations, and in the image $C_2$, almost all of the structure of $S_3$ is ignored; we forget everything but whether a permutation is even or odd.
On the other hand, we have a "copy" of $S_3$ as a subgroup of $S_4$ if we consider all permutations of $S_4$ that leave $4$ fixed. This leads to an "inclusion" homomorphism $\iota: S_3 \to S_4$, sending each permutation to its "copy" in $S_4$. This inclusion homomorphism has only the identity of $S_3$ in its kernel, and is considerably more "respectful" than the sign homomorphism; every bit of information about $S_3$ shows up in the image $\iota(S_3)$.


TITLE: Euler's method approximation
QUESTION [2 upvotes]: I'm supposed to write a program for approximating the value of function $y = y(x)$, which is given as: $$y' = \frac{1+y}{1 + x^2}$$  I also know that $y(0) = 0$. I should approximate the value for $x=2$.
Solving the differential equation gets me that $y = Ke^{\arctan(x)} -1$, where $K\in\mathbb{R}^+$
So I can get the value of $K$ by solving $0 = Ke^{\arctan(0)} - 1$. Therefore $K = 1$.
But I'm not sure if the approximation should use the function in the original shape or the result of the differential equation.
If I use $y = e^{\arctan(x)} - 1$ for the aprroximation, the results are ($h$ is the size of step):
For $h=1$ it is $2.0$
For $h=0.5$ it is $2,02885...$
For $h=0.25$ it is $2,03366...$
The precise value is $y(2) = e^{\arctan(2)} - 1 = 2.0257189050$
Code used to calculate the results:
DWORD ApproximateExampleOne(){
    double stepSize = 0.0;
    double result = 0.0;
    double iterator = 0.0;
    double preciseValue = exp(atan(2.0)) - 1;

    fflush(stdin);
    _tprintf(_T("Approximating equation [y = e^(arctan(x)) - 1] for x = 2 by Euler's method\n\n"));
    _tprintf(_T("Please enter size of the step [#.#]: "));
    _tscanf(_T("%lf"), &stepSize);

    if(stepSize <= 0){
        _tprintf(_T("Size of the step can't be zero or negative! Aborting... \n"));
        return ERROR_INVALID_PARAMETER;
    }

    for(iterator = stepSize; iterator < 2.0; iterator += stepSize){
        result = result + stepSize*(exp(atan(iterator)) - 1);
    }

    _tprintf(_T("\nPrecise result is: %.10f\n"), preciseValue);
    _tprintf(_T("Approximate result is: %.10f\n"), result);
    _tprintf(_T("Approximation error is: %2.0f%%"), (fabs(preciseValue - result)/fabs(preciseValue))*100);

    return ERROR_SUCCESS;
}

REPLY [2 votes]: See, when you're using Euler's method, the assumption there is that you're unable or unwilling to find a symbolic solution to your differential equation, and you just want to estimate the function values of a solution at a few points.
Remember that Euler's method for $y^\prime=f(x,y)$ takes the form
$$\begin{align*}x_{k+1}&=x_k+h\\y_{k+1}&=y_k+hf(x_k,y_k)\end{align*}$$
where $h$ is a predetermined stepsize and $x_0,y_0$ correspond to your initial conditions.
That's what you need to do here: pick a stepsize $h$, let $x_0=y_0=0$ (due to your initial conditions), and then keep running Euler (replacing $f(x,y)$ with whatever's equated to the derivative) up until $x_{k+1}=2$. The corresponding value of $y_{k+1}$ is the Euler estimate, which I presume you'll be asked to compare with the exact solution in this case...


TITLE: If $(x_{2n})_n, (x_{2n+1})_n$ and $(x_{3n})_n$ converges, prove that $(x_n)_n$ also converges.
QUESTION [1 upvotes]: let's say $(x_n)_n$ is a row in $\mathbb{R}$. If $(x_{2n})_n, (x_{2n+1})_n$ and $(x_{3n})_n$ converges, I want to prove that $(x_n)_n$ also converges.
I think it's normal because $(x_n)_n$ exist out of the three rows, but I can't prove it. Can someone help me with this?

REPLY [3 votes]: Suppose that $(x_{2n})$ converges to $R$, $(x_{2n+1})$ converges to $S$ and $(x_{3n})$ converges to $T$.
Consider the sequence $(x_{6n})$: it's a subsequence of $(x_{2n})$, so it converges to $R$. But it's also a subsequence of $(x_{3n})$, so it converges to $T$, showing that $R=T$. Similarly, $(x_{6n+3})$ is a subsequence of both $(x_{2n+1})$ and $(x_{3n})$, showing that $S=T$, and hence $R=S$.
Finally, note that the sequences $(x_{2n})$ and $(x_{2n+1})$ partition the sequence $(x_{n})$ into two disjoint subsequences, both converging to the same limit $R$, implying that $(x_{n})$ converges to $R$.

REPLY [2 votes]: Let $l_1 = \lim\limits_{n \to \infty} x_{2n+1}$, $l_2 = \lim\limits_{n \to \infty} x_{2n}$ and $l_3 = \lim\limits_{n \to \infty} x_{3n}$
It is sufficient to prove that $l_1 = l_2$ which immediately follows from
$$l_1 = \lim\limits_{n \to \infty} x_{2(3n + 1)+1} = \lim\limits_{n \to \infty} x_{3(2n+1)} = l_3 = \lim\limits_{n \to \infty} x_{3(2n)}=\lim\limits_{n \to \infty} x_{2(3n)}=l_2$$


TITLE: Find the value of $a$ and $b$ in $ F(x) = a + b \arcsin x $
QUESTION [1 upvotes]: Given $X$ is a continuous random variable and its probability distribution function is
$$F(x)=
\begin{cases}
0, & x < -1, \\
a+b\arcsin x, & -1 \le x < 1, \\
1, & x \ge 1
\end{cases}$$
I find the value of $b$ in this way.
$$F'(x) = f(x), \int_{-\infty}^{+\infty} f(x)\, dx = 1, \text{ so } b= \frac1\pi$$
How can I find the value of $a$?

REPLY [2 votes]: There is no discrete mass at $1$. So
$$\lim_{x\to 1^-}(a+b\arcsin x)=1.$$
There is also no discrete mass at $-1$. Use that to get another equation.