Python cvxpy.norm() Examples

def __init__(self):
        """
        A large r_scale for a small scale problem will
        ead to numerical problems as parameters become excessively small
        and (it seems) precision is lost in the dynamics.
        """

        self.set_random_initial_state()

        self.x_init = np.concatenate(((self.m_wet,), self.r_I_init, self.v_I_init, self.q_B_I_init, self.w_B_init))
        self.x_final = np.concatenate(((self.m_dry,), self.r_I_final, self.v_I_final, self.q_B_I_final, self.w_B_final))

        self.r_scale = np.linalg.norm(self.r_I_init)
        self.m_scale = self.m_wet

        # slack variable for linear constraint relaxation
        self.s_prime = cvx.Variable((K, 1), nonneg=True)

        # slack variable for lossless convexification
        # self.gamma = cvx.Variable(K, nonneg=True) 

def test_proj_soc(self):
        import cvxpy as cp
        np.random.seed(0)
        n = 100
        for _ in range(15):
            x = np.random.randn(n)
            z = cp.Variable(n)
            objective = cp.Minimize(cp.sum_squares(z - x))
            constraints = [cp.norm(z[1:], 2) <= z[0]]
            prob = cp.Problem(objective, constraints)
            prob.solve(solver="SCS", eps=1e-10)
            p = cone_lib._proj(x, cone_lib.SOC, dual=False)
            np.testing.assert_allclose(
                p, np.array(z.value))
            np.testing.assert_allclose(
                p, cone_lib._proj(x, cone_lib.SOC, dual=True)) 

def running_example():
    print("running example")
    m = 20
    n = 10
    x = cp.Variable((n, 1))
    F = cp.Parameter((m, n))
    g = cp.Parameter((m, 1))
    lambd = cp.Parameter((1, 1), nonneg=True)
    objective_fn = cp.norm(F @ x - g) + lambd * cp.norm(x)
    constraints = [x >= 0]
    problem = cp.Problem(cp.Minimize(objective_fn), constraints)
    assert problem.is_dcp()
    assert problem.is_dpp()
    print("is_dpp: ", problem.is_dpp())

    F_t = torch.randn(m, n, requires_grad=True)
    g_t = torch.randn(m, 1, requires_grad=True)
    lambd_t = torch.rand(1, 1, requires_grad=True)
    layer = CvxpyLayer(problem, parameters=[F, g, lambd], variables=[x])
    x_star, = layer(F_t, g_t, lambd_t)
    x_star.sum().backward()
    print("F_t grad: ", F_t.grad)
    print("g_t grad: ", g_t.grad) 

def __init__(self):
        """
        A large r_scale for a small scale problem will
        ead to numerical problems as parameters become excessively small
        and (it seems) precision is lost in the dynamics.
        """

        # self.set_random_initial_state()

        self.x_init = np.concatenate(((self.m_wet,), self.r_I_init, self.v_I_init, self.q_B_I_init, self.w_B_init))
        self.x_final = np.concatenate(((self.m_dry,), self.r_I_final, self.v_I_final, self.q_B_I_final, self.w_B_final))

        self.r_scale = np.linalg.norm(self.r_I_init)
        self.m_scale = self.m_wet

        # self.nondimensionalize() 

def __init__(self, m, k, n, complex=False):
        if not cvx_available:
            raise RuntimeError('Cannot initialize when cvxpy is not available.')
        # Initialize parameters and variables
        A = cvx.Parameter((m, k), complex=complex)
        B = cvx.Parameter((m, n), complex=complex)
        l = cvx.Parameter(nonneg=True)
        X = cvx.Variable((k, n), complex=complex)
        # Create the problem
        # CVXPY issue:
        #   cvx.norm does not work if axis is not 0.
        # Workaround:
        #   use cvx.norm(X.T, 2, axis=0) instead of cvx.norm(X, 2, axis=1)
        obj_func = 0.5 * cvx.norm(cvx.matmul(A, X) - B, 'fro')**2 + \
                   l * cvx.sum(cvx.norm(X.T, 2, axis=0))
        self._problem = cvx.Problem(cvx.Minimize(obj_func))
        self._A = A
        self._B = B
        self._l = l
        self._X = X 

def test_simple_batch_socp(self):
        set_seed(243)
        n = 5
        m = 1
        batch_size = 4

        P_sqrt = cp.Parameter((n, n), name='P_sqrt')
        q = cp.Parameter((n, 1), name='q')
        A = cp.Parameter((m, n), name='A')
        b = cp.Parameter((m, 1), name='b')

        x = cp.Variable((n, 1), name='x')

        objective = 0.5 * cp.sum_squares(P_sqrt @ x) + q.T @ x
        constraints = [A@x == b, cp.norm(x) <= 1]
        prob = cp.Problem(cp.Minimize(objective), constraints)

        prob_tch = CvxpyLayer(prob, [P_sqrt, q, A, b], [x])

        P_sqrt_tch = torch.randn(batch_size, n, n, requires_grad=True)
        q_tch = torch.randn(batch_size, n, 1, requires_grad=True)
        A_tch = torch.randn(batch_size, m, n, requires_grad=True)
        b_tch = torch.randn(batch_size, m, 1, requires_grad=True)

        torch.autograd.gradcheck(prob_tch, (P_sqrt_tch, q_tch, A_tch, b_tch)) 

def get_constraints(self, X_v, U_v, X_last_p, U_last_p):
        """
        Get model specific constraints.

        :param X_v: cvx variable for current states
        :param U_v: cvx variable for current inputs
        :param X_last_p: cvx parameter for last states
        :param U_last_p: cvx parameter for last inputs
        :return: A list of cvx constraints
        """
        # Boundary conditions:
        constraints = [
            X_v[:, 0] == self.x_init,
            X_v[:, -1] == self.x_final,

            U_v[:, 0] == 0,
            U_v[:, -1] == 0
        ]

        # Input conditions:
        constraints += [
            0 <= U_v[0, :],
            U_v[0, :] <= self.v_max,
            cvx.abs(U_v[1, :]) <= self.w_max,
        ]

        # State conditions:
        constraints += [
            X_v[0:2, :] <= self.upper_bound - self.robot_radius,
            X_v[0:2, :] >= self.lower_bound + self.robot_radius,
        ]

        # linearized obstacles
        for j, obst in enumerate(self.obstacles):
            p = obst[0]
            r = obst[1] + self.robot_radius

            lhs = [(X_last_p[0:2, k] - p) / (cvx.norm((X_last_p[0:2, k] - p)) + 1e-6) * (X_v[0:2, k] - p)
                   for k in range(K)]
            constraints += [r - cvx.vstack(lhs) <= self.s_prime[j]]
        return constraints 

def test_equil(self):
        """Test equilibration.
        """
        from proximal.algorithms.equil import newton_equil
        np.random.seed(1)
        kernel = np.array([1, 1, 1]) / np.sqrt(3)
        kernel_mat = np.ones((3, 3)) / np.sqrt(3)
        x = px.Variable(3)
        wr = np.array([10, 5, 7])
        K = px.mul_elemwise(wr, x)
        K = px.conv(kernel, K)
        wl = np.array([100, 50, 3])
        K = px.mul_elemwise(wl, K)
        K = px.CompGraph(K)

        # Equilibrate
        gamma = 1e-1
        d, e = px.equil(K, 1000, gamma=gamma, M=5)
        tmp = d * wl * kernel_mat * wr * e
        u, v = np.log(d), np.log(e)
        obj_val = np.square(tmp).sum() / 2 - u.sum() - v.sum() + \
            gamma * (np.linalg.norm(v)**2 + np.linalg.norm(u)**2)

        d, e = newton_equil(wl * kernel_mat * wr, gamma, 100)
        tmp = d * wl * kernel_mat * wr * e
        u, v = np.log(d), np.log(e)
        sltn_val = np.square(tmp).sum() / 2 - u.sum() - v.sum() + \
            gamma * (np.linalg.norm(v)**2 + np.linalg.norm(u)**2)
        self.assertAlmostEqual((obj_val - sltn_val) / sltn_val, 0, places=3) 

def _project(self, matrix):
        if self.R >= cvx.norm(matrix, 2).value >= self.r:
            return matrix
        elif cvx.norm(matrix, 2).value == 0:
            result = np.ones(self.shape)
            return self.r*result/cvx.norm(result, 2).value
        elif cvx.norm(matrix, 2).value < self.r:
            return self.r*matrix/cvx.norm(matrix, 2).value
        else:
            return self.R*matrix/cvx.norm(matrix, 2).value 

def test_not_enough_parameters(self):
        x = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        lam2 = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1) + lam2 * cp.sum_squares(x)
        prob = cp.Problem(cp.Minimize(objective))
        with self.assertRaisesRegex(ValueError, "The layer's parameters.*"):
            CvxpyLayer(prob, [lam], [x])  # noqa: F841 

def test_not_enough_parameters_at_call_time(self):
        x = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        lam2 = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1) + lam2 * cp.sum_squares(x)
        prob = cp.Problem(cp.Minimize(objective))
        layer = CvxpyLayer(prob, [lam, lam2], [x])
        with self.assertRaisesRegex(
                ValueError,
                'A tensor must be provided for each CVXPY parameter.*'):
            layer(lam) 

def test_non_dpp(self):
        x = cp.Variable(1)
        y = cp.Variable(1)
        lam = cp.Parameter(1)
        objective = lam * cp.norm(x, 1)
        prob = cp.Problem(cp.Minimize(objective))
        with self.assertRaisesRegex(ValueError, 'Problem must be DPP.'):
            CvxpyLayer(prob, [lam], [x, y])  # noqa: F841 

def test_too_many_variables(self):
        x = cp.Variable(1)
        y = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1)
        prob = cp.Problem(cp.Minimize(objective))
        with self.assertRaisesRegex(ValueError, 'Argument `variables`.*'):
            CvxpyLayer(prob, [lam], [x, y])  # noqa: F841 

def test_least_squares(self):
        set_seed(243)
        m, n = 100, 20

        A = cp.Parameter((m, n))
        b = cp.Parameter(m)
        x = cp.Variable(n)
        obj = cp.sum_squares(A@x - b) + cp.sum_squares(x)
        prob = cp.Problem(cp.Minimize(obj))
        prob_th = CvxpyLayer(prob, [A, b], [x])

        A_th = torch.randn(m, n).double().requires_grad_()
        b_th = torch.randn(m).double().requires_grad_()

        x = prob_th(A_th, b_th, solver_args={"eps": 1e-10})[0]

        def lstsq(
            A,
            b): return torch.solve(
            (A_th.t() @ b_th).unsqueeze(1),
            A_th.t() @ A_th +
            torch.eye(n).double())[0]
        x_lstsq = lstsq(A_th, b_th)

        grad_A_cvxpy, grad_b_cvxpy = grad(x.sum(), [A_th, b_th])
        grad_A_lstsq, grad_b_lstsq = grad(x_lstsq.sum(), [A_th, b_th])

        self.assertAlmostEqual(
            torch.norm(
                grad_A_cvxpy -
                grad_A_lstsq).item(),
            0.0)
        self.assertAlmostEqual(
            torch.norm(
                grad_b_cvxpy -
                grad_b_lstsq).item(),
            0.0) 

def test_not_enough_parameters(self):
        x = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        lam2 = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1) + lam2 * cp.sum_squares(x)
        prob = cp.Problem(cp.Minimize(objective))
        with self.assertRaises(ValueError):
            layer = CvxpyLayer(prob, [lam], [x])  # noqa: F841 

def test_not_enough_parameters_at_call_time(self):
        x = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        lam2 = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1) + lam2 * cp.sum_squares(x)
        prob = cp.Problem(cp.Minimize(objective))
        layer = CvxpyLayer(prob, [lam, lam2], [x])  # noqa: F841
        with self.assertRaisesRegex(
                ValueError,
                'A tensor must be provided for each CVXPY parameter.*'):
            layer(lam) 

def test_too_many_variables(self):
        x = cp.Variable(1)
        y = cp.Variable(1)
        lam = cp.Parameter(1, nonneg=True)
        objective = lam * cp.norm(x, 1)
        prob = cp.Problem(cp.Minimize(objective))
        with self.assertRaises(ValueError):
            layer = CvxpyLayer(prob, [lam], [x, y])  # noqa: F841 

def sparse_lowrank_mse(name_size):
    name, size = name_size
    print(name, size)
    matrix = named_target_matrix(name, size)
    M = matrix
    lambda1 = cp.Parameter(nonneg=True)
    lambda2 = cp.Parameter(nonneg=True)
    L = cp.Variable((size, size))
    S = cp.Variable((size, size))
    prob = cp.Problem(cp.Minimize(cp.sum_squares(M - L - S) / size**2 + lambda1 / size * cp.norm(L, 'nuc') + lambda2 / size**2 * cp.norm(S, 1)))

    result = []
    for _ in range(ntrials):
        l1 = np.exp(np.random.uniform(np.log(1e-2), np.log(1e4)))
        l2 = np.exp(np.random.uniform(np.log(1e-2), np.log(1e4)))
        lambda1.value = l1
        lambda2.value = l2
        try:
            prob.solve()
            nnz = (np.abs(S.value) >= 1e-7).sum()
            singular_values = np.linalg.svd(L.value, compute_uv=False)
            rank = (singular_values >= 1e-7).sum()
            n_params = nnz + 2 * rank * size
            mse = np.sum((matrix - L.value - S.value)**2) / size**2
            result.append((n_params, mse))
        except:
            pass
    budget = 2 * size * np.log2(size)
    if model[name] == 'BPBP':
        budget *= 2
    eligible = [res for res in result if res[0] <= budget]
    if eligible:
        mse = min(m for (n_params, m) in eligible)
    else:
        mse = np.sum(matrix**2) / size**2
    print(name, size, 'done')
    return (name, size, mse) 

def __init__(self, m, k, formulation='penalizedl1', nonnegative=False):
        if not cvx_available:
            raise RuntimeError('Cannot initialize when cvxpy is not available.')
        # Initialize parameters and variables
        A = cvx.Parameter((m, k))
        b = cvx.Parameter((m, 1))
        l = cvx.Parameter(nonneg=True)
        x = cvx.Variable((k, 1))
        # Create the problem
        if formulation == 'penalizedl1':
            obj_func = 0.5 * cvx.sum_squares(cvx.matmul(A, x) - b) + l * cvx.norm1(x)
            constraints = []
        elif formulation == 'constrainedl1':
            obj_func = cvx.sum_squares(cvx.matmul(A, x) - b)
            constraints = [cvx.norm1(x) <= l]
        elif formulation == 'constrainedl2':
            obj_func = cvx.norm1(x)
            constraints = [cvx.norm(cvx.matmul(A, x) - b) <= l]
        else:
            raise ValueError("Unknown formulation '{0}'.".format(formulation))
        if nonnegative:
            constraints.append(x >= 0)
        problem = cvx.Problem(cvx.Minimize(obj_func), constraints)
        self._formulation = formulation
        self._A = A
        self._b = b
        self._l = l
        self._x = x
        self._obj_func = obj_func
        self._constraints = constraints
        self._problem = problem 

def _calc_wts(self, x_i):
        distances = np.array(
            [np.linalg.norm(self.X_[i, :] - x_i) for i in range(self.X_.shape[0])]
        )
        weights = np.exp(-(distances ** 2) / self.sigma)
        if self.span:
            weights = weights * (distances <= np.quantile(distances, q=self.span))
        return weights 

def _solve(self, sensitive, X, y):
        n_obs, n_features = X.shape
        theta = cp.Variable(n_features)
        y_hat = X @ theta

        log_likelihood = cp.sum(
            cp.multiply(y, y_hat)
            - cp.log_sum_exp(
                cp.hstack([np.zeros((n_obs, 1)), cp.reshape(y_hat, (n_obs, 1))]), axis=1
            )
        )
        if self.penalty == "l1":
            log_likelihood -= cp.sum((1 / self.C) * cp.norm(theta[1:]))

        constraints = self.constraints(y_hat, y, sensitive, n_obs)

        problem = cp.Problem(cp.Maximize(log_likelihood), constraints)
        problem.solve(max_iters=self.max_iter)

        if problem.status in ["infeasible", "unbounded"]:
            raise ValueError(f"problem was found to be {problem.status}")

        self.n_iter_ = problem.solver_stats.num_iters

        if self.fit_intercept:
            self.coef_ = theta.value[np.newaxis, 1:]
            self.intercept_ = theta.value[0:1]
        else:
            self.coef_ = theta.value[np.newaxis, :]
            self.intercept_ = np.array([0.0]) 

def block_solve_cvx(r_j, A_j, a_1_j, a_2_j, m, b_j_init, verbose=False):
    # pylint: disable=unused-argument
    """Solve the gelcd optimization problem for a single block with cvx.

    b_j_init and verbose are ignored. b_j_init because cvx doesn't support it.
    verbose because it doesn't go together with tqdm.
    """
    # Convert everything to numpy.
    device = A_j.device
    dtype = A_j.dtype
    r_j = r_j.cpu().numpy()
    A_j = A_j.cpu().numpy()

    # Create the b_j variable.
    b_j = cvx.Variable(A_j.shape[1])

    # Form the objective.
    q_j = r_j - A_j * b_j
    obj_fun = cvx.square(cvx.norm(q_j)) / (2.0 * m)
    obj_fun += a_1_j * cvx.norm(b_j) + (a_2_j / 2.0) * cvx.square(cvx.norm(b_j))

    # Build the optimization problem.
    obj = cvx.Minimize(obj_fun)
    problem = cvx.Problem(obj, constraints=None)

    problem.solve(solver="CVXOPT", verbose=False)
    b_j = np.asarray(b_j.value)
    return torch.from_numpy(b_j).to(device, dtype) 

def _obj(self, b_0, b):
        """Compute the objective function value for the given b_0, b."""
        r = self.y - b_0 - self.X @ b
        g_b = r @ r / (2.0 * self.m)
        b_j_norms = [np.linalg.norm(b[self.groups[j]], ord=2) for j in range(self.p)]
        h_b = self.l_1 * sum(
            np.sqrt(len(self.groups[j])) * b_j_norms[j] for j in range(self.p)
        )
        h_b += self.l_2 * sum(
            np.sqrt(len(self.groups[j])) * (b_j_norms[j] ** 2) for j in range(self.p)
        )
        return g_b + h_b 

def least_squares(m, n):
    """ Create a least squares problem with m datapoints and n dimensions """
    A = np.random.randn(m, n)
    _x = np.random.randn(n)
    b = A.dot(_x)

    x = cp.Variable(n)
    return (x, cp.Problem(cp.Minimize(cp.sum_squares(A * x - b) + cp.norm(x, 2)))) 

def robust_pca(p, suppfrac):
    """ Create a robust PCA problem with a low rank matrix """

    # First, create a rank = "rk" matrix:
    rk = int(round(p * 0.5))
    assert rk <= min(p, p)
    Lstar = np.zeros((p, p))
    for i in range(rk):
        vi = np.random.randn(p, 1)
        mati = vi.dot(vi.T)
        Lstar += mati

    # Then, create a sparse matrix:
    Mstar = np.random.randn(p, p)
    Mstar_vec = Mstar.T.ravel()

    nnz = int(np.floor((1.0 - suppfrac) * p * p))        # Num. nonzeros
    assert nnz <= p * p
    idxes = np.random.randint(0, p * p, nnz)
    Mstar_vec[idxes] = 0
    Mstar = np.reshape(Mstar_vec, (p, p)).T

    # Finally, sum the two matrices "L" and "M":
    X = Lstar + Mstar

    lam = 1.0

    Lhat = cp.Variable(p, p)
    Mhat = cp.Variable(p, p)
    prob = cp.Problem(cp.Minimize(cp.norm(Lhat, "nuc") + cp.sum_squares(Lhat)),
                      [cp.norm(Mhat, 1) <= lam, Lhat + Mhat == X])

    data = prob.get_problem_data(cp.SCS)
    data['beta_from_x'] = cvxpy_beta_from_x(prob, (Lhat, Mhat), data['A'].shape[0])
    return ((Lhat, Mhat), prob, data) 

def test_proj_exp(self):
        import cvxpy as cp
        np.random.seed(0)
        for _ in range(15):
            x = np.random.randn(9)
            var = cp.Variable(9)
            constr = [cp.constraints.ExpCone(var[0], var[1], var[2])]
            constr.append(cp.constraints.ExpCone(var[3], var[4], var[5]))
            constr.append(cp.constraints.ExpCone(var[6], var[7], var[8]))
            obj = cp.Minimize(cp.norm(var[0:3] - x[0:3]) +
                              cp.norm(var[3:6] - x[3:6]) +
                              cp.norm(var[6:9] - x[6:9]))
            prob = cp.Problem(obj, constr)
            prob.solve(solver="SCS", eps=1e-12)
            p = cone_lib._proj(x, cone_lib.EXP, dual=False)
            np.testing.assert_allclose(p, var.value, atol=1e-6)
            # x + Pi_{exp}(-x) = Pi_{exp_dual}(x)
            p_dual = cone_lib._proj(x, cone_lib.EXP_DUAL, dual=False)
            var = cp.Variable(9)
            constr = [cp.constraints.ExpCone(var[0], var[1], var[2])]
            constr.append(cp.constraints.ExpCone(var[3], var[4], var[5]))
            constr.append(cp.constraints.ExpCone(var[6], var[7], var[8]))
            obj = cp.Minimize(cp.norm(var[0:3] + x[0:3]) +
                              cp.norm(var[3:6] + x[3:6]) +
                              cp.norm(var[6:9] + x[6:9]))
            prob = cp.Problem(obj, constr)
            prob.solve(solver="SCS", eps=1e-12)
            np.testing.assert_allclose(
                p_dual, x + var.value, atol=1e-6) 

def dist(self, matrix):
        """Distance from matrix to projection.
        """
        proj_mat = self.project(matrix)
        return cvxpy.norm(cvxpy.vec(matrix - proj_mat), 2).value

    # Project the matrix into the space defined by the non-convex constraint.
    # Returns the updated matrix. 

def relax(self):
        constr = super(Sphere, self).relax()
        return constr + [cvx.norm(self, 2) <= 1] 

def relax(self):
        if self.M is None:
            return []
        else:
            return [cvx.norm(self, 2) <= self.M] 

def init_z(self, random):
        """Initializes the value of the replicant variable.
        """
        if random:
            length = np.random.uniform()
            direction = np.random.randn(self.shape)
            self.z.value = length*direction/norm(direction, 2).value
        else:
            self.z.value = np.zeros(self.shape)