为什么形状 (3,3) 矩阵乘以 (3,) 向量加上 (3,) 向量在 numpy 中形成 (1,3) 形状？

Question

我的代码如下，进行雅可比分解来求解矩阵

def iterative_solver(M, b, Vinit, kmax, epsmax):
    """
    Provides an approximation of M^{-1}b using the proposed iterative method
    ----------   
    parameters:
    M      : matrix in the equation MV = b to solve (numpy array of size N,N)
    b      : vector on the right-hand-side of the equation MV = b (numpy array of size N)
    Vinit  : initial vector of the iterative method (numpy array of size N)
    kmax   : maximum number of iterations, stops if k reaches kmax (integer)
    epsmax : tolerance on the residual, stops if eps reaches epsmax (float)
    
    returns:
    V   : resulting vector at the end of the iterations (numpy array of size N)
    eps : vector composed of the residuals eps at every iteration (numpy array of size k)
    k   : number of iterations performed before the algorithm stops (integer)
    """
    eps = np.inf
    N = len(M)
    d = list()
    for i in range(N):
        d.append(M[i,i])
    D = np.diag(d)
    R = D - M
    Dinv = np.linalg.inv(D)
    V = Vinit
    k =0
    print(f"V: {np.shape(V)} M: {np.shape(M)} b: {np.shape(b)}")
    
    
    while k <= kmax and eps >= epsmax:
        A = np.matmul((-M), V)+ b
        print(f"np.shape(A): {np.shape(A)}")
        V = V + np.matmul(Dinv,A)
        print(f"np.shape(V): {np.shape(V)}")
        #A = np.matmul(R,V)
        #V = np.matmul(Dinv, A + b)
        k += 1
        eps = np.linalg.norm(M*V.T-b)
    return V, eps, k

M = np.mat([[3,1,0],[-1, -5, 1], [0,2,4]])
b = np.array([4,-5,6]).T
Vinit = np.random.rand(3).T
kmax = 100
epsmax = 1e-10

print(f"solved: {iterative_solver(M, b, Vinit, kmax, epsmax)}")

，但是它给出了输出

V: (3,) M: (3, 3) b: (3,) np.shape(A): (1, 3)

这很奇怪，因为形状 (3,3) 矩阵乘以 (3,) 向量加上 (3,) 向量在 numpy 中形成 (1,3) 形状？

Answer 1

形状业务隐藏在您的代码中，但我认为这就是您所指的：

V: (3,) M: (3, 3) b: (3,) 
A = np.matmul((-M), V)+ b

(3,3) 与 (3,) 的 matmul 生成 (3,)。添加 (3,) 应该保留 (3,)

但随后您将 V 更改为 V + np.matmul(Dinv,A)。你打印了它的形状，但我看不到。

In [224]: V = np.arange(3)
     ...: M = np.ones((3, 3))
     ...: b = np.ones(3)
In [225]: V, M, b
Out[225]: 
(array([0, 1, 2]),
 array([[1., 1., 1.],
        [1., 1., 1.],
        [1., 1., 1.]]),
 array([1., 1., 1.]))
In [226]: np.matmul((-M), V) + b
Out[226]: array([-2., -2., -2.])

使用您的代码：

In [23]: iterative_solver(M, b, Vinit, kmax, epsmax)
V: (3,) M: (3, 3) b: (3,)
np.shape(A): (1, 3)
Traceback (most recent call last):
  Input In [23] in <module>
    iterative_solver(M, b, Vinit, kmax, epsmax)
  Input In [4] in iterative_solver
    V = V + np.matmul(Dinv, A)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1 is different from 3)

我尝试过：

In [25]: -M @ Vinit
Out[25]: matrix([[-2.66022252,  2.20102655, -2.90845177]])

我花了一点时间才注意到这是np.matrix！根据定义，它是 2d，因此是 (1,3) 形状。

如果我传递一个 np.array 矩阵：

In [28]: iterative_solver(M.A, b, Vinit, kmax, epsmax)
V: (3,) M: (3, 3) b: (3,)
(3,)
np.shape(A): (3,)
...
Out[28]: (array([1., 1., 1.]), 13.856406460551018, 101)

编辑

---

np.matrix 是 ndarray 的子类。它将自己的 * 和 ** 定义为矩阵运算，而不是按元素运算。它强制执行二维形状。如果需要，它会添加一个前导维度（numpy 默认情况下为“C”顺序）

In [150]: M = np.matrix(np.arange(9).reshape(3, 3))
In [151]: v = np.arange(3)

使用自己的 *：

In [152]: M * v
Traceback (most recent call last):
  Input In [152] in <module>
    M * v
  File /usr/local/lib/python3.8/dist-packages/numpy/matrixlib/defmatrix.py:218 in __mul__
    return N.dot(self, asmatrix(other))
  File <__array_function__ internals>:180 in dot
ValueError: shapes (3,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)

(3,) v 有被‘提升’到矩阵(1,3)。

matmul/@ 是最近添加的内容（np.matrix 已经很旧了，即将淘汰）：

In [153]: M @ v
Out[153]: matrix([[ 5, 14, 23]])
In [154]: M.dot(v)
Out[154]: matrix([[ 5, 14, 23]])

看起来这两个函数都执行乘法，就好像它们是>ndarray，并将结果转换为矩阵。

如果两者都是 np.matrix：

In [155]: V = np.matrix(v)
In [156]: V
Out[156]: matrix([[0, 1, 2]])
In [157]: M * V.T
Out[157]: 
matrix([[ 5],
        [14],
        [23]])

[email protected] ] 做同样的事情。