Doing the inversion directly turns out to be a bad idea since the matrix \( \boldsymbol{X}^T\boldsymbol{X} \) is singular. An alternative approach is to use the singular value decomposition. Using the definition of the Moore-Penrose pseudoinverse we can write the equation for \( \boldsymbol{\beta} \) as $$ \boldsymbol{\beta} = \boldsymbol{X}^{+}\boldsymbol{y}, $$
where the pseudoinverse of \( \boldsymbol{X} \) is given by $$ \boldsymbol{X}^{+} = \frac{\boldsymbol{X}^T}{\boldsymbol{X}^T\boldsymbol{X}}. $$
Using singular value decomposition we can decompose the matrix \( \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^T \), where \( \boldsymbol{U} \) and \( \boldsymbol{V} \) are orthogonal(unitary) matrices and \( \boldsymbol{\Sigma} \) contains the singular values (more details below). where \( X^{+} = V\Sigma^{+} U^T \). This reduces the equation for \( \omega \) to $$ \begin{align} \boldsymbol{\beta} = \boldsymbol{V}\boldsymbol{\Sigma}^{+} \boldsymbol{U}^T \boldsymbol{y}. \tag{26} \end{align} $$
Note that solving this equation by actually doing the pseudoinverse (which is what we will do) is not a good idea as this operation scales as \( \mathcal{O}(n^3) \), where \( n \) is the number of elements in a general matrix. Instead, doing \( QR \)-factorization and solving the linear system as an equation would reduce this down to \( \mathcal{O}(n^2) \) operations.
def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
u, s, v = scl.svd(x)
return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
beta = ols_svd(X_train_own,y_train)
When extracting the \( J \)-matrix we need to make sure that we remove the intercept, as is done here
J = beta[1:].reshape(L, L)
A way of looking at the coefficients in \( J \) is to plot the matrices as images.
fig = plt.figure(figsize=(20, 14))
im = plt.imshow(J, **cmap_args)
plt.title("OLS", fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
cb = fig.colorbar(im)
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
plt.show()
It is interesting to note that OLS considers both \( J_{j, j + 1} = -0.5 \) and \( J_{j, j - 1} = -0.5 \) as valid matrix elements for \( J \). In our discussion below on hyperparameters and Ridge and Lasso regression we will see that this problem can be removed, partly and only with Lasso regression.
In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?