{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"AeroPython\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Álgebra Lineal con NumPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Una vez hemos visto el manejo básico de arrays en Python con NumPy, es hora de pasar a operaciones más interesantes como son las propias del Álgebra Lineal._\n", "\n", "_Los productos escalares y las inversiones de matrices están por todas partes en los programas científicos e ingenieriles, así que vamos a estudiar cómo se realizan en Python._" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Álgebra lineal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Como sabemos, las operaciones del álgebra lineal aparecen con mucha frecuencia a la hora de resolver sistemas de ecuaciones en derivadas parciales y en general al linealizar problemas de todo tipo, y suele ser necesario resolver sistemas con un número enorme de ecuaciones e incógnitas. Gracias a los arrays de NumPy podemos abordar este tipo de cálculos en Python, ya que todas las funciones están escritas en C o Fortran y tenemos la opción de usar bibliotecas optimizadas al límite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "El paquete de álgebra lineal en NumPy se llama `linalg`, así que importando NumPy con la convención habitual podemos acceder a él escribiendo `np.linalg`. Si imprimimos la ayuda del paquete vemos que tenemos funciones para:\n", "\n", "* Funciones básicas (norma de un vector, inversa de una matriz, determinante, traza)\n", "* Resolución de sistemas\n", "* Autovalores y autovectores\n", "* Descomposiciones matriciales (QR, SVD)\n", "* Pseudoinversas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Puede que ya sepas que en la biblioteca `SciPy` se pueden encontrar también funciones de Álgebra Lineal. ¿Cuáles usar? Puedes encontrar la respuesta en este enlace: https://docs.scipy.org/doc/scipy-0.18.1/reference/tutorial/linalg.html#scipy-linalg-vs-numpy-linalg\n", "\n", "Como de momento sólo hemos usado `NumPy` no importaremos las funciones de `SciPy`, aunque como ves, es recomendable hacerlo." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on package numpy.linalg in numpy:\n", "\n", "NAME\n", " numpy.linalg\n", "\n", "DESCRIPTION\n", " Core Linear Algebra Tools\n", " -------------------------\n", " Linear algebra basics:\n", " \n", " - norm Vector or matrix norm\n", " - inv Inverse of a square matrix\n", " - solve Solve a linear system of equations\n", " - det Determinant of a square matrix\n", " - lstsq Solve linear least-squares problem\n", " - pinv Pseudo-inverse (Moore-Penrose) calculated using a singular\n", " value decomposition\n", " - matrix_power Integer power of a square matrix\n", " \n", " Eigenvalues and decompositions:\n", " \n", " - eig Eigenvalues and vectors of a square matrix\n", " - eigh Eigenvalues and eigenvectors of a Hermitian matrix\n", " - eigvals Eigenvalues of a square matrix\n", " - eigvalsh Eigenvalues of a Hermitian matrix\n", " - qr QR decomposition of a matrix\n", " - svd Singular value decomposition of a matrix\n", " - cholesky Cholesky decomposition of a matrix\n", " \n", " Tensor operations:\n", " \n", " - tensorsolve Solve a linear tensor equation\n", " - tensorinv Calculate an inverse of a tensor\n", " \n", " Exceptions:\n", " \n", " - LinAlgError Indicates a failed linear algebra operation\n", "\n", "PACKAGE CONTENTS\n", " _umath_linalg\n", " info\n", " lapack_lite\n", " linalg\n", " setup\n", "\n", "DATA\n", " absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0...\n", " division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192...\n", " print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)...\n", "\n", "FILE\n", " /miniconda/lib/python3.5/site-packages/numpy/linalg/__init__.py\n", "\n", "\n" ] } ], "source": [ "help(np.linalg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recordemos que si queremos usar una función de un paquete pero no queremos escribir la \"ruta\" completa cada vez, podemos usar la sintaxis `from package import func`:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from numpy.linalg import norm, det\n", "norm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "El producto matricial usual (no el que se hace elemento a elemento, sino el del álgebra lineal) se calcula con la misma función que el producto matriz-vector y el producto escalar vector-vector: con la función `dot`, que **no** está en el paquete `linalg` sino directamente en `numpy` y no hace falta importarlo." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Una consideración importante a tener en cuenta es que en NumPy no hace falta ser estricto a la hora de manejar vectores como si fueran matrices columna, siempre que la operación sea consistente. Un vector es una matriz con una sola dimensión: por eso si calculamos su traspuesta no funciona." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "M = np.array([\n", " [1, 2],\n", " [3, 4]\n", "])\n", "v = np.array([1, -1])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1, -1])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v.T" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1, -1])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u = np.dot(M, v)\n", "u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para hacer comparaciones entre arrays de punto flotante se pueden usar las funciones `np.allclose` y `np.isclose`. La primera comprueba si todos los elementos de los arrays son iguales dentro de una tolerancia, y la segunda compara elemento a elemento y devuelve un array de valores `True` y `False`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([-1, -1]), array([ 1, -1]))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u, v" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(u, v)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isclose(0.0, 1e-8, atol=1e-10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__En la versión 3.5 de Python se incorporó un nuevo operador `@` para poder calcular hacer multiplicaciones de matrices de una forma más legible__" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1, -1])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u = M @ v\n", "u" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([23, 67, 66, 49])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat = np.array([[1, 5, 8, 5],\n", " [0, 6, 4, 2],\n", " [9, 3, 1, 6]])\n", "\n", "vec1 = np.array([5, 6, 2])\n", "\n", "vec1 @ mat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si quieres saber más, puedes leer [este artículo](http://pybonacci.org/2016/02/22/el-producto-de-matrices-y-el-nuevo-operador/) en Pybonacci escrito por _Álex Sáez_." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Ejercicios" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1- Hallar el producto de estas dos matrices y su determinante:\n", "\n", "$$\\begin{pmatrix} 1 & 0 & 0 \\\\ 2 & 1 & 1 \\\\ -1 & 0 & 1 \\end{pmatrix} \\begin{pmatrix} 2 & 3 & -1 \\\\ 0 & -2 & 1 \\\\ 0 & 0 & 3 \\end{pmatrix}$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from numpy.linalg import det" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1 0 0]\n", " [ 2 1 1]\n", " [-1 0 1]]\n", "[[ 2 3 -1]\n", " [ 0 -2 1]\n", " [ 0 0 3]]\n" ] } ], "source": [ "A = np.array([\n", " [1, 0, 0],\n", " [2, 1, 1],\n", " [-1, 0, 1]\n", "])\n", "B = np.array([\n", " [2, 3, -1],\n", " [0, -2, 1],\n", " [0, 0, 3]\n", "])\n", "print(A)\n", "print(B)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 2, 3, -1],\n", " [ 4, 4, 2],\n", " [-2, -3, 4]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "C = A @ B\n", "C" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-12.0" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "det(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2- Resolver el siguiente sistema:\n", "\n", "$$ \\begin{pmatrix} 2 & 0 & 0 \\\\ -1 & 1 & 0 \\\\ 3 & 2 & -1 \\end{pmatrix} \\begin{pmatrix} 1 & 1 & 1 \\\\ 0 & 1 & 2 \\\\ 0 & 0 & 1 \\end{pmatrix} \\begin{pmatrix} x \\\\ y \\\\ z \\end{pmatrix} = \\begin{pmatrix} -1 \\\\ 3 \\\\ 0 \\end{pmatrix} $$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 2, 2, 2],\n", " [-1, 0, 1],\n", " [ 3, 5, 6]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M = (np.array([[2, 0, 0],\n", " [-1, 1, 0],\n", " [3, 2, -1]])\n", " @\n", " np.array([[1, 1, 1],\n", " [0, 1, 2],\n", " [0, 0, 1]]))\n", "M" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.5, -4.5, 3.5])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.linalg.solve(M, np.array([-1, 3, 0]))\n", "x" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(M @ x, np.array([-1, 3, 0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3- Hallar la inversa de la matriz $H$ y comprobar que $H H^{-1} = I$ (recuerda la función `np.eye`)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2 4 7 12 21 38]\n", " [ 71 128 265 512 1035 2048]\n", " [ 4109 8206 16399 32784 65553 131090]\n", " [ 262145 524308 1048577 2097174 4194305 8388632]\n", " [ 16777271 33554488 67108921 134217786 268435515 536870972]\n", " [ 1073741855 2147483680 4294967329 8589934626 17179869219 34359738404]]\n" ] } ], "source": [ "# aeropython: preserve\n", "A = np.arange(1, 37).reshape(6,6)\n", "A[1, 1::2] = 0\n", "A[3, ::2] = 1\n", "A[4, :] += 30\n", "B = (2 ** np.arange(36)).reshape((6,6))\n", "H = A + B\n", "print(H)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "456269.58138193813" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.det(H)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "Hinv = np.linalg.inv(H)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[False, False, False, False, False, False],\n", " [False, False, False, False, False, False],\n", " [False, False, True, True, False, True],\n", " [False, False, False, False, False, False],\n", " [ True, False, True, True, False, True],\n", " [False, False, True, True, False, True]], dtype=bool)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isclose(np.dot(Hinv, H), np.eye(6))" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.125 0.25 0.5 1. 1. 4. ]\n", " [-0.094 0.875 -0.25 -0.5 -0.5 -3. ]\n", " [-0.125 0.25 1. 0. 4. 0. ]\n", " [ 0.188 0.25 0.5 2. 1. 6. ]\n", " [ 0. -0.125 0. 0. 0. 0. ]\n", " [-0.031 0.062 0. 0. 1. 1. ]]\n" ] } ], "source": [ "np.set_printoptions(precision=3)\n", "print(np.dot(Hinv, H))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
¡No funciona! Y no solo eso sino que los resultados varían de un ordenador a otro.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4- Comprobar el número de condición de la matriz $H$." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0682988788334074e+17" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.cond(H)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
La matriz está mal condicionada.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Autovalores y autovectores " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La cosa no queda ahí y también se pueden resolver problemas de autovalores y autovectores:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 1., 1., 1.]), array([[ 0.000e+00, 0.000e+00, 4.930e-32],\n", " [ 1.000e+00, -1.000e+00, -1.000e+00],\n", " [ 0.000e+00, 2.220e-16, 2.220e-16]]))" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = np.array([\n", " [1, 0, 0],\n", " [2, 1, 1],\n", " [-1, 0, 1]\n", "])\n", "\n", "np.linalg.eig(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Ya hemos aprendido a efectuar algunas operaciones útiles con NumPy. Estamos en condiciones de empezar a escribir programas más interesantes, pero aún queda lo mejor._\n", "\n", "* [Álgebra lineal en Python con NumPy en Pybonacci](http://pybonacci.org/2012/06/07/algebra-lineal-en-python-con-numpy-i-operaciones-basicas/)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "
\n", "####

¡Síguenos en Twitter!\n", "
\n", "###### Follow @AeroPython \n", "
\n", "###### Este notebook ha sido realizado por: Juan Luis Cano, Mabel Delgado y Álex Sáez \n", "
\n", "##### \"Licencia
Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "_Las siguientes celdas contienen configuración del Notebook_\n", "\n", "_Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como [seguro](http://ipython.org/ipython-doc/dev/notebook/security.html)_\n", "\n", " File > Trusted Notebook" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "/* This template is inspired in the one used by Lorena Barba\n", "in the numerical-mooc repository: https://github.com/numerical-mooc/numerical-mooc\n", "We thank her work and hope you also enjoy the look of the notobooks with this style */\n", "\n", "\n", "\n", "El estilo se ha aplicado =)\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Esta celda da el estilo al notebook\n", "from IPython.core.display import HTML\n", "css_file = '../styles/aeropython.css'\n", "HTML(open(css_file, \"r\").read())" ] } ], "metadata": { "anaconda-cloud": {}, "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.5.4" } }, "nbformat": 4, "nbformat_minor": 1 }