Войти
Физика для игрСтатьи

Решение ЛЦП методом Лемке. Прочитано на gamedev_lecture. (2 стр)

Автор:

[21:29] <fd> Теперь  о представлении
[21:30] <fd> когда эта процедура делается на бумажке, уместно выражать базисные перменные через небазисные
[21:31] <fd> Тогда вместо w - Ax - ez = q мы получим эквивалентную систему, в которой под базисными перменными стоят столбцы единчиной матрицы, а жругие столбци и свободный вектор, соответственно, как-то изменены
[21:32] <fd> Т.е. сисетема w - Ax - ez = q решена относительно базисных перменных
[21:33] <fd> Это удобно по двум причинам. Во-первых, проскольку небазисные пременные у нас 0, transformed_q содержит значения базисных перменных
[21:35] <fd> Во-вторых, увеличиваяя любую из небазисных перменных от 0 ввверх, мы быстро увидим, когда же её рост уппрётся в одно их услвоий неотрицательности - это когда сумма её столбца, назовём его c, в текущем базисе (когда базисные - разрешающие), 'то будет transformed_c, и transformed_q получит отрицательный элемент
[21:37] <fd> Альтернативынй способ узнать то же самое - это исследовать вектор частных transformed_q / transformed_c
[21:38] <fd> мимнимальный элемент этого вектора будет иметь индекс i, это будет как раз индекс блокирующей перменной
[21:39] <fd> Итак, зная текщий ьазис, мы решаем уравнение w - Ax -ez = q относительно базисных переменных, раз
[21:40] <fd> зная входящую перменную (как дополнительную к ушедшей на прошлой итерации, или указанным выше способом дял 1й итерации), мы берём её столбец в трансформированной системе уравнений (решённйо относительно базисныйх перменных), назовём его transformed__c
[21:41] <fd> Берём правую часть из трансформированнйо системы, transformed_q
[21:42] <fd> ищем минимальный элемент среди частных transformed_q / Transfromed_C, и значит перменная i из базиса уходит
[21:43] <fd> Новый базис не содержит больше ушедшей перменной, но содержит пришедшую. Итерация выполена.
[21:44] <fd> Прежде чем покажу код ФайндДепартинг, пару слов
[21:44] <fd> первой слово - лексико-минимум тест
[21:44] <fd> среди частных transformed_q
/ Transfromed_C могут быть равные
[21:45] <fd> тогда нужно из них однозначным образом выбрать один индекс, чтоб исключить зацикливанеи
[21:45] <fd> метод это сделать - лексико минимум тест
[21:48] <fd> суть его в том, что мы находим вектор частных transformed_M[a]
/ transformed_c, т.е. a-го столбца трансформированной системы к столбцу входящей перменной
[21:48] <fd> и ищем минимальное значение среди частных с теми индексами, которые остались после первой проверки
[21:50] <fd> То есть при обычнйо проверки было у нас 3 равных частынх  transformed_q
/ Transfromed_C, значит, для этих 3 i мы считаем  transformed_M[a] / transformed_c и выбираем тот i, где уже такое частное меньше
[21:51] <fd> Я немного завтарлся в обознаечниях.  transformed_M[a]
у меня - iй элемент a-го столбуа именно
[21:51] <fd> Если после такой проверки ещё осталась двусмысленность, берём следующий столбец
[21:52] <fd> т.е. а = 1 ... n
[21:52] <fd> В результате всегда делаем однозначный выбор
[21:52] <fd> Второе, что я хотел сказать
[21:53] <fd> Вообще говоря, на бумажке можно решать систему относительно базисных перменных и в ус н едуть
[21:53] <fd> в компе удобно пермтавлять переменные в базисе физически
[21:54] <fd> Поясню. Так у нас система (I -A -e) (w x z) = q, не трансформированная
[21:56] <fd> сы знаем, какие из перменных (w x z) базисные, какие нет (на каждой итерации, кончено). Мы можем их переставить, так что если b - базисные перменные, n - небазисные, B - столбцы базисных, N - столбцы небазисных, то система будет (B N) * (b n) = q
[21:57] <fd> Тогда, найдя значение B^-1 и домножив матрицу с перестановенными столбцами на на неё, получим (I transformed_N)(b n) = transformed_q
[21:57] <fd> такой подход удобнее..
[21:58] <fd> тогда transformed_C - это столбец из transformed_N, соответсвующий входящей перменной
[21:59] <fd> а матрица tramsformed_M - это как раз B^-1
[21:59] <fd> так.. даю код
[22:00] <fd> http://www.everfall.com/paste/id.php?gl48e69gpe30
[22:01] <fd> строки 5 - 26 - это рассчёт tramsformed_C
[22:02] <fd> это умножение B^-1 на C
[22:03] <fd> Если C - столбец перменной типа Х (y - в коде), то это столбец матрицы A, умноженый на -1, а если столбец переменной w, то это один их столбцов единичной матрицы
[22:03] <fd> перменная, конечно, входящая имеется ввиду
[22:04] <fd> строки 27-33 - это мы праву сторону приводим к текущему базису
[22:05] <fd> 33-60 - эт опервичный поиск минимума среди частных Transformed_C
/ Transformed_Q
[22:10] <fd> вот курю щас в свой код, и по-моему в этом билде ошибочка есть - не вижу я, чтобы выбриался минимальный ПОЛОЖИТЕЛЬНЫЙ элемент
[22:10] <fd> а он должен
[22:11] <fd> если все отрицательные, то неблокирующий рост и бесконченый луч
[22:13] <fd> так, эт онадо на свежу голову проверить, но вроде бы ошибка. Должен выбираться н просто минимальный, а минимальный положительный элемент среди частных
[22:13] <fd> ладно
[22:13] <fd> то, что ниже в этом коде - эт олексико-минимум тест
[22:14] <fd> для тех пременных, для которых осталась двусмысленность
[22:14] <fd> она удаляется
[22:15] <fd> А, нет, фу ты! Нет ошибочки! 43 строка
[22:15] <fd> и никто н ескажет :(
[22:16] <fd> ладно
[22:19] <fd> Теперь пару слов об обращении матрицы B
[22:19] <fd> это мрачно
[22:19] <fd> :)
[22:20] <fd> есть метод для этого
[22:20] <fd> вычитал у Клайна
[22:21] <fd> в матрице B упорядочиваем столбцы, а в векторе базисынх перменных - перменные, особым образом
[22:21] <fd> а именно, чтобы входящая перменная занимала место уходящей
[22:22] <fd> тогда от базиса к базису в матрице B изменяется один столбец
[22:24] <fd> это позволяет использовать формулу шермана-морисона-вудбери http://mathworld.wolfram.com/Sherman-MorrisonFormula.html
[22:25] <fd> вектором u берётся разница нового и старого столбцов B, а v  - вектор-строка с единственной 1 в позиции i, если мы меняем i-й столбец
[22:27] <fd> в чём тут оптимизация
[22:28] <fd> v * A^-1 * v = это просто i-й элемент A^-1 * v, а A^-1 * v - это B^-1 * c, tramsformed_c той перменной, которую мы ввели, чтобы получить этот базис
[22:31] <fd> Првкка: с разницей в 1.0
[22:32] <fd> дело в том, что v - это не c, а с минус тот стобец, который там был, А этот тсолбец, который быЛ, i-й столбец B, при умножении на B^-1 даст i-й столбец единичной матрицы
[22:33] <fd> поэтому v * A^-1 * v = это просто i-й элемент A^-1 * v - 1.0
[22:33] <fd> вот.
[22:34] <fd> ну а v * A^-1 это просто i-й столбец матрицы B^-1 до изменения
[22:34] <fd> щас будет код
[22:36] <fd> http://www.everfall.com/paste/id.php?xutvw0sdu5w8
[22:36] <fd> обращение матрицы с этими соображениями
[22:37] <fd> Возникает вопрос первоначального базиса, матрицы B^-1 lkz ytuj
[22:37] <fd> для него
[22:37] <fd> вопрос имеет ответ :)
[22:37] <fd> матрица B для начального базиса имеет простой вид
[22:38] <fd> ето едичничная матрица (столбцы перменных w - столбцы единичной матрицы), лишь с одним столбцом минус единиц (столбец искусственнйо перменной)
[22:38] <fd> такая матрица - имподентна, т.е. обартна самой себе
[22:38] <fd> поэтому мы тупо формируем B^-1...
[22:39] <fd> Ну, итераци ипродолжаются либо до облома FindDeparting, либо до того, как выходящая перменная не будет z
[22:40] <fd> если при этом правльно скорректировать базис, то в B^-1 * q у нас окажутся значения базисных перменных, что логично
[22:41] <fd> (I transformed_N)(b n) = transformed_q, n = 0, b = transformed_q
[22:41] <fd> восстанавливаем порядок переменных, наслаждаемся
[22:41] <fd> Соображения сходимости
[22:41] <fd> сходится оно не всегда, но часто :)
[22:42] <fd> иногда не сходитс потому, что нет решения, а иногда потмоу, что не смогло ег оустановить
[22:43] <fd> Если Лемке не сошлось с матрицей А - коположительной матрицей, то задача не имеет решения вообще
[22:46] <fd> условие коположительности матрицы: матрица А коположительна-плюс, если для любого неотрицательного вектора x из x.Traspose() * A * x = 0 следует x.Trasposed() * ( A + A.Trasposed()) = 0
[22:46] <fd> речь шла о кополодительных плюс, конечно
[22:47] <fd> если строго коположительна, то завершается и всегда имеет решение
[22:48] <fd> строго коположительна - это когда для любого неотрицательного вектора x : x.Traspose() * A * x = 0
[22:49] <fd> ещё на L-матрицах сходится всегда, когда естиь решение
[22:49] <fd> на L* - всегда есть решение
[22:51] <fd> про L-матрицы найдите где-нибудь..
[22:51] <fd> что важно - эт оположительно-полуопределённые матрицы
[22:54] <fd> на них должно сходиться, если есть решение
[22:54] <fd> сейчас, правда, я не нашёл доказательства, но была такая пьянка
[22:57] <fd> ладно
[22:57] <fd> условия сходимости преварили
[22:57] <fd> теперь что.. сложность..
[22:57] <fd> она ужасна
[22:57] <fd> от o(x^3) до o(X^4) в худшем случае
[22:57] <fd> это делает алгоритм малоприменимым на практике
[22:58] <fd> теперь пару слов о миксд-лцп
[23:01] <fd> плюньте в глаза Клайну, метод лемке невозможно адаптировать к решению миксд-лцп
[23:01] <fd> дело в том, что придётся явно инверить матриу B^-1 на й итерации, да и потом каждый кадр ворочить большую матрицу

Страницы: 1 2 3 Следующая »

30 марта 2006