OpenOpt least norms routine

# OpenOpt routine to solve min {alpha1*||A1 x - b1||_1 + alpha2*||A2 x - b2||^2 + beta1 * ||x||_1 + beta2 * ||x||^2}
# optionally: with specifiable accuracy fTol > 0 (by gsubg and maybe amsg2p, latter requires known good enough fOpt estimation)

# made by Dmitrey

def leastNorms(alpha1 = 0, A1 = None, b1 = None, alpha2 = 0, A2 = None, b2 = None, beta1 = 0, beta2 = 0, x0 = None, fTol = None,  solver = 'gsubg'):
    # solves min alpha1*||A1 x - b1||_1 + alpha2*||A2 x - b2||^2 + beta1 * ||x||_1 + beta2 * ||x||^2 -> min
    # with specifiable accuracy fTol

    assert (alpha1 != 0 and A1 is not None and b1 is not None) or (alpha2 != 0 and A2 is not None and b2 is not None), \
    'at least one of sets (alpha1, A1, b1), (alpha2, A2, b2) must be provided'

    if alpha1 != 0:
        n = A1.shape[1]
    else:
        n = A2.shape[1]

    from openopt import oosolver, NSP
    if type(solver) == str:
        solver = oosolver(solver)
    if solver.__name__ == 'gsubg':
        assert fTol is not None, 'for solver gsubg non-default fTol value must be provided'

    import numpy as np
    try:
        import scipy.sparse as sp
        from scipy.sparse import isspmatrix
    except ImportError:
        isspmatrix = lambda *args: False
        sp = None


    matMultVec =  lambda x, y: np.dot(x, y) if not isspmatrix(x) else x._mul_sparse_matrix(sp.csr_matrix(y.reshape((y.size, 1)))).A.flatten()

    if isspmatrix(A1):
        A1 = A1.tocsc()
        A1T = A1.T.tocsc()
    else:
        A1T = A1.T
    if isspmatrix(A2):
        A2 = A2.tocsc()
        A2T = A2.T.tocsc()
    else:
        A2T = A2.T

    if isspmatrix(b1):
        b1 = b1.A.flatten()
    if isspmatrix(b2):
        b2 = b2.A.flatten()


    def f(x):
        r = 0.0

        # I don't want to use np.linalg norm for better compatibility with PyPy
        if alpha1 != 0:
            tmp = np.abs(matMultVec(A1, x) - b1)
            assert tmp.ndim == 1, 'bug in leastNorms, probably due to incorrect input'
            r += alpha1 * np.sum(tmp)
        if alpha2 != 0:
            tmp = (matMultVec(A2, x) - b2)**2
            assert tmp.ndim == 1, 'bug in leastNorms, probably due to incorrect input'
            r += alpha2 * np.sum(tmp)

        if beta1 != 0:
            r += beta1 * np.sum(np.abs(x))
        if beta2 != 0:
            r += beta2 * np.sqrt(np.sum(x**2))

        return r

    def df(x):
        r = np.zeros(x.size)

        if alpha1 != 0:
            r += alpha1 * matMultVec(A1T, np.sign(matMultVec(A1, x) - b1))
        if alpha2 != 0:
            r += 2 * alpha2 * matMultVec(A2T, matMultVec(A2,x)  - b2)

        if beta1 != 0:
            r += beta1 * np.sign(x)
        if beta2 != 0:
            r += 2 * beta2 * x

        return r

    p = NSP(x0 = np.zeros(n) if x0 is None else x0, f = f, df = df, fTol = fTol)
    #p.checkdf()
    r = p.minimize(solver)
    return r

if __name__ == '__main__':
    n = 10000 # unknowns
    m = 50
    import scipy.sparse as sp

    X = sp.rand(n, 1, density=0.8)

    # Let A1 be matrix of shape (m, n)
    alpha1 = 0.015
    A1 = sp.rand(m, n, density=5.0/n) # in average 5 nonzero elements per each of m columns
    b1 = A1._mul_sparse_matrix(X) + 0.015* sp.rand(m, 1, density=0.8)

    # Let A2 be matrix of shape (m+10, n)
    alpha2 = 0.015
    A2 = sp.rand(m+10, n, density=5.0/n)
    b2 = A2._mul_sparse_matrix(X) + 0.015* sp.rand(m+10, 1, density=0.8)

    beta1 = 0.0005
    beta2 = 0.001

    fTol = 0.01
    x0 = None

    from openopt import oosolver
    # ralg is suitable for medium-scaled problems with nVars (=n) ~ 1000, see http://openopt.org/ralg for details
#    solver = oosolver('ralg', iprint=10)
#    r = leastNorms(alpha1, A1, b1, alpha2, A2, b2, beta1, beta2, x0, fTol, solver)
    # gsubg is suitable for large-scaled problems, see http://openopt.org/gsubg for details
    solver = oosolver('gsubg', iprint=10, maxShoots = 20)
    r = leastNorms(alpha1, A1, b1, alpha2, A2, b2, beta1, beta2, x0, fTol, solver)
    # alternatively you can use leastNorms(alpha2 = alpha2, b1 = b1, alpha1 = alpha1 ...)
'''
Output (may differ due to random variables involved in initial data creation)
------------------------- OpenOpt 0.39 -------------------------
solver: gsubg   problem: unnamed    type: NSP   goal: minimum
 iter   objFunVal
    0  2.018e+00
   10  1.123e-01
   20  9.871e-02
   30  9.267e-02
   40  9.063e-02
   50  8.691e-02
   60  8.544e-02
   70  8.341e-02
   80  8.249e-02
   90  8.216e-02
  100  8.200e-02
Terminated (singular KKT matrix).
  104  8.142e-02
istop: 16 (optimal solution wrt required fTol has been obtained)
Solver:   Time Elapsed = 90.28  CPU Time Elapsed = 85.07
objFunValue: 0.081420241
'''