Schnelle CVX-Löser in Matlab

Question

Ich frage mich, was ist der schnellste konvexe Optimierer in Matlab oder gibt es eine Möglichkeit, aktuelle Löser zu beschleunigen? Ich benutze CVX, aber es dauert ewig, um das Optimierungsproblem zu lösen, das ich habe. Die Optimierung ich habe, ist

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

, wo die Größe von A und b zu lösen, ist sehr groß.

Gibt es eine Möglichkeit, dass ich das mit einem Least-Square-Solver lösen und dann auf die Constraint-Version übertragen kann, um es schneller zu machen?

Answer 1

Ich bin mir nicht sicher, was x.d <= delta bedeutet, aber ich nehme an, es soll x <= delta sein.

Sie können dieses Problem mit der projizierten Gradientenmethode oder einer beschleunigten projizierten Gradientenmethode lösen (dies ist nur eine geringfügige Modifikation der projizierten Gradientenmethode, die "magisch" viel schneller konvergiert). Hier ist ein Python-Code, der zeigt, wie man .5 || minimiert Axe - b ||^2 unterliegt der Einschränkung, dass 0 < = x < = Delta unter Verwendung von FISTA ist, was eine beschleunigte projizierte Gradientenmethode ist. Weitere Details über die projizierte Gradientenmethode und FISTA können zum Beispiel in Boyds manuscript über proximale Algorithmen gefunden werden.

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)