Kürzlich habe ich mit einem Freund ein Buch mit dem Titel "Proteinstruktur und Topologie (Hiraoka)" gelesen. In diesem Buch gab es eine Methode zum Zerlegen einer Ganzzahlmatrix namens Smith-Standardisierung, daher habe ich sie mit Numpy implementiert.
Vor der Erläuterung der Smith-Standardisierung muss die Grundmatrix der Ganzzahlmatrix erläutert werden.
Die Grundmatrix im Vektorraum ist die Matrix zum Durchführen der Grundtransformation der Matrix. (Es sollte im ersten Jahr der Universität gelernt werden.) Es gibt die folgenden drei grundlegenden Transformationen der Matrix.
Jede Transformation kann entweder von links oder von rechts erreicht werden, indem eine ** reguläre ** Matrix angewendet wird (dh eine Matrix, in der eine inverse Matrix existiert). Diese reguläre Matrix wird als Grundmatrix bezeichnet.
Wenn hier die Grundmatrizen, die 1, 2 und 3 entsprechen, als $ P_i (c), Q_ {ij}, R_ {ij} (c) $ ausgedrückt werden, sind sie wie folgt.
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}
Auch die Umkehrung jeder Grundmatrix ist
Betrachten wir nun die Grundmatrix für die Ganzzahlmatrix.
Die Grundmatrix $ Q_ {ij}, R_ {ij} (c) $ im Vektorraum kann als Grundmatrix für die Ganzzahlmatrix übernommen werden, da die Inversmatrix die Ganzzahlmatrix für die Ganzzahl $ c $ ist. Andererseits ist $ P_i (c) $ keine ganzzahlige Matrix, es sei denn, die Konstante $ c $ ist 1 oder -1. Da $ P_i (1) $ eine Einheitsmatrix ist, ist $ P_i (-1) $ eine Grundmatrix einer ganzzahligen Matrix.
Daher sind für die ganze Zahl $ c $ $ P_i (-1), Q_ {ij}, R_ {ij} (c) $ die Grundmatrix der ganzzahligen Matrix.
Die Smith-Standardisierung ist eine Diagonalisierungsmethode für ganzzahlige Matrizen. Die beliebige Ganzzahlmatrix A verwendet die Matrizen $ P und 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)
Kann mit diagonalisiert werden. Hier sind $ P und Q $ die Produkte der Grundmatrizen, und $ c_ {i + 1} $ kann durch $ c_i $ geteilt werden. Da die Matrizen $ P und Q $ die Produkte der Grundmatrizen sind, gibt es eine inverse Matrix von ganzzahligen Matrizen. Die durch das Produkt dieser Grundmatrix diagonalisierte Form wird als Smith-Standardform bezeichnet.
Beispiel:
\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}
Die Smith-Standardform wird durch Wiederholen der folgenden vier Operationen erreicht.
Verschieben Sie das Element $ A [i, j] $ mit dem kleinsten absoluten Wert ungleich Null in die Komponente (1,1). Mit anderen Worten
{\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}
Ist definiert als.
Beispiel:
\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}
Wenn $ q_i = A [i, 0] \ div A [0,0] $ für jedes $ i $ ist (der Rest wird abgeschnitten)
Beispiel:
\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}
Wenn für $ i $ $ q_i = A [0, i] \ div A [0,0] $ ist (der Rest wird abgeschnitten)
Beispiel:
\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}
Sei $ A [i, j] $ ein Element, das nicht durch $ A [1,1] $ teilbar ist, und sei $ q = A [i, j] \ div A [1,1] $ (der Rest wird abgeschnitten).
In diesem Moment,
Beispiel:
\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}
===
Die Smith-Standardform kann durch Wiederholen der obigen vier Transformationen gemäß dem folgenden Algorithmus erzeugt werden.
Dieser Python-Code wird unten vorgestellt.
import numpy as np
#Grundmatrix auf Ganzzahlmatrix
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:]Der kleinste absolute Wert außer 0 oben ist A.[k,k]Konvertieren, um zu verschieben
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]Verwendung einer[k+1:,k]Konvertierung, die auf 0 eingestellt wird.(Wenn Sie es nicht anordnen können, ist das Element A.[k,k]Kleiner)
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ückstand
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]Verwendung einer[k,k+1]Konvertierung, die auf 0 eingestellt wird.(Wenn Sie es nicht anordnen können, ist das Element A.[k,k]Kleiner)
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ückstand
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:]In einem[k,k]Element A, das nicht teilbar ist durch[i,j]EIN[k,k]Umwandlung in Umrechnungen von
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)
Es war sehr schwierig, eine Zusammenfassung im Vergleich zum Python-Code zu schreiben. (Ich habe einige höfliche Erklärungen gebrochen ...)
https://github.com/yuji0001/2020Topology_basic
Reference
Hiroaki Hiraoka, "Proteinstruktur und Topologie, Einführung in die Persistent Homology Group", Kyoritsu Publishing, 2013 Author Yuji Okamoto [email protected]
Recommended Posts