Проверка, имеют ли две матрицы 12x12 одинаковый определитель

Мне дают матрицу Q которая является симметричной, обратимой, положительно определенной и плотной. Мне нужно проверить, если det, где J является матрицей всех единиц.$12 \times 12$$Q$

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

1. Есть ли какой-нибудь прием, который я мог бы использовать, чтобы вычислить определитель быстрее, учитывая, что я знаю их размер? Может быть, грязное расширение для матриц $12 \times12$ которое могло бы работать в этом случае?

2. Есть ли другой эффективный способ проверить равенство $(1)$

Редактировать. Чтобы ответить на комментарии. Мне нужно вычислить все связные не самодополняющие графы $G$ порядка $13$ , чтобы $G$ и $\overline{G}$ имели одинаковое количество остовных деревьев. Мотивация для этого может быть найдена в этом сообщении mathoverflow . Что касается машины, я запускаю ее на 8-ядерном 3.4GHh компьютере параллельно.

Редактировать. Мне удалось сократить ожидаемое время выполнения на 50%, создав программу на C для специального вычисления определителя матрицы . Предложения все еще приветствуются.$12 \times 12$

Так как вы уже используете C ++ и ваши матрицы симметричная положительно определенная, я бы выполнить unpivoted факторизация Q , а также из 12 I - Q - J . Здесь я предполагаю, что 12 I - Q - J также является положительно определенным, в противном случае L D L T потребует поворота для численной устойчивости (также возможно, что, хотя оно не является положительно определенным, поворот не требуется, но вы должны попробуй).$LDL^T$$Q$$12I-Q-J$$12I-Q-J$$LDL^T$

Это быстрее, чем факторизация LU, а также быстрее, чем Cholesky, потому что избегают квадратных корней. Определитель - это просто произведение элементов диагональной матрицы. Код для выполнения факторизации ЛПНП настолько прост, что вы можете написать его менее чем в 50 строках C. На странице Википедии описан алгоритм, и у меня есть несколько простых шаблонных кодов, которые можно сделать здесь по Холецки . Вы можете значительно упростить это и изменить его, чтобы избежать квадратного корня для реализации факторизации L D L $D$$LDL^T$

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

Я бы не пытался разыгрывать трюки, чтобы воспользоваться тем фактом, что содержит Q внутри выражения. Вполне вероятно, что при таких небольших размерах задач эти приемы оказываются медленнее, чем просто выполнение двух отдельных определителей. Конечно, единственный способ проверить эти претензии - это попробовать.$12I-Q-J$$Q$

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

Я скачал пакет Armadillo с Sourceforge и посмотрел документацию. Попробуйте улучшить производительность отдельно вычисления и йе ( 12 I - Q - J ) , где J является ранг одна матрица из всех из них, например , путем установки . В документации отмечается, что это значение по умолчанию для матриц размером до 4 × 4 , поэтому, по пропущению, я полагаю, что это вариант по умолчанию для случая 12 × 12 .$\det(Q)$$\det(12I - Q - J)$$J$det(Q,slow=false)$4\times 4$slow=true$12\times 12$

Что slow=true предположительно делает частичное или полное поворотное в получении формы эшелонированной строки, из которых определитель легко найден. Однако вы заранее знаете, что матрица положительно определена, поэтому поворот не нужен для стабильности (по крайней мере, предположительно для большей части ваших вычислений. Неясно, выбрасывает ли пакет Armadillo исключение, если стержни становятся чрезмерно маленькими, но это должно быть особенность разумного пакета числовой линейной алгебры. РЕДАКТИРОВАТЬ: Я нашел код Armadillo, который реализует в заголовочном файле , используя шаблоны C ++ для существенной функциональности. Похоже, что настройка не влияет на 12 × $Q$detinclude\armadillo_bits\auxlib_meat.hppslow=false$12\times 12$определитель будет выполнен, потому что вычисление «отбрасывается через стену» в LAPACK (или ATLAS) в этой точке без указания того, что поворот не требуется; увидеть det_lapackи его вызовы в этом файле.

Другой момент - следовать их рекомендации по созданию пакета Armadillo, связанного с высокоскоростными заменами для BLAS и LAPACK, если вы действительно их используете; см. гл. 5 файла Armadillo README.TXT для деталей. [Использование выделенной 64-битной версии BLAS или LAPACK также рекомендуется для скорости на современных 64-битных машинах.]

Приведение строки к форме эшелона является по существу гауссовым исключением и имеет арифметическую сложность . Для обеих матриц это в два раза больше, или$\frac{2}{3} n^3 + O(n^2)$. Эти операции вполне могут быть «узким местом» в вашей обработке, но мало надежды на то, что без специальной структуры вQ(или каких-либо известных связей между триллионом тестовых случаев, допускающих амортизацию) работа может быть уменьшена доO(n2).$\frac{4}{3} n^3 + O(n^2)$$Q$$O(n^2)$

Для сравнения, разложение по кофакторам общей матрицы включает n ! операции умножения (и примерно столько же сложений / вычитаний), поэтому при n = 12 сравнение ( 12 ! = 479001600 против $n\times n$$n!$$n=12$$12! = 479001600$) явно способствует исключению по сравнению с кофакторами.$\frac{2}{3} n^3 = 1152$

Другой подход, требующий Работа 3 n3+O(n2)сводитQк трехдиагональной форме с преобразованиями Хаусхолдера, что также переводит12I-Qв трехдиагональную форму. ВычислениеDET(Q)иDET(12I-Q-J)затем может быть сделано вO(п)операций. [Эффект обновления первого ранга-$\frac{4}{3} n^3 + O(n^2)$$Q$$12I - Q$$\det(Q)$$\det(12I - Q - J)$$O(n)$$-J$ во втором определителе можно выразить как скалярный фактор, заданный решением одной трехугольной системы.]

Реализация такого независимого вычисления может быть полезна как проверка результатов успешных (или неудачных) вызовов detфункции Armadillo .

Особый случай: Как следует из Комментария Джерней, предположим, что где J, как и прежде, является матрицей (ранг 1) всех единиц, а D = diag ( d 1 , … , d n ) является неособой (положительной) ) диагональная матрица. Действительно, для предложенного приложения в теории графов это будут целочисленные матрицы. Тогда явная формула для det ( Q ) :$Q = D - J$$J$$D=\text{diag}(d_1,\ldots,d_n)$$\det(Q)$

Схема его доказательства дает возможность проиллюстрировать более широкую применимость, т. Е. Всякий раз, когда имеет известный определитель и система D v = ( 1 ... 1 ) T быстро решается. Начните с факторинга:$D$$Dv = (1\ldots 1)^T$

Теперь снова имеет ранг 1, а именно ( d - 1 1 ... d - 1 n ) T ( 1 ... 1 ) . Обратите внимание, что второй определитель просто:$D^{-1}J$$(d_1^{-1}\ldots d_n^{-1})^T(1\ldots 1)$

где - характеристический многочлен D - 1 Дж . Как матрица ранга 1, f ( x ) должна иметь (как минимум) n - 1 множителей x для учета ее нулевого пространства. «Пропущенное» собственное значение равно ∑ d - 1 i , как видно из вычисления:$f(x)$$D^{-1}J$$f(x)$$n-1$$x$$\sum d_i^{-1}$

Отсюда следует , что характеристический полином , и F ( 1 ) является , как показано выше для DET ( я - Д - 1 Дж ) , 1 - Σ d - 1 я .$f(x) = x^{n-1} (x - \sum d_i^{-1})$$f(1)$$\det(I - D^{-1}J)$$1 - \sum d_i^{-1}$

Также отметим, что если , то 12 I - Q - J = 12 I - D + J - J = 12 I - D , диагональная матрица, определитель которой является просто произведением ее диагональных элементов.$Q = D-J$$12I - Q - J = 12I - D + J - J = 12I - D$

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

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

То есть для матрицы $A$ и векторов $u,v$ :