c - 使用 Lapack&co 求解病态线性方程组

Question

我有以下 11x11 线性方程组 Ax = b:

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};

尽管可以使用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};

我现在想在 C 程序中使用这个矩阵和许多其他类似的矩阵来求解系统。我几乎尝试了所有用于求解线性方程组的 lapack 函数，特别是:

• dgesv
• dsgesv
• dgels
• dgelss
• 德格尔西

但它们都给出了严重错误的结果。

此时，从编程的角度来看，我不希望出现任何拼写错误/错误，因为尝试使用条件良好的矩阵我会得到正确的结果。我想这是概念上的问题，或者也许我必须使用其他工具。我可以做些什么来通过数学库中的一些例程找到正确的解决方案吗？

Answer 1

求解不良线性方程通常很困难。至少您无法使用那些一步式 LAPACK API 来获得具有令人满意的数值误差的答案。

作为一个好的开始，您可以使用截断 SVD 方法来获得数值更稳定的结果。

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.

更有效的方法可能包括在求解之前通过找到预处理矩阵来对矩阵进行良好调节。您需要对原始矩阵的结构有一些了解。您可以在下面的讨论中找到更多想法。

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