Algoritmo de Babai con Python

 En la entrada anterior expliqué y ejemplifiqué el Algoritmo de Babai para encontrar el vértice del n-paralelepípedo más cercano a un punto.

Podemos realizar un script en Python para dicho algoritmo. Utilizando la librería NumPy queda muy sencillo:

 

import numpy as np

def babai(B, w):
    # Calculamos los coeficientes del vector w en la base B
    v = np.linalg.solve(B, w)
    # Devolvemos el producto matricial de B por las coordenadas de v redondeadas
    return B @ np.round(v).astype(int) 

 

By www.python.org - www.python.org, GPL,

 https://commons.wikimedia.org/w/index.php?curid=34991651

 

PD: Podéis encontrar en uno de mis repositorios de Github un script para probar el Algoritmo de Babai y la razón de Hadamard. 

Resolución de sistemas compatibles determinados con wxMaxima (y III): métodos iterativos

En las entradas anteriores de la serie "Resolución de sistemas compatibles determinados con wxMaxima" expliqué como aplicar en wxMaxima algunos métodos directos y algunos métodos basados en factorización para resolver sistemas de ecuaciones lineales (compatibles determinados).

En esta entrada explico cómo aplicar algunos de los métodos denominados iterativos: el método de Jacobi, el método de Gauss-Seidel y el método de relajación.

Podemos decir, simplificando la cuestión, que la principal diferencia entre los métodos iterativos y los métodos vistos en las entradas anteriores es que los métodos iterativos consisten en empezar con un valor inicial que se le da al vector solución x y se aplica un algoritmo a las sucesivas aproximaciones que se van obteniendo. Es decir, la entrada de una aplicación del algoritmo es la salida de la aplicación anterior del algoritmo. Estos algoritmos pueden acercarse más a la solución del sistema en cada iteración (en ese caso se dice que el algoritmo converge) o alejarse de ella (o entrar en un bucle) por lo que no obtendríamos la solución del sistema.

En el caso de que el algoritmo converja a la solución debemos decidir cuándo consideramos que la solución es lo "suficientemente buena". Ese parámetro lo llamamos tolerancia y, de nuevo simplificando la explicación, podemos decir que es el margen de error que vamos a permitir (o una medida de lo lejos que puede estar la aproximación de la solución).

Recordemos que estamos resolviendo el siguiente sistema:

10x+4y+4z+3t+5s=1
  4x+8y+2z+2t+5s=1
  4x+2y+6z+3t+4s=1
  3x+2y+3z+8t+5s=1
  5x+5y+4z+5t+6s=1

Y que en la primera entrada de la serie comprobamos que se trata de un sistema compatible determinado.

Al aplicar en wxMaxima cualquiera de los 3 métodos iterativos que se explican en esta entrada, se definen una serie de parámetros previos (que son comunes a los 3 métodos):

D:zerofor(A);
for i:1 thru matrix_size(A)[2] do D[i,i]:A[i,i];
D;

x0:zerofor(b);

n:matrix_size(A)[2];
tol: 1e-10;
Kmax: 10000;

El parámetro tol es la tolerancia comentada anteriormente y el parámetro Kmax es el número máximo de iteraciones que vamos a permitir al algoritmo (para evitar que entre en un bucle infinito). En nuestro ejemplo permitimos un máximo de 10000 iteraciones. Vamos entonces con los métodos:

Método de Jacobi:

MJ:D;
NJ:MJ-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MJ).(NJ.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

Al ejecutar la instrucción anterior wxMaxima devuelve un error de desbordamiento, lo que significa que, para el sistema que estamos resolviendo, el método de Jacobi no converge.

Método de Gauss-Seidel:

MGS:copymatrix(A);

for i:1 thru n do
 for j:1 thru n do
   if i

MGS;

NGS:MGS-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MGS).(NGS.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

x;

float(%);

Método de relajación:

En este ejemplo escogemos un parámetro de relajación de 1.5 (sobre-relajación).

w:1.5;

MRA:copymatrix(A);

for i:1 thru n do
 for j:1 thru n do
   if i<=j then MRA[i,j]:0;

MRA;

MR:D/w+MRA;

NR: MR-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MR).(NR.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

x;

Con esta entrada doy por finalizada la serie "Resolución de sistemas compatibles determinados con wxMaxima". Espero que pueda ser de utilidad.

Resolución de sistemas compatibles determinados con wxMaxima (II): métodos de factorización

En la primera entrada de la serie expliqué algunos métodos directos para resolver, con wxMaxima, sistemas de ecuaciones lineales que son compatibles determinados.

Vamos a ver ahora, utilizando el mismo ejemplo, otro tipo de métodos para resolver este tipo de sistemas: los métodos de factorización. Podemos resumir dichos métodos diciendo que el objetivo es, dado un sistema Ax=b, encontrar una factorización de A en dos matrices cuyas inversas sean "más fáciles" de calcular (en realidad "más fácil" significa menos coste en términos computacionales). Si A = B · C, entonces el sistema Ax=b nos queda BCx=b y, por consiguiente, podríamos descomponer el sistema Ax=b en dos sistemas By=b, Cx=y (o aplicar el método de la inversa para despejar el vector de soluciones x = C^(-1) · (B^(-1) · b) ).

A continuación explico cómo ejecutar en wxMaxima 3 de estos métodos de factorización: LU, Cholesky y QR.

Recordad que estamos resolviendo el sistema:

10x+4y+4z+3t+5s=1
  4x+8y+2z+2t+5s=1
  4x+2y+6z+3t+4s=1
  3x+2y+3z+8t+5s=1
  5x+5y+4z+5t+6s=1

Y ya comprobamos en la entrada anterior que es un sistema compatible determinado.

Método de factorización LU:

lu_backsub(lu_factor(A),transpose(b));

Método de factorización de Cholesky:

CH:cholesky(A);

 

Generalmente, si A no cumple las condiciones necesarias para poder aplicar el método de Cholesky (que sea Hermítica y definida positiva), wxMaxima devuelve un mensaje de error. En cualquier caso, no está de más comprobar que la matriz devuelta por wxMaxima es la correcta:

is(CH.transpose(CH)=A);

true

invert(transpose(CH)).(invert(CH).b);

ratsimp(%);

Método de factorización QR:

En este caso es necesario cargar primero el paquete "lapack" para poder utilizar la función del procedimiento de factorización QR.

load(lapack);

[Q,R]:dgeqrf(A);

invert(R).(transpose(Q).b);

 

invert(R).(transpose(Q).b);

En la próxima entrada de esta serie trataré algunos de los métodos de resolución denominados iterativos: método de Jacobi, método de Gauss-Siedel y método de relajación.