Juillet 2016

Volume 31, numéro 7

Cet article a fait l'objet d'une traduction automatique.

Série de tests - Inversion de la matrice grâce à C#

Par James McCaffrey

James McCaffreyUne des techniques fondamentales dans les systèmes logiciels (ML) d’apprentissage est une inversion de matrice. Pour des raisons qui ne sont pas clairs pour moi, Microsoft .NET Framework ne semble pas avoir une méthode d’inversion de matrice (ou s’il existe une telle méthode, il est très bien masqué). Dans cet article, je présente et présentent le code pour une méthode d’inversion de matrice qui utilise un algorithme appelé la décomposition de Crout LU.

Permettez-moi de vous être les premiers à reconnaître qu’inversion de matrice n’est pas un sujet très coloré. Mais si vous souhaitez créer des systèmes de ML sans se baser sur les bibliothèques externes, l’utilisation d’une méthode d’inversion de matrice est essentielle car une inversion de matrice est utilisée par des dizaines d’algorithmes ML importants.

Une bonne façon de voir où cet article est menée consiste à examiner le programme de démonstration dans Figure 1.

Démonstration d’Inversion de matrice
Figure 1 matrice Inversion démo

La démonstration commence en configurant et en affichant un 4 x 4 (4 lignes, 4 colonnes) matrice m:

3.0  7.0  2.0  5.0
1.0  8.0  4.0  2.0
2.0  1.0  9.0  3.0
5.0  4.0  7.0  1.0

Ensuite, il calcule l’inverse de la matrice à l’aide d’une méthode définie par le programme et affiche le résultat :

0.097   -0.183   -0.115    0.224
 -0.019    0.146   -0.068    0.010
 -0.087    0.064    0.103   -0.002
  0.204   -0.120    0.123   -0.147

C’est à peu près tout. Mais comme vous le verrez bientôt, inversion de matrice est étonnamment complexe. Ensuite, la démonstration vérifie que l’inverse calculée est correct en multipliant la matrice d’origine heures l’inverse :

1.0  0.0  0.0  0.0
0.0  1.0  0.0  0.0
0.0  0.0  1.0  0.0
0.0  0.0  0.0  1.0

Le résultat de la multiplication est la matrice d’identité (les valeurs 1.0 sur la diagonale, des valeurs 0.0 ailleurs) indiquant le résultat inverse est correct.

Dans les coulisses, la méthode d’inversion de matrice utilise une technique appelée décomposition de matrice. Décomposition intègre une matrice à deux matrices, appelées L (en bas) et U (haut), que lors de la multiplication donner à la matrice d’origine, mais avec certaines lignes réorganisées. La démonstration affiche la décomposition :

5.000    4.000    7.000    1.000
0.200    7.200    2.600    1.800
0.400   -0.083    6.417    2.750
0.600    0.639   -0.602    4.905

La décomposition contienne réellement des matrices le L et U. La matrice L se compose des valeurs dans la partie inférieure gauche de la matrice de %lu combinée, avec les valeurs 1.0 factice sur la diagonale et des valeurs 0.0 dans la partie supérieure :

1.000    0.000    0.000    0.000
0.200    1.000    0.000    0.000
0.400   -0.083    1.000    0.000
0.600    0.639   -0.602    1.000

La matrice U se compose des valeurs dans la partie supérieure droite et la diagonale de la matrice de %lu combinée :

5.000    4.000    7.000    1.000
0.000    7.200    2.600    1.800
0.000    0.000    6.417    2.750
0.000    0.000    0.000    4.905

La démonstration vérifie que la décomposition de %lu est correcte en multipliant les matrices L et U et affiche le résultat :

5.0  4.0  7.0  1.0
1.0  8.0  4.0  2.0
2.0  1.0  9.0  3.0
3.0  7.0  2.0  5.0

Si vous comparez L * U avec la matrice d’origine m, vous verrez que L * U est presque identique m, mais les lignes de L * vous avez été permutés (réorganisé). Les informations de permutation de ligne sont la suivante :

3  1  2  0

Cela signifie que la ligne [0] m est en ligne [3] L * u ; ligne [1] m est en ligne [1] L * u ; ligne [2] m est à la ligne [2] L * u ; et la ligne [3] de m à la ligne [0] * u L.

Présentation d’Inversion de matrice

En arithmétique normale, l’inverse d’un nombre z est un nombre multipliée par z donne 1. Par exemple, si z = 3, l’inverse de z est 1/3 = 0,33 car 3 * (1/3) = 1.

Inversion de matrice étend ce concept. L’inverse d’un nxn (appelé une matrice « carré », car le nombre de lignes est égal au nombre de colonnes) matrice m est un multiple de la matrice par m * mi = I où I est la matrice d’identité (1,0 seconde sur la diagonale, 0.0s ailleurs).

Pas toutes les matrices ont un inverse. En fait, il est une valeur scalaire appelée le déterminant d’une matrice. Si le déterminant d’une matrice est zéro, la matrice n’ont un inverse.

Notez que, pour bien comprendre une inversion de matrice, vous devez comprendre multiplication de matrice. Multiplication de matrice est mieux expliquer par exemple. Examinons l’exemple de Figure 2. La valeur de cellule [r] [c] de la matrice de résultats est le produit des valeurs dans la ligne r de la première matrice et les valeurs dans la colonne c de la deuxième matrice.

Multiplication de matrice
Figure 2 Multiplication des matrices

Lors de la recherche à l’inverse d’une matrice, vous travaillez uniquement avec les matrices carrées, mais la multiplication de matrice peut être appliquée aux matrices de différentes formes. Dans ces situations, les matrices doivent être ce que l'on appelle adaptable. Si la matrice A a forme axn et matrice B a forme nxb, le résultat de la multiplication a forme axb. Le nombre de colonnes dans la première matrice doit égal au nombre de lignes dans la deuxième matrice.

Le programme de démonstration implémente la multiplication de matrice avec une méthode d’assistance et de MatrixProduct MatrixCreate, comme indiqué dans Figure 3. La démonstration utilise une approche en force brute, mais étant donné que le calcul de chaque cellule dans la matrice de résultats est indépendant, multiplication de matrice peut être effectuée en parallèle à l’aide de la méthode Parallel.For à partir de la bibliothèque parallèle de tâches .NET. 

Figure 3 Multiplication des matrices

static double[][] MatrixCreate(int rows, int cols)
{
  double[][] result = new double[rows][];
  for (int i = 0; i < rows; ++i)
    result[i] = new double[cols];
  return result;
}
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");
  double[][] result = MatrixCreate(aRows, bCols);
  for (int i = 0; i < 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;
}

Le programme de démonstration

J’ai codé le programme de démonstration à l’aide de c#, mais vous ne devez avoir aucune difficulté portage du code dans un autre langage, tels que Visual Basic ou Python, si vous le souhaitez. Le code de démonstration est trop long pour être présenté dans son intégralité, mais le code complet est disponible dans le téléchargement qui accompagne cet article. Le code est également disponible à l’adresse quaetrix.com/Matrix/code.html (l’URL respecte la casse).

Pour créer le programme de démonstration, j’ai lancé Visual Studio et créé une nouvelle application de console c# nommée MatrixInverse. Le programme de démonstration n’a aucune dépendance significative de .NET Framework pour n’importe quelle version de Visual Studio fonctionne. Une fois le code du modèle chargé, dans l’Explorateur de solutions fenêtre, j’ai cliqué sur le fichier Program.cs et renommé au MatrixInverseProgram.cs plus descriptif et Visual Studio puis renommer automatiquement la classe programme pour moi.

En haut de la fenêtre de l’éditeur, j’ai supprimé toutes les instructions using qui référence les espaces de noms inutiles, en laissant seulement l’une référence à l’espace de noms système supérieur.

La méthode Main commence en configurant une matrice à inverser :

Console.WriteLine("Begin matrix inverse demo");
double[][] m = MatrixCreate(4, 4);
m[0][0] = 3; m[0][1] = 7; m[0][2] = 2; m[0][3] = 5;
m[1][0] = 1; m[1][1] = 8; m[1][2] = 4; m[1][3] = 2;
m[2][0] = 2; m[2][1] = 1; m[2][2] = 9; m[2][3] = 3;
m[3][0] = 5; m[3][1] = 4; m[3][2] = 7; m[3][3] = 1;

Dans de nombreux cas, vos données sources sont stockées dans un fichier texte, vous devez donc d’écrire une méthode d’assistance pour charger une matrice à partir du fichier. La démonstration utilise une matrice de style tableau de tableaux. Contrairement à la plupart des langages de programmation, c# prend en charge un type de matrice à n dimensions true, mais je préfère utiliser l’approche standard de tableau de tableaux.

Ensuite, principal affiche la matrice m, puis calcule et affiche l’inverse :

Console.WriteLine("Original matrix m is ");
Console.WriteLine(MatrixAsString(m));
double[][] inv = MatrixInverse(m);
Console.WriteLine("Inverse matrix inv is ");
Console.WriteLine(MatrixAsString(inv));

Tout le travail est effectué par la méthode MatrixInverse. Méthode d’assistance MatrixAsString retourne une représentation sous forme de chaîne d’une matrice :

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;
}

Ici, le nombre de décimales (3) et la largeur de la valeur (8) est codées en dur pour plus de simplicité. Une approche plus générale est de passer ces valeurs comme paramètres d’entrée. Ensuite, la démonstration multiplie la matrice d’origine et la matrice inverse afin de vérifier que le résultat est la matrice d’identité :

double[][] prod = MatrixProduct(m, inv);
Console.WriteLine("The product of m * inv is ");
Console.WriteLine(MatrixAsString(prod));

Dans cette version, vous devez vérifier visuellement que le résultat est la matrice d’identité. Une approche plus sophistiquée serait d’écrire une méthode qui accepte une matrice et renvoie la valeur true si la matrice est une matrice d’identité, soumis à quelques petites différences (1.0E-5 est par défaut) dans les valeurs de cellule.

Ensuite, la démonstration illustre certaines des travaux en arrière-plan en décomposant la matrice d’origine :

double[][] lum;
int[] perm;
int toggle = MatrixDecompose(m, out lum, out perm);
Console.WriteLine("The decomposition is");
Console.WriteLine(MatrixAsString(lum));

La signature d’appel de méthode MatrixDecompose peut sembler un peu inhabituelle pour vous. Explicites valeur de retour est + 1 ou -1 en fonction du nombre de permutations ligne comportait (pair ou impair, respectivement). La bascule de retourner la valeur n’est pas utilisée par la démonstration, mais est nécessaire si vous souhaitez calculer le déterminant de la matrice, qui indique si l’inverse d’une matrice existe, comme je l’expliquerai plus loin.

Le lum paramètre out est la décomposition (inférieur à supérieur) de LU combinée. La permission de paramètre de sortie est un tableau de valeurs entières qui codent la façon dont les lignes ont été permutés.

Ensuite, la démonstration extrait les matrices inférieures et supérieures de la matrice de %lu combinée et les affiche :

double[][] lower = ExtractLower(lum);
double[][] upper = ExtractUpper(lum);
Console.WriteLine("The lower part of LUM is");
Console.WriteLine(MatrixAsString(lower));
Console.WriteLine("The upper part of LUM is");
Console.WriteLine(MatrixAsString(upper));

Méthodes d’assistance ExtractLower et ExtractUpper ne sont pas nécessaires pour effectuer une inversion de matrice réellement, mais sont utilisés pour illustrer le fonctionne de la décomposition de matrice.

Le programme de démonstration se termine en affichant les informations de permutation de ligne, en multipliant les matrices de décomposition inférieure et supérieure et affiche le résultat :

Console.WriteLine("The perm[] array is");
ShowVector(perm);
double[][] lowTimesUp = MatrixProduct(lower, upper);
Console.WriteLine("The product of lower * upper is ");
Console.WriteLine(MatrixAsString(lowTimesUp));
Console.WriteLine("End matrix inverse demo");

Méthode d’assistance définie par le programme ShowVector est juste une commodité pour nettoyer la méthode Main. Comme indiqué, le résultat de la partie inférieure * supérieur est la matrice d’origine, sauf que les lignes du résultat sont permutés selon les informations contenues dans le tableau de l’autorisation.

Méthode MatrixInverse

La définition de méthode MatrixInverse commence ainsi :

static double[][] MatrixInverse(double[][] matrix)
{
  int n = matrix.Length;
  double[][] result = MatrixCreate(n, n);
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      result[i][j] = matrix[i][j];
...

La méthode suppose que son paramètre d’entrée, en fait, a une matrice. Cela signifie que vous devez vérifier avant d’appeler, du genre :

int d = MatrixDeterminant(m);
if (d == 0.0)
  Console.WriteLine("No inverse");
else
  double[][] mi = MatrixInverse(m);

Une alternative consiste à placer ce code de vérification des erreurs dans la méthode MatrixInverse, qui crée une copie de la matrice d’entrée. Vous pouvez également effectuer une inversion de matrice en place, ce qui économise de la mémoire mais détruit la matrice d’origine.

Ensuite, MatrixInverse décompose la copie de la matrice d’entrée :

double[][] lum; // Combined lower & upper
int[] perm;
int toggle;
toggle = MatrixDecompose(matrix, out lum, out perm);

Cela peut sembler étrange pour accéder à tous les problèmes de décomposition d’une matrice afin de calculer son inverse, mais croyez-moi, cette approche est beaucoup plus facile qu’inversion directement d’une matrice.

Ensuite, la démonstration calcule l’inverse à l’aide d’une autre méthode d’assistance nommée d’assistance :

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 = Helper(lum, b); //
  for (int j = 0; j < n; ++j)
    result[j][i] = x[j];
}

Ce code est très subtil. Heureusement, vous n’aurez jamais à modifier cette partie du code.

Méthode MatrixInverse se termine en renvoyant l’inverse :

...
  return result;
}

Il doit être clair que méthode MatrixInverse est essentiellement un wrapper autour des méthodes MatrixDecompose et assistance, qui effectuent la plus grande partie du travail. Le code de ces deux méthodes sont dans le téléchargement qui accompagne cet article.                                                    

Le déterminant d’une matrice

Si le déterminant d’une matrice est zéro, la matrice n’a un inverse. Supposons qu’une matrice 3 x 3 est :

1.0  4.0  0.0
3.0  2.0  5.0
7.0  8.0  6.0

Le déterminant de la matrice est la suivante :

+1.0 * [(2.0)(6.0) - (5.0)(8.0)]
-4.0 * [(3.0)(6.0) - (5.0)(7.0)]
+0.0 * [(3.0)(8.0) - (2.0)(7.0)]
= +1.0 * (-28.0) -4.0 * (-17.0) = -28.0 + 68.0 = 40.0

Chaque matrice carrée a un déterminant. Pour les matrices des formes plus de 3 x 3, le calcul déterminant est étonnamment difficile. Toutefois, en fait, si vous décomposez une matrice, vous pouvez utiliser le résultat supérieur-inférieure combiné pour calculer le déterminant facilement en multipliant les éléments du résultat en diagonale. Le programme de démonstration définit une méthode MatrixDeterminant en tant que :

static double MatrixDeterminant(double[][] matrix)
{
  double[][] lum;
  int[] perm;
  int toggle = MatrixDecompose(matrix, out lum,
    out perm);
  double result = toggle;
  for (int i = 0; i < lum.Length; ++i)
    result *= lum[i][i];
  return result;
}

Synthèse

La clé à une inversion de matrice efficace est la décomposition de matrice. Il existe plusieurs algorithmes qui décomposent une matrice. Le code de démonstration utilise une technique appelée algorithme de Crout. Une alternative courante est l’algorithme de Doolittle. J’ai utilisé à préférer les algorithme de Doolittle car il est un peu plus simple que de Crout, mais désormais favoriser algorithme de Crout, car il possède moins de chiffres échec.

J’ai toujours demandé pourquoi le .NET Framework n’a pas une méthode qui calcule l’inverse d’une matrice. Pour être sûr, inversion de matrice est très, très dangereux. J’ai testé le code présenté dans cet article en générant au hasard des matrices carrées de 100 millions de formes entre 2 x 2 et 50 x 50, calcul inverses et par programme en vérifiant que le produit de l’inverse et la matrice d’origine a été la matrice d’identité de la taille appropriée.


Récupération d'urgence. James McCaffrey travaille pour Microsoft Research à Redmond, Wash.  Il a travaillé sur plusieurs produits Microsoft, notamment Internet Explorer et Bing. Dr. Vous pouvez contacter James McCaffrey à jammc@microsoft.com.

Je remercie les experts techniques Microsoft suivants qui cet article : Chris Lee et Kirk Olynyk