Предыдущая лекция

  • Прямые методы для LU
  • Сепараторы графов

Основные темы на сегодня

Концепция итерационных методов для линейных систем: - Итерация Ричардсона и её сходимость - Итерация Чебышева

Итерационные методы

  • Если мы хотим достичь сложности \mathcal{O}(N) при решении разреженных линейных систем, то прямые решатели не подходят.

  • Если мы хотим решить частичную задачу на собственные значения, полное разложение на собственные значения слишком затратно.

  • Для обеих задач мы будем использовать итерационные решатели подпространств Крылова, которые рассматривают матрицу как линейный оператор по принципу черного ящика.

Матрица как черный ящик

  • Теперь у нас совершенно другой взгляд на матрицу: матрица теперь линейный оператор, который действует на вектор,
    и это действие может быть вычислено за \mathcal{O}(N) операций.

  • Это единственная информация, которую мы знаем о матрице: произведение матрицы на вектор (матвек)

  • Можем ли мы решать линейные системы, используя только матвеки?

  • Конечно, мы можем умножать на столбцы единичной матрицы и восстановить полную матрицу, но это не то, что нам нужно.

Итерация Ричардсона

Самая простая идея - это “метод простой итерации” или итерация Ричардсона.

Ax = f, \tau (Ax - f) = 0, x - \tau (Ax - f) = x, x_{k+1} = x_k - \tau (Ax_k - f),

где \tau - это итерационный параметр, который всегда можно выбрать так, чтобы метод сходился.

Связь с ОДУ

  • Итерация Ричардсона имеет глубокую связь с Обыкновенными Дифференциальными Уравнениями (ОДУ).

  • Рассмотрим зависящую от времени задачу

\frac{dy}{dt} + A y = f, \quad y(0) = y_0.

  • Тогда y(t) \rightarrow A^{-1} f при t \rightarrow \infty, и схема Эйлера имеет вид

\frac{y_{k+1} - y_k}{\tau} = -A y_k + f.

что приводит к итерации Ричардсона

y_{k+1} = y_k - \tau(Ay_k -f)

Сходимость метода Ричардсона

  • Пусть x_* - решение; введем ошибку e_k = x_{k} - x_*, тогда

e_{k+1} = (I - \tau A) e_k,

следовательно, если \Vert I - \tau A \Vert < 1 в любой норме, итерация сходится.

  • Для симметричного положительно определенного случая всегда возможно выбрать \tau так, чтобы метод сходился.

  • А как насчет несимметричного случая? Ниже будет представлена демонстрация…

Оптимальный выбор параметра

  • Выбор \tau, который минимизирует \|I - \tau A\|_2 для A = A^* > 0 равен (докажите это!)

\tau_\mathrm{opt} = \frac{2}{\lambda_{\min} + \lambda_{\max}}.

где \lambda_{\min} - минимальное собственное значение, а \lambda_{\max} - максимальное собственное значение матрицы A.

  • Таким образом, чтобы найти оптимальный параметр, нам нужно знать границы спектра матрицы A, и мы можем вычислить их, используя степенной метод.

Число обусловленности и скорость сходимости

Даже при оптимальном выборе параметра, ошибка на следующем шаге удовлетворяет

\|e_{k+1}\|_2 \leq q \|e_k\|_2 , \quad\rightarrow \quad \|e_k\|_2 \leq q^{k} \|e_0\|_2,

где

q = \frac{\lambda_{\max} - \lambda_{\min}}{\lambda_{\max} + \lambda_{\min}} = \frac{\mathrm{cond}(A) - 1}{\mathrm{cond}(A)+1},

\mathrm{cond}(A) = \frac{\lambda_{\max}}{\lambda_{\min}} \quad \text{для} \quad A=A^*>0

является числом обусловленности матрицы A.

Давайте проведем демонстрацию…

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg as spla
import scipy
from scipy.sparse import csc_matrix
n = 500
ex = np.ones(n);
A = sp.sparse.spdiags(np.vstack((-ex,  2*ex, -ex)), [-1, 0, 1], n, n, 'csr'); 
rhs = np.ones(n)
ev1, vec = spla.eigsh(A, k=2, which='LA')
ev2, vec = spla.eigsh(A, k=2, which='SA')
lam_max = ev1[0]
lam_min = ev2[0]

tau_opt = 2.0/(lam_max + lam_min)

fig, ax = plt.subplots()
plt.close(fig)

niters = 1000
x = np.zeros(n)
res_richardson = []
for i in range(niters):
    rr = A.dot(x) - rhs
    x = x - tau_opt * rr
    res_richardson.append(np.linalg.norm(rr))
#Convergence of an ordinary Richardson (with optimal parameter)
plt.semilogy(res_richardson)
plt.xlabel("Number of iterations, $k$", fontsize=20)
plt.ylabel("Residual norm, $\|Ax_k - b\|_2$", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
print("Maximum eigenvalue = {}, minimum eigenvalue = {}".format(lam_max, lam_min))
cond_number = lam_max.real / lam_min.real
print("Condition number = {}".format(cond_number))
#print(np.array(res_richardson)[1:] / np.array(res_richardson)[:-1])
print("Theoretical factor: {}".format((cond_number - 1) / (cond_number + 1)))
Maximum eigenvalue = 3.99984271815512, minimum eigenvalue = 3.932084756989659e-05
Condition number = 101723.207035276
Theoretical factor: 0.9999803389964071

  • Таким образом, для плохо обусловленных матриц ошибка метода простой итерации убывает очень медленно.

  • Это еще одна причина, почему число обусловленности так важно:

    • Помимо оценки погрешности решения, оно также дает оценку количества итераций для итерационных методов.

  • Главный вопрос для итерационного метода - как сделать матрицу лучше обусловленной.

Рассмотрим неэрмитову матрицу A

Возможные случаи поведения итерации Ричардсона: - сходимость - расходимость - почти стабильная траектория

Вопрос: как мы можем определить наш случай до запуска итерационного метода?

# B = np.random.randn(2, 2)
B = np.array([[1, 2], [-1, 0]])
# B = np.array([[0, 1], [-1, 0]])
x_true = np.zeros(2)
f = B.dot(x_true)
eigvals = np.linalg.eigvals(B)
print("Spectrum of the matrix = {}".format(eigvals))

# Run Richardson iteration
x = np.array([0, -1])
tau = 1e-2
conv_x = [x]
r = B.dot(x) - f
conv_r = [np.linalg.norm(r)]
num_iter = 1000
for i in range(num_iter):
    x = x - tau * r
    conv_x.append(x)
    r = B.dot(x) - f
    conv_r.append(np.linalg.norm(r))
Spectrum of the matrix = [0.5+1.32287566j 0.5-1.32287566j]
plt.semilogy(conv_r)
plt.xlabel("Number of iteration, $k$", fontsize=20)
plt.ylabel("Residual norm", fontsize=20)
Text(0, 0.5, 'Residual norm')

plt.scatter([x[0] for x in conv_x], [x[1] for x in conv_x])
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.title("$x_0 = (0, -1)$", fontsize=20)
Text(0.5, 1.0, '$x_0 = (0, -1)$')

Улучшенные итерационные методы

Но прежде чем перейти к предобуславливателям, мы можем использовать улучшенные итерационные методы.

Существует целый зоопарк итерационных методов, но нам нужно знать лишь некоторые из них.

Попытка 1: Метод наискорейшего спуска

  • Предположим, что мы меняем \tau на каждом шаге, т.е.

x_{k+1} = x_k - \tau_k (A x_k - f).

  • Возможный выбор \tau_k - такой, который минимизирует норму текущего остатка

\tau_k = \arg\min_{\tau} \|A(x_k - \tau_k (A x_k - f)) - f\|_2^2.

  • Эта задача может быть решена аналитически (выведите это решение!)

\tau_k = \frac{r_k^{\top}r_k}{r_k^{\top}Ar_k}, \quad r_k = Ax_k - f

  • Этот метод называется метод наискорейшего спуска.

  • Однако, он все еще сходится аналогично итерации Ричардсона.

Попытка 2: Итерация Чебышева

Другой способ найти \tau_k - это рассмотреть

e_{k+1} = (I - \tau_k A) e_k = (I - \tau_k A) (I - \tau_{k-1} A) e_{k-1} = \ldots = p(A) e_0,

где p(A) - это матричный многочлен (простейшая матричная функция)

p(A) = (I - \tau_k A) \ldots (I - \tau_0 A),

и p(0) = 1.

Оптимальный выбор временных шагов

Ошибка записывается как

e_{k+1} = p(A) e_0,

и, следовательно,

\|e_{k+1}\| \leq \|p(A)\| \|e_0\|,

где p(0) = 1 и p(A) - это матричный многочлен.

Для получения лучшего уменьшения ошибки, нам нужно минимизировать

\Vert p(A) \Vert

по всем возможным многочленам p(x) степени k+1 таким, что p(0)=1. Мы будем использовать \|\cdot\|_2.

Многочлены, наименее отклоняющиеся от нуля

Важный частный случай: A = A^* > 0.

Тогда A = U \Lambda U^*,

и

\Vert p(A) \Vert_2 = \Vert U p(\Lambda) U^* \Vert_2 = \Vert p(\Lambda) \Vert_2 = \max_i |p(\lambda_i)| \overset{!}{\leq} \max_{\lambda_\min \leq \lambda {\leq} \lambda_\max} |p(\lambda)|.

Последнее неравенство является единственным приближением. Здесь мы делаем ключевое предположение , что мы не хотим извлекать выгоду из распределения спектра между \lambda_\min и \lambda_\max.

Таким образом, нам нужно найти многочлен такой, что p(0) = 1, который имеет наименьшее возможное отклонение от 0 на [\lambda_\min, \lambda_\max].

Многочлены, наименее отклоняющиеся от нуля (2)

Мы можем выполнить аффинное преобразование интервала [\lambda_\min, \lambda_\max] в интервал [-1, 1]:

\xi = \frac{{\lambda_\max + \lambda_\min - (\lambda_\min-\lambda_\max)x}}{2}, \quad x\in [-1, 1].

Задача тогда сводится к проблеме нахождения многочлена, наименее отклоняющегося от нуля на интервале [-1, 1].

Точное решение: многочлены Чебышёва

Точное решение этой задачи даётся знаменитыми многочленами Чебышёва вида

T_n(x) = \cos (n \arccos x)

Что нужно знать о многочленах Чебышёва

  1. Это многочлен!

  2. Мы можем выразить T_n через T_{n-1} и T_{n-2}:

T_n(x) = 2x T_{n-1}(x) - T_{n-2}(x), \quad T_0(x)=1, \quad T_1(x)=x

  1. |T_n(x)| \leq 1 на x \in [-1, 1].

  2. Он имеет (n+1) точек альтернанса, где достигается максимальное абсолютное значение (это достаточное и необходимое условие для оптимальности) (теорема об альтернансе Чебышёва, доказательство здесь не приводится).

  3. Корни имеют вид

n \arccos x_k = \frac{\pi}{2} + \pi k, \quad \rightarrow\quad x_k = \cos \frac{\pi(2k + 1)}{2n}, \; k = 0, \ldots,n-1

Мы можем их построить…

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x1 = np.linspace(-1, 1, 128)
x2 = np.linspace(-1.1, 1.1, 128)
p = np.polynomial.Chebyshev((0, 0, 0, 0, 0, 0, 0, 0, 0, 1), (-1, 1)) #These are Chebyshev series, a proto of "chebfun system" in MATLAB
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(x1, p(x1))
ax1.set_title('Interval $x\in[-1, 1]$')
ax2.plot(x2, p(x2))
ax2.set_title('Interval $x\in[-1.1, 1.1]$')
Text(0.5, 1.0, 'Interval $x\\in[-1.1, 1.1]$')

Сходимость метода Ричардсона с ускорением Чебышёва

Заметим, что p(x) = (1-\tau_n x)\dots (1-\tau_0 x), следовательно, корни p(x) равны 1/\tau_i, и нам дополнительно нужно отобразить интервал [-1,1] обратно в [\lambda_\min, \lambda_\max]. Это приводит к формуле

\tau_i = \frac{2}{\lambda_\max + \lambda_\min - (\lambda_\max - \lambda_\min)x_i}, \quad x_i = \cos \frac{\pi(2i + 1)}{2n}\quad i=0,\dots,n-1

Сходимость (мы приводим только результат без доказательства) теперь задается выражением

e_{k+1} \leq C q^k e_0, \quad q = \frac{\sqrt{\mathrm{cond}(A)}-1}{\sqrt{\mathrm{cond}(A)}+1},

что лучше, чем в методе Ричардсона.

import numpy as np
import matplotlib.pyplot as plt

n = 64
ex = np.ones(n);

A = sp.sparse.spdiags(np.vstack((-ex,  2*ex, -ex)), [-1, 0, 1], n, n, 'csr'); 
rhs = np.ones(n)
ev1, vec = spla.eigsh(A, k=2, which='LA')
ev2, vec = spla.eigsh(A, k=2, which='SA')
lam_max = ev1[0]
lam_min = ev2[0]

niters = 64
roots = [np.cos((np.pi * (2 * i + 1)) / (2 * niters)) for i in range(niters)]
taus = [(lam_max + lam_min - (lam_min - lam_max) * r) / 2 for r in roots]
x = np.zeros(n)
r = A.dot(x) - rhs
res_cheb = [np.linalg.norm(r)]

print(1/np.array(taus))

for i in range(niters):
    x = x - 1.0 / taus[i] * r
    r = A.dot(x) - rhs
    res_cheb.append(np.linalg.norm(r))
    
plt.semilogy(res_richardson, label="Richardson")
plt.semilogy(res_cheb, label="Chebyshev")
plt.legend(fontsize=20)
plt.xlabel("Number of iterations, $k$", fontsize=20)
plt.ylabel("Residual norm, $\|Ax_k - b\|_2$", fontsize=20)
plt.xticks(fontsize=20)
_ = plt.yticks(fontsize=20)
[2.50622630e-01 2.50924658e-01 2.51530169e-01 2.52442095e-01
 2.53664865e-01 2.55204461e-01 2.57068471e-01 2.59266178e-01
 2.61808653e-01 2.64708878e-01 2.67981890e-01 2.71644951e-01
 2.75717745e-01 2.80222614e-01 2.85184826e-01 2.90632893e-01
 2.96598938e-01 3.03119118e-01 3.10234126e-01 3.17989766e-01
 3.26437629e-01 3.35635888e-01 3.45650225e-01 3.56554925e-01
 3.68434169e-01 3.81383570e-01 3.95511991e-01 4.10943730e-01
 4.27821135e-01 4.46307759e-01 4.66592186e-01 4.88892691e-01
 5.13462953e-01 5.40599097e-01 5.70648426e-01 6.04020340e-01
 6.41200067e-01 6.82766079e-01 7.29412363e-01 7.81977158e-01
 8.41480389e-01 9.09172947e-01 9.86602291e-01 1.07570085e+00
 1.17890673e+00 1.29933105e+00 1.44099335e+00 1.60915915e+00
 1.81083297e+00 2.05549470e+00 2.35622630e+00 2.73148428e+00
 3.20797673e+00 3.82550533e+00 4.64546321e+00 5.76650218e+00
 7.35516474e+00 9.71018736e+00 1.34098577e+01 1.96891165e+01
 3.15513749e+01 5.76947930e+01 1.29218671e+02 3.40581915e+02]

Что случилось с замечательными итерациями Чебышева?

niters = 64
roots = [np.cos((np.pi * (2 * i + 1)) / (2 * niters)) for i in range(niters)]
taus = [(lam_max + lam_min - (lam_min - lam_max) * r) / 2 for r in roots]
x = np.zeros(n)
r = A.dot(x) - rhs
res_cheb_even = [np.linalg.norm(r)]
#print(taus)

# Implementation may be non-optimal if number of iterations is not power of two
def leb_shuffle_2n(n):
    if n == 1:
        return np.array([0,], dtype=int)
    else:
        prev = leb_shuffle_2n(n // 2)
        ans = np.zeros(n, dtype=int)
        ans[::2] = prev
        ans[1::2] = n - 1 - prev
        return ans

good_perm_even = leb_shuffle_2n(niters)
print(good_perm_even, len(good_perm_even))
# good_perm_even = np.random.permutation([i for i in range(niters)])
ts = np.array(taus)[good_perm_even]
plt.figure()
plt.plot(1/ts)
for i in range(niters):
    x = x - 1.0/taus[good_perm_even[i]] * r
    r = A.dot(x) - rhs
    res_cheb_even.append(np.linalg.norm(r))
plt.figure()
    
plt.semilogy(res_richardson, label="Richardson")
plt.semilogy(res_cheb_even, label="Chebyshev")
plt.legend(fontsize=20)
plt.xlabel("Number of iterations, $k$", fontsize=20)
plt.ylabel("Residual norm, $\|Ax_k - b\|_2$", fontsize=20)
plt.xticks(fontsize=20)
_ = plt.yticks(fontsize=20)
[ 0 63 31 32 15 48 16 47  7 56 24 39  8 55 23 40  3 60 28 35 12 51 19 44
  4 59 27 36 11 52 20 43  1 62 30 33 14 49 17 46  6 57 25 38  9 54 22 41
  2 61 29 34 13 50 18 45  5 58 26 37 10 53 21 42] 64

  • Перестановка корней многочлена Чебышева имеет решающее влияние на сходимость
  • Об оптимальной перестановке вы можете прочитать в статье (В. Лебедев, С. Финогенов 1971) (ru, en)

Chebfun проект

  • Проект с открытым исходным кодом для численных вычислений (Python и Matlab интерфейсы)
  • Основан на численных алгоритмах, работающих с кусочно-полиномиальными интерполянтами и многочленами Чебышева
  • Этот проект был инициирован Ником Трефетеном и его студентом Захари Баттлсом в 2002 году, см. статью о проекте chebfun в SISC
  • Инструментарий Chebfun в основном фокусируется на следующих задачах
    • Аппроксимация
    • Квадратуры
    • Обыкновенные дифференциальные уравнения
    • Уравнения в частных производных
    • Поиск корней
    • Глобальная оптимизация в одномерном случае

А что есть кроме Чебышева?

  • Мы сделали важное предположение о спектре: он содержится в интервале на вещественной прямой (и нам нужно знать границы)

  • Если спектр содержится в двух интервалах, и мы знаем границы, мы также можем поставить задачу оптимизации для оптимального многочлена.

Спектр матрицы, содержащийся в нескольких сегментах

  • Для случая двух сегментов наилучший многочлен задается многочленами Золотарева (выраженными в терминах эллиптических функций). Оригинальная работа была опубликована в 1877 году, подробности здесь

  • Для случая более чем двух сегментов наилучший многочлен может быть выражен в терминах гиперэллиптических функций

Как мы можем улучшить это

  • Реализация ускорения Чебышева требует знания спектра.

  • Она хранит только предыдущий вектор x_k и вычисляет новый вектор коррекции

r_k = A x_k - f.

  • Это относится к классу двухчленных итерационных методов, т.е. он аппроксимирует x_{k+1}, используя 2 вектора: x_k и r_k.

  • Оказывается, что если мы храним больше векторов, то мы можем обойтись без оценки спектра (и получить лучшую сходимость на практике)!

Ключевой момент: подпространство Крылова

Метод Чебышева производит приближение вида

x_{k+1} = x_0 + p(A) r_0,

т.е. оно лежит в подпространстве Крылова матрицы, которое определяется как

\mathcal{K}_k(A, r_0) = \mathrm{Span}(r_0, Ar_0, A^2 r_0, \ldots, A^{k-1}r_0 )

Наиболее естественный подход - найти вектор в этом линейном подпространстве, который минимизирует определенную норму ошибки

Идея методов Крылова

Идея заключается в минимизации заданного функционала: - Энергетической нормы ошибки для систем с эрмитовыми положительно-определенными матрицами (метод CG). - Нормы невязки для систем с произвольными матрицами (методы MINRES и GMRES). - Отношения Рэлея для задач на собственные значения (метод Ланцоша).

Чтобы сделать методы практичными, необходимо: 1. Ортогонализировать векторы A^i r_0 подпространства Крылова для обеспечения устойчивости (процесс Ланцоша). 2. Вывести рекуррентные формулы для уменьшения вычислительной сложности.

Мы рассмотрим эти методы подробно на следующей лекции.

Основные выводы

  • Основная идея итерационных методов
  • Итерация Ричардсона: эрмитов и неэрмитов случай
  • Ускорение Чебышева
  • Определение подпространства Крылова

Вопросы?

from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()