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. 

Razón de ortogonalidad de Hadamard con Python

 En una entrada anterior hablé sobre la razón de ortogonalidad de Hadamard:

$$H(B) = \sqrt[n]{\frac{|det(B)|}{\prod_{i=1}^{n}\|\vec{b}_i\|}} $$

Podemos realizar un script en Python para calcular dicha razón, utilizando la librería NumPy:

import numpy as np

def razon_hadamard(B):
    # Calculamos el valor absoluto del determinante de la matriz
    det = np.abs(np.linalg.det(B))
    # Si el determinante es 0 la matriz no es una base y la función devuelve 0
    if not det:
        return 0
    n = len(B)
    # Definimos una variable prod_norm que calcula el producto de la norma de cada vector (columna)
    prod_norm = np.prod(np.linalg.norm(B, axis=0))
    # Devolvemos la razón
    return (det / prod_norm) ** (1/n) 

 

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

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

 

Matrices cuadradas aleatorias con determinante 1 o -1

Objetivo: generar una matriz cuadrada de números enteros aleatorios cuyo determinante sea 1 o -1.

Una primera idea puede ser realizar un bucle en el que generamos una matriz cuadrada de números enteros aleatorios hasta que se satisfaga la condición de que su determinante sea 1 o -1. El problema estaría resuelto, pero el algoritmo no es muy eficiente.

Voy a proponer a continuación otra solución. 

Problema previo: generar una matriz triangular de números enteros aleatorios cuyo determinante sea 1 o -1.

Este problema es muy sencillo teniendo en cuenta que el determinante de una matriz triangular es igual al producto de los elementos de su diagonal.

Únicamente necesitamos rellenar la diagonal con elementos del conjunto {-1, 1}, ceros por encima (o por debajo) de la diagonal y números enteros aleatorios en el resto de elementos de la matriz. 

Ejemplo de función en Python que calcula una matriz triangular de enteros aleatorios (entre -10 y 10, ambos incluidos) utilizando la librería NumPy.

import numpy as np
# Función para crear una matriz triangular de enteros con determinante 1 o -1
# n es la dimensión de la matriz
def matriz_det1(n):
    # Se crea una matriz M_azar con todos sus elementos aleatorios
    M_azar = np.random.randint(-10, 11, size=(n, n))
    # Se elige al azar si se construye una matriz triangular superior o inferior
    if np.random.randint(2):
        # Triangular superior: se cogen los elementos de M_azar por encima de la diagonal
        # con el resto de elementos iguales a 0 y se le suma una matriz diagonal con elementos 1 y -1
        MU = np.triu(M_azar, k=1) + np.diag(np.random.choice([-1, 1], size=n))
    else:
        # Triangular inferior: se cogen los elementos de M_azar por debajo de la diagonal
        # con el resto de elementos iguales a 0 y se le suma una matriz diagonal con elementos 1 y -1
        MU = np.tril(M_azar, k=-1)+ np.diag(np.random.choice([-1, 1], size=n))
    return MU 

Ejemplo de matriz (10x10) creada con dicho algoritmo:

 

Teniendo resuelto el problema previo, ahora podemos abordar el problema original. Las matrices resultantes del algoritmo anterior tienen el "inconveniente" de no ser lo suficientemente aleatorias, dado que son triangulares y tienen demasiados elementos 0 no aleatorios.

Sin embargo, sabemos que: 

- El determinante del producto de matrices es igual al producto de los determinates (|A·B|=|A|·|B|). Por tanto, si multiplicamos matrices cuyo determinante es 1 o -1 el resultado será una matriz cuyo determinante es 1 o -1.

- El producto de matrices triangulares no tiene porqué ser una matriz triangular.

Utilizando esto, podemos idear un algoritmo para el objetivo principal que consista en crear varias matrices triangulares resultantes del algoritmo anterior y multiplicarlas:

# Función para crear una matriz con determinante 1 o -1
# mediante la multiplicación de matrices triangulares
# n es la dimensión de la matriz y m el número de iteraciones
def genU(n, m):
    # Comenzamos con la matriz identidad
    U = np.eye(n, dtype=int)
    # Creamos "m" matrices triangulares con determinante 1 o -1 y las multiplicamos
    for k in range(m):
        U = U @ matriz_det1(n)
    return U

 Ejemplo de matriz creada con dicho algoritmo (n=5, m=10):

 

 

 

Polinomios de Laguerre con wxMaxima y Python

Los polinomios de Laguerre son unos polinomios conocidos principalmente (aunque no únicamente) por ser soluciones de la ecuación diferencial

x · y'' + (1 - x) · y' + n · y = 0 

Su nombre se debe al matemático francés Edmond Nicolas Laguerre, quien publicó muchos artículos, principalmente en las áreas de geometría y análisis.

 Una de las formas de calcular los polinomios de Laguerre viene dada por la siguiente expresión:

 

 También se pueden calcular de manera recursiva:

 ;   ;  

En wxMaxima podemos calcular muy fácilmente los polinomios de Laguerre utilizando la función laguerre().

Ejemplo de cálculo de los primeros 5 polinomios de Laguerre:

 En Python también existen maneras de calcular los polinomios de Laguerre, utilizando funciones de las librerías NumPy, SciPy o SymPy.

A modo de ejemplo, el siguiente código calcula y representa gráficamente los primeros polinomios de Laguerre:

# Cálculo y representación gráfica de los polinomios de Laguerre
# Félix Rodríguez Díaz
# 25 de agosto de 2025

import matplotlib.pyplot as plt
from scipy.special import laguerre as laguerre_sci
from numpy import arange
from numpy.polynomial.laguerre import lag2poly
from numpy.polynomial import Polynomial
from sympy.functions.special.polynomials import laguerre as laguerre_sym
from sympy import symbols, Poly

# Definimos la variable n_polinomios con la cantidad de polinomios que
# queremos calcular (desde grado 0 hasta grado n_polinomios-1)
n_polinomios = 5

# Calculamos los polinomios de Laguerre con SciPy
laguerre_scipy = [laguerre_sci(z) for z in range(n_polinomios)]
# Los mostramos por pantalla
print("Polinomios de Laguerre calculados con SciPy")
for lag_s in laguerre_scipy:
    print(lag_s)

# Calculamos los polinomios de Laguerre con NumPy
# Dado que lag2poly devuelve los coeficientes utilizamos Polynomial
# para crear el correspondiente polinomio
laguerre_numpy = [Polynomial(lag2poly([0] * i + [1])) for i in range(n_polinomios)]
# Los mostramos por pantalla
print("\nPolinomios de Laguerre calculados con NumPy")
for lag_n in laguerre_numpy:
    print(lag_n)
    
# Calculamos los polinomios de Laguerre con SymPy
x_sym = symbols("x")
laguerre_sympy = [laguerre_sym(z, x_sym) for z in range(n_polinomios)]
# Los mostramos por pantalla
print("\nPolinomios de Laguerre calculados con SymPy")
for lag_sym in laguerre_sympy:
    print(lag_sym)
# Si queremos representar gráficamente los polinomios de SymPy con pyplot debemos
# convertir las expresiones simbólicas a expresiones que reconozca numpy
laguerre_sympy_numpy = []
for lag_expr in laguerre_sympy:
    # Obtener coeficientes (del término de mayor grado al menor)
    coeffs = Poly(lag_expr, x_sym).all_coeffs()
    # Convertir a float
    coeffs_float = [float(c) for c in coeffs]
    # Invertir para que NumPy los interprete correctamente (menor a mayor grado)
    coeffs_float.reverse()
    # Crear polinomio de NumPy
    poly_numpy = Polynomial(coeffs_float)
    laguerre_sympy_numpy.append(poly_numpy)

# Representamos gráficamente los polinomios
x = arange(-1, 3.5, 0.01)
fig, ax = plt.subplots()
ax.set_title('Polinomios de Laguerre $L_n$')
for n in range(n_polinomios):
    # Utilizamos por ejemplo los calculados con SciPy
    ax.plot(x, laguerre_scipy[n](x), label=f"$L_{n}$")
    # Si quisiéramos utilizar los calculados con NumPy:
    #ax.plot(x, laguerre_numpy[n](x), label=f"$L_{n}$")
    # Si quisiéramos utilizar los calculados con SymPy:
    #ax.plot(x, laguerre_sympy_numpy[n](x), label=f"$L_{n}$")
plt.legend(loc='best')
plt.grid(True)
plt.show()

 Este es el resultado de ejecutarlo:

Los polinomios asociados/generalizados de Laguerre se utilizan en campos como por ejemplo la mecánica cuántica.

Un problema de datos numéricos con Python y numpy para pensar (solución)

  Esta entrada es la continuación de Un problema de datos numéricos con Python y numpy para pensar (enunciado). Si todavía no lo has leído y quieres pensar por ti mismo el problema estás a tiempo :-)

 Resumiendo la entrada anterior: tenemos unos datos en csv, los importamos y los mostramos por pantalla. Observamos en pantalla varios datos que tienen valor 9. pero al contar la cantidad de datos que son mayores o iguales que 9 el programa devuelve 0.

 Quedó la pregunta abierta para que quien quisiera pensara el posible motivo y en esta entrada explico lo que estaba sucediendo.

Imagen de centrotadi.com

 La importación desde el csv es correcta, así como también la definición del array/matriz. El problema es que al mostrarlo por pantalla los datos se desvirtúan y pierden precisión.

El csv en realidad contiene la siguiente información:

1.999999999999990,7.999999999999990,7.999999999999990,1.999999999999990,8.999999999999989,3.999999999999990,0.999999999999990,1.999999999999990,1.999999999999990,0.999999999999990
6.999999999999990,6.999999999999990,1.999999999999990,7.999999999999990,5.999999999999990,5.999999999999990,3.999999999999990,2.999999999999990,7.999999999999990,8.999999999999989
4.999999999999990,3.999999999999990,5.999999999999990,1.999999999999990,0.999999999999990,5.999999999999990,1.999999999999990,5.999999999999990,5.999999999999990,2.999999999999990
3.999999999999990,6.999999999999990,5.999999999999990,0.999999999999990,0.999999999999990,8.999999999999989,8.999999999999989,4.999999999999990,5.999999999999990,5.999999999999990
1.999999999999990,3.999999999999990,5.999999999999990,7.999999999999990,0.999999999999990,6.999999999999990,5.999999999999990,3.999999999999990,7.999999999999990,6.999999999999990
7.999999999999990,0.999999999999990,5.999999999999990,2.999999999999990,4.999999999999990,6.999999999999990,8.999999999999989,4.999999999999990,3.999999999999990,0.999999999999990
6.999999999999990,2.999999999999990,6.999999999999990,0.999999999999990,0.999999999999990,3.999999999999990,4.999999999999990,8.999999999999989,2.999999999999990,6.999999999999990
5.999999999999990,3.999999999999990,5.999999999999990,0.999999999999990,6.999999999999990,1.999999999999990,8.999999999999989,4.999999999999990,2.999999999999990,5.999999999999990
6.999999999999990,7.999999999999990,2.999999999999990,1.999999999999990,7.999999999999990,2.999999999999990,0.999999999999990,4.999999999999990,1.999999999999990,1.999999999999990
2.999999999999990,3.999999999999990,7.999999999999990,2.999999999999990,6.999999999999990,6.999999999999990,1.999999999999990,0.999999999999990,5.999999999999990,8.999999999999989

 Y podemos observar que en realidad no hay ningún dato mayor o igual que 9, por lo que el programa no estaba fallando. Lo que fallaba era nuestra visualización de los datos.

Esto me sucedió realizando otro programa más complejo y me costó dar con la clave. Para centrarme en el problema real he diseñado este otro programa mucho más sencillo, que se puede utilizar a modo de problema.

 ¡Ojo! Otro fallo que se puede cometer es al ver por pantalla que los datos se muestran como float pero vemos que son enteros trabajar con ellos pasándolos a enteros con int. El resultado de hacer eso cambiaría completamente los datos, dado que al hacer int truncamos el número y por ejemplo 1.999999999999990 pasa a ser 1.

 En este caso, si queremos hacer una conversión a entero, lo correcto es utilizar redondeo (round). Podríamos modificar el código de la siguiente manera:

import numpy as np
from random import choice

# Importamos los datos del archivo CSV
M = np.genfromtxt("datos-ir.csv", delimiter = ",")

print("La matriz importada es: ")
print("M =\n" , M)

print("\nCreamos una matriz donde cada elemento es True si el correspondiente")
print("elemento de la matriz original es mayor o igual que 9 y False en otro caso")

N = (np.round(M,0) >= 9)

print("N =\n", N)
print("\nLa cantidad de True en la nueva matriz es: ", N.sum())

Un problema de datos numéricos con Python y numpy para pensar (enunciado)

Tenemos unos datos numéricos en un archivo CSV. Importamos dichos datos y los mostramos por pantalla:

[[2. 8. 8. 2. 9. 4. 1. 2. 2. 1.]
[7. 7. 2. 8. 6. 6. 4. 3. 8. 9.]
[5. 4. 6. 2. 1. 6. 2. 6. 6. 3.]
[4. 7. 6. 1. 1. 9. 9. 5. 6. 6.]
[2. 4. 6. 8. 1. 7. 6. 4. 8. 7.]
[8. 1. 6. 3. 5. 7. 9. 5. 4. 1.]
[7. 3. 7. 1. 1. 4. 5. 9. 3. 7.]
[6. 4. 6. 1. 7. 2. 9. 5. 3. 6.]
[7. 8. 3. 2. 8. 3. 1. 5. 2. 2.]
[3. 4. 8. 3. 7. 7. 2. 1. 6. 9.]]
 

Vemos que se trata de una matriz de 10x10 con números entre 1 y 9 (incluidos).

Creamos otra matriz de manera que un elemento es True si el correspondiente elemento de la matriz original es mayor o igual que 9 y False en caso contrario. Finalmente contamos la cantidad de True en la nueva matriz.

Es decir, en realidad con el procedimiento anterior estamos contando cuántos elementos de la matriz original son iguales o mayores que 9.

En la pantalla vemos que hay varios 9 en los datos originales, sin embargo la matriz creada no tiene ningún True (ver más abajo la salida del programa). Es decir, el procedimiento concluye que en la matriz original no hay ningún 9. ¿Cómo es esto posible? ¿Dónde está el fallo?

La respuesta la subiré próximamente en otra entrada. Mientras tanto podéis poner vuestras hipótesis en los comentarios.

 

El código del programa es el siguiente

import numpy as np
# Importamos los datos del archivo CSV
M = np.genfromtxt("datos-ir.csv", delimiter = ",")
print("La matriz importada es: ")
print("M =\n" , M)
print("\nCreamos una matriz donde cada elemento es True si el correspondiente")
print("elemento de la matriz original es mayor o igual que 9 y False en otro caso")
N = (M >= 9)
print("N =\n", N)
print("\nLa cantidad de True en la nueva matriz es: ", N.sum())

 

La salida del programa es la siguiente:

La matriz importada es:  
M =
[[2. 8. 8. 2. 9. 4. 1. 2. 2. 1.]
[7. 7. 2. 8. 6. 6. 4. 3. 8. 9.]
[5. 4. 6. 2. 1. 6. 2. 6. 6. 3.]
[4. 7. 6. 1. 1. 9. 9. 5. 6. 6.]
[2. 4. 6. 8. 1. 7. 6. 4. 8. 7.]
[8. 1. 6. 3. 5. 7. 9. 5. 4. 1.]
[7. 3. 7. 1. 1. 4. 5. 9. 3. 7.]
[6. 4. 6. 1. 7. 2. 9. 5. 3. 6.]
[7. 8. 3. 2. 8. 3. 1. 5. 2. 2.]
[3. 4. 8. 3. 7. 7. 2. 1. 6. 9.]]

Creamos una matriz donde cada elemento es True si el correspondiente
elemento de la matriz original es mayor o igual que 9 y False en otro caso
N =
[[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]
[False False False False False False False False False False]]

La cantidad de True en la nueva matriz es:  0