Эффективный способ получить определитель n! х н! матрица в клен

У меня большая матрица, n! x n!, для которого мне нужно взять определитель. Для каждой перестановки n я связываю

  • вектор длины 2n (это легко вычислить)
  • многочлен от 2n переменных (произведение линейных множителей, вычисляемых рекурсивно по n)

Матрица представляет собой оценочную матрицу полиномов на векторах (рассматриваемых как точки). Таким образом, элемент сигма, тау матрицы (индексированный перестановками) представляет собой полином для сигмы, оцененный в векторе для тау.

Пример. Для n=3, если i-й полином равен (x1 - 4)(x3 - 5)(x4 - 4)(x6 - 1), а j-я точка равна (2,2,1,3,5,2), то (i,j)-й элемент матрицы будет (2 - 4)(1 - 5)(3 - 4)(2 - 1) = -8. Здесь n=3, поэтому точки находятся в R^(3!) = R^6, а полиномы имеют 3!=6 переменных.

Моя цель - определить, является ли матрица неособой.

Мой подход сейчас таков:

  • функция point выполняет перестановку и выводит вектор
  • функция poly выполняет перестановку и выводит полином
  • функция nextPerm дает следующую перестановку в лексикографическом порядке

Сокращенная версия псевдокода моего кода такова:

B := [];
P := [];
w := [1,2,...,n];
while w <> NULL do
  B := B append poly(w);
  P := P append point(w);
  w := nextPerm(w);
od;

// BUILD A MATRIX IN MAPLE
M := Matrix(n!, (i,j) -> eval(B[i],P[j]));

// COMPUTE DETERMINANT IN MAPLE
det := LinearAlgebra[Determinant]( M );

// TELL ME IF IT'S NONSINGULAR
if det = 0 then return false;
else return true; fi;

Я работаю в Maple, используя встроенную функцию LinearAlgebra[Determinant], но все остальное — это пользовательская функция, которая использует низкоуровневые функции Maple (например, seq, convert и cat).

Моя проблема в том, что это занимает слишком много времени, то есть я могу набраться терпения до n=7, но чтобы получить n=8, нужны дни. В идеале я хочу иметь возможность добраться до n=10.

У кого-нибудь есть идея, как я могу улучшить время? Я открыт для работы на другом языке, например. Matlab или C, но предпочел бы найти способ ускорить это в Maple.

Я понимаю, что на это может быть трудно ответить без всех кровавых подробностей, но код для каждой функции, например. point и poly уже оптимизированы, поэтому реальный вопрос здесь заключается в том, есть ли более быстрый способ получить определитель путем построения матрицы на лету или что-то в этом роде.

ОБНОВЛЕНИЕ: вот две идеи, с которыми я играл, которые не работают:

  1. Я могу сохранить полиномы (поскольку их вычисление занимает некоторое время, я не хочу переделывать это, если смогу) в вектор длины n!, вычислить точки на лету и подставить эти значения в перестановку. формула определителя:

    формула определителя перестановки

    Проблема здесь в том, что это O(N!) по размеру матрицы, поэтому для моего случая это будет O((n!)!). Когда n=10, (n!)! = 3,628,800! слишком много, чтобы даже думать об этом.

  2. Вычислите определитель, используя разложение LU. К счастью, главная диагональ моей матрицы отлична от нуля, так что это возможно. Поскольку это O(N^3) в размере матрицы, это становится O((n!)^3), что намного ближе к выполнимому. Проблема, однако, в том, что мне нужно хранить всю матрицу, что создает серьезную нагрузку на память, не говоря уже о времени выполнения. Так что это тоже не работает, по крайней мере, без ума. Есть идеи?


algorithm permutation determinants maple
person Daniel    schedule 21.02.2013    source источник
comment
Из какой области взяты коэффициенты ваших полиномов и оценочные баллы? Ваши примеры показывают целые числа - это упрощение или это действительно так?   -  person Erik P.    schedule 05.10.2013

Ответы (2)


arrow_upward
2
arrow_downward

Мне не ясно, является ли ваша проблема пространством или временем. Очевидно, что они торгуются туда и обратно. Если вы только хотите знать, является ли определитель положительным или нет, вам определенно следует использовать LU разложение. Причина в том, что если A = LU с L нижним треугольным и U верхним треугольным, то

det(A) = det(L) det(U) = l_11 * ... * l_nn * u_11 * ... * u_nn

поэтому вам нужно только определить, является ли какая-либо из записей главной диагонали L или U 0.

Для дальнейшего упрощения используйте алгоритм Дулиттла, где l_ii = 1. Если в какой-то момент алгоритм дает сбой, матрица остается единственной, поэтому вы можете остановиться. Вот суть:

for k := 1, 2, ..., n do {
  for j := k, k+1, ..., n do {
    u_kj := a_kj - sum_{s=1...k-1} l_ks u_sj;
  }
  for i = k+1, k+2, ..., n do {
    l_ik := (a_ik - sum_{s=1...k-1} l_is u_sk)/u_kk;
  }
}

Ключевым моментом является то, что вы можете вычислить ith строку U и ith столбец L одновременно, и вам нужно знать только предыдущую строку/столбец, чтобы двигаться вперед. Таким образом, вы параллельно обрабатываете столько, сколько можете, и сохраняете столько, сколько вам нужно. Поскольку вы можете вычислять записи a_ij по мере необходимости, вам потребуется сохранить два вектора длины n при создании еще двух векторов длины n (строки U, столбцы L). Алгоритм занимает n^2 времени. Возможно, вы сможете найти еще несколько трюков, но это зависит от вашего компромисса между пространством и временем.

person PengOne    schedule 05.03.2013
comment
Похоже, мне нужно хранить обе матрицы LU для их вычисления. У меня не хватает памяти для этого. Есть ли другой способ? - person Daniel; 05.03.2013

arrow_upward
-1
arrow_downward

Не уверен, что следил за вашей проблемой; это (или это сводится) к следующему?

У вас есть два вектора из n чисел, назовите их x и c, тогда элемент матрицы является произведением k на (x_k+c_k), причем каждая строка/столбец соответствует разным порядкам x и c?

Если это так, то я считаю, что матрица будет сингулярной всякий раз, когда в x или c будут повторяющиеся значения, поскольку в этом случае матрица будет иметь повторяющиеся строки/столбцы. Попробуйте несколько вычислений Монте-Карло на меньшем n с различными значениями x и c, чтобы увидеть, является ли этот случай в целом несингулярным - весьма вероятно, что если это верно для 6, то это будет верно и для 10.

Что касается грубой силы, ваш метод:

  1. Не запускается
  2. Будет работать намного быстрее (должно быть несколько секунд для n=7), хотя вместо LU вы можете попробовать SVD, который гораздо лучше сообщит вам, насколько хорошо ведет себя ваша матрица.
person Jon    schedule 24.02.2013
comment
Нет. Один вектор имеет точки в R^(n!), например. (2,2,1,3,5,2), а другой имеет многочлены от n! переменные, например (x1 - 4)(x3 - 5)(x4 - 4)(x6 - 1). Матрица представляет собой оценку многочленов в точках, например. запись для примера будет (2 - 4)(1 - 5)(3 - 4)(2 - 1). - person Daniel; 27.02.2013
comment
Монте-Карло вообще не помогает. У меня есть определенный набор точечных векторов и определенный набор полиномиальных векторов, и мне нужно знать, является ли матрица оценки для этих полиномов в этих точках невырожденной. Вы не можете Монте-Карло, по крайней мере, насколько я знаю. - person Daniel; 27.02.2013