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)的矢量
• 第一qazxsw POI（qazxsw POI）是零qazxsw POI的矢量
• y_prev是标量，但为了找到它，我要计算矢量y_0，然后从(355,1)占用最小的项目，并通过一些因素t乘以（通常t*或类似）

我曾尝试：

1. 分析计算t*（通过插入与0 < beta < 1约束下式： 0.9
2. 计算t*使用scipy.sparse.linalg.spsolve从约束： s_next = 0 （实际上t* = (c - A.T * y_prev)/(A.T(AHA.T)^(-1) * b) ，因为y_next因此A.T * y_next = c 不是二次） 然后计算A * A.T * y_next = A * c是这样的： A

这两种方法都给出了预期的结果。

编辑2：

看来，我有计算（AHA.T）的反问题。当我测试它（AHA.T）*（AHA.T）^ - 1我不明白单位矩阵，但事情完全随机的：

A.T

逆本身看起来是这样的：

t*

有没有一种可能，以避免使用逆仍然计算的t_k？

Answer 1

你是正确的，所有3个版本，应计算相同的结果。可能会有一些变化，因为浮点运算既不是联想也不分配：

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

夫妇，与由大量的乘法：

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

打印结果如

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

对于我们进一步调查您的情况，那将是非常有益的，如果你能生产出演示的差异在可运行，可再现的例子。

---

请注意，有强烈不鼓励直接申请NumPy的功能稀疏矩阵“因为NumPy的可能不正确地转换他们的计算，从而导致意外（和不正确的）结果”一

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

。

如果没有，你不能设计出不同的方法，建议你采用任何与NumPy函数之前稀疏矩阵转换为NumPy的阵列。

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

Answer 2

好吧，我已经找到了解决办法。

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

如果我们有一个线性问题warning in the docs，解决的办法是：np.min(s_k_control)

Python有已经计算出该解决方案，例如scipy.sparse.linalg.spsolve的功能。

所以，如果我需要计算s_k_control.min()，我只需要调用Ax = b。

我对于t计算如下：x = A^(-1)*b

这给了我一个稳定的解决方案和20次迭代内我整个算法收敛。