如何从矩阵中找到线性无关的行

32
如何从矩阵中识别线性独立的行?例如,

enter image description here

第4行是独立的。

3
如果我没记错的话,“线性独立”是向量组的一个特征。在这个上下文中,我不确定“识别线性独立的行”是什么意思。输出应该是什么? - cel
抱歉表达不够清晰。在这个例子中,rank 是3,因此删除其中一个依赖行(比如第三行)。实际上,我想要追踪的问题在这里 - SparkAndShine
@sparkandshine 你只考虑整数吗? - Karlo
10个回答

32

通过 ,你可以使用sympy.Matrix.rref来找到线性无关的行:

>>> import sympy 
>>> import numpy as np
>>> mat = np.array([[0,1,0,0],[0,0,1,0],[0,1,1,0],[1,0,0,1]])  # your matrix
>>> _, inds = sympy.Matrix(mat).T.rref()   # to check the rows you need to transpose!
>>> inds
[0, 1, 3]

基本上告诉你第0、1和3行是线性无关的,而第2行不是(它是第0行和第1行的线性组合)。

然后您可以使用切片删除这些行:

>>> mat[inds]
array([[0, 1, 0, 0],
       [0, 0, 1, 0],
       [1, 0, 0, 1]])

这也适用于矩形矩阵(不仅适用于正方形矩阵)。


请问您能帮我理解为什么这个简单的例子不起作用吗? col1 col2 col3 6 0.0789 0.9211 -1 0.0000 1.0000 24 0.0233 0.9767 - Alex
可能是因为您使用了浮点数计算。但如果您询问一个包含[mcve]的新问题,您将获得更准确的答案。 - MSeifert
转置操作中是否真的需要使用 .T 才能使其正常工作? - Christian Singer
据我所知,rref 是基于列而不是行的。因此需要 .T。但我可能错了。我写这个答案已经有一段时间了。 - MSeifert
@ChristianSinger,你确实需要转置,否则会得到不同的结果(0、1、2)。 - EdwardG
这似乎不太准确。例如,矩阵 np.array([[-5,3,2],[4,-6,3],[1,3,-5]]) 返回 (0,1),而同样的矩阵除以10,即 [[-.5,.3,.2],[.4,-.6,.3],[.1,.3,-.5]] 返回 (0,1,2) - Maverick Meerkat

32

首先,您的第三行与第一行和第二行是线性相关的。然而,您的第一列和第四列是线性相关的。

可以使用两种方法:

特征值

如果矩阵的一个特征值为零,则其相应的特征向量是线性相关的。文档eig指出返回的特征值按其重复次数而不一定按顺序重复。但是,假设特征值对应于您的行向量,则一种方法是:

import numpy as np

matrix = np.array(
    [
        [0, 1 ,0 ,0],
        [0, 0, 1, 0],
        [0, 1, 1, 0],
        [1, 0, 0, 1]
    ])

lambdas, V =  np.linalg.eig(matrix.T)
# The linearly dependent row vectors 
print matrix[lambdas == 0,:]

柯西-施瓦茨不等式

为了检验向量的线性相关性并找出哪些向量是线性相关的,你可以使用柯西-施瓦茨不等式。基本上,如果向量的内积等于向量范数的乘积,则这些向量是线性相关的。以下是列向量的一个示例:

import numpy as np

matrix = np.array(
    [
        [0, 1 ,0 ,0],
        [0, 0, 1, 0],
        [0, 1, 1, 0],
        [1, 0, 0, 1]
    ])

print np.linalg.det(matrix)

for i in range(matrix.shape[0]):
    for j in range(matrix.shape[0]):
        if i != j:
            inner_product = np.inner(
                matrix[:,i],
                matrix[:,j]
            )
            norm_i = np.linalg.norm(matrix[:,i])
            norm_j = np.linalg.norm(matrix[:,j])

            print 'I: ', matrix[:,i]
            print 'J: ', matrix[:,j]
            print 'Prod: ', inner_product
            print 'Norm i: ', norm_i
            print 'Norm j: ', norm_j
            if np.abs(inner_product - norm_j * norm_i) < 1E-5:
                print 'Dependent'
            else:
                print 'Independent'

测试行的方法类似。

然后,您可以将其扩展到测试所有向量组合,但我想这种解决方案在规模上会变得糟糕。


2
这段代码只适用于方阵吗?原始问题给出了一个方阵的例子,但你也可以使用矩形矩阵来完成这个任务。 - AnnanFay
@Annan,柯西-施瓦茨方法也适用于矩形矩阵,只需更改for循环中的索引范围即可。 - crlb
2
@hakanc 我认为你的柯西-施瓦茨不等式部分是错误的。考虑矩阵 [[1,0,1], [1,1,0], [0,0,0]],它显然是秩为2的(第三行为0),但你的检查会给出 r1.r2 - r1.r1 * r2.r2 == -1r1.r3 - r1.r1 * r3.r3 == -1r2.r3 - r2.r2 * r3.r3 == -1。你所拥有的检查只能检测一个向量是否是另一个向量的(正)倍数,但向量可以线性相关而没有任何一个共线。 - SirGuy
2
我认为这两种方法比使用高斯消元法找到行阶梯形矩阵的方法要慢。 - R zu
将数据缩放到特定范围,然后应用适当的比例。 - Yahya
显示剩余5条评论

6

我修改了Cauchy-Schwartz不等式的代码,使其在维度上的扩展性更好:输入是矩阵及其维数,输出是一个新的矩形矩阵,其中包含原始矩阵中线性独立的列。这个假设是第一列永远不为空,但可以很容易地推广以实现这种情况。另一个我观察到的事情是 1e-5 似乎是一个“粗略”的阈值,因为在那种情况下,一些特定病理向量被发现是线性相关的:1e-4 不会出现同样的问题。希望这能对你有所帮助:对我来说,找到一个真正有效的提取li向量的程序非常困难,因此我愿意分享我的经验。如果发现错误,请报告给我!

from numpy import dot, zeros
from numpy.linalg import matrix_rank, norm

def find_li_vectors(dim, R):

    r = matrix_rank(R) 
    index = zeros( r ) #this will save the positions of the li columns in the matrix
    counter = 0
    index[0] = 0 #without loss of generality we pick the first column as linearly independent
    j = 0 #therefore the second index is simply 0

    for i in range(R.shape[0]): #loop over the columns
        if i != j: #if the two columns are not the same
            inner_product = dot( R[:,i], R[:,j] ) #compute the scalar product
            norm_i = norm(R[:,i]) #compute norms
            norm_j = norm(R[:,j])

            #inner product and the product of the norms are equal only if the two vectors are parallel
            #therefore we are looking for the ones which exhibit a difference which is bigger than a threshold
            if absolute(inner_product - norm_j * norm_i) > 1e-4:
                counter += 1 #counter is incremented
                index[counter] = i #index is saved
                j = i #j is refreshed
            #do not forget to refresh j: otherwise you would compute only the vectors li with the first column!!

    R_independent = zeros((r, dim))

    i = 0
    #now save everything in a new matrix
    while( i < r ):
        R_independent[i,:] = R[index[i],:] 
        i += 1

    return R_independent

1
对我来说非常好用,谢谢。小问题:你应该从numpy.linalg中导入“norm”,而不是从numpy中导入。 - Edgar H
你说得完全正确,谢谢!我可能在修改代码时没有重新运行它,并忘记了 .linalg。 - Simone Bolognini

4

我知道这个问题已经被问了一段时间,但是这里有一个非常简单的(尽管可能效率不高)解决方案。给定一个数组,以下方法通过逐步添加向量并测试秩是否增加来找到一组线性独立的向量:

from numpy.linalg import matrix_rank

def LI_vecs(dim,M):
    LI=[M[0]]
    for i in range(dim):
        tmp=[]
        for r in LI:
            tmp.append(r)
        tmp.append(M[i])                #set tmp=LI+[M[i]]
        if matrix_rank(tmp)>len(LI):    #test if M[i] is linearly independent from all (row) vectors in LI
            LI.append(M[i])             #note that matrix_rank does not need to take in a square matrix
    return LI                           #return set of linearly independent (row) vectors

#Example
mat=[[1,2,3,4],[4,5,6,7],[5,7,9,11],[2,4,6,8]]
LI_vecs(4,mat)

3
这有多有效?我认为调用matrix_rank会非常昂贵。 - AnnanFay
@Annan matrix_rank是使用奇异值分解计算的,因此我的猜测(根据维基百科上描述的标准算法)是它最多为O(n^3),其中n是矩阵的维数。这应该不会太糟糕,但我不知道在实践中如何。 - Couchy
你的函数有一个未定义的变量 M。请编辑代码以定义 M 或将 mat 映射到 M,因为我认为这是你的意思? - tsando

2
我将这个问题解释为寻找与其他行线性独立的行。 这等价于寻找与其他行线性相关的行。
高斯消元并将小于阈值的数字视为零可以做到这一点。这比查找矩阵的特征值、使用柯西-施瓦茨不等式测试所有行的组合或奇异值分解更快。

See: https://math.stackexchange.com/questions/1297437/using-gauss-elimination-to-check-for-linear-dependence

浮点数问题:

http://numpy-discussion.10968.n7.nabble.com/Reduced-row-echelon-form-td16486.html


1

关于以下讨论:

如何使用Matlab查找矩阵的相关行/列?

from sympy import *
A = Matrix([[1,1,1],[2,2,2],[1,7,5]])
print(A.nullspace())

很明显,第一行和第二行是彼此相乘的。如果我们执行上述代码,我们会得到“[-1/3,-2/3,1]”。零空间中零元素的索引显示独立性。但为什么这里的第三个元素不为零?如果我们将A矩阵与零空间相乘,我们会得到一个零列向量。那么问题出在哪里呢?
我们要找的答案是A的转置的零空间。
B = A.T
print(B.nullspace())

现在我们得到了[-2, 1, 0],这表明第三行是独立的。
这里有两个重要的注意事项:
1. 考虑我们想要检查行依赖性还是列依赖性。 2. 注意,除非矩阵是对称的,否则矩阵的零空间不等于其转置的零空间。

1

您可以使用 SymPy库 的Matrix对象的columnspace()方法,基本上可以找到跨越矩阵列空间的向量。它们自动地是矩阵的线性独立列。

import sympy as sp 
import numpy as np 

M  = sp.Matrix([[0, 1, 0, 0],
                [0, 0, 1, 0],
                [1, 0, 0, 1]])


for i in M.columnspace():
    print(np.array(i))
    print()

# The output is following.

# [[0]
#  [0]
#  [1]]

# [[1]
#  [0]
#  [0]]

# [[0]
#  [1]
#  [0]]

0
一个(可能)更短的解决方案:
def cartesian_product(vec1, vec2):
    return np.transpose(np.array([np.repeat(vec1, len(vec2)), np.tile(vec2, len(vec1))]))

def check_dependencies(matrix_):
    def rows_are_dependent(indices):
        if indices[0] == indices[1]:
            return False
        return np.unique(matrix_[indices[0],] / matrix_[indices[1],]).shape[0] == 1

    return np.any(np.apply_along_axis(rows_are_dependent, 1, cartesian_product(np.arange(matrix_.shape[0]), np.arange(matrix_.shape[0]))))

0

简单地说,使用秩来检查线性相关性/独立性。

如果矩阵的秩小于矩阵的维度(在方阵中)或者秩小于最小维度(dim(matrix)),那么该矩阵是线性相关的。

例如,在您的示例中的方阵

import numpy as np

matrix = np.array(
    [
        [0, 1 ,0 ,0],
        [0, 0, 1, 0],
        [0, 1, 1, 0],
        [1, 0, 0, 1]
    ])

rank = np.linalg.matrix_rank(matrix)

#Check if the matrix is linearly independent
if rank == matrix.shape[0]:
    print("The matrix is linearly independent.")
else:
    print("The matrix is linearly dependent.")

对于非方阵的示例

matrix = np.array(
    [
        [0, 1 ,0 ,0],
        [0, 0, 1, 0],
        [0, 1, 1, 0]
    ])
rank = np.linalg.matrix_rank(matrix)

# Check if the matrix is linearly independent

if rank == min(matrix.shape):
    print("The matrix is linearly independent.")
else:
    print("The matrix is linearly dependent.")

0
这里有一个解决方案,使用numpy和$QR$分解来找到矩阵$A$中线性独立和线性相关的列。$Q$是一个正交矩阵,其列张成$A$的范围。$R$是一个上三角矩阵,它显示如何将$A$的列表示为$Q$的列的线性组合。$R$的非零对角线元素所在的前几列是$A$的线性独立列,而其余列是线性相关的。
例如:
A = np.array([[1, 2, 3]
              [4, 5 , 6]])

在这种情况下,$R$ 是:
 R =  np.array([[-4.12310563, -5.33578375, -6.54846188],
                [ 0.        , -0.72760688, -1.45521375]])

由于对角线元素非零,矩阵$A$的前两列线性无关,第三列线性相关。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接