MAXimal

добавлено: 16 Mar 2011 0:51
редактировано: 16 Mar 2011 0:51

Вычисление определителя методом Краута за O (N3)

Здесь будет рассмотрена модификация метода Краута (Crout), позволяющая вычислить определитель матрицы за O (N3).

 

Собственно алгоритм Краута находит разложение матрицы A в виде A = L U, где L - нижняя, а U - верхняя треугольная матрицы. Без ограничения общности можно считать, что в L все диагональные элементы равны 1. Но, зная эти матрицы, легко вычислить определитель A: он будет равен произведению всех элементов, стоящих на главной диагонали матрицы U.

Имеется теорема, согласно которой любая обратимая матрица обладает LU-разложением, и притом единственным, тогда и только тогда, когда все её главные миноры отличны от нуля. Следует напомнить, что мы рассматриваем только такие разложения, в которых диагональ L состоит только из единиц; иначе же, вообще говоря, разложение не единственно.

 

Пусть A - матрица, N - её размер. Мы найдём элементы матриц L и U.

Сам алгоритм состоит из следующих шагов:

  1. Положим Li i = 1 для i = 1, 2, ..., N
  2. Для каждого j = 1, 2, ..., N выполним:
    1. Для i = 1, 2, ..., j найдём значение Ui j:
      Ui j = Ai j - SUM Li k Uk j ,
      где сумма по всем k = 1, 2, ..., i-1.
    2. Далее, для i = j+1, j+2, ..., N имеем:
      Li j = (Ai j - SUM Li k Uk j) / Uj j ,
      где сумма берётся по всем k = 1, 2, ..., j-1.

Реализация

Код на Java (с использованием дробной длинной арифметики):

static BigInteger det (BigDecimal a [][], int n)
{

	try {

	for (int i=0; i<n; i++)
	{
		boolean nonzero = false;
		for (int j=0; j<n; j++)
			if (a[i][j].compareTo (new BigDecimal (BigInteger.ZERO)) > 0)
				nonzero = true;
		if (!nonzero)
			return BigInteger.ZERO;
	}

	BigDecimal scaling [] = new BigDecimal [n];
	for (int i=0; i<n; i++)
	{
		BigDecimal big = new BigDecimal (BigInteger.ZERO);
		for (int j=0; j<n; j++)
			if (a[i][j].abs().compareTo (big) > 0)
				big = a[i][j].abs();
		scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
			(big, 100, BigDecimal.ROUND_HALF_EVEN);
	}

	int sign = 1;

	for (int j=0; j<n; j++)
	{

		for (int i=0; i<j; i++)
		{
			BigDecimal sum = a[i][j];
			for (int k=0; k<i; k++)
				sum = sum.subtract (a[i][k].multiply (a[k][j]));
			a[i][j] = sum;
		}

		BigDecimal big = new BigDecimal (BigInteger.ZERO);
		int imax = -1;
		for (int i=j; i<n; i++)
		{
			BigDecimal sum = a[i][j];
			for (int k=0; k<j; k++)
				sum = sum.subtract (a[i][k].multiply (a[k][j]));
			a[i][j] = sum;
			BigDecimal cur = sum.abs();
			cur = cur.multiply (scaling[i]);
			if (cur.compareTo (big) >= 0)
			{
				big = cur;
				imax = i;
			}
		}

		if (j != imax)
		{
			
			for (int k=0; k<n; k++)
			{
				BigDecimal t = a[j][k];
				a[j][k] = a[imax][k];
				a[imax][k] = t;
			}
			
			BigDecimal t = scaling[imax];
			scaling[imax] = scaling[j];
			scaling[j] = t;

			sign = -sign;
		}

		if (j != n-1)
			for (int i=j+1; i<n; i++)
				a[i][j] = a[i][j].divide
					(a[j][j], 100, BigDecimal.ROUND_HALF_EVEN);
		
	}

	BigDecimal result = new BigDecimal (1);
	if (sign == -1)
		result = result.negate();
	for (int i=0; i<n; i++)
		result = result.multiply (a[i][i]);

	return result.divide
		(BigDecimal.valueOf(1), 0, BigDecimal.ROUND_HALF_EVEN).toBigInteger();

	}
	catch (Exception e)
	{
		return BigInteger.ZERO;
	}

}