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

Вычислительная математика

УДК 519.612 А. И. Жданов

ОБ ОДНОМ ЧИСЛЕННО УСТОЙЧИВОМ АЛГОРИТМЕ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ НЕПОЛНОГО РАНГА

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

1. Постановка задачи. Рассмотрим произвольную (возможно неполного ранга и несовместную) систему линейных алгебраических уравнений (СЛАУ)

Au = f, A € Rmxn, f € Rm. (1)

Под «решением» СЛАУ (1) в общем случае будет пониматься нормальное псевдорешение и* = A+f, где A+ — псевдообратная матрица.

Пусть U* — вычисленные любым известным прямым или итерационным методом [1,2] на компьютере (в арифметике с плавающей точкой) решения и*. Эти решения различаются по точности. Оценки погрешностей ||U* — и*112 компьютерных решений U* для этих методов зависят от чисел обусловленности матриц A. Спектральное число обусловленности СЛАУ (1) очевидно определяется как К2(A) = oi/оТ, где oi ^ 02 ^ ••• ^ Ok — сингулярные числа матрицы A, k = min(m,n), т = rank (A) ^ k.

Если A — хорошо обусловленная матрица, т. е. её число обусловленности К2 (A) ^ 1/^mach, где ^mach — машинное эпсилон [4], то компьютерные решения U*, вычисленные любыми известными методами [4], по точности сравнимы между собой. В противном случае, известно [1,2], что значительная часть имеющихся рекомендаций по решению плохо обусловленных («2(A) ~ 1/emach) и неполного ранга СЛАУ (т < k) требует нахождения сингулярного разложения (SVD-разложения) матрицы системы или действий, по существу, к этому сводящихся. Это связано с тем, что методы, основанные на SVD-разложении, являются непревзойдёнными (по точности [1]) методами решения этого класса задач. Однако методы, основанные на SVD-разложении, не позволяют эффективно решать СЛАУ большой размерности с разрежёнными матрицами, так как для их применения требуется большой объем оперативной памяти компьютера.

В настоящей работе предлагается новый метод, основанный на регуляризации СЛАУ (1) с использованием эквивалентных расширенных систем. Показано, что этот метод является обратно устойчивым [3] и конкурентно способным (по точности) методам, основанным на SVD-разложении. Кроме того, предлагаемый метод позволяет эффективно использовать математическое обеспечение, созданное для решения больших разрежённых систем на высокопроизводительных параллельных вычислительных системах.

2. Регуляризация на основе расширенных систем. Метод регуляризации А.Н. Тихонова может использоваться для приближенного нахождения нормального псевдорешения и* СЛАУ (1). Именно, в качестве приближения к и* берётся решение регуляризованной нормальной системы

(AT A + a/n)u = ATf (2)

с подходящим значением параметра регуляризации а > 0, где T - знак транспонирования, In — единичная матрица порядка n. Этот подход основан на известном факте [4], что при а ^ 0 решение иа системы (2) стремится к нормальному псевдорешению и* СЛАУ (1). К сожалению, в случае плохо обусловленных или неполного ранга матриц A, при достаточно малых значениях а (а ~ £mach) СЛАУ (2) оказывается очень плохо обусловленной и, как следствие, не может быть решена никаким компьютерным алгоритмом.

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

Регуляризованную нормальную систему уравнений (2) можно записать в виде

Атг — аи = 0,

где г = / — Аи, или

а-1/2Атг — а1/2и = 0. (3)

Объединяя г + Аи = / и (3), получаем систему:

г + Аи = /, а-1/2Аг — а1/2и = 0

или

а1/21т A \ / а-1/2г \ = ( f

AT — а1/21„ ) \ и ) ^ 0

Если обозначить у = а-1/2г и “ = а1/2, то (4) можно записать в виде

(4)

“A? -A, К У) = ( 0 ) - A„x = f, (5)

где

“I? A ^ x = f у \ и 7=/^ f

и

A“ = 1 aT -“In ■ x = и и 0

всегда имеет единственное решение га = (уО;,ит)т, совпадающее с решением регуляризованной нор-

СЛАУ (5) является расширенной системой с квадратной матрицей Аш порядка m + n. Так как а > 0, то и “ > 0. Покажем, что матрица Аш - невырожденная и, следовательно, расширенная СЛАУ енное решение za = (yT uT)T

мальной системы (2)

иа = (ATA + аI„)-1ATf,

и Уа = а-1/2Га, а Га = f — Аиа.

Для этого сформулируем следующую теорему, которая также дает возможность оценить число обусловленности вычислительной задачи (5).

Теорема. Если А = Omxn — нулевая m х n-матрица, то собственными значениями матрицы Аш являются числа и и —и кратностей тип соответственно.

Если А ф ОтХп, то собственными значениями матрицы Аш являются числа ±-^/<т2 + ш2 (г = 1, 2 ,..., т), а также “ и —“ кратностей m — т и n — т соответственно, где т = rank (А) ^ 1.

Доказательство. Так как Аш = AT, то все её собственные значения Aj(A^) € R (i = 1, 2, ..., m + n).

Если А = Omxn, тогда матрица Аш — диагональная, т. е.

Аш = diag (“,...,“, —“,..., —“),

m n

откуда непосредственно следует первое утверждение теоремы.

Пусть А = Omxn и

“Im A z z

\ AT —“In Д V ) =

где z € Rm, v € Rn. Тогда

“z + Av = Az, , ,

ATz — “v = Av. ( )

Если (zT, vT)T = 0 и A = “, то, исключая в (6) z = (A — “)-1Av, получаем

(A — “)-1ATAv — “v = Av

или

Ат Аи = (Л — ш)(Л + ш^ АТА^ = (Л2 — ш2^ (АТА + ш2/га^ = Л2^. (7)

Из (7) видно, что если V = 0, тогда V - собственный вектор и Л2 соответствующее ему собственное

значение матрицы АТА + ш2 1п, т.е. А2 = а2 + ш2 (г = 1, 2, ..., г). Следовательно, числа ±<^/<т2 + ш2

(г = 1, 2, ..., т) являются собственными значениями матрицы Аш.

Если V = 0, а г = 0, то из (7) следует, что шг = Лг и Атг = 0. Так как ёткег Ат = т — т, то ш является собственным значением кратности т — т матрицы Аш.

Если г = 0, а V = 0, то из (7) следует, что —шv = Лv и Av = 0. Так как ёткег А = п — т, то —ш является собственным значением кратности п — т матрицы Аш. □

Так как матрица Аш - симметричная, то из теоремы непосредственно получаем

Следствие. Если т = п или т = п > т, то максимальное сингулярное число матрицы Аш равно

^max(Aw) — \/а2 + u2,

а минимальное сингулярное число сттт(Аш) = ш и спектральное число обусловленности матрицы Аш:

/ 7 \ СГтах(Аш) \/<72 “Ь ^2

К2\Аш) — ------ Т — ------------•

атт(Аш ) ш

Если т = п = т, то <7тах(Аш) = \/а\ +и2, <7тщ(Аш) = д/сг2 + и2, где оп = <тГ7 а

/л ^ а1 + w

К2(Аш) = W -2----------

V an + и

2

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

было устойчивым решением. Таким образом, на основании следствия из теоремы можно получить оценку для числа обусловленности СЛАУ (5), независящую от ранга матрицы А,

К2(Л„) « . (8)

ш

В силу трудоёмкости вычисления ах, учитывая, что ах = ЦА.Ц2 ^ ||А||е, где ||А||е - евклидова матричная норма, вместо неравенства (8) на практике можно воспользоваться неравенством

■ IE +

К2(Аш) ^ К2(АШ) = v

и

3. Решение систем неполного ранга. Если выбрать параметр регуляризации а достаточно малым, то расширенную систему (5) можно использовать для приближенного вычисления численно устойчивых нормальных псевдорешений систем (1) неполного ранга или систем с матрицами A, имеющими неполный численный ранг [2, стр. 75]. В действительности, на практике чаще интересует е-ранг матрицы, который для некоторого малого е определяется по формуле

rank(A, е) = min rank(B).

\\A-B\\2 ^е

Если A € Rmxn — матрица из чисел с плавающей точкой A = fl (A), то разумно считать A численно неполноранговой, если rank (A, е) < min(m,n), где е = emach ЦАЦ2.

Для СЛАУ (1) неполного ранга (rank(A) < k) число обусловленности регуляризованной нормальной системы (2) равно

iT л 2\ _ ^max(ATA + ^2) _ O'2 + U2

K2(AT A + и2) =

■^mln (AT A + и2) и2

а расширенной регуляризованной нормальной системы (РРНС) (Б) —

/ 7 \ Атах(Аш) л/<72 + L02

К2{Аш) — -----— ----------------.

^ln(Aw) и

Таким образом,

K2(AT A + и2) = k2(A^ ).

Если выбрать параметр и из условий: и2 < є^^ и и > є^^, то очевидно, что

fl(ATA + и2/п) =fl(ATA).

Следовательно, при таком значении и компьютерное решение иш РРНС (Б) будет достаточно точным приближением к и*.

Как показано в [4, стр. 210], имеет место оценка

„ „ ^ w2||/||2 w2||/||2 2

Рш—-------------7-------о ^ 5q = const.

ат (ат + и2 а3

Если для решения РРНС (Б) применить любой обратно устойчивый алгоритм [Б], например, основанный на QR-разложении Хаусхолдера [2], то как показано в [3]:

—7= 0(к2(Аш)етжЪ), (9)

где хш - компьютерное решение РРНС (5), ж* = (ут, ит)т, у* = ш-1(/ — Аи*).

4. Тестовые исследования и выводы. В качестве иллюстрации возможностей предлагаемого метода рассмотрим следующий пример. В этом примере задача решалась двумя методами: ЯУБ-методом из пакета МаШЬ и предлагаемым методом на основе РРНС. Для решения РРНС (5) (учитывая, что матрица Аш является симметричной) был использован Ь^Ьт-метод [3], также из пакета МаШЬ. В РРНС (5) параметр ш = 10-15.

Рассмотрим несовместную СЛАУ

/11 1 \ / -94 \

A=

І І І 1 1 1.00000001 V 1 1.0000002 1 J

Є R4x3, / =

106

6.00000003

6.0000004

Є R4. (10)

Точное псевдорешение СЛАУ (10) равно и* = (1, 2, 3)т. Решение вычисленное ЯУБ-методом из пакета МаНаЬ — и = (1.431, 3.668, 1.400)т, ас использованием РРНС (5) — иш = (1.0000000, 1.9999999, 2.9999999)т.

Таким образом, предлагаемый метод на основе РРНС (при соответствующем выборе параметра ш) может даже превосходить по точности ЯУБ-метод. При этом для СЛАУ с т ~ п РРНС превосходит по быстродействию ЯУБ-метод.

*

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Деммель, Дж. Вычислительная линейная алгебра. Теория и алгоритмы [Текст] / Дж. Деммель. —М.: Мир, 2001. — 4З0 с.

2. Голуб, Дж. Матричные вычисления [Текст] / Дж. Голуб, Ч. Ван Лоун. —М.: Мир, 1999. — 548 c.

3. Trefethen, L. N. Numerical Linear Algebra [Text] / L. N. Trefethen, D. Bau. — SIAM, Philadelphia, 1997. — 373 p.

4. Беклемишев, Д. В. Дополнительные главы линейной алгебры [Текст] / Д. В. Беклемишев.—М.: Наука, 1983.— 336 c.

5. Жданов, А. И. Введение в методы решения некорректных задач. Часть 2: учеб. пособие [Текст] / А. И. Жданов. — Самара: Изд-во Самар. гос. аэрокосм. ун-та, 2007. — 79 c.

Самарский государственный аэрокосмический университет Поступила 25.09.2007

им. академика С.П. Королева, г. Самара zhdanov@smr.ru

A. I. Zhdanov

ABOUT ONE NUMERICAL STABLE ALGORITHM FOR SOLVING SYSTEM LINEAR ALGEBRAIC EQUATIONS OF DEFECT RANK

A new method for solving unstable problems that can be reduced to arbitrary systems of linear algebraic equations (which may not be of full rank or may be inconsistent) is examined. This method is based on the reduction of regularization of normal system equations to an equivalent augmented regularization of normal system equations.

Samara State Aerospace University, Samara, Russia zhdanov@smr.ru

Received 25.09.2007