метод наименьших квадратов с ограничением

У меня есть 37 линейных уравнений с 36 переменными в виде матрицы: A x = b. (A имеет 37 строк и 36 столбцов.) У уравнений нет точного решения, поэтому я использовал Matlab чтобы найти ближайший ответ с помощью x = A \ b.

Проблема в том, что у меня также есть условие, все элементы x должны быть положительными: xi > 0 для всех i. x = A \ b дает отрицательные значения для некоторых элементов. Как я могу применить это ограничение?

Вот конкретные значения A и b, с которыми я работаю:

A = [0.83   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.02   0.63    0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.02    -0.2    0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.15   0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0.15    0   0   0   0   0.02    -0.33   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.018   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.17
     1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1];
b = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]';

matlab math least-squares curve-fitting linear-programming
person nasim    schedule 20.05.2015    source источник
comment
Если у вас есть 37 уравнений для 36 неизвестных, это переопределенная система, а это означает, что у нее вряд ли будет несколько решений (и, на самом деле, очень вероятно, что решений не будет).   -  person    schedule 20.05.2015
comment
Я знаю, но нельзя ли хотя бы найти примерное решение с наименьшей ошибкой?   -  person nasim    schedule 20.05.2015
comment
Если вы хотите найти X0, который минимизирует некоторую ошибку, сначала вам нужно определить эту ошибку. :-) Какой тип нормы вы бы использовали для ошибки E = B - A*X0? Евклидово, такси, максимум?   -  person    schedule 20.05.2015
comment
Я даже не знаю этих терминов! (Я читаю о них прямо сейчас) Мои переменные — это вероятности перехода в сети массового обслуживания и обычно находятся в диапазоне от 0,0001 до 0,1. Какой из них вы рекомендуете?   -  person nasim    schedule 20.05.2015
comment
@CST-Link, в заголовке написано метод наименьших квадратов, это определяет функцию стоимости   -  person A. Donda    schedule 20.05.2015
comment
@nasim, можешь опубликовать точные значения A и b? Это упростит отладку возможного ответа. А пока посмотрите мое решение.   -  person A. Donda    schedule 21.05.2015
comment
А не влезет в комментарии :(   -  person nasim    schedule 21.05.2015
comment
@nasim, нет, не комментарии. Отредактируйте вопрос и включите его как код.   -  person A. Donda    schedule 21.05.2015
comment
@nasim, спасибо, но это не имеет смысла. Это A 36x36, а это b 1x38. Пожалуйста, исправьте это.   -  person A. Donda    schedule 21.05.2015

Ответы (3)


arrow_upward
3
arrow_downward

Вы хотите найти приближенное решение x для A x = b по методу наименьших квадратов, т. е. вы хотите минимизировать

||A x - b||2 = xT AT A x + bT b - 2 хT AT б.

Без учета постоянного члена bTb и деления на коэффициент 2 это соответствует форме задача квадратичного программирования, которая заключается в минимизации

1/2 хT Н х + fT х,

если мы выберем H = AT A и f = - AT b.

Соответствующее использование quadprog:

H = A' * A;
f = - A' * b;
x = quadprog(H, f)

Вы также хотите, чтобы элементы x были положительными. Ограничение неотрицательности можно ввести с помощью дополнительных параметров к quadprog, A и b (не путать с вашими матрицами!):

n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

Ограничение положительности не имеет смысла, потому что, если оптимальное решение предполагает, что один или несколько элементов x точно равны 0, то строго положительное решение будет тем лучше, чем меньше соответствующие элементы: 0,01 будет лучше, чем 0,1, 0,001 будет лучше, чем 0,1. быть лучше 0,01 и т. д. и т. д. – естественной границы нет. Если вы хотите убедиться, что решение является полностью положительным, вы должны сами установить конечную границу:

x = quadprog(H, f, -eye(n), zeros(n, 1) + 0.001)

Теперь наименьшее возможное значение элемента x равно 0,001.

Обновление после того, как вопрос был дополнен фактическими данными А и Б: Использование кода

H = A' * A;
f = - A' * b;
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

Я получаю результат:

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

x =
      0.000380906335150292
      3.90638261088393e-05
        0.0111196970167585
        0.0227055107206744
        0.0318402514628274
        0.0371743514880516
      0.000800900221354844
       0.00746652476710186
        0.0180511534370576
        0.0282423767946842
        0.0362606972021829
        0.0417582260990786
       0.00860220929402253
        0.0174105435824309
        0.0265771677458008
        0.0343071472371469
        0.0395176470725881
        0.0419494410289298
        0.0187719294637544
        0.0268976053211278
        0.0336818044612046
        0.0382365751296441
        0.0398823076542831
        0.0391016682549663
        0.0279383031707377
        0.0339393563379992
        0.0377917413001034
        0.0382731422972829
        0.0338557405807941
        0.0217568643500703
        0.0343698083354502
        0.0381554349806972
        0.0392353941260779
        0.0368010570888738
         0.031271868401718
        0.0258232230013864
person A. Donda    schedule 20.05.2015
comment
Предупреждение: Алгоритм Trust-Reflective Region не решает проблемы такого типа, используя алгоритм active-set. Для получения дополнительной справки см. Выбор алгоритма в документации. - person nasim; 21.05.2015
comment
@nasim, я проверил этот код с некоторыми случайно сгенерированными данными, поэтому я не могу сам проверить это предупреждение. Как я комментировал вопрос, лучше бы вы свои матрицы выложили. - person A. Donda; 21.05.2015
comment
Буду рад, но матрица А не влезает в комментарии. Где я могу опубликовать? - person nasim; 21.05.2015
comment
@nasim, щелкните ссылку редактирования под своим вопросом и скопируйте данные в вопрос. - person A. Donda; 21.05.2015
comment
@nasim, спасибо за исправленные данные. Но с этими значениями я получаю правильное решение без предупреждения. Решение полностью положительно даже без указания конечной границы; наименьший x равен 3,9e-5. - person A. Donda; 21.05.2015
comment
Действительно?! Могу ли я увидеть 36 значений, пожалуйста? И могу ли я получить точный фрагмент кода, который вы использовали в качестве ответа? - person nasim; 21.05.2015
comment
Давайте продолжим обсуждение в чате. - person nasim; 21.05.2015

arrow_upward
1
arrow_downward

Это будет проблема линейной оптимизации, когда у вас есть условие x> 0. Лучший алгоритм для решения этого - симплексный алгоритм. Идея состоит в том, что каждое линейное уравнение aix=bi представляет собой линию, а комбинация этих линий дает многоугольник/многогранник, а ответом является одна из вершин этого многогранника/многоугольника. Алгоритм симплекса довольно стандартный, и есть много доступных функций и библиотек, которые могут его вычислить.

person Good Luck    schedule 20.05.2015

arrow_upward
0
arrow_downward

Функции Matlab lsqlin или lsqnonneg можно использовать для решения вашей проблемы. например:

x=lsqnonneg(A,b)

даст вам то, что вы ищете.

person BayesianMonk    schedule 29.08.2018