principle
in this article, let's introduce full rank factorization, which is defined as follows: A ∈ F m × n A\in F^{m\times n} A∈Fm × n. rank(A)=r, if there is a matrix with rank r B ∈ F m × r , C ∈ F r × n B\in F^{m\times r},C\in F^{r\times n} B∈Fm × r,C∈Fr × n. Make A = B C A=BC A=BC, which is called the full rank decomposition of A. From the size of the rank, we can know that, r < m ∣ n r< m|n R < m ∣ n, so the full rank decomposition decomposes A into A "thin high" matrix multiplied by A "fat low" matrix.
for any non-zero matrix A ∈ F m × n A\in F^{m\times n} A∈Fm × n. There is full rank decomposition. This is a strong existence. We can also understand that as long as it is not a zero matrix, the number of its maximum independent groups is at least 1, other column vectors can be linearly expressed by the maximum independent group, and the full rank decomposition can be understood as a process of finding the maximum independent group and giving a linear representation of the original matrix.
Full rank decomposition based on elementary row transformation
first, there are many different solutions to the full rank decomposition. Here, we choose a simple calculation. We can obtain the full rank decomposition of matrix A without introducing inversion (inversion is often an unstable operation).
let's take A look at the examples in the textbook first. Let A be:
A = [ 0 1 0 − 1 5 6 0 2 0 0 0 − 14 2 − 1 2 − 4 0 1 − 2 1 − 2 2 10 25 ] A=\begin{bmatrix} 0 & 1 & 0 & -1 & 5 & 6 \\ 0 & 2 & 0 & 0 & 0 & -14 \\ 2 & -1 & 2 & -4 & 0 & 1\\ -2 & 1 & -2 & 2 & 10 & 25\\ \end{bmatrix} A=⎣⎢⎢⎡002−212−11002−2−10−42500106−14125⎦⎥⎥⎤
after elementary line transformation, it is obtained
A
=
[
1
0
1
0
−
10
−
29
0
1
0
0
0
−
7
0
0
0
1
−
5
−
13
0
0
0
0
0
0
]
A=\begin{bmatrix} 1 & 0 & 1 & 0 & -10 & -29 \\ 0 & 1 & 0 & 0 & 0 & -7 \\ 0 & 0 & 0 & 1 & -5 & -13\\ 0 & 0 & 0 & 0 & 0 & 0\\ \end{bmatrix}
A=⎣⎢⎢⎡1000010010000010−100−50−29−7−130⎦⎥⎥⎤
the maximum independent group is composed of the column vector corresponding to the first non-zero element of the non-zero row, which is set as the B matrix, that is, take out the first, second and fourth columns:
B
=
[
0
1
−
1
0
2
0
2
−
1
−
4
−
2
1
2
]
B=\begin{bmatrix} 0 & 1 & -1 \\ 0 & 2 & 0 \\ 2 & -1 & -4 \\ -2 & 1 & 2 \end{bmatrix}
B=⎣⎢⎢⎡002−212−11−10−42⎦⎥⎥⎤
and the C matrix is the non-zero row of ladder type after simplification:
C
=
[
1
0
1
0
−
10
−
29
0
1
0
0
0
−
7
0
0
0
1
−
5
−
13
]
C=\begin{bmatrix} 1 & 0 & 1 & 0 & -10 & -29 \\ 0 & 1 & 0 & 0 & 0 & -7 \\ 0 & 0 & 0 & 1 & -5 & -13\\ \end{bmatrix}
C=⎣⎡100010100001−100−5−29−7−13⎦⎤
code implementation
@classmethod
def Full_Rank_Fact(self, target, eps=1e-6, test=False):
"""Full Rank Decomposition Solution for 2D-matrix by Junno
Args:
target ([np.darray]): [matrix to be processed]
eps ([float]): numerical threshold.
test (bool, optional): [whether to show testing information]. Defaults to False.
Returns:
[np.darray]: B, C
Last edited: 22-01-02
Author: Junno
"""
assert len(np.shape(target)) == 2
# get row ans col numbers
M, N = np.shape(target)
# A for processing and raw_A for result extraction for B
raw_A = deepcopy(target)
A = deepcopy(target)
if test:
print('original matrix:')
print(target)
for i in range(M):
for j in range(i, N):
if test:
print('During process in row:{}, col:{}'.format(i, j))
# resort rows if needed
if sum(abs(A[i:, j])) > eps:
zero_index = []
non_zero_index = []
for k in range(i, M):
if abs(A[k, j]) > eps:
non_zero_index.append(k)
else:
zero_index.append(k)
non_zero_index = sorted(
non_zero_index, key=lambda x: abs(A[x, j]))
sort_cols = non_zero_index+zero_index
if test:
print("Sort_cols index: {}".format(sort_cols))
A[i:, :] = A[sort_cols, :]
if test:
print('After resorting')
print(A)
# do primary row transformation to eliminate elements belows
prefix = -A[i+1:, j]/A[i, j]
temp = (np.array(prefix).reshape(-1, 1)
)@A[i, :].reshape((1, -1))
A[i+1:, :] += temp
# eliminate elements aboves
for k in range(i):
if A[k, j] != 0:
A[k, :] += -(A[k, j]/A[i, j])*A[i, :]
if test:
print('After primary row transformation:')
print(A)
break
else:
continue
# list to record indexs of principal maximum
principal_indexs = [[], []]
for m in range(M):
for n in range(m, N):
if A[m, n] != 0:
principal_indexs[0].append(n)
principal_indexs[1].append(m)
# normalize
if A[m, n] != 1:
A[m, :] *= 1./A[m, n]
break
# -0 to 0, OCD
A[abs(A) == 0] = 0
if test:
print('final result')
print(A)
print('principal maximum index: {}'.format(principal_indexs))
# choose principal maximum and return B and C
if len(principal_indexs[0]):
B = raw_A[:, principal_indexs[0]]
C = A[principal_indexs[1], :]
return B, C
else:
raise ValueError("Can't find any principle maximum @_@")
example
# test standard example
>>> A=np.array([[0,1,0,-1,5,6],[0,2,0,0,0,-14],[2,-1,2,-4,0,1],[-2,1,-2,2,10,25]]).reshape((4,6)).astype(np.float32)
>>> A
array([[ 0., 1., 0., -1., 5., 6.],
[ 0., 2., 0., 0., 0., -14.],
[ 2., -1., 2., -4., 0., 1.],
[ -2., 1., -2., 2., 10., 25.]], dtype=float32)
>>> B,C=Matrix_Solutions.Full_Rank_Fact(A)
>>> B
array([[ 0., 1., -1.],
[ 0., 2., 0.],
[ 2., -1., -4.],
[-2., 1., 2.]], dtype=float32)
>>> C
array([[ 1., 0., 1., 0., -10., -29.],
[ 0., 1., 0., 0., 0., -7.],
[ 0., 0., 0., 1., -5., -13.]], dtype=float32)
>>> B@C-A
array([[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0.]], dtype=float32)
# test unfull-rank matrix
>>> A=np.array([[1,2,3],[2,4,6],[3,6,9,]]).reshape((3,3)).astype(np.float32)
>>> A
array([[1., 2., 3.],
[2., 4., 6.],
[3., 6., 9.]], dtype=float32)
>>> B,C=Matrix_Solutions.Full_Rank_Fact(A)
>>> B
array([[1.],
[2.],
[3.]], dtype=float32)
>>> C
array([[1., 2., 3.]], dtype=float32)
# test all-zeros matrix
>>> A=np.array([[0,0,0],[0,0,0],[0,0,0,]]).reshape((3,3)).astype(np.float32)
>>> A
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]], dtype=float32)
>>> B,C=Matrix_Solutions.Full_Rank_Fact(A)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "............", line 104, in Full_Rank_Fact
raise ValueError("Can't find any principle maximum @_@")
ValueError: Can't find any principle maximum @_@
Pit filling
come and fill the pit. Remember we were there LU solution of linear equations The function to judge whether the matrix is full rank is used in. Derived, there are row full rank and column full rank, which can be obtained by the above full rank decomposition. The principle is very simple, directly on the code, seconds to understand!
@classmethod
def Check_full_rank(self, target, eps=1e-6):
"""Check_full_rank Function
Args:
target ([np.darray]): input matrix
eps ([float]): numerical threshold.
Returns:
[bool]: [whether full rank]
Last edited: 22-01-02
Author: Junno
"""
if np.sum(target)<eps:
# all-zero
return False
# get row ans col numbers
M, N = np.shape(target)
B, C = self.Full_Rank_Fact(target, False)
return B.shape[1] == min(M, N)
@classmethod
def Get_col_rank(self, target, eps=1e-6):
"""Get_col_rank of target
Args:
target ([np.darray]): input matrix
Returns:
[int]: [rank of cal]
Last edited: 22-01-02
Author: Junno
"""
if np.sum(target)<eps:
# all-zero
return 0
B, C = self.Full_Rank_Fact(target, False)
return B.shape[1]
@classmethod
def Get_row_rank(self, target):
"""Get_row_rank of target
Args:
target ([np.darray]): input matrix
Returns:
[int]: [rank of row]
Last edited: 22-01-02
Author: Junno
"""
# row rank = col rank
return Get_col_rank(target, False)
>>> A=np.array([[1,2,3],[2,4,6],[3,6,9,]]).reshape((3,3)).astype(np.float32)
>>> Matrix_Solutions.Check_full_rank(A)
False
>>> Matrix_Solutions.Get_row_rank(A)
1
>>> Matrix_Solutions.Get_col_rank(A)
1
>>> A=np.zeros((4,4))
>>> Matrix_Solutions.Get_col_rank(A)
0
>>> A=np.eye(4)
>>> A
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]])
>>> Matrix_Solutions.Get_col_rank(A)
4
Unfinished to be continued
continue to fill the pit!
