Matlabの高速CVXソルバー

Question

Matlabで最も高速な凸オプティマイザとは何ですか、現在のソルバを高速化する方法はありますか？私はCVXを使用していますが、それは私が持っている最適化の問題を解決するために永遠に取っています。私が持っている 最適化は、Aの大きさとbが非常に大きい

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

を解決することです。

最小二乗ソルバーでこれを解き、それをより速くするために制約バージョンに転送する方法はありますか？

Answer 1

私はx.d <= deltaの意味は分かりませんが、私はそれがx <= deltaであると仮定します。

この問題は、投影勾配法または加速勾配法（投影勾配法のほんのわずかな変更ですが、「魔法のように」収束するほうがはるかに速くなります）を使用して解決できます。以下は、最小化する方法を示すPythonコードです。 Ax - 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)