# Linear Algebra using SymPy


## Introduction

This notebook is a short tutorial of Linear Algebra calculation using SymPy. For further information refer to SymPy official [tutorial](http://docs.sympy.org/latest/tutorial/index.html).

You can also check the [SymPy in 10 minutes](./SymPy_in_10_minutes.ipynb) tutorial.


```python
from sympy import *
init_session()
```

    IPython console for SymPy 1.0 (Python 2.7.13-64-bit) (ground types: python)
    
    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()
    
    Documentation can be found at http://docs.sympy.org/1.0/


A matrix $A \in \mathbb{R}^{m\times n}$ is a rectangular array of real number with $m$ rows and $n$ columns. To specify  a matrix $A$, we specify the values for its components as a list of lists:


```python
A = Matrix([
        [3, 2, -1, 1],
        [2, -2, 4, -2],
        [-1, S(1)/2, -1, 0]])
display(A)
```


$$\left[\begin{matrix}3 & 2 & -1 & 1\\2 & -2 & 4 & -2\\-1 & \frac{1}{2} & -1 & 0\end{matrix}\right]$$


We can access the matrix elements using square brackets, we can also use it for submatrices


```python
A[0, 1] # row 0, column 1
```


```python
A[0:2, 0:3] # top-left 2x3 submatrix
```




$$\left[\begin{matrix}3 & 2 & -1\\2 & -2 & 4\end{matrix}\right]$$



We can also create some common matrices. Let us create an identity matrix


```python
eye(2)
```




$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$




```python
zeros(2, 3)
```




$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$



We can use algebraic operations like addition $+$, substraction $-$, multiplication $*$, and exponentiation $**$ with ``Matrix`` objects.


```python
B = Matrix([
        [2, -3, -8],
        [-2, -1, 2],
        [1, 0, -3]])
C = Matrix([
        [sin(x), exp(x**2), 1],
        [0, cos(x), 1/x],
        [1, 0, 2]])
```


```python
B + C
```




$$\left[\begin{matrix}\sin{\left (x \right )} + 2 & e^{x^{2}} - 3 & -7\\-2 & \cos{\left (x \right )} - 1 & 2 + \frac{1}{x}\\2 & 0 & -1\end{matrix}\right]$$




```python
B ** 2
```




$$\left[\begin{matrix}2 & -3 & 2\\0 & 7 & 8\\-1 & -3 & 1\end{matrix}\right]$$




```python
C ** 2
```




$$\left[\begin{matrix}\sin^{2}{\left (x \right )} + 1 & e^{x^{2}} \sin{\left (x \right )} + e^{x^{2}} \cos{\left (x \right )} & \sin{\left (x \right )} + 2 + \frac{e^{x^{2}}}{x}\\\frac{1}{x} & \cos^{2}{\left (x \right )} & \frac{1}{x} \cos{\left (x \right )} + \frac{2}{x}\\\sin{\left (x \right )} + 2 & e^{x^{2}} & 5\end{matrix}\right]$$




```python
tan(x) * B ** 5
```




$$\left[\begin{matrix}52 \tan{\left (x \right )} & 27 \tan{\left (x \right )} & - 28 \tan{\left (x \right )}\\- 2 \tan{\left (x \right )} & - \tan{\left (x \right )} & - 78 \tan{\left (x \right )}\\11 \tan{\left (x \right )} & 30 \tan{\left (x \right )} & 57 \tan{\left (x \right )}\end{matrix}\right]$$



And the ``transpose`` of the matrix, that flips the matrix through its main diagonal:


```python
A.transpose() # the same as A.T
```




$$\left[\begin{matrix}3 & 2 & -1\\2 & -2 & \frac{1}{2}\\-1 & 4 & -1\\1 & -2 & 0\end{matrix}\right]$$



## Row operations


```python
M = eye(4)
```


```python
M[1, :] = M[1, :] + 5*M[0, :]
```


```python
M
```




$$\left[\begin{matrix}1 & 0 & 0 & 0\\5 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$$



The notation ``M[1, :]`` refers to entire rows of the matrix. The first argument specifies the 0-based row index, for example the first row of ``M`` is ``M[0, :]``. The code example above implements the row operation $R_2 \leftarrow R_2 + 5R_1$. To scale a row by a constant $c$, use the ``M[1, :] = c*M[1, :]``. To swap rows $1$ and $j$, we can use the Python tuple-assignment syntax ``M[1, :], M[j, :] = M[j, :], M[1, :]``.

## Reduced row echelon form

The Gauss-Jordan elimination procedure is a sequence of row operations that can be performed on any matrix to bring it to its _reduced row echelon form_ (RREF). In Sympy, matrices have a ``rref`` method that compute it:


```python
A.rref()
```




$$\left ( \left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & -2\\0 & 0 & 1 & -2\end{matrix}\right], \quad \left [ 0, \quad 1, \quad 2\right ]\right )$$



It return a tuple, the first value is the RREF of the matrix $A$, and the second tells the location of the leading ones (pivots). If we just want the RREF, we can just get the first entry of the matrix, i.e.


```python
A.rref()[0]
```




$$\left[\begin{matrix}1 & 0 & 0 & 1\\0 & 1 & 0 & -2\\0 & 0 & 1 & -2\end{matrix}\right]$$



## Matrix fundamental spaces

Consider the matrix $A \in \mathbb{R}^{m\times n}$. The fundamental spaces of a matrix are its column space $\mathcal{C}(A)$, its null space $\mathcal{N}(A)$, and its row space $\mathcal{R}(A)$. These vector spaces are importan when we consider the matrix product $A\mathbf{x} = \mathbf{y}$ as a linear transformation $T_A:\mathbb{R}^n\rightarrow \mathbb{R}^n$ of the input vector $\mathbf{x}\in\mathbb{R}^n$ to produce an output vector $\mathbf{y} \in \mathbb{R}^m$.

**Linear transformations** $T_A: \mathbb{R}^n \rightarrow \mathbb{R}^m$ can be represented as $m\times n$ matrices. The fundamental spaces of a matrix $A$ gives us information about the domain and image of the linear transformation $T_A$. The column space $\mathcal{C}(A)$ is the same as the image space $\mathrm{Im}(T_A)$ (the set of all possible outputs). The null space $\mathcal{N}(A)$ is also called kernel $\mathrm{Ker}(T_A)$, and is the set of all input vectors that are mapped to the zero vector. The row space $\mathcal{R}(A)$ is the orthogonal complement of the null space, i.e., the vectors that are mapped to vectors different from zero. Input vectors in the row space of $A$ are in a one-to-one correspondence with the output vectors in the column space of $A$.

Let us see how to compute these spaces, or a base for them!

The non-zero rows in the reduced row echelon form $A$ are a basis for its row space, i.e.


```python
[A.rref()[0][row, :] for row in A.rref()[1]]
```




$$\left [ \left[\begin{matrix}1 & 0 & 0 & 1\end{matrix}\right], \quad \left[\begin{matrix}0 & 1 & 0 & -2\end{matrix}\right], \quad \left[\begin{matrix}0 & 0 & 1 & -2\end{matrix}\right]\right ]$$



The column space of $A$ is the span of the columns of $A$ that contain the pivots.


```python
[A[:, col] for col in A.rref()[1]]
```




$$\left [ \left[\begin{matrix}3\\2\\-1\end{matrix}\right], \quad \left[\begin{matrix}2\\-2\\\frac{1}{2}\end{matrix}\right], \quad \left[\begin{matrix}-1\\4\\-1\end{matrix}\right]\right ]$$



We can also use the ``columnspace`` method


```python
A.columnspace()
```




$$\left [ \left[\begin{matrix}3\\2\\-1\end{matrix}\right], \quad \left[\begin{matrix}2\\-2\\\frac{1}{2}\end{matrix}\right], \quad \left[\begin{matrix}-1\\4\\-1\end{matrix}\right]\right ]$$



Note that we took columns from the original matrix and not from its RREF.

To find (a base for) the null space of $A$ we use the ``nullspace`` method:


```python
A.nullspace()
```




$$\left [ \left[\begin{matrix}-1\\2\\2\\1\end{matrix}\right]\right ]$$



## Determinants

The determinant of a matrix, denoted by $\det(A)$ or $|A|$, isis a useful value that can be computed from the elements of a square matrix. It can be viewed as the scaling factor of the transformation described by the matrix.


```python
M = Matrix([
        [1, 2, 2],
        [4, 5, 6],
        [7, 8, 9]])
```


```python
M.det()
```

## Matrix inverse

For invertible matrices (those with $\det(A)\neq 0$), there is an inverse matrix $A^{-1}$ that have the _inverse_ effect (if we are thinking about linear transformations).


```python
A = Matrix([
        [1, -1, -1],
        [0, 1, 0],
        [1, -2, 1]])
```


```python
A.inv()
```




$$\left[\begin{matrix}\frac{1}{2} & \frac{3}{2} & \frac{1}{2}\\0 & 1 & 0\\- \frac{1}{2} & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]$$




```python
A.inv() * A
```




$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$




```python
A * A.inv()
```




$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$



## Eigenvectors and Eigenvalues

To find the eigenvalues of a matrix, use ``eigenvals``. ``eigenvals`` returns a dictionary of ``eigenvalue:algebraic multiplicity``.


```python
M = Matrix([
        [3, -2,  4, -2],
        [5,  3, -3, -2],
        [5, -2,  2, -2],
        [5, -2, -3,  3]])

M
```




$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$




```python
M.eigenvals()
```

This means that ``M`` has eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.

To find the eigenvectors of a matrix, use ``eigenvects``. ``eigenvects`` returns a list of tuples of the form ``(eigenvalue:algebraic multiplicity, [eigenvectors])``.


```python
M.eigenvects()
```




$$\left [ \left ( -2, \quad 1, \quad \left [ \left[\begin{matrix}0\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\\1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 5, \quad 2, \quad \left [ \left[\begin{matrix}1\\1\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\-1\\0\\1\end{matrix}\right]\right ]\right )\right ]$$



This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors. Because the algebraic and geometric multiplicities are the same for all the eigenvalues, ``M`` is diagonalizable.

To diagonalize a matrix, use diagonalize. diagonalize returns a tuple $(P,D)$, where $D$ is diagonal and $M=PDP^{−1}$.


```python
P, D = M.diagonalize()
```


```python
P
```




$$\left[\begin{matrix}0 & 1 & 1 & 0\\1 & 1 & 1 & -1\\1 & 1 & 1 & 0\\1 & 1 & 0 & 1\end{matrix}\right]$$




```python
D
```




$$\left[\begin{matrix}-2 & 0 & 0 & 0\\0 & 3 & 0 & 0\\0 & 0 & 5 & 0\\0 & 0 & 0 & 5\end{matrix}\right]$$




```python
P * D * P.inv()
```




$$\left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$$




```python
P * D * P.inv() == M
```




    True



Note that since ``eigenvects`` also includes the ``eigenvalues``, you should use it instead of ``eigenvals`` if you also want the ``eigenvectors``. However, as computing the eigenvectors may often be costly, ``eigenvals`` should be preferred if you only wish to find the eigenvalues.

If all you want is the characteristic polynomial, use ``charpoly``. This is more efficient than ``eigenvals``, because sometimes symbolic roots can be expensive to calculate.


```python
lamda = symbols('lamda')
p = M.charpoly(lamda)
factor(p)
```

**Note:** ``lambda`` is a reserved keyword in Python, so to create a Symbol called λ, while using the same names for SymPy Symbols and Python variables, use ``lamda`` (without the b). It will still pretty print as λ.

Non-square matrices don’t have eigenvectors and therefore don’t
have an eigendecomposition. Instead, we can use the singular value
decomposition to break up a non-square matrix A into left singular
vectors, right singular vectors, and a diagonal matrix of singular
values. Use the singular_values method on any matrix to find its
singular values.


```python
A
```




$$\left[\begin{matrix}1 & -1 & -1\\0 & 1 & 0\\1 & -2 & 1\end{matrix}\right]$$




```python
A.singular_values()
```

## References

1. SymPy Development Team (2016). [Sympy Tutorial: Matrices](http://docs.sympy.org/latest/tutorial/matrices.html)
2. Ivan Savov (2016). [Taming math and physics using SymPy](https://minireference.com/static/tutorials/sympy_tutorial.pdf)

The following cell change the style of the notebook.


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>








```python

```

503dab8963b7b367a390b5dbaecd26d88378b788
ipynb
Jupyter Notebook
notebooks/sympy/linear_algebra.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2017-06-29T17:45:20.000Z
2022-02-06T20:14:29.000Z
notebooks/sympy/linear_algebra.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
notebooks/sympy/linear_algebra.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2019-04-22T08:08:56.000Z
2022-01-27T08:15:53.000Z
Qwen/Qwen-72B
1. YES
2. YES