Boolean Quadratic Programming

Question

I am new to optimization and I am trying to understand concepts of semi-definite relaxation (SDR) through examples. It seems my understanding of this topic is not fully clear as I will show in details below.

Let us assume we have the following optimization problem $$\mathbb{OP1}: \min_{\mathbf{a} \in \mathbb{R}^n} \quad \quad \mathbf{a}^{\rm{T}} G \mathbf{a} - 2 \mathbf{y}^{\rm{T}} \mathbf{a} \\ {\rm{s.t.}} \quad a_i^2 = 1, \quad i = 1, ..., n, $$ where $G$ is real symmetric matrix. $\mathbb{OP1}$ is a non-homogeneous Boolean Quadratic Program (BQP). To transform $\mathbb{OP1}$ to a homogeneous BQP, let us introduce a scalar $t \in \mathbb{R}$ and the new optimization problem is written as $$\mathbb{OP2}: \min_{\mathbf{a} \in \mathbb{R}^n, \: t \in \mathbb{R}} \quad \quad \left[ \mathbf{a}^{\rm{T}} t \right ] \begin{bmatrix} G & \mathbf{y}\\ \mathbf{y}^{\rm{T}} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{a}\\ t \end{bmatrix} \\ {\rm{s.t.}} \quad a_i^2 = 1, \quad i = 1, ..., n, \\ t^2 = 1 \hspace{2cm}$$ So, the decision variables and the number of constraints are increased by 1. Let us define $x = \begin{bmatrix} \mathbf{a}\\ t \end{bmatrix}$ and $H = \begin{bmatrix} G & \mathbf{y}\\ \mathbf{y}^{\rm{T}} & 0 \end{bmatrix}$, we can write the following optimization problem $$\mathbb{OP3}: \min_{\mathbf{x} \in \mathbb{R}^{n+1}} \quad \quad \mathbf{x}^{\rm{T}} H \mathbf{x}\\ {\rm{s.t.}} \quad x_i^2 = 1, \quad i = 1, ..., n+1, $$ If we define $X = \mathbf{x} \mathbf{x}^{\rm{T}}$, which means that $X$ is a rank one symmetric positive semi-definite matrix, we can write the following equivalent optimization problem $$\mathbb{OP4}: \min_{X \in \mathbb{S}^{n+1}} \quad \quad {\rm{Trace}}\{H X\}\\ {\rm{s.t.}} \quad X_{ii} = 1, \quad i = 1, ..., n+1, \\ X \succcurlyeq 0, \\ {\rm{Rank}}\{X\} = 1 $$ By relaxing the ${\rm{Rank}}$ 1 constraint, $\mathbb{OP4}$ can be solved by SDR (I used CVX and Matlab). After getting the optimal solution $X^*$, we need to find a feasible solution $\mathbf{x}$ from $X^*$. It seems that one of the efficient methods is to use Gaussian randomization and quantization to generate a number of feasible point and pick the one with the least objective value. For example, we generate $\hat{x} \sim {\rm{sgn}} \{\mathcal{N}(0,X^*)\}$.

The matlab code that solves $\mathbb{OP4}$ is given below

clear all
clc
L = 100;    % number of trials to generate random vector
n = 5;
H = randn(n,n); 
cvx_begin
variable X(n,n) symmetric
minimize(trace(H*X));
subject to
diag(X) == 1;
X == semidefinite(n);
cvx_end
for i = 1:L
    X_hat = X.*randn(1,1);
    eta = sign(X_hat(:,3)); % picked one of the columns
    x_hat(:,i) = eta;
    myrank(i) = rank(eta*eta.');
    objval(i) = trace(eta.'*H*eta);
    constranit_check(i) = sum(eta.^2)/n;
end
[val,loc] = min(objval);
val

My questions are are follows:

1. Are my transformation from $\mathbb{OP1}$ to $\mathbb{OP4}$ correct? or I missed up some of the QCQP and SDR concepts?

2. I am not sure if my understanding of the randomization part is clear as the generated random variables will be a matrix and I need only a vector. That is why I picked column 3 (as an example) of the generated random matrix.

3. There is something wrong in the code as for the 100 realizations I got the same objective function value.

Sorry for the long post. Any help will be appreciated.

Answer 1

The transformation is correct. But the code is wrong. You misunderstand the randomization part. The random variable should be a vector, not a scalar.

I use G-W rounding to obtain a rank-one matrix. I am not sure which method you has used.

Hope my code is helpful to your work.

clear all
clc
numTest = 10; % number of trials to generate random vector 
n = 5;
H = randn(n,n); 

%%
cvx_begin 
variable X(n,n) symmetric semidefinite 
minimize(trace(H*X)); 
subject to diag(X) == 1; 
cvx_end

%% 
ther = 1e-4; 
[U,R] = schur(X); 
eig_value = diag(R); 
rank_x = sum(eig_value>1e-4);

V = sqrt(diag(eig_value(eig_value>ther)))*U(eig_value>ther,:);

for nTest = 1:numTest
 x_hat = sign(diag(randn(n,rank_x)V));
 result(:,nTest) =x_hat;
 objval(nTest) = trace(Hx_hat*x_hat');
 constranit_check(nTest) = mean(x_hat.^2);
end
[val,loc] = min(objval);
val