Récemment, j'ai lu un livre intitulé "Structure et topologie des protéines (Hiraoka)" avec un ami. Dans ce livre, il y avait une méthode de décomposition d'une matrice d'entiers appelée standardisation Smith, donc je l'ai implémentée avec Numpy.
Avant d'expliquer la standardisation de Smith, il est nécessaire d'expliquer la matrice de base de la matrice entière.
La matrice de base dans l'espace vectoriel est la matrice pour effectuer la transformation de base de la matrice. (Il devrait être appris au cours de la première année d'université.). Il y a les trois transformations de base suivantes de la matrice.
Chaque transformation peut être réalisée à partir de la gauche ou de la droite en appliquant une matrice ** régulière ** (c'est-à-dire une matrice dans laquelle une matrice inverse existe). Cette matrice régulière est appelée la matrice de base.
Ici, si les matrices de base correspondant à 1, 2 et 3 sont exprimées comme $ P_i (c), Q_ {ij}, R_ {ij} (c) $, elles sont comme suit.
P_i(c) =\begin{pmatrix}
1&&&&&&\\
&\ddots&&&&&\\
&&1&&&&\\
&&&c&&&\\
&&&&1&&\\
&&&&&\ddots&\\
&&&&&&1\\
\end{pmatrix}
Q_{ij} =\begin{pmatrix}
1&&&&&&\\
&\ddots&&&&&\\
&&0&&1&&\\
&&&\ddots&&&\\
&&1&&0&&\\
&&&&&\ddots&\\
&&&&&&1\\
\end{pmatrix}
R_{ij}(c) =\begin{pmatrix}
1&&&&&&\\
&\ddots&&&&&\\
&&1&&c&&\\
&&&\ddots&&&\\
&&&&1&&\\
&&&&&\ddots&\\
&&&&&&1\\
\end{pmatrix}
De plus, l'inverse de chaque matrice de base est
Maintenant, considérons la matrice de base pour la matrice entière.
La matrice de base $ Q_ {ij}, R_ {ij} (c) $ sur l'espace vectoriel peut être adoptée comme matrice de base pour la matrice entière car la matrice inverse est la matrice entière pour l'entier $ c $. D'autre part, $ P_i (c) $ n'a pas de matrice inverse à moins que la constante $ c $ soit 1 ou -1. Puisque $ P_i (1) $ est une matrice unitaire, $ P_i (-1) $ est une matrice de base d'une matrice entière.
Par conséquent, pour l'entier $ c $, $ P_i (-1), Q_ {ij}, R_ {ij} (c) $ sont la matrice de base de la matrice entière.
La standardisation de Smith est une méthode de diagonalisation pour les matrices entières. La matrice d'entiers arbitraires A utilise les matrices $ P et Q $
B = Q^{-1} A P =
\left(
\begin{array}{c|c}
\begin{matrix}
c_1 & &0 \\
& \ddots &\\
0 & & c_k
\end{matrix}
& 0 \\
\hline
0 &0
\end{array}\right)
Peut être diagonalisé avec. Ici, $ P et Q $ sont les produits des matrices de base, et $ c_ {i + 1} $ peut être divisé par $ c_i $. Puisque les matrices $ P et Q $ sont les produits des matrices de base, il existe une matrice inverse de matrices entières. La forme diagonalisée par le produit de cette matrice de base est appelée la forme standard de Smith.
Exemple:
\begin{pmatrix}
1&0&0\\
7&-7&-4\\
0&2&1
\end{pmatrix}^{-1}
\begin{pmatrix}
7&3&2&1\\
7&6&7&7\\
4&8&2&0
\end{pmatrix}
\begin{pmatrix}
0&0&-6&13\\
0&0&-13&28\\
0&1&65&-138\\
1&-2&-49&101
\end{pmatrix}
=
\begin{pmatrix}
1&0&0&0\\
0&1&0&0\\
0&0&2&0
\end{pmatrix}
Le formulaire standard Smith est obtenu en répétant les quatre opérations suivantes.
Déplacez l'élément $ A [i, j] $ avec la plus petite valeur absolue non nulle vers le composant (1,1). En d'autres termes
{\rm moveMN}(A_t) =
\begin{cases}
Q_{i,0} A_tQ_{0,j},& A[i,j]>0\\
Q_{i,0} A_tQ_{0,j}P_1(-1),& A[i,j]<0
\end{cases}
Est défini comme.
Exemple:
\begin{pmatrix}
7&9&5&8\\
5&1&0&2\\
8&1&8&5
\end{pmatrix}
\rightarrow
\begin{pmatrix}
1&5&0&2\\
9&7&5&8\\
1&8&8&5
\end{pmatrix}
Pour tout $ i $, si $ q_i = A [i, 0] \ div A [0,0] $ (le reste est tronqué)
Exemple:
\begin{pmatrix}
1&5&0&2\\
9&7&5&8\\
1&8&8&5
\end{pmatrix}
\rightarrow
\begin{pmatrix}
1&5&0&2\\
0&-38&5&10\\
0&3&8&3
\end{pmatrix}
Pour tout $ i $, si $ q_i = A [0, i] \ div A [0,0] $ (le reste est tronqué)
Exemple:
\begin{pmatrix}
1&5&0&2\\
0&-38&5&10\\
0&3&8&3
\end{pmatrix}
\rightarrow
\begin{pmatrix}
1&0&0&0\\
0&-38&5&10\\
0&3&8&3
\end{pmatrix}
Soit $ A [i, j] $ un élément non divisible par $ A [1,1] $, et soit $ q = A [i, j] \ div A [1,1] $ (le reste est tronqué).
Cette fois,
Exemple:
\begin{pmatrix}
2&0&0&0\\
0&2&3&4\\
0&2&2&6
\end{pmatrix}
\rightarrow
\begin{pmatrix}
2&0&-2&0\\
2&2&1&4\\
0&2&2&6
\end{pmatrix}
===
Le formulaire standard Smith peut être généré en répétant les quatre transformations ci-dessus selon l'algorithme suivant.
Ce code python est présenté ci-dessous.
import numpy as np
#Matrice de base sur matrice entière
def Eij(i,j,n):
E = np.eye(n)
E[i,i] = 0
E[j,j] = 0
E[i,j] = 1
E[j,i] = 1
return E
def Ei(i,n):
E = np.eye(n)
E[i,i] = -1
return E
def Ec(i,j,c,n):
E = np.eye(n)
E[i,j] = c
return E
# A[k:,k:]La plus petite valeur absolue autre que 0 ci-dessus est A[k,k]Convertir pour passer à
def moveMN(A,k):
tmp_A = A[k:,k:]
a = np.abs(tmp_A[tmp_A != 0]).min()
i = np.where(np.abs(tmp_A) == a)[0][0] + k
j = np.where(np.abs(tmp_A) == a)[1][0]+ k
P = Eij(k,j,A.shape[1])
invQ = Eij(i,k,A.shape[0])
B = invQ.dot(A).dot(P)
if B[k,k]<0:
Pi =Ei(k,A.shape[1])
B = B.dot(Pi)
P = P.dot(Pi)
return invQ.astype(int),B.astype(int),P.astype(int)
#A[k,k]Utilisant un[k+1:,k]Conversion qui s'ajuste à 0.(Si vous ne pouvez pas l'organiser, l'élément est A[k,k]Plus petite)
def rowR(A,k):
B = A.copy()
invQ = np.eye(A.shape[0])
P = np.eye(A.shape[1])
for i in range(k+1,A.shape[0]):
q = A[i,k]//A[k,k]
#Résidu
r = A[i,k]%A[k,k]
invQi = Ec(i,k,-q,A.shape[0])
B = invQi.dot(B)
invQ = invQi.dot(invQ)
return invQ.astype(int),B.astype(int),P.astype(int)
#A[k,k]Utilisant un[k,k+1]Conversion qui s'ajuste à 0.(Si vous ne pouvez pas l'organiser, l'élément est A[k,k]Plus petite)
def colR(A,k):
B = A.copy()
invQ = np.eye(A.shape[0])
P = np.eye(A.shape[1])
for i in range(k+1,A.shape[1]):
q = A[k,i]//A[k,k]
#Résidu
r = A[k,i]%A[k,k]
Pi = Ec(k,i,-q,A.shape[1])
B = B.dot(Pi)
P = P.dot(Pi)
return invQ.astype(int),B.astype(int),P.astype(int)
# A[k+1:,k+1:]Dans un[k,k]Élément A qui n'est pas divisible par[i,j]UNE[k,k]Conversion pour convertir en résidus de
def remR(A,k):
invQ = np.eye(A.shape[0])
P = np.eye(A.shape[1])
# Find i,j
i = np.where(A[k+1:,k+1:]%A[k,k] !=0)[0][0] +k+1
j = np.where(A[k+1:,k+1:]%A[k,k] !=0)[1][0] +k+1
q = A[i,j]//A[k,k]
r = A[i,j]%A[k,k]
invQi = Ec(i,k,q,A.shape[0])
Pi = Ec(k,j,-1,A.shape[1])
B = invQi.dot(A).dot(Pi)
P = P.dot(Pi)
invQ = invQi.dot(invQ)
return invQ.astype(int),B.astype(int),P.astype(int)
# Main Function
def Smith_Normalization(A):
invQ = np.eye(A.shape[0])
P = np.eye(A.shape[1])
A0 = A.copy()
# limit of optimization
N = 1000
for k in range(min(A0.shape)):
# If A0[k:,k:] is zero matrix, then stop calculation
if np.sum(np.abs(A0[k:,k:]))==0:
break
for i in range(N):
if i == N-1 :
print("Error: Time Out")
# minimize A[k,k]
invQi,A1,Pi = moveMN(A0,k)
invQ = invQi.dot(invQ)
P = P.dot(Pi)
# make zero row A[k+1:,k]
invQi,A2,Pi = rowR(A1,k)
invQ = invQi.dot(invQ)
P = P.dot(Pi)
# if row A2[k+1:,k] is zero vector ?
if np.abs(A2[k+1:,k]).sum() ==0:
# make zero col A[k,k+1:]
invQi,A3,Pi = colR(A2,k)
invQ = invQi.dot(invQ)
P = P.dot(Pi)
# if col A3[k+1:,k] is zero vector ?
if np.abs(A3[k,k+1:]).sum() ==0:
# A[k,k]|A[k+1:,k+1:]?
if np.sum(A3[k+1:,k+1:]%A3[k,k]) == 0:
A0 = A3.copy()
break;
else:
# reduce A[k+1:,k+1:]
invQi,A0,Pi = remR(A3,k)
invQ = invQi.dot(invQ)
P = P.dot(Pi)
else:
A0 = A3.copy()
else:
A0 = A2.copy()
B = A0.copy().astype(int)
P = P.astype(int)
invQ = invQ.astype(int)
return invQ,B,P
# Check result
A = np.random.randint(0,10,(4,5))
invQ,B,P = Smith_Normalization(A)
print(A)
print(invQ)
print(B)
print(P)
--Défini la matrice de base de la matrice entière.
Il était très difficile d'écrire un résumé par rapport au code python. (J'ai cassé quelques explications polies ...)
https://github.com/yuji0001/2020Topology_basic
Reference
Hiroaki Hiraoka, "Structure et topologie des protéines, introduction au groupe d'homologie persistante", Kyoritsu Publishing, 2013 Author Yuji Okamoto [email protected]
Recommended Posts