Pivoteo y factorización LU

01 · Mejorar la eliminación

El orden de las filas importa numéricamente

La eliminación de Gauss divide por los pivotes de la diagonal. Un pivote cero detiene el algoritmo básico y uno muy pequeño puede producir multiplicadores enormes, aunque el sistema tenga una solución única.

Si |a[k][k]| es diminuto, pequeñas perturbaciones de los datos o del redondeo pueden amplificarse. Intercambiar filas antes de dividir permite elegir un coeficiente más apropiado.

02 · Motivación

Cuando el primer pivote es cero

Consideremos la matriz aumentada:

No podemos calcular m₂₁=1/0. Sin embargo, intercambiando las filas:

La solución se obtiene inmediatamente: y=2 y x=1. El cero inicial era un problema del orden, no una prueba de singularidad.

El mismo principio se aplica a pivotes no nulos pero demasiado pequeños.

03 · Estrategia principal

Pivoteo parcial

En la etapa k, el pivoteo parcial busca el elemento de mayor valor absoluto en la columna k, desde la fila k hacia abajo:

Si p≠k, intercambia las filas p y k. Después continúa la eliminación normalmente.

1. 1
   Buscar.

   Recorrer la parte activa de la columna y localizar la mayor magnitud.

2. 2
   Comprobar.

   Si incluso el mejor candidato es numéricamente cero, la matriz es singular para la tolerancia usada.

3. 3
   Intercambiar.

   Mover la fila seleccionada a la posición del pivote.

4. 4
   Eliminar.

   Calcular multiplicadores y anular los elementos inferiores.

Esta política añade búsquedas de costo O(n²), inferior al costo cúbico de la eliminación, y suele ofrecer una estabilidad muy buena en la práctica.

04 · Variantes

Pivoteo escalado y pivoteo completo

En el pivoteo escalado se suele definir sᵢ=maxⱼ|aᵢⱼ|. El pivoteo completo requiere registrar también una permutación de variables y conduce a una relación del tipo PAQ=LU.

05 · Registrar intercambios

La matriz de permutación P

Una matriz de permutación se obtiene al reordenar las filas de la identidad. Multiplicar P por A aplica a A esos mismos intercambios.

La matriz anterior intercambia las dos primeras filas. En código no hace falta almacenar todos sus ceros y unos: un vector [1,0,2] indica que la nueva fila 0 proviene de la fila 1, la nueva fila 1 de la fila 0 y la tercera permanece igual.

function permutarVector(permutacion, b) {
  if (permutacion.length !== b.length) {
    throw new Error("Dimensiones incompatibles");
  }
  return permutacion.map(indice => b[indice]);
}

console.log(permutarVector([1, 0, 2], [10, 20, 30]));
// [20, 10, 30]

06 · Separar el trabajo

Qué significa factorizar A=LU

La eliminación puede expresarse como el producto de dos matrices triangulares:

L: triangular inferior

Contiene unos en la diagonal y los multiplicadores utilizados durante la eliminación.

U: triangular superior

Es la matriz que queda después de la eliminación hacia adelante.

Si hubo intercambios de filas, la relación correcta es:

Esta forma separa el trabajo costoso de triangularización de las sustituciones económicas que dependen del término independiente.

07 · Construcción manual

Guardar los multiplicadores en L

Consideremos una matriz que no necesita intercambios para este ejemplo:

Los primeros multiplicadores son m₂₁=2 y m₃₁=−1. Después de la primera columna, el multiplicador de la segunda etapa es m₃₂=−1. El resultado es:

Al multiplicar LU recuperamos exactamente A. Los números que Gauss usa temporalmente como multiplicadores se convierten en información reutilizable.

08 · Dos sustituciones

Resolver Ax=b a partir de LU

Si PA=LU, entonces:

Introducimos un vector intermedio y y resolvemos dos sistemas triangulares:

1. Sustitución progresiva: resolver Ly=Pb desde la primera fila hacia abajo.
2. Sustitución regresiva: resolver Ux=y desde la última fila hacia arriba.

Para la matriz del ejemplo y b=[7,−8,18]ᵀ, obtenemos y=[7,−22,3]ᵀ y finalmente x=[1,2,3]ᵀ.

09 · Algoritmo principal

Factorización LU con pivoteo parcial

Para ahorrar memoria, podemos guardar U en la diagonal y parte superior de una sola matriz, y los multiplicadores de L debajo de la diagonal. Los unos de L quedan implícitos.

function factorizarLU(A, tolerancia = 1e-12) {
  const n = A.length;
  if (n === 0 || A.some(fila => fila.length !== n)) {
    throw new Error("A debe ser una matriz cuadrada no vacía");
  }
  if (A.some(fila => fila.some(valor => !Number.isFinite(valor)))) {
    throw new Error("Los coeficientes deben ser finitos");
  }

  const LU = A.map(fila => [...fila]);
  const P = Array.from({ length: n }, (_, i) => i);
  const escala = Math.max(...A.flat().map(Math.abs));
  if (escala === 0) throw new Error("La matriz es singular");
  let intercambios = 0;

  for (let k = 0; k < n; k++) {
    let filaPivote = k;
    for (let i = k + 1; i < n; i++) {
      if (Math.abs(LU[i][k]) > Math.abs(LU[filaPivote][k])) {
        filaPivote = i;
      }
    }
    if (Math.abs(LU[filaPivote][k]) <= tolerancia * escala) {
      throw new Error("La matriz es singular o casi singular");
    }
    if (filaPivote !== k) {
      [LU[k], LU[filaPivote]] = [LU[filaPivote], LU[k]];
      [P[k], P[filaPivote]] = [P[filaPivote], P[k]];
      intercambios++;
    }

    for (let i = k + 1; i < n; i++) {
      LU[i][k] /= LU[k][k];
      for (let j = k + 1; j < n; j++) {
        LU[i][j] -= LU[i][k] * LU[k][j];
      }
    }
  }
  return { LU, P, intercambios };
}

const ejemploA = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]];
console.log(factorizarLU(ejemploA));

Cuando se intercambian filas completas de la matriz combinada, también se reordenan los multiplicadores guardados en etapas anteriores.

10 · Reutilización

Resolver con una factorización existente

// Versión compacta para que el ejemplo sea ejecutable por separado.
function factorizarLU(A) {
  const LU = A.map(fila => [...fila]);
  const P = Array.from({ length: A.length }, (_, i) => i);

  for (let k = 0; k < A.length; k++) {
    let pivote = k;
    for (let i = k + 1; i < A.length; i++) {
      if (Math.abs(LU[i][k]) > Math.abs(LU[pivote][k])) pivote = i;
    }
    if (Math.abs(LU[pivote][k]) <= 1e-12) {
      throw new Error("Matriz singular o casi singular");
    }
    [LU[k], LU[pivote]] = [LU[pivote], LU[k]];
    [P[k], P[pivote]] = [P[pivote], P[k]];
    for (let i = k + 1; i < A.length; i++) {
      LU[i][k] /= LU[k][k];
      for (let j = k + 1; j < A.length; j++) {
        LU[i][j] -= LU[i][k] * LU[k][j];
      }
    }
  }
  return { LU, P };
}

function resolverLU(factorizacion, b) {
  const { LU, P } = factorizacion;
  const n = LU.length;
  if (b.length !== n || b.some(valor => !Number.isFinite(valor))) {
    throw new Error("Vector b inválido");
  }

  const pb = P.map(indice => b[indice]);
  const y = Array(n).fill(0);
  for (let i = 0; i < n; i++) {
    let suma = 0;
    for (let j = 0; j < i; j++) suma += LU[i][j] * y[j];
    y[i] = pb[i] - suma; // La diagonal de L contiene unos.
  }

  const x = Array(n).fill(0);
  for (let i = n - 1; i >= 0; i--) {
    let suma = 0;
    for (let j = i + 1; j < n; j++) suma += LU[i][j] * x[j];
    x[i] = (y[i] - suma) / LU[i][i];
  }
  return x;
}

const A = [[2, 1, 1], [4, -6, 0], [-2, 7, 2]];
const b = [7, -8, 18];
const factorizacion = factorizarLU(A);
console.log(resolverLU(factorizacion, b)); // [1, 2, 3]

La factorización cuesta O(n³), pero cada nuevo vector b necesita solamente dos sustituciones de costo O(n²).

11 · Aprovechar la factorización

Múltiples términos independientes y determinante

Varios vectores b

En simulaciones temporales, calibraciones o cálculos por escenarios, la matriz A puede permanecer fija mientras cambia b. Factorizar una vez evita repetir toda la eliminación.

Determinante

Como la diagonal de L contiene unos, det(L)=1. Cada intercambio de filas cambia el signo del determinante:

Aquí s es la cantidad de intercambios. Este cálculo aprovecha la factorización sin expandir el determinante por cofactores.

Matriz inversa

Es posible resolver AX=I columna por columna usando la misma LU. Sin embargo, para resolver Ax=b no conviene construir explícitamente A⁻¹: las sustituciones son más directas.

12 · Control de calidad

Verificar PA≈LU y Ax≈b

Una implementación debe comprobar dos niveles diferentes:

• Error de factorización: medir una norma de PA−LU.
• Residuo de solución: medir una norma de b−Ax.

Las formas relativas permiten comparar problemas con escalas distintas. Un error pequeño confirma coherencia computacional, aunque no elimina el efecto de un mal condicionamiento.

13 · Rendimiento

Costo, almacenamiento y decisiones prácticas

En software científico conviene utilizar bibliotecas optimizadas: aprovechan memoria caché, operaciones por bloques, paralelismo y rutinas numéricas ampliamente probadas.

14 · Exploración interactiva

Laboratorio: construir la matriz LU combinada

Elegí una matriz, una estrategia y una etapa. Sobre la diagonal se forma U; debajo se almacenan los multiplicadores de L. El laboratorio usa un vector b=A[1,2,−1]ᵀ para comprobar la solución.

Preguntas para explorar

1. En la matriz general, compará la permutación con y sin pivoteo.
2. En «Pivote inicial cero», observá por qué el algoritmo básico se detiene.
3. En «Filas con escalas distintas», compará pivoteo parcial y escalado.
4. En la matriz singular, identificá la etapa donde desaparece un pivote válido.

15 · Comprensión

Confusiones frecuentes

16 · Práctica guiada

Ejercicios

Ejercicio 1: elegir el pivote

En la columna activa [−0.2, 5, −3]ᵀ, ¿qué fila elige el pivoteo parcial?

Ejercicio 2: interpretar P

¿Qué efecto tiene el vector de permutación [2,0,1] sobre b=[10,20,30]?

Ejercicio 3: resolver con LU

Para las matrices L y U del ejemplo, verificá que b=[7,−8,18] conduce a x=[1,2,3].

Ejercicio 4: determinante

Si diag(U)=[2,−8,1] y no hubo intercambios, calculá det(A).

Ejercicio 5: desempaquetar L y U

Implementá una función que separe una matriz LU combinada.

function separarLU(LU) {
  const n = LU.length;
  const L = Array.from({ length: n }, () => Array(n).fill(0));
  const U = Array.from({ length: n }, () => Array(n).fill(0));

  for (let i = 0; i < n; i++) {
    for (let j = 0; j < n; j++) {
      if (i > j) L[i][j] = LU[i][j];
      else U[i][j] = LU[i][j];
    }
    L[i][i] = 1;
  }
  return { L, U };
}

console.log(separarLU([[2, 1], [2, -8]]));

17 · Síntesis

Ideas para recordar

En el próximo tema estudiaremos los métodos iterativos de Jacobi y Gauss-Seidel, especialmente útiles para sistemas grandes donde construir una factorización completa puede ser costoso.