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

原理

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

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


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

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

​  还是先来看下教材上的例子,设A为:

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=002212110022104250010614125

​  经过初等行变换,得到
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=100001001000001010050297130
​  由非零行首个非零元素对应的列向量组成极大无关组,设为B矩阵,即取出第1,2,4列:
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=002212111042

​  而C矩阵即为化简后阶梯型的非零行
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=100010100001100529713


代码实现

@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

未完待续

  继续填坑!

在这里插入图片描述

  • 1
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值