Il presente articolo è stato tradotto automaticamente.

C#

Scomposizione della matrice

James McCaffrey

Scarica il codice di esempio

Decomposizione di Matrix, una tecnica che si rompe una matrice quadrata numerica in due differenti matrici quadrate, è la base per risolvere in modo efficiente un sistema di equazioni, che a sua volta è la base per invertire una matrice.E inversione di una matrice è una parte di molti importanti algoritmi.Questo articolo presenta e spiega il codice c# che esegue la decomposizione della matrice, inversione di matrice, un sistema di soluzione di equazioni e operazioni correlate.

Certo, decomposizione della matrice non è un argomento appariscente, ma un insieme di metodi di matrice può essere un'importante aggiunta alla libreria di codice personale.I metodi sono spiegati in modo è possibile modificare il codice sorgente per soddisfare le proprie esigenze.Inoltre, alcune delle tecniche utilizzate nei metodi di matrice possono essere riutilizzati in altri scenari di codificante.

Il modo migliore per ottenere un tatto per il tipo di informazioni presentate in questo articolo è quello di dare un'occhiata nella schermata della Figura 1.Il programma demo inizia con la creazione di una matrice quadrata di 4 x 4 e mostrando i suoi valori.Successivamente, la matrice viene scomposto in quello che ha chiamato una matrice LUP.L sta per basso e U sta per upper.La parte P (P sta per permutazione) è una matrice con i valori {3.120} e indica che le righe 0 e 3 sono state scambiate durante il processo di decomposizione.La decomposizione inoltre generato un valore toggle-1, che indica che si è verificato un numero dispari di scambi di riga.Il programma demo Visualizza la scomposizione in due modi: in primo luogo come una matrice di LU combinata e poi separate matrici L ed U.Successivamente, il programma calcola e visualizza l'inversa della matrice originale, utilizzando la matrice LUP dietro le quinte.Il programma demo calcola il determinante della matrice originale, utilizzando nuovamente la decomposizione.Quindi l'inversa della matrice viene utilizzata per risolvere un sistema di equazioni lineari e conclude combinando le matrici L ed U nuovamente dentro la matrice originale.

Figura 1 Matrix decomposizione Demo

Ma perché andare a tutti i problemi di creazione di un metodo di decomposizione di matrice personalizzata e una libreria di metodi correlati?Anche se sono disponibili molti strumenti di matrice autonomo, può essere a volte difficile da integrare in un'applicazione o sistema.E nonostante l'importanza fondamentale della decomposizione di matrix, ci sono alcuni libere, non protetti .NET codice implementazioni disponibili; quelle che esistono sono spesso non spiegate abbastanza dettagliatamente per modificare il codice sorgente per soddisfare gli scenari di codifica.

Questo articolo si presuppone intermediate C# competenze di programmazione e almeno una conoscenza di base della terminologia e operazioni matriciali.Tutto il codice c# è presentato in questo articolo.Il codice è disponibile anche dal sito MSDN download codice a archve.msdn.microsoft.com/mag201212matrix.

Definizione di matrice

Ci sono diversi modi per implementare una matrice in c#.L'approccio tradizionale e quello utilizzato in questo articolo, è quello di utilizzare una matrice di matrici, a volte chiamato una matrice di matrici.Ad esempio, questo codice definisce una tabella con tre righe e due colonne:

double[][] m = new double[3][];
m[0] = new double[2];
m[1] = new double[2];
m[2] = new double[2];
m[2][1] = 5.0; // set row 2, col 1

A differenza della maggior parte dei linguaggi di programmazione c# ha un tipo di matrice multidimensionale integrato, che fornisce un approccio alternativo. Ad esempio:

double[,] m = new double[3,2];
m[2,1] = 5.0;

Un terzo approccio all'implementazione di matrici in C# è quello di utilizzare una singola matrice combinata con la manipolazione di array indice, come questo:

int rows = 3;
int cols = 2;
double[] m = new double[rows * cols];
int i = 2;
int j = 1;
m[i * cols + j] = 5.0;

Indipendentemente dal regime di archiviazione utilizzato, matrici possono essere implementati utilizzando un approccio di metodo statico o un OOP. Ad esempio, potrebbe somigliare un approccio OOP:

public class MyMatrix
{
  int m; // number rows
  int n; // number columns
  data[][]; // the values
  ...
}

Non non c'è nessun singolo scelta migliore per la progettazione a matrice; il miglior design dipende dalla codifica particolare scenario che sta operando e preferenze personali e codificante. Questo articolo utilizza un approccio di metodo statico perché è il più facile da capire ed effettuare il refactoring.

Quando si utilizza una progettazione di matrici di matrici di matrici, perché ogni riga deve essere allocata separatamente, è spesso conveniente definire un metodo di supporto per eseguire l'allocazione di memoria. Ad esempio:

static double[][] MatrixCreate(int rows, int cols)
{
  // creates a matrix initialized to all 0.0s
  // do error checking here?
double[][] result = new double[rows][];
  for (int i = 0; i < rows; ++i)
    result[i] = new double[cols]; // auto init to 0.0
  return result;
}

Il metodo può essere chiamato in questo modo:

double[][] m = MatrixCreate(3,2);
m[2][1] = 5.0;

Questo metodo viene illustrato uno dei vantaggi di creare la propria libreria di metodi di matrice: Se si desidera migliorare le prestazioni, è possibile omettere i parametri di input di controllo errori — a scapito di aumentare il rischio del codice chiamante causando un'eccezione. (Per mantenere questo breve articolo, maggior controllo degli errori è stato rimosso.) Un altro vantaggio è che è possibile personalizzare la vostra libreria per ottimizzare per il tuo scenario esatto. Lo svantaggio principale di creare la tua libreria è che può richiedere più tempo rispetto all'utilizzo di una libreria esistente.

Un altro metodo comodo per aggiungere alla vostra libreria di matrix è uno che può essere utilizzato per visualizzare una matrice come una stringa. Qui è una possibilità:

static string MatrixAsString(double[][] matrix)
{
  string s = "";
  for (int i = 0; i < matrix.Length; ++i)
  {
    for (int j = 0; j < matrix[i].Length; ++j)
      s += matrix[i][j].ToString("F3").PadLeft(8) + " ";
    s += Environment.NewLine;
  }
  return s;
}

Si consiglia di impostare parametri per il numero di decimali da visualizzare, l'imbottitura di larghezza di colonna o di entrambi.

Moltiplicazione di matrici

Il Microsoft .NET Framework 4 e versioni successive forniscono un modo pulito per migliorare significativamente le prestazioni di un metodo di moltiplicazione di matrice. Moltiplicazione di matrici è illustrata Figura 2.

Moltiplicazione di matrici di figura 2

Si noti che il calcolo di ciascun valore di cella del risultato non dipende qualsiasi altri valori di cella nel risultato, quindi ogni calcolo è indipendente e potenzialmente possono essere eseguite in parallelo su un computer con più processori. Qui è un approccio standard alla moltiplicazione di matrici:

static double[][] MatrixProduct(double[][] matrixA, 
  double[][] matrixB)
{
  int aRows = matrixA.Length; int aCols = matrixA[0].Length;
  int bRows = matrixB.Length; int bCols = matrixB[0].Length;
  if (aCols != bRows)
    throw new Exception("Non-conformable matrices in MatrixProduct");
  double[][] result = MatrixCreate(aRows, bCols);
  for (int i = 0; i < aRows; ++i) // each row of A
    for (int j = 0; j < bCols; ++j) // each col of B
      for (int k = 0; k < aCols; ++k)
        result[i][j] += matrixA[i][k] * matrixB[k][j];
  return result;
}

Perché la moltiplicazione è un'operazione Combinalo, le prestazioni possono essere un problema. Ad esempio, se la matrice A ha dimensioni 300 x 200 e matrice B ha dimensioni 200 x 400, calcolo il prodotto della A e B richiede 300 * 200 * 400 = 24.000.000 moltiplicazioni. TPL Task Parallel Library () nello spazio dei nomi System.Threading.Tasks il .NET Framework 4 e più tardi lo rende facile da una semplice versione parallelizzata di moltiplicazione di matrici di codice. Una possibilità è:

static double[][] MatrixProduct(double[][] matrixA, 
  double[][] matrixB)
{
  // error check, compute aRows, aCols, bCols
  double[][] result = MatrixCreate(aRows, bCols);
  Parallel.For(0, aRows, i =>
    {
      for (int j = 0; j < bCols; ++j)
        for (int k = 0; k < aCols; ++k)
          result[i][j] += matrixA[i][k] * matrixB[k][j];
    }
  );
  return result;
}

Questa versione costolette su calcoli di righe. Dietro le quinte, TPL genera tutto il codice dell'impianto idraulico complesso sincronizzazione per eseguire i calcoli su più processori.

Test di coerenza

Un aspetto interessante di una libreria di metodi che sono legati alla vicenda è che spesso è possibile verificarle controllando per vedere se producono risultati coerenti. Ad esempio, si supponga di avere un metodo che crea una matrice casuale:

static double[][] MatrixRandom(int rows, int cols,
  double minVal, double maxVal, int seed)
{
  // return matrix with values between minVal and maxVal
  Random ran = new Random(seed);
  double[][] result = MatrixCreate(rows, cols);
  for (int i = 0; i < rows; ++i)
    for (int j = 0; j < cols; ++j)
      result[i][j] = (maxVal - minVal) * ran.NextDouble() + minVal;
  return result;
}

Inoltre, si supponga che si dispone di una matrice che crea la matrice identità, che è una matrice quadrata con italiano sulla diagonale principale, 0bps altrove:

static double[][] MatrixIdentity(int n)
{
  double[][] result = MatrixCreate(n, n);
  for (int i = 0; i < n; ++i)
    result[i][i] = 1.0;
  return result;
}

E si supponga di che avere un metodo che confronta due matrici per l'uguaglianza:

static bool MatrixAreEqual(double[][] matrixA,
  double[][] matrixB, double epsilon)
{
  // true if all values in A == corresponding values in B
  int aRows = matrixA.Length;
  int bCols = matrixB[0].Length;
  for (int i = 0; i < aRows; ++i) // each row of A and B
    for (int j = 0; j < aCols; ++j) // each col of A and B
      if (Math.Abs(matrixA[i][j] - matrixB[i][j]) > epsilon)
        return false;
  return true;
 }

Si noti che il metodo MatrixAreEqual non confrontare i valori delle celle per uguaglianza esatta perché i valori sono di tipo double. Invece, il metodo controlla per vedere se tutti i valori delle celle sono molto vicino (all'interno di epsilon) a vicenda.

Perché il prodotto di qualsiasi matrice quadrata m e la matrice di identità della stessa dimensione è uguale alla m matrice originale, è possibile testare il metodo del prodotto matrice lungo le linee di questo codice:

double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] i = MatrixIdentity(4);
double[][] mi = MatrixProduct(m, i);
if (MatrixAreEqual(m, mi, 0.00000001) == true)
  Console.WriteLine("Consistent result");
else
  Console.WriteLine("Something is wrong");
Consistency checking lends itself well to random input testing.

Scomposizione della matrice

Decomposizione della matrice è una matrice quadrata M e calcola due nuove matrici quadrate che quando moltiplicati insieme danno la matrice originale M. L'idea è simile al factoring numero ordinario: il numero 6 può essere scomposto in 2 e 3 perché 2 * 3 = 6. In un primo momento che può sembrare che ci sia poco punto in decomposizione di una matrice, ma scopre che decomposizione di matrice rende il compito molto difficile di inversione di matrice abbastanza un po' più semplice. Ci sono molti generi differenti di decomposizione di matrix, e ogni genere può essere calcolate utilizzando diversi algoritmi diversi. La tecnica presentata in questo articolo è chiamata decomposizione LUP e utilizza il metodo di Doolittle con pivoting parziale.

Per capire la decomposizione LUP, è utile conoscere prima la decomposizione LU più semplice, che è stato introdotto il famoso matematico Alan Turing. Si supponga di che disporre di questa matrice 4 x 4 m:

9.000      5.000      3.000      4.000
4.000      8.000      2.000      5.000
3.000      5.000      7.000      1.000
2.000      6.000      0.000      8.000

Una possibile decomposizione LU di M è L =

1.000      0.000      0.000      0.000
0.444      1.000      0.000      0.000
0.333      0.577      1.000      0.000
0.222      0.846     -0.219      1.000

E U =

9.000      5.000      3.000      4.000
0.000      5.778      0.667      3.222
0.000      0.000      5.615     -2.192
0.000      0.000      0.000      3.904

Questo funziona perché L * U = M. Si noti che la matrice L inferiore ha italiano sulla diagonale e 0bps in alto a destra. In altre parole, i valori significativi delle cellule della matrice inferiore sono in basso a sinistra. Allo stesso modo i valori significativi delle cellule della matrice superiore sono sulla diagonale principale e in alto a destra.

Si noti anche che non non c'è nessuna sovrapposizione delle posizioni dei valori delle celle significative di L e U. Così, invece di generare due matrici di risultati, L e U, decomposizione di matrice solitamente memorizza i risultati inferiori e superiori in una singola matrice che detiene L e U per risparmiare spazio di memoria.

LUP matrix decomposizione è una lieve ma importante variazione sulla decomposizione LU. Decomposizione LUP accetta una matrice M e produce matrici L ed U, ma anche una matrice P. Il prodotto di L e U in LUP non è esattamente la matrice originale M ma invece è una versione di M dove alcune delle righe sono state riorganizzate. Il modo in cui le righe sono state riorganizzate viene memorizzato nella matrice P; Questa informazione può essere utilizzata per ricostruire la matrice originale M.

Un cugino vicino alla decomposizione Doolittle presentato in questo articolo è chiamato decomposizione di Crout. La differenza principale tra Doolittle e Crout è che Doolittle luoghi italiano sulla diagonale principale della matrice L e Crout luoghi italiano sulla diagonale principale della matrice U.

Il motivo di decomposizione LUP è usato più spesso di decomposizione LU è sottile. Come vedrà a breve, decomposizione della matrice viene utilizzato per calcolare l'inversa di una matrice. Tuttavia, quando la decomposizione della matrice è utilizzato come supporto per inversione di matrice, si scopre che l'inversione non riuscirà se c'è un valore di 0.0 sulla diagonale principale della matrice LU. Così in decomposizione LUP, quando un 0.0 valore finisce sulla diagonale principale, l'algoritmo scambia due righe per spostare il valore 0.0 fuori diagonale e tiene traccia di quali righe sono state scambiate nella matrice P.

Figura 3 elenca un metodo di decomposizione della matrice.

Figura 3 Metodo di decomposizione Matrix

static double[][] MatrixDecompose(double[][] matrix,
  out int[] perm, out int toggle)
{
  // Doolittle LUP decomposition.
// assumes matrix is square.
int n = matrix.Length; // convenience
  double[][] result = MatrixDuplicate(matrix);
  perm = new int[n];
  for (int i = 0; i < n; ++i) { perm[i] = i; }
  toggle = 1;
  for (int j = 0; j < n - 1; ++j) // each column
  {
    double colMax = Math.Abs(result[j][j]); // largest val in col j
    int pRow = j;
    for (int i = j + 1; i < n; ++i)
    {
      if (result[i][j] > colMax)
      {
        colMax = result[i][j];
        pRow = i;
      }
    }
    if (pRow != j) // swap rows
    {
      double[] rowPtr = result[pRow];
      result[pRow] = result[j];
      result[j] = rowPtr;
      int tmp = perm[pRow]; // and swap perm info
      perm[pRow] = perm[j];
      perm[j] = tmp;
      toggle = -toggle; // row-swap toggle
    }
    if (Math.Abs(result[j][j]) < 1.0E-20)
      return null; // consider a throw
    for (int i = j + 1; i < n; ++i)
    {
      result[i][j] /= result[j][j];
      for (int k = j + 1; k < n; ++k)
        result[i][k] -= result[i][j] * result[j][k];
    }
  } // main j column loop
  return result;
}

Il metodo potrebbe essere chiamato come questo:

double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
int[] perm;
int toggle;
double luMatrix = MatrixDecompose(m, out perm, out toggle);

Il metodo MatrixDecompose accetta come input una matrice quadrata. Il metodo ha tre valori restituiti. Il ritorno di esplicita è una matrice permutata LU. Il metodo restituisce due valori come parametri out. Uno è una matrice di permutazione che tiene le informazioni su come le righe furono permutate. Il secondo parametro è un valore toggle che è + 1 o -1 a seconda se il numero di scambi di riga è stato anche (+ 1) o dispari (-1). Il valore toggle non viene utilizzato per l'inversione di matrice, ma è necessario se la decomposizione della matrice viene utilizzato per calcolare il determinante di una matrice.

Il metodo MatrixDecompose è abbastanza difficile, ma realisticamente, ci sono solo pochi dettagli, che è necessario capire per modificare il codice. La versione qui presentata alloca memoria nuova per la matrice di ritorno LU utilizzando un metodo di supporto MatrixDuplicate:

static double[][] MatrixDuplicate(double[][] matrix)
{
  // assumes matrix is not null.
double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
  for (int i = 0; i < matrix.Length; ++i) // copy the values
    for (int j = 0; j < matrix[i].Length; ++j)
      result[i][j] = matrix[i][j];
  return result;
}

Un approccio alternativo è quello di calcolare il risultato nella matrice di input al fine di risparmiare memoria. Con C# semantica, questo renderebbe il parametro matrice un parametro ref perché è usato per l'input e l'output. Utilizzando questo approccio, la firma del metodo sarebbe:

static void MatrixDecompose(ref double[][] matrix, out int[] perm,
  out int toggle)

O, perché restituisca l'esplicito valore è stato eliminato, è possibile utilizzarlo per la matrice di permutazione o alternare gli scambi. Ad esempio:

static int[] MatrixDecompose(ref double[][] matrix, out int toggle)

Si potrebbe voler eliminare il parametro toggle per semplificare la firma del metodo se non si intende utilizzare la decomposizione della matrice per calcolare un determinante.

Un'altra area di MatrixDecompose, si potrebbe voler modificare è questa affermazione:

if (Math.Abs(result[j][j]) < 1.0E-20)
  return null;

In parole, significa che questo codice: "Se, anche dopo lo scambio di due righe perché c'era un valore di 0.0 sulla diagonale principale, c'è ancora un 0.0 lì, quindi restituire null." Si potrebbe voler modificare il valore di epsilon arbitrario per l'uguaglianza di zero controllo da 1.0 e-20 a qualche altro valore basato sulle caratteristiche di precisione del sistema. E anziché restituire null, è possibile generare un'eccezione; Se il metodo doveva essere chiamati come parte dell'inversione di matrice, l'inversione non riuscirebbe. Infine, se si sta utilizzando la decomposizione della matrice per qualche scopo diverso da inversione di matrice, puoi eliminare del tutto questa affermazione.

Inversione di matrice

La chiave per utilizzare la decomposizione della matrice per invertire una matrice è quello di scrivere un metodo di supporto che risolve un sistema di equazioni. Questo metodo di supporto chiave è presentato Figura 4.

Figura 4 il metodo di HelperSolve

static double[] HelperSolve(double[][] luMatrix, 
  double[] b)
{
  // solve luMatrix * x = b
  int n = luMatrix.Length;
  double[] x = new double[n];
  b.CopyTo(x, 0);
  for (int i = 1; i < n; ++i)
  {
    double sum = x[i];
    for (int j = 0; j < i; ++j)
      sum -= luMatrix[i][j] * x[j];
    x[i] = sum;
  }
  x[n - 1] /= luMatrix[n - 1][n - 1];
  for (int i = n - 2; i >= 0; --i)
  {
    double sum = x[i];
    for (int j = i + 1; j < n; ++j)
      sum -= luMatrix[i][j] * x[j];
    x[i] = sum / luMatrix[i][i];
  }
  return x;
}

Il metodo HelperSolve trova una matrice x che quando moltiplicato per una matrice di LU dà matrice b. Il metodo è abbastanza intelligente, e si può solo pienamente capire dall'analisi tramite alcuni esempi. Ci sono due cicli. Il primo ciclo viene utilizzato sostituzione avanti sulla parte inferiore della matrice LU. Nel secondo ciclo utilizza la sostituzione all'indietro sulla parte superiore della matrice LU. Alcune implementazioni di decomposizione di diversa matrice chiamano loro metodo analogo qualcosa come luBackSub.

Anche se il codice è breve, è difficile, ma non ci dovrebbe essere alcun scenari dove è necessario modificare il codice. Avviso che Helper­Solve accetta una matrice di LU da MatrixDecompose ma non utilizza il parametro out perm. Questo implica che HelperSolve è in realtà un codice di wrapper aggiuntivo Metodo e alle esigenze del supporto per risolvere un sistema di equazioni. Se si il refactoring del codice in questo articolo per un design OOP, si potrebbe voler fare HelperSolve un metodo privato.

Con il metodo HelperSolve nel luogo, il metodo di inversione di matrice può essere implementato, come mostrato Figura 5.

Figura 5 il metodo MatrixInverse

static double[][] MatrixInverse(double[][] matrix)
{
  int n = matrix.Length;
  double[][] result = MatrixDuplicate(matrix);
  int[] perm;
  int toggle;
  double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
  if (lum == null)
    throw new Exception("Unable to compute inverse");
  double[] b = new double[n];
  for (int i = 0; i < n; ++i)
  {
    for (int j = 0; j < n; ++j)
    {
      if (i == perm[j])
        b[j] = 1.0;
      else
        b[j] = 0.0;
    }
    double[] x = HelperSolve(lum, b);
    for (int j = 0; j < n; ++j)
      result[j][i] = x[j];
  }
  return result;
}

Ancora una volta, il codice è ingannevole. L'algoritmo di inversione è basato sull'idea che il prodotto di una matrice M e la sua inversa è la matrice identità. Metodo MatrixInverse risolve essenzialmente un sistema di equazioni Ax = b, dove A è una decomposizione di LU matrice e le costanti b sono 1.0 o 0.0 e corrispondono alla matrice identità. Si noti che MatrixInverse utilizza la matrice perm dalla chiamata a MatrixDecompose.

Chiamare il MatrixInverse metodo potrebbe apparire come questo:

double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] inverse = MatrixInverse(m);
Console.WriteLine(MatrixAsString(inverse));

Per riassumere, un'operazione importante matrice è inversione di matrice, che è abbastanza difficile. Un approccio è quello di scomporre la matrice originale, scrivere un helper risolvere metodo che esegue la sostituzione di avanti e indietro, e poi uso la decomposizione insieme con la matrice di permutazione LU e l'helper risolvere il metodo per trovare l'inverso. Questo approccio può sembrare complicato, ma è solitamente più efficiente e più facile che un inverso della matrice di calcolo direttamente.

Determinante matrice

Con un metodo di decomposizione di matrix in mano, è facile da codificare un metodo per calcolare il determinante di una matrice:

static double MatrixDeterminant(double[][] matrix)
{
  int[] perm;
  int toggle;
  double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
  if (lum == null)
    throw new Exception("Unable to compute MatrixDeterminant");
  double result = toggle;
  for (int i = 0; i < lum.Length; ++i)
    result *= lum[i][i];
  return result;
}

Come si scopre, piuttosto sorprendentemente, il determinante di una matrice è appena più o meno (a seconda del valore toggle) il prodotto dei valori sulla diagonale principale della matrice decomposizione LUP. Si noti che non esiste una conversione implicita dei tipi di valore del cavicchio da int a doppia. Oltre ad aggiungere errore controllo di MatrixDeterminant, si potrebbe voler aggiungere un cortocircuito ritorno in situazioni dove la matrice di input ha dimensioni 1 (poi tornare al valore singolo) o 2 x 2 (quindi restituire un * d - b * c). Chiamare il metodo determinante potrebbe assomigliare a questo:

double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double det = MatrixDeterminant(m);
Console.WriteLine("Determinant = " + det.ToString("F2"));

Risolvere sistemi di equazioni

Il metodo HelperSolve può essere facilmente adattato per risolvere un sistema di equazioni lineari:

static double[] SystemSolve(double[][] A, double[] b)
{
  // Solve Ax = b
  int n = A.Length;
  int[] perm;
  int toggle;
  double[][] luMatrix = MatrixDecompose(A, 
    out perm, out toggle);
  if (luMatrix == null)
    return null; // or throw
  double[] bp = new double[b.Length];
  for (int i = 0; i < n; ++i)
    bp[i] = b[perm[i]];
  double[] x = HelperSolve(luMatrix, bp);
  return x;
}

Ecco il codice che ha prodotto lo screenshot in Figura 1 per risolvere il seguente sistema:

3x0 + 7x1 + 2x2 + 5x3 = 49
 x0 + 8x1 + 4x2 + 2x3 = 30
2x0 +  x1 + 9x2 + 3x3 = 43
5x0 + 4x1 + 7x2 +  x3 = 52

double[][] m = MatrixCreate(4, 4);
m[0][0] = 3.0; m[0][1] = 7.0; m[0][2] = 2.0; m[0][3] = 5.0;
m[1][0] = 1.0; m[1][1] = 8.0; m[1][2] = 4.0; m[1][3] = 2.0;
m[2][0] = 2.0; m[2][1] = 1.0; m[2][2] = 9.0; m[2][3] = 3.0;
m[3][0] = 5.0; m[3][1] = 4.0; m[3][2] = 7.0; m[3][3] = 1.0;
double[] b = new double[] { 49.0, 30.0, 43.0, 52.0 };
double[] x = SystemSolve(m, b);
Console.WriteLine("\nSolution is x = \n" + VectorAsString(x));

Si noti che SystemSolve riorganizza il parametro di input b utilizzando la matrice di perm da MatrixDecompose prima di chiamare HelperSolve.

Capire la matrice di permutazione

Le ultime righe di output nello screenshot in Figura 1 indicano che è possibile moltiplicare le matrici L ed U in modo tale da ottenere la matrice originale. Saper fare questo non vi permetterà di risolvere i problemi pratici di matrix, ma ti aiuterà a capire la parte P di decomposizione LUP. Rigenerare la matrice originale dai suoi componenti L e U può anche essere utile per testare i vostri metodi di libreria matrix per coerenza.

Un modo per rigenerare una matrice originale dopo decomposizione LUP è di moltiplicare L e U e quindi permutare le righe del prodotto utilizzando la matrice P:

// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] lu = MatrixProduct(lower, upper);
double[][] orig = UnPermute(lu, perm);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
  Console.WriteLine("L and U unpermuted using perm array");

Il metodo UnPermute può essere codificato come questo:

static double[][] UnPermute(double[][] luProduct, int[] perm)
{
  double[][] result = MatrixDuplicate(luProduct);
  int[] unperm = new int[perm.Length];
  for (int i = 0; i < perm.Length; ++i)
    unperm[perm[i]] = i; // create un-perm array
  for (int r = 0; r < luProduct.Length; ++r) // each row
    result[r] = luProduct[unperm[r]];
  return result;
}

Un secondo approccio per rigenerare una matrice originale dalla sua decomposizione LUP è quello di convertire la matrice perm in una matrice di perm e quindi moltiplicare la matrice di perm e la matrice di LU combinata:

// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] permMatrix = PermArrayToMatrix(perm);
double[][] orig = MatrixProduct(permMatrix, lu);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
  Console.WriteLine("L and U unpermuted using perm matrix");

Una matrice di perm è una matrice quadrata con uno 1.0 valore in ogni riga e ogni colonna. Il metodo che crea una matrice di perm da una matrice di perm può essere codificato in questo modo:

static double[][] PermArrayToMatrix(int[] perm)
{
  // Doolittle perm array to corresponding matrix
  int n = perm.Length;
  double[][] result = MatrixCreate(n, n);
  for (int i = 0; i < n; ++i)
    result[i][perm[i]] = 1.0;
  return result;
}

Conclusioni

Ci sono molti algoritmi che richiedono la risoluzione di un sistema di equazioni lineari, trovare l'inversa di una matrice o di trovare il determinante di una matrice. Utilizzando la decomposizione della matrice è una tecnica efficace per eseguire tutte queste operazioni. Il codice qui presentato può essere utilizzato in situazioni dove non si desidera dipendenze esterne nella vostra base di codice o hai bisogno la possibilità di personalizzare le operazioni al fine di migliorare le prestazioni o modificare funzionalità. Prendere la pillola rossa!

Dr.James McCaffrey funziona per Volt Information Sciences Inc., dove gestisce la formazione tecnica per gli ingegneri software della Microsoft di Redmond, Washington, campus. Si è occupato di diversi prodotti Microsoft, inclusi Internet Explorer e MSN Search. McCaffrey è l'autore di .NET Test Automation Recipes (Apress, 2006). Può essere raggiunto a jammc@microsoft.com.

Grazie ai seguenti esperti tecnici per la revisione di questo articolo: Paul Koch e Dan Liebling