识别向量是矩形矩阵中另一个向量的倍数

Identify vectors being a multiple of another in rectangular matrix

给定一个 nxm 矩阵 (n > m) 的整数,我想确定是 单个 其他行的倍数的行,所以不是线性组合多个其他行。 我可以将所有行缩放到它们的长度并找到唯一的行,但这很容易出现浮点数问题,并且也不会检测到彼此相反(指向另一个方向)的向量。 有什么想法吗?

例子

A = array([[-1, -1,  0,  0],
       [-1, -1,  0,  1],
       [-1,  0, -1,  0],
       [-1,  0,  0,  0],
       [-1,  0,  0,  1],
       [-1,  0,  1,  1],
       [-1,  1, -1,  0],
       [-1,  1,  0,  0],
       [-1,  1,  1,  0],
       [ 0, -1,  0,  0],
       [ 0, -1,  0,  1],
       [ 0, -1,  1,  0],
       [ 0, -1,  1,  1],
       [ 0,  0, -1,  0],
       [ 0,  0,  0,  1],
       [ 0,  0,  1,  0],
       [ 0,  1, -1,  0],
       [ 0,  1,  0,  0],
       [ 0,  1,  0,  1],
       [ 0,  1,  1,  0],
       [ 0,  1,  1,  1],
       [ 1, -1,  0,  0],
       [ 1, -1,  1,  0],
       [ 1,  0,  0,  0],
       [ 1,  0,  0,  1],
       [ 1,  0,  1,  0],
       [ 1,  0,  1,  1],
       [ 1,  1,  0,  0],
       [ 1,  1,  0,  1],
       [ 1,  1,  1,  0]])

例如第 0 行和 -3 行正好指向相反的方向(1 乘以 -1 使它们相等)。

您可以利用两个归一化线性相关向量的内积给出 1 或 -1 的事实,因此代码可能如下所示:

>>> A_normalized = (A.T/np.linalg.norm(A, axis=-1)).T
>>> M = np.absolute(np.einsum('ix,jx->ij', A_normalized, A_normalized))
>>> i, j = np.where(np.isclose(M, 1))
>>> i, j = i[i < j], j[i < j] # Remove repetitions
>>> print(i, j)

输出:[ 0 2 3 6 7 9 11 13] [27 25 23 22 21 17 16 15]

您可以标准化每行除以其 GCD:

import numpy as np

def normalize(a):
    return a // np.gcd.reduce(a, axis=1, keepdims=True)

并且您可以定义一个距离,将相反的向量视为相等:

def distance(a, b):
    equal = np.all(a == b) or np.all(a == -b)
    return 0 if equal else 1

那么你可以使用标准的聚类方法:

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, fcluster

def cluster(a):
    norm_a = normalize(a)
    distances = pdist(norm_a, metric=distance)
    return fcluster(linkage(distances), t=0.5)

例如:

>>> A = np.array([( 1,  2,  3,  4),
...               ( 0,  2,  4,  8),
...               (-1, -2, -3, -4),
...               ( 0,  1,  2,  4),
...               (-1,  2, -3,  4),
...               ( 2, -4,  6, -8)])
>>> cluster(A)
array([2, 3, 2, 3, 1, 1], dtype=int32)

解释:簇 1 由第 4 行和第 5 行组成,簇 2 由第 0 行和第 2 行组成,簇 3 由第 1 行和第 3 行组成。