Tags:    c#

Pivoting

I visse tilfælde vil vores A matrix være i en form som vores løsningsmetoder kan få problemer med at give en korrekt løsning på - man siger at de er ustabile. Det kan for eksempel opstå hvis element A[0,0] (det allerførste element) er nul. Systemet kan også blive ustabilt af andre årsager.

Så vores første opgave er at arrangere lidt om på A så systemet ikke bliver ustabilt. Denne "flytten-rundt" kaldes på engelsk for pivoting. Det betyder "at dreje sig", eller "at dreje om en tap", og det er altså klart hvorfor at man kaldet denne proces for netop pivoting da vi drejer lidt på matricen.

Når vi laver denne flytten-rundt er målet at befolke diagonalen i matricen med nogle "gode" elementer. Gode elementer er nogle som har en relativt set høj værdi. Så når vi laver "pivoting" går vi vores diagonalelementer igennem og bytter det ud med den række som har den største værdi (man kan også bringe kolonnerne ind i betragtningerne udover rækkerne, men det koster mere rent beregningsmæssigt uden at give en nævneværdig fordel). Vi vil kun bytte en række ud med en række der ligger under det pågældende diagonalelements egen række (ellers omgører vi bare vores tidligere ombytninger i pivoting-processen). Endvidere, for at undgår at en eventuelt skalering påvirker pivot-resultatet vil vi under pivoting betragte hvert element som et tal mellem 0 og 1 ved at skalere det i forhold til rækkens numerisk største element.

I vores skelet for klassen EquationSolver havde vi ikke noget om pivoting - vi vil derfor tilføje en privat metode, kaldet Pivoting, som tager sig at denne rykken-rundt på A. Den er her:

// The array swaps must be properly initialized before calling
// this method. That will be done in the constructor.
private void Pivot()
{
    // Number of rows in the A matrix
    int numberOfRows = a.Length;
    // scaleFactors will contain the needed scaling factors of
    // each row in order to get the row's elements in the
    // 0-1 range.
    double[] scaleFactors = new double[numberOfRows];

    //
    // Detemine the scaling factors
    //
    
    // Number of columns in the A matrix
    int numberOfColumns = a[0].Length;
    for (int i = 0; i < numberOfRows; ++i)
    {
        // We also fill up the swap array with
        // values indicating a clean slate of swap history
        swaps[i][0] = i; 
        swaps[i][1] = i; 
        for (int j = 0; j < numberOfColumns; ++j)
        {
            // We compare the earlier saved scaling factor
            // with the current value from matrix A (a[i][j]),
            // and choose the one with the largest numerical
            // value
            scaleFactors[i] = Math.Max(scaleFactors[i],
                Math.Abs(a[i][j])); 
        }
    }

    //
    // Determine the pivot element
    //

    // Will conatin the row number wiht the pivot element
    int max = 0;

    // Will hold a row when swapping two rows
    double[] tempRow;

    // Will hold an element from the b vector
    double tempBValue;

    // Will hold a temporary value form the scaling array
    double tempScalingValue;
    
    // We start with the left most column
    // and work our way to the next to last
    // column (the last column will then implicitly
    // be set)
    for (int j = 0; j < numberOfColumns - 1; ++j)
    {
        // So far pivot element is in the diagonal's row
        max = j;
        // We only consider rows beneath the current diagnoal element's
        // row, i.e. our starting row number is equal to the current column
        // plus one
        for (int i = j + 1; i < numberOfRows; ++i)
        {
            // We compare the current pivot element, with the
            // element in the current row.
            if (Math.Abs(a[i][j]) / scaleFactors[i] >
                Math.Abs(a[max][j]) / scaleFactors[max])
            {
                // The element in the current row is larger,
                // so save the row number
                max = i;
            }
        }

        //
        // Now we have found the pivot element for the current
        // column, and will do the
        // actual rearranging, if necessary
        //

        // We will only do any rearraing, if we have actually found
        // a pivot element
        if (max != j)
        {
            // Save the row swapping history
            swaps[j][0] = j;
            swaps[j][1] = max;

            // Switch the rows
            tempRow = a[j];
            a[j] = a[max];
            a[max] = tempRow;

            // Don't forget the linear algebra rules!
            // We have manipulated A, and must therefore
            // do a similar manipulation of b
            tempBValue = b[j];
            b[j] = b[max];
            b[max] = tempBValue;

            // Remember to also switch values in the
            // scaling array
            tempScalingValue = scaleFactors[j];
            scaleFactors[j] = scaleFactors[max];
            scaleFactors[max] = tempScalingValue;
        }
    }
}

I vores Pivot metode gør vi brug af vores objekts private felter: matricen a, vektoren b, samt en ny tabel - kaldet swaps - som bruges til at gemme en historie over rækkebytninger som kræves i Gauss-Jordan metoden. Denne tabel er erklæret i vores klasse som de andre felter:

private int[][] swaps; // Needed by Gauss-Jordan

Denne tabel vil vi sørge for bliver initialiseret i klassens konstruktør, som derfor ser sådan ud nu:

public EquationSolver(double[][] a, double[] b, SolveTechnique theSolveTechnique)
{
    this.theSolveTechnique = theSolveTechnique;
    this.a = a;
    this.b = b;            
    swaps = new int[a.Length][];
    for (int i = 0; i < a.Length; ++i)
    {
        swaps[i] = new int[2];
    }            
}

Gauss-Jordan metdoen

Nu er vi langt om længe nået til det som det hele handler om: løsningen af ligningssystemet - men alt vores tidligere arbejde har været strengt nødvendig for overhovedet at være i stand til at finde en løsning.

I denne sektion vil vi kigge på, og implementere, Gauss-Jordan metoden.

Gauss-Jordan metoden (fremefter kaldt GJ) er en ganske god algoritme til at løse et ligningssytem; dens ulempe ligger i at den sammenlignet med de fleste andre metoder, er ganske dyr beregningsmæssigt.

Som nævnt ovenfor er vores overordnede strategi at omdanne vores A matrix til en form som er nem at arbejde med, og i GJ omdanner vi A til enhedsmatricen. En enhedsmatrix er en kvadratisk matrix med værdien 1 i diagonalen og 0 alle andre steder. F.eks. som i denne 4x4 matrix:

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

Vi får omdannet A til enhedsmatricen ved at subtrahere multipla af én række fra en anden. Når A er omdannet til enhedsmatricen, har vi faktisk løsningen:

Ix = A^-1b

Her er I enhedsmatricen, x vores ubekendte, A^1 den inverse af A og b er vores koefficientvektor som nu indeholder løsningen. Som det ses kan GJ altså ikke blot benyttes til at løse ligningssystemet men også til at finde As inverse matrix.

I implementationen vil As inverse blive konstrueret kolonne for kolonne mens A bliver hivet fra hindanden kolonne for kolonne, så vi skal ikke brug mere plads end blot vores to-dimensionalle a variabel. Derudover er det meget vigtigt at bruge pivoting ved GJ. Denne pivoting vil dog forkludre As inverse A^1, forstået på den måde at dets kolonner vil være i den forkerte rækkefølge. Det var derfor vi gemte rækkebytningen i Pivot-metoden i swaps tabelvariablen; med den i hånden kan vi nemlig få den rigtig A^1 frem.

Så i oversigtsform ser vores GJ metode ud som følger:

1. Pivot vores a variabel (og b vektoren)


2. Reducér den "drejede" a variabel til enhedsmatricen


3. Når vi manipulerer a skal vi også huske at manipulerer b


4. Til sidst byt rund på rækkerne i a variablen så vi får den rigtige A^1. Dette er kun nødvendig hvis vi også ønsker at finde A^1

Koden for GJ er som følger:

private void SolveByGaussJordan()
{
    // Start by do a pivoting of a and b
    Pivot();

    // Reduce A row by row by dividing by the diagonal element
    double diagonalElement;
    double tempElement;
    for (int row = 0; row < a.Length; ++row)
    {
        diagonalElement = a[row][row];
        for (int col = 0; col < a[0].Length; ++col)
        {
            a[row][col] /= diagonalElement;
        }
        // We must not forget to manipulate b
        b[row] /= diagonalElement;
        a[row][row] = 1.0 / diagonalElement;

        // Subtract a multiple of the current rows from the
        // other rows
        for (int newRow = 0; newRow < a.Length; ++newRow)
        {
            if (newRow != row)
            {
                tempElement = a[newRow][row];
                for (int newCol = 0; newCol < a[0].Length; ++newCol)
                {
                    a[newRow][newCol] -= tempElement * a[row][newCol];
                }
                b[newRow] -= tempElement * b[row];
                a[newRow][row] = -tempElement * a[row][row];
            }
        }
    }

    // Convert A to A^1 - i.e. undo the wrong columns done by
    // pivoting
    for (int col = a[0].Length - 1; col >= 0; --col)
    {
        // Did a swap occur?
        if (swaps[col][0] != swaps[col][1])
        {
            for (int row = 0; row < a.Length; ++row)
            {
                tempElement = a[row][swaps[col][1]];
                a[row][swaps[col][1]] = a[row][swaps[col][0]];
                a[row][swaps[col][0]] = tempElement;
            }
        }
    }
}

Gauss-elimination

I Gauss-Jordan metoden fandt vi ligningssystemets løsning ved at reducere A matricen til enhedsmatricen - hvilket betød at vi derved også fik givet As inverse matrix.

Hvis vi kun er interesseret i at løse ligningssystemet og er ligeglad med As inverse matrix, så kan vi faktisk nøjes med at reducere A til en "halv enhedsmatrix". En "halv enhedsmatrix" er selvfølgelig noget vrøvl, men betegnelsen dækker over at den nedre del (eller øvre del) af matricen er elementer med værdien 0. Følgende er et eksempel på hvad jeg mener (her har vi kun reduceret den nedre halvdel):

| x00 x01 x02 x03 |
|   0 x11 x12 x13 |
|   0   0 x22 x23 |
|   0   0   0 x33 |

I ovenstående matrix angiver tallene blot elementets indeks, f.eks. betyder x22 "elementet x i matricen X på 3. række og 3. kolonne" (husk vi tæller fra nul).

I Gauss elimination er strategien netop at reducere A "halvt" så at sige - det kan ses af følgende formel hvor matricen b' angiver en modificeret b-vektor:

| x00 x01 x02 x03 | | x0 |   |b0'|
|   0 x11 x12 x13 | | x1 |   |b1'|
|   0   0 x22 x23 | | x2 | = |b2'|
|   0   0   0 x33 | | x3 |   |b3'|

Når vi har fået reduceret A til en halv diagonal matrix, og har b' som i ovenstående formel, finder vi vores løsningsvektor ved at starte med den nederste række. I ovenstående formel kan vi finde løsningen for x3 direkte:

x3 = b3' / x33.

Med løsningen af x3. Kan vi nu løse x2 direkte da vi kender x3; følgende formel isoleres blot for x2:

x2 * x22 + x3 * x23 = b2'

På samme møder finder vi x1 og x0.

Gauss-elimination er en hurtigere metode end Gauss-Jordan metoden, og følgende kode viser en implementation:

private void SolveByGauss()
{
    // Start by do a pivoting of a and b
    Pivot();

    // We will turn a into an upper diagonal matrix. This will
    // be done by subtracting a multipla of the current row from
    // the rows beneath it - we also need to remember to manipulate
    // the b vector accordingly
    for (int i = 0; i < a.Length; ++i)
    {
        b[i] /= a[i][i];

        for (int j = a[0].Length - 1; j >= i; --j)
        {
            a[i][j] /= a[i][i];
        }

        for (int k = i + 1; k < a.Length; ++k)
        {
            b[k] -= a[k][i] * b[i];
            
            for (int m = i + 1; m < a[0].Length; ++m)
            {
                a[k][m] -= a[k][i] * a[i][m];
            }
        }
    }

    // Now we have created out upper diagonal A and likewise
    // manipulated b. So all that's left to do is some back
    // substitution as I explained in the article on
    // udvikleren.dk
    for(int n = a.Length - 2; n >= 0; --n)
    {
        for (int p = n + 1; p < a[0].Length; ++p)
        {
            b[n] -= a[n][p] * b[p];
        }
    }
}

Hvad synes du om denne artikel? Giv din mening til kende ved at stemme via pilene til venstre og/eller lægge en kommentar herunder.

Del også gerne artiklen med dine Facebook venner:  

Karma barometer (30 dage)

IT jobs | Kontakt | Om Udvikleren.dk | Om karma | IT nyheds arkiv | Annoncering

X