АЛГОРИТМЫ

Новости

Рассылка новостей

Форум

AlgoPascal

Редактор блок-схем

Статьи

О сайте

Контакты
Содержание - Системы линейных уравнений - Метод ортогонализации с минимизацией невязки

Решение СЛАУ методом ортогонализации с минимизацией невязки

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

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

Процесс ортогонализации позволяет по системе векторов построить эквивалентную им систему ортогональных векторов. На входе процесс принимает систему векторов A = (a1 , ..., am ) и строит систему векторов B = (b1 , ..., bm ) по следующему правилу:

По завершении процесса полученные вектора нормируются. Система векторов B обладает следующими свойствами:

  1. система B эквивалентна системе A
  2. если bi  = 0 , то ai  линейно выражается через вектора a1 , ..., ai-1 
  3. система векторов B ортонормированная

Если бы матрица системы A состояла из ортонормированных или нулевых столбцов, то i-ая компонента искомого решения с минимальной невязкой получалась бы скалярным умножением i-ого столбца матрицы на столбец свободных членов. При помощи данного алгоритма мы можем привести A к нужному виду, но поскольку операции проводится над столбцами, то мы фактически переходим к новому базису, в котором и найдем решение. Для того же, чтобы получить решение в старом базисе, нам надо осуществить переход к нему. Для этого проведем над единичной матрицей те же линейные операции, которые мы проводили в ходе процесса ортогонализации над матрицей A. Полученная матрица H и есть матрица перехода. Если G - ортогонализированная матрица A, то решение записывается, как:

Данный алгоритм обладает большей функциональностью, чем методы Гаусса, QR- или LU-разложения, но его трудоемкость в четыре-пять раз выше. Также он менее устойчив к накоплению погрешностей при расчетах - системы с "плохими" матрицами решаются им заметно хуже других методов. Для решения классической задачи лучше выбрать другой метод, но в случае, если наличие или единственность решения не гарантированы, он может являться альтернативой ещё более медленному, хотя и более надежному и функциональному SVD-разложению.

Проблема точности

Есть одна тонкость, которую следует обсудить. В формуле ортогонализации есть дробь, в знаменателе которой находится скалярный квадрат вектора bi . Если система векторов A линейно-зависима, то в одной из таких дробей в знаменателе окажется ноль. Решая задачу аналитически, мы просто отбросим нулевой вектор и продолжим процесс без него. Но при численном решении полученный вектор будет очень близок к нолю, но не равен ему из-за погрешностей при расчетах. Поэтому вместо проверки на равенство нолю надо проводить проверку на близость к нолю, и если вектор ближе заданной величины epsilon, то мы считаем его нулевым.

Тем не менее, при таком подходе возможно как ненахождение существующего решения, так и нахождение несуществующего. Рассмотрим это на примере:

Пусть ортогонализируя матрицу A, мы получили два ортогональных вектора g1  и g2 . Правая часть системы равна b. Вектор g2  намного меньше вектора g1 . Вполне возможно, что он неравен нолю из-за возникшей при расчетах погрешности. Если мы не отбросим его (выбрав слишком малое epsilon), то сможем разложить вектор b по базису G и получим ошибочное решение. Т.е. система уравнений несовместна, а точное решение как будто бы существует. Разумеется, в результате мы получим мусор, а не решение.

Возможна и другая ошибка - мы выбрали слишком большое epsilon и отбросили вектор g2 , хотя этот вектор не равен нолю и не близок к нему. Тогда вместо вектора b мы получим только его проекцию на ось g1 , и решение будет ошибочным - мы заключим, что система несовместна, хотя она на самом деле совместна.

Одним из узких мест в работе метода ортогонализации является именно принятие решения о том, выражается ли один вектор через другой или нет. Именно здесь, если задать слишком малое epsilon, и возникают ошибки, искажающие решение. Если epsilon слишком мало, то при делении на него получаются слишком большие числа, которые просто переполняют разрядную сетку, что приводит к гигантской потере точности. К сожалению, нельзя дать никаких общих рекомендаций относительно выбора величины epsilon, так как её оптимальный выбор зависит от входных данных. Здесь всё зависит только от программиста.

Если нашли ошибку в алгоритме - сообщите!


Реализации алгоритма на различных языках:

Реализация алгоритма на C++
Реализация алгоритма на Delphi
Реализация алгоритма на Visual Basic 6

Реализация алгоритма на AlgoPascal:

unit MinimalMisclosureUnit;

interface MinimalMisclosureSolve;

implementation
(*******************************************************************
Функция решения системы линейных уравнений вида:
A[1,1]*X[1] + .. + A[1,N]*X[N] = A[1,N+1]
        ....................
A[N,1]*X[1] + .. + A[N,N]*X[N] = A[N,N+1]
с минимизацией вектора невязки.

Функция   возвращает   True,   если  система уравнений имеет  точное
решение, и False, если система не имеет точного решения.
Вектор X содержит в себе точное решение (или одно из решений)    или
решение  с  минимальной невязкой при отсутствии точного решения.

Epsilon определяет, насколько столбец должен отличаться от ноля, чтобы
считаться неравным нолю и учавствовать в разложении других столбцов в
качестве базисного вектора.
*******************************************************************)
function MinimalMisclosureSolve(
            A   :   array of array of Real;
    const   N   :   Integer;
    out     X   :   array of Real;
    const   Epsilon : Real):Boolean;
var
    I               :   Integer;
    J               :   Integer;
    K               :   Integer;

    L1              :   Real;

    IsBasisVector   :   array of Boolean;

    E               :   array of array of Real;
    C               :   array of Real;
begin
    SetBounds(IsBasisVector, [1,N+1]);
    SetBounds(X, [1,N]);
    SetBounds(C, [1,N]);


    //установим единичную матрицу, добавим фиктивный N+1-ый столбец
    SetBounds(E, [1,N], [1,N+1]);
    for I:=1 to N do
        for J:=1 to N do
            if I=J then
                E[I,J] := 1
            else
                E[I,J] := 0;
                
    //сохраним последний столбец
    for I:=1 to N do
        C[I]:=A[I,N+1];

    //ортонормируем систему стобцов
    for J:=1 to N+1 do
    begin
        for K:=1 to J-1 do
        begin
            if not IsBasisVector[K] then
                Continue;

            L1 := 0;
            for I:=1 to N do
                L1 := L1+A[I,K]*A[I,J];

            for I:=1 to N do
            begin
                A[I,J] := A[I,J]-A[I,K]*L1;
                E[I,J] := E[I,J]-E[I,K]*L1;
            end;
        end;

        L1:=0;
        for I:=1 to N do
            L1 := L1+A[I,J]*A[I,J];
        L1:=Sqrt(L1);
        IsBasisVector[J] := L1>=Epsilon;
        
        if IsBasisVector[J] then
            for I:=1 to N do
            begin
                A[I,J] := A[I,J]/L1;
                E[I,J] := E[I,J]/L1;
            end;
    end;

    //точное ли решение
    Result := not IsBasisVector[N+1];
    
    //вывод решения
    for I:=1 to N do
        X[I] := 0;
    for J:=1 to N do
    begin
        if not IsBasisVector[J] then
            Continue;

        L1:=0;
        for I:=1 to N do
            L1:=L1 + A[I,J]*C[I];
        for I:=1 to N do
            X[I]:=X[I] + L1*E[I,J];
    end;

end;

end.

 


Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU