matrix_operations.qmd

#  Matrix Operations {#sec-matrix-operations}


Addition, subtraction, multiplication, division. These are all things you already know how to do with single numbers. What happens, though, if you want to multiply two different matrices together. Does that simple, 'scalar' operation still translate if you have a $2x3$ matrix and a $3x2$ matrix? If words like matrix and scalar make you break out in a sweat, then this chapter is for you!

Matrix operations, especially multiplication, are critical for understanding core aspects of how modeling actually produces all these cool results that help us discover so many interesting things. Knowing the underlying mechanics of matrix operations helps to demystify several issues that you might run into with your models. It can also help to get the gist of various articles and papers that you might come across. Before we get into any operations, though, let's make sure we are together on some concepts. 


```{r}
#| echo: false
#| label: matrix-vis

base = crossing(
    x = 1:3,
    y = 1:3
)  |> 
    mutate(
        scaler = ifelse(x == 2 & y == 2, 1, 0),
        row_vector = ifelse(y == 2, 1, 0),
        col_vector = ifelse(x == 2, 1, 0),
        matrix = 1
    )

p_scalar = base |> 
    ggplot(aes(x = 1, y = 1)) +
    geom_tile(
        aes(fill = factor(scaler)),
        width=0.8, 
        height=0.8,
        linewidth = 3,
        color = '#fffff8',
        na.rm = FALSE,
        show.legend = FALSE
    ) +
    geom_tile(
        aes(fill = factor(scaler)),
        width=0.7, 
        height=0.7,
        linewidth = 3,
        color = '#fffff8',
        na.rm = FALSE,
        show.legend = FALSE
    ) +
    scale_fill_manual(values = c('gray75', 'white')) +
    labs(
      # title ='Scalar',
    ) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = .5, size = 20)
    )


ggsave('img/app-matrix-scalar.svg', p_scalar, width = 4, height = 3, bg = 'transparent')

p_row = base |> 
    ggplot(aes(x = x, y = y)) +
    # geom_tile(
    #     aes(fill = factor(row_vector)),
    #     width=0.8, 
    #     height=0.8,
    #     linewidth = 3,
    #     color = '#fffff8',
    #     na.rm = FALSE,
    #     show.legend = FALSE
    # ) +
    geom_tile(
        aes(fill = factor(row_vector)),
        width = .7, 
        height = 1,
        linewidth = 3,
        color = '#fffff8',
        na.rm = FALSE,
        show.legend = FALSE
    ) +
    scale_fill_manual(values = c('white', 'gray75')) +
    labs(
      title = 'Row Vector',
    ) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = .5, size = 20)
    )

ggsave('img/app-matrix-row.svg', p_row, width = 8, height = 6, bg = 'transparent')

p_col = base |>
    ggplot(aes(x = x, y = y)) +
    # geom_tile(
    #     aes(fill = factor(col_vector)),
    #     width=0.8, 
    #     height=0.8,
    #     linewidth = 3,
    #     color = '#fffff8',
    #     na.rm = FALSE,
    #     show.legend = FALSE
    # ) +
    geom_tile(
        aes(fill = factor(col_vector)),
        width = 1, 
        height = .7,
        linewidth = 3,
        color = '#fffff8',
        na.rm = FALSE,
        show.legend = FALSE
    ) +
    scale_fill_manual(values = c('white', 'gray75')) +
    labs(
      title = 'Column Vector',
    ) +
    theme_void() +
    theme(
      plot.title = element_text(hjust = .5, size = 20)
    )

ggsave('img/app-matrix-col.svg', p_col, width = 8, height = 8)

# p_row + p_col

ggsave('img/app-matrix-row-col-vector.svg', p_row + p_col, width = 8, height = 6, bg = 'transparent')

p_matrix = base |> 
    ggplot(aes(x = x, y = y)) +
    # geom_tile(
    #     aes(fill = factor(matrix)),
    #     width = 1.05, 
    #     height = 1.05,
    #     linewidth = 3,
    #     color = '#fffff8',
    #     na.rm = FALSE,
    #     show.legend = FALSE
    # ) +
    geom_tile(
        aes(fill = factor(matrix)),
        width = .95, 
        height = .95,
        linewidth = 3,
        color = '#fffff8',
        na.rm = FALSE,
        show.legend = FALSE
    ) +
    scale_fill_manual(values = c( 'gray75')) +
    labs(
      # title = 'Matrix',
    ) +
    theme_void() +
    theme(
      # plot.background = element_rect(fill = 'gray98'),
      plot.title = element_text(hjust = .5, size = 20)
    )

ggsave('img/app-matrix-matrix.svg', p_matrix, width = 8, height = 6, bg = 'transparent')
```


A **scalar** is a single numeric value. It might help if you think about a scalar as a single 'block'. 

![Scalar](img/app-matrix-scalar.svg){width=20%}

:::{.panel-tabset}

##### Code
```{r}
#| eval: false
#| label: matrix-scalar
scalar_example = 1 # scalar value in r or python
```
:::


And just like we can line blocks up on the floor, we can put our scalars together to form a **vector**. A vector is a collection of scalars with a length of *n*. We can also think of a vector as a single row or column of scalars.


![Row and Column Vectors](img/app-matrix-row-col-vector.svg){width=50%}


There are many ways to create a vector in R and Python. Here are a couple.

:::{.panel-tabset}

##### R
```{r}
#| eval: false
#| label: matrix-vector-r
vector_example = 1:6
vector_example = c(1, 2, 3, 4, 5, 6)
vector_example = matrix(1:6, nrow = 1) # or ncol = 1
```

##### Python

```{python}
#| eval: false
#| label: matrix-vector-py
import numpy as np

vector_example = range(5)
vector_example = [1, 2, 3, 4, 5] # as list
row_vector = np.array([1, 2, 3])  # create a row vector
column_vector = np.array([[1], [2], [3]])  # create a column vector
```
:::

Now, we can take a few of our block vectors and stack them into a **matrix**, assuming the vectors are of the same size. A matrix is a 2 dimensional collection of vectors, and it is the fundamental structure for tabular data and beyond.

![Matrix](img/app-matrix-matrix.svg){width=50%}

And here is a matrix of specific values:

$$
\begin{bmatrix}
1 & 2 & 3\\
4 & 5 & 6\\
7 & 8 & 9
\end{bmatrix}
$$

If you think about most tables you've ever seen, you'll see that the simple matrix looks remarkably familiar!

:::{.panel-tabset}

##### R
```{r}
#| eval: false
#| label: matrix-matrix-r

matrix_example = matrix(1:6, nrow = 2, ncol = 3)
dim(matrix_example)
```

##### Python
```{python}
#| eval: false
#| label: matrix-matrix-py

matrix_example = np.array([[1, 2, 3], [4, 5, 6]])
matrix_example.shape
```
:::

A matrix has 2 dimensions, rows and columns, which can be any size. When we talk about the dimensions of a matrix, we always make note of the rows first, followed by the columns. This matrix has 2 rows and 3 columns, therefore, we have a $2x3$, 'two by three' matrix[^majormatrix].


[^majormatrix]: Numpy arrays/matrices are in **row major** order, while R is **column major** order. You'll note how with [numpy]{.pack} we essentially provided two rows to the array function, which automatically created the 2 x 3 matrix. The R matrix is not the same, because by default it fills in the columns. If you add `by_row = TRUE`, you'd then get the same result as the [numpy]{.pack} example. Column major is generally more intuitive for tabular data, because that's how we think of data stored in tables, and why the [pandas]{.pack} package in python is also column major/oriented. However, both R and Python are very flexible, and more generally work in arrays. If you use both it can take a bit to settle with one if you've used the other (especially for 'apply' functions). The [reticulate]{.pack} package has [a vignette](https://cran.r-project.org/web/packages/reticulate/vignettes/arrays.html) that provides a nice overview, while the [rray]{.pack} package actually brings the [numpy]{.pack} approach to the R landscape.



## Addition

Matrix addition, along with subtraction, is the easiest concept when dealing with matrices. There is one rule though: the matrices need to have the same dimensions. From a practical code perspective, if one is a scalar, addition of the scalar will be applied to every element in the matrix.

Let's check out these two matrices:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}  
\ 
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
$$

You probably noticed that we gave each scalar within the matrix a label associated with its row and column position. We can use these to see how we will produce the new matrix: 

Now, we can set this up as an addition problem to produce Matrix C:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}  
+ 
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
A_{11} + B_{11}& A_{12} + B_{12} & A_{13} + B_{13}\\
A_{21} + B_{21}& A_{22} + B_{22} & A_{23} + B_{23}
\end{bmatrix}
}
$$

Now we can pull in the real numbers:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}  
+ 
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
1 + 7  & 2 + 8 & 3 + 9\\
4 + 9 & 5 + 8 & 6 + 7
\end{bmatrix}
}
$$

Giving us Matrix C:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}  
+ 
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
8 & 10 & 12 \\
13 & 13 & 13
\end{bmatrix}
}
$$

First, let's create those matrices in R and Python.

:::{.panel-tabset}

##### R

In R, we can create a matrix a few ways: with the matrix function or by row binding numeric vectors.
```{r}
#| label: matrix-add-subtract-r-setup
matrix_A = rbind(1:3, 4:6)

# The following is an equivalent
# to rbind:
# matrix_A = matrix(
#     c(1:3, 4:6),
#     nrow = 2,
#     ncol = 3,
#     byrow = TRUE
# )

matrix_B = rbind(7:9, 9:7)
```


##### Python

The task is just as easy in Python. We will import `numpy` and then use the `matrix` method to create the matrices:

```{python}
#| label: matrix-add-py-setup
import numpy as np

matrix_A = np.array([[1, 2, 3], [4, 5, 6]])

matrix_B = np.array([[7, 8, 9], [9, 8, 7]])
```

:::


Once we have those matrices created, we can use the standard `+` to add them together:

:::{.panel-tabset}

##### R

```{r}
#| label: matrix-add-r
matrix_A + matrix_B
```

```{r}
#| label: matrix-add-scalar-r

matrix_A + 3
```

##### Python

Just like R, we can use `+` with those matrices.

```{python}
#| label: matrix-add-py
matrix_A + matrix_B
```

```{python}
#| label: matrix-add-scalar-py
matrix_A + 3
```

:::



## Subtraction

Take everything that you just saw with addition and replace it with subtraction! But just like addition, every matrix needs to have the same dimensions. 
<!-- 
Let's see those two matrices again and cast it as subtraction problem:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
-
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
A_{11} - B_{11}& A_{12} - B_{12} & A_{13} - B_{13}\\
A_{21} - B_{21}& A_{22} - B_{22} & A_{23} - B_{23}
\end{bmatrix}
}
$$

And now we can substitute in the real numbers:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
-
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
1 - 7 & 2 - 8 & 3 - 9\\
4 - 9 & 5 - 8 & 6 - 7
\end{bmatrix}
}
$$
-->

Here is the result:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
-
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
-6 & -6 & -6 \\
-5 & -3 & -1
\end{bmatrix}
}
$$

Subtracting matrices in R and Python is the same as addition, just using `-` instead.


:::{.panel-tabset}

##### R

```{r}
#| label: matrix-subtract-r
matrix_A - matrix_B
```

```{r}
#| label: matrix-subtract-scalar-r

matrix_A - 3
```


##### Python


```{python}
#| label: matrix-subtract-py
matrix_A - matrix_B
```

```{python}
#| label: matrix-subtract-scalar-py
matrix_A - 3
```

:::


## Transpose

You might see a matrix denoted as $A^T$ or $A'$. The superscripted $T$ or $'$ for matrix **transpose**. If we transpose a matrix, all we are doing is flipping the rows and columns along the matrix 'main diagonal'. A visual example is much easier:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
->
\stackrel{\mbox{Matrix A transposed} }{
\begin{bmatrix}
1 & 4 \\
2 & 5 \\
3 & 6
\end{bmatrix}
}
$$

:::{.panel-tabset}

##### R

In R, all we need is the `t` function:

```{r}
#| label: matrix-transpose-r
t(matrix_A)
```

##### Python

In Python, we can use numpy's `transpose` method:

```{python}
#| label: matrix-transpose-py
matrix_A.transpose() 
matrix_A.T  # shorthand
```

:::




## Multiplication

Now you probably have some confidence in doing matrix operations. Just as quickly as we built that confidence, it will be crushed when learning about matrix multiplication. 

When dealing with matrix multiplication, we have a huge change to our previous rule. No longer do our dimensions have to be the same! Instead, the matrices need to be *conformable* -- the first matrix needs to have the same number of columns as the number of rows within the second matrix. In other words, the inner dimensions must match. 

Look one more time at these matrices:

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
.
\stackrel{\mbox{Matrix B}}{
\begin{bmatrix}
7_{11} & 8_{12} & 9_{13}\\
9_{21} & 8_{22} & 7_{23}
\end{bmatrix} 
}
$$

Matrix A has dimensions of $2x3$, as does Matrix B. Putting those dimensions side by side -- $2x3 * 2x3$ -- we see that our inner dimensions are 3 and 2 and do not match. 

What if we *transpose* Matrix B?

$$
\stackrel{\mbox{Matrix B}^T}{
\begin{bmatrix}
7_{11} & 9_{12} \\ 
8_{21}& 8_{22}\\
9_{31} & 7_{32}
\end{bmatrix} 
}
$$

Now we have something that works!

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
.
\stackrel{\mbox{Matrix B}^T}{
\begin{bmatrix}
7_{11} & 9_{12} \\ 
8_{21}& 8_{22}\\
9_{31} & 7_{32}
\end{bmatrix} 
}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
. & . \\
. & . \\
\end{bmatrix}
}
$$

Now we have a $2x3 * 3x2$ matrix multiplication problem! The resulting matrix will have the same dimensions as our two matrices' outer dimensions: $2x2$.

Here is how we will get a $2x2$ matrix:

\small

$$
\stackrel{\mbox{Matrix A}}{
\begin{bmatrix}
1_{11} & 2_{12} & 3_{13}\\
4_{21} & 5_{22} & 6_{23}
\end{bmatrix}
}
.
\stackrel{\mbox{Matrix B}^T}{
\begin{bmatrix}
7_{11} & 9_{12} \\ 
8_{21}& 8_{22}\\
9_{31} & 7_{32}
\end{bmatrix} 
}
=
$$


$$
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
(A_{11}*B_{11})+(A_{12}*B_{21})+(A_{13}*B_{31}) & (A_{11}*B_{12})+(A_{12}*B_{22})+(A_{13}*B_{32})\\
(A_{21}*B_{11})+(A_{22}*B_{21})+(A_{23}*B_{31}) & (A_{21}*B_{12})+(A_{22}*B_{22})+(A_{23}*B_{32})
\end{bmatrix} 
}
$$

\normalsize

That might look like a horrible mess and likely isn't easy to commit to memory. Instead, we'd like to show you a way that might make it easier to remember how to multiply matrices. It also gives a nice representation of why your matrices need to be conformable.

We can leave Matrix A exactly where it is, flip Matrix B$^T$, and stack it right on top of Matrix A:

$$
\begin{bmatrix}
9_{b} & 8_{b} & 7_{b} \\
7_{b} & 8_{b} & 9_{b} \\
\\
1_{a} & 2_{a} & 3_{a} \\
4_{a} & 5_{a} & 6_{a}
\end{bmatrix}
$$

Now, we can let those rearranged columns from Matrix B$^T$ 'fall down' through the rows of Matrix A:

$$
\begin{bmatrix}
9_{b} & 8_{b} & 7_{b} \\
\\
1_{a}*7_{b} & 2_{a}*8_{b} & 3_{a}*9_{b}\\
4_{a} & 5_{a} & 6_{a}
\end{bmatrix}
= 
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
50 & .\\
. & .
\end{bmatrix}
}
$$

Adding those products together gives us 50 for $C_{11}$.

Let's move that row down to the next row in the Matrix A, multiply, and sum the products.

$$
\begin{bmatrix}
9_{b} & 8_{b} & 7_{b} \\
\\
1_{a} & 2_{a} & 3_{a}\\
4_{a}*7_{b} & 5_{a}*8_{b} & 6_{a}*9_{b}
\end{bmatrix}
= 
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
50 & .\\
122 & .
\end{bmatrix}
}
$$

We have 122 for $C_{21}$. That first column from Matrix B$^T$ won't be used any more, but now we need to move the second column through Matrix A.

$$
\begin{bmatrix}
1_{a}*9_{b} & 2_{a}*8_{b} & 3_{a}*7_{b}\\
4_{a} & 5_{a} & 6_{a}
\end{bmatrix}
= 
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
50 & 46\\
122 & .
\end{bmatrix}
}
$$

That gives us 46 for $C_{12}$. 

And finally:

$$
\begin{bmatrix}
1_{a} & 2_{a} & 3_{a}\\
4_{a}*9_{b} & 5_{a}*8_{b} & 6_{a}*7_{b}
\end{bmatrix}
=
\stackrel{\mbox{Matrix C}}{
\begin{bmatrix}
50 & 46\\
122 & 118
\end{bmatrix}
}
$$

We have 118 for $C_{22}$.

Now that you know how these work, you can see how easy it is to handle these tasks in R and Python.


:::{.panel-tabset}

##### R

In R, we need to use a fancy operator: `%*%`. This is just R's matrix multiplication operator. We will also use the transpose function: `t`. 

```{r}
#| label: matrix-multiply-r
matrix_A %*% t(matrix_B)
```

##### Python

In Python, we can just use the regular multiplication operator and the transpose method:

```{python}
#| label: matrix-multiply-py
matrix_A @ matrix_B.T
```

:::

You can see that whether we do this by hand, R, or Python, we come up with the same answer! While these small matrices can definitely be done by hand, we will always trust the computer to handle larger matrices. The main thing is to understand the mechanics behind the operation.

:::{.callout-note title='Elementwise Multiplication' collapse='true'}
Matrix multiplication is not the same as **elementwise** multiplication. Elementwise multiplication is when you multiply each element in one matrix by the corresponding element in another matrix. This is done with the `*` operator in R and Python. The matrices must have identical dimensions for this. As with addition and subtraction, if one matrix is a scalar, the operation is automatically applied to every element in the matrix.
:::


## Division

Though addition, subtraction, and multiplication are all pretty straightforward, matrix division is not. In fact, there really isn't such a thing as matrix division, we just use matrix multiplication in a particular way. This is similar to how we can divide two numbers, for example, $a/b$, but we can also multiply by the reciprocal, $a*(1/b)$. In matrix terms this would look something like:


$$
AB^{-1}
$$

While that may also seem straightforward on the surface, **matrix inversion** is not. The basic idea is that we are looking for a matrix that, when multiplied by the original matrix like $B$, gives us the identity matrix. The **identity matrix** is a matrix that has 1s along the diagonal and 0s everywhere else.

$$
\begin{bmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{bmatrix}
$$

Another caveat is that not all matrices have inverses. If the **determinant** of a matrix is 0, then it does not have an inverse. Technically, only square matrices can have inverses, but not all square matrices have inverses. We can, however, get a **pseudo-inverse** for non-square matrices.

:::{.panel-tabset}

##### R

```{r}
#| label: matrix-inverse-r

matrix_B_inv = MASS::ginv(matrix_B)
round(matrix_B %*% matrix_B_inv)
```

##### Python

```{python}
#| label: matrix-inverse-py

matrix_B_inv = np.linalg.pinv(matrix_B)
(matrix_B @ matrix_B_inv).round()
```

:::


More to the point, when would we do this? In the world of modeling, we might use matrix inversion to solve a system of equations. For example, this can be implemented in linear regression, where we are trying to find the coefficients that minimize the error in our model. That problem has an analytical solution that involves matrix inversion.

$$
\beta = (X^TX)^{-1}X^Ty
$$

Let's see this for ourselves. We will create a simple linear regression model and solve for the coefficients using matrix inversion.

:::{.panel-tabset}

##### R
```{r}
#| label: matrix-inverse-r-lm


set.seed(123)
x = rnorm(100)
y = 2*x + rnorm(100)
X = cbind(1, x)

beta = MASS::ginv(t(X) %*% X) %*% t(X) %*% y

tibble(
  ours = beta[,1],
  standard = coef(lm(y ~ x))
)
```


##### Python

```{python}
#| label: matrix

import statsmodels.api as sm
import pandas as pd

np.random.seed(123)
x = np.random.normal(size = 100)
y = 2*x + np.random.normal(size = 100)
X = np.c_[np.ones(100), x]

beta = np.linalg.pinv(X.T @ X) @ X.T @ y
beta

model_sm = sm.OLS(y, X)
results_sm = model_sm.fit()
coefficients_sm = results_sm.params

pd.DataFrame({
    'ours': beta,
    'standard': coefficients_sm
})
```

:::


## Summary

While matrix operations are not something we explicitly do everyday data science, it is always lurking behind the scenes. Having a grasp of the underlying model mechanics helps demystify the modeling process, and can greatly expand a data scientist's abilities when you have to dive into matrix operations for model building. Whether linear regression or deep learning, matrix operations are at the core of almost every model you come across.