Different coding style results different results

[youhua wang]

Hello Everybody,

I want to solve a SDP problem with YALMIP or CVX. I have two versions of MATLAB code which I thought are the same in mathematical viewpoint. But the results of those two codes are different.

Could you please be kind to help me find the reason?

The problem I want to solve can be expressed as

Where Teop(u) is defined as

The first CVX code I wrote is as follows:

%%%%%%%%%%%%%%%%

n = length(y);

cvx_begin

variable u0;

variable u(n-1) complex;

variable Tu(n,n) hermitian;

variable t;

variable xhat(n) complex;

minimize (tau*(n*u0 + t)/2 + 1/2*pow_pos(norm(y - xhat,'fro'),2))

subject to

[toeplitz([u0;conj(u)],[u0;conj(u)]') xhat; xhat' t] == semidefinite(n+1,n+1);

cvx_end

%%%%%%%%%%%%%%%%%%%%

The secod CVX code I wrote is as follows:

%%%%%%%%%%%%%%%%%%%

n = length(y);

cvx_begin

variable u0;

variable u(n-1) complex;

variable Tu(n,n) hermitian;

variable t;

variable xhat(n) complex;

minimize (tau*(n*u0 + t)/2 + 1/2*pow_pos(norm(y - xhat,'fro'),2))

subject to

Tu == toeplitz([u0;conj(u)],[u0;conj(u)]');

[Tu xhat; xhat' t] == semidefinite(n+1,n+1);

cvx_end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I think they are the same. But the results they produced are different largely.

Very similarly, I write two versions of YALMIP code. The results of those two codes are also different.

The first YALMIP code is as follows:

%%%%%%%%%%%%%%%%%%%%%%%%

sdpvar u1;

u2=sdpvar(N-1,1,'full','complex');

u=[u1;u2];

x=sdpvar(N,1,'full','complex');

sdpvar t;

Tu=sdpvar(N,N,'toeplitz','complex');

options=sdpsettings('verbose',1);

 

Constraints=[Tu==toeplitz(u)];

Constraints=[Constraints,[Tu x;x' t]>0];

 

Objective=tau/2*(t+trace(toeplitz(u)))+1/2*norm(y_observerd-x(index),2)^2;

sol=solvesdp(Constraints,Objective,options);

if(sol.problem==0)

    solution=double(x);

else

    disp('Something went wrong!')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The second YALMIP code is as follows:

%%%%%%%%%%%%%%%%%%%

sdpvar u1;

u2=sdpvar(N-1,1,'full','complex');

u=[u1;u2];

x=sdpvar(N,1,'full','complex');

sdpvar t;

Tu=sdpvar(N,N,'toeplitz','complex');

options=sdpsettings('verbose',1);

 

Constraints=[[toeplitz(u) x;x' t]>0];

Objective=tau/2*(t+trace(toeplitz(u)))+1/2*norm(y_observerd-x(index),2)^2;

sol=solvesdp(Constraints,Objective,options);

if(sol.problem==0)

    solution=double(x);

else

    disp('Something went wrong!')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I also attached the complete MATLAB scripts that generate the results as follows:

%%%%%%%%%%%%%%%%

clc

clear

 

%%  simulation setup

f=[0, 0.06,0.15,0.3 0.2]; % normalized frequency

theta=[0, 0,0.1*pi,pi 1.3*pi]; % phase of the individual components

Ampli=[1, 1,0.8,1.2 1.0]; % amplitude of the frequency components

N=32; % number of equal time samples

m=20; % number of undersampling

 

%% measurement processing

M=zeros(N,length(f)); % manifold matrix

for(i=1:length(f))

    for(j=1:N)

        M(j,i)=exp(sqrt(-1)*theta(i))*exp(sqrt(-1)*2*pi*f(i)*(j-1));

    end   

end

y=M*Ampli'; % noiseless measurement

 

 

% We add noise so that the average signal to noise ratio is 20 dB

SNR = 20;

noise_std = norm(y)/sqrt(N)*10^(-SNR/20);

observed = y + noise_std*(randn(N,1) + sqrt(-1)*randn(N,1))/sqrt(2);

 

index=randperm(N);

index=index(1:m);

 

y_observerd=zeros(m,1);

y_observerd=observed(index);

tau=2;

 

%% recover the original data with SDP method

sdpvar u1;

u2=sdpvar(N-1,1,'full','complex');

u=[u1;u2];

x=sdpvar(N,1,'full','complex');

sdpvar t;

Tu=sdpvar(N,N,'toeplitz','complex');

options=sdpsettings('verbose',1);

 

 

% Constraints=[Tu==toeplitz(u)];   % first version of YALMIP code

% Constraints=[Constraints,[Tu x;x' t]>0]; % first version of YALMIP code

 

 

Constraints=[[toeplitz(u) x;x' t]>0]; % seconde version of YALMIP code

Objective=tau/2*(t+trace(toeplitz(u)))+1/2*norm(y_observerd-x(index),2)^2;

sol=solvesdp(Constraints,Objective,options);

if(sol.problem==0)

    solution=double(x);

else

    disp('Something went wrong!')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Another question is how can I get the solution of the dual problem?

[Johan Löfberg]

This

sdpvar(N,N,'toeplitz','complex')
Linear matrix variable 32x32 (full, complex, 64 variables)

does not generate an Hermitian matrix, hence

[Tu x;x' t]>0

defines linear elementwise constraints and not a semidefinite constraint. Not what you want I guess.

To get duals, you use the command dual

[Johan Löfberg]

The problem is that the diagonal of the matrix is allowed to have imaganiry terms, which means it isn't structurally Hermitian. Hence, you can remove these terms and obtain a structurally Hermitian complex Toeplitz, which you then can use directly in a semidefinite constraint

Tu = Tu-sqrt(-1)*diag(imag(diag(Tu)))

[youhua wang]

Thanks Johan,

So the following constraints 

%%%%%%%%%%%%%%%%%%%%%%%%%%

sdpvar u1;

u2=sdpvar(N-1,1,'full','complex');

u=[u1;u2];

x=sdpvar(N,1,'full','complex');

sdpvar t;

Tu=sdpvar(N,N,'toeplitz','complex');

Tu = Tu-sqrt(-1)*diag(imag(diag(Tu)));

Constraints=[Tu==toeplitz(u)];

Constraints=[Constraints,[Tu x;x' t]>0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

is the same as these

%%%%%%%%%%%%%%%%%%%%%

sdpvar u1;

u2=sdpvar(N-1,1,'full','complex');

u=[u1;u2];

x=sdpvar(N,1,'full','complex');

sdpvar t;

Constraints=[[toeplitz(u) x;x' t]>0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Is my understanding right?

Could you please give a more detail about the usage of command dual? Since I use dual(Constraints), a error is occurred which says Dual not applicable on list of constraints.  

[Johan Löfberg]

Depends on what you mean with equivalent. The first model has many more variables thus harder to solve and is a really weird way to describe the second model.

dual can only be applied to a particular constraint
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.Dual

[youhua wang]

Got it! Thanks Johan