在固定的元素集上生成某个等级的"随机"矩阵

Question

我想生成大小为mx n和等级的矩阵r,元素来自指定的有限集,例如{0,1}或{1,2,3,4,5}.我希望它们在某个非常松散的意义上是"随机的",即我希望从算法获得各种可能的输出,其分布模糊地类似于具有指定等级的那组元素上的所有矩阵的分布.

事实上,我实际上并不关心它是否具有等级r,只是它接近于等级矩阵r(由Frobenius规范测量).

当设置在手是实数,我已经执行以下操作,这是完全足够我需要:生成矩阵U尺寸的mX r和V的nX r,从例如普通独立地采样元件(0,2).然后U V'是一个排名的mx n矩阵r(好吧<= r,但我认为这很有r可能).

如果我只是这样做然后舍入到二进制/ 1-5,则等级增加.

通过进行SVD​​并获取第一个r奇异值,也可以得到矩阵的低秩近似.但是,这些值不会出现在所需的集合中,并且对它们进行舍入将再次提高等级.

这个问题是相关的,但接受的答案不是"随机的",另一个答案表明SVD,这在这里不起作用.

我想过的一种可能性是,使r从集线性无关的行或列向量,然后通过这些的线性组合获得其矩阵的其余部分.但是,我不是很清楚如何获得"随机"线性独立向量,或者如何在此之后以准随机方式组合它们.

(并不是说这是超级相关的,但我是在numpy中这样做.)

---

更新:我在评论中尝试了EMS建议的方法,这个简单的实现:

real = np.dot(np.random.normal(0, 1, (10, 3)), np.random.normal(0, 1, (3, 10)))
bin = (real > .5).astype(int)
rank = np.linalg.matrix_rank(bin)
niter = 0

while rank > des_rank:
    cand_changes = np.zeros((21, 5))
    for n in range(20):
        i, j = random.randrange(5), random.randrange(5)
        v = 1 - bin[i,j]
        x = bin.copy()
        x[i, j] = v
        x_rank = np.linalg.matrix_rank(x)
        cand_changes[n,:] = (i, j, v, x_rank, max((rank + 1e-4) - x_rank, 0))
    cand_changes[-1,:] = (0, 0, bin[0,0], rank, 1e-4)

    cdf = np.cumsum(cand_changes[:,-1])
    cdf /= cdf[-1]
    i, j, v, rank, score = cand_changes[np.searchsorted(cdf, random.random()), :]
    bin[i, j] = v
    niter += 1
    if niter % 1000 == 0:
        print(niter, rank)

它适用于小型矩阵但速度很快,例如10x10 - 它似乎陷入了第6或第7级,至少数十万次迭代.

看起来这可能会更好地用更好(即不太平坦)的目标函数,但我不知道那会是什么.

---

我还尝试了一种简单的拒绝方法来构建矩阵:

def fill_matrix(m, n, r, vals):
    assert m >= r and n >= r
    trans = False
    if m > n: # more columns than rows I think is better
        m, n = n, m
        trans = True

    get_vec = lambda: np.array([random.choice(vals) for i in range(n)])

    vecs = []
    n_rejects = 0

    # fill in r linearly independent rows
    while len(vecs) < r:
        v = get_vec()
        if np.linalg.matrix_rank(np.vstack(vecs + [v])) > len(vecs):
            vecs.append(v)
        else:
            n_rejects += 1
    print("have {} independent ({} rejects)".format(r, n_rejects))

    # fill in the rest of the dependent rows
    while len(vecs) < m:
        v = get_vec()
        if np.linalg.matrix_rank(np.vstack(vecs + [v])) > len(vecs):
            n_rejects += 1
            if n_rejects % 1000 == 0:
                print(n_rejects)
        else:
            vecs.append(v)
    print("done ({} total rejects)".format(n_rejects))

    m = np.vstack(vecs)
    return m.T if trans else m

这适用于例如具有任何等级的10x10二进制矩阵,但不适用于0-4矩阵或具有较低等级的更大二进制矩阵.(例如,获得排名15的20x20二进制矩阵花了我42,000个拒绝;排名10的20x20,花了120万.)

这显然是因为第一r行所跨越的空间对于我从中采样的空间的一部分来说太小,例如{0,1}^10在这些情况下.

我们希望第一r行的跨度与有效值集的交集.因此我们可以尝试从跨度中采样并寻找有效值,但由于跨度涉及实值系数,它们永远不会找到我们有效的向量(即使我们规范化,例如第一个组件在有效集中).

也许这可以表示为整数规划问题,或者什么？

Answer 1

我的朋友丹尼尔·约翰逊（Daniel Johnson）在上面发表了评论，他提出了一个想法，但我看到他从未发布过。它不是很充实，但你也许能够适应它。

如果A是 m×r 且B是 r×n 并且两者的秩均为 r，则AB具有秩 r。现在，我们只需选择A仅在给定集合中具有值的值B。AB最简单的情况是S = {0,1,2,...,j}。

一种选择是A使用适当的行/列总和来制作二进制，以保证正确的排名，并且B列总和添加到不超过j（以便产品中的每个术语都在S）中，并选择行总和来导致排名r（或至少鼓励可以将其用作拒绝）。

我只是认为我们可以提出两个独立的采样方案A，B这比尝试立即攻击整个矩阵更简单、更快。不幸的是，我所有的矩阵采样代码都在另一台计算机上。我知道它很容易推广到允许比{0,1}(ie S) 更大的集合中的条目，但我不记得计算是如何随着m*n.