aboutsummaryrefslogtreecommitdiffstats
path: root/src/convex_optimization.py
blob: 1d84db225ffa4d2bd58f9bbdf49aa34b6bcfed2d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
import theano
import cascade_creation
from theano import tensor, function
import numpy as np
import timeout
import cvxopt


@timeout.timeout(20)
def sparse_recovery(M_val, w_val, lbda):
    """
    Solves:
        min lbda * |theta|_1 - sum b^i log(exp(-<M[i], theta_>) -1) - M*theta
        s.t theta_j <= 0
    """
    assert len(M_val) == len(w_val)

    if M_val.dtype == bool:
        M_val = M_val.astype('float32')

    if type(lbda) == int:
        lbda = np.array(lbda)

    theta = tensor.row().T
    theta_ = theta.flatten()

    M = theano.shared(M_val.astype(theano.config.floatX))
    w = theano.shared(w_val.astype(theano.config.floatX))
    lbda = theano.shared(lbda.astype(theano.config.floatX))

    y_tmp = lbda * (theta_).norm(1) - tensor.sum(tensor.dot(M, theta_
        )) + w.dot(tensor.log(tensor.exp(-M.dot(theta_))-1))
    y_diff = tensor.grad(y_tmp, theta_)
    y_hess = theano.gradient.hessian(y_tmp, theta_)
    f = function([theta_], y_hess)

    print(f(-1*np.random.rand(M_val.shape[1]).astype("float32")))

    y = lbda * (theta_).norm(1) - tensor.sum(tensor.dot(M, theta_
        )) - w.dot(tensor.log(tensor.exp(-M.dot(theta_))-1))

    return diff_and_opt(theta, theta_, M, M_val, w, lbda, y)


@timeout.timeout(20)
def type_lasso(M_val, w_val, lbda):
    """
    Solves:
        min - sum_j theta_j + lbda*|e^{M*theta} - (1 - w)|_2
        s.t theta_j <= 0
    """
    assert len(M_val) == len(w_val)

    if M_val.dtype == bool:
        M_val = M_val.astype('float32')

    if type(lbda) == int:
        lbda = np.array(lbda)

    theta = tensor.row().T
    theta_ = theta.flatten()

    M = theano.shared(M_val.astype(theano.config.floatX))
    w = theano.shared(w_val.astype(theano.config.floatX))
    lbda = theano.shared(lbda.astype(theano.config.floatX))

    y = (theta_).norm(1) + lbda * (
        tensor.exp(M.dot(theta_)) - (1 - w)).norm(2)

    return diff_and_opt(theta, theta_, M, M_val, w, lbda, y)


def diff_and_opt(theta, theta_, M, M_val, w, lbda, y):
    z = tensor.row().T
    z_ = z.flatten()

    m, n = M_val.shape

    y_diff = tensor.grad(y, theta_)
    y_hess = z[0] * theano.gradient.hessian(y, theta_)
    f_x = theano.function([theta], [y, y_diff], allow_input_downcast=True)
    f_xz = theano.function([theta, z], [y, y_diff, y_hess],
                           allow_input_downcast=True)

    def F(x=None, z=None):
        if x is None:
            return 0, cvxopt.matrix(.0001, (n,1))
        elif z is None:
            y, y_diff = f_x(x)
            return cvxopt.matrix(float(y), (1, 1)),\
                   cvxopt.matrix(y_diff.astype("float64")).T
        else:
            y, y_diff, y_hess = f_xz(x, z)
            return cvxopt.matrix(float(y), (1, 1)), \
                    cvxopt.matrix(y_diff.astype("float64")).T, \
                    cvxopt.matrix(y_hess.astype("float64"))

    G = cvxopt.spdiag([1 for i in range(n)])
    h = cvxopt.matrix(0.0, (n,1))

    cvxopt.solvers.options['show_progress'] = True
    try:
        theta = cvxopt.solvers.cp(F, G, h)['x']
    except ArithmeticError:
        print("ArithmeticError thrown, change initial point"+\
              " given to the solver")

    return 1 - np.exp(theta), theta


def test():
    """
    unit test
    """
    lbda = 10
    G = cascade_creation.InfluenceGraph(max_proba=.8)
    G.erdos_init(n=100, p = .1)
    A = cascade_creation.generate_cascades(G, .1, 2000)
    M_val, w_val = cascade_creation.icc_matrixvector_for_node(A, 0)
    p_vec, theta = l1obj_l2penalization(M_val, w_val, lbda)
    print(p_vec)

if __name__=="__main__":
    test()