Álgebra Lineal¶
En esta clase abordaremos operaciones de álgebra lineal en Python, en particular el módulo de Álgebra lineal de la biblioteca SciPy (scipy.linalg).
Previos¶
En primer lugar, scipy.linalg contiene todas las funciones de numpy.linalg, en ocasiones agregando funcionalidades, pero posee algunas funciones adicionales. Además podrían existir diferencias de velocidad de cómputo dependiendo de como NumPy fue instalado, por lo que se recomiendo utilizar SciPy para tareas de álgebra lineal.
Como nuevamente no buscamos inventar la rueda de nuevo, esta clase está basada en el tutorial del módulo de álgebra lineal de Scipy (link).
Antes de comenzar, un tema que causa confusiones para matemáticos es el hecho que los objetos 2-D numpy.array (de dos dimensiones) no son matrices como se suelen estudiar en cursos de matemática. Un ejemplo de aquello es que la multiplicación de dos 2-D numpy.array no tienen que coincidir para una multiplicación matricial y además la multiplicación por defecto es elemento a elemento.
En algún momento de la historia de NumPy se optó por una clase especial numpy.matrix pero que al parecer trajo más confusiones consigo e incluso será removida en futuras versiones.
Rutinas Básicas¶
Inversa¶
Sabemos que la inversa de la matriz \(\mathbf{A}\) es la matriz \(\mathbf{A}^{-1}\), tal que
donde \(\mathbf{I}\) es la matriz identidad.
En SciPy, la matriz inversa de un NumPy array A se obtiene utilizando linalg.inv(A) o A.I.
import numpy as np
from scipy import linalg
A = np.array([[1,3,5],[2,5,1],[2,3,8]])
A
array([[1, 3, 5],
[2, 5, 1],
[2, 3, 8]])
linalg.inv(A)
array([[-1.48, 0.36, 0.88],
[ 0.56, 0.08, -0.36],
[ 0.16, -0.12, 0.04]])
A.dot(linalg.inv(A)) # double check
array([[ 1.00000000e+00, -1.11022302e-16, 4.85722573e-17],
[ 3.05311332e-16, 1.00000000e+00, 7.63278329e-17],
[ 2.22044605e-16, -1.11022302e-16, 1.00000000e+00]])
linalg.inv(A).dot(A) # double double check
array([[1.00000000e+00, 3.33066907e-16, 0.00000000e+00],
[0.00000000e+00, 1.00000000e+00, 4.44089210e-16],
[1.38777878e-17, 6.93889390e-18, 1.00000000e+00]])
Attention
Notar que hay un error de precisión pues está relacionado con la precisión de número flotante de NumPy.
Sistema de ecuaciones¶
Es tan sencillo como utilizar el comando linalg.solve que necesita como inputs un matriz y un vector, para luego retornar un vector solución. Por ejemplo para el sistema
se puede escribir en forma matricial
por lo que es posible resolverlo utilizando la inversa
A = np.array([[1, 2], [3, 4]])
A
array([[1, 2],
[3, 4]])
b = np.array([[5], [6]])
b
array([[5],
[6]])
linalg.inv(A).dot(b)
array([[-4. ],
[ 4.5]])
A.dot(linalg.inv(A).dot(b)) - b # check
array([[0.00000000e+00],
[1.77635684e-15]])
Sin embargo linalg.solve ofrece una forma de resolver es sistema mucho más amigable e inclusive más eficiente.
linalg.solve(A, b)
array([[-4. ],
[ 4.5]])
A.dot(linalg.solve(A, b)) - b # check
array([[0.],
[0.]])
Determinante¶
No hay mucho que decir aquí, más que tener en consideración problemas de presición.
A = np.array([[1,2],[3,4]])
A
array([[1, 2],
[3, 4]])
linalg.det(A)
-2.0
Normas¶
SciPy cuenta con una gran variedad de normas que puedes ser calculadas a través de un argumento en linalg.norm. Lo interesante es que esta función puede ser utilizada para elementos 1-D (vectores) o 2-D (matrices), además de tener un argumento opcional order para seleccionar una norma en particular (por defecto es 2).
Sea \(x\) un vector
Sea \(A\) una matrix
donde \(\sigma_{i}\) son los valores singulares de \(\mathbf{A}\).
A=np.array([[1,2],[3,4]])
A
array([[1, 2],
[3, 4]])
linalg.norm(A)
5.477225575051661
linalg.norm(A,'fro') # frobenius norm is the default
5.477225575051661
linalg.norm(A,1) # L1 norm (max column sum)
6.0
linalg.norm(A,-1)
4.0
linalg.norm(A,np.inf) # L inf norm (max row sum)
7.0
Descomposiciones¶
Espectral¶
La descomposición espectral es usual de encontrar en múltiples campos de la ingeniería y la matemática dadas las interpretaciones que se le pueden dar a los valores y vectores propios.
Se define la descomposición espectral para la matriz cuadrada \(\mathbf{A}\) como
donde \(\lambda\) es un escalar y \(\mathbf{v}\) un vector, denominados valor y vector propio respectivamente.
A = np.array([[1, 2], [3, 4]])
la, v = linalg.eig(A)
l1, l2 = la
print(l1, l2) # eigenvalues
(-0.3722813232690143+0j) (5.372281323269014+0j)
print(v[:, 0]) # first eigenvector
[-0.82456484 0.56576746]
print(v[:, 1]) # second eigenvector
[-0.41597356 -0.90937671]
print(np.sum(abs(v**2), axis=0)) # eigenvectors are unitary
[1. 1.]
v1 = np.array(v[:, 0]).T
print(linalg.norm(A.dot(v1) - l1*v1)) # check the computation
5.551115123125783e-17
Valor Singular¶
La descomposición valor singular (SVD) puede ser pensada como una extensión a la espectral cuando las matrices no son cuadradas. Sea \(\mathbf{A}\) de tamaño \(M\times N\), su descomposición valor singular es
aquí \(\mathbf{U}\) y \(\mathbf{V}\) son matrices unitarias de tamaño \(M \times M\) y \(N \times N\) respectivamente, mientras que \(\mathbf{\boldsymbol{\Sigma}}\) es una matriz \(M \times N\) donde los elementos de la diagonal principal son no nulos y se denotan usualmente como valores singulares, mientras que fuera de la diagonal sus valores son cero.
Matriz unitaria: Sea \(\mathbf{D}\) tal que \(\mathbf{D}^H \mathbf{D} = \mathbf{D}\mathbf{D}^H = \mathbf{I}\), es decir \(\mathbf{D}^{-1} = \mathbf{D}^H\).
A = np.array([[1,2,3],[4,5,6]])
A
array([[1, 2, 3],
[4, 5, 6]])
M,N = A.shape
U,s,Vh = linalg.svd(A)
Sig = linalg.diagsvd(s,M,N)
U, Vh = U, Vh
U
array([[-0.3863177 , -0.92236578],
[-0.92236578, 0.3863177 ]])
Sig
array([[9.508032 , 0. , 0. ],
[0. , 0.77286964, 0. ]])
Vh
array([[-0.42866713, -0.56630692, -0.7039467 ],
[ 0.80596391, 0.11238241, -0.58119908],
[ 0.40824829, -0.81649658, 0.40824829]])
U.dot(Sig.dot(Vh)) #check computation
array([[1., 2., 3.],
[4., 5., 6.]])
LU¶
La descomposición LU de una matriz \(\mathbf{A}\) de tamaño \(M\times N\) es
donde \(\mathbf{P}\) es una matriz de permutación de filas para la matriz identidad de \(M\times M\), \(\mathbf{L}\) es triangular inferior \(M\times K\) con \(K=\min\left(M,N\right)\)) con valores unitarios en la diagonal y \(\mathbf{U}\) es triangular superior.
A = np.array(
[
[2, 5, 8, 7],
[5, 2, 2, 8],
[7, 5, 6, 6],
[5, 4, 4, 8]]
)
p, l, u = linalg.lu(A)
l
array([[ 1. , 0. , 0. , 0. ],
[ 0.28571429, 1. , 0. , 0. ],
[ 0.71428571, 0.12 , 1. , 0. ],
[ 0.71428571, -0.44 , -0.46153846, 1. ]])
u
array([[ 7. , 5. , 6. , 6. ],
[ 0. , 3.57142857, 6.28571429, 5.28571429],
[ 0. , 0. , -1.04 , 3.08 ],
[ 0. , 0. , 0. , 7.46153846]])
np.allclose(A - p @ l @ u, np.zeros((4, 4)))
True
Cholesky¶
La descomposición Cholesky es un caso especial de la descomposición LU aplicable a matrices hermitianas definidas positivas. Cuando \(\mathbf{A}=\mathbf{A}^{H}\) y \(\mathbf{x}^{H}\mathbf{Ax}\geq 0\) para todo \(\mathbf{x}\),
donde \(\mathbf{L}\) es triangular inferior y \(\mathbf{U}\) es triangular superior.
Matriz hermitiana: Sea \(\mathbf{D}\) tal que \(\mathbf{D}^H = \mathbf{D}\).
A = np.array([[1,-2j],[2j,5]])
L = linalg.cholesky(A, lower=True)
L
array([[1.+0.j, 0.+0.j],
[0.+2.j, 1.+0.j]])
L @ L.T.conj()
array([[1.+0.j, 0.-2.j],
[0.+2.j, 5.+0.j]])
np.allclose(A, L @ L.T.conj())
True
QR¶
La descomposición QR es aplicable a cualquier matriz \(M\times N\), obteniendo una matriz unitaria \(\mathbf{Q}\) de tamaño \(M\times M\) y una matriz trapezoidal superior \(\mathbf{R}\) de \(M\times N\) tal que
A = np.random.randn(9, 6)
q, r = linalg.qr(A)
np.allclose(A, np.dot(q, r))
True
q.shape, r.shape
((9, 9), (9, 6))