Descomposició LU

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

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]

${\displaystyle 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:

${\displaystyle {\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 ${\displaystyle a_{11}=l_{11}u_{11}}$. Si ${\displaystyle a_{11}=0}$, llavors o bé ${\displaystyle l_{11}=0}$ o bé ${\displaystyle 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 esmenar 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:

${\displaystyle 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:

${\displaystyle 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

${\displaystyle 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

Descomponem la següent matriu 2×2:

${\displaystyle {\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:

${\displaystyle \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ó:

${\displaystyle \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

${\displaystyle {\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

Matrius quadrades

Tota matriu quadrada ${\displaystyle A}$ admet una descomposició LUP.[4] Si ${\displaystyle 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 ${\displaystyle A}$ és una matriu singular de rang ${\displaystyle k}$, llavors adment una descomposició LU si els ${\displaystyle 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 que totes les entrades de la diagonal de ${\displaystyle L}$ (o de la de ${\displaystyle U}$) valguin 1.

Matrius simètriques definides positives

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

${\displaystyle A=LL^{*}.\,}$

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

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

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 2n^3 / 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

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, ${\displaystyle D_{1}=A_{1,1}}$ i per ${\displaystyle i=2,\ldots ,n}$, ${\displaystyle D_{i}}$ és el quocient de la i-sima submatriu principal entre la (i-1)-sima submatriu principal.

Algorisme de Doolittle

Donada una matriu N × N

${\displaystyle A=(a_{n,n})}$

definim

${\displaystyle 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

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

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

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

Definim

${\displaystyle 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ó

${\displaystyle 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 ${\displaystyle L=L_{1}^{-1}\ldots L_{N-1}^{-1}}$. Com que la inversa d'una matriu triangular inferior L_n é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

${\displaystyle 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 ${\displaystyle A=LU}$.

És obvi que, per tal que aquest algorisme funcioni, hom necessita que els elements ${\displaystyle a_{n,n}^{(n-1)}}$ siguin diferents de 0 en cada pas (vegeu la definició de ${\displaystyle 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 ${\displaystyle P^{-1}A=LU}$.

Algorismes de Crout i LUP

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 ${\displaystyle A}$ té una entrada no-nul·la en la primera fila, llavors prenem una matriu permutació ${\displaystyle P_{1}}$ tal que ${\displaystyle AP_{1}}$ tingui una entrada no-nul·la com a element superior esquerre. Altrament, agafem ${\displaystyle P_{1}}$ que sigui la matriu identitat. Sigui ${\displaystyle A_{1}=AP_{1}}$.
2. Sigui ${\displaystyle A_{2}}$ la matriu que s'obté a partir de ${\displaystyle A_{1}}$ tot eliminant-ne la primera fila i la primera columna. Descomponem ${\displaystyle A_{2}=L_{2}U_{2}P_{2}}$ de forma recursiva. Construïm ${\displaystyle L}$ a partir de ${\displaystyle L_{2}}$ afegint-hi una fila amb zeros a la part superior, i després afegint-hi la primera columna de ${\displaystyle A_{1}}$ a l'esquerra.
3. Construïm ${\displaystyle U_{3}}$ a partir de ${\displaystyle 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 ${\displaystyle P_{3}}$ a partir de ${\displaystyle P_{2}}$ de forma semblant, i definim ${\displaystyle A_{3}=A_{1}/P_{3}=AP_{1}/P_{3}}$. Sigui ${\displaystyle P}$ la inversa de ${\displaystyle P_{1}/P_{3}}$.
4. En aquest punt, ${\displaystyle A_{3}}$ és igual a ${\displaystyle LU_{3}}$, excepte (possiblement) la primera fila. Si la primera fila de ${\displaystyle A}$ és zero, llavors ${\displaystyle A_{3}=LU_{3}}$, ja que totes dues tenen la primera fila a zeros, i llavors es compleix que ${\displaystyle A=LU_{3}P}$, com desitjàvem. Altrament, ${\displaystyle A_{3}}$ i ${\displaystyle LU_{3}}$ tenen la mateixa entrada no-nul·la a l'extrem superior esquerre, i ${\displaystyle A_{3}=LU_{3}U_{1}}$ per alguna matriu quadrada triangular superior ${\displaystyle U_{1}}$ amb valors 1 a la diagonal (${\displaystyle U_{1}}$ «neteja» entrades de ${\displaystyle LU_{3}}$ i afegeix entrades de ${\displaystyle A_{3}}$ a través de l'extrem superior esquerre). En aquest moment, ${\displaystyle A=LU_{3}U_{1}P}$ és una descomposició de la forma desitjada.

Complexitat

Si dues matrius d'ordre n es poden multiplicar en temps M(n), on M(n)≥n^a 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(n^2,376) basat en l'algorisme de Coppersmith–Winograd.

Decomposició per matrius disperses

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.

Reducció LU

Vegeu també: Reducció LU

La reducció LU és un algorisme relacionat amb la descomposició LU. Aquest terme s'acostuma a utilitzar en el context de supercomputació i més específicament en computació paral·lela. En aquest context, la reducció LU s'usa com un algorisme de benchmarking, és a dir, proporciona una mesura comparativa de la velocitat per diferents computadors. La reducció LU és una versió especial paral·lelitzada d'un algorisme per la descomposició LU (vegeu-ne un exemple[13]). La versió paral·lelitzada distribueix la càrrega de treball sobre una fila de la matriu cap a un sol processador, i després sincronitza el resultat amb la matriu sencera.[14][15]

Aplicacions

Resolució d'equacions lineals

Donat un sistema d'equacions lineals en forma matricial

${\displaystyle Ax=b,\,}$

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

${\displaystyle LUx=Pb.\,}$

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

1. Primer, resolem l'equació ${\displaystyle Ly=Pb}$ per y.
2. Segon, resolem l'equació ${\displaystyle Ux=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 ${\displaystyle {\frac {2}{3}}n^{3}}$ operacions de coma flotant si la matriu ${\displaystyle A}$ té grandrària ${\displaystyle n}$. Això ho fa dues vegades més ràpid que els algorismes basats en una descomposició QR, que té un cost aproximat de ${\displaystyle {\frac {4}{3}}n^{3}}$ operacions de coma flotant, si usem transformacions de Householder. Per aquest motiu, hom prefereix usar la descomposició LU.[16]

Càlcul de la inversa d'una matriu

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

${\displaystyle AX=LUX=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.[17]

Càlcul del determinant

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

${\displaystyle \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

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

1. MATLAB function reference. «lu». Arxivat de l'original el 13 de maig 2014. [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. Okunev & Johnson (1997), Corol·lari 3
5. Trefethen & Bau (1997), p. 166
6. Trefethen & Bau (1997), p. 161
7. 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. J. Guitart, X. Martorell, J. Torres, E. Ayguadé: Improving Java Multithreading Facilities: the Java Nanos Environment, Research Report UPC-DAC-2001-8, Departament d'Arquitectura de Computadors, Universitat Politècnica de Catalunya, març 2001, [1]^[Enllaç no actiu]
14. (eds.), José M.L.M. Palma, Jack Dongarra, Vicente Hernández. Vector and parallel processing - VECPAR 2000 4th International Conference, Porto, Portugal, June 21-23, 2000 : selected papers and invited talks. [Elektronische Ressource]. Berlin [etc.]: Springer, 2001, p. 128-141. ISBN 978-3-540-41999-0 [Consulta: 1r gener 2019].  Arxivat 2006-06-21 a Wayback Machine.
15. Oliver, José; Guitart, Jordi; Ayguadé, Eduard; Navarro, Nacho; Torres, Jordi «Strategies for the efficient exploitation of loop-level parallelism in Java». Concurrency and Computation: Practice and Experience, 13, 8-9, 01-07-2001, pàg. 663–680. Arxivat de l'original el 21 d’octubre 2007. DOI: 10.1002/cpe.573 [Consulta: 9 agost 2013]. Arxivat 21 October 2007^[Date mismatch] a Wayback Machine.
16. Trefethen & Bau (1997), p. 152
17. Golub & Van Loan (1996), p. 121

Bibliografia

• Bunch, James R.; Hopcroft, John «Triangular factorization and inversion by fast matrix multiplication» (en anglès). Mathematics of Computation, 28, 1974, pàg. 231-236. ISSN: 0025-5718.
• al.], Thomas H. Cormen ... [et. Introduction to algorithms (en anglès). 2a ed.. Cambridge, Mass.: MIT Press, 2001. ISBN 978-0-262-03293-3. 
• Golub, Gene H.; Van Loan, Charles F. Matrix computations (en anglès). 3a ed.. Baltimore, Md. [u.a.]: Johns Hopkins Univ. Press, 1996. ISBN 978-0-8018-5414-9. 
• Johnson, Roger A. Horn; Charles R.. «Secció 3.5». A: Matrix analysis (en anglès). 19a ed.. Cambridge [u.a.]: Cambridge Univ. Press, 2005. ISBN 0-521-38632-2. 
• Householder, Alston S. «The Theory of Matrices in Numerical Analysis» (en anglès). Mathematical Reviews. Dover Publications [Nova York], 1975.
• Okunev, Pavel; Johnson, Charles R. «Necessary and Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix» (en anglès). University of Connecticut, 1997.
• Poole, David. Linear algebra : a modern introduction (en anglès). 2a edició. Australia: Brooks/Cole-Thomson Learning, 2006. ISBN 0-534-99845-3. 
• al.], William H. Press ... [et; Teukolsky, S.A.; Vetterling, W.T.; Finnery, B.P.. «Secció 2.3». A: Numerical recipes : the art of scientific computing (en anglès). 3a edició. Cambridge, UK: Cambridge University Press, 2007. ISBN 978-0-521-88068-8. ^[Enllaç no actiu]
• III, / Lloyd N. Trefethen, David Bau. Numerical linear algebra. (en anglès). 3a ed.. Philadelphia: Society for Industrial and Applied Mathematics, 1996. ISBN 978-0-89871-361-9. 

Vegeu també

• Factorització de Cholesky
• Descomposició de matrius
• Descomposició QR

Enllaços externs

• LU decomposition a MathWorld.
• LU decomposition a Math-Linux.
• LU matrix factorization. MATLAB reference.
• LU Decomposition by Ed Pegg, Jr., The Wolfram Demonstrations Project, 2007.