Revista de Matema´tica: Teor´ıa y Aplicaciones 2(2): 35–48 (1995)
algoritmo de karmarkar y matrices ralas
Juan Fe´lix Avila Herrera1
Resumen
Este es el segundo de una serie de dos art´ıculos en los que se estudia el me´todo de
Karmarkar. Se muestra co´mo utilizar la teor´ıa de matrices ralas para obtener una im-
plementacio´n eficiente del proceso de Karmarkar, presentado en el primer art´ıculo. En
la fase I del proceso de Karmarkar, se pone en evidencia la forma como se incrementa
el taman˜o de la matriz de restricciones tecnolo´gicas. La nueva matriz sin embargo,
posee una estructura peculiar bastante favorable, debido a la presencia de bloques de
ceros que la hacen parte de una familia de matrices bastante conocida, a saber las
matrices ralas. Se discute en este trabajo algunas te´cnicas para manejar este tipo de
matrices, y finalmente el autor propone una variante del me´todo de Karmarkar que
aprovecha dicha situacio´n.
Abstract
This is the second of a series of two articles in wich we study the Karmarkar’s
method. In this article we are going to show how can we use sparse matrix theory
to get an effcient implementation of the Karmarkar’s process presented in the first
article. In phase I of the Karmarkar’s process, it was evident how the size of the
technological matrix increased. However, the new matrix has a special structure in
wich we observed the presence of zero’s blocks that make it a sparse matrix. We will
discuss here some techniques to be used with this kind of matrix. Finaly we propose
a Kamarkar’s variant that takes advantage of this situation.
1. Resumen del me´todo de Karmarkar
A continuacio´n se mencionan algunos de los aspectos ma´s importantes sobre el me´todo
de Karmarkar explicado en [2]. Este me´todo permite resolver programas lineales en tiempo
polinomial, mejorando en este sentido el consabido algoritmo del s´ımplex. Es recomendable
en todo caso consultar [2] para tener una mejor visio´n de lo que se expondra´ seguidamente.
El me´todo de Karmarkar consiste en tres fases que se resumen a continuacio´n
1. (Fase I: Convirtiendo un programa lineal a la forma de Karmarkar) En esta fase se
convierte un programa lineal dado en la forma esta´ndar
Maximizar { z = cx : Ax ≥ b, x ≥ 0 } , (1)
1Escuela de Informa´tica, Universidad Nacional
35
36 j. f. avila
a la forma esta´ndar de Karmarkar (FEK):
Minimizar : z = cx, sujeta a: Ax = 0, Unx = 1, x ≥ 0, (2)
donde A es m× n de rango m, n ≥ 2, c y A esta´n compuestos por nu´meros enteros,
Un es un vector con n unos, y cumple las siguientes condiciones (S1) el punto
x0 =
(
1
n , · · · , 1n
)
, es factible, y (S2) el valor objetivo o´ptimo del problema es cero.
Al convertir (1) a la (FEK) obtenemos:
Minimizar : z = x′2(m+n)+1
sujeta a:

H ′x′ = 0,
2(m+n+1)∑
i=1
x′i = 1,
x′i ≥ 0, para 1 ≤ i ≤ 2(m+ n+ 1).
(3)
en donde
H ′ :=
 Am×n − Im Zm×(m+n) αm×1 −bm×1Zn×(m+n) ATn×m In βn×1 −cn×1
c1×n Z1×m −b1×m Z1×n γ 0
 ,
con
α = b+ Um −AUn,
β = c− Un −AT Um,
γ =
∑m
i=1 bi −
∑n
i=1 ci
Como se observo´ en [2] H ′ es una matriz rala.
2. (Fase 2: Resumen del Algoritmo de Karmarkar)
El siguiente es el resumen del algoritmo de Karmarkar que aplica solo para programas
lineales que esta´n bajo la (FEK).
a) Iniciacio´n: (k = 0)
1) Calcule el radio r = 1/
√
n(n− 1).
2) Calcule L = d1 + log(1 + |cjma´x|) + log(|detma´x |)e. Aqu´ı d·e, significa la
parte entera superior y log(·) denota logaritmo en base 2, adema´s |cjma´x| =
ma´x
1≤j≤n
{ |cj | } , y la cantidad |detma´x | es
|detma´x | = ma´x { |det(B)| : B es una base para el problema (2) } .
La constante L se llama la longitud de entrada del problema (2). Una forma
de estimar L es utilizar L˜ en donde:
L˜ := d1 + log(1 + |cjma´x|) + log(1 +m) +
m∑
i=1
n∑
j=1
log(1 + |aij |)e, (4)
y se satisface que L ≤ L˜.
algoritmo de karmarkar y matrices ralas 37
3) Calcule α = (n− 1)/3n.
4) Coloque x0 = (1/n, · · · , 1/n).
b) Elegir Solucio´n: Si cxk < 2−L, la solucio´n actual xk es factible y satisfacto-
ria. Parar el algoritmo.
c) Paso principal: Se definen las siguientes cantidades matriciales:
1) La matriz Dk = diag {xk1, · · · , xkn }, en donde xk = (xk1, · · · , xkn).
2) La matriz P =
[
ADk
Un
]
.
3) El vector c¯ = cDk.
Se calcula entonces x∗ como x∗ = x0−αr cp||cp|| , con cp =
[
I − P T (PP T )−1P ] c¯T .
De esta forma se obtiene xk+1 = (Dkx∗)/(UnDkx∗). Incremente k en uno y
regrese al paso Elegir Solucio´n.
3. (Fase 3: Rutina de redondeo optimal) La solucio´n xk obtenida en la Fase 2, no
necesariamente es un punto extremo. Es por esto que se debe “redondear” xk a una
solucio´n o´ptima, de suerte que tenga tan buen valor objetivo como xk. A continuacio´n
se resume la forma de lograr esto.
a) Elegir Solucio´n: Si al sustituir xk en el programa lineal (2), n restricciones
linealmente independientes esta´n siendo satisfechas en forma activa, se concluye
entonces que xk es una solucio´n ba´sica o´ptima. Parar.
b) Calcular direccio´n: Si xk no es una solucio´n de punto extremo, existen
varias restricciones (digamos l < n) del problema (2) que xk satisface en forma
activa. Denotemos con A˜ la matriz formada a partir de esas restricciones, de la
siguiente forma:
1) Si una restriccio´n de Ax = b se satisface en forma activa con xk, se incluye
la fila correspondiente de A en A˜.
2) Si Unxk = 1, el vector Un = (1, 1, · · · , 1) se incluye en A˜.
3) Si la entrada i-e´sima de xk es nula, se agrega ei = (0, 0, · · · , 1, · · · , 0) (con
un u´nico 1 en la i-e´sima posicio´n) a A˜.
La matriz A˜ resultante es entonces de taman˜o l×n con l < n. Considere ahora
el sistema A˜x = 0. y sea d una solucio´n no trivial de este sistema.
c) Actualizar solucio´n: Se obtiene una nueva solucio´n x∗k, movie´ndose a lo
largo de la direccio´n d si cd < 0 y a lo largo de −d en otro caso2, hasta que
algunas restricciones bloquee´n cualquier movimiento adicional por considera-
ciones de factibilidad. As´ı:
1) si cd < 0, y todas las entradas de d son no negativas, se concluye que la
solucio´n es no acotada. En el caso contrario x∗k = xk + λ · d, en donde
λ = mı´n
1≤i≤n
{−xki/di : di < 0 } , en donde xk = (xk1, xk2, · · · , xkn) y d =
(d1, d2, · · · , dn).
2En virtud de la fase I, basta con considerar el signo de d[n− 1].
38 j. f. avila
2) si cd ≥ 0, y todas las entradas de d son no positivas, se concluye que la
solucio´n es no acotada. En el caso contrario
x∗k = xk − λ · d,
en donde
λ = mı´n
1≤i≤n
{ xki/di : di > 0 } .
Se regresa entonces al paso Elegir Solucio´n.
2. Matrices ralas y eliminacio´n Gaussiana
Cuando una matriz posee “muchas” entradas nulas, se advierte la necesidad de ca-
racterizarlas y buscar estructuras de datos que permitan manejar esta clase especial de
matrices de una forma apropiada.
Definicio´n 1 La densidad δ(A) de una matriz A se define como el cociente entre el nu´mero
de entradas no nulas y el taman˜o total de la matriz.
Definicio´n 2 Se dice que una matriz A es -rala para 0 ≤  ≤ 1, si su densidad δ(A)
satisface δ(A) ≤ .
De esta forma se dira´ que una matriz es rala si es -rala para  cercano a cero. Ma´s
concretamente, se aceptara´ como rala una matriz A (por lo menos en el presente trabajo),
si δ(A) ≤ 1/3.
Se utilizara´n dos estructuras de datos distintas para representar una matriz rala segu´n
convenga. Estas ideas fueron tomadas de [7].
1. (EU): disen˜ada para matrices ralas de caracter´ısticas generales, y que se empleara´ para
almacenar la matriz U en su descomposicio´n L–U. La estructura (EU) la conforman
una lista de vectores { A, Indr, Lonr, Locr }, donde
a) A contiene las entradas no nulas de U , listadas por filas no necesariamente
ordenadas.
b) Indr[i] indica la columna en la que se halla la i-e´sima entrada no nula de U
almacenada en A[i]. Por conveniencia se supondra´ que Indr[i] = 0 si la entrada
A[i] esta´ “vac´ıa” o disponible.
c) Lonr[i] indica la longitud (en te´rmino de entradas no nulas) de la fila i.
d) Locr[i] indica la entrada de A en donde empieza la fila i. Si la i-e´sima fila de
U es nula se coloca Locr[i] en cero.
Es importante notar que los elementos no nulos de la i-e´sima fila son Lonr[i] en
total, y que se hallan ubicados en las posiciones
A[ Locr[i]], · · ·A[ Locr[i] + Lonr[i]− 1],
algoritmo de karmarkar y matrices ralas 39
del vector A.
(Ejemplo) Una representacio´n (EU) de la matriz
U =
 0 2 31 9 4
8 4 7

esta´ dada en la Fig. 1. La entradas 3 y 4 del vector A se consideran “vac´ıas” o
disponibles.
U →

A = 2 3 5 7 8 4 7 1 9 4
Indr = 2 3 0 0 1 2 3 1 2 3
Lonr = 2 3 3
Locr = 1 8 5
Figura 1: Ejemplo de una representacio´n (EU).
2. (EL): disen˜ada para almacenar matrices ralas del tipo triangular inferior, com-
plementando de esta forma, la tarea asignada a la (EU) en la descomposicio´n
L–U. La estructura (EL) esta´ constituida por una lista de vectores de la forma
{ A, Indc, Indr}, que representan una matriz L como sigue:
a) A contiene la entradas no nulas ubicadas debajo de la diagonal de L.
b) La entrada A[i] se halla en la posicio´n ( Indc[i], Indr[i]) de L.
3. Factorizacio´n L–U de una matriz rala
Se tienen ya las estructuras necesarias para almacenar la descomposicio´n L–U de una
matriz rala A de taman˜o n × n. Este proceso consiste en factorizar la matriz A como
A = LU , en donde U es una matriz triangular superior, mientras que L esta´ expresada
como un producto de transposiciones 3 y matrices de tipo especial. Formalizamos estas
ideas en la siguiente definicio´n.
Definicio´n 3 Sean A, L y U matrices de taman˜o n × n. L y U son una descomposicio´n
o factorizacio´n L–U de A si A = LU , y satisfacen que:
a) U es una matriz triangular superior y
3Una transposicio´n es una permutacio´n de dos filas de la matriz identidad.
40 j. f. avila
b) L = (P1M1)(P2M2)(P3M3) · · · (PrMr), en donde r ≤ n(n − 1)/2, Pi es una
transposicio´n y cada factor Mk tiene la forma
Mk = I − µkeikeTjk , k = 1, 2, · · · , r, (5)
donde eik y ejk son vectores unitarios, ik 6= jk, y a los escalares µk se les llama
multiplicadores.
Con el fin de mantener un balance entre estabilidad nume´rica y la preservacio´n de la es-
tructura rala se exigira´ que los multiplicadores dados en (5) este´n acotados de acuerdo con
|µk| ≤ µ¯, para un µ¯ ≥ 1. T´ıpicamente un valor de µ¯ = 10 logra este propo´sito. Recue´rdese
del a´lgebra lineal esta´ndar que una potencia de µ¯ aparece en la cota de crecimiento sobre
los elementos de la factorizacio´n L–U.
El algoritmo que se construira´ para obtener la descomposicio´n L–U de una matriz A
de taman˜o 4 n× n, esta´ basado en la siguiente factorizacio´n:[
α vT
β wT
]
=
[
1 0
−µ 1
] [
α vT
0 w¯T
]
, (6)
en donde µ = −β/α (α 6= 0, β 6= 0) y w¯ = w + µv. La ecuacio´n (6), permite hallar la
descomposicio´n L–U de la matriz A, factorizando progresivamente dos de sus filas. Los
vectores (α, vT ) y (β,wT ) son del mismo taman˜o con una longitud menor o igual a n.
Normalmente (β,wT ), por ejemplo, esta´ formado por entradas de una determinada fila de
A, tomada a partir de su primer elemento no nulo β.
La idea para factorizar la matriz A es sencilla, y se expone someramente a continuacio´n:
1. Inicie L := I, U := A y coloque P1 = P2 = · · · = Pr = I.
2. Al pretender formar un escalo´n 5 en la fila 2 se tienen varias posibilidades:
a) Si A[2, 1] 6= 0 y A[1, 1] 6= 0, se emplea la factorizacio´n (6) para “hacer” un cero
en A[2, 1]. Defina entonces M1 = I −µe2eT1 y reemplace la segunda fila de A, a
saber (β,wT ) por (0, w¯T ).
b) Si A[2, 1] = 0, def´ınase entonces M1 = (I − 0e2eT1 ) = I.
c) Si A[2, 1] 6= 0 pero A[1, 1] = 0, se efectu´a un intercambio entre las filas 1 y 2,
que se registra permutando las filas 1 y 2 de P1. Se llega, por tanto, al caso (b)
anterior.
Colo´quese entonces L := L · (P1M1).
4La ideas son tambie´n va´lidas para matrices rectangulares.
5Entie´ndase por hacer un escalo´n en la fila i, al proceso de hacer ceros las entradas
A[i, 1], · · · , A[i, i− 1],
utilizando para ello operaciones elementales sobre las filas precedentes.
algoritmo de karmarkar y matrices ralas 41
3. Para formar un escalo´n en la fila 3 se utilizan las filas 1 y 2, siguiendo las mismas
ideas presentadas en el punto interior. Este proceso arroja dos nuevos factores M2 y
M3. Se define entonces L := L · (P2M2) · (P3M3).
4. Se procede de igual forma con las filas restantes. Observe que el nu´mero de factores
Mk que se obtiene es 1 + 2 + 3 + · · ·+ (n− 1) = n(n− 1)/2.
A continuacio´n se describen las 10 subrutinas principales que intervienen en la solucio´n
del sistema Ax = b, mediante la factorizacio´n L–U. Los detalles de cada uno de estos
procedimientos se hallan en [1].
1. Comprimir: se encarga de compactar el vector A de la estructura (EU) para
procurar un uso apropiado del espacio. Esta rutina debe utilizarse con moderacio´n.
Si un porcentaje significativo 6 del espacio reservado en memoria principal ha sido
usado, la utilizacio´n de Comprimir es justificable, de otro modo, sencillamente se
estar´ıa adrede retardando el algoritmo que lo invoque.
2. CRF (correr renglo´n al final): se encarga de trasladar “f´ısicamente” el i-e´simo
renglo´n de A (en la estructura (EU)) hasta el final del vector A. Esta rutina es
necesaria puesto que al efectuar operaciones elementales sobre un renglo´n, e´ste puede
incrementar su taman˜o (en te´rmino de entradas no nulas), siendo necesario trasladar-
lo hasta el final del vector A donde puede crecer “ilimitadamente 7”.
3. Apuntar: Para la i-e´sima fila de la matriz A (representada mediante la estructura
(EU)), la subrutina Apuntar se encarga de anotar en un vector n-dimensional v,
las posiciones de A en las que se encuentran sus entradas no nulas, i.e. se carga v
de la siguiente forma:
v[j] =
{
l si A[i, j] = A[l], (A[l] 6= 0, l > 0);
0 en otro caso.
Al tratar de formar un escalo´n en la i-e´sima fila de A, en la eliminacio´n Gaussiana,
el vector v permite escudrin˜ar dicha fila en busca de la pro´xima entrada no nula.
4. Valor U: Este procedimiento se encarga de retornar la posicio´n U [i, j] de una matriz
U , dada en la forma (EU).
5. Valor L: Este procedimiento se encarga de retornar la posicio´n L[i, j] de una matriz
L, dada en la forma (EL).
6. Agregue: Se encarga de agregar un nuevo factor Mk a L. Las transposiciones Pk se
manejan mediante un vector R que se explicara´ despue´s.
6Un 90% puede considerarse razonable.
7En realidad el crecimiento esta´ limitado por la cantidad de memoria reservada o disponible para A,
ya sea que se use programacio´n esta´tica o dina´mica.
42 j. f. avila
7. Escalo´n: Se encarga de formar un escalo´n en una fila cualquiera de la matriz U . Tal
fila se llama la espiga. Escalo´n es la rutina clave de la descomposicio´n L–U.
Es pertinente decir algunas palabras sobre las permutaciones. Cuando se trata de
formar un escalo´n en una fila de una matriz, algunas veces se tropieza con 2 tipos
de dificultades:
a) La cantidad µ = −β/α, no puede ser calculado pues α = 0.
b) El test de estabilidad falla, i.e. |µ| > 10.
En tales casos es necesario intercambiar la fila espiga y la que se estaba usando
para formar un escalo´n en una entrada de e´sta. Tal operacio´n es registrada por dos
vectores que llamaremos P y R. El n-vector P es una permutacio´n, e´sto quiere decir
que si P [i] = j, entonces se intercambiaron las filas i y j. R, por otro lado es un
r-vector (con r = n(n− 1)/2) que se encarga de registrar las transposiciones. Cada
entrada R[i] tiene dos componentes que se llamaran respectivamente R[i].f y R[i].g.
Antes de invocar el procedimiento Escalo´n se inicia el vector R como
R =
([
1
1
]
,
[
2
2
]
, · · · ,
[
r
r
])
. (7)
Si en el ca´lculo de el i-e´simo factor de L, se deben intercambiar las filas j y k, se
modifica R como sigue R[i].f := j, y R[i].g := k. Se puede utilizar un tercer compo-
nente R[i].p, para almacenar solamente la entradas significativas del R propuesto en
(7).
8. Factorizar: Este procedimiento obtiene la factorizacio´n L–U de una matriz, al
aplicar (n − 1) veces la subrutina Escalo´n. Si se producen permutaciones de fi-
las, e´stas se anotan (usando P y R) con el fin de resolver cualquier sistema que
utilice esta descomposicio´n. El lector debe notar que en el caso particular de la ma-
triz P˜ , dada en (14), se puede empezar la eliminacio´n de Gauss a partir de la fila
1 + 2(m+ n+ 1).
9. Evaluar: Si se define L˜ como la matriz triangular inferior dada por 8
L˜ =

1 0 0 · · · 0
µ1 1 0 · · · 0
µ2 µ3 1 · · · 0
...
...
µs µs+1 µs+2 · · · 1
 ,
Evaluar se encarga de regresar U [i, j] si i ≤ j, y L˜[i, j] en caso contrario.
10. Resolver: Esta rutina se encarga de resolver un sistema de la forma Ax = b, uti-
lizando su descomposicio´n L–U. Si la matriz A en cuestio´n es singular, se produce
8Es fa´cil ver que s = 1 + (n− 1)(n− 2)/2.
algoritmo de karmarkar y matrices ralas 43
un mensaje de error y se detiene el proceso. De acuerdo con la Def. 3
A = LU = (P1M1)(P2M2)(P3M3) · · · (PrMr) · U,
en donde r = n(n − 1)/2. De esta forma para resolver el sistema original Ax = b,
sera´ suficiente resolver primero Lz = b y despue´s Ux = z. El sistema Lz = b, se
puede escribir expl´ıcitamente como:
(P1M1)(P2M2)(P3M3) · · · (PrMr)z = b. (8)
El sistema (8) se puede resolver en r ≤ n(n− 1)/2 pasos. Veamos el primero para
comprobar esta aseveracio´n. No´tese que (8) se puede poner como
(P2M2)(P3M3) · · · (PrMr)z =M−11 P−11 b.
Para calcular M−11 P
−1
1 b se debe efectuar primero la transposicio´n P
−1
1 = P1, sobre
el vector b. Esto transforma b en un nuevo vector b1 = P1b. Se nota ahora que
M−11 P
−1
1 b = (I + µ1e2e
T
1 )b
1 = b1 + µ1e2eT1 b
1,
o bien
M−11 P
−1
1 b = b
1 + (µ1b11)e2 =

b11
b12 + µ1b
1
1
...
b1n
 .
Es fa´cil generalizar este paso. En el paso k-e´simo (Mk = (I + µkeseTt )) se debe
entonces efectuar primero la transposicio´n Pk, obteniendo as´ı un nuevo vector bk.
Se actualiza posteriormente la s-e´sima entrada de bk mediante bks := b
k
s + µkb
k
t . La
solucio´n de Ux = z se obtiene mediante la conocida te´cnica de sustitucio´n hacia
atra´s o tambie´n llamada sustitucio´n regresiva, consulte por ejemplo [5].
4. Una implementacio´n eficiente del algoritmo de Karmarkar
En la Fase I del proceso de Karmarkar (consulte [2]) se puso en evidencia co´mo la
matriz tecnolo´gica H ′ del programa lineal transformado, ten´ıa una estructura bastante
rala. De hecho si δ(A) es la densidad de la matriz Am×n del problema original
Maximizar { cx : Ax ≥ b, x ≥ 0 } ,
es posible estimar la densidad de la matriz H ′. El resultado se presenta en la siguiente
proposicio´n.
Proposicio´n 1 La densidad de la matriz tecnolo´gica H ′, del programa lineal transforma-
do (presentado en el primer art´ıculo de esta serie de tres), satisface
δ(H ′) ≤ 2mnδ(A) + 4m+ 4n+ 1
2(m+ n+ 1)2
. (9)
44 j. f. avila
Demostracio´n: El resultado se obtiene simplemente contando el nu´mero de entradas no
nulas en la matriz H ′ dada en (3).
En en efecto, no´tese que si δ(A) = k/mn, en donde k es el nu´mero de entradas no
cero de A, entonces k = mnδ(A), de esta forma k′ (el total de entradas no nulas de H ′)
satisface
k′ ≤ mnδ(A) +m+m+m+mnδ(A) + n+ n+ n+ n+m+ 1,
puesto que en el peor caso b, c, α, β son vectores sin entradas nulas y γ 6= 0. Dado que H ′
es de taman˜o (m + n + 1) × 2(m + n+ 1), se concluye fa´cilmente que la desigualdad (9)
es correcta. 
Es fa´cil verificar nume´ricamente como disminuye la densidad de H ′ conforme aumenta
el taman˜o de la matriz A. Esto hace muy conveniente el uso de estructuras de datos (como
(EU) y (EL)) que sean sensibles a tal feno´meno. Se debe notar tambie´n que si m = n
entonces
δ(H ′) ≤ 2n
2δ(A) + 8n+ 1
2(2n+ 1)2
,
por tanto
l´ım
n→+∞ δ(H
′) ≤ l´ım
n→+∞
2n2δ(A)
8n2
,
de esta forma, en el peor caso (δ(A) = 1) se tiene:
l´ım
n→+∞ δ(H
′) ≤ 0,25.
De hecho, au´n en el caso m < n, la desigualdad de Cauchy 2mn ≤ m2 + n2 da 4mn ≤
(m+ n)2 < (m+ n+ 1)2, as´ı que
δ(H ′) ≤ 2mnδ(A) + 4m+ 4n+ 1
2(m+ n+ 1)2
< 14δ(A) +
2
m+ n+ 1
.
De esta forma para, m y n “grandes”, δ(H ′) es aproximadamente 1/4.
El ana´lisis anterior nos deja concluir que utilizando estructuras de datos como (EU)
y (EL) se ahorra un porcentaje significativo de memoria principal, permitiendo de esta
forma, manejar problemas de mayor taman˜o. Hay sin embargo un paralogismo sobre el que
se debe advertir. El hecho de que la densidad disminuya conforme el taman˜o de problema
lineal aumente, no garantiza de ninguna manera que cualquier programa lineal sea mane-
jable, pues la matriz H ′ puede ser virtualmente “gigantesca”, agotando inexorablemente
la memoria principal. La ventaja de utilizar estructuras como (EU) y (EL) reside en que
no se almacena complementamente la matriz rala, sino sus entradas no nulas. Sin embargo
el 25% de una matriz “gingantesca” puede tal vez no caber en memoria principal.
Se demuestra as´ı que, el enfoque sobre el proceso de Karmarkar adoptado por esta in-
vestigacio´n teo´ricamente rinde frutos. En [1] se muestra como estas cavilaciones funcionan
bien en la pra´ctica, adema´s de hacer ana´lisis de tiempos de ejecucio´n
algoritmo de karmarkar y matrices ralas 45
Se debe observar de la fase II del proceso de Karmarkar, que el esfuerzo computacional
esta´ dominado por el ca´lculo del vector cp (la proyeccio´n ortogonal de c) dado por
cp = [I − P T (PP T )−1P ]c¯T ,
(consu´ltese el primer art´ıculo de esta serie). Debemos entonces, encontrar la inversa de
PP T o bien resolver un sistema lineal de la forma 9
(PP T )x = P c¯. (10)
Sustituyendo 10 P =
[
AD
Un
]
, se obtiene
PP T =
[
AD2AT AD Un
(AD Un)T Un Un
]
=
[
AD2AT 0
0 n
]
. (11)
(Es importante hacer e´nfasis aqu´ı en que la matriz A que se menciona en P , es la matriz
H ′ de la ecuacio´n 3).
No´tese de (11) que la matriz PP T no es necesariamente rala, au´n cuando A y conse-
cuentemente P lo sean. Un manejo despreocupado en la solucio´n de (10), echar´ıa por tierra
el ahorro de espacio en memoria principal que se ha venido pregonando, adema´s de que
obligar´ıa a disen˜ar rutinas para efectuar operaciones elementales con matrices utilizando
las estructuras (EU) y (EL).
Se debe notar que el sistema (10) es equivalente a[
I P T
P 0
] [
y
x
]
=
[
c¯
0
]
. (12)
En efecto, no´tese que (12) es equivalente a{
y + P Tx = c¯,
P y = 0,
de esta forma
y = c¯− P Tx, (13)
y por tanto al multiplicar ambos miembros de (13) por P , se obtiene 0 = Py = P c¯−PP Tx,
que es justamente (10). Es justo aclarar que esta idea (ecuacio´n (12)), fue tomada de [6].
Es pertinente hacer algunos comentarios sobre este enfoque. Lo primero es que se evita
disen˜ar algoritmos para efectuar productos matriciales utilizando las estructuras de datos
(EL) y (EU). Esto reduce el tiempo computacional. Por otro lado, la matriz
P˜ :=
[
I P T
P 0
]
, (14)
hereda la estructura rala de P , y e´sto es algo deseable. De hecho si δ(A) es la densidad de la
matriz Am×n del problema original Maximizar { cx : Ax ≥ b, x ≥ 0 } , es posible estimar
la densidad del la matriz P˜ . El resultado se resume en la siguiente proposicio´n.
9Obs´ervese que esta es la “ecuacio´n normal” para el problema de mı´nimos cuadrados: hallar x que
minimiza ||P Tx− c¯T ||. Ve´ase [8].
10Se usara´ por comodidad D en lugar de Dk.
46 j. f. avila
Proposicio´n 2 La densidad de la matriz P˜ dada en (14) satisface
δ(P˜ ) ≤ 4mnδ(A) + 14m+ 14n+ 8
(3m+ 3n+ 4)2
. (15)
Demostracio´n: Se debe notar primero que
δ(P ) ≤ 2mnδ(A) + 4m+ 4n+ 1 + 2(m+ n+ 1)
2(m+ n+ 1)(m + n+ 2)
=
2mnδ(A) + 6m+ 6nA+ 3
2(m+ n+ 1)(m+ n+ 2)
,
entonces dado que P˜ =
[
I P T
P 0
]
, se tiene que:
δ(P˜ ) ≤ 2[2mnδ(A) + 6m+ 6n+ 3] + 2(m+ n+ 1)
(3m+ 3n+ 4)2
,
y esto es justamente lo que se quer´ıa demostrar. 
De nuevo se puede analizar nume´ricamente el comportamiento de δ(P˜ ) conforme au-
menta el taman˜o de la matriz A. Se notara´ lo conveniente que resulta el uso de las es-
tructuras (EL) y (EU) para manejar esta matriz. Obse´rvese, adema´s que si m = n,
entonces
l´ım
n→+∞ δ(P˜ ) ≤ 4/36 = 1/9 ≈ 0,111,
y de nuevo, el estimado 4mn < (m+n+1)2, permite concluir que δ(P˜ ) es aproximadamente
1/9 para m y n “grandes”.
La principal desventaja de este enfoque es el incremento en el uso de memoria principal
al reemplazar P por P˜ . Pareciera estar en contraposicio´n con nuestra lucha a “ultranza”
por el ahorro de este recurso. Sin embargo, dado que P˜ es rala, e´ste es un precio razonable
que se puede pagar, y que redunda en una mejora de la velocidad del algoritmo final.
No´tese por ejemplo, que si la matriz tecnolo´gica del problema original es de taman˜o
50 × 50, la matriz P˜ correspondiente es de taman˜o 304 × 304, aproximadamente unas 37
veces ma´s grande que la original. Sin embargo, dado que en este caso P˜ tiene densidad
aproximada de 0.123, se puede ahorrar un 87.7% del espacio en memoria principal. De
esta forma al utilizar estructura de datos como (EU), la matriz P˜ , ocupa solamente unas
4.5 veces ma´s espacio que la matriz original.
La “variante” del algoritmo de Karmarkar, que se propone en esta investigacio´n,
sera´ llamada (K & MR) (Karmarkar y matrices ralas), y puede ser explicada ahora
con ma´s propiedad. Las tres fases del me´todo de Karmarkar (consulte [2]) permanecen
ide´nticas. La modificacio´n estriba en el uso de estructuras de datos sensibles a las carac-
ter´ısticas de las matrices ralas presentes a lo largo del proceso. En efecto, el punto a´lgido
de esta concepcio´n se situ´a en la solucio´n del sistema (10), utilizando un enfoque (como el
adoptado) que pretende un ahorro considerable de espacio en memoria principal, sin que
e´sto vaya en detrimento significativo de la rapidez del algoritmo final.
algoritmo de karmarkar y matrices ralas 47
(K & MR) es una implementacio´n eficiente del algoritmo de Karmarkar, en el sentido
de que permite manejar problemas ma´s grandes que los tolerados por una implantacio´n
rasa (usando arreglos bidimensionales) de este algoritmo. La consigna que ha marcado
la directriz en el desarrollo de esta investigacio´n, ha sido construir una versio´n del algo-
ritmo Karmarkar aplicable a programa lineales “grandes”, administrando eficientemente
la memoria principal del computador. Esto no quiere decir que no se haya cuidado la
eficiencia de los procedimientos involucrados, no obstante, cuando se debio´ escoger entre
velocidad y memoria RAM, casi siempre se sacrifico´ velocidad. Un aspecto en que se ev-
idencia perfectamente este lineamiento, es en el manejo de la matriz tecnolo´gica original.
En efecto, dicha matriz se carga directamente de disco duro, cada vez que sea necesaria y,
como se explica seguidamente, esto se da en cada iteracio´n requerida en las fases II y III.
Se debe reconocer que d´ıa con d´ıa las limitaciones de memoria RAM van disminuyendo y
esto puede ser una cr´ıtica al enfoque adoptado aqu´ı, sin embargo las ideas plasmadas en
este trabajo permitira´n utilizar con ma´s eficiencia cualquier cota de memoria RAM que
se proponga.
La matriz P˜ , se construye usando como piedra angular, la matriz tecnolo´gica original
A. Dado que se debe resolver el sistema (10), la invocacio´n del procedimiento que obtiene
la factorizacio´n L–U de P˜ , destruye esta matriz, por lo que se hace necesario construirla
despue´s de cada invocacio´n de esta rutina. Al analizar la matriz P˜ , se observa que A,
esta´ “presente” 4 veces en e´sta, no obstante basta con cargarla una sola vez del disco
duro. Este razonamiento no es tan obvio como suena, cuando se usa la estructura de
dados (EU). El procedimiento para cargar P˜ utilizando (EU) puede ser el siguiente:
1. Colo´quese la matriz Am×n entre las filas 2m+2n+3 y 3m+2n+2 de P˜ , empezando
en al columna 1 y terminado (obviamente) en la columna n. Llamaremos a este
bloque (principal), (BP).
Es importante dejar 3 campos disponibles despue´s de colocar cada fila de A en P˜ ,
dado que la filas, comprendidas entre la 2m + 2n + 3 y la 3m + 2n + 2 no han
sido construidas del todo. En efecto, existen tres entradas (por cada una de estas
filas) que pueden ser no cero: la que determina Im (que obviamente es no nula), la
correspondiente en αm×1 y la determinada por −bm×1.
2. Calcu´lesen los vectores α y β utilizando la matriz P˜ construida hasta ahora. Dado
que A esta´ contenida en P˜ , este ca´lculo requiere simplemente acceder las entradas
apropiadas de P˜ .
3. Comple´tense las filas comprendidas entre las 2m+2n+3 y 3m+ 2n+2, utilizando
los vectores α y b.
4. Mediante el uso del bloque (BP) de P˜ , constru´yase el resto de la matriz.
5. Para terminar de formar P˜ , se debe multiplicar la i-e´sima columna del bloque for-
mado por las filas comprendidas entre la 2m+2n+3 y 3m+ 2n+2 y las columnas
ubicadas entre la 1 y la 2(m + n+ 1), por la i-e´sima entrada de xk (la solucio´n de
turno).
48 j. f. avila
Es importante notar que a la hora de resolver el sistema en donde se involucra P˜ , se
puede empezar la eliminacio´n de Gauss a partir de la fila 2m+ 2n+ 3, dado que las filas
precedentes no necesitan ser procesadas.
En la fase III del me´todo de Karmarkar, tambie´n es necesario resolver sistemas de
ecuaciones para purificar la solucio´n obtenida en la fase II. En este caso, sin embargo, la
matriz (que llamaremos R) de cada sistema, es ma´s sencilla que P˜ , y de hecho tambie´n
ma´s pequen˜a. Dado que R tambie´n utiliza la matriz tecnolo´gica original A, es necesario
cargarla directamente desde el disco duro, esto se puede hacer coloca´ndola de una vez en
las primeras m filas y las primeras n columnas de R. Dado que R no es rectangular, se le
puede completar hasta obtener una matriz cuadrada, utilizando renglones adicionales ge-
nerados aleatoriamente, pero cuidando que tengan “suficientes” entradas nulas, de suerte
que la estructura (EU) siga siendo apropiada.
(K & MR) tiene una virtud en la que se ha hecho poco e´nfasis, a saber la obtencio´n
concomitante de la solucio´n para los problemas dual y primario. En efecto, la fase I del
me´todo de Karmarkar que busca transformar un problema dado en la forma cano´nica
a la forma esta´ndar de Karmarkar, indirectamente permite la obtencio´n de la solucio´n
del problema dual. Si tanto el problema primal como el dual son o´ptimos, se sabe que
ambos objetivos o´ptimos son iguales (esto es en virtud de las condiciones de optimalidad
de Kuhn–Tucker, consu´ltese [4]). Para obtener entonces la solucio´n o´ptima del problema
dual, bastara´ con utilizar las posiciones ubicadas entre la m + n + 1 y la 2m + n de la
solucio´n o´ptima para el problema transformado (presentado en el segundo art´ıculo de esta
serie). En [1] el lector encontrara´ los detalles de (K & MR)., adema´s ah´ı se efectuan
experimentos tendientes a valorar las posibilidades de (K & MR), compara´ndolo con el
me´todo del s´ımplex.
Otro aspecto que se debe mencionar es la utilizacio´n del procesamiento en paralelo
para acelerar la obtencio´n de soluciones. Ser´ıa agradable que algu´n lector pueda estar
interesado en continuar con esta l´ınea de investigacio´n.
Referencias
[1] J. F. Avila Herrera (1993) Una implementacio´n eficiente del algoritmo de Karmarkar .
Tesis de Maestr´ıa ITCR, Cartago, C.R..
[2] J. F. Avila Herrera (1995) “Me´todo de Karmarkar”, Revista de Matema´ticas 1
[3] M. S. Bazaraa y J. J. Jarvis (1984) Programacio´n Lineal y Flujo de Redes, 1era
Edicio´n. Limusa, Me´xico.
[4] M. S. Bazaraa, J. J. Jarvis and H. Sherali (1990) Linear Programming and Network
Flows, 2nd Edition. John Wiley & Sons, N.Y..
[5] R. L. Burden y J. D. Douglas (1985) Ana´lisis Nume´rico. Grupo Editorial Iberoame´ri-
ca, Me´xico.
[6] I. E. Duff and J. K. Reid (1989) Direct Methods for Sparse Matrices. Oxford Science
Publications, Oxford.
algoritmo de karmarkar y matrices ralas 49
[7] P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright (1986) Mantaining LU
Factors of a General Sparse Matrix . Department of Operations Research, Stanford
University, California.
[8] G. H. Golub and C. F. Van Loan (1989) Matrix Computations, 2nd Edition. The
Johns Hopkins University Press, USA.
[9] N. Karmarkar (1984) “A new polynomial-time algorithm for Linear Programming”,
Combinatorica 4
[10] M. Nun˜ez y V. Kong (1992) El problema de la Mezcla de Alimentos. Reporte Te´cnico,
ITCR.