Решение переопределенной системы ограничений

Question

У меня есть n переменные действительного числа (не знаю, на самом деле все равно), давайте назовем их X[n].У меня также есть m >> n отношения между ними, назовем их R[m], в форме:

X[i] = alpha*X[j], alpha - ненулевое положительное действительное число, i и j различны, нопара (i, j) не обязательно уникальна (т. е. между одними и теми же переменными может быть два отношения с разным альфа-фактором)

Я пытаюсь найти набор alpha параметров, которые решаютпереопределенная система в некотором смысле наименьших квадратов.Идеальным решением было бы минимизировать возведенную в квадрат сумму разностей между каждым параметром уравнения и его выбранным значением, но меня устраивает следующее приближение:

Если я превращу m уравнений в переопределенную систему из n неизвестныхлюбой псевдообратный числовой решатель даст мне очевидное решение (все нули).Итак, в настоящее время я добавляю еще одно уравнение в микс x[0] = 1 (на самом деле подойдет любая константа) и решает сгенерированную систему в смысле наименьших квадратов, используя псевдообрат Мура-Пенроуза.Хотя это пытается минимизировать сумму (x[0] - 1)^2 и квадратную сумму x[i] - alpha*x[j], я считаю это хорошим и численно устойчивым приближением к моей проблеме.Вот пример:

a = 1
a = 2*b
b = 3*c
a = 5*c

в октаве:

A = [
  1  0  0;
  1 -2  0;
  0  1 -3;
  1  0 -5;
]

B = [1; 0; 0; 0]

C = pinv(A) * B or better yet:
C = pinv(A)(:,1)

, который дает значения для a, b, c: [0.99383; 0.51235; 0.19136], который дает мнеследующие (разумные) отношения:

a = 1.9398*b
b = 2.6774*c
a = 5.1935*c

Так что сейчас мне нужно реализовать это в C / C ++ / Java, и у меня есть следующие вопросы:

Есть ли более быстрый способ решениямоя проблема, или я на правильном пути с созданием переопределенной системы и вычислением псевдообратного?

Мое текущее решение требует разложения по сингулярному значению и трех умножений матриц, что немного с учетом mможет быть 5000 или даже 10000. Существуют ли более быстрые способы вычисления псевдообратного (на самом деле, мне нужен только первый столбец этого столбца, а не вся матрица, учитывая, что B равен нулю, кроме первой строки), учитывая разреженность матрицы(каждая строка содержит ровно два ненулевых значения, одно из которых всегда одно, а другое всегда отрицательное)

Какие математические библиотеки вы бы предложили использоватьза это?Лапак в порядке?

Я также открыт для любых других предложений, при условии, что они численно устойчивы и асимптотически бывают быстрыми (скажем, k*n^2, где k может быть большим).

Answer 1

Ваша проблема некорректна.Если вы рассматриваете проблему как функцию от n переменных (наименьшего квадрата разности), функция имеет ровно ОДНУ глобальную минимальную гиперплоскость.

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

Если вам нужна параматеризация гиперплоскости решения, вы можете получить ее из псевдообратного уравнения Мура-Пенроуза http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse и проверить раздел о получении всех решений.

(Обратите внимание, что я использовал слово «гиперплоскость» технически некорректным образом. Я имею в виду «плоское» неограниченное подмножество вашего пространства параметров ... Линия, плоскость, что-то, что может быть параматеризовано одним или несколькими векторами. Почему-то не могу найти общее существительное для таких объектов)

Answer 2

Подход SVD численно очень стабилен, но не очень быстр.Если вы используете SVD, то LAPACK - хорошая библиотека для использования.Если это просто однократное вычисление, то оно, вероятно, достаточно быстрое.

Если вам нужен существенно более быстрый алгоритм, вам, возможно, придется пожертвовать стабильностью.Одной из возможностей будет использование QR-факторизации.Вы должны будете прочитать это, чтобы увидеть детали, но часть рассуждений выглядит следующим образом.Если AP = QR (где P - матрица перестановок, Q - ортогональная матрица, а R - треугольная матрица) - это экономичное QR-разложение A, тогда уравнение AX = B становится QRP ^ {- 1} X = Bи решение X = PR ^ {- 1} Q ^ T B. Следующий октавный код иллюстрирует это, используя те же A и B, что и в вашем коде.

[Q,R,P] = qr(A,0)
C(P) = R \ (Q' * B)

Хорошая вещь об этом состоит в том, чтоВы можете использовать разреженность A, выполняя разреженную QR-декомпозицию.В справке Octave есть некоторое объяснение функции qr, но она не сработала для меня сразу.

Еще быстрее (но также и менее стабильно) использовать нормальные уравнения: если AX = B, то A ^TAX = A ^ T B. Матрица A ^ TA является квадратной матрицей (надеюсь) полного ранга, поэтому вы можете использовать любой решатель для линейных уравнений.Октавный код:

C = (A' * A) \ (A' * B)

Опять же, в этом подходе можно использовать разреженность.Есть много методов и библиотек для решения разреженных линейных систем;популярным кажется UMFPACK .

Добавлено позже: Я недостаточно знаю об этом поле для количественной оценки.Целые книги были написаны на этом.Возможно, QR примерно в 3 или 5 раз быстрее SVD, а нормальные уравнения в два раза быстрее.Влияние на числовую стабильность зависит от вашей матрицы A. Разреженные алгоритмы могут быть намного быстрее (скажем, с коэффициентом m), но их вычислительные затраты и численная стабильность очень сильно зависят от проблемы способами, которые иногда не совсем понятны.

В вашем случае использования я бы порекомендовал попытаться вычислить решение с помощью SVD, посмотреть, сколько времени это займет, и если это приемлемо, то просто используйте это (я думаю, это будет около минуты для n =1000 и m = 10000).Если вы хотите изучить его дальше, попробуйте также QR и нормальные уравнения и посмотрите, насколько они быстрее и насколько точны;если они дают примерно то же решение, что и SVD, то вы можете быть уверены, что они достаточно точны для ваших целей.Только если все это слишком медленно и вы готовы потратить некоторое время на это, посмотрите на редкие алгоритмы.