原理

这一篇我们来介绍满秩分解（Full-Rank-Factorization），其定义如下：设 $A\in F^{m\times n}$ ，rank(A)=r，若存在秩为r的矩阵 $B\in F^{m\times r},C\in F^{r\times n}$ ，使得 $A = B C$ ，称该分解为A的满秩分解。从秩的大小来看，我们可以知道， $r < m ∣ n$ ，所以满秩分解把A分解成了一个“瘦高”矩阵乘上一个“胖矮”矩阵。

对于任何非零的矩阵 $A\in F^{m\times n}$ ，都存在满秩分解。 这是一个很强的存在性，我们也可以这么理解，只要不是零矩阵，其极大无关组的个数至少是1，其他列向量可以被该极大无关组线性表出，满秩分解可以理解为是找极大无关组并给出线性表示原矩阵的过程。

基于初等行变换的满秩分解

首先满秩分解可以有很多不同的求法，这里我们选择计算较为简单，不需要引入求逆（求逆往往是不稳定的操作），就可以得到矩阵A的满秩分解。

还是先来看下教材上的例子，设A为：

$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=\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}$
由非零行首个非零元素对应的列向量组成极大无关组，设为B矩阵，即取出第1，2，4列：
$B=\begin{bmatrix} 0 & 1 & -1 \\ 0 & 2 & 0 \\ 2 & -1 & -4 \\ -2 & 1 & 2 \end{bmatrix}$

而C矩阵即为化简后阶梯型的非零行：
$C=\begin{bmatrix} 1 & 0 & 1 & 0 & -10 & -29 \\ 0 & 1 & 0 & 0 & 0 & -7 \\ 0 & 0 & 0 & 1 & -5 & -13\\ \end{bmatrix}$

代码实现

@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, 强迫症
    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 @_@

填坑

来填坑啦，还记得我们在LU解线性方程组中使用到需要判断矩阵是否满秩的函数，衍生一下还有行满秩和列满秩，都可以用上面的满秩分解来得到。原理很简单，直接上代码，秒懂！

@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

未完待续

继续填坑！

在这里插入图片描述

矩阵论代码实践之满秩分解

胖瘦通吃——满秩分解

原理

基于初等行变换的满秩分解

代码实现

example

填坑

未完待续