how to reproduce mldivide solver for a linear system with the exact method?

1 次查看（过去 30 天）
JT Yang
JT Yang 2020-4-20
编辑： JT Yang 2020-4-21
Hello, I'm solving a poorly conditioned linear system A*x = b which cames from a discretized poisson like equation. However, the A matrix has a high condition number.
I've been reading some literatures suggesting to use iterative methods such as cgs and gmres, but either could not get a convergence. But I could solve the linear system with matlab mrdivide. Is there anyways to understand the exact method mrdivide used for solving the system and reproduce the same result by the exact method used?
The example linear system is as followed
A = [0.111904781196914 0.0168898905984570 0 0.0168898905984570 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0 0 0 0 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715;
0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0;
0 0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715;
0.0168898905984570 0 0.0168898905984570 0.111904781196914 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0 0 0 0 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570;
0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0.111904781196914 0.0168898905984570 0 0.0168898905984570 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0 0 0 0;
-0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0 0;
0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0;
-0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0.0168898905984570 0 0.0168898905984570 0.111904781196914 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0 0 0 0;
0 0 0 0 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0.111904781196914 0.0168898905984570 0 0.0168898905984570 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715;
0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0;
0 0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715;
0 0 0 0 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0.0168898905984570 0 0.0168898905984570 0.111904781196914 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570;
0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0 0 0 0 0.0168898905984570 -0.0110863047007715 0 -0.0110863047007715 0.111904781196914 0.0168898905984570 0 0.0168898905984570;
-0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570 0;
0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0 0 0 0 -0.0110863047007715 0.0168898905984570 -0.0110863047007715 0 0.0168898905984570 0.111904781196914 0.0168898905984570;
-0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0 0 0 0 -0.0110863047007715 0 -0.0110863047007715 0.0168898905984570 0.0168898905984570 0 0.0168898905984570 0.111904781196914]
b = [-0.0186760600000000 -0.101709940000000 0.0286474500000000 0.0465896100000000 0.0531106000000000 -0.523516250000000 -0.0214476100000000 -0.171042970000000 0.0408699700000000 0.0735594600000000 -0.0196324800000000 0.00217344000000000 0.0158687800000000 0.149093200000000 -0.0893736500000000 0.535486420000000]

请先登录，再回答此问题。

回答（1 个）

Ameer Hamza
Ameer Hamza 2020-4-20
编辑：Ameer Hamza 2020-4-20
I guess your question is related to mldivide, not mrdivide. This shows the top-level details of the algorithm: https://www.mathworks.com/help/matlab/ref/mldivide.html#bt42oms_head. However, as far as I know, the low-level details of the algorithm are owned by MathWorks and not publicly available. They have optimized several aspects of the algorithm, based on the type of the input matrix. However, In general, it can be considered to be solving the following optimization problem
For example, see this
A = [1 3 3 4; 8 6 7 2; 9 5 2 1; 5 4 8 7];
b = [10; 20; 30; 40];
x1 = A\b;
x2 = fmincon(@(x) sum((A*x-b).^2), rand(4,1));
Result:
>> x1
x1 =
5.7594
-5.1085
-0.9467
5.6016
>> x2
x2 =
5.7594
-5.1085
-0.9467
5.6016
But for your matrices A and B, because of high condition number, mldrive must be using some other technique. The solutions are not the same for both methods. I think the solution given by fmincon() seems more reasonable.
A = % 16x16 matrix
B = % 1x16 matrix
x1 = A\b.';
x2 = fmincon(@(x) sum((A*x-b.').^2), rand(16,1));
Result
>> x1
x1 =
1.0e+14 *
-2.4735
2.4735
-2.4735
2.4735
2.4735
-2.4735
2.4735
-2.4735
-2.4735
2.4735
-2.4735
2.4735
2.4735
-2.4735
2.4735
-2.4735
>> x2
x2 =
-0.2313
-0.4032
0.4525
-0.1399
1.5872
-5.1760
0.8600
-1.8094
0.1117
1.2050
-0.1681
-0.4949
-0.8393
1.6084
-1.8407
5.2781
The results given by mldivide are of the order .Also, the solution given by them seem to fit the equation quite well
>> norm(A*x1-b.')
ans =
0.0122
>> norm(A*x2-b.')
ans =
0.0106

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by