Mostrando entradas con la etiqueta Métodos númericos. Mostrar todas las entradas
Mostrando entradas con la etiqueta Métodos númericos. Mostrar todas las entradas

7 ene 2020

El algoritmo del gradiente descendente

El gradiente descendente (GD) es un algoritmo de optimización genérico, capaz de encontrar soluciones óptimas para una amplia gama de problemas. La idea del gradiente descendente es ajustar los parámetros de forma iterativa para minimizar una función.

Concretamente, se tiene de una función diferenciable convexa1, $f:\Omega\subset\mathbb{R}^n \to \mathbb{R}$, el algoritmo GD permite encontrar un $w$ en $\Omega$ tal que $f(w)$ es un mínimo, en otras palabras, GD se utiliza para determinar los elementos del siguiente conjunto2:

\begin{equation}\tag{1} w \in \operatorname*{argmin\,\,}_{ w\in \Omega} f(w). \end{equation}

Para determinar los valores de $w$ que optimiza la función $f(w)$, GD hace uso de una serie de iteracciones que se hacen de acuerdo con la siguiente regla de actualización:

\begin{equation}\tag{2} w_{t+1} = w_{t} -\eta_{t} \nabla f(w_{t}), \end{equation}

que usualmente se inicializa en cero y cada iteración, como se puede observar, se hace en la dirección negativa del gradiente. Recuerde que el gradiente se define como el vector:

\begin{equation}\nonumber \nabla f(w)=\left(\frac{\partial f}{\partial x_{1}}(w),\dots, \frac{\partial f}{\partial x_{n}}(w)\right). \end{equation}

Un hiperparámetro importante en GD es el $\eta_{t} > 0$, denominado como tasa de aprendizaje. Si la tasa de aprendizaje es demasiado pequeña, entonces el algoritmo tendrá que pasar por muchas iteraciones para converger, lo que llevará mucho tiempo. Por otro lado, si la tasa de aprendizaje es demasiado alta, es posible que se salte el mínimo global y termine en otro lado, posiblemente incluso más alto que antes. Esto podría hacer que el algoritmo diverja, con valores cada vez mayores, sin encontrar una buena solución.

¿Cómo funciona este algoritmo?

Este algoritmo se define a partir de las dos características esenciales que tiene el gradiente, las cuales se mencionan a continuación:

  1. El gradiente es perpendicular a las curvas de nivel de $f$, de manera que para cualquier dirección $v\in \mathbb{R}^{n}$ ortogonal a $\nabla f(p)$, es una dirección de cambio nulo. Esto se observa facilmente al parametrizar la curva $S_{k}=\{p\in \Omega: f(p)=k\}$ mediante una función $\alpha:I\subset \mathbb{R}\to S_{k}$ tal que $\alpha(0)=p$, pues al calcular el producto punto de $\nabla f(p)$ con la velocidad de $\alpha$ en $p$ se obtiene que la tasa de cambio es:

    \begin{equation}\nonumber|df_{p}(\alpha'(0))|=|\nabla f(p)\cdot \alpha' (0)| = 0, \end{equation}
    es decir, la tasa de cambio en la dirección de $\alpha'(0)$ es cero.



  1. El gradiente indica la dirección ascendente de la tasa de máximo cambio de $f$ en el punto $p$. La tasa máxima se calcula como $||\nabla f(p)||$. La razón de esto, se aprecia cuando se considera un vector $v\in \mathbb{R}$ tal que $||v||=1$, de manera que para este vector la tasa de cambio es:

    \begin{equation}\nonumber |df_{p}(v)|=||\nabla f(p)||||v|||cos\theta|\leq ||\nabla f(p)|| \end{equation}
    Dicha magnitud es máxima cuando $\theta = 2n\pi$ con $n\in \mathbb{Z}$, es decir, para que $|df_{p}(v)|$ sea máxima, los vectores $\nabla f(p)$ y $v$ deben ser paralelos, de esta manera, la función $f$ crece más rápidamente en la dirección del vector $\nabla f(p)$ y decrece más rápidamente en la dirección de $-\nabla f(p)$, en efecto, si $v=\frac{\nabla f(p)}{||\nabla f(p)||}$, entonces $df_{p}(v)=||\nabla f(p)||$.

La iteración definida en la Ecuación 2 permite construir una sucesión de puntos $\{w_t\}_{t\in [m]_{\mathbb{N}_0}}$ de tal manera que $f(w_{t+1})<f(w_{t})$. Este hecho se puede evidenciar mediante el polinomio de Taylor, en efecto si se considera un punto inicial $w_{0}$, entonces para el primer termino de la expansión de Taylor alrededor de $w_{0}$ se tiene:

\begin{equation}\nonumber f(w_{1})-f(w_{0})\approx \langle w_{1}-w_{0}, \nabla f(w_{0})\rangle=-\eta ||\nabla f(w_{o})||^2 \end{equation}

Por consiguiente:

\begin{equation}\nonumber f(w_{1})-f(w_{0}))=-\eta ||\nabla f(w_{o})||^2 + o(\eta) \end{equation}

de tal manera que para un $\eta$ adecuado se puede garantizar que $f(w_{1})<f(w_{0})$. El razanomiento se puede repetir para $w_{t+1}$ y $w_{t}$, de tal manera que $w_{t+1}$ es una mejora de $w_{t}$.

De acuerdo a las consideraciones anteriores, el algoritmo del gradiente descendente se define de la siguiente manera:


Algoritmo del gradiente descendente (AGD)

Input: $w_0$, $m$, $\eta$, $\nabla f(w)$ 
1.     for $k=0$ to $m$ do: 
2.          $w\leftarrow w - \eta \nabla f(w)$ 
3.     end 
return: $w$

donde $w_0$ es la condición inicial del algoritmo, $m$ es el número máximo de iteraciones, $\eta$ es la tasa de aprendizaje y $\nabla f(w)$ es la función gradiente de $f$. Observe que la regla de actualización se definió apartir de la Ecuación 2. Es importante notar que finalmente $w_{t+1}$ es tal que:

\begin{equation}\nonumber w_{t+1}\in \operatorname*{argmin\,\,}_{ w\in \Omega} \frac{1}{2\eta_{t}}||w-w_{t}+\eta_{t}\nabla f(w_{t})||. \end{equation}

Es fácil comprobar que el problema de optimización anterior se puede reescribir como:

\begin{equation} \operatorname*{argmin\,\,}_{ w\in \Omega} \frac{1}{2\eta_{t}}||w-w_{t}+\eta_{t}\nabla f(w_{t})||=\operatorname*{argmin\,\,}_{ w\in \Omega} \left(f(w_{t}) + \langle \nabla f(w_t), w-w_t \rangle + \frac{1}{2\eta_{t}}||w-w_{t}||\right) \end{equation}

Por lo tanto, el $w_{t+1}$ es obtenido para minimizar la linealización de la función $f$ alrededor del punto $w_{t}$, manteniendolo lo suficientemente aproximado a este el punto $w_{t}$.

¿Cuándo parar las iteraciones?

Es importante aclarar, que a pesar de que el algoritmo ejecute el máximo de iteraciones, el resultado arrojado no necesariamente es una buena aproximación a el elemento minimizador de la función $f$. Por lo tanto es necesario definir un criterio que permita decidir si el resultado obteniendo es adecuado o no.

Como primer criterio de parada del algoritmo que se le puede ocurrir al lector es detener las iteraciones cuando $||\nabla f(x_{t})||=0$, pero esto no es práctico debido a diferentes factores que influyen, como el comportamiento de la coma flotante y la elección adecuada de la tasa de aprendizaje. Usualmente en la práctica se suele definir un parametro de tolerancia $\epsilon>0$ que junto con alguno de los siguientes criterios se usan para detener las iteraciones:

  • Condición sobre el gradiente:

    \begin{equation}\nonumber ||\nabla f(x_{t})||<\epsilon\end{equation}
  • Condición sobre las diferencias sucesivas relativas de la función objetivo:

    \begin{equation}\nonumber \frac{|f(x_{t+1})-f(x_{t})|}{|f(x_{t})|}<\epsilon,\end{equation}
    si el denominador es muy pequeño, es conveniente remplazarlo por $\max\{1, |f(x)|\}$.
  • Condición sobre las diferencias sucesivas relativas de la variable independiente:

    \begin{equation}\nonumber \frac{||x_{t+1}-x_{t}||}{||x_{t}||}<\epsilon,\end{equation}
    si el denominador es muy pequeño, es conveniente remplazarlo por $\max\{1, ||x_{t}||\}$.

Algunas consideraciones sobre las tasas de aprendizaje.

  • La principal desventaja del AGD se encuentra en el ajuste adecuado de la tasa de aprendizaje $\eta$. Si $\eta$ toma un valor muy pequeño, es necesario un gran número de iteraciones para que el proceso converga; si por otro lado $\eta$ es muy grande, entonces puede ocurrir que el proceso no converga.

  • La tasa de aprendizaje $\eta$ es determinada por la minimización exacta de:

    \begin{equation}\nonumber\eta_{t} \in \operatorname*{argmin\,\,}_{\eta>0}f(x_{t}-\eta\nabla f(x_{t})). \end{equation}
    Esto se usa principalmente para problemas de caracter cuadrático y en donde el cálculo de $\eta$ es económico pero una evaluación de gradiente costosa; de lo contrario, no vale la pena el esfuerzo de resolver este subproblema exactamente.

Ejemplo: Gradiente descente para una forma cuadrática

Asuma que $Q$ es simétrica y definida positiva ($x^{\top}Q x>0$ para cualquier $x\neq 0$). Considere la forma cuadrática:

\begin{equation}\nonumber f(x)=\frac{1}{2}x^{\top}Qx - b^{\top}x \end{equation}

el lector puede comprobar que su gradiente es:

\begin{equation}\nonumber \nabla f(x)=Qx-b. \end{equation}

Así la secuencia de $\{x_{t}\}_{t\in [m]}$ que inicia en cualquier $x_{0}$ viene dada por:

\begin{equation}\nonumber x_{t+1}=x_{t}-\eta_{k}(Qx-b) \end{equation}

con $g_{t}:=\nabla f(x_t)$ se define:

\begin{equation}\nonumber \eta_{t}=\frac{g_{t}^{\top}g_{t}}{g_{t}^{\top}Qg_{t}}. \end{equation}

El lector puede validar que con el valor de $\eta_{t}$ definido anteriormente se tiene:

\begin{equation}\nonumber \eta_{t}\in \operatorname*{argmin\,\,}_{\eta>0}f(x_{t}-\eta\nabla f(x_{t})). \end{equation}

¿Cómo implementarlo en python?

Para las consideraciones del ejemplo anterior, es fácil definir una clase en Python para todas las formas cuadraticas, junto con tres métodos principales que permiten evaluar la forma cuadrática, calcular $\eta$ y el gradiente en un punto $x$. Esto sería algo así:

In [1]:
import numpy as np
class QuadraticForm:
def __init__(self, Q, b):
"""
Inputs:
Q: Positive definite symmetric matrix.
b: Rn vector.
"""
self.Q = np.array(Q)
self.b = np.array(b).reshape(-1, 1)
def evaluate(self, x):
"""
Method to evaluate the quadratic form.
Inputs:
x: Rn vector.
Ouput:
Value of the quadratic form in x.
"""
x = np.array(x).reshape(-1, 1)
Q = self.Q
b = self.b
return (1/2) * x.T.dot(Q.dot(x)) - b.T.dot(x)
def eta(self, x):
"""
Method to evaluate eta.
Inputs:
x: Rn vector.
Output:
Value of eta in x.
"""
gradient_x = self.gradient(x)
numerator = gradient_x.T.dot(gradient_x)
denominator = gradient_x.T.dot(Q.dot(gradient_x))
return numerator / denominator
def gradient(self, x):
"""
Method to evaluate the gradient.
Inputs:
x: Rn vector.
Output:
Value of the gradient in x.
"""
x = np.array(x).reshape(-1, 1)
return Q.dot(x) - b

El AGD se implementaría de la siguiente forma:

In [2]:
def gradient_descent(initial_x, eta, epsilon, function):
"""
Gradient descent Algorithm.
Inputs:
initial_x: Rn vector.
eta: learning rate.
epsilon: precision.
function: Quadractic Form.
Output:
Point where the function reaches the minimum.
"""
k = 0
x = np.array(initial_x).reshape(-1, 1)
while True:
gradient_x = function.gradient(x)
if np.linalg.norm(gradient_x) < epsilon:
print('stop: {}'.format(k))
break
if isinstance(eta, (int, float)):
x = x - eta * gradient_x
else:
x = x - eta(x) * gradient_x
k += 1
return x

Para ejemplicar el funcionamiento del código anterior, se considera la matriz simétrica y positiva definida:

\begin{equation}\nonumber Q = \left(\begin{array}{cc} 1 & 0.5 \\ 0.5 & 3 \end{array}\right) \end{equation}

y vector $b$ dado por:

\begin{equation}\nonumber b = \left(\begin{array}{c} 3 \\ 0.5 \end{array}\right) \end{equation}

De esta manera la forma cuadrática en Python y usando código anterior quedaría así:

In [3]:
Q = np.array([[1, 0.5], [0.5, 3]])
b = np.array([3, 0.5]).reshape(-1, 1)
In [4]:
quadratic_form = QuadraticForm(Q, b) 

Para encontrar el $x$ minimizador de la función se tiene dos opciones, una es definir el parámetro $\eta$ manualmente, y la otra es usar el método eta de la clase QuadracticForm. A continuación se hará uso de los dos casos.

Para empezar se ejecuta el algoritmo en algún punto $x$:

In [5]:
initial_x = np.array([105.5, 105.8])

Luego ejecuta el AGD con $\eta = 0.0001$ y $\epsilon=0.0000001$:

In [6]:
gradient_descent(initial_x, eta=0.0001, epsilon=0.0000001, function=quadratic_form)
stop: 230300
Out[6]:
array([[ 3.18181829],
[-0.36363639]])

Como el lector podrá notar el algoritmo tardó 230300 iteraciones para obtener un candidato al mínimo con la precisión deseada. A continuación se ejecuta usando el método eta de la forma cuadrática:

In [7]:
gradient_descent(initial_x, eta=quadratic_form.eta, epsilon=0.0000001, function=quadratic_form)
stop: 15
Out[7]:
array([[ 3.1818182 ],
[-0.36363637]])

En esta ocasión el algoritmo tardó 15 iteraciones. Esto representa una mejora en tiempo de ejecución bastante considerable con respecto a el experimento anterior.

¿Cómo se puede asegurar que este es el mínimo de la función? Para este caso en particular, el mínimo ocurre cuando $Qx = b$, por lo tanto es suficiente con resolver este sistema lineal. Usando los métodos de la librería numpy se puede resolver rápidamente así:

In [8]:
np.linalg.solve(Q, b)
Out[8]:
array([[ 3.18181818],
[-0.36363636]])

Como se puede observar, el resultado es muy aproximado al valor que se obtuvo al ejecutar AGD, por lo que el algoritmo funciona bastante bien.

Conclusiones

  • Se logró aprender que el AGD es un algoritmo iterativo empleado principalmente para resolver problemas de optimización.
  • La principal desventaja del AGD se encuentra en el ajuste adecuado de la tasa de aprendizaje $\eta$. Si $\eta$ toma un valor muy pequeño, es necesario un gran número de iteraciones para que el proceso converga; si por otro lado $\eta$ es muy grande, entonces puede ocurrir que el proceso no converga.
  • Finalmente te invito a leer este otro artículo (Jugando con el gradiente descendente y Python) en donde se análiza el comportamiento general de los algoritmos por gradiente descendente.

No olvides comentar y suscribirte al blog para que estés enterando de los posts que voy a ir subiendo semana a semana. También sientete en libertad de seguirme en LinkedIn, Twitter, Github e Instagram.

Notas

[1]. En el contexto de este artículo, se dira que una función $f:\Omega \subset \mathbb{R}^{m}\to \mathbb{R}$ es diferenciable si tiene derivada continua en $\Omega$ y es convexa, si para todo $w, z\in \Omega$ y $\alpha \in [0, 1]$ se cumple la siguiente condición:

\begin{equation}\nonumber f(\alpha w + (1-\alpha)z)\leq \alpha f(w) + (1-\alpha)f(z). \end{equation}

Note que la condición es valida para todo $m\geq 1$.



[2]. Si $y_{\min}$ es el mínimo global de una función $f(x)$, entonces se define el operador $\operatorname*{argmin\,\,}$ como siguiente conjunto: $$\operatorname*{argmin\,\,}_{ x \in X} f(x) = \{x\in X: f(x) = y_\min\}.$$

Referencias

Contacto

  • Participa de la canal de Nerve a través de Discord.
  • Se quieres conocer más acerca de este tema me puedes contactar a través de Classgap.

22 may 2018

Los cuatro espacios fundamentales de una matriz

Cada que se considera una matriz $A\in \mathbb{F}^{n,m}$, donde $\mathbb{F}$ puede ser $\mathbb{R}$ o $\mathbb{C}$, esta induce naturalmente dos aplicaciones lineales: $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$ y $A^{\top}:\mathbb{F}^{n}\to \mathbb{F}^{m}$; las cuales permiten determinar cuatro subespacios vectoriales, dos subespacios vectoriales en $\mathbb{F}^{n}$ y dos subespacios vectoriales en $\mathbb{F}^{m}$, conocidos usualmente como espacio nulo, espacio fila, espacio nulo izquierdo y espacio columna de $A$ respectivamente; por lo tanto, el objetivo en esta ocasión es estudiar estos conjuntos y las relaciones entre ellos.

Para comenzar, veamos la definición de espacio nulo y de rango de una matriz $A$: el espacio nulo de la aplicación $A:\mathbb{F}^m\to \mathbb{F}^n$, es el conjunto $N(A)\subseteq \mathbb{F}^{m}$ dado por \begin{equation} N(A)=\{x\in \mathbb{F}^{n}:Ax=0\} \end{equation}; y el espacio columna o rango, es el conjunto $R(A)\subseteq \mathbb{F}^{n}$ dado por \begin{equation} R(A)=\{y\in \mathbb{F}^{n}:y=Ax, \forall\,x\in \mathbb{F}^{m}\}. \end{equation}

En la Figura (1) se muestra el espacio nulo y el espacio columna asociados con la matriz $A$. En este gráfico, se puede apreciar que para la función $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$ el espacio nulo es un subconjunto de $\mathbb{F}^{m}$; mientras que, el espacio columna es un subconjunto de $\mathbb{F}^{n}$. Recuerda que la ecuación homogénea $Ax=0$, siempre tiene como solución trivial $x=0$. Esta propiedad implica que, $O_{_{\mathbb{F}^{m}}}\in \mathbb{F}^{m}$ también pertenece a $N(A)$ y $O_{_{\mathbb{F}^{n}}}\in \mathbb{F}^{n}$ también pertenece a $R(A)$.

Figura 1. El espacio nulo y el rango de la aplicación $A:\mathbb{F}^m\to \mathbb{F}^n$.

El lector puede observar que, realmente $A$ tiene «cuatro subespacios vectoriales» asociados que son $N(A)$, $R(A)$, $N(A^{\top})$ y $R(A^\top)$. Observe que, el espacio nulo y el espacio columna son subespacios de espacios diferentes; más precisamente, una matriz $A\in \mathbb{F}^{m,n}$ define las aplicaciones $A:\mathbb{F}^{m}\to \mathbb{F}^{n}$, $A^{\top}:\mathbb{F}^{n}\to \mathbb{F}^{m}$ y los conjuntos $N(A)$, $R(A^\top)\subseteq \mathbb{F}^{m}$ y $N(A^\top), R(A)\subseteq \mathbb{F}^{n}$.

En efecto, para ver que $N(A)$, $R(A)$, $N(A^{\top})$ y $R(A^\top)$ son subespacios vectoriales, recordemos que un subconjunto $U\subset \mathbb{F}^{m}$ es cerrado bajo combinaciones lineales, si para cada par de elementos $x,y\in U$ y todos los escalares $a, b\in \mathbb{F}$ se tiene que $ax+by\in U$.

Los conjuntos $N(A)$ y $R(A)$ son cerrados bajo combinaciones lineales; dado que, la función inducida por $A$ es una transformación lineal. En efecto, considere dos elementos arbitrarios $x_{_1}, x_{_2}\in N(A)$, entonces $Ax_{_1}=0$ y $Ax_{_2}=0$. Adicionalmente, para cualquier par de elementos $a,b\in \mathbb{F}$ se tiene: \begin{equation} A(ax_{_1}+bx_{_2})=aAx_{_1}+bAx_{_2}=0 \end{equation} de modo que, $(ax_{_1}+bx_{_2})\in N(A)$. Por consiguiente, $N(A)\subseteq \mathbb{F}^m$ es cerrado bajo combinaciones lineales. Análogamente, considere dos elementos $y_{_1}, y_{_2}\in R(A)$; esto es, existe $x_{_1}, x_{_2}\in \mathbb{F}^{m}$ tal que, $y_{_1}=Ax_{_1}$ y $y_{_2}=Ax_{_2}$. Entonces, para cualquier $a, b\in \mathbb{F}$ se tiene: \begin{equation} (ay_{_1}+by_{_2})=aAx_{_1}+bAx_{_2}=A(ax_{_1}+bx_{_2}) \Rightarrow (ay_{_1}+by_{_2})\in R(A). \end{equation} Por lo tanto, $(A)\subset \mathbb{K}^{n}$ es cerrado bajo combinaciones lineales.

En consecuencia, si $A=[A_{_{:1}},\dots, A_{_{:n}}]$, cualquier elemento $y\in R(A)$ puede ser expresado como $y=Ax$ para algún $x\in \mathbb{F}^{n}$; esto es, \begin{equation} y=Ax =[A_{_{:1}},\dots, A_{_{:n}}]\left[\begin{array}{c} x_{_1}\\ \vdots \\ x_{_n} \end{array}\right]=A_{_{:1}}x_{_1}+\cdots + A_{_{:n}}x_{_n}\in gen(\left\{A_{_{:1}}x_{_1},\dots, A_{_{:n}}x_{_n}\right\}). \end{equation} Es por esta razón que, el espacio columna y el rango son el mismo espacio. El lector puede verificar que cualquier elemento de la forma $y=A_{_{:1}}x_{_1}+\cdots + A_{_{:n}}x_{_n}$ es también un elemento de $R(A)$; por lo tanto, $y=Ax$.

El espacio nulo y el rango de una matriz cuadrada, caracterizan si la matriz es invertible o no; esto se establece en el siguiente resultado:

Teorema 1. Dada una matriz $A\in \mathbb{F}^{n,n}$, entonces las siguientes afirmaciones son equivalentes:
  • La matriz $A^{-1}$ existe;
  • $N(A)=\{0\}$;
  • $R(A)=\mathbb{F}^{n}$.

Veamos ahora algunas de las relaciones que existen entre los cuatro espacios asociados a un par de matrices $A$ y $B$, tales que, la matriz $B$ se pueda obtener mediante operaciones de fila sobre la matriz $A$. Asumimos entonces la siguiente notación: $A\stackrel{fila}{\longleftrightarrow}B$, para indicar que la matriz $A$ se puede transformar en la matriz $B$ mediante operaciones de Gauss sobre las filas de $A$. En tal caso, se tiene el siguiente teorema:

Teorema 2. Sean las matrices $A, B \in \mathbb{R}^{m,n}$, entonces:
  • $A\stackrel{fila}{\longleftrightarrow}B\Longleftrightarrow N(A)=N(B)$;
  • $A\stackrel{fila}{\longleftrightarrow}B\Longleftrightarrow R(A^\top)=R(B^\top)$.

Es fácil ver que este último resultado es verdadero; dado que, las operaciones de Gauss no cambian las soluciones de sistemas lineales; por lo cual, el sistema lineal $Ax = 0$ tiene exactamente las mismas soluciones de $x$ como el sistema lineal $Bx = 0$; es decir, $N (A) = N (B)$. La segunda propiedad también es cierta; ya que, las operaciones de Gauss en las filas de $A$, son equivalentes a las operaciones de Gauss sobre las columnas de $A^{\top}$. Ahora bien, es fácil ver que cada una de las operaciones de Gauss sobre las columnas de $A^{\top}$ no cambian el $R(A^{\top})$; por lo tanto, $R(A^{\top})=R(B^{\top})$.

Sin embargo, no podemos pensar que el párrafo anterior es suficiente para hacer una demostración del teorema. Veamos una presentación más detallada de las ideas anteriores. Una forma de hacerlo, es utilizar la multiplicación de matrices para expresar la propiedad en la cual las operaciones de Gauss no cambian las soluciones de sistemas lineales. Con esto, se puede probar lo siguiente: Si se tienen las matrices $A$ y $B$ de orden $n\times m$ relacionadas por las operaciones de Gauss en sus filas, entonces existe una matriz $G$ de orden $m × m$ invertible $GA = B$. La prueba de esta propiedad es simple; debido a que cada una de las operaciones de Gauss se asocia con una matriz invertible, $E$, llamada una matriz de Gauss elemental. Cada matriz de Gauss elemental es invertible, ya que cada operación de Gauss siempre se puede revertir. El resultado de varias operaciones de Gauss sobre una matriz $A$, es el producto de las matrices de Gauss elementales apropiadas en el mismo orden que se realizan las operaciones de Gauss. Si se obtiene la matriz $B$ de la matriz $A$, haciendo operaciones de Gauss dadas por matrices $E_{_i}$, para $i = 1,\dots k$, en ese orden, podemos expresar el resultado del método de Gauss de la siguiente manera: \begin{equation} E_{_k}\cdots E_{_1}A=B\hspace{0.5cm} G=E_{_k}\cdots E_{_1} \Longrightarrow GA =B. \end{equation} Donde cada matriz elemental de Gauss es una matriz invertible; y por lo tanto, $G$ también es invertible.

Considere dos matrices $A$ y $B$ de orden $m\times n$ que están relacionadas por operaciones de Gauss en sus filas; siendo así, existe una matriz $G$ de orden $m\times m$, tal que, $GA = B$. Esta observación es la clave para mostrar que $N(A) = N(B)$, ya que dado cualquier elemento $x \in N(A)$ \begin{equation} Ax=0\Longleftrightarrow GAx=0, \end{equation} donde la equivalencia se sigue del hecho de que $G$ es invertible. Entonces es simple (¡Lo simple es una provocación para que tú lo verifiques!) ver que: \begin{equation} 0=GAx=Bx \Longleftrightarrow x\in N(B). \end{equation} Así pues, se tiene que $N(A)=N(B)$. Ahora veamos la afirmación opuesta: Si $N(A)=N(B)$, significa que sus formas escalonadas reducidas $E_{_A}, E_{_B}$ son las mismas; es decir, $E_{_A}=E_{_B}$; esto significa que existen operaciones de Gauss en las filas $A$ que la transforman en la matriz $B$

Ahora mostraremos que $R(A^{\top})= R(B^{\top})$. Considere un elemento $x\in R(A^{\top})$, por lo tanto, existe un elemento $y\in \mathbb{F}^{m}$, tal que: \begin{equation} x=A^{\top}y = A^{\top}G^{\top}(G^{\top})^{-1}y=(GA)^{\top}\bar{y}=B^{\top}\bar{y},\hspace{0.5cm} \bar{y}=(G^{\top})^{-1}\bar{y} \end{equation} Veamos que dado un $x\in R(A^{\top})$, se puede decir que, $x\in R(B^{\top})$; esto es, $R(A^{\top})\subset R(B^{\top})$. La implicación opuesta se prueba de igual forma: Considere $x\in R(B^{\top})$, por ende existe $\bar{y}\in \mathbb{F}^{m}$, tal que: \begin{equation} x=B^{\top}\bar{y}=B^{\top}(G^{\top})^{-1}G^{\top}\bar{y}=(G^{-1}B)^{\top}y=A^{\top}y,\hspace{0.5cm} y=G^{\top}\bar{y}. \end{equation} De modo que, se ha mostrado que para cualquier $x\in R(B^{\top})$, por lo cual $x\in R(A^{\top})$; es decir, $R(B^{\top})\subset R(A^{\top})$; y en consecuencia, $R(A^{\top})=R(B^{\top})$. Ahora sí se asume que $R(A^{\top})=R(B^{\top})$; esto quiere decir que, cada fila de $A$ es una combinación lineal de las filas de $B$; que a su vez, significa que existen operaciones de Gauss sobre las filas de $A$ que transforman a $A$ en $B$; y por lo tanto, con eso se establece el resultado final del teorema 2.

Un argumento similar también dice que, $E_{_A^{\top}}=E_{_{B^{\top}}}$ sí y sólo sí $A^{\top}\stackrel{fila}{\longleftrightarrow} B^{\top}$. También se puede concluir que, $E_{_{A^{\top}}}=E_{_{B^{\top}}}$ es equivalente a $R(A)=R(B)$ y esto también es equivalente a $N(A^{\top})=N(B^{\top})$.

Otro resultado con respecto a la transpuesta de $A$ dice:

Teorema 3.Para cada matriz $A\in \mathbb{F}^{n,m}$ se cumple que $\dim R(A)=\dim R(A^{\top})$.
Sabiendo que una matriz se dice que es de rango completo, sí y sólo sí, $\dim R(A)=\min(m,n)$. Entonces podemos enunciar el siguiente teorema donde se establecen algunas de las relaciones principales de los cuatro espacios fundamentales:
Teorema 4. Si una matriz $A\in \mathbb{F}^{n,m}$ es de rango completo, entonces:
  • Si $m=n$, por lo tanto: \[\dim R(A)=\dim R(A^{\top})=n=m\Leftrightarrow \{0_{_{\mathbb{F}^{m}}}\}=N(A)=N(A^{\top})\subset \mathbb{F}^{m};\]
  • Si $n < m$, de forma que: \[\dim R(A)=\dim R(A^{\top})=n < m \Leftrightarrow \left\{\begin{array}{l} \{0_{_{\mathbb{F}^{m}}}\}\varsubsetneq N(A)\subset \mathbb{F}^{m},\\ \{0_{_{\mathbb{F}^{n}}}\}=N(A^{\top})\subset \mathbb{F}^{n};\end{array}\right.\]
  • Si $m>n$, siendo así: \[ \dim R(A)=\dim R(A^{\top})= m < n \Leftrightarrow \left\{\begin{array}{l} \{0_{_{\mathbb{F}^{m}}}\}= N(A)\subset \mathbb{F}^{m},\\ \{0_{_{\mathbb{F}^{n}}}\}\varsubsetneq N(A^{\top})\subset \mathbb{F}^{n}; \end{array}\right.\]

Recordemos que, el rango de una matriz $A$ es el número de pivotes columna de la matriz escalonada reducida $E_{_A}$ y que coincide con la dimensión de $R(A)$. Si una matriz $A$ de orden $m\times n$ tiene $rang(A)=n$, esto quiere decir dos cosas: Primero, $n\leq m$; y segundo, que cada columna de $E_{_A}$ tiene un pivote; es decir, que no hay variables libres en la solución de la ecuación $Ax=0$, y así $x=0$ es la única solución. Debido a esto, $N(A)=0$. Por otro lado, al estudiar $N(A^{\top})$ es necesario considerar dos casos: $n=m$ o $n < m$. Si $n=m$, entonces las matrices son cuadradas; debido a esto, se puede decir que no hay variables libres en la ecuación $A^{\top}y=0$; y por lo tanto, se concluye que $N(A^{\top})=0_{_{\mathbb{F}^{n}}}$. Por otra parte, si $n < m$ entonces hay variables libres en la solución de la ecuación $A^{\top}y=0$; y por ende, $\{0_{_{\mathbb{F}^{n}}}\}\varsubsetneq N(A^{\top})$. Si se tiene una matriz $A$ de orden $m\times n$ con $rang(A)=m$, note que $rang(A)=rang(A^{\top})$; esto quiere decir dos cosas: Primero que $m \leq n$; y segundo que cada columna de $E_{_A^{\top}}$ tiene un pivote. Esta última afirmación muestra que $A^{\top}$ es de rango completo; es decir, $N(A^{\top})=0$. Ya se han analizado todos los casos en que $n=m$; sólo falta analizar qué pasa cuando $m < n$. En este caso, hay variables libres en la solución de la ecuación $Ax=0$; por lo tanto, $\{0_{_{\mathbb{F}^{m}}}\}\varsubsetneq N(A)$.

Figura 2. Relaciones entre los cuatro espacios fundamentales de $A:\mathbb{F}^m\to\mathbb{F}^n$.

Para terminar, enunciamos el siguiente resultado que relaciona las dimensiones de los espacios nulos y el rango de una transformación lineal en espacios vectoriales de dimensión finita. Este resultado se suele llamar Teorema de la Nulidad y el Rango, donde la nulidad de una transformación lineal es la dimensión de su espacio nulo, y el rango es la dimensión de su espacio de columna. Este resultado también se le conoce como el Teorema de la Dimensión.

Teorema 4. Para cada transformación lineal $T:V\to W$ entre espacios vectoriales $V$ y $W$ se cumple que: \begin{equation} \dim N(T)+\dim R(T) = \dim V. \end{equation}

Si consideramos el producto punto de dos vectores $u, v\in \mathbb{F}^{m}$ como el producto matricial definido por $u\cdot v =u^{\top}v$; podemos observar lo siguiente: \begin{equation} Ax=\left[\begin{array}{c} A_{_{1:}}x\\ \vdots \\ A_{_{m:}}x \end{array}\right] \end{equation} Note que las filas $A$ son las columnas de $A^{\top}$, por lo tanto, $A^{\top}_{_{i:}}\cdot x= A_{_{i*}}x$ Si $x\in N(A)$, y por ende, se puede concluir que: $N(A)$ es el complemento ortogonal de $R(A^{\top})$,y por consiguiente, $N(A)\cap R(A^{\top})=\emptyset$. En otro orden de ideas, si la dimensión de $\dim R(A)=r$, entonces la dimensión del espacio nulo $\dim N(A) = m-r$; y $\dim R(A^{\top})=r$, luego $\dim N(A^{\top})=n-r$. Note que, $\dim N(A)\oplus R(A^{\top}) = m$ y $\dim N(A^{\top})\oplus R(A) = n$; es decir que, $ N(A)\oplus R(A^{\top})\cong \mathbb{F}^{m}$ y $N(A^{\top})\oplus R(A) \cong \mathbb{F}^{n}$. Estas relaciones se resumen en la Figura 2.

Bueno, ya nos hemos extendido lo suficiente por esta ocasión. Esperamos que aprovechen mucho este post.

Referencias

21 may 2018

Números de coma o punto flotante y sistemas lineales

En esta ocasión, estudiaremos algunos aspectos generales de los números de punto flotante y algunas dificultades que se presentan cuando se realizan cálculos elementales en la resolución de sistemas lineales por eliminación Gaussiana y Gauss-Jordan. Esperamos que lo disfruten mucho.

Los números de coma o punto flotante, son un conjunto finito de números racionales que se emplean para representar números reales empleando computadoras. Existen diferentes tipos de números de punto flotante; todos ellos, se caracterizan porque tienen un número finito de dígitos cuando se escriben en una base particular. La necesidad de estos números es precisamente porque las computadoras sólo pueden representar los números reales con un conjunto finito de dígitos.

Definición. Un número racional \(x\) es un número de punto flotante en base $b\in \mathbb{N}$, de precisión $p\in \mathbb{N}$, con rango de exponencial $N\in \mathbb{Z}$, sí y sólo sí, existen enteros $d_{_i}$, para $i=1,\dots, b-1$ y $d_{_1}\neq 0$ tal que, $x$ tiene la forma: \begin{equation}\label{eq:01} x=\pm 0.d_{_1}\cdots d_{_p}\times b^{n},\hbox{ con } -N\leq n\leq - N. \end{equation} Se denota por $\mathbb{F}_{p, b, N}$ el conjunto de todos los números de punto flotante con precisión $p$, base $b$ y rango exponencial $N$.

Lo primero que se puede observar es que el conjunto $\mathbb{F}_{p, b, N}$ es un conjunto finito, y que no hay una distribución homogénea de sus elementos. Por ejemplo, si consideramos el conjunto de números flotantes de precisión $1$, base $10$ y rango exponencial $10$, como el lector podrá notar, es fácil listar los elementos positivos de este conjunto. En efecto sus elementos son: \[\{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09\}=\{0.i\times 10^{-1}\}_{i=1}^{9},\] \[\{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9\}=\{0.i\times 10^{0}\}_{i=1}^{9},\] \[\{1, 2, 3, 4, 5, 6, 7, 8, 9\}=\{0.i\times 10^{0}\}_{i=1}^{9}.\]

Como se puede observar, los elementos de $\mathbb{F}_{1, 10, 1}$ no se encuentran homogéneamente distribuidos sobre el intervalo $[0,10]$. Esta distribución irregular afecta notablemente la posibilidad de hacer cálculos de adición y multiplicación entre números pequeños y números grandes, o solamente entre números grandes. Observe por ejemplo que, cuando sumamos $0.01$ con $9$ obtenemos $9.01$, que no pertenece a $\mathbb{F}_{1, 10, 1}$; igualmente, se puede verificar que el producto de $8$ por $9$ tampoco pertenece.


Figura 1. Distribución no homogénea de los elementos de $\mathbb{F}_{1, 10, 1}$.

En general, los conjuntos de números flotantes $\mathbb{F}_{p, b, N}\subset\mathbb{R}$ no son cerrados bajo la suma o la multiplicación. Por lo tanto, al tratar de representar números reales en algún conjunto $\mathbb{F}_{p, b, N}\subset\mathbb{R}$, habrán ciertos cálculos aritméticos que no serán permitidos. Una manera de realizar cálculos con números reales en $\mathbb{F}_{p, b, N}\subset\mathbb{R}$ es, primero proyectar los números reales en los números de punto flotante y luego realizar el cálculo; ya que, el resultado podría no estar en el conjunto $\mathbb{F}_{_{p, b, N}}\subset\mathbb{R}$, uno debe proyectar de nuevo el resultado en $\mathbb{F}_{p, b, N}\subset\mathbb{R}$. La acción de proyectar un número real en el conjunto de números flotantes, es llamada redondear un número real. Existen muchas formas de hacer esto. A continuación, presentamos la forma más común de hacerlo.

Definición. Sea $X_{_N}=\max \mathbb{F}_{p, b, N}$; por esto, la función de redondeo se define como una aplicación $fl:\mathbb{R}\cap [-X_{_N},X_{_N}]\to \mathbb{F}_{p, b, N}$ definida como: Dado un número real $x\in \mathbb{R}\cap [-X_{_N},X_{_N}]$, con $x=\pm 0.d_{_1}\cdots d_{_p}d_{_{p+1}}\cdots \times b^{n}$ y $-N\leq n\leq N$, se tiene que: \begin{equation} fl(x)=\left\{\begin{array}{cl} \pm 0.d_{_1}\cdots d_{_p}\times b^n & \hbox{ si }d_{_{p+1}}<\frac{b}{2}, \\ \pm (0.d_{_1}\cdots d_{_p}+b^{-p})\times b^n & \hbox{ si }d_{_{p+1}}\geq \frac{b}{2}.\end{array}\right. \end{equation}

Por ejemplo, si consideramos el conjunto de números flotantes $\mathbb{F}_{3,10,3}$, entonces algunas proyecciones de los números reales serían $fl(0.2103\times 10^{3})=0.210\times 10^3$, $fl(0.21037\times 10^{3})=0.210\times 10^3$, $fl(0.2105\times 10^{3})=0.211\times 10^3$ y $fl(0.2107\times 10^{3})=0.211\times 10^3$. Observe cómo la función $fl$ no es una función inyectiva, y por lo tanto, diferentes números reales pueden ser representados con el mismo número flotante.

Otra consecuencia que cabe mencionar, es que la aritmética de los número reales cambia notablemente; es decir, si $x,y\in \mathbb{R}$, siendo así, la suma de reales sobre un conjunto $\mathbb{F}_{p, b, N}$ queda definida por: \begin{equation} x+_{f}y=fl(fl(x)+fl(y)) \end{equation} y la multiplicación por: \begin{equation} x\cdot_{f}y=fl(fl(x)\cdot fl(y)) \end{equation} Aquí las operaciones $+_{f}$ y $\cdot_{f}$ son la suma y la multiplicación en el conjunto de coma flotante $\mathbb{F}_{p, b, N}$. Éstas operaciones son diferentes de la suma ($+$) y la multiplicación ($\cdot$) usual de números reales. En efecto, se tiene la siguiente proposición:

Proposición. Sea $\mathbb{F}_{p, b, N}$, entonces siempre existen $x,y\in \mathbb{R}$ tales que, \begin{equation} fl(fl(x)+fl(y))\neq fl(x+y),\hspace{0.5cm} fl(fl(x)\cdot fl(y))\neq fl(xy).\end{equation}

Estamos seguros de que el lector puede encontrar algunos ejemplos sencillos que verifican la proposición anterior.

Como se ha visto, las operaciones aritméticas usuales no son posibles, en general, en los conjuntos $\mathbb{F}_{p, b, N}$; debido a que, estos conjuntos no son cerrados bajo estas operaciones. Por lo que, se han definido de una manera un poco artificial las operaciones $+_{f}$ y $\cdot_{f}$. A continuación, veremos algunos efectos de estas operaciones cuando se trata de resolver sistemas de ecuaciones lineales en algún $\mathbb{F}_{p, b, N}$.

Considere el conjunto de números flotantes $\mathbb{F}_{3,10,3}$ para resolver el siguiente sistema de ecuaciones lineales: \[ \begin{array}{rcc} 5x_{_1}+x_{_2} &=& 6,\\ 9.34x_{_1}+1.57x_{_2} &=& 11. \end{array} \] Observe que, al resolver este sistema $\mathbb{R}$, sin las limitaciones de $\mathbb{F}_{3,10,3}$, encontramos que las soluciones son $x_{_1}=x_{_2}=1$. Veamos ahora qué ocurre cuando empleamos $\mathbb{F}_{3,10,3}$ con las operaciones $+_{f}$ y $\cdot_{f}$ al realizar Gauss-Jordan sobre este sistema de ecuaciones. Inicialmente tenemos el sistema matricial: \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\]; para cambiar la segunda fila, podemos hacer las siguientes operaciones sobre las filas: empleamos la suma y la multiplicación de redondeo en $\mathbb{F}_{3,10,3}$. \[\hat{a}_{_{2i}}=fl(a_{_{21}}-fl\left(fl(a_{_{1i}})\cdot fl\left(\frac{9.53}{5}\right)\right)\hbox{ con } i=1,2,3\], donde $\hat{a}_{_{2i}}$ son las entradas de la nueva fila $2$, y $a_{_{1i}}$ son las entradas de la fila $1$ en el sistema inicial. Al hacer estas operaciones se tiene que: \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\rightarrow \left[\begin{array}{cc|c} 5 & 1 & 6\\ 0.02 & -0.32 & -0.3\end{array}\right].\] En este caso no es posible continuar con el método de Gauss-Jordan al menos que, $\hat{a}_{_{21}}=0$; lo cual no es posible. Si introducimos la modificación de que $\hat{a}_{_{21}}=0$, se tiene para nuestro ejemplo que, \[\left[\begin{array}{cc|c} 5 & 1 & 6\\ 9.43 & 1.57 & 11\end{array}\right]\rightarrow \left[\begin{array}{cc|c} 5 & 1 & 6\\ 0 & -0.32 & -0.3\end{array}\right].\] Es importante notar que, aquí no se ha redondeado el error, pues el valor $0.02$ pertenece a $\mathbb{F}_{3,10,3}$. Lo que se ha hecho es una modificación al sistema para poder continuar con el proceso de Gauss-Jordan. Si completamos el proceso bajo las operaciones $+_{f}$ y $\cdot_{f}$ de $\mathbb{F}_{3,10,3}$, se encuentra \[\left[\begin{array}{cc|c} 1 & 0 & 1.01\\ 0 & 1 & 0.938\end{array}\right]\rightarrow \begin{array}{l} x_{_1}=0.101\times 10,\\ x_{_2}=0.938.\end{array}\] Note que la solución bajo $\mathbb{F}_{3,10,3}$, difiere de la solución exacta $x_{_1}=x_{_2}=1$. El error se produce por los redondeos y por la modificación hecha para poder completar el procedimiento de Gauss-Jordan.

En general, los errores por redondeo son muy importantes cuando se hacen sumas entre números pequeños con números grandes o cuando se divide por números muy pequeños. Por ejemplo, considere $x=0.100\times 10^4$ y $y=0.400\times 10$; estos números pertenecen al conjunto $\mathbb{F}_{3,10,4}$; y observe que, $fl(x)=x$ y $fl(y)=y$. Por lo tanto, cuando se hace la suma, se obtiene: \[fl(x+y)=fl(1000+4)=fl(0.1004\times 10^3)=1\times 10^3=x.\] Es decir, la información de $y$ se pierde completamente.

Hay tres estrategias conocidas como el pivoteo parcial, pivoteo parcial escalado y pivoteo completo que buscan evitar la división por números grandes; para así, reducir un poco los errores por redondeo. Considere una matriz $A$ de tamaño $m\times n$ y veamos en qué consisten estos dos métodos:

  • Pivoteo Parcial: En cada paso $k$ de la eliminación Gaussiana, se elige en calidad de pivote a la entrada con índice $(k,p)$ tal que: \begin{equation} |A_{_{kp}}|=\max_{p\leq i \leq m}|A_{_{ip}}|, \end{equation}; es decir, la mayor entrada de la $p$ - ésima columna empezando con $(p,p)$ hasta $(m,p)$. Por ejemplo, supongamos que vamos a realizar el paso $p$ de la eliminación Gaussiana, entonces nuestra matriz tendrá una forma parecida a esto: \[\left[\begin{array}{cccccc|c} * & * & * & & * & & * \\ 0 & * & * & & * & & * \\ 0 & 0 & * & & * & & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots \\ 0 & 0 & 0 & & A_{_{pp}} & & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots \\ 0 & 0 & 0 & & A_{_{mp}} & & * \\ \end{array}\right]_{m\times n}\] Y suponga que nuestro máximo es $|A_{_{kp}}|=\max_{p\leq i \leq m}|A_{_{ip}}|$, luego la entrada $A_{kp}$ se usará como pivote. Esto significa que se aplicará el intercambio de las filas $R_{_p}\leftrightarrow R_{_k}$, y luego se aplican las operaciones elementales para eliminar las entradas desde $(p+1,p)$ hasta $(m,p)$.
  • Pivoteo Parcial Escalado: Suponga que estamos en el paso $p$. En cada renglón $i$, con $p\leq i \leq m$, se calcula el valor máximo absoluto en la parte principal de la matriz: \[s_{_i}=\max_{p\leq j\leq m} |A_{_{ij}}|.\] Suponga que $s_{_i}>0$ para todo $i\in \{k,\dots, m\}$. Si se elige el menor entero $q$ con \[\frac{|A_{_{qk}}|}{s_{_q}}=\max_{p\leq i\leq m}\frac{|A_{_{ip}}|}{s_{_i}}.\] En otras palabras, sea $q$ el menor de los índices $i$ en los cuales la expresión $\frac{|A_{_{ip}}|}{s_{_i}}$ alcanza el máximo. Si $q\neq p$, entonces se intercambian los renglones $p$ y $q$, y luego se eliminan las entradas por debajo de $(p,p)$.
  • Pivoteo Completo: En el $k$ - ésimo paso de la eliminación Gaussina, se buscan los índices $p,q\in \{k,\dots, m\}$ tales que, \[|A_{_{rs}}|=\max_{\substack{p\leq i\leq m \\ p\leq j\leq n}}|A_{_{ij}}|.\] En otras palabras, se busca el máximo entre los números $|A_{_{ij}}|$ con $p\leq i \leq m$ y $p\leq j \leq m$. En este caso, antes de efectuar el paso $p$ de la eliminación Gaussiana, la matriz tendrá una forma como ésta: \[\left[\begin{array}{ccccccc|c} * & * & * & & * & & * & *\\ 0 & * & * & & * & & * & * \\ 0 & 0 & * & & * & & * & *\\ \vdots & \vdots & \vdots & & \vdots & & \vdots & \vdots\\ 0 & 0 & 0 & & A_{_{pp}} & \cdots & A_{_{pn}} & * \\ \vdots & \vdots & \vdots & & \vdots & & \vdots & *\\ 0 & 0 & 0 & & A_{_{mp}} & \cdots & A_{_{mn}} & *\\ \end{array}\right]_{m\times n}\] Si $(r,s)\neq (p,p)$, entonces se intercambian los renglones y las columnas de tal manera que la entrada $A_{_{rs}}$ se pone en la posición $(p,p)$, y luego se aplican las operaciones elementales para anular las entradas por debajo de $(p,p)$.
Esperamos que esta entrada haya sido de su agrado, nos vemos en otra ocasión.

Referencias

13 ene 2017

Red neuronal artificial simple



El objetivo en esta ocasión es introducirnos en el mundo de las redes neuronales. Para empezar comentaremos de una manera general y quizás un poco atrevida el concepto de neurona artificial, para finalmente construir un perceptrón digital simple empleando python 3.

¿Qué es un sistema neuronal artificial?

La idea de los sistemas neuronales artificiales fue inspirada por los sistemas neuronales biológicos. De acuerdo con Ramón y Cajal (1888), un sistema neuronal biológico está compuesto por una red de células individuales, ampliamente interconectadas entre sí. Estas células son denominadas neuronas y son como pequeños procesadores de información. Estos procesadores se encuentra compuesto principalmente por:
  1. Las dendritas, son el canal receptor de la información.
  2. El soma, es el órgano encargado de procesar la información.
  3. El axón, es el canal de emisión de información a otras neuronas.
El cerebro humano tiene cerca 90.000.000.000 neuronas, además cada neurona recibe información de aproximadamente 10.000 neuronas y envía impulsos a cientos de ellas. También hay neuronas que reciben información directamente de exterior. Es importante observar que el cerebro se modela durante el desarrollo del ser vivo, por lo tanto, algunas cualidades no son innatas, sino adquiridas por la influencia de la información que del medio externo recibe.

El objetivo de las redes será construir un gran conjunto de neuronas artificiales para simular un comportamiento similar al del cerebro humano.

El modelo estándar de una neurona artificial

Según los principios descritos por Rumelhart y McClelland (1986). Una neurona artificial estandar tiene los siguientes componentes:

  1. Las dendritas artificiales. Las dentritas artificiales son un conjunto de parámetros $x_{_i}(t)$ con los cuales se codifica la información de un problema que se quiere resolver. Las variables de entrada y salida pueden ser binarias (digitales) o continuas (analógicas), dependiendo del modelo y la aplicación. Por ejemplo en un perceptrón multicapa (MLP), por lo general las salidas son señales digitales representadas por $1$ y $-1$, en el caso de las salidas analógicas, la señal se da en un cierto intervalo.
  2. Los pesos sinápticos $w_{_{ij}}(t)$ de la neurona $i$ son variables relacionadas a la sinapsis o conexión entre neuronas, los cuales representan la intensidad de interacción entre la neurona presináptica $j$ y la postsináptica $i$. Dada una entrada positiva (puede ser una señal proveniente de una neurona), si el peso es positivo tenderá a excitar a la neurona postsináptica, si el peso es negativo tenderá a inhibirla.
  3. La función de potencial, permite obtener a partir de las entradas (dendritas) y los pesos sinápticos, el valor de potencial postsináptico $h_{_i}$ de la neurona $i$ en función de los pesos y entradas \begin{equation} h_{_i}(t)=\sigma_{_i}(w_{_{ij}}(t), x_{_j}(t)). \end{equation} La función más habitual es de tipo lineal, y se basa en la suma ponderada de las entradas con los pesos sinápticos, es decir, \begin{equation} h_{_i}(t) = \sum_{j}w_{_{ij}}(t)x_{_j}(t). \end{equation} Habitualmente se agrega al conjunto de pesos de la neurona un parámetro adicional $\theta_{_i}$, que se denomina umbral de excitación, el cual se acostumbra a restar al potencial postsináptico. Es decir: \begin{equation} h_{_i}(t) = \sum_{j}w_{_{ij}}(t)x_{_j}(t)-\theta_{_i}. \end{equation} Si se tiene un número finito de dendritas y hacemos que los índices $i$ y $j$ comiencen en cero, y denotamos por $w_{_{i0}}=\theta_{_i}$ y $x_{_0}=-1$, la función de potencial lineal se puede expresar como \begin{equation} h_{_i}(t) = \sum_{j=0}^{n}w_{_{ij}}(t)x_{_j}(t)=\pmb{w}^{\top}_{_i}(t)\cdot \pmb{x}(t), \end{equation} con $\pmb{w}_{_i}(t) = (w_{_{i0}}(t),\dots,w_{_{in}}(t))$ y $\pmb{x}(t)=(x_{_0}(t),\dots,x_{_n}(t))$.
  4. La función de activación $f_{_i}$ de la neurona $i$ proporciona el estado de activación actual $a_{_i}(t)$ a partir del potencial postsináptico $h_{_i}(t)$ y del propio estado de activación anterior, $a_{_{i}}(t-1)$, es decir, \begin{equation} a_{_i}(t)=f_{_i}(a_{_i}(t-1), h_{_i}(t)). \end{equation} Sin embargo, en muchos modelos de redes artificiales se considera que el estado actual de la neurona no depende de su estado anterior, sino unicamente del actual, por lo tanto, \begin{equation} a_{_i}(t)=f_{_i}(h_{_i}(t)). \end{equation}
  5. La función de emisión. Esta función proporciona la salida global $y_{_i}(t)$ y es el componente principal del axón artificial de la neurona $i$ en función de su estado de activación actual. Muy frecuentemente la función de emisión es simplemente la identidad $F(x) = x$ de tal modo que el estado de activación de la neurona se considera como la propia señal de la neurona.

Gráficamente, una neuronal artificial se puede representar de la siguiente forma:




Figura 1. Esquema de una neurona artificial.

Los pesos sinápticos, la función de potencial y la función de activación son los componentes que definen el soma artificial de la neurona.

Ejemplo de una red neuronal simple: el perceptrón

En los 50's Frank Rosenblatt propuso una red neuronal denominada perceptrón digital simple. Éste consiste de una o varias neuronas, donde la función de activación para cada neurona es \begin{equation} y = F\left(\sum_{i=1}^{n}w_{_{i}}x_{_i}+\theta(t)\right). \end{equation} Generalmente la función de activación $F$ puede ser lineal, y se dice por lo tanto que la conexión es de lineal, aunque puede ser no lineal. Para los objetivos ilustrativos, vamos a considerar la siguiente función: \begin{equation} F(s) = \left\{\begin{array}{cc} 1 & \mbox{ si } s > 0 \\ -1 & \mbox{ si } s \leq 0 \end{array}\right. \end{equation} Para simplificar la exposición de las ideas, suponga que las señal emitida por la neurona sera $1$ o $-1$, aunque puede ser cualquier otra cosa. Este tipo de neurona se puede usar para tareas de clasificación, es decir, se puede usar para decir si una patrón de entrada pertenece a alguna de las clases definidas por los valores $1$ o $-1$. Si el potencial de entrada es positivo, entonces el patrón se le asignará la etiqueta $1$, sino por el contrario es cero o menor que cero se le asignará la etiqueta de $-1$. Observe que los patrones de entrada siempre se pueden identificar con algún vector de $\mathbb{R}^{n}$, de manera que esta neurona separa a $\mathbb{R}^{n}$ en dos clases mediante un hiperplano dado por la ecuación: \begin{equation} x_{_n} = -\frac{w_{_1}}{w_{_n}}x_{_1}-\frac{w_{_2}}{w_{_n}}x_{_1}-\cdots-\frac{w_{_{n-1}}}{w_{_n}}x_{_1}+\frac{\theta}{w_{_n}}. \end{equation} La anterior función se denomina, función discriminante.

En el caso en que el espacio de patrones de entrada se pueda identificar con $\mathbb{R}^{2}$, la situación se puede representar gráficamente, en este caso, el hiperplano que define las dos clases es la linea recta dada por \begin{equation} w_{_1}x_{_1}+w_{_2}x_{_2}-\theta = 0, \end{equation} la cual se puede escribir como \begin{equation} x_{_2}=\frac{w_{_1}}{w_{_2}}x_{_1}+\frac{\theta}{w_{_2}} = 0, \end{equation} observe que el cociente $\frac{w_{_1}}{w_{_2}}$ determina la pendiente de la recta y $\frac{\theta}{w_{_2}}$ su bias. Note también que el vector $(w_{_1}, w_{_2})$ es siempre perpendicular a recta.


Figura 2. Función discriminante de un perceptrón simple.

Suponga ahora que se tiene un conjunto de datos $S\subset\mathbb{R}^{2}$, y un vector $\pmb{x}\in S$ para el cual se desea obtener la señal $\hat{y}(x)$. Como se ha dicho anteriormente $\hat{y}(x)$ es usualmente es un vector donde cada entrada es $+1$ o $-1$. ¿Cómo aprende el perceptrón a clasificar adecuadamente? Para eso el perceptrón sigue la siguiente rutina de aprendizaje:

  1. Iniciar con un conjunto aleatorio de pesos sinápticos.
  2. Seleccionar un patrón de entrada $x\in S$.
  3. Si $y(x) \neq \hat{y}(x)$, entonces los pesos se modifican de acuerdo a la regla: \begin{equation} \Delta w_{_i}= \hat{y}(x)x_{_i}; \end{equation}.
  4. Volver al paso 2.

El lector podrá verificar que este procedimiento es muy similar a la regla de aprendizaje de Hebb, la única diferencia es que cuando la neurona responde correctamente, los pesos sinápticos no son modificados. Por otro otro lado, $\theta$ como es considerado es el peso sináptico $w_{_0}$ que siempre recibe por la dendrita $x_{0}$ el valor de $-1$. Para el caso de $\theta$, la regla de aprendizaje viene dada por: \begin{equation} \Delta \theta = \left\{\begin{array}{cc} 0 & \mbox{ si el perceptron responde correctamente} \\ \hat{y}(x) & \mbox{ si el perceptron responde incorrectamente} \end{array}\right. \end{equation} Por ahora no entraremos en más detalles teóricos y vamos a ver como tener nuestro propio perceptrón con Python 3.

¿Cómo construir un perceptrón con Python?

El perceptron que vamos a construir consiste de una capa de $n$ neuronas artificiales, cada una para reconocer un único patrón. El objetivo de la operación del perceptrón es aprender una tranformación dada de la forma $\hat{y}:\{1,-1\}^{m}\to \{1,-1\}^{m}$ usando un conjunto de muestras donde cada elemtento es de la forma $(\pmb{x}, \pmb{y})$ donde $\pmb{x},\pmb{y}\in \{1,-1\}^{m}$ y además se le indican cuales son los vectores $\pmb{x}, \pmb{y}$, para esto se usará la función $tanh\,\theta$, de la siguiente forma: \begin{equation} f(t)=tanh(\pmb{w}\cdot \pmb{x}), \end{equation}

Tomaremos como bias para el criterio de desición el valor $\theta = 0.9999999999$ y además vamos a considerar que $\pmb{w}=\pmb{y}$, la razón de esta consideración es debido a que la función $tanh\,\theta$ nos permite decir que tan diferentes son dos vectores, así, si cuando se da un valor de entrada $\pmb{x}$ y la neurona artificial nos entrega el valor correcto de $\pmb{y}$ es porque el patrón $\pmb{x}$ está muy cercado al vector $w$, es decir, que el valor de $tanh,\theta$ es muy cercano a uno. No se preocupen por esto, en otra ocasión explicaré como hacer el entrenamiento de la neurona. Por ahora solo veamos como construir un ejemplo completamente funcional.

Lo primero que haremos en importar las librerías necesarias para nuestro algoritmo.

import numpy as np
import math
import pickle
from itertools import product
import tkinter as tk
from typing import Callable

Se define la clase Neuron con los siguientes métodos:

1. El constructor de la neurona:

    def __init__(self, dendrite: np.array, sensitivity:
                 int=2, dilatation: float=0.5) -> object:

        sensitivity = '9' * int(sensitivity)

        # Attributes of class.
        self.sensitivity = 1 - 100 / int(sensitivity)
        self.dilatation = 1/dilatation
        self.n_weight = len(dendrite)
        self.synaptic_weight = np.ones(self.n_weight)

        # We define the memory of neuron.
        pickle.dump([], open('memory.mem', 'wb'))
        self.memory = pickle.load(open('memory.mem', 'rb'))

2. El soma de la neurona:

    def soma(self, dendrite: np.array, potential_function: Callable=np.dot,
        potential = potential_function(dendrite, self.synaptic_weight)

        potential = potential / self.dilatation
        order = active_function(potential)
        return order

3. Toda neurona necesita una método de aprendizaje. Esto se programa a continuación:

    def learn(self, dendrite: np.array) -> None:

        self.memory.append(dendrite)
        self.memory.reverse()
        print('Signal learned')

3. Y también deba saber olvidar:

    def forget(self) -> None:

        self.memory = []
        pickle.dump(self.memory, open("memory.mem", 'wb'))
        print('Memory deleted')

4. Y por último se programa el axón o función de activación:

    def axon(self, dendrite: np.array) -> np.array:

        candidate_signals = {}

        # Se identifican todas las señales posibles.
        for mem in self.memory:
            self.synaptic_weight = mem
            weight = self.soma(dendrite)
            if weight > self.sensitivity:
                candidate_signals[weight] = mem

        if not candidate_signals:
            predict_signal = None
        else:
            # Se selecciona la mejor señal.
            predict_weight = max(candidate_signals.keys())
            print(predict_weight)
            predict_signal = candidate_signals[predict_weight]
        return predict_signal

Con la anterior se ha programado el cerebro de nuestro percetrón. Veamos ahora como implementarlo, para esto se debe construir una clase Perceptron con los siguientes métodos:

1. El constructor de la clase Perceptrón con una instancia de la clase Neuron. Aquí también se implementa la interfaz gráfica:

    def __init__(self, sqrt_n_receptors: int=5, sensitivity: int=17, focus:
                 float=1) -> object:

        self.sqrt_n_receptors = sqrt_n_receptors
        self.sensitivity = sensitivity
        self.focus = focus

        # Atributos de la clase.
        self.signal = np.array([-1 for _ in range(self.sqrt_n_receptors ** 2)])
        self.neuron = Neuron(self.signal, self.sensitivity, self.focus)

        # Ventana principal con botones.
        main_window = tk.Tk()
        main_window.title('Perceptron')
        main_window.resizable(width=False, height=False)

        # Botones de la ventana principal.
        n_row = range(sqrt_n_receptors)
        self.buttons = [[None]*sqrt_n_receptors for _ in n_row]
        self.buttons = np.array(self.buttons)
        self.values = np.zeros((self.sqrt_n_receptors, self.sqrt_n_receptors))
        self.coordinates = {}

        kwargs = dict(text=' ', bg='gray', relief='flat', width=2, height=2)

        for row, col in product(n_row, n_row):
            self.buttons[row, col] = tk.Button(main_window, **kwargs)
            self.buttons[row, col].grid(row=row, column=col, padx=2, pady=2)
            self.coordinates[self.buttons[row, col]] = [row, col]

        # Se detectan los eventos de cada uno de los botones.
        for button in self.buttons.flat:
            button.bind("", self.button_pressed)

        # Ventana secundaria.
        second_window = tk.Tk()
        second_window.title('')
        second_window.geometry("110x110")
        second_window.resizable(width=False, height=False)

        # Botones de la ventana secundaria.
        kwargs_memorize = dict(text='learn', command=self.memorize, width=455)
        button_memorize = tk.Button(second_window, **kwargs_memorize)
        kwargs_analyze = dict(text='Analyze', command=self.id_signal, width=455)
        button_analyze = tk.Button(second_window, **kwargs_analyze)
        kwargs_reset = dict(text='Reset', command=self.reset, width=455)
        button_reset = tk.Button(second_window, **kwargs_reset)
        kwargs_forget = dict(text='Forget', command=self.forget, width=455)
        button_forget = tk.Button(second_window, **kwargs_forget)

        # La ventana secundaria se construye con un pack.
        button_memorize.pack()
        button_analyze.pack()
        button_reset.pack()
        button_forget.pack()

        main_window.mainloop()
        second_window.mainloop()

2. El siguiente método define que ocurre cuando se presiona algún botón de la ventana principal.

    def button_pressed(self, event: Callable) -> None:
        # Se identifican las coordenadas del evento.
        row, col = self.coordinates[event.widget]
        n = self.sqrt_n_receptors * row + col

        if self.signal[n] == -1:
            self.signal[n] = 1
            self.buttons[row][col]['bg'] = '#5BADFF'
        else:
            self.signal[n] = -1
            self.buttons[row][col]['bg'] = 'gray'

3. Método para reconocer la señal:

    def id_signal(self) -> None:

        signal = self.neuron.axon(self.signal)
        try:
            n = len(signal)
        except TypeError:
            print('I do not know is this')
        else:
            for i in range(n):
                if signal[i] == 1:
                    row = i // self.sqrt_n_receptors
                    col = i % self.sqrt_n_receptors
                    self.buttons[row][col]['bg'] = '#01D826'

4. Método para resetear el estado de los botones en la ventana principal.

    def reset(self) -> None:

        self.signal = -1 * np.ones(len(self.signal))
        size = range(self.sqrt_n_receptors)
        for row, col in product(size, size):
            self.buttons[row][col]['bg'] = 'gray'
                    self.buttons[row][col]['bg'] = '#01D826'

5. Método para aprender nuevas señales.

    def memorize(self) -> None:

        self.neuron.learn(self.signal)
        self.reset()

5. Método para olvidar todas las señales.

    def forget(self) -> None:

        self.neuron.forget()

Finalmente instanciamos la clase del perceptrón.

if __name__ == "__main__":
    Perceptron()

El lector debe notar, cuando tenga funcionando el perceptrón, que cada que ingresa un nuevo patrón de aprendizaje, se está entrenando una nueva neurona. En otras palabras, para cada patrón se tiene un vector de pesos que permite caracterizar cada nueva neurona entrenada. Estos pesos, se almacenan en memoria.mem. Así, que el lector, puede imaginar el funcionamiento del perceptrón como se muestra en la siguiente figura:



Figura 3. Estructura global del perceptrón.


Otra observación importante, es que el perceptron aprende adecuadamente los pesos sinápticos en un tipo finito. Teóricamente, esto es:

Teorema. Se tiene un perceptrón con un conjunto adecuado de pesos sinápticos $w^*$ para el resultado $\hat{y}(x) = y$. Entonces el perceptrón converge en un tiempo finito, sin importar quién sea $w^*$ inicial.
En efecto, si consideramos que $w^*$ es una solución adecuada, entonces $||w^*|| = 1$ (esto si se considera como criterio de desición a la función $sgn$, hacer esta consideración no representa ninguna perdida de generalidad). Ahora si calculamos $|w^*\cdot x|$, entonces tenemos dos posibilidades o que el resultado sea cero, o que existe un $\delta > 0$ tal que $|w^*\cdot x|>0$ para la entrada $x$. Si ahora se considera \begin{equation} \cos \alpha = \frac{w\cdot w^*}{||w||}, \end{equation} entonces de acuerdo a las reglas de aprendizaje del perceptrón se tiene que $\Delta w = \hat{y}x$, y por lo tanto la modificación a los pesos sería $w' = w +\Delta w$. De esto se sigue que: $$w'\cdot w^* = w\cdot w^*+\hat{y}w^*\cdot x = w\cdot w^*+sgn(w^*\cdot x)w^*\cdot x > w\cdot w^* +\delta$$ por otro lado se tiene: $$||w'||^2=w^2+2\hat{y}w\cdot x + x^2 < w^2 + x^2 = w^2+ M$$ dado que $\hat{y}=-sgn(w\cdot x)$. Después de estás modificaciones, entonces es puede concluir que: $$\cos\alpha > \frac{w^*\cdot w + \delta}{\sqrt{w^2+tM}}.$$ De esta última expresión se concluye que el tiempo de convergencia debe ser finito, dado que $cos \alpha \leq 1$. Con algunas modificaciones, se puede considerar como el tiempo máximo a $t_{máx}=\frac{M}{\delta^2}$.

Conclusiones

  1. Lo que se aprendió hoy fue que una neurona artificial tiene gran parecido a la neurona biológica tanto en su estructura como en su funcionalidad. Así como la unión de las neuronas dan origen al cerebro, las neuronas artificiales dan lugar a arreglos de neuronas llamadas redes neurales, que intentan emular el funcionamiento del cerebro humano, tarea que todavía no se consigue. La neurona artificial al igual que la neurona biológica, maneja diferentes tipos de señales, como se mencionó con antelación estas pueden ser continuas o digitales.
  2. El codigo completo del perceptron lo encuentras aquí.

Referencias