I tried Smith standardizing an integer matrix with Numpy

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.

Elementary matrix

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.

  1. Multiply a row or column by a constant (excluding 0).
  2. Swap two rows or columns.
  3. Add a constant multiple of one row or column to another.

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 1.~ P_i(c)^{-1} = P_i(1/c) 2.~Q_{ij}^{-1} = Q_{ij} 3.~R_{ij}(c)^{-1} = R_{ij}(-c) Therefore, the regularity of the elementary matrix is derived. The basic theory of linear algebra such as the sweep method is derived using this basic matrix.

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 standard type

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.

  1. $ A_ {t + 1} = {\ rm moveMN} (A_t) $: (1,1) A transformation that moves the smallest value to an element.

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}
  1. $ A_ {t + 1} = {\ rm rowR} (A_t) $: Conversion to make the leftmost column excluding $ A [1,1] $ smaller.

For any $ i $, if $ q_i = A [i, 0] \ div A [0,0] $ (the remainder is truncated), $ {\rm rowR}(A_t) := \prod_{i>1} R_{i,0}(-q_i) A_t $ Is defined as.

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}
  1. $ A_ {t + 1} = {\ rm colR} (A_t) $: Transform to make the top row smaller except $ A [1,1] $.

For any $ i $, if $ q_i = A [0, i] \ div A [0,0] $ (the remainder is truncated) $ {\rm colR}(A_t) := \prod_{i>1} A_t R_{0,i}(-q_i) $ Is defined as.

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}
  1. $ A_ {t + 1} = {\ rm remR} (A_t) $: A conversion that excludes elements that are not divisible by $ A [1,1] $ in $ A [2 :, 2:] $.

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, $ {\rm remR}(A_t) := R_{0,j}(q) A_t R_{i,0}(-1) $ Is defined as.

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.

  1. A = moveMN(A)
  2. A = rowR(A)
  3. If A [2 :, 1]! = 0, return to 1.
  4. A = colR(A)
  5. If A [1,2:]! = 0, return to 1.
  6. If all the elements of A [2 :, 2:] are not divisible by A [1,1], then A = remR (A) and 1. Return to.
  7. Set A = A [2 :, 2:] and return to 1.

This python code is introduced below.

Python code

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)

Summary

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

Code details

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

I tried Smith standardizing an integer matrix with Numpy
I tried sending an email with python.
I tried non-negative matrix factorization (NMF) with TensorFlow
I tried to detect an object with M2Det!
I tried sending an email with SendGrid + Python
Matrix concatenation with Numpy
I tried to implement an artificial perceptron with python
I tried to make an OCR application with PySimpleGUI
I tried to find an alternating series with tensorflow
I tried fp-growth with python
I tried scraping with Python
I wrote GP with numpy
I tried Learning-to-Rank with Elasticsearch!
Try matrix operation with NumPy
I tried clustering with PyCaret
I tried gRPC with Python
I tried scraping with python
I made an N-dimensional matrix operation library Matft with Swift
I tried sending an email from Amazon SES with Python
I tried to create an article in Wiki.js with SQLAlchemy
I tried trimming efficiently with OpenCV
I tried summarizing sentences with summpy
I tried machine learning with liblinear
I tried web scraping with python.
I tried moving food with SinGAN
I sent an SMS with Python
I tried implementing DeepPose with PyTorch
I tried sending an email from the Sakura server with flask-mail
I tried face detection with MTCNN
I tried running prolog with python 3.8.2.
I tried SMTP communication with Python
I tried sentence generation with GPT-2
I tried learning LightGBM with Yellowbrick
I tried face recognition with OpenCV
I tried to make an image similarity function with Python + OpenCV
I tried to implement deep learning that is not deep with only NumPy
I tried to make an open / close sensor (Twitter cooperation) with TWE-Lite-2525A
I get an error with import pandas.
I tried multiple regression analysis with polynomial regression
I tried using Amazon SQS with django-celery
I tried to implement Autoencoder with TensorFlow
I tried linebot with flask (anaconda) + heroku
I tried to get started with Hy
I tried scraping Yahoo News with Python
I tried using Selenium with Headless chrome
I tried factor analysis with Titanic data!
I tried learning with Kaggle's Titanic (kaggle②)
When I get an error with PyInstaller
I tried non-photorealistic rendering with Python + opencv
I tried a functional language with Python
I tried batch normalization with PyTorch (+ note)
I tried recursion with Python ② (Fibonacci sequence)
I tried implementing DeepPose with PyTorch PartⅡ
I tried to implement CVAE with PyTorch
I tried playing with the image with Pillow
I made a life game with Numpy
I tried to solve TSP with QAOA
I tried simple image recognition with Jupyter
I tried CNN fine tuning with Resnet
I tried natural language processing with transformers.
#I tried something like Vlookup with Python # 2