I have a matrix A and a column vector b as follows:
- A
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
Then I conduct the left division in Matlab as
A = importdata('A.txt');
b = importdata('b.txt');
A\b
which gives output
-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
I also want to do the same thing in R and I try
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))
which gives
[,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
One can find that two results in parts are different. Can we make R output the same with the Matlab output?
CodePudding user response：
tldr; Ax = b has infinitely many solutions.
This is more of a linear algebra question than a coding question. Since rank(A) and rank(A|b) (A|b being the augmented matrix with b appended as a column vector) are equal and less than the number of columns/rows of A (i.e. A is not full rank), there exist infinitely many solutions.
We can check in R
qr(A)$rank
#[1] 8
qr(cbind(A, b))$rank
#[1] 8
ncol(A)
#[1] 18
See e.g. here, here on the Mathematics Stack Exchange.
To verify we can calculate Ax using both solutions (x1 from Matlab, and x2 from R).
all(round(A %*% x1, 3) == round(b, 3))
#[1] TRUE
all(round(A %*% x2, 3) == round(b, 3))
#[1] TRUE
Sample data
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)