A method of extended normal equations for tikhonov’s regulatization problems with differentiation operator

Full Text

1. Постановка задачи. Рассмотрим задачу решения системы линейных алгебраических уравнений (СЛАУ) Au = f, A ∈ Rm×n , f ∈ Rm . (1) Если система вида (1) является результатом дискретизации интегрального уравнения Фредгольма первого рода [1-4], то в этом случае известно, что она является плохо обусловленной. Для решения таких задач используется метод регуляризации Тихонова [5-10]. Регуляризованное решение СЛАУ (1) определяется выражением u∗ = arg min n u∈R Au - f 2 2 + α Lu 2 2 , (2) где α > 0 - параметр регуляризации, а L ∈ R(n-p)×n , p = 1, 2, . . . - дискретный аналог оператора дифференцирования для случая p = 1 (оператор дифференцирования первого порядка)   1 -1 0 . . . 0 0 0 1 -1 . . . 0 0  L= ∈ R(n-1)×n . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 . . . 1 -1 Если p = 2 (оператор дифференцирования второго порядка), то   1 -2 1 0 . . . 0 0 0 0 1 -2 1 . . . 0 0 0 ∈ R(n-2)×n . L= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 . . . 1 -2 1 В случае L = E, где E - единичная матрица порядка n, мы имеем стандартную задачу регуляризации Тихонова. Если L ∈ R(n-p)×n , то матрица L не имеет обратной. Известные приемы основаны на преобразовании (2) к стандартной задаче регуляризации Тихонова [3, 11]. Однако в этом случае требуется вычислять псевдообратную матрицу L+ . Возможно вычислить решение в (2), как и для любой задачи наименьших квадратов, с помощью следующей регуляризованной системы нормальных уравнений: A A + αL L u = A f. (3) Однако такой подход невыгоден, как будет показано в данной работе, в силу большого числа обусловленности матрицы данной системы. 133 Ж д а н о в А. И., М и х а й л о в И. А. В данной работе предлагается подход к решению исходной задачи (2) на основе её приведения к задаче решения эквивалентной нормальной расширенной регуляризованной системы. Этот подход не требует дополнительных преобразований исходных данных [12]. 2. Регуляризация на основе расширенных систем. Регуляризованную нормальную систему уравнений (3) можно записать в виде A r - αL Lu = 0, где r = f - Au, или α-1/2 A r - α1/2 L Lu = 0. (4) Объединяя r + Au = f и (4), получаем систему r + Au = f, α или -1/2 A r-α 1/2 L Lu = 0 α1/2 Em A 1/2 L L A -α α1/2 r u = f . 0 (5) Если обозначить y = α-1/2 r и ω = α1/2 , то (5) можно записать в виде ωEm A A -ωL L y u = f 0 ˜ ˜ ⇔ Aω x = f , (6) где ˜ Aω = ωEm A , A -ωL L x= y , u ˜ f= f . 0 Систему расширенных нормальных уравнений (6) удобно использовать на практике, так как матрица L L довольно просто вычисляется аналитически. Для оператора p-того порядка она имеет 2p + 1 диагональную структуру, так для оператора первого порядка матрица L L имеет следующий вид:   1 -1 0 0 ... 0 0 0 -1 2 -1 0 . . . 0 0 0   0 0  0 -1 2 -1 . . . 0 L L=  ∈ Rn×n , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 . . . -1 2 -1 0 0 0 0 . . . 0 -1 1 а для оператора второго порядка - следующий:   1 -2 1 0 0 ... 0 0 0 0 -2 5 -4 1 0 ... 0 0 0 0    1 -4 6 -4 1 . . . 0 0 0 0   1 -4 6 -4 . . . 0 0 0 0 0 L L=  ∈ Rn×n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 . . . -4 6 -4 1    0 0 0 0 0 . . . 1 -4 5 -2 0 0 0 0 0 ... 0 1 -2 1 134 Метод расширенных нормальных уравнений для задач регуляризации . . . Обозначим L L(p = 1) = L1 , L L(p = 2) = L2 , тогда в аналитическом виде данные матрицы будут иметь вид  L1 (i, i) = 1, L1 (i, i + 1) = -1, если i = 1;  L (i, i - 1) = -1, L1 (i, i) = 2, L1 (i, i + 1) = -1, если i = 2 . . . n - 1;  1 L (i, i - 1) = -1, L (i, i) = 1, если i = n, 1 1  L1 (i, i) = 1, L1 (i, i + 1) = -2, L1 (i, i + 2) = 1,    L1 (i, i - 1) = -2, L1 (i, i) = 5,     L (i, i + 1) = -4, L (i, i + 2) = 1,  1  1  L (i, i - 2) = 1, L (i, i - 1) = -4, L (i, i) = 6, 1 1 1  L1 (i, i + 1) = -4, L1 (i, i + 2) = 1,   L (i, i - 2) = 1, L (i, i - 1) = -4,  1  1    L (i, i) = 5, L (i, i + 1) = -2,  1 1   L (i, i - 2) = 1, L (i, i - 1) = -2, L (i, i) = 1, 1 1 1 если i = 1; если i = 2; если i = 3 . . . n - 2; если i = n - 1; если i = n. Таким образом, разреженная структура матрицы L L позволяет достаточно просто решать систему (6), например, с помощью метода исключения Гаусса. ˜ Найдём верхнюю границу для числа обусловленности матрицы Aω . Соб˜ω определяются из уравнения ственные значения матрицы A ωEm A A -ωL L z v =λ z , v где z ∈ Rm , v ∈ Rn . Тогда ωz + Av = λz, A z - ωL Lv = λv. Выражая z = (λ - ω)-1 Av, получаем (λ-ω)-1 A Av-ωL Lv = λv⇐⇒ A A + (ω 2 - λω)(L L - E) v = (λ2 -ω 2 )v. ˜ Точно вычислить спектр матрицы Aω невозможно, поэтому получим оценки сверху и снизу значения |λ|. Для некоторых встречающихся на практике ˜ частных случаев матрицы Aω возможно и точное нахождение спектра матрицы, а также её собственных векторов [13]. Обозначим Bω = A A + (ω 2 - λω)(L L - E) ˜ и найдём верхнюю границу |λ| для матрицы Aω . Для этого определим верхнюю границу для λmax матрицы Bω , так как этому значению будет соответствовать максимум выражения λ2 - ω 2 и, соответственно, максимум для |λ|. 135 Ж д а н о в А. И., М и х а й л о в И. А. По минимаксной теореме Куранта-Фишера [14] rT A A + (ω 2 - λω)(L L - E) r = 0=r∈R r r r A Ar rT (ω 2 - λω)(L L - E)r = max n + 0=r∈R r r r r T (ω 2 - λω)(L L - E)r r A Ar r max + max n = 0=r∈Rn 0=r∈R r r r r λmax (Bω ) = max n = λmax (A A) + λmax (ω 2 - λω)(L L - E) . Рассмотрим два случая. 1. Если λ > 0, тогда (ω 2 - λω) < 0 и λmax (ω 2 - λω)(L L - E) = (ω 2 - λω) σn (L L) - 1 , где σn обозначает n-ное сингулярное число матрицы. Так как L L - вырожденная матрица, σn (L L) = 0 и λmax (ω 2 - λω)(L L - E) = λω - ω 2 . Следовательно, λmax (Bω ) = λ2 - ω 2 2 λω - ω 2 + σ1 (A), где σ1 - первое (максимальное) сингулярное число матрицы; σ1 σ2 ... 2 Решая неравенство λ2 - λω - σ1 (A) |λ| ω+ σn . 0, окончательно получаем 2 ω 2 + 4σ1 (A) . 2 (7) 2. Если λ < 0, то (ω 2 - λω) > 0 и 2 λmax (ω 2 - λω)(L L - E) = (ω 2 - λω) σ1 (L) - 1 . Получаем λmax (Bω ) = λ2 - ω 2 2 2 (ω 2 - λω) σ1 (L) - 1 + σ1 (A). Решая это неравенство, окончательно находим |λ| 136 2 ω σ1 (L) - 1 + 2 ω 2 σ1 (L) + 1 2 2 2 + 4σ1 (A) . (8) Метод расширенных нормальных уравнений для задач регуляризации . . . Видно, что в (8) граница для |λ| больше, чем в (7), поэтому окончательно |λ| 2 ω σ1 (L) - 1 + 2 ω 2 σ1 (L) + 1 2 2 + 4σ1 (A) 2 . ˜ Найдём нижнюю границу |λ| для матрицы Aω . Для этого определим нижнюю границу для λmin матрицы (Bω ), так как этому значению будет соответствовать минимум выражения λ2 - ω 2 и, соответственно, минимум для |λ|. По минимаксной теореме Куранта-Фишера [14] rT A A + (ω 2 - λω)(L L - E) r = 0=r∈R r r r A Ar rT (ω 2 - λω)(L L - E)r = min n + 0=r∈R r r r r T (ω 2 - λω)(L L - E)r r A Ar r min n + min n = 0=r∈R 0=r∈R r r r r λmin (Bω ) = min n = λmin (A A) + λmin (ω 2 - λω)(L L - E) . При |λ| < ω ω 2 - λω > 0, поэтому λmin (ω 2 - λω)(L L - E) = (ω 2 - λω) σmin (L L) - 1 = λω - ω 2 , и λmin (Bω ) = λ2 - ω 2 2 λω - ω 2 + σn (A). Решая это неравенство, получаем 2 ω 2 + 4σn (A) - ω . 2 ˜ ˜ ˜ Так как Aω - симметричная матрица, σi (Aω ) = |λi (Aω )| ∀i = 1,2, . . . , n. Таким образом, получаем верхнюю границу для спектрального числа обу˜ словленности матрицы Aω : |λ| 2 ω σ1 (L) - 1 + ˜ σ1 (Aω ) ˜ µ(Aω ) = ˜ σn (Aω ) 2 ω 2 σ1 (L) + 1 2 2 + 4σ1 (A) 2 ω 2 + 4σn (A) - ω . (9) Таким же образом можно получить верхнюю границу для спектрального числа обусловленности матрицы в (3). Обозначим A1 = A A + ω 2 L L. Тогда верхняя граница для спектрального числа обусловленности матрицы A1 µ(A1 ) 2 2 σ1 (A) + ω 2 σ1 (L) . 2 σn (A) (10) Вычислим отношение границ для спектральных чисел обусловленности ˜ матриц Aω и A1 . С учётом (9) и (10) получаем µ(A1 ) ≈ ˜ µ(Aω ) 2 2 σ1 (A) + ω 2 σ1 (L) 2 σn (A) ω 2 σ1 (L) -1 + 2 ω 2 + 4σn (A) - ω ω2 2 σ1 (L) +1 2 + . 2 4σ1 (A) 137 Ж д а н о в А. И., М и х а й л о в И. А. Рассмотрим случай, когда выполняются условия σn (A) ω, ω σ1 (L) σ1 (A), тогда предыдущее выражение принимает вид σ1 (A) и µ(A1 ) σ1 (A) ≈ . ˜ ω µ(Aω ) (11) Данная оценка хотя и получена при определенных условиях, к тому же не для отношения самих чисел обусловленности, а только их верхних границ, тем не менее, как будет видно из тестовых исследований, хорошо аппрокси˜ мирует отношение спектральных чисел обусловленности матриц A1 и Aω . Чтобы показать, насколько велика оценка (11), выберем параметр регуляризации следующим образом [5]: 2 δ · σ1 (A) , f 2+δ √ где δ - погрешность задачи. Так как ω = α, имеем α= µ(A1 ) σ1 (A) ≈ = ˜ ω µ(Aω ) σ1 (A) 2 δ·σ1 (A) f 2 +δ = f 2 δ +δ ≈ f 2 . δ Полученное отношение в реальных задачах велико, так как норма вектора правой части матричного уравнения намного больше погрешности задачи; если бы эти величины были сопоставимы, то задачу не имело бы смысла решать из-за слишком большого шума в исходных данных. Из полученных результатов можно определить, при каком условии верх˜ ˜ ний предел для σ1 (Aω ), а также нижний предел для σn (Aω ) достигаются. Так как в выражении rT A A + (ω 2 - λω)(L L - E) r = 0=r∈R r r r A Ar rT (ω 2 - λω)(L L - E)r = max n + 0=r∈R r r r r T (ω 2 - λω)(L L - E)r r A Ar r max + max n = 0=r∈Rn 0=r∈R r r r r λmax (Bω ) = max n = λmax (A A) + λmax (ω 2 - λω)(L L - E) максимум достигается, когда r является собственным вектором, соответствующим наибольшему собственному значению, неравенство превращается в равенство при совпадении собственных векторов, соответствующих наибольшим собственным значениям матриц A A и (ω 2 - λω)(L L - E) , а так как собственные векторы этих матриц равны правым сингулярным векторам мат˜ риц A и L, верхний предел для σ1 (Aω ) достигается, когда совпадают первые правые сингулярные векторы матриц A и L. Аналогично, нижний предел ˜ для σn (Aω ) достигается, когда совпадают n-ные правые сингулярные векторы матриц A и L. 138 Метод расширенных нормальных уравнений для задач регуляризации . . . Условия достижения верхнего предела для σ1 (A1 ), а также нижнего пре˜ дела для σn (A1 ) аналогичны условиям для матрицы Aω , что подтверждает правильность оценки (11). 3. Тестовые исследования и выводы. В ходе тестовых исследований генерировались плохо обусловленные матрицы из R10×6 . Первые пять столбцов матрицы задавались случайно, а шестой являлся линейной комбинацией первых пяти столбцов и дополнительного столбца погрешностей, норма которого была в 100 раз меньше, чем у первых пяти столбцов. Из миллиона сгенерированных таким образом тестовых матриц максимальное число обусловлен˜ ности матрицы Aω составило 6 % от полученной в данной работе верхней границы (9). Это говорит о том, что хотя эта граница теоретически достижима, на практике число обусловленности реальных задач много меньше данной ˜ границы. Отношение спектральных чисел обусловленности матриц A1 и Aω в тестовых исследованиях в среднем составило 95, при этом оценка (11) была в среднем больше истинного значения в 1.23 раза, что говорит о том, что эта оценка хорошо аппроксимирует отношение спектральных чисел обусловлен˜ ности матриц A1 и Aω . Таким образом можно сделать вывод, что метод расширенных нормальных уравнений можно использовать для решения задачи (2), а метод (3) невыгодно использовать для вычисления регуляризованного решения из-за большого значения числа обусловленности данной системы.

About the authors

Alexandr I Zhdanov

Samara State Technical University

Email: zhdanovaleksan@yandex.ru
(Dr. Phys. & Math. Sci.; zhdanovaleksan@yandex.ru), Dean, Faculty of the Distance and Additional Education 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ivan A Mikhaylov

Samara State Technical University

Email: mikhaylovivan90@mail.ru
(mikhaylovivan90@mail.ru; Corresponding Author), Postgraduate Student, Dept. of Higher Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files