Использование Eigen для решения плотного, ограниченного наименьшего квадрата
Мне нужно решить классическую задачу вида Ax = b для вектора x размером 4. A имеет порядок ~ 500 точек данных и, следовательно, является плотной матрицей 500×4.
В настоящее время я могу решить эту проблему с помощью нормальных уравнений, описанных Вот и это прекрасно работает, однако я хотел бы ограничить один из моих параметров в х, чтобы никогда не быть выше определенного значения.
Есть ли хороший способ сделать это программно с Eigen?
Решение
Вы можете попробовать мой решатель квадратического программирования на основе Eigen там. Вам все равно придется составить нормальное уравнение.
Другие решения
Вот некоторые питон-демо (используя NumPy, который не так далеко от Eigen), показывая ускоренный алгоритм проецируемого градиента для твоего рода проблемы. Обычно этот подход используется для крупномасштабных задач (где другие алгоритмы, включающие информацию второго порядка, могут не справиться), но его также приятно реализовать.
Эта версия, которая является небольшой модификацией моего старого кода, не является самым простым подходом, который можно использовать, поскольку мы используем:
- ускорение / импульс (более быстрая итерация)
- поиск строки (избавляет нас от проблемы настройки размера шага)
Вы можете удалить поиск строки и настроить размер шага. Импульс тоже не нужен.
Поскольку я сейчас не очень много занимаюсь C ++, я не думаю, что перенесу это на Eigen. Но я уверен, что если вы захотите перенести его, это не так сложно. Eigen не должен слишком отличаться от NumPy.
Я не измерял производительность, но результаты рассчитываются мгновенно (перцептивно).
Редактировать: некоторые ненаучные моменты времени (больший импульс; меньший допуск, чем в следующем коде):
A=(500,4):
Solve non-bounded with scipys lsq_linear
used: 0.004898870004675975
cost : 244.58267993
Solve bounded (0.001, 0.05) with scipys lsq_linear
used: 0.005605718416479959
cost : 246.990611225
Solve bounded (0.001, 0.05) with accelerated projected gradient
early-stopping @ it: 3
used: 0.002282825315435914
cost: 246.990611225
A=(50000, 500):
Solve non-bounded with scipys lsq_linear
used: 4.118898701951786 secs
cost : 24843.7115776
Solve bounded (0.001, 0.05) with scipys lsq_linear
used: 14.727660030288007 secs
cost : 25025.0328661
Solve bounded (0.001, 0.05) with accelerated projected gradient
early-stopping @ it: 14
used: 5.319953458329618 secs
cost: 25025.0330754
Основная идея состоит в том, чтобы использовать некоторый алгоритм, подобный градиентному спуску, и проецировать на наши ограничения после каждого шага градиента. Этот подход очень эффективен, если этот прогноз может быть выполнен эффективно. Коробка-ограничение-проекции просты!
Страница 4 в этом PDF показывает проекцию рамочного ограничения.
Мы просто обрезаем наш вектор-решение на lower_bound, upper_bound, Отсечение в numpy описывается как: При заданном интервале значения за пределами интервала обрезаются до краев интервала. Например, если задан интервал [0, 1], значения меньше 0 становятся 0, а значения больше 1 становятся 1.
Это итеративный алгоритм, приближающий решение, и я думаю, что каждый используемый алгоритм будет повторяющимся.
Код
import numpy as np
from scipy.optimize import lsq_linear
np.random.seed(1)
A = np.random.normal(size=(500, 4))
b = np.random.normal(size=500)
""" Solve Ax=b
----------
"""print('Solve non-bounded with scipys lsq_linear')
sol = lsq_linear(A, b)
print('Ax=b sol: ', sol['x'])
print('cost : ', sol['cost'])
print()
""" Solve Ax=b with box-constraints
-------------------------------
"""print('Solve bounded (0.001, 0.05) with scipys lsq_linear')
sol = lsq_linear(A, b, bounds=(0.001, 0.05))
print('Ax=b constrained sol: ', sol['x'])
print('cost : ', sol['cost'])
print()
""" Solve Ax=b with box-constraints using a projected gradient algorithm
--------------------------------------------------------------------
"""def solve_pg(A, b, bounds=(-np.inf, np.inf), momentum=0.9, maxiter=1000):
""" remarks:
algorithm: accelerated projected gradient
projection: proj onto box constraints
line-search: armijo-rule along projection-arc (Bertsekas book)
stopping-criterion: naive
gradient-calculation: precomputes AtA
"""
lb = np.empty(A.shape[1])
ub = np.empty(A.shape[1])
if len(bounds) == 2:
# apply lb & ub to all variables
lb = bounds[0]
ub = bounds[1]
else:
# assume dimensions are ok
lb = np.array(bounds[0])
ub = np.array(bounds[1])
M, N = A.shape
x = np.zeros(N)
AtA = A.T.dot(A)
Atb = A.T.dot(b)
stop_count = 0
def gradient(x):
return AtA.dot(x) - Atb
def obj(x):
return 0.5 * np.linalg.norm(A.dot(x) - b)**2
it = 0
while True:
grad = gradient(x)
# line search
alpha = 1
beta = 0.5
sigma=1e-2
old_obj = obj(x)
while True:
new_x = x - alpha * grad
new_obj = obj(new_x)
if old_obj - new_obj >= sigma * grad.dot(x - new_x):
break
else:
alpha *= beta
x_old = x[:]
x = x - alpha*grad
# projection
np.clip(x, lb, ub, out=x) # Projection onto box constraints
# see SO-text
# in-place clipping
y = x + momentum * (x - x_old)
if np.abs(old_obj - obj(x)) < 1e-2:
stop_count += 1
else:
stop_count = 0
if stop_count == 3:
print('early-stopping @ it: ', it)
return x
it += 1
if it == maxiter:
return x
print('Solve bounded (0.001, 0.05) with accelerated projected gradient')
sol = solve_pg(A, b, bounds=(0.001, 0.05))
print(sol)
print('cost: ', 0.5 * (np.square(np.linalg.norm(A.dot(sol) - b))))
Выход
Solve non-bounded with scipys lsq_linear
Ax=b sol: [ 0.06627173 -0.06104991 -0.07010355 0.04024075]
cost : 244.58267993
Solve bounded (0.001, 0.05) with scipys lsq_linear
Ax=b constrained sol: [ 0.05 0.001 0.001 0.03902291]
cost : 246.990611225
Solve bounded (0.001, 0.05) with accelerated projected gradient
early-stopping @ it: 3
[ 0.05 0.001 0.001 0.03902229]
cost: 246.990611225