Factorización LU

Ya saben, el algoritmo que permite descomponer una matriz como producto de dos matrices triangulares (superior e inferior). Eso sí, hace falta que la matriz original cumpla ciertos requisitos, claro.

Hasta donde he podido averiguar fue A. M. Turing quien propuso por primera vez esta factorización en su artículo Rounding-off errors in matrix processes. [Ref. 1]

Concretamente, en el punto 3. de dicho artículo: theorem on triangular resolution.

El teorema en cuestión afirma lo siguiente:

If the principal minors of the matrix A are non-singular, then there is a unique unit lower triangular matrix L, a unique diagonal matrix D, with non-zero diagonal element d_{j} and a unique unit upper triangular matrix U such that A = LDU.

Más claro, la matriz A puede descomponerse en el producto de estas otras tres:

\begin{pmatrix} 1 & 0 & 0 & \cdots & 0\\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & 1 \end{pmatrix} \begin{pmatrix} d_{11} & 0 & 0 & \cdots & 0\\ 0 & d_{22} & 0 & \cdots & 0 \\ 0 & 0 & d_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & d_{nn} \end{pmatrix}  \begin{pmatrix} 1 & u_{12} & u_{13} & \cdots & u_{1n}\\ 0 & 1 & u_{23} & \cdots & u_{2n} \\ 0 & 0 & 1 & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{pmatrix}

A continuación, Turing pasa a demostrarnos dicho Teorema. Vamos con ello.

The kth diagonal element of D will be denoted by d_{k}

Por simplicidad de notación, Turing denomina d_{k} al elemento d_{kk} de la matriz diagonal D.

The l_{k} coefficient of the equation A = LDU gives us l_{11}d_{1}u_{1k} = a_{1k}

Una forma rápida de llegar a esta ecuación es darnos cuenta de que LD es el producto de dos matrices triangulares inferiores (una matriz diagonal es a la vez una matriz triangular inferior). Por tanto su producto es también una matriz triangular inferior.

Por lo que la primera \, fila \, de \, la \, matriz \, LD = [l_{11}d_{1} \, 0 \, 0 \, ... \, 0]

Si ahora multiplicamos esta primera \, fila \, de \, LD por la matriz U (que es triangular superior)

primera \, fila \, de \, LDU = [l_{11}d_{1} \, 0 \, 0 \, ... \, 0] \begin{pmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n}\\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{pmatrix}

Es decir,

primera \, fila \, de \, LDU = [l_{11}d_{1}u_{11} \, \, \, l_{11}d_{1}u_{12} \, \, \, l_{11}d_{1}u_{13} \, ... \, l_{11}d_{1}u_{1n}]

Recordando que A = LDU, hemos llegado a la expresión

a_{1k} = l_{11}d_{1}u_{1k}

and since l_{11} = u_{11} = 1 this determines d_{1} to be a_{11} and u_{1k} to be a_{1k}/d_{1}, these choices satisfy the equations in question.

Es decir, dada la matriz A, sabemos que d_{1} = a_{11} y que la primera fila de la matriz U es u_{1k} = a_{1k}/d_{1} = a_{1k}/a_{11}

Suppose now that we have found values of l_{ij}, u_{jk} with j < i_{0} (that is, we have found the first i_{0} - 1 rows of L and columns of U) and the first i_{0} - 1 diagonal elements d_{k},  so that the equations arising from the first i_{0} - 1 rows of the equation A = LDU are satisfied; and suppose further that these choices are unique and d_{k} \not= 0

Aquí hay que recordar que la matriz L es triangular inferior unitaria, por lo que l_{{(i_{0}-1)}{(i_{0}-1)}}= 1 y l_{{(i_{0}-1)}{(k)}} = 0 con k \textgreater i_{0} - 1.

Ídem, pero para columnas, en la matriz triangular superior unitaria U.

It will be shown how the next row of L and the next column of U, and the next diagonal element d_{i_{0}} \not= 0 are to be chosen so as to satisfy the equations arising from the next row of A = LDU, and that the choice is unique.

Asumiendo lo anterior, Turing va a demostrar a continuación cómo escoger la siguiente fila (la fila i_{o}) de la matriz L y el siguiente elemento de la matriz D (es decir, d_{i_{0}}) para que satisfagan la ecuación A = LDU. Además demostrará que dicha elección es única.

The equations to be satisfied in fact state

$latex l_{i_{0}i_{0}}d_{i_{0}}u_{i_{0}k} = a_{i_{0}k} – \sum\limits_{j \textless i_{0}}{l_{i_{0}j}d_{j}u_{jk}} \, \, \, \, (k \ge i_{0})$

$latex l_{i_{0}k}d_{k}u_{kk} = a_{i_{0}k} – \sum\limits_{j \textless k}{l_{i_{0}j}d_{j}u_{jk}} \, \, \, \, (k \textless i_{0})$

Para los perezosos, les dejo aquí cómo se han calculado las expresiones anteriores.

Recordemos que A = LDU. Llamemos B = LD

Por definición, b_{ij} = \sum\limits_{k} {l_{ik}d_{kj}}

Como la matriz D es una matriz diagonal, d_{kj} = 0 para k \not= j

Por lo tanto: b_{ij} = l_{ij}d_{j} (ec. 1)

Ahora, la matriz A = BU. Aplicando nuevamente la definición de producto de matrices:

a_{pq} = \sum\limits_{r}{b_{pr}u_{rq}}

Sustituyendo la ec. 1,

a_{pq} = \sum\limits_{r}{l_{pr}d_{r}u_{rq}}

Una simple sustitución de los índices nos lleva a la fórmula para la fila i_{o}

a_{i_{0}k} = \sum\limits_{j}{l_{i_{0}j}d_{j}u_{jk}} (ec. 2)

Las dos ecuaciones que Turing nos muestra en su artículo se obtienen de la ec. 2, teniendo en cuenta la triangularidad de las matrices L y U, así:

Primer caso: por ser L una matriz triangular inferior, l_{ij} = 0 para j > i

a_{i_{0}k} = {\sum\limits_{j \le {i_{0}}}}{l_{i_{0}j}d_{j}u_{jk}} = l_{i_{0}i_{0}}d_{i_{0}}u_{i_{0}k} + {\sum\limits_{j < {i_{0}}}{l_{i_{0}j}d_{j}u_{jk}}}

Segundo caso: por ser U una matriz triangular superior, u_{ij} = 0 para i > j

a_{i_{0}k} = {\sum\limits_{j \le k}}{l_{i_{0}j}d_{j}u_{jk}} = l_{i_{0}k}d_{k}u_{kk} + {\sum\limits_{j < {k}}{l_{i_{0}j}d_{j}u_{jk}}}

The right-hand sides of these equations are entirely in terms of quantities already determined.

Al conocer todos los términos de la parte derecha de las ecuaciones podemos resolverlas. Turing analiza los tres casos: k = i_{0}, k > i_{0} y k < i_{0} y demuestra que existe una solución que además es única.

When k = i_{0} the first equation is satisfied and can only be satisfied by putting d_{i_{0}} = right-hand side, determining d_{i_{0}}

En efecto, si hacemos k = i_{0} tenemos

$latex l_{i_{0}i_{0}}d_{i_{0}}u_{i_{0}i_{0}} = a_{i_{0}i_{0}} – \sum\limits_{j \textless i_{0}}{l_{i_{0}j}d_{j}u_{ji_{0}}}$

Recordando que l_{i_{0}i_{0}} = u_{i_{0}i_{0}} = 1, obtenemos

$latex d_{i_{0}} = a_{i_{0}i_{0}} – \sum\limits_{j \textless i_{0}}{l_{i_{0}j}d_{j}u_{ji_{0}}}$

The equations for k > i_{0} can then be satisfied by one and only one set of values of u_{i_{0}k}, provided d_{i_{0}} \neq 0

En efecto, la ecuación para k > i_{0} es

$latex l_{i_{0}i_{0}}d_{i_{0}}u_{i_{0}k} = a_{i_{0}k} – \sum\limits_{j \textless i_{0}}{l_{i_{0}j}d_{j}u_{jk}}$

Sabemos que l_{i_{0}i_{0}} = 1 y puesto que ya hemos calculado d_{i_{0}}, entonces

$latex u_{i_{0}k} = (a_{i_{0}k} – \sum\limits_{j \textless i_{0}}{l_{i_{0}j}d_{j}u_{jk}}) / d_{i_{0}}$

The equations for k < i_{0} can also be satisfied by one and only one set of values of l_{i_{o}k}, since each d_{k} is different from 0.

En efecto, para el caso k < i_{0} se aplica la ecuación

$latex l_{i_{0}k}d_{k}u_{kk} = a_{i_{0}k} – \sum\limits_{j \textless k}{l_{i_{0}j}d_{j}u_{jk}} \, \, \, \, (k \textless i_{0})$

Pero como u_{kk} = 1,

l_{i_{0}k}d_{k} = a_{i_{0}k} - {\sum\limits_{j < {k}}{l_{i_{0}j}d_{j}u_{jk}}}

Despejando,

l_{i_{0}k} = (a_{i_{0}k} - {\sum\limits_{j < {k}}{l_{i_{0}j}d_{j}u_{jk}}}) / d_{k}

donde sabemos que cada d_{k} \not= 0, y por tanto el conjunto de valores l_{i_{0}k} es único.

En todo lo anterior quedaba por comprobar que d_{i_{0}} \not= 0, pero esto ha de ser así, ya que

The new diagonal element d_{i_{0}} is not 0 because the i_{0}th principal minor of A is equal to the product of the first i_{0} diagonal elements d_{k}.

***************************

Nota: habrán observado que en un caso se ha demostrado que los valores de la matriz U son únicos y en otro caso se demuestra que los valores de L son únicos. Es obvio que si los valores de U son únicos, también lo son los de L y viceversa.

****************************

COMENTARIOS

1. Turing empleó la factorización LDU como punto de partida para analizar (entre otros) el método de eliminación de Gauss. Sin embargo, este modo de analizar el problema es muy fructífero y ha dado lugar a muchas otras factorizaciones útiles.

“The point of view taken here in which Gaussian elimination is regarded as being accomplished by premultiplications by elementary lower triangular matrices originated with Turing (1948) and has been extensively exploited by Wilkinson (AEP, for example). This approach has the advantage of generalizing readily; any other computationally convenient class of matrices may be used to introduce zeros.” [Ref. 2]

“The process of replacing the rows of a matrix by a linear combination of other rows may be regarded as left-multiplication of the matrix by another matrix, this second matrix having coefficients which describe the linear combinations required.” [Ref. 3]

2. En el capítulo 3.1 del libro Introduction to matrix computations [Ref. 4] pueden ver una demostración más moderna (y sin tantos subíndices) del teorema aquí expuesto.

3. Finalmente, y con respecto a la inversión de matrices, les recomiendo que lean este artículo Don’t Invert that matrix, del blog de John D. Cook.

REFERENCIAS

[Ref. 1] Rounding-off Errors in Matrix Processes, Quart. J. Mech. Appl. Math. 1, pp 287-308 (1948)

[Ref. 2] Turing’s Paper on Rounding-off Errors in Matrix Processes.

[Ref. 3] Blog de Alexandre Borovic. En su entrada Alan Turing and Linear Algebra.

[Ref. 4] Introduction to matrix computation. Gilbert W. Stewart.

Nota: en este enlace pueden encontrar el capítulo 3.

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s