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}

\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}

\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}

\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}

\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}

\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}

\texttt{(A.T * c).T} = \begin{pmatrix}2&4\\6&9\\15&20\end{pmatrix}

\texttt{(A.T / c).T} = \begin{pmatrix}0.5&1\\0.66...&1\\0.75&0.8\end{pmatrix}

q = \begin{pmatrix}1\\2\\2\\3\end{pmatrix}

G = \begin{pmatrix}1&2&2&3\\3&4&4&5\\5&6&6&7 \end{pmatrix}

\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}

\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}

\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}

\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}

\texttt{m = np.array([3, 7],dtype=np.float64)} = \begin{pmatrix}3&7\end{pmatrix}

\texttt{np.dot(T.reshape((T.shape[0],-1)), (T/m).reshape((T.shape[0],-1)).T)}