通过 sympy,你可以使用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]])
这也适用于矩形矩阵(不仅适用于正方形矩阵)。
rref
是基于列而不是行的。因此需要 .T
。但我可能错了。我写这个答案已经有一段时间了。 - MSeifertnp.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首先,您的第三行与第一行和第二行是线性相关的。然而,您的第一列和第四列是线性相关的。
可以使用两种方法:
特征值
如果矩阵的一个特征值为零,则其相应的特征向量是线性相关的。文档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'
测试行的方法类似。
然后,您可以将其扩展到测试所有向量组合,但我想这种解决方案在规模上会变得糟糕。
[[1,0,1], [1,1,0], [0,0,0]]
,它显然是秩为2的(第三行为0),但你的检查会给出 r1.r2 - r1.r1 * r2.r2 == -1
、r1.r3 - r1.r1 * r3.r3 == -1
和 r2.r3 - r2.r2 * r3.r3 == -1
。你所拥有的检查只能检测一个向量是否是另一个向量的(正)倍数,但向量可以线性相关而没有任何一个共线。 - SirGuy我修改了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
我知道这个问题已经被问了一段时间,但是这里有一个非常简单的(尽管可能效率不高)解决方案。给定一个数组,以下方法通过逐步添加向量并测试秩是否增加来找到一组线性独立的向量:
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)
matrix_rank
会非常昂贵。 - AnnanFayM
。请编辑代码以定义 M 或将 mat
映射到 M
,因为我认为这是你的意思? - tsandohttp://numpy-discussion.10968.n7.nabble.com/Reduced-row-echelon-form-td16486.html
关于以下讨论:
from sympy import *
A = Matrix([[1,1,1],[2,2,2],[1,7,5]])
print(A.nullspace())
B = A.T
print(B.nullspace())
[-2, 1, 0]
,这表明第三行是独立的。您可以使用 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]]
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]))))
简单地说,使用秩来检查线性相关性/独立性。
如果矩阵的秩小于矩阵的维度(在方阵中)或者秩小于最小维度(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.")
A = np.array([[1, 2, 3]
[4, 5 , 6]])
R = np.array([[-4.12310563, -5.33578375, -6.54846188],
[ 0. , -0.72760688, -1.45521375]])