Mini-Lecture 6-29 - Using the new coordinates to visualize the behavior of a linear model.ipynb

Kernel: SageMath 9.2

import numpy as np
from ipywidgets import IntSlider, SelectionSlider, Checkbox, Text, \
    interactive_output, VBox, HBox, HTMLMath, Layout

Learning goals:

• For a discrete-time linear model, once you look at the state space using the $R,S$ coordinates*, the system behaves exactly like a diagonal (decoupled) linear model. That is: $$\begin{bmatrix} X_{t+1} \\ Y_{t+1} \end{bmatrix} = \begin{pmatrix} \dotsm & \dotsm \\ \dotsm & \dotsm \end{pmatrix} \begin{bmatrix} X_t \\ Y_t \end{bmatrix} \implies \begin{bmatrix} R_{t+1} \\ S_{t+1} \end{bmatrix} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \begin{bmatrix} R_t \\ S_t \end{bmatrix}$$

• This process is called diagonalizing the matrix, or diagonalizing the linear model.

• This means everything we know about diagonal matrices (in particular, the concept of the dominant eigenvalue) now carries over to any discrete-time linear model!

* The $R,S$ coordinates are defined using the eigenvectors of the model as a basis, and the $\lambda_1$ and $\lambda_2$ on the diagonal of the diagonalized matrix are the eigenvalues of the model.

def plot_graph_paper(T, **options):
    color1, color2 = options.get("colors", ("red", "green"))
    v1_color = colors[color1].darker(0.25)
    v2_color = colors[color2].darker(0.25)
    grid1_color = colors[color1].lighter(0.75)
    grid2_color = colors[color2].lighter(0.75)
    axis1_label, axis2_label = options.get("axes_labels", ("$R$", "$S$"))
    show_axes = options.get("show_axes", True)
    gridlines1, gridlines2 = options.get("gridlines", (6, 6))
    v1 = T.column(0)
    v2 = T.column(1)
    xymax = min(g / a for g, a in zip((gridlines1, gridlines2), ~T*vector((1,1))) if a != 0)
    xmax, ymax = options.get("xyrange", (xymax, xymax))
    p = Graphics()
    if options.get("show_grid", True):
        for n in range(-gridlines1, gridlines1 + 1):
            if n == 0 and show_axes:
                continue
            p += line((v1*n - gridlines2*v2, v1*n + gridlines2*v2), color=grid2_color, thickness=0.5)
        for n in range(-gridlines2, gridlines2 + 1):
            if n == 0 and show_axes:
                continue
            p += line((v2*n - gridlines1*v1, v2*n + gridlines1*v1), color=grid1_color, thickness=0.5)
    if show_axes:
        p += line((-gridlines1*v1, gridlines1*v1), color=color1, thickness=1.5)
        p += line((-gridlines2*v2, gridlines2*v2), color=color2, thickness=1.5)
        p += text(axis1_label, 1.1*min(abs(xmax / v1[0]), abs(ymax / v1[1]))*v1, color=color1, fontsize=16)
        p += text(axis2_label, 1.1*min(abs(xmax / v2[0]), abs(ymax / v2[1]))*v2, color=color2, fontsize=16)
    if options.get("show_vectors", True):
        p += plot(v1, color=v1_color, thickness=2, zorder=5)
        p += plot(v2, color=v2_color, thickness=2, zorder=5)
    p.set_axes_range(xmin=-xmax, xmax=xmax, ymin=-ymax, ymax=ymax)
    return p

def initchoice(a, vmax):
    if abs(a) < 1:
        return randint(round(0.5*vmax + 1), vmax)
    if abs(a) > 1:
        return randint(1, round(0.5*vmax - 1))
    return randint(1, vmax - 1)

def nondiagonal_interactive(T, eigenvalues, **options):
    initial_xy = options.get("initial_state", None)
    state_vars = options.get("state_vars", ("X", "Y"))
    rsmax = 10
    tmax = 50
    lambda1, lambda2 = eigenvalues
    D = diagonal_matrix(RDF, eigenvalues)
    M = matrix(RDF, T*D*~T).round(4)
    if initial_xy is None:
        initial_rs = vector([initchoice(lambda1, rsmax), initchoice(lambda2, rsmax)])
        initial_xy = vector([round(x,3) for x in T * initial_rs])
    else:
        initial_rs = vector([round(x, 3) for x in ~T * vector(RDF, initial_xy)])
    solution_rs = solution_r = solution_s = solution_xy = []

    dont_update = False
    def update(rt, st, rst, xyt, show_eigen, show_grid, axis1max, axis2max):
        if dont_update or len(solution_rs) == 0:
            return
        rsmax = ceil(axis1max / abs(np.array(T).flatten()).min())
        p1  = plot_graph_paper(T, show_vectors=show_eigen, show_axes=show_eigen, 
                               show_grid=show_grid, gridlines=(rsmax,rsmax), 
                               xyrange=(axis1max,axis1max))
        p1 += list_plot(solution_xy[:xyt+1], size=30, color="black", zorder=10)
        p1 += list_plot(solution_xy[:xyt+1], plotjoined=True, color="lightblue", zorder=5)
        if show_eigen:
            p1 += list_plot(solution_r[:rt+1], size=30, color="red")
            p1 += list_plot(solution_s[:st+1], size=30, color="green")
        p2  = list_plot(solution_rs[:rst+1], size=30, color="black", zorder=10)
        p2 += list_plot(solution_rs[:rst+1], plotjoined=True, color="lightblue", zorder=5)
        p2 += list_plot(solution_rs[:rt+1] * (1,0), size=30, color="red")
        p2 += list_plot(solution_rs[:st+1] * (0,1), size=30, color="green")
        if show_grid:
            p1 += line((solution_r[xyt], solution_xy[xyt], solution_s[xyt]), color="gray")
            p1 += point(solution_r[xyt], size=50, color="red")
            p1 += point(solution_s[xyt], size=50, color="green")
            p2 += line((solution_rs[rst] * (1,0), solution_rs[rst], solution_rs[rst] * (0,1)), 
                       color="gray")
            p2 += point(solution_rs[rst] * (1,0), size=50, color="red")
            p2 += point(solution_rs[xyt] * (0,1), size=50, color="green")
        p1.set_axes_range(xmin=-axis1max, xmax=axis1max, ymin=-axis1max, ymax=axis1max)
        p2.set_axes_range(xmin=-axis2max, xmax=axis2max, ymin=-axis2max, ymax=axis2max)
        p1.axes_labels(["$%s$" % label for label in state_vars])
        p2.axes_labels(("$R$", "$S$"))
        p1.set_aspect_ratio(1)
        p2.set_aspect_ratio(1)
        both = multi_graphics([(p1, (0, 0.25, 0.57, 0.57)), (p2, (0.43, 0.25, 0.57, 0.57))])
        both.show(figsize=12)

    def change_initial_rs(change):
        nonlocal dont_update, solution_rs, solution_r, solution_s, solution_xy
        initial_rs = change["new"].strip("()[]").split(",")
        try:
            initial_rs = vector([float(x) for x in initial_rs])
        except:
            return
        solution_rs = np.array([D^t * initial_rs for t in range(tmax+1)], dtype=float)
        solution_r = [initial_rs[0] * lambda1^t * T.column(0) for t in range(tmax+1)]
        solution_s = [initial_rs[1] * lambda2^t * T.column(1) for t in range(tmax+1)]
        solution_xy = [a + b for a, b in zip(solution_r, solution_s)]
        dont_update = True
        rt.value = 0
        st.value = 0
        rst.value = 0
        xyt.value = 1
        dont_update = False
        xyt.value = 0

    def change_initial_xy(change):
        initial_xy = change["new"].strip("()[]").split(",")
        try:
            initial_xy = vector([float(x) for x in initial_xy])
        except:
            return
        state_rs.value = str(vector([round(x, 3) for x in ~T * initial_xy]))

    p2b = str.maketrans("()", "[]")
    label = ((r"$ \begin{bmatrix} %s_{t+1} \\ %s_{t+1} \end{bmatrix} = " % state_vars) + 
             latex(M) + (r"\begin{bmatrix} %s_t \\ %s_t \end{bmatrix} \qquad " % state_vars) + 
             r"\vec{v}_1 = " + latex(T[:,0]).translate(p2b) + r" \qquad " + 
             (r"\text{with } \lambda_1 = %.2f \qquad" % lambda1) + 
             r"\vec{v}_2 = " + latex(T[:,1]).translate(p2b) + r" \qquad " + 
             (r"\text{with } \lambda_2 = %.2f $" % lambda2) )
    equation = HTMLMath(label)
    layout = Layout(width="250px", margin="0px 0px 0px 0px")
    xyt = IntSlider(min=0, max=tmax, layout=layout, description="Sol'n (left):")
    rst = IntSlider(min=0, max=tmax, layout=layout, description="Sol'n (right):")
    rt  = IntSlider(min=0, max=tmax, layout=layout, description="R steps:")
    st  = IntSlider(min=0, max=tmax, layout=layout, description="S steps:")
    show_eigen = Checkbox(value=False, layout=layout, description="Show eigenlines")
    show_grid = Checkbox(value=False, layout=layout, description="Show gridlines")
    zoom1 = SelectionSlider(options=[1, 2, 5, 10, 20, 50, 100, 200, 500, 1000], value=10, 
                            layout=layout, description="Zoom (left):")
    zoom2 = SelectionSlider(options=[1, 2, 5, 10, 20, 50, 100, 200, 500, 1000], value=10, 
                            layout=layout, description="Zoom (right):")
    state_rs = Text(continuous_update=False, layout=layout, description="Initial (R,S):")
    state_xy = Text(value=str(initial_xy), continuous_update=False, 
                    layout=layout, description="Initial (%s,%s):" % state_vars)
    state_rs.observe(change_initial_rs, "value")
    state_xy.observe(change_initial_xy, "value")
    state_rs.value = str(initial_rs)
    output = interactive_output(update, dict(rt=rt, st=st, rst=rst, xyt=xyt, 
                                show_eigen=show_eigen, show_grid=show_grid, 
                                axis1max=zoom1, axis2max=zoom2))
    controls1 = HBox((state_xy, state_rs, show_eigen, show_grid))
    controls2 = HBox((rt, st))
    controls3 = HBox((xyt, zoom1, rst, zoom2))
    display(VBox((equation, controls1, controls2, controls3, output)))

v1 = vector((2, 1))
v2 = vector((3, -1))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (1.03, 0.28), initial_state=(2,6), state_vars=("J", "A"))

v1 = vector((3, 4))
v2 = vector((7, -3))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (1.03, 0.29), initial_state=(20,2), state_vars=("J", "A"))

M = matrix(RDF, [[0.36, 0.35], [0.24, 0.82]])
M0 = matrix(ZZ, 100*M)
show(M0)
for evalue, evectors, mult in M0.eigenvectors_right():
    print(evalue, evectors[0])

96 (1, 12/7)
22 (1, -2/5)

v1 = vector(RDF, (7, 12))*0.25
v2 = vector(RDF, (5, -2))*0.5
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (0.96, 0.22), initial_state=(17,8), state_vars=("J", "A"))

v1 = vector((1, 1))
v2 = vector((1, -2))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (1.3, 1.12), initial_state=(0,3))

v1 = vector((-1, 1))
v2 = vector((2, 1))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (-0.9, 1.2), initial_state=(1,5))

v1 = vector((-1, 1))
v2 = vector((1, 2))
T = column_matrix((v1, v2))
nondiagonal_interactive(T, (-0.7, 0.89), initial_state=(10,5))

Conclusions:

1. For a discrete-time linear model, if you use the eigenvectors of the model to define a new $R,S$ coordinate system, then once you look at the state space using these new coordinates, the system behaves exactly like a diagonal (decoupled) linear model: $$\begin{bmatrix} X_{t+1} \\ Y_{t+1} \end{bmatrix} = \begin{pmatrix} \dotsm & \dotsm \\ \dotsm & \dotsm \end{pmatrix} \begin{bmatrix} X_t \\ Y_t \end{bmatrix} \implies \begin{bmatrix} R_{t+1} \\ S_{t+1} \end{bmatrix} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \begin{bmatrix} R_t \\ S_t \end{bmatrix}$$ Furthermore, the entries on the diagonal are the eigenvalues of the model.

2. This process is called diagonalizing the matrix, or diagonalizing the linear model.

3. This means we can apply what we know about the dominant eigenvalue to any discrete-time linear model: the long-term behavior will be exponential growth or decay, depending on the dominant eigenvalue (the one with largest absolute value, and whether that's $>1$ or $<1$). This long-term growth or decay will be parallel to the dominant eigenvector line (the one corresponding to the dominant eigenvalue).