2A. Reporte escrito. Experimentos y análisis de estructuras de datos.

Author

José Alberto Villegas Díaz Disciplina.

Lista de cambios:

Sin cambios. Esta entrega se calificó con la nota más alta.

Introducción

En este reporte, se lleva a cabo un análisis exhaustivo de dos algoritmos de álgebra lineal aplicados sobre matrices: la multiplicación de matrices y la eliminación gaussiana.El objetivo principal de este reporte es implementar y analizar el desempeño de ambos algoritmos al trabajar con matrices aleatorias de tamaño \(nxn\) para valores de \(n=100,300, 1000\). Para evaluar de manera rigurosa su eficiencia, se compararán los resultados en términos del número de operaciones requeridas y el tiempo real de ejecución para cada algoritmo. La elección de estos parámetros permitirá observar cómo cada uno de estos algoritmos escala con el tamaño de las matrices, proporcionando así una visión clara de sus ventajas y limitaciones bajo diferentes condiciones.

Algoritmo para multiplicación de matrices aleatorias de tamaño \(n×n\)

En primera instancia crearemos dos funciones, una para la generación de matrices aleatorias de tamaño \(nxn\) y otra implementando un algoritmo convencional para la multiplicación de las mismas.

La definición dada por Grossman & Flores (2012), de producto de matrices establece:

Sea \(A = (a_{ij})\) una matriz de dimensión \(m \times n\), y sea \(B = (b_{ij})\) una matriz de dimensión \(n \times p\). Entonces el producto de \(A\) y \(B\) es una matriz \(C = (c_{ij})\) de dimensión \(m \times p\), en donde

\[ c_{ij} = (\text{renglón } i \text{ de } A) \cdot (\text{columna } j \text{ de } B) \]

Es decir, el elemento \(c_{ij}\) de \(AB\) es el producto punto del renglón \(i\) de \(A\) y la columna \(j\) de \(B\). Si esto se extiende, se obtiene

\[ c_{ij} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{in}b_{nj} \]

Si el número de columnas de \(A\) es igual al número de renglones de \(B\), entonces se dice que \(A\) y \(B\) son compatibles bajo la multiplicación.

Bajo ese esquema teórico se realizara el algoritmo.

Definición y prueba del algoritmo

# Función para matrices aleatorias de tamaño n x n
function generar_matriz_aleatoria(n)
    return rand(n, n)  # Matriz con números aleatorios entre 0 y 1
end

# Función en la que se implementa un algoritmo convencional para el producto de matrices
function multiplicar_matrices(A, B)
    n = size(A, 1)  #A y B son de tamaño n x n
    C = zeros(n, n)  # Matriz resultado en la que se asignarán los valores de cij
    operaciones = 0  # Contador de operaciones

    # Algoritmo de multiplicación de matrices convencional
    for i in 1:n
        for j in 1:n
            suma = 0  # Variable para la suma en la posición (i, j)
            for k in 1:n
                suma += A[i, k] * B[k, j]  # Multiplicación y suma
                operaciones += 1  # Contamos cada multiplicación
            end
            C[i, j] = suma  # Asignamos el valor a la matriz resultado
        end
    end

    return C, operaciones  # Retorna la matriz resultado y el número de operaciones
end

multiplicar_matrices (generic function with 1 method)

Una vez definidas las funciones podemos realizar una validación simple para matrices pequeñas de orden 3x3.

# Validación de uso con matrices 3x3
n = 3  # Tamaño de la matriz
A = generar_matriz_aleatoria(n)
B = generar_matriz_aleatoria(n)

# Multiplicación de matrices
C, operaciones = multiplicar_matrices(A, B)

# Mostrar el resultado

println("Matriz A:")
println(A)
println("Matriz B:")
println(B)
println("Matriz C (Resultado):")
println(C)

# Mostrar el número de operaciones realizadas
println("Número de operaciones realizadas: $operaciones")

Matriz A:
[0.48776183234824144 0.6124433790334284 0.028635546875128703; 0.8370848622572715 0.7473292661917925 0.9877981495667337; 0.9735150429578445 0.19481858305396282 0.057932071864141954]
Matriz B:
[0.3442737221849198 0.8344937410238792 0.5451937807016894; 0.766291104724915 0.08561097205382118 0.6577453667078036; 0.9520486858842626 0.7301190823571463 0.9768559577402947]
Matriz C (Resultado):
[0.664495929835295 0.48037342841918884 0.696729316958151; 1.8012900205112459 1.483731941698228 1.9128623305140917; 0.539597547560282 0.8713681295588075 0.7154866567372042]
Número de operaciones realizadas: 27

Los algoritmos mostraron un correcto funcionamiento por lo que se procederá a evaluar para los distintos valores de \(n\). Nótese que para esta validación se imprimieron las matrices, esto no es viable en los siguientes casos debido al tamaño de las entradas.

Evaluación para \(n=100\)

# Validación de uso con matrices 100x100
n = 100  # Tamaño de la matriz
A = generar_matriz_aleatoria(n)
B = generar_matriz_aleatoria(n)

# Multiplicación de matrices
C, operaciones = multiplicar_matrices(A, B)
# Mostrar el número de operaciones realizadas
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1000000

Se puede observar que el número de operaciones realizadas cuando \(n=100\) es de 1,000,000 es decir \(100^3\)

using BenchmarkTools
@btime multiplicar_matrices(A, B)

  825.900 μs (4 allocations: 78.23 KiB)

([28.132338392688823 23.908113166359033 … 25.024304177684872 25.7096938977009; 26.667022471404753 23.660024037205346 … 24.639590229379817 26.44517290558687; … ; 26.002885424023155 22.99216990611572 … 23.492611418137535 24.916168129409794; 27.243874373771543 25.87820795626854 … 26.581182371428145 26.91368518219324], 1000000)

Evaluación para \(n=300\)

# Validación de uso con matrices 100x100
n = 300  # Tamaño de la matriz
A = generar_matriz_aleatoria(n)
B = generar_matriz_aleatoria(n)

# Multiplicación de matrices
C, operaciones = multiplicar_matrices(A, B)
# Mostrar el número de operaciones realizadas
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 27000000

Se puede observar que el número de operaciones realizadas cuando \(n=300\) es de 27,000,000 es decir \(300^3\)

using BenchmarkTools
@btime multiplicar_matrices(A, B)

  23.998 ms (4 allocations: 703.23 KiB)

([80.53091461158125 77.91382521022068 … 78.23327276120801 77.15153895643374; 71.37048085506102 70.54930640968472 … 67.77446808097446 69.81725780702683; … ; 78.54949084586757 74.26211081012234 … 74.54586378001189 75.42220879028378; 72.79998868719915 72.28643983521658 … 71.78528089018043 70.06249995014716], 27000000)

Podemos observar que el tiempo en el que incrementó fue apróximadamente de un orden más respecto de \(n=100\) y el número de asignaciones permanece constante

Evaluación para \(n=1000\)

# Validación de uso con matrices 1000x1000
n = 1000  # Tamaño de la matriz
A = generar_matriz_aleatoria(n)
B = generar_matriz_aleatoria(n)

# Multiplicación de matrices
C, operaciones = multiplicar_matrices(A, B)
# Mostrar el número de operaciones realizadas
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1000000000

Se puede observar que el número de operaciones realizadas cuando \(n=1000\) es de 1,000,000,000 es decir \(1000^3\)

using BenchmarkTools
@btime multiplicar_matrices(A, B)

  1.331 s (4 allocations: 7.63 MiB)

([252.28375919010392 256.163478843609 … 255.81051340050507 252.5790668154193; 245.61903803038666 244.96190533494607 … 244.73860377052003 243.5577681299076; … ; 249.72986247348513 252.9162720651867 … 252.232903123629 244.19782743234276; 264.62505711595907 257.55400576599106 … 260.0656411956669 253.25510157969072], 1000000000)

Podemos observar que el tiempo en el que incrementó fue apróximadamente de un orden más respecto de \(n=300\) y el número de asignaciones permanece constante

Experimento para matrices dispersas

using SparseArrays

A = sprand(1000, 1000, 0.01)  # Matriz 1000x1000 con 1% de elementos no ceros
B = sprand(1000, 1000, 0.01)
# Multiplicación de matrices
C, operaciones = multiplicar_matrices(A, B)
# Mostrar el número de operaciones realizadas
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1000000000

using BenchmarkTools
@btime multiplicar_matrices(A, B)

  21.248 s (4 allocations: 7.63 MiB)

([0.0 0.5179738303041588 … 0.0 0.0; 0.0 0.0 … 0.0 0.28833088395964407; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], 1000000000)

Se observa que para \(n=1000\) el tiempo de ejecución incrementa en un orden. Las asignaciones son constantes

Algortimo para Eliminación gaussiana / Gauss-Jordan para matrices aleatorias de tamaño \(n×n\)

Definición y prueba del algoritmo

En algebra lineal una de las aplicaciones del método de Gauss Jordan es su uso para encontrar la inversa de una matriz invertible \(nxn\), para lo cual el primer paso es agregar una matriz identidad a la derecha de la matriz inverible en forma de matriz aumentada y resolver el sistema por el método de Gauss Jordan. De acuerdo a Grossman & Flores (2012), el método para realizar Gauss-Jordan con pivoteo consiste en:

Escribir el sistema en la forma de matriz aumentada. De la primer columna con componentes diferentes de cero (denominada columna pivote), seleccione la componente con el valor absoluto. Esta componente se denomina pivote:
Reacomodar los renglones para mover el pivote hasta arriba:
Dividir el primer renglón entre el pivote:
Sumar múltiplos del primer renglón a los otros renglones para hacer cero todas las componentes de la columna pivote:
Tapar el primer renglón y realice los pasos 1 al 4 en la submatriz que resulta:
Continuar de esta manera hasta que la matriz esté en la forma escalonada por renglones.
Resolver el sistema

Bajo este esquema se implementará el algoritmo.

Definición y prueba del algoritmo

using LinearAlgebra 

# Crear una matriz aleatoria invertible de tamaño n x n
function generar_matriz_invertible(n)
    while true
        A = rand(n, n)  # Generamos una matriz aleatoria
        if det(A) != 0   # Verificamos que sea invertible
            return A
        end
    end
end

function gauss_jordan(A)
    n = size(A, 1)  # Tamaño de la matriz (nxn)
    operaciones = 0  # Contador de operaciones

    # Convertimos A en una matriz aumentada con la identidad (para calcular la inversa)
    A = hcat(A, Matrix(I, n, n))  # Agregamos la matriz identidad a la derecha para expresarla en forma de matriz aumentada

    # Eliminación hacia adelante
    for i in 1:n
        # Pivoteo parcial (si es necesario)
        max = argmax(abs.(A[i:n, i])) + (i - 1)  # Encuentra el máximo en la columna
        if max != i
            A[i, :], A[max, :] = A[max, :], A[i, :]  # Intercambia filas
            operaciones += n  # Contamos el intercambio de filas
        end

        # Normalizar el pivote a 1
        pivote = A[i, i]
        A[i, :] /= pivote
        operaciones += n  # División de toda la fila

        # Eliminación de los elementos debajo y arriba del pivote
        for j in 1:n
            if j != i  # Evitar la fila del pivote
                factor = A[j, i]
                A[j, :] -= factor * A[i, :]
                operaciones += n  # Contamos la operación de eliminación
            end
        end
    end

    # La parte derecha de A ahora es la matriz inversa
    return A[:, n+1:end], operaciones
end

gauss_jordan (generic function with 1 method)

Una vez definidas las funciones podemos realizar una validación simple para una matriz pequeña de orden 3x3.

# Prueba con una matriz 3x3
n = 3
A = generar_matriz_invertible(n)

# Aplicamos Gauss-Jordan para encontrar la inversa
A_inv, operaciones = gauss_jordan(A)

# Mostrar resultados
println("Matriz A:")
println(A)
println("Matriz Inversa A^-1:")
println(A_inv)
println("Número de operaciones realizadas: $operaciones")

Matriz A:
[0.7176888293553598 0.1159503346927897 0.7824477172960047; 0.26816290254571196 0.5436422839447754 0.09663661895942965; 0.5776198687754796 0.502034634542894 0.6842905799728082]
Matriz Inversa A^-1:
[4.201278869398719 4.071118358154523 -5.378855536411074; -1.6582303726782912 0.5084526591919578 1.8242884828509491; -2.32978959049701 -3.809521529362638 4.663337360923403]
Número de operaciones realizadas: 27

Se corrobora que el algoritmo funciona correctamente. También desde este momento podemos observar que el número de operaciones realizadas para \(n=3\) es de 27, es decir \(3^3\)

Evaluación para \(n=100\)

# Experimento para n=100
n = 100
A = generar_matriz_invertible(n)
A_inv, operaciones = gauss_jordan(A)
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1009600

Se puede observar que el número de operaciones realizadas cuando \(n=100\) es de apróximadamente 1,000,000 es decir \(100^3\)

using BenchmarkTools
@btime gauss_jordan(A)

  7.225 ms (80394 allocations: 63.79 MiB)

([0.28970373469019167 -0.17420484581562626 … 0.27307472669636884 -0.3588290552142628; 0.11215087930590079 -0.25156162543991023 … 0.14969481586464503 -0.4030287746902638; … ; -0.0754963833616388 -0.10535739685169299 … 0.078313017746731 -0.22571431989707574; -0.04624036046053944 0.34960007812129645 … 0.02804706715744275 -0.30342600533361846], 1009600)

Evaluación para \(n=300\)

# Experimento para n=300
n = 300
A = generar_matriz_invertible(n)
A_inv, operaciones = gauss_jordan(A)
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 27087900

Se puede observar que el número de operaciones realizadas cuando \(n=300\) es de apróximadamente 27,000,000 es decir \(100^3\)

using BenchmarkTools
@btime gauss_jordan(A)

  332.074 ms (1081266 allocations: 1.63 GiB)

([0.2446861370288894 0.2603838908496686 … 0.9280233428073322 0.09375809661532152; 0.1298458317229163 0.12987395132258728 … 0.35016998210347056 0.01506863357009279; … ; 0.13790527314543724 0.35112067320686813 … 0.026894330076987827 -0.13166694676382443; 0.07039322841938805 -0.06504029318774748 … 0.5015589283627557 0.039208192096636965], 27087900)

El tiempo de ejecución aumenta aproximadamente en dos órdenes respecto a \(n=100\) lo mismo sucede para las asignaciones

Evaluación para \(n=1000\)

# Experimento para n=1000
n = 1000
A = generar_matriz_invertible(n)
A_inv, operaciones = gauss_jordan(A)
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1000994000

Se puede observar que el número de operaciones realizadas cuando \(n=1000\) es de apróximadamente 1,000,000,000 es decir \(1000^3\)

using BenchmarkTools
@btime gauss_jordan(A)

  30.021 s (12005472 allocations: 59.87 GiB)

([0.047490859301975295 -0.060854470450899374 … -0.12558169672686376 0.1434304526853882; -0.1200956782253656 -0.428885355426993 … -0.013085833224851684 0.30326859892872826; … ; 1.7067673603290574 1.0086591100969922 … -0.831195709124633 -0.3635885587092932; 1.237672201837348 0.8858751875164227 … -0.7664065225726228 -0.16137824678541648], 1000994000)

El tiempo de ejecución aumenta aproximadamente en dos órdenes respecto a \(n=300\) las asignaciones incrementan en un orden

Experimento para matrices dispersas

Bajor la misma lógica que se evaluó este algoritmo se generará una matriz dispersa tamaño \(nxn\) que sea invertible.

using SparseArrays, LinearAlgebra

function generar_matriz_invertible_dispersada(n::Int, densidad::Float64)
    # Primero generamos una matriz dispersa aleatoria de tamaño n x n
    A = sprandn(n, n, 0.01)

    # Hacemos la matriz simétrica A = A + A^T para asegurar que sea simétrica
    A = A + transpose(A)

    # Añadimos una constante al diagonal para asegurar que sea positiva definida
    # Esto ayuda a evitar que tenga autovalores no positivos
    for i in 1:n
        A[i,i] += n  # Incrementamos la diagonal para garantizar la positividad
    end

    # A es ahora simétrica y positiva definida, y como es dispersa, debería ser invertible
    return A
end

generar_matriz_invertible_dispersada (generic function with 1 method)

# Experimento para n=1000
n = 1000
A = generar_matriz_invertible(n)
A_inv, operaciones = gauss_jordan(A)
println("Número de operaciones realizadas: $operaciones")

Número de operaciones realizadas: 1000987000

Observarmos que el número de operaciones continua siendo aproximadamente \(n^3\)

using BenchmarkTools
@btime gauss_jordan(A)

  31.942 s (12005430 allocations: 59.87 GiB)

([-0.286088270150775 -0.04415745714555952 … -0.02341198529697472 -0.060345876344513556; -0.6395347509959397 0.20227489477122187 … 0.015799245413789784 -0.0076372936972852745; … ; -0.10946769430792393 0.008548864678396426 … -0.007695531472638406 -0.08864286022847348; -0.21194227741102356 0.08071600512590583 … 0.009539117934955094 -0.04134071560159264], 1000987000)

No se notó una diferencia significativa en tiempo de ejecución ni en asignaciones

Discusión y conclusiones

Del experimento podemos concluir que el orden de crecimiento para el algoritmo del método tradicional de multiplicación de matrices como para el algoritmo de eliminación gaussiana es de \(O(n^3)\). Este resultado es previsible ya que está bien documentado. Sin embargo, se debe apuntar que para los algoritmos de multiplicación de matrices el mínimo exponente para multiplicar matrices \(nxn\) oscila entre 2 y 3 (Anderson & Barman, 2009) y al día de hoy existen múltiples algoritmos, de un orden menor al tradicional, siendo uno de los más conocidos el algoritmo de Strassen (1969) de orden \(O(n^{2.81})\), el cual se basa en el método de “Divide y vencerás” (Cormen et al. 2022). Sin embargo estos algoritmos van más allá del alcance de este reporte.

Si bien ambos algoritmos poseen un orden de crecimiento igual, el algoritmo de eliminación gaussiana, plantea una problemática adicional que es el número de asignaciones, estas son muy altas lo que se interpreta como que Julia está creando múltiples copias de matrices, mientras que en el algoritmo de multiplicación de matrices se aprecia un número de asignaciones constante. Esta discrepacia en los algoritmos afecta directamente los tiempos de ejecución los costos en memoria para el algoritmo de eliminación gaussiana.

Finalmente el experimento haciendo uso de matrices dispersas demostró que el tiempo de ejecución fue mayor para una misma entrada \(n=1000\), para el algoritmo de multiplicación de matrices, mientras que para el algoritmo de eliminación gaussiana no demostró diferencias significativas.

¿Cuál es el impacto de acceder los elementos contiguos en memoria de una matriz?

El acceso a los elementos contiguos en memoria de una matriz podría tener un impacto significativo en los tiempos de ejecución de los algoritmos.Cuando accedes a elementos contiguos en memoria (por ejemplo, en un orden fila-por-fila o columna-por-columna), se puede aprovechar la memoria caché, para acceder a estos y que no tengan que ser cargados desde la memoria principal, que es más lenta.

¿Qué cambiarías si utilizas matrices dispersas? ¿Cuáles serían los costos?

Como se pudo corroborar en los experimentos el uso de matrices dispersas se comportó distinto para cada algoritmo, para el de multiplicación incrementó los tiempos de ejecución, mientras que para el de eliminación gaussiana no tuvo diferencias considerables para \(n=1000\)

Referencias

Anderson, M., & Barman, S. (2009, diciembre 6). The Coppersmith-Winograd Matrix Multiplication Algorithm.

Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford (2022). Introduction to Algorithms (2nd ed.). MIT Press.

Grossman, S. I., & Flores, J. J. (2012). Álgebra lineal (7a ed.). McGraw-Hill.