Об одном оптимизационном алгоритме, вычисляющем наибольшее сингулярное число вещественной матрицы Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 519.6 А. Н. Борзых

ОБ ОДНОМ ОПТИМИЗАЦИОННОМ АЛГОРИТМЕ,

ВЫЧИСЛЯЮЩЕМ НАИБОЛЬШЕЕ СИНГУЛЯРНОЕ ЧИСЛО ВЕЩЕСТВЕННОЙ МАТРИЦЫ

Введение

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

1. Применение степенного метода к матрице ЛЛТ (см., например, [1, 2]). Возможно как явное перемножение матриц Л и ЛТ с последующим применением степенного метода к результирующей матрице ЛЛТ, так и «неявное перемножение», когда на каждом шаге текущий вектор сперва умножается на ЛТ, а затем на Л. Также возможно применение степенного метода к матрице (А АО, что является строго эквивалентным его применению к матрице ЛЛТ.

2. Приведение заданной вещественной матрицы к двухдиагональной форме с последующим применением метода бисекции (см. [1]).

3. Применение метода релаксации отношения Рэлея к матрице ЛЛТ (с явным или неявным перемножением матриц) или к матрице (А АТ ). Возможно использование либо покоординатного спуска (когда на каждом шаге изменяется лишь одна компонента текущего вектора), либо градиентного спуска (см. [2]).

Степенной метод при близких старших сингулярных числах может иметь сколь угодно медленную сходимость, что делает его крайне «ненадежным» и редко используемым на практике. Приведение матрицы к двухдиагональной форме имеет вычислительную сложность 8/3п3 флопов и может быть «разумным» лишь для матриц относительно небольших размерностей. Метод релаксации отношения Рэлея, можно полагать, наиболее оптимален как в плане его «надежности», так и в плане применимости к матрицам любой размерности.

В работе [3] был предложен новый оптимизационный алгоритм, вычисляющий наибольшее собственное число заданной симметричной матрицы. В данной работе предла-

© А. Н. Борзых, 2008

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

Схема алгоритма

1. Если Др > Ад, где Др = шах(рк) — тт(рк), Ад = тах(дк) — тт(дк), рк =

к к к к

£аЬ (строковые суммы матрицы), Як = агк (столбцевые суммы матрицы), то

г г

выполняем шаг 1.1 (преобразование строк); в противном случае выполняем шаг 1.2 (преобразование столбцов).

(1.1) Преобразуем в матрице А две строки:

А = Ta A,

(1.1)

( 1

где та =

1

s = sin (а), c = cos (a), i —номер строки с минимальной строковой суммой, j — номер) строки с максимальной строковой суммой, а вычисляется из условия s (А) max, где s (А) = Выведем явную формулу для а.

п “г,з

4 / 4 / i,j Сумма элементов строк i и j после преобразования имеет вид Pi (а) = cos (а) • pi (A) - sin (а) • pj (A) ;

Pj ^A^ = sin (а) • pi (A) + cos (а) • pj (A).

Получим производную s' ^A^ по а:

.'(Л) = 1(Р;(Л)+Р;(Д)) =

= — (— sin (а) ■ pi (А) — cos (а) • pj (А) + cos (а) ■ pi (А) — sin (а) • pj (А)).

, (Pi - Pj\ ,

«2 = arctg -------------- + 7Г.

VPi + Pj/

Решая в' = 0, получаем

, (рг — рЛ

а\ = ап^ -----------

\рг + рз)

Одно из этих решений обеспечивает глобальный максимум в (А), другое глобальный минимум в (А).

c

—s

i

1

s

c

(1.2) Преобразуем в матрице A два столбца:

A = AT”, (1.2)

где i — номер столбца с минимальной столбцевой суммой, j — номер) столбца с максимальной столбцевой суммой, а вычисляется из условия s ^A^ —U max.

2. Пусть е — требуемая точность вычислений. Если Д ^А^ < е, где Д ^A^ =

= max ^Др ^A^ , Д^ ^A^ ^, то прекратить вычисления и в качестве приближенного значения наибольшего сингулярного числа исходной матрицы А взять s ^А^. В противном случае положить A = A4 и перейти к пункту 1.

Замечание 1. В представленном алгоритме а вычисляется из следующего условия s ^A^j —U max. Легко проверить, что при таком выборе а оказывается справедливым

pi ^A^ = pj ^A^, если выполняются преобразования строк, и qi ^A^ = qj ^A^, если выполняются преобразования столбцов. То есть данный выбор а не только максимизирует сумму всех элементов текущей матрицы, но и «уравнивает» суммы элементов двух строк или двух столбцов.

Доказательство линейной сходимости предложенного алгоритма

Утверждение 1. Преобразования (1.1) и (1.2) не меняют сингулярных чисел матрицы A.

Доказательство. Пусть А — сингулярное число A, т. е. А2 собственное число AAT. Пусть A = T1AT2, где Ti, T2 —произвольные ортогональные матрицы. Заметим, что AAT и AAT суть подобные матрицы:

AAT = T1AT2T2tAttT = Ti (AAT) tT .

Следовательно, А2 есть собственное число AAT, т. е. А — сингулярное число A.

Следствие. Для всякого к справедливо ||Afc|| = ||Ao||, где Ak —матрица A, получаемая на шаге алгоритма к.

Утверждение 2. Имеет место неравенство

s (A) <||A||.

Доказательство. Оценивая квадрат нормы матрицы через отношение Рэлея, получаем

12 п лТ /,ц^ ((АтА)е,е^ (Ае,Ае) (Ае, Ае) (е, е) ^ 1 u л ^ 2

(е>е) (е>е) (е,е)2

Ш = 11^ All > " ^ 7 ^ ’ > - Ue, е)Г = s2 (А),

где e = (1, 1 . . . , 1)T.

(e,e)

T

|s2 - Л2| < 2||A|| • Д, где А — одно из сингулярных чисел матрицы А.

Доказательство .

Определим вектора р, 4:

р = Ae — s (A) e, 4 = ATe — s (A) e,

то есть

Ae = s (A) e + p, ATe = s (A) e + 4.

Отметим, что

||р||то = max |pi| = max |pi — s (A) | < Др < Д,

|І4І|ТО = max |qfj| = max |q* — s (A) | < Д< < Д,

\\p\\ < VnA, ||(?|| < л/пА.

Вычислим ATAe:

ATAe = AT (s • e + p) = s2 • e + s • 4 + ATp = s2e + f,

где f = s • 4 + ATp.

Оценим ||f ||:

\\f\\<\s\-\\q\\ + \\AT\\-\\p\\<2M\A\\^

Применим к симметричной матрице ATA теорему об оценке невязки ([4, §4.5], [5, §7]). Согласно теореме

V*.3 Л2 є <т (АТА) : А2 - „| < ~ Н|.

||х||

Возьмем x = e, ^ = s2, тогда

|л2 2| ||ATAe — s2e|| ||f|| л

A2 - s2 < i!-------0----" = Ш < 2 A • Д.

||e|| y/n

Утверждение 4. Имеет место предельное соотношение

* i А \ k *

3s : s (Afc) —> s .

Доказательство. По определению a (А) = ~^2pk = ~^2чк- На каждом шаге

k k

алгоритма р* ^Aj + p (Aj > p* (A) + pj (A) или q* (Aj + qj ^A^ > q* (A) + qj (A).

Следовательно, s ^A^ > s (A), т. е. s — монотонно-возрастающая функция. С другой стороны, s (A) < IIA|, т. е. s —ограниченная сверху функция. Значит, существует s* : s (Ak) k—^ s*.

в(л) - 5 (А) > СхД2,

гдв С1 = 4ТВД-

Доказательство. Рассмотрим один шаг алгоритма. Будем считать, что на данном шаге выполняется преобразование строк (доказательство для столбцов проводится аналогично). Преобразование (1.1) сохраняет неизменным сумму квадратов строковых сумм: р2 + р2 = Р + Р/. С другой стороны, преобразование (1.1) обеспечивает равенство сумм элементов строк г, _?’: р = р^. Решая систему из этих двух уравнений, получаем

¡P2i+Pj

Pi = Pj = тогда

S(A) -s(A) = (Я-+Я.)-(».+»■)= fc+Pi? -b>i+Pi? > (Pi+Pjf -(Pi+Pjf y W ( j [Рі+Рз) [Рі+Рз) (рг+рэ) + (рг+рэ) - 2IIAIU + 2IIAIU -

1

>

(pi + pj)2 - (pi + pj)2

4^||А|| L l J 4а/п||А||

Утверждение 6. Имеет место предельное соотношение

A(Afc) 0.

Доказательство .

По утверждению 4 существует s*: s (Ak) —2? s*. Следовательно, s ^A^ — s (A) 2 0. Тогда по утверждению 5 Д 2 0.

Теорема 1. Алгоритм сходится к одному из сингулярных чисел:

s (Ak) —? A.

Доказательство. По утверждению 3 |s2 — A2| < 2||A|| • Д, где A — одно из сингулярных чисел матрицы А. По утверждению 6 Д 2 0. Следовательно, s2 2 A2, s 2 A.

Лемма 1. Пусть векторы х, у: ||х|| = ||у||. Пусть в —острый угол между векторами х, у,

а (х,У)

COS# = -г—г,-г—77,

Тогда существует ортогональная матрица V, такая что

их = у; \\и - /|| = 2 • віп ^

где I — единичная матрица.

Доказательство. См. [3].

x

в < sin 2в

при в £ [0, |].

Доказательство тривиально.

Введем определения.

1. Если x — собственный вектор A, отвечающий собственному числу А, то (x, А) называем собственной парой A.

2. Если l, r — левый и правый сингулярные вектора A, отвечающие сингулярному числу А, то (l, r, А) называем собственной тройкой A.

Утверждение 8. Пусть (l*, r*, А^) — собственная тройка A*. Пусть A = A* — F, || F У < У Ay. Пусть А1,...,А„ —сингулярные числа A. Тогда у матрицы A имеется собственная тройка (l, r, Ат), удовлетворяющая условиям

|« - Am|< 3||A||||F

1й„29,<МИ; 1»т2»2<МИ,

2 & (А, Ат) 2 и;(А, Ат)

где = Z (1, 1*) — угол между векторами 1 и I*, #2 = ^ (г, г*) — угол между векторами г и г*, ш (A, Am) = min |A2 - Amm|-i=m

Доказательство. Установим связь между матрицей AAT и матрицей A*A*T:

AAT = (A* - F) (A* - F)T = A* A*T + E, где E = -FA*T - A*FT + FFT- Для нормы E справедлива оценка ||E|| < 3||A||||F||.

Симметричную матрицу AAT в отношении симметричной матрицы A*A*T будем считать возмущенной:

AAT = A* A*T + E.

Отметим, что (l*,Am2) есть собственная пара A*A*T (так как 1* есть левый сингулярный вектор A*, отвечающий сингулярному числу Am). Тогда из теории возмущений собственных чисел и собственных векторов симметричных матриц (см. [1, с. 214]) следует, что существует собственная пара (1, A^m) матрицы AAT, удовлетворяющая условиям

|Am2 - Am|<||E||< 3||A||||FII;

Idna., < . ,1? < 3|I'4|II|F»

2 тт(А^ - А2) тт(Ат - А2) ’

г=т г=т

где 01 = / (1, Г), А2 — собственные числа ААТ.

Доказательство для правых собственных векторов аналогично.

Утверждение 9. Пусть А* — сингулярные числа матрицы А. Пусть в ^ А& (по теореме 1). Начиная с некоторого шага алгоритма, матрица А имеет левый сингулярный вектор 1 и правый сингулярный вектор г, удовлетворяющие усло-

виям

/ (1,е) < С2Д, / (г,е) < С2Д,

где = и{АМ)1 Ш(А’ Л*) ~Хк\'

Доказательство. Определим матрицу Е следующим образом:

/и = 2в (А) - р1 (А) - 41 (А) ;

/¿1 = в (А) — р®(А), г > 1; /и = в (А) — 9®(А), г > 1;

/у =0, г, з > 1.

Отметим следующие свойства Е:

||Е|| < пД;

П

Р1(Е) = 2в (А) - р1(А) - 91(А) + £ (в (А) - 9®(А)) = в (А) -р1 (А) ;

¿=2

р® (Е) = в (А) - р® (А), г > 1;

П

91 (Е) = 2в (А) -Р1 (А) - 91 (А) + £ (в (А) -р®(А)) = в (А) - 91 (А);

¿=2

9® (Е) = в (А) - 9® (А), г > 1.

В силу того, что Д ^ 0, можно утверждать, что с некоторого шага алгоритма

НЕ|| < ||А||.

Определим матрицу А* следующим образом:

А* = А + Е

Отметим следующие свойства А*:

р® (А*) = р® (А + Е) = р® (А) + р® (Е) = р® (А) + в (А) - р® (А) = в (А);

9® (А*) = 9® (А + Е) = 9®(А) + 9® (Е) = 9® (А) + в (А) - 9® (А) = в (А) ,

то есть все строковые и столбцовые суммы А* равны в (А). Тогда е одновременно и правый, и левый сингулярный вектор А*:

А*А*Те = в (А)2 е, А*тА*е = в (А)2 е.

Матрицу А по отношению к матрице А* будем считать возмущенной:

А = А* - Е, ||Е|| < пД.

Пусть (е, е, в2) —есть собственная тройка А*. Тогда из утверждения 8 следует, что существует собственная тройка (1, г, Ат) матрицы А, удовлетворяющая следующим условиям:

1) |в2 - Ат|< 3||А||||Е|| < 3п||А|| • Д;

2) і 8ІП20! < • А, где вх = / (е, 0;

3) ^ 8Іп202 < ^'(1'лГ) < • А’ гДе 02 = ^ (е, г).

Так как Д ^ 0, из условия 1 вытекает, что в2 ^ Ат, в ^ Ат. Однако по условиям утверждения 9 в ^ А&. Следовательно к = т, А& есть Ат.

Из условий 2 и 3 вытекает, что 01,02 ^ 0, и можно утверждать, что, начиная с некоторого шага алгоритма, 01, 02, € [0, |-]. Применяя утверждение 7, получаем

1 Зп||А|| 1 Зп||А||

2 1 - ш(А,Хк) ' ’ 2 2 - ш(А,Хк) ' •

Тогда

/(,’'> = ^Т1ТТД = С"Д; ¿М = <>2<-?Ш-.ь = с2ь.

ш (А, Ай) ш (А, Ай)

Утверждение 10. Пусть и — ортогональная матрица, ||£/ — /|| = х < Тогда и и иТ представляются в виде

и = I + 5 + Р, ит = I - 5 + Рт,

где 5 —кососимметричная, ||51| < 2х, ||Р|| < 2х2.

Доказательство. См. [3].

Утверждение 11. Начиная с некоторого шага алгоритма,

А2 (А) > Сз (^і - в (А)),

где С3 = 6С22||А||, и>1 — точное значение одного из сингулярных чисел матрицы А.

Доказательство. По утверждению 9, начиная с некоторого шага алгоритма, матрица А имеет левый сингулярный вектор / и правый сингулярный вектор г: А (/, е) < С2А, А (г, е) < С2А. Будем считать, что ||/|| = ||г| = ||е|.

По лемме 1 существуют ортогональные матрицы Ь и Д:

е = Дг,

\\Ц-1\\ = 2*т(Щ^-) <2*п(Щ < С2Д;

= Ь/,

¿|£>)<2ЙП(^Л

\Ь - ІII = 2віп ( —^— ) < 2віп ( < С'гА.

В силу того, что Д ^ 0, можно утвержать, что ||Д - 1|| ^ 0, ||Ь - 1|| ^ 0. Будем считать, что

1|й-Л1<^ \\ь-1\\<1-.

Представим A в ее сингулярном разложении: A = UWVT, где U — матрица левых сингулярных векторов (один из столбцов — вектор 1), V — матрица правых сингулярных векторов (один из столбцов — вектор г), W — диагональная матрица сингулярных чисел.

Определим матрицу A следующим образом:

A = LART = LUWVT RT = LUW (RV )T = U0WV0T,

где Uo = LU, причем один из столбцов Uo есть вектор e (т. к. 1 есть столбец U; L1 = e), Vo = RV, причем один из столбцов V0 есть вектор e (т.к. г есть столбец V; Rr = e).

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

Ae = U0WV0Te = U0W (1, 0, - - - 0)T = U0 (wb0, - - - 0)T = w1e,

где wi —одно из собственных чисел A (совпадает с сингулярным числом). Умножив слева данное уравнение на eT, получаем

T л T

e Ae = e wie-

Тогда точное значение сингулярного числа матрицы A (и, следовательно, A) определяется равенством

W\ = —eTÄe = s .

По утверждению 10 возможно представление

LT = I + S1 + Pi, ||Si||< 2С2Д, ||Pi|| < 2С|Д2;

R = I + S2 + P2, HS2II < 2С2Д, HP2H < 2С22Д2,

где Si, S2 —кососимметричные.

Отметим, что

s fs'iA') = — eTS\Äe = —-eTS\e = 0, s fAsO = 0.

V / n n V /

Оценим wi - s (A):

wi - s (A) = s ^A^ - s (A) = s ^A^ - s Ar^ - s ^(I + Si + Pi) A (I + S2 + P2^ =

= s ^A^ - s (,4) - s (^A) - s (As2) -s (s^As^ - s (PlAR - s (ltAp2 j <

< IIS1H • ||A|| • ys2y + ||Pi|| • ||Ar|| + ||LTA|| • yp2y < бс|д2||Ац < 6С2||А||Д2.

Теорема 2. Алгоритм имеет линейную скорость сходимости.

Доказательство. По утверждению 5

s(A) - s (A) > С1Д2-

По утверждению 11

А2 (А) > Сз (^і - в (А)).

Тогда

в ^А^ — в (А) > С1С3 (^і — в (А)).

то есть

в (А - в (А)

У У >С71С73.

аді - в (А)

Отметим, что и>і - в (А) —погрешность на текущем шаге, и>і - в ^А^ —погрешность на следующем шаге. Для их отношения справедлива оценка

и>і - в ^А^ и>і - в - в (А) + в (А) в - в (А)

--------, д ' = -------------;— -------- = 1 — ----------— < 1 — Сі Сч.

аді - в (А) ^і - в (А) ^і - в (А)

Таким образом, показано, что сумма всех элементов матрицы, максимизируемая на каждом шаге алгоритма, достигает предельного значения с линейной скоростью сходимости. Предельное значение есть пАі, где Аі, строго говоря, одно из сингулярных чисел А (требуется, чтобы оно было простым). Может ли алгоритм сойтись не к максимальному сингулярному числу?

Матрицу А, получаемую на к-ом шаге, можно выразить через начальную матрицу А:

Ад = Тк ... ТіА5і... ,

где Ті... Тд, 5і ... , —некие ортогональные матрицы, порождаемые алгоритмом.

Введем обозначения: Хд = Тд... Ті, Уд = 5і ... . Тогда

Ад = X АУД,

в(Ак) = - (Аке,е) = - (АУке,ХІе) = (Аук,хк), п п

где хк = ук = -^=Уке— некие векторы, порождаемые алгоритмом, причем

11хй У = 1Ы1 = і.

Следовательно, можем говорить, что алгоритм ищет максимум функционала в* (ж, у) = (Ау, ж) по всем х, у: ||ж|| = ||у| = 1. При каких условиях функционал в* (ж, у) имеет максимумы? Есть у него локальные максимумы, отличные от глобальных?

Если у фиксировано, то максимум по ж обеспечивается при ж = Ау/||Ау|| (т. е. когда вектора ж и Ау коллинеарные). Поэтому

т«х(8*(»у!/» = тах (.,* (р^,»)) = тах (^¡^) = тах(Ц.ЗД),

т. е. мы можем утверждать, что алгоритм ищем максимум функционала в** (у) = ||Ау|| по всем у: ||у|| = 1. Про данный функционал хорошо известно, что он имеет лишь одну точку максимума, значение в которой есть максимальное сингулярное число. Локальных максимумов, отличных от глобального, в** не имеет, но имеет точки минимакса.

Таким образом, показано, что представленный алгоритм может интерпретироваться как релаксационный алгоритм, вычисляющий максимум функционала в**. Из теории

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

Теоретически может, но про сходимость подобного рода можно утверждать, что она является неустойчивой. Действительно, малейшее возмущение точки минимакса (например, из-за ошибок округления) может приводить как к увеличению целевой функции, так и уменьшению. Если возмущение приводит к увеличению целевой функции, то начинается новая сходимость к другой точке минимакса или точке локального максимума.

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

Замечание 2.

Конечно можно придумать исключительные примеры, для которых алгоритм «застрянет» на не максимальном сингулярном числе. Например, если исходная матрица еще до начала выполнения алгоритма имеет равные строковые и столбцовые суммы, то алгоритм остановится на первом же шаге (хотя «найденное» сингулярное число не обязательно максимальное). В качестве практического решения, предотвращающего поиск «нежелательного» сингулярного числа (т. е. не максимального), можно предложить выполнение нескольких «случайных» преобразований подобия и последующий повторный запуск алгоритма. Если по результатам двух или более запусков получается один и тот же результат, то с большой уверенностью можно утверждать, что вычислено именно максимальное сингулярное число.

Замечание 3.

Очевидно, что предлагаемый алгоритм можно рассматривать как релаксационный алгоритм. Целевой функцией выступает сумма всех элементов текущей матрицы, а пространством, в котором ведется оптимизация, — пространство ортогональных матриц:

Е (и, и2) = ^(иТЛи2)у = еТиТАи2е ^ тах. м

В теории методов оптимизации известно множество теорем, представляющих достаточные условия линейной скорости сходимости релаксационных методов. Тем не менее, применение какой-либо из существующих теорем в отношении рассматриваемой задачи оказалось чрезвычайно трудоемким (быть может, невозможным) по следующим причинам:

1) точка экстремума является вырожденной;

2) множество, на котором ищется решение, является нелинейным.

Заключение

В статье предложен новый алгоритм, вычисляющий максимальное сингулярное число заданной вещественной матрицы. Показана линейная сходимость предложенного алгоритма. Теоретическое сравнение с ранее известными алгоритмами (сравнение констант, определяющих скорость сходимости), а также проведение численных экспериментов планируется сделать предметом дальнейшей работы и последующих публикаций.

1. Деммель Дж. Вычислительная линейная алгебра / Пер. с англ. Х. Д. Икрамова. М.: Мир, 2001.

2. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. СПб.: Лань, 2002.

3. Борзых А. Н. О сходимости одного оптимизационного алгоритма расчета наибольшего собственного значения // Записки научных семинаров ПОМИ. 2007. Т. 346. С. 5-20.

4. Парлетт Б. Симметричная проблема собственных значений / Пер. с англ. Х. Д. Икрамова, Ю. А. Кузнецова. М.: Мир, 1983.

5. Икрамов Х. Д. Несимметричная проблема собственных значений. М.: Наука, 1991.

Статья поступила в редакцию 18 мая 2008 г.