Журнал вычислительной математики и математической физики, 2019, T. 59, № 10, стр. 1641-1647

ВВЕДЕНИЕ

В настоящей работе рассматривается задача исключения для системы линейных алгебраических уравнений (СЛАУ), которая в простейшей ее постановке (см., например, [1]) состоит в том, что нужно вычислить линейную форму от решения СЛАУ; само решение при этом не требуется. Настоящая работа примыкает к работам [2]–[4]. Существенное отличие от этих работ состоит в том, что рассматривается в общем случае переопределенная СЛАУ, которая может быть несовместной. При этом ее решение понимается в смысле теории наименьших квадратов. Таким образом, в работе задача исключения формулируется и решается для СЛАУ, полученной из исходной системы после 1-й трансформации Гаусса. В работе формулируются условия, при которых значение линейной формы не зависит от выбора решения трансформированной системы, если она неоднозначно разрешима.

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

Кроме того, этот метод имеет значительное преимущество перед другими способами решения рассматриваемой задачи исключения в том случае, если решается возмущенная задача (например, исходная совместная СЛАУ становится неразрешимой в результате внесения в исходные данные достаточно малых погрешностей). В этом случае, следуя идеям регуляризации А.Н. Тихонова (см. [9]) в работе [10] для метода Крейга предложен критерий прекращения итераций, играющий роль регуляризатора и позволяющий получить приближенное решение, близкое, что очень важно, к решению исходной невозмущенной задачи, а не той, которая решается непосредственно.

Основные теоретические положения метода продемонстрированы в работе на простом модельном примере.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим систему линейных агебраических уравнений вида

где $A \in {{C}^{{m \times n}}}$ и $b \in {{C}^{m}}$ заданы, $x \in {{C}^{n}}$ неизвестно. При $m > n$ система является переопределенной. Эта система может не иметь решения. В таком случае будем использовать метод наименьших квадратов для нахождения вектора, дающего наименьшую невязку уравнений системы (1.1). Для этого рассмотрим вспомогательную СЛАУ, полученную с помощью первой трансформации Гаусса (см., например, [1]) которая всегда имеет решение, возможно не единственное.

Здесь и далее $A{\kern 1pt} {\text{*}}$ обозначает сопряженную матрицу (в вещественном случае – просто транспонированную).

Каждое решение системы (1.2) (и только эти решения) дает наименьшее из возможных значение скалярного квадрата вектора невязки, т.е. минимум функционала $(Ax - b,\;Ax - b)$.

Для пары векторов в ${{C}^{m}}$ и ${{C}^{n}}$ будем использовать стандартное скалярное произведение $(u,{v}) = {v}{\kern 1pt} {\text{*}}u$. Как обычно, модуль вектора означает $\left| u \right| = \sqrt {(u,u)} $.

Если система (1.1) несовместна, то векторы, минимизирующие указанный функционал, называются ее псевдорешениями.

Нам понадобится следующая известная (см., например, [4])

Лемма. Если СЛАУ (1.1) имеет хотя бы одно решение, то СЛАУ (1.1) и (1.2) эквивалентны, т.е. совокупности решений этих систем совпадают.

В работе рассматривается следующая задача, называемая задачей исключения, формулируемая для системы (1.2).

Пусть задан вектор $f \in {{C}^{n}}$. Требуется вычислить значение линейной формы

от решения $x$ СЛАУ (1.2). При этом само решение $x$ нам не нужно.

Будем рассматривать такие векторы $f \in {{C}^{n}}$, для которых значение $\sigma $ не зависит от того, какое именно из решений системы (1.2) выбрано в случае ее неоднозначной разрешимости. Для того, чтобы $\sigma $ не зависело от этого выбора, необходимо и достаточно, чтобы вектор $f$ был ортогонален $\ker A{\kern 1pt} {\text{*}}A$ – ядру оператора с матрицей $A{\kern 1pt} {\text{*}}A$, т.е. чтобы для всех решений $z \in {{C}^{n}}$ СЛАУ

выполнялось соотношение Тогда, как известно, (см., например, [11]) $f \in \operatorname{Im} A{\kern 1pt} {\text{*}}A$, т.е. имеет место представление $f = A{\kern 1pt} {\text{*}}A{v}$.

С другой стороны, согласно приведенной лемме СЛАУ (1.4) эквивалентна СЛАУ

поэтому вектор $f$ должен быть ортогонален $\ker A$ – ядру оператора с матрицей $A$.

Таким образом, значение $\sigma $ не зависит от выбора решения $x$ СЛАУ (1.2) тогда и только тогда, когда вектор $f$ имеет вид

где $u \in {{C}^{m}}$ – некоторый вектор, представимый, как следует из вышеизложенного, в виде для некоторого ${v} \in {{С}^{n}}$.

Действительно, в этом случае $\sigma $ может быть вычислено следующим образом:

Очевидно, это значение не зависит от выбора решения $x.$

2. МЕТОДЫ РЕШЕНИЯ

Рассмотрим некоторые методы решения рассматриваемой задачи. Итак, мы имеем переопределенную ($m > n$) СЛАУ (1.1), возможно, несовместную. В этом случае она решается методом наименьших квадратов, так что задача исключения формулируется для системы (1.2). В этом состоит существенное отличие поставленной в работе задачи исключения от аналогичной задачи и методов ее решения, рассмотренных ранее в работах [2]–[4].

Естественный подход к нахождению $\sigma $ состоит в том, что вычисляется матрица $A{\kern 1pt} {\text{*}}A$ и решается СЛАУ (1.2), особенно, если матрица исходной системы (1.1) полного ранга, т.е. $r = {\text{rank}}\,A = n$. В этом случае система (1.2) имеет единственное решение. Используя это решение в явном виде, находится $\sigma $ по формуле (1.8). При таком подходе приходится иметь дело с матрицей $A{\kern 1pt} {\text{*}}A$, вычисление которой затратно при больших размерах задачи.

В настоящей работе мы уделяем внимание методам решения сформулированной задачи, не требующим использования матрицы $A{\kern 1pt} {\text{*}}A$, а оперирующим только с матрицей $A$.

Один из таких методов – это применение к решению СЛАУ (1.1) метода сопряженных градиентов после первой трансформации Гаусса (см. [1]). При этом сама трансформация не делается (т.е. матрица $A{\kern 1pt} {\text{*}}A$ не вычисляется). Этот метод реализуем при любом выборе начального приближения и при $r = n$ дает единственное решение СЛАУ (1.2), которое является для СЛАУ (1.1) наилучшим в смысле теории наименьших квадратов. Затем используется формула (1.8). Однако при $r < n$ метод сходится к решению, зависящему от начального приближения, и поскольку на практике условие (1.5) для вектора $f$ не проверяется, значение $\sigma $ также может зависеть от выбора этого приближения, т.е. от выбора решения задачи.

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

где $u \in {{C}^{m}}$ – искомый вектор. Допустим, что матрица $\,A$ может быть неполного ранга, т.е. $r \leqslant n$. Система (2.1) является недоопределенной. Если вектор $f$ выбран указанным выше образом, то она разрешима и имеет $(m - r)$-мерное линейное многообразие решений. Выделим из него решение, ортогональное подпространству решений соответствующей однородной системы, т.е. решение, минимальное по норме. (Известно, что в случае разрешимости системы (2.1) такое решение существует и единственно.) Используем для этого второе преобразование Гаусса (см. [1]). А именно, будем искать векторы $u$, представимые в виде (1.7). Относительно вектора ${v} \in {{C}^{n}}$ получается задача Уравнение (2.2) имеет и притом единственное решение, которое позволяет получить интересующее нас решение $u$ задачи (2.1) по формуле (1.7).

Однако можно находить решение $u$ непосредственно из системы (2.1) без ее предварительного преобразования. Будем для этого использовать метод Крейга, т.е. метод сопряженных градиентов после второй трансформации Гаусса (см. [1]). Сама эта трансформация в методе не делается.

Поскольку на практике, как уже было сказано, условия (1.4), (1.5) не проверяются, то условием разрешимости СЛАУ (2.1) является условие реализуемости метода Крейга, т.е. условие успешного завершения процесса решения. В работе [6] показано, что в случае разрешимости задачи (2.1) метод Крейга дает последовательность элементов, сходящуюся к искомому минимальному по норме решению.

При использовании метода Крейга мы получаем приближенное решение уравнения (2.1), представимое в виде

где ${{u}_{0}}$ – некоторое начальное приближение, $y \in {{С}^{n}}$ – некоторый вектор. А тогда искомое значение $\sigma $ вычисляется следующим образом (см. (1.8)): Отсюда видно, что если $x$ является решением не только СЛАУ (1.2), но и СЛАУ (1.1) , или если ${{u}_{0}}$ представимо в виде ${{u}_{0}} = A{{{v}}_{0}}$, в частности, ${{u}_{0}} = 0$, то $\sigma $ вычисляется по формуле (1.8). Действительно, в этом случае Поскольку система (1.1) может быть несовместной, то из вышеизложенного следует

Вывод. Если метод Крейга для решения СЛАУ (2.1) реализуем, то значение $\sigma $ не зависит от выбора решения СЛАУ (1.2) и при ${{u}_{0}} = A{{{v}}_{0}}$ (в частности, ${{u}_{0}} = 0$) мы получим то решение СЛАУ (2.1), которое дает возможность вычислить значение $\sigma $ по формуле (1.8).

3. ФОРМУЛЫ МЕТОДА

Рассмотрим формулы метода решения системы (2.1), т.е. метода Крейга. Отметим, что этот метод относится к методам сопряженных направлений, общим свойством которых является построение ортогональной в некоторой выбранной метрике системы векторов, по которой раскладывается искомое решение. В конечномерном пространстве такие методы являются конечными. При наличии ошибок округления, особенно в плохо обусловленных задачах больших размеров, процесс построения указанной системы векторов становится бесконечным, т.е. эти методы становятся итерационными. В этом случае важной характеристикой метода является его устойчивость относительно накопления погрешностей в процессе итераций. В [7], [8] предложены модификации для некоторых методов такого типа. Модифицированные методы при точных вычислениях эквивалентны исходным, а при вычислениях с погрешностью гарантируют устойчивость итерационного процесса, в котором неизбежные ошибки округления, в отличие от исходных методов, не приводят к расходимости и разбалтыванию.

Приведем формулы модифицированного метода Крейга применительно к системе (2.1):

где ${{u}_{0}}$ – начальное приближение, ${{r}_{0}}$ – начальная невязка;

для $k \geqslant \,1$ имеем

а коэффициенты ${{\alpha }_{k}}$ вычисляются по следующей формуле, которая составляет суть модификации из [7], [8]:

Поскольку само решение в рассматриваемой задаче нам не нужно, то в приведенном алгоритме метода вычисление решения можно исключить и получить рекуррентную формулу для нахождения $\sigma $. Это тем более удачно, так как обычно значение функционала – величина, более устойчиво зависящая от исходных данных, чем решение линейной системы. Типичный пример – это СЛАУ, полученная аппроксимацией интегрального уравнения Фредгольма I рода, а функционал – умножение на гладкую функцию. При этом погрешность приближенного решения – это быстроколеблющаяся функция, функционал от которой может быть мал. Тогда последовательность приближенных значений функционала сходится намного быстрее, чем последовательность приближений к решению этого уравнения (см. расчеты в [4]).

Итак, в результате несложных преобразований получим следующие расчетные формулы:

для $k \geqslant \,1$ имеем

4. О СВОЙСТВАХ ПРЕДЛАГАЕМОГО МЕТОДА

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

Метод Крейга использовался в достаточно многочисленных расчетах, проведенных для ряда различных плохо обусловленных и некорректных задач больших размеров (см., например, [5]–[8], [10]). При этом успешное применение этого метода было обусловлено, в частности, такими его свойствами.

Пусть, например, $A$ – плохо обусловленная матрица и такая, что при разложении искомого решения СЛАУ (2.1) по собственным векторам матрицы $AA{\kern 1pt} {\text{*}}$ основной компонентой являются те слагаемые, которые соответствуют верхней части спектра этой матрицы. Случай, когда вектор-решение задачи (2.1) “почти” лежит в подпространстве сравнительно небольшой размерности, порожденном старшими собственными векторами, является типичным: такая ситуация возникает, например, если рассматриваемая СЛАУ является аппроксимацией уравнения Фредгольма I рода. Для таких задач при использовании метода Крейга, как показано в [4], достаточная точность достигается за сравнительно небольшое число итераций, существенно меньшее, чем размерность задачи.

Преимущество предлагаемого метода перед указанными выше способами решения рассматриваемой задачи исключения проявляется также в том случае, если решается возмущенная задача, т.е. исходные данные известны с некоторой достаточно малой погрешностью. В этом случае, как сказано в [9], целесообразно согласовывать число фактически проводимых итераций с точностью исходных данных. Следуя этой идее, в работе [10] предложены эффективные критерии прекращения итераций по Крейгу, учитывающие величину возмущения и играющие роль регуляризатора. Как показано в этой работе и подтверждено проведенными расчетами, эти критерии за относительно малое число итераций, существенно меньшее размеров задачи, позволяют получить приближение, достаточно близкое к точному решению исходной невозмущенной задачи, в то время как решение непосредственно рассматриваемой возмущенной задачи существенно от него отличается.

Приведем формулы метода Крейга и критериев прекращения итераций для возмущенной задачи (2.1).

Итак, вместо (2.1) решается задача

где тильдой обозначены возмущенные входные величины. Формулы модифицированного метода Крейга (отметим, что в работе [10], как более ранней, по сравнению с [7], [8], этот метод использован в немодифицированном виде) имеют вид: для $k \geqslant \,1$ имеем

Один из двух предложенных в [10] критериев останова итераций имеет следующий вид. Пусть нам известна оценка

где ${{u}_{*}}$ – искомое минимальное по норме решение исходной задачи (2.1). Тогда итерации метода Крейга выполняются до тех пор, пока

Вариант этого критерия, не требующий знания оценки точного решения задачи, имеет следующий вид. Пусть известно, что $\left| {A{\kern 1pt} {\text{*}} - \tilde {A}{\kern 1pt} {\text{*}}} \right| \leqslant Q$, $\left| {f - \tilde {f}} \right| \leqslant R$. Введем величину ${{\delta }_{k}} = Q\left| {{{u}_{k}}} \right| + R$ и число $D \geqslant 2$. Тогда условие продолжения процесса имеет вид

В работе [10] доказано, что описанный итерационный процесс завершается за конечное число шагов и полученное приближенное решение при $\delta \to 0$ в критерии (4.1) или при $Q + R \to 0$ в критерии (4.2) сходится к точному решению невозмущенной задачи.

5. ИЛЛЮСТРАТИВНЫЙ ПРИМЕР ИСПОЛЬЗОВАНИЯ МЕТОДА

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

Рассмотрим систему уравнений вида (1.1):

Это – переопределенная система, которая решения не имеет. Система (1.2), полученная для этого случая, имеет вид Ранг матрицы этой системы, т.е. $\operatorname{rank} A{\kern 1pt} {\text{*}}A$ равен единице, так что эта система имеет однопараметрическое семейство решений вида $x = {{({3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2} - 2C,\;C)}^{{\text{т}}}}$, где $C$ – произвольная постоянная. Пространство $\ker A = \ker A{\kern 1pt} {\text{*}}A$ составляют векторы $z = {{({{z}_{1}},{{z}_{2}})}^{{\text{т}}}}$, удовлетворяющие уравнению ${{z}_{1}} + 2{{z}_{2}} = 0$, т.е. векторы вида $z = C{{(2, - 1)}^{{\text{т}}}}$. Векторы $f\, \bot \,\ker A$ удовлетворяют уравнению $2{{f}_{1}} - {{f}_{2}} = 0$. Возьмем, например, $f = (1,2)$. Для такого вектора $f$ СЛАУ (2.1) эквивалентна уравнению В работе было сказано, что метод Крейга для этой системы сходится к минимальному по норме решению, и, следовательно, ортогональному решениям соответствующей однородной системы. Именно для этого решения искомое значение $\sigma $, как показано выше, вычисляется по формуле (1.8). Для рассматриваемой модельной задачи найдем это решение аналитически, используя условие ортогональности $u\, \bot \,\ker A{\kern 1pt} {\text{*}}$ и условие $z \in \ker A{\kern 1pt} {\text{*}}$: $A{\kern 1pt} {\text{*}}z = {{z}_{1}} + {{z}_{2}} + 2{{z}_{3}} = 0$. Нетрудно видеть, что $\ker A{\kern 1pt} {\text{*}}$ является линейной оболочкой векторов ${{(1, - 1,0)}^{{\text{т}}}}$ и ${{(2,0, - 1)}^{{\text{т}}}}.$ Тогда для этого примера уравнение (2.1) и условия ортогональности составляют следующую систему: Из нее получим интересующее нас единственное решение, которое имеет вид $u = \frac{1}{6}{{\left( {1,1,2} \right)}^{{\text{т}}}}.$ Для него по формуле (1.8) вычислим значение $\sigma = (b,u) = \frac{3}{2}$.

Отметим, что по формуле $\sigma = (x,f) = \left( {3{\text{/}}2 - 2C} \right) + 2C = 3{\text{/}}2$ получается то же самое значение независимо от взятого решения системы (1.2), поскольку $f$ выбрано соответствующим образом. Это является иллюстрацией приведенного выше утверждения, что если задача (1.2) имеет неединственное решение, то независимость значения $\sigma $ от выбора этого решения эквивалентна разрешимости уравнения (2.1).

Таким образом, на простом примере показаны основные теоретические положения предложенного в работе метода решения задачи исключения для СЛАУ, полученной после 1-го преобразования Гаусса.