# symmetricTridiagonal.mpl # # A collection of routines to investigate Bohemians # that are complex symmetric and have zero diagonal # # The two purposes of this collection are: # 1. To explore such Bohemians (much is unknown) # 2. To illustrate the use of the Computer Algebra # System Maple. Maple is a registered trademark # of Waterloo Maple Inc. # # ----------------------------------------------------------------------------- # (c) Robert M. Corless 2019 (MIT Release licence) # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in all # copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. # ------------------------------------------------------------------------------- # # Code intended for a "masterclass" at the Isaac # Newton Institute during the programme # Complex Analysis: Tools, techniques, and applications. # # Accompanying material: # IntroWithBohemians.maple # which links to # --- MatrixImagewithevalhf.mw # --- QuestionAboutEigenvalues.mw # --- polyImagewithEvalhf.mw # --- FirstDay.mw # --- this code (symmetricTridiagonal.mpl) # # I thank the Isaac Newton Institute for Mathematical # Sciences for support and hospitality during this # programme when work on this project was undertaken. # This work was supported by EPSRC Grant # EP/R014604/1 # # I am also grateful for the advice of Juergen Gerhard, # David Linder, John May, and Erik Postma of Maplesoft. # # Previous work on Bohemians, which this builds on, # was done in collaboration with Eunice Y. S. Chan, # Steven E. Thornton, Laureano Gonzalez-Vega, # J. Rafael Sendra, and Juana Sendra. Additional # work with Cristian Ardelean and others is in progress. # # See bohemianmatrices.com and the Bohemian Matrices # github repository for more information. # makeMatrix --- make either a symbolic # matrix for a generic entry # in the family, or a particular # matrix in the family # A := makeMatrix( m, [u1, u2, ..., u[m-1]]); # can be called with the result of an iterator # or with a vector or with a list # # demonstrates argument pass-through, loops, # Matrix creation, args, nargs # # Could use sparse data structure but the matrices # we use are so small it does not matter here. # makeMatrix := proc( m::posint, U ) local A, i; # This routine knows what Bohemian # family we are constructing. # We pass in any extra arguments. # This one is complex symmetric, # has zero main diagonal, and is # tridiagonal. Therefore it only # needs to specify a vector of m-1 # elements U[i] from the population. A := Matrix(m, m, 0, args[3..nargs]): for i to m-1 do A[i,i+1] := U[i]; # U could be a symbol A[i+1,i] := U[i]; end do: return A end proc: # getAllEigenvals --- get all eigenvalues numerically # from the Bohemian family # specified by m (dimension) # and P (population) # # Example: # m := 9; # P := {1, I}; # eigvals := getAllEigenvals( m, P ): # # demonstrates parameter type-checking # Iterators package # nested procedures and scoping # datatypes float[8] and complex[8] getAllEigenvals := proc( m::posint, P::list, $ ) local A, eigvals, ev, i, j, k, matrixentries, N, numberMats, u; description "Compute eigenvalues of all complex symmetric tridiagonal Bohemians with zero diagonal"; # Generate all nops(P)^(m-1) possible vectors for the U values matrixentries := Iterator:-CartesianProduct(seq(P, k=1..m-1)): N := (m-1); # Could sample randomly. Here we will do exhaustive computation numberMats := nops(P)^N; eigvals := Vector(m*numberMats,datatype=complex[8]); i := 0; for u in matrixentries do A := makeMatrix(m,evalf[15](u),datatype=complex[8]); ev := LinearAlgebra:-Eigenvalues(A); for j to m do eigvals[i+j] := ev[j]; end do; i := i+m; end do: return eigvals end proc: # # ROUTINES FOR PLOTTING AND CREATING IMAGES # # basicPlot, getImageData, colourizeImageData # # basicPlot --- plot of complex eigenvalues # (just a call to plots[complexplot] # with some options specified) # # demonstrates argument pass-through basicPlot := proc( rtvec ) plots[complexplot]( [seq(rtvec)], style = POINT, symbol = SOLIDCIRCLE, symbolsize = 2, colour = BLACK, axes = BOXED, size = [0.5, "square"], args[2..nargs]); end proc: # getImageData --- get data for a density plot # of eigenvalues in the complex plane # # Meant to be used when the number of eigenvalues is # large enough that plotting them all just gives a # blackened plot. # # Allows for the creation of colourized images # (see below) # # getImageData should be called with evalhf # (will be sometimes fifty times faster) # # Demonstrates evalhf, ImageTools package # # N := 9; # a := lower x limit, e.g. -2; # b := upper x limit, e.g. 2; # c := lower y limit, e.g. -2; # d := upper y limit, e.g. 2; # nrow := 800; # ncol := 800; # more if you want higher res # imageData := evalhf( getImageData( rtvec, N, a, b, c, d, nrow, ncol )): getImageData := proc( rtvec::Vector, N::posint, a::numeric, b::numeric, c::numeric, d::numeric, nrow::posint, ncol::posint ) local dx, dy, i, imageData, j, k, x, y, z; imageData := Array(1..nrow,1..ncol,0.,datatype=float[8]); dx := (b-a)/ncol; # a <= x <= b dy := (d-c)/nrow; # c <= y <= d for k to N do z := rtvec[k]; x := Re(z); y := Im(z); i := floor((y-c)/dy); j := floor((x-a)/dx); if i > 0 and i <= nrow and j>0 and j<= ncol then imageData[i,j] := imageData[i,j]+1.0; end if; end do; return imageData; end proc: # colourizeImage --- use predefined table to # colourize the intensity # Replace intensities with RGB colour chosen according to some # previous scheme (recorded in a table called in this routine "f") colourizeImage := proc( nrow, ncol, imageData::Array, f::Array) local colourizedImageData, i, j, ell; colourizedImageData := Array(1..nrow,1..ncol,1..3); for i to nrow do for j to ncol do for ell to 3 do colourizedImageData[i,j,ell] := f[floor(imageData[i,j]), ell]; od; od; od; return colourizedImageData end proc; colorizeImage := colourizeImage: # # ROUTINES FOR EXACT CHARACTERISTIC POLYNOMIALS # # makePolyFun, makePolyFun2, getAllCharPolys, # getSomeCharPolys, getRoots # # makePolyFun --- generate a function to compute the # characteristic polynomial of the # Bohemian matrix with entries u # uses symbolic computation, code generation # and optimization to brutally make the function. # Surprisingly effective, and quite generalizable. # makePolyFun := proc( m::posint, $ ) local A, F, p, U, u, x; # We first build a symbolic matrix, and use it to generate # a symbolic characteristic polynomial, and turn it # into a function that generates the vector of coefficients # given the entries of the matrix. # A := makeMatrix( m, U ); # U is a symbol # We use that matrix only once. p := LinearAlgebra:-CharacteristicPolynomial(A, x): u := PolynomialTools:-CoefficientVector(p, x); # We make a procedure out of the coefficient vector F := codegen[makeproc](u, parameters = [U]); return codegen[optimize](F, args[2..nargs]); # Optimize == reuse subexpressions end proc: # getSomeCharPolys --- sample characteristic polynomials # at random from the Bohemian family # # Demonstrates the RandomTools package # keyword parameters # default values for parameters getSomeCharPolys := proc( Nsamp::posint, m::posint, P::list, F::procedure, {getReps::boolean:=false}, $ ) local countPolys, i, matrixentries, Matrixrep, N, p, plist, R, r, u; description "Sample and count characteristic polynomial frequencies of all complex symmetric tridiagonal Bohemians with zero diagonal"; # Generate all nops(P)^(m-1) possible vectors for the U values matrixentries := Iterator:-CartesianProduct(seq(P, k=1..m-1)): N := (m-1); # Could do exhaustive computation. Here we will sample randomly. R := RandomTools:-MersenneTwister:-NewGenerator(range=1..nops(P)^N): if Nsamp > nops(P)^N/10 then warning("Sampling more than one tenth of the total. Consider exhaustive computation instead"); end if; countPolys := table(sparse): if getReps then Matrixrep := table(sparse): end if; for i to Nsamp do r := R(); u := Unrank( matrixentries, r ); # Pick an entry with uniform probability p := F(u); #Charpoly of A # The purpose of the following is to give a unique key # into the table; p is a Vector and its name is not # a good key, though its values are (say, as a list) plist := [seq( Normalizer(p[k]), k=1..m+1 )]; # Thx Juergen Gerhard for the "sparse" trick below countPolys[plist] := 1 + countPolys[plist]; if getReps then if Matrixrep[plist]=0 then Matrixrep[plist] := [seq(u)]; end if; end if; end do: if getReps then return eval(countPolys), eval(Matrixrep) else return eval(countPolys) end if; end proc: # getRoots --- compute roots of the collection # of characteristic polynomials # Demonstrates fsolve, modular arithmetic # uses factors to get multiplicities exactly # uses the output of getAllPolys or getSomePolys. # # Demonstrates table indices, and doubly-indexed # for-loops # getRoots := proc( m::posint, countPolys::table, $) local coefs, counts, degq, f, facs, i, k, npolys, p, polys, rtsq, rtvec, x; description "Compute all roots of all polynomials in the table indices."; polys := [indices(countPolys)]; npolys := nops( polys ); # m*npolys roots (some of which will be duplicated) rtvec := Vector(m*npolys,0,datatype=complex[8]): i := 0; # counter for polynomial roots for coefs in polys do p := add( coefs[1][1+k]*x^k, k=0..m ); facs := factors(p); # [u, [...]] for f in facs[2] do try rtsq := fsolve( f[1], x, 'complex'); if degree(f[1],x)=1 then to f[2] do i := i+1; rtvec[i] := rtsq; end do; else for k to degree(f[1],x) do to f[2] do i := i+1; rtvec[i] := rtsq[k]; # multiple roots adjacent in list end do; end do; end if; catch: error("fsolve is having trouble with", f, rtsq ); end try; end do; end do: if i<>m*npolys then warning("Missed some roots somehow", i, m*npolys); end if; return rtvec; end proc: getRootsOld := proc( m::posint, countPolys::table, $) local coefs, counts, degq, i, k, npolys, offset, q, rtsq, rootsof, rtvec, tmp, x; description "Compute all roots of all polynomials in the table indices."; npolys := nops( [indices(countPolys)] ); # m*npolys roots (some of which will be duplicated) rtvec := Vector(m*npolys): rootsof := Vector(m); i := 0; # counter for polynomial roots degq := floor(m/2); offset := m mod 2; # even m means no zero root # odd m means always a zero root # In either case half the coefficients are zero # so we take advantage of this. for coefs,counts in eval(countPolys) do q := add( coefs[1+offset + 2*k]*x^k, k=0..degq ); rtsq := fsolve( factor(q), x, 'complex' ); if offset=1 then rootsof[1] := 0.; end if; for k to degq do tmp := sqrt(rtsq[k]); rootsof[offset+2*k-1] := tmp; rootsof[offset+2*k] := -tmp; od; for k to m do rtvec[i+k] := rootsof[k]; od; i := i+m; end do: return rtvec; end proc: # makePolyFun2 --- an alternative way to # get the polynomials # via the recurrence relation # # # The following artisanal, hand-crafted routine # to compute the characteristic polynomial from # the recurrence relation # P[n] = lambda*P[n-1] - U[n-1]^2*P[n-2] # runs only twice as slowly (!) as the brutal thing built # automatically above. # # The next step in speed improvement is to make # everything a float[8] and compile the code... # but that obviates any benefit of exact computation. # # Demonstrates in-place array work with the # hackware package; code generation; scoping; # sequence syntax with $ makePolyFun2 := proc( m::posint, $ ) if m=1 then return proc( U ) Vector(2,[0,1]) end proc elif m=2 then return proc( U ) Vector(3,[-U[1]^2,0,1]) end proc else return proc( U ) local deg, j, p1, p2, p3, tmp, v1, v2, v3; v1 := Vector(m+1, [0,1,0$m-1] ); v2 := Vector(m+1, [-U[1]^2, 0, 1, 0$m-2] ); v3 := Vector(m+1, 0 ); p1 := addressof(v1); p2 := addressof(v2); p3 := addressof(v3); for deg from 3 to m do # Build vector P = lambda*p2 - U[deg-1]^2*p1 tmp := -U[deg-1]^2; pointto(p3)[1] := tmp*pointto(p1)[1]; for j from 2 to deg-1 do pointto(p3)[j] := pointto(p2)[j-1] + tmp*pointto(p1)[j]; end do; pointto(p3)[deg] := 0; pointto(p3)[deg+1] := 1; # Now cycle the pointers tmp := p1; p1 := p2; p2 := p3; p3 := tmp; end do; return pointto(p2) end proc end if end proc: # getAllCharPolys --- get all characteristic polys # from the family. # # Demonstrates Iterators package. # keyword parameters # parameter defaults # Normalizer # (makes 0 expressions have the 0 representation) # # charPolyFun := makePolyFun(m); # countPolys := getAllCharPolys(m,P,charPolyFun); # or # (countPolys,matrixReps) := getAllCharPolys(m,P,charPolyFun,getReps=true); getAllCharPolys := proc( m::posint, P::list, F::procedure, {getReps::boolean:=false}, {getAllReps::boolean:=false}, $ ) local countPolys, k, Matrixrep, matrixentries, p, plist, u; description "Count characteristic polynomial frequencies of all complex symmetric tridiagonal Bohemians with zero diagonal"; # Generate all nops(P)^(m-1) possible vectors for the U values matrixentries := Iterator:-CartesianProduct(seq(P, k=1..m-1)): # Could sample randomly. Here we will do exhaustive computation countPolys := table(sparse): if getReps or getAllReps then Matrixrep := table(); # Disallow getReps if getAllReps is chosen if getAllReps and getReps then error("Choose only one of getReps or getAllReps.") end if end if; for u in matrixentries do p := F(u); #Charpoly of A # The purpose of the following is to give a unique key # into the table; p is a Vector and its name is not # a good key, though its values are (say, as a list) plist := [seq(Normalizer(p[k]),k=1..m+1)]; # Thx Juergen Gerhard for the "sparse" trick below countPolys[plist] := countPolys[plist]+1; if getReps and not assigned(Matrixrep[plist]) then # Return the parameters for a matrix representatio # of this particular characteristic polynomial Matrixrep[plist] := [seq(u)]; elif getAllReps then if not assigned(Matrixrep[plist]) then Matrixrep[plist] := [[seq(u)]]; else # Much more expensive than just getting one rep Matrixrep[plist] := [[seq(u)],op(Matrixrep[plist])]; end if; end if; end do: if getReps or getAllReps then return eval(countPolys), eval(Matrixrep) else return eval(countPolys) end if; end proc: # # Missing (left as exercises for the student) # # 1. getSomeEigenvalues # # sample the matrices using the RandomTools # package, as we did for getSomePolys above # but for matrices. # # 2. Count multiple roots # # Getting all (or some) characteristic polynomials # allows us to factor them, and to count how many # polynomials have multiple roots (and thus matrices # with multiple eigenvalues). # # 3. Count distinct roots in the family # Another thing that can be done is to compute a # "GCD-Free Basis", the roots of which are the # distinct eigenvalues of the family. Some matrices # in the family share eigenvalues, of course (this # is not the same as multiple eigenvalues) # # 4. A faster way to compute roots from polynomials # is to compute eigenvalues of the matrix representative. # This would give the benefits of the compression to # polynomials and keep the speed of computation. It # sacrifices the high-quality solution of fsolve, # but that is expensive and may not be necessary. # # 5. Circulant matrices with zero diagonal # # Such matrices also take m-1 parameters. # It should be simple to modify this package # to work with that family. The first pass # could ignore the connection with the DFT/FFT # but any serious exploration should use that. # (My student Cristian Ardelean has worked on # this, and in the symbolic case the results # are quite interesting). #