Аналитическое решение для итеративной (асимптотической) нормализации матрицы строк и столбцов

#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 десятичной точки точности на итерацию.