不同的结果分布

Question

我用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))

tuk是标量，其他变量都是稀疏矩阵（csc\u矩阵）

如果我不是完全错的话，由于点积（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？你知道吗

Answer 1

好的，我找到解决办法了。你知道吗

正确的方法是完全避免计算逆。你知道吗

如果我们有一个线性问题Ax = b，解决方法是：x = A^(-1)*b

例如，Python已经有了计算这个解决方案的函数scipy.sparse.linalg公司斯索尔夫先生。你知道吗

如果我需要计算(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次迭代中收敛。你知道吗

Answer 2

所有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

对我们进一步调查你的情况，如果你 可以产生一个可运行的，可复制的例子来证明这种差异。你知道吗

---

请注意，有一个warning in the docs强烈反对将NumPy函数直接应用于稀疏矩阵，“因为NumPy可能无法正确地将它们转换为计算，从而导致意外（和不正确的）结果”。可用时使用稀疏矩阵方法。因此，使用s_k_control.min()代替np.min(s_k_control)。你知道吗

如果不可用，并且无法设计其他方法，建议在应用任何NumPy函数之前将稀疏矩阵转换为NumPy数组。你知道吗

在您的案例中，这个问题似乎不是问题的原因，但需要注意。你知道吗