Find maximum of quadratically constrained quadratic problem using MATLAB

Hey,

I am trying to solve for maximized SNR at receiver of a two hop network. First hop channel coefficients are denoted by N*N matrix g and second hop channel coefficients are denoted by N*N matrix h. We have an N size array of phase tuning varibales using which we have to maximize frobenious norm of h*diag(exp(1j.*phi))*g.

with constraints that each diagnoal element should be unit modulas.

I have written the code as below. MATLAB wasn't allowing to combine complex data with optimization variable so I separted the real and imaginary part.

N=4;
h=sqrt(0.5)*(randn(N,N)+1j*randn(N,N));
h1=real(h);h2=imag(h);
g=sqrt(0.5)*(randn(N,N)+1j*randn(N,N));
g1=real(g);g2=imag(g);

k=optimproblem;
p=optimvar('p',N);
phi1=diag(cos(p));
phi2=diag(sin(p));
k.ObjectiveSense='maximize';
k.Objective=sum(h1*phi1*g1-h2*phi2*g1-h1*phi2*g2-h2*phi1*g2,'all')^2 +sum(h1*phi1*g2-h2*phi2*g2-h1*phi2*g1+h2*phi1*g1,'all')^2;
k.Constraints.intercept1 = (0<=p);
k.Constraints.intercept2 = (p<=2*pi);
sol0.p=zeros(1,N);
sol=solve(k,sol0);

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

cNew=sol.p;

It seems to be a non convex QCQP, can someone please help in solving this. Above code is using fmincon but it is not giving correct results.

[Matt J]

Above code is using fmincon but it is not giving correct results.

As indicated by what?

[Shankul Saini]

Thanks @Matt J for replying, I know the results for simple case. If you take h a 1*N vector and g a N*1 vector. We can trivially see what should be the value of phi, it should be cancelling the phase of corresponding g transpose and h, i.e. it should be sum of phase of h and g transpose. But even for that case results using fmincon aren't matching that's why i wrote the that results aren't matching.

[Shankul Saini]

yeah it does seem like it. Just one this why you have taken @phi * minus * there, because fmincon will minimize and we need to maximize that's why.. and the other thing is, this is also giving incorrect results. If you simply take h=[1,-2] and g=[1;2] then answer for phi should be [0,pi] but it is giving something else only.

Like the maximum value of obj function should be 25 but it is giving 9.

[Matt J]

The expression h*diag(exp(1j.*phi))*g is not scalar-valued if h and g are both NxN matrices. So, it is not clear what it even means to "maximize" it.

[Shankul Saini]

yeah @Matt J you are right, but I have to maximize the frobenious norm of this matrix.

i.e. trace((h*diag(exp(1j.*phi))*g)*(h*diag(exp(1j.*phi))*g)')

[Bruno Luong]

Your original code with the objective function expressed with real/imag of phi doesn't seem to return the frobenious norm

N=4;
h=sqrt(0.5)*(randn(N,N)+1j*randn(N,N));
h1=real(h);h2=imag(h);
g=sqrt(0.5)*(randn(N,N)+1j*randn(N,N));
g1=real(g);g2=imag(g);
p=2*pi*rand(N,1);
phi1=diag(cos(p));
phi2=diag(sin(p));

sum(h1*phi1*g1-h2*phi2*g1-h1*phi2*g2-h2*phi1*g2,'all')^2 +sum(h1*phi1*g2-h2*phi2*g2-h1*phi2*g1+h2*phi1*g1,'all')^2

ans = 0.9829

M = h*diag(exp(1i*p))*g;
norm(M,'fro')^2

ans = 67.4239

trace(M*M')

ans = 67.4239

[Shankul Saini]

yeah correctly pointed out @Bruno Luong, let me look into this why this isn't the same and where i was doing mistake. But thank you very much for helping.

sum((h1*phi1*g1-h2*phi2*g1-h1*phi2*g2-h2*phi1*g2).^2,"all") + sum((h1*phi1*g2-h2*phi2*g2+h1*phi2*g1+h2*phi1*g1).^2,'all')

Unrecognized function or variable 'h1'.

this was the line supposed to be there @Bruno Luong, thanks again for pointing out.

Answer 1

[Torsten]

Is this the simplified version of the problem you are trying to solve ?

It seems fmincon gives the correct solution.

N = 4;
h = sqrt(0.5)*(randn(1,N)+1i*randn(1,N));
g = sqrt(0.5)*(randn(N,1)+1i*randn(N,1));
phi0 = zeros(N,1);
lb = zeros(N,1);
ub = 2*pi*ones(N,1);
obj = @(phi)-(h*diag(exp(1i*phi))*g)*(h*diag(exp(1i*phi))*g)';
[sol,fval] = fmincon(obj,phi0,[],[],[],[],lb,ub)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = 4×1

2.2056 2.7240 4.1670 3.1435

fval = -1.3543

phi1 = -(phase(h)+phase(g.'));
obj(phi1)

ans = -1.3543

[Shankul Saini]

ok yeah @Torsten for the above programme it is working absolutely fine, but can you please take value of N=2, h as [1,-2] and g=[1;2]. Then can you please tell why is it not working.

For simplicity since many days i was working with latter case and it is not giving correct results.

and again thank you very much @Torsten for looking into it.

[Torsten]

The reason you don't get the global maximum is that your problem is not convex.

So you can get stuck in local optima if the initial guess for phi is bad.

[Shankul Saini]

so does this mean, this solution isn't trustworthy is i m taking it in loop running 100's of times to get optimized phi every time?

cause we can't change initial guess every time..

[Torsten]

You could try Bruno's code or MATLAB's MultiStart.

I don't know if a different MATLAB solver could overcome this problem.

I'm quite surprised that fmincon stops at a local maximum instead of a local minimum (the objective function for your case from above is -17 + 8*cos(phi1-phi2) which has a local maximum for phi1 = phi2 instead of a local minimum).

N = 2;
h = [1,-2] ;
g = [1;2] ;
phi0 = zeros(N,1);
lb = zeros(N,1);
ub = 2*pi*ones(N,1);
obj = @(phi)-(h*diag(exp(1i*phi))*g)*(h*diag(exp(1i*phi))*g)';
[sol,fval] = fmincon(obj,phi0,[],[],[],[],lb,ub)

Initial point is a local minimum that satisfies the constraints.

Optimization completed because at the initial point, the objective function is non-decreasing 
in feasible directions to within the value of the optimality tolerance, and 
constraints are satisfied to within the value of the constraint tolerance.

sol = 2×1

0.9900 0.9900

fval = -9

[Shankul Saini]

Okay yeah @torsten i will try one of those, but another thing is how to solve the original problem where h and g were matrices. How would objective function look like in that case ? And regrading @bruno code, does that get applies on my scenrio i am not very sure about that?

[Torsten]

how to solve the original problem where h and g were matrices.

You should know better than I do what your objective is in this case. Isn't it

obj = @(phi) -norm(h*diag(exp(1i*phi))*g,'fro')^2;

or equivalently

obj = @(phi) -trace((h*diag(exp(1i*phi))*g)*(h*diag(exp(1i*phi))*g)');

?

You wrote above about maximizing h*diag(exp(1j.*phi))*g. But if h and g are matrices, it's difficult to tell what you mean by "maximizing a matrix".

[Shankul Saini]

yeah correct I was wrong there, I should had written frobenious norm of h*diag(exp(1j.*phi))*g.

But thanks for helping @Torsten.

[Shankul Saini]

clear all
clc
N = 4;
h = sqrt(0.5)*(randn(2,N)+1i*randn(2,N));
g = sqrt(0.5)*(randn(N,2)+1i*randn(N,2));
phi0 = zeros(N,1);
lb = zeros(N,1);
ub = 2*pi*ones(N,1);
obj = @(phi) -trace((h*diag(exp(1i*phi))*g)*(h*diag(exp(1i*phi))*g)');
[sol,fval] = fmincon(obj,phi0,[],[],[],[],lb,ub)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = 4×1

4.1887 3.8264 1.9413 2.5359

fval = -10.0484

Hey @Torsten, when i m changing the size of h and g as my need my matlab is giving error. Then i pasted the code here and surprisingly it ran successfully here. Also when i tried code in matlab online simulator it ran there also. Can you please tell the problem which is occuring in my pc matlab by seeing the following error.

[Bruno Luong]

Try to change to

obj = @(phi)  norm( h*diag(exp(1i*phi))*g,'fro')^2;

Might be the Trace expression returns some small non-zero imaginary part in the output.

I would also remove the constraint 0 <= phi = 2*pi, and let the minimizer do the unconstrainted and take the modulo 2*pi after getting the result. The constraints just make life difficult for fmincon for no apparent reason.

[Shankul Saini]

ok thanks, using norm solved the error. and also yeah taking modulo instead of constraints will ease out the code correct. Thank you.

[Shankul Saini]

N = 1;
h = sqrt(0.5)*(randn(1,N)+1i*randn(1,N));
g = sqrt(0.5)*(randn(N,1)+1i*randn(N,1));
phi0 = zeros(N,1);
lb = zeros(N,1);
ub = 2*pi*ones(N,1);
obj = @(phi)-((h*diag(exp(1i*phi))*g));
[sol,fval] = fmincon(obj,phi0,[],[],[],[],lb,ub)

Error using barrier
Objective function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 886)
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

phi1 = -(angle(h)+angle(g.'));

If in the above code i am trying to maximize ((h*diag(exp(1i*phi))*g)) where all h,g are complex scalar. Why matlab is giving such error, is it because objective function is complex, or some other reason?

[Bruno Luong]

The complex has not ordering, so MATLAB throw the error telling that what you ask to minimize has no sense.

Any reason why you not longer want norm(...'fro') as you told earlier?

[Shankul Saini]

yeah @Bruno Luong, because I wanted to see if the point "phi1" is equal to the optimization point "sol", although objective function is giving same value on both points in every case.

So I tried for simplest case, but when i m taking norm it is getting satisfied at inital point only because it just takes the magnitude part and makes angle part 1 so but phi1 is giving different value.

N = 1;
h = sqrt(0.5)*(randn(1,N)+1i*randn(1,N));
g = sqrt(0.5)*(randn(N,1)+1i*randn(N,1));
phi0 = zeros(N,1);
lb = zeros(N,1);
ub = 2*pi*ones(N,1);
obj = @(phi)-norm((h*diag(exp(1i*phi))*g));
[sol,fval] = fmincon(obj,phi0,[],[],[],[],lb,ub)

Initial point is a local minimum that satisfies the constraints.

Optimization completed because at the initial point, the objective function is non-decreasing 
in feasible directions to within the value of the optimality tolerance, and 
constraints are satisfied to within the value of the constraint tolerance.

sol = 0.9900

fval = -0.3117

phi1 = -(angle(h)+angle(g.'))

phi1 = -2.1901

obj(phi1)

ans = -0.3117

So now here is you see, value of obj(phi1) and fval are same but this is not the case with sol and phi1.

and thank you for pointing out that complex numbers have no ordering, now can you please tell me some way to workout my problem..

[Bruno Luong]

The issue is NOT fmincon but it seems the formulation of the simplified problem, which is constant (=abs(h*g)) regardless the value of phi on the admissible set, so it does not give what you think is good.

The formulation and expectation, and intuition of what is right is on your side, not MATLAB so noone can help you. You told us that the fro norm is the SNR of your system. Obviously either that is wrong or either phi has no a right control of your SNR. Only you can make sense of it (or not).

[Shankul Saini]

Yeah @bruno luong, you are correct and thank you very much for the patient explanation.

Answer 2

[Bruno Luong]

The problem of least-square minimization under quadratic constraints is known to might have many local minima. So if you use fmincon or such without care on first guess initialization it could converge to local minima.

You could use my FEX https://fr.mathworks.com/matlabcentral/fileexchange/27596-least-square-with-2-norm-constraint?s_tid=srchtitle

where one of the method use reformulation as quadratic-eigenvalue could overcome such numerical issue.

[Shankul Saini]

yeah your solution could help but in your FEX constraint is on norm of vector but in my part it is on every value of the vector.

[Bruno Luong]

OK then you do not have quadratic constrained problem as the subject missleadingly suggests

https://en.wikipedia.org/wiki/Quadratically_constrained_quadratic_program

[Shankul Saini]

@Bruno Luong yeah definietly i may be wrong, i do not posses much background in optimization theory but I saw from this research paper that my problem is also a non convex QCQP .

[Torsten]

You've put the constraints in the objective function.

This way, you get an unconstraint optimization problem with general objective.