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

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]

Answer 1

[Ameer Hamza]

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