% SDP-based method for integer least squares (with concave quadratic cuts)
%
% Finds an SDP-based lower bound and an upper bound, along with a feasible 
% point that achieves the upper bound, to the following integer problem:
%
% minimize    f(x) = ||Ax - b||_2^2
% subject to  x in Z^n,
%
% where A^T A is positive definite. Utilizes concave quadratic cuts hold
% for every integer points, namely
%
%   (a^T x - b)(a^T x - (b + 1)) >= 0
%
% for integer-valued a and b.
%
% Input arguments:
%   A: m-by-n real matrix
%   b: m-by-1 real vector
%   K: number of samples used in the randomized algorithm, set to K = 3n by
%      default
%
% Outputs:
%   lb1: an SDP-based lower bound on the optimal value of the problem
%   lb2: same SDP-based bound with concave quadratic cuts
%   ub: an upper bound on the optimal value of the problem
%   xhat: an integer point that achieves f(xhat) = ub
%
% Notes:
%   Requires CVX. If MOSEK solver is available, use the code as it is. If 
%   not, comment out line 53, or change the solver to one that's available
%   to the system.

% Copyright (c) 2015. Jaehyun Park and Stephen Boyd

function [lb1, lb2, ub, xhat] = ILS_concave_cuts(A, b, K)

    % Input sanity check
    assert(nargin >= 2);
    P = A'*A; q = -A'*b; r = b'*b;
    n = size(P, 1);
    if nargin < 3
        K = 3*n;
    end
    chol(P); % Checks positive definiteness of P

    % Change of coordinates
    xcts = -P\q;
    xflr = floor(xcts);
    [q, r] = NEW_ORIGIN(xflr, P, q, r); % Translate by -xflr

    % Use CVX to solve initial SDP
    M = [P, q; q', r];
    cvx_solver mosek % Needs CVX professional license
    cvx_begin quiet
        variables X(n, n) x(n)
        expression Z
        Z = [X, x; x', 1];
        minimize sum(sum(M.*Z))
        subject to
            Z == semidefinite(n+1);
            diag(X) >= x;
    cvx_end
    lb1 = cvx_optval;
    X = (X + X')/2; % Force symmetry

    % Add concave quadratic cuts
    As = []; Bs = [];
    for i = 1:n
        for j = i+1:n
            a = zeros(n, 1);
            a(i) = 1; a(j) = 1;
            ax = x(i) + x(j);
            aXa = X(i, i) + 2*X(i, j) + X(j, j);
            b = floor(ax);
            if aXa - (2*b+1)*ax + b*(b+1) < 0 % found cut
               As = [As, a]; Bs = [Bs; b];
            end

            a(i) = 1; a(j) = -1;
            ax = x(i) - x(j);
            aXa = X(i, i) - 2*X(i, j) + X(j, j);
            b = floor(ax);
            if aXa - (2*b+1)*ax + b*(b+1) < 0 % found cut
               As = [As, a]; Bs = [Bs; b];
            end
        end
    end

    % Use CVX to solve SDP
    cvx_begin quiet
        variables X(n, n) x(n)
        expression Z
        Z = [X, x; x', 1];
        minimize sum(sum(M.*Z))
        subject to
            Z == semidefinite(n+1);
            diag(X) >= x;
            % Constraints corresponding to concave quadratic cuts
            diag(As'*X*As)-(2*Bs+1).*(As'*x)+Bs.*(Bs+1) >= 0;
    cvx_end
    lb2 = cvx_optval;

    % Randomized algorithm
    mu = x; Sigma = X - x*x';
    % Using eigenvalue decomposition because chol() complains about
    % near-zero negative eigenvalues
    [V, D] = eig(Sigma);
    A = V*sqrt(max(D, 0));

    ub = r; xhat = zeros(n, 1);
    for k = 1:K
        x = round(mulrandn_cached(mu, A));
        [x, val] = ONEOPT(x, P, q);
        if ub > val + r
            % Undo the change of coordinates
            ub = val + r; xhat = x + xflr;
        end
    end
end



% Subroutine for change of coordinates

function [q2, r2] = NEW_ORIGIN(x0, P, q, r)
    q2 = P*x0 + q;
    r2 = x0'*P*x0 + 2*q'*x0 + r;
end


% Subroutine for multivariate normal distribution
%
% Samples from a multivariate normal distribution with mean mu
% and square root matrix of covariance, A (covariance = A*A').
%
% Reference:
% http://en.wikipedia.org/wiki/Multivariate_normal_distribution

function x = mulrandn_cached(mu, A)
    n = size(mu, 1);
    z = randn(n, 1);
    x = mu + A*z;
end


% 1-opt greedy descent subroutine
%
% Starting from x, greedily chooses and optimizes over a single coordinate 
% to reduce the objective x^T P x + 2 q^T x the most. Stops when no single 
% coordinate change can reduce the objective. Returns the 1-opt point and 
% the function value.

function [x, val] = ONEOPT(x, P, q)
    g = 2*(P*x+q);
    v = diag(P);
    iters = 0;
    while true
        iters = iters + 1;
        if v >= abs(g)
            break;
        end
        c = round(-g./(2*v));
        diffs = (c.^2).*v + c.*g;
        [~, i] = min(diffs);
        x(i) = x(i) + c(i);
        g = g + 2*c(i)*P(:, i);
    end
    val = x'*P*x + 2*q'*x;
end
