Décomposition LU en Python

introduction

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

Matrice triangulaire

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 $.

Décomposition LU

Aperçu

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.)

Comment trouver

Maintenant, je vais vous expliquer comment le trouver concrètement, mais avant cela, je vais vous expliquer la petite matrice.

Petite procession

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.

Comment trouver le sujet principal

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.

la mise en oeuvre

Implémentation d'économie de mémoire

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.

Sujet de mise en œuvre

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())

Sélection de pivot partiel

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()

Supplément

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.

référence

Calcul numérique pour l'ingénierie

finalement

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

Décomposition LU en Python
Introduction à l'algèbre linéaire avec Python: Décomposition A = LU
Quadtree en Python --2
CURL en Python
Métaprogrammation avec Python
Python 3.3 avec Anaconda
Géocodage en python
SendKeys en Python
Méta-analyse en Python
Unittest en Python
Époque en Python
Discord en Python
Allemand en Python
DCI en Python
tri rapide en python
nCr en python
N-Gram en Python
Programmation avec Python
Plink en Python
Constante en Python
FizzBuzz en Python
Sqlite en Python
Étape AIC en Python
LINE-Bot [0] en Python
CSV en Python
Assemblage inversé avec Python
Réflexion en Python
Constante en Python
nCr en Python.
format en python
Scons en Python 3
Puyopuyo en python
python dans virtualenv
PPAP en Python
Quad-tree en Python
Réflexion en Python
Chimie avec Python
Hashable en Python
DirectLiNGAM en Python
LiNGAM en Python
Aplatir en Python
Aplatir en python
Liste triée en Python
AtCoder # 36 quotidien avec Python
Texte de cluster en Python
AtCoder # 2 tous les jours avec Python
Daily AtCoder # 32 en Python
Daily AtCoder # 6 en Python
Modifier les polices en Python
Motif singleton en Python
Opérations sur les fichiers en Python
Lire DXF avec python
Daily AtCoder # 53 en Python
Séquence de touches en Python
Utilisez config.ini avec Python
Daily AtCoder # 33 en Python
Résoudre ABC168D en Python
Distribution logistique en Python