One goal of this section is to introduce you to the joys of numerical linear algebra with sage.

Suppose $A$ is the adjacency matrix of $X$. We define $H(t)=H_X(t)$ by \[ H =\exp(iAt) \] We are concerned with the absolute values of the entries $H(t)$, and in particular want to know when we can have \[ |H(\tau)_{u,v}| =1 \] for some time $\tau$ and some vertices $u$ and $v$. Since $A$ is symmetric it has a spectral decomposition \[ A =\sum_\th \th E_\th \] and consequently \[ H(t) = \sum_\th \exp(i\th t)E_\th \]

To use this formula we need an orthogonal basis for each eigenspace. The difficulty if that if we work in floating point, then we have to decide which eigenvalues are actually equal before we determine an eigenspace.

We are going to use scipy and numpy. It seems that numpy is basically an array package, which provides such things as $k$-dimensional arrays. (Computer scientists seem to think these are useful, I have no idea why.) Scipy builds on numpy, providing access to standard linear algebra routines (eigenthings, svd, etc.) There is on-line documentation for scipy and numpy, but it is not very good. (For example it illustrates how solve linear equations by inverting the matrix of coefficients.) The syntax for scipy and numpy is not consistent with sage, or with python.

Anyway, here is the opening incantation:

import scipy
from scipy import linalg as lin

We will go through some simple computations, so you can get a feel for things.

P = graphs.PetersenGraph()
AP = P.am()

Now with

la, v = lin.eigh(AP)
la

array([-2., -2., -2., -2.,  1.,  1.,  1.,  1.,  1.,  3.])

v

array([[ 0.        ,  0.        ,  0.        ,  0.63245553,  0.        ,
         0.        ,  0.        ,  0.        ,  0.70710678, -0.31622777],
       [-0.18497308, -0.36304065,  0.23708367, -0.42163702, -0.1783604 ,
         0.38972155,  0.08372533, -0.50372532,  0.23570226, -0.31622777],
       [-0.00842549,  0.62345912,  0.01079912,  0.10540926, -0.06962648,
         0.15213545,  0.61464175, -0.19663905, -0.23570226, -0.31622777],
       [-0.22096121, -0.4419388 , -0.38046699,  0.10540926, -0.1672403 ,
         0.26202423,  0.26545821,  0.52664064, -0.23570226, -0.31622777],
       [ 0.41435978,  0.18152033,  0.1325842 , -0.42163702,  0.41994923,
         0.22885117, -0.04186266,  0.46255999,  0.23570226, -0.31622777],
       [-0.2293867 ,  0.18152033, -0.36966787, -0.42163702, -0.24158883,
        -0.61857272, -0.04186266,  0.04116533,  0.23570226, -0.31622777],
       [ 0.37837165,  0.10262218, -0.48496647,  0.10540926, -0.10873392,
         0.2375861 , -0.53091642, -0.30708626, -0.23570226, -0.31622777],
       [ 0.42278527, -0.4419388 ,  0.12178508,  0.10540926,  0.27597422,
        -0.49961034,  0.26545821, -0.21955438, -0.23570226, -0.31622777],
       [ 0.03598813,  0.07889815,  0.61755067,  0.10540926, -0.51756305,
        -0.11896239, -0.30732088,  0.2607197 , -0.23570226, -0.31622777],
       [-0.60775834,  0.07889815,  0.11529859,  0.10540926,  0.58718953,
        -0.03317306, -0.30732088, -0.06408065, -0.23570226, -0.31622777]])

we get the eigenvalues of $AP$ (in $\la$) and the eigenvectors (in \(v\)). The routine linalg.eigh takes a Hermitian matrix and returns real eigenvalues. Eigenvectors corresponding to distinct eigenvalues are orthogonal. Note that the eigenvectors are returned by linalg.eigh as numpy arrays, and there is a argument for wrapping it in a routine that immediately transforms them into sage vectors.

We intend to make use of the fact that if the vector $\seq z1m$ are an orthogonal basis for a subspace then the matrix representing projection onto the subspace is \[ z_1z_1^T+\cdots+z_mz_m^T \] The problem is that our eigenvalues have been computed in floating point, and so we need a routine with arguments $\la$ and $\eps$ which divides the eigenvalues into equivalence classes, where two eigenvalues are equivalent if their difference is less than $\eps$. The output of the routine is a dictionary keyed on distinct eigenvalues and the entry with given eigenvalue is the set of indices of the eigenvalues equivalent to the key.

Given this it is straightforward to write the code that takes this dictionary and returns a list of pairs consisting of an eigenvalue and its associated projection.