Scipy другой результат распределительного закона - PullRequest
Новая
       37

Scipy другой результат распределительного закона

0 голосов
/ 02 февраля 2019

Я использую scipy для вычисления алгоритма двойного аффинного масштабирования.

В заключительной части шага итерации я получаю разные результаты при вычислении контрольного значения sk. 

s_k_control = c - (AT.dot(y_k) + AT.dot(t_k * (AHAT_inv.dot(b))))
print("SK_control 1:",np.min(s_k_control))

s_k_control = c - AT.dot(y_k + (t_k * (AHAT_inv.dot(b))))
print("SK_control 2:",np.min(s_k_control))

y_k_1 = y_k + t_k * (AHAT_inv.dot(b))

s_k_control = c - AT.dot(y_k_1)
print("SK_control 3:",np.min(s_k_control))

t_k isвсе скалярные и другие переменные являются разреженными матрицами (csc_matrix)

Если я не совсем не прав, из-за закона распределения точечного произведения ( wiki ) приведенный выше код должен возвращать тот же результат ввсе три случая.

Вместо этого я получаю следующий результат: 

SK_control 1: 0.026123046875
SK_control 2: 0.0
SK_control 3: 0.0

Что я могу сделать, чтобы вычислить y_k_1 таким образом, чтобы последующее вычисление sk давало тот же результат, что и первый элемент управления?

Редактировать:

Вот исходная проблема:

Существует ограничение: s_next = c - A.T * y_next >= 0

Использовать ограничение и следующиеформула для расчета размера шага t: 

y_next = y_prev + t(AHA.T)^(-1)*b
  • A - разреженная матрица формы (355,729)
  • c - вектор единиц (729,1)
  • b является вектором единиц (355,1)
  • первый y_prev (y_0) - это вектор нулей (355,1)
  • t - скаляр, но чтобы найти его, мне нужно вычислить вектор t*, а затем взять наименьший элемент из t*и умножить его на некоторый коэффициент 0 < beta < 1 (обычно 0.9 или аналогичный)

Что я пробовал:

  1. Аналитически вычислить t* (повставка формулы в ограничение с помощью s_next = 0: 

    t* = (c - A.T * y_prev)/(A.T(AHA.T)^(-1) * b)

  2. вычисления y_next с использованием scipy.sparse.linalg.spsolve из ограничения: 

    A.T * y_next = c

    (на самом деле A * A.T * y_next = A * c, потому что A и, следовательно, A.T не являются квадратичными)

    , а затем вычисляет t* следующим образом: 

    t* = (y_next - y_prev)/(AHA.T)^(-1) * b

Ни то, ни другоеМетод дает ожидаемый результат.

РЕДАКТИРОВАТЬ 2:

Кажется, у меня возникли проблемы с вычислением обратного (AHA.T).Когда я проверяю это (AHA.T) * (AHA.T) ^ - 1, я не получаю единичную матрицу, но что-то совершенно случайное: 

AHAT = (A * H * A.T).todense()
print(AHAT)
[[9. 0. 0. ... 0. 0. 0.]
 [0. 9. 0. ... 0. 0. 0.]
 [0. 0. 9. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 1. 0. 0.]
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 0. 1.]]
print(np.dot(np.linalg.inv(AHAT), AHAT))
[[ 6.43630605e+16 -3.53205583e+17  1.82309332e+16 ... -2.47507371e+15
   1.93886558e+15 -4.17941428e+15]
 [ 7.72005634e+16 -1.32187302e+17 -1.82278681e+16 ... -1.51094942e+16
  -1.67465411e+16  4.12101169e+15]
 [ 8.58099974e+14  1.89665457e+16 -1.23638446e+16 ... -3.81892219e+15
  -2.21686073e+15  2.61939698e+15]
 ...
 [ 4.44089210e-16 -5.32907052e-15  1.72084569e-15 ...  1.00000000e+00
   1.66533454e-16  1.31838984e-16]
 [ 6.66133815e-16 -1.77635684e-15  1.99840144e-15 ... -8.18789481e-16
   1.00000000e+00 -7.97972799e-16]
 [-1.11022302e-15  3.10862447e-15 -4.44089210e-16 ...  9.99200722e-16
   3.88578059e-16  1.00000000e+00]]

Инверсия сама по себе выглядит так: 

[[-6.10708114e+14 -4.24172270e+16 -1.62348045e+14 ... -1.80059454e-01
   7.58665399e-02  9.93203316e-01]
 [-2.81790056e+15 -8.26584741e+15  3.08108915e+14 ... -1.06861647e+02
  -1.96226676e-01  7.66381784e-01]
 [-4.36162847e+13  5.27325574e+15 -1.20358871e+15 ... -7.79860964e+00
  -3.24595030e-01 -1.50920847e-01]
 ...
 [ 9.20066618e-02  1.39924661e+00 -5.81619213e-02 ...  1.52230844e+00
   1.14720794e-02 -3.70994069e-02]
 [-2.31455053e-01  2.33160131e+00 -3.65460727e-02 ...  1.14720794e-02
   1.52976177e+00 -1.95029578e-02]
 [ 1.44223299e-01 -1.10202460e+00 -7.57449990e-02 ... -3.70994069e-02
  -1.95029578e-02  1.52462702e+00]]

Есть ли возможность избежать использования обратного и при этом вычислить t_k?

Ответы [ 2 ]

0 голосов
/ 03 февраля 2019

Хорошо, я нашел решение.

Правильный путь - полностью избежать вычисления обратного.

Если у нас есть линейная задача Ax = b, решение: x = A^(-1)*b

В Python уже есть функции для расчета этого решения, например, scipy.sparse.linalg.spsolve.

Так что, если мне нужно вычислить (AHA.T)^(-1) * b, мне просто нужно вызвать spsolve(A * H * A.T, b).

Мой расчет для t выглядит следующим образом: t_k = (c - A.T * y_k)/(A.T * spsolve(A * H * A.T, b))

Это дает мне стабильное решение, и весь мой алгоритм сходится за 20 итераций.

0 голосов
/ 02 февраля 2019

Вы правы, что все 3 версии должны вычислять один и тот же результат.Могут быть некоторые различия, поскольку арифметика с плавающей точкой не является ни ассоциативной, ни дистрибутивной: 

In [147]: ((0.1+0.2)+0.3) != (0.1+(0.2+0.3))
Out[147]: True

In [153]: 0.3*(0.1+0.2) != 0.3*0.1 + 0.3*0.2
Out[153]: True

Соедините это с умножением на большое число: 

In [164]: 1e15 * ((0.1+0.2)+0.3) - 1e15 * (0.1+(0.2+0.3))
Out[164]: 0.125

, и расхождение может стать значительным.

Но в типичном случае ваш код работает должным образом: 

import numpy as np
import scipy.sparse as sparse
# np.random.seed(2019)

K, M, N, P = 100, 200, 300, 400
AT = sparse.random(K, M, density=0.001, format='csc')
y_k = sparse.random(M, P, density=0.001, format='csc')
t_k = np.exp(1)
AHAT_inv = sparse.random(M, N, density=0.001, format='csc')
b = sparse.random(N, P, density=0.0001, format='csc')
c = sparse.random(K, P, density=0.001, format='csc')

s_k_control = c - (AT.dot(y_k) + AT.dot(t_k * (AHAT_inv.dot(b))))
print("SK_control 1:", s_k_control.min())

s_k_control = c - AT.dot(y_k + (t_k * (AHAT_inv.dot(b))))
print("SK_control 2:", s_k_control.min())

y_k_1 = y_k + t_k * (AHAT_inv.dot(b))

s_k_control = c - AT.dot(y_k_1)
print("SK_control 3:", s_k_control.min())

печатает такой результат, как 

SK_control 1: -0.6701900742964602
SK_control 2: -0.6701900742964602
SK_control 3: -0.6701900742964602

Для дальнейшего изучения вашей ситуации это будетБудьте очень полезны, если вы можете создать работоспособный воспроизводимый пример, демонстрирующий несоответствие.

Обратите внимание, что в документах есть предупреждение *1019*, что настоятельно не рекомендует применение функций NumPy непосредственно к разреженным матрицам: «, поскольку NumPy может неправильно преобразовывать их для вычислений, что приводит к неожиданным (и неправильным) результатам ».Вместо этого используйте методы разреженной матрицы.Поэтому вместо np.min(s_k_control) используйте s_k_control.min().

Если недоступно, и вы не можете придумать другой метод, рекомендуется преобразовать разреженную матрицу в массив NumPy перед применением любого NumPy.функция.

Эта проблема, по-видимому, не является причиной проблемы в вашем случае, но об этом нужно знать.

...