генерация псевдослучайной положительно определенной матрицы
Я хотел протестировать простой код Холецкого, написанный на C ++. Итак, я генерирую случайную нижнюю треугольную букву L и умножаю ее на транспонирование для генерации А.
A = L * Lt;
Но мой код не учитывает фактор А. Поэтому я попробовал это в Matlab:
N=200; L=tril(rand(N, N)); A=L*L'; [lc,p]=chol(A,'lower'); p
Это выводит ненулевое p, что означает, что Matlab также не может учесть A. Я предполагаю, что случайность генерирует матрицы с недостатком ранга. Я прав?
Обновить:
Я забыл упомянуть, что следующий код Matlab работает, как указано Малифом ниже:
N=200; L=rand(N, N); A=L*L'; [lc,p]=chol(A,'lower'); p
Разница в том, что L является первым треугольником в первом коде, а не во втором. Почему это должно иметь значение?
Я также попробовал следующее с scipy после прочтения Простой алгоритм генерации положительно-полуопределенных матриц:
from scipy import random, linalg
A = random.rand(100, 100)
B = A*A.transpose()
linalg.cholesky(B)
Но это ошибки с:
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 66, in cholesky
c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True)
File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 24, in _cholesky
raise LinAlgError("%d-th leading minor not positive definite" % info)
numpy.linalg.linalg.LinAlgError: 2-th leading minor not positive definite
Я не понимаю, почему это происходит со Сципи. Есть идеи?
Спасибо,
Nilesh.
Решение
Проблема не в холесской факторизации. Проблема со случайной матрицей L,
rand(N,N) гораздо лучше обусловлен, чем tril(rand(N,N)), Чтобы увидеть это, сравните cond(rand(N,N)) в cond(tril(rand(N,N))), Я получил что-то вроде 1e3 для первого и 1e19 для второго, так что число кондиционирования второй матрицы намного выше, и вычисления будут менее стабильны в численном отношении.
Это приведет к получению некоторых небольших отрицательных собственных значений в плохо обусловленном случае — чтобы увидеть этот взгляд на собственные значения, используя eig()некоторые маленькие будут отрицательными.
Так что я бы предложил использовать rand(N,N) генерировать численно устойчивую случайную матрицу.
Кстати, если вас интересует теория, почему это происходит, вы можете посмотреть на эту статью:
Другие решения
Как было сказано ранее, собственные значения треугольной матрицы лежат на диагонали. Следовательно, делая
L=tril(rand(n))
Вы убедились, что eig (L) дает только положительные значения. Вы можете улучшить число условий L * L ‘, добавив к диагонали достаточно большое положительное число, например
L=L+n*eye(n)
и L * L ‘положительно определен и хорошо обусловлен:
> cond(L*L')
ans =
1.8400
Чтобы сгенерировать случайную положительно определенную матрицу в MATLAB, ваш код должен выглядеть так:
N=200;
L=rand(N, N);
A=L*transpose(L);
[lc,p]=chol(A,'lower');
eig(A)
p
И вы должны действительно иметь собственные значения больше нуля и p быть ноль.
Вы спрашиваете о нижнем треугольном регистре. Посмотрим, что происходит и почему возникают проблемы. Это часто хорошая вещь, чтобы посмотреть на тестовый пример.
Для простой матрицы 5×5
L = tril(rand(5))
L =
0.72194 0 0 0 0
0.027804 0.78422 0 0 0
0.26607 0.097189 0.77554 0 0
0.96157 0.71437 0.98738 0.66828 0
0.024571 0.046486 0.94515 0.38009 0.087634
eig(L)
ans =
0.087634
0.66828
0.77554
0.78422
0.72194
Конечно, собственные значения треугольной матрицы — это просто диагональные элементы. Поскольку элементы, сгенерированные rand, всегда находятся между 0 и 1, в среднем они будут примерно 1/2. Возможно, рассмотрение распределения определителя L поможет. Лучше рассмотреть распределение log (det (L)). Так как определитель будет просто произведением диагональных элементов, log является суммой логарифмов диагональных элементов. (Да, я знаю, что детерминант — плохая мера сингулярности, но распределение log (det (L)) легко вычисляется, и мне лень думать о распределении числа условия.)
Ах, но отрицательный логарифм равномерной случайной величины является экспоненциальной переменной, в данном случае экспоненциальной с лямбда = 1. Сумма логарифмов набора из n равномерных случайных чисел из интервала (0,1) будет Центральная предельная теорема Гаусса. Среднее значение этой суммы будет -n. Следовательно, определитель нижней треугольной матрицы nxn, порожденной такой схемой, будет exp (-n). Когда n равно 200, MATLAB говорит мне, что
exp(-200)
ans =
1.3839e-87
Таким образом, для матрицы любого заметного размера мы можем видеть, что она будет плохо обусловлена. Хуже того, когда вы формируете произведение L * L ‘, оно обычно будет численно единичным. Те же аргументы применяются к условному номеру. Таким образом, даже для матрицы 20×20 видно, что число условий такой нижней треугольной матрицы довольно велико. Тогда, когда мы сформируем матрицу L * L ‘, условие будет возведено в квадрат, как и ожидалось.
L = tril(rand(20));
cond(L)
ans =
1.9066e+07
cond(L*L')
ans =
3.6325e+14
Посмотрите, насколько лучше вещи для полной матрицы.
A = rand(20);
cond(A)
ans =
253.74
cond(A*A')
ans =
64384