Il existe de nombreux articles sur la décomposition LU, mais aucun article ne mentionnait la fonctionnalité de décomposition LU selon laquelle les calculs peuvent être effectués avec une économie de mémoire, j'ai donc essayé d'expliquer la décomposition LU un peu plus attentivement, en faisant à nouveau attention à l'économie de mémoire. pense.
Après avoir étudié React, j'ai créé un site qui peut facilement effectuer une décomposition LU. Cependant, pour faciliter le rendu, l'implémentation d'économie de mémoire décrite ci-dessous n'est pas implémentée. Expérience de décomposition LU
Avant de discuter de la décomposition LU, parlons des matrices triangulaires. Le but de la décomposition LU est ** de représenter les coefficients d'équations simultanées uniquement dans une matrice triangulaire **. J'expliquerai la raison. Tout d'abord, considérons une équation simultanée simple.
\left\{
\begin{array}{rrrrrrr}
2x_1 & & & & &=& 2 \\
x_1 &+& 4x_2 & & &=& 9 \\
4x_1 &-& 3x_2 &+& 3x_3 &=& -5
\end{array}
\right.
La solution de cette équation simultanée peut être facilement obtenue en la résolvant séquentiellement.
\begin{align}
x_1 &= \frac{2}{2} = 1 \\
x_2 &= \frac{9 - x_1}{4} = \frac{9-1}{4} = 2 \\
x_3 &= \frac{-5 - 4x_1 + 3x_2}{3} = \frac{-5 - 4 + 6}{3} = -1
\end{align}
Ici, si cette équation simultanée est exprimée par une matrice,
\left(
\begin{matrix}
2 & 0 & 0 \\
1 & 4 & 0 \\
4 & -3 & 3
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right) \tag{1}
On dirait. Vous pouvez voir que les coefficients des équations simultanées sont représentés par une ** matrice triangulaire inférieure ** dont les éléments au-dessus de la composante diagonale sont $ 0 $. Dans ce cas, lorsque vous trouvez $ x_i (i \ in 2, 3) $, vous connaissez déjà $ x_ {i-1} (i \ in 1,2) $, vous pouvez donc les résoudre dans l'ordre. Cette méthode est appelée ** substitution avant **. Ensuite, considérez les équations simultanées suivantes.
\left\{
\begin{array}{rrrrrrr}
3x_1 &-& 3x_2 &+& 4x_3 &=& -5 \\
& & 4x_2 &+& x_3 &=& 9 \\
& & & & 2x_3 &=& 2
\end{array}
\right. \\
\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
-5 \\
9 \\
2
\end{matrix}
\right)
Dans ce cas, les coefficients des équations simultanées sont représentés par une ** matrice triangulaire supérieure ** où les composantes sous la composante diagonale sont $ 0 $. Dans ce cas, quand on trouve $ x_i (i \ in 2, 1) $, $ x_ {i + 1} (i \ in 3, 2) $ est connu, donc l'ordre inverse dans le cas de la matrice triangulaire inférieure Peut être résolu. Cette méthode est appelée ** substitution arrière **. De cette manière, lorsque les coefficients d'équations simultanées sont représentés par une matrice triangulaire, il peut être résolu avec un algorithme simple. Ici, je vais définir les mots. L'expression $ (1) $ est exprimée comme suit.
A \boldsymbol{x} = \boldsymbol{b} \\
A =
\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right), \
\boldsymbol{x} =
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right), \
\boldsymbol{b} =
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right)
Ici, $ A $ est appelé une matrice système, $ \ boldsymbol {x} $ est appelé un vecteur de solution, et $ \ boldsymbol {b} $ est appelé un vecteur de droite. Dans ce qui suit, cette notation est généralisée, où $ A $ est la matrice $ n \ times n $ et $ \ boldsymbol {x}, \ \ boldsymbol {b} $ est le vecteur de dimension $ n $.
Comme mentionné ci-dessus, si les coefficients d'équations simultanées peuvent être représentés par une matrice triangulaire, c'est presque comme résoudre les équations. Considérons maintenant d'exprimer la matrice de coefficients $ A $ comme le produit de la matrice triangulaire inférieure $ L $ et de la matrice triangulaire supérieure $ U $.
LU \boldsymbol{x} = \boldsymbol{b} \\
Si $ U \ boldsymbol {x} = \ boldsymbol {y} $, alors $ \ boldsymbol {y} $ peut être obtenu par substitution avant, et en utilisant le résultat, $ \ boldsymbol {x} $ peut être obtenu par substitution arrière. Vous pouvez voir que c'est obligatoire.
\begin{align}
L \boldsymbol{y} &= \boldsymbol{b} \\
U \boldsymbol{x} &= \boldsymbol{y}
\end{align}
Alors, comment trouvez-vous $ L $ et $ U $? Avant d'entrer dans le sujet principal, concentrons-nous sur le degré de liberté. Il y a respectivement $ n ^ 2 $ éléments pour $ A $ et $ n (n + 1) / 2 $ éléments pour $ L et \ U $, pour un total de $ n ^ 2 + n $ éléments. .. Par conséquent, après décomposition, le degré de liberté est $ n $ plus grand, mais cela signifie que les composantes diagonales de $ L $ sont toutes $ 1 $, donc le degré de liberté devient $ n ^ 2 $, qui devient $ A $. D'un autre côté, $ L et \ U $ peuvent être déterminés de manière unique. (À propos, peu importe que la composante diagonale de $ L $ soit de 1 $ ou la composante diagonale de $ U $ soit de 1 $. Cela dépend de la personne.)
Maintenant, je vais vous expliquer comment le trouver concrètement, mais avant cela, je vais vous expliquer la petite matrice.
Considérons la matrice suivante $ M $.
M = \left(
\begin{matrix}
m_{11} & m_{12} &|& m_{13} & m_{14} & m_{15} & m_{16} \\
m_{21} & m_{22} &|& m_{23} & m_{24} & m_{25} & m_{26} \\
- & - & - & - & - & - & - \\
m_{31} & m_{32} &|& m_{33} & m_{34} & m_{35} & m_{36} \\
m_{41} & m_{42} &|& m_{43} & m_{44} & m_{45} & m_{46} \\
m_{51} & m_{52} &|& m_{53} & m_{54} & m_{55} & m_{56} \\
m_{61} & m_{62} &|& m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)
Séparez cette matrice par une ligne pointillée et donnez un nom à chaque bloc (** petite matrice **).
M_{11} =
\left(
\begin{matrix}
m_{11} & m_{12} \\
m_{21} & m_{22}
\end{matrix}
\right)
, \
M_{12} =
\left(
\begin{matrix}
m_{13} & m_{14} & m_{15} & m_{16} \\
m_{23} & m_{24} & m_{25} & m_{26}
\end{matrix}
\right) \\
M_{21} =
\left(
\begin{matrix}
m_{31} & m_{32} \\
m_{41} & m_{42} \\
m_{51} & m_{52} \\
m_{61} & m_{62}
\end{matrix}
\right)
, \
M_{22} =
\left(
\begin{matrix}
m_{33} & m_{34} & m_{35} & m_{36} \\
m_{43} & m_{44} & m_{45} & m_{46} \\
m_{53} & m_{54} & m_{55} & m_{56} \\
m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)
Alors la matrice $ M $ est exprimée comme suit:
M =
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)
Considérons maintenant la matrice $ 6 \ times 6 $ $ N $ et essayez de la séparer de la même manière que $ M $.
N =
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right)
Et le produit de $ M $ et $ N $ peut en fait être exprimé comme:
\begin{align}
MN &=
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
N_{11}N_{11} + N_{12}N_{21} & N_{11}N_{12} + N_{12}N_{22} \\
N_{21}N_{11} + N_{22}N_{21} & N_{21}N_{12} + N_{22}N_{22}
\end{matrix}
\right) \\
\end{align}
En d'autres termes, vous pouvez traiter une petite matrice comme un élément d'une matrice. (Vous pouvez le voir en le calculant.) Cependant, chaque sous-matrice doit être divisée en $ M $ et $ N $ afin qu'elle puisse être calculée de manière cohérente.
L'introduction a été allongée, mais voici les étapes pour enfin trouver les matrices $ L $ et $ U $. Commencez par diviser $ L $ et $ U $ en matrices plus petites, comme indiqué ci-dessous. Où $ O $ est une matrice avec tous les éléments $ 0 $.
L =
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right), \
U =
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right)
$ L_ {22} et U_ {22} $ sont les triangles inférieur et supérieur de $ (n-1) \ times (n-1) $, respectivement. De même, divisez la matrice de coefficients $ A $ de la même manière.
A =
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)
Puisque $ A = LU $,
\begin{align}
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)
&=
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
u_{11} & U_{12} \\
L_{21}u_{11} & L_{21}U_{12} + L_{22}U_{22}
\end{matrix}
\right)
\end{align}
Est exprimé comme. De cette relation
u_{11} = a_{11}, \ U_{12} = A_{12}
Obtenir Ici, en supposant $ a_ {11} \ neq 0 $,
L_{21} = \frac{A_{21}}{a_{11}}
Et $ L_ {21} $ sont également déterminés. Cela vous donnera la première colonne de $ L $ et la première ligne de $ U $. Eh bien, je pense que c'est probablement facile à comprendre jusqu'à présent. La question est de savoir comment gérer $ A_ {22} = L_ {21} U_ {12} + L_ {22} U_ {22} $. Ici, comme mentionné précédemment, $ L_ {21} $ et $ U_ {12} $ ont déjà été obtenus, donc le premier terme sur le côté droit a déjà été obtenu. Ensuite, concernant le deuxième terme sur le côté droit, comment le trouver Comme mentionné au début du sujet principal, $ L_ {22} $ est la matrice triangulaire inférieure et $ U_ {22} $ est la matrice triangulaire supérieure. Par conséquent, en transposant le premier terme,
\hat{A}_{22} = A_{22} - L_{21}U_{12}
étant donné que,
A^\prime_{22} = L_{22}U_{22}
Peut être exprimé comme. Cela signifie que la $ (n-1) \ times (n-1) $ matrice $ A ^ \ prime_ {22} $ est décomposée en une matrice triangulaire inférieure et une matrice triangulaire supérieure. Par conséquent, en appliquant la même procédure que ci-dessus pour ce $ A ^ \ prime_ {22} $, le problème de la décomposition LU de la matrice de $ (n-2) \ times (n-2) $ est réduit. De plus, il peut être résolu séquentiellement comme le problème de la décomposition LU de la matrice de la matrice $ (n-3) \ times (n-3) $. Une fois que vous avez compris cela, tout ce que vous avez à faire est de le mettre dans votre code.
Maintenant, j'écrirai le code pour trouver les matrices $ L $ et $ U $, mais si je ne pense à rien, le tableau bidimensionnel pour la matrice de coefficients $ A $ et le tableau bidimensionnel pour la matrice triangulaire inférieure $ L $ , Vous voudrez peut-être avoir trois variables dans le tableau à deux dimensions pour la matrice triangulaire supérieure $ U $. En fait, dans l'article qui présente la décomposition LU (par exemple, j'ai écrit la méthode de décomposition LU par programme), les variables pour chaque matrice sont définies. Faire. Cependant, bien que ce soit une méthode destructive, si vous prenez la méthode d'écrasement de $ A $, une seule variable est requise. Plus précisément, envisagez d'utiliser le fait que la composante diagonale de $ L $ est de 1 $ et de gérer $ L $ et $ U $ dans une matrice $ V $ comme suit.
L =
\left(
\begin{matrix}
1 & 0 & 0 \\
l_{11} & 1 & 0 \\
l_{21} & l_{22} & 1 \\
\end{matrix}
\right), \
U =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{22} \\
0 & 0 & u_{33} \\
\end{matrix}
\right) \\
V =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
l_{11} & u_{22} & u_{22} \\
l_{21} & l_{22} & u_{33} \\
\end{matrix}
\right)
En considérant un tel $ V $ et en écrasant $ V $ sur $ A $, la mémoire requise pour sauvegarder le tableau n'est que de $ A $, et la mémoire peut être sauvegardée.
L'introduction est devenue assez volumineuse, mais je présenterai un exemple d'implémentation en Python. Comme je l'ai expliqué en détail, la mise en œuvre elle-même est étonnamment facile.
import numpy as np
def main():
n = int(input())
A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])
for i in range(n):
for j in range(i + 1, n):
A[j][i] /= A[i][i]
for k in range(i + 1, n):
A[j][k] -= A[j][i] * A[i][k]
return A
if __name__ == "__main__":
print(main())
Un commentaire sera ajouté. Je ne mettrai que le code.
import numpy as np
def main():
n = int(input())
A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])
P = np.identity(n)
for i in range(n):
maxidx = np.argmax(abs(A[i:n, i])) + i
A[i, :], A[maxidx, :] = A[maxidx, :], A[i, :].copy()
P[i, :], P[maxidx, :] = P[maxidx, :], P[i, :].copy()
for j in range(i + 1, n):
A[j][i] /= A[i][i]
for k in range(i + 1, n):
A[j][k] -= A[j][i] * A[i][k]
return A, P
if __name__ == "__main__":
A, P = main()
En fait, le calcul peut ne pas être stable tel quel. C'est parce que nous supposons $ a_ {ii} \ neq 0 $. Pour contourner ce problème, vous devez effectuer une ** sélection pivot partielle ** qui permute les lignes. J'ajouterai cette explication quand j'aurai le temps.
Calcul numérique pour l'ingénierie
Si vous avez des questions, n'hésitez pas à demander! Je ne connais pas très bien les mathématiques, donc je ne pourrai peut-être pas m'en occuper.
Recommended Posts