Метод решения некоторой вычислительной задачи (например, решения линейной системы) называют точным, если точное решение задачи может быть получено за конечное число операций (без учёта вычислительной погрешности).
Метод решения вычислительной задачи называют итерационным, если он заключается в вычислении последовательности приближений к точному ответу xk ⟶ x при k → ∞.
На практике применяют и те, и другие методы: учитывая наличие вычислительной погрешности, точные методы не всегда дают лучший результат, чем итерационные.
Есть задачи, точного (прямого) метода решения которых не существует — например, задача трёх тел решается только итерационно.
(про них подробнее — на 4 курсе)
Для системы Ax = b решение определяется с помощью итерации:
A = L* + U
L*x(k + 1) = b − Ux(k)
x(0) выбирается произвольно, метод сходится к решению если
A — симметричная, положительно определённая или с диагональным преобладанием.
Пара из собственного значения λk и собственного вектора vk матрицы A — такие комплексное число и вектор, что:
A · vk = λkvk
Любая матрица из 𝐌n имеет n собственных значений (не обязательно различных).
По определению, (A − λkE)·vk = 0, т.е. столбцы матрицы A − λkE линейно зависимы с коэффициентами vk. Поэтому собственные значения можно найти как корни λk многочлена n-й степени:
$ \big| A - λE \big| = 0 $
Теорема. Для матриц из пространства 𝐌n, n ⩾ 5, не существует точного метода нахождения всех собственных значений.
Идея доказательства:
⟹ нужно использовать итерационные методы.
«Где искать собственные значения?»
Обозначим:
Теорема. (Гершгорин) Для любой комплексной матрицы A ∈ 𝐌n собственные значения (на комплексной плоскости) лежат в кругах с центрами в диагональных элементах и радиусами — сумма модулей внедиагональных элементов строк, соответственно:
$$ σ(A) ⊂ \bigcup\limits_{i=1}^{n} \{ z : | z - a_{ii} | ⩽ R'_i(A) \} $$
Если же k кругов образуют связную область (например, диагональные элементы — действительные числа aii ∈ ℝ), то в такой области находится в точности k собственных значений матрицы.
Для матрицы A ∈ 𝐌n, диагонализируемой с использованием невырожденного преобразования C (A = CΛC−1, Λ = diag(λ1, λ2, …, λn)), ошибка при погрешности во входных данных E выражается следующим образом:
$$ | \hat{λ} - λ | ⩽ κ_∞(C)‖E‖_∞, $$
где $\hat{λ}$ — собственное значение матрицы с погрешностью A + E, λ — собственное значение исходной матрицы A.
Пусть A ∈ 𝐌n имеет полную систему ортонормированных собственных векторов ei. Собственные числа матрицы:
|λ1|>|λ2|⩾|λ3|⩾…⩾|λn|,
т.е. максимальное по модулю собственное значение — не кратное.
Тогда для любого начального вектора x(0), не перпендикулярного собственному вектору e1, следующий итерационный процесс сходится λ1, а вектора x(k) сходятся к вектору, коллинеарному собственному вектору e1:
$$ x^{(k+1)} = Ax^{(k)}, \qquad λ^{(k)}_1 = \dfrac{(x^{(k+1)}, x^{(k)})}{(x^{(k)}, x^{(k)})} $$
Критерий остановки — ‖x(k + 1) − x(k)‖<ε
Матрица B подобна A, если существует такая C, что:
A = C · B · C−1
Подобные матрицы имеют один и тот же набор собственных значений.
Матрица B унитарно подобна A, если C — унитарная, т.е. $C^{-1} = C^* $ (в вещественном случае — ортогональная)
Матрица A симметрична
Строится последовательность A = A0, A1, A2, …
ортогонально подобных матриц Ak = OkAk − 1Ok*,
т.ч. Σ(Ak)<Σ(Ak − 1) и limΣ(Ak)→0.
В методе Якоби Ok = Tij(φk), φk подбирается так, чтобы обнулить aij.
Для данных i, j положим x = −2aij, y = aii − ajj, если aii ≠ ajj (в~противном~случае~положим~$φ = \frac{π}{4}$).
Элементы матрицы Tij(φk) вычисляются по формулам:
$$ \cos φ = \sqrt{\dfrac{1}{2} \left( 1 + \dfrac{y}{\sqrt{x^2 + y^2}} \right)} \qquad \sin φ = \dfrac{\sgn(xy) | x | }{2 \cos φ \sqrt{x^2 + y^2}} $$
$(i, j)^{(k)} = \arg \max\limits_{l≠m} |a^{(k-1)}_{l,m} | $
$(i, j) = \left((1, 2), (1,3), \dots, (1,n), (2,3), (2,4), \dots, (n-1, n), {\mathbf\circlearrowright} \right)$
$i = \arg \max\limits_l \sum_{j≠l} |a^{(k-1)}_{lj}|^2$
$j = \arg \max\limits_m |a^{(k-1)}_{im}|$
Ряд способов нахождения собственных значений связан с декомпозицией матрицы на каждом шаге итеративного алгоритма.
В общем случае сложность декомпозиции матрицы — n3, поэтому такие алгоритмы никогда не применяют к матрицам произвольного вида.
Унитарным подобием можно привести произвольную матрицу к почти треугольному виду:
$$
\begin{pmatrix}
* & * & * & … & * & * \\
* & * & * & … & * & * \\
& * & * & … & * & * \\
& & * & \ddots & \vdots & \vdots \\
& & & … & * & *
\end{pmatrix}
$$
Для симметричной матрицы приведение к почти треугольному виду унитарным подобием методами вращений и отражений даёт трёхдиагональную матрицу.
Метод вращений основан на матрицах вращения Tij, которые преобразуют «поддиагональный» вектор ai…n, i − 1(i − 1) в (‖ai + 1…n, i(i − 1)‖,0, …, 0). Такие матрицы можно получить аналогично обычному методу вращений (с той поправкой, что норма собирается не в диагональном, а в поддиагональном элементе; диагональный элемент в определении углов вращения не участвует).
$$\large A' = \underbrace{(T_{i,n}·T_{i,n-1}·…·T_{i,i+1}·A)}_{\small\text{Почти как для приведения матрицы к треугольному виду}} · T^t_{i,i+1} · T^t_{i,i+2} · … · T^t_{i,n-1}·T^t_{i,n} $$
Умножение на матрицу вращения Ti, j справа изменяет i, j столбцы матрицы, а слева — i, j строки матрицы.
(поэтому к треугольному виду привести унитарным подобием произвольную матрицу нельзя)
Метод отражений основан на матрицах отражения $U(x_i) = I - 2 x_i x^{* }_i $.
В отличие от матриц вращения, они являются самосопряжёнными (симметричными в ℝ).
Умножение строки на U(xi) справа тоже выражается через скалярное произведение:
yt · U(xi)=yt − 2xi*(y, xi)
Вектор xi строится как для приведения матрицы к треугольному виду методом отражений с единственным отличием — он строится для части столбца матрицы под диагональю, не включая диагональный элемент.
Шаг приведения к почти треугольному виду методом отражений:
$$\large A' = \underbrace{(U(x_i) · A)} · U(x_i) $$
I. Матрицы приводятся к почти треугольному или трёхдиагональному виду ортогональным подобием (чтобы сохранить собственные значения).
Построение разложения оптимизируется с учётом известной структуры матрицы (заведомо нулевые элементы исключаются из формул).
В книге Богачёва показано, что полученная матрица ортогонально подобна исходной, и последовательность Ak сходится к диагональной матрице, элементами которой являются собственные значения матрицы A.
При реализации QR, LR-алгоритмах и методе Холецкого обязательно использовать один из методов ускорения сходимости (на ваш выбор), в противном случае для многих матриц метод оказывается неприменим.
Используется два метода ускорения сходимости — исчерпывание матрицы и сдвиги, каждый даёт свой критерий остановки.
Пусть ε — требуемая точность и на некотором шаге k один из элементов под диагональю ai + 1, i такой, что |ai + 1, i|<ε‖Ak‖∞ (‖A‖∞ — максимальная строчная норма).
Тогда можно заменить его на 0 и запускать метод незавимимо для двух подматриц (a11…aii и ai + 1, i + 1…ann):
$$
\begin{pmatrix}
a_{11} & a_{12} & … & a_{1,i} & a_{1, i+1} & … & a_{1, n-1} & a_{1, n} \\
a_{21} & a_{22} & … & a_{2,i} & a_{2, i+1} & … & a_{2, n-1} & a_{2, n} \\
& a_{32} & \ddots & \vdots & \vdots & … & \vdots & \vdots \\
& & \ddots & a_{i,i} & a{i, i+1} & … & a_{i, n-1} & a_{i, n} \\
& & & 0 & a_{i+1, i+1} & … & … & … \\
& & & & a_{i+2, i+1} & … & … & … \\
& & & & & & \ddots & \vdots \\
& & & & & & a_{n-1, n-1} & a_{n-1, n} \\
& & & & & & a_{n, n-1} & a_{n, n}
\end{pmatrix}
$$
Для уменьшения вычислительных затрат норму A вычисляют только на первом шаге, далее она изменяется незначительно.
Для матрицы размером 2 собственные значения находятся в явном виде как решение квадратного уравнения.
На каждом шаге k определяется значение сдвига sk (в нашем случае можно взять sk = ann).
Разложение выполняется для матрицы Ak − skI (значение сдвига вычитается из диагональных элементов).
После выполнения шага метода, к матрица Ak + 1 прибавляется skI.
Итерация продолжается пока модуль последнего элемента под диагональю
не станет меньше ε‖A‖∞:
|an, n − 1|<ε‖Ak‖∞
В этом случае считаем, что ann — это одно из искомых собственных значений и переходим к рассмотрению меньшей подматрицы a11…an − 1, n − 1.