Генерация «случайной» матрицы определенного ранга по фиксированному набору элементов
Я хотел бы создать матрицы размераm
Иксn
и званиеr
с элементами, поступающими из указанного конечного набора, например,{0,1}
или же{1,2,3,4,5}
, Я хочу, чтобы они были "случайными" в некотором очень слабом смысле этого слова, т. е. я хочу получить множество возможных выходных данных из алгоритма с распределением, неопределенно подобным распределению всех матриц по этому набору элементов с указанным рангом.
На самом деле, меня не волнует, имеет ли он званиеr
просто, что этоclose матрице рангаr
(измеряется по норме Фробениуса).
Когда под рукой находятся реальные вещи, я делаю следующее, что вполне соответствует моим потребностям: генерировать матрицыU
размераm
Иксr
а такжеV
изn
Иксr
с элементами, независимо выбранными из, например, Нормальный (0, 2). затемU V'
являетсяm
Иксn
матрица рангаr
(Что ж,<= r
но я думаю, что этоr
с большой вероятностью).
Если я просто сделаю это, а затем округлю до двоичного / 1-5, тем не менее, ранг увеличивается.
Также возможно получить более низкое приближение к матрице, выполнив SVD и взяв первоеr
особые значения. Эти значения, однако, не будут лежать в желаемом наборе, и округление их снова повысит ранг.
Этот вопрос связан, но принятый ответ не является «случайным»; и другой ответ предполагает SVD, который здесь не работает, как отмечено.
Одна возможность, о которой я подумал, состоит в том, чтобыr
линейно независимые векторы строки или столбца из набора, а затем получить остальную часть матрицы с помощью линейных комбинаций из них. Однако я не совсем понимаю, как получить «случайный» результат. линейно независимые векторы или как их объединить квазислучайным образом после этого.
(Не то чтобы это было супер-актуально, но я делаю это просто.)
Update: Я попробовал подход, предложенный EMS в комментариях, с этой простой реализацией:
<code>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) </code>
Это работает быстро для маленьких матриц, но разваливается, например, 10x10 - кажется, застрял на 6 или 7 уровне, по крайней мере, на сотни тысяч итераций.
Кажется, что это могло бы работать лучше с лучшей (то есть менее плоской) целевой функцией, но я не знаю, что это будет.
Я также попробовал простой метод отклонения для построения матрицы:
<code>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 </code>
Это работает нормально, например, Двоичные матрицы 10x10 с любым рангом, но не для 0-4 матриц или гораздо больших двоичных файлов с более низким рангом. (Например, получение двоичной матрицы 20х20 ранга 15 заняло у меня 42 000 отказов; при 20х20 ранга 10 потребовалось 1,2 миллиона).
Это ясно, потому что пространство, охватываемое первымr
ряды слишком малы, часть пространства, из которого я делаю выборку, например,{0,1}^10
, в этих случаях.
Мы хотим пересечение пролета первогоr
строки с набором допустимых значений.
Таким образом, мы могли бы попробовать произвести выборку из диапазона и найти допустимые значения, но поскольку диапазон включает в себя действительные коэффициенты, которые никогда не найдут нас действительными векторами (даже если мы нормализуем так, что, например, первый компонент находится в допустимом наборе) ,
Может быть, это можно сформулировать как задачу целочисленного программирования или как-то так?