{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"hide_cell": true
},
"source": [
"Make me look good. Click on the cell below and press Ctrl+Enter."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide_cell": true
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"HTML(open('css/custom.css', 'r').read())"
]
},
{
"cell_type": "markdown",
"metadata": {
"hide_cell": true
},
"source": [
"SM286D · Introduction to Applied Mathematics with Python · Spring 2020 · Uhan
\n",
"\n",
"Lesson 11.
\n",
"\n",
"Matrix algebra with Python
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## This lesson..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Matrix algebra with NumPy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What is NumPy?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- __NumPy__ is the fundamental scientific computing package for Python.\n",
"\n",
"- Numerous scientific and mathematical packages for Python use NumPy, including __pandas__ (data analysis and wrangling) and __scikit-learn__ (machine learning).\n",
"\n",
"- Let's start by importing NumPy:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## NumPy arrays"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- The main data structure in NumPy is the __array__. Arrays are typically used to represent matrices.\n",
"\n",
"- We can create a one-dimensional array like this:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[2 3 4]\n",
"\n"
]
}
],
"source": [
"# Create one-dimensional array\n",
"a = np.array([2, 3, 4])\n",
"\n",
"# Print the array\n",
"print(a)\n",
"\n",
"# Check the type of the array\n",
"print(type(a))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Similarly, we can create a two-dimensional array like this:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1 2]\n",
" [3 4]]\n",
"\n"
]
}
],
"source": [
"# Create two-dimensional array\n",
"b = np.array([[1, 2], [3, 4]])\n",
"\n",
"# Print the array\n",
"print(b)\n",
"\n",
"# Check the type of the array\n",
"print(type(b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- We can access and change values of an array like this:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3\n",
"[2 3 0]\n"
]
}
],
"source": [
"# Print element of a\n",
"print(a[1])\n",
"\n",
"# Change element of a\n",
"a[2] = 0\n",
"print(a)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3\n",
"[[ 1 100]\n",
" [ 3 4]]\n"
]
}
],
"source": [
"# Print element of b\n",
"print(b[1, 0])\n",
"\n",
"# Change element of b\n",
"b[0, 1] = 100\n",
"print(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Some basic Numpy array methods and attributes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- NumPy arrays have a number of useful methods and attributes.\n",
"\n",
"- Learn about some of the more common of these methods by answering the questions below.\n",
"\n",
"- If you're stuck or have doubts, [try searching the documentation.](https://docs.scipy.org/doc/numpy/index.html) \n",
" - Wading through the documentation can be hard at first but you'll get better at reading documentation with practice. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 1\n",
"\n",
"Run the code below. What does `np.arange` do? "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1 2 3]\n"
]
}
],
"source": [
"c = np.arange(1, 4)\n",
"print(c)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `np.arange(start, stop)` produces a one-dimensional array, with values starting at `start`, increasing by 1, and ending at `stop - 1`.\n",
"\n",
"- [Documentation for `np.arange`](https://numpy.org/devdocs/reference/generated/numpy.arange.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 2\n",
"\n",
"Run the code below. What does `np.eye` do?"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0. 0. 0.]\n",
" [0. 1. 0. 0.]\n",
" [0. 0. 1. 0.]\n",
" [0. 0. 0. 1.]]\n"
]
}
],
"source": [
"I = np.eye(4)\n",
"print(I)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `np.eye(n)` produces a two-dimensional array with `n` rows and columns, with ones on the diagonal and zeros elsewhere (in other words, the identity matrix of size `n`).\n",
"\n",
"- [Documentation for `np.eye`](https://numpy.org/devdocs/reference/generated/numpy.eye.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 3\n",
"\n",
"Run the code below. What does `np.zeros()` do?"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]]\n"
]
}
],
"source": [
"d = np.zeros((8,2))\n",
"print(d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `np.zeros((a, b))` produces a two-dimensional array of zeros with `a` rows and `b` columns.\n",
"\n",
"- [Documentation for `np.zeros`](https://numpy.org/devdocs/reference/generated/numpy.zeros.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 4\n",
"\n",
"Run the code below. `.transpose()` is a NumPy array method. What does `.transpose()` do? Does `.transpose()` change the array permanently?"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0. 0. 0. 0. 0. 0. 0. 0.]\n",
" [0. 0. 0. 0. 0. 0. 0. 0.]]\n",
"[[0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]\n",
" [0. 0.]]\n"
]
}
],
"source": [
"print(d.transpose())\n",
"print(d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `.transpose()` takes the transpose of the associated array. It does not change the array permanently.\n",
"\n",
"- [Documentation for `.transpose()`](https://numpy.org/devdocs/reference/generated/numpy.ndarray.transpose.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 5\n",
"\n",
"Run the code below. `.shape` is a NumPy array attribute. What information does `.shape` give?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(8, 2)\n",
"(2, 8)\n"
]
}
],
"source": [
"print(d.shape)\n",
"print(d.transpose().shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- For a two-dimensional array, `.shape` contains a tuple `(a, b)` where `a` is the number of rows, and `b` is the number of columns.\n",
"\n",
"- [Documentation for `.shape`](https://numpy.org/devdocs/reference/generated/numpy.ndarray.shape.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 6\n",
"\n",
"Run the code below. `.reshape()` is a NumPy array method. What does the `.reshape()` method do? Does `.reshape()` change the array permanently?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]\n",
"[[ 1 2 3 4 5 6 7 8]\n",
" [ 9 10 11 12 13 14 15 16]]\n",
"[[ 1 2 3 4]\n",
" [ 5 6 7 8]\n",
" [ 9 10 11 12]\n",
" [13 14 15 16]]\n"
]
}
],
"source": [
"e = np.arange(1, 17)\n",
"print(e)\n",
"\n",
"f = e.reshape(2, 8)\n",
"print(f)\n",
"\n",
"g = e.reshape(4, 4)\n",
"print(g)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `.reshape(a, b)` takes the entries of the associated array, and puts them back, row-by-row, into an array with `a` rows and `b` columns.\n",
"\n",
"- [Documentation for `.reshape`](https://numpy.org/devdocs/reference/generated/numpy.ndarray.reshape.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 7\n",
"\n",
"Run the code below. What does `np.linspace()` do?"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 1. 3.25 5.5 7.75 10. ]\n"
]
}
],
"source": [
"h = np.linspace(1, 10, 5)\n",
"print(h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `np.linspace(start, stop, num)` returns an array of `num` evenly-spaced values between `start` and `stop`.\n",
"\n",
"- [Documentation for `np.linspace`](https://numpy.org/devdocs/reference/generated/numpy.linspace.html)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Matrix multiplication with NumPy arrays"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- As you might expect, you can perform various matrix algebra tasks using NumPy arrays.\n",
"\n",
"- Let's start with matrix multiplication."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Let\n",
"\n",
"$$\n",
"A = \\begin{bmatrix}\n",
"1 & 2 & 3\\\\\n",
"4 & 5 & 6\n",
"\\end{bmatrix}\n",
"\\qquad\n",
"B = \\begin{bmatrix}\n",
"2\\\\\n",
"0\\\\\n",
"1\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"\n",
"- Let's define these matrices as NumPy arrays:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1 2 3]\n",
" [4 5 6]]\n",
"[[2]\n",
" [0]\n",
" [1]]\n"
]
}
],
"source": [
"A = np.array([[1, 2, 3], [4, 5, 6]])\n",
"B = np.array([[2], [0], [1]])\n",
"\n",
"print(A)\n",
"print(B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Here's one way to get $AB$:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 5]\n",
" [14]]\n"
]
}
],
"source": [
"print(A @ B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Here's another way to get $AB$:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 5]\n",
" [14]]\n"
]
}
],
"source": [
"print(A.dot(B))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Here's yet another way to get $AB$:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 5]\n",
" [14]]\n"
]
}
],
"source": [
"print(np.matmul(A, B))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Here's one way __not__ to get $AB$:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "operands could not be broadcast together with shapes (2,3) (3,1) ",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mB\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (2,3) (3,1) "
]
}
],
"source": [
"print(A * B)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Note that `*` does __not__ perform matrix multiplication!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Continue your exploration of how to use NumPy arrays by answering the questions below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 8\n",
"\n",
"Run the code cell below. What does `*` do to an array when used with a scalar value?"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[2. 0. 0.]\n",
" [0. 2. 0.]\n",
" [0. 0. 2.]]\n"
]
}
],
"source": [
"twoi = 2 * np.eye(3)\n",
"print(twoi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- `*` acts as a scalar multiplication operator when used with a scalar and an array."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 9\n",
"\n",
"Run the code cell below. How do you add arrays with the same dimensions?"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[3. 0. 0.]\n",
" [0. 3. 0.]\n",
" [0. 0. 3.]]\n"
]
}
],
"source": [
"threei = twoi + np.eye(3)\n",
"print(threei)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Write your answer here. Double-click to edit.__\n",
"\n",
"- The `+` operator adds two arrays together."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Determinants and inverses"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Next, let's explore how to compute determinants and inverses with NumPy arrays."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Let $C = \\begin{bmatrix} 1 & 2 \\\\ 4 & 5 \\end{bmatrix}$:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"C = np.array([[1, 2], [4, 5]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- We can find $\\det(C)$ like this:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-2.9999999999999996\n"
]
}
],
"source": [
"print(np.linalg.det(C))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- We can find $C^{-1}$ like this:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-1.66666667 0.66666667]\n",
" [ 1.33333333 -0.33333333]]\n"
]
}
],
"source": [
"print(np.linalg.inv(C))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 10\n",
"\n",
"Write code to compute $C C^{-1}$ using NumPy. Before you run your code, what do you expect to get? Now run your code. Did you get what you expected?"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1., 0.],\n",
" [0., 1.]])"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C @ np.linalg.inv(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving a system of linear equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Consider the system $Cx = d$, where\n",
"\n",
"$$\n",
"C = \\begin{bmatrix} 1 & 2 \\\\ 4 & 5 \\end{bmatrix} \\qquad d = \\begin{bmatrix} 5 \\\\ 14 \\end{bmatrix}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 11\n",
"\n",
"Solve for $x$ by computing $C^{-1} d$."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1.]\n",
" [2.]]\n"
]
}
],
"source": [
"d = np.array([[5], [14]])\n",
"x = np.linalg.inv(C) @ d\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- You can also use `np.linalg.solve` like this:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1.]\n",
" [2.]]\n"
]
}
],
"source": [
"# Solve Cx = d\n",
"x = np.linalg.solve(C, d)\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Eigenvalues and eigenvectors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Next, let's compute the eigenvalues and eigenvectors of $C$:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The eigenvalues of C are [-0.46410162 6.46410162].\n",
"\n",
"The eigenvectors of C are the column vectors of\n",
"[[-0.80689822 -0.34372377]\n",
" [ 0.59069049 -0.9390708 ]]\n"
]
}
],
"source": [
"# Note that np.linalg.eig returns 2 objects\n",
"w, V = np.linalg.eig(C)\n",
"\n",
"print(f\"The eigenvalues of C are {w}.\\n\")\n",
"print(f\"The eigenvectors of C are the column vectors of\")\n",
"print(V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Note that `np.linalg.eig` returns a tuple of 2 objects: `w` and `V`\n",
" - The first object `w` is a one-dimensional array of the eigenvalues. The eigenvalues are repeated according to their multiplicity.\n",
" - The second object `V` is a two-dimensional array of the eigenvectors: the $i$th column of `V` corresponds to the $i$th eigenvalue in `w`.\n",
"\n",
"- [Here is the documentation for `numpy.linalg.eig`.](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html?highlight=eig)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- To get a column of a two-dimensional array, we can slice the array like this:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The eigenvector corresponding to -0.4641016151377544 is\n",
"[-0.80689822 0.59069049]\n",
"\n",
"The eigenvector corresponding to 6.464101615137754 is\n",
"[-0.34372377 -0.9390708 ]\n"
]
}
],
"source": [
"# Get the first column of the eigenvector matrix\n",
"v0 = V[:, 0]\n",
"print(f\"The eigenvector corresponding to {w[0]} is\")\n",
"print(v0)\n",
"\n",
"# Get the second column of the eigenvector matrix\n",
"v1 = V[:, 1]\n",
"print(f\"\\nThe eigenvector corresponding to {w[1]} is\")\n",
"print(v1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 12\n",
"\n",
"Check if `v0` defined in the above code cell is in fact an eigenvector of $C$: compute $C v_0$ and $w_0 v_0$, and check if they are equal."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"C v0 is\n",
"[ 0.37448277 -0.27414041]\n",
"\n",
"w0 v0 is\n",
"[ 0.37448277 -0.27414041]\n"
]
}
],
"source": [
"print(\"C v0 is\")\n",
"print(C @ v0)\n",
"\n",
"print(\"\\nw0 v0 is\")\n",
"print(w[0] * v0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 13\n",
"\n",
"Find the eigenvalues and eigenvectors of the matrix $\\displaystyle{\\begin{bmatrix} 0 & 0 & 2 \\\\ -2 & 1 & 4 \\\\ -1 & 0 & 3 \\end{bmatrix}}$. Can you find eigenvectors with integer entries? "
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The eigenvalues of this matrix are\n",
"[1. 1. 2.]\n",
"\n",
"The eigenvectors of this matrix are the column vectors of\n",
"[[ 0. -0.89442719 -0.40824829]\n",
" [ 1. 0. -0.81649658]\n",
" [ 0. -0.4472136 -0.40824829]]\n"
]
}
],
"source": [
"# Define the matrix as a NumPy array\n",
"E = np.array([[0, 0, 2], [-2, 1, 4], [-1, 0, 3]])\n",
"\n",
"# Find the eigenvalues and eigenvectors of E\n",
"w, V = np.linalg.eig(E)\n",
"\n",
"# Print the eigenvalues\n",
"print(\"The eigenvalues of this matrix are\")\n",
"print(w)\n",
"\n",
"# Print the eigenvectors\n",
"print(\"\\nThe eigenvectors of this matrix are the column vectors of\")\n",
"print(V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- The next code cell factors the matrix $C$ into a product of the form $V D V^{-1}$, where\n",
" - $V$ is a matrix whose columns are the eigenvectors $v_0$ and $v_1$ of $C$ (i.e. the second output from `np.linalg.eig(C)`)\n",
" - $D$ is a diagonal matrix whose entries are the eigenvalues $w_0$ and $w_1$ of $C$ (i.e. the first output from `np.linalg.eig(C)`)\n",
"\n",
"- In other words, we can factor $C$ like this:\n",
"\\begin{equation*}\n",
"C = \\begin{bmatrix} \\mid & \\mid \\\\ v_0 & v_1 \\\\ \\mid & \\mid \\end{bmatrix}\\; \\begin{bmatrix} w_0 & 0 \\\\ 0 & w_1 \\end{bmatrix} \\; \n",
"\\begin{bmatrix} \\mid & \\mid \\\\ v_0 & v_1 \\\\ \\mid & \\mid \\end{bmatrix}^{-1}.\n",
"\\end{equation*}"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-0.46410162 0. ]\n",
" [ 0. 6.46410162]]\n",
"[[1. 2.]\n",
" [4. 5.]]\n",
"[[1 2]\n",
" [4 5]]\n"
]
}
],
"source": [
"# Find the eigenvalues and eigenvectors of C\n",
"w, V = np.linalg.eig(C)\n",
"\n",
"# np.diag takes a one-dimensional array and makes it into a diagonal matrix\n",
"D = np.diag(w)\n",
"print(D)\n",
"\n",
"# C should be equal to V * D * V^-1\n",
"print(V @ D @ np.linalg.inv(V))\n",
"print(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 14\n",
"\n",
"Factor the matrix from Question 13 as a product of the form $VDV^{-1}$, where $D$ is a diagonal matrix."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0. 0.]\n",
" [0. 1. 0.]\n",
" [0. 0. 2.]]\n",
"[[ 0. 0. 2.]\n",
" [-2. 1. 4.]\n",
" [-1. 0. 3.]]\n",
"[[ 0 0 2]\n",
" [-2 1 4]\n",
" [-1 0 3]]\n"
]
}
],
"source": [
"# Find the eigenvalues and eigenvectors of E, the matrix from Question\n",
"w, V = np.linalg.eig(E)\n",
"\n",
"# np.diag takes a one-dimensional array and makes it into a diagonal matrix\n",
"D = np.diag(w)\n",
"print(D)\n",
"\n",
"# C should be equal to V * D * V^-1\n",
"print(V @ D @ np.linalg.inv(V))\n",
"print(E)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Where is the RREF?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__WARNING. The material in this section is much more advanced than the rest of the material in this lesson. We recommend skipping this section on a first pass, only coming back to it if you have extra time.__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now you might be interested in finding the row reduced echelon form (RREF) for a matrix; after all, that was the main tool in SM261 and we used it extensively to perform all sorts of computations. However, there is no such command in NumPy and there is a good reason for its absence. It turns out that row reduction is numerically unstable: small changes in the entries of a matrix can cause qualitative changes in its row reduction. \n",
"\n",
"Fortunately, almost every sensible application of row reduction can be accomplished by a numerically stable operation called the __singular value decomposition (SVD)__. The SVD of a $m \\times n$ matrix $C$ is a factorization, \n",
"\n",
"$$\n",
"C = U\\Sigma V^T,\n",
"$$ \n",
"\n",
"where $U$ and $V$ are orthonormal matrices (the dot product of the columns of these matrices are zero if the columns are different and 1 if the columns are equal) of sizes $m \\times m$ and $n \\times n$, respectively, and $\\Sigma$ is a block matrix whose upper left corner is a square diagonal matrix with nonzero entries on the diagonal and all other blocks are zero. The number of nonzero entries of $\\Sigma$ is the rank $r$ of $C$. To be explicit, \n",
"\n",
"$$\n",
"\\Sigma = \\begin{bmatrix} \\Sigma & 0 \\\\ 0 & 0 \\end{bmatrix},\n",
"$$ \n",
"\n",
"where the matrix $\\Sigma$ is an $r\\times r$ diagonal matrix, the zero to its right is a $r \\times (n-r)$ matrix, and the zero in the bottom corner is an $(m-r) \\times (n-r)$ matrix. \n",
"\n",
"As an example, consider the matrix \n",
"\n",
"$$\n",
"C = \\begin{bmatrix} \\phantom{-}1 & -1 \\\\ -2 & \\phantom{-}2 \\\\ \\phantom{-}2 & -2 \\end{bmatrix}.\n",
"$$ \n",
"\n",
"The matrix $C$ has rank 1 and has SVD given by:\n",
"\n",
"\\begin{equation} \n",
"\\label{svdeq} \n",
"C = \\begin{bmatrix} \n",
"-1/3 & \\phantom{-}2/3 & -2/3 \\\\\n",
"\\phantom{-}2/3 & \\phantom{-}2/3 & \\phantom{-}1/3 \\\\ \n",
"-2/3 & \\phantom{-}1/3 & \\phantom{-}2/3\n",
"\\end{bmatrix} \\, \n",
"\\begin{bmatrix} \n",
"3 \\sqrt{2} & 0 \\\\ \n",
"0 & 0 \\\\ \n",
"0 & 0\n",
"\\end{bmatrix} \\, \n",
"\\begin{bmatrix} \n",
"-1/\\sqrt{2} & \\phantom{-}1/\\sqrt{2} \\\\ \n",
"\\phantom{-}1/\\sqrt{2} & \\phantom{-}1/\\sqrt{2} \n",
"\\end{bmatrix}.\n",
"\\end{equation} \n",
"\n",
"Note that $C$ is a $3 \\times 2$ matrix, $U$ is $3 \\times 3$ and $V$ is $2 \\times 2$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can recover four critical subspaces of $C$ by writing the SVD factorization in the following form\n",
"\n",
"$$\n",
"C = \\begin{bmatrix} U_r & U_{r,\\perp} \\end{bmatrix} \\begin{bmatrix} S_r & 0 \\\\ 0 & 0 \\end{bmatrix} \\begin{bmatrix} V_r^T \\\\ V_{r,\\perp}^T \\end{bmatrix},\n",
"$$\n",
"\n",
"where $S_r$ is an $r \\times r$ diagonal matrix with nonzero entries on the diagonal. Then \n",
"\n",
"$$ \n",
"\\begin{array}{rll} \\text{range}(C) & = & \\text{range}(U_r), \\\\\n",
"\\text{null}(C) & = & \\text{range}(V_{r,\\perp}), \\\\\n",
"\\text{range}(C^T) & = & \\text{range}(V_r), \\\\\n",
"\\text{null}(C^T) & = & \\text{range}(U_{r,\\perp}). \\end{array}\n",
"$$ \n",
"\n",
"Recalling that the range of a matrix is the span of its columns, this gives us a way to describe each of the four key subspaces associated to a matrix $C$. Indeed, since the matrices $U$ and $V$ are orthonormal, their columns and rows form a basis for the space that they span. \n",
"\n",
"Let's see how some of these ideas work for our SVD of $C$. There is just one nonzero element in $\\Sigma$, so $C$ must have rank 1. Then $U_r$ should be $3 \\times 1$ ($U_r$ has $r$ columns). So a basis for the range of $C$ is the set consisting of the first column vector in $U$: \n",
"\n",
"$$\n",
"\\displaystyle{\\left\\{\\begin{bmatrix} -1/3 \\\\ 2/3 \\\\ -2/3 \\end{bmatrix} \\right\\}}.\n",
"$$\n",
"\n",
"Check that this really is the range of $C$! The null space of $C$ is the range of the matrix $V_{r,\\perp}$. The matrix $V_r^T$ is an $r \\times n$ matrix (has $r$ __rows__), meaning that the matrix $V_{r,\\perp}^T$ has size $(n-r)$ rows. That means that $V_{r,\\perp}$ has $n-r$ columns, each the transpose of one of the last $n-r$ rows of $V_{r,\\perp}^T$. Now \n",
"\n",
"$$\n",
"\\begin{bmatrix} \n",
"V_{r}^T \\\\ \n",
"V_{r,\\perp}^T \n",
"\\end{bmatrix} \n",
"= \n",
"\\begin{bmatrix} \n",
"-1/\\sqrt{2} & 1/\\sqrt{2} \\\\ \n",
"\\phantom{-}1/\\sqrt{2} & 1/\\sqrt{2} \n",
"\\end{bmatrix}\n",
"$$ \n",
"\n",
"so $V_{r,\\perp}^T = \\begin{bmatrix} 1/\\sqrt{2} & \\phantom{-}1/\\sqrt{2} \\end{bmatrix}$ and $V_{r,\\perp} = \\begin{bmatrix} 1/\\sqrt{2} \\\\ 1/\\sqrt{2} \\end{bmatrix}$. A basis for the null space of $C$ is \n",
"\n",
"$$\n",
"\\left\\{ \\begin{bmatrix} 1/\\sqrt{2} \\\\ 1/\\sqrt{2} \\end{bmatrix} \\right\\}.\n",
"$$\n",
"\n",
"Check that this is true! \n",
"\n",
"In the next code cell we compute the Singular Value Decomposition (SVD) of the matrix $C$. Note that on line 14 we learn that `vh`, the third output of `np.linalg.svd`, is already transposed: `vh` = $V^T$."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix U is\n",
"[[-0.33333333 0.66666667 -0.66666667]\n",
" [ 0.66666667 0.66666667 0.33333333]\n",
" [-0.66666667 0.33333333 0.66666667]]\n",
"\n",
"The matrix Sigma is\n",
"[[4.24264069 0. ]\n",
" [0. 0. ]\n",
" [0. 0. ]]\n",
"\n",
"The matrix V^T is\n",
"[[-0.70710678 0.70710678]\n",
" [ 0.70710678 0.70710678]]\n",
"\n",
"Multiplying these together, we should get back C:\n",
"[[ 1. -1.]\n",
" [-2. 2.]\n",
" [ 2. -2.]]\n"
]
}
],
"source": [
"# Define C\n",
"C = np.array([[1, -1], [-2, 2], [2, -2]])\n",
"\n",
"# Get dimensions of C\n",
"m = C.shape[0]\n",
"n = C.shape[1]\n",
"\n",
"# Find the SVD of C\n",
"u, s, vh = np.linalg.svd(C)\n",
"\n",
"print(\"The matrix U is\")\n",
"print(u)\n",
"\n",
"# s returned by np.linalg.svd is a one-dimensional array containing the diagonal entries of sigma\n",
"n_s = s.shape[0] # number of diagonal entries\n",
"sigma = np.zeros((m, n)) # initialize sigma to be a matrix of zeros\n",
"sigma[:n_s,:n_s] = np.diag(s) # fill entries of sigma\n",
"\n",
"print(\"\\nThe matrix Sigma is\")\n",
"print(sigma)\n",
"\n",
"print(\"\\nThe matrix V^T is\")\n",
"print(vh)\n",
"\n",
"# Check the factorization: this should equal C\n",
"print(\"\\nMultiplying these together, we should get back C:\")\n",
"print(u @ sigma @ vh)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 15\n",
"\n",
"(a) Compute the SVD of the matrix \n",
"\n",
"$$\n",
"Q = \\begin{bmatrix} 2 & 1 & 0 &1 & 0 \\\\\n",
" 1 & 0 & 1 & 1 & 1 \\\\\n",
" -1 & 1 & 1 & 0 & 1 \\\\\n",
" 2 & 2 & 2 & 2 & 2 \\end{bmatrix}.\n",
"$$\n",
"\n",
"(b) Verify that the rank $r$ of the matrix $Q$ is 3 using its SVD. \n",
"\n",
"(c) Find the shapes of the matrices $U$, $S$ and $V^T$.\n",
"\n",
"__The next three parts of this question are considered challenge material.__\n",
"\n",
"(d) The first $r$ columns of $U$ make up $U_r$ and the remaining columns make up $U_{r,\\perp}$. The first $r$ rows of $V^T$ make up $V_{r}^T$ and the remaining rows make up $V_{r,\\perp}^T$. Find these submatrices.\n",
"\n",
"(e) Find a basis for the null space of the original matrix $Q$. Does this basis have the correct number of vectors? Do your vectors live in the correct ambient space? \n",
"\n",
"(f) Find a basis for the range of the original matrix $Q$. Does this basis have the correct number of vectors? Do your vectors live in the correct ambient space?"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix U is\n",
"[[-0.36206122 0.6246702 0.47822465 -0.5 ]\n",
" [-0.34561768 0.00735298 -0.79403674 -0.5 ]\n",
" [-0.14568802 -0.7687844 0.37114088 -0.5 ]\n",
" [-0.85336692 -0.13676122 0.05532879 0.5 ]]\n",
"\n",
"The matrix Sigma is\n",
"[[5.22657769e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
" 0.00000000e+00]\n",
" [0.00000000e+00 2.34830781e+00 0.00000000e+00 0.00000000e+00\n",
" 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 1.08089598e+00 0.00000000e+00\n",
" 0.00000000e+00]\n",
" [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.32285481e-16\n",
" 0.00000000e+00]]\n",
"\n",
"The matrix V^T is\n",
"[[-5.03347715e-01 -4.23696576e-01 -4.20550437e-01 -4.61949076e-01\n",
" -4.20550437e-01]\n",
" [ 7.46050124e-01 -1.77845780e-01 -4.40723250e-01 1.52663437e-01\n",
" -4.40723250e-01]\n",
" [-9.07309675e-02 8.88173457e-01 -2.88869869e-01 -1.89800418e-01\n",
" -2.88869869e-01]\n",
" [ 3.45876835e-02 6.24500451e-17 7.22070511e-01 -6.91753671e-02\n",
" -6.87482827e-01]\n",
" [ 4.24996322e-01 0.00000000e+00 1.55140977e-01 -8.49992645e-01\n",
" 2.69855345e-01]]\n",
"\n",
"We recover the matrix Q\n",
"[[ 2.00000000e+00 1.00000000e+00 9.10181558e-17 1.00000000e+00\n",
" 2.22044605e-16]\n",
" [ 1.00000000e+00 1.11022302e-15 1.00000000e+00 1.00000000e+00\n",
" 1.00000000e+00]\n",
" [-1.00000000e+00 1.00000000e+00 1.00000000e+00 3.09886780e-16\n",
" 1.00000000e+00]\n",
" [ 2.00000000e+00 2.00000000e+00 2.00000000e+00 2.00000000e+00\n",
" 2.00000000e+00]]\n"
]
}
],
"source": [
"# Define the matrix\n",
"Q = np.array([[2, 1, 0, 1, 0], [1, 0, 1, 1, 1], [-1, 1, 1, 0, 1], [2, 2, 2, 2, 2]])\n",
"\n",
"# Get dimensions of Q\n",
"m = Q.shape[0]\n",
"n = Q.shape[1]\n",
"\n",
"# (a) Find SVD of Q\n",
"u, s, vh = np.linalg.svd(Q)\n",
"\n",
"print(\"The matrix U is\")\n",
"print(u)\n",
"\n",
"n_s = s.shape[0]\n",
"sigma = np.zeros((m, n))\n",
"sigma[:n_s,:n_s] = np.diag(s)\n",
"\n",
"print(\"\\nThe matrix Sigma is\")\n",
"print(sigma)\n",
"\n",
"print(\"\\nThe matrix V^T is\")\n",
"print(vh)\n",
"\n",
"print(\"\\nWe recover the matrix Q\")\n",
"print(u @ sigma @ vh)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) The matrix $\\Sigma$ (the array `s`) has 3 nonzero entries, so the rank of $Q$ is 3. (Note the last entry of `s` is on the order of magnitude of $10^{-16}$: it is effectively zero.)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Shape of U: (4, 4)\n",
"Shape of Sigma: (4, 5)\n",
"Shape of V^T: (5, 5)\n"
]
}
],
"source": [
"# (c) Shapes of SVD factors\n",
"print(f\"Shape of U: {u.shape}\")\n",
"print(f\"Shape of Sigma: {sigma.shape}\")\n",
"print(f\"Shape of V^T: {vh.shape}\")"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix U_r is\n",
"[[-0.36206122 0.6246702 0.47822465]\n",
" [-0.34561768 0.00735298 -0.79403674]\n",
" [-0.14568802 -0.7687844 0.37114088]\n",
" [-0.85336692 -0.13676122 0.05532879]]\n",
"\n",
"The matrix U_{r, perp} is\n",
"[[-0.5]\n",
" [-0.5]\n",
" [-0.5]\n",
" [ 0.5]]\n",
"\n",
"The matrix V_r^T is\n",
"[[-0.50334771 -0.42369658 -0.42055044 -0.46194908 -0.42055044]\n",
" [ 0.74605012 -0.17784578 -0.44072325 0.15266344 -0.44072325]\n",
" [-0.09073097 0.88817346 -0.28886987 -0.18980042 -0.28886987]]\n",
"\n",
"The matrix V_{r, perp}^T is\n",
"[[ 3.45876835e-02 6.24500451e-17 7.22070511e-01 -6.91753671e-02\n",
" -6.87482827e-01]\n",
" [ 4.24996322e-01 0.00000000e+00 1.55140977e-01 -8.49992645e-01\n",
" 2.69855345e-01]]\n"
]
}
],
"source": [
"# (d) Submatrices\n",
"r = 3 # rank of Q\n",
"ur = u[:, :r]\n",
"urperp = u[:, r:]\n",
"vhr = vh[:r, :]\n",
"vhrperp = vh[r:, :]\n",
"\n",
"print(\"The matrix U_r is\")\n",
"print(ur)\n",
"\n",
"print(\"\\nThe matrix U_{r, perp} is\")\n",
"print(urperp)\n",
"\n",
"print(\"\\nThe matrix V_r^T is\")\n",
"print(vhr)\n",
"\n",
"print(\"\\nThe matrix V_{r, perp}^T is\")\n",
"print(vhrperp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) From the discussion above, we know that $\\text{null}(Q) = \\text{range}(V_{r,\\perp})$. We computed $V_{r,\\perp}^T$ in part (d), which has 2 rows. So, $\\text{null}(Q)$ is the span of these 2 row vectors, which live in $\\mathbb{R}^5$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(f) From the discussion above, we know that $\\text{range}(Q) = \\text{range}(U_r)$. We computed $U_r$ in part (d), which has 3 columns. So, $\\text{range}(U_r)$ is the span of these 3 column vectors, which live in $\\mathbb{R}^4$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}