# Miroslav and the separating wand

In a huge graph, it is hard to determine a bottleneck by just looking at it, or doing any type of search. And by a bottleneck, I mean part of graph that is not very much connected to the other part. In different setting this could mean different things, for example, if your graph is a communication network, then a bottleneck mean the information doesn’t get to all of the nodes very easily, there will be traffic around the bottleneck.

One of the celebrated results in algebraic graph theory is by Miroslav Fiedler, a Czech mathematician [1, 2], about the algebraic connectivity of a graph. It simply says the second smallest eigenvalue of the Laplacian matrix of a graph determines how connected the graph is. In particular, if it is 0, the graph is disconnected. Furthermore, it shows the eigenvector corresponding to this eigenvalue distinguishes the vertices in the two connected components. In the case that the graph is connected, this eigenvector will find bottlenecks of the graph. Below I first generate a random graph that has a visible bottleneck, then I’ll use Fiedler’s result to find it computationally.

First, let’s generate a random graph with a bottleneck. In order to do so, I’ll generate two random graphs, then I look at their disjoint union and put few random edges between them.

def random_bottleneck_graph(n1,n2,p1,p2,p):
#n1 = number of vertices of first part
#n2 = number of vertices of second part
#p1 = probability of edges of first part
#p2 = probability of edges of second part
#p  = probability of edges between two parts

import random

#generate two parts randomly and put them together
H1 = graphs.RandomGNP(n1,p1) #generate a random GNP graph for first part
H2 = graphs.RandomGNP(n2,p2) #generate a random GNP graph for second part
G = H1.disjoint_union(H2) #form the disjoint union of the two parts

#add edges between two parts
for i in range(n1): #pick each vertex in H1
for j in range(n2): #pick each vertex in H2
a = random.random() #generate a random number between 0 and 1
if a > 1-p: #with probability p add an edge between the two vertices of H1 and H2
G.add_edge([(0,i),(1,j)])
return(G)


And here is a sample run:

#initialization
n1 = 40 #number of vertices of first part
n2 = 60 #number of vertices of second part
p1 = .7 #probability of edges of first part
p2 = .6 #probability of edges of second part
p = .01 #probability of edges between two parts

#run
G = random_bottleneck_graph(n1,n2,p1,p2,p)


Of course when I look at the graph, sage who has generated it knows that the graph has two main parts, and even the vertices are numbered in a way that anyone can recognize it. The left picture is generated by the above ‘show’ command, and the right picture is using the following code.

def show_random_bottleneck_graph(G,layout=None):
EdgeColors={'red':[], 'blue':[], 'black':[]}

for e in G.edges():
if not e[0][0] == e[1][0]:
EdgeColors['black'].append(e)
elif e[0][0] == 0:
EdgeColors['red'].append(e)
else:
EdgeColors['blue'].append(e)

G.show(vertex_labels=False,figsize=[10,10],layout=layout,edge_colors=EdgeColors)


So, in order to really get a random graph with a bottleneck, I take the above graph and try to loose the information that I had from it.

def relable_graph(G,part=None):
# input: a graph G
# output: if part is not given, it gives a random labeling of nodes of graph with integers
#         if part is given, then it partitions the vertices of G labeling each vertex with pairs (i,j) where i in the number of partition the vertex is in part, and j is an integer.

V = G.vertices()
E = G.edges()

H = Graph({})

map = {}
if part == None:
new_part = None
import random
random.shuffle(V)

for i in range(len(V)):
map[V[i]] = i
else:
new_part = [ [] for p in part]
count = 0
for i in range(len(part)):
for w in part[i]:
map[w] = (i,w)
new_part[i].append((i,w))
count += 1

for v,w,l in E:
H.add_edge([map[v],map[w],l])

return(H,new_part)

###############################
#run
H = relable_graph(G)[0]
H.show(layout='circular',vertex_labels=False, figsize=[10,10])

Now if I look at the graph I won’t be able to recognize the bottleneck:

Now that I truly have random graph, which I know has a visible bottleneck, let’s see if we can find it. Below I use ‘networkx’ to compute the Fiedler vector, an eigenvector corresponding to the second smallest eigenvalue of the Laplacian matrix. Then I’ll partition the vertices according to the sign of the corresponding entry on this vector:

def partition_with_fiedler(G):
# input: a graph G
# algorithm: using networkx finds the fiedler vector of G
# a partitioning of vertices into two pieces according to positive and negative entries of the Fiedler vector

import networkx
H = G.networkx_graph()
f = networkx.fiedler_vector(H, method='lanczos')

V = H.nodes()
P = [[],[]]
for i in range(len(V)):
if f[i] >= 0:
P[0].append(V[i])
else:
P[1].append(V[i])

return(P)


Let’s see how it is partitioned:

P = partition_with_fiedler(H)
H.show(partition=P, layout='circular', vertex_labels=True, figsize=[10,10])

Well, not much can be seen yet. Now let me put the red vertices in one side and the blue vertices in the other side:

def position_bottleneck(G):
P = partition_with_fiedler(G)

V = G.vertices()
new_pos = {}

for v in V:
if v in P[0]:
new_pos[v] = [6*(1+random.random()),6*(2*random.random()-1)]
else:
new_pos[v] = [-6*(1+random.random()),6*(2*random.random()-1)]

return(new_pos)
pos = position_bottleneck(H)
H.show(partition=P, pos=pos, vertex_labels=True, figsize=[10,10])

Now it’s easy to see the bottleneck of the graph! We can even recover the original picture as follows:

(K,newP) = relable_graph(H,part=P)
K.show(partition=newP, layout='circular', vertex_labels=False, figsize=[10,10])

And here is the representation of it using the colored edges:

show_random_bottleneck_graph(K,layout='circular') #suggested layouts are "circular" and "spring"

UPDATE:

Since I posted this, I’ve been thinking about multiway Fiedler clustering. Here is what I got so far. First I define a list N of number of vertices in each part, and a matrix A which is basically a probability matrix where entry (i,j) is the probability of an edge between part i and part j, hence it is symmetric. If you choose diagonal entries close to 1 and off-diagonal entries close to zero, then you’ll get a graph with multiple communities. I’ll shuffle the vertices around after creating the multi-community graph.

def multi_bottleneck_graph(N,A):
#N = list of size n of number of vertices in each part
#A = an nxn matrix where entry (i,j) shows the probability of edges between part i and part j of vertices, and n is the number of parts.

import random

n = len(N)
G = Graph({})

#add vertices for each part
for i in range(n):
for j in range(N[i]):
G.add_vertex((i,j))

V = G.vertices()

#add edges
for v in V:
for w in V:
if not v == w:
a = random.random() #generate a random number between 0 and 1
p = A[v[0],w[0]]
if a > 1-p: #with probability p add an edge between the two vertices of H1 and H2
G.add_edge([v,w])

return(G)

def relable_graph(G,part=None):
# input: a graph G
# output: if part is not given, it gives a random labeling of nodes of graph with integers
#         if part is given, then it partitions the vertices of G labeling each vertex with pairs (i,j) where i in the number of partition the vertex is in part, and j is an integer.

V = G.vertices()
E = G.edges()

H = Graph({})

map = {}
if part == None:
new_part = None
import random
random.shuffle(V)

for i in range(len(V)):
map[V[i]] = i
else:
new_part = [ [] for p in part]
count = 0
for i in range(len(part)):
for w in part[i]:
map[w] = (i,w)
new_part[i].append((i,w))
count += 1

for v,w,l in E:
H.add_edge([map[v],map[w],l])

return(H,new_part)

def shuffled_multi_bottleneck_graph(N,A):
G = multi_bottleneck_graph(N,A)
H = relable_graph(G)[0]
return(H)

# run
N = [10,15,20]
A = matrix([[.9,.01,.01],[.01,.9,.01],[.01,.01,.9]])
G = shuffled_multi_bottleneck_graph(N,A)
G.show(layout='circular',vertex_labels=False, vertex_size=30, figsize=[10,10])


Then I break the graph into two pieces using Fiedler vector, compare the two pieces according their algebraic connectivities (second smallest eigenvalue of the Laplacian) and break the smaller one into two pieces and repeat this process, until I get n communities.

def partition_with_fiedler(G):
# input: a graph G
# algorithm: using networkx finds the fiedler vector of G
# a partitioning of vertices into two pieces according to positive and negative entries of the Fiedler vector

import networkx
H = G.networkx_graph()
f = networkx.fiedler_vector(H, method='lanczos')

V = H.nodes()
P = [[],[]]
for i in range(len(V)):
if f[i] >= 0:
P[0].append(V[i])
else:
P[1].append(V[i])

return(P)

def multi_fiedler(G,n=2):
if n == 2:
P = partition_with_fiedler(G)
return(P)

else:
import networkx

P = partition_with_fiedler(G)
H = [ G.subgraph(vertices=P[i]) for i in range(2) ]
C = [ (networkx.algebraic_connectivity(H[i].networkx_graph(), method='lanczos')/H[i].order(),i) for i in range(2) ]
l = min(C)[1]

while n > 2:
n = n-1

temP = partition_with_fiedler(H[l])
nul = P.pop(l)
P.append(temP[0])
P.append(temP[1])

temH = [ G.subgraph(vertices=temP[i]) for i in range(2) ]
nul = H.pop(l)
H.append(temH[0])
H.append(temH[1])

temC = [ (networkx.algebraic_connectivity(temH[i].networkx_graph(), method='lanczos')/temH[i].order(),i) for i in range(2) ]
nul = C.pop(l)
C.append(temC[0])
C.append(temC[1])
l = min(C)[1]

return(P)
# run
G.show(partition=K, layout='circular', figsize=[10,10])


Now I can separate same colour vertices to see the communities in the graph:

def show_communities(G,n):
#G: a graph
#n: number of communities > 1
if n == 1:
P = [G.vertices()]
elif n > 1:
P = multi_fiedler(G,n)
else:
raise ValueError('n > 0 is the number of communities.')
L,newP = relable_graph(G,part=P)
L.show(partition=newP, layout='circular', vertex_labels=False, figsize=[10,10])

# run
show_communities(G,3)


There is only one important thing in using networkx. It seems like the default method (tracemin) for algebraic_connectivity and for fiedler_vector has a bug. So, I used the method=’lanczos’ option there.

You might be interested to see what happens if you ask for more or less than 3 communities. Here are a couple of them:

The next question is how to decide about number of communities.

1. M. Fiedler, Algebraic connectivity of graphs. Czechoslovak Math. J. 23(98):298–305 (1973).
2. M. Fiedler, A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory, Czech. Math. J. 25:619–637, (1975).
Advertisements

# Step-by-step reduction

One of the things that I always tell my students is to check their solutions when they are done solving a problem. That by itself can mean several things, depending on what the problems is. Of course after solving a system of linear equations, one would plug in the solutions into the original system of equations to see if they are satisfied. The harder part is to come up with a way to check if you’ve done the row-reduction correctly or not. One can easily use Sage to see if the reduced row echelon form is computed correctly. If A is your matrix, just input A.rref() to find the reduced row echelon form of it.

But still sometimes in the case that the answer is wrong, one might get frustrated by not trying to figure out in which step they’ve made an error, or sometimes they might just get stuck. Considering that row-reduction is a time consuming process, it is plausible to assume that some one can easily give up after a few tries. I have found this nice code in Sage that shows step-by-step Gauss-Jordan elimination process, and actually tells you what to do in each step, from here. Then I did a little bit of customization on it. Here is the final version:

# Naive Gaussian reduction
def gauss_method(MATRIX,rescale_leading_entry='Last'):
"""Describe the reduction to echelon form of the given matrix of rationals.

MATRIX  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])
rescale_leading_entry='First' make the leading entries to 1's while doing Gaussisan ellimination
rescale_leading_entry='Last' (Default) make the leading entries to 1's while reducing

Returns: reduced form.  Side effect: prints steps of reduction.

"""
M = copy(MATRIX)
num_rows=M.nrows()
num_cols=M.ncols()
show(M.apply_map(lambda t:t.full_simplify()))

col = 0   # all cols before this are already done
for row in range(0,num_rows):
# ?Need to swap in a nonzero entry from below
while (col < num_cols
and M[row][col] == 0):
for i in M.nonzero_positions_in_column(col):
if i > row:
print " swap row",row+1,"with row",i+1
M.swap_rows(row,i)
show(M.apply_map(lambda t:t.full_simplify()))
break
else:
col += 1

if col >= num_cols:
break

# Now guaranteed M[row][col] != 0
if (rescale_leading_entry == 'First'
and M[row][col] != 1):
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
show(M.apply_map(lambda t:t.full_simplify()))

for changed_row in range(row+1,num_rows):
if M[changed_row][col] != 0:
factor = -1*M[changed_row][col]/M[row][col]
print " take",factor,"times row",row+1,"plus row",changed_row+1
M.add_multiple_of_row(changed_row,row,factor)
show(M.apply_map(lambda t:t.full_simplify()))
col +=1

print "Above is a row echelon form, let's keep cruising to get the reduced row echelon form:\n"

for i in range(num_rows):
row = num_rows-i-1
if M[row] != 0:
for col in range(num_cols):
if M[row,col] != 0:
if M[row,col] != 1:
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
show(M.apply_map(lambda t:t.full_simplify()))
break

for changed_row in range(row):
factor = -1 * M[row-1-changed_row,col]
if factor != 0:
print " take", factor,"times row", row+1, "plus row", row-1-changed_row+1
M.add_multiple_of_row(row-1-changed_row,row,factor)
show(M.apply_map(lambda t:t.full_simplify()))
return(M.apply_map(lambda t:t.full_simplify()))


And here is a sample run:

M = matrix(SR,[[3,-1,4,6],[0,1,8,0],[-2,1,0,-4]])
gauss_method(M,rescale_leading_entry='First')


And the output looks like this:

Click here to run it in the sage cell server. Note that the code is implemented to work over rational field since it is being used for pedagogical reasons.

In the comments Jason Grout has added an interactive implementation of the first part of the code:

# Naive Gaussian reduction
@interact
def gauss_method(M=random_matrix(QQ,4,algorithm='echelonizable',rank=3),rescale_leading_entry=False):
"""Describe the reduction to echelon form of the given matrix of rationals.

M  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])
rescale_leading_entry=False  boolean  make the leading entries to 1's

Returns: None.  Side effect: M is reduced.  Note: this is echelon form,
not reduced echelon form; this routine does not end the same way as does
M.echelon_form().

"""
num_rows=M.nrows()
num_cols=M.ncols()
print M

col = 0   # all cols before this are already done
for row in range(0,num_rows):
# ?Need to swap in a nonzero entry from below
while (col < num_cols                and M[row][col] == 0):              for i in M.nonzero_positions_in_column(col):                 if i > row:
print " swap row",row+1,"with row",i+1
M.swap_rows(row,i)
print M
break
else:
col += 1

if col >= num_cols:
break

# Now guaranteed M[row][col] != 0
if (rescale_leading_entry
and M[row][col] != 1):
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
print M
change_flag=False
for changed_row in range(row+1,num_rows):
if M[changed_row][col] != 0:
change_flag=True
factor=-1*M[changed_row][col]/M[row][col]
print " take",factor,"times row",row+1,"plus row",changed_row+1
M.add_multiple_of_row(changed_row,row,factor)
if change_flag:
print M
col +=1


One of the uses of this can be to find the inverse of a matrix step by step. For example:

var('a b c d')
A = matrix([[a,b],[c,d]])
R = gauss_method(A.augment(identity_matrix(A.ncols()),subdivide=True))

will give you: Or you can run it for a generic matrix of size 3 (for some reason it didn’t factor i, so I used letter k instead:)

var('a b c d e f g h i j k l m n o p q r s t u v w x y z')
A = matrix([[a,b,c],[d,e,f],[g,h,k]])
R = gauss_method(A.augment(identity_matrix(A.ncols()),subdivide=True))

The steps are too long, so I’m not going to include a snapshot, if you are interested look at the steps here. But you can check out only the final result by

view(R.subdivision(0,1))

which will give you And to check if the denominators and the determinant of the matrix have any relations with each other you can multiply by the determinant and simplify:

view((det(A)*R.subdivision(0,1)).apply_map(lambda t:t.full_simplify()))

to get

Do you think this can help students?

# Generating not-so-much-random matrices

I’ve been using sage for a while now. I’ve also been teaching linear algebra for a few years. One of the problems in teaching linear algebra is coming up with a handful of examples of matrices that can nicely be row-reduced. One way to do this is to start with your “nice” row-reduced matrix and take a few (or many) “nice” elementary matrices and there is your example. So, one might want to code this in sage. I just recently found out that actually sage has such a function built in it. For example,

random_matrix(ZZ,4,5,algorithm='echelonizable',rank=3, upper_bound=7)


will create a random $4 \times 5$ matrix with rank $3$, integer entries (that’s the ZZ), and entries no bigger than $7$, which can “nicely” be turned into rref. So, here it is in sage cell server. I think this is a good thing to share with students so they can practice as many times as they want.

For the same reasons when generating a random matrix to find its inverse, we usually want a matrix with small-ish integer entries such that the inverse also has relatively small integer entries. It is easy to see that for a matrix $A$ with integer entries to have an inverse with integer entries, it is necessary that its determinant is $\pm 1$. It turns out that it is the sufficient condition too, surprisingly, thanks to Cramer’s rule! Though it might lose some pedagogical purposes, like relating the inverse to the determinant, but here is a code that generates random matrices with integer entries where their inverses also have integer entries:

random_matrix(ZZ,3,algorithm='unimodular',upper_bound=6)

To play around with it in sage cell server click here.

The sage documentation says:

Warning: Matrices generated are not uniformly distributed. For
unimodular matrices over finite field this function does not even
generate all of them: for example "Matrix.random(GF(3), 2,
algorithm='unimodular')" never generates "[[2,0],[0,2]]". This
function is made for teaching purposes.

And

Random matrices in echelon form.  The "algorithm='echelon_form'"
keyword, along with a requested number of non-zero rows
("num_pivots") will return a random matrix in echelon form.  When
the base ring is "QQ" the result has integer entries.  Other exact
rings may be also specified.

sage: A=random_matrix(QQ, 4, 8, algorithm='echelon_form', num_pivots=3); A # random
[ 1 -5  0 -2  0  1  1 -2]
[ 0  0  1 -5  0 -3 -1  0]
[ 0  0  0  0  1  2 -2  1]
[ 0  0  0  0  0  0  0  0]
sage: A.base_ring()
Rational Field
sage: (A.nrows(), A.ncols())
(4, 8)
sage: A in sage.matrix.matrix_space.MatrixSpace(ZZ, 4, 8)
True
sage: A.rank()
3
sage: A==A.rref()
True

For more, see the documentation of the "random_rref_matrix()"
function.  In the notebook or at the Sage command-line, first
execute the following to make this further documentation available:

from sage.matrix.constructor import random_rref_matrix

Random matrices with predictable echelon forms.  The
"algorithm='echelonizable'" keyword, along with a requested rank
("rank") and optional size control ("upper_bound") will return a
random matrix in echelon form.  When the base ring is "ZZ" or "QQ"
the result has integer entries, whose magnitudes can be limited by
the value of "upper_bound", and the echelon form of the matrix also
has integer entries.  Other exact rings may be also specified, but
there is no notion of controlling the size.  Square matrices of
full rank generated by this function always have determinant one,
and can be constructed with the "unimodular" keyword.

sage: A=random_matrix(QQ, 4, 8, algorithm='echelonizable', rank=3, upper_bound=60); A # random
sage: A.base_ring()
Rational Field
sage: (A.nrows(), A.ncols())
(4, 8)
sage: A in sage.matrix.matrix_space.MatrixSpace(ZZ, 4, 8)
True
sage: A.rank()
3
sage: all([abs(x)<60 for x in A.list()])
True
sage: A.rref() in sage.matrix.matrix_space.MatrixSpace(ZZ, 4, 8)
True

For more, see the documentation of the
"random_echelonizable_matrix()" function.  In the notebook or at
the Sage command-line, first execute the following to make this
further documentation available:

from sage.matrix.constructor import random_echelonizable_matrix
The "x" and "y" keywords can be used to distribute entries
uniformly. When both are used "x" is the minimum and "y" is one
greater than the maximum.

sage: random_matrix(ZZ, 4, 8, x=70, y=100)
[81 82 70 81 78 71 79 94]
[80 98 89 87 91 94 94 77]
[86 89 85 92 95 94 72 89]
[78 80 89 82 94 72 90 92]
If only "x" is given, then it is used as the upper bound of a range
starting at 0.

sage: random_matrix(ZZ, 5, 5, x=25)
[20 16  8  3  8]
[ 8  2  2 14  5]
[18 18 10 20 11]
[19 16 17 15  7]
[ 0 24  3 17 24]
To control the number of nonzero entries, use the "density" keyword
at a value strictly below the default of 1.0.  The "density"
keyword is used to compute the number of entries that will be
nonzero, but the same entry may be selected more than once.  So the
value provided will be an upper bound for the density of the
created matrix.  Note that for a square matrix it is only necessary
to set a single dimension.

sage: random_matrix(ZZ, 5, x=-10, y=10, density=0.75)
[-6  1  0  0  0]
[ 9  0  0  4  1]
[-6  0  0 -8  0]
[ 0  4  0  6  0]
[ 1 -9  0  0 -8]
For a matrix with low density it may be advisable to insist on a
sparse representation, as this representation is not selected
automatically.

sage: random_matrix(ZZ, 5, 5, density=0.3, sparse=True)
[ 4  0  0  0 -1]
[ 0  0  0  0 -7]
[ 0  0  2  0  0]
[ 0  0  1  0 -4]
[ 0  0  0  0  0]
Random rational matrices:
sage: random_matrix(QQ, 2, 8, num_bound=20, den_bound=4)
[ -1/2 6 13 -12 -2/3 -1/4 5 5]
[ -9/2 5/3 19 15/2 19/2 20/3 -13/4 0]
Random matrices over other rings.  Several classes of matrices have
specialized "randomize()" methods.  You can locate these with the
Sage command:

search_def('randomize')

# Gershgorin Disks

Gershgorin circle theorem roughly states that the eigenvalues of an $n \times n$ matrix lie inside $n$ circles where $i$-th circle is centred at $A_{i,i}$ and its radius is the sum of the absolute values of the off-diagonal entries of the $i$-th row of $A$. Applying this to $A^{\top}$ implies that the radius of this circles shall be the smaller of the sums for rows and columns. The following Sage code draws the circles for a given matrix.

def Gershgorin(A,evals=False):
# A is a square matrix
# Output is the Gershgorin disks of A

from sage.plot.circle import Circle
from sage.plot.colors import rainbow

n = A.ncols()
Colors = rainbow(n, 'rgbtuple')
E = point([])

B = A.transpose()
R = [(A[i-1,i-1], min(sum( [abs(a) for a in (A[i-1:i]).list()] )-abs(A[i-1,i-1]),sum( [abs(a) for a in (B[i-1:i]).list()] )-abs(B[i-1,i-1]))) for i in range(1,n+1)]
C = [ circle((real(R[i][0]),imaginary(R[i][0])), R[i][1], color="black") for i in range(n)]
if evals == True:
E = point(A.eigenvalues(), color="black", size=30)
CF = [ circle((real(R[i][0]),imaginary(R[i][0])), R[i][1], rgbcolor=Colors[i], fill=True, alpha=1/n) for i in range(n)]

(sum(C)+sum(CF)+E).show()


And here is a sample output:

A =  random_matrix(CDF,4)
Gershgorin(A)

Gershgorin disks for a randomly generated $4\times 4$ complex matrix

If you want to see the actual eigenvalues, call the function as below:

Gershgorin(A,evals=True)


# Linear approximation and Taylor polynomials

For a given function $f(x)$ around a point $a$, if you zoom in close enough on the graph of the function, you can’t tell the difference between the curve and the tangent line to the graph of the function at point $a$. This is the basic idea for linear approximations. To explore that statement here I wrote a little sage interactive code that will let you play around with it:

@interact
def _(
f = ('$f(x)=$', input_box(x^3 + 6*x^2 - 15*x + 7,width=20)),
a = ('around the point $x = a$', input_box(0,width=10)),
zoom = ('zoom', slider(.1,12,default=1))):

g(x) = f(a) + diff(f,x)(a)*(x-a)
P = plot(f,(x,a-10*2^(-zoom), a+10*2^(-zoom)),  color="blue", axes=True, thickness=4) # for horizontal tangent lines add 'aspect_ratio=1' to see something meaningful, otherwise sage will just change the ration so that the graph is not just a horizontal line.
Q = plot(g,(x,a-10*2^(-zoom), a+10*2^(-zoom)),  color="red", axes=True, thickness=2)
show(P + Q)


Click here to play with it in the sage cell server. Here is an output:

The linear approximations have errors though. One way to get better approximations is by taking into account the higher order derivatives of the function at the given point. The following code will calculate and draw the n-th degree Taylor polynomial of the function so that you can compare them:

@interact
def _(
f = ('$f(x)=$', input_box(sin(x),width=20)),
n = ('Degree of the Taylor polynomial, $n=$', slider([0..9],default=3)),
a = ('around the point $x = a$', input_box(0,width=10)),
xzoom = ('range of $x$', range_slider(-12,12,default=(-12,12))),
yzoom = ('range of $y$', range_slider(-20,20,default=(-4,4)))):

g = f.taylor(x, a, n)
P = plot(f,(x,xzoom[0], xzoom[1]), ymin = yzoom[0], ymax = yzoom[1], color="blue", axes=True, aspect_ratio=1, thickness=3)
Q = plot(g,(x,xzoom[0], xzoom[1]), ymin = yzoom[0], ymax = yzoom[1], color="red", axes=True, aspect_ratio=1)
R = text(r"$%s \approx %s$" %(f,g),(0,yzoom[1]+1),rgbcolor=(1,0,0))
show(P + Q + R)


Click here to see it in the sage cell server. Here is an output:

The only problem with that code is that it doesn’t render the latex formula correctly. That is, the two digit exponents are shown wrong. So, for now I’ve limited $n$ to 9. If you have any suggestions on how to fix it please see the discussion here.

# What to do with a long list?

I mean a list in SAGE. The answer is to

save(mylist, 'mylist')


and then later

mynewlist = load('mylist')


This will save your list into a file called mylist.sobj and then loads it. The story is I calculated a list through a long calculation to do some other calculations with it. So I saved the list in a .txt file and in the file it looks like this:

[[[-1/10, 1.02421477975960],
[-99/1000, 1.02369404664481],
[-49/500, 1.02317986236459],
[-97/1000, 1.02267219319285],
...


But then that I needed to read it, I had to open the text file and copy and paste it into SAGE. If the list is short, I think it’s an easy way to do it, but if the list is long (a few gigabytes maybe?) you don’t want to open that text file, forget about copying it! So, one way to read it was to strip all the \n‘s and get rid of the brackets and all the other things in a nice way and then read it as a csv file, which I could save it as a csv in the first place with:

import csv
with open('mylist.csv', 'w') as f1:
writefile = csv.writer(f1)
writefile.writerows(mylist)


And then read it with

import csv
with open('mylist.csv','rU') as f1:
mynewlist =load('myfile')


The problem is I still had to go through the list and change strings to numbers. Something I really don’t want to do. So, the save and load file is the best thing.

save(mylist, 'mylist')
mynewlist = load('mylist')


Resource: AskSageMath community.

# A full matrix with given eigenvalues

## Induction: comprehension and precision

In Theorem 3.1 of “The combinatorial inverse eigenvalue problems: complete graphs and small graphs with strict inequality, W. Barrett, A. Lazenby, N. Malloy, C. Nelson, W. Sexton, Electronic Journal of Linear Algebra, 656-672 (26) 2013” the authors show in an inductive process that $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n$ are eigenvalues of a real symmetric matrix $A$ with no zero off-diagonal entries if and only if $\lambda_1 > \lambda_n$. Even though they don’t mention it, the proof constructs a matrix whose diagonal entries are also nonzero, hence the result matrix is a full matrix. The induction is straightforward and uses a clever starting case, and then extends it using orthogonal similarities. A cool idea in the construction is that the diagonal entries of a real symmetric matrix are bounded by its smallest and largest eigenvalues (Corollary 2.6). Furthermore, if the matrix is irreducible, then the diagonal entries of the matrix are strictly bounded by the smallest and largest eigenvalues of the matrix (Lemma 2.7).

The following is a SAGE code to produce a full matrix for the given set of eigenvalues following the steps in the proof of Theorem 3.1 of the mentioned paper.

###########################################################################
# this is an auxiliary function which computes the direct sum of two given matrices
###########################################################################
def direct_sum(A,B):
ar = A.nrows()
ac = A.ncols()
br = B.nrows()
bc = B.ncols()
M = matrix.block([[A, zero_matrix(ar,bc)],[zero_matrix(br,ac) ,B]])
return(M)
############################################################
# this is the main function which construct the matrix
############################################################
def nonwhere_zero_matrix_with_spectrum(L):
# L is a list of desired eigenvalues with at least two of them distinct
L.sort()
# This checks to see if at least two of them are distinct
if L[0] == L[-1]:
print('Not all eigenvalues can be equal. Such matrix does not exist')
return()
# We do an induction. The base of the iduction is a matrix with two eigenvalues that are the two smallest distinct numbers in the list L
# Start with the first two numbers in the list:
a = L[0]
b = L[1]
count = 0
# Search until you find the second smallest number in the list that is not equal to the first number in the list
while true:
if L[count] > L[0]:
b = L[count]
count += 1 # this becomes the size of B
break
count += 1
# Construct the matrix for the base:
B = (b-a) / (count) * ones_matrix(count) + a * identity_matrix(count)
# Here c and d will give me an orthogonal matrix Q. I'm letting them to be 1, but one could change them as desired, as long as they are nonzero we'll get what we need
c = 1
d = 1
# This is the inductive step. Each time add a new row and column with a diagonal entry equal to a the next eigenvalue in the list, and then multiply by the orthogonal matrix to make that row and column nowherezero
while count < len(L):
B = direct_sum(matrix([L[count]]), B)
Q = direct_sum(1/sqrt(c^2 + d^2) * matrix([[c,-d],[d,c]]),identity_matrix(count-1))
B = (Q.transpose()) * B * Q
count += 1
# The desired matrix is ready
return(B)
#########################################################################
# Usage
#########################################################################
A = nonwhere_zero_matrix_with_spectrum([1,1,1,3,3,6])
A.n(16)

To run the code in sage cell server click here, and to follow the code on github click here.

## Jacobian method: limitations vs. liberty

On the other hand, if you want to have control over the off-diagonal entries, and if the eigenvalues of the matrix are distinct, then you can use Theorem 4.2 from our paper “Construction of matrices with a given graph and prescribed interlaced spectral data, Keivan Hassani Monfared and Bryan L. Shader, Linear Algebra and Its Applications; 4348–4358 (438) 2013“. An implementation of the process in the proof of the theorem in SAGE:

#########################################################################
# Build variables, and the matrix corresponding to it
#########################################################################

def build_variables(n):
names = [ [[] for i in range(n)] for j in range(n) ]

for i in range(n):
for j in range(i+1):
names[i][j] = (SR('x_' + str(j) + '_' + str(i)))
names[j][i] = (SR('x_' + str(j) + '_' + str(i)))

return(names)

#########################################################################
# Define the function f that maps a matrix to the coefficients of its characteristic polynomial
#########################################################################
def CharPoly(Mat):
X = matrix(Mat)
n = X.ncols()
C_X = X.characteristic_polynomial()
Y = []
for i in range(n):
Y.append(C_X[i])
return(Y)

############################################################################
# This solves that lambda SIEP
############################################################################
def lambda_siep(G,L,iter=100,epsilon = .1):
# G is any graph on n vertices
# L is the list of n desired distinct eigenvalues
# m is the number of itterations of the Newton's method
# epsilon: the off-diagonal entries will be equal to epsilon

n = G.order()
my_variables = build_variables(n)

R = PolynomialRing(CC,[my_variables[i][j] for i in range(n) for j in range(n)])
R.gens()
R.inject_variables()

X = [ [ R.gens()[n*i+j] for i in range(n) ] for j in range(n) ]

Y = matrix(CharPoly(X)) - matrix(CharPoly(diagonal_matrix(L)))

J = matrix(R,n)
for i in range(n):
for j in range(n):
J[i,j] = derivative(Y[0][i],my_variables[j][j])

B = diagonal_matrix(L) + epsilon * G.adjacency_matrix()

count = 0
while count < iter:

T = [ B[i,j] for i in range(n) for j in range(n)]

C = (J(T)).solve_right(Y(T).transpose())
LC = list(C)

B = B - diagonal_matrix([LC[i][0] for i in range(n)])

count = count + 1

return(B)

##########################################################################
# This shows the output matrix, its eigenvalues and the eigenvlaues of A(i), and its graph
##########################################################################

def check_output_lambda_siep(A,precision=8):
# A is a matrix which is the output of lambda_siep()
# i is the one that also is entered in lambda_siep()
# precision is an integer that shows how many digits do I want to be printed at the end, and I set the default to be 8

eigA = A.eigenvalues()
EigA = []
for e in eigA:
EigA = EigA + [e.n(precision)]
print('A is:')
print(A.n(precision))
print(' ')
print('Eigenvalues of A are: %s') %(EigA)

AdjA = matrix(A.ncols())
for i in range(A.ncols()):
for j in range(A.ncols()):
if i != j:
if A[i,j] != 0:
AdjA[i,j] = 1
FinalGraph = Graph(AdjA)

print(' ' )
print('And the graph of A is:')
FinalGraph.show()

# Usage:
G = graphs.CompleteGraph(5)
L = [1,2,3,4,5]

A = lambda_siep(G, L, iter=1000, epsilon=.1)
check_output_lambda_siep(A,precision=32)


The proof relies on the Implicit Function Theorem, and that if the eigenvalues that you start with are away from zero, the diagonals of the matrix stay away from zero after adjustment. This way, you can set the off-diagonal entries of the matrix to be any small numbers. In particular you can set them all be equal to each other (that is what the above code does for simplicity), or choose them all to be positive numbers.

To run the code in sage cell server click here, and to follow it on github click here.

Assuming the eigenvalues are all positive, the induction method gets really close to producing an entrywise positive matrix, but the rotations actually ensure that some of the entries become negative! On the other hand the Jacobian method would result in an entrywise positive matrix if the perturbations of the off-diagonal entries are small enough that the adjustments of the diagonal entries wont make them zero or negative. To be fair, even if you start with a zero eigenvalue (you cannot have more than one, because the Jacobian method requires the matrix to have distinct eigenvalues), since the diagonal entries are strictly bounded by the smallest and largest eigenvalues, we have the following theorem:

Theorem. Given real numbers $\lambda_1 > \lambda_2 > \ldots > \lambda_n \geq 0$, there is an entrywise positive matrix with the given eigenvalues.