Быстрые решатели CVX в Matlab

Мне интересно, какой самый быстрый выпуклый оптимизатор в Matlab или есть ли способ ускорить текущие решатели? Я использую CVX, но решение моей проблемы с оптимизацией занимает целую вечность. У меня есть оптимизация, чтобы решить

minimize norm(Ax-b, 2)
subject to
x >= 0
and x d <= delta

где размеры A и b очень велики.

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


matlab optimization least-squares cvx convex-optimization
person Erin    schedule 05.12.2016    source источник
comment
norm(Ax,b) выглядит странно для меня. Вы имеете в виду norm(Ax-b,2)?   -  person Erwin Kalvelagen    schedule 06.12.2016
comment
Что означает x.d?   -  person littleO    schedule 20.12.2016
comment
d — другой вектор. В идеале я хочу, чтобы x и d были ортогональны, что контролируется значением дельты.   -  person Erin    schedule 29.12.2016

Ответы (1)


arrow_upward
0
arrow_downward

Я не уверен, что означает x.d <= delta, но я просто предполагаю, что это должно быть x <= delta.

Вы можете решить эту проблему, используя метод прогнозируемого градиента или метод ускоренного прогнозируемого градиента (который является лишь небольшой модификацией метода прогнозируемого градиента, который «волшебным образом» сходится намного быстрее). Вот некоторый код Python, который показывает, как минимизировать .5|| Ax - b ||^2 при условии, что 0 ‹= x ‹= дельта с использованием FISTA, который представляет собой ускоренный метод прогнозируемого градиента. Более подробную информацию о методе прогнозируемого градиента и FISTA можно найти, например, в рукописи Бойда. по проксимальным алгоритмам.

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)
person littleO    schedule 20.12.2016