How to solve numerical range¶

In previous tutorial, we have defined the numerical range of a matrix $A$ as the set of all possible values of $x^*Ax$ for all unit vectors $x$. In this tutorial, we will focus on different methods to solve the numerical range of a matrix.

In [1]:

import numpy as np
import matplotlib.pyplot as plt

try:
    import numqi
except ImportError:
    %pip install numqi
    import numqi

np_rng = np.random.default_rng()

import numpy as np import matplotlib.pyplot as plt try: import numqi except ImportError: %pip install numqi import numqi np_rng = np.random.default_rng()

Root finding method¶

In paper "Bounding Real Tensor Optimizations via the Numerical Range" arxiv-link appendix A. For a given angle $\theta$, the boundary can be find by solving the root of a scalar function. (TODO add more details)

In [2]:

N0 = 5
matA = np_rng.normal(size=(N0,N0)) + 1j*np_rng.normal(size=(N0,N0))
alpha_list = np_rng.uniform(0, 2*np.pi, size=5)

boundary = numqi.matrix_space.get_matrix_numerical_range(matA, num_point=200)
max_list = []
min_list = []
for alpha_i in alpha_list:
    min_list.append(numqi.matrix_space.get_matrix_numerical_range_along_direction(matA, alpha_i, kind='min')[0])
    max_list.append(numqi.matrix_space.get_matrix_numerical_range_along_direction(matA, alpha_i, kind='max')[0])
min_list = np.array(min_list)
max_list = np.array(max_list)

fig,ax = plt.subplots()
ax.plot(boundary.real, boundary.imag)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
for ind0,alpha_i in enumerate(alpha_list):
    tmp0 = np.sqrt(max(abs(xlim[0]), abs(xlim[1]))**2 + max(abs(ylim[0]), abs(ylim[1]))**2)
    ax.plot([0,tmp0*np.cos(alpha_i)], [0, tmp0*np.sin(alpha_i)], 'k')
    ax.plot([0,tmp0*np.cos(alpha_i+np.pi)], [0, tmp0*np.sin(alpha_i+np.pi)], 'k:')
tmp0 = max_list*np.exp(1j*alpha_list)
ax.plot(tmp0.real, tmp0.imag, 'rx', markersize=10)
tmp0 = min_list*np.exp(1j*alpha_list)
ax.plot(tmp0.real, tmp0.imag, 'r+', markersize=10)
ax.set_xlim(*xlim)
ax.set_ylim(*ylim)
fig.tight_layout()

N0 = 5 matA = np_rng.normal(size=(N0,N0)) + 1j*np_rng.normal(size=(N0,N0)) alpha_list = np_rng.uniform(0, 2*np.pi, size=5) boundary = numqi.matrix_space.get_matrix_numerical_range(matA, num_point=200) max_list = [] min_list = [] for alpha_i in alpha_list: min_list.append(numqi.matrix_space.get_matrix_numerical_range_along_direction(matA, alpha_i, kind='min')[0]) max_list.append(numqi.matrix_space.get_matrix_numerical_range_along_direction(matA, alpha_i, kind='max')[0]) min_list = np.array(min_list) max_list = np.array(max_list) fig,ax = plt.subplots() ax.plot(boundary.real, boundary.imag) xlim = ax.get_xlim() ylim = ax.get_ylim() for ind0,alpha_i in enumerate(alpha_list): tmp0 = np.sqrt(max(abs(xlim[0]), abs(xlim[1]))**2 + max(abs(ylim[0]), abs(ylim[1]))**2) ax.plot([0,tmp0*np.cos(alpha_i)], [0, tmp0*np.sin(alpha_i)], 'k') ax.plot([0,tmp0*np.cos(alpha_i+np.pi)], [0, tmp0*np.sin(alpha_i+np.pi)], 'k:') tmp0 = max_list*np.exp(1j*alpha_list) ax.plot(tmp0.real, tmp0.imag, 'rx', markersize=10) tmp0 = min_list*np.exp(1j*alpha_list) ax.plot(tmp0.real, tmp0.imag, 'r+', markersize=10) ax.set_xlim(*xlim) ax.set_ylim(*ylim) fig.tight_layout()

Above, the blue boundary points are obtained by performing calculating for all angle (see tutorial/numerical-range for more details). The black lines are some desired angle at which to calculate the boundary points. The red points are the boundary points obtained by solving the root of a scalar function. One can see that this method also gives the correct boundary points. And if the number of desired angle is not too much, one can anticipate that this method will be more efficient.