Descomposició LU

De Viquipèdia
Dreceres ràpides: navegació, cerca

En àlgebra lineal la descomposició LU (també anomenada factorització LU o LR[1][2]) descompon una matriu com a producte d'una matriu triangular inferior i una matriu triangular superior. Sovint, aquest producte inclou també una matriu permutació. La descomposició LU es pot interpretar com una forma matricial del mètode de reducció de Gauss. En computació, és usual emprar la descomposició LU per resoldre sistemes d'equacions lineals; també és un procés clau en el càlcul de la inversa d'una matriu, i en el càlcul del determinant. El mètode de descomposició LU fou desenvolupat pel matemàtic Alan Turing.[3]

Definicions[modifica | modifica el codi]

Descomposició LDU d'una matriu de Walsh

Sigui A una matriu quadrada. Una descomposició LU és una factorització de A, amb reordenacions o permutacions de les seves files i/o columnes, en dos factors, una matriu triangular inferior L (de l'anglès lower, inferior) i una matriu triangular superior U (de l'anglès upper, superior),[nota 1]

 A = LU, \,

En la matriu triangular inferior, tots els elements per sobre de la diagonal són 0; anàlogament, en la matriu triangular superior, tots els elements per sota de la diagonal són 0. Per exemple, per una matriu A 3×3, la seva descomposició LU té la forma:


        \begin{bmatrix}
           a_{11} & a_{12} & a_{13} \\
           a_{21} & a_{22} & a_{23} \\
           a_{31} & a_{32} & a_{33} \\
        \end{bmatrix} =
      \begin{bmatrix}
           l_{11} & 0 & 0 \\
           l_{21} & l_{22} & 0 \\
           l_{31} & l_{32} & l_{33} \\
        \end{bmatrix}
        \begin{bmatrix}
           u_{11} & u_{12} & u_{13} \\
           0 & u_{22} & u_{23} \\
           0 & 0 & u_{33} \\
        \end{bmatrix}.

Aquesta factorització pot no ser possible si no es realitzen les ordenacions o permutacions adequades sobre la matriu original. Per exemple, és fàcil comprovar (calculant explícitament el producte matricial, entrada per entrada) que a_{11} = l_{11} u_{11}. Si a_{11} = 0, llavors o bé l_{11}=0 o bé u_{11}=0 (o tots dos), la qual cosa implicaria que o bé L o bé U és singular. Això no és possible si A no és singular. Hom pot subsanar aquest problema mitjançant reordenant les files de A de tal manera que el primer element de la matriu permutada sigui no-nul. El mateix problema pot sorgir en posteriors passos del procediment, que es pot resoldre de la mateixa forma; vegeu el procediment bàsic més endavant.

Hom pot demostrar que una permutació adequada de les files (o de les columnes) és suficient per la descomposició LU. L'anomenada descomposició LU amb pivot parcial fa referència a la descomposició LU que només té permutacions de files:

 PA = LU, \,

on L i U són les matrius triangulars inferior i superior respectivament, i P és una matriu permutació que, quan multiplica per l'esquerra a A, reordena les files de A. Hom pot veure que qualsevol matriu quadrada es pot descompondre d'aquesta manera,[4] i la descomposició és numèricament estable.[5] Això fa que la descomposició LUP sigui una tècnica útil a la pràctica.

En una descomposició LU amb pivot complet intervenen permutacions de files i de columnes:

 PAQ = LU, \,

on L, U i P són com abans, i Q és una matriu permutació que reordena les columnes de A.[6]

Una descomposició LDU és una descomposició de la forma

 A = LDU, \,

on D és una matriu diagonal, i L i U són matrius triangulars unitàries, la qual cosa vol dir que totes les entrades de les diagonals de L i de U valen 1.

Abans hem suposat que A és una matriu quadrada, però aquestes descomposicions es poden generalitzar per matrius rectangulars. En aquest cas L i P són matrius quadrades amb el mateix nombre de files que A, mentre que U té exactament el mateix nombre de files i de columnes que A. Aquí, s'ha d'interpretar triangular superior com que té entrades 0 per sota de la diagonal principal, que comença a l'entrada superior esquerra.

Exemple[modifica | modifica el codi]

Descomponem la següent matriu 2×2:


        \begin{bmatrix}
           4 & 3 \\
           6 & 3 \\
        \end{bmatrix} =
      \begin{bmatrix}
           l_{11} & 0 \\
           l_{21} & l_{22} \\
        \end{bmatrix}
        \begin{bmatrix}
           u_{11} & u_{12} \\
           0 & u_{22} \\
        \end{bmatrix}.

Una manera de trobar la descomposició LU d'aquesta matriu senzilla seria resoldre el sistema d'equacions lineals. Si multipliquem explícitament L per U, obtenim:

\left\{
\begin{array}{lcll}
l_{11} \cdot u_{11} & + & 0 \cdot 0          & = 4 \\
l_{11} \cdot u_{12} & + & 0 \cdot u_{22}     & = 3 \\
l_{21} \cdot u_{11} & + & l_{22} \cdot 0      & = 6 \\
l_{21} \cdot u_{12} & + & l_{22} \cdot u_{22} & = 3
\end{array}
\right.{.}

Aquest sistema d'equacions és indeterminat. En aquest cas, dos elements no-nuls qualssevol de les matrius L i U poden actuar com a paràmetres de la solució, i poden prendre qualsevol valor no-nul. Per tant, per trobar la descomposició LU única, és necessari fixar alguna restricció sobre les matrius L i U. Per exemple, podem requerir que la matriu triangular inferior L sigui unitària (és a dir, que totes les entrades de la seva diagonal principal valguin 1). Aleshores, el sistema d'equacions té la següent solució:

\left\{
\begin{array}{lcl}
l_{21} & = & 1,5 \\
u_{11} & = & 4 \\
u_{12} & = & 3 \\
u_{22} & = & -1,5 \\
\end{array}
\right.{.}

Si substituïm aquests valors en la descomposició LU anterior, obtenim


        \begin{bmatrix}
           4 & 3 \\
           6 & 3 \\
        \end{bmatrix} =
      \begin{bmatrix}
           1 & 0 \\
           1,5 & 1 \\
        \end{bmatrix}
        \begin{bmatrix}
           4 & 3 \\
           0 & -1,5 \\
        \end{bmatrix}.

Existència i unicitat[modifica | modifica el codi]

Matrius quadrades[modifica | modifica el codi]

Tota matriu quadrada  A admet una descomposició LUP.[4] Si  A és invertible, llavors admet una descomposició LU (o LDU) si i només si tots els seus menors principals són no-nuls.[7] Si  A és una matriu singular de rang  k , llavors adment una descomposició LU si els  k primers menors principals són no-nuls (el recíproc no és cert).[8]

Si una matriu quadrada invertible té una descomposició LDU, llavors és única.[7] En aquest cas, la descomposició LU també és única si afegim la hipòtesi de què totes les entrades de la diagonal de  L (o de la de  U ) valguin 1.

Matrius simètriques definides positives[modifica | modifica el codi]

Si A és una matriu simètrica (o hermítica, si A és complexa) definida positiva, podem aconseguir que U sigui la transposada conjugada de L. És a dir, podem escriure A com

 A = L L^{*}. \,

Aquesta descomposició s'anomena descomposició de Cholesky. La descomposició de Cholesky sempre existeix i és única. Addicionalment, el càlcul de la descomposició de Cholesky és més eficient i numèricament més estable que el càlcul d'altres descomposicions LU.

Cas general[modifica | modifica el codi]

Per una matriu (no necessàriament invertible) sobre un cos qualsevol, hom coneix de manera precisa les condicions necessàries i suficients per la seva factorització LU. Aquestes condicions s'expressen en termes dels rangs de certes submatrius. L'algorisme de reducció de Gauss per obtenir la descomposició LU es pot estendre a aquest cas més general.[9]

Algorismes[modifica | modifica el codi]

La descomposició LU és bàsicament una forma modificada del mètode de reducció de Gauss. Transformem la matriu A en una matriu triangular superior U, mitjançant l'eliminació de les entrades per sota de la diagonal principal. L'algorisme de Doolittle realitza aquesta eliminació columna per columna començant des de l'esquerra, multiplicant A per l'esquerra per matrius triangulars inferiors atòmiques. Això proporciona una matriu triangular inferior unitària i una matriu triangular superior. L'algorisme de Crout és lleugerament diferent, i construeix una matriu triangular inferior i una matriu triangular superior unitària.

El càlcul de la descomposició LU mitjançant qualsevol d'aquests algorismes necessita 2n3 / 3 operacions en coma flotant, si ignorem els termes d'ordre inferior. Si fem servir el mètode del pivot parcial, només s'hi afegeix un terme quadràtic; però això no és cert pel mètode de pivot complet.[10]

Fórmula tancada[modifica | modifica el codi]

Quan existeix una factorització LDU única, existeix una fórmula tancada (explícita) pels elements de L, D i U en termes dels quocients dels determinants de certes submatrius de la matriu original A.[11] En particular, D_1 = A_{1,1} i per i = 2, \ldots, n, D_i és el quocient de la i-sima submatriu principal entre la (i-1)-sima submatriu principal.

Algorisme de Doolittle[modifica | modifica el codi]

Donada una matriu N × N


A= (a_{n,n})

definim

 A^{(0)} = A

i llavors iterem n = 1,...,N-1 de la següent manera.

Eliminem els elements per sota de la diagonal principal a la n-sima columna de A(n-1) tot sumant a la i-sima fila d'aquesta matriu la n-sima columna multiplicada per

l_{i,n} = -\frac{a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}

per i = n+1,\ldots,N. Això es pot realitzar multiplicant A(n-1) per l'esquerra per la matriu triangular inferior


L_n =
\begin{pmatrix}
     1 &        &           &         &         & 0 \\
       & \ddots &           &         &         &   \\
       &        &         1 &         &         &   \\
       &        & l_{n+1,n} &  \ddots &         &   \\
       &        &    \vdots &         &  \ddots &   \\
     0 &        &   l_{N,n} &         &         & 1 \\
\end{pmatrix}.

Definim

 A^{(n)} = L_n A^{(n-1)}.

Després de N-1 passos, hem eliminat tots els elements per sota de la diagonal principal, de tal manera que tenim una matriu triangular superior A(N-1). Hem trobat, doncs, la descomposició


A = L_{1}^{-1} L_{1} A^{(0)}
= L_{1}^{-1} A^{(1)} = L_{1}^{-1} L_{2}^{-1} L_{2} A^{(1)} = 
L_{1}^{-1}L_{2}^{-1} A^{(2)} =\ldots = L_{1}^{-1} \ldots L_{N-1}^{-1} A^{(N-1)}.

Denotem la matriu triangular superior A(N-1) per U, i L=L_{1}^{-1} \ldots L_{N-1}^{-1}. Com que la inversa d'una matriu triangular inferior Ln és també una matriu triangular inferior, i la multiplicació de dues matrius triangulars inferiors també és una matriu triangular inferior, tenim que L és una matriu triangular inferior. És mes, podem veure que


L =
\begin{pmatrix}
       1 &        &            &        &            & 0 \\
-l_{2,1} & \ddots &            &        &            &   \\
         &        &          1 &        &            &   \\
  \vdots &        & -l_{n+1,n} & \ddots &            &   \\
         &        &     \vdots &        &       1    &   \\
-l_{N,1} &        & -l_{N,n}   &        & -l_{N,N-1} & 1 \\
\end{pmatrix}.

Així obtenim A=LU.

És obvi que, per tal que aquest algorisme funcioni, hom necessita que els elements a_{n,n}^{(n-1)} siguin diferents de 0 en cada pas (vegeu la definició de l_{i,n}). Si això no es compleix en algun moment, es necessita intercanviar la n-sima fila per una altra fila que estigui per sota, abans de continuar. Aquest és el motiu pel qual la descomposició LU, en general, té la forma P^{-1}A = L U .

Algorismes de Crout i LUP[modifica | modifica el codi]

L'algorisme de Cormen et al. per la descomposició LUP és una generalització de la descomposició matricial de Crout. Es pot descriure de la següent manera:

  1. Si A té una entrada no-nul·la en la primera fila, llavors prenem una matriu permutació P_1 tal que A P_1 tingui una entrada no-nul·la com a element superior esquerre. Altrament, agafem P_1 que sigui la matriu identitat. Sigui A_1 = A P_1.
  2. Sigui A_2 la matriu que s'obté a partir de A_1 tot eliminant-ne la primera fila i la primera columna. Descomponem A_2 = L_2 U_2 P_2 de forma recursiva. Construïm L a partir de L_2 afegint-hi una fila amb zeros a la part superior, i després afegint-hi la primera columna de A_1 a l'esquerra.
  3. Construïm U_3 a partir de U_2 afegint-hi una fila amb zeros a la part superior i una columna amb zeros a l'esquerra, i llavors substituint l'element superior esquerre (que de moment val 0) per 1. Construïm P_3 a partir de P_2 de forma semblant, i definim A_3 = A_1 / P_3 = A P_1 / P_3. Sigui P la inversa de P_1 / P_3.
  4. En aquest punt, A_3 és igual a L U_3, excepte (possiblement) la primera fila. Si la primera fila de A és zero, llavors A_3 =  L U_3, ja que totes dues tenen la primera fila a zeros, i llavors es compleix que A = L U_3 P, com desitjàvem. Altrament, A_3 i L U_3 tenen la mateixa entrada no-nul·la a l'extrem superior esquerre, i A_3 = L U_3 U_1 per alguna matriu quadrada triangular superior U_1 amb valors 1 a la diagonal (U_1 «neteja» entrades de L U_3 i afegeix entrades de A_3 a través de l'extrem superior esquerre). En aquest moment, A = L U_3 U_1 P és una descomposició de la forma desitjada.

Complexitat[modifica | modifica el codi]

Si dues matrius d'ordre n es poden multiplicar en temps M(n), on M(n)≥na per algun a>2, llavors la descomposició LU es pot calcular en temps O(M(n)).[12] Això significa, per exemple, que existeix un algorisme d'ordre O(n2,376) basat en l'algorisme de Coppersmith–Winograd.

Decomposició per matrius disperses[modifica | modifica el codi]

Existeixen algorismes específics per factoritzar matrius disperses grans. Aquests algorismes intenten trobar factors dispersos L i U. En general, el cost de càlcul està determinat pel nombre d'entrades no-nul·les, en comptes de per la grandària de la matriu.

Aquests algorismes usen l'intercanvi de files i columnes per evitar l'aparició d'entrades que canvien de 0 a un valor diferent durant l'execució de l'algorisme. Mitjançant la teoria de grafs, hom pot analitzar aquestes reordenacions de files i columnes.

Aplicacions[modifica | modifica el codi]

Resolució d'equacions lineals[modifica | modifica el codi]

Donat un sistema d'equacions lineals en forma matricial

A x = b, \,

volem resoldre l'equació per x, si coneixem A i b. Suposem que hem obtingut la descomposició LUP de A tal que PA = LU, (o la descomposició LU si existeix, i llavors P = I); aleshores podem reescriure l'equació com

L U x = P b. \,

En aquest cas, hom troba la solució en dos passos lògics:

  1. Primer, resolem l'equació  L y = P b per y.
  2. Segon, resolem l'equació  U x = y per x.

Notem que, en tots dos passos, estem treballant amb matrius triangulars (L i U) que es poden resoldre directament per substitució enrere i endavant, sense necessitat de fer servir el mètode de reducció de Gauss.

Hom pot emprar el procediment anterior de forma repetida per resoldre l'equació diverses vegades per diferents b. En aquest cas, és més ràpid (i més convenient) realitzar una descomposició LU de la matriu A una sola vegada, i llavors resoldre les matrius triangulars pels diferents b, en comptes d'usar la reducció de Gauss cada cop. Es pot interpretar que les matrius L i U han «codificat» el procés de reducció de Gauss.

El cost de resoldre un sistema d'equacions lineals és aproximadament  \frac{2}{3} n^3 operacions de coma flotant si la matriu  A té grandrària  n . Això ho fa dues vegades més ràpid que els algorismes basats en una descomposició QR, que té un cost aproximat de \frac{4}{3} n^3 operacions de coma flotant, si usem transformacions de Householder. Per aquest motiu, hom prefereix usar la descomposició LU.[13]

Càlcul de la inversa d'una matriu[modifica | modifica el codi]

En la resolució de sistemes d'equacions lineals, normalment tenim un vector b de longitud igual a l'alçada de la matriu A. En comptes d'un vector b, podem tenir una matriu B de dimensió n per p, amb la qual cosa estem intentant trobar una matriu X (també de dimensió n per p) tal que


A X = L U X = B.

Podem usar el mateix algorisme que hem vist per resoldre cada columna de la matriu X. Si ara suposem que B és la matriu identitat de dimensió n, tenim que la matriu X resultant és la inversa de A.[14]

Càlcul del determinant[modifica | modifica el codi]

Donada la descomposició LUP A = P^{-1} L U d'una matriu quadrada A, es pot calcular directament el determinant de A de la següent manera:

\det(A) = \det(P^{-1}) \det(L) \det(U) = (-1)^S \left( \prod_{i=1}^n l_{ii} \right)  \left( \prod_{i=1}^n u_{ii} \right).

La segona equació és una conseqüència del fet que el determinan d'una matriu triangular és simplement el producte dels elements de la seva diagonal, i que el determinant d'una matriu permutació és igual a (−1)S, on S és el nombre d'intercanvis de files de la descomposició.

El mateix procediment és vàlid per la descomposició LU; només cal considerar P com la matriu identitat.

Notes[modifica | modifica el codi]

  1. La nomenclatura LR també prové de l'anglès, en aquest cas de left, esquerra, i right, dreta, car una matriu triangular inferior L té els seus elements no nuls a l'esquerra de la diagonal principal (inclosa); anàlogament, una matriu triangular superior U o R té els seus elements no nuls a la dreta de la diagonal principal (inclosa).

Referències[modifica | modifica el codi]

  1. MATLAB function reference. «lu». [Consulta: 8 setembre 2013].
  2. Dasgupta, Bhaskar. Applied mathematical methods. India: Dorling Kindersley (India) Pvt. Ltd., licensees of Pearson Education in South Asia, 2006. ISBN 978-81-317-0068-6 [Consulta: 8 setembre 2013]. 
  3. Poole (2006)
  4. 4,0 4,1 Okunev & Johnson (1997), Corol·lari 3
  5. Trefethen & Bau (1997), p. 166
  6. Trefethen & Bau (1997), p. 161
  7. 7,0 7,1 Horn & Johnson (1985), Corol·lari 3.5.5
  8. Horn & Johnson (1985), Teorema 3.5.2
  9. Okunev & Johnson (1997)
  10. Golub & Van Loan (1996), p. 112, 119
  11. Householder (1975)
  12. Bunch & Hopcroft (1974)
  13. Trefethen & Bau (1997), p. 152
  14. Golub & Van Loan (1996), p. 121

Bibliografia[modifica | modifica el codi]

Vegeu també[modifica | modifica el codi]

Enllaços externs[modifica | modifica el codi]

Referències

Codis font en llenguatges de programació

Recursos en línia