-4

I have the following 11x11 linear system of equations Ax = b with:

A = {
{1.0000000000000000, 8.0000000000000000, 6.0000000000000000, 12.0000000000000000, 24.0000000000000000, 24.0000000000000000, 8.0000000000000000, 6.0000000000000000, 24.0000000000000000, 24.0000000000000000, 24.0000000000000000},
{4.5999999999999996, 41.8531411531233601, 33.0479488942856037, 87.8349057232554173, 149.3783917109033439, 195.3689938163366833, 121.0451669808013690, 48.8422484540841708, 223.6406089026404516, 851.8470736603384239, 269.3015780207464900},
{21.1599999999999966, 218.9606780479085160, 182.0278210198854936, 642.9142219510971472, 929.7459962556697519, 1590.3768227003254196, 1831.4915561762611560, 397.5942056750813549, 2083.9634145976574473, 30235.1432043200838962, 3021.8058301860087340},
{97.3359999999999701, 1145.5240206653393216, 1002.6076877338904296, 4705.8591727678940515, 5786.8317341801457587, 12946.2633183243797248, 27711.6501551604087581, 3236.5658295810949312, 19419.1186238102454809, 1073154.9275125553831458, 33907.3782725576675148},
{447.7455999999998539, 5992.9723163999815370, 5522.3546042079124163, 34444.8913989153879811, 36017.8173980603314703, 105387.4349242659372976, 419295.1650431178859435, 26346.8587310664843244, 180954.3130575636751018, 38090161.8577392920851707, 380471.2698060897528194},
{0.0000000000000000, 34.2801357124991952, 168.4702728821191613, 2101.6181209908259007, 1236.1435394200643714, 6631.0420254749351443, 38374.2674650820554234, 4069.0485156323466072, 28291.8793721561523853, 7044717.1197200166061521, 60211.4334496619121637},
{2059.6297599999993508, 31353.0895356311411888, 30417.0821226643129194, 252121.9823892920394428, 224178.4848274685500655, 857893.2134182706940919, 6344206.6583608603104949, 214473.3033545676735230, 1686197.1981563565786928, 1351958038.0734937191009521, 4269229.7229307144880295},
{0.0000000000000000, 179.3414198404317403, 927.9328280691040618, 15382.9524602928686363, 7693.8805767663707229, 53979.1670196200575447, 580627.4516345988959074, 33123.5797620395824197, 263633.8804078772664070, 250042569.2999326586723328, 675626.4184535464737564},
{0.0000000000000000, 938.2502198978935439, 5111.0461132262771571, 112596.6815912620077142, 47887.4794405465727323, 439410.6478194649680518, 8785268.3545934017747641, 269638.3520710353623144, 2456635.0642409822903574, 8874917956.1941699981689453, 7581135.8600852200761437},
{0.0000000000000000, 938.2502198978935439, 0.0000000000000000, 56298.3407956310038571, 23943.7397202732863661, 319571.3802323381532915, 8785268.3545934017747641, 0.0000000000000000, 269630.6777825467870571, 3293783983.7421655654907227, 1735440.7390556528698653},
{0.0000000000000000, 70.9608494071368625, 1546.2151390406352220, 34063.2210755480555235, 13279.8613116998949408, 129911.1650312914862297, 2657756.2850107550621033, 183537.2854802548536099, 1654054.3836708476301283, 5487391301.6329326629638672, 5049794.3807012736797333}
};


b = {1, 6.167551546217714, 39.66265463865314, 267.9960092725794, 1918.2310370808632, 137.49061855461255, 14662.396462231256, 1216.4598834815756, 11424.520672986631, 3808.17355766221, 6082.299417407878};

The matrix is clearly ill-conditioned, although the correct solution can be found with mathematica:

x = {0.0775277, 0.00771443, 0.087553, 0.0208838, 8.47931*1e-7, 0.00197285, 0.0000611365, 0.00187375, 0.000283606, 3.82771*1e-9, 0.000788588};

I now want to solve the system using this and many other similar matrices inside a C program. I have tried almost every lapack function for solving a linear system of equations, in particular:

  • dgesv
  • dsgesv
  • dgels
  • dgelss
  • dgelsy

but they all give severely wrong results.

At this point I don't expect to have any typo / mistake from a programming point of view, since trying with well-conditioned matrices I get correct results. I guess it's something conceptually or maybe I have to use other tools.. Is there anything I can do to find get the correct solution with some routine from mathematical libraries?

kangshiyin
  • 9,681
  • 1
  • 17
  • 29
user1584773
  • 699
  • 7
  • 19
  • 2
    Posting the values in this case is completely unnecessary, instead you should post how you **tried** to use *lapack* to solve the problem. Also, if you know that the result is wrong it has to be because you can solve this with a different method, why do you want it in [tag:c]? – Iharob Al Asimi Jun 14 '16 at 23:34
  • 1
    You are asking a math question not a programming question. – kangshiyin Jun 15 '16 at 05:03

1 Answers1

1

Solving ill-coditioned linear equations is generally hard. At least you could not use those one-step LAPACK APIs to get an answer with satisfied numerical error.

As a good start, you could use truncated SVD method to get a more numerically stable result.

https://en.m.wikipedia.org/wiki/Linear_least_squares_(mathematics)

This method is the most computationally intensive, but is particularly useful if the normal equations matrix, XTX, is very ill-conditioned (i.e. if its condition number multiplied by the machine's relative round-off error is appreciably large). In that case, including the smallest singular values in the inversion merely adds numerical noise to the solution. This can be cured with the truncated SVD approach, giving a more stable and exact answer, by explicitly setting to zero all singular values below a certain threshold and so ignoring them, a process closely related to factor analysis.

More effective methods may involve making the matrix well-conditioned before solving by finding a pre conditioning matrix. You need to have some understanding on the structure of the original matrix. You could find some more ideas in the following discussion.

https://www.researchgate.net/post/How_can_I_solve_an_ill-conditioned_linear_system_of_equations

kangshiyin
  • 9,681
  • 1
  • 17
  • 29