Метод спряженого градієнта

У математиці метод спря́женого градієнта є алгоритмом чисельного рішення окремих систем лінійних рівнянь, а саме тих, чия матриця симетрична і позитивно-визначена. Метод спряженого градієнта часто реалізовується як ітераційний алгоритм, застосовний до розріджених систем, які занадто великі, щоб обробляти їх шляхом прямої реалізації або інших прямих методів, таких як декомпозиція Холеського. Великі розріджені системи часто виникають при чисельному вирішенні часткових диференціальних рівнянь або задачах оптимізації.

Порівняння збіжності градієнтного спуску з оптимальним розміром кроку (зеленим) та кон'югованим вектором (червоним кольором) для мінімізації квадратичної функції, пов'язаної із заданою лінійною системою. Спряжений градієнт, припускаючи точну арифметику, сходиться не більше n кроків, де n - розмір матриці системи (тут n = 2).

Метод спряженого градієнта також може бути використаний для вирішення необмежених задач оптимізації, таких як мінімізація енергії . Його в основному розробили Магнус Гестенес та Едуард Стіфель [1] які запрограмували його на Z4 .

Метод двобічного спряженого градієнта забезпечує узагальнення до несиметричних матриць. Різні методи нелінійного спряженого градієнта шукають мінімуми нелінійних рівнянь.

Опис задачі, котру вирішують сполучені градієнти

Припустимо, ми хочемо розв’язати систему лінійних рівнянь

\mathbf {A} \mathbf {x} =\mathbf {b}

для вектора x, де відома n × n матриця A симетрична (тобто A ^T = A ), позитивно-визначена (тобто x ^T Ax > 0 для всіх ненульових векторів x в R ⁿ ), і реальна, і b також відомо. Позначимо унікальний розв'язок цієї системи через $\mathbf {x} _{*}$ .

Прямий метод

Ми припускаємо, що два ненульові вектори u і v є сполученими (щодо А ), якщо

\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.

Оскільки A симетрична і позитивно-визначена, ліва частина визначає внутрішній добуток

\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .

Два вектори є сполученими тоді і лише тоді, коли вони ортогональні щодо цього внутрішнього добутку. Будучи сполученим - це симетричне відношення: якщо u є спряженим на v, то v є спряженим на u . Припустимо, що

P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}

являє собою сукупність n взаємно сполучених векторів (щодо А ). Тоді $P$ становить основу для $\mathbb {R} ^{n}$ , і ми можемо висловити рішення $x *$ of $\mathbf {Ax} =\mathbf {b}$ виходячи з цього:

\mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}.

На основі цього розширення ми обчислюємо:

\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.

Ліву частину множимо на $\mathbf {p} _{k}^{\mathsf {T}}$ :

\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i},

підставляючи $\mathbf {Ax_{*}} =\mathbf {b}$ і $\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }$ :

\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} },

потім $\mathbf {u} ^{\mathsf {T}}\mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle$ і використання $\forall i\neq k:\langle \mathbf {p} _{k},\mathbf {p} _{i}\rangle _{\mathbf {A} }=0$ врожайність

\langle \mathbf {p} _{k},\mathbf {b} \rangle =\alpha _{k}\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} },

що означає

\alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.

Це дає наступний метод розв’язання рівняння $Ax = b$ : знайти послідовність n спрямованих напрямків, а потім обчислити коефіцієнти $α k$ .

Як ітеративний метод

Якщо ми обережно оберемо сполучені вектори p _k, то, можливо, нам не знадобляться всі, щоб отримати гарне наближення до рішення $x *$ . Отже, ми хочемо розглянути метод спряженого градієнта як ітераційний метод. Це також дозволяє приблизно вирішити системи, де n настільки велике, що прямий метод зайняв би занадто багато часу.

Позначимо початкове припущення для $x *$ через $x 0$ (можна без втрати загальності вважати, що $x 0 = 0$ , інакше розглянемо систему Az = b - Ax ₀ ). Починаючи з x₀ ми шукаємо вирішення і в кожній ітерації ми повинні мати метрику, котра зможєе сказати нам чи ми ближче до вирішення $x *$ , нам це невідомо). Ця метрика випливає з того, що рішення $x *$ також є унікальним мінімізатором наступної квадратичної функції

f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.

Існування унікального мінімізатора очевидно, оскільки його друга похідна задана симетричною позитивно-визначеною матрицею

\nabla ^{2}f(\mathbf {x} )=\mathbf {A} \,,

і що мінімалізатор (виокристовує Df(x) = 0) вирішує початкову задачу очевидно з її першої похідної

\nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.

Це говорить про те, щоб перший базовий вектор p ₀ був від'ємним градієнтом f при x = x ₀ . Градієнт f дорівнює $Ax - b$ . Починаючи з початкової здогадки x ₀, це означає, що беремо p ₀ = b - Ax ₀ . Інші вектори в основі будуть спряжені з градієнтом, звідси і назва метод спряженого градієнта . Зауважимо, що p ₀ також є залишковим, передбаченим цим початковим кроком алгоритму.

Нехай r _k - залишок на k- му кроці:

\mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.

Як було зазначено вище, r _k - від'ємний градієнт f при x = x _k, тому метод спуску градієнтом потребує руху в напрямку r _k . Тут, однак, ми наполягаємо, щоб напрямки p _k були сполучені один з одним. Практичний спосіб забезпечити це - вимагаючи, щоб наступний напрямок пошуку був побудований з поточного залишкового та всіх попередніх напрямків пошуку. [2] Це дає такий вираз:

\mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}

(див. малюнок у верхній частині статті про вплив обмеження спряженості на збіжність). Слідуючи цьому напрямку, наступне оптимальне місце задається

\mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}

з

\alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},

де остання рівність випливає з визначення r _k . Вираз для $\alpha _{k}$ може бути отримано, якщо підміняти вираз x _{k +1} на f і мінімізувати його wrt $\alpha _{k}$

{\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}

Отриманий алгоритм

Наведений вище алгоритм дає найбільш просте пояснення методу спряженого градієнта. Здається, алгоритм, як заявлено, вимагає зберігання всіх попередніх напрямків пошуку та векторів залишків, а також багатьох матричних векторних множень, і, таким чином, може бути обчислювально дорогим. Однак більш детальний аналіз алгоритму показує, що r _i є ортогональним до r _j, тобто $\mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0$ , для i ≠ j. І p _i - A-ортогональна до p _j, тобто $\mathbf {p} _{i}^{\mathsf {T}}A\mathbf {p} _{j}=0$ , для i ≠ j. Це можна вважати, що в міру просування алгоритму p _i і r _i охоплюють той самий підпростір Крилова. Якщо r _я утворює ортогональну основу відносно стандартного внутрішнього добутку, а p _i утворює ортогональну основу відносно внутрішнього добутку, індукованого А. Тому x _k можна розглядати як проекцію x на підпростір Крилова.

Алгоритм детально описаний нижче для розв’язання Ax = b, де A - реальна, симетрична, позитивно-визначена матриця. Вхідний вектор x ₀ може бути приблизним початковим рішенням або 0 . Це інша рецептура точної процедури, описаної вище.

{\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&{\hbox{if }}\mathbf {r} _{0}{\text{ is sufficiently small, then return }}\mathbf {x} _{0}{\text{ as the result}}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if }}\mathbf {r} _{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}

Це найбільш часто використовуваний алгоритм. Така ж формула для $β k$ також використовується в нелінійному методі градієнта Флетчера-Рівза.

Розрахунок альфа та бета-версії

В алгоритмі $α k$ вибирається таким, що $\mathbf {r} _{k+1}$ є ортогональним до r _k . Знаменник спрощено від

\alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}

з тих пір $\mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}$ . $β k$ вибирається таким, що $\mathbf {p} _{k+1}$ сполучається з p _k . Спочатку $β k$ є

\beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}

використовуючи

\mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}

і рівнозначно

$\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),$

чисельник $β k$ переписується як

\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}

оскільки $\mathbf {r} _{k+1}$ і r _k є ортогональними за конструкцією. Знаменник переписується як

\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}

використовуючи, що напрямки пошуку p _k кон'югуються і знову, що залишки є ортогональними. Це дає $β$ в алгоритмі після скасування $α k$ .

Приклад коду в MATLAB / GNU Octave

function x = conjgrad(A, b, x)
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-10
              break;
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end

Числовий приклад

Розглянемо лінійну систему Ax = b, задану через

\mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},

ми виконаємо два етапи методу спряженого градієнта, починаючи з початкової здогадки

\mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}

щоб знайти приблизне рішення для системи.

Рішення

Для довідки правильне рішення

\mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}

Наш перший крок - обчислити залишковий вектор r _0, пов'язаний з x ₀ . Цей залишок обчислюється за формулою r ₀ = b - Ax ₀, а в нашому випадку дорівнює

\mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}.

Оскільки це перша ітерація, ми будемо використовувати залишковий вектор r ₀ як наш початковий напрямок пошуку p ₀ ; метод вибору p _k зміниться в подальших ітераціях.

Тепер обчислимо скалярний $α 0$ використовуючи відношення

\alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}.

Тепер ми можемо обчислити х _1, використовуючи формулу

\mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.

Цей результат завершує першу ітерацію, результатом якої є "покращене" приблизне рішення для системи, x ₁ . Тепер ми можемо перейти і обчислити наступний залишковий вектор r ₁ за формулою

\mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.

Наступним нашим кроком у процесі є обчислення скалярного $β 0$ яке згодом буде використано для визначення наступного напрямку пошуку p ₁ .

\beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.

Тепер, використовуючи цей скаляр $β 0$ , ми можемо обчислити наступний напрямок пошуку p _1, використовуючи відношення

\mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.

Тепер ми обчислюємо скалярний $α 1$ використовуючи нещодавно придбаний p _1, використовуючи той самий метод, що і для $α 0$ .

\alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.

Нарешті, ми знаходимо х _2, використовуючи той самий метод, що і для знаходження х ₁ .

\mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.

Результат, x ₂, є "кращим" наближенням до рішення системи, ніж x ₁ і x ₀ . Якби точна арифметика повинна використовуватися в цьому прикладі замість обмеженої точності, то точне рішення теоретично було б досягнуте після n = 2 ітерацій ( n - це порядок системи).

Властивості збіжності

Метод спряженого градієнта теоретично можна розглядати як прямий метод, оскільки він дає точне рішення після кінцевого числа ітерацій, що не перевищує розмір матриці, за відсутності помилки округлення . Однак метод градієнта спряжених нестабільний щодо навіть невеликих збурень, наприклад, більшість напрямків на практиці не є сполученими, і точного рішення так і не отримати. На щастя, метод спряженого градієнта може бути використаний як ітераційний метод, оскільки він забезпечує монотонно поліпшення наближень $\mathbf {x} _{k}$ до точного рішення, яке може досягти необхідного допуску після відносно невеликої (порівняно з розміром проблеми) кількості ітерацій. Поліпшення, як правило, лінійне і його швидкість визначається числом умови $\kappa (A)$ системної матриці $A$ : тим більше $\kappa (A)$ є, чим повільніше поліпшення. [3]

Якщо $\kappa (A)$ велика, попередня умова використовується для заміни вихідної системи $\mathbf {Ax} -\mathbf {b} =0$ з $\mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0$ такий як $\kappa (\mathbf {M} ^{-1}\mathbf {A} )$ менше, ніж $\kappa (\mathbf {A} )$ , Дивіться нижче.

Теорема конвергенції

Визначте підмножину многочленів як

{\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\

де $\Pi _{k}$ - це множина многочленів максимального ступеня $k$ .

Дозволяти $\left(\mathbf {x} _{k}\right)_{k}$ бути ітераційним наближенням точного рішення $\mathbf {x} _{*}$ , і визначити помилки як $\mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}$ . Тепер швидкість конвергенції можна приблизно оцінити як [4]

{\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}

де $\sigma (\mathbf {A} )$ позначає спектр, і $\kappa (\mathbf {A} )$ позначає номер умови .

Зауважте, важлива межа, коли $\kappa (\mathbf {A} )$ схиляється до $\infty$

{\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.

Ця межа показує більш швидкий коефіцієнт конвергенції порівняно з ітераційними методами Якобі або Гаусса-Сейделя, які масштабуються як $\approx 1-{\frac {2}{\kappa (\mathbf {A} )}}$ .

Метод попередньо обумовленого градієнта

У більшості випадків попередня підготовка необхідна для забезпечення швидкої конвергенції методу градієнта спряжених. Метод попередньо обумовленого градієнта має такий вигляд:

\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}

\mathbf {z} _{0}:=\mathbf {M} ^{-1}\mathbf {r} _{0}

\mathbf {p} _{0}:=\mathbf {z} _{0}

k:=0\,

repeat

\alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}

\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}

\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}

if r_k+1 is sufficiently small then exit loop end if

\mathbf {z} _{k+1}:=\mathbf {M} ^{-1}\mathbf {r} _{k+1}

\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}

\mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}

k:=k+1\,

end repeat

The result is x_k+1

Вищевказаний склад еквівалентний застосуванню методу градієнта спряженого без попереднього обумовлення системи ^[1]

\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b}

де

\mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .

Матриця попереднього кондиціонера M повинна бути симетричною-позитивно визначеною і фіксованою, тобто не може змінюватися від ітерації до ітерації. Якщо будь-яке з цих припущень щодо попереднього кондиціонера порушено, поведінка методу попередньо обумовленого градієнта може стати непередбачуваним.

Прикладом часто використовуваного попереднього кондиціонера є неповна факторизація Холеського .

Метод гнучких попередньо обумовлених градієнтів

У важкозахисних програмах застосовуються складні попередні кондиціонери, що може призвести до змінної попередньої кондиціонування, що змінюється між ітераціями. Навіть якщо попередній кондиціонер є симетричним позитивно-визначеним на кожній ітерації, той факт, що він може змінитися, робить аргументи вище недійсними, а на практичних тестах призводить до значного уповільнення конвергенції алгоритму, представленого вище. Використовуючи формулу Поляка-Ріб'єра

\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\left(\mathbf {r} _{k+1}-\mathbf {r} _{k}\right)}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}

замість формули Флетчер-Рівз

\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}

може різко покращити конвергенцію в цьому випадку. [5] Цей варіант попередньо обумовленого методу градієнта кон'югату можна назвати [6] гнучким, оскільки він дозволяє змінювати попередню умову. Також показано, що гнучка версія [7] є надійною, навіть якщо попередній кондиціонер не є симетричним позитивним значенням (SPD).

Реалізація гнучкої версії вимагає зберігання додаткового вектора. Для фіксованого попереднього кондиціонера SPD, $\mathbf {z} _{k+1}^{\mathsf {T}}\mathbf {r} _{k}=0,$ тому обидві формули для $β k$ еквівалентні в точній арифметиці, тобто без похибки округлення .

Математичне пояснення кращої поведінки конвергенції методу за формулою Поляка-Ріб'єра полягає в тому, що метод в цьому випадку є локально оптимальним, зокрема, він не зближується повільніше, ніж локально оптимальний метод найбільш крутого спуску. [8]

Приклад коду в MATLAB / GNU Octave

function [x, k] = cgp(x0, A, C, b, mit, stol, bbA, bbC)
% Synopsis:
% x0: initial point
% A: Matrix A of the system Ax=b
% C: Preconditioning Matrix can be left or right
% mit: Maximum number of iterations
% stol: residue norm tolerance
% bbA: Black Box that computes the matrix-vector product for A * u
% bbC: Black Box that computes:
%      for left-side preconditioner : ha = C \ ra
%      for right-side preconditioner: ha = C * ra
% x: Estimated solution point
% k: Number of iterations done 
%
% Example:
% tic;[x, t] = cgp(x0, S, speye(1), b, 3000, 10^-8, @(Z, o) Z*o, @(Z, o) o);toc
% Elapsed time is 0.550190 seconds.
%
% Reference:
%  Métodos iterativos tipo Krylov para sistema lineales
%  B. Molina y M. Raydan - {{ISBN|908-261-078-X}}
        if nargin < 8, error('Not enough input arguments. Try help.'); end;
        if isempty(A), error('Input matrix A must not be empty.'); end;
        if isempty(C), error('Input preconditioner matrix C must not be empty.'); end;
        x = x0;
        ha = 0;
        hp = 0;
        hpp = 0;
        ra = 0;
        rp = 0;
        rpp = 0;
        u = 0;
        k = 0;

        ra = b - bbA(A, x0); % <--- ra = b - A * x0;
        while norm(ra, inf) > stol
                ha = bbC(C, ra); % <--- ha = C \ ra;
                k = k + 1;
                if (k == mit), warning('GCP:MAXIT', 'mit reached, no conversion.'); return; end;
                hpp = hp;
                rpp = rp;
                hp = ha;
                rp = ra;
                t = rp' * hp;
                if k == 1
                        u = hp;
                else
                        u = hp + (t / (rpp' * hpp)) * u;
                end;
                Au = bbA(A, u); % <--- Au = A * u;
                a = t / (u' * Au);
                x = x + a * u;
                ra = rp - a * Au;
        end;

Місцево оптимальний метод найбільш стрімкого спуску

І в оригінальному, і в попередньо обумовленому методах градієнта кон'югату потрібно лише встановити $\beta _{k}:=0$ щоб зробити їх локально оптимальними, використовуючи пошук лінії, найкрутіші методи спуску . При цій підстановці вектори $p$ завжди такі ж, як вектори $z$ , тому немає необхідності зберігати вектори $p$ . Таким чином, кожна ітерація цих найбільш стрімких методів спуску є дещо дешевшою порівняно з методом спряженого градієнта. Однак останні сходяться швидше, якщо не застосовується (високо) змінна та / або попередній кондиціонер, який не є SPD, див. Вище.

Виведення методу

Метод спряженого градієнта може бути отриманий з кількох різних точок зору, включаючи спеціалізацію методу спряженого спрямування для оптимізації та варіацію ітерації Арнольді / Ланцоса для проблем власного значення. Незважаючи на розбіжність у підходах, ці виводи поділяють загальну тему - доказуючи ортогональність залишків та сукупність напрямків пошуку. Ці дві властивості мають вирішальне значення для розробки добре відомого стислого способу.

Спряження градієнта на нормальних рівняннях

Кон'югат градиентного метод може бути застосований до довільного п матриця з розмірністю м матриці, застосовуючи його до нормальним рівнянням ^Т А і права частина вектора А ^Т Ь, так як ^Т А є симетричною позитивно-полуопределена матрицею для будь-якого А. Результат - це спряжений градієнт у звичайних рівняннях (CGNR).

A ^T Ax = A ^T b

Як ітераційний метод не потрібно явно формувати A ^T A в пам'яті, а лише виконувати матричний вектор і транспонувати множення матричного вектора. Отже, CGNR особливо корисний, коли A є розрідженою матрицею, оскільки ці операції зазвичай є надзвичайно ефективними. Однак недоліком формування нормальних рівнянь є те, що число умови κ ( A ^T A ) дорівнює κ ² ( A ), тому швидкість конвергенції CGNR може бути повільною і якість приблизного рішення може бути чутливою до округлення помилки. Пошук хорошого попереднього кондиціонера часто є важливою частиною використання методу CGNR.

Запропоновано кілька алгоритмів (наприклад, CGLS, LSQR). Нібито алгоритм LSQR має найкращу числову стійкість, коли A погано обумовлений, тобто A має велике число умов .

Дивись також

Метод двобічного градієнта (BiCG)
Спосіб кон'югації залишків
Пропаганда вірувань Гаусса
Ітеративний метод: Лінійні системи
Крилова підпростір
Метод нелінійного спряженого градієнта
Підготовка
Рідке множення матричного вектора

Примітки

Straeter, T. A. On the Extension of the Davidon-Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems. NASA.
The conjugation constraint is an orthonormal-type constraint and hence the algorithm bears resemblance to Gram-Schmidt orthonormalization.
Saad, Yousef (2003). Iterative methods for sparse linear systems (вид. 2nd). Philadelphia, Pa.: Society for Industrial and Applied Mathematics. с. 195. ISBN 978-0-89871-534-7.
Hackbusch, W. (21 червня 2016). Iterative solution of large sparse systems of equations (вид. 2nd). Switzerland: Springer. ISBN 9783319284835. OCLC 952572240.
Golub, Gene H.; Ye, Qiang (1999). Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration. SIAM Journal on Scientific Computing 21 (4): 1305. doi:10.1137/S1064827597323415. Проігноровано невідомий параметр |citeseerx= (довідка)
Notay, Yvan (2000). Flexible Conjugate Gradients. SIAM Journal on Scientific Computing 22 (4): 1444–1460. doi:10.1137/S1064827599362314. Проігноровано невідомий параметр |citeseerx= (довідка)
Henricus Bouwmeester, Andrew Dougherty, Andrew V Knyazev. Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods. Procedia Computer Science, Volume 51, Pages 276-285, Elsevier, 2015. https://doi.org/10.1016/j.procs.2015.05.241
Knyazev, Andrew V.; Lashuk, Ilya (2008). Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning. SIAM Journal on Matrix Analysis and Applications 29 (4): 1267. arXiv:math/0605767. doi:10.1137/060675290.

Список літератури

Спосіб спряженого градієнта спочатку був запропонований в

Hestenes, Magnus R.; Stiefel, Eduard (December 1952). Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards 49 (6): 409. doi:10.6028/jres.049.044.

Описи методу можна знайти в наступних підручниках:

Atkinson, Kendell A. (1988). Section 8.9. An introduction to numerical analysis (вид. 2nd). John Wiley and Sons. ISBN 978-0-471-50023-0.
Avriel, Mordecai (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. ISBN 978-0-486-43227-4.
Golub, Gene H.; Van Loan, Charles F. (15 жовтня 1996). Chapter 10. Matrix computations (вид. 3rd). Johns Hopkins University Press. ISBN 978-0-8018-5414-9.
Saad, Yousef (1 квітня 2003). Chapter 6. Iterative methods for sparse linear systems (вид. 2nd). SIAM. ISBN 978-0-89871-534-7.

зовнішні посилання

Hazewinkel, Michiel, ed. (2001) [1994], "Conjugate gradients, method of", Encyclopedia of Mathematics, Springer Science+Business Media B.V. / Kluwer Academic Publishers, ISBN 978-1-55608-010-4

This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[1] Straeter, T. A. On the Extension of the Davidon-Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems. NASA.

[2] The conjugation constraint is an orthonormal-type constraint and hence the algorithm bears resemblance to Gram-Schmidt orthonormalization.

[saad1996iterative-3] Saad, Yousef (2003). Iterative methods for sparse linear systems (вид. 2nd). Philadelphia, Pa.: Society for Industrial and Applied Mathematics. с. 195. ISBN 978-0-89871-534-7.

[4] Hackbusch, W. (21 червня 2016). Iterative solution of large sparse systems of equations (вид. 2nd). Switzerland: Springer. ISBN 9783319284835. OCLC 952572240.

[5] Golub, Gene H.; Ye, Qiang (1999). Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration. SIAM Journal on Scientific Computing 21 (4): 1305. doi:10.1137/S1064827597323415. Проігноровано невідомий параметр |citeseerx= (довідка)

[6] Notay, Yvan (2000). Flexible Conjugate Gradients. SIAM Journal on Scientific Computing 22 (4): 1444–1460. doi:10.1137/S1064827599362314. Проігноровано невідомий параметр |citeseerx= (довідка)

[7] Henricus Bouwmeester, Andrew Dougherty, Andrew V Knyazev. Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods. Procedia Computer Science, Volume 51, Pages 276-285, Elsevier, 2015. https://doi.org/10.1016/j.procs.2015.05.241

[8] Knyazev, Andrew V.; Lashuk, Ilya (2008). Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning. SIAM Journal on Matrix Analysis and Applications 29 (4): 1267. arXiv:math/0605767. doi:10.1137/060675290.