Генерация случайных неособых целочисленных матриц

Question

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

a _{i, j} (i, j: 1..n) / & forall; (i, j) a _{i, j} & isin; ℤ и 0 & le; a _{i, j} & le; k и Det [a] & ne; 0

но a _{i, j} также должны быть случайными после равномерного распределения в [0, k].

В своем нынешнем воплощении проблема имеет n & cong; 300, k & cong; 100.

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

Проблема в том, что эта проверка для матриц 300x300 занимает что-то около 2 секунд, и я не могу себе этого позволить.

Конечно, я могу построить строки, выбрав случайную первую строку, а затем создав последовательные ортогональные строки, но я не уверен, как гарантировать, что эти строки будут иметь свои элементы после равномерного распределения в [0, k].

Я ищу решение в Mathematica, но также приветствуется более быстрый алгоритм генерации матриц.

NB> Условие U [0, k] означает, что для набора матриц каждая позиция (i, j) в наборе должна следовать равномерному распределению.

Answer 1

Вместо этого вы можете использовать MatrixRank.На моей машине это примерно в n / 10 раз быстрее для больших nxn целочисленных матриц.

Answer 2

Как в Matlab, так и в Octave детерминант и факторизация LU матрицы 500x500 в основном происходят мгновенно.Есть ли у Mathematica режим, в котором он может вызывать LAPACK или какую-то похожую библиотеку?Возможно, вам придется аннотировать, что ваши массивы должны рассматриваться как числа с плавающей запятой, а не символически;это может сделать это намного быстрее.Для сравнения: LU на матрице 5000x5000 занимает 8,66 секунды в моей системе с использованием Octave;500x500 должно быть примерно в 1000 раз быстрее.

Answer 3

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

k = 100; n = 500;
mat = RandomInteger[100, {n, n}];

AbsoluteTiming[Det[mat] == 0]

Out [57] = {6.8640000, False}

AbsoluteTiming[Det[N@mat] == 0.] (*warning light!!*)

Out [58] = {0.0230000, False}

AbsoluteTiming[MatrixRank[N@mat] != n]

Out [59] = {0.1710000, False}

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

mat2 = mat;
mat2[[-1]] = Total[Most[mat]];

AbsoluteTiming[Det[mat2] == 0]

Out [70] = {9.4750000, True}

AbsoluteTiming[Det[N@mat2] == 0.]

Out [69] = {0.0470000, False}

AbsoluteTiming[MatrixRank[N@mat2] != n]

Out [71] = {0.1440000, True}

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

Продолжая приведенные выше примеры:

AbsoluteTiming[Det[mat, Modulus -> Prime[1000]]]

Out [77] = {0,6320000, 4054}

AbsoluteTiming[Det[mat2, Modulus -> Prime[1000]]]

Out [78] = {0,6470000, 0}

Это медленно, но быстрее, чем работа над рациональными. Что бы это ни стоило, в большинстве случаев я был бы достаточно уверен в результатах более быстрого тестирования через MatrixRank [N [matrix]].

Даниэль Лихтблау Wolfram Research

Answer 4

Вот расширение комментария, который я сделал немного.Я согласен с Дэном в том, что очень маловероятно, что числовая версия выдаст ложное срабатывание.Тем не менее, вы можете избежать этого сценария, изучив единичные значения и вернув False надежно, если наименьшее единичное значение больше некоторой погрешности.(Надо полагать, найти доказуемый допуск может быть немного сложно.) Если наименьшее единственное значение неудобно мало, вы можете применить Det к целочисленной матрице.

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

singularQ[M_?MatrixQ] := If[
  Last[SingularValueList[N[M], Tolerance -> 0]] > 1/10.0^8,
 False, Det[M] == 0];

Вот 200 матриц, которые соответствуют вашему описанию.Одна посередине была сфальсифицирована, чтобы быть единственной.

SeedRandom[1];
matrices = RandomInteger[{0, 100}, {200, 300, 300}];
matrices[[100, -1]] = Sum[RandomInteger[{0, 10}]*matrices[[100, k]],
 {k, 1, Length[matrices[[100]]] - 1}];

Теперь давайте найдем индексы всех сингулярных матриц, наблюдая, как мы идем.

Flatten@Monitor[Table[
  If[singularQ[matrices[[k]]], k, {}],
  {k, 1, Length[matrices]}], k]