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

Question

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

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

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

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

Answer 1

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

Вы можете решить эту проблему, используя метод проецируемого градиента или метод ускоренного проецируемого градиента (который является лишь небольшой модификацией метода проецируемого градиента, который "волшебным образом" сходится гораздо быстрее). Вот код Python, который показывает, как минимизировать.5|| Ax - b ||^2 при условии ограничения 0 <= x <= delta с использованием 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)