In [2]:
import matplotlib.pyplot as plt
%matplotlib inline
from skimage import io
from skimage.viewer import ImageViewer
from numpy import linalg
import numpy as np
In [25]:
matrix = io.imread('face.png', as_grey=True)
print(matrix.shape)
In [4]:
#### La clase image viewer te permite ver imagenes
#viewer = ImageViewer(img)
#viewer.show()
plt.imshow(matrix, cmap='gray')
plt.show()
Ahora la descomponemos en SVD.
In [26]:
matrix.shape
Out[26]:
In [27]:
u, s, vt = linalg.svd(matrix, full_matrices=False)
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)
In [28]:
S=np.diag(s)
m_check = np.matmul(np.matmul(u,S),vt)
print(m_check.shape)
Sea m la matrix en blanco y negro correspondiente a una imagen
A continuación se grafica la comprobación
In [29]:
plt.imshow(m_check, cmap='gray')
plt.show()
Ahora usaremos la SVD para aproximar un rango menor de m. K será 100.
In [30]:
k=10
In [31]:
u_r = u[:,:k]
S_r = S[:k,:k]
vt_r = vt[:k,:]
matrix_k = np.matmul(u_r, np.matmul(S_r,vt_r))
plt.imshow(matrix_k, cmap='gray')
plt.show()
Esta descomposición es fundamental cuando se trata de comprimir imágenes. La matriz aproximación K mostrada previamente es una compresión de la matriz original.
In [32]:
# A es una matriz cualquiera.
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Au, As, Avt = linalg.svd(A, full_matrices=True)
AS = np.diag(As)
In [33]:
Ainv = np.linalg.inv(A)
Au_inv = np.linalg.inv(Au)
As_inv = np.linalg.inv(AS)
Avt_inv = np.linalg.inv(Avt)
print(Au_inv,'\n', As_inv,'\n', Avt_inv)
Para Au y Avt, las transpestas son las inversas, en cuanto a S, por ser diagonal, los reciprocos de su diagonal, son la diagonal principal de la inversa de S, siempre y cuando distitnos de cero.
In [34]:
AS = np.diag(np.array([1/x if x != 0 else 0 for x in As]))
AS
Out[34]:
In [35]:
# Verificación
print(np.dot(Au, Au_inv))
print(np.dot(As, As_inv))
print(np.dot(Avt, Avt_inv))
In [ ]:
Por que la definición es muy clara. Si la SVD de $A = U\sum V^{T}$
La pseudoinversa de A, llamemosla $A^{+} = V\sum^{-1}U^{T} $
Multiplicar de esta forma sería un error $U^{T}\sum^{-1}V$
In [36]:
print(np.dot(Avt_inv,np.dot( As_inv,Au_inv)) )
# Ahora todo va a coincidir
print(np.linalg.pinv(A))
print(Ainv)
Programar una función que dada cualquier matriz devuelva la pseudainversa usando la descomposición SVD.
In [37]:
def get_pseudoinverse(matrix):
"""
Get pseudo inverse of a matrix
"""
Au, As, Avt = linalg.svd(matrix, full_matrices=False)
Au_inv = Au.transpose()
As_inv = np.diag(np.array([1/x if x != 0 else 0 for x in As]))
Avt_inv = Avt.transpose()
pseudo_matrix = np.dot(Avt_inv,np.dot( As_inv, Au_inv))
return pseudo_matrix
In [38]:
#Verificamos
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Ainv = get_pseudoinverse(A)
iden = np.dot(A,Ainv)
iden
Out[38]:
In [39]:
def solve_system(coefficient_matrix, image_vector):
pseudoinverse = get_pseudoinverse(coefficient_matrix)
domain_vector = np.dot(pseudoinverse,image_vector)
return domain_vector
In [40]:
### Prueba de calidad
B = np.array([[-10,9],[10,5]])
b = np.array([-9,-5])
print(B,'\n',b)
x = solve_system(B,b)
print(x)
np.dot(B,x)
###### Funciona!!!!!!!!
Out[40]:
In [41]:
A = np.array([[1,1],[0,0]])
A
Out[41]:
La imagen de A es el siguiente conjunto $ b \hspace{3mm} x \hspace{3mm} \left(\begin{array}{ccc} 1 \\ 0 \end{array} \right)$ con $b \in \Bbb R$
In [42]:
b = np.array([1,1])
Evidentemente b no esta en la imagen La solución no trivial al sistema homogeneo es $ x \hspace{3mm} = \hspace{3mm} \left(\begin{array}{ccc} a \\ -a \end{array} \right)$ con $a \in \Bbb R$
In [43]:
solve_system(A,b)
Out[43]:
Por la pseudoinversa se encuentra una solución que no es otra cosa más que el vector $b$ proyectado en el espacio generado por la matriz A.
In [44]:
np.linalg.solve(A,b)
La matriz A es singular por definición no tiene inversa. Sin embargo, por el método de la pseudo inversa podemos aproximarla solución aún cuando no esté en el espacio de la imagen de x o, equivalentemente, el espacio generado por $A$ el método que diseñamos, obtiene la pseudoinversa y aproxima una solución proyectando el vector $b$ sobre el espacio generado por la matriz A.
Si usaramos el tradicional de python o R, de obten la inversa, alzaría una excepción del tipo Error, la matriz es singular.
Ahora $ A = \left( \begin{array}{ccc} 1 & 1 \\ 0 & 1e-32 \end{array} \right)$
In [45]:
A = np.array([[1,1], [0,1e-32]])
A
Out[45]:
In [48]:
np.linalg.det(A)
Out[48]:
El determinante es distinto de cero, la inversa existe.
In [49]:
b = np.array([2,0])
c = np.array([1,1])
In [ ]:
In [50]:
solution_pseudo = solve_system(A,b)
solution_inv = np.linalg.solve(A,b)
print(np.dot(A,solution_pseudo))
print(np.dot(A,solution_inv))
In [51]:
solution_pseudo_c = solve_system(A,c)
solution_inv_c = np.linalg.solve(A,c)
print(solution_pseudo_c)
print(solution_inv_c)
#print(np.dot(A,solution_pseudo_c))
print(np.dot(A,solution_inv_c))
El método de la inversa sí encontró una solución, sin embargo, no es el vector (1,1)
In [52]:
print(solution_pseudo_c)
print(np.dot(A,solution_pseudo_c))
In [ ]:
In [55]:
import pandas as pd
In [56]:
base =pd.read_csv('study_vs_sat.csv')
base.head(5)
Out[56]:
Sea $ y $ sat_score y la queremos explicar mediante una relación lineal con $ x $ study_hours, entonces queremos encontrar $ \alpha , \beta $ tal que para cada observación aproximamos $ y_i$ con $ \hat{y}_i = \alpha_i + \beta x_i $
La función de error que deseamos minimizar: $ E = \sum\limits_{i=1}^n (y_i - \hat{y}_i )^2 = \sum\limits_{i=1}^n (y_i - \alpha_i + \beta x_i )^2 $
Encontrar $ \alpha $ y $ \beta $ tal que minimizen E.
Las condiciones de primer orden implican obtener las derivadas parciales con el vector gradiente e igualarlas 0.
Vector Gradiente $$ \nabla F(\alpha,\beta) = (\frac{\partial E}{\partial \alpha},\frac{\partial E}{\partial \beta}) $$
Condiciones de Primer Orden
$$ \frac{\partial E}{\partial \alpha} = -2 * \sum\limits_{i=1}^n (y_i - \alpha_i + \beta x_i ) = 0 $$
$$ \frac{\partial E}{\partial \beta} = -2 * \sum\limits_{i=1}^n (y_i - \alpha_i + \beta x_i )*x_i = 0 $$
La siguiente función recibe los parámetros alfa, beta y un numpy array de una dimensión con las horas de estudio y regresa un numpy array con predicciones de la forma $ y_i = \alpha + \beta x_i $
In [53]:
def predict_score( alpha, beta, study_hours):
predictions = np.array([alpha + beta * study_time for study_time in study_hours])
return predictions
Definan un numpy array X de dos columnas, la primera con unos en todas sus entradas y la segunda con la variable study_hours. Observen que X[alpha,beta] nos devuelve alpha + betastudy_hours_i en cada entrada y que entonces el problema se vuelve sat_score ~ X*[alpha,beta]
In [57]:
base =pd.read_csv('study_vs_sat.csv')
base.head(5)
study_hours = base.iloc[:,0]
study_hours = study_hours.values
alpha, beta = 2,3
X_array = np.empty( shape=[20,2], dtype=int)
for i in range(20):
X_array[i] = [1, study_hours[i]]
X_array * [alpha, beta]
Out[57]:
Calculen la pseudoinversa X^+ de X y computen (X^+)*sat_score para obtener alpha y beta soluciones
In [58]:
#La pseudoinversa de X es la siguiente
X_pseudoinv = get_pseudoinverse(X_array)
X_pseudoinv
Out[58]:
In [59]:
sat_score = base['sat_score'].values
solution_vector = np.dot(X_pseudoinv, sat_score)
print('Alfa óptima con pseudoinversa es :',solution_vector[0])
print('\n'*1)
print('Beta óptima con pseudoinversa es :', solution_vector[1])
Comparen la solución anterior con la de la fórmula directa de solución exacta (alpha,beta)=(X^tX)^(-1)X^t*sat_score
Resolvamos por partes, queremos $((X^{t}X))^{-1}X^{t}$ * study_hours
Sea $ J = ((X^{t}X))^{-1}$
In [60]:
J = np.linalg.inv(np.dot(X_array.transpose(), X_array))
K = np.dot(J, X_array.transpose())
sol_vector = np.dot(K, sat_score)
print('Alfa óptima con inversa tradicional es :',sol_vector[0])
print('\n'*1)
print('Beta óptima con inversa tradicional es :', sol_vector[1])
print('La solución es la misma que con la pseudoinversa')
In [ ]: