Numpy Matrix Stuff

Define some small matrices and vectors to do numpy operations on

\texttt{A = np.array([[1,2],[2,3],[3,4]])} = \begin{pmatrix}1&2\\2&3\\3&4\end{pmatrix}
\texttt{b = np.array([2,3])} = \begin{pmatrix}2&3\end{pmatrix}
\texttt{c = np.array([2,3,4])} = \begin{pmatrix}2&3&5\end{pmatrix}

Component-wise multiplication

Some examples of the * operator
\texttt{b * b} = \begin{pmatrix}2&3\end{pmatrix} * \begin{pmatrix}2&3\end{pmatrix} = \texttt{array([4, 9])} = \begin{pmatrix}4&9\end{pmatrix}
It doesn't make any difference the shape of b
\texttt{b.T * b} = \begin{pmatrix}2\\3\end{pmatrix} * \begin{pmatrix}2&3\end{pmatrix} = \texttt{array([4, 9])} = \begin{pmatrix}4&9\end{pmatrix}
Now for 2d stuff
\texttt{A * b} = \begin{pmatrix}1&2\\2&3\\3&4\end{pmatrix} * \begin{pmatrix}2&3\end{pmatrix} = \texttt{array([[2, 6],[4, 9],[6, 12]])} = \begin{pmatrix}2&6\\4&9\\6&12\end{pmatrix}
Again the shape of b doesn't matter
\texttt{A * b.T} = \begin{pmatrix}1&2\\2&3\\3&4\end{pmatrix} * \begin{pmatrix}2\\3\end{pmatrix} = \texttt{array([[2,6],[4,9],[6,12]])} = \begin{pmatrix}2&6\\4&9\\6&12\end{pmatrix}
But the shape of A does matter. If you want to multiply each row of A by something, then you need to transpose it
\texttt{A.T * c}=\begin{pmatrix}1&2&3\\2&3&4\end{pmatrix} * \begin{pmatrix}2&3&5\end{pmatrix} = \texttt{[[2,6,15],[4,9,20]]} = \begin{pmatrix}2&6&15\\4&9&20\end{pmatrix}
and of course you will probably want to transpose that back
\texttt{(A.T * c).T} = \begin{pmatrix}2&4\\6&9\\15&20\end{pmatrix}
and just like before it doesn't matter the shape of c. Note that division is just the same
\texttt{(A.T / c).T} = \begin{pmatrix}0.5&1\\0.66...&1\\0.75&0.8\end{pmatrix}

Matrix multiplication of matrix and position vector

In most of the literature the positions are expanded out into a vector. For instance suppose we have two particles in 2d: one is (1,2), and the other is (2,3). then the position vector is given as
q = \begin{pmatrix}1\\2\\2\\3\end{pmatrix}
And if we have a matrix, G, as in the Rattle method, like so
G = \begin{pmatrix}1&2&2&3\\3&4&4&5\\5&6&6&7 \end{pmatrix}
and want to evaluate Gq we do this
\texttt{np.dot(G, q)} = \begin{pmatrix}1&2&2&3\\3&4&4&5\\5&6&6&7 \end{pmatrix} \begin{pmatrix}1\\2\\2\\3\end{pmatrix} = \texttt{array([18, 34, 50])} = \begin{pmatrix}18\\34\\50\end{pmatrix}

But this involves expanding out the positions into a vector, it is much nicer if they can be in a matrix. Then in ones code you can refer to position i as pos[i] instead of pos[3*i:3*i+3]. So now the position is given by

Q = \begin{pmatrix}1&2\\2&3\end{pmatrix}
This means the matrix G is changed so that each element of the matrix is actually a row vector. In other words G is a tensor, call it T. Eg
\texttt{T = np.array([[[1,2],[2,3]],[[3,4],[4,5]],[[5,6],[6,7]]])} = \begin{pmatrix} \begin{pmatrix}1&2\end{pmatrix} & \begin{pmatrix}2&3\end{pmatrix} \\ \begin{pmatrix}3&4\end{pmatrix} & \begin{pmatrix}4&5\end{pmatrix} \\ \begin{pmatrix}5&6\end{pmatrix} & \begin{pmatrix}6&7\end{pmatrix} \end{pmatrix}
Now to work out Gq, I would like to do a np.dot(T,Q) with a few transposes and so on, but it doesn't work. So you have to do a reshape, which gets us back to the original G above. I suppose the reshape in numpy is very fast so it doesn't matter much.
\texttt{np.dot(T.reshape((T.shape[0],-1)), Q.flatten()).T} = \begin{pmatrix}1&2&2&3\\3&4&4&5\\5&6&6&7\end{pmatrix} \begin{pmatrix}1\\2\\2\\3\end{pmatrix}

Another thing in the Rattle method is GM^{-1}G^T. Again if G and M are just regular matrices then you can use np.dot(np.dot(G, np.linalg.inv(M)), G.T)

\texttt{M = np.diag([3, 3, 7, 7])} = \begin{pmatrix}3&0&0&0\\0&3&0&0\\0&0&7&0\\0&0&0&7\end{pmatrix}
then GM^{-1}G^T is just
\texttt{np.dot(np.dot(G, np.linalg.inv(M)), G.T)} = \begin{pmatrix}3.52380952&6.95238095&10.38095238\\6.95238095&14.19047619&21.42857143\\10.38095238&21.42857143&32.4761904\\\end{pmatrix}
But if G is actually a tensor, T like above, and M is actually just a vector of masses you have to do something slightly different. Again we use the reshape trick. We can deal with the mass by a simple element-wise division, we don't have to do a full blown inverse multiplication.
\texttt{m = np.array([3, 7],dtype=np.float64)} = \begin{pmatrix}3&7\end{pmatrix}
then GM^{-1}G^T is just
\texttt{np.dot(T.reshape((T.shape[0],-1)), (T/m).reshape((T.shape[0],-1)).T)}
and the result is identical to above.