#python #arrays #math #matrix #linear-algebra
Вопрос:
Рассмотрим квадратную матрицу X
, обладающую этим свойством X.sum(axis=0) = X.sum(axis=1) = x1
, для некоторого вектора x1
. Теперь подумайте, что X
трансформируется Y
вот так:
Y = np.exp(np.log(X) R)
где R
— случайная матрица, обладающая свойством R.mean(axis=0) = R.mean(axis=1) = 0
.
Я хотел бы масштабировать строки и столбцы Y
так, чтобы: Y.sum(axis=0) = Y.sum(axis=1) = x1
. Я обнаружил на практике, что этого можно достичь таким образом:
for i in range(N):
Y *= (x1 / Y.sum(axis=0))[None,:]
Y *= (x1 / Y.sum(axis=1))[:,None]
который асимптотически сходится с большим N
. Однако это, очевидно, неэффективно с вычислительной точки зрения, и я подозреваю, что существует аналитическое решение. Есть ли? Спасибо,
МВЕ
import numpy as np
_ = None
# define x1 amp; X
x1 = np.array([.5,.3,.2])
X = x1[:,_] * x1[_,:]
# define R
R = np.random.random(X.shape)
R = R - R.mean(axis=0)[_,:]
R = R - R.mean(axis=1)[:,_]
# transformation
Y = np.exp(np.log(X) R)
# normalizing Y -> this is what I hope can be replaced with analytical solution
for i in range(1000):
Y *= (x1 / Y.sum(axis=0))[_,:]
Y *= (x1 / Y.sum(axis=1))[:,_]
# end normalizing, show it worked
print(Y.sum(axis=0))
print(Y.sum(axis=1))
Комментарии:
1. Я перекрестно разместил этот вопрос в math.stackexchange и получил ответ там , в комментариях. Похоже, аналитического решения не существует, и этот итерационный алгоритм хорошо известен («Синкхорн-Кнопп»). В моем приложении конвергенция была на самом деле достаточно быстрой: около 1 десятичной точки точности на итерацию.