Python meets linear algebra

Statistics and linear algebra seem like a pair that should always be married to each other. This is hardly the case however. If you are like me and you fancy remembering the equation for say the correlation between two variables without memorizing the formula, then you should probably pick a linear algebra text by Gilbert Strang or head over to his MIT vidoes on youtube which can be found here. Trust me this will be very illuminating.

I will give some examples in case you are still not convinced. We have two variables

$$\begin{array}{c}x=\left[\begin{array}{c}1\\ 2\\ 3\\ 4\\ 5\\ 6\end{array}\right]y=\left[\begin{array}{c}7\\ 8\\ 9\\ 10\\ 11\\ 12\end{array}\right]\end{array}$$

Pardon my lack of creativity but these two vectors will do. If you wanted to calculate the correlation between these two vectors(variables), from statistics we would need to do $\frac{Cov(x,y)}{{\sigma }_x{\sigma }_y}$$\frac{Cov(x,y)}{\sigma_x \sigma_y}$​ . Step in linear algebra and the correlation between two vectors(variables) is basically the angle between the two variables given by $cos\theta =\frac{\overrightarrow{x}.\overrightarrow{y}}{|x||y}$$cos\theta = \frac{\vec{x}.\vec{y}}{|x||y}$ . Comparing these two equations it becomes apparent that the covariance of two vectors is their dot product whilst their standard deviation(after centering by taking away the mean of the vectors from each component) is basically the length of the vector. This for me is a pretty tidy way to view statistical objects from a linear algebra perspective as I think it is easier to understand what the dot product and length of two vectors are and then relate it to their statistical equivalent.

I think the most interesting of these objects is the regression of y on x. Linear algebra tells us that if y is not in the column space of x then there is no way to arrive at a solution however, we can arrive at an approximation by projecting y orthogonally on x. This is a pretty neat discovery as one can easily imagine two lines x and y, going in different directions, then from the shadow y makes on x, we can create a new vector $\overrightarrow{p}$$\vec{p}$ that sits on $\overrightarrow{x}$$\vec{x}$ and is orthogonal to the vector $\overrightarrow{y}-\overrightarrow{p}$$\vec{y} - \vec{p}$ . The image this creates in our head can be translated in different dimensions even if we can't represent these dimensions with pen and paper. This connection between statistics and linear algebra is beyond the scope of this article but if I have managed to pique your interest, then I think it is worthwhile to pick up a text like Geometry of multivariate statistics by Thomas D. Wickens and explore this subject further.

I have always liked to make images move around on my screen, however, until I understood better transition matrices and change of basis I never had a good grasp of how things worked under the hood. Do you feel anxious about transition matrices and change of basis? Don't worry they are just fancy names of what happens when you multiply matrices. The matrix multiplication $A\overrightarrow{x}=\overrightarrow{b}$$A\vec{x} = \vec{b}$ is all we probably need to know to understand this article. A is a transition matrix that acts on $\overrightarrow{x}$$\vec{x}$ to yield $\overrightarrow{b}$$\vec{b}$.

I want to explain how to generate a transition matrix, that flips a vector in $R^2$$R^{2}$ say $\begin{array}{c}x=\left[\begin{array}{c}1\\ 2\end{array}\right]\end{array}$$x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}$ By ${90}^{\mathrm{°}}$$90^{\degree}$ in the anticlockwise direction. This should yeild the vector $\left[\begin{array}{c}-2\\ 1\end{array}\right]$$\begin{bmatrix} -2 \\ 1 \end{bmatrix}$. Firstly I will like to talk about a basis of a space. This is basically the minimal vectors that span a particular space. For example we are in $R^2$$R^{2}$ and we want two vectors that span $R^2$$R^{2}$ we can easily go with the standard unit vectors $\left[\begin{array}{c}1\\ 0\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix}$ Or $\left[\begin{array}{c}1\\ 2\end{array}\right]\left[\begin{array}{c}-2\\ 1\end{array}\right]$$\begin{bmatrix} 1 \\ 2 \end{bmatrix} \begin{bmatrix} -2 \\ 1 \end{bmatrix}$ with a linear combination of these two vectors(and many more) we can get to any point on the plane that is generated by $R^2$$R^{2}$ for $R^3$$R^{3}$ we can use the standard unit vectors in $R^3$$R^{3}$ -> $\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$ the point of this whole discussion is to point out that to get a transition matrix, we only need to observe what happens to for example the standard unit vectors of $R^2$$R^{2}$ . This will help us generate our matrix. In the next section, it might be a bit difficult to do this and that is where change of basis comes to play but let us get on with it. When our transition matrix A acts on the vector $\left[\begin{array}{c}1\\ 0\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ to get a ${90}^{\mathrm{°}}$$90^{\degree}$ rotation the x axis becomes 0 and the y axis becomes 1 that is $\left[\begin{array}{c}0\\ 1\end{array}\right]$$\begin{bmatrix} 0 \\ 1 \end{bmatrix}$ When the transition matrix acts on $\left[\begin{array}{c}0\\ 1\end{array}\right]$$\begin{bmatrix} 0 \\ 1 \end{bmatrix}$ the cordinates change to -1 and 0 that is $\left[\begin{array}{c}-1\\ 0\end{array}\right]$$\begin{bmatrix} -1 \\ 0 \end{bmatrix}$. These two vectors give us our transition matrix $\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$\begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$. Check for yourself that if you multiply this matrix by $\left[\begin{array}{c}1\\ 2\end{array}\right]$$\begin{bmatrix} 1 \\ 2 \end{bmatrix}$ that we get $\left[\begin{array}{c}-2\\ 1\end{array}\right]$$\begin{bmatrix} -2 \\ 1 \end{bmatrix}$ .

Lets try to extend this result with a little bit of trigonometry. We know that $cos90=0$$cos90 = 0$, $sin90$$sin90$ = 1 so our transition matrix could possibly look like $\left[\begin{array}{cc}cos90& -sin90\\ sin90& cos90\end{array}\right]$$\begin{bmatrix} cos 90 & -sin 90 \\ sin 90 & cos 90 \end{bmatrix}$ in fact we can generalize this matrix for any angle $\theta$$\theta$ to become $\left[\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right]$$\begin{bmatrix} cos \theta & -sin\theta \\ sin\theta & cos\theta \end{bmatrix}$. We can implement this in python to verify the results.

x
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np

# The angles we want to test with
thetas = [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360]
plt.ioff()
for theta in thetas:
    # The transition matrix generalized with sines and cosines
    matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
    fig = plt.figure()
    ax = fig.add_subplot(111, aspect='equal')
    ax.set_xlim(-100, 100)
    ax.set_ylim(-100, 100)

    # Using numpy to multiply the matrix and the vector to get new cordinates
    vec1 = np.matmul(matrix, np.array([-25, -25]))
    vec2 = np.matmul(matrix, np.array([-25, 25]))
    vec3 = np.matmul(matrix, np.array([25, 25]))
    vec4 = np.matmul(matrix, np.array([25, -25]))
    vec5 = np.matmul(matrix, np.array([-10,24]))
    vec6 = np.matmul(matrix, np.array([-10,40]))
    vec7 = np.matmul(matrix, np.array([10,40]))
    vec8 = np.matmul(matrix, np.array([10,24]))
    

    x = [vec1[0], vec2[0], vec3[0], vec4[0]]
    y = [vec1[1], vec2[1], vec3[1], vec4[1]]
    ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))
    
    x = [vec5[0], vec6[0], vec7[0], vec8[0]]
    y = [vec5[1], vec6[1], vec7[1], vec8[1]]
    ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))
    plt.draw()
    plt.pause(0.4)
    plt.clf()

If you run the above code you will get two rectangles that look like a head and a body moving around by the angle specified by theta.

This was pretty tidy, however what happens when we want to say for example reflect about the line $\left[\begin{array}{c}1\\ 2\end{array}\right]$$\begin{bmatrix} 1 \\ 2 \end{bmatrix}$ this is not as straighforward as reflecting about ${90}^{\mathrm{°}}$$90^{\degree}$ try it. However, what we can do is to switch basis. If you have some knowledge about eigen values and eigen vectors this will come in handy even though I will explain a trick that will help us arrive at the answer. What we need to do is to find a complete set of new eigen vectors that span $R^2$$R^{2}$ we already have one. To find the other one we need a vector that is orthogonal to it. That is they form a ${90}^{\mathrm{°}}$$90^{\degree}$ angle. to achieve this we need a new vector such that if we obtain the dot product of our vector with it, we get the zero vector. This vector is $\left[\begin{array}{c}2\\ -1\end{array}\right]$$\begin{bmatrix} 2 \\ -1 \end{bmatrix}$ we can easily check that $1\ast 2+2\ast (-1)=0$$1* 2 + 2 * (-1) = 0$. What we need to do now is use $\left[\begin{array}{c}1\\ 2\end{array}\right]$$\begin{bmatrix} 1 \\ 2 \end{bmatrix}$ and $\left[\begin{array}{c}2\\ -1\end{array}\right]$$\begin{bmatrix} 2 \\ -1 \end{bmatrix}$ as our new cordinate system and instead of reflecting about the traditional cordinate system we know, we will reflect about this cordinate system. Imagine a world where the standard cordinates for y and x are the lines $\begin{array}{c}c\left[\begin{array}{c}1\\ 2\end{array}\right]\end{array}$$c\begin{bmatrix} 1 \\ 2 \end{bmatrix}$ and $\begin{array}{c}c\left[\begin{array}{c}2\\ -1\end{array}\right]\end{array}$$c\begin{bmatrix} 2 \\ -1 \end{bmatrix}$ where c is any real number as against the lines $\begin{array}{c}c\left[\begin{array}{c}1\\ 0\end{array}\right]\end{array}$$c\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ and $\begin{array}{c}c\left[\begin{array}{c}0\\ 1\end{array}\right]\end{array}$$c\begin{bmatrix} 0 \\ 1 \end{bmatrix}$this takes a bit of a stretch but once you have it, we can then think about what flipping over the y axis which is now the line $\left[\begin{array}{c}1\\ 2\end{array}\right]$$\begin{bmatrix} 1 \\ 2 \end{bmatrix}$ if we apply our previous mechanics by checking what happens to the vectors $\left[\begin{array}{c}1\\ 0\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix}$ we would gain an insight on the transition matrix in our new cordinate system. The vector $\left[\begin{array}{c}0\\ 1\end{array}\right]$$\begin{bmatrix} 0 \\ 1 \end{bmatrix}$ should remain the same since we are flipping about the y axis in our new system. However, the vector $\left[\begin{array}{c}1\\ 0\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ will become $\left[\begin{array}{c}-1\\ 0\end{array}\right]$$\begin{bmatrix} -1 \\ 0 \end{bmatrix}$. We get the transition matrix as $\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]$$\begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}$. The numbers -1 and 1 are basically the eigen values of our matrix $\left[\begin{array}{cc}1& 2\\ 2& -1\end{array}\right]$$\begin{bmatrix} 1 & 2 \\ 2 & -1 \end{bmatrix}$ try to compute it and the matrix $\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]$$\begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}$ is basically placing the eigen values in the diagonal of 2 by 2 identity matrix.

Since we humans haven't yet developed the sense of understanding this new cordinate system, if you wanted to show your new results to your friends, you will need to make some adjustments. Linear algebra comes to the rescue as we can switch from this new cordinate system to a more traditional one. There is a pretty neat technique for doing this called diagonilazation of a matrix. The formula is given by $A=SE{S}^{-1}$$A = SES^{-1}$ where S is the eigen vector matrix and E is a matrix that contains the eigen values in the diagonals. If we plug in our values we get A as $\left[\begin{array}{cc}-3/5& 4/5\\ 4/5& 3/5\end{array}\right]$$\begin{bmatrix} -3/5 & 4/5 \\ 4/5 & 3/5 \end{bmatrix}$ in standard cordinates. So if we wanted to find the transformed cordinates for say $\left[\begin{array}{c}1\\ 0\end{array}\right]$$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ this should yield $\left[\begin{array}{c}-3/5\\ 4/5\end{array}\right]$$\begin{bmatrix} -3/5 \\ 4/5 \end{bmatrix}$ . We can check our results with python by inspecting the outcome of the script below

x
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np

theta = 0
matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal')

ax.set_xlim(-100, 100)
ax.set_ylim(-100, 100)

vec1 = np.matmul(matrix, np.array([-25, -25]))
vec2 = np.matmul(matrix, np.array([-25, 25]))
vec3 = np.matmul(matrix, np.array([25, 25]))
vec4 = np.matmul(matrix, np.array([25, -25]))
vec5 = np.matmul(matrix, np.array([-10,24]))
vec6 = np.matmul(matrix, np.array([-10,40]))
vec7 = np.matmul(matrix, np.array([10,40]))
vec8 = np.matmul(matrix, np.array([10,24]))

x = [vec1[0], vec2[0], vec3[0], vec4[0]]
y = [vec1[1], vec2[1], vec3[1], vec4[1]]

ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))

x = [vec5[0], vec6[0], vec7[0], vec8[0]]
y = [vec5[1], vec6[1], vec7[1], vec8[1]]
ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))

# Set the line x = [1, 2]
x = [-100, 100]
y = [-200, 200]
ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=False))
plt.show()


matrix = np.array([[-3/5, 4/5], [4/5, 3/5]])

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal')
ax.set_xlim(-100, 100)
ax.set_ylim(-100, 100)

vec1 = np.matmul(matrix, np.array([-25, -25]))
vec2 = np.matmul(matrix, np.array([-25, 25]))
vec3 = np.matmul(matrix, np.array([25, 25]))
vec4 = np.matmul(matrix, np.array([25, -25]))
vec5 = np.matmul(matrix, np.array([-10,24]))
vec6 = np.matmul(matrix, np.array([-10,40]))
vec7 = np.matmul(matrix, np.array([10,40]))
vec8 = np.matmul(matrix, np.array([10,24]))


x = [vec1[0], vec2[0], vec3[0], vec4[0]]
y = [vec1[1], vec2[1], vec3[1], vec4[1]]

ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))

x = [vec5[0], vec6[0], vec7[0], vec8[0]]
y = [vec5[1], vec6[1], vec7[1], vec8[1]]
ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=True))

# Set the line x = [1, 2]
x = [-100, 100]
y = [-200, 200]
ax.add_patch(patches.Polygon(xy=list(zip(x,y)), fill=False))
plt.show()

The results can be extended into 3 dimensions with a little more work.

Thanks.