{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving Linear Systems" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.linalg as la" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solve $Ax = b$ where $A$ is full rank square matrix\n", "\n", "Solve $Ax = b$ for all the $b$ vectors" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(123)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[2, 2, 6],\n", " [1, 3, 9],\n", " [6, 1, 0]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = np.random.randint(0, 10, (3,3))\n", "A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check if A is invertible." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-12.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.det(A)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 9, 0, 0, 9, 3, 4, 0, 0, 4],\n", " [1, 7, 3, 2, 4, 7, 2, 4, 8, 0],\n", " [7, 9, 3, 4, 6, 1, 5, 6, 2, 1]])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = np.random.randint(0, 10, (3, 10))\n", "b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Direct matrix inversion \n", "\n", "Not recommended. Less numerically stable and slower when there are multiple $b$ to solve for." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.39 ms, sys: 0 ns, total: 1.39 ms\n", "Wall time: 1.12 ms\n" ] } ], "source": [ "%%time\n", "\n", "x = np.empty((3, 10))\n", "for i in range(10):\n", " x[:, i] = la.inv(A) @ b[:, i]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.25 , 3.25 , -1.5 , -1. ,\n", " 4.75 , -1.25 , 2. , -2. ,\n", " -4. , 3. ],\n", " [ 5.5 , -10.5 , 12. , 10. ,\n", " -22.5 , 8.5 , -7. , 18. ,\n", " 26. , -17. ],\n", " [ -1.75 , 3.91666667, -3.5 , -3. ,\n", " 7.41666667, -1.91666667, 2.33333333, -5.33333333,\n", " -7.33333333, 5.33333333]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using `solve`\n", "\n", "Recommended method. " ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 741 µs, sys: 1.13 ms, total: 1.88 ms\n", "Wall time: 779 µs\n" ] }, { "data": { "text/plain": [ "array([[ 0.25 , 3.25 , -1.5 , -1. ,\n", " 4.75 , -1.25 , 2. , -2. ,\n", " -4. , 3. ],\n", " [ 5.5 , -10.5 , 12. , 10. ,\n", " -22.5 , 8.5 , -7. , 18. ,\n", " 26. , -17. ],\n", " [ -1.75 , 3.91666667, -3.5 , -3. ,\n", " 7.41666667, -1.91666667, 2.33333333, -5.33333333,\n", " -7.33333333, 5.33333333]])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "\n", "la.solve(A, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using PLU decomposition\n", "\n", "The permuation matrix $P$ indicates the row switches. Solving requires\n", "\n", "- forward substitution $Ly = P^Tb$\n", "- backward substitution $Ux = y$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 842 µs, sys: 1.3 ms, total: 2.15 ms\n", "Wall time: 795 µs\n" ] } ], "source": [ "%%time\n", "\n", "P, L, U = la.lu(A)\n", "y = la.solve_triangular(L, P.T@b, lower=True)\n", "x = la.solve_triangular(U, y, lower=False)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.25 , 3.25 , -1.5 , -1. ,\n", " 4.75 , -1.25 , 2. , -2. ,\n", " -4. , 3. ],\n", " [ 5.5 , -10.5 , 12. , 10. ,\n", " -22.5 , 8.5 , -7. , 18. ,\n", " 26. , -17. ],\n", " [ -1.75 , 3.91666667, -3.5 , -3. ,\n", " 7.41666667, -1.91666667, 2.33333333, -5.33333333,\n", " -7.33333333, 5.33333333]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving least squares equations\n", "\n", "We solve $X\\beta = y$ where the dimensions of $X$ are such that $m > n$. Here we interpreet the solution as finding the best coefficinets for fittting the eequation \n", "\n", "$y = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$\n", "\n", "Note that we are looking for the projection of $y$ onto the column space (image) of the column vectors of $X$. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "X = np.c_[np.ones((10, 1)), np.random.randint(0, 10, (10, 2))]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 8., 3.],\n", " [1., 5., 0.],\n", " [1., 2., 6.],\n", " [1., 2., 4.],\n", " [1., 4., 6.],\n", " [1., 3., 0.],\n", " [1., 6., 4.],\n", " [1., 7., 6.],\n", " [1., 7., 1.],\n", " [1., 5., 7.]])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "y = np.random.randint(0, 10, (10, 1))" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[9],\n", " [2],\n", " [4],\n", " [8],\n", " [1],\n", " [2],\n", " [1],\n", " [1],\n", " [3],\n", " [5]])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using `lstsq` method\n", "\n", "Recommended." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 643 µs, sys: 0 ns, total: 643 µs\n", "Wall time: 623 µs\n" ] } ], "source": [ "%%time\n", "\n", "β, res, rank, s = la.lstsq(X, y)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 3.51230193],\n", " [-0.02659447],\n", " [ 0.05892189]])" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "β" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using the normal equations" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 3.51230193],\n", " [-0.02659447],\n", " [ 0.05892189]])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.solve(X.T@X, X.T@y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using optimization" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "def res(β, X, y):\n", " β = β.reshape((-1,1))\n", " return (y - X@β).T @ (y - X@β)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "β0 = np.zeros(3)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "import scipy.optimize as opt" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 76.138865\n", " Iterations: 169\n", " Function evaluations: 304\n" ] }, { "data": { "text/plain": [ "array([ 3.51234729, -0.02660293, 0.05892217])" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opt.fmin(res, β0, args=(X, y))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }