Коротко — выбор методов, проверка устойчивости и предобработки зависят от размера и структуры задачи (малые плотные матрицы vs большие разреженные) и от того, есть ли шум в $b$. Ниже с конкретикой.
Рекомендуемые численные методы
- Прямые плотные:
- LU с частичным/полным выбором главного элемента (Gaussian elimination with pivoting). Частичное пикетирование обычно обратимо и обладает свойствами назад-устойчивости.
- QR-разложение (например, через Householder) — лучше для плохо обусловленных систем, более устойчиво к потере значений, даёт решение через $x=R^{-1}{Q}^{T}b$.
- SVD (полное разложение $A=U\Sigma V^T$) — наиболее надёжно для почти вырожденных: даёт псевдообратную $A^{+}=V{\Sigma }^{+}{U}^T$ и позволяет оценить ранг/отрезать малые сингулярные числа (т.к. даёт минимальное 2‑нормное решение).
- Регуляризация (при шумном $b$ или при близкой к сингулярности):
- Тихоновская регуляризация: ${\hat{x}}_{\lambda }=(A^{T}A+\lambda I{)}^{-1}{A}^{T}b$. Выбор $\lambda$ через L‑кривую, GCV или кросс‑валидацию.
- Усечение сингулярных чисел (truncated SVD): отбросить ${\sigma }_i$ ниже порога.
- Итеративные методы (большие разреженные):
- Для SPD: Conjugate Gradient + предобуславливатель (IC — incomplete Cholesky, diagonal scaling).
- Для общего: GMRES, BiCGStab, MINRES с предобуславливателями (ILU, approximate inverse, мультиуровневые).
- Всегда применять предобуславливание и масштабирование для улучшения сходимости.
Как оценить устойчивость/качество решения
- Число обусловленности:
- $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert =\frac{{\sigma }_{\max }(A)}{{\sigma }_{\min }(A)}$. Большое $\kappa$ → сильная чувствительность.
- Остаток:
- $r=b-A\hat{x}$. Малый $\Vert r\Vert$ не гарантирует малую погрешность, если $\kappa (A)$ велико.
- Оценка погрешности (приближённо):
- $\frac{\Vert \hat{x}-x\Vert }{\Vert x\Vert }\lesssim \kappa (A)\frac{\Vert r\Vert }{\Vert b\Vert }$.
- Нормы и назад/вперёд‑ошибки:
- Для назад‑устойчивого метода вычисленное $\hat{x}$ удовлетворяет $(A+\Delta A)\hat{x}=b$ с $\frac{\Vert \Delta A\Vert }{\Vert A\Vert }=O({\epsilon }_{\text{mach}})$. Тогда вперёд‑ошибка масштабируется примерно как $\kappa (A){\epsilon }_{\text{mach}}$.
- Норматив обратной ошибки (backward error) можно оценить:
- $\eta \approx \frac{\Vert r\Vert }{\Vert A\Vert \Vert \hat{x}\Vert +\Vert b\Vert }$.
- Оценка ранга/малых сингулярных чисел:
- Эффективный ранг = число ${\sigma }_i$ таких, что $\frac{{\sigma }_i}{{\sigma }_{\max }}>\text{tol}$ (обычно $\text{tol}$ — несколько ${\epsilon }_{\text{mach}}$ или зависит от уровня шума).
Предобработки и улучшение численной устойчивости
- Масштабирование/эквилинг:
- Диагональное масштабирование: найдите диагональные $D_r,D_c$ и решайте $(D_{r}A{D}_c)y=D_{r}b$, затем $x=D_{c}y$. Уменьшает разброс по порядкам величин. Алгоритм Руиса (iterative equilibration) часто эффективен.
- Нормирование столбцов/строк (деление на их нормы) для улучшения поведения QR/итеративных методов.
- Ранж‑разрушающие разложения:
- Rank‑revealing QR (RRQR) — быстрый тест на линейную зависимость столбцов.
- SVD — наиболее информативен для определения малых сингулярных чисел и построения псевдообратной.
- Удаление/регуляризация плохо обусловленных направлений:
- Усечение малых ${\sigma }_i$ или Тихоновская регуляризация.
- Предобуславливатели для итеративных методов:
- Простые: Jacobi (диагональ), SSOR.
- Сложные: ILU, incomplete Cholesky, approximate inverse, AMG.
- Пивотирование при LU и QR для защиты от разложения по маленьким элементам.
- Итеративное уточнение (iterative refinement):
- После LU/QR вычислить $r=b-A\hat{x}$, решить систему для ошибки и обновить $\hat{x}$. Лучше выполнять вычисление остатка в удвоенной точности, если доступно.
Практические рекомендации
- Малые плотные матрицы (н ≤ несколько тысяч): используйте SVD для диагностики/регуляризации; QR для решения, если SVD слишком дорого.
- Большие разреженные: итеративные методы с хорошим предобуславливателем и эквилингом; если задача почти сингулярна, добавляйте регуляризацию.
- Если данные шумные — предпочитайте регуляризацию (Тихонов, усечение).
- Оценивайте $\kappa (A)$ (через SVD или оценитель обусловленности) и смотрите на распределение сингулярных чисел; маленькие ${\sigma }_{\min }$ — сигнал к регуляризации или сокращению ранга.
- При подозрении на потерю значащих цифр — используйте повышенную точность для остатка/уточнения.
Короткая шпаргалка формул
- $\kappa (A)=\frac{{\sigma }_{\max }}{{\sigma }_{\min }}$.
- Остаток: $r=b-A\hat{x}$.
- Тихонов: ${\hat{x}}_{\lambda }=(A^{T}A+\lambda I{)}^{-1}{A}^{T}b$.
- Псевдообратная (SVD): $A^{+}=V{\Sigma }^{+}{U}^T$.
Если нужно — могу предложить конкретный алгоритм/код (LU+итеративное уточнение, либо SVD/Tikhonov) для вашей матрицы с учётом размера, плотности и уровня шума.