Recently, I read a book called "Protein Structure and Topology (Hiraoka)" with a friend. In that book, there was an integer matrix factorization method called Smith standardization, so I implemented it with Numpy.
Before explaining Smith standardization, we need to explain the elementary matrix of integer matrices.
The basic matrix on the vector space is the matrix for performing the basic transformation of the matrix. (You should be learning in the first year of college.). There are the following three basic transformations of a matrix.
Each transformation can be achieved from either the left or right by applying a ** regular ** matrix (that is, a matrix with an inverse matrix). This regular matrix is called the elementary matrix.
Here, if the elementary matrices corresponding to 1, 2, and 3 are expressed as $ P_i (c), Q_ {ij}, R_ {ij} (c) $, they are as follows.
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}
Also, the inverse matrix of each elementary matrix is
Now, let's consider the elementary matrix for an integer matrix.
The basic matrix $ Q_ {ij}, R_ {ij} (c) $ on the vector space can be adopted as the basic matrix for the integer matrix because the inverse matrix is the integer matrix for the integer $ c $. On the other hand, $ P_i (c) $ is not an integer matrix unless the constant $ c $ is 1 or -1. Since $ P_i (1) $ is an identity matrix, $ P_i (-1) $ is the basic matrix of an integer matrix.
Therefore, for the integer $ c $, $ P_i (-1), Q_ {ij}, R_ {ij} (c) $ are the basic matrix of the integer matrix.
Smith standardization is a diagonalization method for integer matrices. Arbitrary integer matrix A uses matrices $ P and 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)
Can be diagonalized with. Here, $ P and Q $ are the products of the elementary matrices, and $ c_ {i + 1} $ can be divided by $ c_i $. Since the matrices $ P and Q $ are the products of the elementary matrices, there is an inverse matrix of integer matrices. The form diagonalized by the product of this elementary matrix is called the Smith standard form.
Example:
\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}
The Smith standard form is achieved by repeating the following four operations.
Move the element $ A [i, j] $ with the smallest nonzero absolute value to the (1,1) component. In other words
{\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}
Is defined as.
Example:
\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}
For any $ i $, if $ q_i = A [i, 0] \ div A [0,0] $ (the remainder is truncated),
Example:
\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}
For any $ i $, if $ q_i = A [0, i] \ div A [0,0] $ (the remainder is truncated)
Example:
\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}
Let $ A [i, j] $ be an element that is not divisible by $ A [1,1] $, and let $ q = A [i, j] \ div A [1,1] $ (the remainder is truncated).
This time,
Example:
\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}
===
The Smith standard form can be generated by repeating the above four transformations according to the following algorithm.
This python code is introduced below.
import numpy as np
#Elementary matrix on an integer matrix
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:]The smallest absolute value other than 0 above is A[k,k]Convert to move to
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]Using, A[k+1:,k]Conversion that adjusts to 0.(If you can't arrange it, the element is A[k,k]Smaller)
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]
#Residue
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]Using, A[k,k+1]Conversion that adjusts to 0.(If you can't arrange it, the element is A[k,k]Smaller)
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]
#Residue
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 A[k,k]Element A that is not divisible by[i,j]A[k,k]Conversion to convert to residuals of
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)
--Defined the elementary matrix of an integer matrix. --Introduced Smith standard form and its conversion method. --I created a python code to convert to Smith standard form.
It was very difficult to write a summary compared to the python code. (I've broken some polite explanations ...)
https://github.com/yuji0001/2020Topology_basic
Reference
Hiroaki Hiraoka, "Protein Structure and Topology, Introduction to Persistent Homology Group", Kyoritsu Shuppan, 2013 Author Yuji Okamoto [email protected]
Recommended Posts