Введение

Во многих областях науки и техники, при решении различных задач, возникает необходимость расчета систем линейных уравнений, которые иногда являются сильно разреженными. Поскольку обычно такие системы уравнений имеют большой порядок, то найти их решение (в сколь-нибудь обозримые сроки) возможно только с применением вычислительной техники. Для решения подобных задач можно использовать прямые методы (Гаусса, Крамера), и приближенные (Зейделя, градиентного спуска и др.). Прямые методы дают точное решение (ошибка может быть только из-за выхода величин при преобразованиях за пределы точности типов для используемой платформы), но обычно требуют большее количество итераций. Приближенные методы обычно работают быстрее, но могут требовать особых начальных условий или какой-либо подготовки исходных данных, и решение не всегда сходится. В данной статье рассматривается реализация метода Гаусса, модифицированная для решения разреженных систем линейных уравнений. Данный метод не является «математическим», а скорее «алгоритмическим» и показывает, как путем простых приемов можно ускорить расчет в сотни раз.

Условия применения

Рассматриваемый метод подходит для любых невырожденных (имеющих одно решение) систем линейных уравнений. В дальнейшем под словом матрица будет пониматься матрица коэффициентов линейных уравнений. Единственное условие - метод имеет смысл только для сильно разреженных матриц больших размерностей N>1000, где не менее двух третей коэффициентов равны 0. Для обычных матриц метод не даст выигрыша ни в скорости, ни в экономии памяти.

Описание метода

Метод Гаусса заключается в обнулении нижней (или верхней) полуматрицы a[i,k], где строка i=2..N, столбец k=1..i-1. Используя известное правило, что решение системы не изменится, если к одному из ее уравнений добавить другое, умноженное на произвольный коэффициент, строки матрицы поочередно, начиная с первой, домножаются на специальные коэффициенты и суммируются со всеми ниже лежащими строками. Коэффициенты подбираются таким образом, чтобы каждая строка обнуляла соответствующий ей (i=k) ниже лежащий столбец (только его поддиагональную часть). Так первая строка обнуляет первый столбец A[2..N,1], вторая – второй A[3..N,2] и т.д. до предпоследней строки. Коэффициенты рассчитываются так: например для того чтобы обнулить первый столбец нужно первую строку домножить на k=-A[2,1]/A[1,1] и прибавить ко второй, затем на k=-A[3,1]/A[1,1] и прибавить к третей и т.д. После операция повторяется со второй строкой (она домножается на соответствующие коэффициенты и складывается со всеми нижележащими строками) и т.д. до предпоследней. В результате получается матрица, под главной диагональю которой остаются одни нули. Теперь в последнем уравнении остался только один коэффициент – соответственно можно легко вычислить значение последней неизвестной: X[N]=B[N]/A[N,N], где B – вектор свободных членов уравнений, который кстати, также участвует в преобразованиях при обнулении матрицы A. В предпоследнем уравнении осталось два коэффициента, и последняя неизвестная уже известна, можно вычислить предпоследнюю: X[N-1]=(B[N-1]-A[N-1,N]*X[N])/A[N-1,N-1]. Таким образом, используя обратный ход, можно вычислить все неизвестные. В процессе преобразования матрицы не должны обнуляться коэффициенты главной диагонали. Если такое происходит, то матрица не была факторизована или неправильно составлены уравнения. Или же данная система имеет не одно решение. Более подробно о методе Гаусса можно узнать здесь.

Приведенный выше алгоритм легко реализуем. Количество итераций известно заранее и зависит только от порядка матрицы. Если его запрограммировать в классическом стиле, то он легко будет справляться с небольшими матрицами N~1000. Однако матрицы иногда бывают порядка 100 000 и более. На расчет такой матрицы классическим методом Гаусса уйдут годы, поскольку количество итераций будет огромным. Если данную матрицу хранить в обычном двумерном массиве, то при использовании типа double она будет занимать в памяти 100 000 х 100 000 х 8 = 800 000 000 байт = 800 Мб что совсем не мало. А если учесть что размер матрицы возрастает пропорционально квадрату порядка матрицы, то для матрицы порядка 400 000 это будет в 16 раз больше, т.е. 12.8 Гб!

На основании выше приведенных рассуждений можно выяснить два направления оптимизации расчета:

• Более компактное хранение информации.
• Уменьшение количества итераций.

Для разреженных матриц возможно и то, и другое:

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

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

• для главной диагонали это одна координата d, которая равна соответствующему индексу элемента главной диагонали
• для верхней и нижней полуматрицы координаты d,j, где j – смещение от главной диагонали вправо или вниз (как показано на рисунке)

Диаграмма классов

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

Реализация метода

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

Реализация интерфейса строк и столбцов:

public interface IRow
{
      // записывает значение value в ячейку с индексом num
      void setValue(int num, double value);

      // добавляет value к существующему значению в ячейке num
      void addValue(int num, double value);

      // возвращает значение ячейки num
      double getValue(int num);

      // возвращает все ненулевые ячейки строки/столбца: 
      // индексы ячеек - в indexes, значения в values
      void getValues(ref int[] indexes, ref double[] values);
}

Реализация интерфейса матрицы:

public interface IMatrix
{
	  // устанавливает значение value в ячейку с координатами [row,col];
	  // row - номер строки матрицы
	  // col - номер столбца матрицы
	  void setValue(int row, int col, double value);

      // добавляет значение value к ячейке [row,col]
      void addValue(int row, int col, double value);

      // возвращает значение ячейки [row,col]
      double getValue(int row, int col);

      // возвращает порядок матрицы
      int getN();

      // возвращает ненулевые значения и индексы ячеек строки d,
      // которые находятся правее главной диагонали
      void getJRow(int d, ref int[] indexes, ref double[] values);

       // возвращает ненулевые значения и индексы ячеек столбца d, 
      // которые находятся ниже главной диагонали
      void getJCol(int d, ref int[] indexes, ref double[] values);
}

Абстрактный класс матрицы, который содержит поля и методы, идентичные в различных реализациях классов матриц:

public abstract class Matrix : IMatrix
{
      // порядок матрицы
      private int N;
	  public Matrix(int N)
      {
            this.N = N;
      }

      // возвращает порядок матрицы
      public int getN()
      {
            return N;
      }
      // методы интерфеса IMatrix имплементируемые в реализациях
      public abstract void setValue(int row, int col, double value);
      public abstract void addValue(int row, int col, double value);
      public abstract double getValue(int row, int col);
      public abstract void getJRow(int d, ref int[] indexes, ref double[] values);
      public abstract void getJCol(int d, ref int[] indexes, ref double[] values);

      // сервисные методы для преобразования координат внутри матрицы
      // возвращает 0 при x<=0, x при x>0
      protected static int sigma(int x)
      {
            return (x>0)?x:0;
      }
      // преобразование row,col в координаты полуматриц (d,j)
      // для верхней полуматрицы: d=0..(n-1), j=-(n-d-1)..-1
      // для нижней полуматрицы:  d=0..(n-1), j=1..(n-d-1)
      protected static int getD(int row, int col)
      {
            return row-sigma(row-col);
      }
      protected static int getJ(int row, int col)
      {
            return col-row;
      }

      // обратное преобразование координат полуматриц (d,j) в row,col
      protected static int getRow(int d, int j)
      {
            return d + sigma(-j);
      }
      protected static int getCol(int d, int j)
      {
            return d + sigma(j);
      }
}

На данном этапе нужно решить в каком формате хранить значения ячеек матрицы. Для очень сильно разреженных матриц, когда в каждой строке/столбце матрицы ненулевых значений совсем мало, данные лучше хранить в обычных массивах, так как поиск в небольшом массиве работает быстро. Если же матрица не сильно разрежена, то лучше использовать хэш-таблицу. Ниже приведены реализации интерфейса IRow для обоих случаев:

Реализация интерфейса строки c использованием массива для хранения значений

public class Array : IRow
{
      // массив индексов (номера ненулевых ячеек)
      private ArrayList indexes = new ArrayList();

      // массив значений ячеек
      private ArrayList values = new ArrayList();

      public void setValue(int num, double value)
      {
            int i = indexes.IndexOf(num);

            // если записывается 0-вое значение, то удаляем данную ячейку
            if (value==0)
            {
                  if (i!=-1)
                  {
                        indexes.RemoveAt(i);
                        values.RemoveAt(i);
                  }
                  return;
            }

            // если значение не 0-вое, то перезаписываем или добавляем ячейку
            if (i!=-1)
            {
                  values[i] = value;
            }
            else
            {
                  indexes.Add(num);
                  values.Add(value);
            }
      }

     public void addValue(int num, double value)
      {
            int i = indexes.IndexOf(num);
            if (i!=-1)
            {
                  values[i] = (double)values[i] + value;
            }
            else
            {
                  indexes.Add(num);
                  values.Add(value);
            }
      }


      public double getValue(int num)
      {
            int i = indexes.IndexOf(num);
            if (i!=-1) return (double)values[i];
            return 0;
      }

      public void getValues(ref int[] indexes, ref double[] values)
      {
         indexes = (int[])this.indexes.ToArray(System.Type.GetType("System.Int32"));
         values = (double[])this.values.ToArray(System.Type.GetType("System.Double"));
      }
}

IРеализация интерфейса строки c использованием хэш-таблицы для хранения значений:

public class  Hash : IRow
{
      // хэш-таблица в формате номер_ненулевой_ячейки = значение_ячейки
      private Hashtable hash = new Hashtable();

      public void setValue(int num, double value)
      {
            // если значение 0-вое, то удаляем данную ячейку
            if (value==0)
            {
                  if (hash.ContainsKey(num)) hash.Remove(num);
                  return;
            }
            // если значение не 0-вое, то перезаписываем или добавляем ячейку
            hash[num] = value;
      }

      public void addValue(int num, double value)
      {
            if (hash.ContainsKey(num)) hash[num] = (double)hash[num] + value;
              else hash[num] = value;
      }

      public double getValue(int num)
      {
            if (hash.ContainsKey(num)) return (double)hash[num];
            return 0;
      }

      public void getValues(ref int[] indexes, ref double[] values)
      {
            System.Collections.ICollection keys = hash.Keys;
            int count = keys.Count;
            indexes = new int[count];
            values = new double[count];
            int i=0;
            for (IEnumerator it = keys.GetEnumerator(); it.MoveNext();)
            {
                  indexes[i] = (int)it.Current;
                  values[i++] = (double)hash[it.Current];
            }
      }
}

Полная реализация класса "матрица":

public class HashMatrix : Matrix
{
      // тип хранения данных: массив или хэш-таблица
      public const byte ARRAY_TYPE = 0;
      public const byte HASH_TYPE = 1;

      // jp - строки верхней полуматрицы
      // jm - столбцы нижней полуматрицы
      // dd - главная диагональ
      private IRow[] jp,jm;
      private double[] dd;

      // конструктор матрицы
      // type - тип хранения данных ARRAY_TYPE/HASH_TYPE
      // N - порядок матрицы
      public HashMatrix(byte type, int N) : base(N)
      {
            dd = new double[N];
            if (type==HASH_TYPE)
            {
                  jp = new Hash[N];
                  jm = new Hash[N];
                  for (int i=0; i<N; i++)
                        {
                             jp[i] = new Hash();
                             jm[i] = new Hash();
                        }
            }
            else
            {
                  jp = new Array[N];
                  jm = new Array[N];
                  for (int i=0; i<N; i++)
                  {
                        jp[i] = new Array();
                        jm[i] = new Array();
                  }
            }
      }

      public override void setValue(int row, int col, double value)
      {
            if (row==col)
            {
                  dd[row] = value;
                  return;
            }
            int d = getD(row,col);
            int j = getJ(row,col);
            if (j>0) jp[d].setValue(j,value);
				else jm[d].setValue(-j,value);
      }

      public override void addValue(int row, int col, double value)
      {
            if (value==0) return;
            if (row==col)
			{
                  dd[row] += value;
                  return;
            }
            int  d = getD(row,col);
            int  j = getJ(row,col);
            if  (j>0) jp[d].addValue(j, value);
            else  jm[d].addValue(-j,value);
      }

      public override double  getValue(int  row, int  col)
      {
            if  (row==col) return  dd[row];
            int  d = getD(row,col);
            int  j = getJ(row,col);
            return  (j>0)?jp[d].getValue(j):jm[d].getValue(-j);
      }

      public override void  getJRow(int  d, ref int [] indexes, ref double [] values)
      {
            jp[d].getValues(ref  indexes, ref  values);
            for  (int  i=0; i<indexes.Length; i++) indexes[i] = getCol(d,indexes[i]);
      }
      public override void  getJCol(int  d, ref int [] indexes, ref double [] values)
      {
            jm[d].getValues(ref indexes, ref values);
            for  (int  i=0; i<indexes.Length; i++) indexes[i] = getRow(d,-indexes[i]);
      }
}

И, наконец, реализация метода Гаусса:

public class  Gauss
{
      // матрица
      private IMatrix matrix;

      // конструктор, принимает созданную матрицу коэффициентов
      public  Gauss(IMatrix matrix)
      {
            this .matrix = matrix;
      }

      // главный метод, возвращающий решение, принимает вектор свободных членов
      public double [] calculate(double [] B)
      {
            // обнуляем нижнюю полуматрицу, перебирая сверху вниз все строки
            // и складывая каждую со всеми нижележащими
            for  (int  row=0; row<(matrix.getN()-1); row++)
            {
                  // get all non-zero values of zeroing column
                  int [] colIndexes = new int [0];
                  double [] colValues = new double [0];
                  matrix.getJCol(row,ref  colIndexes,ref  colValues);

                  // получаем все ненулевые значения обнуляемого столбца
                  // получаем индексы и значения ячеек строки, правее главной диагонали
                  int [] rowIndexes = new int [0];
                  double [] rowValues = new double [0];
                  matrix.getJRow(row, ref  rowIndexes, ref  rowValues);

                  // получаем элемент главной диагонали, которым будем обнулять столбец
                  double  dd = matrix.getValue(row,row);
                  for  (int  i=0; i// высчитываем коэффициент 
                     double  k = colValues[i]/dd;

                     // k подобран таким образом чтобы обнулить ячейку столбца,
                     matrix.setValue(colIndexes[i],row,0);

                     // складываем строки
                     for  (int  ii=0; ii<rowIndexes.Length; ii++)
                     {
                        matrix.addValue(colIndexes[i],rowIndexes[ii],-rowValues[ii]*k);
                     }

                     // складываем соответствующие свободные члены
                     B[colIndexes[i]] -= B[row]*k;
                  }
            }
            // создаем вектор неизвестных
            double [] x = new double [matrix.getN()];

            // используя обратный ход находим неизвестные
            for  (int  row = (matrix.getN()-1); row>=0; row--)
            {
               double  e = 0;
               int [] indexes = new int [0];
               double [] values = new double [0];
               matrix.getJRow(row,ref  indexes,ref  values);
               for  (int  i=0; i<indexes.Length; i++) e += x[indexes[i]]*values[i];
               x[row] = (B[row]-e)/matrix.getValue(row,row);
            }
            return  x;
      }
}

Заключение

В статье рассмотрен метод решения систем линейных уравнений по алгоритму Гаусса, и приведен код реализующий этот алгоритм для C#. Хотя данный код и является полным, для использования в пользовательских приложениях он нуждается в следующей доработке:

• Как минимум нужно сделать обработку исключительных ситуаций, например таких как деление на ноль или выход индекса за пределы массива
• Реализовать расчетный алгоритм в отдельном потоке команд. Если вызывать метод calculate из основного UI потока, это остановит прием сообщений окружения («повесит программу») до момента завершения расчета
• Если необходимо рассчитать несколько систем линейных уравнений, которые отличаются только правой частью (свободными членами), то естественно нет необходимости каждый раз обнулять нижнюю полуматрицу. Для этого нужно доработать класс Gauss таким образом, чтобы он, обнуляя матрицу, сохранял историю изменений. После этого для расчета неизвестных нужно будет преобразовать только вектор свободных членов, передаваемых в метод calculate в соответствии с историей. Историю можно хранить в отдельной матрице, или в нижней полуматрице основной матрицы, которая все равно не используется при обратном ходе
• Если данный алгоритм планируется реализовать на платформе, где возможна прямая работа с указателями (Delphi, C++), то лучше по возможности использовать их, это позволит в некоторых местах выполнять меньшее количество преобразований
• Для некоторых матриц необходима доработка алгоритма факторизации, для исключения попадания нулей на главную диагональ
• В реализациях интерфейса IRow, при создании экземпляров массивов и хэш-таблиц, в явном виде не задавалось свойство capacity. На практике, чтобы оптимально использовать быстродействие и экономить память, нужно подобрать компромисное значение capacity. Так с увеличением capacity увеличивается быстродействие и расход памяти, с уменьшением – наоборот. Оптимальное значение capacity зависит от структуры матрицы и не превышает максимального количества ненулевых элементов в строке или столбце матрицы

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

Литература

1. Эндрю Троелсен. С# и платформа .NET. Библиотека программиста. — СПб.: Питер, 2004.



2. Техническая документация MSDN.