Eigenvector for Circulant Matrices

Define the \(n\times n\) matrix \(A\) by

\[ A=\begin{bmatrix} -2&1&0&0&\cdots&0&0&1\\ 1&-2&1&0&\cdots&0&0&0\\ 0&1&-2&1&\cdots&0&0&0\\ 0&0&1&2&\cdots&0&0&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\\ 0&0&0&0&\cdots&-2&1&0\\ 0&0&0&0&\cdots&1&-2&1\\ 1&0&0&0&\cdots&0&1&-2 \end{bmatrix} \] In this project you will explore (and explain) the behavior of the iterative map

\[\begin{equation} x^{(k+1)}=\left(\mathcal{I}+\eta A\right)x^{(k)} \end{equation}\]

for various initial vectors \(x^{(0)}\) and for \(\eta=0.1\)

a)

Let \(n=32\). Create the \(n\) dimensional initial vector \(x^{(0)}=[1,0,0,\cdots,0]^t\) and plot \(x^{(0)},x^{(1)},x^{(2)},\cdots\) Repeat this process for a selection of unitial vector until you can qualitatively describe the behavior of the map \(x\mapsto\left(\mathcal{I}+\eta A\right)x\)

• Let’s perform the mapping to the canonical vector \(e_1=[1,0,0,\cdots,0]^t\)

X<-rep(0,32)
X[1]<-1

#X is our initial vector (1,0,...,0) of n=32 entries


Prod<-function(k,n=32,mu=0.2,x=X){
  #Performs the map x\mapsto\left(\mathcal{I}+\eta A\right)x k times over the 
  #vector x, again we take n=32, mu=0.2 and x=X=x0
  A=A(n)
  #create the A matrix of 32x32 entries
  I=diag(n)
  #take the 32x32 identity matrix
  D=I+mu*A
  #construct the matrix we will be applying to X k times
  DK=D%^%k
  #instead of looping, we multiply D by itself k times
  xk=DK%*%x
  #which is equivalent to performing the map k times over X
  return(xk)
}

tr<-1:20
b<-lapply(tr,Prod)
#We perform the map 20 times and store the results of each iteration to plot the
#entries of those vectors.
c<-data.frame(do.call(rbind, b))
C1<-rep(matrix(1:32,ncol = 1),20)
Delta<-rbind(matrix(rep(1,32),ncol=1),matrix(rep(2,32),ncol=1),
             matrix(rep(3,32),ncol=1),matrix(rep(4,32),ncol=1),
             matrix(rep(5,32),ncol=1),matrix(rep(6,32),ncol=1),
             matrix(rep(7,32),ncol=1),matrix(rep(8,32),ncol=1),
             matrix(rep(9,32),ncol=1),matrix(rep(10,32),ncol=1),
             matrix(rep(11,32),ncol=1),matrix(rep(12,32),ncol=1),
             matrix(rep(13,32),ncol=1),matrix(rep(14,32),ncol=1),
             matrix(rep(15,32),ncol=1),matrix(rep(16,32),ncol=1),
             matrix(rep(17,32),ncol=1),matrix(rep(18,32),ncol=1),
             matrix(rep(19,32),ncol=1),matrix(rep(20,32),ncol=1))
Data<-data.frame(cbind(C1,c$do.call.rbind..b.,Delta))
colnames(Data)<-c('Entry','Value','Iteration')
PLOT1<-ggplot(data = Data,aes(x=Entry,y=Value,size=Iteration,color=Iteration))+
  geom_point(aes(frame = Iteration))+ggtitle('Applying (Id+muA)x n-times')

## Warning: Ignoring unknown aesthetics: frame

• Now, we perform the mapping to a random vector of dimension \(32\)

• As \(k\) gets larger, the entries of the initial vector \(x^{(0)}\) diffuse, meaning that the distribution of the entries will get stabilized: the influence of the matrix \(A\) over the noise of the vector \(x^{(0)}\) will make the high frequency oscillations of \(x^{(0)}\) dissapear quickly (after not that many iterations) and the low oscillations of \(x^{(0)}\) dissapear slowly (after many iterations).

• Finally, at the \(\infty-\text{th}\) iteration \(x^{(\infty)}\) will be an stationary distribution for map \(\left(\mathcal{I}+\eta A\right)\), in other words:

\[ x^{(\infty)}=\left(\mathcal{I}+\eta A\right)x^{(\infty)}\implies\\ Ax^{(\infty)}=0 \]

• In other words, this iterative map approximates an eigen vector of \(A\).

b)

For a fixed value of \(n\), and for each \(k\in\{0,\dots,n-1\}\), let \(v_k\) be the vector

\[ v_k=\frac{1}{\sqrt{n}}\begin{bmatrix}e^{\frac{(2\pi i)0k}{n}}\\e^{\frac{(2\pi i)1k}{n}}\\e^{\frac{(2\pi i)2k}{n}}\\\vdots\\e^{\frac{(2\pi i)(n-1)k}{n}}\\\end{bmatrix} \]

Show algebraically, and verify computationally, that the vectors \(v_k\) are eignevectors of \(A\). find the corresponding eigenvalues \(\lambda_k\)

• Notice that the matrix \(A\) has the following form:

\[ A=\begin{bmatrix} -2&1&0&0&\cdots&0&0&1\\ 1&-2&1&0&\cdots&0&0&0\\ 0&1&-2&1&\cdots&0&0&0\\ 0&0&1&2&\cdots&0&0&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\\ 0&0&0&0&\cdots&-2&1&0\\ 0&0&0&0&\cdots&1&-2&1\\ 0&0&0&0&\cdots&0&1&-2 \end{bmatrix}= \begin{bmatrix} a_0&a_1&a_2&\cdots&a_{n-1}\\ a_{n-1}&a_0&a_1&\cdots&a_{n-2}\\ a_{n-2}&a_{n-1}&a_0&\cdots&a_{n-3}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ a_1&a_2&a_3&\cdots&a_{0} \end{bmatrix} \]

• Suppose \(x=[x_0,x_1,\dots,x_{n-1}]\in\mathbb{C}^{n}\) and \(\lambda\in\mathbb{C}\) such that:

\[ Ax=\lambda x \]

• Therefore

\[ \forall_{m=0,1,\dots,n-1}~~~\lambda x_m=\sum_{k=0}^{n-1-m}a_{k}x_{k+m}+\sum_{k=n-m}^{n-1}a_{k}x_{n-m} \] - Suppose \(x_k=\rho^{k}\in\mathbb{C}\), then

\[ \lambda=\sum_{k=0}^{n-1-m}a_{k}\rho^k+\rho^{-n}\sum_{k=n-m}^{n-1}a_{k}\rho^k \]

• If we choose \(\rho^{-n}=1\), i.e. \(\rho\) is one of the \(n\) distinct complex \(n^{\text{th}}\) roots of the unity, then

\[ \lambda=\sum_{k=0}^{n-1}a_{k}\rho^k \]

• Thereby

\[ x=\frac{1}{\sqrt{n}}\begin{bmatrix}1\\\rho\\\rho^2\\\vdots\\\rho^{n-1}\end{bmatrix} \]

• Choosing \(\rho_m\) (meaning that \(x^{m}\) is the \(m^{\text{th}}\) eigenvector of \(A\)) as the \(m^{\text{th}}\) complex root of the unity such that \(\rho_m^{n}=1\), i.e. 

\[ \rho_m=e^{\frac{(-2\pi i)m}{n}} \]

• We will arrive at the following:

\[ \lambda_m=\sum_{k=0}^{n-1}a_{k}e^{\frac{(-2\pi i)mk}{n}} \]

• And

\[ x^{m}=\frac{1}{\sqrt{n}}\begin{bmatrix}1\\ e^{\frac{(-2\pi i)m}{n}}\\ e^{\frac{(-2\pi i)m2}{n}}\\\vdots\\ e^{\frac{(-2\pi i)m(n-1)}{n}}\end{bmatrix} \]

c)

Plot enough of the eigenvectors so you get a sense for what they look like.

d)

Without any further computation, and by only considering eigenvectors and their corresponding eigenvalues, describe the evolution ofthe iterative map in equation \(x^{(k+1)}=\left(\mathcal{I}+\eta A\right)x^{(k)}\) for an initial vector of the form

\[ x^{(0)}=\left[\cos\left(\frac{2\pi j}{n}\right)+5\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]\in\mathbb{R}^n \]

• Notice that

\[ x^{(0)}=\left[\cos\left(\frac{2\pi j}{n}\right)+5\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]=\\\left[\cos\left(\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]+5\left[\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]= \]

\[ \left[\mathcal{Re}\left(e^{\frac{-2\pi i}{n}j}\right)~,~j=0,\dots,n-1\right]+5\left[\mathcal{Re}\left(e^{\frac{-2\pi i}{n}4j}\right)~,~j=0,\dots,n-1\right]= \]

\[ \mathcal{Re}\left(v_1\right)+5\cdot\mathcal{Re}\left(v_4\right) \]

• Thereby

\[ A x^{(0)}=A\left(\mathcal{Re}\left(v_1\right)+5\cdot\mathcal{Re}\left(v_4\right)\right)=A\mathcal{Re}\left(v_1\right)+5A\mathcal{Re}\left(v_4\right)= \]

\[ \mathcal{Re}\left(\lambda_1v_1\right)+5\mathcal{Re}\left(\lambda_4v_4\right)= \] \[ \mathcal{Re}\left(\lambda_1\right)\cdot\left[\cos\left(\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]+5\mathcal{Re}\left(\lambda_4\right)\cdot\left[\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right] \]

• Therefore, when performing the iterative map at the \(n^{th}\) step, we get:

\[ x^{(k+1)}=\left(\mathcal{I}+\eta A\right)x^{(k)}=\left(\mathcal{I}+\eta A\right)^{k}x^{(0)}= \]

\[ \left(1+\eta\mathcal{Re}(\lambda_1)\right)^{k}\cdot\left[\cos\left(\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]+\\5\left(1+\eta\mathcal{Re}(\lambda_4)\right)^{k}\cdot\left[\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]= \]

\[ \left(1-\eta\cdot0.03842944\right)^{k}\cdot\left[\cos\left(\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]+\\5\left(1-\eta\cdot0.5857864\right)^{k}\cdot\left[\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right] \]

• Since \(\eta=0.1\), as \(k\rightarrow\infty\)

\[ \left(1-\eta\cdot0.03842944\right)^{k}\rightarrow 0\\5\left(1-\eta\cdot0.5857864\right)^{k}\rightarrow 0 \]

• Then

\[ x^{(k+1)}\rightarrow \hat{0} \]

• Depending on the eigenvalues of \(A\), we can determine the effects of the iterative map over any initial distribution \(x^{(0)}\)

e)

Repeat these steps for the \(n\times n\) matrix \(B\), defined by

\[ B=\begin{bmatrix} 0&1&0&0&\cdots&0&0&-1\\ -1&0&1&0&\cdots&0&0&0\\ 0&-1&0&1&\cdots&0&0&0\\ 0&0&-1&0&\cdots&0&0&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\\ 0&0&0&0&\cdots&0&1&0\\ 0&0&0&0&\cdots&-1&0&1\\ 1&0&0&0&\cdots&0&-1&0 \end{bmatrix} \]

B<-function(n){
  #Creates the nxn Circulant Matrix, takes n as the dimention of the matrix
  d1<-rep(0,n) #main diagonal
  d2<-rep(1,n-1) 
  d3<-rep(-1,n-1)
  M<-triDiag(d1,d2,d3) #creates a tridiagonal matrix using d1 and d2 as the main 
  # and upper/lower diagonals
  M[n,1]=1 #making sure this is a Circulant Matrix
  M[1,n]=-1 #making sure this is a Circulant Matrix
  return(M)
}

X0<-function(n){
  #Define x^{(0)}
  Z=rep(0,n)
  for (i in 1:n){
    Z[i]=cos(2*pi*(i-1)/n)+5*cos(4*2*pi*(i-1)/n)
  }
  return(Z)
}

Xo=X0(32)

• Let’s perform the mapping to the vector \(x^{(0)}=\left[\cos\left(\frac{2\pi j}{n}\right)+5\cos\left(4\cdot\frac{2\pi j}{n}\right)~,~j=0,\dots,n-1\right]\in\mathbb{R}^n\)

• Notice that \(\left(\mathcal{I}+\eta B\right)\) is ‘moving’ the entries of the distribution from ‘right’ to ‘left’, in this case it is because the matrix \(B\) is associated to a ‘Transportation Equation’, i.e. how a scalar quantity is transported in a space.