首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >Matlab中的快速CVX求解器

Matlab中的快速CVX求解器
EN

Stack Overflow用户
提问于 2016-12-06 05:59:58
回答 1查看 698关注 0票数 0

我想知道Matlab中最快的凸优化程序是什么,或者有什么方法可以加速当前的求解器?我正在使用CVX,但它永远解决不了我的优化问题。我要做的优化就是要解决

代码语言:javascript
复制
minimize norm(Ax-b, 2)
subject to
x >= 0
and x d <= delta

其中a和b的大小非常大。

有没有什么方法可以用最小二乘求解器来解决这个问题,然后把它转移到约束版本中,使它更快?

EN

回答 1

Stack Overflow用户

发布于 2016-12-20 08:45:21

我不确定x.d <= delta是什么意思,但我只是假设它应该是x <= delta

您可以使用投影梯度法或加速投影梯度法来解决此问题(这只是对投影梯度法的轻微修改,它“神奇地”收敛得更快)。下面是一些python代码,它展示了如何使用FISTA最小化.5|| Ax -b|^2,它受0 <= x <=增量的约束,这是一种加速的投影梯度方法。关于投影梯度法和FISTA的更多细节可以在博伊德关于近邻算法的manuscript中找到。

代码语言:javascript
复制
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)
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/40984131

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档