0
我在想什麼是Matlab中最快的凸優化器,或者有什麼方法來加速當前求解器?我正在使用CVX,但它永久解決了我的優化問題。 我有優化是解決Matlab中的快速CVX求解器
minimize norm(Ax-b, 2)
subject to
x >= 0
and x d <= delta
其中A的大小,B是非常大的。
有沒有什麼辦法可以通過最小二乘法求解,然後將其轉移到約束版本以使其更快?
我在想什麼是Matlab中最快的凸優化器,或者有什麼方法來加速當前求解器?我正在使用CVX,但它永久解決了我的優化問題。 我有優化是解決Matlab中的快速CVX求解器
minimize norm(Ax-b, 2)
subject to
x >= 0
and x d <= delta
其中A的大小,B是非常大的。
有沒有什麼辦法可以通過最小二乘法求解,然後將其轉移到約束版本以使其更快?
我不知道什麼x.d <= delta
手段,但我就認爲它應該是x <= delta
。
可以使用投影梯度方法或加速投影梯度方法(這是投影梯度方法,其中「神奇」收斂速度更快的只是輕微的修改)解決了這個問題。這裏有一些python代碼展示瞭如何最小化.5 ||斧 - B ||^2受約束0 < = X < =使用FISTA,這是一個加速的投影梯度方法增量。關於投影梯度方法和FISTA的更多細節可以在例如在近端算法Boyd的manuscript找到。
import numpy as np
import matplotlib.pyplot as plt
def fista(gradf,proxg,evalf,evalg,x0,params):
# This code does FISTA with line search
maxIter = params['maxIter']
t = params['stepSize'] # Initial step size
showTrigger = params['showTrigger']
increaseFactor = 1.25
decreaseFactor = .5
costs = np.zeros((maxIter,1))
xkm1 = np.copy(x0)
vkm1 = np.copy(x0)
for k in np.arange(1,maxIter+1,dtype = np.double):
costs[k-1] = evalf(xkm1) + evalg(xkm1)
if k % showTrigger == 0:
print "Iteration: " + str(k) + " cost: " + str(costs[k-1])
t = increaseFactor*t
acceptFlag = False
while acceptFlag == False:
if k == 1:
theta = 1
else:
a = tkm1
b = t*(thetakm1**2)
c = -t*(thetakm1**2)
theta = (-b + np.sqrt(b**2 - 4*a*c))/(2*a)
y = (1 - theta)*xkm1 + theta*vkm1
(gradf_y,fy) = gradf(y)
x = proxg(y - t*gradf_y,t)
fx = evalf(x)
if fx <= fy + np.vdot(gradf_y,x - y) + (.5/t)*np.sum((x - y)**2):
acceptFlag = True
else:
t = decreaseFactor*t
tkm1 = t
thetakm1 = theta
vkm1 = xkm1 + (1/theta)*(x - xkm1)
xkm1 = x
return (xkm1,costs)
if __name__ == '__main__':
delta = 5.0
numRows = 300
numCols = 50
A = np.random.randn(numRows,numCols)
ATrans = np.transpose(A)
xTrue = delta*np.random.rand(numCols,1)
b = np.dot(A,xTrue)
noise = .1*np.random.randn(numRows,1)
b = b + noise
def evalf(x):
AxMinusb = np.dot(A, x) - b
val = .5 * np.sum(AxMinusb ** 2)
return val
def gradf(x):
AxMinusb = np.dot(A, x) - b
grad = np.dot(ATrans, AxMinusb)
val = .5 * np.sum(AxMinusb ** 2)
return (grad, val)
def evalg(x):
return 0.0
def proxg(x,t):
return np.maximum(np.minimum(x,delta),0.0)
x0 = np.zeros((numCols,1))
params = {'maxIter': 500, 'stepSize': 1.0, 'showTrigger': 5}
(x,costs) = fista(gradf,proxg,evalf,evalg,x0,params)
plt.figure()
plt.plot(x)
plt.plot(xTrue)
plt.figure()
plt.semilogy(costs)
'norm(Ax,b)'對我來說看起來很奇怪。你的意思是「標準(Ax-b,2)」嗎? –
'x.d'是什麼意思? – littleO
d是另一個向量。理想情況下,我希望x和d是正交的,它是由delta值控制的。 – Erin