книги / Регулярные методы решения некорректно поставленных задач
..pdfПусть теперь N = PX\ P 2 - сингулярное разложение матрицыN, т.е. Pi и Р2 —ортогональные, Л - ’’диагональная” матрица. Используя тот факт, чтоха удовлетворяет уравнению
(NTN + otE)xa = N Tf,
и выполнив несложные преобразования в (5), получаем оконча тельные формулы
£ |
/ * \ 211/2 |
■ |
, ч г 4" / |
V I 1/2 |
■ |
*(«)= i.Urnv |
|
|
|||
|
i t |
|
|
|
(6) |
F(a)= 2 |
|
|
|
|
|
0i\f + 1 |
|
|
|
|
|
i = 1 |
|
|
|
|
где Л,- есть /-е сингулярное число матрицы N при /< и, X,
/ > есть /-я компонента вектора эТ/.
п, g( Р(
Представление (6) делает очевидным утверждение леммы о монотонности функций R s, ys и F s. Сформулированные в лемме свойства выпуклости следуют из того, что знак второй производной функции
Г |
2 |
а, |
J;(°0= |
1(Ъ(а. + а )к _ |
|
|
/ = |
( г д е bh Cj неотрицательны,к > 0 при а > 0, l> - \jk) совпадает со знаком / . Лемма доказана.
В дальнейшем значения s > —1, s Ф 0 удобно называть допусти мыми. Заметим еще, что функции R s, ys и F s легко доопределить
по непрерывности справа в точке а = 0 |
с сохранением гладкости. |
Рассмотрим теперь уравнения |
|
R s{a) = el, ys(a) = esy, F \a)= е*. |
(7) |
Непосредственным следствием леммы и условия разрешимости уравнений (2 )-(4 ) является
Т е о р е м а 99. Метод касательных Ньютона для определения корней уравнений (7) сходится для допустимых значений s при любом начальном приближении а > 0.
3. Выберем теперь значения степеней s, дающие наибольшее приращение параметра а на каждом шаге (начиная со второго) при
решении (7) |
методом касательных Ньютона. Для приращения Аа |
||
получаем |
|
|
|
Да = |
е* - |
Vs |
1 - i e J y Y |
— - — = |
(8) |
||
|
S)’s~ У |
У |
|
где у - |
одна из функций R, у или F, еу обозначает ер , еу или |
||
соответственно. |
|
14* |
211 |
Исследуем Act как функцию s. Используя лемму 53 при s = 1, получаем, что множитель ( -у / у' ) в (8) положителен, а положитель ная величина q = еу/у меньше единицы, по крайней мере, начиная
со второго шага. Отсюда следует, что при любом s Ф 0 |
выражение |
|
в квадратных скобках положительно. Так как функция |
(1 - |
qs)/s |
является убывающей функцией s, то ее наибольшее значение |
для |
s > —1 достигается при s = —1. Заметим, что при s < — 1 безуслов ная сходимость метода Ньютона для уравнений (7) не гарантирует ся, так что значение — 1 обеспечивает наибольшее приращение Да для любой из рассматриваемых функций на одном шаге метода Ньютона, т.е. наибольшую скорость сходимости метода.
Далее оценим скорость сходимости метода Ньютона при реше нии уравнений (7). При этом ограничимся случаем s = — 1 как наилучшим в указанном смысле.
Уравнения (7) запишем в виде |
|
Y(°0 = У 1(а) ~ еу1 = 0. |
(9) |
Здесь ^(а) имеет тот же смысл, что и в |
(8) . Как известно, скорость |
сходимости метода касательных Ньютона определяется формулой
~ |
+ 1 “ Му(Oty —0tn) , |
|
||
где an ,a„ + 1 - |
соответственно |
п-е и п + 1-е приближения к корню |
||
а^ уравнения (9), а |
|
|
||
„ |
У"(Ь») |
|
|
|
Mv = ------ ;----- , |
|
|
||
J |
2 У ' Ю |
|
|
|
ще при установленных выше знаках Y* и Y " имеем |
||||
^ |
^ |
^ |
+ 1 ^ |
П ^ 2. |
При otn , достаточно близких к значению корня соответствующе* го уравнения, справедливы следующие оценки:
MR < |
Л К ) |
|
зх> |
|
|
■Л($„) |
2 [ 1 - К |
- “«)Хтах ]3’ |
|||
|
|||||
Му |
: -2 |
|
|
|
|
|
2^min(^ |
(®7 ®n)Aminl |
|||
|
|
||||
|
|
|
2 |
|
|
MF < |
[1 |
(**«£ |
Xmax |
||
|
&п) ^ т ах] |
||||
где Хтах = max | \ |, Xmin = |
min |
| X,- |, |
как и выше, суть сингу- |
//:X/ФО
212
лярные числа матрицы N. Таким образом, асимптотически
3 |
, |
3 |
„ |
(10) |
MR ^ “ |
^шах» |
Му ^ ~ |
\ n i n > Мр ^.Хj2 |
|
|
|
|
max• |
|
4. Переходим к описанию схемы вычислений при определении корней уравнений (7). Как видно из (7), для реализации одного шага метода Ньютона требуются значения у(ап) иу ' (ап) •Так как решение задачи минимизации (1) эквивалентно решению системы линейных уравнений
(аС +А ТА) иа =А г/, |
(11) |
то значение функции у (а) |
при а = а„ можно вычислить после опре |
деления решения иа системы (И ) по одной из формул, определяю щих функцию у (а ).
Для вычисления значений у ’ (а) потребуются значения производ
ных Р ' (а ), у ' (а) |
или <р' (а ). Используя определения этих функций, |
|
получаем |
|
|
р ' (а) = - а (Mot>£^а) |
|
|
|
р(а) |
|
(ма | |
) |
( 12) |
7'(<*) = |
|
|
7 (а) |
|
|
</(а) = (ма, Сма), |
|
|
где вектор и’а = dujdot можно найти, решая систему |
|
|
(аС + ЛТА) иа = -Сма |
(13) |
с той же матрицей, что и для вычисления иа , но с другой правой частью. Выражение (13) получается прямым дифференцировани ем (11). Производные для функций R (OL) и F(a) легко вывести, используя соотношения (12), а именно
R' (а) = рз (ц р> Сир ) |
У (« ) = |
7(a) |
Р(Р) |
|
|
|
|
(14) |
Р'(а) = - Р2 (ир,Сир),
где принято обозначение р = 1/а. Векторы ир и Мр находятся из уравнений
(pC + A TA)up * A Tf, (рС + А тА)и'р = -Сир . |
(15) |
Для сокращения объема вычислительной работы можно выпол нить преобразования, предложенные В.В. Воеводиным. Матрицу С
методом квадратного корня разложим в произведение КТК = С, где К — верхняя треугольная матрица. Затем матрицу N = АК~1 приведем к двухдиагональному виду: N= QDR, где Q VLR —ортого-
213
нальные матрицы, D — верхняя двухдиагональная. Полагая z a =
= RKua, приведем (И ) |
к системе с трехдиагональной матрицей |
|
(аЕ + D TD)za = D T>p, |
у = QTf |
(16) |
Решение иа исходной системы (И ) |
определяется через решение |
z a преобразованной системы (16) из соотношения
Киа = R Tza.
Легко получить выражения для функций p(a)=\\Dza -<р\\т,
У (л) = II za ||и, ^ (а) = р2 (а) + осу2(а)
идля их производных
р' (а) = -о. (^а>*а)
р(а )
'/ \ |
(га»га-) |
|
г, \ / |
\ |
У (<*) = ---- 7— |
> *(<*) = (Za , ZQ) , |
|||
|
7 (<*) |
|
|
|
где z^ |
определяется путем решения системы |
|||
(a£ + £)r Z))z; = - z a . |
|
(18) |
||
Тогда имеем |
|
|
|
|
|
(У * р ) |
7' (а) = |
(z'g, za) |
|
* '(а ) = А>3 |
|
У (а) |
||
|
Р (Р) |
|
|
F '(a ) = - P 2(2p, zp ),
где р = l/o, zp и Zp определяются последовательным решением трехдиагональных систем уравнений
(рЕ + £>гО) Zp = £ > V (р£ + D r £>) z; = -Zp |
(19) |
с одной и той же матрицей коэффициентов.
Расчетные формулы (17), (18) и (19) требуют значительно меньшего объема вычислений, чем формулы (2 ) - (4 ), (14) и (15).
5. Мы рассматривали параллельно три способа определения па раметра регуляризации, которые для удобства назовем р-, <р- и 7-ме- тодами. Однако эти методы не одинаково эффективны с вычисли тельной точки зрения.
При «^-методе система уравнений решается один раз на каждом шаге метода Ньютона, а в других методах - дважды на шаге.
Далее, приближение к корням соответствующих уравнений про исходит слева, по крайней мере, начиная со второго шага. Для р- и ^-методов это благоприятный фактор, так как обусловленность систем уравнений (15) или (19) будет наилучшей при а = 0; с рос-
214
том а обусловленность ухудшается. Для 7-метода указанное обсто ятельство является неблагоприятным, так как при а = 0 систе мы (11) и (13) или (16) и (18) наиболее плохо обусловлены (воз можна даже вырожденность этих систем). Это может привести к большой потере точности или даже к переполнению разрядной сетки при определении иа или za, если а окажется достаточно близ ким к нулю в начальный момент или в результате выполнения ша га метода Ньютона.
Неравенства (10) показывают, что и в отношении скорости схо димости 7-метод находится в менее благоприятных условиях, чем остальные. Это связано с тем, что величина X^fin существенно за висит от поведения сингулярных значений матрицы N, в то время как величина Х^ах всегда ограничена числом || jVr jV||.
Заметим, что в р- и ^-методах в качестве начального приближе ния к корню можно всегда взять а = 0. Как легко вывести,
R (0+) = H/IU = II * IL. F(0+) = ||/и», = II * и*,.
Rl |
{C~i A Tf, A Tf) _ |
\\NTf\\% |
_ \\DT*Wl |
|
\\f\\m |
WfWm |
IM U ’ |
F' ( 0+) = -dC-'lA Tf, A Tf ) = -II N Tf\\l = -II D T$ ||’ .
Описанный здесь p-метод выбора параметра регуляризации реа лизован на фортране (Численный анализ на ФОРТРАНЕ, вып. 6, 7, под ред. В.В. Воеводина. - М.: МГУ, 1974), при этом использова лись матричные разложения, приводящие исходную регуляризованную систему уравнений к системе (16) с трехдиагональной матри цей. Регуляризованное решение иа восстанавливается лишь после определения корня ар уравнения (2).
В качестве начального приближения к корню всегда берется а = 0. Выполненные расчеты показали неожиданно быструю сходимость p-метода. Так, нормальное решение системы уравнений 8-го поряд ка с вырожденной матрицей было определено за две итерации при относительной точности вычисления корня, равной 10’ 3. Это об стоятельство наводит на мысль, что при задании невысокой точ ности вычисления корня (2) и решении системы (15) с обычной точностью нет необходимости в выполнении соответствующих мат ричных разложений, так как они требуют значительной затраты ма шинного времени, окупаемой лишь при многократном решении регуляризованной системы уравнений. В рассматриваемом случае регуляризованное решение, удовлетворяющее (2), может быть получено за вполне разумное время непосредственно на основе формул (15).
6. Пусть L — матрица размера s X п такая, что квадратичная форма \\Аи ||2 + \\Lu\\2 положительно определена. Это предполо-
215
жение обеспечивает, очевидно, выполнение условия дополнитель ности. Пусть иа реализует минимальное значение функционала
Фа [м] = \\Аи - / | | 2 + a \ \ L u - g \ W u e R ni
где g Е Rs - заданный вектор. Таким образом, мы рассматриваем
конечномерный аналог основной задачи. Полагая а = 1/Л, обозна |
||
чим через их |
вектор, минимизирующий квадратичную |
форму |
X \ \ A u - f \ \ 2 |
+ \ \Lu - g\ \2, u e R n. |
(20) |
Очевидно, их = м1/л. Обозначим p0 Q0 = II А и х - / | | . Тогда р0 (X) -
- |
Р(1/Х) =Р(«). |
100. Если квадратичная форма \\Аи ||2 |
+ || Lu ||2, |
|||
|
Т е о р е м а ' |
|||||
и |
Е |
положительно |
определена, то при всех |
X > 0 функ |
||
ция Ро (X) дважды непрерывно дифференцируема |
и |
Ро (X) < 0, |
||||
Ро"(Х) ^ |
0, г.е. функция Ро (X) выпукла вниз при X > 0. |
|
||||
|
Д о к а з а т е л ь с т в о . Решение их задачи (20) определяется |
|||||
из уравнений Эйлера |
|
|
|
|||
|
L T (Lux - g ) + \ A T (Aux - / ) = 0, |
|
|
|||
вектор |
dux/dX - |
из уравнения |
|
|
||
|
|
|
du х |
|
|
|
|
(LTL + \ A TA ) ------- = - А т (Аих - f ) , |
|
|
|||
|
|
|
dX |
|
|
|
а вектор вторых производных d2ux/dX2 —из уравнения |
|
|||||
|
гг, |
^ |
-------- |
= -2А тА ------- |
|
|
|
(LJL +XAJA) |
dX2 |
dX |
|
|
|
|
|
|
|
|
Используя определение р0 (X) и приведенные уравнения, получаем
Р о ( Х ) Р о ( Х ) = ( л ^ - |
, |
= |
Лг 04их - / ) ) |
= |
|||
= - (( LJL + \ А ТА У 1 А т{Аик - |
/), А 1 {Аик - |
f))< 0 |
(21) |
||||
в силу положительности матрицы |
(LTL + ХАТА ) |
, т.е. p,J(X) < 0 . |
|||||
Далее, еще раз дифференцируя функцию (21), получаем |
|
||||||
Ро (X)Ро (X) + Ро2 (X) = |
(ЛнЛ - / ) |
) + |
|
||||
+ ( /1------- , |
А ------ |
) |
|
|
|
|
|
\ |
dX |
dX |
|
|
|
|
216
Используя неравенство Коши—Буняковского, имеем
|
|
dux |
. |
\ 2 |
|| |
|
dux II 2 |
Ро2 (Х) = Ро(Х) |
( |
d \ |
А и к - / ) |
< |
/4 |
<*Х |
|
Поэтому |
d 2ux |
|
\ |
|
|
|
|
/ |
|
= |
|
|
|||
Ро (X) Ро (Х)>^--^ |
, A T (Aux - f ) J |
|
|
||||
= 2 ((LTL + ХЛ7^ ) ' 1 / l 7^ |
(1 г1 + ХЛ7^ ) - 1 у47 (/4мх - / ) , |
||||||
/47 (Лмх - Л ) > 0 , |
Х>0 . |
|
|
|
|
|
Теорема доказана.
З а м е ч а н и е . Матрица С = L TL может и не быть положитель но определенной. Если С > 0, то в теореме можно считать X > 0.
7. Пусть С > 0 и для определенности вектор g = 0. Выполняя разложения, указанные в п. 4, можно записать
Ро (X) = ||0z* -<pll,
где z x определяется из системы, аналогичной (16),
(E + \ D TD)zK= \ D T<p9 |
(22) |
матрица которой трехдиагональна. Нетрудно видеть, что z x - z Xfx =
—za , ft —1 /X.
Использование (22) для вычисления р0 (X) связано с явным оп ределением решения z x этого уравнения. На самом деле для даль нейшего нам достаточно знать лишь вектор ех = Dzx — <р. Ниже предлагается простой метод, позволяющий непосредственно вычис лять только этот вектор. Это позволяет обеспечивать вычисление с большей точностью значений как функции р0 (X), так и ее произ водной, необходимых при численном решении уравнения
Ро (X) = А. |
(23) |
Именно, умножая (22) на матрицу р , видим, что ех удовлет воряет уравнению
(E + \ DDT)ex = -ip.
Легко |
видеть, |
что |
решение |
г Хщуравнения (22) |
следующим |
образом |
связано' |
с |
е*: z x = |
—\D Tex. Заметим |
также, что |
II (Е + \DDT) “1 || < 1 при всех X > 0.
Последовательные приближения к корню Хд уравнения 1/Ро (X)А 1/А,
равносильного (23), определяются согласно формуле
Ро (Х„) (Д - |
ро (Х„)) |
Vi+i = ХЯ +• |
, и = 0 , 1 , . . . , |
^Po(Xw) |
|
217
где начальное значение Х0, как правило, выбирается так, чтобы О ^Х 0 ^ Х д .
Остановимся на методе выбора начального приближения. Оче видно, в данном случае всегда можно брать Х0 = 0. Однако при незначительных дополнительных вычислительных затратах можно указать более точное начальное приближение к корню.
Идея заключается в построении простой функции г (X), являю
щейся минорантой для функции р0 (X) при X > |
0. Тогда в качест |
||||
ве Х0 естественно взять решение уравнения г(Х) |
= Д. Имеем, оче |
||||
видно, |
|
|
|
|
|
Ро(А)р0 (X) = - ( ( £ + XDDTy l DDTex,ex)> |
|
||||
\D\ |
Po (X). |
|
|
||
> - |
|
|
|||
1 + А ||Л ||2 |
|
\D\ |
|
|
|
Следовательно, Po (X) > |
Po (X), p-0 (0) = || ip ||. Решая |
||||
1 + A || D | |
|||||
|
|
|
|
||
полученное дифференциальное неравенство, получаем |
|||||
Ро (А)>т(А) = До(0)(1 |
+ А || /) II2) ' 1. |
|
|||
Искомое значение Х0 теперь определяется по формулам |
|||||
0 < Х о = Ро (0) - |
А |
II * II - А < А |
|
||
ДНЯ II2 |
ДНЯ II2 |
Д> |
|
||
|
|
все величины в которых легко вычислить.
8.Остановимся еще на некоторых моментах численной реализа
ции. Замечая, что матрица
г* л
D D= ... ,
0
А
где D — верхняя двухдиагональная матрица порядка п, получаем
|
Л |
|
0 |
|
|
DDT | |
|
||
DDT = --------- 1- - . |
|
|||
|
0 |
| |
о |
|
Обозначим |
|
|
|
|
ех = (ег , e „ + i,. . . , е т)т, |
<р = ($Г, <р„+\, • • • ,<рт)Т, |
|||
dex |
_ /jie |
|
den+1 |
dem |
dX |
\c?X |
’ |
dX |
dX ) ’ |
где e и |
суть «-мерные векторы. Нетрудно видеть, что |
|||
|
def |
= 0, i > n + 1, |
||
ei = -Vi, ~ТГ |
218
а векторы е и de/dX определя ются из ’’урезанных” систем уравнений
(£, + Х Ш г )е = -*>,
алг |
de |
= |
|
СЕ + \D D T) — |
|
||
|
d \ |
|
|
|
А |
А |
|
А А„, А |
€ ^ |
ip |
(24) |
= - DDTe =-------- |
Преимущества, достигаемые при переходе к уравнениям (24), очевидны, особенно при т > п. Решение (24) можно осуществить методом прогонки. Вычисление р0 (X) иро(Х) осуществляется сог
ласно формулам |
л |
|
Ро(К)= 2 |
<pj + (е,е), ро (Л)Ро (Х )=(-^- |
.«)• |
Решение z K системы уравнений (22) находится согласно формуле
z A= - XD'e
или непосредственно из системы уравнений
{E + \ D TD) z k =\DT$
при найденном значении параметра X = Хд, определяемом уравне нием (23).
9. Применение модифицированного (обобщенного) критерия выбора параметра регуляризации связано с необходимостью ре шать уравнение
р (a) = hy (а) + 7Я, |
(25) |
где h > 0, т > 0 —некоторые фиксированные величины.
Опишем процесс нахождения корня а уравнения (25) при усло
вии его существования. Зададимся некоторым а 0 > |
0 и заменим |
|
функцию 7 (а) уравнением прямойу г (а) = 7 (а о) + 7 |
СЛо) (ос —ос0)- |
|
Найдем значение |
из условия пересечения графика функции |
|
р(а) с прямой у =у 1(а), т.е. как корень уравнения |
|
|
Р (а) =^1 (а). |
|
(26) |
Из геометрических соображений (рис. 4) видно, что такое зна чение Oi существует и единственно. Для численного решения (26) естественно воспользоваться свойством выпуклости вниз функции
р0 (Х) — р ( 1/Х). Полагая а = 1/Х в |
(26), будем иметь уравнение |
||
Ро (X )- 7 ^ |
= m+h (у (а0) - |
аоУ'(<*<>)), |
(27) |
X |
|
|
|
219
левая часть которого является выпуклой вниз, так к а к 7#(а0) < 0. Поэтому корень Л* уравнения (27) можно эффективно найти, например, методом касательных Ньютона. В качестве нулевого приближения в процессе определения \ \ естественно взять корень уравнения р0 (X) = т. Положим ах = 1/Хi . Очевидно, <*х является корнем уравнения (26).
После определения |
описанный процесс можно повторить. А |
|||
именно, линеаризуем 7 (a) в окрестности |
точки а |
= а*, записав |
||
уравнение касательной у |
= у 2 (а) к |
7 (a) |
в точке |
(a i, 7 (^1) ) , и |
ищем корень а = а 2 уравнения р(а) |
- у 2 (а), предварительно по |
ложив a = 1/Х. В качестве нулевого приближения к корню Х2 урав нения Ро (X) =у2 ( 1/Х) целесообразно взять X i . Ясно, что а2 = 1/Х2.
Аналогично находятся а3, а4, . . . |
Из геометрических соображений, |
||||||
отражающих свойства функций р(а) |
и 7(a), |
с очевидностью сле |
|||||
дует, |
что |
последовательность {а^} |
сходится |
к |
а , причем <хк < |
||
< а |
(к > 1). Можно |
ожидать, |
что |
указанный |
способ решения |
||
уравнения |
(25) не |
потребует |
существенных |
вычислительных, |
|||
затрат. |
|
|
|
|
|
|
Рассмотрения данного параграфа были связаны с конечно мерными операторами А и L. Однако сформулированные резуль таты имеют место и в общем случае. Мы не будем их приво дить из-за излишней громоздкости соответствующих формули ровок.
§ 27. Эвристические методы выбора параметра регуляризации
1. |
Пусть в условиях основной задачи G = Н и L = Е, |
- ре |
шение задачи |
|
|
inf |
Фа [щи*] =Фа [м^1);м -], |
|
“ е я |
|
(1) |
Фа [и;и*] =\\Аи -f\\}? + a\\u - и * \\гн , а > 0 .
Предполагаем также, что оператор А ограничен и определен на всем Н.
'Наряду |
с задачей |
(1), определяющей регуляризованные реше |
|
ния |
рассмотрим задачу минимизации того же функционала |
||
при и* |
т.е. задачу минимизации |
|
|
Ф Л и -.и ^ Ь |
2 |
(2) |
|
Решение задачи (2) |
обозначим через и ^ и назовем { и ™} |
вто |
ричным регуляризованным семейством. Аналогично определяет ся р-ичное регуляризованное семейство.
В целях упрощения дальнейшего анализа выразим р-ичное ре гуляризованное семейство через первичное. Имеем в силу урав-
220