Operaciones sobre matrices

Comparación de algoritmos de álgebra lineal en Julia, evaluando rendimiento en multiplicación, reducción y cálculo de inversa de matrices.

using BenchmarkPlots, BenchmarkTools, DataFrames, Random, StatsBase, StatsPlots
rng = Xoshiro(20250217)

Xoshiro(0x36f7c6e102ef02a0, 0x23bced4a6d0428c2, 0xb526ca0691c1c909, 0x25316b9df93bb889, 0x304a7232f36675d2)

1 Introducción

En este documento se implementan en Julia algoritmos de multiplicación de matrices, eliminación gaussiana y eliminación de Gauss-Jordan, todos ellos descritos en Cormen et al. (2022). A continuación, se comparan sus desempeños con respecto a matrices aleatorias de tamaño \(n \times n\) para \(n=100,\, 300,\, 1000\), mediante el uso de BenchmarkTools.

A partir del análisis de los algoritmos y de los resultados obtenidos en los experimentos, se concluye que el algoritmo de multiplicación de matrices implementado satisface la expectativa del tiempo de ejecución teórica, como se explica en la Sección 2. En esta misma sección, se señalan posibles alternativas para mejorar el desempeño de la función para matrices dispersas, al aprovechar la cantidad de ceros disponibles.¹ Asimismo, se discute el impacto del acceso a bloques de memoria continuos y discontinuos que habría tenido una implementación alternativa del algoritmo de multiplicación de matrices, disponible en Téllez Avila (2025, sec. 3.4).

A continuación, la Sección 3 implementa un algoritmo para hacer la reducción Gaussiana de una matriz. Los resultados obtenidos indican que para valores bajos de \(n\) se obtuvieron los rendimientos teóricamente esperados de implementar el algoritmo, pero esto no sucedió para \(n = 1000\), en donde se obtuvo un tiempo de ejecución muy deteriorado. Además, se observa que hay efectos diferenciados en términos del consumo de tiempo y memoria. Se concluye que es necesaria una mayor investigación para determinar la causa de este comportamiento.

Finalmente, la Sección 4 discute la relación entre invertir una matriz y realizar la reducción de Gauss-Jordan e implementa un algoritmo para encontrar la inversa de una matriz a partir de la ejecución iterada de un procedimiento que utiliza el resultado del algoritmo discutido en la Sección 3. En este caso, se encontraron rendimientos más deteriorados que en dicha sección, por lo que también se requiere una mayor investigación para determinar la causa del deterioro.

2 Multiplicación de matrices

Cormen et al. (2022, 374) consideran el Algoritmo 1 que calcula \(C = C + A B\) para tres matrices \(A = (a_{ij}),\, B = (b_{ij}),\, C = (c_{ij})\) con dimensiones \(p \times q,\, q \times r, \, p \times r\), respectivamente. Este algoritmo sirve para el cálculo de \(A B\), para lo cual únicamente se requiere inicializar \(C\) como una matriz de ceros. Este código se implementó en Julia mediante la función matmul!().

\begin{algorithm} \caption{MULTIPLICAR-MATRIZ-RECTANGULAR($A, B, C, p, q, r$)}\begin{algorithmic} \For{$i=1$ a $p$} \For{$j = 1$ a $r$} \For{$k = 1$ a $q$} \State $c_{ij} = c_{ij} + a_{ik} \cdot b_{kj}$ \EndFor \EndFor \EndFor \end{algorithmic} \end{algorithm}

function matmul!(A::Matrix, B::Matrix, C::Union{Nothing, Matrix} = nothing)
    p, q = size(A)
    r = size(B)[2]
    @assert q == size(B)[1]
    if isnothing(C)
        C = zeros(Float64, (p, r) )
    else
        @assert size(C) == (p, r)
        C = copy(C)
        C = convert(Matrix{Float64}, C)
    end
    for i in 1:p
        for j in 1:r
            for k in 1:q
                C[i, j] = C[i, j] + A[i, k] * B[k, j]
            end
        end
    end
    return C
end

matmul! (generic function with 2 methods)

Al respecto, se realizaron experimentos para ejecutar la función matmul!() en multiplicaciones de matrices cuadradas aleatorias de \(n \times n, \, n = 100, 300, 1000\). En la Tabla 1 se muestra la mediana del tiempo de ejecución obtenido para cada uno de estos experimentos y en la Figura 1 se grafican las distribuciones de dicho tiempo.

segundos = 5*60

experimentos = [100, 300, 1000]
lab_matmul = BenchmarkGroup()
for i in experimentos
    lab_matmul[i] =  @benchmarkable matmul!(A, B) setup = (A = $(rand(rng, i, i)); B = $(rand(rng, i, i)))
end
tune!(lab_matmul)
lab_matmul = run(lab_matmul, seconds = segundos)
DataFrame(median(lab_matmul), ["n", "Tiempo de ejecución"])

Tabla 1: Tiempo mediano de ejecución de matmul!(), para \(n = 100, 300, 1000\)

3×2 DataFrame

Row	n	Tiempo de ejecución
	Int64	TrialEst…
1	300	TrialEstimate(78.425 ms)
2	1000	TrialEstimate(3.102 s)
3	100	TrialEstimate(2.363 ms)

plot(lab_matmul, yaxis=:log10, xlabel = "n")

Figura 1: Distribución del tiempo de ejecución de matmul!(), para \(n = 100,\, 300,\, 1000\)

Cormen et al. (2022, 374) señalan que el Algoritmo 1 está dominado por el término \(npr\). Por lo tanto, el costo de ejecución esperado depende de \(npr = n^3\) para matrices cuadradas; es decir:

• La ejecución del algoritmo para \(n = 300\) debería ser unas 27 veces más demorada que para el caso de \(n = 100\) (la razón entre el número de operaciones de ambos está dada por \(\frac{300^3}{100^3} = 27\)).
• La ejecución del algoritmo para \(n = 1000\) debería ser unas 37 veces más demorada que para el caso de \(n = 300\).

A continuación se compara el rendimiento obtenido para estos casos, a partir de lo que se observa que las diferencias teóricamente esperadas se encuentran en un rango cercano al de las realmente obtenidas.

for i in eachindex(experimentos)
    if i == firstindex(experimentos)
        continue
    end

    println("Comparación de n = $(experimentos[i]) y n = $(experimentos[i - 1])")
    display(
        judge( median( lab_matmul[experimentos[i]] ),
               median( lab_matmul[experimentos[i - 1]] ))
        )
end

Comparación de n = 300 y n = 100

BenchmarkTools.TrialJudgement: 
  time:   +3219.42% => regression (5.00% tolerance)
  memory: +799.12% => regression (1.00% tolerance)

Comparación de n = 1000 y n = 300

BenchmarkTools.TrialJudgement: 
  time:   +3855.18% => regression (5.00% tolerance)
  memory: +1010.99% => regression (1.00% tolerance)

2.1 Multiplicación de matrices dispersas

De una revisión del Algoritmo 1, no pasa desapercibido que los cambios en \(c_{ij}, \, i = 1, \dots, p,\, j = 1, \dots, r\) ocurren únicamente cuando \(a_{ik} \neq 0\) y \(b_{kj} \neq 0\), como se observa en la línea 4 de dicho pseudocódigo. Se observa que esta misma línea es la que domina el tiempo de ejecución, pues se encuentra en el interior de los tres iteradores de la función; por lo tanto, evitar su ejecución en los casos en que sea posible ayuda a mejorar el rendimiento general del algoritmo.

En particular, para el caso en que \(A\) o \(B\) son matrices dispersas (es decir, que tienen una proporción alta de ceros) es posible realizar una modificación al algoritmo para que verifique si \(a_{ik} \neq 0\) y \(b_{kj} \neq 0\) y, únicamente si ambas condiciones se satisfacen, realice la multiplicación correspondiente. Esta modificación se muestra en el Algoritmo 2, en donde:

• Se intercambió el orden del segundo y tercer iterador, de forma que se puede verificar si \(a_{ik} \neq 0\) antes de iterar sobre la k-ésima fila de \(B\) (línea 3). De esta forma, en la iteración de la línea 4 no se realiza \(n^2\) veces, sino únicamente la cantidad de valores diferentes a cero de \(A\) (denotada \(n_1^A\)).
• Se adicionó la verificación de que el elemento \(b_{kj} \neq 0\) (línea 5) antes de realizar la operación costosa de multiplicado de la línea 6.

Con base en estas modificaciones, la iteración de la línea 4 únicamente se realiza \(n_1^A \cdot n\) veces, en lugar de \(n^3\) como con el algoritmo original.

\begin{algorithm} \caption{MULTIPLICAR-MATRIZ-RECTANGULAR-DISPERSA($A, B, C, p, q, r$)}\begin{algorithmic} \For{$i=1$ a $p$} \For{$k = 1$ a $q$} \If{$a_{ik} \neq 0$} \For{$j = 1$ a $r$} \If{$b_{kj} \neq 0$} \State $c_{ij} = c_{ij} + a_{ik} \cdot b_{kj}$ \EndIf \EndFor \EndIf \EndFor \EndFor \end{algorithmic} \end{algorithm}

Este código se implementó en Julia mediante la adición de un nuevo método a la función matmul!().

function matmul!(A::Matrix, B::Matrix, dispersa::Bool, C::Union{Nothing, Matrix} = nothing)
    p, q = size(A)
    r = size(B)[2]
    @assert q == size(B)[1]
    if isnothing(C)
        C = zeros(Float64, (p, r) )
    else
        @assert size(C) == (p, r)
        C = copy(C)
        C = convert(Matrix{Float64}, C)
    end
    if dispersa
        for i in 1:p
            for k in 1:q
                if A[i, k] != 0
                    for j in 1:r
                        if B[k, j] != 0
                            C[i, j] = C[i, j] + A[i, k] * B[k, j]
                        end
                    end
                end
            end
        end
    end
    return C
end

matmul! (generic function with 4 methods)

Al respecto, se realizaron experimentos para ejecutar la función matmul!() en multiplicaciones de matrices cuadradas aleatorias dispersas, que únicamente tienen el 10% de los valores diferentes de 0, de \(n \times n, \, n = 100, 300, 1000\), tanto con el método señalado en el Algoritmo 1 como con el señalado en el Algoritmo 2. En la Tabla 2 se muestra la mediana del tiempo de ejecución obtenido para cada uno de estos experimentos.

# 1. Función para generar matrices aleatorias dispersas
function disp(n, p)
    A = zeros(n, n)
    listado = sample([(num1, num2) for num1 in 1:n for num2 in 1:n], 
                     round(Int64, p * n^2))
    valores = rand(rng, length(listado))
    for i in eachindex(listado)
        A[listado[i][1], listado[i][2]] = valores[i]
    end
    return A
end

# 2. Ejecución de experimentos
lab_matmul_d = BenchmarkGroup()
for i in experimentos
    lab_matmul_d[1][i] = @benchmarkable matmul!(A, B) setup = (A = disp($i, 0.1); B = disp($i, 0.1))
    lab_matmul_d[2][i] = @benchmarkable matmul!(A, B, true) setup = (A = disp($i, 0.1); B = disp($i, 0.1))
end
tune!(lab_matmul_d)
lab_matmul_d = run(lab_matmul_d, seconds = segundos)

# 3. Mostrar los resultados
innerjoin(DataFrame(median(lab_matmul_d)[1], ["n", "Tiempo (Original)"]),
          DataFrame(median(lab_matmul_d)[2], ["n", "Tiempo (Dispersa)"]),
          on = :n)

Tabla 2: Tiempo mediano de ejecución de matmul!() con matrices dispersas, para \(n = 100,, 300,, 1000\)

3×3 DataFrame

Row	n	Tiempo (Original)	Tiempo (Dispersa)
	Int64	TrialEst…	TrialEst…
1	300	TrialEstimate(79.393 ms)	TrialEstimate(7.357 ms)
2	1000	TrialEstimate(3.129 s)	TrialEstimate(371.363 ms)
3	100	TrialEstimate(2.356 ms)	TrialEstimate(191.833 μs)

Como se señaló previamente, la ejecución de la iteración de la línea 4 del Algoritmo 2 se realiza \(n_1^A \cdot n\) veces; en lugar de las \(n^3\) veces que se ejecuta la de la línea 3 del Algoritmo 1. Por lo tanto, dado que en los experimentos realizados la matriz \(A\) tiene el 10% de sus elementos diferentes a cero (\(n_1^A = 0.1 n^2\)) se espera que los métodos implementados en el algoritmo original tengan una duración de ejecución 10 veces superior a la de los implementados en el modificado. A continuación se compara el rendimiento obtenido para estos casos, a partir de lo que se observa que no se rechazan las hipótesis nulas de que se obtuvieron estas diferencias para \(n = 100,\, 300\); pero sí se obtuvieron unas diferencias un poco menos grandes para \(n = 1000\).

for i in eachindex(experimentos)
    println("Comparación de n = $(experimentos[i])")
    display(
        judge( median(lab_matmul_d[1][experimentos[i]]),
               median(lab_matmul_d[2][experimentos[i]]) )
        )
end

Comparación de n = 100

BenchmarkTools.TrialJudgement: 
  time:   +1128.24% => regression (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

Comparación de n = 300

BenchmarkTools.TrialJudgement: 
  time:   +979.09% => regression (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

Comparación de n = 1000

BenchmarkTools.TrialJudgement: 
  time:   +742.57% => regression (5.00% tolerance)
  memory: +0.00% => invariant (1.00% tolerance)

2.2 Acceso a bloques de memoria

Tanto el Algoritmo 1 como el Algoritmo 2 acceden a los elementos de las matrices \(A, \, B,\, C\) de forma individual, por lo que no hay un impacto significativo de la forma en que están representadas dichas matrices en la memoria. Sin embargo, Téllez Avila (2025, sec. 3.4) propone un algoritmo que calcula \(c_{ij} = A_i B^j\), donde \(A_i\) es la i-ésima fila de \(A\) y \(B^j\) es la j-ésima columna de \(B\). Por lo tanto, la velocidad del acceso a los elementos depende del orden en que estos se encuentran representados: si la representación es por columnas (column major), como en el caso de Julia o Fortran, entonces es más rápida la operación de extraer el vector \(B^j\) que aquella del vector \(A_i\); si la representación es por filas (row major), como en el caso de numpy en Python y C, entonces sucede lo contrario. Al respecto, Téllez Avila (2025) señala que

La representación precisa en memoria es significativa en el desempeño de operaciones matriciales como pueden ser el producto entre matrices o la inversión de las mismas. La manera como se acceden los datos es crucial en el diseño de los algoritmos.

3 Reducción gaussiana

Cormen et al. (2022, 824-32) explican el procedimiento para obtener las descomposiciones \(LU\) y \(LUP\) de una matriz \(A\) a través del método de eliminación gaussiana. A través de la descomposición \(LUP\) de \(A\) se obtienen una matriz triangular inferior \(L\), una matriz triangular superior \(U\) y una matriz de permutación \(P\) tales que \(PA = LU\). Para ello, Cormen et al. (2022) plantean el Algoritmo 3, que genera las matrices señaladas. Este código se implementó en Julia mediante la función descomposicionlup().

\begin{algorithm} \caption{DESCOMPOSICIÓN-LUP($A, n$)}\begin{algorithmic} \State Sea $\pi[1:n]$ un nuevo arreglo \For{$i = 1$ a $n$} \State $\pi[i] = i$ \Comment{Inicializar $\pi$ como la permutación identidad} \EndFor \For{$k = 1$ a $n$} \State $p = 0$ \For{$i = k$ a $n$} \Comment{Encontrar el mayor valor absoluto en la columna $k$} \If{$|a_{ik}| > p$} \State $p = |a_{ik}|$ \State $k' = i$ \Comment{Número de fila del mayor encontrado hasta ahora} \EndIf \EndFor \If{$p = 0$} \State Error: 'Matriz singular' \EndIf \State Intercambiar $a_{ki}$ con $a_{k'i}$ \For{$i = 1$ a $n$} \State Intercambiar $a_{ki}$ con $a_{k'i}$ \Comment{Intercambiar las filas $k$ y $k'$} \EndFor \For{$i = k + 1$ a $n$} \State $a_{ik} = a_{ik} / a_{kk}$ \For{$j = k + 1$ a $n$} \State $a_{ij} = a_{ij} - a_{ik}a_{kj}$ \Comment{Calcular $L$ y $U$ en la propia $A$} \EndFor \EndFor \EndFor \end{algorithmic} \end{algorithm}

function descomposicionlup(A)
    A = copy(A)
    A = convert(Matrix{Float64}, A)

    # Implementa el código igual que en Cormen
    n = size(A)[1]
    @assert n == size(A)[2]
    permutacion = [x for x in 1:n]
    for k in 1:n
        p = 0
        k1 = nothing
        for i in k:n
            if abs(A[i, k]) > p
                p = abs(A[i, k])
                k1 = i
            end
        end
        @assert p != 0
        if k1 != nothing
            permutacion[k], permutacion[k1] = permutacion[k1], permutacion[k]
            for i in 1:n
                A[k, i], A[k1, i] = A[k1, i], A[k, i]
            end
        end
        for i in (k + 1):n
            A[i, k] = A[i, k] / A[k, k]
            for j in (k + 1):n
                A[i, j] = A[i, j] - A[i, k] * A[k, j]
            end
        end
    end

    # Entrega el resultado
    L = zeros(Float64, (n, n))
    U = zeros(Float64, (n, n))
    for i in 1:n
        L[i, i] = 1
        for j in 1:n
            if i > j
                L[i, j] = A[i, j]
            else
                U[i, j] = A[i, j]
            end
        end
    end
    return L, U, permutacion
end

descomposicionlup (generic function with 1 method)

Al respecto, se realizaron experimentos para ejecutar la función descomposicionlup() en matrices cuadradas aleatorias de \(n \times n, \, n = 100, 300, 1000\). En la Tabla 3 se muestra la mediana del tiempo de ejecución obtenido para cada uno de estos experimentos y en la Figura 2 se grafican las distribuciones de dicho tiempo.

lab_lup = BenchmarkGroup()
for i in experimentos
    lab_lup[i] =  @benchmarkable descomposicionlup(A) setup = (A = rand(rng, $i, $i) )
end
tune!(lab_lup)
lab_lup = run(lab_lup, seconds = segundos)
DataFrame(median(lab_lup), ["n", "Tiempo de ejecución"])

Tabla 3: Tiempo mediano de ejecución de descomposicionlup(), para \(n = 100, 300, 1000\)

3×2 DataFrame

Row	n	Tiempo de ejecución
	Int64	TrialEst…
1	300	TrialEstimate(11.390 ms)
2	1000	TrialEstimate(687.338 ms)
3	100	TrialEstimate(355.298 μs)

plot(lab_lup, yaxis=:log10, xlabel = "n")

Figura 2: Distribución del tiempo de ejecución de descomposicionlup(), para \(n = 100,\, 300,\, 1000\)

Cormen et al. (2022, 832) señalan que el Algoritmo 3 tiene un tiempo de ejecución de \(\Theta(n^3)\). Por lo tanto, teóricamente se esperaría una comparación de rendimientos entre experimentos similar a la que se señaló en la Sección 2. No obstante, aunque la comparación de rendimientos de los experimentos con \(n=100\) y \(n=300\) sí se encuentra cercana a esta expectativa, la comparación con el de \(n=1000\) difiere en casi el 70%, como se muestra a continuación.

for i in eachindex(experimentos)
    if i == firstindex(experimentos)
        continue
    end

    println("Comparación de n = $(experimentos[i]) y n = $(experimentos[i - 1])")
    display(
        judge( mean( lab_lup[experimentos[i]] ),
               mean( lab_lup[experimentos[i - 1]] ))
        )
end

Comparación de n = 300 y n = 100

BenchmarkTools.TrialJudgement: 
  time:   +2590.42% => regression (5.00% tolerance)
  memory: +796.59% => regression (1.00% tolerance)

Comparación de n = 1000 y n = 300

BenchmarkTools.TrialJudgement: 
  time:   +5761.58% => regression (5.00% tolerance)
  memory: +1010.06% => regression (1.00% tolerance)

Este comportamiento diferente al esperado se podría explicar porque hay una sustitución no lineal entre la memoria y el tiempo de ejecución para cada uno de los algoritmos. Así, aunque se esperaría que matmul!() y descomposicionlup() tuvieran un comportamiento similar para cada \(n\), en realidad mejora el rendimiento con respecto al tiempo pero empeora con respecto al uso de memoria, al comparar la segunda con la primera función. En cada caso, se rechaza la hipótesis nula de que las diferencias observadas sean iguales a cero.

for i in experimentos
    println("Comparación de matmul!() y descomposicionlup() para n = $i")
    display( judge( mean( lab_lup[i] ), mean( lab_matmul[i] )) )
end

Comparación de matmul!() y descomposicionlup() para n = 100

BenchmarkTools.TrialJudgement: 
  time:   -84.67% => improvement (5.00% tolerance)
  memory: +201.20% => regression (1.00% tolerance)

Comparación de matmul!() y descomposicionlup() para n = 300

BenchmarkTools.TrialJudgement: 
  time:   -84.39% => improvement (5.00% tolerance)
  memory: +200.35% => regression (1.00% tolerance)

Comparación de matmul!() y descomposicionlup() para n = 1000

BenchmarkTools.TrialJudgement: 
  time:   -77.77% => improvement (5.00% tolerance)
  memory: +200.10% => regression (1.00% tolerance)

Este comportamiento de los rendimientos requiere una mayor investigación, pues se observa que no basta con contabilizar el número de operaciones para estimar el rendimiento esperado de los algoritmos de una forma unívoca. Lo anterior, ya que en el sistema en que estos se implementan puede haber sustitución entre tiempo de ejecución y uso de memoria de una forma que no es lineal. Luego, es necesario profundizar en el estudio del uso de diversos recursos escasos para la ejecución de los algoritmos y su relación con \(n\).

4 Reducción de Gauss-Jordan

Grossman (2012, 14-16) señala la diferencia entre el método de eliminación Gaussiana y el método de eliminación de Gauss-Jordan. El primero reduce la matriz \(A\) a su forma escalonada por renglones, que es la que se halló mediante el procedimiento descrito en la Sección 3. Por su parte, el segundo reduce la matriz \(A\) a su forma escalonada reducida por renglones, donde se adiciona la restricción de que cada elemento de la columna en donde hay un pivote tiene todos los otros elementos iguales a cero. Dado que la matriz \(A\) que se está analizando es cuadrada y no singular (véase las líneas 13 a 15 del Algoritmo 3), este procedimiento se reduce a convertir \(A\) en \(I_n\); es decir, en encontrar la inversa de la matriz \(A\).

Cormen et al. (2022, 833-34) señalan que se puede pensar en la resolución del problema para hallar la inversa de \(A\) (la matrix \(X\) que resuelve \(AX = I_n\)) como un problema de \(n\) ecuaciones de la forma \(AX_i = e_i\), en donde \(X_i\) es una columna de la inversa de \(A\) y \(e_i\) es la i-ésima columna de \(I_n\).

Para resolver cada una de estas ecuaciones, se puede utilizar el Algoritmo 4 que utiliza la descomposición LUP de \(A\) y sustituciones hacia adelante y hacia atrás para resolver el sistema de ecuaciones \(Ax = b\) (Cormen et al. 2022, 823-24).

\begin{algorithm} \caption{SOLUCIÓN-LUP($L, U, \pi, b, n$)}\begin{algorithmic} \State Sean $x$ y $y$ nuevos vectores de longitud $n$ \For{$i = 1$ a $n$} \State $y_i = b_{\pi[i]} - \sum_{j = 1}^{i-1} l_{ij} y_j$ \EndFor \For{$i = n$ decreciendo a 1} \State $x_i = \left( y_i - \sum_{j = i + 1}^n u_{ij}x_j \right) / u_{ii}$ \EndFor \State Entregar $x$ \end{algorithmic} \end{algorithm}

Al respecto, se implementó la función solucion_lup() en Julia, que resuelve un sistema \(Ax = b\) a partir de la descomposición LUP de \(A\) (implementada mediante descomposicionlup()) y el Algoritmo 4 de Cormen et al. (2022, 824).

function solucion_lup(A, b)
    n = size(A)[1]
    L, U, permutacion = descomposicionlup(A)
    x = Vector{Float64}(undef, n)
    y = Vector{Float64}(undef, n)
    for i in 1:n
        suma = 0
        for j in 1:(i-1)
            suma += L[i, j] * y[j]
        end
        y[i] = b[permutacion[i]] - suma
    end
    for i in n:-1:1
        suma = 0
        for j in (i+1):n
            suma += U[i, j] * x[j]
        end
        x[i] = (y[i] - suma) / U[i, i]
    end
    return x
end

solucion_lup (generic function with 1 method)

Finalmente, se implementó el procedimiento para resolver \(n\) sistemas lineales de la forma \(AX_i = e_i\) con base en el procedimiento de solucion_lup(), a través de una función denominada inversa().

function inversa(A)
    n = size(A)[1]
    A_inv = Matrix{Float64}(undef, (n, n))
    for i in 1:n
        ei = zeros(n)
        ei[i] = 1
        x = solucion_lup(A, ei)
        for j in 1:n
            A_inv[j, i] = x[j]
        end
    end
    return A_inv
end

inversa (generic function with 1 method)

Al respecto, se realizaron experimentos para ejecutar la función inversa() en matrices cuadradas aleatorias de \(n \times n, \, n = 100, 300, 1000\). En la Tabla 4 se muestra la mediana del tiempo de ejecución obtenido para cada uno de estos experimentos y en la Figura 3 se grafican las distribuciones de dicho tiempo.

lab_inv = BenchmarkGroup()
for i in experimentos
    lab_inv[i] =  @benchmarkable inversa(A) setup = (A = rand(rng, $i, $i))
end
tune!(lab_inv)
lab_inv = run(lab_inv, seconds = segundos)
DataFrame(median(lab_inv), ["n", "Tiempo de ejecución"])

Tabla 4: Tiempo mediano de ejecución de inversa(), para \(n = 100, 300, 1000\)

3×2 DataFrame

Row	n	Tiempo de ejecución
	Int64	TrialEst…
1	300	TrialEstimate(3.509 s)
2	1000	TrialEstimate(604.235 s)
3	100	TrialEstimate(36.127 ms)

plot(lab_inv, yaxis=:log10, xlabel = "n")

Figura 3: Distribución del tiempo de ejecución de inversa(), para \(n = 100,\, 300,\, 1000\)

Un análisis de la función inversa() permite observar que hay dos niveles de ciclos anidamiento (para un orden de crecimiento de \(\Theta(n^2)\)) y un llamado a la función solucion_lup(); por lo que el orden de crecimiento es el mayor de ambos. Por su parte, un análisis de la función solucion_lup() también permite concluir que hay dos niveles de ciclos de anidamiento y un llamado a la función descomposicionlup(). Finalmente, como se explicó en la Sección 3, Cormen et al. (2022, 832) señalan que el Algoritmo 3 implementado a través de la función descomposicionlup() tiene un tiempo de ejecución de \(\Theta(n^3)\).

Por lo tanto, teóricamente se esperaría que la función inversa() tuviera un rendimiento con un orden de crecimiento \(\Theta(n^3)\), que es el mayor de los señalados en el párrafo anterior. Este, en efecto, es el rendimiento esperado para dicha función, de acuerdo con Cormen et al. (2022, 834).

Luego, se esperaría una comparación de rendimientos entre experimentos similar a la que se señaló en la Sección 2. No obstante, las comparaciones de rendimiento se deterioran con una velocidad mucho mayor a la esperada, a medida que crece \(n\), como se muestra a continuación.

for i in eachindex(experimentos)
    if i == firstindex(experimentos)
        continue
    end

    println("Comparación de n = $(experimentos[i]) y n = $(experimentos[i - 1])")
    display(
        judge( mean( lab_inv[experimentos[i]] ),
               mean( lab_inv[experimentos[i - 1]] ))
        )
end

Comparación de n = 300 y n = 100

BenchmarkTools.TrialJudgement: 
  time:   +8869.74% => regression (5.00% tolerance)
  memory: +2562.85% => regression (1.00% tolerance)

Comparación de n = 1000 y n = 300

BenchmarkTools.TrialJudgement: 
  time:   +17102.85% => regression (5.00% tolerance)
  memory: +3588.31% => regression (1.00% tolerance)

Este comportamiento diferente al esperado puede tener la misma explicación que el señalado en la Sección 3, toda vez que dicha solución sirve como parte del procesamiento utilizado para la construcción de la matriz inversa. Por lo tanto, al igual que en la discusión desarrollada en esa sección, se requiere una mayor investigación sobre este comportamiento.

5 Referencias

Cormen, Thomas H., Charles E. Leiserson, Ronald L. Rivest, y Clifford Stein. 2022. Introduction to algorithms. Cuarta Edición. Inglaterra: The MIT Press.

Grossman, Stanley I. 2012. Algebra lineal. México: McGraw Hill.

Téllez Avila, Eric Sadit. 2025. Curso Introductorio al Análisis de Algoritmos con Julia. https://sadit.github.io/ALGO-IR/.

Notas

1. En Téllez Avila (2025, sec. 6) se plantean algoritmos de intersección que pueden ser utilizados para realizar multiplicaciones de matrices dispersas, al representarlas en un arreglo ordenado de pares (fila, valor) o (columna, valor) que solo contiene los valores diferentes de cero. Al intersecar los conjuntos de índices en estas representaciones, se puede realizar la multiplicación matricial únicamente cuando ambos valores son diferentes de cero.