Matlab与r的左除不同

Question

我有一个矩阵A和一个列向量b，如下所示:

0.4585 0.9135 1.1685  1.4235  1.6785  1.9335  2.1885  2.4435  2.6985  2.9535  3.2085  3.4635  3.7185  3.9735  4.2285  4.4835   4.7385  4.9935   
0.9135 1.0685 1.6235  1.9785  2.3335  2.6885  3.0435  3.3985  3.7535  4.1085  4.4635  4.8185  5.1735  5.5285  5.8835  6.2385   6.5935  6.9485
1.1685 1.6235 1.8785  2.5335  2.9885  3.4435  3.8985  4.3535  4.8085  5.2635  5.7185  6.1735  6.6285  7.0835  7.5385  7.9935   8.4485  8.9035
1.4235 1.9785 2.5335  2.8885  3.6435  4.1985  4.7535  5.3085  5.8635  6.4185  6.9735  7.5285  8.0835  8.6385  9.1935  9.7485  10.3035 10.8585
1.6785 2.3335 2.9885  3.6435  4.0985  4.9535  5.6085  6.2635  6.9185  7.5735  8.2285  8.8835  9.5385 10.1935 10.8485 11.5035  12.1585 12.8135
1.9335 2.6885 3.4435  4.1985  4.9535  5.5085  6.4635  7.2185  7.9735  8.7285  9.4835 10.2385 10.9935 11.7485 12.5035 13.2585  14.0135 14.7685
2.1885 3.0435 3.8985  4.7535  5.6085  6.4635  7.3185  8.1735  9.0285  9.8835 10.7385 11.5935 12.4485 13.3035 14.1585 15.0135  15.8685 16.7235
2.4435 3.3985 4.3535  5.3085  6.2635  7.2185  8.1735  9.1285 10.0835 11.0385 11.9935 12.9485 13.9035 14.8585 15.8135 16.7685  17.7235 18.6785
2.6985 3.7535 4.8085  5.8635  6.9185  7.9735  9.0285 10.0835 11.1385 12.1935 13.2485 14.3035 15.3585 16.4135 17.4685 18.5235  19.5785 20.6335
2.9535 4.1085 5.2635  6.4185  7.5735  8.7285  9.8835 11.0385 12.1935 13.3485 14.5035 15.6585 16.8135 17.9685 19.1235 20.2785  21.4335 22.5885
3.2085 4.4635 5.7185  6.9735  8.2285  9.4835 10.7385 11.9935 13.2485 14.5035 15.7585 17.0135 18.2685 19.5235 20.7785 22.0335  23.2885 24.5435
3.4635 4.8185 6.1735  7.5285  8.8835 10.2385 11.5935 12.9485 14.3035 15.6585 17.0135 18.3685 19.7235 21.0785 22.4335 23.7885  25.1435 26.4985
3.7185 5.1735 6.6285  8.0835  9.5385 10.9935 12.4485 13.9035 15.3585 16.8135 18.2685 19.7235 21.1785 22.6335 24.0885 25.5435  26.9985 28.4535
3.9735 5.5285 7.0835  8.6385 10.1935 11.7485 13.3035 14.8585 16.4135 17.9685 19.5235 21.0785 22.6335 24.1885 25.7435 27.2985  28.8535 30.4085
4.2285 5.8835 7.5385  9.1935 10.8485 12.5035 14.1585 15.8135 17.4685 19.1235 20.7785 22.4335 24.0885 25.7435 27.3985 29.0535  30.7085 32.3635
4.4835 6.2385 7.9935  9.7485 11.5035 13.2585 15.0135 16.7685 18.5235 20.2785 22.0335 23.7885 25.5435 27.2985 29.0535 30.8085  32.5635 34.3185
4.7385 6.5935 8.4485 10.3035 12.1585 14.0135 15.8685 17.7235 19.5785 21.4335 23.2885 25.1435 26.9985 28.8535 30.7085 32.5635  34.4185 36.2735
4.9935 6.9485 8.9035 10.8585 12.8135 14.7685 16.7235 18.6785 20.6335 22.5885 24.5435 26.4985 28.4535 30.4085 32.3635 34.3185  36.2735 38.2285

b

5.97
7.07
8.17
9.27
10.37
11.47
9.57
10.67
11.77
12.87
13.97
15.07
16.17
17.27
18.37
19.47
20.57
21.67

然后在Matlab中用

进行左除

A = importdata('A.txt');
b = importdata('b.txt');
Ab

输出

-15.0000
-15.0000
-15.0000
-15.0000
-15.0000
-15.0000
75.5708
95.1859
135.8990
-66.8754
-85.3158
-2.5355
-23.0411
25.1833
-61.6936
40.9536
-56.1559
32.8247

我也想在R中做同样的事情我试着

A = tseries::read.matrix('A.txt')
b = tseries::read.matrix('b.txt')
# matrix(solve(A, b)) #Error in solve.default(A, b) 
matrix(solve(A, b, tol=9e-20))

,

[,1]
[1,]  -15.000000
[2,]  -15.000000
[3,]  -15.000000
[4,]  -15.000000
[5,]  -15.000000
[6,]  -15.000000
[7,]   22.749715
[8,]   93.524478
[9,]  132.015938
[10,]  -74.173920
[11,] -128.108988
[12,]  105.734883
[13,]   82.000007
[14,]   26.968011
[15,] -112.856319
[16,]    5.686058
[17,]  -23.565352
[18,]  -19.974511

可以发现两个部分的结果是不同的。我们能让R的输出和Matlab的输出一样吗?

Answer 1

tldr;Ax = b有无穷多个解

---

这更像是一个线性代数问题而不是编码问题。由于秩(A)和秩(A|b) (A|b是加增矩阵b作为列向量)等于且小于A的列数/行数(即A不是满秩)，因此存在无穷多个解。

我们可以检查R

qr(A)$rank
#[1] 8
qr(cbind(A, b))$rank
#[1] 8
ncol(A)
#[1] 18

参见这里，这里的数学堆栈交换。

---

为了验证我们可以使用两种解(x1来自Matlab,x2来自R)计算Ax。

all(round(A %*% x1, 3) == round(b, 3))
#[1] TRUE
all(round(A %*% x2, 3) == round(b, 3))
#[1] TRUE

示例数据

x1 <- c(-15, -15, -15, -15, -15, -15, 75.5708, 95.1859, 135.899, -66.8754, 
-85.3158, -2.5355, -23.0411, 25.1833, -61.6936, 40.9536, -56.1559, 
32.8247)
x2 <- c(-15, -15, -15, -15, -15, -15, 22.749715, 93.524478, 132.015938, 
-74.17392, -128.108988, 105.734883, 82.000007, 26.968011, -112.856319, 
5.686058, -23.565352, -19.974511)
A <- structure(
c(0.4585, 0.9135, 1.1685, 1.4235, 1.6785, 1.9335, 2.1885, 
2.4435, 2.6985, 2.9535, 3.2085, 3.4635, 3.7185, 3.9735, 4.2285, 
4.4835, 4.7385, 4.9935, 0.9135, 1.0685, 1.6235, 1.9785, 2.3335, 
2.6885, 3.0435, 3.3985, 3.7535, 4.1085, 4.4635, 4.8185, 5.1735, 
5.5285, 5.8835, 6.2385, 6.5935, 6.9485, 1.1685, 1.6235, 1.8785, 
2.5335, 2.9885, 3.4435, 3.8985, 4.3535, 4.8085, 5.2635, 5.7185, 
6.1735, 6.6285, 7.0835, 7.5385, 7.9935, 8.4485, 8.9035, 1.4235, 
1.9785, 2.5335, 2.8885, 3.6435, 4.1985, 4.7535, 5.3085, 5.8635, 
6.4185, 6.9735, 7.5285, 8.0835, 8.6385, 9.1935, 9.7485, 10.3035, 
10.8585, 1.6785, 2.3335, 2.9885, 3.6435, 4.0985, 4.9535, 5.6085, 
6.2635, 6.9185, 7.5735, 8.2285, 8.8835, 9.5385, 10.1935, 10.8485, 
11.5035, 12.1585, 12.8135, 1.9335, 2.6885, 3.4435, 4.1985, 4.9535, 
5.5085, 6.4635, 7.2185, 7.9735, 8.7285, 9.4835, 10.2385, 10.9935, 
11.7485, 12.5035, 13.2585, 14.0135, 14.7685, 2.1885, 3.0435, 
3.8985, 4.7535, 5.6085, 6.4635, 7.3185, 8.1735, 9.0285, 9.8835, 
10.7385, 11.5935, 12.4485, 13.3035, 14.1585, 15.0135, 15.8685, 
16.7235, 2.4435, 3.3985, 4.3535, 5.3085, 6.2635, 7.2185, 8.1735, 
9.1285, 10.0835, 11.0385, 11.9935, 12.9485, 13.9035, 14.8585, 
15.8135, 16.7685, 17.7235, 18.6785, 2.6985, 3.7535, 4.8085, 5.8635, 
6.9185, 7.9735, 9.0285, 10.0835, 11.1385, 12.1935, 13.2485, 14.3035, 
15.3585, 16.4135, 17.4685, 18.5235, 19.5785, 20.6335, 2.9535, 
4.1085, 5.2635, 6.4185, 7.5735, 8.7285, 9.8835, 11.0385, 12.1935, 
13.3485, 14.5035, 15.6585, 16.8135, 17.9685, 19.1235, 20.2785, 
21.4335, 22.5885, 3.2085, 4.4635, 5.7185, 6.9735, 8.2285, 9.4835, 
10.7385, 11.9935, 13.2485, 14.5035, 15.7585, 17.0135, 18.2685, 
19.5235, 20.7785, 22.0335, 23.2885, 24.5435, 3.4635, 4.8185, 
6.1735, 7.5285, 8.8835, 10.2385, 11.5935, 12.9485, 14.3035, 15.6585, 
17.0135, 18.3685, 19.7235, 21.0785, 22.4335, 23.7885, 25.1435, 
26.4985, 3.7185, 5.1735, 6.6285, 8.0835, 9.5385, 10.9935, 12.4485, 
13.9035, 15.3585, 16.8135, 18.2685, 19.7235, 21.1785, 22.6335, 
24.0885, 25.5435, 26.9985, 28.4535, 3.9735, 5.5285, 7.0835, 8.6385, 
10.1935, 11.7485, 13.3035, 14.8585, 16.4135, 17.9685, 19.5235, 
21.0785, 22.6335, 24.1885, 25.7435, 27.2985, 28.8535, 30.4085, 
4.2285, 5.8835, 7.5385, 9.1935, 10.8485, 12.5035, 14.1585, 15.8135, 
17.4685, 19.1235, 20.7785, 22.4335, 24.0885, 25.7435, 27.3985, 
29.0535, 30.7085, 32.3635, 4.4835, 6.2385, 7.9935, 9.7485, 11.5035, 
13.2585, 15.0135, 16.7685, 18.5235, 20.2785, 22.0335, 23.7885, 
25.5435, 27.2985, 29.0535, 30.8085, 32.5635, 34.3185, 4.7385, 
6.5935, 8.4485, 10.3035, 12.1585, 14.0135, 15.8685, 17.7235, 
19.5785, 21.4335, 23.2885, 25.1435, 26.9985, 28.8535, 30.7085, 
32.5635, 34.4185, 36.2735, 4.9935, 6.9485, 8.9035, 10.8585, 12.8135, 
14.7685, 16.7235, 18.6785, 20.6335, 22.5885, 24.5435, 26.4985, 
28.4535, 30.4085, 32.3635, 34.3185, 36.2735, 38.2285), 
.Dim = c(18L, 18L), .Dimnames = list(NULL, NULL))
b <- c(5.97, 7.07, 8.17, 9.27, 10.37, 11.47, 9.57, 10.67, 11.77, 12.87, 
13.97, 15.07, 16.17, 17.27, 18.37, 19.47, 20.57, 21.67)