Numpy与mldivide，"\\" matlab算子的比较

Question

在matlab中，A \ B会给出一个特殊解，而numpy.linalg.lstsq则不会。

A = [1 2 0; 0 4 3];
b = [8; 18];
c_mldivide = A \ b
c_mldivide =

                 0
                 4
  0.66666666666667
 c_lstsq = np.linalg.lstsq([[1 ,2, 0],[0, 4, 3]],[[8],[18]])
 print c_lstsq
 c_lstsq = (array([[ 0.91803279],
                   [ 3.54098361],
                   [ 1.27868852]]), array([], dtype=float64), 2, array([ 5.27316304,1.48113184]))

1. mldivide在matlab中如何提供特殊解？
2. 这个解对于实现计算精度有用吗？
3. 为什么这个解很特殊，你如何在numpy中实现它？

- Schrodinger

---

“特殊解”是什么意思？Python的解（[0.918 3.541 1.279]）也是正确的解。您有3个未知数的2个方程，因此没有唯一的解。解为任何实数s的[-1 9/2 0]+s*[3/2 -3/4 1]。将s=2/3用于Matlab解，将s=1.27868852用于Python解。 - David

Octave中的A\b解法与numpy相同。MATLAB文档建议使用计算成本更高的方法pinv(A)*B。在Octave中，这也会产生相同的结果。numpy也有pinv，结果相同。 - hpaulj

这个问题没有最小二乘解——存在无限数量的精确解。如果你的代码中有其他内容依赖于获得“正确”的解决方案，那么你需要指定另一个条件或者你的算法存在问题。 - David

@David 当存在无限数量的精确解时，lstsq 返回具有最小范数的解；即，具有变量平方和最小的解。这是一种可以描述为“最小二乘”的东西，尽管“最小范数”是一个更可取的术语，以避免混淆。 - user3717023

@Bookend 我不完全同意。最小二乘解意味着在最小二乘意义下，近似解中的残差（误差）被最小化，即真实解与近似解之间的距离被最小化。选择范数最小的解是确定要选择哪个解的一种方法，就像Matlab试图尽可能使向量的许多元素为零一样。因此，我认为您的方法提供了一个最小二乘解，但使用了不同的选择返回哪个解的方式。 - David

Answer 1

对于您这样的欠定系统（秩小于变量数），“mldivide”返回尽可能多的零值解。哪个变量将被设置为零是任意选择的。

相比之下，“lstsq”方法在这种情况下返回最小范数的解：也就是说，在无限精确解族中，它会选择具有变量平方和最小的解之一。

因此，Matlab的“特殊”解有些随意：在此问题中可以将三个变量中的任何一个设为零。NumPy给出的解实际上更加特殊：存在唯一的最小范数解。

哪种解决方案更适合您的目的取决于您的目的是什么。解决方案的非唯一性通常是重新思考方程式方法的原因。但既然您要求，这里提供了生成类似Matlab解的NumPy代码。

import numpy as np
from itertools import combinations
A = np.matrix([[1 ,2, 0],[0, 4, 3]])
b = np.matrix([[8],[18]])

num_vars = A.shape[1]
rank = np.linalg.matrix_rank(A)
if rank == num_vars:              
    sol = np.linalg.lstsq(A, b)[0]    # not under-determined
else:
    for nz in combinations(range(num_vars), rank):    # the variables not set to zero
        try: 
            sol = np.zeros((num_vars, 1))  
            sol[nz, :] = np.asarray(np.linalg.solve(A[:, nz], b))
            print(sol)
        except np.linalg.LinAlgError:     
            pass                    # picked bad variables, can't solve

以您的示例为例，它输出了三个“特殊”的解决方案，其中最后一个是Matlab选择的。

[[-1. ]
 [ 4.5]
 [ 0. ]]

[[ 8.]
 [ 0.]
 [ 6.]]

[[ 0.        ]
 [ 4.        ]
 [ 0.66666667]] 

- user3717023

---

似乎 np.linalg.solve(mA, vB) 只适用于方阵。是否有适用于任何大小矩阵的求解器（当然，该矩阵与 vB 具有相同的行数）？ - Royi