scipy 不同的结果分配规律

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

t_k 是标量，其他变量都是稀疏矩阵（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. 使用来自约束的 scipy.sparse.linalg.spsolve 计算 y_next：

   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

您是正确的，所有 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 可能无法正确转换它们以进行计算，从而导致意外（且不正确）结果”。如果可用，请使用稀疏矩阵方法。因此，不要使用np.min(s_k_control)，而是使用s_k_control.min()。

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

在您的情况下，这个问题似乎不是问题的原因，但需要注意。

Answer 2

好的，我找到了解决办法。

正确的方法是完全避免计算逆。

如果我们有一个线性问题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 次迭代内收敛。