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

Abstract


This article is devoted to a new method of ill-conditioned linear algebraic systems solving with the help of differentiation operator. These problems appear while solving the first kind integral Fredholm equations. The most difficult thing about this method is that differential operator discrete analogue matrix is rank deficiency matrix. The generalized singular value decomposition methods are used to solve those problems. The approach has high computational complexity. This also leads to additional computational error. Our method is based on the original regularized problem transformation into equivalent augmented regularized normal equation system using differential operator discrete analogue. The problem of spectrum matrix investigation of augmented regularized normal equation system with rank deficiency differential operator discrete analogue matrix is very relevant nowadays. Accurate eigenvalue spectrum research for this problem is impossible. That is why we estimated spectrum matrix bounds. Our estimation is based on a wellknown Courant-Fisher theorem. It is shown that estimated spectrum matrix bounds are rather accurate. The comparison between the proposed method and standard method based on the solving of normal system of equations is done. As shown in the paper, the condition number of normal method matrix is bigger than the condition number of augmented normal equations method matrix. In conclusion test problems description is given which proves our theoretical background.

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
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation
(Dr. Phys. & Math. Sci.; zhdanovaleksan@yandex.ru), Dean, Faculty of the Distance and Additional Education

Ivan A Mikhaylov

Samara State Technical University

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

References

  1. Abdelmalek N. N. A program for the solution of ill-posed linear systems arising from the discretization of the Fredholm integral equation of the first kind // Computer Physics Communications, 1990. vol. 58, no. 3. pp. 285-292. doi: 10.1016/0010-4655(90)90064-8.
  2. Delves L. M., Mohamed J. L. Computational Methods for Integral Equations. Cambridge: Cambridge University Press, 1985. 376+xii pp. doi: 10.1017/CBO9780511569609.
  3. Hansen P. C. REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems // Numerical Algorithms, 1994. vol. 6, no. 1. pp. 1-35. doi: 10. 1007/BF02149761.
  4. Bouhamidi A., Jbilou K., Reichel L., Sadok H. An extrapolated TSVD method for linear discrete ill-posed problems with Kronecker structure // Linear Algebra and Its Applications, 2011. vol. 434, no. 7. pp. 1677-1688. doi: 10.1016/j.laa.2010.06.001.
  5. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979. 286 с.
  6. Phillips D. L. A technique for the numerical solution of certain integral equations of the first kind // JACM, 1962. vol. 9, no. 1. pp. 84-97. doi: 10.1145/321105.321114.
  7. Björck Å., Eldén L. Methods in numerical algebra for ill-posed problems: Technical Report LiTH-MAT-R33-1979. Linköping, Sweden, 1979. 267 pp.
  8. Wing G. M. A Primer on Integral Equations of the First Kind / Other Titles in Applied Mathematics. Los Alamos, New Mexico: Los Alamos National Laboratory, 1991. 141+xiv pp. doi: 10.1137/1.9781611971675.
  9. Bauer F., Lukas M. A. Comparingparameter choice methods for regularization of ill-posed problems // Mathematics and Computers in Simulation, 2011. vol. 81, no. 9. pp. 1795-1841. doi: 10.1016/j.matcom.2011.01.016.
  10. Liu C.-S. A dynamical Tikhonov regularization for solving ill-posed linear algebraic systems // Acta Applicandae Mathematicae, 2013. vol. 123, no. 1. pp. 285-307. doi: 10.1007/s10440-012-9766-3.
  11. Hansen P. C. Regularization Tools version 4.0 for Matlab 7.3 // Numer. Algor., 2007. vol. 46, no. 2. pp. 189-194. doi: 10.1007/s11075-007-9136-9.
  12. Жданов А. И. Метод расширенных регуляризованных нормальных уравнений // Ж. вычисл. матем. и матем. физ., 2012. Т. 52, № 2. С. 205-208.
  13. Stor N. J., Slapničar I., Barlow J. L. Accurate eigenvalue decomposition of real symmetric arrowhead matrices and applications // Linear Algebra and its Application, 2015. vol. 464, no. 1. pp. 62-89, arXiv: 1302.7203 [math.NA]. doi: 10.1016/j.laa.2013.10.007.
  14. Demmel J. W. Applied Numerical Linear Algebra / Other Titles in Applied Mathematics. Berkeley: University of California, 1997. 416+xi pp. doi: 10.1137/1.9781611971446.

Statistics

Views

Abstract - 19

PDF (Russian) - 4

Cited-By


PlumX

Dimensions

Refbacks

  • There are currently no refbacks.

Copyright (c) 2014 Samara State Technical University

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies