In class, we saw that $o = [1,1,\ldots,1,1]$ is an eigenvector of $M^T$ with eigenvalue $\lambda = 1$ for any Markov matrix $M$.
(a) If $x_k$ is an eigenvector of $M$ ($M x_k = \lambda_k x_k$) for any other eigenvalue $\lambda_k \ne 1$ of $M$, show that we must have $o^T x_k = 0$: it must be orthogonal to $o$. (Hint: use $o^T = o^T M$.)
(b) Check your result from (a) numerically for a random $5 \times 5$ Markov matrix M = rand(5,5); M = M ./ sum(M, dims=1), with eigenvalues eigvals(M) and eigenvectors X = eigvecs(M). (Do using LinearAlgebra to get eigvecs and eigvals.)
(Note: if you have a long vector v, Julia only shows a few elements by default, but you can show all the elements with @show v. You can also look at the absolute values of the elements with abs.(v), which can be easier to read than complex numbers in checking that entries are small.)
(c) If we expand an arbitrary $x$ in an eigenvector basis $x = c_1 x_1 + \cdots + c_m x_m$, letting $x_m$ be a steady-state eigenvector ($\lambda_m = 1$) and supposing all of the other eigenvalues are $\ne 1$, show that $o^T x$ gives us a simple formula for $c_m = \_\_\_\_\_\_\_\_$.
(d) Hence, if all other eigenvalues have magnitude $<1$, then $M^n x \to \_\_\_\_\_\_\_\_$ (simple formula in $o,x,x_m$) as $n \to \infty$. Check this formula against M^100 * [1,2,3,4,5] for your M from (b).
(a) As $M$ is Markov, we have $o^T = o^TM$. Let us multiply both sides of $o^T = o^TM$ against the vector $x_k$, where $x_k$ is an eigenvector such that $Mx_k = \lambda_k x_k$ and $\lambda_k \neq 1$. We have $$o^Tx_k = o^TMx_k$$ $$\implies o^Tx_k = \lambda_k o^Tx_k$$ $$\implies (1 - \lambda_k) o^Tx_k = 0$$ $$\implies o^T x_k = 0,$$ where the last implication follows since $\lambda_k \neq 1$.
(b) The code is below, and indeed the dot products $o^T X$ are $\approx 0$ (up to roundoff errors) except for the (last) eigenvector corresponding to $\lambda = 1$.
using LinearAlgebra;
# generate Markov matrix
M = rand(5, 5);
M = M ./ sum(M, dims=1);
# compute eigenvalues and eigenvectors
λ = eigvals(M);
X = eigvecs(M);
# check that o^T x_k is 0 for eigenvectors whose corresponding eigenvalue is not 1
o = [1, 1, 1, 1, 1];
result = o' * X;
@show λ
@show abs.(result)
λ = ComplexF64[-0.025239843953041324 - 0.17137257909892123im, -0.025239843953041324 + 0.17137257909892123im, 0.0569700076493306 + 0.0im, 0.1222868283145958 + 0.0im, 1.0000000000000002 + 0.0im] abs.(result) = [4.726604209672303e-16 4.726604209672303e-16 4.440892098500626e-16 2.3592239273284576e-16 2.154503435713057]
1×5 Matrix{Float64}:
4.7266e-16 4.7266e-16 4.44089e-16 2.35922e-16 2.1545
(c) Taking the dot product of $o$ with $x$ gives $$o^Tx = o^T(c_1x_1 + c_2x_2 + \ldots + c_mx_m)$$ $$= c_1o^Tx_1 + c_2o^Tx_2 + \ldots + c_mo^Tx_m.$$ From (a), we know that $o^Tx_i = 0$ for $i \neq m$. Thus, we have $$o^Tx = c_mo^Tx_m \implies c_m = \frac{o^Tx}{o^Tx_m}.$$
(d) We have $M^nx = c_1 \lambda_1^n x_1 + \ldots + c_m \lambda_m^n x_m$. As $n\to \infty$, $\lambda_i^n \to 0$ for $i \neq m$. So, $$M^nx \to c_m \lambda_m^n x_m = \frac{o^Tx}{o^Tx_m} x_m$$ where we use the formula for $c_m$ from (c). The code below observes this convergence.
x = [1, 2, 3, 4, 5];
# calculate true result
xTrue = M^100 * x;
# calculate expected result from part (c)'s work
xExpected = ((o' * x) / (o' * X[:,5])) .* X[:,5];
# compare the two
@show xTrue
@show xExpected
@show abs.(xTrue - xExpected)
xTrue = [1.9338510947374465, 3.8824772333448783, 3.9405657589028347, 3.0603805344264408, 2.182725378588311] xExpected = ComplexF64[1.9338510947374572 + 0.0im, 3.8824772333449022 + 0.0im, 3.9405657589028564 + 0.0im, 3.0603805344264576 + 0.0im, 2.1827253785883243 + 0.0im] abs.(xTrue - xExpected) = [1.0658141036401503e-14, 2.398081733190338e-14, 2.1760371282653068e-14, 1.687538997430238e-14, 1.3322676295501878e-14]
5-element Vector{Float64}:
1.0658141036401503e-14
2.398081733190338e-14
2.1760371282653068e-14
1.687538997430238e-14
1.3322676295501878e-14
From Strang, section 6.2. Consider $A = \begin{pmatrix} 0.6 & 0.4 \\ 0.4 & 0.6 \end{pmatrix}$ and $B = \begin{pmatrix} 0.6 & 0.9 \\ 0.1 & 0.6 \end{pmatrix}$. For this problem you keep in mind the diagonalization of matrices like $A$ and $B$.
(a) Which of $A^n$ or $B^n$ (or both, or neither) go $\to \begin{pmatrix}0 & 0 \\ 0 & 0 \end{pmatrix}$ as $n \to \infty$? Double-check your answer by looking at $A^{100}$ and $B^{100}$ in Julia — these are approaching matrices of rank ________ and ________, respectively.
(b) For what values of the real scalar $\mu$ is $\sqrt{A - \mu I}$ a real matrix? Check your answer by trying sqrt(A - μ*Ι) for a few values of μ in Julia (do using LinearAlgebra to get I).
(a) $A$ is a positive Markov matrix, so it must have one eigenvalue $\lambda_1 = 1$ and another eigenvalue with magnitude $|\lambda_2| < 1$. So, with no calculation, we can see that we must have $$A^n = X\begin{pmatrix}\lambda_1 &\\&\lambda_2\end{pmatrix}^nX^{-1} = X\begin{pmatrix}\lambda_1^n&\\&\lambda_2^n\end{pmatrix}X^{-1} \to X\begin{pmatrix}1&\\&1\end{pmatrix}X^{-1}.$$ As one of the diagonal entries is zero, $A^n$ is approaching a matrix of $\boxed{\text{rank } 1}$ (and so not going to $\begin{pmatrix}0&0\\0&0\end{pmatrix}$)
We can compute the eigenvalues explicitly, of course. From the quadratic formula for $\det(A-\lambda I) = \lambda^2 - \mathrm{trace}(A) \lambda + \det A$, we get $\lambda = 0.6 \pm \sqrt{0.6^2 - (0.6^2 - 0.4^2)} = 0.6 \pm 0.4 = 1 \textrm{ or } 0.2$. We can also check this in Julia below. So the diagonalization of $A$ is $X\begin{pmatrix}1&\\&0.2\end{pmatrix}X^{-1}$ where $X$ is the matrix of $A$'s eigenvectors. As $n \to \infty$, we have $$A^n = X\begin{pmatrix}1&\\&0.2\end{pmatrix}^nX^{-1} = X\begin{pmatrix}1^n&\\&0.2^n\end{pmatrix}X^{-1} \to X\begin{pmatrix}1&\\&0\end{pmatrix}X^{-1}.$$ As one of the diagonal entries is zero, $A^n$ is approaching a matrix of $\boxed{\text{rank } 1}$ (and so not going to $\begin{pmatrix}0&0\\0&0\end{pmatrix}$).
$B$ is not a Markov matrix, so we'll just calculate the eigenvalues explicitly. Via the quadratic formula, we get $\lambda = 0.6 \pm \sqrt{0.6^2 - (0.6^2 - 0.9 \times 0.1)} = 0.6 \pm \sqrt{0.09} = 0.6 \pm 0.3 = 0.3 \textrm{ or } 0.9$ So the diagonalization of $B$ is $X\begin{pmatrix}0.3&\\&0.9\end{pmatrix}X^{-1}$ where $X$ is the matrix of $B$'s eigenvectors. As $n \to \infty$, we have $$B^n = X\begin{pmatrix}0.3&\\&0.9\end{pmatrix}^nX^{-1} = X\begin{pmatrix}0.3^n&\\&0.9^n\end{pmatrix}X^{-1} \to X\begin{pmatrix}0&\\&0\end{pmatrix}X^{-1}.$$ As both of the diagonal entries are zero, $B^n$ is approaching a matrix of $\boxed{\text{rank } 0}$, and so must be going to $\begin{pmatrix}0&0\\0&0\end{pmatrix}$.
The code below confirms this.
A = [0.6 0.4; 0.4 0.6];
B = [0.6 0.9; 0.1 0.6];
using LinearAlgebra
# compute eigenvalues
λ_A = eigvals(A);
λ_B = eigvals(B);
@show λ_A
@show λ_B
# check high powers
@show A^100;
@show B^100;
λ_A = [0.19999999999999996, 1.0] λ_B = [0.29999999999999993, 0.8999999999999999] A ^ 100 = [0.5 0.5; 0.5 0.5] B ^ 100 = [1.3280699443793745e-5 3.984209833138124e-5; 4.426899814597916e-6 1.3280699443793748e-5]
(b) Note that the eigenvalues of $A - \mu I$ are $0.2 - \mu$ and $1.0 - \mu$, with the same corresponding eigenvectors as $A$. Then $A - \mu I = X\begin{pmatrix}0.2 - \mu&\\&1.0 - \mu\end{pmatrix}X^{-1}$ and in particular, $$\sqrt{A - \mu I} = X\begin{pmatrix}0.2 - \mu&\\&1.0 - \mu\end{pmatrix}^{\frac{1}{2}}X^{-1} = X\begin{pmatrix}\sqrt{0.2 - \mu}&\\&\sqrt{1.0 - \mu}\end{pmatrix}X^{-1}.$$ Note that $X$ and $X^{-1}$ are real matrices. For $\sqrt{A - \mu I}$ to be a real matrix, we thus need $\sqrt{0.2 - \mu}$ and $\sqrt{1.0 - \mu}$ to be real, or equivalently, $\mu \le 0.2$ and $\mu \le 1.0$. So the answer is $\boxed{\mu \le 0.2}$.
The code below confirms this.
# test μ from 0 to 1 in 0.1 increments
for μ in LinRange(0.0, 1.0, 11)
X = sqrt(A - μ * I);
println("μ: ", μ);
println(X);
end
μ: 0.0 [0.7236067977499787 0.276393202250021; 0.276393202250021 0.7236067977499787] μ: 0.1 [0.6324555320336757 0.3162277660168379; 0.3162277660168379 0.6324555320336757] μ: 0.2 [0.4472135954999578 0.4472135954999578; 0.4472135954999578 0.4472135954999578] μ: 0.3 ComplexF64[0.4183300132670377 + 0.15811388300841897im 0.4183300132670377 - 0.15811388300841897im; 0.4183300132670377 - 0.15811388300841897im 0.4183300132670377 + 0.15811388300841897im] μ: 0.4 ComplexF64[0.38729833462074165 + 0.223606797749979im 0.38729833462074165 - 0.223606797749979im; 0.38729833462074165 - 0.223606797749979im 0.38729833462074165 + 0.223606797749979im] μ: 0.5 ComplexF64[0.35355339059327373 + 0.27386127875258304im 0.35355339059327373 - 0.27386127875258304im; 0.35355339059327373 - 0.27386127875258304im 0.35355339059327373 + 0.27386127875258304im] μ: 0.6 ComplexF64[0.3162277660168379 + 0.3162277660168379im 0.3162277660168379 - 0.3162277660168379im; 0.3162277660168379 - 0.3162277660168379im 0.3162277660168379 + 0.3162277660168379im] μ: 0.7 ComplexF64[0.27386127875258304 + 0.35355339059327373im 0.27386127875258304 - 0.35355339059327373im; 0.27386127875258304 - 0.35355339059327373im 0.27386127875258304 + 0.35355339059327373im] μ: 0.8 ComplexF64[0.22360679774997888 + 0.38729833462074165im 0.22360679774997888 - 0.38729833462074165im; 0.22360679774997888 - 0.38729833462074165im 0.22360679774997888 + 0.38729833462074165im] μ: 0.9 ComplexF64[0.15811388300841892 + 0.4183300132670377im 0.15811388300841892 - 0.4183300132670377im; 0.15811388300841892 - 0.4183300132670377im 0.15811388300841892 + 0.4183300132670377im] μ: 1.0 ComplexF64[0.0 + 0.4472135954999578im 0.0 - 0.4472135954999578im; 0.0 - 0.4472135954999578im 0.0 + 0.4472135954999578im]