• I: Introduction to Numerical Range
• II: The Algorithms
  • W: Produces a graph of the W(M) by connecting points found on the boundary of W(M)
  • WLines: Produces a graph of W(M) bounded by support lines
  • NumRadius: Gives the greatest distance from the origin to a point in W(M).
• III. Procedure Code (With Comments)
  • W(M::'Matrix'(square))
  • WLines(M::'Matrix'(square))
    • OpNorm(M::Matrix)
    • WLines(M::'Matrix'(square))
  • NumRadius(M::'Matrix'(square), n::posint)
• IV. Examples
  • A. Normal Matrix
  • B. Elliptic Numerical Range
  • C. 3 x 3 Example
  • D. Hermitian Example
  • E. Infinite-Dimensional Applications
    • Approximating the Numerical Range and the Norm of a Composition Operator
      C[phi];
• V. Conclusions
• VI. References
• VII. Acknowledgment
• VIII. Disclaimer

Numerical Ra nge

2001 Paul Bourdon and Tina Harbilas
Washington and Lee University
USA

ABSTRACT: This worksheet contains procedures that compute the numerical range and

numerical radius of a square matrix.

I: Introduction to Numerical Range

The numerical range of the n x n matrix M denoted W( M ) is given by

W( M ) = { < M v , v > : || v || =1 }

where < , > denotes the usual inner product on . Thus,

< a , b > =

where a = ( , , ... , ), b = ( , . ... , ), and the bar represents complex conjugation.

The numerical range of a matrix M is sometimes called the field of values of M .

Observe that the numerical range of the matrix M is a closed bounded subset of the complex plane containing all of the eigenvalues of M . The Toeplitz-Hausdorff Theorem states the numerical range of any matrix is a convex subset of the plane.

The numerical range of a 2 x 2 matrix is an ellipse; for example, the numerical range of is an ellipse with foci and having minor axis length of .

The numerical range of a Hermitian matrix, that is a matrix M such that M equals its conjugate transpose M *, contains only real numbers. More exactly, the numerical range of a Hermitian matrix M is the line segment extending from the smallest eigenvalue of M to the largest eigenvalue of M .

Related to the numerical range is the numerical radius , which is the greatest distance between something in the numerical range and the origin, or

= max { | w | : w in W( M ) }

where W( M ) is the numerical range of the matrix M .


Gustafson and Rao's monograph Numerical Range ([2]) is a good general reference.

II: The Algorithms

This worksheet contains two procedures that find numerical ranges and a third procedure that finds numerical radii. The algorithms are similar in all three procedures. The numerical range procedures work with the plots and LinearAlgebra packages, taking a square Matrix M as the input and producing a graph of the numerical range of M as the output. The numerical-radius procedure, which also calls on the LinearAlgebra package, takes a square Matrix M and returns the numerical radius of M .

Note that the WLines algorithm is not ideal for finding the numerical range of a Hermitian matrix. See Example IV.D below. Additional discussion of these algorithms may be found in references [2] and [3].

W: Produces a graph of the W( M ) by connecting points found on the boundary of W( M )

Given a square matrix M , Maple finds the largest eigenvalue of ( M + M *)/2. This eigenvalue, which will be real because ( M + M *)/2 is Hermitian, represents max{Re( w ) : w in W( M )}. Maple finds a corresponding eigenvector v of length one for the maximum eigenvalue. By the definition of numerical range,

= < M v , v >

belongs to W( M ); moreover, it must be a point on the boundary of W( M ) because there is no number in W( M ) having larger real part.

For a positive integer n (default value 50), Maple repeats the process described in the preceding paragraph with M replaced by its rotate , where takes values for j =1.. n . Connecting all the calculated points, , j = 0.. n , produces a picture of the numerical range. The larger the value of n , the more accurate the representation of the numerical range. Note that the procedure W draws the boundary of the numerical range; the numerical range consists of the boundary and its interior.

WLines: Produces a graph of W(M) bounded by support lines

Given a square matrix M and a positive integer n (default value 36), Maple finds the largest eigenvalue of

L[ ] := ( M + M *) / 2

for j = 1.. n , where = . The largest eigenvalue of L[ ] represents the largest real part of any number in the rotated numerical range W( M ). For each j from 1 to n , Maple draws a vertical line segment through the largest eigenvalue of L[ ] and then rotates the segment by to its appropriate position. Every point in W( M ) will either lie on the rotated segment or to one side of it. The result is a picture of the numerical range bounded by n support lines. The larger the value of n , the more accurate the representation of the numerical range.

NumRadius: Gives the greatest distance from the origin to a point in W(M).

Given a square matrix M and a positive integer n , Maple finds the set of largest eigenvalues of the matrices

L[ ] := ( M + M *) / 2

for j = 1.. n , where = . The maximum eigenvalue of this set is returned as the numerical radius. Larger values of n typically result in a more accurate determination of numerical radius.

III. Procedure Code (With Comments)

> restart;

> with(LinearAlgebra):

> with(plots):

Warning, the name changecoords has been redefined

Note: The following procedures were developed under Maple 6.

W(M::'Matrix'(square))

Calling Sequence:

W( M , n , plot options)

Parameters:

M ::'Matrix'(square) - Procedure returns a plot of the boundary of the numerical range of M

n ::posint - (optional) n is the number of boundary points plotted; default value is 50

Plot Options - (optional) Any of Maple's plot options such as color=blue; default plot settings match those of Maple 2d plots.

Discussion:

For example, W( M , 75, scaling=constrained, color=magenta) produces a magenta plot of 75 boundary points of the numerical range of M (connected by line segments) with units along the real and imaginary axes having the same length. The command W( M , scaling=constrained, color=magenta) would produce approximately the same output except that 50 boundary points would be plotted.

> W:=proc(M::'Matrix'(square))

#DECLARE LOCAL VARIABLES
local n,p,r,T,DM,R,i,j,E,H,G,B,k,J,V,W,points;

#PROVIDE FOR OPTIONAL ARGUMENTS
#n IS NUMBER OF BOUNDARY POINTS PLOTTED (DEFAULT VALUE 50)
#PLOT OPTIONS START WITH args[p]
if nargs > 1 then
if type(args[2],posint) then
n:=args[2];
p:=3;
else
n:=50;
if (type(args[2], equation)) then
p:=2;
else
WARNING("Second argument should either be a positive integer or a plot option");
p:=3;
fi;
fi;
else
n:=50;
fi;

#FORCE MAPLE TO USE FLOATS
M[1,1]:=evalf(M[1,1]);

#NEXT LOOP IMPLEMENTS ALGORITHM DESCRIBED IN SECTION II
for r from 0 to n do
T:=MatrixAdd((evalf(exp(2*Pi*I*r/n))*M), (HermitianTranspose(evalf(exp(2*Pi*I*r/n))*M)))/2;
DM:=RowDimension(T);
R:=Matrix(1..DM,1..DM,shape=hermitian,datatype=complex(float));
for i from 1 to DM do

#T[i,i] SHOULD BE REAL
#GUARD AGAINST ROUNDING ERRORS BY TAKING REAL PART
R[i,i]:=Re(T[i,i]);

od;
for i from 1 to DM do
for j from i+1 to DM do
R[i,j]:=T[i,j];
od;
od;
E:=Eigenvectors(R)[1]; #R IS HERMITIAN SO E-VALUES REAL
H:=Eigenvectors(R)[2];
G:=Dimension(E);
B:=max(seq(E[k], k=1..G));

#GET INDEX J FOR MAX E-VALUE B
for k from 1 to G do
if testfloat(E[k], B, 5)=true then
J:=k;
fi;
od;

#V IS E-VECTOR FOR MAX E-VALUE B
V:=Column(H,J);

#W[r] IS BOUNDARY POINT CORRESPONDING TO r-TH ROTATE
W[r]:=evalf(Multiply(HermitianTranspose(V), Multiply(M,V)));
od;

#RETURN PLOT OF BOUNDARY POINTS, CONNECT=TRUE
points:= [seq([Re(W[r]), Im(W[r])], r=0..n)];
if nargs>1 then
RETURN(display(pointplot(points, connect=true, seq(args[k],k=p..nargs))));
else
RETURN(display(pointplot(points, connect=true)));
fi;
end proc:

WLines(M::'Matrix'(square))

Calling Sequence:

WLines( M , n , plot options)

Parameters:

M ::'Matrix'(square) - Procedure returns a plot of the numerical range of M bounded by suport lines.

n ::posint - (optional) n is the number of support lines plotted; default value is 36

Plot Options - (optional) Any of Maple's plot options such as color=blue; default plot settings match those of Maple 2d plots.

Discussion:

For example, WLines( M , 50, scaling=constrained, color=magenta) produces a plot of 50 magenta support lines surrounding the numerical range of M , with units along the real and imaginary axes having the same length. The command WLines( M , scaling=constrained, color=magenta) would produce approximately the same output except that 36 support lines would be plotted.

The WLines procedure relies on the following OpNorm procedure to determine lengths of segments of support lines plotted. The OpNorm procedure must be executed/compiled before running WLines.

OpNorm(M::Matrix)

Calling Sequence:

OpNorm( M )

Parameters:

M ::Matrix - Procedure returns the operator norm of the matrix M

Discussion:

The operator norm of the matrix M , is defined by

max{|| M v ||: || v || = 1 },

where || || denotes the standard norm on (giving the length of the vector that is its argument). It's easy to prove the operator norm of M is given by

.

OpNorm(M) computes the operator norm of the matrix M by using the preceding formula. Maple's Norm( M , 2) computes the operator norm of M using the same formula; however, Maple's Norm( , 2) procedure does not take advantage of the fact that the product matrix

HermitianTranspose( M ) M

is Hermitian.

Maple's Eigenvalues routine converges more reliably when its argument is a Hermitian matrix. Thus OpNorm( M ) produces accurate results in some cases where Maple's Norm( M , 2) does not. See Example IV.E below.

> OpNorm:=proc(M::Matrix)

#DECLARE LOCAL VARIABLES
local Dim,S,W,i,j,E,N;

#DIMENSION OF M*M WILL BE COLUMNDIMENSION(M)
Dim:=ColumnDimension(M);

#SET S TO M*M
S:=Matrix(1..Dim,1..Dim,shape=hermitian);
W:=evalf(Multiply(HermitianTranspose(M),M));
for i from 1 to Dim do S[i,i]:=Re(W[i,i]) od;
for i from 1 to Dim do
for j from i+1 to Dim do
S[i,j]:=W[i,j];
od;
od;

#GET OPERATOR NORM = SQRT(MAX(EIGENVALUES(S))
E:=Eigenvalues(S);
N:=max(seq(E[i],i=1..Dim));

#ALL E-VALUES OF S SHOULD BE NON-NONGETATIVE
#GUARD AGAINST ERRORS
if N < 0 then
error "Operator Norm could not be calculated"
else
RETURN(sqrt(N))
fi;
end proc:

WLines(M::'Matrix'(square))

> WLines:=proc(M::'Matrix'(square))

#DECLARE LOCAL VARIABLES
local n,p,N,r,T,DM,R,i,j,E,G,A,F;

#PROVIDE FOR OPTIONAL ARGUMENTS
#n IS NUMBER OF SUPPORT LINES PLOTTED (DEFAULT VALUE 36)
#PLOT OPTIONS START WITH args[p]
if nargs > 1 then
if type(args[2],posint) then
n:=args[2];
p:=3;
else
n:=36;
if (type(args[2], equation)) then
p:=2;
else
WARNING("Second option should either be a postive integer or a plot option");
p:=3;
fi;
fi;
else
n:=36
fi;

#FORCE MAPLE TO USE FLOATS
M[1,1]:=evalf(M[1,1]);

#N = OpNorm(M) DETERMINES LENGTH OF SUPPORT LINES
N:=OpNorm(M);
if type(N, float) = false then
error "Make sure OpNorm procedure has been loaded."
fi;

#NEXT LOOP IMPLEMENTS ALGORITHM DESCRIBED IN SECTION II
for r from 1 to n do
T:=MatrixAdd(evalf(exp(2*Pi*I*r/n))*M,(HermitianTranspose(evalf(exp(2*Pi*I*r/n))*M)))/2;
DM:=RowDimension(T);
R:=Matrix(1..DM,1..DM, shape=hermitian);
for i from 1 to DM do R[i,i]:=Re(T[i,i]); od;
for i from 1 to DM do
for j from i+1 to DM do
R[i, j]:=T[i,j];
od;
od;
E:=Eigenvalues(R);#R IS HERMITIAN SO E-VALUES ARE REAL
G:=Dimension(E);
A:=max(seq(E[k],k=1..G));

#F[r] IS r-TH SUPPORT LINE
F[r]:=pointplot([[Re(exp(-2*Pi*I*r/n)*(A+N*I)), Im(exp(-2*Pi*I*r/n)*(A+N*I))],\ [Re(exp(-2*Pi*I*r/n)*(A-N*I)), Im(exp(-2*Pi*I*r/n)*(A-N*I))]], connect=true);
od;

#RETURN PLOT OF NUMERICAL RANGE BOUNDED BY SUPPORT LINES
if nargs > 1 then
RETURN(display(seq(F[k], k=1..n),seq(args[k], k=p..nargs)));
else
RETURN(display(seq(F[k], k=1..n)));
fi;
end proc:

NumRadius(M::'Matrix'(square), n::posint)

Calling Sequence:

Numradius( M , n )

Parameters:

M ::'Matrix'(square) - Procedure returns the numerical radius of M

n :: posint - n is the number of points in the boundary of the numerical range of M whose distances to the origin are computed.

Discussion:

The maximum of the n distances computed is returned as an approximation for the numerical radius of M . Larger values of n typically produce more accurate results; recommended starting value for n is 36. For example, NumRadius( M , 36) produces an approximation of the numerical radius of M using 36 boundary points. Note well the second argument is required.

> NumRadius:=proc(M::'Matrix'(square), n::posint)

#DECLARE LOCAL VARIABLES
local r,T,DM,R,i,j,E,J,A;

#FORCE MAPLE TO USE FLOATS
M[1,1]:=evalf(M[1,1]);

#NEXT LOOP IMPLEMENTS ALGORITHM DESCRIBED IN SECTION II
for r from 1 to n do
T:=MatrixAdd((evalf(exp(2*Pi*I*r/n))*M), (HermitianTranspose(evalf(exp(2*Pi*I*r/n))*M)))/2;
DM:=RowDimension(T);
R:=Matrix(1..DM,1..DM,shape=hermitian);
for i from 1 to DM do R[i,i]:=Re(T[i,i]); od;
for i from 1 to DM do
for j from i+1 to DM do
R[i,j]:=T[i,j];
od;
od;
E:=Eigenvalues(R);
J:=Dimension(E);
A[r]:=max(seq(Re(E[k]),k=1..J));
od;

#RETURN APPROXIMATE NUMERICAL RADIUS
RETURN(max(seq(A[k], k=1..n)));
end proc:

IV. Examples

Note: At this point, Maple's LinearAlgebra and plots packages have been loaded as well as the numerical -range and numerical-radius procedures of Section III.

A. Normal Matrix

The numerical range of a normal matrix is the convex hull of the eigenvalues of the matrix.

> M1:=Matrix([[0, 0, 1, 0, 0], [1, 0, 0, 0, 0], [0, 0, 0, 0, 1], [0, 1, 0, 0, 0], [0, 0, 0, 1, 0]]);

> W(M1);

> WLines(M1, 30, scaling=constrained, color=magenta);

> NumRadius(M1,36); #Exact Value is 1

B. Elliptic Numerical Range

The numerical matrix of a 2 x 2 matrix is an ellipse.

> M2:=Matrix([[-1,1],[0,1]]);

> WLines(M2, scaling=constrained);

> NumRadius(M2,36);

C. 3 x 3 Example

The numerical range of the following matrix clearly contains 3I as well as the elliptic numerical range W( M2 ) of the matrix M2 in the preceding example. A little more thought reveals that the numerical range is the convex hull of 3I and W( M2 ).

> M3:=Matrix([[3*I,0,0],[0,-1,1],[0,0,1]]);

> W(M3, scaling=constrained);

> NumRadius(M3,36);

D. Hermitian Example

The numerical range of a Hermitian matrix is a line segment joining the smallest eigenvalue of the matrix to the largest.

> M4:=Matrix([[3,-I],[I,1]]);

> WLines(M4,24);

> W(M4);

> W(M4, axes=framed, color=green);

> NumRadius(M4,1);

E. Infinite-Dimensional Applications

Suppose that T is a continuous linear operator on the Hilbert space H (complete, inner-product space) and that H has denumerable orthonormal basis B = { , ...}. The definitions of numerical range, numerical radius, and operator norm apply, without modification, to T .

The matrix for T relative to B has entries

= < T , >

i = 1..infinity, j = 1..infinity. Let be the orthogonal projection of H on the linear span of = { , ... , }. Let be the linear transformation from to defined by = , and let be the matrix for relative to . It's easy to see that W( ) is contained in W( T ) and not difficult to prove that W( ) converges to the closure of W( T ) in the Hausdorff metric. Also, OpNorm( ) converges to OpNorm( T ).

Composition Operators on the Hardy Space

The Hardy space consists of those analytic functions on the unit disk U in the complex plane whose Taylor coefficients in the expansion about the origin are square summable. The Hardy space is a Hilbert space with inner product

< , > = .

Observe that B = {1, z , , , ...} is an orthonormal basis for .

Suppose that is a holomorphic function mapping the open unit disk in the complex plane into itself. For example, might be given by = . Any such mapping induces a continuos linear operator on via composition: f = f( ). (Continuity is not obvious, but follows from a theorem by J.E. Littlewood. For references concerning composition operators and their numerical ranges, the reader may consult [1].) The matrix for relative to B has i - j entry given by the Maple command

coeftayl( , z = 0, i -1).

Approximating the Numerical Range and the Norm of a Composition Operator

Define the "symbol" of the composition operator.

> phi:=z -> 2*(3-3*z-I*sqrt(3)+2*I*sqrt(3)*z)/(9-6*z-I*sqrt(3)+2*I*sqrt(3)*z);

The preceding function is a conformal automorphism of the unit disk fixing .

The next execution group defines the 60 x 60 submatrix for relative to the basis { : n = 0..59}. Many digits are required to get accurate entries.

> Digits:=40:
EntryIJ:=(i,j) -> evalf(coeftayl(phi(z)^(j-1), z=0, i-1)):
M:=Matrix(60,EntryIJ):

Now generate an increasing sequence of approximations for W( ).

> for i from 1 to 6 do
P[i]:=W(SubMatrix(M,1..i*10,1..i*10)):
od:

display(seq(P[i],i=1..6),scaling=constrained,title=`10 by 10 to 60 by 60 approximations in steps of 10`);

Now we turn to the norm approximation. The operator norm of is known to be .

> NormCphi:=evalf(sqrt((1+abs(phi(0)))/(1-abs(phi(0)))));

> for i from 1 to 6 do
OpNorm(SubMatrix(M,1..i*10,1..i*10));
od;

The sequence above is increasing toward NormCphi, consistent with theory.

> for i from 1 to 6 do
Norm(SubMatrix(M,1..i*10,1..i*10),2);
od;

Note that the final value return by Norm( , 2) is inaccurate.

V. Conclusions

Maple's LinearAlgebra package together with its plotting routines facilitates development of procedures for computing numerical ranges and radii of square matrices. As illustrated by the sequence of examples above, output of these procedures is entirely consistent with theorems concerning numerical range. The procedures can be used for investigation of numerical ranges in both the finite-dimensional and infinite-dimensional settings.

VI. References

[1] Bourdon, Paul S. and Shapiro, Joel H., "The Numerical Ranges of Automorphic Composition Operators", Journal of Mathematical Analysis and Applications 251 (2000), 839-854.

[2] Gustafson, Karl E. and Rao, Duggirala K.M., Numerical Range , Springer-Verlag, New York, 1997.

[3] Johnson, Charles R., "Numerical Determination of the Field of Values of a General Complex Matrix", SIAM Journal on Numerical Analysis 15 (1978), 595-602.

VII. Acknowledgment

The work of the authors was supported in part by a grant from the National Science Foundation of the United States through its RUI/REU program. (Grant number: 0100290)

VIII. Disclaimer

While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.