Разложение матрицы

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

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

Лучший способ прочувствовать то, о чем пойдет речь в этой статье, — взглянуть на экранный снимок на рис. 1. Демонстрационная программа начинает с создания квадратной матрицы 4 × 4 и отображения ее значений. Затем матрица раскладывается в так называемую матрицу LUP (lower, upper, permutation) (нижняя часть, верхняя часть и часть, относящаяся к перестановке). Последняя часть представляет собой массив со значениями {3,1,2,0} и указывает, что строки 0 и 3 поменялись местами в процессе разложения. В этом процессе также было сгенерировано значение-переключатель (toggle value), равное –1, сообщающее, что выполнено нечетное количество перестановок строк. Эта программа демонстрирует разложение двумя способами: сначала в комбинированную матрицу LU, а затем в отдельные матрицы L и U. Далее программа вычисляет и отображает обращенную по отношению к исходной матрицу, используя «за кулисами» матрицу LUP. Демонстрационная программа вычисляет определитель исходной матрицы, вновь применяя разложение. Далее она использует обратную матрицу для решения системы линейных уравнений и завершает свою работу объединением матриц L и U в исходную матрицу.

Рис. 1. Демонстрация разложения матрицы

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

В этой статье предполагается, что вы владеете программированием на C# хотя бы на среднем уровне и по крайней мере на базовом уровне понимаете матричные операции и сопутствующую терминологию. Все ключевые части кода на C# приведены в самой статье. Полный исходный код можно скачать по ссылке archive.msdn.microsoft.com/mag201207matrix.

Определение матрицы

Существует несколько способов реализации матрицы на C#. Традиционный подход (он используется в этой статье) — применение массива массивов, который иногда называют массивом с неровными границами (jagged array). Например, следующий код определяет матрицу с тремя строками и двумя столбцами:

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; // присваиваем значение строке 2, столбцу 1

В отличие от большинства языков программирования в C# есть встроенный тип многомерного массива, который предоставляет альтернативный подход к матрицам, например:

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

Третий подход к реализации матриц на C# — использование одного массива в сочетании с манипуляциями над индексами массива:

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

Независимо от примененной схемы хранения матрицы можно реализовать, используя либо ООП (объектно-ориентированное программирование), либо подход со статическими методами. Так, подход с ООП мог бы быть примерно таким:

public class MyMatrix
{
  int m; // количество строк
  int n; // количество столбцов
  data[][]; // значения
  ...
}

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

При использовании для матриц массива массивов (поскольку каждая строка должна создаваться отдельно) зачастую удобно определять вспомогательный метод, выполняющий выделение памяти. Например:

static double[][] MatrixCreate(int rows, int cols)
{
  // Создаем матрицу, полностью инициализированную
  // значениями 0.0. Проверка входных параметров опущена.
  double[][] result = new double[rows][];
  for (int i = 0; i < rows; ++i)
    result[i] = new double[cols]; // автоинициализация в 0.0
  return result;
}

Этот метод можно вызывать так:

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

Данный метод демонстрирует одно из преимуществ создания собственной библиотеки матричных методов: если вы хотите повысить производительность, то можете опустить проверку входных параметров на ошибки ценой большего риска получения исключения. (Чтобы сохранить лаконичность этой статьи, большинство проверок на ошибки было удалено.) Другое преимущество в том, что вы можете модифицировать свою библиотеку для оптимизации именно вашего сценария. Основной недостаток — создание собственной библиотеки отнимает больше времени, чем использование существующей.

Другой удобный метод, который стоит добавить в вашу библиотеку для работы с матрицами, — такой, который позволяет отображать матрицу как строку. Вот один из возможных вариантов:

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

Возможно, вы захотите параметризовать количество отображаемых десятичных разрядов, дополнение ширины столбцов или и то, и другое.

Перемножение матриц

Microsoft .NET Framework 4 и более поздние версии предоставляют хороший способ значительно повысить производительность метода перемножения матриц. Перемножение матриц иллюстрируется на рис. 2.

Рис. 2. Перемножение матриц

Заметьте, что вычисление значения каждой ячейки не зависит от значений других ячеек, поэтому каждое вычисление независимо, и потенциально их можно выполнять параллельно на компьютере с несколькими процессорами (или ядрами). Вот стандартный подход к перемножению матриц:

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) // каждая строка A
    for (int j = 0; j < bCols; ++j) // каждый столбец B
      for (int k = 0; k < aCols; ++k)
        result[i][j] += matrixA[i][k] * matrixB[k][j];
  return result;
}

Task Parallel Library упрощает кодирование простой распараллеленной версии перемножения матриц.

Поскольку перемножение матриц является операцией порядка O(n^3), производительность может быть проблемой. Так, если размер матрицы A равен 300 × 200, а матрицы B — 200 × 400, вычисление результата произведения A и B требует 300 * 200 * 400 = 24 000 000 операций умножения. Task Parallel Library (TPL) в пространстве имен System.Threading.Tasks в .NET Framework 4 и более поздних версиях упрощает кодирование простой распараллеленной версии перемножения матриц. Один из возможных вариантов выглядит так:

static double[][] MatrixProduct(double[][] matrixA,
  double[][] matrixB)
{
  // Проверка ошибок, вычисление 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;
}

В этой версии вычисления делятся по строкам. «За кулисами» TPL генерирует весь сложный код инфраструктуры синхронизации для параллельных вычислений на нескольких процессорах.

Проверка согласованности

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

static double[][] MatrixRandom(int rows, int cols,
  double minVal, double maxVal, int seed)
{
  // Возвращаем матрицу со значениями
  // в диапазоне от minVal до 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;
}

В дополнение предположим, что ваша матрица создает матрицу тождественности, или единичную матрицу (identity matrix), т. е. квадратную матрицу со значениями 1.0 по основной диагонали и значениями 0.0 в остальных ячейках:

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

И у вас есть метод, сравнивающий две матрицы на тождественность:

static bool MatrixAreEqual(double[][] matrixA,
  double[][] matrixB, double epsilon)
{
  // True, если все значения в A равны
  // соответствующим значениям в B
  int aRows = matrixA.Length;
  int bCols = matrixB[0].Length;
  for (int i = 0; i < aRows; ++i) // каждая строка A и B
    for (int j = 0; j < aCols; ++j) // каждый столбец A и B
      if (Math.Abs(matrixA[i][j] - matrixB[i][j]) > epsilon)
        return false;
  return true;
}

Заметьте, что метод MatrixAreEqual не сравнивает значения ячеек на точное равенство, так как эти значения имеют тип double. Вместо этого метод проверяет, очень ли близки друг к другу значения ячеек (в пределах epsilon).

Поскольку произведение любой квадратной матрицы m и единичной матрицы той же размерности равно исходной матрице m, вы можете протестировать метод MatrixProduct наряду со следующими строками кода:

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");
Проверка согласованности хорошо поддается тестированию случайным вводом.

Разложение матрицы

При разложении берется квадратная матрица M и вычисляются две новые квадратные матрицы, которые при перемножении дают исходную матрицу M. Идея аналогичная обычному разложению чисел на множители: число 6 можно разложить на множители 2 и 3, потому что 2 * 3 = 6. Поначалу может показаться, что в разложении матриц нет особого смысла, но оказывается, что такое разложение существенно упрощает очень трудную задачу обращения матрицы. Если много видов разложения матриц, и каждый из них может вычисляться по нескольким алгоритмам. Методика, представленная в этой статье, называется разложением LUP и использует метод Дулитла (Doolittle's method) с выбором ведущего элемента столбца (partial pivoting).

Чтобы понять разложение LUP, полезно сначала разобраться в простом разложении LU, которое было введено знаменитым математиком Аланом Тьюрингом (Alan Turing). Допустим, у вам имеется вот такая матрица M размером 4 × 4:

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

Одно из возможных разложений LU для 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

и 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

Это работает, так как L * U = M. Заметьте, что нижняя матрица L имеет значения 1.0 по диагонали и 0.0 вверху справа. Иначе говоря, значимые значения ячеек нижней матрицы находятся внизу слева. Аналогично значимые значения ячеек верхней матрицы расположены на основной диагонали и вверху справа.

Также обратите внимание на отсутствие перекрытия местонахождений значимых значений ячеек в L и U. Таким образом, вместо генерации двух матриц результатов, L и U, алгоритм разложения матрицы обычно сохраняет как нижние, так и верхние результаты в одной матрице, которая хранит L и U для экономии памяти.

Разложение LUP матрицы — небольшая, но важная вариация разложения LU. При разложении LUP берется матрица matrix M и создаются матрицы L и U, а также массив P. Произведение L и U в LUP не точно соответствует исходной матрице M, а является версией M, где некоторые строки переупорядочены. Информация о том, как эти строки были переупорядочены, хранится в массиве P; эту информацию можно использовать для реконструкции исходной матрицы M.

Близкий вариант разложения Дулитла, представленного в этой статье, называется разложением Краута (Crout). Главное различие между методами Дулитла и Краута в том, что Дулитл помещает значения 1.0 на основную диагональ матрицы L, а Краут — на основную диагональ матрицы U.

Причина, по которой разложение LUP используется чаще, чем разложение LU, весьма тонкая. Как вы вскоре увидите, разложение используется для обращения матрицы. Однако, когда разложение матрицы применяется как вспомогательный метод для ее обращения, оказывается, что инверсия заканчивается неудачей, если по основной диагонали матрицы LU имеется значение 0.0. Поэтому в разложении LUP, когда на основной диагонали появляется значение 0.0, алгоритм меняет местами две строки, чтобы сместить значение 0.0 с диагонали, и отслеживает, какие строки были переставлены, в массиве P.

Метод разложения матрицы показан на рис. 3.

Рис. 3. Метод разложения матрицы

static double[][] MatrixDecompose(double[][] matrix,
  out int[] perm, out int toggle)
{
  // Разложение LUP Дулитла. Предполагается,
  // что матрица квадратная.
  int n = matrix.Length; // для удобства
  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) // каждый столбец
  {
    double colMax = Math.Abs(result[j][j]); // Наибольшее значение в столбце 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) // перестановка строк
    {
      double[] rowPtr = result[pRow];
      result[pRow] = result[j];
      result[j] = rowPtr;
      int tmp = perm[pRow]; // Меняем информацию о перестановке
      perm[pRow] = perm[j];
      perm[j] = tmp;
      toggle = -toggle; // переключатель перестановки строк
    }
    if (Math.Abs(result[j][j]) < 1.0E-20)
      return null;
    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];
    }
  } // основной цикл по столбцу j
  return result;
}

Этот метод можно было бы вызвать так:

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

Метод MatrixDecompose принимает квадратную матрицу и возвращает три значения. Явно возвращается переставленная матрица LU. Остальные два значения возвращаются через выходные параметры. Один из них является массивом перестановки (permutation array), где хранится информация о том, как переставлены строки. Второй выходной параметр — переключатель, принимающий значение либо +1, либо –1 в зависимости от того, было количество перестановок четным (+1) или нечетным (–1). Значение переключателя не используется для обращения матрицы, но необходимо, если разложение применяется для вычисления определителя матрицы.

Метод MatrixDecompose весьма замысловатый, но на самом деле для модификации кода достаточно понимать в нем лишь несколько деталей. Версия, представленная здесь, выделяет новую память под возвращаемую матрицу LU через вспомогательный метод MatrixDuplicate:

static double[][] MatrixDuplicate(double[][] matrix)
{
  // Предполагается, что матрица не нулевая
  double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
  for (int i = 0; i < matrix.Length; ++i) // Копирование значений
   for (int j = 0; j < matrix[i].Length; ++j)
      result[i][j] = matrix[i][j];
  return result;
}

Альтернативный подход — помещение результата во входную матрицу для экономии памяти. Семантика C# требует сделать параметр matrix передаваемым по ссылке (ref), потому что он используется и для ввода, и для вывода. При таком подходе сигнатура метода выглядела бы следующим образом:

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

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

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

Возможно, вы предпочтете исключить параметр toggle, чтобы упростить сигнатуру метода, если вы не собираетесь использовать разложение для вычисления определителя матрицы.

Другая часть MatrixDecompose, которую вам, может быть, захочется изменить, — это выражение:

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

Этот код означает следующее: если даже после обмена двух строк из-за того, что на основной диагонали было значение 0.0, там все равно остается 0.0, то вернуть null». Возможно, вы захотите сменить произвольное значение epsilon с 1.0E-20 на какое-то другое значение в зависимости от уровня точности, используемого в вашей системе. И вместо возврата null можно генерировать исключение; если бы этот метод вызывался в процессе обращения матрицы, данный процесс заканчивался бы неудачей. Наконец, если вы используете разложение матрицы для целей, отличных от ее обращения, то можно вообще исключить это выражение.

Обращение матрицы

Чтобы использовать разложение для обращения матрицы, крайне важно написать вспомогательный метод, решающий систему уравнений. Этот метод представлен на рис. 4.

Рис. 4. Метод HelperSolve

static double[] HelperSolve(double[][] luMatrix,
 double[] b)
{
  // Решаем 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;
}

Метод HelperSolve находит массив x, который при умножении на матрицу LU дает массив b. Этот метод весьма непрост, и вы можете полностью понять его, только проследив его работу на нескольких примерах. В нем два цикла. Первый цикл использует прямую подстановку (forward substitution) в нижней части матрицы LU, а второй цикл — обратную подстановку (backward substitution) в верхней части той же матрицы. В некоторых других реализациях разложения матрицы аналогичный метод называют, например, luBackSub.

Хотя этот код короткий, он весьма сложен, но никаких сценариев, где вам могло бы понадобиться его изменение, нет. Заметьте, что HelperSolve принимает матрицу LU от MatrixDecompose, но не использует выходной параметр perm. Это подразумевает, что HelperSolve, по сути, является вспомогательным методом и нуждается в дополнительной оболочке кода для решения системы уравнений. При рефакторинге кода из этой статьи под объектно-ориентированное программирование вы, вероятно, предпочтете сделать метод HelperSolve закрытым.

Создав метод HelperSolve, можно реализовать метод обращения матрицы, как показано на рис. 5.

Рис. 5. Метод 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;
}

И вновь сложный код. Алгоритм обращения основан на том, что результат перемножения матрицы M и ее обращения является единичной матрицей. Метод MatrixInverse в основном отвечает за решение системы уравнений Ax = b, где A — матрица разложения LU , а константы b равны либо 1.0, либо 0.0 и соответствуют единичной матрице. Заметьте, что MatrixInverse использует массив perm, возвращаемый после вызова MatrixDecompose.

Вызов метода MatrixInverse мог бы выглядеть так:

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

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

Определитель матрицы

Когда у вас есть метод разложения матрицы, написать метод для вычисления определителя матрицы нетрудно:

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

Как это ни удивительно, оказывается, что определитель матрицы — это просто результат перемножения значений на основной диагонали матрицы разложения LUP со знаком «плюс» или «минус» в зависимости от значения toggle. Обратите внимание на неявное преобразование типа значения toggle из int в double. Помимо добавления проверки ошибок в MatrixDeterminant, вы должны добавить прямой возврат в ситуациях, где размер входной матрицы равен 1 (возвращается единственное значение) или 2 × 2 (возвращается a*d – b*c). Вызов этого метода мог бы выглядеть так:

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

Решение систем уравнений

Метод HelperSolve можно легко адаптировать для решения системы линейных уравнений:

static double[] SystemSolve(double[][] A, double[] b)
{
  // Решаем Ax = b
  int n = A.Length;
  int[] perm;
  int toggle;
  double[][] luMatrix = MatrixDecompose(
    A, out perm, out toggle);
  if (luMatrix == null)
    return null; // или исключение
  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;
}

Вот код, который дал вывод, показанный на рис. 1, для решения следующей системы:

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

Заметьте, что SystemSolve переупорядочивает свой входной параметр b, используя массив perm из MatrixDecompose до вызова HelperSolve.

Что такое массив перестановок

Последние несколько строк вывода на рис. 1 указывают, что матрицы L и U можно перемножить так, чтобы получить исходную матрицу. Знание того, как это делается, не поможет вам в решении практических задач операций над матрицами, но позволит разобраться, что представляет собой часть P в разложении LUP. Восстановление исходной матрицы из ее компонентов L и U также пригодится для тестирования ваших библиотечных методов работы с матрицами на согласованность.

Один из способов восстановления исходной матрицы после разложения LUP — перемножение L и U с последующей перестановкой строк результата на основе массива P:

// Создаем матрицу m.
// Вызываем MatrixDecompose, возвращаем lu и perm.
// Извлекаем нижнюю часть lu как матрицу 'lower'.
// Извлекаем верхнюю часть lu как матрицу '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");

Метод UnPermute можно закодировать так:

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; // создаем массив unperm
  for (int r = 0; r < luProduct.Length; ++r) //по каждой строке
    result[r] = luProduct[unperm[r]];
  return result;
}

Второй подход — преобразование массива perm в матрицу perm с последующим перемножением матрицы perm и комбинированной матрицы LU:

// Создаем матрицу m.
// Вызываем MatrixDecompose, возвращаем lu и perm.
// Извлекаем нижнюю часть lu как матрицу 'lower'.
// Извлекаем верхнюю часть lu как матрицу '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");

Матрица perm является квадратной с одним значением 1.0 в каждой строке и каждом столбце. Метод, который создает матрицу perm из массива perm, можно написать следующим образом:

static double[][] PermArrayToMatrix(int[] perm)
{
  // Массив perm преобразуем
  // в соответствующую матрицу по Дулитлу
  int n = perm.Length;
  double[][] result = MatrixCreate(n, n);
  for (int i = 0; i < n; ++i)
    result[i][perm[i]] = 1.0;
  return result;
}

Заключение

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