当A是奇异平方和全(稠密)矩阵时，MATLAB A\B是如何计算的

Question

我正在研究像A*X = B这样的线性方程。我得到了一个类似A = magic(4)和B = [1;3;2;4]的矩阵。我使用了几种方法来解决这个问题，还包括MATLABAB。

令人恼火的是，MATLABAB的结果与NumPy、SciPy和LAPACK(LAPACKE_dgetrf、LAPACKE_dgels、LAPACKE_dgelsd)的结果不同。

结果如下。

MATLABAB(Ubuntu MATLAB 2021a和MATLAB Online)

>> A = magic(4)

A =
16     2     3    13
5    11    10     8
9     7     6    12
4    14    15     1
>> B = [1;3;2;4]
B =
1
3
2
4
>> A  B
RCOND =  4.625929e-18。  
ans =
-0.0098
0.0735
0.1985
0.0319
>> lsqr(A, B)
ans =
0.0110
0.1360
0.1360
0.0110
>> x = pinv(A)*B
x =
0.0110
0.1360
0.1360
0.0110
>> [L, U, P] = lu(A) %lu decomposition
L =
1.0000         0         0         0
0.2500    1.0000         0         0
0.5625    0.4352    1.0000         0
0.3125    0.7685    1.0000    1.0000

U =
16.0000    2.0000    3.0000   13.0000
0   13.5000   14.2500   -2.2500
0         0   -1.8889    5.6667
0         0         0    0.0000

P =
1     0     0     0
0     0     0     1
0     0     1     0
0     1     0     0
>> y = L  (P * B) % LU tril solver
y =
1.0000
3.7500
-0.1944
0.0000
>> x = U  y % LU triu solver
x =
-0.0098
0.0735
0.1985
0.0319

Pythonnumpy.linalg.solve(a, b)

import numpy as np
a = np.array([[16., 2., 3., 13.], [5., 11., 10., 8.], [9., 7., 6., 12.], [4., 14., 15., 1.]]);
b = np.array([[1.],[3.],[2.],[4]]);
x = np.linalg.solve(a, b)
'''
x = 
array([[ 0.10799632],
[ 0.42693015],
[-0.15487132],
[-0.0859375 ]])
'''

Pythonscipy.linalg.solve(a, b)

from scipy import linalg
x1 = linalg.solve(a, b)
'''
x1 = 
array([[ 0.10799632],
[ 0.42693015],
[-0.15487132],
[-0.0859375 ]])
'''

Octave在线

octave:1> A = magic(4)
A =
16    2    3   13
5   11   10    8
9    7    6   12
4   14   15    1
octave:2> B = [1;3;2;4]
B =
1
3
2
4
octave:3> AB
warning: matrix singular to machine precision, rcond = 1.30614e-17
ans =
0.011029
0.136029
0.136029
0.011029

LapackLAPACKE_degelsd

a = 
Matrix(Row = 4, Col = 4, Major = ColMajor)
16             2             3            13
5            11            10             8
9             7             6            12
4            14            15             1
B = 
Vector(Size = 4)
1
3
2
4
info = 0
X = 
Vector(Size = 4)
0.107996
0.42693
-0.154871
-0.0859375

我的问题是MATLABA B到底计算了什么。我知道pinv(A) * B计算最小范数最小二乘解，np.linalg.solve和scipy.linalg.solve也是如此。

MATLAB参考页面上写着：

• 如果A是一个平方矩阵，那么AB大致等于inv(A)*B，但MATLAB对AB的处理方式不同，而且更稳健。

• 如果A的秩小于A中的列数，则x = AB不一定是最小范数解。可以使用x = lsqminnorm(A,B)或x = pinv(A)*B计算最小范数最小二乘解。

我知道当A几乎为奇异时，解不是唯一的。有许多可行的解决方案(如numpy、倍频程和LAPACK。它们通过QR、SVD(pinv)最小化||AX - B||的范数来解决这个问题)。然而，是什么让matlab得到了它的解决方案，它的魔力是什么？

所以，我想问，当A是奇异平方矩阵时，AB的运算到底是什么，AB的结果是什么。

p.S.

我知道，当线性方程的系数矩阵是奇异的或近似奇异的，也就是病态的时，并没有唯一的解。根据(mldive)的参考页面，x = AB不一定是最小范数解。我想知道AB的基本细节。如果它不一定是最小范数解，那么呢？正如我所知，scipy、numpy、LAPACK和Octave得到了一个最小范数解。这个问题还有另一个我不知道的公式吗？

Answer 1

有几种方法可以尝试解决这个问题，但我建议的一种方法是，当A条件不好时，尝试解决Ax = B是没有意义的：

以1D为例。什么是x而不是0.x = b？

我们知道，如果b = 0，这个问题有无限多个解，如果b~=0，这个问题没有解。MATLAB对此怎么说：

0
ans =
NaN

01
ans =
Inf

不是很有说服力。lsqminnorm:怎么样

lsqminnorm(0,0)
ans =
0
lsqminnorm(0,1)
ans =
0

看到了吗？只有当存在无穷多个解时，lsqminnorm才会给你一个有效的答案。当没有时，它只是试图最小化|| 0.x - b ||。这有道理吗？您对它返回x = 0作为0.x = 1的解决方案感到满意吗？

回到你的问题上来。让我们为B:尝试一个不同的值

B2 = [1;1;2;1];
out2 = AB2;
A*out2
out_minnorm = lsqminnorm(A,B2);
A*out_minnorm
out_pinv = pinv(A)*B2;
A*out_pinv
ans =
0
-1.0000
0
-2.6250

ans =
1.1500
1.4500
1.5500
0.8500

ans =
1.1500
1.4500
1.5500
0.8500

x的这些值中有任何一个可以接受吗？尽管最后两个结果是最小范数解，但这并没有多大意义。

所以，我对你的问题的看法是检查rcond是否接近机器精度，因为在我看来，不返回任何结果比返回没有意义的东西要好：

if rcond(A) < 1e-12
error('A is singular, the problem seems to be ill-conditioned');
else
out = BA;
end