この内容は，Numerical Determination of the Field of Values of a General Complex Matrixを読んで実装したものです．
恐らくこの論文のアルゴリズム内の$\theta_j$は，$\theta_j = (j-1) / (2 \pi / k)$ではなく，$\theta_j = (j-1) * (2 \pi / k)$だと思います．
下記の内容に間違っている記述がある場合は，申し訳ありません．是非とも教えて頂けると嬉しいです．

線形方程式$Ax=b$のconvergence boundsを考えます．
一般には$A$の固有値で抑えたり，固有値分布で評価します．
またfield-of-valuesを用いて抑える手法もあります．
field-of-valuesは他にnumerical range，Wertovorrat，Hausdorff domain，range of valuesとも呼ばれるようです．
(Wikiではnumerical rangeですね)

field-of-valuesの定義は次式です． これは長さ$1$の単位ベクトルのRayleigh quotientのrangeと見れます． $$ F(A) = {x^* A x, x^*x = 1, x \in \mathbb{C}^n } $$

field-of-valuesは次の性質を持ちます．

  1. $F(A)$ is convex, closed and bounded.
  2. $\sigma(A) \subseteq F(A)$. ($\sigma(A)$ is set of eigenvalues.)
  3. $F(U^* A U) = F(A)$ for any unitary matrix.
  4. $F(A + zI) = F(A) + z$ and $F(zA) = zF(A)$ for any complex number $z$.
  5. If $A$ is hermite, then $F(A)$ is a segment of the real line $[\lambda_{min}, \lambda_{max}]$.
  6. If $A$ is normal, then $F(A)$ is $\text{Convex}(\sigma(A))$.

また，$A$のHermitian part($H(A)$)の最大固有値は，$F(A)$の境界に乗ります．
そのため，この最大固有値に対応する固有ベクトル$x$を用いて，$x^* A x$を計算すれば，$F(A)$の境界点が求まります．

この動作を$F(A)$を角度$\theta$で回転を繰り返し，境界点を求めていきます．
$F(A)$の$\theta$回転は$e^{i\theta}F(A)$ですが，これは$F(e^{i \theta}A)$なので(性質4)，$H(e^{i \theta}A)$の最大固有値，最大固有値に対応する固有ベクトルが分かれば良くなります．

求めた$F(A)$の境界点で凸多角形を作り，それを$F_{in}(A)$とします．
$F_{in}(A)$は必ず$F(A)$に内包されます．

また$F(A)$がconvexである条件を用いると(性質1)，回転角度$\theta_j$で求めた境界点$u$に対して，$F(A)$は直線$u+i\theta_{j}$を跨がないはずなので，$F(A)$を内包する凸多角形$F_{out}(A)$が作れます．

必ず$F(A)$は，$F_{in}(A)$と$F_{out}(A)$に挟まれているので，$F(A)$が存在する範囲を求めることができて
また，回転角度$\theta_j$を無限に小さくすれば$F(A)=F_{in}(A)=F_{out}(A)$になります．

問題行列としてGRCAR(10, 5)で，$[0, 2\pi]$区間の分割数を$3 \sim 20$まで変えていった場合のプロットを乗せます．
青点はGRCARの固有値です．
赤点と赤点線で表される凸多角形は$F_{in}(A)$です．
青点と青点線で表される凸多角形は$F_{out}(A)$です．
つまり$F(A)$は赤と青の点線の間にあります．

$[0, 2\pi]$区間の分割数の終了条件は，次の面積比と，基準値$\epsilon$を決めます．
$$ \frac{\text{area}(F_{out}(A)) - \text{area}(F_{in}(A))}{\text{area}(F_{out}(A))} \leq \epsilon $$

また収束性解析では，$F(A)$の中で原点と最も距離が小さいものを$\nu(F(A))$として，$\nu(F(A))$と$\nu(F(A^{-1}))$を用います． $$ \frac{\|r_{n}\|}{\|r_{0}\|} \leq \left(1 - \nu(F(A)) \nu(F(A^{-1}))\right)^{\frac{n}{2}} $$

今回の例では$F(A)$に$(0, 0)$が入ってしまって，上手く見積もりとして用いられません．
その場合は多項式$\|p(A)\|$を用いて見積もるようです．

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as LA

def H(B):
    B_star = np.conjugate(B.T)
    return 0.5 * (B + B_star)

def lambda_M(B):
    eig, vec = LA.eig(B)
    index = np.argmax(eig)
    return eig[index], vec[:, index]


def field_of_values(A, number_of_points, eps):
    n = len(A)

    # =========== 1. if A is hermitian
    is_hermite = True
    for i in range(n):
        for j in range(i+1, n):
            if abs(A[i, j] - A[j, i]) <= eps:
                continue
            is_hermite = False

    if is_hermite:
        eig, _ = LA.eig(A)
        poly = []
        poly.append(np.min(eig))
        poly.append(np.max(eig))
        return np.array(poly)

    # =========== 2. if A is normal
    is_normal = True
    AAT = np.dot(A, A.T)
    ATA = np.dot(A.T, A)
    
    for i in range(n):
        for j in range(n):
            if abs(AAT[i, j] - ATA[i, j]) <= eps:
                continue
            is_normal = False
    
    if is_normal:
        eig, _ = LA.eig(A)
        return eig

    # =========== 3.
    k = number_of_points  # Choose a k-mesh
    theta = []
    for j in range(0, k):
        theta.append((2*np.pi/k) * j)
    
    ps = []
    lambdas = []
    for j in range(len(theta)):
        # ===== Determine Fin: p_theta
        A_theta = np.exp(1j * theta[j]) * A
        HA = H(A_theta)
        lambda_theta, x_theta = lambda_M(HA)
        p_theta = np.dot(np.dot(np.conjugate(x_theta), A), x_theta)

        lambdas.append(lambda_theta)
        ps.append(p_theta)

    qs = []
    for j in range(len(theta)):
        # ===== Determine Fout: q_theta
        theta_one = theta[j]
        theta_two = theta[(j+1) % len(theta)]
        delta = theta_two - theta_one
        lambda_one = lambdas[j]
        lambda_two = lambdas[(j+1) % len(lambdas)]

        tmp = lambda_one + ((lambda_one * np.cos(delta) - lambda_two) / np.sin(delta)) * 1j
        q_theta = np.exp(-1j * theta[j]) * tmp
        qs.append(q_theta)
    
    return np.array(ps)

参考文献