Tarea 2: Álgebra Lineal y Descomposición SVD

Teoría de Algebra Lineal y Optimización

1. ¿Por qué una matriz equivale a una transformación lineal entre espacios vectoriales?

Una multiplicación de una matriz, implica escalamientos, reducciones de dimensiones y rotaciones que se pueden observar en el espacio geométrico. Todas estas transformaciones son lineales en tanto se pueden representar como combinaciones lineales de la forma $\alpha*x + \beta$, donde $\alpha$ y $\beta$ son parámetros que determina la matriz.

2. ¿Cuál es el efecto de transformación lineal de una matriz diagonal y el de una matriz ortogonal?

Una matriz diagonal va a rescalar cada una de las columnas del vector o matriz al que se le aplique la transformación lineal, mientras que una matriz ortogonal va a generar una transformación isométrica que puede ser cualquiera de tres tipos: rotación, traslación o generar el reflejo.

3. ¿Qué es la descomposición en valores singulares de una matriz (Single Value Descomposition)?

Es la factorización de una matriz en la multiplicación de tres matrices: la de eigenvectores, los componentes singulares y la transpuesta de eigenvectores. La descomposicion se lleva a cabo de la siguiente manera: $ A = U \Sigma V^{T} $ En donde U = matriz de eigenvectores, $\Sigma$ = matriz de valores singulares o raices cuadradas de eigenvalores de $A^{T}A$ y $V^{T}$ = matriz transpuesta de eigenvectores por derecha. La descomposición SVD tiene muchas aplicaciones en las áreas de principal component analytsis, dimensionality reduction, compresión de imágenes, etcétera.

4. ¿Qué es diagonalizar una matriz y que representan los eigenvectores?

Es la factorización de la matriz en las tres matrices básicas mencionadas arriba (eigenvectores, eigenvalores e inversa de eigenvectores). Esto se observa de la siguiente manera: $A = PDP^{-1}$ donde P = matriz de eigenvectores, D = Matriz diagonal de eigenvalores y $P^{-1}$ = matriz inversa de eigenvectores. Los eigenvectores de una transformación lineal son vectores que cumplen que el resultado de multiplicarlos por la matriz de transformación equivale a multiplicarlos por un escalar que llamamos 'eigenvalor' ($ Ax = \lambda x $). Esto implica que a estos vectores al aplicarles la transformación lineal (multiplicarlos por la matriz de transformación) no cambian su dirección solamente cambian su longitud o su sentido (en caso de que haya eigenvalor con signo negativo).

5. ¿Intuitivamente que son los eigenvectores?

Se pueden interpretar como ejes cartesianos de la transformación lineal.

6. ¿Cómo interpretas la descomposición en valores singulares como una composición de tres tipos de transformaciones lineales simples?

Las tres matrices que conforman la SVD se pueden ver como transformaciones simples que son una rotación inicial, un escalamiento a lo largo de los ejes principales (valores singulares) y una rotación final.

7. ¿Qué relación hay entre la descomposición en valores singulares y la diagonalización?

La descomposición svd es una generalización de la diagonalización de una matriz. Cuando una matriz no es cuadrada, no es diagonalizable, sin embargo si se puede descomponer mediante svd. Así mismo si se utiliza la descomposición svd de una matriz, es posible resolver sistemas de ecuaciones lineales que no tengan una única solución sino que la solución devuelta es el ajuste a mínimos cuadrados, (obviamente si tienen solución, regresa la solución al sistema).

8. ¿Cómo se usa la descomposición en valores singulares para dar una aproximación de rango menor a una matriz?

En la descomposición svd de una matriz mxn, se obtiene una descomposición de tres matrices cuya multiplicación resulta en la matriz original completa, sin embargo dadas las propiedades de la descomposición, se pueden utilizar solamente las componentes principales de mayor relevancia desechando columnas de la matriz $U_{mxr}$, filas o columnas de la matriz de componentes principales $\Sigma _{rxr} $ y filas de la matriz $V_{rxn}^{T}$ La multiplicación de estas resulta en el tamaño mxn de la matriz original pero con un determinado error y a medida que se consideran más componentes principales, columnas y vectores de las matrices U,S,VT , mejor será la reconstrucción de la matriz original. Esta descomposición es muy útil para comprimir información y análisis de componentes principales (PCA).

9. Describe el método de minimización por descenso gradiente

Es un método iterativo de primer orden de minimización de una función dado su gradiente. Conociendo la función a optimizar, se calcula el gradiente (vector de derivadas parciales), posteriormente se considera un punto aleatorio de inicio, se substituye la coordenada en el vector gradiente, se analiza en cuál de las componentes del vector gradiente (x, y, o z) es el valor más negativo y se da un pequeño incremento en esa dirección de manera que se acerca al mínimo local, y así sucesivamente hasta llegar a un punto de convergencia (mínimo local). Entre sus aplicaciones están: Encontrar mínimos locales de funciones, Solución de sistemas de ecuaciones lineales y solución de sistemas de ecuaciones no lineales.

10. Menciona 4 ejemplos de problemas de optimización(dos con restricciones y dos sin restricciones) que te parezcan interesantes como científico de datos

Con restricciones: Modelos de crecimiento económico, Modelos de impuestos óptimos. Sin restricciones: Maximización intertemporal de extracción de recursos de una mina, minimización la varianza.

Aplicaciones en Python

Ejercicio 1

Recibir el path de un archivo de una imagen y convertirlo en una matriz numérica que represente la versión en blanco y negro de la imagen


In [24]:
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np

#url = sys.argv[1]
url = 'pikachu.png'

img = Image.open(url)

imggray = img.convert('LA')

Realizar y verificar la descomposición SVD


In [25]:
imggrayArray = np.array(list(imggray.getdata(band=0)), float)
imggrayArray.shape = (imggray.size[1], imggray.size[0])
imggrayArray = np.matrix(imggrayArray)

plt.imshow(imggray)

plt.show()

u, s, v = np.linalg.svd(imggrayArray)

print("U: ")
print(u)
print("S: ")
print(s)
print("V: ")
print(v)


U: 
[[-0.03210234 -0.01230276  0.01447713 ...,  0.00031968  0.00701922
  -0.19292895]
 [-0.03204528 -0.01400614  0.01570612 ...,  0.00941665  0.05242116
   0.08919738]
 [-0.03201338 -0.01494325  0.01640772 ..., -0.01038472  0.06065611
  -0.08362323]
 ..., 
 [-0.03173323 -0.01063853  0.00752548 ..., -0.07724167 -0.07910892
  -0.01163805]
 [-0.03188642 -0.01072858  0.00994695 ...,  0.10842794 -0.01935532
   0.07194286]
 [-0.03206187 -0.01065193  0.0120968  ...,  0.0013839   0.12562333
   0.06909291]]
S: 
[  2.79287511e+05   3.15299224e+04   1.92331515e+04 ...,   3.09418230e-11
   1.97960667e-11   1.69294955e-11]
V: 
[[ -3.22810737e-02  -3.22810737e-02  -3.22810737e-02 ...,  -3.22810737e-02
   -3.22810737e-02  -3.22810737e-02]
 [  5.05425041e-03   5.05425041e-03   5.05425041e-03 ...,   5.05425041e-03
    5.05425041e-03   5.05425041e-03]
 [  2.29661334e-04   2.29661334e-04   2.29661334e-04 ...,   2.29661334e-04
    2.29661334e-04   2.29661334e-04]
 ..., 
 [  0.00000000e+00   2.31701671e-18  -4.75710117e-19 ...,   3.82000361e-18
    1.04609550e-17  -3.11632844e-18]
 [  0.00000000e+00   1.75372480e-01  -2.07157595e-02 ...,  -6.07112677e-03
   -1.80198425e-01  -1.71406168e-01]
 [ -9.86038849e-01  -8.99936678e-02   1.54318692e-02 ...,   1.76419364e-02
    1.04520228e-02   1.04848300e-02]]

Para alguna imagen de su elección, elegir distintos valores de aproximación a la imagen original


In [26]:
for i in range(1,60,15):
    reconstimg = np.matrix(u[:, :i]) * np.diag(s[:i]) * np.matrix(v[:i,:])
    plt.imshow(reconstimg, cmap='gray')
    plt.show()


Ejercicio 2


In [5]:
def lin_solve_pseudo(A,b): 
    pseudoinv = pseudoinverse(A)
    return np.matmul(pseudoinv,b)

def pseudoinverse(A): 
    u,s,v = np.linalg.svd(A)            
    diagonal = np.diag(s)
    if v.shape[0] > diagonal.shape[1]:       
        print("Agregando columnas a la sigma")
        vector = np.array([[0 for x in range(v.shape[0] - diagonal.shape[1])] for y in range(diagonal.shape[0])])                      
        diagonal = np.concatenate((diagonal, vector), axis=1)

    elif u.shape[1] > diagonal.shape[0]:
        print("Agregando renglones a la sigma")
        vector = np.array([[0 for x in range(diagonal.shape[0])] for y in range(u.shape[1]-diagonal.shape[0])])
        diagonal = np.concatenate((diagonal, vector), axis=0)
        
    for a in range(diagonal.shape[0]):
        for b in range(diagonal.shape[1]):
            if diagonal[a][b] != 0:
                diagonal[a][b] = 1/diagonal[a][b]
                
    resultante = np.dot(np.transpose(v),np.transpose(diagonal))
    resultante = np.dot(resultante,np.transpose(u))
    return resultante

Programar una función que dada cualquier matriz devuelva la pseudoinversa usando la descomposición SVD. Hacer otra función que resuelva cualquier sistema de ecuaciones de la forma Ax=b usando esta pseudoinversa.


In [6]:
A = np.array([[1,1,1],[1,1,3],[2,4,4]])
b = np.array([[18,30,68]])

solve = lin_solve_pseudo(A,np.transpose(b))
print(solve)


[[ 2.]
 [10.]
 [ 6.]]

2. Jugar con el sistema Ax=b donde A=[[1,1],[0,0]] y b puede tomar distintos valores.

(a) Observar que pasa si b esta en la imagen de A (contestar cuál es la imagen) y si no está (ej. b = [1,1]).


In [7]:
print("La imagen de A es cualquier vector de dos coordenadas en donde la segunda componente siempre sea cero")
print("Vector b en imagen de A")
A = np.array([[1,1],[0,0]])
b = np.array([[12,0]])
solve = lin_solve_pseudo(A, np.transpose(b))
print(solve)
print("Cuando b esta en la imagen, la funcion lin_solve_pseudo devuelve la solucion unica a su sistema")
print("Vector b no en imagen de A")
b = np.array([[12,8]])
solve = lin_solve_pseudo(A, np.transpose(b))
print(solve)
print("Cuando b no esta en la imagen, devuelve la solucion mas aproximada a su sistema")


La imagen de A es cualquier vector de dos coordenadas en donde la segunda componente siempre sea cero
Vector b en imagen de A
[[ 6.]
 [ 6.]]
Cuando b esta en la imagen, la funcion lin_solve_pseudo devuelve la solucion unica a su sistema
Vector b no en imagen de A
[[ 6.]
 [ 6.]]
Cuando b no esta en la imagen, devuelve la solucion mas aproximada a su sistema

(b) Contestar, ¿la solución resultante es única? Si hay más de una solución, investigar que carateriza a la solución devuelta.


In [8]:
print("Vector b no en imagen de A")
b = np.array([[12,8]])
solve = lin_solve_pseudo(A, np.transpose(b))
print(solve)
print("Cuando b no esta en la imagen, devuelve la solucion mas aproximada a su sistema")


Vector b no en imagen de A
[[ 6.]
 [ 6.]]
Cuando b no esta en la imagen, devuelve la solucion mas aproximada a su sistema

(c) Repetir cambiando A=[[1,1],[0,1e-32]], ¿En este caso la solucíon es única? ¿Cambia el valor devuelto de x en cada posible valor de b del punto anterior?


In [9]:
A = np.array([[1,1],[0,1e-32]])
b = np.array([[12,9]])
solve = lin_solve_pseudo(A, np.transpose(b))
print(solve)
cadena = """En este caso, la solucion devuelta siempre es el valor de la segunda coordenada del vector b por e+32\
             y es el valor de ambas incognitas, solo que con signos contrarios ej(x1=-9.0e+32, x2=9.0e+32) \
             esto debido a que cualquier numero entre un numero muy pequenio tiende a infinito, de manera que la \
             coordenada dos del vector tiene mucho peso con referencia a la coordenada uno del vector
         """
print(cadena)


[[ -9.00000000e+32]
 [  9.00000000e+32]]
En este caso, la solucion devuelta siempre es el valor de la segunda coordenada del vector b por e+32             y es el valor de ambas incognitas, solo que con signos contrarios ej(x1=-9.0e+32, x2=9.0e+32)              esto debido a que cualquier numero entre un numero muy pequenio tiende a infinito, de manera que la              coordenada dos del vector tiene mucho peso con referencia a la coordenada uno del vector
         

Ejercicio3

Leer el archivo study_vs_sat.csv y almacenearlo como un data frame de pandas.


In [11]:
import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv("./study_vs_sat.csv", sep=',')
print(data)


    study_hours  sat_score
0             4        390
1             9        580
2            10        650
3            14        730
4             4        410
5             7        530
6            12        600
7            22        790
8             1        350
9             3        400
10            8        590
11           11        640
12            5        450
13            6        520
14           10        690
15           11        690
16           16        770
17           13        700
18           13        730
19           10        640

Plantear como un problema de optimización que intente hacer una aproximación de la forma sat_score ~ alpha + beta*study_hours minimizando la suma de los errores de predicción al cuadrado, ¿Cuál es el gradiente de la función que se quiere optimizar (hint: las variables que queremos optimizar son alpha y beta)?


In [12]:
hrs_studio = np.array(data["study_hours"])
sat_score = np.array(data["sat_score"])

A = np.vstack([hrs_studio, np.ones(len(hrs_studio))]).T

m,c = np.linalg.lstsq(A,sat_score)[0]

print("Beta y alfa: ")
print(m,c)


Beta y alfa: 
25.3264677779 353.164879499

Programar una función que reciba valores de alpha, beta y el vector study_hours y devuelva un vector array de numpy de predicciones alpha + beta*study_hours_i, con un valor por cada individuo


In [13]:
def predict(alfa, beta, study_hours):
    study_hours_i=[]
    for a in range(len(study_hours)):
        study_hours_i.append(alfa + beta*np.array(study_hours[a]))
    return study_hours_i

print("prediccion")
print(predict(353.165, 25.326, hrs_studio))


prediccion
[454.46900000000005, 581.09900000000005, 606.42499999999995, 707.72900000000004, 454.46900000000005, 530.447, 657.077, 910.33699999999999, 378.49100000000004, 429.14300000000003, 555.77300000000002, 631.75099999999998, 479.79500000000002, 505.12100000000004, 606.42499999999995, 631.75099999999998, 758.38100000000009, 682.40300000000002, 682.40300000000002, 606.42499999999995]

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 + beta study_hours_i en cada entrada y que entonces el problema se vuelve sat_score ~ X*[alpha,beta]


In [14]:
unos = np.ones((len(hrs_studio),1))
hrs_studio = [hrs_studio]
hrs_studio = np.transpose(hrs_studio)
x = np.hstack((unos, hrs_studio))
print("La prediccion es: ")
print(np.matmul(x,np.array([[353.165],[25.326]])))


La prediccion es: 
[[ 454.469]
 [ 581.099]
 [ 606.425]
 [ 707.729]
 [ 454.469]
 [ 530.447]
 [ 657.077]
 [ 910.337]
 [ 378.491]
 [ 429.143]
 [ 555.773]
 [ 631.751]
 [ 479.795]
 [ 505.121]
 [ 606.425]
 [ 631.751]
 [ 758.381]
 [ 682.403]
 [ 682.403]
 [ 606.425]]

Calculen la pseudoinversa X^+ de X y computen (X^+)*sat_score para obtener alpha y beta soluciones.


In [15]:
X_pseudo = pseudoinverse(x)
print("Las alfas y betas son: ")
print(np.matmul(X_pseudo,sat_score))


Agregando renglones a la sigma
Las alfas y betas son: 
[ 353.1648795    25.32646778]

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


In [16]:
def comparacion(X, sat_score):
    x_transpose = np.transpose(X)    
    return np.matmul(np.linalg.inv(np.matmul(x_transpose,X)), np.matmul(x_transpose,sat_score))

Usen la libreria matplotlib para visualizar las predicciones con alpha y beta solución contra los valores reales de sat_score.


In [27]:
plt.plot(hrs_studio, sat_score, 'x', label='Datos', markersize=20)
plt.plot(hrs_studio, m*hrs_studio + c, 'r', label='Línea de regresión')
plt.legend()
plt.show()


Programen el método de descenso gradiente para obtener alpha y beta por vía de un método numérico de optimización. Experimenten con distintos learning rates (tamaños de paso).

No hice el ejercicio muy avanzado.